Fix iterating over entries of parallel PETSc matrices.

Previously, it was not possible to iterate over the local range because
one would have to call
  matrix.end(row)
where 'row' is the last locally owned row, and that triggered an
assertion because this end iterator is also the begin iterator of
the next row -- which is not locally owned any more.

Fix this.

diff --git a/include/deal.II/lac/petsc_matrix_base.h b/include/deal.II/lac/petsc_matrix_base.h
index 2be755ce3c3116f9225f19251d4054f538798ab6..6ab9d445f6be26dab8e2fd4c5df7d4c16b124d9c 100644 (file)
--- a/include/deal.II/lac/petsc_matrix_base.h
+++ b/include/deal.II/lac/petsc_matrix_base.h
@@ -1111,15 +1111,16 @@ namespace PETScWrappers
 
       ++accessor.a_index;
 
-      // if at end of line: do one step, then
-      // cycle until we find a row with a
-      // nonzero number of entries
+      // if at end of line: do one step, then cycle until we find a
+      // row with a nonzero number of entries
       if (accessor.a_index >= accessor.colnum_cache->size())
         {
           accessor.a_index = 0;
           ++accessor.a_row;
 
           while ((accessor.a_row < accessor.matrix->m())
+                 &&
+                 (accessor.a_row < accessor.matrix->local_range().second)
                  &&
                  (accessor.matrix->row_length(accessor.a_row) == 0))
             ++accessor.a_row;
@@ -1492,7 +1493,9 @@ namespace PETScWrappers
   MatrixBase::const_iterator
   MatrixBase::begin(const size_type r) const
   {
-    Assert (r < m(), ExcIndexRange(r, 0, m()));
+    Assert (in_local_range(r),
+            ExcIndexRange(r, local_range().first, local_range().second));
+
     if (row_length(r) > 0)
       return const_iterator(this, r, 0);
     else
@@ -1504,13 +1507,21 @@ namespace PETScWrappers
   MatrixBase::const_iterator
   MatrixBase::end(const size_type r) const
   {
-    Assert (r < m(), ExcIndexRange(r, 0, m()));
-
-    // place the iterator on the first entry
-    // past this line, or at the end of the
-    // matrix
+    Assert (in_local_range(r),
+            ExcIndexRange(r, local_range().first, local_range().second));
+
+    // place the iterator on the first entry past this line, or at the
+    // end of the matrix
+    //
+    // in the parallel case, we need to put it on the first entry of
+    // the first row after the locally owned range. this of course
+    // doesn't exist, but we can nevertheless create such an
+    // iterator. we need to check whether 'i' is past the locally
+    // owned range of rows first, before we ask for the length of the
+    // row since the latter query leads to an exception in case we ask
+    // for a row that is not locally owned
     for (size_type i=r+1; i<m(); ++i)
-      if (row_length(i) > 0)
+      if (i == local_range().second || (row_length(i) > 0))
         return const_iterator(this, i, 0);
 
     // if there is no such line, then take the

diff --git a/source/lac/petsc_matrix_base.cc b/source/lac/petsc_matrix_base.cc
index 65bc5e910768fb8ae221d0db9ecf6bde2ba87f31..80395aaaa88e220da8b1c8b7bad535a57faff053 100644 (file)
--- a/source/lac/petsc_matrix_base.cc
+++ b/source/lac/petsc_matrix_base.cc
@@ -34,11 +34,10 @@ namespace PETScWrappers
     MatrixBase::const_iterator::Accessor::
     visit_present_row ()
     {
-      // if we are asked to visit the
-      // past-the-end line, then simply
-      // release all our caches and go on
-      // with life
-      if (this->a_row == matrix->m())
+      // if we are asked to visit the past-the-end line (or a line that is not
+      // stored on the current processor), then simply release all our caches
+      // and go on with life
+      if (matrix->in_local_range(this->a_row) == false)
         {
           colnum_cache.reset ();
           value_cache.reset ();