Chombo + EB + MF  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: 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