]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new Bicgstab
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 8 Sep 1998 15:16:22 +0000 (15:16 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 8 Sep 1998 15:16:22 +0000 (15:16 +0000)
Iterative methods perform reinit on auxiliary vectors
PrimitiveVectorMemory does not take a size

git-svn-id: https://svn.dealii.org/trunk@590 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/solver_bicgstab.h [new file with mode: 0644]
deal.II/lac/include/lac/solver_cg.h

diff --git a/deal.II/lac/include/lac/solver_bicgstab.h b/deal.II/lac/include/lac/solver_bicgstab.h
new file mode 100644 (file)
index 0000000..5441fd5
--- /dev/null
@@ -0,0 +1,183 @@
+/*----------------------------   solver_bicgstab.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __solver_bicgstab_H
+#define __solver_bicgstab_H
+/*----------------------------   solver_bicgstab.h     ---------------------------*/
+
+#include <base/logstream.h>
+#include <lac/solver.h>
+#include <lac/solver_control.h>
+
+/**
+ * Bicgstab algorithm by van der Vorst.
+ */
+template<class Matrix, class Vector>
+class SolverBicgstab : public Solver<Matrix,Vector>
+{
+public:
+  SolverBicgstab(SolverControl &cn, VectorMemory<Vector> &mem) :
+                   Solver<Matrix,Vector>(cn,mem) {};
+
+                                  /**
+                                   * Solve primal problem only.
+                                   */
+  virtual ReturnState solve (const Matrix &A,
+                            Vector       &x,
+                            const Vector &b);
+
+  protected:
+                                    /**
+                                     * Computation of the stopping criterion.
+                                     */
+    virtual double criterion();
+
+                                  /**
+                                   * Auxiliary vectors.
+                                   */
+  Vector *Vx, *Vr, *Vrbar, *Vp, *Vy, *Vz, *Vs, *Vt, *Vv;
+  const Vector *Vb;
+  
+                                  /**
+                                   * Pointer to the system matrix.
+                                   */
+  const Matrix *MA;
+  
+                                  /**
+                                   * Auxiliary values.
+                                   */
+  double alpha, beta, omega, rho, rhobar;
+  
+                                  /**
+                                   * Current iteration step.
+                                   */
+  unsigned step;
+  
+                                  /**
+                                   * Residual.
+                                   */
+  double res;
+  
+private:
+                                  /**
+                                   * Everything before the iteration loop.
+                                   */
+  SolverControl::State start();
+
+                                  /**
+                                   * The iteration loop itself.
+                                   */
+  ReturnState iterate();
+  
+};
+
+template<class Matrix, class Vector> double
+SolverBicgstab<Matrix, Vector>::criterion()
+{
+  res = MA->residual(*Vt, *Vx, *Vb);
+  return res;
+}
+
+
+template < class Matrix, class Vector > SolverControl::State
+SolverBicgstab<Matrix, Vector>::start()
+{
+  res = MA->residual(*Vr, *Vx, *Vb);
+  Vp->reinit(*Vx);
+  Vv->reinit(*Vx);
+  Vrbar->equ(1.,*Vr);
+  SolverControl::State state = control().check(step, res);
+  return state;
+}
+
+template<class Matrix, class Vector> Solver<Matrix,Vector>::ReturnState
+SolverBicgstab<Matrix, Vector>::iterate()
+{
+  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;
+  
+  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);
+    
+    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)
+{
+  deallog.push("Bicgstab");
+  Vr    = memory.alloc(); Vr->reinit(x);
+  Vrbar = memory.alloc(); Vrbar->reinit(x);
+  Vp    = memory.alloc();
+  Vy    = memory.alloc(); Vy->reinit(x);
+  Vz    = memory.alloc(); Vz->reinit(x);
+  Vs    = memory.alloc(); Vs->reinit(x);
+  Vt    = memory.alloc(); Vt->reinit(x);
+  Vv    = memory.alloc();
+
+  MA = &A;
+  Vx = &x;
+  Vb = &b;
+
+  step = 0;
+
+  ReturnState state;
+  
+  do 
+  {
+    deallog << "Go!" << endl;
+    if (start() == SolverControl::success) break;  
+    state = iterate();
+  }
+  while (state == breakdown);
+
+  memory.free(Vr);
+  memory.free(Vrbar);
+  memory.free(Vp);
+  memory.free(Vy);
+  memory.free(Vz);
+  memory.free(Vs);
+  memory.free(Vt);
+  memory.free(Vv);
+  
+  deallog.pop();
+  return state;
+}
+
+
+/*----------------------------   solver_bicgstab.h     ---------------------------*/
+/* end of #ifndef __solver_bicgstab_H */
+#endif
+/*----------------------------   solver_bicgstab.h     ---------------------------*/
index 548b5489f4a8d72be4281f3694fa2bc6425bdb82..e7e8cf6b60669b242116f1d24bde6249e23eb2b2 100644 (file)
@@ -83,10 +83,10 @@ SolverPCG<Matrix,Vector>::solve (const Matrix &A,
   SolverControl::State conv=SolverControl::iterate;
   
                                   // Memory allocation
-  Vr  = memory.alloc(); Vector& g  = *Vr;
-  Vp  = memory.alloc(); Vector& h  = *Vp;
-  Vz  = memory.alloc(); Vector& d  = *Vz;
-  VAp = memory.alloc(); Vector& Ad = *VAp; 
+  Vr  = memory.alloc(); Vector& g  = *Vr; g.reinit(x);
+  Vp  = memory.alloc(); Vector& h  = *Vp; h.reinit(x);
+  Vz  = memory.alloc(); Vector& d  = *Vz; d.reinit(x);
+  VAp = memory.alloc(); Vector& Ad = *VAp; Ad.reinit(x);
 
                                   // Implementation taken from the DEAL
                                   // library

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.