MFEM v2.0
gmres.cpp
Go to the documentation of this file.
```00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.googlecode.com.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011
00012 #include <iostream>
00013 #include <iomanip>
00014 #include <math.h>
00015 #include "vector.hpp"
00016 #include "matrix.hpp"
00017 #include "densemat.hpp"
00018
00019 //*****************************************************************
00020 // Iterative template routine -- GMRES
00021 //
00022 // GMRES solves the unsymmetric linear system Ax = b using the
00023 // Generalized Minimum Residual method
00024 //
00025 // GMRES follows the algorithm described on p. 20 of the
00026 // SIAM Templates book.
00027 //
00028 // The return value indicates convergence within max_iter (input)
00029 // iterations (0), or no convergence within max_iter iterations (1).
00030 //
00031 // Upon successful return, output arguments have the following values:
00032 //
00033 //        x  --  initial approximation / approximate solution to Ax = b
00034 // max_iter  --  the number of iterations performed before the
00035 //               tolerance was reached
00036 //      tol  --  the residual after the final iteration
00037 //
00038 //*****************************************************************
00039
00040 inline void GeneratePlaneRotation (double &dx, double &dy, double &cs, double &sn)
00041 {
00042    if (dy == 0.0) {
00043       cs = 1.0;
00044       sn = 0.0;
00045    } else if (fabs(dy) > fabs(dx)) {
00046       double temp = dx / dy;
00047       sn = 1.0 / sqrt( 1.0 + temp*temp );
00048       cs = temp * sn;
00049    } else {
00050       double temp = dy / dx;
00051       cs = 1.0 / sqrt( 1.0 + temp*temp );
00052       sn = temp * cs;
00053    }
00054 }
00055
00056 inline void ApplyPlaneRotation (double &dx, double &dy, double &cs, double &sn)
00057 {
00058    double temp  =  cs * dx + sn * dy;
00059    dy = -sn * dx + cs * dy;
00060    dx = temp;
00061 }
00062
00063 inline void Update (Vector &x, int k, DenseMatrix &h, Vector &s, Array<Vector*> &v)
00064 {
00065    Vector y(s.Size());  y = s;
00066
00067    // Backsolve:
00068    for (int i = k; i >= 0; i--) {
00069       y(i) /= h(i,i);
00070       for (int j = i - 1; j >= 0; j--)
00071          y(j) -= h(j,i) * y(i);
00072    }
00073
00074    for (int j = 0; j <= k; j++)
00075       x.Add(y(j),*v[j]);
00076 }
00077
00078 inline double norm (Vector & u)
00079 {
00080    return sqrt(u*u);
00081 }
00082
00083 int
00084 GMRES(const Operator &A, Vector &x, const Vector &b,
00085       const Operator &M, int &max_iter,
00086       int m, double &tol, double &atol, int printit)
00087 {
00088    int n = A.Size();
00089
00090    DenseMatrix H(m+1,m);
00091    Vector s(m+1), cs(m+1), sn(m+1);
00092    Vector w(n), av(n);
00093
00094    double resid;
00095    int i, j, k;
00096
00097    M.Mult(b,w);
00098    double normb = norm(w); // normb = ||M b||
00099    if (normb == 0.0)
00100       normb = 1;
00101
00102    Vector r(n);
00103    A.Mult(x, r);
00104    subtract(b,r,w);
00105    M.Mult(w, r);           // r = M (b - A x)
00106    double beta = norm(r);  // beta = ||r||
00107
00108    resid = beta / normb;
00109
00110    if (resid * resid <= tol) {
00111       tol = resid * resid;
00112       max_iter = 0;
00113       return 0;
00114    }
00115
00116    if (printit)
00117       cout << "   Pass : " << setw(2) << 1
00118            << "   Iteration : " << setw(3) << 0
00119            << "  (r, r) = " << beta*beta << endl;
00120
00121    tol *= (normb*normb);
00122    tol = (atol > tol) ? atol : tol;
00123
00124    Array<Vector *> v(m+1);
00125    for (i= 0; i<=m; i++)
00126       v[i] = NULL;
00127
00128    j = 1;
00129    while (j <= max_iter)
00130    {
00131       if (v[0] == NULL)  v[0] = new Vector(n);
00132       (*v[0]) = 0.0;
00133       v[0] -> Add (1.0/beta, r);   // v[0] = r / ||r||
00134       s = 0.0; s(0) = beta;
00135
00136       for (i = 0; i < m && j <= max_iter; i++) {
00137          A.Mult((*v[i]),av);
00138          M.Mult(av,w);              // w = M A v[i]
00139
00140          for (k = 0; k <= i; k++) {
00141             H(k,i) = w * (*v[k]);    // H(k,i) = w * v[k]
00142             w.Add(-H(k,i), (*v[k])); // w -= H(k,i) * v[k]
00143          }
00144
00145          H(i+1,i)  = norm(w);       // H(i+1,i) = ||w||
00146          if (v[i+1] == NULL)  v[i+1] = new Vector(n);
00147          (*v[i+1]) = 0.0;
00148          v[i+1] -> Add (1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
00149
00150          for (k = 0; k < i; k++)
00151             ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
00152
00153          GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
00154          ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
00155          ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
00156
00157          resid = fabs(s(i+1));
00158          if (printit)
00159             cout << "   Pass : " << setw(2) << j
00160                  << "   Iteration : " << setw(3) << i+1
00161                  << "  (r, r) = " << resid*resid << endl;
00162
00163          if ( resid*resid < tol) {
00164             Update(x, i, H, s, v);
00165             tol = resid * resid;
00166             max_iter = j;
00167             for (i= 0; i<=m; i++)
00168                if (v[i])  delete v[i];
00169             return 0;
00170          }
00171       }
00172
00173       if (printit)
00174          cout << "Restarting..." << endl;
00175
00176       Update(x, i-1, H, s, v);
00177
00178       A.Mult(x, r);
00179       subtract(b,r,w);
00180       M.Mult(w, r);           // r = M (b - A x)
00181       beta = norm(r);         // beta = ||r||
00182       if ( resid*resid < tol) {
00183          tol = resid * resid;
00184          max_iter = j;
00185          for (i= 0; i<=m; i++)
00186             if (v[i])  delete v[i];
00187          return 0;
00188       }
00189
00190       j++;
00191    }
00192
00193    tol = resid * resid;
00194    for (i= 0; i<=m; i++)
00195       if (v[i])  delete v[i];
00196    return 1;
00197 }
00198
00199 int
00200 aGMRES(const Operator &A, Vector &x, const Vector &b,
00201        const Operator &M, int &max_iter,
00202        int m_max, int m_min, int m_step, double cf,
00203        double &tol, double &atol, int printit)
00204 {
00205    int n = A.Size();
00206
00207    int m = m_max;
00208
00209    DenseMatrix H(m+1,m);
00210    Vector s(m+1), cs(m+1), sn(m+1);
00211    Vector w(n), av(n);
00212
00213    double r1, resid;
00214    int i, j, k;
00215
00216    M.Mult(b,w);
00217    double normb = norm(w); // normb = ||M b||
00218    if (normb == 0.0)
00219       normb = 1;
00220
00221    Vector r(n);
00222    A.Mult(x, r);
00223    subtract(b,r,w);
00224    M.Mult(w, r);           // r = M (b - A x)
00225    double beta = norm(r);  // beta = ||r||
00226
00227    resid = beta / normb;
00228
00229    if (resid * resid <= tol) {
00230       tol = resid * resid;
00231       max_iter = 0;
00232       return 0;
00233    }
00234
00235    if (printit)
00236       cout << "   Pass : " << setw(2) << 1
00237            << "   Iteration : " << setw(3) << 0
00238            << "  (r, r) = " << beta*beta << endl;
00239
00240    tol *= (normb*normb);
00241    tol = (atol > tol) ? atol : tol;
00242
00243    m = m_max;
00244    Array<Vector *> v(m+1);
00245    for (i= 0; i<=m; i++) {
00246       v[i] = new Vector(n);
00247       (*v[i]) = 0.0;
00248    }
00249
00250    j = 1;
00251    while (j <= max_iter) {
00252       (*v[0]) = 0.0;
00253       v[0] -> Add (1.0/beta, r);   // v[0] = r / ||r||
00254       s = 0.0; s(0) = beta;
00255
00256       r1 = beta;
00257
00258       for (i = 0; i < m && j <= max_iter; i++) {
00259          A.Mult((*v[i]),av);
00260          M.Mult(av,w);              // w = M A v[i]
00261
00262          for (k = 0; k <= i; k++) {
00263             H(k,i) = w * (*v[k]);    // H(k,i) = w * v[k]
00264             w.Add(-H(k,i), (*v[k])); // w -= H(k,i) * v[k]
00265          }
00266
00267          H(i+1,i)  = norm(w);       // H(i+1,i) = ||w||
00268          (*v[i+1]) = 0.0;
00269          v[i+1] -> Add (1.0/H(i+1,i), w); // v[i+1] = w / H(i+1,i)
00270
00271          for (k = 0; k < i; k++)
00272             ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
00273
00274          GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
00275          ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
00276          ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
00277
00278          resid = fabs(s(i+1));
00279          if (printit)
00280             cout << "   Pass : " << setw(2) << j
00281                  << "   Iteration : " << setw(3) << i+1
00282                  << "  (r, r) = " << resid*resid << endl;
00283
00284          if ( resid*resid < tol) {
00285             Update(x, i, H, s, v);
00286             tol = resid * resid;
00287             max_iter = j;
00288             for (i= 0; i<=m; i++)
00289                delete v[i];
00290             return 0;
00291          }
00292       }
00293
00294       if (printit)
00295          cout << "Restarting..." << endl;
00296
00297       Update(x, i-1, H, s, v);
00298
00299       A.Mult(x, r);
00300       subtract(b,r,w);
00301       M.Mult(w, r);           // r = M (b - A x)
00302       beta = norm(r);         // beta = ||r||
00303       if ( resid*resid < tol) {
00304          tol = resid * resid;
00305          max_iter = j;
00306          for (i= 0; i<=m; i++)
00307             delete v[i];
00308          return 0;
00309       }
00310
00311       if (beta/r1 > cf)
00312       {
00313          if (m - m_step >= m_min)
00314             m -= m_step;
00315          else
00316             m = m_max;
00317       }
00318
00319       j++;
00320    }
00321
00322    tol = resid * resid;
00323    for (i= 0; i<=m; i++)
00324       delete v[i];
00325    return 1;
00326 }
```