From a91021a79119bd8263c09ad54e93118dee89279b Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Mon, 18 Sep 2000 19:32:52 +0000 Subject: [PATCH] Relaxation parameter in PreconditionBlock git-svn-id: https://svn.dealii.org/trunk@3333 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/precondition_block.h | 131 +++++++++++--- .../lac/precondition_block.templates.h | 171 +++++++++++++----- deal.II/lac/source/precondition_block.cc | 46 +++++ 3 files changed, 279 insertions(+), 69 deletions(-) diff --git a/deal.II/lac/include/lac/precondition_block.h b/deal.II/lac/include/lac/precondition_block.h index 49c3778fd8..509bab60fa 100644 --- a/deal.II/lac/include/lac/precondition_block.h +++ b/deal.II/lac/include/lac/precondition_block.h @@ -96,9 +96,14 @@ class PreconditionBlock * second step, the inverses of * the diagonal blocks may be * computed. + * + * Additionally, a relaxation + * parameter for derived classes + * may be provided. */ void initialize (const SparseMatrix& A, - const unsigned int block_size); + const unsigned int block_size, + const double relaxation = 1.); /** * Deletes the inverse diagonal @@ -192,6 +197,11 @@ class PreconditionBlock * Exception */ DeclException0 (ExcMatrixNotSquare); + + /** + * Exception + */ + DeclException0 (ExcDiagonalsNotStored); protected: /** @@ -200,6 +210,12 @@ class PreconditionBlock */ const FullMatrix& inverse (unsigned int i) const; + /** + * Access to the diagonal + * blocks. + */ + const FullMatrix& diagonal (unsigned int i) const; + /** * Size of the blocks. Each * diagonal block is assumed to @@ -220,6 +236,18 @@ class PreconditionBlock * derived classes. */ SmartPointer > A; + /** + * Relaxation parameter to be + * used by derived classes. + */ + double relaxation; + + + /** + * Flag for storing the diagonal + * blocks of the matrix. + */ + bool store_diagonals; private: /** @@ -227,14 +255,22 @@ class PreconditionBlock * matrices of the diagonal * blocks matrices as * @p{FullMatrix} - * matrices. For BlockSOR as - * preconditioning using + * matrices. Using * @p{inverse_type=float} saves * memory in comparison with * @p{inverse_type=double}. */ - vector > _inverse; + 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. + */ + vector > var_diagonal; + /** * Flag for diagonal compression. * @see set_same_diagonal() @@ -302,29 +338,19 @@ class PreconditionBlockJacobi : public Subscriptor, /** * Block SOR preconditioning. * - * The diagonal blocks and the elements - * (of arbitray structure) below the diagonal blocks are used - * in the @p{operator ()} function of this class. + * The functions @p{vmult} and @p{Tvmult} execute a (transposed) + * block-SOR step, based on the blocks in @ref{PreconditionBlock}. The + * elements outside the diagonal blocks may be distributed + * arbitrarily. * * See @ref{PreconditionBlock} for requirements on the matrix. * @author Ralf Hartmann, Guido Kanschat, 1999, 2000 */ template class PreconditionBlockSOR : public Subscriptor, - private PreconditionBlock + protected PreconditionBlock { public: - /** - * Constructor. Takes the damping - * Parameter as - */ - PreconditionBlockSOR (const number omega=1.); - - /** - * Destructor. - */ - ~PreconditionBlockSOR(); - /** * Make initialization function * publicly available. @@ -346,11 +372,6 @@ class PreconditionBlockSOR : public Subscriptor, */ PreconditionBlock::invert_diagblocks; - /** - * Set the relaxation parameter. - */ - void set_omega(number omega); - /** * Execute block SOR * preconditioning. @@ -375,12 +396,68 @@ class PreconditionBlockSOR : public Subscriptor, */ template void Tvmult (Vector&, const Vector&) const; +}; - private: + +/** + * Block SSOR preconditioning. + * + * The functions @p{vmult} and @p{Tvmult} execute a block-SSOR step, + * based on the implementation in @ref{PreconditionBlockSOR}. This + * class requires storage of the diagonal blocks and their inverses. + * + * See @ref{PreconditionBlock} for requirements on the matrix. + * @author Ralf Hartmann, Guido Kanschat, 1999, 2000 + */ +template +class PreconditionBlockSSOR : private PreconditionBlockSOR +{ + public: /** - * Relaxation parameter. + * Constructor. + */ + PreconditionBlockSSOR (); + /** + * Make initialization function + * publicly available. + */ + PreconditionBlock::initialize; + + /** + * Make function of base class public again. + */ + PreconditionBlock::clear; + + /** + * Make function of base class public again. + */ + PreconditionBlock::set_same_diagonal; + + /** + * Make function of base class public again. */ - number omega; + PreconditionBlock::invert_diagblocks; + + /** + * Execute block SSOR + * preconditioning. + * + * This function will + * automatically use the inverse + * matrices if they exist, if not + * then BlockSOR will waste much + * time inverting the diagonal + * block matrices in each + * preconditioning step. + */ + template + void vmult (Vector&, const Vector&) const; + + /** + * Same as @ref{vmult} + */ + template + void Tvmult (Vector&, const Vector&) const; }; diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index 6d53a57a2a..42aafca092 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -26,6 +26,7 @@ template PreconditionBlock::PreconditionBlock (): blocksize(0), A(0), + store_diagonals(false), same_diagonal(false) {}; @@ -33,16 +34,20 @@ PreconditionBlock::PreconditionBlock (): template PreconditionBlock::~PreconditionBlock () { - if (_inverse.size()!=0) - _inverse.erase(_inverse.begin(), _inverse.end()); + 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 (_inverse.size()!=0) - _inverse.erase(_inverse.begin(), _inverse.end()); + 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()); blocksize = 0; same_diagonal = false; A = 0; @@ -51,7 +56,8 @@ void PreconditionBlock::clear () template void PreconditionBlock::initialize (const SparseMatrix &M, - unsigned int bsize) + unsigned int bsize, + const double relax) { clear(); Assert (M.m() == M.n(), ExcMatrixNotSquare()); @@ -59,6 +65,7 @@ void PreconditionBlock::initialize (const SparseMatrix0, ExcIndexRange(bsize, 1, M.m())); Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m())); blocksize=bsize; + relaxation = relax; } @@ -67,10 +74,24 @@ const FullMatrix& PreconditionBlock::inverse(unsigned int i) const { if (same_diagonal) - return _inverse[0]; + return var_inverse[0]; - Assert (i < _inverse.size(), ExcIndexRange(i,0,_inverse.size())); - return _inverse[i]; + 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]; } @@ -85,7 +106,7 @@ template void PreconditionBlock::set_same_diagonal() { - Assert(_inverse.size()==0, ExcInverseMatricesAlreadyExist()); + Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); same_diagonal = true; } @@ -94,7 +115,7 @@ template bool PreconditionBlock::inverses_ready() const { - return (_inverse.size() != 0); + return (var_inverse.size() != 0); } @@ -105,7 +126,7 @@ void PreconditionBlock::invert_diagblocks() Assert (blocksize!=0, ExcNotInitialized()); const SparseMatrix &M=*A; - Assert (_inverse.size()==0, ExcInverseMatricesAlreadyExist()); + Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); const unsigned int n_cells = M.m()/blocksize; @@ -117,13 +138,17 @@ void PreconditionBlock::invert_diagblocks() // Invert only the first block // This is a copy of the code in the // 'else' part, stripped of the outer loop - _inverse.resize(1); + if (store_diagonals) + var_diagonal.resize(1); + var_inverse.resize(1); for (unsigned int row_cell=0; row_cell::invert_diagblocks() // row, column are the global numbering // of the unkowns. - // set the @p{_inverse} array to the right + // set the @p{var_inverse} array to the right // size. we could do it like this: - // _inverse = vector<>(n_cells,FullMatrix<>()) + // var_inverse = vector<>(n_cells,FullMatrix<>()) // but this would involve copying many // FullMatrix objects. // // the following is a neat trick which // avoids copying + if (store_diagonals) + { + vector > tmp(n_cells, + FullMatrix(blocksize)); + var_diagonal.swap (tmp); + } + vector > tmp(n_cells, FullMatrix(blocksize)); - _inverse.swap (tmp); + var_inverse.swap (tmp); M_cell.clear (); @@ -156,7 +188,9 @@ void PreconditionBlock::invert_diagblocks() column_cell begin_diag_block+=blocksize; } + dst.scale(relaxation); } @@ -254,27 +289,6 @@ void PreconditionBlockJacobi /*--------------------- PreconditionBlockSOR -----------------------*/ -template -PreconditionBlockSOR::PreconditionBlockSOR(const number omega): - omega(omega) -{} - - -template -PreconditionBlockSOR::~PreconditionBlockSOR() -{} - - - -template -void -PreconditionBlockSOR::set_omega(number om) -{ - omega = om; -} - - - template template void PreconditionBlockSOR::vmult (Vector &dst, @@ -339,7 +353,7 @@ void PreconditionBlockSOR::vmult (Vector &ds M_cell.backward(x_cell,b_cell); // distribute x_cell to dst for (row=cell*blocksize, row_cell=0; row_cell::vmult (Vector &ds inverse(cell).vmult(x_cell, b_cell); // distribute x_cell to dst for (row=cell*blocksize, row_cell=0; row_cell::Tvmult (Vector &d M_cell.backward(x_cell,b_cell); // distribute x_cell to dst for (row=cell*blocksize, row_cell=0; row_cell::Tvmult (Vector &d inverse(cell).vmult(x_cell, b_cell); // distribute x_cell to dst for (row=cell*blocksize, row_cell=0; row_cell +PreconditionBlockSSOR::PreconditionBlockSSOR () +{ + store_diagonals = 1; +} + + +template +template +void PreconditionBlockSSOR::vmult (Vector &dst, + const Vector &src) const +{ + Vector help; + help.reinit(dst); + + PreconditionBlockSOR::vmult(help, src); + + Vector cell_src(blocksize); + Vector cell_dst(blocksize); + + const unsigned int n_cells = A->m()/blocksize; + + // Multiply with diagonal blocks + for (unsigned int cell=0; cell::Tvmult(dst, help); +} + +template +template +void PreconditionBlockSSOR::Tvmult (Vector &dst, + const Vector &src) const +{ + Vector help; + help.reinit(dst); + + PreconditionBlockSOR::Tvmult(help, src); + + Vector cell_src(blocksize); + Vector cell_dst(blocksize); + + const unsigned int n_cells = A->m()/blocksize; + + // Multiply with diagonal blocks + for (unsigned int cell=0; cell::vmult(dst, help); +} #endif diff --git a/deal.II/lac/source/precondition_block.cc b/deal.II/lac/source/precondition_block.cc index eed558300d..0427e542e8 100644 --- a/deal.II/lac/source/precondition_block.cc +++ b/deal.II/lac/source/precondition_block.cc @@ -112,3 +112,49 @@ template void PreconditionBlockSOR::Tvmult ( Vector &, const Vector &) const; +/*--------------------- PreconditionBlockSSOR -----------------------*/ + + +// explicit instantiations for "float" PreconditionBlock +template class PreconditionBlockSSOR; + +template void PreconditionBlockSSOR::vmult ( + Vector &, const Vector &) const; +template void PreconditionBlockSSOR::vmult ( + Vector &, const Vector &) const; +template void PreconditionBlockSSOR::Tvmult ( + Vector &, const Vector &) const; +template void PreconditionBlockSSOR::Tvmult ( + Vector &, const Vector &) const; + + +// the instantiation for class PreconditionBlockSSOR is skipped +// because it does not make sense to have inverse block matrices with +// higher precision than the matrix itself + + +// explicit instantiations for "double" PreconditionBlockSSOR +template class PreconditionBlockSSOR; + + +template void PreconditionBlockSSOR::vmult ( + Vector &, const Vector &) const; +template void PreconditionBlockSSOR::vmult ( + Vector &, const Vector &) const; +template void PreconditionBlockSSOR::Tvmult ( + Vector &, const Vector &) const; +template void PreconditionBlockSSOR::Tvmult ( + Vector &, const Vector &) const; + +template class PreconditionBlockSSOR; + +template void PreconditionBlockSSOR::vmult ( + Vector &, const Vector &) const; +template void PreconditionBlockSSOR::vmult ( + Vector &, const Vector &) const; +template void PreconditionBlockSSOR::Tvmult ( + Vector &, const Vector &) const; +template void PreconditionBlockSSOR::Tvmult ( + Vector &, const Vector &) const; + + -- 2.39.5