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 36a674c85fa2cd7acd1fea11e9f3f8076a6429e6..02d7e4bdcae9edeefcd03a7c8fa4f3b9c8efaa54 100644 (file)
--- 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 78752e46c2a1297886b04d95e6369a56a58daf53..e3498494fa0fa1bb2c256786155daf61f3190427 100644 (file)
--- 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 4031c13e7bf8b83de7d1177d4c4704643c6efce8..0f4ee74df04ceded9aac4891c43dc8fb7a2539ec 100644 (file)
--- 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> >;