]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add minres
authorThomas Richter <thomas.richter@iwr.uni-heidelberg.de>
Wed, 8 Mar 2000 14:03:49 +0000 (14:03 +0000)
committerThomas Richter <thomas.richter@iwr.uni-heidelberg.de>
Wed, 8 Mar 2000 14:03:49 +0000 (14:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@2567 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/lac/include/lac/solver_minres.h b/deal.II/lac/include/lac/solver_minres.h
new file mode 100644 (file)
index 0000000..2e50c30
--- /dev/null
@@ -0,0 +1,310 @@
+//----------------------------  solver_minres.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  solver_minres.h  -----------------------
+#ifndef __deal2__solver_minres_h
+#define __deal2__solver_minres_h
+
+
+#include <lac/solver.h>
+#include <lac/solver_control.h>
+#include <base/logstream.h>
+#include <cmath>
+
+
+
+/**
+ * Preconditioned MinRes method.
+ *
+ * Like all other solver classes, this class has a local structure called
+ * #AdditionalData# which is used to pass additional parameters to the
+ * solver, like damping parameters or the number of temporary vectors. We
+ * use this additional structure instead of passing these values directly
+ * to the constructor because this makes the use of the #SolverSelector# and
+ * other classes much easier and guarantees that these will continue to
+ * work even if number or type of the additional parameters for a certain
+ * solver changes.
+ *
+ * However, since the MinRes method does not need additional data, the respective
+ * structure is empty and does not offer any functionality. The constructor
+ * has a default argument, so you may call it without the additional
+ * parameter.
+ *
+ * The preconditioner has to be positive definite and symmetric
+ *
+ * @author Thomas Richter, 2000
+ */
+template <class Matrix = SparseMatrix<double>, class Vector = Vector<double> >
+class SolverMinRes : public Solver<Matrix,Vector>
+{
+  public:
+                                    /**
+                                     * Standardized data struct to
+                                     * pipe additional data to the
+                                     * solver. This solver does not
+                                     * need additional data yet.
+                                     */
+    struct AdditionalData {};
+
+                                    /**
+                                     * Constructor.
+                                     */
+    SolverMinRes (SolverControl &cn,
+                 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,
+          const Preconditioner& precondition);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcPreconditionerNotDefinite);
+    
+
+  protected:
+                                    /**
+                                     * Implementation of the computation of
+                                     * the norm of the residual.
+                                     */
+    virtual long double criterion();
+    
+                                    /**
+                                     * Temporary vectors, allocated through
+                                     * the #VectorMemory# object at the start
+                                     * of the actual solution process and
+                                     * deallocated at the end.
+                                     */
+    Vector *Vu0, *Vu1, *Vu2;
+    Vector *Vm0, *Vm1, *Vm2;   
+    Vector *Vv;
+    
+                                    /**
+                                     * Within the iteration loop, the
+                                     * square of the residual vector is
+                                     * stored in this variable. The
+                                     * function #criterion# uses this
+                                     * variable to compute the convergence
+                                     * value, which in this class is the
+                                     * norm of the residual vector and thus
+                                     * the square root of the #res2# value.
+                                     */
+    long double res2;
+};
+
+
+/*------------------------- Implementation ----------------------------*/
+
+
+template<class Matrix, class Vector>
+SolverMinRes<Matrix,Vector>::SolverMinRes(SolverControl &cn,
+                                         VectorMemory<Vector> &mem,
+                                         const AdditionalData &) :
+               Solver<Matrix,Vector>(cn,mem) {};
+
+
+template<class Matrix, class Vector>
+long double
+SolverMinRes<Matrix,Vector>::criterion()
+{
+  return res2;
+};
+
+
+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)
+{
+  SolverControl::State conv=SolverControl::iterate;
+
+  deallog.push("minres");
+
+
+  unsigned int VS = b.size();
+
+  
+                                  // Memory allocation
+  Vu0  = memory.alloc();
+  Vu1  = memory.alloc();
+  Vu2  = memory.alloc();
+  Vv   = memory.alloc();
+  Vm0  = memory.alloc();
+  Vm1  = memory.alloc();
+  Vm2  = memory.alloc();
+                                  // define some aliases for simpler access
+  typedef Vector vecref;
+  vecref u[3] = {*Vu0, *Vu1, *Vu2};
+  vecref m[3] = {*Vm0, *Vm1, *Vm2};
+  vecref v    = *Vv;
+                                  // resize the vectors, but do not set
+                                  // the values since they'd be overwritten
+                                  // soon anyway.
+  u[0].reinit(VS,true);
+  u[1].reinit(VS,true);
+  u[2].reinit(VS,true);
+  m[0].reinit(VS,true);
+  m[1].reinit(VS,true);
+  m[2].reinit(VS,true);
+  v.reinit(VS,true);
+
+                                  // some values needed
+  vector<double> delta(3);
+  vector<double> f(2);
+  vector<double> e(2); 
+
+  double r_l2 = 0;
+  double r0   = 0;
+  double tau = 0;
+  double c    = 0;
+  double gamma = 0;
+  double s = 0;
+  double d_ = 0;
+  double d = 0;  
+
+                                  // The iteration step.
+  int j = 1;
+  
+
+                                  // Start of the solving process
+  
+                                  // The algorithm is taken from
+                                  // Astrid Battermann, Master thesis
+                                  // with some changes
+
+  A.vmult(m[0],x);
+  u[1] = b;
+  u[1].add(-1.,m[0]);
+                                  // Precondition is applied.
+                                  // The preconditioner has to be
+                                  // positiv definite and symmetric
+
+                                  // M v = u[1]
+  precondition (v,u[1]);
+  
+  delta[1] = v * u[1];
+  Assert (delta[1]>=0, ExcPreconditionerNotDefinite());
+  
+  r0 = sqrt(delta[1]);
+  r_l2 = r0;
+  
+  
+  u[0].reinit(VS,0);
+  delta[0] = 1.;
+  m[0].reinit(VS,0);
+  m[1].reinit(VS,0);
+  m[2].reinit(VS,0);
+                                  
+  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);
+
+      A.vmult(u[2],v);
+      u[2].add(-sqrt(delta[1]/delta[0]),u[0]);
+
+      gamma = u[2] * v;
+      u[2].add (-gamma / sqrt(delta[1]), u[1]);
+      m[0] = v;
+      
+                                      // precondition: solve M v = u[2]
+                                      // Preconditioner has to be positiv
+                                      // definite and symmetric.
+      precondition(v,u[2]);
+      delta[2] = v * u[2];
+
+      Assert (delta[2]>=0, ExcPreconditionerNotDefinite());
+
+      if (j==1)
+       {
+         d_ = gamma;
+         e[1] = sqrt(delta[2]);
+       }
+      if (j>1)
+       {
+         d_ = s * e[0] - c * gamma;
+         e[0] = c * e[0] + s * gamma;
+         f[1] = s * sqrt(delta[2]);
+         e[1] = -c * sqrt(delta[2]);
+       }
+
+      d = sqrt (d_*d_ + delta[2]);
+      
+      if (j>1) tau *= s / c;
+      c = d_ / d;
+      tau *= c;
+      
+      s = sqrt(delta[2]) / d;
+
+      if (j==1)
+       tau = r0 * c;
+
+      m[0].add (-e[0], m[1]);
+      if (j>1)
+       m[0].add (-f[0],m[2]);
+      m[0].scale(1./d);
+      x.add (tau, m[0]);
+      r_l2 *= fabs(s);
+
+      conv = control().check(j,r_l2);
+      
+                                      // next iteration step
+      ++j;
+                                      // All vectors have to be shifted
+                                      // one iteration step.
+                                      // This should be changed one time.
+      m[2] = m[1];
+      m[1] = m[0];
+
+      u[0] = u[1];
+      u[1] = u[2];
+      f[0] = f[1];
+      e[0] = e[1];
+      delta[0] = delta[1];
+      delta[1] = delta[2];
+    }
+
+                                  // Deallocation of Memory
+  memory.free(Vu0);
+  memory.free(Vu1);
+  memory.free(Vu2);
+  memory.free(Vv); 
+  memory.free(Vm0);
+  memory.free(Vm1);
+  memory.free(Vm2);
+                                  // Output
+  deallog.pop ();
+  
+  if (conv == SolverControl::failure)
+    return exceeded;
+  
+  return success;
+};
+
+
+#endif
+

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.