]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Check the requirements on vectors and matrices in the linear solvers. Minor clean...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 14 Aug 2001 09:39:53 +0000 (09:39 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 14 Aug 2001 09:39:53 +0000 (09:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@4880 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/solver.h
deal.II/lac/include/lac/solver_bicgstab.h
deal.II/lac/include/lac/solver_cg.h
deal.II/lac/include/lac/solver_gmres.h
deal.II/lac/include/lac/solver_minres.h
deal.II/lac/include/lac/solver_qmrs.h
deal.II/lac/include/lac/solver_richardson.h

index f638719d1d3930b0b810f9391448db7db7330942..1a8b33a7ee9b5f1d46b099ab793df951185722b1 100644 (file)
@@ -19,84 +19,86 @@ class SolverControl;
 
 
 /**
- * Base class for iterative solvers.
- *
-//TODO:[?] * HAS TO BE UPDATED!
- *
- * This class defines possible
+ * Base class for iterative solvers.  This class defines possible
  * return states of linear solvers and provides interfaces to a memory
  * pool and the control object.
  *
  *
  * @sect3{Requirements for template classes}
  *
- * The class is templated to allow for different matrix and vector
- * classes, since iterative solvers do not rely on any special structure
- * of matrices or the format of storage. However, there are some common
- * requirements a matrix or vector type must fulfil to qualify as an
- * applicable type for the solvers in this hierarchy. These requirements
- * are listed following. The classes do not declare any concrete
- * class, they are rather intended to form a `signature' which a concrete
- * class has to conform to.
+ * Since iterative solvers do not rely on any special structure of
+ * matrices or the format of storage, but only require that matrices
+ * and vector define certain operations such as matrix-vector
+ * products, or scalar products between vectors, this class as well as
+ * the derived classes implementing concrete linear solvers are
+ * templated on the types of matrices and vectors. However, there are
+ * some common requirements a matrix or vector type must fulfill to
+ * qualify as an applicable type for the solvers in this
+ * hierarchy. These requirements are listed following. The listed
+ * classes are not any concrete class, they are rather intended to
+ * form a `signature' which a concrete class has to conform to. Note
+ * that the matrix and vector classes within this library of course
+ * conform to this interface.
  *
  * @begin{verbatim}
  * class Matrix
  * {
  *   public:
- *                        // Application to a Vector
- *     void vmult (Vector& dst, const Vector& src) const;
- *
- *                        // Application of a preconditioner to
- *                        // a Vector, i.e. $dst=\tilde A^(-1) src$,
- *                        // where $\tilde A^(-1)$ is an approximation
- *                        // to the inverse if the matrix stored in
- *                        // this object.
- *     void precondition (Vector& dst, const Vector& src) const;
+ *                        // Application of matrix to vector src.
+ *                        // write result into dst
+ *     void vmult (Vector &dst, const Vector &src) const;
  * 
  *                        // Application of transpose to a Vector.
- *                        // Only used by special iterative methods.
- *     void T_vmult (Vector& dst, const Vector& src) const;
- *
- *                        // Application of a transposed preconditioner
- *                        // to a Vector. Only used by special
- *                        // iterative methods
- *    
- *     void T_precondition (Vector& dst, const Vector& src) const;
+ *                        // Only used by certain iterative methods.
+ *     void Tvmult (Vector &dst, const Vector &src) const;
  * };
  *
  *
  * class Vector
  * {
  *   public:
+ *                        // resize and/or clear vector. note
+ *                        // that the second argument must have
+ *                        // a default value equal to false
+ *     void reinit (const unsigned int size,
+ *                  bool  leave_elements_uninitialized = false);
+ *
  *                        // scalar product
- *     double operator * (const Vectorv) const;
+ *     double operator * (const Vector &v) const;
  *
  *                        // addition of vectors
  *                        // $y = y + x$.
- *     void add (const Vector& x);
+ *     void add (const Vector &x);
+ *
  *                        // $y = y + ax$.
- *     void add (double a, const Vector& x);
+ *     void add (const double  a,
+ *               const Vector &x);
  *
- *                        // scaled addition of vectors
- *                        // $y = ay + x$.
- *     void sadd (double a,
- *                const Vector& x);
  *                        // $y = ay + bx$.
- *     void sadd (double a,
- *                double b, const Vector& x);
- *                        // $y = ay + bx + cz$.
- *     void sadd (double a,
- *                double b, const Vector& x,
- *                double c, const Vector& z);
+ *     void sadd (const double  a,
+ *                const double  b,
+ *                const Vector &x);
  * 
  *                        // $y = ax$.
- *     void equ (double a, const Vector& x);
- *                        // $y = ax + bz$.
- *     void equ (double a, const Vector& x,
- *               double b, const Vector& z);
+ *     void equ (const double  a,
+ *               const Vector &x);
+ *
+ *                        // scale the elements of the vector
+ *                        // by a fixed value
+ *     void scale (const double a);
+ *
+ *                        // return the l2 norm of the vector
+ *     double l2_norm () const;
  * };
  * @end{verbatim}
  *
+ * In addition, for some solvers there has to be a global function
+ * @p{swap(vector &a, vector &b)} that exchanges the values of the two vectors.
+ *
+ * The preconditioners used must have the same interface as matrices,
+ * i.e. in particular they have to provide a member function @p{vmult}
+ * which denotes the application of the preconditioner.
+ *
  *
  * @sect3{AdditionalData}
  *
index 9bff8e3e00b9bf82a3bd50be2212fb289d9d04af..534720c70867b504f9c1e6116e30b561417f9f34 100644 (file)
@@ -22,6 +22,9 @@
 /**
  * Bicgstab algorithm by van der Vorst.
  *
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
  * Like all other solver classes, this class has a local structure called
  * @p{AdditionalData} which is used to pass additional parameters to the
  * solver, like damping parameters or the number of temporary vectors. We
index efacba01da1c9699ab363c36b318d97cb5d64f6f..98c50a0f05ebb7d0b6362f65de157af783f14869 100644 (file)
@@ -25,6 +25,9 @@
 /**
  * Preconditioned cg method.
  *
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
  * Like all other solver classes, this class has a local structure called
  * @p{AdditionalData} which is used to pass additional parameters to the
  * solver, like damping parameters or the number of temporary vectors. We
@@ -85,14 +88,15 @@ class SolverCG : public Subscriptor, private Solver<VECTOR>
     virtual ~SolverCG ();
 
                                     /**
-                                     * Solver method.
+                                     * Solve the linear system $Ax=b$
+                                     * for x.
                                      */
     template<class MATRIX, class PRECONDITIONER>
     void
-    solve (const MATRIX &A,
-          VECTOR       &x,
-          const VECTOR &b,
-          const PRECONDITIONERprecondition);
+    solve (const MATRIX         &A,
+          VECTOR               &x,
+          const VECTOR         &b,
+          const PRECONDITIONER &precondition);
 
   protected:
                                     /**
@@ -189,10 +193,10 @@ SolverCG<VECTOR>::print_vectors(const unsigned int,
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
 void
-SolverCG<VECTOR>::solve (const MATRIX &A,
-                        VECTOR       &x,
-                        const VECTOR &b,
-                        const PRECONDITIONERprecondition)
+SolverCG<VECTOR>::solve (const MATRIX         &A,
+                        VECTOR               &x,
+                        const VECTOR         &b,
+                        const PRECONDITIONER &precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 
index b5a4103754d84e67131eb13fcf95ca0df924b0a2..ef6f3f6de2144715ebe5345e10b5631b9019d575 100644 (file)
@@ -69,6 +69,9 @@
  * of temporary vectors as commented upon above. By default, the number
  * of these vectors is set to 30.
  *
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
  * @author Wolfgang Bangerth
  */
 template <class VECTOR = Vector<double> >
@@ -107,13 +110,15 @@ class SolverGMRES : public Subscriptor, private Solver<VECTOR>
                 const AdditionalData &data=AdditionalData());
 
                                     /**
-                                     * Solver method.
+                                     * Solve the linear system $Ax=b$
+                                     * for x.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    void solve (const MATRIX &A,
-               VECTOR       &x,
-               const VECTOR &b,
-               const PRECONDITIONER& precondition);
+    void
+    solve (const MATRIX         &A,
+          VECTOR               &x,
+          const VECTOR         &b,
+          const PRECONDITIONER &precondition);
 
     DeclException1 (ExcTooFewTmpVectors,
                    int,
@@ -205,10 +210,10 @@ SolverGMRES<VECTOR>::givens_rotation (Vector<double> &h,
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
 void
-SolverGMRES<VECTOR>::solve (const MATRIXA,
-                                  VECTOR      & x,
-                                  const VECTOR& b,
-                                  const PRECONDITIONER& precondition)
+SolverGMRES<VECTOR>::solve (const MATRIX         &A,
+                           VECTOR               &x,
+                           const VECTOR         &b,
+                           const PRECONDITIONER &precondition)
 {
                                   // this code was written a very
                                   // long time ago by people not
index 1103da70d6b14d636aca054fc96129d1151e594b..711c7ba1d480e25aa978f103ab1349474aafa4a6 100644 (file)
@@ -24,6 +24,9 @@
 /**
  * Preconditioned MinRes method.
  *
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
  * Like all other solver classes, this class has a local structure called
  * @p{AdditionalData} which is used to pass additional parameters to the
  * solver, like damping parameters or the number of temporary vectors. We
@@ -74,14 +77,15 @@ class SolverMinRes : public Subscriptor, private Solver<VECTOR>
     virtual ~SolverMinRes ();
     
                                     /**
-                                     * Solver method.
+                                     * Solve the linear system $Ax=b$
+                                     * for x.
                                      */
     template<class MATRIX, class PRECONDITIONER>
     void
-    solve (const MATRIX &A,
-          VECTOR       &x,
-          const VECTOR &b,
-          const PRECONDITIONERprecondition);
+    solve (const MATRIX         &A,
+          VECTOR               &x,
+          const VECTOR         &b,
+          const PRECONDITIONER &precondition);
 
                                     /**
                                      * Exception
@@ -172,10 +176,10 @@ SolverMinRes<VECTOR>::print_vectors(const unsigned int,
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
 typename Solver<VECTOR>::ReturnState 
-SolverMinRes<VECTOR>::solve (const MATRIX &A,
-                            VECTOR       &x,
-                            const VECTOR &b,
-                            const PRECONDITIONERprecondition)
+SolverMinRes<VECTOR>::solve (const MATRIX         &A,
+                            VECTOR               &x,
+                            const VECTOR         &b,
+                            const PRECONDITIONER &precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 
@@ -246,21 +250,20 @@ SolverMinRes<VECTOR>::solve (const MATRIX &A,
   r_l2 = r0;
   
   
-  u[0].reinit(VS,0);
+  u[0].reinit(VS);
   delta[0] = 1.;
-  m[0].reinit(VS,0);
-  m[1].reinit(VS,0);
-  m[2].reinit(VS,0);
+  m[0].reinit(VS);
+  m[1].reinit(VS);
+  m[2].reinit(VS);
                                   
   conv = control().check(0,r_l2);
   
   while (conv==SolverControl::iterate)
-    {
-      
+    {      
       if (delta[1]!=0)
        v.scale(1./sqrt(delta[1]));
       else
-       v.reinit(VS,0);
+       v.reinit(VS);
 
       A.vmult(u[2],v);
       u[2].add (-sqrt(delta[1]/delta[0]), u[0]);
index b4c8788594c98341f437628a91fc561fb763ee6c..2bca2aec6c0b797e8cb18ed469e56e4c4e9c2aec 100644 (file)
@@ -32,6 +32,9 @@
  * preconditioner is used: left preconditioning seems to require the
  * inverse.
  *
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
  * Like all other solver classes, this class has a local structure called
  * @p{AdditionalData} which is used to pass additional parameters to the
  * solver, like damping parameters or the number of temporary vectors. We
@@ -106,14 +109,15 @@ class SolverQMRS : public Subscriptor, private Solver<VECTOR>
              const AdditionalData &data=AdditionalData());
 
                                     /**
-                                     * Solver method.
+                                     * Solve the linear system $Ax=b$
+                                     * for x.
                                      */
     template<class MATRIX, class PRECONDITIONER>
     void
-    solve (const MATRIX &A,
-          VECTOR       &x,
-          const VECTOR &b,
-          const PRECONDITIONERprecondition);
+    solve (const MATRIX         &A,
+          VECTOR               &x,
+          const VECTOR         &b,
+          const PRECONDITIONER &precondition);
 
                                     /**
                                      * Interface for derived class.
@@ -218,10 +222,10 @@ SolverQMRS<VECTOR>::print_vectors(const unsigned int,
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
 void
-SolverQMRS<VECTOR>::solve (const MATRIX &A,
-                          VECTOR       &x,
-                          const VECTOR &b,
-                          const PRECONDITIONERprecondition)
+SolverQMRS<VECTOR>::solve (const MATRIX         &A,
+                          VECTOR               &x,
+                          const VECTOR         &b,
+                          const PRECONDITIONER &precondition)
 {
   deallog.push("QMRS");
   
@@ -277,8 +281,8 @@ SolverQMRS<VECTOR>::solve (const MATRIX &A,
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
 bool
-SolverQMRS<VECTOR>::iterate(const MATRIXA,
-                           const PRECONDITIONERprecondition)
+SolverQMRS<VECTOR>::iterate(const MATRIX         &A,
+                           const PRECONDITIONER &precondition)
 {
 /* Remark: the matrix A in the article is the preconditioned matrix.
  * Therefore, we have to precondition x before we compute the first residual.
@@ -319,12 +323,9 @@ SolverQMRS<VECTOR>::iterate(const MATRIX& A,
   precondition.vmult(q,p);
  
   tau = v.norm_sqr();
-  //deallog << "tau:" << tau << std::endl;
   rho = q*v;
-  //deallog << "rho:" << rho << std::endl;
-
-
-while (state == SolverControl::iterate)
+  
+  while (state == SolverControl::iterate)
     {
       step++; it++;
                                       // Step 1
@@ -337,7 +338,6 @@ while (state == SolverControl::iterate)
        return true;
                                       // Step 3
       alpha = rho/sigma;
-      //deallog << "alpha:" << alpha << std::endl;
 
       v.add(-alpha,t);
                                       // Step 4
@@ -346,10 +346,6 @@ while (state == SolverControl::iterate)
       psi = 1./(1.+theta);
       tau *= theta*psi;
 
-      //deallog << "psi:" << psi << std::endl;
-      //deallog << "theta:" << theta << std::endl;
-      //deallog << "tau:" << tau << std::endl;
-
       d.sadd(psi*theta_old, psi*alpha, p);
       x.add(d);
 
index 71b87816d6dc9aef7fa2e58ac1fe8b01df288606..60682892c0e9a8e68d7a1bd11def134a2fb7bed5 100644 (file)
@@ -22,6 +22,9 @@
  * Implementation of the richardson iteration method. The stopping criterion
  * is the norm of the residual.
  *
+ * For the requirements on matrices and vectors in order to work with
+ * this class, see the documentation of the @ref{Solver} base class.
+ *
  * Like all other solver classes, this class has a local structure called
  * @p{AdditionalData} which is used to pass additional parameters to the
  * solver, like damping parameters or the number of temporary vectors. We
@@ -76,22 +79,25 @@ class SolverRichardson : public Subscriptor, private Solver<VECTOR>
     virtual ~SolverRichardson ();
     
                                     /**
-                                     * Solve $Ax=b$ for $x$.
+                                     * Solve the linear system $Ax=b$
+                                     * for x.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    void solve (const MATRIX &A,
-               VECTOR       &x,
-               const VECTOR &b,
-               const PRECONDITIONER& precondition);
+    void
+    solve (const MATRIX         &A,
+          VECTOR               &x,
+          const VECTOR         &b,
+          const PRECONDITIONER &precondition);
 
                                     /**
                                      * Solve $A^Tx=b$ for $x$.
                                      */
     template<class MATRIX, class PRECONDITIONER>
-    void Tsolve (const MATRIX &A,
-                VECTOR       &x,
-                const VECTOR &b,
-                const PRECONDITIONER& precondition);
+    void
+    Tsolve (const MATRIX         &A,
+           VECTOR               &x,
+           const VECTOR         &b,
+           const PRECONDITIONER &precondition);
 
                                     /**
                                      * Set the damping-coefficient.
@@ -169,10 +175,10 @@ SolverRichardson<VECTOR>::~SolverRichardson()
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
 void
-SolverRichardson<VECTOR>::solve (const MATRIX &A,
-                                VECTOR       &x,
-                                const VECTOR &b,
-                                const PRECONDITIONERprecondition)
+SolverRichardson<VECTOR>::solve (const MATRIX         &A,
+                                VECTOR               &x,
+                                const VECTOR         &b,
+                                const PRECONDITIONER &precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 
@@ -218,10 +224,10 @@ SolverRichardson<VECTOR>::solve (const MATRIX &A,
 template<class VECTOR>
 template<class MATRIX, class PRECONDITIONER>
 void
-SolverRichardson<VECTOR>::Tsolve (const MATRIX &A,
-                                 VECTOR       &x,
-                                 const VECTOR &b,
-                                 const PRECONDITIONERprecondition)
+SolverRichardson<VECTOR>::Tsolve (const MATRIX         &A,
+                                 VECTOR               &x,
+                                 const VECTOR         &b,
+                                 const PRECONDITIONER &precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 

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.