From 791d6173901cd513017fdf1430c8747004c48aeb Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 22 Mar 2005 01:03:04 +0000 Subject: [PATCH] final? BlockMatrixArray git-svn-id: https://svn.dealii.org/trunk@10198 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.html | 13 +- deal.II/lac/include/lac/block_matrix_array.h | 434 ++++--------------- deal.II/lac/source/block_matrix_array.cc | 364 ++++++++++++++++ 3 files changed, 456 insertions(+), 355 deletions(-) create mode 100644 deal.II/lac/source/block_matrix_array.cc diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index 161eff2d8b..cd94ccb2f1 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -42,14 +42,17 @@ inconvenience this causes.
    -
  1. Changed: BlockMatrixArray - received a second template argument number, equal to the - same argument of the used

    Changed: The template argument of BlockMatrixArray was changed to + number, equal to the same argument of the used BlockVector. Furthermore, its constructor requires an additional argument of type VectorMemory<Vector<number> > - providing space for auxiliary vectors. -
    (GK 2005/03/14) + providing space for auxiliary vectors. Since the entries are now of + type PointerMatrixBase, even matrices + with blocks of different types can be constructed. +
    (GK 2005/03/21)

  2. Changed: The GeometryInfo:: #include -#include #include +#include #include #include #include #include +#include #ifdef HAVE_STD_STRINGSTREAM # include @@ -40,15 +41,16 @@ template class Vector; /** - * Block matrix composed of different single matrices. + * Block matrix composed of different single matrices; these matrices + * may even be of different types. * * Given a set of arbitrary matrices Ai, this class * implements a block matrix with block entries of the form * Mjk = sjkAi. Each * Ai may be used several times with different * prefix. The matrices are not copied into the BlockMatrixArray - * object, but rather references to them will be stored along with - * factors and transposition flags. + * object, but rather a PointerMatrix referencing each of them will be + * stored along with factors and transposition flags. * * Non-zero entries are registered by the function enter(), zero * entries are not stored at all. Using enter() with the same location @@ -69,7 +71,9 @@ template class Vector; * matrix-vector multiplication functions vmult(), Tvmult(), * vmult_add() and Tvmult_add() used in this class, but with arguments * of type Vector<number> instead of - * BlockVector<number>. SparseMatrix is a possible entry type. + * BlockVector<number>. Every matrix which can be used by + * PointerMatrix is allowed, in particular SparseMatrix is a possible + * entry type. * *

    Example program

    * We document the relevant parts of examples/doxygen/block_matrix_array.cc. @@ -118,7 +122,7 @@ template class Vector; * * @author Guido Kanschat, 2000 - 2005 */ -template +template class BlockMatrixArray : public Subscriptor { public: @@ -156,6 +160,7 @@ class BlockMatrixArray : public Subscriptor * ExcDimensionMismatch in one of * the multiplication functions. */ + template void enter (const MATRIX &matrix, const unsigned int row, const unsigned int col, @@ -263,10 +268,33 @@ class BlockMatrixArray : public Subscriptor * Constructor initializing all * data fields. */ + template Entry (const MATRIX& matrix, unsigned row, unsigned int col, double prefix, bool transpose); - + /** + * Copy constructor + * invalidating the old + * object. Since it is only + * used for entering + * temporary objects into a + * vector, this is ok. + * + * For a deep copy, we would + * need a reproduction + * operator in + * PointerMatixBase. + */ + Entry(const Entry&); + + /** + * Destructor, where we + * delete the PointerMatrix + * created by the + * constructor. + */ + ~Entry(); + /** * Row number in the block * matrix. @@ -295,7 +323,7 @@ class BlockMatrixArray : public Subscriptor /** * The matrix block itself. */ - SmartPointer matrix; + PointerMatrixBase >* matrix; }; /** @@ -382,9 +410,9 @@ class BlockMatrixArray : public Subscriptor * * @author Guido Kanschat, 2001, 2005 */ -template +template class BlockTrianglePrecondition - : private BlockMatrixArray + : private BlockMatrixArray { public: /** @@ -403,6 +431,15 @@ class BlockTrianglePrecondition */ void reinit(const unsigned int n_block_rows); + /** + * Add a new block. Calls BlockMatrixArray::enter(). + */ + template + void enter (const MATRIX &matrix, + const unsigned int row, + const unsigned int col, + const double prefix = 1., + const bool transpose = false); /** * Preconditioning. */ @@ -432,25 +469,25 @@ class BlockTrianglePrecondition /** * Make function of base class available. */ - BlockMatrixArray::enter; + BlockMatrixArray::enter; /** * Make function of base class available. */ - BlockMatrixArray::print_latex; + BlockMatrixArray::print_latex; /** * Make function of base class available. */ - BlockMatrixArray::n_block_rows; + BlockMatrixArray::n_block_rows; /** * Make function of base class available. */ - BlockMatrixArray::n_block_cols; - BlockMatrixArray::clear; - BlockMatrixArray::Subscriptor::subscribe; - BlockMatrixArray::Subscriptor::unsubscribe; + BlockMatrixArray::n_block_cols; + BlockMatrixArray::clear; + BlockMatrixArray::Subscriptor::subscribe; + BlockMatrixArray::Subscriptor::unsubscribe; /** @addtogroup Exceptions * @{ */ @@ -484,9 +521,10 @@ class BlockTrianglePrecondition ///@if NoDoc //----------------------------------------------------------------------// -template +template +template inline -BlockMatrixArray::Entry::Entry (const MATRIX& matrix, +BlockMatrixArray::Entry::Entry (const MATRIX& matrix, unsigned row, unsigned int col, double prefix, bool transpose) : @@ -494,40 +532,16 @@ BlockMatrixArray::Entry::Entry (const MATRIX& matrix, col (col), prefix (prefix), transpose (transpose), - matrix (&matrix, typeid(*this).name()) -{} - - - -template -inline -BlockMatrixArray::BlockMatrixArray ( - const unsigned int n_block_rows, - const unsigned int n_block_cols, - VectorMemory >& mem) - : block_rows (n_block_rows), - block_cols (n_block_cols), - mem(&mem, typeid(*this).name()) + matrix (new PointerMatrix >(&matrix)) {} -template +template +template inline void -BlockMatrixArray::reinit ( - const unsigned int n_block_rows, - const unsigned int n_block_cols) -{ - clear(); - block_rows = n_block_rows; - block_cols = n_block_cols; -} - -template -inline -void -BlockMatrixArray::enter (const MATRIX& matrix, +BlockMatrixArray::enter (const MATRIX& matrix, unsigned row, unsigned int col, double prefix, bool transpose) { @@ -537,189 +551,31 @@ BlockMatrixArray::enter (const MATRIX& matrix, } - -template -inline -void -BlockMatrixArray::clear () -{ - entries.clear(); -} - - -template -inline -void -BlockMatrixArray::vmult_add (BlockVector& dst, - const BlockVector& src) const +template +struct BlockMatrixArrayPointerMatrixLess { - Assert (dst.n_blocks() == block_rows, - ExcDimensionMismatch(dst.n_blocks(), block_rows)); - Assert (src.n_blocks() == block_cols, - ExcDimensionMismatch(src.n_blocks(), block_cols)); - - Vector* p_aux = mem->alloc(); - Vector& aux = *p_aux; - - typename std::vector::const_iterator m = entries.begin(); - typename std::vector::const_iterator end = entries.end(); - - for (; m != end ; ++m) - { - aux.reinit(dst.block(m->row)); - if (m->transpose) - m->matrix->Tvmult(aux, src.block(m->col)); - else - m->matrix->vmult(aux, src.block(m->col)); - dst.block(m->row).add (m->prefix, aux); - } - mem->free(p_aux); -} - - - - -template -inline -void -BlockMatrixArray::vmult (BlockVector& dst, - const BlockVector& src) const -{ - dst = 0.; - vmult_add (dst, src); -} - - - - -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 std::vector::const_iterator m = entries.begin(); - typename std::vector::const_iterator end = entries.end(); - - Vector* p_aux = mem->alloc(); - Vector& aux = *p_aux; - - for (; m != end ; ++m) - { - aux.reinit(dst.block(m->col)); - if (m->transpose) - m->matrix->vmult(aux, src.block(m->row)); - else - m->matrix->Tvmult(aux, src.block(m->row)); - dst.block(m->col).add (m->prefix, aux); - } - mem->free(p_aux); -} - - - -template -inline -void -BlockMatrixArray::Tvmult (BlockVector& dst, - const BlockVector& src) const -{ - dst = 0.; - Tvmult_add (dst, src); -} - - - - -template -inline -number -BlockMatrixArray::matrix_scalar_product ( - const BlockVector& u, - const BlockVector& v) const -{ - Assert (u.n_blocks() == block_rows, - ExcDimensionMismatch(u.n_blocks(), block_rows)); - Assert (v.n_blocks() == block_cols, - ExcDimensionMismatch(v.n_blocks(), block_cols)); - - Vector* p_aux = mem->alloc(); - Vector& aux = *p_aux; - - typename std::vector::const_iterator m; - typename std::vector::const_iterator end = entries.end(); - - number result = 0.; + bool operator () (const PointerMatrixBase >*a, + const PointerMatrixBase >* b) const + { + return *a < *b; + } +}; - for (unsigned int i=0;irow != i) - continue; - if (m->transpose) - m->matrix->Tvmult_add(aux, v.block(m->col)); - else - m->matrix->vmult(aux, v.block(m->col)); - } - result += u.block(i)*aux; - } - mem->free(p_aux); - - return result; -} - - - -template -inline -number -BlockMatrixArray::matrix_norm_square ( - const BlockVector& u) const -{ - return matrix_scalar_product(u,u); -} - - - -template -inline -unsigned int -BlockMatrixArray::n_block_rows () const -{ - return block_rows; -} - - - -template -inline -unsigned int -BlockMatrixArray::n_block_cols () const -{ - return block_cols; -} - - -template +template template inline void -BlockMatrixArray::print_latex (STREAM& out) const +BlockMatrixArray::print_latex (STREAM& out) const { out << "\\begin{array}{" << std::string(n_block_cols(), 'c') << "}" << std::endl; Table<2,std::string> array(n_block_rows(), n_block_cols()); - typedef std::map NameMap; + + typedef std::map >*, + std::string, BlockMatrixArrayPointerMatrixLess > NameMap; NameMap matrix_names; typename std::vector::const_iterator m = entries.begin(); @@ -732,7 +588,7 @@ BlockMatrixArray::print_latex (STREAM& out) const { std::pair x = matrix_names.insert( - std::pair (m->matrix, + std::pair >*, std::string> (m->matrix, std::string("M"))); #ifdef HAVE_STD_STRINGSTREAM std::ostringstream stream; @@ -779,144 +635,22 @@ BlockMatrixArray::print_latex (STREAM& out) const out << "\\end{array}" << std::endl; } -//----------------------------------------------------------------------// - -template -BlockTrianglePrecondition::BlockTrianglePrecondition( - unsigned int block_rows, - VectorMemory >& mem, - bool backward) - : - BlockMatrixArray (block_rows, block_rows, mem), - backward(backward) -{} - - -template -inline -void -BlockTrianglePrecondition::reinit ( - const unsigned int n) -{ - BlockMatrixArray::reinit(n,n); -} - -template -inline -void -BlockTrianglePrecondition::do_row ( - BlockVector& dst, - unsigned int row_num) const -{ - typename std::vector::Entry>::const_iterator - m = this->entries.begin(); - typename std::vector::Entry>::const_iterator - end = this->entries.end(); - typename std::vector::Entry>::const_iterator - diagonal = end; - - Vector* p_aux = this->mem->alloc(); - Vector& aux = *p_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 == end, 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); - - this->mem->free(p_aux); -} - - - -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 -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 +template inline void -BlockTrianglePrecondition::Tvmult ( - BlockVector&, - const BlockVector&) const +BlockTrianglePrecondition::enter ( + const MATRIX &matrix, + const unsigned int row, + const unsigned int col, + const double prefix, + const bool transpose) { - Assert (false, ExcNotImplemented()); + BlockMatrixArray::enter(matrix, row, col, prefix, transpose); } -template -inline -void -BlockTrianglePrecondition::Tvmult_add ( - BlockVector&, - const BlockVector&) const -{ - Assert (false, ExcNotImplemented()); -} - ///@endif #endif diff --git a/deal.II/lac/source/block_matrix_array.cc b/deal.II/lac/source/block_matrix_array.cc new file mode 100644 index 0000000000..4295cfbe08 --- /dev/null +++ b/deal.II/lac/source/block_matrix_array.cc @@ -0,0 +1,364 @@ +//------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998 - 2005 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//------------------------------------------------------------------- + +#include +#include +#include + + +template +BlockMatrixArray::Entry::Entry (const Entry& e) + : + row(e.row), + col(e.col), + prefix(e.prefix), + transpose(e.transpose), + matrix(e.matrix) +{ + Entry& e2 = const_cast(e); + e2.matrix = 0; +} + + + +template +BlockMatrixArray::Entry::~Entry () +{ + if (matrix) + delete matrix; +} + + + +template +BlockMatrixArray::BlockMatrixArray ( + const unsigned int n_block_rows, + const unsigned int n_block_cols, + VectorMemory >& mem) + : block_rows (n_block_rows), + block_cols (n_block_cols), + mem(&mem, typeid(*this).name()) +{} + + + +template +void +BlockMatrixArray::reinit ( + const unsigned int n_block_rows, + const unsigned int n_block_cols) +{ + clear(); + block_rows = n_block_rows; + block_cols = n_block_cols; +} + + + +template +void +BlockMatrixArray::clear () +{ + entries.clear(); +} + + +template +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)); + + Vector* p_aux = mem->alloc(); + Vector& aux = *p_aux; + + typename std::vector::const_iterator m = entries.begin(); + typename std::vector::const_iterator end = entries.end(); + + for (; m != end ; ++m) + { + aux.reinit(dst.block(m->row)); + if (m->transpose) + m->matrix->Tvmult(aux, src.block(m->col)); + else + m->matrix->vmult(aux, src.block(m->col)); + dst.block(m->row).add (m->prefix, aux); + } + mem->free(p_aux); +} + + + + +template +void +BlockMatrixArray::vmult (BlockVector& dst, + const BlockVector& src) const +{ + dst = 0.; + vmult_add (dst, src); +} + + + + +template +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 std::vector::const_iterator m = entries.begin(); + typename std::vector::const_iterator end = entries.end(); + + Vector* p_aux = mem->alloc(); + Vector& aux = *p_aux; + + for (; m != end ; ++m) + { + aux.reinit(dst.block(m->col)); + if (m->transpose) + m->matrix->vmult(aux, src.block(m->row)); + else + m->matrix->Tvmult(aux, src.block(m->row)); + dst.block(m->col).add (m->prefix, aux); + } + mem->free(p_aux); +} + + + +template +void +BlockMatrixArray::Tvmult (BlockVector& dst, + const BlockVector& src) const +{ + dst = 0.; + Tvmult_add (dst, src); +} + + + + +template +number +BlockMatrixArray::matrix_scalar_product ( + const BlockVector& u, + const BlockVector& v) const +{ + Assert (u.n_blocks() == block_rows, + ExcDimensionMismatch(u.n_blocks(), block_rows)); + Assert (v.n_blocks() == block_cols, + ExcDimensionMismatch(v.n_blocks(), block_cols)); + + Vector* p_aux = mem->alloc(); + Vector& aux = *p_aux; + + typename std::vector::const_iterator m; + typename std::vector::const_iterator end = entries.end(); + + number result = 0.; + + for (unsigned int i=0;irow != i) + continue; + if (m->transpose) + m->matrix->Tvmult_add(aux, v.block(m->col)); + else + m->matrix->vmult(aux, v.block(m->col)); + } + result += u.block(i)*aux; + } + mem->free(p_aux); + + return result; +} + + + +template +number +BlockMatrixArray::matrix_norm_square ( + const BlockVector& u) const +{ + return matrix_scalar_product(u,u); +} + + + +template +unsigned int +BlockMatrixArray::n_block_rows () const +{ + return block_rows; +} + + + +template +unsigned int +BlockMatrixArray::n_block_cols () const +{ + return block_cols; +} + + + +//----------------------------------------------------------------------// + +template +BlockTrianglePrecondition::BlockTrianglePrecondition( + unsigned int block_rows, + VectorMemory >& mem, + bool backward) + : + BlockMatrixArray (block_rows, block_rows, mem), + backward(backward) +{} + + +template +void +BlockTrianglePrecondition::reinit ( + const unsigned int n) +{ + BlockMatrixArray::reinit(n,n); +} + +template +void +BlockTrianglePrecondition::do_row ( + BlockVector& dst, + unsigned int row_num) const +{ + typename std::vector::Entry>::const_iterator + m = this->entries.begin(); + typename std::vector::Entry>::const_iterator + end = this->entries.end(); + typename std::vector::Entry>::const_iterator + diagonal = end; + + Vector* p_aux = this->mem->alloc(); + Vector& aux = *p_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 == end, 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); + + this->mem->free(p_aux); +} + + + +template +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 +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 +void +BlockTrianglePrecondition::Tvmult ( + BlockVector&, + const BlockVector&) const +{ + Assert (false, ExcNotImplemented()); +} + + +template +void +BlockTrianglePrecondition::Tvmult_add ( + BlockVector&, + const BlockVector&) const +{ + Assert (false, ExcNotImplemented()); +} + +template class BlockMatrixArray; +template class BlockMatrixArray; +template class BlockTrianglePrecondition; +template class BlockTrianglePrecondition; -- 2.39.5