From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
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 <base/subscriptor.h>
 #include <lac/vector_memory.h>
 #include <lac/pointer_matrix.h>
+#include <lac/solver_richardson.h>
 
 template<typename number> class Vector;
 template<typename number> 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 <tt>dst</tt>. Functions
+ * vmult_add() and Tvmult_add() start the iteration with a zero
+ * vector.
+ *
+ * @ref Instantiations: some (Vector<float>, Vector<double>, BlockVector<float>, BlockVector<double>)
+ * @author Guido Kanschat, 2005
+ */
+template<class VECTOR>
+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<VECTOR>& mem);
+    
+				     /**
+				      * Initialization
+				      * function. Provide a solver
+				      * object, a matrix, and another
+				      * preconditioner for this.
+				      */
+    template <class MATRIX, class PRECONDITION>
+    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<VECTOR>& mem;
+    
+				     /**
+				      * The solver object.
+				      */
+    mutable SolverRichardson<VECTOR> solver;
+    
+				     /**
+				      * The matrix in use.
+				      */
+    PointerMatrixBase<VECTOR>* matrix;
+    
+				     /**
+				      * The preconditioner to use.
+				      */
+    PointerMatrixBase<VECTOR>* precondition;
+};
+
+
+
+
 /*@}*/
 //---------------------------------------------------------------------------
 
@@ -403,5 +495,16 @@ MeanValueFilter::Tvmult_add(VECTOR&, const VECTOR&) const
   Assert(false, ExcNotImplemented());
 }
 
+//-----------------------------------------------------------------------//
+
+template <class VECTOR>
+template <class MATRIX, class PRECONDITION>
+inline void
+InverseMatrixRichardson<VECTOR>::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<number>& dst,
       dst.block(i).add(src.block(i));
 }
 
+
+//----------------------------------------------------------------------//
+
+
+template <class VECTOR>
+InverseMatrixRichardson<VECTOR>::InverseMatrixRichardson(
+  SolverControl& c,
+  VectorMemory<VECTOR>& m)
+		:
+		mem(m),
+		solver(c,m),
+		matrix(0),
+		precondition(0)
+{}
+
+
+
+
+template <class VECTOR>
+void
+InverseMatrixRichardson<VECTOR>::vmult(VECTOR& dst, const VECTOR& src) const
+{
+  Assert (matrix != 0, ExcNotInitialized());
+  Assert (precondition != 0, ExcNotInitialized());
+  solver.solve(*matrix, dst, src, *precondition);
+}
+
+
+
+template <class VECTOR>
+void
+InverseMatrixRichardson<VECTOR>::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 <class VECTOR>
+void
+InverseMatrixRichardson<VECTOR>::Tvmult(VECTOR& dst, const VECTOR& src) const
+{
+  Assert (matrix != 0, ExcNotInitialized());
+  Assert (precondition != 0, ExcNotInitialized());
+  solver.Tsolve(*matrix, dst, src, *precondition);
+}
+
+
+
+template <class VECTOR>
+void
+InverseMatrixRichardson<VECTOR>::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<float>&,
 template void MeanValueFilter::vmult_add(BlockVector<double>&,
 					 const BlockVector<double>&) const;
 
+template class InverseMatrixRichardson<Vector<float> >;
+template class InverseMatrixRichardson<Vector<double> >;
+template class InverseMatrixRichardson<BlockVector<float> >;
+template class InverseMatrixRichardson<BlockVector<double> >;