template <typename number> class FullMatrix;
template <typename number> class Vector;
+template<class MATRIX, typename inverse_type>
+class PreconditionBlockJacobi;
+
/*! @addtogroup Preconditioners
*@{
*/
* yields a preconditioner.
*/
void set_same_diagonal ();
+
+ /**
+ * Does the matrix use only one
+ * diagonal block?
+ */
+ bool same_diagonal () const;
/**
* Stores the inverse of the
*/
unsigned int block_size () const;
+ /**
+ * The number of blocks of the
+ * matrix.
+ */
+ unsigned int n_blocks() const;
+
/**
* Determine an estimate for the
* memory consumption (in bytes)
*/
bool store_diagonals;
+ /**
+ * Number of blocks.
+ */
+ unsigned int nblocks;
private:
+
/**
* Storage of the inverse
* matrices of the diagonal
* Flag for diagonal compression.
* @ref set_same_diagonal()
*/
- bool same_diagonal;
+ bool var_same_diagonal;
};
/**
* Block Jacobi preconditioning.
* See @ref{PreconditionBlock} for requirements on the matrix.
- * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
+ * @author Ralf Hartmann, Guido Kanschat, 1999, 2000, 2003
*/
template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
class PreconditionBlockJacobi : public virtual Subscriptor,
typedef typename MATRIX::value_type number;
public:
+ /**
+ * STL conforming iterator.
+ */
+ class const_iterator
+ {
+ private:
+ /**
+ * Accessor class for iterators
+ */
+ class Accessor
+ {
+ public:
+ /**
+ * Constructor. Since we use
+ * accessors only for read
+ * access, a const matrix
+ * pointer is sufficient.
+ */
+ Accessor (const PreconditionBlockJacobi<MATRIX, inverse_type> *matrix,
+ const unsigned int row);
+
+ /**
+ * Row number of the element
+ * represented by this
+ * object.
+ */
+ unsigned int row() const;
+
+ /**
+ * Column number of the
+ * element represented by
+ * this object.
+ */
+ unsigned int column() const;
+
+ /**
+ * Value of this matrix entry.
+ */
+ inverse_type value() const;
+
+ protected:
+ /**
+ * The matrix accessed.
+ */
+ const PreconditionBlockJacobi<MATRIX, inverse_type>* matrix;
+
+ /**
+ * Save block size here
+ * for further reference.
+ */
+ unsigned int bs;
+
+ /**
+ * Current block number.
+ */
+ unsigned int a_block;
+
+ /**
+ * Current block for extracting data.
+ *
+ * This is either
+ * @p{a_block} or zero,
+ * depending on
+ * @p{same_diagonal}.
+ */
+ unsigned int u_block;
+
+ /**
+ * Iterator inside block.
+ */
+ typename FullMatrix<inverse_type>::const_iterator b_iterator;
+
+ /**
+ * End of current block.
+ */
+ typename FullMatrix<inverse_type>::const_iterator b_end;
+
+ /**
+ * Make enclosing class a
+ * friend.
+ */
+ friend class const_iterator;
+ };
+
+ public:
+ /**
+ * Constructor.
+ */
+ const_iterator(const PreconditionBlockJacobi<MATRIX, inverse_type>* matrix,
+ const unsigned int row);
+
+ /**
+ * Prefix increment.
+ */
+ const_iterator& operator++ ();
+
+ /**
+ * Postfix increment.
+ */
+ const_iterator& operator++ (int);
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor& operator* () const;
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor* operator-> () const;
+
+ /**
+ * Comparison. True, if
+ * both iterators point to
+ * the same matrix
+ * position.
+ */
+ bool operator == (const const_iterator&) const;
+ /**
+ * Inverse of @p{==}.
+ */
+ bool operator != (const const_iterator&) const;
+
+ /**
+ * Comparison
+ * operator. Result is true
+ * if either the first row
+ * number is smaller or if
+ * the row numbers are
+ * equal and the first
+ * index is smaller.
+ */
+ bool operator < (const const_iterator&) const;
+
+ private:
+ /**
+ * Store an object of the
+ * accessor class.
+ */
+ Accessor accessor;
+ };
+
/**
* Make initialization function
* publicly available.
template <typename number2>
void Tvmult_add (Vector<number2>&, const Vector<number2>&) const;
+ /**
+ * STL-like iterator with the
+ * first entry.
+ */
+ const_iterator begin () const;
+
+ /**
+ * Final iterator.
+ */
+ const_iterator end () const;
+
+ /**
+ * STL-like iterator with the
+ * first entry of row @p{r}.
+ */
+ const_iterator begin (const unsigned int r) const;
+
+ /**
+ * Final iterator of row @p{r}.
+ */
+ const_iterator end (const unsigned int r) const;
+
+
private:
/**
* Actual implementation of the
}
+template<class MATRIX, typename inverse_type>
+inline unsigned int
+PreconditionBlock<MATRIX, inverse_type>::n_blocks () const
+{
+ return nblocks;
+}
+
+
template<class MATRIX, typename inverse_type>
inline inverse_type
PreconditionBlock<MATRIX, inverse_type>::el (
return B(ib, jb);
}
+//----------------------------------------------------------------------//
+
+template<class MATRIX, typename inverse_type>
+inline
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::
+Accessor (const PreconditionBlockJacobi<MATRIX, inverse_type>* matrix,
+ const unsigned int row)
+ :
+ matrix(matrix),
+ b_iterator(&matrix->inverse(0), 0, 0),
+ b_end(&matrix->inverse(0), 0, 0)
+{
+ bs = matrix->block_size();
+ a_block = row / bs;
+
+ // This is the end accessor, which
+ // does not hava a valid block.
+ if (a_block == matrix->n_blocks())
+ return;
+
+ const unsigned int r = row % bs;
+
+ u_block = (matrix->same_diagonal()) ? 0 : a_block;
+
+ b_iterator = matrix->inverse(u_block).begin(r);
+ b_end = matrix->inverse(u_block).end();
+
+ Assert (a_block < matrix->n_blocks(),
+ ExcIndexRange(a_block, 0, matrix->n_blocks()));
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+unsigned int
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::row() const
+{
+ Assert (a_block < matrix->n_blocks(),
+ ExcIteratorPastEnd());
+
+ return bs * a_block + b_iterator->row();
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+unsigned int
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::column() const
+{
+ Assert (a_block < matrix->n_blocks(),
+ ExcIteratorPastEnd());
+
+ return bs * a_block + b_iterator->column();
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+inverse_type
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::value() const
+{
+ Assert (a_block < matrix->n_blocks(),
+ ExcIteratorPastEnd());
+
+ return b_iterator->value();
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
+const_iterator(const PreconditionBlockJacobi<MATRIX, inverse_type> *matrix,
+ const unsigned int row)
+ :
+ accessor(matrix, row)
+{}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator &
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::operator++ ()
+{
+ Assert (*this != accessor.matrix->end(), ExcIteratorPastEnd());
+
+ ++accessor.b_iterator;
+ if (accessor.b_iterator == accessor.b_end)
+ {
+ ++accessor.a_block;
+ if (!accessor.matrix->same_diagonal())
+ ++accessor.u_block;
+
+ if (accessor.a_block < accessor.matrix->n_blocks())
+ {
+ accessor.b_iterator = accessor.matrix->inverse(accessor.a_block).begin();
+ accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
+ }
+ }
+ return *this;
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+const typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor &
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::operator* () const
+{
+ return accessor;
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+const typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor *
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::operator-> () const
+{
+ return &accessor;
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+bool
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
+operator == (const const_iterator& other) const
+{
+ if (accessor.a_block == accessor.matrix->n_blocks() &&
+ accessor.a_block == other.accessor.a_block)
+ return true;
+
+ if (accessor.a_block != other.accessor.a_block)
+ return false;
+
+ return (accessor.row() == other.accessor.row() &&
+ accessor.column() == other.accessor.column());
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+bool
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
+operator != (const const_iterator& other) const
+{
+ return ! (*this == other);
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+bool
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
+operator < (const const_iterator& other) const
+{
+ return (accessor.row() < other.accessor.row() ||
+ (accessor.row() == other.accessor.row() &&
+ accessor.column() < other.accessor.column()));
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
+PreconditionBlockJacobi<MATRIX, inverse_type>::begin () const
+{
+ return const_iterator(this, 0);
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
+PreconditionBlockJacobi<MATRIX, inverse_type>::end () const
+{
+ return const_iterator(this, n_blocks() * block_size());
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
+PreconditionBlockJacobi<MATRIX, inverse_type>::begin (
+ const unsigned int r) const
+{
+ Assert (r<A->m(), ExcIndexRange(r,0,A->m()));
+ return const_iterator(this, r);
+}
+
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
+PreconditionBlockJacobi<MATRIX, inverse_type>::end (
+ const unsigned int r) const
+{
+ Assert (r<A->m(), ExcIndexRange(r,0,A->m()));
+ return const_iterator(this, r+1);
+}
+
#endif
blocksize(0),
A(0),
store_diagonals(false),
- same_diagonal(false)
+ var_same_diagonal(false)
{}
if (var_diagonal.size()!=0)
var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
blocksize = 0;
- same_diagonal = false;
+ var_same_diagonal = false;
A = 0;
}
Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
blocksize=bsize;
relaxation = parameters.relaxation;
- same_diagonal = parameters.same_diagonal;
+ var_same_diagonal = parameters.same_diagonal;
+ nblocks = A->m()/bsize;
+
if (parameters.invert_diagonal)
invert_diagblocks();
}
const FullMatrix<inverse_type>&
PreconditionBlock<MATRIX,inverse_type>::inverse(unsigned int i) const
{
- if (same_diagonal)
+ if (same_diagonal())
return var_inverse[0];
Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
{
Assert(store_diagonals, ExcDiagonalsNotStored());
- if (same_diagonal)
+ if (same_diagonal())
return var_diagonal[0];
Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size()));
PreconditionBlock<MATRIX,inverse_type>::set_same_diagonal()
{
Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
- same_diagonal = true;
+ var_same_diagonal = true;
+}
+
+
+template <class MATRIX, typename inverse_type>
+bool
+PreconditionBlock<MATRIX,inverse_type>::same_diagonal() const
+{
+ return var_same_diagonal;
}
const MATRIX &M=*A;
Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
- const unsigned int n_cells = M.m()/blocksize;
-
FullMatrix<inverse_type> M_cell(blocksize);
- if (same_diagonal)
+ if (same_diagonal())
{
deallog << "PreconditionBlock uses only one diagonal block" << std::endl;
// Invert only the first block
// set the @p{var_inverse} array to the right
// size. we could do it like this:
- // var_inverse = vector<>(n_cells,FullMatrix<>())
+ // var_inverse = vector<>(nblocks,FullMatrix<>())
// but this would involve copying many
// FullMatrix objects.
//
if (store_diagonals)
{
std::vector<FullMatrix<inverse_type> >
- tmp(n_cells, FullMatrix<inverse_type>(blocksize));
+ tmp(nblocks, FullMatrix<inverse_type>(blocksize));
var_diagonal.swap (tmp);
}
if (true)
{
std::vector<FullMatrix<inverse_type> >
- tmp(n_cells, FullMatrix<inverse_type>(blocksize));
+ tmp(nblocks, FullMatrix<inverse_type>(blocksize));
var_inverse.swap (tmp);
// make sure the tmp object
// goes out of scope as
M_cell.clear ();
- for (unsigned int cell=0; cell<n_cells; ++cell)
+ for (unsigned int cell=0; cell<nblocks; ++cell)
{
const unsigned int cell_start = cell*blocksize;
for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
Assert(this->A!=0, ExcNotInitialized());
const MATRIX &M=*this->A;
- const unsigned int n_cells=M.m()/this->blocksize;
Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
if (!this->inverses_ready())
{
FullMatrix<number> M_cell(this->blocksize);
- for (unsigned int cell=0; cell<n_cells; ++cell)
+ for (unsigned int cell=0; cell<nblocks; ++cell)
{
for (row=cell*this->blocksize, row_cell=0;
row_cell<this->blocksize;
}
}
else
- for (unsigned int cell=0; cell<n_cells; ++cell)
+ for (unsigned int cell=0; cell<nblocks; ++cell)
{
for (row=cell*this->blocksize, row_cell=0;
row_cell<this->blocksize;
Assert (this->A!=0, ExcNotInitialized());
const MATRIX &M=*this->A;
- const unsigned int n_cells=M.m()/this->blocksize;
Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
if (!this->inverses_ready())
{
FullMatrix<number> M_cell(this->blocksize);
- for (unsigned int cell=0; cell<n_cells; ++cell)
+ for (unsigned int cell=0; cell<nblocks; ++cell)
{
for (row = block_start, row_cell = 0;
row_cell < this->blocksize;
}
}
else
- for (unsigned int cell=0; cell<n_cells; ++cell)
+ for (unsigned int cell=0; cell<nblocks; ++cell)
{
for (row = block_start, row_cell = 0;
row_cell<this->blocksize;
Assert (this->A!=0, ExcNotInitialized());
const MATRIX &M=*this->A;
- const unsigned int n_cells=M.m()/this->blocksize;
Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
// row, column are the global numbering
// of the unkowns.
unsigned int row, row_cell, block_start = 0;
- unsigned int block_end=this->blocksize *n_cells;
+ unsigned int block_end=this->blocksize *nblocks;
number2 b_cell_row;
if (!this->inverses_ready())
{
FullMatrix<number> M_cell(this->blocksize);
- for (int icell=n_cells-1; icell>=0 ; --icell)
+ for (int icell=nblocks-1; icell>=0 ; --icell)
{
unsigned int cell = (unsigned int) icell;
// Collect upper triangle
}
}
else
- for (int icell=n_cells-1; icell >=0 ; --icell)
+ for (int icell=nblocks-1; icell >=0 ; --icell)
{
unsigned int cell = (unsigned int) icell;
for (row=cell*this->blocksize, row_cell=0;
Vector<inverse_type> cell_src(this->blocksize);
Vector<inverse_type> cell_dst(this->blocksize);
- const unsigned int n_cells = this->A->m()/this->blocksize;
-
// Multiply with diagonal blocks
- for (unsigned int cell=0; cell<n_cells; ++cell)
+ for (unsigned int cell=0; cell<nblocks; ++cell)
{
unsigned int row = cell*this->blocksize;
Vector<inverse_type> cell_src(this->blocksize);
Vector<inverse_type> cell_dst(this->blocksize);
- const unsigned int n_cells = this->A->m()/this->blocksize;
-
// Multiply with diagonal blocks
- for (unsigned int cell=0; cell<n_cells; ++cell)
+ for (unsigned int cell=0; cell<nblocks; ++cell)
{
unsigned int row = cell*this->blocksize;