From c42af443d5791553e66812bbab05fd818bf3e895 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 20 Mar 2018 18:45:24 -0600 Subject: [PATCH] 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. --- include/deal.II/lac/petsc_matrix_base.h | 31 +++++++++++++++++-------- source/lac/petsc_matrix_base.cc | 9 ++++--- 2 files changed, 25 insertions(+), 15 deletions(-) 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 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 (); -- 2.39.5