dune-istl  2.2.0
solvers.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set ts=4 sw=2 et sts=2:
3 
4 #ifndef DUNE_SOLVERS_HH
5 #define DUNE_SOLVERS_HH
6 
7 #include<cmath>
8 #include<complex>
9 #include<iostream>
10 #include<iomanip>
11 #include<string>
12 
13 #include "istlexception.hh"
14 #include "operators.hh"
15 #include "preconditioners.hh"
16 #include "scalarproducts.hh"
17 #include <dune/common/timer.hh>
18 #include <dune/common/ftraits.hh>
19 #include <dune/common/static_assert.hh>
20 
21 namespace Dune {
47  {
50  {
51  clear();
52  }
53 
55  void clear ()
56  {
57  iterations = 0;
58  reduction = 0;
59  converged = false;
60  conv_rate = 1;
61  elapsed = 0;
62  }
63 
66 
68  double reduction;
69 
71  bool converged;
72 
74  double conv_rate;
75 
77  double elapsed;
78  };
79 
80 
81  //=====================================================================
93  template<class X, class Y>
95  public:
97  typedef X domain_type;
98 
100  typedef Y range_type;
101 
103  typedef typename X::field_type field_type;
104 
114  virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
115 
126  virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
127 
129  virtual ~InverseOperator () {}
130 
131  protected:
132  // spacing values
133  enum { iterationSpacing = 5 , normSpacing = 16 };
134 
136  void printHeader(std::ostream& s) const
137  {
138  s << std::setw(iterationSpacing) << " Iter";
139  s << std::setw(normSpacing) << "Defect";
140  s << std::setw(normSpacing) << "Rate" << std::endl;
141  }
142 
144  template <class DataType>
145  void printOutput(std::ostream& s,
146  const double iter,
147  const DataType& norm,
148  const DataType& norm_old) const
149  {
150  const DataType rate = norm/norm_old;
151  s << std::setw(iterationSpacing) << iter << " ";
152  s << std::setw(normSpacing) << norm << " ";
153  s << std::setw(normSpacing) << rate << std::endl;
154  }
155 
157  template <class DataType>
158  void printOutput(std::ostream& s,
159  const double iter,
160  const DataType& norm) const
161  {
162  s << std::setw(iterationSpacing) << iter << " ";
163  s << std::setw(normSpacing) << norm << std::endl;
164  }
165  };
166 
167 
168  //=====================================================================
169  // Implementation of this interface
170  //=====================================================================
171 
180  template<class X>
181  class LoopSolver : public InverseOperator<X,X> {
182  public:
184  typedef X domain_type;
186  typedef X range_type;
188  typedef typename X::field_type field_type;
189 
209  template<class L, class P>
210  LoopSolver (L& op, P& prec,
211  double reduction, int maxit, int verbose) :
212  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
213  {
214  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
215  "L and P have to have the same category!");
216  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
217  "L has to be sequential!");
218  }
219 
240  template<class L, class S, class P>
241  LoopSolver (L& op, S& sp, P& prec,
242  double reduction, int maxit, int verbose) :
243  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
244  {
245  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
246  "L and P must have the same category!");
247  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
248  "L and S must have the same category!");
249  }
250 
251 
253  virtual void apply (X& x, X& b, InverseOperatorResult& res)
254  {
255  // clear solver statistics
256  res.clear();
257 
258  // start a timer
259  Timer watch;
260 
261  // prepare preconditioner
262  _prec.pre(x,b);
263 
264  // overwrite b with defect
265  _op.applyscaleadd(-1,x,b);
266 
267  // compute norm, \todo parallelization
268  double def0 = _sp.norm(b);
269 
270  // printing
271  if (_verbose>0)
272  {
273  std::cout << "=== LoopSolver" << std::endl;
274  if (_verbose>1)
275  {
276  this->printHeader(std::cout);
277  this->printOutput(std::cout,0,def0);
278  }
279  }
280 
281  // allocate correction vector
282  X v(x);
283 
284  // iteration loop
285  int i=1; double def=def0;
286  for ( ; i<=_maxit; i++ )
287  {
288  v = 0; // clear correction
289  _prec.apply(v,b); // apply preconditioner
290  x += v; // update solution
291  _op.applyscaleadd(-1,v,b); // update defect
292  double defnew=_sp.norm(b); // comp defect norm
293  if (_verbose>1) // print
294  this->printOutput(std::cout,i,defnew,def);
295  //std::cout << i << " " << defnew << " " << defnew/def << std::endl;
296  def = defnew; // update norm
297  if (def<def0*_reduction || def<1E-30) // convergence check
298  {
299  res.converged = true;
300  break;
301  }
302  }
303 
304  // print
305  if (_verbose==1)
306  this->printOutput(std::cout,i,def);
307 
308  // postprocess preconditioner
309  _prec.post(x);
310 
311  // fill statistics
312  res.iterations = i;
313  res.reduction = def/def0;
314  res.conv_rate = pow(res.reduction,1.0/i);
315  res.elapsed = watch.elapsed();
316 
317  // final print
318  if (_verbose>0)
319  {
320  std::cout << "=== rate=" << res.conv_rate
321  << ", T=" << res.elapsed
322  << ", TIT=" << res.elapsed/i
323  << ", IT=" << i << std::endl;
324  }
325  }
326 
328  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
329  {
330  std::swap(_reduction,reduction);
331  (*this).apply(x,b,res);
332  std::swap(_reduction,reduction);
333  }
334 
335  private:
337  LinearOperator<X,X>& _op;
338  Preconditioner<X,X>& _prec;
339  ScalarProduct<X>& _sp;
340  double _reduction;
341  int _maxit;
342  int _verbose;
343  };
344 
345 
346  // all these solvers are taken from the SUMO library
348  template<class X>
349  class GradientSolver : public InverseOperator<X,X> {
350  public:
352  typedef X domain_type;
354  typedef X range_type;
356  typedef typename X::field_type field_type;
357 
363  template<class L, class P>
364  GradientSolver (L& op, P& prec,
365  double reduction, int maxit, int verbose) :
366  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
367  {
368  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
369  "L and P have to have the same category!");
370  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
371  "L has to be sequential!");
372  }
378  template<class L, class S, class P>
379  GradientSolver (L& op, S& sp, P& prec,
380  double reduction, int maxit, int verbose) :
381  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
382  {
383  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
384  "L and P have to have the same category!");
385  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
386  "L and S have to have the same category!");
387  }
388 
394  virtual void apply (X& x, X& b, InverseOperatorResult& res)
395  {
396  res.clear(); // clear solver statistics
397  Timer watch; // start a timer
398  _prec.pre(x,b); // prepare preconditioner
399  _op.applyscaleadd(-1,x,b); // overwrite b with defect
400 
401  X p(x); // create local vectors
402  X q(b);
403 
404  double def0 = _sp.norm(b);// compute norm
405 
406  if (_verbose>0) // printing
407  {
408  std::cout << "=== GradientSolver" << std::endl;
409  if (_verbose>1)
410  {
411  this->printHeader(std::cout);
412  this->printOutput(std::cout,0,def0);
413  }
414  }
415 
416  int i=1; double def=def0; // loop variables
417  field_type lambda;
418  for ( ; i<=_maxit; i++ )
419  {
420  p = 0; // clear correction
421  _prec.apply(p,b); // apply preconditioner
422  _op.apply(p,q); // q=Ap
423  lambda = _sp.dot(p,b)/_sp.dot(q,p);// minimization
424  x.axpy(lambda,p); // update solution
425  b.axpy(-lambda,q); // update defect
426 
427  double defnew=_sp.norm(b);// comp defect norm
428  if (_verbose>1) // print
429  this->printOutput(std::cout,i,defnew,def);
430 
431  def = defnew; // update norm
432  if (def<def0*_reduction || def<1E-30) // convergence check
433  {
434  res.converged = true;
435  break;
436  }
437  }
438 
439  if (_verbose==1) // printing for non verbose
440  this->printOutput(std::cout,i,def);
441 
442  _prec.post(x); // postprocess preconditioner
443  res.iterations = i; // fill statistics
444  res.reduction = def/def0;
445  res.conv_rate = pow(res.reduction,1.0/i);
446  res.elapsed = watch.elapsed();
447  if (_verbose>0) // final print
448  std::cout << "=== rate=" << res.conv_rate
449  << ", T=" << res.elapsed
450  << ", TIT=" << res.elapsed/i
451  << ", IT=" << i << std::endl;
452  }
453 
459  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
460  {
461  std::swap(_reduction,reduction);
462  (*this).apply(x,b,res);
463  std::swap(_reduction,reduction);
464  }
465 
466  private:
468  LinearOperator<X,X>& _op;
469  Preconditioner<X,X>& _prec;
470  ScalarProduct<X>& _sp;
471  double _reduction;
472  int _maxit;
473  int _verbose;
474  };
475 
476 
477 
479  template<class X>
480  class CGSolver : public InverseOperator<X,X> {
481  public:
483  typedef X domain_type;
485  typedef X range_type;
487  typedef typename X::field_type field_type;
488 
494  template<class L, class P>
495  CGSolver (L& op, P& prec, double reduction, int maxit, int verbose) :
496  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
497  {
498  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
499  "L and P must have the same category!");
500  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
501  "L must be sequential!");
502  }
508  template<class L, class S, class P>
509  CGSolver (L& op, S& sp, P& prec, double reduction, int maxit, int verbose) :
510  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
511  {
512  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
513  "L and P must have the same category!");
514  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
515  "L and S must have the same category!");
516  }
517 
523  virtual void apply (X& x, X& b, InverseOperatorResult& res)
524  {
525  res.clear(); // clear solver statistics
526  Timer watch; // start a timer
527  _prec.pre(x,b); // prepare preconditioner
528  _op.applyscaleadd(-1,x,b); // overwrite b with defect
529 
530  X p(x); // the search direction
531  X q(x); // a temporary vector
532 
533  double def0 = _sp.norm(b);// compute norm
534  if (def0<1E-30) // convergence check
535  {
536  res.converged = true;
537  res.iterations = 0; // fill statistics
538  res.reduction = 0;
539  res.conv_rate = 0;
540  res.elapsed=0;
541  if (_verbose>0) // final print
542  std::cout << "=== rate=" << res.conv_rate
543  << ", T=" << res.elapsed << ", TIT=" << res.elapsed
544  << ", IT=0" << std::endl;
545  return;
546  }
547 
548  if (_verbose>0) // printing
549  {
550  std::cout << "=== CGSolver" << std::endl;
551  if (_verbose>1) {
552  this->printHeader(std::cout);
553  this->printOutput(std::cout,0,def0);
554  }
555  }
556 
557  // some local variables
558  double def=def0; // loop variables
559  field_type rho,rholast,lambda,alpha,beta;
560 
561  // determine initial search direction
562  p = 0; // clear correction
563  _prec.apply(p,b); // apply preconditioner
564  rholast = _sp.dot(p,b); // orthogonalization
565 
566  // the loop
567  int i=1;
568  for ( ; i<=_maxit; i++ )
569  {
570  // minimize in given search direction p
571  _op.apply(p,q); // q=Ap
572  alpha = _sp.dot(p,q); // scalar product
573  lambda = rholast/alpha; // minimization
574  x.axpy(lambda,p); // update solution
575  b.axpy(-lambda,q); // update defect
576 
577  // convergence test
578  double defnew=_sp.norm(b);// comp defect norm
579 
580  if (_verbose>1) // print
581  this->printOutput(std::cout,i,defnew,def);
582 
583  def = defnew; // update norm
584  if (def<def0*_reduction || def<1E-30) // convergence check
585  {
586  res.converged = true;
587  break;
588  }
589 
590  // determine new search direction
591  q = 0; // clear correction
592  _prec.apply(q,b); // apply preconditioner
593  rho = _sp.dot(q,b); // orthogonalization
594  beta = rho/rholast; // scaling factor
595  p *= beta; // scale old search direction
596  p += q; // orthogonalization with correction
597  rholast = rho; // remember rho for recurrence
598  }
599 
600  if (_verbose==1) // printing for non verbose
601  this->printOutput(std::cout,i,def);
602 
603  _prec.post(x); // postprocess preconditioner
604  res.iterations = i; // fill statistics
605  res.reduction = def/def0;
606  res.conv_rate = pow(res.reduction,1.0/i);
607  res.elapsed = watch.elapsed();
608 
609  if (_verbose>0) // final print
610  {
611  std::cout << "=== rate=" << res.conv_rate
612  << ", T=" << res.elapsed
613  << ", TIT=" << res.elapsed/i
614  << ", IT=" << i << std::endl;
615  }
616  }
617 
623  virtual void apply (X& x, X& b, double reduction,
625  {
626  std::swap(_reduction,reduction);
627  (*this).apply(x,b,res);
628  std::swap(_reduction,reduction);
629  }
630 
631  private:
633  LinearOperator<X,X>& _op;
634  Preconditioner<X,X>& _prec;
635  ScalarProduct<X>& _sp;
636  double _reduction;
637  int _maxit;
638  int _verbose;
639  };
640 
641 
642  // Ronald Kriemanns BiCG-STAB implementation from Sumo
644  template<class X>
645  class BiCGSTABSolver : public InverseOperator<X,X> {
646  public:
648  typedef X domain_type;
650  typedef X range_type;
652  typedef typename X::field_type field_type;
653 
659  template<class L, class P>
660  BiCGSTABSolver (L& op, P& prec,
661  double reduction, int maxit, int verbose) :
662  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
663  {
664  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category), "L and P must be of the same category!");
665  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential), "L must be sequential!");
666  }
672  template<class L, class S, class P>
673  BiCGSTABSolver (L& op, S& sp, P& prec,
674  double reduction, int maxit, int verbose) :
675  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
676  {
677  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
678  "L and P must have the same category!");
679  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
680  "L and S must have the same category!");
681  }
682 
688  virtual void apply (X& x, X& b, InverseOperatorResult& res)
689  {
690  const double EPSILON=1e-80;
691 
692  double it;
693  field_type rho, rho_new, alpha, beta, h, omega;
694  field_type norm, norm_old, norm_0;
695 
696  //
697  // get vectors and matrix
698  //
699  X& r=b;
700  X p(x);
701  X v(x);
702  X t(x);
703  X y(x);
704  X rt(x);
705 
706  //
707  // begin iteration
708  //
709 
710  // r = r - Ax; rt = r
711  res.clear(); // clear solver statistics
712  Timer watch; // start a timer
713  _prec.pre(x,r); // prepare preconditioner
714  _op.applyscaleadd(-1,x,r); // overwrite b with defect
715 
716  rt=r;
717 
718  norm = norm_old = norm_0 = _sp.norm(r);
719 
720  p=0;
721  v=0;
722 
723  rho = 1;
724  alpha = 1;
725  omega = 1;
726 
727  if (_verbose>0) // printing
728  {
729  std::cout << "=== BiCGSTABSolver" << std::endl;
730  if (_verbose>1)
731  {
732  this->printHeader(std::cout);
733  this->printOutput(std::cout,0,norm_0);
734  //std::cout << " Iter Defect Rate" << std::endl;
735  //std::cout << " 0" << std::setw(14) << norm_0 << std::endl;
736  }
737  }
738 
739  if ( norm < (_reduction * norm_0) || norm<1E-30)
740  {
741  res.converged = 1;
742  _prec.post(x); // postprocess preconditioner
743  res.iterations = 0; // fill statistics
744  res.reduction = 0;
745  res.conv_rate = 0;
746  res.elapsed = watch.elapsed();
747  return;
748  }
749 
750  //
751  // iteration
752  //
753 
754  for (it = 0.5; it < _maxit; it+=.5)
755  {
756  //
757  // preprocess, set vecsizes etc.
758  //
759 
760  // rho_new = < rt , r >
761  rho_new = _sp.dot(rt,r);
762 
763  // look if breakdown occured
764  if (std::abs(rho) <= EPSILON)
765  DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho "
766  << rho << " <= EPSILON " << EPSILON
767  << " after " << it << " iterations");
768  if (std::abs(omega) <= EPSILON)
769  DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - omega "
770  << omega << " <= EPSILON " << EPSILON
771  << " after " << it << " iterations");
772 
773 
774  if (it<1)
775  p = r;
776  else
777  {
778  beta = ( rho_new / rho ) * ( alpha / omega );
779  p.axpy(-omega,v); // p = r + beta (p - omega*v)
780  p *= beta;
781  p += r;
782  }
783 
784  // y = W^-1 * p
785  y = 0;
786  _prec.apply(y,p); // apply preconditioner
787 
788  // v = A * y
789  _op.apply(y,v);
790 
791  // alpha = rho_new / < rt, v >
792  h = _sp.dot(rt,v);
793 
794  if ( std::abs(h) < EPSILON )
795  DUNE_THROW(ISTLError,"h=0 in BiCGSTAB");
796 
797  alpha = rho_new / h;
798 
799  // apply first correction to x
800  // x <- x + alpha y
801  x.axpy(alpha,y);
802 
803  // r = r - alpha*v
804  r.axpy(-alpha,v);
805 
806  //
807  // test stop criteria
808  //
809 
810  norm = _sp.norm(r);
811 
812  if (_verbose>1) // print
813  {
814  this->printOutput(std::cout,it,norm,norm_old);
815  }
816 
817  if ( norm < (_reduction * norm_0) )
818  {
819  res.converged = 1;
820  break;
821  }
822  it+=.5;
823 
824  norm_old = norm;
825 
826  // y = W^-1 * r
827  y = 0;
828  _prec.apply(y,r);
829 
830  // t = A * y
831  _op.apply(y,t);
832 
833  // omega = < t, r > / < t, t >
834  omega = _sp.dot(t,r)/_sp.dot(t,t);
835 
836  // apply second correction to x
837  // x <- x + omega y
838  x.axpy(omega,y);
839 
840  // r = s - omega*t (remember : r = s)
841  r.axpy(-omega,t);
842 
843  rho = rho_new;
844 
845  //
846  // test stop criteria
847  //
848 
849  norm = _sp.norm(r);
850 
851  if (_verbose > 1) // print
852  {
853  this->printOutput(std::cout,it,norm,norm_old);
854  }
855 
856  if ( norm < (_reduction * norm_0) || norm<1E-30)
857  {
858  res.converged = 1;
859  break;
860  }
861 
862  norm_old = norm;
863  } // end for
864 
865  if (_verbose==1) // printing for non verbose
866  this->printOutput(std::cout,it,norm);
867 
868  _prec.post(x); // postprocess preconditioner
869  res.iterations = static_cast<int>(std::ceil(it)); // fill statistics
870  res.reduction = norm/norm_0;
871  res.conv_rate = pow(res.reduction,1.0/it);
872  res.elapsed = watch.elapsed();
873  if (_verbose>0) // final print
874  std::cout << "=== rate=" << res.conv_rate
875  << ", T=" << res.elapsed
876  << ", TIT=" << res.elapsed/it
877  << ", IT=" << it << std::endl;
878  }
879 
885  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
886  {
887  std::swap(_reduction,reduction);
888  (*this).apply(x,b,res);
889  std::swap(_reduction,reduction);
890  }
891 
892  private:
894  LinearOperator<X,X>& _op;
895  Preconditioner<X,X>& _prec;
896  ScalarProduct<X>& _sp;
897  double _reduction;
898  int _maxit;
899  int _verbose;
900  };
901 
908  template<class X>
909  class MINRESSolver : public InverseOperator<X,X> {
910  public:
912  typedef X domain_type;
914  typedef X range_type;
916  typedef typename X::field_type field_type;
918  typedef typename FieldTraits<field_type>::real_type real_type;
919 
925  template<class L, class P>
926  MINRESSolver (L& op, P& prec, double reduction, int maxit, int verbose) :
927  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
928  {
929  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
930  "L and P must have the same category!");
931  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
932  "L must be sequential!");
933  }
939  template<class L, class S, class P>
940  MINRESSolver (L& op, S& sp, P& prec, double reduction, int maxit, int verbose) :
941  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
942  {
943  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
944  "L and P must have the same category!");
945  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
946  "L and S must have the same category!");
947  }
948 
954  virtual void apply (X& x, X& b, InverseOperatorResult& res)
955  {
956  res.clear(); // clear solver statistics
957  Timer watch; // start a timer
958  _prec.pre(x,b); // prepare preconditioner
959  _op.applyscaleadd(-1,x,b); // overwrite b with defect/residual
960 
961  real_type def0 = _sp.norm(b); // compute residual norm
962 
963  if (def0<1E-30) // convergence check
964  {
965  res.converged = true;
966  res.iterations = 0; // fill statistics
967  res.reduction = 0;
968  res.conv_rate = 0;
969  res.elapsed=0;
970  if (_verbose>0) // final print
971  std::cout << "=== rate=" << res.conv_rate << ", T=" << res.elapsed << ", TIT=" << res.elapsed << ", IT=0" << std::endl;
972  return;
973  }
974 
975  if (_verbose>0) // printing
976  {
977  std::cout << "=== MINRESSolver" << std::endl;
978  if (_verbose>1) {
979  this->printHeader(std::cout);
980  this->printOutput(std::cout,0,def0);
981  }
982  }
983 
984  // some local variables
985  real_type def=def0; // the defect/residual norm
986  field_type alpha, // recurrence coefficients as computed in the Lanczos alg making up the matrix T
987  c[2]={0.0, 0.0}, // diagonal entry of Givens rotation
988  s[2]={0.0, 0.0}; // off-diagonal entries of Givens rotation
989  real_type beta;
990 
991  field_type T[3]={0.0, 0.0, 0.0}; // recurrence coefficients (column k of Matrix T)
992 
993  X z(b.size()), // some temporary vectors
994  dummy(b.size());
995 
996  field_type xi[2]={1.0, 0.0};
997 
998  // initialize
999  z = 0.0; // clear correction
1000 
1001  _prec.apply(z,b); // apply preconditioner z=M^-1*b
1002 
1003  beta = std::sqrt(std::abs(_sp.dot(z,b)));
1004  real_type beta0 = beta;
1005 
1006  X p[3]; // the search directions
1007  X q[3]; // Orthonormal basis vectors (in unpreconditioned case)
1008 
1009  q[0].resize(b.size());
1010  q[1].resize(b.size());
1011  q[2].resize(b.size());
1012  q[0] = 0.0;
1013  q[1] = b;
1014  q[1] /= beta;
1015  q[2] = 0.0;
1016 
1017  p[0].resize(b.size());
1018  p[1].resize(b.size());
1019  p[2].resize(b.size());
1020  p[0] = 0.0;
1021  p[1] = 0.0;
1022  p[2] = 0.0;
1023 
1024 
1025  z /= beta; // this is w_current
1026 
1027  // the loop
1028  int i=1;
1029  for ( ; i<=_maxit; i++)
1030  {
1031  dummy = z; // remember z_old for the computation of the search direction p in the next iteration
1032 
1033  int i1 = i%3,
1034  i0 = (i1+2)%3,
1035  i2 = (i1+1)%3;
1036 
1037  // Symmetrically Preconditioned Lanczos (Greenbaum p.121)
1038  _op.apply(z,q[i2]); // q[i2] = Az
1039  q[i2].axpy(-beta, q[i0]);
1040  alpha = _sp.dot(q[i2],z);
1041  q[i2].axpy(-alpha, q[i1]);
1042 
1043  z=0.0;
1044  _prec.apply(z,q[i2]);
1045 
1046  beta = std::sqrt(std::abs(_sp.dot(q[i2],z)));
1047 
1048  q[i2] /= beta;
1049  z /= beta;
1050 
1051  // QR Factorization of recurrence coefficient matrix
1052  // apply previous Givens rotations to last column of T
1053  T[1] = T[2];
1054  if (i>2)
1055  {
1056  T[0] = s[i%2]*T[1];
1057  T[1] = c[i%2]*T[1];
1058  }
1059  if (i>1)
1060  {
1061  T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
1062  T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
1063  }
1064  else
1065  T[2] = alpha;
1066 
1067  // recompute c, s -> current Givens rotation \TODO use BLAS-routine drotg instead for greater robustness
1068 // cblas_drotg (a, b, c, s);
1069  c[i%2] = 1.0/std::sqrt(T[2]*T[2] + beta*beta);
1070  s[i%2] = beta*c[i%2];
1071  c[i%2] *= T[2];
1072 
1073  // apply current Givens rotation to T eliminating the last entry...
1074  T[2] = c[i%2]*T[2] + s[i%2]*beta;
1075 
1076  // ...and to xi, the right hand side of the least squares problem min_y||beta*xi-T*y||
1077  xi[i%2] = -s[i%2]*xi[(i+1)%2];
1078  xi[(i+1)%2] *= c[i%2];
1079 
1080  // compute correction direction
1081  p[i2] = dummy;
1082  p[i2].axpy(-T[1],p[i1]);
1083  p[i2].axpy(-T[0],p[i0]);
1084  p[i2] /= T[2];
1085 
1086  // apply correction/update solution
1087  x.axpy(beta0*xi[(i+1)%2], p[i2]);
1088 
1089  // remember beta_old
1090  T[2] = beta;
1091 
1092  // update residual - not necessary if in the preconditioned case we are content with the residual norm of the
1093  // preconditioned system as convergence test
1094 // _op.apply(p[i2],dummy);
1095 // b.axpy(-beta0*xi[(i+1)%2],dummy);
1096 
1097 // convergence test
1098  real_type defnew = std::abs(beta0*xi[i%2]); // the last entry the QR-transformed least squares RHS is the new residual norm
1099 
1100  if (_verbose>1) // print
1101  this->printOutput(std::cout,i,defnew,def);
1102 
1103  def = defnew; // update norm
1104  if (def<def0*_reduction || def<1E-30 || i==_maxit) // convergence check
1105  {
1106  res.converged = true;
1107  break;
1108  }
1109  }
1110 
1111  if (_verbose==1) // printing for non verbose
1112  this->printOutput(std::cout,i,def);
1113 
1114  _prec.post(x); // postprocess preconditioner
1115  res.iterations = i; // fill statistics
1116  res.reduction = def/def0;
1117  res.conv_rate = pow(res.reduction,1.0/i);
1118  res.elapsed = watch.elapsed();
1119 
1120  if (_verbose>0) // final print
1121  {
1122  std::cout << "=== rate=" << res.conv_rate
1123  << ", T=" << res.elapsed
1124  << ", TIT=" << res.elapsed/i
1125  << ", IT=" << i << std::endl;
1126  }
1127 
1128  }
1129 
1135  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
1136  {
1137  std::swap(_reduction,reduction);
1138  (*this).apply(x,b,res);
1139  std::swap(_reduction,reduction);
1140  }
1141 
1142  private:
1143  SeqScalarProduct<X> ssp;
1144  LinearOperator<X,X>& _op;
1145  Preconditioner<X,X>& _prec;
1146  ScalarProduct<X>& _sp;
1147  double _reduction;
1148  int _maxit;
1149  int _verbose;
1150  };
1151 
1163  template<class X, class Y=X, class F = Y>
1165  {
1166  public:
1168  typedef X domain_type;
1170  typedef Y range_type;
1172  typedef typename X::field_type field_type;
1174  typedef F basis_type;
1175 
1183  template<class L, class P>
1184  RestartedGMResSolver (L& op, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) :
1185  _A_(op), _M(prec),
1186  ssp(), _sp(ssp), _restart(restart),
1187  _reduction(reduction), _maxit(maxit), _verbose(verbose),
1188  _recalc_defect(recalc_defect)
1189  {
1190  dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1191  "P and L must be the same category!");
1192  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
1193  "L must be sequential!");
1194  }
1195 
1203  template<class L, class S, class P>
1204  RestartedGMResSolver (L& op, S& sp, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) :
1205  _A_(op), _M(prec),
1206  _sp(sp), _restart(restart),
1207  _reduction(reduction), _maxit(maxit), _verbose(verbose),
1208  _recalc_defect(recalc_defect)
1209  {
1210  dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1211  "P and L must have the same category!");
1212  dune_static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1213  "P and S must have the same category!");
1214  }
1215 
1217  virtual void apply (X& x, X& b, InverseOperatorResult& res)
1218  {
1219  apply(x,b,_reduction,res);
1220  };
1221 
1227  virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res)
1228  {
1229  int m = _restart;
1230  field_type norm;
1231  field_type norm_old = 0.0;
1232  field_type norm_0;
1233  field_type beta;
1234  int i, j = 1, k;
1235  std::vector<field_type> s(m+1), cs(m), sn(m);
1236  // helper vector
1237  X w(b);
1238  std::vector< std::vector<field_type> > H(m+1,s);
1239  std::vector<F> v(m+1,b);
1240 
1241  // start timer
1242  Timer watch; // start a timer
1243 
1244  // clear solver statistics
1245  res.clear();
1246  _M.pre(x,b);
1247  if (_recalc_defect)
1248  {
1249  // norm_0 = norm(M^-1 b)
1250  w = 0.0; _M.apply(w,b); // w = M^-1 b
1251  norm_0 = _sp.norm(w);
1252  // r = _M.solve(b - A * x);
1253  w = b;
1254  _A_.applyscaleadd(-1,x, /* => */ w); // w = b - Ax;
1255  v[0] = 0.0; _M.apply(v[0],w); // r = M^-1 w
1256  beta = _sp.norm(v[0]);
1257  }
1258  else
1259  {
1260  // norm_0 = norm(M^-1 b)
1261  w = 0.0; _M.apply(w,b); // w = M^-1 b
1262  // r = _M.solve(b - A * x);
1263  _A_.applyscaleadd(-1,x, /* => */ b); // b = b - Ax;
1264  norm_0 = _sp.norm(b);
1265  v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b
1266  beta = _sp.norm(v[0]);
1267  }
1268 
1269  // avoid division by zero
1270  if (norm_0 == 0.0)
1271  norm_0 = 1.0;
1272  norm = norm_old = _sp.norm(v[0]);
1273 
1274  // print header
1275  if (_verbose > 0)
1276  {
1277  std::cout << "=== RestartedGMResSolver" << std::endl;
1278  if (_verbose > 1)
1279  {
1280  this->printHeader(std::cout);
1281  this->printOutput(std::cout,0,norm_0);
1282  this->printOutput(std::cout,0,norm, norm_0);
1283  }
1284  }
1285 
1286  // check convergence
1287  if (norm <= reduction * norm_0) {
1288  _M.post(x); // postprocess preconditioner
1289  res.converged = true;
1290  if (_verbose > 0) // final print
1291  print_result(res);
1292  return;
1293  }
1294 
1295  while (j <= _maxit && res.converged != true) {
1296  v[0] *= (1.0 / beta);
1297  for (i=1; i<=m; i++) s[i] = 0.0;
1298  s[0] = beta;
1299 
1300  for (i = 0; i < m && j <= _maxit && res.converged != true; i++, j++) {
1301  w = 0.0;
1302  v[i+1] = 0.0; // use v[i+1] as temporary vector
1303  _A_.apply(v[i], /* => */ v[i+1]);
1304  _M.apply(w, v[i+1]);
1305  for (k = 0; k <= i; k++) {
1306  H[k][i] = _sp.dot(w, v[k]);
1307  // w -= H[k][i] * v[k];
1308  w.axpy(-H[k][i], v[k]);
1309  }
1310  H[i+1][i] = _sp.norm(w);
1311  if (H[i+1][i] == 0.0)
1312  DUNE_THROW(ISTLError,"breakdown in GMRes - |w| "
1313  << w << " == 0.0 after " << j << " iterations");
1314  // v[i+1] = w * (1.0 / H[i+1][i]);
1315  v[i+1] = w; v[i+1] *= (1.0 / H[i+1][i]);
1316 
1317  for (k = 0; k < i; k++)
1318  applyPlaneRotation(H[k][i], H[k+1][i], cs[k], sn[k]);
1319 
1320  generatePlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
1321  applyPlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
1322  applyPlaneRotation(s[i], s[i+1], cs[i], sn[i]);
1323 
1324  norm = std::abs(s[i+1]);
1325 
1326  if (_verbose > 1) // print
1327  {
1328  this->printOutput(std::cout,j,norm,norm_old);
1329  }
1330 
1331  norm_old = norm;
1332 
1333  if (norm < reduction * norm_0) {
1334  res.converged = true;
1335  }
1336  }
1337 
1338  if (_recalc_defect)
1339  {
1340  // update x
1341  update(x, i - 1, H, s, v);
1342 
1343  // update residuum
1344  // r = M^-1 (b - A * x);
1345  w = b; _A_.applyscaleadd(-1,x, /* => */ w);
1346  _M.apply(v[0], w);
1347  beta = _sp.norm(v[0]);
1348  norm = beta;
1349  }
1350  else
1351  {
1352  // calc update vector
1353  w = 0;
1354  update(w, i - 1, H, s, v);
1355 
1356  // update x
1357  x += w;
1358 
1359  // r = M^-1 (b - A * x);
1360  // update defect
1361  _A_.applyscaleadd(-1,w, /* => */ b);
1362  // r = M^-1 (b - A * x);
1363  v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b
1364  beta = _sp.norm(v[0]);
1365  norm = beta;
1366 
1367  res.converged = false;
1368  }
1369 
1370  if (_verbose > 1) // print
1371  {
1372  this->printOutput(std::cout,j,norm,norm_old);
1373  }
1374 
1375  norm_old = norm;
1376 
1377  if (norm < reduction * norm_0) {
1378  // fill statistics
1379  res.converged = true;
1380  }
1381 
1382  if (res.converged != true && _verbose > 0)
1383  std::cout << "=== GMRes::restart\n";
1384  }
1385 
1386  _M.post(x); // postprocess preconditioner
1387 
1388  res.iterations = j;
1389  res.reduction = norm / norm_0;
1390  res.conv_rate = pow(res.reduction,1.0/j);
1391  res.elapsed = watch.elapsed();
1392 
1393  if (_verbose>0)
1394  print_result(res);
1395  }
1396  private:
1397 
1398  void
1399  print_result (const InverseOperatorResult & res) const
1400  {
1401  int j = res.iterations>0?res.iterations:1;
1402  std::cout << "=== rate=" << res.conv_rate
1403  << ", T=" << res.elapsed
1404  << ", TIT=" << res.elapsed/j
1405  << ", IT=" << res.iterations
1406  << std::endl;
1407  }
1408 
1409  static void
1410  update(X &x, int k,
1411  std::vector< std::vector<field_type> > & h,
1412  std::vector<field_type> & s, std::vector<F> v)
1413  {
1414  std::vector<field_type> y(s);
1415 
1416  // Backsolve:
1417  for (int i = k; i >= 0; i--) {
1418  y[i] /= h[i][i];
1419  for (int j = i - 1; j >= 0; j--)
1420  y[j] -= h[j][i] * y[i];
1421  }
1422 
1423  for (int j = 0; j <= k; j++)
1424  // x += v[j] * y[j];
1425  x.axpy(y[j],v[j]);
1426  }
1427 
1428  void
1429  generatePlaneRotation(field_type &dx, field_type &dy, field_type &cs, field_type &sn)
1430  {
1431  if (dy == 0.0) {
1432  cs = 1.0;
1433  sn = 0.0;
1434  } else if (std::abs(dy) > std::abs(dx)) {
1435  field_type temp = dx / dy;
1436  sn = 1.0 / std::sqrt( 1.0 + temp*temp );
1437  cs = temp * sn;
1438  } else {
1439  field_type temp = dy / dx;
1440  cs = 1.0 / std::sqrt( 1.0 + temp*temp );
1441  sn = temp * cs;
1442  }
1443  }
1444 
1445 
1446  void
1447  applyPlaneRotation(field_type &dx, field_type &dy, field_type &cs, field_type &sn)
1448  {
1449  field_type temp = cs * dx + sn * dy;
1450  dy = -sn * dx + cs * dy;
1451  dx = temp;
1452  }
1453 
1454  LinearOperator<X,X>& _A_;
1455  Preconditioner<X,X>& _M;
1456  SeqScalarProduct<X> ssp;
1457  ScalarProduct<X>& _sp;
1458  int _restart;
1459  double _reduction;
1460  int _maxit;
1461  int _verbose;
1462  bool _recalc_defect;
1463  };
1464 
1467 } // end namespace
1468 
1469 #endif