From: guido Date: Tue, 19 Jul 2005 11:17:36 +0000 (+0000) Subject: new class InverseMatrixRichardson X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=47c4c2dc91bd6eb1581420ff273c86b9132336ad;p=dealii-svn.git new class InverseMatrixRichardson git-svn-id: https://svn.dealii.org/trunk@11178 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/matrix_lib.h b/deal.II/lac/include/lac/matrix_lib.h index 36a674c85f..02d7e4bdca 100644 --- a/deal.II/lac/include/lac/matrix_lib.h +++ b/deal.II/lac/include/lac/matrix_lib.h @@ -16,6 +16,7 @@ #include #include #include +#include template class Vector; template class BlockVector; @@ -292,6 +293,97 @@ class MeanValueFilter : public Subscriptor unsigned int component; }; + + +/** + * Inverse matrix computed approximately by using the SolverRichardson + * iterative solver. In particular, the function + * SolverRichardson::Tsolve() allows for the implementation of + * transpose matrix vector products. + * + * The functions vmult() and Tvmult() appoximate the inverse + * iteratively starting with the vector dst. Functions + * vmult_add() and Tvmult_add() start the iteration with a zero + * vector. + * + * @ref Instantiations: some (Vector, Vector, BlockVector, BlockVector) + * @author Guido Kanschat, 2005 + */ +template +class InverseMatrixRichardson +{ + public: + /** + * Constructor, initializing the + * solver with a control and + * memory object. The inverted + * matrix and the preconditioner + * are added in initialize(). + */ + InverseMatrixRichardson (SolverControl& control, + VectorMemory& mem); + + /** + * Initialization + * function. Provide a solver + * object, a matrix, and another + * preconditioner for this. + */ + template + void initialize (const MATRIX&, + const PRECONDITION&); + + /** + * Access to the SolverControl + * object used by the solver. + */ + SolverControl& control() const; + /** + * Execute solver. + */ + void vmult (VECTOR&, const VECTOR&) const; + + /** + * Execute solver. + */ + void vmult_add (VECTOR&, const VECTOR&) const; + + /** + * Execute transpose solver. + */ + void Tvmult (VECTOR&, const VECTOR&) const; + + /** + * Execute transpose solver. + */ + void Tvmult_add (VECTOR&, const VECTOR&) const; + + private: + /** + * Access to the provided + * VectorMemory object. + */ + mutable VectorMemory& mem; + + /** + * The solver object. + */ + mutable SolverRichardson solver; + + /** + * The matrix in use. + */ + PointerMatrixBase* matrix; + + /** + * The preconditioner to use. + */ + PointerMatrixBase* precondition; +}; + + + + /*@}*/ //--------------------------------------------------------------------------- @@ -403,5 +495,16 @@ MeanValueFilter::Tvmult_add(VECTOR&, const VECTOR&) const Assert(false, ExcNotImplemented()); } +//-----------------------------------------------------------------------// + +template +template +inline void +InverseMatrixRichardson::initialize (const MATRIX& m, const PRECONDITION& p) +{ + matrix = &m; + precondition = &p; +} + #endif diff --git a/deal.II/lac/include/lac/matrix_lib.templates.h b/deal.II/lac/include/lac/matrix_lib.templates.h index 78752e46c2..e3498494fa 100644 --- a/deal.II/lac/include/lac/matrix_lib.templates.h +++ b/deal.II/lac/include/lac/matrix_lib.templates.h @@ -116,4 +116,75 @@ MeanValueFilter::vmult_add(BlockVector& dst, dst.block(i).add(src.block(i)); } + +//----------------------------------------------------------------------// + + +template +InverseMatrixRichardson::InverseMatrixRichardson( + SolverControl& c, + VectorMemory& m) + : + mem(m), + solver(c,m), + matrix(0), + precondition(0) +{} + + + + +template +void +InverseMatrixRichardson::vmult(VECTOR& dst, const VECTOR& src) const +{ + Assert (matrix != 0, ExcNotInitialized()); + Assert (precondition != 0, ExcNotInitialized()); + solver.solve(*matrix, dst, src, *precondition); +} + + + +template +void +InverseMatrixRichardson::vmult_add(VECTOR& dst, const VECTOR& src) const +{ + Assert (matrix != 0, ExcNotInitialized()); + Assert (precondition != 0, ExcNotInitialized()); + VECTOR* aux = mem.alloc(); + aux->reinit(dst); + solver.solve(*matrix, *aux, src, *precondition); + dst += *aux; + mem.free(aux); +} + + + +template +void +InverseMatrixRichardson::Tvmult(VECTOR& dst, const VECTOR& src) const +{ + Assert (matrix != 0, ExcNotInitialized()); + Assert (precondition != 0, ExcNotInitialized()); + solver.Tsolve(*matrix, dst, src, *precondition); +} + + + +template +void +InverseMatrixRichardson::Tvmult_add(VECTOR& dst, const VECTOR& src) const +{ + Assert (matrix != 0, ExcNotInitialized()); + Assert (precondition != 0, ExcNotInitialized()); + VECTOR* aux = mem.alloc(); + aux->reinit(dst); + solver.Tsolve(*matrix, *aux, src, *precondition); + dst += *aux; + mem.free(aux); +} + + + + #endif diff --git a/deal.II/lac/source/matrix_lib.cc b/deal.II/lac/source/matrix_lib.cc index 4031c13e7b..0f4ee74df0 100644 --- a/deal.II/lac/source/matrix_lib.cc +++ b/deal.II/lac/source/matrix_lib.cc @@ -118,3 +118,7 @@ template void MeanValueFilter::vmult_add(BlockVector&, template void MeanValueFilter::vmult_add(BlockVector&, const BlockVector&) const; +template class InverseMatrixRichardson >; +template class InverseMatrixRichardson >; +template class InverseMatrixRichardson >; +template class InverseMatrixRichardson >;