From: wolf Date: Tue, 3 May 2005 03:02:23 +0000 (+0000) Subject: Add row_length(). Make sure that matrix iterators also work if there are empty rows... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=64d18875075e9bba932bc444d444700669d876c9;p=dealii-svn.git Add row_length(). Make sure that matrix iterators also work if there are empty rows in the matrices. git-svn-id: https://svn.dealii.org/trunk@10603 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/petsc_matrix_base.h b/deal.II/lac/include/lac/petsc_matrix_base.h index df58d602c7..ace5bf7eff 100644 --- a/deal.II/lac/include/lac/petsc_matrix_base.h +++ b/deal.II/lac/include/lac/petsc_matrix_base.h @@ -491,6 +491,12 @@ namespace PETScWrappers */ unsigned int n_nonzero_elements () const; + + /** + * Number of entries in a specific row. + */ + unsigned int row_length (const unsigned int row) const; + /** * Return the l1-norm of the matrix, that is * $|M|_1=max_{all columns j}\sum_{all @@ -669,12 +675,27 @@ namespace PETScWrappers /** * STL-like iterator with the * first entry of row @p r. + * + * Note that if the given row is empty, + * i.e. does not contain any nonzero + * entries, then the iterator returned by + * this function equals + * end(r). Note also that the + * iterator may not be dereferencable in + * that case. */ const_iterator begin (const unsigned int r) const; /** - * Final iterator of row - * @p r. + * Final iterator of row r. It + * points to the first element past the + * end of line @p r, or past the end of + * the entire sparsity pattern. + * + * Note that the end iterator is not + * necessarily dereferencable. This is in + * particular the case if it is the end + * iterator for the last row of a matrix. */ const_iterator end (const unsigned int r) const; @@ -764,9 +785,6 @@ namespace PETScWrappers a_index(index) { visit_present_row (); - Assert ((row == matrix->m()) || - (index < colnum_cache->size()), - ExcInvalidIndexWithinRow (index, row)); } @@ -824,10 +842,25 @@ namespace PETScWrappers Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd()); ++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 (accessor.a_index >= accessor.colnum_cache->size()) { accessor.a_index = 0; - accessor.a_row++; + ++accessor.a_row; + + while (accessor.a_index >= accessor.matrix->row_length(accessor.a_row)) + { + ++accessor.a_row; + + // if we happened to find the end + // of the matrix, then stop here + if (accessor.a_row == accessor.matrix->m()) + break; + } + accessor.visit_present_row(); } return *this; @@ -839,16 +872,7 @@ namespace PETScWrappers const_iterator::operator++ (int) { const const_iterator old_state = *this; - - Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd()); - - ++accessor.a_index; - if (accessor.a_index >= accessor.colnum_cache->size()) - { - accessor.a_index = 0; - accessor.a_row++; - accessor.visit_present_row(); - } + ++(*this); return old_state; } @@ -932,7 +956,10 @@ namespace PETScWrappers MatrixBase::begin(const unsigned int r) const { Assert (r < m(), ExcIndexRange(r, 0, m())); - return const_iterator(this, r, 0); + if (row_length(r) > 0) + return const_iterator(this, r, 0); + else + return end (r); } @@ -941,7 +968,17 @@ namespace PETScWrappers MatrixBase::end(const unsigned int r) const { Assert (r < m(), ExcIndexRange(r, 0, m())); - return const_iterator(this, r+1, 0); + + // place the iterator on the first entry + // past this line, or at the end of the + // matrix + for (unsigned int i=r+1; i 0) + return const_iterator(this, i, 0); + + // if there is no such line, then take the + // end iterator of the matrix + return end(); } diff --git a/deal.II/lac/source/petsc_matrix_base.cc b/deal.II/lac/source/petsc_matrix_base.cc index 7e5f0c9c03..e09f59fd22 100644 --- a/deal.II/lac/source/petsc_matrix_base.cc +++ b/deal.II/lac/source/petsc_matrix_base.cc @@ -63,18 +63,15 @@ namespace PETScWrappers AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr)); // copy it into our caches if the line - // isn't empty - if (ncols != 0) - { - colnum_cache.reset (new std::vector (colnums, - colnums+ncols)); - value_cache.reset (new std::vector (values, values+ncols)); - } - else - { - colnum_cache.reset (); - value_cache.reset (); - } + // isn't empty. if it is, then we've + // done something wrong, since we + // shouldn't have initialized an + // iterator for an empty line (what + // would it point to?) + Assert (ncols != 0, ExcInternalError()); + colnum_cache.reset (new std::vector (colnums, + colnums+ncols)); + value_cache.reset (new std::vector (values, values+ncols)); // and finally restore the matrix ierr = MatRestoreRow(*matrix, this->a_row, &ncols, &colnums, &values); @@ -309,6 +306,48 @@ namespace PETScWrappers + unsigned int + MatrixBase:: + row_length (const unsigned int row) const + { +//TODO: this function will probably only work if compress() was called on the +//matrix previously. however, we can't do this here, since it would impose +//global communication and one would have to make sure that this function is +//called the same number of times from all processors, something that is +//unreasonable. there should simply be a way in PETSc to query the number of +//entries in a row bypassing the call to compress(), but I can't find one + Assert (row < m(), ExcInternalError()); + + // get a representation of the present + // row + int ncols; + +#if (PETSC_VERSION_MAJOR <= 2) && \ + ((PETSC_VERSION_MINOR < 2) || \ + ((PETSC_VERSION_MINOR == 2) && (PETSC_VERSION_SUBMINOR == 0))) + int *colnums; + PetscScalar *values; +#else + const int *colnums; + const PetscScalar *values; +#endif + +//TODO: this is probably horribly inefficient; we should lobby for a way to +//query this information from PETSc + int ierr; + ierr = MatGetRow(*this, row, &ncols, &colnums, &values); + Assert (ierr == 0, MatrixBase::ExcPETScError(ierr)); + + // then restore the matrix and return the + // number of columns in this row as + // queried previously + ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values); + AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr)); + + return ncols; + } + + PetscReal MatrixBase::l1_norm () const {