template <typename> class Vector;
template <typename> class BlockVector;
+template <typename> class BlockSparseMatrix;
/*! @addtogroup Matrix1
*@{
*/
+
+namespace internal
+{
+
/**
- * Blocked sparse matrix. The behaviour of objects of this type is
- * almost as for the SparseMatrix objects, with most of the
- * functions being implemented in both classes. The main difference is
- * that the matrix represented by this object is composed of an array
- * of sparse matrices (i.e. of type SparseMatrix<number>) and all
- * accesses to the elements of this object are relayed to accesses of
- * the base matrices.
- *
- * In addition to the usual matrix access and linear algebra
- * functions, there are functions block() which allow access to the
- * different blocks of the matrix. This may, for example, be of help
- * when you want to implement Schur complement methods, or block
- * preconditioners, where each block belongs to a specific component
- * of the equation you are presently discretizing.
- *
- * Note that the numbers of blocks and rows are implicitly determined
- * by the sparsity pattern objects used.
- *
- * Objects of this type are frequently used when a system of differential
- * equations has solutions with variables that fall into different
- * classes. For example, solutions of the Stokes or Navier-Stokes equations
- * have @p dim velocity components and one pressure component. In this case,
- * it may make sense to consider the linear system of equations as a system of
- * 2x2 blocks, and one can construct preconditioners or solvers based on this
- * 2x2 block structure. This class can help you in these cases, as it allows
- * to view the matrix alternatively as one big matrix, or as a number of
- * individual blocks.
- *
- * @ref Instantiations: some (<tt>@<float@> @<double@></tt>)
+ * Namespace in which iterators in block matrices are implemented.
*
- * @author Wolfgang Bangerth, 2000
+ * @author Wolfgang Bangerth, 2004
*/
-template <typename number>
-class BlockSparseMatrix : public Subscriptor
-{
- public:
- /**
- * Type of matrix entries. In analogy to
- * the STL container classes.
- */
- typedef number value_type;
-
- class const_iterator;
+ namespace BlockMatrixIterators
+ {
+ template <typename> class ConstIterator;
+
/**
* Accessor class for iterators
*/
+ template <typename number>
class Accessor
{
public:
+ /**
+ * Typedef the value type of the
+ * matrix we point into.
+ */
+ typedef
+ typename BlockSparseMatrix<number>::value_type
+ value_type;
+
/**
* Constructor. Since we use
* accessors only for read
*/
Accessor (const BlockSparseMatrix<number> *m,
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
/**
* Value of this matrix entry.
*/
- number value() const;
+ value_type value() const;
/**
* Block row of the
* Iterator of the underlying matrix
* class.
*/
- typename SparseMatrix<number>::const_iterator base_iterator;
+ typename
+ BlockSparseMatrix<number>::BlockType::const_iterator
+ base_iterator;
/**
* Number of block where row lies in.
*/
unsigned int a_index;
- friend class const_iterator;
+ /**
+ * Let the iterator class be a
+ * friend.
+ */
+ template <typename>
+ friend class ConstIterator;
};
/**
* STL conforming iterator.
*/
- class const_iterator : private Accessor
+ template <typename number>
+ class ConstIterator : private Accessor<number>
{
public:
/**
* Constructor.
*/
- const_iterator(const BlockSparseMatrix<number>*,
- unsigned int row,
- unsigned short index);
+ ConstIterator(const BlockSparseMatrix<number>*,
+ const unsigned int row,
+ const unsigned int index);
/**
* Prefix increment.
*/
- const_iterator& operator++ ();
+ ConstIterator& operator++ ();
/**
* Postfix increment.
*/
- const_iterator& operator++ (int);
+ ConstIterator& operator++ (int);
/**
* Dereferencing operator.
*/
- const Accessor& operator* () const;
+ const Accessor<number> & operator* () const;
/**
* Dereferencing operator.
*/
- const Accessor* operator-> () const;
+ const Accessor<number> * operator-> () const;
/**
* Comparison. True, if
* the same matrix
* position.
*/
- bool operator == (const const_iterator&) const;
+ bool operator == (const ConstIterator&) const;
/**
* Inverse of operator==().
*/
- bool operator != (const const_iterator&) const;
+ bool operator != (const ConstIterator&) const;
/**
* Comparison
* equal and the first
* index is smaller.
*/
- bool operator < (const const_iterator&) const;
+ bool operator < (const ConstIterator&) const;
};
+ }
+}
+
+
+/**
+ * Blocked sparse matrix. The behaviour of objects of this type is
+ * almost as for the SparseMatrix objects, with most of the
+ * functions being implemented in both classes. The main difference is
+ * that the matrix represented by this object is composed of an array
+ * of sparse matrices (i.e. of type SparseMatrix<number>) and all
+ * accesses to the elements of this object are relayed to accesses of
+ * the base matrices.
+ *
+ * In addition to the usual matrix access and linear algebra
+ * functions, there are functions block() which allow access to the
+ * different blocks of the matrix. This may, for example, be of help
+ * when you want to implement Schur complement methods, or block
+ * preconditioners, where each block belongs to a specific component
+ * of the equation you are presently discretizing.
+ *
+ * Note that the numbers of blocks and rows are implicitly determined
+ * by the sparsity pattern objects used.
+ *
+ * Objects of this type are frequently used when a system of differential
+ * equations has solutions with variables that fall into different
+ * classes. For example, solutions of the Stokes or Navier-Stokes equations
+ * have @p dim velocity components and one pressure component. In this case,
+ * it may make sense to consider the linear system of equations as a system of
+ * 2x2 blocks, and one can construct preconditioners or solvers based on this
+ * 2x2 block structure. This class can help you in these cases, as it allows
+ * to view the matrix alternatively as one big matrix, or as a number of
+ * individual blocks.
+ *
+ * @ref Instantiations: some (<tt>@<float@> @<double@></tt>)
+ *
+ * @author Wolfgang Bangerth, 2000
+ */
+template <typename number>
+class BlockSparseMatrix : public Subscriptor
+{
+ public:
+ /**
+ * Typedef the type of the underlying
+ * matrix.
+ */
+ typedef SparseMatrix<number> BlockType;
+
+ /**
+ * Type of matrix entries. In analogy to
+ * the STL container classes.
+ */
+ typedef typename BlockType::value_type value_type;
+ typedef value_type *pointer;
+ typedef const value_type *const_pointer;
+ typedef internal::BlockMatrixIterators::ConstIterator<number> iterator;
+ typedef internal::BlockMatrixIterators::ConstIterator<number> const_iterator;
+ typedef value_type &reference;
+ typedef const value_type &const_reference;
+ typedef size_t size_type;
/**
* Constructor; initializes the
* Access the block with the
* given coordinates.
*/
- SparseMatrix<number> &
+ BlockType &
block (const unsigned int row,
const unsigned int column);
* given coordinates. Version for
* constant objects.
*/
- const SparseMatrix<number> &
+ const BlockType &
block (const unsigned int row,
const unsigned int column) const;
* STL-like iterator with the
* first entry of row <tt>r</tt>.
*/
- const_iterator begin (unsigned int r) const;
+ const_iterator begin (const unsigned int r) const;
/**
* Final iterator of row <tt>r</tt>.
*/
- const_iterator end (unsigned int r) const;
+ const_iterator end (const unsigned int r) const;
/**
* Determine an estimate for the
/**
* Array of sub-matrices.
*/
- Table<2,SmartPointer<SparseMatrix<number> > > sub_objects;
+ Table<2,SmartPointer<BlockType> > sub_objects;
+ /**
+ * Make the iterator class a
+ * friend. We have to work around
+ * a compiler bug here again.
+ */
+#ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG
+ template <typename>
+ friend class internal::BlockMatrixIterators::Accessor;
+
+ template <typename>
+ friend class internal::BlockMatrixIterators::ConstIterator;
+#else
+ typedef internal::BlockMatrixIterators::Accessor<number> Accessor;
friend class Accessor;
+
friend class const_iterator;
+#endif
};
/*@}*/
/* ------------------------- Template functions ---------------------- */
-template <typename number>
-inline
-BlockSparseMatrix<number>::Accessor::
-Accessor (const BlockSparseMatrix<number> *matrix,
- const unsigned int r,
- const unsigned short i)
- :
- matrix(matrix),
- base_iterator(matrix->block(0,0).begin()),
- row_block(0),
- row_start(0),
- col_block(0),
- col_start(0),
- a_index(0)
-{
- Assert (i==0, ExcNotImplemented());
- if (r < matrix->m())
- {
- std::pair<unsigned int,unsigned int> indices
- = matrix->sparsity_pattern->get_row_indices().global_to_local(r);
- row_block = indices.first;
- base_iterator = matrix->block(indices.first, 0).begin(indices.second);
- row_start = matrix->sparsity_pattern
- ->get_row_indices().local_to_global(row_block, 0);
- }
- else
+namespace internal
+{
+ namespace BlockMatrixIterators
+ {
+ template <typename number>
+ inline
+ Accessor<number>::
+ Accessor (const BlockSparseMatrix<number> *matrix,
+ const unsigned int r,
+ const unsigned int i)
+ :
+ matrix(matrix),
+ base_iterator(matrix->block(0,0).begin()),
+ row_block(0),
+ row_start(0),
+ col_block(0),
+ col_start(0),
+ a_index(0)
{
- row_block = matrix->n_block_rows();
- base_iterator = matrix->block(0, 0).begin();
+ Assert (i==0, ExcNotImplemented());
+
+ if (r < matrix->m())
+ {
+ std::pair<unsigned int,unsigned int> indices
+ = matrix->sparsity_pattern->get_row_indices().global_to_local(r);
+ row_block = indices.first;
+ base_iterator = matrix->block(indices.first, 0).begin(indices.second);
+ row_start = matrix->sparsity_pattern
+ ->get_row_indices().local_to_global(row_block, 0);
+ }
+ else
+ {
+ row_block = matrix->n_block_rows();
+ base_iterator = matrix->block(0, 0).begin();
+ }
}
-}
-template <typename number>
-inline
-unsigned int
-BlockSparseMatrix<number>::Accessor::row() const
-{
- return row_start + base_iterator->row();
-}
+ template <typename number>
+ inline
+ unsigned int
+ Accessor<number>::row() const
+ {
+ return row_start + base_iterator->row();
+ }
-template <typename number>
-inline
-short unsigned int
-BlockSparseMatrix<number>::Accessor::index() const
-{
- return a_index;
-}
+ template <typename number>
+ inline
+ unsigned int
+ Accessor<number>::index() const
+ {
+ return a_index;
+ }
-template <typename number>
-inline
-unsigned int
-BlockSparseMatrix<number>::Accessor::column() const
-{
- return col_start + base_iterator->column();
-}
+ template <typename number>
+ inline
+ unsigned int
+ Accessor<number>::column() const
+ {
+ return col_start + base_iterator->column();
+ }
-template <typename number>
-inline
-unsigned int
-BlockSparseMatrix<number>::Accessor::block_row() const
-{
- return row_block;
-}
+ template <typename number>
+ inline
+ unsigned int
+ Accessor<number>::block_row() const
+ {
+ return row_block;
+ }
-template <typename number>
-inline
-unsigned int
-BlockSparseMatrix<number>::Accessor::block_column() const
-{
- return col_block;
-}
+ template <typename number>
+ inline
+ unsigned int
+ Accessor<number>::block_column() const
+ {
+ return col_block;
+ }
-template <typename number>
-inline
-number
-BlockSparseMatrix<number>::Accessor::value () const
-{
- return base_iterator->value();
-}
+ template <typename number>
+ inline
+ typename Accessor<number>::value_type
+ Accessor<number>::value () const
+ {
+ return base_iterator->value();
+ }
//----------------------------------------------------------------------//
-template <typename number>
-inline
-BlockSparseMatrix<number>::const_iterator::
-const_iterator(const BlockSparseMatrix<number>* m,
- unsigned int r,
- unsigned short i)
- :
- BlockSparseMatrix<number>::Accessor(m, r, i)
-{}
+ template <typename number>
+ inline
+ ConstIterator<number>::
+ ConstIterator (const BlockSparseMatrix<number>* m,
+ const unsigned int r,
+ const unsigned int i)
+ :
+ Accessor<number>(m, r, i)
+ {}
-template <typename number>
-inline
-typename BlockSparseMatrix<number>::const_iterator&
-BlockSparseMatrix<number>::const_iterator::operator++ ()
-{
- Assert (this->row_block<this->matrix->n_block_rows(), ExcIteratorPastEnd());
+ template <typename number>
+ inline
+ ConstIterator<number> &
+ ConstIterator<number>::operator++ ()
+ {
+ Assert (this->row_block<this->matrix->n_block_rows(), ExcIteratorPastEnd());
- // Remember current row inside block
- unsigned int local_row = this->base_iterator->row();
+ // Remember current row inside block
+ unsigned int local_row = this->base_iterator->row();
- // Advance inside block
- ++this->base_iterator;
- ++this->a_index;
+ // Advance inside block
+ ++this->base_iterator;
+ ++this->a_index;
- // If end of row inside block,
- // advance to next block
- if (this->base_iterator ==
- this->matrix->block(this->row_block, this->col_block).end(local_row))
- {
- if (this->col_block < this->matrix->n_block_cols()-1)
- {
- // Advance to next block in
- // row
+ // If end of row inside block,
+ // advance to next block
+ if (this->base_iterator ==
+ this->matrix->block(this->row_block, this->col_block).end(local_row))
+ {
+ if (this->col_block < this->matrix->n_block_cols()-1)
+ {
+ // Advance to next block in
+ // row
//TODO: what if this row in that block is empty?
- ++this->col_block;
- this->col_start = this->matrix->sparsity_pattern
- ->get_column_indices().local_to_global(this->col_block, 0);
- }
- else
- {
- // Advance to first block
- // in next row
- this->col_block = 0;
- this->col_start = 0;
- this->a_index = 0;
- ++local_row;
- if (local_row >= this->matrix->block(this->row_block,0).m())
- {
- // If final row in
- // block, go to next
- // block row
- local_row = 0;
- ++this->row_block;
- if (this->row_block < this->matrix->n_block_rows())
- this->row_start = this->matrix->sparsity_pattern
- ->get_row_indices().local_to_global(this->row_block, 0);
- }
- }
- // Finally, set base_iterator
- // to start of row determined
- // above
- if (this->row_block < this->matrix->n_block_rows())
- this->base_iterator = this->matrix->block(this->row_block, this->col_block).begin(local_row);
- else
- // Set base_iterator to a
- // defined state for
- // comparison. This is the
- // end() state.
+ ++this->col_block;
+ this->col_start = this->matrix->sparsity_pattern
+ ->get_column_indices().local_to_global(this->col_block, 0);
+ }
+ else
+ {
+ // Advance to first block
+ // in next row
+ this->col_block = 0;
+ this->col_start = 0;
+ this->a_index = 0;
+ ++local_row;
+ if (local_row >= this->matrix->block(this->row_block,0).m())
+ {
+ // If final row in
+ // block, go to next
+ // block row
+ local_row = 0;
+ ++this->row_block;
+ if (this->row_block < this->matrix->n_block_rows())
+ this->row_start = this->matrix->sparsity_pattern
+ ->get_row_indices().local_to_global(this->row_block, 0);
+ }
+ }
+ // Finally, set base_iterator
+ // to start of row determined
+ // above
+ if (this->row_block < this->matrix->n_block_rows())
+ this->base_iterator = this->matrix->block(this->row_block, this->col_block).begin(local_row);
+ else
+ // Set base_iterator to a
+ // defined state for
+ // comparison. This is the
+ // end() state.
//TODO: this is a particularly bad choice, since it is actually a valid iterator. it should rather be something like the end iterator of the last block!
- this->base_iterator = this->matrix->block(0, 0).begin();
+ this->base_iterator = this->matrix->block(0, 0).begin();
+ }
+ return *this;
}
- return *this;
-}
// template <typename number>
// inline
// const_iterator&
-// BlockSparseMatrix<number>::const_iterator::operator++ (int)
+// ConstIterator::operator++ (int)
// {
// Assert (false, ExcNotImplemented());
// }
-template <typename number>
-inline
-const typename BlockSparseMatrix<number>::Accessor&
-BlockSparseMatrix<number>::const_iterator::operator* () const
-{
- return *this;
-}
+ template <typename number>
+ inline
+ const Accessor<number> &
+ ConstIterator<number>::operator* () const
+ {
+ return *this;
+ }
-template <typename number>
-inline
-const typename BlockSparseMatrix<number>::Accessor*
-BlockSparseMatrix<number>::const_iterator::operator-> () const
-{
- return this;
-}
+ template <typename number>
+ inline
+ const Accessor<number> *
+ ConstIterator<number>::operator-> () const
+ {
+ return this;
+ }
-template <typename number>
-inline
-bool
-BlockSparseMatrix<number>::const_iterator::
-operator == (const const_iterator& i) const
-{
- if (this->matrix != i->matrix)
- return false;
+ template <typename number>
+ inline
+ bool
+ ConstIterator<number>::
+ operator == (const ConstIterator& i) const
+ {
+ if (this->matrix != i->matrix)
+ return false;
- if (this->row_block == i->row_block
- && this->col_block == i->col_block
- && this->base_iterator == i->base_iterator)
- return true;
- return false;
-}
+ if (this->row_block == i->row_block
+ && this->col_block == i->col_block
+ && this->base_iterator == i->base_iterator)
+ return true;
+ return false;
+ }
-template <typename number>
-inline
-bool
-BlockSparseMatrix<number>::const_iterator::
-operator != (const const_iterator& i) const
-{
- return !(*this == i);
-}
+ template <typename number>
+ inline
+ bool
+ ConstIterator<number>::
+ operator != (const ConstIterator& i) const
+ {
+ return !(*this == i);
+ }
-template <typename number>
-inline
-bool
-BlockSparseMatrix<number>::const_iterator::
-operator < (const const_iterator& i) const
-{
- if (this->row_block < i->row_block)
- return true;
- if (this->row_block == i->row_block)
+ template <typename number>
+ inline
+ bool
+ ConstIterator<number>::
+ operator < (const ConstIterator& i) const
{
- if (this->base_iterator->row() < i->base_iterator->row())
- return true;
- if (this->base_iterator->row() == i->base_iterator->row())
- {
- if (this->a_index < i->a_index)
- return true;
- }
+ if (this->row_block < i->row_block)
+ return true;
+ if (this->row_block == i->row_block)
+ {
+ if (this->base_iterator->row() < i->base_iterator->row())
+ return true;
+ if (this->base_iterator->row() == i->base_iterator->row())
+ {
+ if (this->a_index < i->a_index)
+ return true;
+ }
+ }
+ return false;
}
- return false;
+
+ }
}
//----------------------------------------------------------------------//
template <typename number>
inline
-SparseMatrix<number> &
+typename BlockSparseMatrix<number>::BlockType &
BlockSparseMatrix<number>::block (const unsigned int row,
const unsigned int column)
{
template <typename number>
inline
-const SparseMatrix<number> &
+const typename BlockSparseMatrix<number>::BlockType &
BlockSparseMatrix<number>::block (const unsigned int row,
const unsigned int column) const
{
template <typename number>
inline
typename BlockSparseMatrix<number>::const_iterator
-BlockSparseMatrix<number>::begin (unsigned int r) const
+BlockSparseMatrix<number>::begin (const unsigned int r) const
{
Assert (r<m(), ExcIndexRange(r,0,m()));
return const_iterator(this, r, 0);
template <typename number>
inline
typename BlockSparseMatrix<number>::const_iterator
-BlockSparseMatrix<number>::end (unsigned int r) const
+BlockSparseMatrix<number>::end (const unsigned int r) const
{
Assert (r<m(), ExcIndexRange(r,0,m()));
return const_iterator(this, r+1, 0);