From 9489fd2ee9b5b062282f99dd135033bfb4e22a58 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 15 Mar 2005 20:59:44 +0000 Subject: [PATCH] matrix scalar product and matrix norm, VectorMemory git-svn-id: https://svn.dealii.org/trunk@10161 0785d39b-7218-0410-832d-ea1e28bc413d --- .../examples/doxygen/block_matrix_array.cc | 11 +- deal.II/lac/include/lac/block_matrix_array.h | 236 ++++++++++++------ 2 files changed, 166 insertions(+), 81 deletions(-) diff --git a/deal.II/examples/doxygen/block_matrix_array.cc b/deal.II/examples/doxygen/block_matrix_array.cc index a69b774a3f..af584a033b 100644 --- a/deal.II/examples/doxygen/block_matrix_array.cc +++ b/deal.II/examples/doxygen/block_matrix_array.cc @@ -66,7 +66,9 @@ int main () B2.fill(B2data); C.fill(Cdata); - BlockMatrixArray > matrix(2,2); + GrowingVectorMemory > simple_mem; + + BlockMatrixArray, double> matrix(2, 2, simple_mem); matrix.enter(A,0,0,2.); matrix.enter(B1,0,1,-1.); @@ -97,12 +99,17 @@ int main () x.add(-1., result); deallog << "Error " << x.l2_norm() << std::endl; + deallog << "Error A-norm " + << std::sqrt(matrix.matrix_norm_square(x)) + << std::endl; + FullMatrix Ainv(4,4); Ainv.invert(A); FullMatrix Cinv(2,2); Cinv.invert(C); - BlockTrianglePrecondition > precondition(2,2); + BlockTrianglePrecondition, double> + precondition(2, simple_mem); precondition.enter(Ainv,0,0,.5); precondition.enter(Cinv,1,1); diff --git a/deal.II/lac/include/lac/block_matrix_array.h b/deal.II/lac/include/lac/block_matrix_array.h index 7faf15845c..8c0112183e 100644 --- a/deal.II/lac/include/lac/block_matrix_array.h +++ b/deal.II/lac/include/lac/block_matrix_array.h @@ -18,6 +18,8 @@ #include #include +#include + #include #include #include @@ -28,7 +30,6 @@ # include #endif -//TODO:[GK] Obtain aux vector from VectorMemory template class BlockVector; template class Vector; @@ -65,11 +66,10 @@ template class Vector; *

Requirements on MATRIX

* * The template argument MATRIX is a class providing the - * matrix-vector multiplication functions vmult, - * Tvmult, vmult_add and Tvmult_add used in - * this class, but with arguments of type VECTOR instead of - * BlockVector. SparseMatrix is a possible entry - * type. + * 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. * *

Example program

