OpenMEEG
Loading...
Searching...
No Matches
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
16namespace 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