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