OpenMEEG
gmres.h
Go to the documentation of this file.
1 // Project Name: OpenMEEG (http://openmeeg.github.io)
2 // © INRIA and ENPC under the French open source license CeCILL-B.
3 // See full copyright notice in the file LICENSE.txt
4 // If you make a copy of this file, you must either:
5 // - provide also LICENSE.txt and modify this header to refer to it.
6 // - replace this header by the LICENSE.txt content.
7 
8 #pragma once
9 
10 #include "vector.h"
11 #include "matrix.h"
12 #include "sparse_matrix.h"
13 
14 #include <OpenMEEG_Export.h>
15 
16 namespace OpenMEEG {
17 
18  // ===================================
19  // = Define a Jacobi preconditionner =
20  // ===================================
21  template <typename M>
22  class Jacobi {
23  public:
24  Jacobi (const M& m): J(m.nlin()) {
25  for ( unsigned i = 0; i < m.nlin(); ++i) {
26  J(i, i) = 1.0 / m(i,i);
27  }
28  }
29 
30  Vector operator()(const Vector& g) const {
31  return J*g;
32  }
33 
34  ~Jacobi () {};
35  private:
36  SparseMatrix J; // diagonal
37  };
38 
39  // =========================
40  // = Define a GMRes solver =
41  // =========================
42 
43  void GeneratePlaneRotation(double &dx, double &dy, double &cs, double &sn)
44  {
45  if (dy == 0.0) {
46  cs = 1.0;
47  sn = 0.0;
48  } else if (std::abs(dy) > std::abs(dx)) {
49  double temp = dx / dy;
50  sn = 1.0 / sqrt( 1.0 + temp*temp );
51  cs = temp * sn;
52  } else {
53  double temp = dy / dx;
54  cs = 1.0 / sqrt( 1.0 + temp*temp );
55  sn = temp * cs;
56  }
57  }
58 
59  void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
60  {
61  double temp = cs * dx + sn * dy;
62  dy = -sn * dx + cs * dy;
63  dx = temp;
64  }
65 
66  template<class T>
67  void Update(Vector &x, int k, T &h, Vector &s, Vector v[])
68  {
69  Vector y(s);
70  // Backsolve:
71  for (int i = k; i >= 0; i--) {
72  y(i) /= h(i,i);
73  for (int j = i - 1; j >= 0; j--)
74  y(j) -= h(j,i) * y(i);
75  }
76  for (int j = 0; j <= k; j++) {
77  x += v[j] * y(j);
78  }
79  }
80 
81  // code taken from http://www.netlib.org/templates/cpp/gmres.h and modified
82  template<class T,class P> // T should be a linear operator, and P a preconditionner
83  unsigned GMRes(const T& A, const P& M, Vector &x, const Vector& b, int max_iter, double tol,unsigned m) {
84 
85  // m is the size of the Krylov subspace, if m<A.nlin(), it is a restarted GMRes (for saving memory)
86  Matrix H(m+1,m);
87  x.set(0.0);
88 
89  double resid;
90  int i, j = 1, k;
91  Vector s(m+1), cs(m+1), sn(m+1), w;
92 
93  double normb = (M(b)).norm();//(M*b).norm()
94  Vector r = M(b-A*x);//M.solve(b - A * x);
95  double beta = r.norm();
96 
97  if (normb == 0.0)
98  normb = 1;
99 
100  if ((resid = r.norm() / normb) <= tol) {
101  tol = resid;
102  max_iter = 0;
103  return 0;
104  }
105  Vector *v = new Vector[m+1];
106 
107  while (j <= max_iter) {
108  v[0] = r * (1.0 / beta);
109  s.set(0.0);
110  s(0) = beta;
111 
112  for (i = 0; i < m && j <= max_iter; i++, j++) {
113  w = M(A*v[i]); //M.solve(A * v[i]);
114  for (k = 0; k <= i; k++) {
115  H(k, i) = w*v[k];
116  w -= H(k, i) * v[k];
117  }
118  H(i+1, i) = w.norm();
119  v[i+1] = (w / H(i+1, i));
120 
121  for (k = 0; k < i; k++)
122  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
123 
124  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
125  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
126  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
127 
128  if ((resid = std::abs(s(i+1)) / normb) < tol) {
129  Update(x, i, H, s, v);
130  tol = resid;
131  max_iter = j;
132  // std::cout<<max_iter <<std::endl;
133  delete [] v;
134  return 0;
135  }
136  }
137  Update(x, i - 1, H, s, v);
138  r = M(b-A*x);//M.solve(b - A * x);
139  beta = r.norm();
140  if ((resid = beta / normb) < tol) {
141  tol = resid;
142  max_iter = j;
143  // std::cout<<max_iter <<std::endl;
144  delete [] v;
145  return 0;
146  }
147  }
148 
149  tol = resid;
150  delete [] v;
151  return 1;
152  }
153 }
Vector operator()(const Vector &g) const
Definition: gmres.h:30
Jacobi(const M &m)
Definition: gmres.h:24
Matrix class Matrix class.
Definition: matrix.h:28
void set(const double x)
double norm() const
Definition: vector.h:193
void GeneratePlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: gmres.h:43
unsigned GMRes(const T &A, const P &M, Vector &x, const Vector &b, int max_iter, double tol, unsigned m)
Definition: gmres.h:83
void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
Definition: gmres.h:59
void Update(Vector &x, int k, T &h, Vector &s, Vector v[])
Definition: gmres.h:67