From c6af1369bd510373f32e442d2e88e60ec5887a53 Mon Sep 17 00:00:00 2001 From: kanschat Date: Tue, 26 Oct 2010 19:56:35 +0000 Subject: [PATCH] allow for QR method in PreconditionBlockBase git-svn-id: https://svn.dealii.org/trunk@22508 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/lac/householder.h | 33 ++- .../include/deal.II/lac/precondition_block.h | 6 +- .../lac/precondition_block.templates.h | 44 +++- .../deal.II/lac/precondition_block_base.h | 204 +++++++++++++++--- .../include/deal.II/lac/relaxation_block.h | 5 + .../deal.II/lac/relaxation_block.templates.h | 25 ++- 6 files changed, 270 insertions(+), 47 deletions(-) diff --git a/deal.II/include/deal.II/lac/householder.h b/deal.II/include/deal.II/lac/householder.h index c9a90496d6..e41be1b9db 100644 --- a/deal.II/include/deal.II/lac/householder.h +++ b/deal.II/include/deal.II/lac/householder.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 2005, 2006, 2007, 2008, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -101,6 +101,18 @@ class Householder : private FullMatrix double least_squares (BlockVector &dst, const BlockVector &src) const; + /** + * A wrapper to least_squares(), + * implementing the standard + * MATRIX interface. + */ + template + void vmult (VECTOR &dst, const VECTOR &src) const; + + template + void Tvmult (VECTOR &dst, const VECTOR &src) const; + + private: /** * Storage for the diagonal @@ -280,6 +292,25 @@ Householder::least_squares (BlockVector& dst, } +template +template +void +Householder::vmult (VECTOR &dst, const VECTOR &src) const +{ + least_squares (dst, src); +} + + +template +template +void +Householder::Tvmult (VECTOR &, const VECTOR &) const +{ + Assert(false, ExcNotImplemented()); +} + + + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/include/deal.II/lac/precondition_block.h b/deal.II/include/deal.II/lac/precondition_block.h index 631e718bc1..5d75ce1ae3 100644 --- a/deal.II/include/deal.II/lac/precondition_block.h +++ b/deal.II/include/deal.II/lac/precondition_block.h @@ -325,7 +325,11 @@ class PreconditionBlock * @{ */ /** - * Exception + * For non-overlapping block + * preconditioners, the block + * size must divide the matrix + * size. If not, this exception + * gets thrown. */ DeclException2 (ExcWrongBlockSize, int, int, diff --git a/deal.II/include/deal.II/lac/precondition_block.templates.h b/deal.II/include/deal.II/lac/precondition_block.templates.h index 0fb19a7c70..76eb4cacd3 100644 --- a/deal.II/include/deal.II/lac/precondition_block.templates.h +++ b/deal.II/include/deal.II/lac/precondition_block.templates.h @@ -130,7 +130,18 @@ void PreconditionBlock::invert_permuted_diagblocks( } if (this->store_diagonals()) this->diagonal(0) = M_cell; - this->inverse(0).invert(M_cell); + switch (this->inversion) + { + case PreconditionBlockBase::gauss_jordan: + this->inverse(0).invert(M_cell); + break; + case PreconditionBlockBase::householder: + this->inverse_householder(0).initialize(M_cell); + break; + default: + Assert(false, ExcNotImplemented()); + + } } else { @@ -170,7 +181,18 @@ void PreconditionBlock::invert_permuted_diagblocks( if (this->store_diagonals()) this->diagonal(cell) = M_cell; - this->inverse(cell).invert(M_cell); + switch (this->inversion) + { + case PreconditionBlockBase::gauss_jordan: + this->inverse(cell).invert(M_cell); + break; + case PreconditionBlockBase::householder: + this->inverse_householder(cell).initialize(M_cell); + break; + default: + Assert(false, ExcNotImplemented()); + + } } } this->inverses_computed(true); @@ -250,9 +272,9 @@ void PreconditionBlock::forward_step ( if (this->inverses_ready()) { if (transpose_diagonal) - this->inverse(cell).Tvmult(x_cell, b_cell); + this->inverse_Tvmult(cell, x_cell, b_cell); else - this->inverse(cell).vmult(x_cell, b_cell); + this->inverse_vmult(cell, x_cell, b_cell); } else { @@ -347,9 +369,9 @@ void PreconditionBlock::backward_step ( if (this->inverses_ready()) { if (transpose_diagonal) - this->inverse(cell).Tvmult(x_cell, b_cell); + this->inverse_Tvmult(cell, x_cell, b_cell); else - this->inverse(cell).vmult(x_cell, b_cell); + this->inverse_vmult(cell, x_cell, b_cell); } else { @@ -547,7 +569,7 @@ void PreconditionBlockJacobi { b_cell(row_cell)=src(row); } - this->inverse(cell).vmult(x_cell, b_cell); + this->inverse_vmult(cell, x_cell, b_cell); // distribute x_cell to dst for (row=cell*this->blocksize, row_cell=0; row_cellblocksize; @@ -737,9 +759,9 @@ void PreconditionBlockSOR::forward ( if (this->inverses_ready()) { if (transpose_diagonal) - this->inverse(cell).Tvmult(x_cell, b_cell); + this->inverse_Tvmult(cell, x_cell, b_cell); else - this->inverse(cell).vmult(x_cell, b_cell); + this->inverse_vmult(cell, x_cell, b_cell); } else { @@ -851,9 +873,9 @@ void PreconditionBlockSOR::backward ( if (this->inverses_ready()) { if (transpose_diagonal) - this->inverse(cell).Tvmult(x_cell, b_cell); + this->inverse_Tvmult(cell, x_cell, b_cell); else - this->inverse(cell).vmult(x_cell, b_cell); + this->inverse_vmult(cell, x_cell, b_cell); } else { diff --git a/deal.II/include/deal.II/lac/precondition_block_base.h b/deal.II/include/deal.II/lac/precondition_block_base.h index b7962bd7fb..d8bc501f1b 100644 --- a/deal.II/include/deal.II/lac/precondition_block_base.h +++ b/deal.II/include/deal.II/lac/precondition_block_base.h @@ -19,6 +19,7 @@ #include #include #include +#include #include @@ -50,11 +51,32 @@ template class PreconditionBlockBase { public: + /** + * Choose a method for inverting + * the blocks, and thus the data + * type for the inverses. + */ + enum Inversion + { + /** + * Use the standard + * Gauss-Jacobi method + * implemented in FullMatrix::inverse(). + */ + gauss_jordan, + /** + * Use QR decomposition of + * the Householder class. + */ + householder + }; + /** * Constructor initializing * default values. */ - PreconditionBlockBase(bool store_diagonals = false); + PreconditionBlockBase(bool store_diagonals = false, + Inversion method = gauss_jordan); /** * The virtual destructor @@ -78,7 +100,8 @@ class PreconditionBlockBase * then only one block will be * stored. */ - void reinit(unsigned int nblocks, unsigned int blocksize, bool compress); + void reinit(unsigned int nblocks, unsigned int blocksize, bool compress, + Inversion method = gauss_jordan); /** * Tell the class that inverses @@ -137,13 +160,33 @@ class PreconditionBlockBase * are stored. */ number el(unsigned int i, unsigned int j) const; + + /** + * Multiply with the inverse + * block at position i. + */ + template + void inverse_vmult(unsigned int i, Vector& dst, const Vector& src) const; + + /** + * Multiply with the transposed inverse + * block at position i. + */ + template + void inverse_Tvmult(unsigned int i, Vector& dst, const Vector& src) const; /** * Access to the inverse diagonal - * blocks. + * blocks if Inversion is #gauss_jordan. */ FullMatrix& inverse (unsigned int i); + /** + * Access to the inverse diagonal + * blocks if Inversion is #householder. + */ + Householder& inverse_householder (unsigned int i); + /** * Access to the inverse diagonal * blocks. @@ -177,6 +220,12 @@ class PreconditionBlockBase */ DeclException0 (ExcDiagonalsNotStored); + protected: + /** + * The method used for inverting blocks. + */ + Inversion inversion; + private: /** * The number of (inverse) @@ -185,17 +234,31 @@ class PreconditionBlockBase */ unsigned int n_diagonal_blocks; - /** + /** * Storage of the inverse * matrices of the diagonal * blocks matrices as * FullMatrix - * matrices. Using + * matrices, if Inversion + * #gauss_jordan is used. Using * number=float saves * memory in comparison with * number=double. */ - std::vector > var_inverse; + std::vector > var_inverse_full; + + /** + * Storage of the inverse + * matrices of the diagonal + * blocks matrices as + * Householder + * matrices if Inversion + * #householder is used. Using + * number=float saves + * memory in comparison with + * number=double. + */ + std::vector > var_inverse_householder; /** * Storage of the original diagonal blocks. @@ -229,8 +292,10 @@ class PreconditionBlockBase template inline -PreconditionBlockBase::PreconditionBlockBase(bool store) +PreconditionBlockBase::PreconditionBlockBase( + bool store, Inversion method) : + inversion(method), n_diagonal_blocks(0), var_store_diagonals(store), var_same_diagonal(false), @@ -243,8 +308,10 @@ inline void PreconditionBlockBase::clear() { - if (var_inverse.size()!=0) - var_inverse.erase(var_inverse.begin(), var_inverse.end()); + if (var_inverse_full.size()!=0) + var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end()); + if (var_inverse_full.size()!=0) + var_inverse_householder.erase(var_inverse_householder.begin(), var_inverse_householder.end()); if (var_diagonal.size()!=0) var_diagonal.erase(var_diagonal.begin(), var_diagonal.end()); var_same_diagonal = false; @@ -255,16 +322,29 @@ PreconditionBlockBase::clear() template inline void -PreconditionBlockBase::reinit(unsigned int n, unsigned int b, bool compress) +PreconditionBlockBase::reinit(unsigned int n, unsigned int b, bool compress, +Inversion method) { + inversion = method; var_same_diagonal = compress; var_inverses_ready = false; n_diagonal_blocks = n; if (compress) { - var_inverse.resize(1); - var_inverse[0].reinit(b,b); + switch(inversion) + { + case gauss_jordan: + var_inverse_full.resize(1); + var_inverse_full[0].reinit(b,b); + break; + case householder: + var_inverse_householder.resize(1); + break; + default: + Assert(false, ExcNotImplemented()); + } + if (store_diagonals()) { var_diagonal.resize(1); @@ -287,15 +367,21 @@ PreconditionBlockBase::reinit(unsigned int n, unsigned int b, bool compr tmp(n, FullMatrix(b)); var_diagonal.swap (tmp); } - - if (true) + + switch(inversion) { - std::vector > - tmp(n, FullMatrix(b)); - var_inverse.swap (tmp); - // make sure the tmp object - // goes out of scope as - // soon as possible + case gauss_jordan: + { + std::vector > + tmp(n, FullMatrix(b)); + var_inverse_full.swap (tmp); + break; + } + case householder: + var_inverse_householder.resize(n); + break; + default: + Assert(false, ExcNotImplemented()); } } } @@ -309,6 +395,55 @@ PreconditionBlockBase::size() const return n_diagonal_blocks; } +template +template +inline +void +PreconditionBlockBase::inverse_vmult( + unsigned int i, Vector& dst, const Vector& src) const +{ + const unsigned int ii = same_diagonal() ? 0U : i; + + switch (inversion) + { + case gauss_jordan: + AssertIndexRange (ii, var_inverse_full.size()); + var_inverse_full[ii].vmult(dst, src); + break; + case householder: + AssertIndexRange (ii, var_inverse_householder.size()); + var_inverse_householder[ii].vmult(dst, src); + break; + default: + Assert(false, ExcNotImplemented()); + } +} + + +template +template +inline +void +PreconditionBlockBase::inverse_Tvmult( + unsigned int i, Vector& dst, const Vector& src) const +{ + const unsigned int ii = same_diagonal() ? 0U : i; + + switch (inversion) + { + case gauss_jordan: + AssertIndexRange (ii, var_inverse_full.size()); + var_inverse_full[ii].Tvmult(dst, src); + break; + case householder: + AssertIndexRange (ii, var_inverse_householder.size()); + var_inverse_householder[ii].Tvmult(dst, src); + break; + default: + Assert(false, ExcNotImplemented()); + } +} + template inline @@ -316,10 +451,10 @@ const FullMatrix& PreconditionBlockBase::inverse(unsigned int i) const { if (same_diagonal()) - return var_inverse[0]; + return var_inverse_full[0]; - Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size())); - return var_inverse[i]; + Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size())); + return var_inverse_full[i]; } @@ -344,10 +479,23 @@ FullMatrix& PreconditionBlockBase::inverse(unsigned int i) { if (same_diagonal()) - return var_inverse[0]; + return var_inverse_full[0]; + + Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size())); + return var_inverse_full[i]; +} + + +template +inline +Householder& +PreconditionBlockBase::inverse_householder(unsigned int i) +{ + if (same_diagonal()) + return var_inverse_householder[0]; - Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size())); - return var_inverse[i]; + AssertIndexRange (i, var_inverse_householder.size()); + return var_inverse_householder[i]; } @@ -404,8 +552,8 @@ unsigned int PreconditionBlockBase::memory_consumption () const { unsigned int mem = sizeof(*this); - for (unsigned int i=0; i::Inversion inversion; }; /** diff --git a/deal.II/include/deal.II/lac/relaxation_block.templates.h b/deal.II/include/deal.II/lac/relaxation_block.templates.h index 4206ff170b..003ea60b6d 100644 --- a/deal.II/include/deal.II/lac/relaxation_block.templates.h +++ b/deal.II/include/deal.II/lac/relaxation_block.templates.h @@ -24,7 +24,8 @@ RelaxationBlock::AdditionalData::AdditionalData () : relaxation(1.), invert_diagonal(true), - same_diagonal(false) + same_diagonal(false), + inversion(PreconditionBlockBase::gauss_jordan) {} @@ -39,7 +40,8 @@ RelaxationBlock::AdditionalData::AdditionalData ( block_list(&bl), relaxation(relaxation), invert_diagonal(invert_diagonal), - same_diagonal(same_diagonal) + same_diagonal(same_diagonal), + inversion(PreconditionBlockBase::gauss_jordan) {} @@ -57,7 +59,8 @@ RelaxationBlock::initialize ( A = &M; additional_data = parameters; - this->reinit(additional_data.block_list->size(), 0, additional_data.same_diagonal); + this->reinit(additional_data.block_list->size(), 0, additional_data.same_diagonal, + additional_data.inversion); if (parameters.invert_diagonal) invert_diagblocks(); @@ -120,8 +123,18 @@ RelaxationBlock::invert_diagblocks () this->diagonal(block).reinit(bs, bs); this->diagonal(block) = M_cell; } - this->inverse(block).reinit(bs, bs); - this->inverse(block).invert(M_cell); + switch(this->inversion) + { + case PreconditionBlockBase::gauss_jordan: + this->inverse(block).reinit(bs, bs); + this->inverse(block).invert(M_cell); + break; + case PreconditionBlockBase::householder: + this->inverse_householder(block).initialize(M_cell); + break; + default: + Assert(false, ExcNotImplemented()); + } } } this->inverses_computed(true); @@ -162,7 +175,7 @@ RelaxationBlock::do_step ( b_cell(row_cell) -= entry->value() * prev(entry->column()); } // Apply inverse diagonal - this->inverse(block).vmult(x_cell, b_cell); + this->inverse_vmult(block, x_cell, b_cell); // Store in result vector row=additional_data.block_list->begin(block); for (unsigned int row_cell=0; row_cell