*/
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
/**
* 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
+ * <tt>end(r)</tt>. 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 <tt>r</tt>. 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;
a_index(index)
{
visit_present_row ();
- Assert ((row == matrix->m()) ||
- (index < colnum_cache->size()),
- ExcInvalidIndexWithinRow (index, row));
}
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;
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;
}
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);
}
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<m(); ++i)
+ if (row_length(i) > 0)
+ return const_iterator(this, i, 0);
+
+ // if there is no such line, then take the
+ // end iterator of the matrix
+ return end();
}
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<unsigned int> (colnums,
- colnums+ncols));
- value_cache.reset (new std::vector<PetscScalar> (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<unsigned int> (colnums,
+ colnums+ncols));
+ value_cache.reset (new std::vector<PetscScalar> (values, values+ncols));
// and finally restore the matrix
ierr = MatRestoreRow(*matrix, this->a_row, &ncols, &colnums, &values);
+ 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
{