From 77da56935464177c028625d8750774124aa4ed12 Mon Sep 17 00:00:00 2001 From: guido Date: Tue, 13 Feb 2001 15:09:07 +0000 Subject: [PATCH] Block matrix based on an array of matrix pointers git-svn-id: https://svn.dealii.org/trunk@3923 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/2001/c-3-1.html | 9 + deal.II/lac/include/lac/block_matrix_array.h | 323 +++++++++++++++++++ 2 files changed, 332 insertions(+) create mode 100644 deal.II/lac/include/lac/block_matrix_array.h diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 57fc7a7e2a..d04cf26100 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -79,6 +79,15 @@ documentation, etc.

lac

    +

  1. New: class BlockMatrixArray<MATRIX> implements + a block matrix based on an array of matrix pointers. Since this + array may be incomplete, no direct access to entries is + allowed. Matrix blocks may be scaled and transposed. +
    + (GK 2001/02/13) +

    +

  2. New: There is now some support to include and use routines from the +#include +#include + +template class BlockVector; + + +/** + * Block matrix composed of different single matrices. + * + * Given a set of arbitrary matrices @p{A_i}, this class implements a + * block matrix with block entries of the form @p{M_{jk} = s_{jk}A_i}. + * Each @p{A_i} may be used several times with different prefix. + * + * Non-zero entries are registered by the function @p{enter}, zero + * entries are not stored at all. Using @p{enter} with the same + * location @p{(i,j)} several times will add the corresponding + * matrices in matrix-vector multiplications. + * + * @subsection{Requirements} + * + * The template argument @p{MATRIX} is a class providing the the + * matrix-vector multiplication functions @p{vmult} etc. defined in + * this class, but with arguments of type @p{VECTOR} instead of + * @p{BlockVector}. @ref{SparseMatrix} is a possible entry + * type. + * + * @author Guido Kanschat, 2000, 2001 */ +template +class BlockMatrixArray : + public Subscriptor +{ +public: + /** + * Constructor fixing the + * dimensions. + */ + BlockMatrixArray (const unsigned int n_block_rows, + const unsigned int n_block_cols); + + /** + * Add a block matrix entry. + */ + void enter (const MATRIX& matrix, + unsigned row, + unsigned int col, + double prefix = 1., + bool transpose = false); + + /** + * Number of block-entries per + * column. + */ + unsigned int n_block_rows () const; + + /** + * Number of block-entries per + * row. + */ + unsigned int n_block_cols () const; + + /** + * Matrix-vector multiplication. + */ + 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; + /** + * Transposed matrix-vector + * multiplication. + */ + 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; + +private: + /** + * Internal data structure. + * + * For each entry of a + * @p{BlockMatrixArray}, its + * position, matrix, prefix and + * optional transposition must be + * stored. This structure + * encapsulates all of them. + * + * @author Guido Kanschat, 2000, 2001 + */ + class Entry + { + public: + /** + * Constructor initializing all + * data fields. + */ + Entry (const MATRIX& matrix, + unsigned row, unsigned int col, + double prefix, bool transpose); + + /** + * Row number in the block + * matrix. + */ + unsigned int row; + + /** + * Column number in the block + * matrix. + */ + unsigned int col; + + /** + * Factor in front of the matrix + * block. + */ + double prefix; + + /** + * Indicates that matrix block + * must be transposed for + * multiplication. + */ + bool transpose; + + /** + * The matrix block itself. + */ + SmartPointer matrix; + }; + + /** + * Array of block entries in the + * matrix. + */ + vector entries; + + /** + * Number of blocks per column. + */ + unsigned int block_rows; + /** + * number of blocks per row. + */ + unsigned int block_cols; +}; + + +template +inline +BlockMatrixArray::Entry::Entry (const MATRIX& matrix, + unsigned row, unsigned int col, + double prefix, bool transpose) + : + row (row), + col (col), + prefix (prefix), + transpose (transpose), + matrix (&matrix) +{} + + + +template +inline +BlockMatrixArray::BlockMatrixArray (const unsigned int n_block_rows, + const unsigned int n_block_cols) + : block_rows (n_block_rows), + block_cols (n_block_cols) +{} + + + +template +inline +void +BlockMatrixArray::enter (const MATRIX& matrix, + unsigned row, unsigned int col, + double prefix, bool transpose) +{ + entries.push_back(Entry(matrix, row, col, prefix, transpose)); +} + + + +template +template +inline +void +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)); + + typename vector::const_iterator m = entries.begin(); + typename vector::const_iterator end = entries.end(); + + for (; m != end ; ++m) + { + Assert (m->prefix==1., ExcNotImplemented()); + if (m->transpose) + m->matrix->Tvmult_add(dst.block(m->row), + src.block(m->col)); + else + m->matrix->vmult_add(dst.block(m->row), + src.block(m->col)); + } +} + + + + +template +template +inline +void +BlockMatrixArray::vmult (BlockVector& dst, + const BlockVector& src) const +{ + dst.equ = 0.; + vmult_add (dst, src); +} + + + + +template +template +inline +void +BlockMatrixArray::Tvmult_add (BlockVector& dst, + const BlockVector& src) const +{ + Assert (dst.n_blocks() == block_cols, + ExcDimensionMismatch(dst.n_blocks(), block_cols)); + Assert (src.n_blocks() == block_rows, + ExcDimensionMismatch(src.n_blocks(), block_rows)); + + typename vector::const_iterator m = entries.begin(); + typename vector::const_iterator end = entries.end(); + + for (; m != end ; ++m) + { + Assert (m->prefix==1., ExcNotImplemented()); + if (m->transpose) + m->matrix->vmult_add(dst.block(m->col), + src.block(m->row)); + else + m->matrix->Tvmult_add(dst.block(m->col), + src.block(m->row)); + } +} + + + +template +template +inline +void +BlockMatrixArray::Tvmult (BlockVector& dst, + const BlockVector& src) const +{ + dst = 0.; + Tvmult_add (dst, src); +} + + + + +template +inline +unsigned int +BlockMatrixArray::n_block_rows () const +{ + return block_rows; +} + + + +template +inline +unsigned int +BlockMatrixArray::n_block_cols () const +{ + return block_cols; +} + + + + +#endif -- 2.39.5