]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Solver classes without matrix template
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 25 Apr 2000 19:08:00 +0000 (19:08 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 25 Apr 2000 19:08:00 +0000 (19:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@2770 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/eigen.h
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
deal.II/lac/include/lac/solver_selector.h
tests/lac/mg.cc

index e12dcdea7270703295b6e0a0970a451fdfd8e4f1..e9489ccae5fe0a11cbcbf9453da4f05bdea88adb 100644 (file)
@@ -89,9 +89,8 @@ class ShiftedMatrix
  *
  * @author Guido Kanschat, 2000
  */
-template <class MATRIX = SparseMatrix<double>,
-          class VECTOR = Vector<double> >
-class EigenPower : private Solver<MATRIX,VECTOR>
+template <class VECTOR = Vector<double> >
+class EigenPower : private Solver<VECTOR>
 {
   public:
                                     /**
@@ -139,7 +138,8 @@ class EigenPower : private Solver<MATRIX,VECTOR>
                                      * corresponding eigenvector,
                                      * normalized with respect to the l2-norm.
                                      */
-    typename Solver<MATRIX,VECTOR>::ReturnState
+    template <class MATRIX>
+    typename Solver<VECTOR>::ReturnState
     solve (double       &value,
           const MATRIX &A,
           VECTOR       &x);
@@ -158,21 +158,38 @@ class EigenPower : private Solver<MATRIX,VECTOR>
  *
  * This class implements an adaptive version of the inverse iteration by Wieland.
  *
+ * There are two choices for the stopping criterion: by default, the
+ * norm of the residual $A x - l x$ is computed. Since this might not
+ * converge to zero for non-symmetric matrices with non-trivial Jordan
+ * blocks, it can be replaced by checking the difference of successive
+ * eigenvalues. Use @p{AdditionalData::use_residual} for switching
+ * these options.
+ *
  * @author Guido Kanschat, 2000
  */
-template <class MATRIX = SparseMatrix<double>,
-          class VECTOR = Vector<double> >
-class EigenInverse : private Solver<MATRIX, VECTOR>
+template <class VECTOR = Vector<double> >
+class EigenInverse : private Solver<VECTOR>
 {
   public:
                                     /**
                                      * Standardized data struct to
                                      * pipe additional data to the
-                                     * solver. This solver does not
-                                     * need additional data yet.
+                                     * solver.
                                      */
     struct AdditionalData
-    {};
+    {
+                                        /**
+                                         * Flag for the stopping criterion.
+                                         */
+       bool use_residual;
+                                        /**
+                                         * Constructor.
+                                         */
+       AdditionalData (bool use_residual = true):
+                       use_residual(use_residual)
+         {}
+       
+    };
     
                                     /**
                                      * Constructor.
@@ -200,7 +217,8 @@ class EigenInverse : private Solver<MATRIX, VECTOR>
                                      * normalized with respect to the
                                      * l2-norm.
                                      */
-    typename Solver<MATRIX,VECTOR>::ReturnState
+    template <class MATRIX>
+    typename Solver<VECTOR>::ReturnState
     solve (double       &value,
           const MATRIX &A,
           VECTOR       &x);
@@ -267,25 +285,26 @@ ShiftedMatrix<MATRIX>::residual (VECTOR& dst,
 //----------------------------------------------------------------------
 
 
-template <class MATRIX, class VECTOR>
-EigenPower<MATRIX, VECTOR>::EigenPower (SolverControl &cn,
-                                       VectorMemory<VECTOR> &mem,
-                                       const AdditionalData &data):
-               Solver<MATRIX, VECTOR>(cn, mem),
+template <class VECTOR>
+EigenPower<VECTOR>::EigenPower (SolverControl &cn,
+                               VectorMemory<VECTOR> &mem,
+                               const AdditionalData &data):
+               Solver<VECTOR>(cn, mem),
                additional_data(data)
 {}
 
 
-template <class MATRIX, class VECTOR>
-EigenPower<MATRIX, VECTOR>::~EigenPower ()
+template <class VECTOR>
+EigenPower<VECTOR>::~EigenPower ()
 {}
 
 
-template <class MATRIX, class VECTOR>
-typename Solver<MATRIX,VECTOR>::ReturnState
-EigenPower<MATRIX, VECTOR>::solve (double       &value,
-                                  const MATRIX &A,
-                                  VECTOR       &x)
+template <class VECTOR>
+template <class MATRIX>
+typename Solver<VECTOR>::ReturnState
+EigenPower<VECTOR>::solve (double       &value,
+                          const MATRIX &A,
+                          VECTOR       &x)
 {
   SolverControl::State conv=SolverControl::iterate;
 
@@ -350,25 +369,26 @@ EigenPower<MATRIX, VECTOR>::solve (double       &value,
 
 //----------------------------------------------------------------------//
 
-template <class MATRIX, class VECTOR>
-EigenInverse<MATRIX, VECTOR>::EigenInverse (SolverControl &cn,
+template <class VECTOR>
+EigenInverse<VECTOR>::EigenInverse (SolverControl &cn,
                                            VectorMemory<VECTOR> &mem,
                                            const AdditionalData &data):
-               Solver<MATRIX, VECTOR>(cn, mem),
+               Solver<VECTOR>(cn, mem),
                additional_data(data)
 {}
 
 
-template <class MATRIX, class VECTOR>
-EigenInverse<MATRIX, VECTOR>::~EigenInverse ()
+template <class VECTOR>
+EigenInverse<VECTOR>::~EigenInverse ()
 {}
 
 
-template <class MATRIX, class VECTOR>
-typename Solver<MATRIX,VECTOR>::ReturnState
-EigenInverse<MATRIX, VECTOR>::solve (double       &value,
-                                    const MATRIX &A,
-                                    VECTOR       &x)
+template <class VECTOR>
+template <class MATRIX>
+typename Solver<VECTOR>::ReturnState
+EigenInverse<VECTOR>::solve (double       &value,
+                            const MATRIX &A,
+                            VECTOR       &x)
 {
   deallog.push("Wieland");
 
@@ -380,7 +400,7 @@ EigenInverse<MATRIX, VECTOR>::solve (double       &value,
                                   // Define solver
   ReductionControl inner_control (100, 1.e-16, 1.e-8, false, false);
   PreconditionIdentity prec;
-  SolverQMRS<ShiftedMatrix <MATRIX>, VECTOR>
+  SolverCG<VECTOR>
     solver(inner_control, memory);
 
                                   // Next step for recomputing the shift
@@ -433,11 +453,15 @@ EigenInverse<MATRIX, VECTOR>::solve (double       &value,
                                       // Update normalized eigenvector
       x.equ (1./length, y);
                                       // Compute residual
-      y.equ (value, x);
-      double res = A.residual (r,x,y);
-
-                                      // Check the residual
-      conv = control().check (iter, res);
+      if (additional_data.use_residual)
+       {
+         y.equ (value, x);
+         double res = A.residual (r,x,y);
+                                          // Check the residual
+         conv = control().check (iter, res);
+       } else {
+         conv = control().check (iter, fabs(1.-old_value/value));
+       }
       old_value = value;
     }
 
index 4c9fc9964c58482845f81a6699affe7669c2e321..88156952e9633d260ac0c09a0258df884d54077b 100644 (file)
  *
  * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997, 1998, 1999
  */
-template <class Matrix = SparseMatrix<double>,
-          class Vector = Vector<double> >
+template <class Vector = Vector<double> >
 class Solver
 {
   public:
@@ -179,17 +178,18 @@ class Solver
 
 /*-------------------------------- Inline functions ------------------------*/
 
-template <class Matrix, class Vector>
+template <class Vector>
 inline
-SolverControl & Solver<Matrix,Vector>::control() const
+SolverControl &
+Solver<Vector>::control() const
 {
   return cntrl;
 };
 
 
-template<class Matrix, class Vector>
+template<class Vector>
 inline
-Solver<Matrix, Vector>::Solver(SolverControl &cn, VectorMemory<Vector> &mem)
+Solver<Vector>::Solver(SolverControl &cn, VectorMemory<Vector> &mem)
                : cntrl(cn),
                  memory(mem)
 {};
index 3e6020b3a1890dee92df0e839d6559d97c35c54e..c1838ed788004bbec674bbcfb030287bbe2607d0 100644 (file)
@@ -36,8 +36,8 @@
  * has a default argument, so you may call it without the additional
  * parameter.
  */
-template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverBicgstab : public Solver<Matrix,Vector>
+template <class VECTOR = Vector<double> >
+class SolverBicgstab : private Solver<VECTOR>
 {
   public:
                                     /**
@@ -52,7 +52,7 @@ class SolverBicgstab : public Solver<Matrix,Vector>
                                      * Constructor.
                                      */
     SolverBicgstab (SolverControl &cn,
-                   VectorMemory<Vector> &mem,
+                   VectorMemory<VECTOR> &mem,
                    const AdditionalData &data=AdditionalData());
 
                                     /**
@@ -63,18 +63,19 @@ class SolverBicgstab : public Solver<Matrix,Vector>
                                     /**
                                      * Solve primal problem only.
                                      */
-    template<class Preconditioner>
-    typename Solver<Matrix,Vector>::ReturnState
-    solve (const Matrix &A,
-          Vector       &x,
-          const Vector &b,
+    template<class MATRIX, class Preconditioner>
+    typename Solver<VECTOR>::ReturnState
+    solve (const MATRIX &A,
+          VECTOR       &x,
+          const VECTOR &b,
           const Preconditioner& precondition);
 
   protected:
                                     /**
                                      * Computation of the stopping criterion.
                                      */
-    virtual double criterion();
+    template <class MATRIX>
+    double criterion (const MATRIX& A, const VECTOR& x, const VECTOR& b);
 
                                     /**
                                      * Interface for derived class.
@@ -86,55 +87,50 @@ class SolverBicgstab : public Solver<Matrix,Vector>
                                      * convergence history.
                                      */
     virtual void print_vectors(const unsigned int step,
-                              const Vector& x,
-                              const Vector& r,
-                              const Vector& d) const;
+                              const VECTOR& x,
+                              const VECTOR& r,
+                              const VECTOR& d) const;
 
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector *Vx;
+    VECTOR *Vx;
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector *Vr;
+    VECTOR *Vr;
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector *Vrbar;
+    VECTOR *Vrbar;
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector *Vp;
+    VECTOR *Vp;
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector *Vy;
+    VECTOR *Vy;
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector *Vz;
+    VECTOR *Vz;
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector *Vs;
+    VECTOR *Vs;
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector *Vt;
+    VECTOR *Vt;
                                     /**
                                      * Auxiliary vector.
                                      */
-    Vector *Vv;
+    VECTOR *Vv;
                                     /**
                                      * Right hand side vector.
                                      */
-    const Vector *Vb;
-  
-                                    /**
-                                     * Pointer to the system matrix.
-                                     */
-    const Matrix *MA;
+    const VECTOR *Vb;
   
                                     /**
                                      * Auxiliary value.
@@ -171,50 +167,53 @@ class SolverBicgstab : public Solver<Matrix,Vector>
                                     /**
                                      * Everything before the iteration loop.
                                      */
-    SolverControl::State start();
+    template <class MATRIX>
+    SolverControl::State start(const MATRIX& A);
 
                                     /**
                                      * The iteration loop itself.
                                      */
-    template<class Preconditioner>
-    typename Solver<Matrix,Vector>::ReturnState 
-    iterate(const Preconditioner& precondition);
+    template<class MATRIX, class PRECONDITIONER>
+    typename Solver<VECTOR>::ReturnState 
+    iterate(const MATRIX& A, const PRECONDITIONER& precondition);
   
 };
 
 /*-------------------------Inline functions -------------------------------*/
 
-template<class Matrix, class Vector>
-SolverBicgstab<Matrix, Vector>::SolverBicgstab (SolverControl &cn,
-                                              VectorMemory<Vector> &mem,
-                                              const AdditionalData &)
+template<class VECTOR>
+SolverBicgstab<VECTOR>::SolverBicgstab (SolverControl &cn,
+                                       VectorMemory<VECTOR> &mem,
+                                       const AdditionalData &)
                :
-               Solver<Matrix,Vector>(cn,mem)
+               Solver<VECTOR>(cn,mem)
 {}
 
 
 
-template<class Matrix, class Vector>
-SolverBicgstab<Matrix, Vector>::~SolverBicgstab ()
+template<class VECTOR>
+SolverBicgstab<VECTOR>::~SolverBicgstab ()
 {}
 
 
 
-template<class Matrix, class Vector>
+template <class VECTOR>
+template <class MATRIX>
 double
-SolverBicgstab<Matrix, Vector>::criterion()
+SolverBicgstab<VECTOR>::criterion (const MATRIX& A, const VECTOR& x, const VECTOR& b)
 {
-  res = MA->residual(*Vt, *Vx, *Vb);
+  res = A.residual(*Vt, x, b);
   return res;
 }
 
 
 
-template < class Matrix, class Vector >
+template <class VECTOR >
+template <class MATRIX>
 SolverControl::State
-SolverBicgstab<Matrix, Vector>::start()
+SolverBicgstab<VECTOR>::start(const MATRIX& A)
 {
-  res = MA->residual(*Vr, *Vx, *Vb);
+  res = A.residual(*Vr, *Vx, *Vb);
   Vp->reinit(*Vx);
   Vv->reinit(*Vx);
   Vrbar->equ(1.,*Vr);
@@ -224,32 +223,33 @@ SolverBicgstab<Matrix, Vector>::start()
 
 
 
-template<class Matrix, class Vector>
+template<class VECTOR>
 void
-SolverBicgstab<Matrix,Vector>::print_vectors(const unsigned int,
-                                            const Vector&,
-                                            const Vector&,
-                                            const Vector&) const
+SolverBicgstab<VECTOR>::print_vectors(const unsigned int,
+                                     const VECTOR&,
+                                     const VECTOR&,
+                                     const VECTOR&) const
 {}
 
 
 
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverBicgstab<Matrix, Vector>::iterate(const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class PRECONDITIONER>
+typename Solver<VECTOR>::ReturnState
+SolverBicgstab<VECTOR>::iterate(const MATRIX& A,
+                               const PRECONDITIONER& precondition)
 {
   SolverControl::State state = SolverControl::iterate;
   alpha = omega = rho = 1.;
 
-  Vector& r = *Vr;
-  Vector& rbar = *Vrbar;
-  Vector& p = *Vp;
-  Vector& y = *Vy;
-  Vector& z = *Vz;
-  Vector& s = *Vs;
-  Vector& t = *Vt;
-  Vector& v = *Vv;
+  VECTOR& r = *Vr;
+  VECTOR& rbar = *Vrbar;
+  VECTOR& p = *Vp;
+  VECTOR& y = *Vy;
+  VECTOR& z = *Vz;
+  VECTOR& s = *Vs;
+  VECTOR& t = *Vt;
+  VECTOR& v = *Vv;
   
   do
     {
@@ -258,7 +258,7 @@ SolverBicgstab<Matrix, Vector>::iterate(const Preconditioner& precondition)
       rho    = rhobar;
       p.sadd(beta, 1., r, -beta*omega, v);
       precondition(y,p);
-      MA->vmult(v,y);
+      A.vmult(v,y);
       rhobar = rbar * v;
 
       alpha = rho/rhobar;
@@ -269,13 +269,13 @@ SolverBicgstab<Matrix, Vector>::iterate(const Preconditioner& precondition)
     
       s.equ(1., r, -alpha, v);
       precondition(z,s);
-      MA->vmult(t,z);
+      A.vmult(t,z);
       rhobar = t*s;
       omega = rhobar/(t*t);
       Vx->add(alpha, y, omega, z);
       r.equ(1., s, -omega, t);
 
-      state = control().check(++step, criterion());
+      state = control().check(++step, criterion(A, *Vx, *Vb));
       print_vectors(step, *Vx, r, y);
     }
   while (state == SolverControl::iterate);
@@ -284,13 +284,13 @@ SolverBicgstab<Matrix, Vector>::iterate(const Preconditioner& precondition)
 }
 
 
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState
-SolverBicgstab<Matrix, Vector>::solve(const Matrix &A,
-                                     Vector       &x,
-                                     const Vector &b,
-                                     const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class PRECONDITIONER>
+typename Solver<VECTOR>::ReturnState
+SolverBicgstab<VECTOR>::solve(const MATRIX &A,
+                             VECTOR       &x,
+                             const VECTOR &b,
+                             const PRECONDITIONER& precondition)
 {
   deallog.push("Bicgstab");
   Vr    = memory.alloc(); Vr->reinit(x);
@@ -302,7 +302,6 @@ SolverBicgstab<Matrix, Vector>::solve(const Matrix &A,
   Vt    = memory.alloc(); Vt->reinit(x);
   Vv    = memory.alloc();
 
-  MA = &A;
   Vx = &x;
   Vb = &b;
 
@@ -314,8 +313,8 @@ SolverBicgstab<Matrix, Vector>::solve(const Matrix &A,
     {
       if (step)
        deallog << "Restart step " << step << endl;
-      if (start() == SolverControl::success) break;  
-      state = iterate(precondition);
+      if (start(A) == SolverControl::success) break;  
+      state = iterate(A, precondition);
     }
   while (state == breakdown);
 
index 5a30403b45d958af0493002f9afd0f27226cbe52..81c49d474b9829ac18a19a323f4f92c55fd2b698 100644 (file)
@@ -43,8 +43,8 @@
  *
  * @author Original implementation by G. Kanschat, R. Becker and F.-T. Suttmeier, reworking and  documentation by Wolfgang Bangerth
  */
-template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverCG : public Solver<Matrix,Vector>
+template <class VECTOR = Vector<double> >
+class SolverCG : private Solver<VECTOR>
 {
   public:
                                     /**
@@ -59,7 +59,7 @@ class SolverCG : public Solver<Matrix,Vector>
                                      * Constructor.
                                      */
     SolverCG (SolverControl &cn,
-             VectorMemory<Vector> &mem,
+             VectorMemory<VECTOR> &mem,
              const AdditionalData &data=AdditionalData());
 
                                     /**
@@ -70,11 +70,11 @@ class SolverCG : public Solver<Matrix,Vector>
                                     /**
                                      * Solver method.
                                      */
-    template<class Preconditioner>
-    typename Solver<Matrix,Vector>::ReturnState
-    solve (const Matrix &A,
-          Vector       &x,
-          const Vector &b,
+    template<class MATRIX, class Preconditioner>
+    typename Solver<VECTOR>::ReturnState
+    solve (const MATRIX &A,
+          VECTOR       &x,
+          const VECTOR &b,
           const Preconditioner& precondition);
 
   protected:
@@ -96,9 +96,9 @@ class SolverCG : public Solver<Matrix,Vector>
                                      * convergence history.
                                      */
     virtual void print_vectors(const unsigned int step,
-                              const Vector& x,
-                              const Vector& r,
-                              const Vector& d) const;
+                              const VECTOR& x,
+                              const VECTOR& r,
+                              const VECTOR& d) const;
 
                                     /**
                                      * Temporary vectors, allocated through
@@ -106,10 +106,10 @@ class SolverCG : public Solver<Matrix,Vector>
                                      * of the actual solution process and
                                      * deallocated at the end.
                                      */
-    Vector *Vr;
-    Vector *Vp;
-    Vector *Vz;
-    Vector *VAp;
+    VECTOR *Vr;
+    VECTOR *Vp;
+    VECTOR *Vz;
+    VECTOR *VAp;
     
                                     /**
                                      * Within the iteration loop, the
@@ -128,48 +128,48 @@ class SolverCG : public Solver<Matrix,Vector>
 /*------------------------- Implementation ----------------------------*/
 
 
-template<class Matrix, class Vector>
-SolverCG<Matrix,Vector>::SolverCG(SolverControl &cn,
-                                 VectorMemory<Vector> &mem,
-                                 const AdditionalData &)
+template<class VECTOR>
+SolverCG<VECTOR>::SolverCG(SolverControl &cn,
+                          VectorMemory<VECTOR> &mem,
+                          const AdditionalData &)
                :
-               Solver<Matrix,Vector>(cn,mem)
+               Solver<VECTOR>(cn,mem)
 {}
 
 
 
-template<class Matrix, class Vector>
-SolverCG<Matrix,Vector>::~SolverCG ()
+template<class VECTOR>
+SolverCG<VECTOR>::~SolverCG ()
 {}
 
 
 
-template<class Matrix, class Vector>
+template<class VECTOR>
 double
-SolverCG<Matrix,Vector>::criterion()
+SolverCG<VECTOR>::criterion()
 {
   return sqrt(res2);
 }
 
 
 
-template<class Matrix, class Vector>
+template<class VECTOR>
 void
-SolverCG<Matrix,Vector>::print_vectors(const unsigned int,
-                                      const Vector&,
-                                      const Vector&,
-                                      const Vector&) const
+SolverCG<VECTOR>::print_vectors(const unsigned int,
+                               const VECTOR&,
+                               const VECTOR&,
+                               const VECTOR&) const
 {}
 
 
 
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState 
-SolverCG<Matrix,Vector>::solve (const Matrix &A,
-                               Vector       &x,
-                               const Vector &b,
-                               const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class Preconditioner>
+typename Solver<VECTOR>::ReturnState 
+SolverCG<VECTOR>::solve (const MATRIX &A,
+                        VECTOR       &x,
+                        const VECTOR &b,
+                        const Preconditioner& precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 
@@ -181,10 +181,10 @@ SolverCG<Matrix,Vector>::solve (const Matrix &A,
   Vz  = memory.alloc();
   VAp = memory.alloc();
                                   // define some aliases for simpler access
-  Vector& g  = *Vr; 
-  Vector& h  = *Vp; 
-  Vector& d  = *Vz; 
-  Vector& Ad = *VAp;
+  VECTOR& g  = *Vr; 
+  VECTOR& h  = *Vp; 
+  VECTOR& d  = *Vz; 
+  VECTOR& Ad = *VAp;
                                   // resize the vectors, but do not set
                                   // the values since they'd be overwritten
                                   // soon anyway.
index d20e0bdb5f0e8f7a59a3638037aab6187d8e3a2c..2dbf439fe807be7ebcb5f5f73b3563dc933bb42a 100644 (file)
@@ -71,9 +71,8 @@
  *
  * @author Wolfgang Bangerth
  */
-template <class MATRIX = SparseMatrix<double>,
-          class VECTOR = Vector<double> >
-class SolverGMRES : public Solver<MATRIX, VECTOR>
+template <class VECTOR = Vector<double> >
+class SolverGMRES : private Solver<VECTOR>
 {
   public:
                                     /**
@@ -109,11 +108,11 @@ class SolverGMRES : public Solver<MATRIX, VECTOR>
                                     /**
                                      * Solver method.
                                      */
-    template<class Preconditioner>
-    typename Solver<MATRIX,VECTOR>::ReturnState solve (const MATRIX &A,
+    template<class MATRIX, class PRECONDITIONER>
+    typename Solver<VECTOR>::ReturnState solve (const MATRIX &A,
                                                       VECTOR       &x,
                                                       const VECTOR &b,
-                                                      const Preconditioner& precondition);
+                                                      const PRECONDITIONER& precondition);
 
     DeclException1 (ExcTooFewTmpVectors,
                    int,
@@ -149,11 +148,11 @@ class SolverGMRES : public Solver<MATRIX, VECTOR>
 /* --------------------- Inline and template functions ------------------- */
 
 
-template <class MATRIX, class VECTOR>
-SolverGMRES<MATRIX,VECTOR>::SolverGMRES (SolverControl        &cn,
+template <class VECTOR>
+SolverGMRES<VECTOR>::SolverGMRES (SolverControl        &cn,
                                         VectorMemory<VECTOR> &mem,
                                         const AdditionalData &data) :
-               Solver<MATRIX,VECTOR> (cn,mem),
+               Solver<VECTOR> (cn,mem),
                additional_data(data)
 {
   Assert (additional_data.max_n_tmp_vectors >= 10, 
@@ -161,14 +160,14 @@ SolverGMRES<MATRIX,VECTOR>::SolverGMRES (SolverControl        &cn,
 };
 
 
-template <class MATRIX, class VECTOR>
+template <class VECTOR>
 inline
 void
-SolverGMRES<MATRIX,VECTOR>::givens_rotation (Vector<double> &h,
-                                            Vector<double> &b,
-                                            Vector<double> &ci,
-                                            Vector<double> &si, 
-                                            int     col) const
+SolverGMRES<VECTOR>::givens_rotation (Vector<double> &h,
+                                     Vector<double> &b,
+                                     Vector<double> &ci,
+                                     Vector<double> &si, 
+                                     int     col) const
 {
   for (int i=0 ; i<col ; i++)
     {
@@ -188,13 +187,13 @@ SolverGMRES<MATRIX,VECTOR>::givens_rotation (Vector<double> &h,
 }
 
 
-template<class MATRIX, class VECTOR>
-template<class Preconditioner>
-typename Solver<MATRIX,VECTOR>::ReturnState
-SolverGMRES<MATRIX,VECTOR>::solve (const MATRIX& A,
+template<class VECTOR>
+template<class MATRIX, class PRECONDITIONER>
+typename Solver<VECTOR>::ReturnState
+SolverGMRES<VECTOR>::solve (const MATRIX& A,
                                   VECTOR      & x,
                                   const VECTOR& b,
-                                  const Preconditioner& precondition)
+                                  const PRECONDITIONER& precondition)
 {
                                   // this code was written by the fathers of
                                   // DEAL. I take absolutely no guarantees
@@ -400,9 +399,9 @@ unsigned int dim = 0;
 };
 
 
-template<class MATRIX, class VECTOR>
+template<class VECTOR>
 double
-SolverGMRES<MATRIX,VECTOR>::criterion () 
+SolverGMRES<VECTOR>::criterion () 
 {
                                   // dummy implementation. this function is
                                   // not needed for the present implementation
index 0ef4b63d1b9db7742913be830a81943a054b1507..6bf5e26a356b576d42010512266f72ce3d53e9e0 100644 (file)
@@ -47,8 +47,8 @@
  *
  * @author Thomas Richter, 2000
  */
-template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverMinRes : public Solver<Matrix,Vector>
+template <class VECTOR = Vector<double> >
+class SolverMinRes : private Solver<VECTOR>
 {
   public:
                                     /**
@@ -63,7 +63,7 @@ class SolverMinRes : public Solver<Matrix,Vector>
                                      * Constructor.
                                      */
     SolverMinRes (SolverControl &cn,
-                 VectorMemory<Vector> &mem,
+                 VectorMemory<VECTOR> &mem,
                  const AdditionalData &data=AdditionalData());
 
                                     /**
@@ -74,11 +74,11 @@ class SolverMinRes : public Solver<Matrix,Vector>
                                     /**
                                      * Solver method.
                                      */
-    template<class Preconditioner>
-    typename Solver<Matrix,Vector>::ReturnState
-    solve (const Matrix &A,
-          Vector       &x,
-          const Vector &b,
+    template<class MATRIX, class Preconditioner>
+    typename Solver<VECTOR>::ReturnState
+    solve (const MATRIX &A,
+          VECTOR       &x,
+          const VECTOR &b,
           const Preconditioner& precondition);
 
                                     /**
@@ -103,9 +103,9 @@ class SolverMinRes : public Solver<Matrix,Vector>
                                      * convergence history.
                                      */
     virtual void print_vectors(const unsigned int step,
-                              const Vector& x,
-                              const Vector& r,
-                              const Vector& d) const;
+                              const VECTOR& x,
+                              const VECTOR& r,
+                              const VECTOR& d) const;
 
                                     /**
                                      * Temporary vectors, allocated through
@@ -113,9 +113,9 @@ class SolverMinRes : public Solver<Matrix,Vector>
                                      * of the actual solution process and
                                      * deallocated at the end.
                                      */
-    Vector *Vu0, *Vu1, *Vu2;
-    Vector *Vm0, *Vm1, *Vm2;   
-    Vector *Vv;
+    VECTOR *Vu0, *Vu1, *Vu2;
+    VECTOR *Vm0, *Vm1, *Vm2;   
+    VECTOR *Vv;
     
                                     /**
                                      * Within the iteration loop, the
@@ -134,46 +134,46 @@ class SolverMinRes : public Solver<Matrix,Vector>
 /*------------------------- Implementation ----------------------------*/
 
 
-template<class Matrix, class Vector>
-SolverMinRes<Matrix,Vector>::SolverMinRes (SolverControl &cn,
-                                         VectorMemory<Vector> &mem,
-                                         const AdditionalData &)
+template<class VECTOR>
+SolverMinRes<VECTOR>::SolverMinRes (SolverControl &cn,
+                                   VectorMemory<VECTOR> &mem,
+                                   const AdditionalData &)
                :
-               Solver<Matrix,Vector>(cn,mem)
+               Solver<VECTOR>(cn,mem)
 {}
 
 
-template<class Matrix, class Vector>
-SolverMinRes<Matrix,Vector>::~SolverMinRes ()
+template<class VECTOR>
+SolverMinRes<VECTOR>::~SolverMinRes ()
 {}
 
 
 
-template<class Matrix, class Vector>
+template<class VECTOR>
 double
-SolverMinRes<Matrix,Vector>::criterion()
+SolverMinRes<VECTOR>::criterion()
 {
   return res2;
 };
 
 
-template<class Matrix, class Vector>
+template<class VECTOR>
 void
-SolverMinRes<Matrix,Vector>::print_vectors(const unsigned int,
-                                          const Vector&,
-                                          const Vector&,
-                                          const Vector&) const
+SolverMinRes<VECTOR>::print_vectors(const unsigned int,
+                                   const VECTOR&,
+                                   const VECTOR&,
+                                   const VECTOR&) const
 {}
 
 
 
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState 
-SolverMinRes<Matrix,Vector>::solve (const Matrix &A,
-                                   Vector       &x,
-                                   const Vector &b,
-                                   const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class Preconditioner>
+typename Solver<VECTOR>::ReturnState 
+SolverMinRes<VECTOR>::solve (const MATRIX &A,
+                            VECTOR       &x,
+                            const VECTOR &b,
+                            const Preconditioner& precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 
@@ -192,7 +192,7 @@ SolverMinRes<Matrix,Vector>::solve (const Matrix &A,
   Vm1  = memory.alloc();
   Vm2  = memory.alloc();
                                   // define some aliases for simpler access
-  typedef Vector vecref;
+  typedef VECTOR vecref;
   vecref u[3] = {*Vu0, *Vu1, *Vu2};
   vecref m[3] = {*Vm0, *Vm1, *Vm2};
   vecref v    = *Vv;
index e41631fdd78d9a063419a3eb083e5d6aa4dacd4b..c6f181d570eb9e9fa936508f917d60ab96f2eb05 100644 (file)
@@ -48,8 +48,8 @@
  *
  * @author Guido Kanschat, 1999
  */
-template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverQMRS : public Solver<Matrix,Vector>
+template <class VECTOR = Vector<double> >
+class SolverQMRS : private Solver<VECTOR>
 {
   public:
                                     /**
@@ -96,17 +96,17 @@ class SolverQMRS : public Solver<Matrix,Vector>
                                      * Constructor.
                                      */
     SolverQMRS (SolverControl &cn,
-             VectorMemory<Vector> &mem,
+             VectorMemory<VECTOR> &mem,
              const AdditionalData &data=AdditionalData());
 
                                     /**
                                      * Solver method.
                                      */
-    template<class Preconditioner>
-    typename Solver<Matrix,Vector>::ReturnState
-    solve (const Matrix &A,
-                      Vector       &x,
-                      const Vector &b,
+    template<class MATRIX, class Preconditioner>
+    typename Solver<VECTOR>::ReturnState
+    solve (const MATRIX &A,
+                      VECTOR       &x,
+                      const VECTOR &b,
                       const Preconditioner& precondition);
 
                                     /**
@@ -119,9 +119,9 @@ class SolverQMRS : public Solver<Matrix,Vector>
                                      * convergence history.
                                      */
     virtual void print_vectors(const unsigned int step,
-                              const Vector& x,
-                              const Vector& r,
-                              const Vector& d) const;
+                              const VECTOR& x,
+                              const VECTOR& r,
+                              const VECTOR& d) const;
    protected:
                                     /**
                                      * Implementation of the computation of
@@ -135,24 +135,20 @@ class SolverQMRS : public Solver<Matrix,Vector>
                                      * of the actual solution process and
                                      * deallocated at the end.
                                      */
-    Vector *Vv;
-    Vector *Vp;
-    Vector *Vq;
-    Vector *Vt;
-    Vector *Vd;
+    VECTOR *Vv;
+    VECTOR *Vp;
+    VECTOR *Vq;
+    VECTOR *Vt;
+    VECTOR *Vd;
                                     /**
                                      * Iteration vector.
                                      */
-    Vector *Vx;
+    VECTOR *Vx;
                                     /**
                                      * RHS vector.
                                      */
-    const Vector *Vb;
+    const VECTOR *Vb;
     
-                                    /**
-                                     * Pointer to the matrix to be inverted.
-                                     */
-    const Matrix* MA;
                                     /**
                                      * Within the iteration loop, the
                                      * square of the residual vector is
@@ -172,9 +168,9 @@ class SolverQMRS : public Solver<Matrix,Vector>
                                     /**
                                      * The iteration loop itself.
                                      */
-    template<class Preconditioner>
-    typename Solver<Matrix,Vector>::ReturnState 
-    iterate(const Preconditioner& precondition);
+    template<class MATRIX, class Preconditioner>
+    typename Solver<VECTOR>::ReturnState 
+    iterate(const MATRIX& A, const Preconditioner& precondition);
                                     /**
                                      * The current iteration step.
                                      */
@@ -185,41 +181,41 @@ class SolverQMRS : public Solver<Matrix,Vector>
 /*------------------------- Implementation ----------------------------*/
 
 
-template<class Matrix, class Vector>
-SolverQMRS<Matrix,Vector>::SolverQMRS(SolverControl &cn,
-                                 VectorMemory<Vector> &mem,
-                                 const AdditionalData &data) :
-               Solver<Matrix,Vector>(cn,mem),
+template<class VECTOR>
+SolverQMRS<VECTOR>::SolverQMRS(SolverControl &cn,
+                              VectorMemory<VECTOR> &mem,
+                              const AdditionalData &data) :
+               Solver<VECTOR>(cn,mem),
                additional_data(data)
 {};
 
 
-template<class Matrix, class Vector>
+template<class VECTOR>
 double
-SolverQMRS<Matrix,Vector>::criterion()
+SolverQMRS<VECTOR>::criterion()
 {
   return sqrt(res2);
 }
 
 
 
-template<class Matrix, class Vector>
+template<class VECTOR>
 void
-SolverQMRS<Matrix,Vector>::print_vectors(const unsigned int,
-                                              const Vector&,
-                                              const Vector&,
-                                              const Vector&) const
+SolverQMRS<VECTOR>::print_vectors(const unsigned int,
+                                 const VECTOR&,
+                                 const VECTOR&,
+                                 const VECTOR&) const
 {}
 
 
 
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState 
-SolverQMRS<Matrix,Vector>::solve (const Matrix &A,
-                                 Vector       &x,
-                                 const Vector &b,
-                                 const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class Preconditioner>
+typename Solver<VECTOR>::ReturnState 
+SolverQMRS<VECTOR>::solve (const MATRIX &A,
+                          VECTOR       &x,
+                          const VECTOR &b,
+                          const Preconditioner& precondition)
 {
   deallog.push("QMRS");
   
@@ -230,7 +226,6 @@ SolverQMRS<Matrix,Vector>::solve (const Matrix &A,
   Vt  = memory.alloc();
   Vd  = memory.alloc();
 
-  MA = &A;
   Vx = &x;
   Vb = &b;
                                   // resize the vectors, but do not set
@@ -249,7 +244,7 @@ SolverQMRS<Matrix,Vector>::solve (const Matrix &A,
     {
       if (step)
        deallog << "Restart step " << step << endl;
-      state = iterate(precondition);
+      state = iterate(A, precondition);
     }
   while (state == breakdown);
     
@@ -268,10 +263,11 @@ SolverQMRS<Matrix,Vector>::solve (const Matrix &A,
   return state;
 };
 
-template<class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState 
-SolverQMRS<Matrix,Vector>::iterate(const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class Preconditioner>
+typename Solver<VECTOR>::ReturnState 
+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.
@@ -282,15 +278,13 @@ SolverQMRS<Matrix,Vector>::iterate(const Preconditioner& precondition)
   SolverControl::State state = SolverControl::iterate;
 
                                   // define some aliases for simpler access
-  Vector& v  = *Vv; 
-  Vector& p  = *Vp; 
-  Vector& q  = *Vq; 
-  Vector& t  = *Vt;
-  Vector& d  = *Vd;
-  Vector& x  = *Vx;
-  const Vector& b = *Vb;
-  
-  const Matrix& A  = *MA;
+  VECTOR& v  = *Vv; 
+  VECTOR& p  = *Vp; 
+  VECTOR& q  = *Vq; 
+  VECTOR& t  = *Vt;
+  VECTOR& d  = *Vd;
+  VECTOR& x  = *Vx;
+  const VECTOR& b = *Vb;
   
   int  it=0;
  
index e9ae82f5148384c2fe7d1331c4ac61c92a986e1b..eb5a6332e955a86fe87804d4823aef9b68185499 100644 (file)
@@ -37,8 +37,8 @@
  *
  * @author Ralf Hartmann
  */
-template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
-class SolverRichardson : public Solver<Matrix, Vector>
+template <class VECTOR = Vector<double> >
+class SolverRichardson : private Solver<VECTOR>
 {
   public:
                                     /**
@@ -67,16 +67,16 @@ class SolverRichardson : public Solver<Matrix, Vector>
                                      * Constructor.
                                      */
     SolverRichardson (SolverControl &cn,
-                     VectorMemory<Vector> &mem,
+                     VectorMemory<VECTOR> &mem,
                      const AdditionalData &data=AdditionalData());
 
                                     /**
                                      * Solve $Ax=b$ for $x$.
                                      */
-    template<class Preconditioner>
-    typename Solver<Matrix,Vector>::ReturnState solve (const Matrix &A,
-                                                      Vector       &x,
-                                                      const Vector &b,
+    template<class MATRIX, class Preconditioner>
+    typename Solver<VECTOR>::ReturnState solve (const MATRIX &A,
+                                                      VECTOR       &x,
+                                                      const VECTOR &b,
                                                       const Preconditioner& precondition);
 
                                     /**
@@ -95,9 +95,9 @@ class SolverRichardson : public Solver<Matrix, Vector>
                                      * convergence history.
                                      */
     virtual void print_vectors(const unsigned int step,
-                              const Vector& x,
-                              const Vector& r,
-                              const Vector& d) const;
+                              const VECTOR& x,
+                              const VECTOR& r,
+                              const VECTOR& d) const;
     
   protected:
                                     /**
@@ -112,8 +112,8 @@ class SolverRichardson : public Solver<Matrix, Vector>
                                      * of the actual solution process and
                                      * deallocated at the end.
                                      */
-    Vector *Vr;
-    Vector *Vd;
+    VECTOR *Vr;
+    VECTOR *Vd;
 
                                     /**
                                      * Damping parameter.
@@ -137,27 +137,27 @@ class SolverRichardson : public Solver<Matrix, Vector>
 /*----------------- Implementation of the Richardson Method ------------------*/
 
 
-template<class Matrix, class Vector>
-SolverRichardson<Matrix,Vector>::SolverRichardson(SolverControl &cn,
-                                                 VectorMemory<Vector> &mem,
-                                                 const AdditionalData &data):
-               Solver<Matrix,Vector> (cn,mem),
+template<class VECTOR>
+SolverRichardson<VECTOR>::SolverRichardson(SolverControl &cn,
+                                          VectorMemory<VECTOR> &mem,
+                                          const AdditionalData &data):
+               Solver<VECTOR> (cn,mem),
                additional_data(data)  {};
 
 
-template<class Matrix,class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState 
-SolverRichardson<Matrix,Vector>::solve (const Matrix &A,
-                                       Vector       &x,
-                                       const Vector &b,
-                                       const Preconditioner& precondition)
+template<class VECTOR>
+template<class MATRIX, class Preconditioner>
+typename Solver<VECTOR>::ReturnState 
+SolverRichardson<VECTOR>::solve (const MATRIX &A,
+                                VECTOR       &x,
+                                const VECTOR &b,
+                                const Preconditioner& precondition)
 {
   SolverControl::State conv=SolverControl::iterate;
 
                                   // Memory allocation
-  Vr  = memory.alloc(); Vector& r  = *Vr; r.reinit(x);
-  Vd  = memory.alloc(); Vector& d  = *Vd; d.reinit(x);
+  Vr  = memory.alloc(); VECTOR& r  = *Vr; r.reinit(x);
+  Vd  = memory.alloc(); VECTOR& d  = *Vd; d.reinit(x);
 
   deallog.push("Richardson");
 
@@ -188,27 +188,27 @@ SolverRichardson<Matrix,Vector>::solve (const Matrix &A,
 }
 
 
-template<class Matrix, class Vector>
+template<class VECTOR>
 void
-SolverRichardson<Matrix,Vector>::print_vectors(const unsigned int,
-                                              const Vector&,
-                                              const Vector&,
-                                              const Vector&) const
+SolverRichardson<VECTOR>::print_vectors(const unsigned int,
+                                       const VECTOR&,
+                                       const VECTOR&,
+                                       const VECTOR&) const
 {}
 
 
 
-template<class Matrix,class Vector>
+template<class VECTOR>
 inline double
-SolverRichardson<Matrix,Vector>::criterion()
+SolverRichardson<VECTOR>::criterion()
 {
   return res;
 }
 
 
-template<class Matrix,class Vector>
+template<class VECTOR>
 inline void
-SolverRichardson<Matrix,Vector>::set_omega(double om)
+SolverRichardson<VECTOR>::set_omega(double om)
 {
   additional_data.omega=om;
 }
index aa47cfbfc65beeffa631aef7a4fd0af4dcfe8113..1122e2a704674747451bf76b22ff0b921eb24718 100644 (file)
@@ -85,8 +85,7 @@
  * 
  * @author Ralf Hartmann, 1999
  */
-template <class Matrix = SparseMatrix<double>,
-          class Vector = Vector<double> >
+template <class Vector = Vector<double> >
 class SolverSelector
 {
   public:
@@ -112,8 +111,8 @@ class SolverSelector
                                      * was specified in the constructor.
                                      * 
                                      */
-    template<class Preconditioner>
-    typename Solver<Matrix,Vector>::ReturnState solve(const Matrix &A,
+    template<class Matrix, class Preconditioner>
+    typename Solver<Vector>::ReturnState solve(const Matrix &A,
                                                      Vector &x,
                                                      const Vector &b,
                                                      const Preconditioner &precond) const;
@@ -122,28 +121,28 @@ class SolverSelector
                                      * Set the additional data. For more
                                      * info see the #Solver# class.
                                      */
-    void set_data(const typename SolverRichardson<Matrix,Vector>
+    void set_data(const typename SolverRichardson<Vector>
                  ::AdditionalData &data);
 
                                     /**
                                      * Set the additional data. For more
                                      * info see the #Solver# class.
                                      */
-    void set_data(const typename SolverCG<Matrix,Vector>
+    void set_data(const typename SolverCG<Vector>
                  ::AdditionalData &data);
 
                                     /**
                                      * Set the additional data. For more
                                      * info see the #Solver# class.
                                      */
-    void set_data(const typename SolverBicgstab<Matrix,Vector>
+    void set_data(const typename SolverBicgstab<Vector>
                  ::AdditionalData &data); 
 
                                     /**
                                      * Set the additional data. For more
                                      * info see the #Solver# class.
                                      */
-    void set_data(const typename SolverGMRES<Matrix,Vector>
+    void set_data(const typename SolverGMRES<Vector>
                  ::AdditionalData &data);
 
                                     /**
@@ -183,116 +182,116 @@ class SolverSelector
                                     /**
                                      * Stores the additional data.
                                      */
-    typename SolverRichardson<Matrix,Vector>::AdditionalData richardson_data;
+    typename SolverRichardson<Vector>::AdditionalData richardson_data;
 
                                     /**
                                      * Stores the additional data.
                                      */
-    typename SolverCG<Matrix,Vector>::AdditionalData cg_data;
+    typename SolverCG<Vector>::AdditionalData cg_data;
 
                                     /**
                                      * Stores the additional data.
                                      */
-    typename SolverBicgstab<Matrix,Vector>::AdditionalData bicgstab_data;
+    typename SolverBicgstab<Vector>::AdditionalData bicgstab_data;
 
                                     /**
                                      * Stores the additional data.
                                      */
-    typename SolverGMRES<Matrix,Vector>::AdditionalData gmres_data;
+    typename SolverGMRES<Vector>::AdditionalData gmres_data;
 };
 
 
 /* --------------------- Inline and template functions ------------------- */
 
 
-template <class Matrix, class Vector>
-SolverSelector<Matrix, Vector>::SolverSelector(string solver_name,
-                                              SolverControl &control,
-                                              VectorMemory<Vector> &vector_memory) :
+template <class Vector>
+SolverSelector<Vector>::SolverSelector(string solver_name,
+                                      SolverControl &control,
+                                      VectorMemory<Vector> &vector_memory) :
                solver_name(solver_name),
                control(&control),
                vector_memory(&vector_memory)
 {};
 
 
-template <class Matrix, class Vector>
-SolverSelector<Matrix, Vector>::~SolverSelector()
+template <class Vector>
+SolverSelector<Vector>::~SolverSelector()
 {};
 
 
-template <class Matrix, class Vector>
-template<class Preconditioner>
-typename Solver<Matrix,Vector>::ReturnState 
-SolverSelector<Matrix, Vector>::solve(const Matrix &A,
-                                     Vector &x,
-                                     const Vector &b,
-                                     const Preconditioner &precond) const
+template <class Vector>
+template<class Matrix, class Preconditioner>
+typename Solver<Vector>::ReturnState 
+SolverSelector<Vector>::solve(const Matrix &A,
+                             Vector &x,
+                             const Vector &b,
+                             const Preconditioner &precond) const
 {
   if (solver_name=="richardson")
     {
-      SolverRichardson<Matrix,Vector> solver(*control,*vector_memory,
+      SolverRichardson<Vector> solver(*control,*vector_memory,
                                             richardson_data);
       return solver.solve(A,x,b,precond);
     }       
   else if (solver_name=="cg")
     {
-      SolverCG<Matrix,Vector> solver(*control,*vector_memory,
+      SolverCG<Vector> solver(*control,*vector_memory,
                                     cg_data);
       return solver.solve(A,x,b,precond);
     }
   else if (solver_name=="bicgstab")
     {
-      SolverBicgstab<Matrix,Vector> solver(*control,*vector_memory,
+      SolverBicgstab<Vector> solver(*control,*vector_memory,
                                           bicgstab_data);
       return solver.solve(A,x,b,precond);
     }
   else if (solver_name=="gmres")
     {
-      SolverGMRES<Matrix,Vector> solver(*control,*vector_memory,
+      SolverGMRES<Vector> solver(*control,*vector_memory,
                                        gmres_data);
       return solver.solve(A,x,b,precond);
     }
   else
     Assert(false,ExcSolverDoesNotExist(solver_name));
 
-  return Solver<Matrix,Vector>::breakdown;
+  return Solver<Vector>::breakdown;
 };
 
 
-template <class Matrix, class Vector>
-string SolverSelector<Matrix, Vector>::get_solver_names()
+template <class Vector>
+string SolverSelector<Vector>::get_solver_names()
 {
   return "richardson|cg|bicgstab|gmres";
 };
 
 
-template <class Matrix, class Vector>
-void SolverSelector<Matrix, Vector>::set_data(
-  const typename SolverGMRES<Matrix,Vector>::AdditionalData &data)
+template <class Vector>
+void SolverSelector<Vector>::set_data(
+  const typename SolverGMRES<Vector>::AdditionalData &data)
 { 
   gmres_data=data; 
 }
 
 
-template <class Matrix, class Vector>
-void SolverSelector<Matrix, Vector>::set_data(
-  const typename SolverRichardson<Matrix,Vector>::AdditionalData &data)
+template <class Vector>
+void SolverSelector<Vector>::set_data(
+  const typename SolverRichardson<Vector>::AdditionalData &data)
 { 
   richardson_data=data; 
 }
 
 
-template <class Matrix, class Vector>
-void SolverSelector<Matrix, Vector>::set_data(
-  const typename SolverCG<Matrix,Vector>::AdditionalData &data) 
+template <class Vector>
+void SolverSelector<Vector>::set_data(
+  const typename SolverCG<Vector>::AdditionalData &data) 
 { 
   cg_data=data; 
 }
 
 
-template <class Matrix, class Vector>
-void SolverSelector<Matrix, Vector>::set_data(
-  const typename SolverBicgstab<Matrix,Vector>::AdditionalData &data) 
+template <class Vector>
+void SolverSelector<Vector>::set_data(
+  const typename SolverBicgstab<Vector>::AdditionalData &data) 
 { 
   bicgstab_data=data; 
 };
index 9c8e3541a72fcd21cb45cc28b4ff1efc5ad60ab8..19a864ceb900b95300eeea037478800b6ada83a9 100644 (file)
@@ -155,8 +155,8 @@ for (unsigned int level = 0; level <= maxlevel; ++level)
       PreconditionMG<FDMG, Vector<TYPE> >
        precondition(multigrid, smoother, smoother, coarsegrid);
 
-//      SolverRichardson<SparseMatrix<double> , Vector<TYPE> > solver(control, mem);
-      SolverCG<SparseMatrix<double> , Vector<TYPE> > solver(control, mem);
+//      SolverRichardson<Vector<TYPE> > solver(control, mem);
+      SolverCG<Vector<TYPE> > solver(control, mem);
 
       Vector<TYPE> u(dim);
       Vector<TYPE> f(dim);

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.