]> https://gitweb.dealii.org/ - dealii.git/commitdiff
use library class instead of hand woven thingy
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 2 Jul 2010 20:30:37 +0000 (20:30 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 2 Jul 2010 20:30:37 +0000 (20:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@21452 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-20/doc/intro.dox
deal.II/examples/step-20/doc/results.dox
deal.II/examples/step-20/step-20.cc

index 0f992b73621218bcae5d838ae918ae0e57ec6bc0..91c7b3ca77ae3127560d9f6ee4b340f16d30350d 100644 (file)
@@ -397,46 +397,7 @@ matrix-vector products with it. We can do so in steps: to compute $Sv$, we
   positive definite and symmetric mass matrix $M$;
  <li> form $z=By$ to obtain $Sv=z$.
 </ol>
-We will implement a class that does that in the program. Before showing its
-code, let us first note that we need to multiply with $M^{-1}$ in several
-places here: in multiplying with the Schur complement $S$, forming the right
-hand side of the first equation, and solving in the second equation. From a
-coding viewpoint, it is therefore appropriate to relegate such a recurring
-operation to a class of its own. We call it <code>InverseMatrix</code>. As far as
-linear solvers are concerned, this class will have all operations that solvers
-need, which in fact includes only the ability to perform matrix-vector
-products; we form them by using a CG solve (this of course requires that the
-matrix passed to this class satisfies the requirements of the CG
-solvers). Here are the relevant parts of the code that implements this:
-
-@code
-class InverseMatrix
-{
-  public:
-    InverseMatrix (const SparseMatrix<double> &m);
-
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const SparseMatrix<double> > matrix;
-    // ...
-};
-
-
-void InverseMatrix::vmult (Vector<double>       &dst,
-                           const Vector<double> &src) const
-{
-  SolverControl solver_control (src.size(), 1e-8*src.l2_norm());
-  SolverCG<>    cg (solver_control, vector_memory);
-
-  cg.solve (*matrix, dst, src, PreconditionIdentity());        
-}
-@endcode
-
-Once created, objects of this class can act as matrices: they perform
-matrix-vector multiplications. How this is actually done is irrelevant to the
-outside world.
+This is accomplished by using the class IterativeInverse.
 
 Using this class, we can then write a class that implements the Schur
 complement in much the same way: to act as a matrix, it only needs to offer a
index bf1d2e814e7b46da0b75be723b4d2f33fb4feb72..4ea4ae7634c8af3f897db769443986a1dfd4e214 100644 (file)
@@ -13,7 +13,7 @@ examples/step-20> make run
 Number of active cells: 64
 Total number of cells: 85
 Number of degrees of freedom: 208 (144+64)
-10 CG Schur complement iterations to obtain convergence.
+15 CG Schur complement iterations to obtain convergence.
 Errors: ||e_p||_L2 = 0.178055,   ||e_u||_L2 = 0.0433435
 @endcode
 
index 36f38faf298224b0f8a3afdd0586a29125eec7d2..febeb0602f06d83387a3efc2e1cf8bf9f0df948b 100644 (file)
@@ -3,7 +3,7 @@
 
 /*    $Id$       */
 /*                                                                */
-/*    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     */
 #include <lac/block_sparse_matrix.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
+                                // For our Schur complement solver,
+                                // we need two new objects. One is a
+                                // matrix object which acts as the
+                                // inverse of a matrix by calling an
+                                // iterative solver.
+#include <lac/iterative_inverse.h>
+
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <grid/tria_accessor.h>
@@ -719,121 +726,6 @@ void MixedLaplaceProblem<dim>::assemble_system ()
                                  // rather only comment on
                                  // implementational aspects.
 
-                                 // @sect4{The <code>InverseMatrix</code> class template}
-
-                                 // The first component of our linear
-                                 // solver scheme was the creation of
-                                 // a class that acts like the inverse
-                                 // of a matrix, i.e. which has a
-                                 // <code>vmult</code> function that multiplies
-                                 // a vector with an inverse matrix by
-                                 // solving a linear system.
-                                 //
-                                 // While most of the code below
-                                 // should be obvious given the
-                                 // purpose of this class, two
-                                 // comments are in order. First, the
-                                 // class is derived from the
-                                 // <code>Subscriptor</code> class so that we
-                                 // can use the <code>SmartPointer</code> class
-                                 // with inverse matrix objects. The
-                                 // use of the <code>Subscriptor</code> class
-                                 // has been explained before in
-                                 // step-7 and step-20. The present
-                                 // class also sits on the receiving
-                                 // end of this
-                                 // <code>Subscriptor</code>/<code>SmartPointer</code>
-                                 // pair: it holds its pointer to the
-                                 // matrix it is supposed to be the
-                                 // inverse of through a
-                                 // <code>SmartPointer</code> to make sure that
-                                 // this matrix is not destroyed while
-                                 // we still have a pointer to it.
-                                 //
-                                 // Secondly, we realize that we will
-                                 // probably perform many
-                                 // matrix-vector products with
-                                 // inverse matrix objects. Now, every
-                                 // time we do so, we have to call the
-                                 // CG solver to solve a linear
-                                 // system. To work, the CG solver
-                                 // needs to allocate four temporary
-                                 // vectors that it will release again
-                                 // at the end of its operation. What
-                                 // this means is that through
-                                 // repeated calls to the <code>vmult</code>
-                                 // function of this class we have to
-                                 // allocate and release vectors over
-                                 // and over again.
-                                 //
-                                 // The natural question is then:
-                                 // Wouldn't it be nice if we could
-                                 // avoid this, and allocate vectors
-                                 // only once? In fact, deal.II offers
-                                 // a way to do exactly this and we
-                                 // don't even have to do anything
-                                 // special about it (so this comment
-                                 // is purely educational). What all
-                                 // the linear solvers do is not to
-                                 // allocate memory using
-                                 // <code>new</code> and
-                                 // <code>delete</code>, but rather to
-                                 // allocate them from an object
-                                 // derived from the
-                                 // <code>VectorMemory</code> class
-                                 // (see the module on Vector memory
-                                 // management in the API reference
-                                 // manual). By default, the linear
-                                 // solvers use a derived class
-                                 // <code>GrowingVectorMemory</code>
-                                 // that, every time a vector is
-                                 // requested, allocates one from a
-                                 // pool that is shared by all
-                                 // <code>GrowingVectorMemory</code>
-                                 // objects.
-template <class Matrix>
-class InverseMatrix : public Subscriptor
-{
-  public:
-    InverseMatrix (const Matrix &m);
-
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const Matrix> matrix;
-};
-
-
-template <class Matrix>
-InverseMatrix<Matrix>::InverseMatrix (const Matrix &m)
-                :
-                matrix (&m)
-{}
-
-
-                                 // Here now is the function that
-                                 // implements multiplication with the
-                                 // inverse matrix by calling a CG
-                                 // solver. Note that we set the solution
-                                 // vector to zero before starting the
-                                 // solve, since we do not want to use
-                                 // the possible previous and unknown
-                                 // content of that variable as
-                                 // starting vector for the linear
-                                 // solve:
-template <class Matrix>
-void InverseMatrix<Matrix>::vmult (Vector<double>       &dst,
-                                   const Vector<double> &src) const
-{
-  SolverControl solver_control (src.size(), 1e-8*src.l2_norm());
-  SolverCG<>    cg (solver_control);
-
-  dst = 0;
-  
-  cg.solve (*matrix, dst, src, PreconditionIdentity());        
-}
-
 
                                  // @sect4{The <code>SchurComplement</code> class template}
 
@@ -882,21 +774,21 @@ class SchurComplement : public Subscriptor
 {
   public:
     SchurComplement (const BlockSparseMatrix<double> &A,
-                     const InverseMatrix<SparseMatrix<double> > &Minv);
+                     const IterativeInverse<Vector<double> > &Minv);
 
     void vmult (Vector<double>       &dst,
                 const Vector<double> &src) const;
 
   private:
     const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-    const SmartPointer<const InverseMatrix<SparseMatrix<double> > > m_inverse;
+    const SmartPointer<const IterativeInverse<Vector<double> > > m_inverse;
     
     mutable Vector<double> tmp1, tmp2;
 };
 
 
 SchurComplement::SchurComplement (const BlockSparseMatrix<double> &A,
-                                  const InverseMatrix<SparseMatrix<double> > &Minv)
+                                  const IterativeInverse<Vector<double> > &Minv)
                 :
                 system_matrix (&A),
                 m_inverse (&Minv),
