]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
First cleaned-up version.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Mar 1999 10:50:44 +0000 (10:50 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Mar 1999 10:50:44 +0000 (10:50 +0000)
git-svn-id: https://svn.dealii.org/trunk@1055 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/solver_gmres.h

index 183107e714e6eef17d1602b360565aa93e7edcf7..8ed068e45181dc2567b23164fb35ca989ff34502 100644 (file)
@@ -106,115 +106,6 @@ SolverPGMRES<Matrix,Vector>::givens_rotation (Vector& h, Vector& b,
 
 
 
-/*
-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
@@ -228,10 +119,10 @@ SolverPGMRES<Matrix,Vector>::solve (const Matrix& A,
                                   // 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)
+                                  // but whoever wrote this code in the first
+                                  // place should get stoned, IMHO! (WB)
 
-  int kmax = n_tmp_vectors-1;
+  const unsigned int kmax = n_tmp_vectors-1;
                                   // allocate an array of n_tmp_vectors
                                   // temporary vectors from the VectorMemory
                                   // object
@@ -246,71 +137,82 @@ SolverPGMRES<Matrix,Vector>::solve (const Matrix& A,
 //  int k0   = info.usediter();
   int k0 = 0;
   
-
+                                  // matrix used for the orthogonalization
+                                  // process later
   FullMatrix<double> H(kmax+1, kmax);
+
+                                  // some additional vectors, also used
+                                  // in the orthogonalization
   ::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];
+
+  unsigned int dim;
+
+  SolverControl::State reached = SolverControl::iterate;
+  
+                                  // switch to determine whether we want a
+                                  // left or a right preconditioner. at
+                                  // present, left is default, but both
+                                  // ways are implemented
+  const bool left_precondition = true;
+
+                                  // define two aliases
+  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);
-  }
+    {
+      A.residual(p,x,b);
+      A.precondition(v,p);
+    } else {
+      A.residual(v,x,b);
+    };
  
-  double rho = sqrt(v*v);
+  double rho = v.l2_norm();
   gamma(0) = rho;
-  
-  v.equ(1./rho,v);
 
-  for (k=0 ; k<kmax-1 && (!reached) ; k++)
+  v.scale (1./rho);
+
+  unsigned int k;
+  for (k=0 ; k<kmax-1 && (reached==SolverControl::iterate) ; 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);
-    }
+      {
+       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);
+    for (unsigned int i=0 ; i<dim ; ++i)
+      {
+       h(i) = vv * *tmp_vectors[i];
+       vv.add(-h(i),*tmp_vectors[i]);
+      };
+    
+    double s = vv.l2_norm();
     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);
+    for (unsigned i=0; i<dim; ++i)
+      {
+       double htmp = vv * *tmp_vectors[i];
+       h(i) += htmp;
+       vv.add(-htmp,*tmp_vectors[i]);
+      };
+    
+    s = vv.l2_norm();
     h(k+1) = s;
     
-    vv.equ(1./s, vv);
+    vv.scale(1./s);
     
     /*  Transformation into triagonal structure  */
     
@@ -318,15 +220,13 @@ SolverPGMRES<Matrix,Vector>::solve (const Matrix& A,
 
     /*  append vector on matrix  */
 
-    for (i=0 ; i<dim ; i++)
+    for (unsigned int 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);
   }
   
@@ -335,28 +235,23 @@ SolverPGMRES<Matrix,Vector>::solve (const Matrix& A,
   h.reinit(dim);
   FullMatrix<double> H1(dim+1,dim);
 
-  for (i=0 ; i<dim+1 ; i++) {
-    for (int j=0 ; j<dim ; j++) {
+  for (unsigned int i=0; i<dim+1; ++i)
+    for (unsigned 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++) {
+    for (unsigned int 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);
-  }
+    {
+      p = 0.;
+      for (unsigned int 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

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.