Chombo + EB  3.2
GMRESSolver.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 // MFA, March 2011
12 
13 #ifndef _GMRESSOLVER_H_
14 #define _GMRESSOLVER_H_
15 
16 #include "LinearSolver.H"
17 #include "parstream.H"
18 #include "CH_Timer.H"
19 #include "NamespaceHeader.H"
20 
21 ///
22 /**
23  Krylov solver using the GMRES algorithm.
24  Important parmeters are:
25  1) relative tolerance: m_reps(1.e-12)
26  2) restart length: m_restrtLen(30)
27  */
28 template <class T>
29 class GMRESSolver : public LinearSolver<T>
30 {
31 public:
32 
33  GMRESSolver();
34 
35  virtual ~GMRESSolver();
36 
37  void clearData();
38 
39  virtual void setHomogeneous(bool a_homogeneous)
40  {
41  m_homogeneous = a_homogeneous;
42  }
43 
44  ///
45  virtual void setConvergenceMetrics(Real a_metric,
46  Real a_tolerance);
47 
48  ///
49  /**
50  define the solver. a_op is the linear operator.
51  a_homogeneous is whether the solver uses homogeneous boundary
52  conditions.
53  */
54  virtual void define(LinearOp<T>* a_op, bool a_homogeneous=true);
55 
56  ///solve the equation.
57  virtual void solve(T& a_phi, const T& a_rhs);
58 
59  ///
60  /**
61  public member data: whether the solver is restricted to
62  homogeneous boundary conditions
63  */
65 
66  ///
67  /**
68  public member data: operator to solve.
69  */
71 
72  ///
73  /**
74  private member data: restart length, need to allocate stuff
75  */
76 private:
78 public:
79  int restartLen()const
80  {
81  return m_restrtLen;
82  }
83  void setRestartLen(int mm);
84  ///
85  /**
86  public member data: max iterations (eg, >= m_restrtLen)
87  */
88  int m_imax;
89  ///
90  /**
91  public member data: how much screen out put the user wants.
92  set = 0 for no output.
93  */
95 
96  ///
97  /**
98  public member data: solver tolerance
99  */
101 
102  ///
103  /**
104  public member data: relative solver tolerance
105  */
107 
108  ///
109  /**
110  public member data:
111  set = -1 if solver exited for an unknown reason
112  set = 1 if solver converged to tolerance
113  set = 2 if rho = 0
114  set = 3 if max number of restarts was reached
115  */
117 
118  ///
119  /**
120  public member data: what the algorithm should consider "close to zero"
121  */
123 
125 
126 private:
127  void allocate();
128  void CycleGMRES( T &xx,
129  int &reason, int &itcount, Real &rnorm0,
130  const bool avoidnorms = false);
131 
132  void ResidualGMRES( T &a_vv, const T &a_xx,
133  const T &a_bb, T &a_temp );
134 
135  void BuildGMRESSoln( Real nrs[], T &a_xx, const int it,
136  const T vv_0[] );
137 
138  void UpdateGMRESHessenberg( const int it, bool hapend, Real &res );
139 
140  //void ModifiedGramSchmidtOrthogonalization( const int it );
141  //void UnmodifiedGramSchmidtOrthogonalization( const int it );
142  void TwoUnmodifiedGramSchmidtOrthogonalization( const int it );
143 
144  void ApplyAB( T &a_dest, const T &a_xx, T &a_temp ) const;
145 
146  // data
150 };
151 
152 // *******************************************************
153 // GMRESSolver Implementation
154 // *******************************************************
155 
156 template <class T>
158 {
159  int max_k = m_restrtLen;
160  int hh = (max_k + 2) * (max_k + 1);
161  int hes = (max_k + 1) * (max_k + 1);
162  int rs = (max_k + 2);
163  int cc = (max_k + 1);
164  int size = (hh + hes + rs + 2*cc + 1);
165 
166  m_data = new Real[size];
167  m_hh = m_data; // hh
168  m_hes = m_data + hh; // hes_
169  m_d = m_hes + hes; // rs_
170  m_ee = m_d + rs; // cc_
171  m_dd = m_ee + cc; // ss_
172 
173  m_work_arr = 0;
174 }
175 
176 template <class T>
178 {
179  delete [] m_data;
180 }
181 
182 template <class T>
184 {
185  CH_assert(mm>0);
186  clearData();
187  m_restrtLen = mm;
188  allocate();
189 }
190 
191 template <class T>
193  :m_homogeneous(false),
194  m_op(NULL),
195  m_restrtLen(30),
196  m_imax(1000),
197  m_verbosity(3),
198  m_eps(1.0E-50),
199  m_reps(1.0E-12),
200  m_exitStatus(-1),
201  m_small(1.0E-30),
202  m_normType(2)
203 {
204  allocate();
205 }
206 
207 template <class T>
209 {
210  m_op = NULL;
211  clearData();
212 }
213 
214 template <class T>
215 void GMRESSolver<T>::define(LinearOp<T>* a_operator, bool a_homogeneous/*=true*/)
216 {
217  m_homogeneous = a_homogeneous;
218  m_op = a_operator;
219 }
220 
221 ////////////////////////////////////////////////////////////////////////
222 // GMRES
223 ///////////////////////////////////////////////////////////////////////
224 #define HH(a,b) (m_hh + (b)*(m_restrtLen+2) + (a))
225 #define HES(a,b) (m_hes + (b)*(m_restrtLen+1) + (a))
226 #define CC(a) (m_ee + (a))
227 #define SS(a) (m_dd + (a))
228 #define GRS(a) (m_d + (a))
229 
230 /* vector names */
231 #define VEC_OFFSET 2
232 #define VEC_TEMP_RHS m_work_arr[0]
233 #define VEC_TEMP_LHS m_work_arr[1]
234 #define VEC_VV(i) m_work_arr[VEC_OFFSET + i]
235 
236 ////////////////////////////////////////////////////////////////////////
237 // GMRESSolver<T>::solve
238 ///////////////////////////////////////////////////////////////////////
239 template <class T>
240 void GMRESSolver<T>::solve( T& a_xx, const T& a_bb )
241 {
242  CH_TIMERS("GMRESSolver::solve");
243 
244  CH_TIMER("GMRESSolver::solve::Initialize",timeInitialize);
245  CH_TIMER("GMRESSolver::solve::MainLoop",timeMainLoop);
246  CH_TIMER("GMRESSolver::solve::Cleanup",timeCleanup);
247 
248  CH_START(timeInitialize);
249 
250  if (m_verbosity >= 3)
251  {
252  pout() << "GMRESSolver::solve" << endl;
253  }
254 
255  const int nwork = VEC_OFFSET + m_restrtLen + 1; // flex = VEC_OFFSET + 2*(m_restrtLen + 1);
256  m_work_arr = new T[nwork];
257  m_op->create(VEC_TEMP_RHS, a_bb);
258  m_op->create(VEC_TEMP_LHS, a_xx);
259  for (int i=VEC_OFFSET;i<nwork;i++)
260  {
261  m_op->create(m_work_arr[i], a_bb);
262  }
263 
264  int it;
265  Real rnorm0 = 0.0;
266  T &vv_0 = VEC_VV(0);
267 
268  /* Compute the initial (preconditioned) residual (into 'vv_0')*/
269  m_op->assign( vv_0, a_bb );
270  m_op->setToZero( a_xx );
271 
272  CH_STOP(timeInitialize);
273 
274  ///
275  /** m_exitStatus
276  set = -1 if solver exited for an unknown reason
277  set = 1 if solver converged to tolerance
278  set = 2 if rho = 0
279  set = 3 if max number of restarts was reached
280  */
281 
282  CH_START(timeMainLoop);
283 
284  // doit
285  it = 0; m_exitStatus = -1;
286  CycleGMRES( a_xx, m_exitStatus, it, rnorm0 );
287  // loop for restarts
288  while ( m_exitStatus==-1 && it < m_imax )
289  {
290  if (m_verbosity >= 3)
291  {
292  pout() << "*";
293  }
294  ResidualGMRES( vv_0, a_xx, a_bb, VEC_TEMP_RHS );
295  CycleGMRES( a_xx, m_exitStatus, it, rnorm0 );
296  }
297  if (m_exitStatus==-1 && it >= m_imax)
298  {
299  m_exitStatus = 3;
300  }
301 
302  CH_STOP(timeMainLoop);
303 
304  CH_START(timeCleanup);
305 
306  // clean up
307  for (int i=0;i<nwork;i++)
308  {
309  m_op->clear(m_work_arr[i]);
310  }
311  delete [] m_work_arr; m_work_arr = 0;
312 
313  if (m_verbosity >= 3)
314  {
315  pout() << "GMRESSolver::solve done, status = " << m_exitStatus << endl;
316  }
317 
318  CH_STOP(timeCleanup);
319 }
320 
321 #define CONVERGED(r0,r) (r<r0*m_reps || r<m_eps)
322 
323 template <class T>
325  int &a_reason, int &a_itcount, Real &a_rnorm0,
326  const bool a_avoidnorms/*=false*/ )
327 {
328  CH_assert(m_op != 0);
329  Real res_norm,res,hapbnd,tt;
330  int it;
331  bool hapend = false;
332  T &vv_0 = VEC_VV(0);
333 
334  /* scale VEC_VV (the initial residual) */
335  //ierr = vv_0->Normalize( &res_norm ); CHKERRQ(ierr);
336  res_norm = m_op->norm( vv_0, m_normType );
337  res = res_norm;
338  *GRS(0) = res_norm;
339 
340  /* check for the convergence */
341  if ( res == 0. )
342  {
343  if (m_verbosity >= 3)
344  {
345  pout() << "GMRESSolver::solve zero residual!!!" << endl;
346  }
347  a_reason = 1; // 1 == converged,
348  return;
349  }
350  // normilize
351  tt = 1./res;
352  m_op->scale( vv_0, tt);
353 
354  if ( a_itcount == 0 ) a_rnorm0 = res;
355  bool conv = CONVERGED( a_rnorm0, res );
356  a_reason = conv ? 1 : -1;
357  it = 0;
358  while ( a_reason == -1 && it < m_restrtLen && a_itcount < m_imax )
359  {
360  if ( m_verbosity>=4 && (it!=0 || a_itcount==0) )
361  {
362  pout() << a_itcount << ") GMRES residual = " << res << endl;
363  }
364 
365  const T &vv_it = VEC_VV(it);
366  T &vv_it1 = VEC_VV(it+1);
367 
368  // apply AB
369  {
370  T &Mb = VEC_TEMP_LHS;
371  ApplyAB( vv_it1, vv_it, Mb );
372  }
373  /* update hessenberg matrix and do Gram-Schmidt */
374  TwoUnmodifiedGramSchmidtOrthogonalization(it);
375 
376  /* vv(i+1) . vv(i+1) */
377  //vv_it1->Normalize( &tt ); CHKERRQ(ierr);
378  tt = m_op->norm( vv_it1, m_normType );
379  /* check for the happy breakdown */
380  //CH_assert( *GRS(it) != 0.0 );
381  //hapbnd = Abs( tt / *GRS(it) );
382  hapbnd = 1.e-99; // hard wired happy tol!!!
383  if (tt < hapbnd)
384  {
385  if ( m_verbosity>=3 )
386  {
387  pout() << "Detected happy breakdown, "<<it<<") current hapbnd = "<< tt << endl;
388  }
389  hapend = true;
390  }
391  else
392  {
393  m_op->scale( vv_it1, 1./tt);
394  }
395 
396  /* save the magnitude */
397  *HH(it+1,it) = tt; *HES(it+1,it) = tt;
398 
399  UpdateGMRESHessenberg( it, hapend, res );
400 
401  // increment
402  it++; (a_itcount)++;
403  conv = CONVERGED( a_rnorm0, res );
404  a_reason = conv ? 1 : -1;
405  /* Catch error in happy breakdown and signal convergence and break from loop */
406  if ( hapend )
407  {
408  if ( !conv && m_verbosity>=3)
409  {
410  pout() << "You reached the happy break down, but convergence was not indicated. Residual norm=" << res << endl;
411  }
412  break;
413  }
414  }
415  /* Monitor if we know that we will not return for a restart */
416  if ( (a_reason!=1 || a_itcount>=m_imax) && m_verbosity>=4 )
417  {
418  pout() << a_itcount << ") GMRES residual = " << res << endl;
419  }
420 
421  /*
422  Down here we have to solve for the "best" coefficients of the Krylov
423  columns, add the solution values together, and possibly unwind the
424  preconditioning from the solution
425  */
426  /* Form the solution (or the solution so far) */
427  BuildGMRESSoln( GRS(0), a_xx, it-1, &VEC_VV(0) );
428 }
429 
430 template <class T>
431 void GMRESSolver<T>::ResidualGMRES( T &a_vv, const T &a_xx,
432  const T &a_bb, T &a_temp_rhs )
433 {
434  CH_assert(m_op != 0);
435  m_op->applyOp( a_temp_rhs, a_xx, true); /* t <- Ax */
436  m_op->assign( a_vv, a_bb );
437  m_op->incr( a_vv, a_temp_rhs, -1.0); /* b - A(B)x */
438 }
439 template <class T>
440 void GMRESSolver<T>::BuildGMRESSoln( Real nrs[], T &a_xx, const int it,
441  const T vv_0[] )
442 {
443  Real tt;
444  int ii,k,j;
445 
446  /* Solve for solution vector that minimizes the residual */
447 
448  /* If it is < 0, no gmres steps have been performed */
449  if (it < 0)
450  {
451  return;
452  }
453  if (*HH(it,it) == 0.0)
454  {
455  pout() << "HH(it,it) is identically zero!!! GRS(it) = " << *GRS(it) << endl;
456  }
457  if (*HH(it,it) != 0.0)
458  {
459  nrs[it] = *GRS(it) / *HH(it,it);
460  }
461  else
462  {
463  nrs[it] = 0.0;
464  }
465 
466  for (ii=1; ii<=it; ii++)
467  {
468  k = it - ii;
469  tt = *GRS(k);
470  for (j=k+1; j<=it; j++) tt = tt - *HH(k,j) * nrs[j];
471  nrs[k] = tt / *HH(k,k);
472  }
473 
474  /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
475  T &temp = VEC_TEMP_RHS;
476  m_op->setToZero(temp);
477  //temp->MAXPY( it+1, nrs, vv_0 );
478  for (ii=0; ii<it+1; ii++)
479  {
480  m_op->incr(temp, vv_0[ii], nrs[ii]);
481  }
482  /* unwind pc */
483  /*If we preconditioned on the right, we need to solve for the correction to
484  the unpreconditioned problem */
485  T &temp_matop = VEC_TEMP_LHS;
486  //ierr = pc->Apply( temp, temp_matop );CHKERRQ(ierr);
487  m_op->preCond( temp_matop, temp );
488  m_op->incr( a_xx, temp_matop, 1.0 );
489 }
490 /* GMRESSolver::UpdateGMRESHessenberg *****************************************
491  *
492  * INPUT:
493  * - it:
494  * - hapend: good breakdown?
495  * - res: residual (out)
496  *
497  * SIDE EFFECTS:
498  * - sets 'nwork_' and allocs 'work_'
499  *
500  * RETURN:
501  * - PETSc error code
502  */
503 template <class T>
504 void GMRESSolver<T>::UpdateGMRESHessenberg( const int it, bool hapend, Real &res )
505 {
506  Real *hh,*cc,*ss,tt;
507  int j;
508 
509  hh = HH(0,it);
510  cc = CC(0);
511  ss = SS(0);
512 
513  /* Apply all the previously computed plane rotations to the new column
514  of the Hessenberg matrix */
515  for (j=1; j<=it; j++)
516  {
517  tt = *hh;
518  *hh = *cc * tt + *ss * *(hh+1);
519  hh++;
520  *hh = *cc++ * *hh - (*ss++ * tt);
521  }
522 
523  /*
524  compute the new plane rotation, and apply it to:
525  1) the right-hand-side of the Hessenberg system
526  2) the new column of the Hessenberg matrix
527  thus obtaining the updated value of the residual
528  */
529  if ( !hapend )
530  {
531  tt = sqrt( *hh * *hh + *(hh+1) * *(hh+1) );
532  if (tt == 0.0)
533  {
534  pout() << "Your matrix or preconditioner is the null operator\n";
535  }
536  *cc = *hh / tt;
537  *ss = *(hh+1) / tt;
538  *GRS(it+1) = - (*ss * *GRS(it));
539  *GRS(it) = *cc * *GRS(it);
540  *hh = *cc * *hh + *ss * *(hh+1);
541  res = Abs( *GRS(it+1) );
542  }
543  else
544  {
545  /* happy breakdown: HH(it+1, it) = 0, therfore we don't need to apply
546  another rotation matrix (so RH doesn't change). The new residual is
547  always the new sine term times the residual from last time (GRS(it)),
548  but now the new sine rotation would be zero...so the residual should
549  be zero...so we will multiply "zero" by the last residual. This might
550  not be exactly what we want to do here -could just return "zero". */
551 
552  res = 0.0;
553  }
554 }
555 /*
556  This is the basic orthogonalization routine using modified Gram-Schmidt.
557 */
558 // template <class T>
559 // void GMRESSolver<T>::ModifiedGramSchmidtOrthogonalization( const int it )
560 // {
561 
562 // int ierr,j;
563 // double *hh,*hes,tmp;
564 // PromVector_base *vv_1 = VEC_VV(it+1), *vv_j;
565 // //const PromMatrix_base * const Amat = alloced_ ? A_ : alias_->A_;
566 
567 // PetscFunctionBegin;
568 // /* update Hessenberg matrix and do Gram-Schmidt */
569 // hh = HH(0,it);
570 // hes = HES(0,it);
571 // for (j=0; j<=it; j++) {
572 // /* (vv(it+1), vv(j)) */
573 // vv_j = VEC_VV(j); ierr = vv_1->Dot( vv_j, hh );CHKERRQ(ierr);
574 // *hes++ = *hh;
575 // /* vv(it+1) <- vv(it+1) - hh[it+1][j] vv(j) */
576 // tmp = - (*hh++);
577 // ierr = vv_1->AXPY( tmp, vv_j );CHKERRQ(ierr);
578 // }
579 // PetscFunctionReturn(0);
580 // }
581 
582 /*
583  This version uses UNMODIFIED Gram-Schmidt. It is NOT always recommended,
584  but it can give MUCH better performance than the default modified form
585  when running in a parallel environment.
586  */
587 
588 // int PromSolver::UnmodifiedGramSchmidtOrthogonalization( const int it )
589 // {
590 // int j,ierr;
591 // double *hh,*hes;
592 // PromVector_base *vv_1 = VEC_VV(it+1);
593 // const PromVector_base *const*vv_0 = &(VEC_VV(0));
594 // //const PromMatrix_base * const Amat = alloced_ ? A_ : alias_->A_;
595 
596 // PetscFunctionBegin;
597 // /* update Hessenberg matrix and do unmodified Gram-Schmidt */
598 // hh = HH(0,it);
599 // hes = HES(0,it);
600 // /*
601 // This is really a matrix-vector product, with the matrix stored
602 // as pointer to rows
603 // */
604 // ierr = vv_1->MDot( it+1, vv_0, hes );CHKERRQ(ierr);
605 
606 // /*
607 // This is really a matrix-vector product:
608 // [h[0],h[1],...]*[ v[0]; v[1]; ...] subtracted from v[it+1].
609 // */
610 
611 // for (j=0; j<=it; j++) hh[j] = -hes[j];
612 // ierr = vv_1->MAXPY( it+1, hh, vv_0 );CHKERRQ(ierr);
613 // for (j=0; j<=it; j++) hh[j] = -hh[j];
614 // PetscFunctionReturn(0);
615 // }
616 
617 /*
618  uses 1 iteration of iterative refinement of UNMODIFIED Gram-Schmidt.
619  It can give better performance when running in a parallel
620  environment and in some cases even in a sequential environment (because
621  MAXPY has more data reuse).
622 
623  Care is taken to accumulate the updated HH/HES values.
624  */
625 template <class T>
627 {
628  Real *hh,*hes,*lhh = 0;
629  T &vv_1 = VEC_VV(it+1);
630  const T *vv_0 = &(VEC_VV(0));
631 
632  /* Don't allocate small arrays */
633  lhh = new Real[it+1];
634 
635  /* update Hessenberg matrix and do unmodified Gram-Schmidt */
636  hh = HH(0,it);
637  hes = HES(0,it);
638 
639  /* Clear hh and hes since we will accumulate values into them */
640  for (int j=0; j<=it; j++)
641  {
642  hh[j] = 0.0;
643  hes[j] = 0.0;
644  }
645 
646  for ( int ncnt = 0 ; ncnt < 2 ; ncnt++ )
647  {
648  /*
649  This is really a matrix-vector product, with the matrix stored
650  as pointer to rows
651  */
652  //ierr = vv_1->MDot( it+1, vv_0, lhh ); CHKERRQ(ierr); /* <v,vnew> */
653 // for (int j=0; j<=it; j++) // need a vector dot product!!!
654 // {
655 // lhh[j] = m_op->dotProduct(vv_1, vv_0[j]);
656 // }
657  m_op->mDotProduct(vv_1, it+1, vv_0, lhh);
658 
659  /*
660  This is really a matrix vector product:
661  [h[0],h[1],...]*[ v[0]; v[1]; ...] subtracted from v[it+1].
662  */
663  for (int j=0; j<=it; j++) lhh[j] = - lhh[j];
664  //ierr = vv_1->MAXPY( it+1, lhh, vv_0 ); CHKERRQ(ierr);
665  for (int j=0; j<=it; j++)
666  {
667  m_op->incr(vv_1, vv_0[j], lhh[j]);
668  }
669  for (int j=0; j<=it; j++)
670  {
671  hh[j] -= lhh[j]; /* hh += <v,vnew> */
672  hes[j] += lhh[j]; /* hes += - <v,vnew> */
673  }
674  }
675 
676  delete [] lhh;
677 }
678 
679 /* PromSolver::PromPCApplyBAorAB ******************************************
680  *
681  */
682 template <class T>
683 void GMRESSolver<T>::ApplyAB( T &a_dest, const T &a_xx, T &a_tmp_lhs ) const
684 {
685  m_op->preCond( a_tmp_lhs, a_xx );
686  m_op->applyOp( a_dest, a_tmp_lhs, true);
687 }
688 
689 template <class T>
691  Real a_tolerance)
692 {
693  m_eps = a_tolerance;
694 }
695 
696 #include "NamespaceFooter.H"
697 #endif /*_GMRESSOLVER_H_*/
std::ostream & pout()
Use this in place of std::cout for program output.
#define CH_TIMERS(name)
Definition: CH_Timer.H:133
Definition: LinearSolver.H:28
#define SS(a)
Definition: GMRESSolver.H:227
LinearOp< T > * m_op
Definition: GMRESSolver.H:70
Real * m_ee
Definition: GMRESSolver.H:148
virtual void setConvergenceMetrics(Real a_metric, Real a_tolerance)
Set a convergence metric, along with solver tolerance, if desired.
Definition: GMRESSolver.H:690
#define CH_assert(cond)
Definition: CHArray.H:37
void allocate()
Definition: GMRESSolver.H:157
int m_exitStatus
Definition: GMRESSolver.H:116
void ResidualGMRES(T &a_vv, const T &a_xx, const T &a_bb, T &a_temp)
Definition: GMRESSolver.H:431
Real * m_hes
Definition: GMRESSolver.H:148
Definition: GMRESSolver.H:29
T * m_work_arr
Definition: GMRESSolver.H:149
Real * m_data
Definition: GMRESSolver.H:147
#define CH_START(tpointer)
Definition: CH_Timer.H:145
#define HH(a, b)
Definition: GMRESSolver.H:224
Real * m_d
Definition: GMRESSolver.H:148
#define CONVERGED(r0, r)
Definition: GMRESSolver.H:321
void setRestartLen(int mm)
Definition: GMRESSolver.H:183
virtual void solve(T &a_phi, const T &a_rhs)
solve the equation.
Definition: GMRESSolver.H:240
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:63
#define VEC_VV(i)
Definition: GMRESSolver.H:234
Real m_eps
Definition: GMRESSolver.H:100
void clearData()
Definition: GMRESSolver.H:177
double Real
Definition: REAL.H:33
T Abs(const T &a_a)
Definition: Misc.H:53
virtual ~GMRESSolver()
Definition: GMRESSolver.H:208
virtual void setHomogeneous(bool a_homogeneous)
Definition: GMRESSolver.H:39
#define HES(a, b)
Definition: GMRESSolver.H:225
void TwoUnmodifiedGramSchmidtOrthogonalization(const int it)
Definition: GMRESSolver.H:626
#define VEC_TEMP_LHS
Definition: GMRESSolver.H:233
Real * m_dd
Definition: GMRESSolver.H:148
Real m_small
Definition: GMRESSolver.H:122
Real * m_hh
Definition: GMRESSolver.H:148
int m_verbosity
Definition: GMRESSolver.H:94
virtual void define(LinearOp< T > *a_op, bool a_homogeneous=true)
Definition: GMRESSolver.H:215
Definition: LinearSolver.H:156
#define GRS(a)
Definition: GMRESSolver.H:228
bool m_homogeneous
Definition: GMRESSolver.H:64
Real m_reps
Definition: GMRESSolver.H:106
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
int m_imax
Definition: GMRESSolver.H:88
GMRESSolver()
Definition: GMRESSolver.H:192
int m_normType
Definition: GMRESSolver.H:124
int restartLen() const
Definition: GMRESSolver.H:79
#define VEC_OFFSET
Definition: GMRESSolver.H:231
void CycleGMRES(T &xx, int &reason, int &itcount, Real &rnorm0, const bool avoidnorms=false)
Definition: GMRESSolver.H:324
void BuildGMRESSoln(Real nrs[], T &a_xx, const int it, const T vv_0[])
Definition: GMRESSolver.H:440
void UpdateGMRESHessenberg(const int it, bool hapend, Real &res)
Definition: GMRESSolver.H:504
#define CC(a)
Definition: GMRESSolver.H:226
int m_restrtLen
Definition: GMRESSolver.H:77
void ApplyAB(T &a_dest, const T &a_xx, T &a_temp) const
Definition: GMRESSolver.H:683
#define VEC_TEMP_RHS
Definition: GMRESSolver.H:232