From 348e79a5b4af19f485cb55b65e030a2ba98ea73f Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 9 Nov 2010 15:01:03 +0000 Subject: [PATCH] add option of inverting blocks with SVD git-svn-id: https://svn.dealii.org/trunk@22643 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/lac/matrix_block.h | 1 + .../deal.II/lac/precondition_block_base.h | 38 ++++++++++++++++ .../include/deal.II/lac/relaxation_block.h | 44 +++++++++++++++++++ .../deal.II/lac/relaxation_block.templates.h | 8 ++-- 4 files changed, 88 insertions(+), 3 deletions(-) diff --git a/deal.II/include/deal.II/lac/matrix_block.h b/deal.II/include/deal.II/lac/matrix_block.h index 1dba79fd2a..ca8d59c04a 100644 --- a/deal.II/include/deal.II/lac/matrix_block.h +++ b/deal.II/include/deal.II/lac/matrix_block.h @@ -465,6 +465,7 @@ class MatrixBlockVector NamedData >::subscribe; NamedData >::unsubscribe; NamedData >::size; + NamedData >::name; }; 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 0b15a13d43..b3c3a6a004 100644 --- a/deal.II/include/deal.II/lac/precondition_block_base.h +++ b/deal.II/include/deal.II/lac/precondition_block_base.h @@ -205,6 +205,18 @@ class PreconditionBlockBase */ const FullMatrix& inverse (unsigned int i) const; + /** + * Access to the inverse diagonal + * blocks if Inversion is #householder. + */ + const Householder& inverse_householder (unsigned int i) const; + + /** + * Access to the inverse diagonal + * blocks if Inversion is #householder. + */ + const LAPACKFullMatrix& inverse_svd (unsigned int i) const; + /** * Access to the diagonal * blocks. @@ -507,6 +519,32 @@ PreconditionBlockBase::inverse(unsigned int i) const } +template +inline +const Householder& +PreconditionBlockBase::inverse_householder(unsigned int i) const +{ + if (same_diagonal()) + return var_inverse_householder[0]; + + AssertIndexRange (i, var_inverse_householder.size()); + return var_inverse_householder[i]; +} + + +template +inline +const LAPACKFullMatrix& +PreconditionBlockBase::inverse_svd(unsigned int i) const +{ + if (same_diagonal()) + return var_inverse_svd[0]; + + AssertIndexRange (i, var_inverse_svd.size()); + return var_inverse_svd[i]; +} + + template inline const FullMatrix& diff --git a/deal.II/include/deal.II/lac/relaxation_block.h b/deal.II/include/deal.II/lac/relaxation_block.h index 8052b27c13..d561195593 100644 --- a/deal.II/include/deal.II/lac/relaxation_block.h +++ b/deal.II/include/deal.II/lac/relaxation_block.h @@ -106,6 +106,17 @@ class RelaxationBlock : * method for the blocks. */ typename PreconditionBlockBase::Inversion inversion; + + /** + * The if #inversion is SVD, + * the threshold below which + * a singular value will be + * considered zero and thus + * not inverted. This + * parameter is used in the + * call to LAPACKFullMatrix::compute_inverse_svd(). + */ + double threshold; }; /** @@ -267,6 +278,22 @@ class RelaxationBlockSOR : public virtual Subscriptor, * Make function of base class public again. */ using RelaxationBlock::empty; + /** + * Make function of base class public again. + */ + using RelaxationBlock::size; + /** + * Make function of base class public again. + */ + using RelaxationBlock::inverse; + /** + * Make function of base class public again. + */ + using RelaxationBlock::inverse_householder; + /** + * Make function of base class public again. + */ + using RelaxationBlock::inverse_svd; /** * Perform one step of the SOR * iteration. @@ -335,6 +362,23 @@ class RelaxationBlockSSOR : public virtual Subscriptor, * Make function of base class public again. */ using RelaxationBlock::empty; + + /** + * Make function of base class public again. + */ + using RelaxationBlock::size; + /** + * Make function of base class public again. + */ + using RelaxationBlock::inverse; + /** + * Make function of base class public again. + */ + using RelaxationBlock::inverse_householder; + /** + * Make function of base class public again. + */ + using RelaxationBlock::inverse_svd; /** * Perform one step of the SOR * iteration. 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 6298b4bfb3..dbc3a56980 100644 --- a/deal.II/include/deal.II/lac/relaxation_block.templates.h +++ b/deal.II/include/deal.II/lac/relaxation_block.templates.h @@ -25,7 +25,8 @@ RelaxationBlock::AdditionalData::AdditionalData () relaxation(1.), invert_diagonal(true), same_diagonal(false), - inversion(PreconditionBlockBase::gauss_jordan) + inversion(PreconditionBlockBase::gauss_jordan), + threshold(0.) {} @@ -41,7 +42,8 @@ RelaxationBlock::AdditionalData::AdditionalData ( relaxation(relaxation), invert_diagonal(invert_diagonal), same_diagonal(same_diagonal), - inversion(PreconditionBlockBase::gauss_jordan) + inversion(PreconditionBlockBase::gauss_jordan), + threshold(0.) {} @@ -135,7 +137,7 @@ RelaxationBlock::invert_diagblocks () case PreconditionBlockBase::svd: this->inverse_svd(block).reinit(bs, bs); this->inverse_svd(block) = M_cell; - this->inverse_svd(block).compute_inverse_svd(0.); + this->inverse_svd(block).compute_inverse_svd(additional_data.threshold); break; default: Assert(false, ExcNotImplemented()); -- 2.39.5