*/
Accessor (const MatrixBase *matrix,
const unsigned int row,
- const unsigned short index);
+ const unsigned int index);
/**
* Row number of the element
* represented by this
* object.
*/
- unsigned short index() const;
+ unsigned int index() const;
/**
* Column number of the
/**
* The matrix accessed.
*/
- const MatrixBase *matrix;
+ mutable MatrixBase *matrix;
/**
* Current row number.
* into the matrix @p matrix for the
* given row and the index within it.
*/
- const_iterator (const MatrixBase *matrix,
- const unsigned int row,
- const unsigned short index);
+ const_iterator (const MatrixBase *matrix,
+ const unsigned int row,
+ const unsigned int index);
/**
* Prefix increment.
PetscScalar el (const unsigned int i,
const unsigned int j) const;
+ /**
+ * Return the main diagonal
+ * element in the <i>i</i>th
+ * row. This function throws an
+ * error if the matrix is not
+ * quadratic.
+ *
+ * Since we do not have direct access
+ * to the underlying data structure,
+ * this function is no faster than the
+ * elementwise access using the el()
+ * function. However, we provide this
+ * function for compatibility with the
+ * SparseMatrix class.
+ */
+ PetscScalar diag_element (const unsigned int i) const;
+
/**
* Return the number of rows in this
* matrix.
* Exception
*/
DeclException0 (ExcSourceEqualsDestination);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcMatrixNotSquare);
protected:
/**
inline
const_iterator::Accessor::
- Accessor (const MatrixBase *matrix,
- const unsigned int row,
- const unsigned short index)
+ Accessor (const MatrixBase *matrix,
+ const unsigned int row,
+ const unsigned int index)
:
- matrix(matrix),
+ matrix(const_cast<MatrixBase*>(matrix)),
a_row(row),
a_index(index)
{
visit_present_row ();
- Assert (index < colnum_cache->size(),
+ Assert ((row == matrix->m()) ||
+ (index < colnum_cache->size()),
ExcInvalidIndexWithinRow (index, row));
}
inline
- unsigned short
+ unsigned int
const_iterator::Accessor::index() const
{
Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
inline
const_iterator::
- const_iterator(const MatrixBase *matrix,
- const unsigned int row,
- const unsigned short index)
+ const_iterator(const MatrixBase *matrix,
+ const unsigned int row,
+ const unsigned int index)
:
accessor(matrix, row, index)
{}
const_iterator::
operator == (const const_iterator& other) const
{
- return (accessor.row() == other.accessor.row() &&
- accessor.index() == other.accessor.index());
+ return (accessor.a_row == other.accessor.a_row &&
+ accessor.a_index == other.accessor.a_index);
}
MatrixBase::const_iterator::Accessor::
visit_present_row ()
{
- const int a_row_int = this->a_row; // need signed int for PETSc?
+ // 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())
+ {
+ colnum_cache.reset ();
+ value_cache.reset ();
+
+ return;
+ }
+
+ // otherwise first flush PETSc caches
+ matrix->compress ();
+
+ const signed int a_row_int = this->a_row;
// get a representation of the present
// row
int ncols;
int *colnums;
PetscScalar *values;
+
int ierr;
ierr = MatGetRow(*matrix, a_row_int, &ncols, &colnums, &values);
AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr));
// copy it into our caches
- colnum_cache =
- boost::shared_ptr<std::vector<unsigned int> >
- ( new std::vector<unsigned int> (colnums, colnums+ncols));
- value_cache =
- boost::shared_ptr<std::vector<PetscScalar> >
- ( new std::vector<PetscScalar> (values, values+ncols));
+ 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, a_row_int, &ncols, &colnums, &values);
}
+
+ PetscScalar
+ MatrixBase::diag_element (const unsigned int i) const
+ {
+ Assert (m() == n(), ExcMatrixNotSquare());
+
+ // this doesn't seem to work any
+ // different than any other element
+ return el(i,i);
+ }
+
+
+
void
MatrixBase::compress ()
{