]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Check in an initial working version of the preconditioned GMRES solver.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Mar 1999 10:14:50 +0000 (10:14 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Mar 1999 10:14:50 +0000 (10:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@1051 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h
new file mode 100644 (file)
index 0000000..183107e
--- /dev/null
@@ -0,0 +1,394 @@
+/*----------------------------   solver_pgmres.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __solver_pgmres_H
+#define __solver_pgmres_H
+/*----------------------------   solver_pgmres.h     ---------------------------*/
+
+#include <lac/solver.h>
+#include <lac/solver_control.h>
+#include <lac/fullmatrix.h>
+#include <vector>
+
+
+
+/**
+ * Implementation of the Restarted Preconditioned Direct Generalized
+ * Minimal Residual Method. The stopping criterion is the norm of the
+ * residual.
+ *
+ * @author Original implementation by the DEAL authors, adapted by Wolfgang Bangerth
+ */
+template<class Matrix, class Vector>
+class SolverPGMRES : public Solver<Matrix, Vector> {
+  public:
+                                    /**
+                                     * Constructor.
+                                     */
+    SolverPGMRES (SolverControl        &cn,
+                 VectorMemory<Vector> &mem,
+                 const unsigned int    n_tmp_vectors) :
+                   Solver<Matrix,Vector> (cn,mem),
+                   n_tmp_vectors (n_tmp_vectors)
+      {};
+    
+                                    /**
+                                     * Solver method.
+                                     */
+    virtual ReturnState solve (const Matrix &A,
+                              Vector       &x,
+                              const Vector &b);
+
+  protected:
+    const unsigned int n_tmp_vectors;
+    
+                                    /**
+                                     * Implementation of the computation of
+                                     * the norm of the residual.
+                                     */
+    virtual double criterion();
+    
+                                    /**
+                                     * Transformation of an upper
+                                     * Hessenberg matrix into
+                                     * tridiagonal structure by givens
+                                     * rotation of the last column
+                                     */
+    void givens_rotation (Vector &h,  Vector &b,
+                         Vector &ci, Vector &si,
+                         int col) const;
+};
+
+
+
+
+/* ------------------------- Inline functions ----------------------------- */
+
+template <class Matrix, class Vector>
+inline
+void
+SolverPGMRES<Matrix,Vector>::givens_rotation (Vector& h, Vector& b,
+                                             Vector& ci, Vector& si, 
+                                             int col) const
+{
+  for (int i=0 ; i<col ; i++)
+    {
+      const double s = si(i);
+      const double c = ci(i);
+      const double dummy = h(i);
+      h(i)   =  c*dummy + s*h(i+1);
+      h(i+1) = -s*dummy + c*h(i+1);
+    };
+  
+  const double r = 1./sqrt(h(col)*h(col) + h(col+1)*h(col+1));
+  si(col) = h(col+1) *r;
+  ci(col) = h(col)   *r;
+  h(col)  =  ci(col)*h(col) + si(col)*h(col+1);
+  b(col+1)= -si(col)*b(col);
+  b(col) *=  ci(col);
+}
+
+
+// // restarted method
+// template<class Matrix, class Vector>
+// inline int
+// SolverPGMRES<Matrix,Vector>::solve (Matrix& A, Vector& x, Vector& b)
+// {
+//   int reached, kmax = mem.n()-1;
+//   for(int j=0;;j++)
+//   {
+//     info.usediter() = j*(kmax-1);
+//     reached = dgmres(A,x,b,mem,info);
+//     if(reached) break;
+//   }
+//   if (reached<0) return 1;
+//   return 0;
+// }
+
+
+
+/*
+template<class Matrix, class Vector>
+inline
+Solver<Matrix,Vector>::ReturnState
+SolverPGMRES<Matrix,Vector>::solve (const Matrix& A,
+                                   Vector& x,
+                                   const Vector& b)
+{
+                                  // this code was written by the fathers of
+                                  // DEAL. I take absolutely no guarantees
+                                  // for any failures or airplane-explosions
+                                  // or nuclear wars or whatever resulting
+                                  // from this code. I tried to clean a bit,
+                                  // but whoever wrote this code should get
+                                  // stone, IMHO! (WB)
+
+  int kmax = n_tmp_vectors;
+  FullMatrix<double> H(kmax+1, kmax), H1(kmax+1, kmax);
+  
+  ::Vector<double> y(kmax), b0(kmax+1);
+  int i,k;
+
+  SolverControl::State conv=SolverControl::iterate;
+  
+  double rho,beta;
+
+                                  // allocate an array of n_tmp_vectors
+                                  // temporary vectors from the VectorMemory
+                                  // object
+  vector<Vector*> tmp_vectors (n_tmp_vectors, 0);
+  for (unsigned int tmp=0; tmp<n_tmp_vectors; ++tmp)
+    {
+      tmp_vectors[tmp] = memory.alloc();
+      tmp_vectors[tmp]->reinit (x.size());
+    };
+
+  
+  A.residual(*tmp_vectors[0],x,b);
+  
+  rho = tmp_vectors[0]->l2_norm();
+  beta = rho;
+  
+  tmp_vectors[0]->scale (1./rho);
+  
+  for (k=0 ; k<kmax-1 && (conv==SolverControl::iterate) ; k++)
+  {
+    A.vmult(*tmp_vectors[k+1], *tmp_vectors[k]);
+    
+    H.reinit(k+2,k+1);
+    if (k>0)  H.fill(H1);
+    
+    for (i=0 ; i<=k ; i++)
+    {
+      H(i,k) = *tmp_vectors[k+1] * *tmp_vectors[i];
+      tmp_vectors[k+1]->add(-H(i,k),*tmp_vectors[i]);
+    }
+    
+    double s = tmp_vectors[k+1]->l2_norm();
+    H(k+1,k) = s;
+    
+    // Re-orthogonalization
+    
+    //printf("\n");
+    for (i=0 ; i<=k ; i++)
+    {
+      double htmp = *tmp_vectors[k+1] * *tmp_vectors[i];
+      //printf(" %e ",htmp);
+      H(i,k) += htmp;
+      tmp_vectors[k+1]->add(-htmp,*tmp_vectors[i]);
+    }
+    //printf("\n");
+    
+    s = tmp_vectors[k+1]->l2_norm();
+    H(k+1,k) = s;
+    
+    tmp_vectors[k+1]->scale(1./s);
+    
+    // Least - Squares
+    
+    y.reinit(k+1);
+    b0.reinit(k+2);
+    b0(0) = beta;
+    H1 = H;
+    rho = H.least_squares(y,b0);
+    conv = control().check(k,rho);
+  }
+
+                                  // this will miserably fail if the
+                                  // loop above was left before k=kmax-1!
+  for (i=0 ; i<kmax ; i++)
+//  for (i=0 ; i<k ; i++)
+    x.add(y(i), *tmp_vectors[i]);
+
+
+                                  // free the allocated memory before
+                                  // leaving
+  for (unsigned int tmp=0; tmp<n_tmp_vectors; ++tmp)
+    memory.free (tmp_vectors[tmp]);
+
+
+  if (conv == SolverControl::failure)
+    return exceeded;
+  else
+    return success;
+}
+*/
+
+
+
+
+template<class Matrix, class Vector>
+inline
+Solver<Matrix,Vector>::ReturnState
+SolverPGMRES<Matrix,Vector>::solve (const Matrix& A,
+                                   Vector& x,
+                                   const Vector& b)
+{
+                                  // this code was written by the fathers of
+                                  // DEAL. I take absolutely no guarantees
+                                  // for any failures or airplane-explosions
+                                  // or nuclear wars or whatever resulting
+                                  // from this code. I tried to clean a bit,
+                                  // but whoever wrote this code should get
+                                  // stone, IMHO! (WB)
+
+  int kmax = n_tmp_vectors-1;
+                                  // allocate an array of n_tmp_vectors
+                                  // temporary vectors from the VectorMemory
+                                  // object
+  vector<Vector*> tmp_vectors (n_tmp_vectors, 0);
+  for (unsigned int tmp=0; tmp<n_tmp_vectors; ++tmp)
+    {
+      tmp_vectors[tmp] = memory.alloc();
+      tmp_vectors[tmp]->reinit (x.size());
+    };
+
+// WB  
+//  int k0   = info.usediter();
+  int k0 = 0;
+  
+
+  FullMatrix<double> H(kmax+1, kmax);
+  ::Vector<double> gamma(kmax+1), ci(kmax), si(kmax), h(kmax);
+  int i,k,reached=0,dim;
+  int left_precondition = 1;
+
+  Vector& v = *tmp_vectors[0];
+  Vector& p = *tmp_vectors[kmax];
+
+  if (left_precondition)
+  {
+    A.residual(p,x,b);
+    A.precondition(v,p);
+  }
+  else
+  {
+    A.residual(v,x,b);
+  }
+  double rho = sqrt(v*v);
+  gamma(0) = rho;
+  
+  v.equ(1./rho,v);
+
+  for (k=0 ; k<kmax-1 && (!reached) ; k++)
+  {
+    Vector& vv = *tmp_vectors[k+1];
+    
+    if (left_precondition)
+    {
+      A.vmult(p, *tmp_vectors[k]);
+      A.precondition(vv,p);
+    }
+    else
+    {
+      A.precondition(p,*tmp_vectors[k]);
+      A.vmult(vv,p);
+    }
+
+// WB why is this here?    
+//    double s0 = sqrt(vv*vv);
+    dim = k+1;
+
+    /* Orthogonalization */
+    
+    for (i=0 ; i<dim ; i++)
+    {
+      h(i) = vv * *tmp_vectors[i];
+      vv.add(-h(i),*tmp_vectors[i]);
+    }
+    double s = sqrt(vv*vv);
+    h(k+1) = s;
+    
+    /* Re-orthogonalization */
+    
+    for (i=0 ; i<dim ; i++)
+    {
+      double htmp = vv * *tmp_vectors[i];
+      h(i) += htmp;
+      vv.add(-htmp,*tmp_vectors[i]);
+    }
+    s = sqrt(vv*vv);
+    h(k+1) = s;
+    
+    vv.equ(1./s, vv);
+    
+    /*  Transformation into triagonal structure  */
+    
+    givens_rotation(h,gamma,ci,si,k);
+
+    /*  append vector on matrix  */
+
+    for (i=0 ; i<dim ; i++)
+      H(i,k) = h(i);
+
+    /*  residual  */
+
+    rho = fabs(gamma(dim));
+    
+// WB    
+//    reached = info.check_residual(k0+k,rho);
+    reached = control().check (k0+k, rho);
+  }
+  
+  /*  Calculate solution  */  
+
+  h.reinit(dim);
+  FullMatrix<double> H1(dim+1,dim);
+
+  for (i=0 ; i<dim+1 ; i++) {
+    for (int j=0 ; j<dim ; j++) {
+      H1(i,j) = H(i,j);
+    }
+  }
+
+  H1.backward(h,gamma);
+
+  if (left_precondition)
+  {
+    for (i=0 ; i<dim ; i++) {
+      x.add(h(i), *tmp_vectors[i]);
+    }
+  }
+  else
+  {
+    p = 0.;
+    for (i=0 ; i<dim ; i++)
+      p.add(h(i), *tmp_vectors[i]);
+    A.precondition(v,p);
+    x.add(1.,v);
+  }
+
+                                  // free the allocated memory before
+                                  // leaving
+  for (unsigned int tmp=0; tmp<n_tmp_vectors; ++tmp)
+    memory.free (tmp_vectors[tmp]);
+
+// WB  
+//  return reached;
+  if (reached)
+    return success;
+  else
+    return exceeded;
+};
+
+
+  
+
+
+template<class Matrix, class Vector>
+double
+SolverPGMRES<Matrix,Vector>::criterion () 
+{
+                                  // dummy implementation. this function is
+                                  // not needed for the present implementation
+                                  // of gmres
+  Assert (false, ExcInternalError());
+  return 0;
+};
+
+
+
+/*----------------------------   solver_pgmres.h     ---------------------------*/
+/* end of #ifndef __solver_pgmres_H */
+#endif
+/*----------------------------   solver_pgmres.h     ---------------------------*/

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.