From f9689e8cd7c52c6623fda7cf75a27f3001013c35 Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 7 Jan 2013 14:56:43 +0000
Subject: [PATCH] Fix implementation of SparseMatrix::iterator::operator- for
 the general case. Add a test.

git-svn-id: https://svn.dealii.org/trunk@27961 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/lac/sparse_matrix.h   |  14 +--
 tests/bits/sparse_matrix_iterator_13.cc       | 107 ++++++++++++++++++
 .../sparse_matrix_iterator_13/cmp/generic     |   2 +
 3 files changed, 116 insertions(+), 7 deletions(-)
 create mode 100644 tests/bits/sparse_matrix_iterator_13.cc
 create mode 100644 tests/bits/sparse_matrix_iterator_13/cmp/generic

diff --git a/deal.II/include/deal.II/lac/sparse_matrix.h b/deal.II/include/deal.II/lac/sparse_matrix.h
index f787951f9e..96b39d28ff 100644
--- a/deal.II/include/deal.II/lac/sparse_matrix.h
+++ b/deal.II/include/deal.II/lac/sparse_matrix.h
@@ -2240,13 +2240,13 @@ namespace SparseMatrixIterators
       return ((*this)->index() - other->index());
     else
       {
-//TODO: this shouldn't be so hard to implement. it could either be done as
-// std::difference(*this, other), but for that we lack a bunch of typedefs
-// in the iterator class; it is also inefficient since it has linear complexity
-// in the distance. alternatively, one should be able to just add up the
-// entries in all of the rows of the matrix between *this and other
-        Assert (false, ExcNotImplemented());
-        return 0;
+        const SparsityPattern &sparsity = accessor.get_matrix().get_sparsity_pattern();
+        const unsigned int this_position
+        = sparsity.get_rowstart_indices()[(*this)->row()] + (*this)->index();
+        const unsigned int other_position
+        = sparsity.get_rowstart_indices()[other->row()] + other->index();
+
+        return (this_position - other_position);
       }
   }
 
diff --git a/tests/bits/sparse_matrix_iterator_13.cc b/tests/bits/sparse_matrix_iterator_13.cc
new file mode 100644
index 0000000000..94661d23ba
--- /dev/null
+++ b/tests/bits/sparse_matrix_iterator_13.cc
@@ -0,0 +1,107 @@
+//----------------------------  sparse_matrix_iterator_13.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparse_matrix_iterator_13.cc  ---------------------------
+
+// test SparseMatrix::iterator::operator-
+
+#include "../tests.h"
+#include <deal.II/lac/sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test ()
+{
+  SparsityPattern sp (5,5,3);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        sp.add (i,j);
+  sp.compress ();
+
+  SparseMatrix<double> m(sp);
+
+  for (unsigned int row=0; row<sp.n_rows(); ++row)
+    Assert (m.begin(row)-m.begin(row) == 0,
+	    ExcInternalError());
+
+  for (unsigned int row=0; row<sp.n_rows(); ++row)
+    Assert (m.end(row)-m.begin(row) == (int)sp.row_length(row),
+	    ExcInternalError());
+  for (unsigned int row=0; row<sp.n_rows(); ++row)
+    Assert (m.begin(row)-m.end(row) == -(int)sp.row_length(row),
+	    ExcInternalError());
+
+  {
+    unsigned int counter = 0;
+    for (unsigned int row=0; row<sp.n_rows(); ++row)
+      {
+	Assert (m.begin(row)-m.begin(0) == (int)counter,
+		ExcInternalError());
+	Assert (m.begin(0)-m.begin(row) == -(int)counter,
+		ExcInternalError());
+	counter += sp.row_length(row);
+      }
+  }
+
+  Assert (m.begin() - m.begin(0) == 0, ExcInternalError());
+  Assert (m.begin(0) - m.begin() == 0, ExcInternalError());
+  Assert (m.end(sp.n_rows()-1) - m.end() == 0, ExcInternalError());
+  Assert (m.end() - m.end(sp.n_rows()-1) == 0, ExcInternalError());
+  Assert (m.end() - m.begin() == (int)sp.n_nonzero_elements(),
+	  ExcInternalError());
+  Assert (m.begin() - m.end() == -(int)sp.n_nonzero_elements(),
+	  ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("sparse_matrix_iterator_13/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test ();
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      deallog << "Exception on processing: " << std::endl
+		<< exc.what() << std::endl
+		<< "Aborting!" << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      deallog << "Unknown exception!" << std::endl
+		<< "Aborting!" << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      return 1;
+    };
+}
diff --git a/tests/bits/sparse_matrix_iterator_13/cmp/generic b/tests/bits/sparse_matrix_iterator_13/cmp/generic
new file mode 100644
index 0000000000..0fd8fc12f0
--- /dev/null
+++ b/tests/bits/sparse_matrix_iterator_13/cmp/generic
@@ -0,0 +1,2 @@
+
+DEAL::OK
-- 
2.39.5