ShyLU  Version of the Day
shylu_internal_gmres_tools.h
1 #ifndef IQR_GMRES_TOOLS_H
2 #define IQR_GMRES_TOOLS_H
3 
4 #include <iostream>
5 #include <shylu_internal_gmres.h>
6 
7 #include <Epetra_Operator.h>
8 
9 /*
10  This typename provides some tools to handle with gmres.
11  In particular it stores informations of a computed gmres which are useful
12  for the followings gmres, such as initial guess or preconditioner.
13  It is even possible to have nested right preconditioners using
14  the pointer to shylu_internal_gmres_tools *P.
15  Data stored are:
16  v : array of basis vectors (for "AP^-1 y" and "Ax")
17  w : for future use: same as v but basis of "x" and "P^-1y".
18  for the moment v and w are he same.
19  Q,R : QR-factorization of the hessenberg matrix
20 */
21 
22 namespace IQR
23 {
24 
25 template < typename Map, typename MultiVector, typename LocalMatrix, typename LocalVector >
27 {
28 public:
29  // Constructor and destructor
30 
31  // flex != 0 means that we will use flexible gmres and
32  // we have to use w
33  GMRESManager(const Map& map, const int &rest, const int flex = 0, const bool scaling = true);
34  ~GMRESManager();
35 
36  // Public methods
37  int initialGuess(const MultiVector &b,
38  MultiVector& X) const;
39  // solve min_{x\inV} || Ax -b||
40  int ApplyInverse(const MultiVector &b, MultiVector& X) const;
41  // delete P and restart the gmres.
42  int start();
43 
44  // Public data
45  int restart;
46  int m;
47  int isFlex;
48  bool doScaling;
49  bool isFirst;
50 
51  // Q
52  LocalVector cs;
53  LocalVector sn;
54  // R , with abuse of notation
55  LocalMatrix H;
56  // v basis of AP^-1 y = Ax and of y
57  MultiVector *v;
58  // w basis of x = P^-1 y
59  MultiVector *w;
60  GMRESManager *P; // preconditioner(s)
61  Epetra_Operator* P2;
62 
63  // Related to parallel operation
64  const Map& map_;
65 };
66 
67 template < typename Map, typename MultiVector, typename LocalMatrix, typename LocalVector >
68 GMRESManager< Map, MultiVector, LocalMatrix, LocalVector >::GMRESManager(const Map& map, const int &rest, const int flex, const bool scaling)
69  : restart(rest),
70  m(-1),
71  isFlex(flex),
72  doScaling(scaling),
73  cs(rest + 1, 0.0),
74  sn(rest + 1, 0.0),
75  H(rest + 1, LocalVector(rest, 0.0)),
76  P(0),
77  map_(map)
78 {
79  if (rest < 0) {
80  std::cout << "shylu_internal_gmres_tools: restart must be positive" << std::endl;
81  exit(1);
82  }
83 
84  // After you get it working, try removing these...
85  v = new MultiVector(map_, rest + 1, true);
86 
87  if (isFlex == 0) {
88  w = v;
89  } else {
90  w = new MultiVector(map_, rest + 1, true);
91  }
92 }
93 
94 template < typename Map, typename MultiVector, typename LocalMatrix, typename LocalVector >
96 {
97  delete v;
98  if (isFlex != 0) {
99  delete w;
100  }
101  //delete P;
102 }
103 
104 template < typename Map, typename MultiVector, typename LocalMatrix, typename LocalVector >
106  MultiVector& X) const
107 {
108  int mm(m);
109 
110  if (m > restart) {
111  std::cout << "GMRES_TOOLS.initialGuess: mm>m, not expected using m"
112  << std::endl;
113  mm = restart;
114  }
115 
116  if (P != 0 && isFlex == 0) {// if there is a (right) preconditioner, use it
117  this->ApplyInverse(b, X);
118  // x= P->initialGuess(x);
119  return 0;
120  }
121 
122  if (mm <= -1 ) {
123  X = b;
124  return 0;
125  }
126 
127  //std::cout << "GMRES_TOOLS.initialGuess: guess with m = " << mm
128  // << std::endl;
129  int i;
130 
131  LocalVector bp(mm + 1, 0.0);
132  // bp = Q_{n+1}^T b and x = x - bp Q_{n+1}
133  for (i = 0; i < mm + 1; i++) {
134  b.Dot(*(*v)(i), &bp[i]);
135  }
136 
137  // applying Q: bp = \Omega_n bp. At this time, bp = \Omega_n Q_{n+1}^T b
138  for (i = 0; i < mm; i++) {
139  ApplyPlaneRotation(bp[i], bp[i + 1], cs[i], sn[i]);
140  }
141 
142  // solving R_{n} y_{n} = bp and setting x = x + Q_{n} y_{n}
143  Update(X, mm-1, H, bp, *w);
144  //Update(y, i, G.H, , *(G.v));
145 
146  return 0;
147 }
148 
149 
150 template < typename Map, typename MultiVector, typename LocalMatrix, typename LocalVector >
152 {
153  m = -1;
154 
155  if (P != 0 ) {
156  delete P;
157  P = 0;
158  }
159 
160  return 0;
161 }
162 
163 
164 template < typename Map, typename MultiVector, typename LocalMatrix, typename LocalVector >
165 int GMRESManager< Map, MultiVector,
166  LocalMatrix, LocalVector >::ApplyInverse(const MultiVector &b,
167  MultiVector& X) const
168 {
169  int mm(m);
170 
171  if (m > restart) {
172  std::cout << "GMRES_TOOLS.solve: mm>m, not expected using m"
173  << std::endl;
174  mm = restart;
175  }
176 
177  if (mm <= -1 ) {
178  std::cout << "GMRES_TOOLS.solve: no preconditioner" << std::endl;
179  X = b;
180  return 1;
181  }
182 
183  //std::cout << "GMRES_TOOLS.solve: prec with m = " << mm << std::endl;
184 // int myPID = map_.Comm().MyPID();
185 // if (! myPID) std::cout << "Apply inverse " << id << "\n";
186 
187  int i;
188  LocalVector bp(mm + 1, 0.0);
189  MultiVector x(map_, 1, false);
190  x = b;
191 
192  // bp = Q_{n+1}^T b and x = x - bp Q_{n+1}
193  for (i = 0; i < mm + 1; i++) { // using " V is an orthogonal basis
194  MultiVector* vi = (*v)(i);
195  b.Dot(*vi, &bp[i]);
196  //std::cout << "bp: " << bp[i] << std::endl;
197  // x-= bp(i) * w[i]; // computing the orthogonal component of b
198  x.Update(-bp[i], *vi, 1.0);
199  }
200 
201  // applying Q: bp = \Omega_n bp. At this time, bp = \Omega_n Q_{n+1}^T b
202  for (i = 0; i < mm; i++) {
203  ApplyPlaneRotation(bp[i], bp[i + 1], cs[i], sn[i]);
204  }
205 
206  typedef typename LocalVector::value_type Scalar;
207  Scalar scaling(0);
208  for (i = 0; i < mm; i++)
209  {
210  scaling += H[i][i];
211  }
212  scaling = static_cast<Scalar>(mm) / scaling;
213 
214  LocalVector bq(mm + 1, 0.0);
215  bq[mm] = bp[mm];
216 
217  // bq = \omega_n^T (0...0 bp(n+1))^T
218  for (i = mm-1; i >= 0; i--) {
219  ApplyPlaneRotation(bq[i], bq[i + 1], cs[i], -sn[i]);
220 // bq[i] = - sn[i] * bq[i + 1];
221 // bq[i + 1] = cs[i] * bq[i + 1];
222  }
223 
224  // x = x + Q_n bq
225  for (i = 0; i <= mm; i++) {
226  MultiVector* vi = (*v)(i);
227  // x+= bq(i)*w[i];
228  x.Update(bq[i], *vi, 1.0);
229  }
230 
231  if (doScaling) {
232  x.Scale(scaling);
233  }
234 
235  // solving R_{n} y_{n} = bp and setting x = x + Q_{n} y_{n}
236  Update(x, mm-1, H, bp, *w); // summing to x the solution of Rx=bp
237 
238  if (isFirst) {
239  if (P2 !=0 && isFlex == 0) {
240  P2->ApplyInverse(x, X);
241  } else {
242  X = x;
243  }
244  } else {
245  if (P != 0 && isFlex == 0) {
246  P->ApplyInverse(x, X);
247  } else {
248  X = x;
249  }
250  }
251  return 0;
252 }
253 
254 } // namespace IQR
255 
256 #endif // IQR_GMRES_TOOLS_H