From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Wed, 21 Mar 2018 00:45:24 +0000 (-0600)
Subject: Fix iterating over entries of parallel PETSc matrices.
X-Git-Tag: v9.0.0-rc1~271^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c42af443d5791553e66812bbab05fd818bf3e895;p=dealii.git

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 2be755ce3c..6ab9d445f6 100644
--- 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 65bc5e9107..80395aaaa8 100644
--- 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 ();