* We document the relevant parts of examples/doxygen/block_matrix_array.cc. @@ -84,6 +84,15 @@ template class Vector; * @skip main * @until C.fill * + * The BlockMatrixArray needs a VectorMemory<Vector<number> + * > object to allocate auxiliary vectors. number must + * equal the second template argument of BlockMatrixArray and also the + * number type of the BlockVector used later. We use the + * GrowingVectorMemory type, since it remembers the vector and avoids + * reallocating. + * + * @ line Growing + * * Now, we are ready to build a 2x2 BlockMatrixArray. * @line Block * First, we enter the matrix A multiplied by 2 in the upper left block @@ -109,7 +118,7 @@ template class Vector; * * @author Guido Kanschat, 2000 - 2005 */ -template +template class BlockMatrixArray : public Subscriptor { public: @@ -118,7 +127,8 @@ class BlockMatrixArray : public Subscriptor * dimensions. */ BlockMatrixArray (const unsigned int n_block_rows, - const unsigned int n_block_cols); + const unsigned int n_block_cols, + VectorMemory >& mem); /** * Add a block matrix entry. The @@ -166,7 +176,6 @@ class BlockMatrixArray : public Subscriptor /** * Matrix-vector multiplication. */ - template void vmult (BlockVector& dst, const BlockVector& src) const; @@ -174,7 +183,6 @@ class BlockMatrixArray : public Subscriptor * Matrix-vector multiplication * adding to dst. */ - template void vmult_add (BlockVector& dst, const BlockVector& src) const; @@ -182,7 +190,6 @@ class BlockMatrixArray : public Subscriptor * Transposed matrix-vector * multiplication. */ - template void Tvmult (BlockVector& dst, const BlockVector& src) const; @@ -191,10 +198,24 @@ class BlockMatrixArray : public Subscriptor * multiplication adding to * dst. */ - template void Tvmult_add (BlockVector& dst, const BlockVector& src) const; + /** + * Matrix scalar product between + * two vectors (at least for a + * symmetric matrix). + */ + number matrix_scalar_product (const BlockVector& u, + const BlockVector& v) const; + + /** + * Square of the matrix norm of a + * vector (at least for a + * symmetric matrix). + */ + number matrix_norm_square (const BlockVector& u) const; + /** * Print the block structure as a * LaTeX-array. This output will @@ -285,6 +306,12 @@ class BlockMatrixArray : public Subscriptor * number of blocks per row. */ unsigned int block_cols; + protected: + /** + * The memory used for auxiliary + * vectors. + */ + mutable SmartPointer > > mem; }; /*@}*/ @@ -348,8 +375,9 @@ class BlockMatrixArray : public Subscriptor * * @author Guido Kanschat, 2001, 2005 */ -template -class BlockTrianglePrecondition : private BlockMatrixArray +template +class BlockTrianglePrecondition + : private BlockMatrixArray { public: /** @@ -359,13 +387,13 @@ class BlockTrianglePrecondition : private BlockMatrixArray * insertion instead of forward. */ BlockTrianglePrecondition(unsigned int block_rows, + VectorMemory >& mem, bool backward = false); /** * Preconditioning. */ - template void vmult (BlockVector& dst, const BlockVector& src) const; @@ -373,14 +401,12 @@ class BlockTrianglePrecondition : private BlockMatrixArray * Preconditioning * adding to dst. */ - template void vmult_add (BlockVector& dst, const BlockVector& src) const; /** * Transposed preconditioning */ - template void Tvmult (BlockVector& dst, const BlockVector& src) const; @@ -388,32 +414,31 @@ class BlockTrianglePrecondition : private BlockMatrixArray * Transposed preconditioning * adding to dst. */ - template void Tvmult_add (BlockVector& dst, const BlockVector& src) const; /** * 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 * @{ */ @@ -433,7 +458,6 @@ class BlockTrianglePrecondition : private BlockMatrixArray * entry of the diagonal element * for one row. */ - template void do_row (BlockVector& dst, unsigned int row_num) const; @@ -448,9 +472,9 @@ class BlockTrianglePrecondition : private BlockMatrixArray ///@if NoDoc //----------------------------------------------------------------------// -template +template inline -BlockMatrixArray::Entry::Entry (const MATRIX& matrix, +BlockMatrixArray::Entry::Entry (const MATRIX& matrix, unsigned row, unsigned int col, double prefix, bool transpose) : @@ -463,20 +487,23 @@ BlockMatrixArray::Entry::Entry (const MATRIX& matrix, -template +template inline -BlockMatrixArray::BlockMatrixArray (const unsigned int n_block_rows, - const unsigned int n_block_cols) +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) + block_cols (n_block_cols), + mem(&mem) {} -template +template inline void -BlockMatrixArray::enter (const MATRIX& matrix, +BlockMatrixArray::enter (const MATRIX& matrix, unsigned row, unsigned int col, double prefix, bool transpose) { @@ -487,20 +514,19 @@ BlockMatrixArray::enter (const MATRIX& matrix, -template +template inline void -BlockMatrixArray::clear () +BlockMatrixArray::clear () { entries.clear(); } -template -template +template inline void -BlockMatrixArray::vmult_add (BlockVector& dst, +BlockMatrixArray::vmult_add (BlockVector& dst, const BlockVector& src) const { Assert (dst.n_blocks() == block_rows, @@ -508,7 +534,8 @@ BlockMatrixArray::vmult_add (BlockVector& dst, Assert (src.n_blocks() == block_cols, ExcDimensionMismatch(src.n_blocks(), block_cols)); - static Vector aux; + 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(); @@ -522,16 +549,16 @@ BlockMatrixArray::vmult_add (BlockVector& dst, m->matrix->vmult(aux, src.block(m->col)); dst.block(m->row).add (m->prefix, aux); } + mem->free(p_aux); } -template -template +template inline void -BlockMatrixArray::vmult (BlockVector& dst, +BlockMatrixArray::vmult (BlockVector& dst, const BlockVector& src) const { dst = 0.; @@ -541,11 +568,10 @@ BlockMatrixArray::vmult (BlockVector& dst, -template -template +template inline void -BlockMatrixArray::Tvmult_add (BlockVector& dst, +BlockMatrixArray::Tvmult_add (BlockVector& dst, const BlockVector& src) const { Assert (dst.n_blocks() == block_cols, @@ -556,7 +582,8 @@ BlockMatrixArray::Tvmult_add (BlockVector& dst, typename std::vector::const_iterator m = entries.begin(); typename std::vector::const_iterator end = entries.end(); - static Vector aux; + Vector* p_aux = mem->alloc(); + Vector& aux = *p_aux; for (; m != end ; ++m) { @@ -567,15 +594,15 @@ BlockMatrixArray::Tvmult_add (BlockVector& dst, m->matrix->Tvmult(aux, src.block(m->row)); dst.block(m->col).add (m->prefix, aux); } + mem->free(p_aux); } -template -template +template inline void -BlockMatrixArray::Tvmult (BlockVector& dst, +BlockMatrixArray::Tvmult (BlockVector& dst, const BlockVector& src) const { dst = 0.; @@ -585,31 +612,83 @@ BlockMatrixArray::Tvmult (BlockVector& dst, -template +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.; + + 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 +BlockMatrixArray::n_block_rows () const { return block_rows; } -template +template inline unsigned int -BlockMatrixArray::n_block_cols () const +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') @@ -678,32 +757,34 @@ BlockMatrixArray::print_latex (STREAM& out) const //----------------------------------------------------------------------// -template -BlockTrianglePrecondition::BlockTrianglePrecondition( +template +BlockTrianglePrecondition::BlockTrianglePrecondition( unsigned int block_rows, + VectorMemory >& mem, bool backward) : - BlockMatrixArray (block_rows, block_rows), - backward(backward) + BlockMatrixArray (block_rows, block_rows, mem), + backward(backward) {} -template -template +template inline void -BlockTrianglePrecondition::do_row ( +BlockTrianglePrecondition::do_row ( BlockVector& dst, unsigned int row_num) const { - typename std::vector::Entry>::const_iterator + typename std::vector::Entry>::const_iterator m = this->entries.begin(); - typename std::vector::Entry>::const_iterator + typename std::vector::Entry>::const_iterator end = this->entries.end(); - typename std::vector::Entry>::const_iterator + typename std::vector::Entry>::const_iterator diagonal = end; - static Vector aux; + Vector* p_aux = mem->alloc(); + Vector& aux = *p_aux; + aux.reinit(dst.block(row_num), true); for (; m != end ; ++m) @@ -732,15 +813,15 @@ BlockTrianglePrecondition::do_row ( diagonal->matrix->vmult(aux, dst.block(row_num)); dst.block(row_num).equ (diagonal->prefix, aux); + mem->free(p_aux); } -template -template +template inline void -BlockTrianglePrecondition::vmult_add ( +BlockTrianglePrecondition::vmult_add ( BlockVector& dst, const BlockVector& src) const { @@ -757,11 +838,10 @@ BlockTrianglePrecondition::vmult_add ( -template -template +template inline void -BlockTrianglePrecondition::vmult ( +BlockTrianglePrecondition::vmult ( BlockVector& dst, const BlockVector& src) const { @@ -783,11 +863,10 @@ BlockTrianglePrecondition::vmult ( } -template -template +template inline void -BlockTrianglePrecondition::Tvmult ( +BlockTrianglePrecondition::Tvmult ( BlockVector&, const BlockVector&) const { @@ -795,11 +874,10 @@ BlockTrianglePrecondition::Tvmult ( } -template -template +template inline void -BlockTrianglePrecondition::Tvmult_add ( +BlockTrianglePrecondition::Tvmult_add ( BlockVector&, const BlockVector&) const { -- 2.39.5