From: Guido Kanschat Date: Tue, 20 Jul 2010 05:54:46 +0000 (+0000) Subject: introduce a base class for PreconditionBlock which only stores matrices X-Git-Tag: v8.0.0~5788 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a637886715ea3eb636238a675b4f969377d7e22b;p=dealii.git introduce a base class for PreconditionBlock which only stores matrices git-svn-id: https://svn.dealii.org/trunk@21542 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/precondition_block.h b/deal.II/lac/include/lac/precondition_block.h index 4ac9f84287..100be2e1d8 100644 --- a/deal.II/lac/include/lac/precondition_block.h +++ b/deal.II/lac/include/lac/precondition_block.h @@ -18,14 +18,12 @@ #include #include #include +#include #include DEAL_II_NAMESPACE_OPEN -template class FullMatrix; -template class Vector; - template class PreconditionBlockJacobi; @@ -35,15 +33,14 @@ class PreconditionBlockJacobi; /** - * Base class for @p PreconditionBlockJacobi, @p PreconditionBlockSOR, - * ... This class assumes the MATRIX consisting of - * invertible blocks of @p blocksize on the diagonal and provides the - * inversion of the diagonal blocks of the matrix. NOT only block - * diagonal matrices are allowed but all matrices of arbitrary - * structure with the minimal property of having invertible blocks on - * the diagonal! Still the matrix must have access to single matrix - * entries. therefore, BlockMatrixArray is not a possible matrix - * class. + * Base class for actual block preconditioners. This class assumes the + * MATRIX consisting of invertible blocks of @p blocksize on + * the diagonal and provides the inversion of the diagonal blocks of + * the matrix. NOT only block diagonal matrices are allowed but all + * matrices of arbitrary structure with the minimal property of having + * invertible blocks on the diagonal! Still the matrix must have + * access to single matrix entries. therefore, BlockMatrixArray is not + * a possible matrix class. * * This block matrix structure is given e.g. for the DG method for the * transport equation. For a downstream numbering the matrices even @@ -72,10 +69,13 @@ class PreconditionBlockJacobi; * choice. * * @see @ref GlossBlockLA "Block (linear algebra)" - * @author Ralf Hartmann, Guido Kanschat, 1999, 2000 + * @author Ralf Hartmann, Guido Kanschat + * @date 1999, 2000, 2010 */ template -class PreconditionBlock : public virtual Subscriptor +class PreconditionBlock + : public virtual Subscriptor, + protected PreconditionBlockBase { private: /** @@ -132,7 +132,7 @@ class PreconditionBlock : public virtual Subscriptor /** * Constructor. */ - PreconditionBlock(); + PreconditionBlock(bool store_diagonals = false); /** * Destructor. @@ -225,27 +225,6 @@ class PreconditionBlock : public virtual Subscriptor value_type el(unsigned int i, unsigned int j) const; - /** - * Use only the inverse of the - * first diagonal block to save - * memory and computation time. - * - * Possible applications: - * computing on a cartesian grid, - * all diagonal blocks are the - * same or all diagonal blocks - * are at least similar and - * inversion of one of them still - * yields a preconditioner. - */ - void set_same_diagonal (); - - /** - * Does the matrix use only one - * diagonal block? - */ - bool same_diagonal () const; - /** * Stores the inverse of the * diagonal blocks in @@ -342,12 +321,6 @@ class PreconditionBlock : public virtual Subscriptor */ unsigned int memory_consumption () const; - /** - * Determine, whether inverses - * have been computed. - */ - bool inverses_ready () const; - /** @addtogroup Exceptions * @{ */ @@ -360,39 +333,13 @@ class PreconditionBlock : public virtual Subscriptor << " and the size of the matrix " << arg2 << " do not match."); - /** - * Exception - */ - DeclException2 (ExcWrongNumberOfInverses, - int, int, - << "There are " << arg1 - << " inverse matrices but " << arg2 - << " cells."); - /** * Exception */ DeclException0 (ExcInverseMatricesAlreadyExist); - /** - * Exception - */ - DeclException0 (ExcDiagonalsNotStored); - //@} - /** - * Access to the inverse diagonal - * blocks. - */ - const FullMatrix& inverse (unsigned int i) const; - - /** - * Access to the diagonal - * blocks. - */ - const FullMatrix& diagonal (unsigned int i) const; - protected: /** * Size of the blocks. Each @@ -419,14 +366,7 @@ class PreconditionBlock : public virtual Subscriptor * used by derived classes. */ double relaxation; - - - /** - * Flag for storing the diagonal - * blocks of the matrix. - */ - bool store_diagonals; - + /** * Number of blocks. */ @@ -446,30 +386,6 @@ class PreconditionBlock : public virtual Subscriptor * Flag for diagonal compression. * @ref set_same_diagonal() */ - private: - - /** - * Storage of the inverse - * matrices of the diagonal - * blocks matrices as - * FullMatrix - * matrices. Using - * inverse_type=float saves - * memory in comparison with - * inverse_type=double. - */ - std::vector > var_inverse; - - /** - * Storage of the original diagonal blocks. - * These are only filled if @p store_diagonals - * is @p true. - * - * Used by the blocked SSOR method. - */ - std::vector > var_diagonal; - - bool var_same_diagonal; }; @@ -808,6 +724,11 @@ class PreconditionBlockSOR : public virtual Subscriptor, protected PreconditionBlock { public: + /** + * Default constructor. + */ + PreconditionBlockSOR(); + /** * Define number type of matrix. */ @@ -944,6 +865,12 @@ class PreconditionBlockSOR : public virtual Subscriptor, void Tstep (Vector& dst, const Vector& rhs) const; protected: + /** + * Constructor to be used by + * PreconditionBlockSSOR. + */ + PreconditionBlockSOR(bool store); + /** * Implementation of the forward * substitution loop called by @@ -1138,7 +1065,7 @@ PreconditionBlock::el ( const unsigned int bs = blocksize; const unsigned int nb = i/bs; - const FullMatrix& B = inverse(nb); + const FullMatrix& B = this->inverse(nb); const unsigned int ib = i % bs; const unsigned int jb = j % bs; diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index e86b7e04c5..4ad3844a75 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -27,33 +27,23 @@ DEAL_II_NAMESPACE_OPEN template -PreconditionBlock::PreconditionBlock (): - blocksize(0), - A(0, typeid(*this).name()), - store_diagonals(false), - var_same_diagonal(false) +PreconditionBlock::PreconditionBlock (bool store) + : PreconditionBlockBase(store), + blocksize(0), + A(0, typeid(*this).name()) {} template PreconditionBlock::~PreconditionBlock () -{ - if (var_inverse.size()!=0) - var_inverse.erase(var_inverse.begin(), var_inverse.end()); - if (var_diagonal.size()!=0) - var_diagonal.erase(var_diagonal.begin(), var_diagonal.end()); -} +{} template void PreconditionBlock::clear () { - if (var_inverse.size()!=0) - var_inverse.erase(var_inverse.begin(), var_inverse.end()); - if (var_diagonal.size()!=0) - var_diagonal.erase(var_diagonal.begin(), var_diagonal.end()); + PreconditionBlockBase::clear(); blocksize = 0; - var_same_diagonal = false; A = 0; } @@ -72,13 +62,14 @@ void PreconditionBlock::initialize ( Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m())); blocksize=bsize; relaxation = parameters.relaxation; - var_same_diagonal = parameters.same_diagonal; nblocks = A->m()/bsize; + this->reinit(nblocks, blocksize, parameters.same_diagonal); if (parameters.invert_diagonal) invert_diagblocks(); } + template void PreconditionBlock::initialize ( const MATRIX &M, @@ -96,8 +87,8 @@ void PreconditionBlock::initialize ( Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m())); blocksize=bsize; relaxation = parameters.relaxation; - var_same_diagonal = parameters.same_diagonal; nblocks = A->m()/bsize; + this->reinit(nblocks, blocksize, parameters.same_diagonal); if (parameters.invert_diagonal) invert_permuted_diagblocks(permutation, inverse_permutation); @@ -113,21 +104,14 @@ void PreconditionBlock::invert_permuted_diagblocks( Assert (blocksize!=0, ExcNotInitialized()); const MATRIX &M=*A; - Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); + Assert (this->inverses_ready()==0, ExcInverseMatricesAlreadyExist()); FullMatrix M_cell(blocksize); - if (same_diagonal()) + if (this->same_diagonal()) { deallog << "PreconditionBlock uses only one diagonal block" << std::endl; - // Invert only the first block - // This is a copy of the code in the - // 'else' part, stripped of the outer loop - if (store_diagonals) - var_diagonal.resize(1); - var_inverse.resize(1); - var_inverse[0].reinit(blocksize, blocksize); - + for (unsigned int row_cell=0; row_cell::invert_permuted_diagblocks( ++entry; } } - if (store_diagonals) - var_diagonal[0] = M_cell; - var_inverse[0].invert(M_cell); + if (this->store_diagonals()) + this->diagonal(0) = M_cell; + this->inverse(0).invert(M_cell); } else { @@ -152,32 +136,6 @@ void PreconditionBlock::invert_permuted_diagblocks( // blocks. // row, column are the global numbering // of the unkowns. - - // set the var_inverse array to the right - // size. we could do it like this: - // var_inverse = vector<>(nblocks,FullMatrix<>()) - // but this would involve copying many - // FullMatrix objects. - // - // the following is a neat trick which - // avoids copying - if (store_diagonals) - { - std::vector > - tmp(nblocks, FullMatrix(blocksize)); - var_diagonal.swap (tmp); - } - - if (true) - { - std::vector > - tmp(nblocks, FullMatrix(blocksize)); - var_inverse.swap (tmp); - // make sure the tmp object - // goes out of scope as - // soon as possible - }; - M_cell = 0; for (unsigned int cell=0; cell::invert_permuted_diagblocks( M_cell(row_cell, column_cell) = entry->value(); } } - - if (store_diagonals) - var_diagonal[cell] = M_cell; - var_inverse[cell].invert(M_cell); + + if (this->store_diagonals()) + this->diagonal(cell) = M_cell; + this->inverse(cell).invert(M_cell); } } + this->inverses_computed(true); } @@ -405,32 +364,6 @@ void PreconditionBlock::backward_step ( } -template -const FullMatrix& -PreconditionBlock::inverse(unsigned int i) const -{ - if (same_diagonal()) - return var_inverse[0]; - - Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size())); - return var_inverse[i]; -} - - -template -const FullMatrix& -PreconditionBlock::diagonal(unsigned int i) const -{ - Assert(store_diagonals, ExcDiagonalsNotStored()); - - if (same_diagonal()) - return var_diagonal[0]; - - Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size())); - return var_diagonal[i]; -} - - template unsigned int PreconditionBlock::block_size() const { @@ -438,31 +371,6 @@ unsigned int PreconditionBlock::block_size() const } -template -void -PreconditionBlock::set_same_diagonal() -{ - Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); - var_same_diagonal = true; -} - - -template -bool -PreconditionBlock::same_diagonal() const -{ - return var_same_diagonal; -} - - -template -bool -PreconditionBlock::inverses_ready() const -{ - return (var_inverse.size() != 0); -} - - template void PreconditionBlock::invert_diagblocks() { @@ -470,21 +378,13 @@ void PreconditionBlock::invert_diagblocks() Assert (blocksize!=0, ExcNotInitialized()); const MATRIX &M=*A; - Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); - + Assert (this->inverses_ready()==0, ExcInverseMatricesAlreadyExist()); + FullMatrix M_cell(blocksize); - if (same_diagonal()) + if (this->same_diagonal()) { deallog << "PreconditionBlock uses only one diagonal block" << std::endl; - // Invert only the first block - // This is a copy of the code in the - // 'else' part, stripped of the outer loop - if (store_diagonals) - var_diagonal.resize(1); - var_inverse.resize(1); - var_inverse[0].reinit(blocksize, blocksize); - for (unsigned int row_cell=0; row_cell::invert_diagblocks() ++entry; } } - if (store_diagonals) - var_diagonal[0] = M_cell; - var_inverse[0].invert(M_cell); + if (this->store_diagonals()) + this->diagonal(0) = M_cell; + this->inverse(0).invert(M_cell); } else { - // cell_row, cell_column are the - // numbering of the blocks (cells). - // row_cell, column_cell are the local - // numbering of the unknowns in the - // blocks. - // row, column are the global numbering - // of the unkowns. - - // set the @p var_inverse array to the right - // size. we could do it like this: - // var_inverse = vector<>(nblocks,FullMatrix<>()) - // but this would involve copying many - // FullMatrix objects. - // - // the following is a neat trick which - // avoids copying - if (store_diagonals) - { - std::vector > - tmp(nblocks, FullMatrix(blocksize)); - var_diagonal.swap (tmp); - } - - if (true) - { - std::vector > - tmp(nblocks, FullMatrix(blocksize)); - var_inverse.swap (tmp); - // make sure the tmp object - // goes out of scope as - // soon as possible - }; - M_cell = 0; for (unsigned int cell=0; cell::invert_diagblocks() } } - if (store_diagonals) - var_diagonal[cell] = M_cell; - var_inverse[cell].invert(M_cell); + if (this->store_diagonals()) + this->diagonal(cell) = M_cell; + this->inverse(cell).invert(M_cell); } } + this->inverses_computed(true); } @@ -587,11 +455,8 @@ template unsigned int PreconditionBlock::memory_consumption () const { - unsigned int mem = sizeof(*this); - for (unsigned int i=0; i); + mem += PreconditionBlockBase::memory_consumption(); return mem; } @@ -768,6 +633,18 @@ void PreconditionBlockJacobi /*--------------------- PreconditionBlockSOR -----------------------*/ +template +PreconditionBlockSOR::PreconditionBlockSOR () + : PreconditionBlock (false) + +{} + +template +PreconditionBlockSOR::PreconditionBlockSOR (bool store) + : PreconditionBlock (store) + +{} + template template void PreconditionBlockSOR::forward ( @@ -1059,9 +936,9 @@ void PreconditionBlockSOR template PreconditionBlockSSOR::PreconditionBlockSSOR () -{ - this->store_diagonals = 1; -} + : PreconditionBlockSOR (true) + +{} template diff --git a/deal.II/lac/include/lac/precondition_block_base.h b/deal.II/lac/include/lac/precondition_block_base.h new file mode 100644 index 0000000000..e06473088a --- /dev/null +++ b/deal.II/lac/include/lac/precondition_block_base.h @@ -0,0 +1,412 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__precondition_block_base_h +#define __deal2__precondition_block_base_h + + +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +template class FullMatrix; +template class Vector; + +/** + * A class storing the inverse diagonal blocks for block + * preconditioners and block relaxation methods. + * + * This class does the book keeping for preconditioners and relaxation + * methods based on inverting blocks on the diagonal of a matrix. + * It allows us to either store all diagonal blocks and their + * inverses or the same block for each entry, and it keeps track of + * the choice. Thus, after initializing it and filling the inverse + * diagonal blocks correctly, a derived class can use inverse() with + * an integer argument referring to the block number. + * + * Additionally, it allows the storage of the original diagonal + * blocks, not only the inverses. These are for instance used in the + * intermediate step of the SSOR preconditioner. + * + * @author Guido Kanschat + * @date 2010 + */ +template +class PreconditionBlockBase +{ + public: + /** + * Constructor initializing + * default values. + */ + PreconditionBlockBase(bool store_diagonals); + + /** + * Deletes the inverse diagonal + * block matrices if existent hence + * leaves the class in the state + * that it had directly after + * calling the constructor. + */ + void clear(); + + /** + * Resize to this number of + * diagonal blocks with the given + * block size. If + * compress is true, + * then only one block will be + * stored. + */ + void reinit(unsigned int nblocks, unsigned int blocksize, bool compress); + + /** + * Tell the class that inverses + * are computed. + */ + void inverses_computed(bool are_they); + + /** + * Use only the inverse of the + * first diagonal block to save + * memory and computation time. + * + * Possible applications: + * computing on a cartesian grid, + * all diagonal blocks are the + * same or all diagonal blocks + * are at least similar and + * inversion of one of them still + * yields a preconditioner. + */ + void set_same_diagonal (); + + /** + * Does the matrix use only one + * diagonal block? + */ + bool same_diagonal () const; + + /** + * Check, whether diagonal blocks + * (not their inverses) + * should be stored. + */ + bool store_diagonals() const; + + /** + * Return true, if inverses are + * ready for use. + */ + bool inverses_ready () const; + + /** + * Checks whether the object is empty. + */ + bool empty () const; + + /** + * The number of blocks. + */ + unsigned int size() const; + + /** + * Read-only access to entries. + * This function is only possible + * if the inverse diagonal blocks + * are stored. + */ + number el(unsigned int i, unsigned int j) const; + + /** + * Access to the inverse diagonal + * blocks. + */ + FullMatrix& inverse (unsigned int i); + + /** + * Access to the inverse diagonal + * blocks. + */ + const FullMatrix& inverse (unsigned int i) const; + + /** + * Access to the diagonal + * blocks. + */ + FullMatrix& diagonal (unsigned int i); + + /** + * Access to the diagonal + * blocks. + */ + const FullMatrix& diagonal (unsigned int i) const; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. + */ + unsigned int memory_consumption () const; + + /** + * You are trying to access a + * diagonal block (not its + * inverse), but you decided not + * to store the diagonal blocks. + */ + DeclException0 (ExcDiagonalsNotStored); + + private: + /** + * The number of (inverse) + * diagonal blocks, if only one + * is stored. + */ + unsigned int n_diagonal_blocks; + + /** + * Storage of the inverse + * matrices of the diagonal + * blocks matrices as + * FullMatrix + * matrices. Using + * number=float saves + * memory in comparison with + * number=double. + */ + std::vector > var_inverse; + + /** + * Storage of the original diagonal blocks. + * + * Used by the blocked SSOR method. + */ + std::vector > var_diagonal; + + + /** + * This is true, if the field + * #var_diagonal is to be used. + */ + bool var_store_diagonals; + + /** + * This is true, if only one + * inverse is stored. + */ + bool var_same_diagonal; + + /** + * The inverse matrices are + * usable. Set by the parent + * class via inverses_computed(). + */ + bool var_inverses_ready; +}; + +//----------------------------------------------------------------------// + +template +inline +PreconditionBlockBase::PreconditionBlockBase(bool store) + : + n_diagonal_blocks(0), + var_store_diagonals(store), + var_same_diagonal(false), + var_inverses_ready(false) +{} + + +template +inline +void +PreconditionBlockBase::clear() +{ + if (var_inverse.size()!=0) + var_inverse.erase(var_inverse.begin(), var_inverse.end()); + if (var_diagonal.size()!=0) + var_diagonal.erase(var_diagonal.begin(), var_diagonal.end()); + var_same_diagonal = false; + var_inverses_ready = false; + n_diagonal_blocks = 0; +} + +template +inline +void +PreconditionBlockBase::reinit(unsigned int n, unsigned int b, bool compress) +{ + var_same_diagonal = compress; + var_inverses_ready = false; + n_diagonal_blocks = n; + + if (compress) + { + var_inverse.resize(1); + var_inverse[0].reinit(b,b); + if (store_diagonals()) + { + var_diagonal.resize(1); + var_diagonal[0].reinit(b,b); + } + } + else + { + // set the arrays to the right + // size. we could do it like this: + // var_inverse = vector<>(nblocks,FullMatrix<>()) + // but this would involve copying many + // FullMatrix objects. + // + // the following is a neat trick which + // avoids copying + if (store_diagonals()) + { + std::vector > + tmp(n, FullMatrix(b)); + var_diagonal.swap (tmp); + } + + if (true) + { + std::vector > + tmp(n, FullMatrix(b)); + var_inverse.swap (tmp); + // make sure the tmp object + // goes out of scope as + // soon as possible + } + } +} + + +template +inline +unsigned int +PreconditionBlockBase::size() const +{ + return n_diagonal_blocks; +} + + +template +inline +const FullMatrix& +PreconditionBlockBase::inverse(unsigned int i) const +{ + if (same_diagonal()) + return var_inverse[0]; + + Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size())); + return var_inverse[i]; +} + + +template +inline +const FullMatrix& +PreconditionBlockBase::diagonal(unsigned int i) const +{ + Assert(store_diagonals(), ExcDiagonalsNotStored()); + + if (same_diagonal()) + return var_diagonal[0]; + + Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size())); + return var_diagonal[i]; +} + + +template +inline +FullMatrix& +PreconditionBlockBase::inverse(unsigned int i) +{ + if (same_diagonal()) + return var_inverse[0]; + + Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size())); + return var_inverse[i]; +} + + +template +inline +FullMatrix& +PreconditionBlockBase::diagonal(unsigned int i) +{ + Assert(store_diagonals(), ExcDiagonalsNotStored()); + + if (same_diagonal()) + return var_diagonal[0]; + + Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size())); + return var_diagonal[i]; +} + + +template +inline bool +PreconditionBlockBase::same_diagonal() const +{ + return var_same_diagonal; +} + + +template +inline bool +PreconditionBlockBase::store_diagonals() const +{ + return var_store_diagonals; +} + + +template +inline void +PreconditionBlockBase::inverses_computed(bool x) +{ + var_inverses_ready = x; +} + + +template +inline bool +PreconditionBlockBase::inverses_ready() const +{ + return var_inverses_ready; +} + + +template +inline +unsigned int +PreconditionBlockBase::memory_consumption () const +{ + unsigned int mem = sizeof(*this); + for (unsigned int i=0; i