]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
minor changes
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 12 Mar 1999 09:28:02 +0000 (09:28 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 12 Mar 1999 09:28:02 +0000 (09:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@999 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/ivector.h
deal.II/lac/include/lac/solver_bicgstab.h
deal.II/lac/source/solver_control.cc

index ce6ac8aa946ca320d565eb4329951587844ba2a3..cc2d22ee48bba0beba28dad5444d2f9291291c06 100644 (file)
@@ -158,12 +158,7 @@ protected:
 };
 
 
-
-
-
-
-
-
+/*-------------------------Inline functions -------------------------------*/
 
 
 inline unsigned int iVector::n() const
index 5441fd59c4457353bd78cb33b51e21d101ef83d3..f8b7724e34e36fa3342a8263f15f75bcc05c7cd6 100644 (file)
  * Bicgstab algorithm by van der Vorst.
  */
 template<class Matrix, class Vector>
-class SolverBicgstab : public Solver<Matrix,Vector>
+class SolverBicgstab //<Matrix,Vector>
+  : public Solver<Matrix,Vector>
 {
-public:
-  SolverBicgstab(SolverControl &cn, VectorMemory<Vector> &mem) :
-                   Solver<Matrix,Vector>(cn,mem) {};
+  public:
+                                    /**
+                                     * Constructor.
+                                     */
+    SolverBicgstab(SolverControl &cn, VectorMemory<Vector> &mem);
 
-                                  /**
-                                   * Solve primal problem only.
-                                   */
-  virtual ReturnState solve (const Matrix &A,
-                            Vector       &x,
-                            const Vector &b);
+                                    /**
+                                     * Solve primal problem only.
+                                     */
+    virtual ReturnState solve (const Matrix &A,
+                              Vector       &x,
+                              const Vector &b);
 
   protected:
                                     /**
@@ -31,45 +34,106 @@ public:
                                      */
     virtual double criterion();
 
-                                  /**
-                                   * Auxiliary vectors.
-                                   */
-  Vector *Vx, *Vr, *Vrbar, *Vp, *Vy, *Vz, *Vs, *Vt, *Vv;
-  const Vector *Vb;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    Vector *Vx;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    Vector *Vr;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    Vector *Vrbar;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    Vector *Vp;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    Vector *Vy;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    Vector *Vz;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    Vector *Vs;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    Vector *Vt;
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    Vector *Vv;
+                                    /**
+                                     * Right hand side vector.
+                                     */
+    const Vector *Vb;
   
-                                  /**
-                                   * Pointer to the system matrix.
-                                   */
-  const Matrix *MA;
+                                    /**
+                                     * Pointer to the system matrix.
+                                     */
+    const Matrix *MA;
   
-                                  /**
-                                   * Auxiliary values.
-                                   */
-  double alpha, beta, omega, rho, rhobar;
+                                    /**
+                                     * Auxiliary value.
+                                     */
+    double alpha;
+                                    /**
+                                     * Auxiliary value.
+                                     */
+    double beta;
+                                    /**
+                                     * Auxiliary value.
+                                     */
+    double omega;
+                                    /**
+                                     * Auxiliary value.
+                                     */
+    double rho;
+                                    /**
+                                     * Auxiliary value.
+                                     */
+    double rhobar;
   
-                                  /**
-                                   * Current iteration step.
-                                   */
-  unsigned step;
+                                    /**
+                                     * Current iteration step.
+                                     */
+    unsigned int step;
   
-                                  /**
-                                   * Residual.
-                                   */
-  double res;
+                                    /**
+                                     * Residual.
+                                     */
+    double res;
   
-private:
-                                  /**
-                                   * Everything before the iteration loop.
-                                   */
-  SolverControl::State start();
-
-                                  /**
-                                   * The iteration loop itself.
-                                   */
-  ReturnState iterate();
+  private:
+                                    /**
+                                     * Everything before the iteration loop.
+                                     */
+    SolverControl::State start();
+
+                                    /**
+                                     * The iteration loop itself.
+                                     */
+    ReturnState iterate();
   
 };
 
+/*-------------------------Inline functions -------------------------------*/
+
+template<class Matrix, class Vector>
+SolverBicgstab<Matrix, Vector>::SolverBicgstab(SolverControl &cn,
+                                              VectorMemory<Vector> &mem)
+               :
+               Solver<Matrix,Vector>(cn,mem)
+{}
+
+
 template<class Matrix, class Vector> double
 SolverBicgstab<Matrix, Vector>::criterion()
 {
@@ -105,37 +169,38 @@ SolverBicgstab<Matrix, Vector>::iterate()
   Vector& v = *Vv;
   
   do
-  {
-    rhobar = r*rbar;
-    beta   = rhobar * alpha / (rho * omega);
-    rho    = rhobar;
-    p.sadd(beta, 1., r, -beta*omega, v);
-    MA->precondition(y,p);
-    MA->vmult(v,y);
-    rhobar = rbar * v;
-
-    if (fabs(rhobar) < 1.e-19) return ReturnState(breakdown);
+    {
+      rhobar = r*rbar;
+      beta   = rhobar * alpha / (rho * omega);
+      rho    = rhobar;
+      p.sadd(beta, 1., r, -beta*omega, v);
+      MA->precondition(y,p);
+      MA->vmult(v,y);
+      rhobar = rbar * v;
+
+      if (fabs(rhobar) < 1.e-19) return ReturnState(breakdown);
     
-    alpha = rho/rhobar;
-    s.equ(1., r, -alpha, v);
-    MA->precondition(z,s);
-    MA->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());
-  }
+      alpha = rho/rhobar;
+      s.equ(1., r, -alpha, v);
+      MA->precondition(z,s);
+      MA->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());
+    }
   while (state == SolverControl::iterate);
   if (state == SolverControl::success) return success;
   return exceeded;
 }
 
+
 template<class Matrix, class Vector> Solver<Matrix,Vector>::ReturnState
 SolverBicgstab<Matrix, Vector>::solve(const Matrix &A,
-                     Vector       &x,
-                     const Vector &b)
+                                     Vector       &x,
+                                     const Vector &b)
 {
   deallog.push("Bicgstab");
   Vr    = memory.alloc(); Vr->reinit(x);
@@ -156,11 +221,11 @@ SolverBicgstab<Matrix, Vector>::solve(const Matrix &A,
   ReturnState state;
   
   do 
-  {
-    deallog << "Go!" << endl;
-    if (start() == SolverControl::success) break;  
-    state = iterate();
-  }
+    {
+      deallog << "Go!" << endl;
+      if (start() == SolverControl::success) break;  
+      state = iterate();
+    }
   while (state == breakdown);
 
   memory.free(Vr);
index edd8592fe7d7db62310943dd6ece2a99e56cf410..a38ebc5ce72d04253096c87c89672a8586fa7905 100644 (file)
@@ -29,8 +29,18 @@ SolverControl::check (const unsigned int step,
   
   lstep  = step;
   lvalue = check_value;
-  if (step>=maxsteps) return failure;
-  if (check_value <= tol) return success;
+  if (step>=maxsteps)
+    {
+      deallog << "Failure step " << step
+             << " value " << check_value << endl;
+      return failure;
+    }
+  if (check_value <= tol)
+    {
+      deallog << "Convergence step " << step
+             << " value " << check_value << endl;      
+      return success;
+    }
   return 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.