ShyLU  Version of the Day
shylu_internal_gmres.h
1 //*****************************************************************
2 // Iterative template routine -- GMRES
3 //
4 // GMRES solves the unsymmetric linear system Ax = b using the
5 // Generalized Minimum Residual method
6 //
7 // GMRES follows the algorithm described on p. 20 of the
8 // SIAM Templates book.
9 //
10 // The return value indicates convergence within max_iter (input)
11 // iterations (0), or no convergence within max_iter iterations (1).
12 //
13 // Upon successful return, output arguments have the following values:
14 //
15 // x -- approximate solution to Ax = b
16 // max_iter -- the number of iterations performed before the
17 // tolerance was reached
18 // tol -- the residual after the final iteration
19 //
20 //*****************************************************************
21 
22 #ifndef IQR_GMRES_H
23 #define IQR_GMRES_H
24 
25 #include <cmath>
26 #include <iostream>
27 
28 #include <shylu_internal_gmres_tools.h>
29 
30 namespace IQR
31 {
32 
34 {
35  void ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y)
36  {
37  Y = X;
38  }
39 };
40 
42 template <typename Scalar>
43 void GeneratePlaneRotation(const Scalar &dx, const Scalar &dy, Scalar &cs,
44  Scalar &sn)
45 {
46  if (dy == 0.0) {
47  cs = 1.0;
48  sn = 0.0;
49  } else if (std::abs(dy) > std::abs(dx)) {
50  Scalar temp = dx / dy;
51  sn = 1.0 / std::sqrt( 1.0 + temp * temp );
52  cs = temp * sn;
53  } else {
54  Scalar temp = dy / dx;
55  cs = 1.0 / std::sqrt( 1.0 + temp * temp );
56  sn = temp * cs;
57  }
58 }
59 
62 template <typename Scalar>
63 void ApplyPlaneRotation(Scalar &dx, Scalar &dy, const Scalar &cs,
64  const Scalar &sn)
65 {
66  Scalar temp = cs * dx + sn * dy;
67  dy = -sn * dx + cs * dy;
68  dx = temp;
69 }
70 
71 
73 template < typename LocalMatrix, typename LocalVector, typename MultiVector >
74 void Update(MultiVector &x, const int k, const LocalMatrix &h,
75  const LocalVector &s, const MultiVector &v)
76 {
77  LocalVector y(s);
78 
79  // Backsolve:
80  for (int i = k; i >= 0; i--) {
81  y[i] /= h[i][i];
82  for (int j = i - 1; j >= 0; j--) {
83  y[j] -= h[j][i] * y[i];
84  }
85  }
86 
87  for (int j = 0; j <= k; j++) {
88  x.Update(y[j], *v(j), 1.0);
89  }
90 }
91 
92 template < typename Operator, typename MultiVector, typename LeftPrec,
93  typename RightPrec, typename GMRESManager, typename LocalVector,
94  typename Scalar>
95 int GMRES(const Operator &A, MultiVector &x, const MultiVector &b,
96  LeftPrec *L, RightPrec *M, GMRESManager &G, int &max_iter,
97  Scalar &tol)
98 {
99  // Storing a reference to the parallel map
100  //auto& b.Map() = b.Map();
101  //int myPID = b.Map().Comm().MyPID();
102 
103  Scalar resid;
104  int i(0), j(1), k(0);
105  // The following initial guess was wrong! Indeed it was altering the whole QR factorization
106  // initial guess from previous solves : compute x
107  //M->ApplyInverse(b, x);
108 
109  LocalVector s(G.restart + 1);
110  MultiVector w(b.Map(), 1, true);
111 
112  Scalar normb;
113  L->ApplyInverse(b, w);
114  w.Norm2(&normb);
115 
116  MultiVector t(b.Map(), 1, true);
117  A.Apply(x, t);
118  w.Update(1.0, b, -1.0, t, 0.0);
119 
120  MultiVector r(b.Map(), 1, true);
121  L->ApplyInverse(w, r);
122  Scalar beta;
123  r.Norm2(&beta);
124 
125  if (normb == 0.0) {
126  normb = 1;
127  }
128 
129  if ((resid = beta / normb) <= tol) { // qui ho migliorato
130  tol = resid;
131  max_iter = 0;
132  return 0;
133  }
134 
135  while (j <= max_iter) {
136  MultiVector* v0 = (*G.v)(0);
137  v0->Update(1.0 / beta, r, 0.0);
138  s.assign(G.restart + 1, 0.0);
139  s[0] = beta;
140 
141  for (i = 0; i < G.restart && j <= max_iter; i++, j++) {
142  M->ApplyInverse(*((*G.v)(i)), t);
143  A.Apply(t, r);
144  L->ApplyInverse(r, w);
145 
146  for (k = 0; k <= i; k++) {
147  MultiVector* vk = (*G.v)(k);
148  w.Dot(*vk, &(G.H[k][i]));
149  w.Update(-G.H[k][i], *vk, 1.0);
150  }
151  w.Norm2(&(G.H[i + 1][i]));
152  MultiVector* vi1 = (*G.v)(i + 1);
153  // Set (*G.v)(i + 1) to w/||w||
154  vi1->Scale(1.0 / G.H[i + 1][i], w);
155 
156  for (k = 0; k < i; k++) {
157  ApplyPlaneRotation(G.H[k][i], G.H[k + 1][i], G.cs[k], G.sn[k]);
158  }
159 
160  // Generate i-th Given rotation Gi. Note Qn= G0 * ... * Gn
161  GeneratePlaneRotation(G.H[i][i], G.H[i + 1][i], G.cs[i], G.sn[i]);
162  // Apply Gi
163  ApplyPlaneRotation(G.H[i][i], G.H[i + 1][i], G.cs[i], G.sn[i]);
164  ApplyPlaneRotation(s[i], s[i + 1], G.cs[i], G.sn[i]);
165 
166 // if (! myPID) std::cout << "iter: " << j << ", residual: " << resid << std::endl;
167 
168  if ((resid = abs(s[i + 1]) / normb) < tol) {
169  MultiVector y(b.Map(), 1, true);
170  Update(y, i, G.H, s, *(G.v));
171  M->ApplyInverse(y, t);
172  x.Update(1.0, t, 1.0);
173  tol = resid;
174  max_iter = j;
175  G.m = i;
176  return 2;
177  }
178  }
179 
180  MultiVector y(b.Map(), 1, true);
181  Update(y, i - 1, G.H, s, *(G.v));
182  M->ApplyInverse(y, t);
183  x.Update(1.0, t, 1.0);
184  A.Apply(x, t);
185  w.Update(1.0, b, -1.0, t, 0.0);
186  L->ApplyInverse(w, r);
187  r.Norm2(&beta);
188 
189  if ((resid = beta / normb) < tol) {
190  tol = resid;
191  max_iter = j;
192  G.m = i;
193  return 3;
194  }
195  }
196 
197  tol = resid;
198  G.m = i;
199  return 1;
200 }
201 
202 } // namespace IQR
203 
204 #endif // IQR_GMRES_H
KOKKOS_INLINE_FUNCTION RealType abs(const complex< RealType > &x)