dune-istl  2.4
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 <cmath>
8 #include <complex>
9 #include <iostream>
10 #include <iomanip>
11 #include <memory>
12 #include <string>
13 #include <vector>
14 
15 #include "istlexception.hh"
16 #include "operators.hh"
17 #include "scalarproducts.hh"
18 #include "solver.hh"
19 #include "preconditioner.hh"
20 #include <dune/common/array.hh>
21 #include <dune/common/deprecated.hh>
22 #include <dune/common/timer.hh>
23 #include <dune/common/ftraits.hh>
24 #include <dune/common/typetraits.hh>
25 
26 namespace Dune {
44  //=====================================================================
45  // Implementation of this interface
46  //=====================================================================
47 
56  template<class X>
57  class LoopSolver : public InverseOperator<X,X> {
58  public:
60  typedef X domain_type;
62  typedef X range_type;
64  typedef typename X::field_type field_type;
66  typedef typename FieldTraits<field_type>::real_type real_type;
67 
87  template<class L, class P>
88  LoopSolver (L& op, P& prec,
89  real_type reduction, int maxit, int verbose) :
90  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
91  {
92  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
93  "L and P have to have the same category!");
94  static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
95  "L has to be sequential!");
96  }
97 
118  template<class L, class S, class P>
119  LoopSolver (L& op, S& sp, P& prec,
120  real_type reduction, int maxit, int verbose) :
121  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
122  {
123  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
124  "L and P must have the same category!");
125  static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
126  "L and S must have the same category!");
127  }
128 
129 
131  virtual void apply (X& x, X& b, InverseOperatorResult& res)
132  {
133  // clear solver statistics
134  res.clear();
135 
136  // start a timer
137  Timer watch;
138 
139  // prepare preconditioner
140  _prec.pre(x,b);
141 
142  // overwrite b with defect
143  _op.applyscaleadd(-1,x,b);
144 
145  // compute norm, \todo parallelization
146  real_type def0 = _sp.norm(b);
147 
148  // printing
149  if (_verbose>0)
150  {
151  std::cout << "=== LoopSolver" << std::endl;
152  if (_verbose>1)
153  {
154  this->printHeader(std::cout);
155  this->printOutput(std::cout,real_type(0),def0);
156  }
157  }
158 
159  // allocate correction vector
160  X v(x);
161 
162  // iteration loop
163  int i=1; real_type def=def0;
164  for ( ; i<=_maxit; i++ )
165  {
166  v = 0; // clear correction
167  _prec.apply(v,b); // apply preconditioner
168  x += v; // update solution
169  _op.applyscaleadd(-1,v,b); // update defect
170  real_type defnew=_sp.norm(b); // comp defect norm
171  if (_verbose>1) // print
172  this->printOutput(std::cout,real_type(i),defnew,def);
173  //std::cout << i << " " << defnew << " " << defnew/def << std::endl;
174  def = defnew; // update norm
175  if (def<def0*_reduction || def<1E-30) // convergence check
176  {
177  res.converged = true;
178  break;
179  }
180  }
181 
182  //correct i which is wrong if convergence was not achieved.
183  i=std::min(_maxit,i);
184 
185  // print
186  if (_verbose==1)
187  this->printOutput(std::cout,real_type(i),def);
188 
189  // postprocess preconditioner
190  _prec.post(x);
191 
192  // fill statistics
193  res.iterations = i;
194  res.reduction = def/def0;
195  res.conv_rate = pow(res.reduction,1.0/i);
196  res.elapsed = watch.elapsed();
197 
198  // final print
199  if (_verbose>0)
200  {
201  std::cout << "=== rate=" << res.conv_rate
202  << ", T=" << res.elapsed
203  << ", TIT=" << res.elapsed/i
204  << ", IT=" << i << std::endl;
205  }
206  }
207 
209  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
210  {
211  real_type saved_reduction = _reduction;
212  _reduction = reduction;
213  (*this).apply(x,b,res);
214  _reduction = saved_reduction;
215  }
216 
217  private:
219  LinearOperator<X,X>& _op;
220  Preconditioner<X,X>& _prec;
221  ScalarProduct<X>& _sp;
222  real_type _reduction;
223  int _maxit;
224  int _verbose;
225  };
226 
227 
228  // all these solvers are taken from the SUMO library
230  template<class X>
231  class GradientSolver : public InverseOperator<X,X> {
232  public:
234  typedef X domain_type;
236  typedef X range_type;
238  typedef typename X::field_type field_type;
240  typedef typename FieldTraits<field_type>::real_type real_type;
241 
242 
248  template<class L, class P>
249  GradientSolver (L& op, P& prec,
250  real_type reduction, int maxit, int verbose) :
251  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
252  {
253  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
254  "L and P have to have the same category!");
255  static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
256  "L has to be sequential!");
257  }
263  template<class L, class S, class P>
264  GradientSolver (L& op, S& sp, P& prec,
265  real_type reduction, int maxit, int verbose) :
266  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
267  {
268  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
269  "L and P have to have the same category!");
270  static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
271  "L and S have to have the same category!");
272  }
273 
279  virtual void apply (X& x, X& b, InverseOperatorResult& res)
280  {
281  res.clear(); // clear solver statistics
282  Timer watch; // start a timer
283  _prec.pre(x,b); // prepare preconditioner
284  _op.applyscaleadd(-1,x,b); // overwrite b with defect
285 
286  X p(x); // create local vectors
287  X q(b);
288 
289  real_type def0 = _sp.norm(b); // compute norm
290 
291  if (_verbose>0) // printing
292  {
293  std::cout << "=== GradientSolver" << std::endl;
294  if (_verbose>1)
295  {
296  this->printHeader(std::cout);
297  this->printOutput(std::cout,real_type(0),def0);
298  }
299  }
300 
301  int i=1; real_type def=def0; // loop variables
302  field_type lambda;
303  for ( ; i<=_maxit; i++ )
304  {
305  p = 0; // clear correction
306  _prec.apply(p,b); // apply preconditioner
307  _op.apply(p,q); // q=Ap
308  lambda = _sp.dot(p,b)/_sp.dot(q,p); // minimization
309  x.axpy(lambda,p); // update solution
310  b.axpy(-lambda,q); // update defect
311 
312  real_type defnew=_sp.norm(b); // comp defect norm
313  if (_verbose>1) // print
314  this->printOutput(std::cout,real_type(i),defnew,def);
315 
316  def = defnew; // update norm
317  if (def<def0*_reduction || def<1E-30) // convergence check
318  {
319  res.converged = true;
320  break;
321  }
322  }
323 
324  //correct i which is wrong if convergence was not achieved.
325  i=std::min(_maxit,i);
326 
327  if (_verbose==1) // printing for non verbose
328  this->printOutput(std::cout,real_type(i),def);
329 
330  _prec.post(x); // postprocess preconditioner
331  res.iterations = i; // fill statistics
332  res.reduction = static_cast<double>(def/def0);
333  res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
334  res.elapsed = watch.elapsed();
335  if (_verbose>0) // final print
336  std::cout << "=== rate=" << res.conv_rate
337  << ", T=" << res.elapsed
338  << ", TIT=" << res.elapsed/i
339  << ", IT=" << i << std::endl;
340  }
341 
347  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
348  {
349  real_type saved_reduction = _reduction;
350  _reduction = reduction;
351  (*this).apply(x,b,res);
352  _reduction = saved_reduction;
353  }
354 
355  private:
357  LinearOperator<X,X>& _op;
358  Preconditioner<X,X>& _prec;
359  ScalarProduct<X>& _sp;
360  real_type _reduction;
361  int _maxit;
362  int _verbose;
363  };
364 
365 
366 
368  template<class X>
369  class CGSolver : public InverseOperator<X,X> {
370  public:
372  typedef X domain_type;
374  typedef X range_type;
376  typedef typename X::field_type field_type;
378  typedef typename FieldTraits<field_type>::real_type real_type;
379 
385  template<class L, class P>
386  CGSolver (L& op, P& prec, real_type reduction, int maxit, int verbose) :
387  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
388  {
389  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
390  "L and P must have the same category!");
391  static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
392  "L must be sequential!");
393  }
399  template<class L, class S, class P>
400  CGSolver (L& op, S& sp, P& prec, real_type reduction, int maxit, int verbose) :
401  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
402  {
403  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
404  "L and P must have the same category!");
405  static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
406  "L and S must have the same category!");
407  }
408 
414  virtual void apply (X& x, X& b, InverseOperatorResult& res)
415  {
416  res.clear(); // clear solver statistics
417  Timer watch; // start a timer
418  _prec.pre(x,b); // prepare preconditioner
419  _op.applyscaleadd(-1,x,b); // overwrite b with defect
420 
421  X p(x); // the search direction
422  X q(x); // a temporary vector
423 
424  real_type def0 = _sp.norm(b); // compute norm
425  if (def0<1E-30) // convergence check
426  {
427  res.converged = true;
428  res.iterations = 0; // fill statistics
429  res.reduction = 0;
430  res.conv_rate = 0;
431  res.elapsed=0;
432  if (_verbose>0) // final print
433  std::cout << "=== rate=" << res.conv_rate
434  << ", T=" << res.elapsed << ", TIT=" << res.elapsed
435  << ", IT=0" << std::endl;
436  return;
437  }
438 
439  if (_verbose>0) // printing
440  {
441  std::cout << "=== CGSolver" << std::endl;
442  if (_verbose>1) {
443  this->printHeader(std::cout);
444  this->printOutput(std::cout,real_type(0),def0);
445  }
446  }
447 
448  // some local variables
449  real_type def=def0; // loop variables
450  field_type rho,rholast,lambda,alpha,beta;
451 
452  // determine initial search direction
453  p = 0; // clear correction
454  _prec.apply(p,b); // apply preconditioner
455  rholast = _sp.dot(p,b); // orthogonalization
456 
457  // the loop
458  int i=1;
459  for ( ; i<=_maxit; i++ )
460  {
461  // minimize in given search direction p
462  _op.apply(p,q); // q=Ap
463  alpha = _sp.dot(p,q); // scalar product
464  lambda = rholast/alpha; // minimization
465  x.axpy(lambda,p); // update solution
466  b.axpy(-lambda,q); // update defect
467 
468  // convergence test
469  real_type defnew=_sp.norm(b); // comp defect norm
470 
471  if (_verbose>1) // print
472  this->printOutput(std::cout,real_type(i),defnew,def);
473 
474  def = defnew; // update norm
475  if (def<def0*_reduction || def<1E-30) // convergence check
476  {
477  res.converged = true;
478  break;
479  }
480 
481  // determine new search direction
482  q = 0; // clear correction
483  _prec.apply(q,b); // apply preconditioner
484  rho = _sp.dot(q,b); // orthogonalization
485  beta = rho/rholast; // scaling factor
486  p *= beta; // scale old search direction
487  p += q; // orthogonalization with correction
488  rholast = rho; // remember rho for recurrence
489  }
490 
491  //correct i which is wrong if convergence was not achieved.
492  i=std::min(_maxit,i);
493 
494  if (_verbose==1) // printing for non verbose
495  this->printOutput(std::cout,real_type(i),def);
496 
497  _prec.post(x); // postprocess preconditioner
498  res.iterations = i; // fill statistics
499  res.reduction = static_cast<double>(def/def0);
500  res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
501  res.elapsed = watch.elapsed();
502 
503  if (_verbose>0) // final print
504  {
505  std::cout << "=== rate=" << res.conv_rate
506  << ", T=" << res.elapsed
507  << ", TIT=" << res.elapsed/i
508  << ", IT=" << i << std::endl;
509  }
510  }
511 
517  virtual void apply (X& x, X& b, double reduction,
519  {
520  real_type saved_reduction = _reduction;
521  _reduction = reduction;
522  (*this).apply(x,b,res);
523  _reduction = saved_reduction;
524  }
525 
526  private:
528  LinearOperator<X,X>& _op;
529  Preconditioner<X,X>& _prec;
530  ScalarProduct<X>& _sp;
531  real_type _reduction;
532  int _maxit;
533  int _verbose;
534  };
535 
536 
537  // Ronald Kriemanns BiCG-STAB implementation from Sumo
539  template<class X>
540  class BiCGSTABSolver : public InverseOperator<X,X> {
541  public:
543  typedef X domain_type;
545  typedef X range_type;
547  typedef typename X::field_type field_type;
549  typedef typename FieldTraits<field_type>::real_type real_type;
550 
556  template<class L, class P>
557  BiCGSTABSolver (L& op, P& prec,
558  real_type reduction, int maxit, int verbose) :
559  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
560  {
561  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
562  "L and P must be of the same category!");
563  static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
564  "L must be sequential!");
565  }
571  template<class L, class S, class P>
572  BiCGSTABSolver (L& op, S& sp, P& prec,
573  real_type reduction, int maxit, int verbose) :
574  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
575  {
576  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
577  "L and P must have the same category!");
578  static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
579  "L and S must have the same category!");
580  }
581 
587  virtual void apply (X& x, X& b, InverseOperatorResult& res)
588  {
589  const real_type EPSILON=1e-80;
590  double it;
591  field_type rho, rho_new, alpha, beta, h, omega;
592  real_type norm, norm_old, norm_0;
593 
594  //
595  // get vectors and matrix
596  //
597  X& r=b;
598  X p(x);
599  X v(x);
600  X t(x);
601  X y(x);
602  X rt(x);
603 
604  //
605  // begin iteration
606  //
607 
608  // r = r - Ax; rt = r
609  res.clear(); // clear solver statistics
610  Timer watch; // start a timer
611  _prec.pre(x,r); // prepare preconditioner
612  _op.applyscaleadd(-1,x,r); // overwrite b with defect
613 
614  rt=r;
615 
616  norm = norm_old = norm_0 = _sp.norm(r);
617 
618  p=0;
619  v=0;
620 
621  rho = 1;
622  alpha = 1;
623  omega = 1;
624 
625  if (_verbose>0) // printing
626  {
627  std::cout << "=== BiCGSTABSolver" << std::endl;
628  if (_verbose>1)
629  {
630  this->printHeader(std::cout);
631  this->printOutput(std::cout,real_type(0),norm_0);
632  //std::cout << " Iter Defect Rate" << std::endl;
633  //std::cout << " 0" << std::setw(14) << norm_0 << std::endl;
634  }
635  }
636 
637  if ( norm < (_reduction * norm_0) || norm<1E-30)
638  {
639  res.converged = 1;
640  _prec.post(x); // postprocess preconditioner
641  res.iterations = 0; // fill statistics
642  res.reduction = 0;
643  res.conv_rate = 0;
644  res.elapsed = watch.elapsed();
645  return;
646  }
647 
648  //
649  // iteration
650  //
651 
652  for (it = 0.5; it < _maxit; it+=.5)
653  {
654  //
655  // preprocess, set vecsizes etc.
656  //
657 
658  // rho_new = < rt , r >
659  rho_new = _sp.dot(rt,r);
660 
661  // look if breakdown occured
662  if (std::abs(rho) <= EPSILON)
663  DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho "
664  << rho << " <= EPSILON " << EPSILON
665  << " after " << it << " iterations");
666  if (std::abs(omega) <= EPSILON)
667  DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - omega "
668  << omega << " <= EPSILON " << EPSILON
669  << " after " << it << " iterations");
670 
671 
672  if (it<1)
673  p = r;
674  else
675  {
676  beta = ( rho_new / rho ) * ( alpha / omega );
677  p.axpy(-omega,v); // p = r + beta (p - omega*v)
678  p *= beta;
679  p += r;
680  }
681 
682  // y = W^-1 * p
683  y = 0;
684  _prec.apply(y,p); // apply preconditioner
685 
686  // v = A * y
687  _op.apply(y,v);
688 
689  // alpha = rho_new / < rt, v >
690  h = _sp.dot(rt,v);
691 
692  if ( std::abs(h) < EPSILON )
693  DUNE_THROW(ISTLError,"h=0 in BiCGSTAB");
694 
695  alpha = rho_new / h;
696 
697  // apply first correction to x
698  // x <- x + alpha y
699  x.axpy(alpha,y);
700 
701  // r = r - alpha*v
702  r.axpy(-alpha,v);
703 
704  //
705  // test stop criteria
706  //
707 
708  norm = _sp.norm(r);
709 
710  if (_verbose>1) // print
711  {
712  this->printOutput(std::cout,real_type(it),norm,norm_old);
713  }
714 
715  if ( norm < (_reduction * norm_0) )
716  {
717  res.converged = 1;
718  break;
719  }
720  it+=.5;
721 
722  norm_old = norm;
723 
724  // y = W^-1 * r
725  y = 0;
726  _prec.apply(y,r);
727 
728  // t = A * y
729  _op.apply(y,t);
730 
731  // omega = < t, r > / < t, t >
732  omega = _sp.dot(t,r)/_sp.dot(t,t);
733 
734  // apply second correction to x
735  // x <- x + omega y
736  x.axpy(omega,y);
737 
738  // r = s - omega*t (remember : r = s)
739  r.axpy(-omega,t);
740 
741  rho = rho_new;
742 
743  //
744  // test stop criteria
745  //
746 
747  norm = _sp.norm(r);
748 
749  if (_verbose > 1) // print
750  {
751  this->printOutput(std::cout,real_type(it),norm,norm_old);
752  }
753 
754  if ( norm < (_reduction * norm_0) || norm<1E-30)
755  {
756  res.converged = 1;
757  break;
758  }
759 
760  norm_old = norm;
761  } // end for
762 
763  //correct i which is wrong if convergence was not achieved.
764  it=std::min(static_cast<double>(_maxit),it);
765 
766  if (_verbose==1) // printing for non verbose
767  this->printOutput(std::cout,real_type(it),norm);
768 
769  _prec.post(x); // postprocess preconditioner
770  res.iterations = static_cast<int>(std::ceil(it)); // fill statistics
771  res.reduction = static_cast<double>(norm/norm_0);
772  res.conv_rate = static_cast<double>(pow(res.reduction,1.0/it));
773  res.elapsed = watch.elapsed();
774  if (_verbose>0) // final print
775  std::cout << "=== rate=" << res.conv_rate
776  << ", T=" << res.elapsed
777  << ", TIT=" << res.elapsed/it
778  << ", IT=" << it << std::endl;
779  }
780 
786  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
787  {
788  real_type saved_reduction = _reduction;
789  _reduction = reduction;
790  (*this).apply(x,b,res);
791  _reduction = saved_reduction;
792  }
793 
794  private:
796  LinearOperator<X,X>& _op;
797  Preconditioner<X,X>& _prec;
798  ScalarProduct<X>& _sp;
799  real_type _reduction;
800  int _maxit;
801  int _verbose;
802  };
803 
810  template<class X>
811  class MINRESSolver : public InverseOperator<X,X> {
812  public:
814  typedef X domain_type;
816  typedef X range_type;
818  typedef typename X::field_type field_type;
820  typedef typename FieldTraits<field_type>::real_type real_type;
821 
827  template<class L, class P>
828  MINRESSolver (L& op, P& prec, real_type reduction, int maxit, int verbose) :
829  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
830  {
831  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
832  "L and P must have the same category!");
833  static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
834  "L must be sequential!");
835  }
841  template<class L, class S, class P>
842  MINRESSolver (L& op, S& sp, P& prec, real_type reduction, int maxit, int verbose) :
843  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
844  {
845  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
846  "L and P must have the same category!");
847  static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
848  "L and S must have the same category!");
849  }
850 
856  virtual void apply (X& x, X& b, InverseOperatorResult& res)
857  {
858  // clear solver statistics
859  res.clear();
860  // start a timer
861  Dune::Timer watch;
862  watch.reset();
863  // prepare preconditioner
864  _prec.pre(x,b);
865  // overwrite rhs with defect
866  _op.applyscaleadd(-1,x,b);
867 
868  // compute residual norm
869  real_type def0 = _sp.norm(b);
870 
871  // printing
872  if(_verbose > 0) {
873  std::cout << "=== MINRESSolver" << std::endl;
874  if(_verbose > 1) {
875  this->printHeader(std::cout);
876  this->printOutput(std::cout,real_type(0),def0);
877  }
878  }
879 
880  // check for convergence
881  if(def0 < 1e-30 ) {
882  res.converged = true;
883  res.iterations = 0;
884  res.reduction = 0;
885  res.conv_rate = 0;
886  res.elapsed = 0.0;
887  // final print
888  if(_verbose > 0)
889  std::cout << "=== rate=" << res.conv_rate
890  << ", T=" << res.elapsed
891  << ", TIT=" << res.elapsed
892  << ", IT=" << res.iterations
893  << std::endl;
894  return;
895  }
896 
897  // the defect norm
898  real_type def = def0;
899  // recurrence coefficients as computed in Lanczos algorithm
900  field_type alpha, beta;
901  // diagonal entries of givens rotation
902  Dune::array<real_type,2> c{{0.0,0.0}};
903  // off-diagonal entries of givens rotation
904  Dune::array<field_type,2> s{{0.0,0.0}};
905 
906  // recurrence coefficients (column k of tridiag matrix T_k)
907  Dune::array<field_type,3> T{{0.0,0.0,0.0}};
908 
909  // the rhs vector of the min problem
910  Dune::array<field_type,2> xi{{1.0,0.0}};
911 
912  // some temporary vectors
913  X z(b), dummy(b);
914 
915  // initialize and clear correction
916  z = 0.0;
917  _prec.apply(z,b);
918 
919  // beta is real and positive in exact arithmetic
920  // since it is the norm of the basis vectors (in unpreconditioned case)
921  beta = std::sqrt(_sp.dot(b,z));
922  field_type beta0 = beta;
923 
924  // the search directions
925  Dune::array<X,3> p{{b,b,b}};
926  p[0] = 0.0;
927  p[1] = 0.0;
928  p[2] = 0.0;
929 
930  // orthonormal basis vectors (in unpreconditioned case)
931  Dune::array<X,3> q{{b,b,b}};
932  q[0] = 0.0;
933  q[1] *= 1.0/beta;
934  q[2] = 0.0;
935 
936  z *= 1.0/beta;
937 
938  // the loop
939  int i = 1;
940  for( ; i<=_maxit; i++) {
941 
942  dummy = z;
943  int i1 = i%3,
944  i0 = (i1+2)%3,
945  i2 = (i1+1)%3;
946 
947  // symmetrically preconditioned Lanczos algorithm (see Greenbaum p.121)
948  _op.apply(z,q[i2]); // q[i2] = Az
949  q[i2].axpy(-beta,q[i0]);
950  // alpha is real since it is the diagonal entry of the hermitian tridiagonal matrix
951  // from the Lanczos Algorithm
952  // so the order in the scalar product doesn't matter even for the complex case
953  alpha = _sp.dot(z,q[i2]);
954  q[i2].axpy(-alpha,q[i1]);
955 
956  z = 0.0;
957  _prec.apply(z,q[i2]);
958 
959  // beta is real and positive in exact arithmetic
960  // since it is the norm of the basis vectors (in unpreconditioned case)
961  beta = std::sqrt(_sp.dot(q[i2],z));
962 
963  q[i2] *= 1.0/beta;
964  z *= 1.0/beta;
965 
966  // QR Factorization of recurrence coefficient matrix
967  // apply previous givens rotations to last column of T
968  T[1] = T[2];
969  if(i>2) {
970  T[0] = s[i%2]*T[1];
971  T[1] = c[i%2]*T[1];
972  }
973  if(i>1) {
974  T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
975  T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
976  }
977  else
978  T[2] = alpha;
979 
980  // update QR factorization
981  generateGivensRotation(T[2],beta,c[i%2],s[i%2]);
982  // to last column of T_k
983  T[2] = c[i%2]*T[2] + s[i%2]*beta;
984  // and to the rhs xi of the min problem
985  xi[i%2] = -s[i%2]*xi[(i+1)%2];
986  xi[(i+1)%2] *= c[i%2];
987 
988  // compute correction direction
989  p[i2] = dummy;
990  p[i2].axpy(-T[1],p[i1]);
991  p[i2].axpy(-T[0],p[i0]);
992  p[i2] *= 1.0/T[2];
993 
994  // apply correction/update solution
995  x.axpy(beta0*xi[(i+1)%2],p[i2]);
996 
997  // remember beta_old
998  T[2] = beta;
999 
1000  // check for convergence
1001  // the last entry in the rhs of the min-problem is the residual
1002  real_type defnew = std::abs(beta0*xi[i%2]);
1003 
1004  if(_verbose > 1)
1005  this->printOutput(std::cout,real_type(i),defnew,def);
1006 
1007  def = defnew;
1008  if(def < def0*_reduction || def < 1e-30 || i == _maxit ) {
1009  res.converged = true;
1010  break;
1011  }
1012  } // end for
1013 
1014  if(_verbose == 1)
1015  this->printOutput(std::cout,real_type(i),def);
1016 
1017  // postprocess preconditioner
1018  _prec.post(x);
1019  // fill statistics
1020  res.iterations = i;
1021  res.reduction = static_cast<double>(def/def0);
1022  res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
1023  res.elapsed = watch.elapsed();
1024 
1025  // final print
1026  if(_verbose > 0) {
1027  std::cout << "=== rate=" << res.conv_rate
1028  << ", T=" << res.elapsed
1029  << ", TIT=" << res.elapsed/i
1030  << ", IT=" << i << std::endl;
1031  }
1032  }
1033 
1039  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
1040  {
1041  real_type saved_reduction = _reduction;
1042  _reduction = reduction;
1043  (*this).apply(x,b,res);
1044  _reduction = saved_reduction;
1045  }
1046 
1047  private:
1048 
1049  void generateGivensRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1050  {
1051  real_type norm_dx = std::abs(dx);
1052  real_type norm_dy = std::abs(dy);
1053  if(norm_dy < 1e-15) {
1054  cs = 1.0;
1055  sn = 0.0;
1056  } else if(norm_dx < 1e-15) {
1057  cs = 0.0;
1058  sn = 1.0;
1059  } else if(norm_dy > norm_dx) {
1060  real_type temp = norm_dx/norm_dy;
1061  cs = 1.0/std::sqrt(1.0 + temp*temp);
1062  sn = cs;
1063  cs *= temp;
1064  sn *= dx/norm_dx;
1065  // dy is real in exact arithmetic
1066  // so we don't need to conjugate here
1067  sn *= dy/norm_dy;
1068  } else {
1069  real_type temp = norm_dy/norm_dx;
1070  cs = 1.0/std::sqrt(1.0 + temp*temp);
1071  sn = cs;
1072  sn *= dy/dx;
1073  // dy and dx is real in exact arithmetic
1074  // so we don't have to conjugate both of them
1075  }
1076  }
1077 
1078  SeqScalarProduct<X> ssp;
1079  LinearOperator<X,X>& _op;
1080  Preconditioner<X,X>& _prec;
1081  ScalarProduct<X>& _sp;
1082  real_type _reduction;
1083  int _maxit;
1084  int _verbose;
1085  };
1086 
1100  template<class X, class Y=X, class F = Y>
1102  {
1103  public:
1105  typedef X domain_type;
1107  typedef Y range_type;
1109  typedef typename X::field_type field_type;
1111  typedef typename FieldTraits<field_type>::real_type real_type;
1113  typedef F basis_type;
1114 
1115  template<class L, class P>
1116  DUNE_DEPRECATED_MSG("recalc_defect is a unused parameter! Use RestartedGMResSolver(L& op, P& prec, real_type reduction, int restart, int maxit, int verbose) instead")
1117  RestartedGMResSolver (L& op, P& prec, real_type reduction, int restart, int maxit, int verbose, bool recalc_defect)
1118  : _A(op)
1119  , _W(prec)
1120  , ssp()
1121  , _sp(ssp)
1122  , _restart(restart)
1123  , _reduction(reduction)
1124  , _maxit(maxit)
1125  , _verbose(verbose)
1126  {
1127  static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1128  "P and L must be the same category!");
1129  static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
1130  "L must be sequential!");
1131  }
1132 
1133 
1140  template<class L, class P>
1141  RestartedGMResSolver (L& op, P& prec, real_type reduction, int restart, int maxit, int verbose) :
1142  _A(op), _W(prec),
1143  ssp(), _sp(ssp), _restart(restart),
1144  _reduction(reduction), _maxit(maxit), _verbose(verbose)
1145  {
1146  static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1147  "P and L must be the same category!");
1148  static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
1149  "L must be sequential!");
1150  }
1151 
1152  template<class L, class S, class P>
1153  DUNE_DEPRECATED_MSG("recalc_defect is a unused parameter! Use RestartedGMResSolver(L& op, S& sp, P& prec, real_type reduction, int restart, int maxit, int verbose) instead")
1154  RestartedGMResSolver(L& op, S& sp, P& prec, real_type reduction, int restart, int maxit, int verbose, bool recalc_defect)
1155  : _A(op)
1156  , _W(prec)
1157  , _sp(sp)
1158  , _restart(restart)
1159  , _reduction(reduction)
1160  , _maxit(maxit)
1161  , _verbose(verbose)
1162  {
1163  static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1164  " P and L must have the same category!");
1165  static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1166  "P and S must have the same category!");
1167  }
1168 
1175  template<class L, class S, class P>
1176  RestartedGMResSolver (L& op, S& sp, P& prec, real_type reduction, int restart, int maxit, int verbose) :
1177  _A(op), _W(prec),
1178  _sp(sp), _restart(restart),
1179  _reduction(reduction), _maxit(maxit), _verbose(verbose)
1180  {
1181  static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1182  "P and L must have the same category!");
1183  static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1184  "P and S must have the same category!");
1185  }
1186 
1188  virtual void apply (X& x, Y& b, InverseOperatorResult& res)
1189  {
1190  apply(x,b,_reduction,res);
1191  }
1192 
1198  virtual void apply (X& x, Y& b, real_type reduction, InverseOperatorResult& res)
1199  {
1200  const real_type EPSILON = 1e-80;
1201  const int m = _restart;
1202  real_type norm, norm_old = 0.0, norm_0;
1203  int j = 1;
1204  std::vector<field_type> s(m+1), sn(m);
1205  std::vector<real_type> cs(m);
1206  // need copy of rhs if GMRes has to be restarted
1207  Y b2(b);
1208  // helper vector
1209  Y w(b);
1210  std::vector< std::vector<field_type> > H(m+1,s);
1211  std::vector<F> v(m+1,b);
1212 
1213  // start timer
1214  Dune::Timer watch;
1215  watch.reset();
1216 
1217  // clear solver statistics and set res.converged to false
1218  res.clear();
1219  _W.pre(x,b);
1220 
1221  // calculate defect and overwrite rhs with it
1222  _A.applyscaleadd(-1.0,x,b); // b -= Ax
1223  // calculate preconditioned defect
1224  v[0] = 0.0; _W.apply(v[0],b); // r = W^-1 b
1225  norm_0 = _sp.norm(v[0]);
1226  norm = norm_0;
1227  norm_old = norm;
1228 
1229  // print header
1230  if(_verbose > 0)
1231  {
1232  std::cout << "=== RestartedGMResSolver" << std::endl;
1233  if(_verbose > 1) {
1234  this->printHeader(std::cout);
1235  this->printOutput(std::cout,real_type(0),norm_0);
1236  }
1237  }
1238 
1239  if(norm_0 < EPSILON) {
1240  _W.post(x);
1241  res.converged = true;
1242  if(_verbose > 0) // final print
1243  print_result(res);
1244  }
1245 
1246  while(j <= _maxit && res.converged != true) {
1247 
1248  int i = 0;
1249  v[0] *= 1.0/norm;
1250  s[0] = norm;
1251  for(i=1; i<m+1; i++)
1252  s[i] = 0.0;
1253 
1254  for(i=0; i < m && j <= _maxit && res.converged != true; i++, j++) {
1255  w = 0.0;
1256  // use v[i+1] as temporary vector
1257  v[i+1] = 0.0;
1258  // do Arnoldi algorithm
1259  _A.apply(v[i],v[i+1]);
1260  _W.apply(w,v[i+1]);
1261  for(int k=0; k<i+1; k++) {
1262  // notice that _sp.dot(v[k],w) = v[k]\adjoint w
1263  // so one has to pay attention to the order
1264  // the in scalar product for the complex case
1265  // doing the modified Gram-Schmidt algorithm
1266  H[k][i] = _sp.dot(v[k],w);
1267  // w -= H[k][i] * v[k]
1268  w.axpy(-H[k][i],v[k]);
1269  }
1270  H[i+1][i] = _sp.norm(w);
1271  if(std::abs(H[i+1][i]) < EPSILON)
1272  DUNE_THROW(ISTLError,
1273  "breakdown in GMRes - |w| == 0.0 after " << j << " iterations");
1274 
1275  // normalize new vector
1276  v[i+1] = w; v[i+1] *= 1.0/H[i+1][i];
1277 
1278  // update QR factorization
1279  for(int k=0; k<i; k++)
1280  applyPlaneRotation(H[k][i],H[k+1][i],cs[k],sn[k]);
1281 
1282  // compute new givens rotation
1283  generatePlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1284  // finish updating QR factorization
1285  applyPlaneRotation(H[i][i],H[i+1][i],cs[i],sn[i]);
1286  applyPlaneRotation(s[i],s[i+1],cs[i],sn[i]);
1287 
1288  // norm of the defect is the last component the vector s
1289  norm = std::abs(s[i+1]);
1290 
1291  // print current iteration statistics
1292  if(_verbose > 1) {
1293  this->printOutput(std::cout,real_type(j),norm,norm_old);
1294  }
1295 
1296  norm_old = norm;
1297 
1298  // check convergence
1299  if(norm < reduction * norm_0)
1300  res.converged = true;
1301 
1302  } // end for
1303 
1304  // calculate update vector
1305  w = 0.0;
1306  update(w,i,H,s,v);
1307  // and current iterate
1308  x += w;
1309 
1310  // restart GMRes if convergence was not achieved,
1311  // i.e. linear defect has not reached desired reduction
1312  // and if j < _maxit
1313  if( res.converged != true && j <= _maxit ) {
1314 
1315  if(_verbose > 0)
1316  std::cout << "=== GMRes::restart" << std::endl;
1317  // get saved rhs
1318  b = b2;
1319  // calculate new defect
1320  _A.applyscaleadd(-1.0,x,b); // b -= Ax;
1321  // calculate preconditioned defect
1322  v[0] = 0.0;
1323  _W.apply(v[0],b);
1324  norm = _sp.norm(v[0]);
1325  norm_old = norm;
1326  }
1327 
1328  } //end while
1329 
1330  // postprocess preconditioner
1331  _W.post(x);
1332 
1333  // save solver statistics
1334  res.iterations = j-1; // it has to be j-1!!!
1335  res.reduction = static_cast<double>(norm/norm_0);
1336  res.conv_rate = static_cast<double>(pow(res.reduction,1.0/(j-1)));
1337  res.elapsed = watch.elapsed();
1338 
1339  if(_verbose>0)
1340  print_result(res);
1341 
1342  }
1343 
1344  private :
1345 
1346  void print_result(const InverseOperatorResult& res) const {
1347  int k = res.iterations>0 ? res.iterations : 1;
1348  std::cout << "=== rate=" << res.conv_rate
1349  << ", T=" << res.elapsed
1350  << ", TIT=" << res.elapsed/k
1351  << ", IT=" << res.iterations
1352  << std::endl;
1353  }
1354 
1355  void update(X& w, int i,
1356  const std::vector<std::vector<field_type> >& H,
1357  const std::vector<field_type>& s,
1358  const std::vector<X>& v) {
1359  // solution vector of the upper triangular system
1360  std::vector<field_type> y(s);
1361 
1362  // backsolve
1363  for(int a=i-1; a>=0; a--) {
1364  field_type rhs(s[a]);
1365  for(int b=a+1; b<i; b++)
1366  rhs -= H[a][b]*y[b];
1367  y[a] = rhs/H[a][a];
1368 
1369  // compute update on the fly
1370  // w += y[a]*v[a]
1371  w.axpy(y[a],v[a]);
1372  }
1373  }
1374 
1375  template<typename T>
1376  typename enable_if<is_same<field_type,real_type>::value,T>::type conjugate(const T& t) {
1377  return t;
1378  }
1379 
1380  template<typename T>
1381  typename enable_if<!is_same<field_type,real_type>::value,T>::type conjugate(const T& t) {
1382  return conj(t);
1383  }
1384 
1385  void
1386  generatePlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1387  {
1388  real_type norm_dx = std::abs(dx);
1389  real_type norm_dy = std::abs(dy);
1390  if(norm_dy < 1e-15) {
1391  cs = 1.0;
1392  sn = 0.0;
1393  } else if(norm_dx < 1e-15) {
1394  cs = 0.0;
1395  sn = 1.0;
1396  } else if(norm_dy > norm_dx) {
1397  real_type temp = norm_dx/norm_dy;
1398  cs = 1.0/std::sqrt(1.0 + temp*temp);
1399  sn = cs;
1400  cs *= temp;
1401  sn *= dx/norm_dx;
1402  sn *= conjugate(dy)/norm_dy;
1403  } else {
1404  real_type temp = norm_dy/norm_dx;
1405  cs = 1.0/std::sqrt(1.0 + temp*temp);
1406  sn = cs;
1407  sn *= conjugate(dy/dx);
1408  }
1409  }
1410 
1411 
1412  void
1413  applyPlaneRotation(field_type &dx, field_type &dy, real_type &cs, field_type &sn)
1414  {
1415  field_type temp = cs * dx + sn * dy;
1416  dy = -conjugate(sn) * dx + cs * dy;
1417  dx = temp;
1418  }
1419 
1420  LinearOperator<X,Y>& _A;
1421  Preconditioner<X,Y>& _W;
1422  SeqScalarProduct<X> ssp;
1423  ScalarProduct<X>& _sp;
1424  int _restart;
1425  real_type _reduction;
1426  int _maxit;
1427  int _verbose;
1428  };
1429 
1430 
1444  template<class X>
1446  {
1447  public:
1449  typedef X domain_type;
1451  typedef X range_type;
1453  typedef typename X::field_type field_type;
1455  typedef typename FieldTraits<field_type>::real_type real_type;
1456 
1464  template<class L, class P>
1465  GeneralizedPCGSolver (L& op, P& prec, real_type reduction, int maxit, int verbose,
1466  int restart=10) :
1467  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit),
1468  _verbose(verbose), _restart(std::min(maxit,restart))
1469  {
1470  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
1471  "L and P have to have the same category!");
1472  static_assert(static_cast<int>(L::category) ==
1473  static_cast<int>(SolverCategory::sequential),
1474  "L has to be sequential!");
1475  }
1483  template<class L, class P, class S>
1484  GeneralizedPCGSolver (L& op, S& sp, P& prec,
1485  real_type reduction, int maxit, int verbose, int restart=10) :
1486  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose),
1487  _restart(std::min(maxit,restart))
1488  {
1489  static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
1490  "L and P must have the same category!");
1491  static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
1492  "L and S must have the same category!");
1493  }
1499  virtual void apply (X& x, X& b, InverseOperatorResult& res)
1500  {
1501  res.clear(); // clear solver statistics
1502  Timer watch; // start a timer
1503  _prec.pre(x,b); // prepare preconditioner
1504  _op.applyscaleadd(-1,x,b); // overwrite b with defect
1505 
1506  std::vector<std::shared_ptr<X> > p(_restart);
1507  std::vector<typename X::field_type> pp(_restart);
1508  X q(x); // a temporary vector
1509  X prec_res(x); // a temporary vector for preconditioner output
1510 
1511  p[0].reset(new X(x));
1512 
1513  real_type def0 = _sp.norm(b); // compute norm
1514  if (def0<1E-30) // convergence check
1515  {
1516  res.converged = true;
1517  res.iterations = 0; // fill statistics
1518  res.reduction = 0;
1519  res.conv_rate = 0;
1520  res.elapsed=0;
1521  if (_verbose>0) // final print
1522  std::cout << "=== rate=" << res.conv_rate
1523  << ", T=" << res.elapsed << ", TIT=" << res.elapsed
1524  << ", IT=0" << std::endl;
1525  return;
1526  }
1527 
1528  if (_verbose>0) // printing
1529  {
1530  std::cout << "=== GeneralizedPCGSolver" << std::endl;
1531  if (_verbose>1) {
1532  this->printHeader(std::cout);
1533  this->printOutput(std::cout,real_type(0),def0);
1534  }
1535  }
1536  // some local variables
1537  real_type def=def0; // loop variables
1538  field_type rho, lambda;
1539 
1540  int i=0;
1541  int ii=0;
1542  // determine initial search direction
1543  *(p[0]) = 0; // clear correction
1544  _prec.apply(*(p[0]),b); // apply preconditioner
1545  rho = _sp.dot(*(p[0]),b); // orthogonalization
1546  _op.apply(*(p[0]),q); // q=Ap
1547  pp[0] = _sp.dot(*(p[0]),q); // scalar product
1548  lambda = rho/pp[0]; // minimization
1549  x.axpy(lambda,*(p[0])); // update solution
1550  b.axpy(-lambda,q); // update defect
1551 
1552  // convergence test
1553  real_type defnew=_sp.norm(b); // comp defect norm
1554  if (_verbose>1) // print
1555  this->printOutput(std::cout,real_type(++i),defnew,def);
1556  def = defnew; // update norm
1557  if (def<def0*_reduction || def<1E-30) // convergence check
1558  {
1559  res.converged = true;
1560  if (_verbose>0) // final print
1561  {
1562  std::cout << "=== rate=" << res.conv_rate
1563  << ", T=" << res.elapsed
1564  << ", TIT=" << res.elapsed
1565  << ", IT=" << 1 << std::endl;
1566  }
1567  return;
1568  }
1569 
1570  while(i<_maxit) {
1571  // the loop
1572  int end=std::min(_restart, _maxit-i+1);
1573  for (ii=1; ii<end; ++ii )
1574  {
1575  //std::cout<<" ii="<<ii<<" i="<<i<<std::endl;
1576  // compute next conjugate direction
1577  prec_res = 0; // clear correction
1578  _prec.apply(prec_res,b); // apply preconditioner
1579 
1580  p[ii].reset(new X(prec_res));
1581  _op.apply(prec_res, q);
1582 
1583  for(int j=0; j<ii; ++j) {
1584  rho =_sp.dot(q,*(p[j]))/pp[j];
1585  p[ii]->axpy(-rho, *(p[j]));
1586  }
1587 
1588  // minimize in given search direction
1589  _op.apply(*(p[ii]),q); // q=Ap
1590  pp[ii] = _sp.dot(*(p[ii]),q); // scalar product
1591  rho = _sp.dot(*(p[ii]),b); // orthogonalization
1592  lambda = rho/pp[ii]; // minimization
1593  x.axpy(lambda,*(p[ii])); // update solution
1594  b.axpy(-lambda,q); // update defect
1595 
1596  // convergence test
1597  real_type defnew=_sp.norm(b); // comp defect norm
1598 
1599  if (_verbose>1) // print
1600  this->printOutput(std::cout,real_type(++i),defnew,def);
1601 
1602  def = defnew; // update norm
1603  if (def<def0*_reduction || def<1E-30) // convergence check
1604  {
1605  res.converged = true;
1606  break;
1607  }
1608  }
1609  if(res.converged)
1610  break;
1611  if(end==_restart) {
1612  *(p[0])=*(p[_restart-1]);
1613  pp[0]=pp[_restart-1];
1614  }
1615  }
1616 
1617  // postprocess preconditioner
1618  _prec.post(x);
1619 
1620  // fill statistics
1621  res.iterations = i;
1622  res.reduction = def/def0;
1623  res.conv_rate = pow(res.reduction,1.0/i);
1624  res.elapsed = watch.elapsed();
1625 
1626  if (_verbose>0) // final print
1627  {
1628  std::cout << "=== rate=" << res.conv_rate
1629  << ", T=" << res.elapsed
1630  << ", TIT=" << res.elapsed/i
1631  << ", IT=" << i+1 << std::endl;
1632  }
1633  }
1634 
1640  virtual void apply (X& x, X& b, double reduction,
1641  InverseOperatorResult& res)
1642  {
1643  real_type saved_reduction = _reduction;
1644  _reduction = reduction;
1645  (*this).apply(x,b,res);
1646  _reduction = saved_reduction;
1647  }
1648  private:
1649  SeqScalarProduct<X> ssp;
1650  LinearOperator<X,X>& _op;
1651  Preconditioner<X,X>& _prec;
1652  ScalarProduct<X>& _sp;
1653  real_type _reduction;
1654  int _maxit;
1655  int _verbose;
1656  int _restart;
1657  };
1658 
1661 } // end namespace
1662 
1663 #endif
CGSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: solvers.hh:400
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: solvers.hh:549
STL namespace.
Y range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1107
virtual void apply(X &x, Y &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:1188
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:43
Preconditioned loop solver.
Definition: solvers.hh:57
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1109
Definition: basearray.hh:19
MINRESSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:842
RestartedGMResSolver(L &op, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1141
double conv_rate
Convergence rate (average reduction per step)
Definition: solver.hh:59
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: solvers.hh:378
GeneralizedPCGSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1484
Define base class for scalar product and norm.
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:1451
derive error class from the base class in common
Definition: istlexception.hh:16
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1499
RestartedGMResSolver(L &op, S &sp, P &prec, real_type reduction, int restart, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:1176
implements the Generalized Minimal Residual (GMRes) method
Definition: solvers.hh:1101
gradient method
Definition: solvers.hh:231
F basis_type
The field type of the basis vectors.
Definition: solvers.hh:1113
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1039
LoopSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up loop solver.
Definition: solvers.hh:119
BiCGSTABSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:557
double elapsed
Elapsed time in seconds.
Definition: solver.hh:62
int iterations
Number of iterations.
Definition: solver.hh:50
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:1449
virtual void post(X &x)=0
Clean up.
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:818
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:543
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:587
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:517
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:374
X::field_type field_type
The field type of the operator that we do the inverse for.
Definition: solvers.hh:238
Define general, extensible interface for operators. The available implementation wraps a matrix...
Category for sequential solvers.
Definition: solvercategory.hh:21
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:1640
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: solvers.hh:131
Minimal Residual Method (MINRES)
Definition: solvers.hh:811
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const =0
apply operator to x, scale and add:
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:540
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: solvers.hh:240
Define general, extensible interface for inverse operators.
X::field_type field_type
The field type of the operator that we do the inverse for.
Definition: solvers.hh:64
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: solvers.hh:1455
MINRESSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up MINRES solver.
Definition: solvers.hh:828
virtual void pre(X &x, Y &b)=0
Prepare the preconditioner.
double reduction
Reduction achieved: .
Definition: solver.hh:53
X domain_type
The domain type of the operator that we do the inverse for.
Definition: solvers.hh:60
void printHeader(std::ostream &s) const
helper function for printing header of solver output
Definition: solver.hh:121
Abstract base class for all solvers.
Definition: solver.hh:79
GradientSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:249
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:814
LoopSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up Loop solver.
Definition: solvers.hh:88
virtual void apply(const X &x, Y &y) const =0
apply operator to x: The input vector is consistent and the output must also be consistent on the in...
GradientSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:264
conjugate gradient method
Definition: solvers.hh:369
virtual void apply(X &x, Y &b, real_type reduction, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:1198
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:376
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: solvers.hh:820
X range_type
The range type of the operator that we do the inverse for.
Definition: solvers.hh:62
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:1105
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:279
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: solvers.hh:1111
X domain_type
The domain type of the operator to be inverted.
Definition: solvers.hh:372
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:816
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:347
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: solvers.hh:786
void printOutput(std::ostream &s, const DataType &iter, const DataType &norm, const DataType &norm_old) const
helper function for printing solver output
Definition: solver.hh:130
Definition: example.cc:34
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:856
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: solvers.hh:414
X range_type
The range type of the operator to be inverted.
Definition: solvers.hh:545
X domain_type
The domain type of the operator that we do the inverse for.
Definition: solvers.hh:234
X range_type
The range type of the operator that we do the inverse for.
Definition: solvers.hh:236
GeneralizedPCGSolver(L &op, P &prec, real_type reduction, int maxit, int verbose, int restart=10)
Set up nonlinear preconditioned conjugate gradient solver.
Definition: solvers.hh:1465
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:1453
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
Default implementation for the scalar case.
Definition: scalarproducts.hh:94
Generalized preconditioned conjugate gradient solver.
Definition: solvers.hh:1445
void clear()
Resets all data.
Definition: solver.hh:40
Statistics about the application of an inverse operator.
Definition: solver.hh:31
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: solvers.hh:66
BiCGSTABSolver(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up solver.
Definition: solvers.hh:572
CGSolver(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: solvers.hh:386
X::field_type field_type
The field type of the operator to be inverted.
Definition: solvers.hh:547
virtual void apply(X &v, const Y &d)=0
Apply one step of the preconditioner to the system A(v)=d.
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: solvers.hh:209