@@ -920,7 +812,7 @@ void SchurComplement::vmult (Vector<double>       &dst,
                                  // and preconditioner system is the
                                  // class that approximates the Schur
                                  // complement so we can form a
-                                 // <code>InverseMatrix@<ApproximateSchurComplement@></code>
+                                 // an InverseIterate
                                  // object that approximates the
                                  // inverse of the Schur
                                  // complement. It follows the same
@@ -932,6 +824,12 @@ void SchurComplement::vmult (Vector<double>       &dst,
                                  // step. Consequently, the class also
                                  // does not have to store a pointer
                                  // to an inverse mass matrix object.
+                                //
+                                // Since InverseIterate follows the
+                                // standard convention for matrices,
+                                // we need to provide a
+                                // <tt>Tvmult</tt> function here as
+                                // well.
 class ApproximateSchurComplement : public Subscriptor
 {
   public:
@@ -939,6 +837,8 @@ class ApproximateSchurComplement : public Subscriptor
 
     void vmult (Vector<double>       &dst,
                 const Vector<double> &src) const;
+    void Tvmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
 
   private:
     const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
@@ -964,6 +864,15 @@ void ApproximateSchurComplement::vmult (Vector<double>       &dst,
 }
 
 
+void ApproximateSchurComplement::Tvmult (Vector<double>       &dst,
+                                        const Vector<double> &src) const
+{
+  system_matrix->block(1,0).Tvmult (dst, tmp2);
+  system_matrix->block(0,0).precondition_Jacobi (tmp2, tmp1);
+  system_matrix->block(0,1).Tvmult (tmp1, src);
+}
+
+
 
                                  // @sect4{MixedLaplace::solve}
 
@@ -985,8 +894,14 @@ void ApproximateSchurComplement::vmult (Vector<double>       &dst,
 template <int dim>
 void MixedLaplaceProblem<dim>::solve () 
 {
-  const InverseMatrix<SparseMatrix<double> >
-    m_inverse (system_matrix.block(0,0));
+  PreconditionIdentity identity;
+  IterativeInverse<Vector<double> >
+    m_inverse;
+  m_inverse.initialize(system_matrix.block(0,0), identity);
+  m_inverse.solver.select("cg");
+  ReductionControl inner_control(1000, 0., 1.e-13);
+  m_inverse.solver.control = inner_control;
+  
   Vector<double> tmp (solution.block(0).size());
 
                                    // Now on to the first
@@ -1016,8 +931,11 @@ void MixedLaplaceProblem<dim>::solve ()
     ApproximateSchurComplement
       approximate_schur_complement (system_matrix);
       
-    InverseMatrix<ApproximateSchurComplement>
-      preconditioner (approximate_schur_complement);
+    IterativeInverse<Vector<double> >
+      preconditioner;
+    preconditioner.initialize(approximate_schur_complement, identity);
+    preconditioner.solver.select("cg");
+    preconditioner.solver.control = inner_control;
 
     
     SolverControl solver_control (solution.block(1).size(),

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.