From: guido Date: Tue, 13 Nov 2001 22:28:56 +0000 (+0000) Subject: BlockTrianglePrecondition X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fba86afdf10dbacd3a09b7e2332f627edcb05dad;p=dealii-svn.git BlockTrianglePrecondition git-svn-id: https://svn.dealii.org/trunk@5192 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/block_matrix_array.h b/deal.II/lac/include/lac/block_matrix_array.h index 1ce794b928..4b3fe35011 100644 --- a/deal.II/lac/include/lac/block_matrix_array.h +++ b/deal.II/lac/include/lac/block_matrix_array.h @@ -30,6 +30,7 @@ # define STRINGSTREAM std::ostrstream #endif +//TODO:[GK] Obtain aux vector from VectorMemory template class BlockVector; @@ -91,34 +92,34 @@ class BlockMatrixArray : public Subscriptor /** * Matrix-vector multiplication. */ - template - void vmult (BlockVector& dst, - const BlockVector& src) const; + template + void vmult (BlockVector& dst, + const BlockVector& src) const; /** * Matrix-vector multiplication * adding to @p{dst}. */ - template - void vmult_add (BlockVector& dst, - const BlockVector& src) const; + template + void vmult_add (BlockVector& dst, + const BlockVector& src) const; /** * Transposed matrix-vector * multiplication. */ - template - void Tvmult (BlockVector& dst, - const BlockVector& src) const; + template + void Tvmult (BlockVector& dst, + const BlockVector& src) const; /** * Transposed matrix-vector * multiplication adding to * @p{dst}. */ - template - void Tvmult_add (BlockVector& dst, - const BlockVector& src) const; + template + void Tvmult_add (BlockVector& dst, + const BlockVector& src) const; /** * Print the block structure as a @@ -126,7 +127,7 @@ class BlockMatrixArray : public Subscriptor */ void print_latex (ostream& out) const; - private: + protected: /** * Internal data structure. * @@ -187,6 +188,7 @@ class BlockMatrixArray : public Subscriptor */ vector entries; + private: /** * Number of blocks per column. */ @@ -198,6 +200,122 @@ class BlockMatrixArray : public Subscriptor }; +/** + * Inversion of a block-triangular matrix. + * + * In this block matrix, the inverses of the diagonal blocks are + * stored together with the off-diagonal blocks of a block + * matrix. Then, forward or backward insertion is performed + * block-wise. The diagonal blocks are NOT inverted for this purpose! + * + * While block indices may be duplicated (see @see{BlockMatrixArray}) + * to add blocks, this is not allowed on the diagonal. A short + * computation reveals why. + * + * Like for all preconditioners, the preconditioning operation is + * performed by the @p{vmult} member function. + * + * The implementation may be a little clumsy, but it should be + * sufficient as long as the block sizes are much larger than the + * number of blocks. + * + * @author Guido Kanschat, 2001 + */ +template +class BlockTrianglePrecondition : private BlockMatrixArray +{ + public: + /** + * Constructor. This matrix must be + * block-quadratic. The additional + * parameter allows for backward + * insertion instead of forward. + */ + BlockTrianglePrecondition(unsigned int block_rows, + bool backward = false); + + + /** + * Preconditioning. + */ + template + void vmult (BlockVector& dst, + const BlockVector& src) const; + + /** + * Preconditioning + * adding to @p{dst}. + */ + template + void vmult_add (BlockVector& dst, + const BlockVector& src) const; + + /** + * Transposed preconditioning + */ + template + void Tvmult (BlockVector& dst, + const BlockVector& src) const; + + /** + * Transposed preconditioning + * adding to @p{dst}. + */ + template + void Tvmult_add (BlockVector& dst, + const BlockVector& src) const; + + /** + * Make function of base class available. + */ + BlockMatrixArray::enter; + + /** + * Make function of base class available. + */ + BlockMatrixArray::print_latex; + + /** + * Make function of base class available. + */ + BlockMatrixArray::n_block_rows; + + /** + * Make function of base class available. + */ + BlockMatrixArray::n_block_cols; + BlockMatrixArray::Subscriptor::subscribe; + BlockMatrixArray::Subscriptor::unsubscribe; + + + /** + * Multiple diagonal element. + */ + DeclException1(ExcMultipleDiagonal, + unsigned int, + << "Inverse diagonal entries may not be added in block " + << arg1); + + private: + /** + * Add all off-diagonal + * contributions and return the + * entry of the diagonal element + * for one row. + */ + template + void do_row (BlockVector& dst, + unsigned int row_num) const; + + /** + * Flag for backward insertion. + */ + bool backward; +}; + + +//----------------------------------------------------------------------// + template inline BlockMatrixArray::Entry::Entry (const MATRIX& matrix, @@ -236,18 +354,18 @@ BlockMatrixArray::enter (const MATRIX& matrix, template -template +template inline void -BlockMatrixArray::vmult_add (BlockVector& dst, - const BlockVector& src) const +BlockMatrixArray::vmult_add (BlockVector& dst, + const BlockVector& src) const { Assert (dst.n_blocks() == block_rows, ExcDimensionMismatch(dst.n_blocks(), block_rows)); Assert (src.n_blocks() == block_cols, ExcDimensionMismatch(src.n_blocks(), block_cols)); - static Vector aux; + static Vector aux; typename vector::const_iterator m = entries.begin(); typename vector::const_iterator end = entries.end(); @@ -279,11 +397,11 @@ BlockMatrixArray::vmult_add (BlockVector& dst, template -template +template inline void -BlockMatrixArray::vmult (BlockVector& dst, - const BlockVector& src) const +BlockMatrixArray::vmult (BlockVector& dst, + const BlockVector& src) const { dst = 0.; vmult_add (dst, src); @@ -293,11 +411,11 @@ BlockMatrixArray::vmult (BlockVector& dst, template -template +template inline void -BlockMatrixArray::Tvmult_add (BlockVector& dst, - const BlockVector& src) const +BlockMatrixArray::Tvmult_add (BlockVector& dst, + const BlockVector& src) const { Assert (dst.n_blocks() == block_cols, ExcDimensionMismatch(dst.n_blocks(), block_cols)); @@ -307,7 +425,7 @@ BlockMatrixArray::Tvmult_add (BlockVector& dst, typename vector::const_iterator m = entries.begin(); typename vector::const_iterator end = entries.end(); - static Vector aux; + static Vector aux; for (; m != end ; ++m) { @@ -335,11 +453,11 @@ BlockMatrixArray::Tvmult_add (BlockVector& dst, template -template +template inline void -BlockMatrixArray::Tvmult (BlockVector& dst, - const BlockVector& src) const +BlockMatrixArray::Tvmult (BlockVector& dst, + const BlockVector& src) const { dst = 0.; Tvmult_add (dst, src); @@ -394,6 +512,7 @@ BlockMatrixArray::print_latex (ostream& out) const pair (m->matrix, string("M"))); STRINGSTREAM stream; stream << matrix_number++; + stream << ends; x.first->second += stream.str(); } @@ -403,6 +522,7 @@ BlockMatrixArray::print_latex (ostream& out) const stream << matrix_names.find(m->matrix)->second; if (m->transpose) stream << "^T"; + stream << ends; array(m->row, m->col) += stream.str(); } for (unsigned int i=0;i::print_latex (ostream& out) const cout << "\\end{array}" << endl; } +//----------------------------------------------------------------------// + +template +BlockTrianglePrecondition::BlockTrianglePrecondition( + unsigned int block_rows, + bool backward) + : + BlockMatrixArray (block_rows, block_rows), + backward(backward) +{} + + +template +template +inline +void +BlockTrianglePrecondition::do_row ( + BlockVector& dst, + unsigned int row_num) const +{ + const typename BlockMatrixArray::Entry* diagonal = 0; + + typename vector::Entry>::const_iterator + m = entries.begin(); + typename vector::Entry>::const_iterator + end = entries.end(); + + static Vector aux; + aux.reinit(dst.block(row_num), true); + + for (; m != end ; ++m) + { + const unsigned int i=m->row; + if (i != row_num) + continue; + const unsigned int j=m->col; + if (((j > i) && !backward) || ((j < i) && backward)) + continue; + if (j == i) + { + Assert (diagonal == 0, ExcMultipleDiagonal(j)); + diagonal = m; + } else { + if (m->transpose) + m->matrix->Tvmult(aux, dst.block(j)); + else + m->matrix->vmult(aux, dst.block(j)); + dst.block(i).add (-1 * m->prefix, aux); + } + } + if (diagonal->transpose) + diagonal->matrix->Tvmult(aux, dst.block(row_num)); + else + diagonal->matrix->vmult(aux, dst.block(row_num)); + dst.block(row_num).equ (diagonal->prefix, aux); + +} + + + +template +template +inline +void +BlockTrianglePrecondition::vmult_add ( + BlockVector& dst, + const BlockVector& src) const +{ + Assert (dst.n_blocks() == n_block_rows(), + ExcDimensionMismatch(dst.n_blocks(), n_block_rows())); + Assert (src.n_blocks() == n_block_cols(), + ExcDimensionMismatch(src.n_blocks(), n_block_cols())); + + BlockVector& aux; + aux.reinit(dst); + vmult(aux, src); + dst.add(aux); +} + + + +template +template +inline +void +BlockTrianglePrecondition::vmult ( + BlockVector& dst, + const BlockVector& src) const +{ + Assert (dst.n_blocks() == n_block_rows(), + ExcDimensionMismatch(dst.n_blocks(), n_block_rows())); + Assert (src.n_blocks() == n_block_cols(), + ExcDimensionMismatch(src.n_blocks(), n_block_cols())); + + dst.equ(1., src); + + if (backward) + { + for (unsigned int i=n_block_rows(); i>0;) + do_row(dst, --i); + } else { + for (unsigned int i=0; i +template +inline +void +BlockTrianglePrecondition::Tvmult ( + BlockVector&, + const BlockVector&) const +{ + Assert (false, ExcNotImplemented()); +} + + +template +template +inline +void +BlockTrianglePrecondition::Tvmult_add ( + BlockVector&, + const BlockVector&) const +{ + Assert (false, ExcNotImplemented()); +} + + #endif