dune-istl  2.6-git
solvers.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_ISTL_SOLVERS_HH
5 #define DUNE_ISTL_SOLVERS_HH
6 
7 #include <array>
8 #include <cmath>
9 #include <complex>
10 #include <iostream>
11 #include <memory>
12 #include <type_traits>
13 #include <vector>
14 
15 #include <dune/common/conditional.hh>
16 #include <dune/common/exceptions.hh>
17 #include <dune/common/hybridutilities.hh>
18 #include <dune/common/rangeutilities.hh>
19 #include <dune/common/std/type_traits.hh>
20 #include <dune/common/timer.hh>
21 
22 #include <dune/istl/allocator.hh>
23 #include <dune/istl/bcrsmatrix.hh>
26 #include <dune/istl/operators.hh>
29 #include <dune/istl/solver.hh>
30 
31 namespace Dune {
43  //=====================================================================
44  // Implementation of this interface
45  //=====================================================================
46 
55  template<class X>
56  class LoopSolver : public IterativeSolver<X,X> {
57  public:
58  using typename IterativeSolver<X,X>::domain_type;
59  using typename IterativeSolver<X,X>::range_type;
60  using typename IterativeSolver<X,X>::field_type;
61  using typename IterativeSolver<X,X>::real_type;
62 
63  // copy base class constructors
65 
66  // don't shadow four-argument version of apply defined in the base class
68 
70  virtual void apply (X& x, X& b, InverseOperatorResult& res)
71  {
72  // clear solver statistics
73  res.clear();
74 
75  // start a timer
76  Timer watch;
77 
78  // prepare preconditioner
79  _prec->pre(x,b);
80 
81  // overwrite b with defect
82  _op->applyscaleadd(-1,x,b);
83 
84  // compute norm, \todo parallelization
85  real_type def0 = _sp->norm(b);
86 
87  // printing
88  if (_verbose>0)
89  {
90  std::cout << "=== LoopSolver" << std::endl;
91  if (_verbose>1)
92  {
93  this->printHeader(std::cout);
94  this->printOutput(std::cout,0,def0);
95  }
96  }
97 
98  // allocate correction vector
99  X v(x);
100 
101  // iteration loop
102  int i=1; real_type def=def0;
103  for ( ; i<=_maxit; i++ )
104  {
105  v = 0; // clear correction
106  _prec->apply(v,b); // apply preconditioner
107  x += v; // update solution
108  _op->applyscaleadd(-1,v,b); // update defect
109  real_type defnew=_sp->norm(b); // comp defect norm
110  if (_verbose>1) // print
111  this->printOutput(std::cout,i,defnew,def);
112  //std::cout << i << " " << defnew << " " << defnew/def << std::endl;
113  def = defnew; // update norm
114  if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
115  {
116  res.converged = true;
117  break;
118  }
119  }
120 
121  //correct i which is wrong if convergence was not achieved.
122  i=std::min(_maxit,i);
123 
124  // print
125  if (_verbose==1)
126  this->printOutput(std::cout,i,def);
127 
128  // postprocess preconditioner
129  _prec->post(x);
130 
131  // fill statistics
132  res.iterations = i;
133  res.reduction = static_cast<double>(max_value(def/def0));
134  res.conv_rate = pow(res.reduction,1.0/i);
135  res.elapsed = watch.elapsed();
136 
137  // final print
138  if (_verbose>0)
139  {
140  std::cout << "=== rate=" << res.conv_rate
141  << ", T=" << res.elapsed
142  << ", TIT=" << res.elapsed/i
143  << ", IT=" << i << std::endl;
144  }
145  }
146 
147  protected:
154  };
155 
156 
157  // all these solvers are taken from the SUMO library
159  template<class X>
160  class GradientSolver : public IterativeSolver<X,X> {
161  public:
162  using typename IterativeSolver<X,X>::domain_type;
163  using typename IterativeSolver<X,X>::range_type;
164  using typename IterativeSolver<X,X>::field_type;
165  using typename IterativeSolver<X,X>::real_type;
166 
167  // copy base class constructors
169 
170  // don't shadow four-argument version of apply defined in the base class
172 
178  virtual void apply (X& x, X& b, InverseOperatorResult& res)
179  {
180  res.clear(); // clear solver statistics
181  Timer watch; // start a timer
182  _prec->pre(x,b); // prepare preconditioner
183  _op->applyscaleadd(-1,x,b); // overwrite b with defect
184 
185  X p(x); // create local vectors
186  X q(b);
187 
188  real_type def0 = _sp->norm(b); // compute norm
189 
190  if (_verbose>0) // printing
191  {
192  std::cout << "=== GradientSolver" << std::endl;
193  if (_verbose>1)
194  {
195  this->printHeader(std::cout);
196  this->printOutput(std::cout,0,def0);
197  }
198  }
199 
200  int i=1; real_type def=def0; // loop variables
201  field_type lambda;
202  for ( ; i<=_maxit; i++ )
203  {
204  p = 0; // clear correction
205  _prec->apply(p,b); // apply preconditioner
206  _op->apply(p,q); // q=Ap
207  lambda = _sp->dot(p,b)/_sp->dot(q,p); // minimization
208  x.axpy(lambda,p); // update solution
209  b.axpy(-lambda,q); // update defect
210 
211  real_type defnew=_sp->norm(b); // comp defect norm
212  if (_verbose>1) // print
213  this->printOutput(std::cout,i,defnew,def);
214 
215  def = defnew; // update norm
216  if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
217  {
218  res.converged = true;
219  break;
220  }
221  }
222 
223  //correct i which is wrong if convergence was not achieved.
224  i=std::min(_maxit,i);
225 
226  if (_verbose==1) // printing for non verbose
227  this->printOutput(std::cout,i,def);
228 
229  _prec->post(x); // postprocess preconditioner
230  res.iterations = i; // fill statistics
231  res.reduction = static_cast<double>(max_value(def/def0));
232  res.conv_rate = pow(res.reduction,1.0/i);
233  res.elapsed = watch.elapsed();
234  if (_verbose>0) // final print
235  std::cout << "=== rate=" << res.conv_rate
236  << ", T=" << res.elapsed
237  << ", TIT=" << res.elapsed/i
238  << ", IT=" << i << std::endl;
239  }
240 
241  protected:
248  };
249 
250 
251 
253  template<class X>
254  class CGSolver : public IterativeSolver<X,X> {
255  public:
256  using typename IterativeSolver<X,X>::domain_type;
257  using typename IterativeSolver<X,X>::range_type;
258  using typename IterativeSolver<X,X>::field_type;
259  using typename IterativeSolver<X,X>::real_type;
260 
261  // copy base class constructors
263 
264  private:
266 
267  protected:
268 
269  using enableConditionEstimate_t = Dune::Std::bool_constant<(std::is_same<field_type,float>::value || std::is_same<field_type,double>::value)>;
270 
271  public:
272 
273  // don't shadow four-argument version of apply defined in the base class
275 
284  scalar_real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, prec, reduction, maxit, verbose),
285  condition_estimate_(condition_estimate)
286  {
287  if (condition_estimate && !(enableConditionEstimate_t{})) {
288  condition_estimate_ = false;
289  std::cerr << "WARNING: Condition estimate was disabled. It is only available for double and float field types!" << std::endl;
290  }
291  }
292 
301  scalar_real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
302  condition_estimate_(condition_estimate)
303  {
304  if (condition_estimate && !(std::is_same<field_type,float>::value || std::is_same<field_type,double>::value)) {
305  condition_estimate_ = false;
306  std::cerr << "WARNING: Condition estimate was disabled. It is only available for double and float field types!" << std::endl;
307  }
308  }
309 
310 
322  virtual void apply (X& x, X& b, InverseOperatorResult& res)
323  {
324  using std::isfinite;
325 
326  res.clear(); // clear solver statistics
327  Timer watch; // start a timer
328  _prec->pre(x,b); // prepare preconditioner
329  _op->applyscaleadd(-1,x,b); // overwrite b with defect
330 
331  X p(x); // the search direction
332  X q(x); // a temporary vector
333 
334  real_type def0 = _sp->norm(b); // compute norm
335 
336  if (!all_true(isfinite(def0))) // check for inf or NaN
337  {
338  if (_verbose>0)
339  std::cout << "=== CGSolver: abort due to infinite or NaN initial defect"
340  << std::endl;
341  DUNE_THROW(SolverAbort, "CGSolver: initial defect=" << def0
342  << " is infinite or NaN");
343  }
344 
345  if (max_value(def0)<1E-30) // convergence check
346  {
347  res.converged = true;
348  res.iterations = 0; // fill statistics
349  res.reduction = 0;
350  res.conv_rate = 0;
351  res.elapsed=0;
352  if (_verbose>0) // final print
353  std::cout << "=== rate=" << res.conv_rate
354  << ", T=" << res.elapsed << ", TIT=" << res.elapsed
355  << ", IT=0" << std::endl;
356  return;
357  }
358 
359  if (_verbose>0) // printing
360  {
361  std::cout << "=== CGSolver" << std::endl;
362  if (_verbose>1) {
363  this->printHeader(std::cout);
364  this->printOutput(std::cout,0,def0);
365  }
366  }
367 
368  // Remember lambda and beta values for condition estimate
369  std::vector<real_type> lambdas(0);
370  std::vector<real_type> betas(0);
371 
372  // some local variables
373  real_type def=def0; // loop variables
374  field_type rho,rholast,lambda,alpha,beta;
375 
376  // determine initial search direction
377  p = 0; // clear correction
378  _prec->apply(p,b); // apply preconditioner
379  rholast = _sp->dot(p,b); // orthogonalization
380 
381  // the loop
382  int i=1;
383  for ( ; i<=_maxit; i++ )
384  {
385  // minimize in given search direction p
386  _op->apply(p,q); // q=Ap
387  alpha = _sp->dot(p,q); // scalar product
388  lambda = rholast/alpha; // minimization
389  Hybrid::ifElse(enableConditionEstimate_t{}, [&](auto id) {
390  if (condition_estimate_)
391  lambdas.push_back(std::real(id(lambda)));
392  });
393  x.axpy(lambda,p); // update solution
394  b.axpy(-lambda,q); // update defect
395 
396  // convergence test
397  real_type defnew=_sp->norm(b); // comp defect norm
398 
399  if (_verbose>1) // print
400  this->printOutput(std::cout,i,defnew,def);
401 
402  def = defnew; // update norm
403  if (!all_true(isfinite(def))) // check for inf or NaN
404  {
405  if (_verbose>0)
406  std::cout << "=== CGSolver: abort due to infinite or NaN defect"
407  << std::endl;
408  DUNE_THROW(SolverAbort,
409  "CGSolver: defect=" << def << " is infinite or NaN");
410  }
411 
412  if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
413  {
414  res.converged = true;
415  break;
416  }
417 
418  // determine new search direction
419  q = 0; // clear correction
420  _prec->apply(q,b); // apply preconditioner
421  rho = _sp->dot(q,b); // orthogonalization
422  beta = rho/rholast; // scaling factor
423  Hybrid::ifElse(enableConditionEstimate_t{}, [&](auto id) {
424  if (condition_estimate_)
425  betas.push_back(std::real(id(beta)));
426  });
427  p *= beta; // scale old search direction
428  p += q; // orthogonalization with correction
429  rholast = rho; // remember rho for recurrence
430  }
431 
432  //correct i which is wrong if convergence was not achieved.
433  i=std::min(_maxit,i);
434 
435  if (_verbose==1) // printing for non verbose
436  this->printOutput(std::cout,i,def);
437 
438  _prec->post(x); // postprocess preconditioner
439  res.iterations = i; // fill statistics
440  res.reduction = static_cast<double>(max_value(max_value(def/def0)));
441  res.conv_rate = pow(res.reduction,1.0/i);
442  res.elapsed = watch.elapsed();
443 
444  if (_verbose>0) // final print
445  {
446  std::cout << "=== rate=" << res.conv_rate
447  << ", T=" << res.elapsed
448  << ", TIT=" << res.elapsed/i
449  << ", IT=" << i << std::endl;
450  }
451 
452 
453  if (condition_estimate_) {
454 #if HAVE_ARPACKPP
455  Hybrid::ifElse(enableConditionEstimate_t{}, [&](auto id) {
456 
457  // Build T matrix which has extreme eigenvalues approximating
458  // those of the original system
459  // (see Y. Saad, Iterative methods for sparse linear systems)
460 
461  COND_MAT T(i, i, COND_MAT::row_wise);
462 
463  for (auto row = T.createbegin(); row != T.createend(); ++row) {
464  if (row.index() > 0)
465  row.insert(row.index()-1);
466  row.insert(row.index());
467  if (row.index() < T.N() - 1)
468  row.insert(row.index()+1);
469  }
470  for (int row = 0; row < i; ++row) {
471  if (row > 0) {
472  T[row][row-1] = std::sqrt(id(betas[row-1])) / lambdas[row-1];
473  }
474 
475  T[row][row] = 1.0 / id(lambdas[row]);
476  if (row > 0) {
477  T[row][row] += betas[row-1] / lambdas[row-1];
478  }
479 
480  if (row < i - 1) {
481  T[row][row+1] = std::sqrt(id(betas[row])) / lambdas[row];
482  }
483  }
484 
485  // Compute largest and smallest eigenvalue of T matrix and return as estimate
487 
488  real_type eps = 0.0;
489  COND_VEC eigv;
490  real_type min_eigv, max_eigv;
491  arpack.computeSymMinMagnitude (eps, eigv, min_eigv);
492  arpack.computeSymMaxMagnitude (eps, eigv, max_eigv);
493 
494  res.condition_estimate = id(max_eigv / min_eigv);
495 
496  if (this->_verbose > 0) {
497  std::cout << "Min eigv estimate: " << min_eigv << std::endl;
498  std::cout << "Max eigv estimate: " << max_eigv << std::endl;
499  std::cout << "Condition estimate: " << max_eigv / min_eigv << std::endl;
500  }
501  });
502 #else
503  std::cerr << "WARNING: Condition estimate was requested. This requires ARPACK, but ARPACK was not found!" << std::endl;
504 #endif
505  }
506  }
507 
508  private:
509  bool condition_estimate_ = false;
510 
511  // Matrix and vector types used for condition estimate
514 
515  protected:
522  };
523 
524 
525  // Ronald Kriemanns BiCG-STAB implementation from Sumo
527  template<class X>
528  class BiCGSTABSolver : public IterativeSolver<X,X> {
529  public:
530  using typename IterativeSolver<X,X>::domain_type;
531  using typename IterativeSolver<X,X>::range_type;
532  using typename IterativeSolver<X,X>::field_type;
533  using typename IterativeSolver<X,X>::real_type;
534 
535  // copy base class constructors
537 
538  // don't shadow four-argument version of apply defined in the base class
540 
548  virtual void apply (X& x, X& b, InverseOperatorResult& res)
549  {
550  using std::abs;
551  const real_type EPSILON=1e-80;
552  using std::abs;
553  double it;
554  field_type rho, rho_new, alpha, beta, h, omega;
555  real_type norm, norm_old, norm_0;
556 
557  //
558  // get vectors and matrix
559  //
560  X& r=b;
561  X p(x);
562  X v(x);
563  X t(x);
564  X y(x);
565  X rt(x);
566 
567  //
568  // begin iteration
569  //
570 
571  // r = r - Ax; rt = r
572  res.clear(); // clear solver statistics
573  Timer watch; // start a timer
574  _prec->pre(x,r); // prepare preconditioner
575  _op->applyscaleadd(-1,x,r); // overwrite b with defect
576 
577  rt=r;
578 
579  norm = norm_old = norm_0 = _sp->norm(r);
580 
581  p=0;
582  v=0;
583 
584  rho = 1;
585  alpha = 1;
586  omega = 1;
587 
588  if (_verbose>0) // printing
589  {
590  std::cout << "=== BiCGSTABSolver" << std::endl;
591  if (_verbose>1)
592  {
593  this->printHeader(std::cout);
594  this->printOutput(std::cout,0,norm_0);
595  //std::cout << " Iter Defect Rate" << std::endl;
596  //std::cout << " 0" << std::setw(14) << norm_0 << std::endl;
597  }
598  }
599 
600  if ( all_true(norm < (_reduction * norm_0)) || max_value(norm)<1E-30)
601  {
602  res.converged = 1;
603  _prec->post(x); // postprocess preconditioner
604  res.iterations = 0; // fill statistics
605  res.reduction = 0;
606  res.conv_rate = 0;
607  res.elapsed = watch.elapsed();
608  return;
609  }
610 
611  //
612  // iteration
613  //
614 
615  for (it = 0.5; it < _maxit; it+=.5)
616  {
617  //
618  // preprocess, set vecsizes etc.
619  //
620 
621  // rho_new = < rt , r >
622  rho_new = _sp->dot(rt,r);
623 
624  // look if breakdown occurred
625  if (all_true(abs(rho) <= EPSILON))
626  DUNE_THROW(SolverAbort,"breakdown in BiCGSTAB - rho "
627  << rho << " <= EPSILON " << max_value(EPSILON)
628  << " after " << it << " iterations");
629  if (all_true(abs(omega) <= EPSILON))
630  DUNE_THROW(SolverAbort,"breakdown in BiCGSTAB - omega "
631  << omega << " <= EPSILON " << max_value(EPSILON)
632  << " after " << it << " iterations");
633 
634 
635  if (it<1)
636  p = r;
637  else
638  {
639  beta = ( rho_new / rho ) * ( alpha / omega );
640  p.axpy(-omega,v); // p = r + beta (p - omega*v)
641  p *= beta;
642  p += r;
643  }
644 
645  // y = W^-1 * p
646  y = 0;
647  _prec->apply(y,p); // apply preconditioner
648 
649  // v = A * y
650  _op->apply(y,v);
651 
652  // alpha = rho_new / < rt, v >
653  h = _sp->dot(rt,v);
654 
655  if ( all_true(abs(h) < EPSILON) )
656  DUNE_THROW(SolverAbort,"abs(h) < EPSILON in BiCGSTAB - abs(h) "
657  << abs(h) << " < EPSILON " << max_value(EPSILON)
658  << " after " << it << " iterations");
659 
660  alpha = rho_new / h;
661 
662  // apply first correction to x
663  // x <- x + alpha y
664  x.axpy(alpha,y);
665 
666  // r = r - alpha*v
667  r.axpy(-alpha,v);
668 
669  //
670  // test stop criteria
671  //
672 
673  norm = _sp->norm(r);
674 
675  if (_verbose>1) // print
676  {
677  this->printOutput(std::cout,it,norm,norm_old);
678  }
679 
680  if ( all_true(norm < (_reduction * norm_0)) )
681  {
682  res.converged = 1;
683  break;
684  }
685  it+=.5;
686 
687  norm_old = norm;
688 
689  // y = W^-1 * r
690  y = 0;
691  _prec->apply(y,r);
692 
693  // t = A * y
694  _op->apply(y,t);
695 
696  // omega = < t, r > / < t, t >
697  omega = _sp->dot(t,r)/_sp->dot(t,t);
698 
699  // apply second correction to x
700  // x <- x + omega y
701  x.axpy(omega,y);
702 
703  // r = s - omega*t (remember : r = s)
704  r.axpy(-omega,t);
705 
706  rho = rho_new;
707 
708  //
709  // test stop criteria
710  //
711 
712  norm = _sp->norm(r);
713 
714  if (_verbose > 1) // print
715  {
716  this->printOutput(std::cout,it,norm,norm_old);
717  }
718 
719  if ( all_true(norm < (_reduction * norm_0)) || max_value(norm)<1E-30)
720  {
721  res.converged = 1;
722  break;
723  }
724 
725  norm_old = norm;
726  } // end for
727 
728  //correct i which is wrong if convergence was not achieved.
729  it=std::min(static_cast<double>(_maxit),it);
730 
731  if (_verbose==1) // printing for non verbose
732  this->printOutput(std::cout,it,norm);
733 
734  _prec->post(x); // postprocess preconditioner
735  res.iterations = static_cast<int>(std::ceil(it)); // fill statistics
736  res.reduction = static_cast<double>(max_value(norm/norm_0));
737  res.conv_rate = pow(res.reduction,1.0/it);
738  res.elapsed = watch.elapsed();
739  if (_verbose>0) // final print
740  std::cout << "=== rate=" << res.conv_rate
741  << ", T=" << res.elapsed
742  << ", TIT=" << res.elapsed/it
743  << ", IT=" << it << std::endl;
744  }
745 
746  protected:
753  };
754 
761  template<class X>
762  class MINRESSolver : public IterativeSolver<X,X> {
763  public:
764  using typename IterativeSolver<X,X>::domain_type;
765  using typename IterativeSolver<X,X>::range_type;
766  using typename IterativeSolver<X,X>::field_type;
767  using typename IterativeSolver<X,X>::real_type;
768 
769  // copy base class constructors
771 
772  // don't shadow four-argument version of apply defined in the base class
774 
780  virtual void apply (X& x, X& b, InverseOperatorResult& res)
781  {
782  using std::sqrt;
783  using std::abs;
784  // clear solver statistics
785  res.clear();
786  // start a timer
787  Dune::Timer watch;
788  watch.reset();
789  // prepare preconditioner
790  _prec->pre(x,b);
791  // overwrite rhs with defect
792  _op->applyscaleadd(-1,x,b);
793 
794  // compute residual norm
795  real_type def0 = _sp->norm(b);
796 
797  // printing
798  if(_verbose > 0) {
799  std::cout << "=== MINRESSolver" << std::endl;
800  if(_verbose > 1) {
801  this->printHeader(std::cout);
802  this->printOutput(std::cout,0,def0);
803  }
804  }
805 
806  // check for convergence
807  if( max_value(def0) < 1e-30 ) {
808  res.converged = true;
809  res.iterations = 0;
810  res.reduction = 0;
811  res.conv_rate = 0;
812  res.elapsed = 0.0;
813  // final print
814  if(_verbose > 0)
815  std::cout << "=== rate=" << res.conv_rate
816  << ", T=" << res.elapsed
817  << ", TIT=" << res.elapsed
818  << ", IT=" << res.iterations
819  << std::endl;
820  return;
821  }
822 
823  // the defect norm
824  real_type def = def0;
825  // recurrence coefficients as computed in Lanczos algorithm
826  field_type alpha, beta;
827  // diagonal entries of givens rotation
828  std::array<real_type,2> c{{0.0,0.0}};
829  // off-diagonal entries of givens rotation
830  std::array<field_type,2> s{{0.0,0.0}};
831 
832  // recurrence coefficients (column k of tridiag matrix T_k)
833  std::array<field_type,3> T{{0.0,0.0,0.0}};
834 
835  // the rhs vector of the min problem
836  std::array<field_type,2> xi{{1.0,0.0}};
837 
838  // some temporary vectors
839  X z(b), dummy(b);
840 
841  // initialize and clear correction
842  z = 0.0;
843  _prec->apply(z,b);
844 
845  // beta is real and positive in exact arithmetic
846  // since it is the norm of the basis vectors (in unpreconditioned case)
847  beta = sqrt(_sp->dot(b,z));
848  field_type beta0 = beta;
849 
850  // the search directions
851  std::array<X,3> p{{b,b,b}};
852  p[0] = 0.0;
853  p[1] = 0.0;
854  p[2] = 0.0;
855 
856  // orthonormal basis vectors (in unpreconditioned case)
857  std::array<X,3> q{{b,b,b}};
858  q[0] = 0.0;
859  q[1] *= real_type(1.0)/beta;
860  q[2] = 0.0;
861 
862  z *= real_type(1.0)/beta;
863 
864  // the loop
865  int i = 1;
866  for( ; i<=_maxit; i++) {
867 
868  dummy = z;
869  int i1 = i%3,
870  i0 = (i1+2)%3,
871  i2 = (i1+1)%3;
872 
873  // symmetrically preconditioned Lanczos algorithm (see Greenbaum p.121)
874  _op->apply(z,q[i2]); // q[i2] = Az
875  q[i2].axpy(-beta,q[i0]);
876  // alpha is real since it is the diagonal entry of the hermitian tridiagonal matrix
877  // from the Lanczos Algorithm
878  // so the order in the scalar product doesn't matter even for the complex case
879  alpha = _sp->dot(z,q[i2]);
880  q[i2].axpy(-alpha,q[i1]);
881 
882  z = 0.0;
883  _prec->apply(z,q[i2]);
884 
885  // beta is real and positive in exact arithmetic
886  // since it is the norm of the basis vectors (in unpreconditioned case)
887  beta = sqrt(_sp->dot(q[i2],z));
888 
889  q[i2] *= real_type(1.0)/beta;
890  z *= real_type(1.0)/beta;
891 
892  // QR Factorization of recurrence coefficient matrix
893  // apply previous givens rotations to last column of T
894  T[1] = T[2];
895  if(i>2) {
896  T[0] = s[i%2]*T[1];
897  T[1] = c[i%2]*T[1];
898  }
899  if(i>1) {
900  T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
901  T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
902  }
903  else
904  T[2] = alpha;
905 
906  // update QR factorization
907  generateGivensRotation(T[2],beta,c[i%2],s[i%2]);
908  // to last column of T_k
909  T[2] = c[i%2]*T[2] + s[i%2]*beta;
910  // and to the rhs xi of the min problem
911  xi[i%2] = -s[i%2]*xi[(i+1)%2];
912  xi[(i+1)%2] *= c[i%2];
913 
914  // compute correction direction
915  p[i2] = dummy;
916  p[i2].axpy(-T[1],p[i1]);
917  p[i2].axpy(-T[0],p[i0]);
918  p[i2] *= real_type(1.0)/T[2];
919 
920  // apply correction/update solution
921  x.axpy(beta0*xi[(i+1)%2],p[i2]);
922 
923  // remember beta_old
924  T[2] = beta;
925 
926  // check for convergence
927  // the last entry in the rhs of the min-problem is the residual
928  real_type defnew = abs(beta0*xi[i%2]);
929 
930  if(_verbose > 1)
931  this->printOutput(std::cout,i,defnew,def);
932 
933  def = defnew;
934  if(all_true(def < def0*_reduction)
935  || max_value(def) < 1e-30 || i == _maxit ) {
936  res.converged = true;
937  break;
938  }
939  } // end for
940 
941  if(_verbose == 1)
942  this->printOutput(std::cout,i,def);
943 
944  // postprocess preconditioner
945  _prec->post(x);
946  // fill statistics
947  res.iterations = i;
948  res.reduction = static_cast<double>(max_value(def/def0));
949  res.conv_rate = pow(res.reduction,1.0/i);
950  res.elapsed = watch.elapsed();
951 
952  // final print
953  if(_verbose > 0) {
954  std::cout << "=== rate=" << res.conv_rate
955  << ", T=" << res.elapsed
956  << ", TIT=" << res.elapsed/i
957  << ", IT=" << i << std::endl;
958  }
959  }
960 
961  private:
962 
963  // helper function to extract the real value of a real or complex number
964  inline
965  real_type to_real(const real_type & v)
966  {
967  return v;
968  }
969 
970  inline
971  real_type to_real(const std::complex<real_type> & v)
972  {
973  return v.real();
974  }
975 
976  void generateGivensRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
977  {
978  using std::sqrt;
979  using std::abs;
980  using std::max;
981  using std::min;
982  const real_type eps = 1e-15;
983  real_type norm_dx = abs(dx);
984  real_type norm_dy = abs(dy);
985  real_type norm_max = max(norm_dx, norm_dy);
986  real_type norm_min = min(norm_dx, norm_dy);
987  real_type temp = norm_min/norm_max;
988  // we rewrite the code in a vectorizable fashion
989  cs = cond(norm_dy < eps,
990  real_type(1.0),
991  cond(norm_dx < eps,
992  real_type(0.0),
993  cond(norm_dy > norm_dx,
994  real_type(1.0)/sqrt(real_type(1.0) + temp*temp)*temp,
995  real_type(1.0)/sqrt(real_type(1.0) + temp*temp)
996  )));
997  sn = cond(norm_dy < eps,
998  field_type(0.0),
999  cond(norm_dx < eps,
1000  field_type(1.0),
1001  cond(norm_dy > norm_dx,
1002  // dy and dx are real in exact arithmetic
1003  // thus dx*dy is real so we can explicitly enforce it
1004  field_type(1.0)/sqrt(real_type(1.0) + temp*temp)*dx*dy/norm_dx/norm_dy,
1005  // dy and dx is real in exact arithmetic
1006  // so we don't have to conjugate both of them
1007  field_type(1.0)/sqrt(real_type(1.0) + temp*temp)*dy/dx
1008  )));
1009  }
1010 
1011  protected:
1018  };
1019 
1033  template<class X, class Y=X, class F = Y>
1035  {
1036  public:
1037  using typename IterativeSolver<X,Y>::domain_type;
1038  using typename IterativeSolver<X,Y>::range_type;
1039  using typename IterativeSolver<X,Y>::field_type;
1040  using typename IterativeSolver<X,Y>::real_type;
1041 
1042  private:
1044 
1046  using fAlloc = ReboundAllocatorType<X,field_type>;
1048  using rAlloc = ReboundAllocatorType<X,real_type>;
1049 
1050  public:
1051 
1058  RestartedGMResSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
1059  IterativeSolver<X,Y>::IterativeSolver(op,prec,reduction,maxit,verbose),
1060  _restart(restart)
1061  {}
1062 
1069  RestartedGMResSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
1070  IterativeSolver<X,Y>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
1071  _restart(restart)
1072  {}
1073 
1082  virtual void apply (X& x, Y& b, InverseOperatorResult& res)
1083  {
1084  apply(x,b,max_value(_reduction),res);
1085  }
1086 
1095  virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res)
1096  {
1097  using std::abs;
1098  const real_type EPSILON = 1e-80;
1099  const int m = _restart;
1100  real_type norm, norm_old = 0.0, norm_0;
1101  int j = 1;
1102  std::vector<field_type,fAlloc> s(m+1), sn(m);
1103  std::vector<real_type,rAlloc> cs(m);
1104  // need copy of rhs if GMRes has to be restarted
1105  Y b2(b);
1106  // helper vector
1107  Y w(b);
1108  std::vector< std::vector<field_type,fAlloc> > H(m+1,s);
1109  std::vector<F> v(m+1,b);
1110 
1111  // start timer
1112  Dune::Timer watch;
1113  watch.reset();
1114 
1115  // clear solver statistics and set res.converged to false
1116  res.clear();
1117  _prec->pre(x,b);
1118 
1119  // calculate defect and overwrite rhs with it
1120  _op->applyscaleadd(-1.0,x,b); // b -= Ax
1121  // calculate preconditioned defect
1122  v[0] = 0.0; _prec->apply(v[0],b); // r = W^-1 b
1123  norm_0 = _sp->norm(v[0]);
1124  norm = norm_0;
1125  norm_old = norm;
1126 
1127  // print header
1128  if(_verbose > 0)
1129  {
1130  std::cout << "=== RestartedGMResSolver" << std::endl;
1131  if(_verbose > 1) {
1132  this->printHeader(std::cout);
1133  this->printOutput(std::cout,0,norm_0);
1134  }
1135  }
1136 
1137  if(all_true(norm_0 < EPSILON)) {
1138  _prec->post(x);
1139  res.converged = true;
1140  if(_verbose > 0) // final print
1141  print_result(res);
1142  }
1143 
1144  while(j <= _maxit && res.converged != true) {
1145 
1146  int i = 0;
1147  v[0] *= real_type(1.0)/norm;
1148  s[0] = norm;
1149  for(i=1; i<m+1; i++)
1150  s[i] = 0.0;
1151 
1152  for(i=0; i < m && j <= _maxit && res.converged != true; i++, j++) {
1153  w = 0.0;
1154  // use v[i+1] as temporary vector
1155  v[i+1] = 0.0;
1156  // do Arnoldi algorithm
1157  _op->apply(v[i],v[i+1]);
1158  _prec->apply(w,v[i+1]);
1159  for(int k=0; k<i+1; k++) {
1160  // notice that _sp->dot(v[k],w) = v[k]\adjoint w
1161  // so one has to pay attention to the order
1162  // in the scalar product for the complex case
1163  // doing the modified Gram-Schmidt algorithm
1164  H[k][i] = _sp->dot(v[k],w);
1165  // w -= H[k][i] * v[k]
1166  w.axpy(-H[k][i],v[k]);
1167  }
1168  H[i+1][i] = _sp->norm(w);
1169  if(all_true(abs(H[i+1][i]) < EPSILON))
1170  DUNE_THROW(SolverAbort,
1171  "breakdown in GMRes - |w| == 0.0 after " << j << " iterations");
1172 
1173  // normalize new vector
1174  v[i+1] = w; v[i+1] *= real_type(1.0)/H[i+1][i];
1175 
1176  // update QR factorization
1177  for(int k=0; k<i; k++)
1178  applyPlaneRotation(H[k][i],H[k+1][i],cs[k],sn[k]);
1179 
1180  // compute new givens rotation
1181  generatePlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1182  // finish updating QR factorization
1183  applyPlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1184  applyPlaneRotation(s[i],s[i+1],cs[i],sn[i]);
1185 
1186  // norm of the defect is the last component the vector s
1187  norm = abs(s[i+1]);
1188 
1189  // print current iteration statistics
1190  if(_verbose > 1) {
1191  this->printOutput(std::cout,j,norm,norm_old);
1192  }
1193 
1194  norm_old = norm;
1195 
1196  // check convergence
1197  if(all_true(norm < real_type(reduction) * norm_0))
1198  res.converged = true;
1199 
1200  } // end for
1201 
1202  // calculate update vector
1203  w = 0.0;
1204  update(w,i,H,s,v);
1205  // and current iterate
1206  x += w;
1207 
1208  // restart GMRes if convergence was not achieved,
1209  // i.e. linear defect has not reached desired reduction
1210  // and if j < _maxit
1211  if( res.converged != true && j <= _maxit ) {
1212 
1213  if(_verbose > 0)
1214  std::cout << "=== GMRes::restart" << std::endl;
1215  // get saved rhs
1216  b = b2;
1217  // calculate new defect
1218  _op->applyscaleadd(-1.0,x,b); // b -= Ax;
1219  // calculate preconditioned defect
1220  v[0] = 0.0;
1221  _prec->apply(v[0],b);
1222  norm = _sp->norm(v[0]);
1223  norm_old = norm;
1224  }
1225 
1226  } //end while
1227 
1228  // postprocess preconditioner
1229  _prec->post(x);
1230 
1231  // save solver statistics
1232  res.iterations = j-1; // it has to be j-1!!!
1233  res.reduction = static_cast<double>(max_value(norm/norm_0));
1234  res.conv_rate = pow(res.reduction,1.0/(j-1));
1235  res.elapsed = watch.elapsed();
1236 
1237  if(_verbose>0)
1238  print_result(res);
1239 
1240  }
1241 
1242  private :
1243 
1244  void print_result(const InverseOperatorResult& res) const {
1245  int k = res.iterations>0 ? res.iterations : 1;
1246  std::cout << "=== rate=" << res.conv_rate
1247  << ", T=" << res.elapsed
1248  << ", TIT=" << res.elapsed/k
1249  << ", IT=" << res.iterations
1250  << std::endl;
1251  }
1252 
1253  void update(X& w, int i,
1254  const std::vector<std::vector<field_type,fAlloc> >& H,
1255  const std::vector<field_type,fAlloc>& s,
1256  const std::vector<X>& v) {
1257  // solution vector of the upper triangular system
1258  std::vector<field_type,fAlloc> y(s);
1259 
1260  // backsolve
1261  for(int a=i-1; a>=0; a--) {
1262  field_type rhs(s[a]);
1263  for(int b=a+1; b<i; b++)
1264  rhs -= H[a][b]*y[b];
1265  y[a] = rhs/H[a][a];
1266 
1267  // compute update on the fly
1268  // w += y[a]*v[a]
1269  w.axpy(y[a],v[a]);
1270  }
1271  }
1272 
1273  template<typename T>
1274  typename std::enable_if<std::is_same<field_type,real_type>::value,T>::type conjugate(const T& t) {
1275  return t;
1276  }
1277 
1278  template<typename T>
1279  typename std::enable_if<!std::is_same<field_type,real_type>::value,T>::type conjugate(const T& t) {
1280  using std::conj;
1281  return conj(t);
1282  }
1283 
1284  // helper function to extract the real value of a real or complex number
1285  inline
1286  real_type to_real(const real_type & v)
1287  {
1288  return v;
1289  }
1290 
1291  inline
1292  real_type to_real(const std::complex<real_type> & v)
1293  {
1294  return v.real();
1295  }
1296 
1297  void
1298  generatePlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1299  {
1300  using std::sqrt;
1301  using std::abs;
1302  using std::max;
1303  using std::min;
1304  const real_type eps = 1e-15;
1305  real_type norm_dx = abs(dx);
1306  real_type norm_dy = abs(dy);
1307  real_type norm_max = max(norm_dx, norm_dy);
1308  real_type norm_min = min(norm_dx, norm_dy);
1309  real_type temp = norm_min/norm_max;
1310  // we rewrite the code in a vectorizable fashion
1311  cs = cond(norm_dy < eps,
1312  real_type(1.0),
1313  cond(norm_dx < eps,
1314  real_type(0.0),
1315  cond(norm_dy > norm_dx,
1316  real_type(1.0)/sqrt(real_type(1.0) + temp*temp)*temp,
1317  real_type(1.0)/sqrt(real_type(1.0) + temp*temp)
1318  )));
1319  sn = cond(norm_dy < eps,
1320  field_type(0.0),
1321  cond(norm_dx < eps,
1322  field_type(1.0),
1323  cond(norm_dy > norm_dx,
1324  field_type(1.0)/sqrt(real_type(1.0) + temp*temp)*dx*conjugate(dy)/norm_dx/norm_dy,
1325  field_type(1.0)/sqrt(real_type(1.0) + temp*temp)*conjugate(dy/dx)
1326  )));
1327  }
1328 
1329 
1330  void
1331  applyPlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1332  {
1333  field_type temp = cs * dx + sn * dy;
1334  dy = -conjugate(sn) * dx + cs * dy;
1335  dx = temp;
1336  }
1337 
1344  int _restart;
1345  };
1346 
1347 
1361  template<class X>
1363  {
1364  public:
1365  using typename IterativeSolver<X,X>::domain_type;
1366  using typename IterativeSolver<X,X>::range_type;
1367  using typename IterativeSolver<X,X>::field_type;
1368  using typename IterativeSolver<X,X>::real_type;
1369 
1370  private:
1372 
1374  using fAlloc = ReboundAllocatorType<X,field_type>;
1375 
1376  public:
1377 
1378  // don't shadow four-argument version of apply defined in the base class
1380 
1387  GeneralizedPCGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
1388  IterativeSolver<X,X>::IterativeSolver(op,prec,reduction,maxit,verbose),
1389  _restart(restart)
1390  {}
1391 
1399  GeneralizedPCGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
1400  IterativeSolver<X,X>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
1401  _restart(restart)
1402  {}
1403 
1409  virtual void apply (X& x, X& b, InverseOperatorResult& res)
1410  {
1411  res.clear(); // clear solver statistics
1412  Timer watch; // start a timer
1413  _prec->pre(x,b); // prepare preconditioner
1414  _op->applyscaleadd(-1,x,b); // overwrite b with defect
1415 
1416  std::vector<std::shared_ptr<X> > p(_restart);
1417  std::vector<field_type,fAlloc> pp(_restart);
1418  X q(x); // a temporary vector
1419  X prec_res(x); // a temporary vector for preconditioner output
1420 
1421  p[0].reset(new X(x));
1422 
1423  real_type def0 = _sp->norm(b); // compute norm
1424  if ( max_value(def0) < 1E-30 ) // convergence check
1425  {
1426  res.converged = true;
1427  res.iterations = 0; // fill statistics
1428  res.reduction = 0;
1429  res.conv_rate = 0;
1430  res.elapsed=0;
1431  if (_verbose>0) // final print
1432  std::cout << "=== rate=" << res.conv_rate
1433  << ", T=" << res.elapsed << ", TIT=" << res.elapsed
1434  << ", IT=0" << std::endl;
1435  return;
1436  }
1437 
1438  if (_verbose>0) // printing
1439  {
1440  std::cout << "=== GeneralizedPCGSolver" << std::endl;
1441  if (_verbose>1) {
1442  this->printHeader(std::cout);
1443  this->printOutput(std::cout,0,def0);
1444  }
1445  }
1446  // some local variables
1447  real_type def=def0; // loop variables
1448  field_type rho, lambda;
1449 
1450  int i=0;
1451  int ii=0;
1452  // determine initial search direction
1453  *(p[0]) = 0; // clear correction
1454  _prec->apply(*(p[0]),b); // apply preconditioner
1455  rho = _sp->dot(*(p[0]),b); // orthogonalization
1456  _op->apply(*(p[0]),q); // q=Ap
1457  pp[0] = _sp->dot(*(p[0]),q); // scalar product
1458  lambda = rho/pp[0]; // minimization
1459  x.axpy(lambda,*(p[0])); // update solution
1460  b.axpy(-lambda,q); // update defect
1461 
1462  // convergence test
1463  real_type defnew=_sp->norm(b); // comp defect norm
1464  ++i;
1465  if (_verbose>1) // print
1466  this->printOutput(std::cout,i,defnew,def);
1467  def = defnew; // update norm
1468  if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
1469  {
1470  res.converged = true;
1471  if (_verbose>0) // final print
1472  {
1473  std::cout << "=== rate=" << res.conv_rate
1474  << ", T=" << res.elapsed
1475  << ", TIT=" << res.elapsed
1476  << ", IT=" << 1 << std::endl;
1477  }
1478  return;
1479  }
1480 
1481  while(i<_maxit) {
1482  // the loop
1483  int end=std::min(_restart, _maxit-i+1);
1484  for (ii=1; ii<end; ++ii )
1485  {
1486  //std::cout<<" ii="<<ii<<" i="<<i<<std::endl;
1487  // compute next conjugate direction
1488  prec_res = 0; // clear correction
1489  _prec->apply(prec_res,b); // apply preconditioner
1490 
1491  p[ii].reset(new X(prec_res));
1492  _op->apply(prec_res, q);
1493 
1494  for(int j=0; j<ii; ++j) {
1495  rho =_sp->dot(q,*(p[j]))/pp[j];
1496  p[ii]->axpy(-rho, *(p[j]));
1497  }
1498 
1499  // minimize in given search direction
1500  _op->apply(*(p[ii]),q); // q=Ap
1501  pp[ii] = _sp->dot(*(p[ii]),q); // scalar product
1502  rho = _sp->dot(*(p[ii]),b); // orthogonalization
1503  lambda = rho/pp[ii]; // minimization
1504  x.axpy(lambda,*(p[ii])); // update solution
1505  b.axpy(-lambda,q); // update defect
1506 
1507  // convergence test
1508  defnew=_sp->norm(b); // comp defect norm
1509 
1510  ++i;
1511  if (_verbose>1) // print
1512  this->printOutput(std::cout,i,defnew,def);
1513 
1514  def = defnew; // update norm
1515  if (all_true(def<def0*_reduction) || max_value(def)<1E-30) // convergence check
1516  {
1517  res.converged = true;
1518  break;
1519  }
1520  }
1521  if(res.converged)
1522  break;
1523  if(end==_restart) {
1524  *(p[0])=*(p[_restart-1]);
1525  pp[0]=pp[_restart-1];
1526  }
1527  }
1528 
1529  // postprocess preconditioner
1530  _prec->post(x);
1531 
1532  // fill statistics
1533  res.iterations = i;
1534  res.reduction = static_cast<double>(max_value(def/def0));
1535  res.conv_rate = pow(res.reduction,1.0/i);
1536  res.elapsed = watch.elapsed();
1537 
1538  if (_verbose>0) // final print
1539  {
1540  std::cout << "=== rate=" << res.conv_rate
1541  << ", T=" << res.elapsed
1542  << ", TIT=" << res.elapsed/i
1543  << ", IT=" << i+1 << std::endl;
1544  }
1545  }
1546 
1547  private:
1554  int _restart;
1555  };
1556 
1559 } // end namespace
1560 
1561 #endif
Base class for all implementations of iterative solvers.
Definition: solver.hh:195
GeneralizedPCGSolver(LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1387
X::field_type field_type
The field type of the operator.
Definition: solver.hh:100
Definition: allocator.hh:7
double reduction
Reduction achieved: .
Definition: solver.hh:62
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:422
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:548
typename AllocatorTraits< T >::type::template rebind< X >::other ReboundAllocatorType
Definition: allocator.hh:33
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:780
gradient method
Definition: solvers.hh:160
Define general, extensible interface for inverse operators.
virtual void apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1082
SimdScalar< real_type > scalar_real_type
scalar type underlying the field_type
Definition: solver.hh:106
FieldTraits< field_type >::real_type real_type
The real type of the field type (is the same if using real numbers, but differs for std::complex) ...
Definition: solver.hh:103
std::shared_ptr< ScalarProduct< X > > _sp
Definition: solver.hh:302
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1056
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1894
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:155
A linear operator.
Definition: operators.hh:64
Generalized preconditioned conjugate gradient solver.
Definition: solvers.hh:1362
Define base class for scalar product and norm.
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:528
RestartedGMResSolver(LinearOperator< X, Y > &op, Preconditioner< X, Y > &prec, scalar_real_type reduction, int restart, int maxit, int verbose)
Set up RestartedGMResSolver solver.
Definition: solvers.hh:1058
std::shared_ptr< LinearOperator< X, X > > _op
Definition: solver.hh:300
scalar_real_type _reduction
Definition: solver.hh:303
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:178
Dune::Std::bool_constant<(std::is_same< field_type, float >::value||std::is_same< field_type, double >::value)> enableConditionEstimate_t
Definition: solvers.hh:269
double elapsed
Elapsed time in seconds.
Definition: solver.hh:74
Implementation of the BCRSMatrix class.
Preconditioned loop solver.
Definition: solvers.hh:56
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1409
Define general, extensible interface for operators. The available implementation wraps a matrix...
Thrown when a solver aborts due to some problem.
Definition: istlexception.hh:43
Wrapper to use a range of ARPACK++ eigenvalue solvers.
Definition: arpackpp.hh:241
bool converged
True if convergence criterion has been met.
Definition: solver.hh:65
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:1034
conjugate gradient method
Definition: solvers.hh:254
CGSolver(LinearOperator< X, X > &op, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
Constructor to initialize a CG solver.
Definition: solvers.hh:283
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:70
void computeSymMaxMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition: arpackpp.hh:286
int _maxit
Definition: solver.hh:304
Minimal Residual Method (MINRES)
Definition: solvers.hh:762
CGSolver(LinearOperator< X, X > &op, ScalarProduct< X > &sp, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
Constructor to initialize a CG solver.
Definition: solvers.hh:300
X domain_type
Type of the domain of the operator to be inverted.
Definition: solver.hh:94
A vector of blocks with memory management.
Definition: bvector.hh:316
void clear()
Resets all data.
Definition: solver.hh:49
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:322
void computeSymMinMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition: arpackpp.hh:388
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:46
virtual void apply(X &x, Y &b, double reduction, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1095
Y range_type
Type of the range of the operator to be inverted.
Definition: solver.hh:97
GeneralizedPCGSolver(LinearOperator< X, X > &op, ScalarProduct< X > &sp, Preconditioner< X, X > &prec, scalar_real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1399
Statistics about the application of an inverse operator.
Definition: solver.hh:40
int _verbose
Definition: solver.hh:305
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1062
RestartedGMResSolver(LinearOperator< X, Y > &op, ScalarProduct< X > &sp, Preconditioner< X, Y > &prec, scalar_real_type reduction, int restart, int maxit, int verbose)
Set up RestartedGMResSolver solver.
Definition: solvers.hh:1069
int iterations
Number of iterations.
Definition: solver.hh:59
std::shared_ptr< Preconditioner< X, X > > _prec
Definition: solver.hh:301
void printOutput(std::ostream &s, const CountType &iter, const DataType &norm, const DataType &norm_old) const
helper function for printing solver output
Definition: solver.hh:164
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:68
double condition_estimate
Estimate of condition number.
Definition: solver.hh:71