dune-istl  2.2.1
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;
654  typedef typename FieldTraits<field_type>::real_type real_type;
655 
661  template<class L, class P>
662  BiCGSTABSolver (L& op, P& prec,
663  double reduction, int maxit, int verbose) :
664  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
665  {
666  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category), "L and P must be of the same category!");
667  dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential), "L must be sequential!");
668  }
674  template<class L, class S, class P>
675  BiCGSTABSolver (L& op, S& sp, P& prec,
676  double reduction, int maxit, int verbose) :
677  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
678  {
679  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
680  "L and P must have the same category!");
681  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
682  "L and S must have the same category!");
683  }
684 
690  virtual void apply (X& x, X& b, InverseOperatorResult& res)
691  {
692  const double EPSILON=1e-80;
693  double it;
694  field_type rho, rho_new, alpha, beta, h, omega;
695  real_type norm, norm_old, norm_0;
696 
697  //
698  // get vectors and matrix
699  //
700  X& r=b;
701  X p(x);
702  X v(x);
703  X t(x);
704  X y(x);
705  X rt(x);
706 
707  //
708  // begin iteration
709  //
710 
711  // r = r - Ax; rt = r
712  res.clear(); // clear solver statistics
713  Timer watch; // start a timer
714  _prec.pre(x,r); // prepare preconditioner
715  _op.applyscaleadd(-1,x,r); // overwrite b with defect
716 
717  rt=r;
718 
719  norm = norm_old = norm_0 = _sp.norm(r);
720 
721  p=0;
722  v=0;
723 
724  rho = 1;
725  alpha = 1;
726  omega = 1;
727 
728  if (_verbose>0) // printing
729  {
730  std::cout << "=== BiCGSTABSolver" << std::endl;
731  if (_verbose>1)
732  {
733  this->printHeader(std::cout);
734  this->printOutput(std::cout,0,norm_0);
735  //std::cout << " Iter Defect Rate" << std::endl;
736  //std::cout << " 0" << std::setw(14) << norm_0 << std::endl;
737  }
738  }
739 
740  if ( norm < (_reduction * norm_0) || norm<1E-30)
741  {
742  res.converged = 1;
743  _prec.post(x); // postprocess preconditioner
744  res.iterations = 0; // fill statistics
745  res.reduction = 0;
746  res.conv_rate = 0;
747  res.elapsed = watch.elapsed();
748  return;
749  }
750 
751  //
752  // iteration
753  //
754 
755  for (it = 0.5; it < _maxit; it+=.5)
756  {
757  //
758  // preprocess, set vecsizes etc.
759  //
760 
761  // rho_new = < rt , r >
762  rho_new = _sp.dot(rt,r);
763 
764  // look if breakdown occured
765  if (std::abs(rho) <= EPSILON)
766  DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho "
767  << rho << " <= EPSILON " << EPSILON
768  << " after " << it << " iterations");
769  if (std::abs(omega) <= EPSILON)
770  DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - omega "
771  << omega << " <= EPSILON " << EPSILON
772  << " after " << it << " iterations");
773 
774 
775  if (it<1)
776  p = r;
777  else
778  {
779  beta = ( rho_new / rho ) * ( alpha / omega );
780  p.axpy(-omega,v); // p = r + beta (p - omega*v)
781  p *= beta;
782  p += r;
783  }
784 
785  // y = W^-1 * p
786  y = 0;
787  _prec.apply(y,p); // apply preconditioner
788 
789  // v = A * y
790  _op.apply(y,v);
791 
792  // alpha = rho_new / < rt, v >
793  h = _sp.dot(rt,v);
794 
795  if ( std::abs(h) < EPSILON )
796  DUNE_THROW(ISTLError,"h=0 in BiCGSTAB");
797 
798  alpha = rho_new / h;
799 
800  // apply first correction to x
801  // x <- x + alpha y
802  x.axpy(alpha,y);
803 
804  // r = r - alpha*v
805  r.axpy(-alpha,v);
806 
807  //
808  // test stop criteria
809  //
810 
811  norm = _sp.norm(r);
812 
813  if (_verbose>1) // print
814  {
815  this->printOutput(std::cout,it,norm,norm_old);
816  }
817 
818  if ( norm < (_reduction * norm_0) )
819  {
820  res.converged = 1;
821  break;
822  }
823  it+=.5;
824 
825  norm_old = norm;
826 
827  // y = W^-1 * r
828  y = 0;
829  _prec.apply(y,r);
830 
831  // t = A * y
832  _op.apply(y,t);
833 
834  // omega = < t, r > / < t, t >
835  omega = _sp.dot(t,r)/_sp.dot(t,t);
836 
837  // apply second correction to x
838  // x <- x + omega y
839  x.axpy(omega,y);
840 
841  // r = s - omega*t (remember : r = s)
842  r.axpy(-omega,t);
843 
844  rho = rho_new;
845 
846  //
847  // test stop criteria
848  //
849 
850  norm = _sp.norm(r);
851 
852  if (_verbose > 1) // print
853  {
854  this->printOutput(std::cout,it,norm,norm_old);
855  }
856 
857  if ( norm < (_reduction * norm_0) || norm<1E-30)
858  {
859  res.converged = 1;
860  break;
861  }
862 
863  norm_old = norm;
864  } // end for
865 
866  if (_verbose==1) // printing for non verbose
867  this->printOutput(std::cout,it,norm);
868 
869  _prec.post(x); // postprocess preconditioner
870  res.iterations = static_cast<int>(std::ceil(it)); // fill statistics
871  res.reduction = norm/norm_0;
872  res.conv_rate = pow(res.reduction,1.0/it);
873  res.elapsed = watch.elapsed();
874  if (_verbose>0) // final print
875  std::cout << "=== rate=" << res.conv_rate
876  << ", T=" << res.elapsed
877  << ", TIT=" << res.elapsed/it
878  << ", IT=" << it << std::endl;
879  }
880 
886  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
887  {
888  std::swap(_reduction,reduction);
889  (*this).apply(x,b,res);
890  std::swap(_reduction,reduction);
891  }
892 
893  private:
895  LinearOperator<X,X>& _op;
896  Preconditioner<X,X>& _prec;
897  ScalarProduct<X>& _sp;
898  double _reduction;
899  int _maxit;
900  int _verbose;
901  };
902 
909  template<class X>
910  class MINRESSolver : public InverseOperator<X,X> {
911  public:
913  typedef X domain_type;
915  typedef X range_type;
917  typedef typename X::field_type field_type;
919  typedef typename FieldTraits<field_type>::real_type real_type;
920 
926  template<class L, class P>
927  MINRESSolver (L& op, P& prec, double reduction, int maxit, int verbose) :
928  ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
929  {
930  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
931  "L and P must have the same category!");
932  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
933  "L must be sequential!");
934  }
940  template<class L, class S, class P>
941  MINRESSolver (L& op, S& sp, P& prec, double reduction, int maxit, int verbose) :
942  _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
943  {
944  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category),
945  "L and P must have the same category!");
946  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category),
947  "L and S must have the same category!");
948  }
949 
955  virtual void apply (X& x, X& b, InverseOperatorResult& res)
956  {
957  res.clear(); // clear solver statistics
958  Timer watch; // start a timer
959  _prec.pre(x,b); // prepare preconditioner
960  _op.applyscaleadd(-1,x,b); // overwrite b with defect/residual
961 
962  real_type def0 = _sp.norm(b); // compute residual norm
963 
964  if (def0<1E-30) // convergence check
965  {
966  res.converged = true;
967  res.iterations = 0; // fill statistics
968  res.reduction = 0;
969  res.conv_rate = 0;
970  res.elapsed=0;
971  if (_verbose>0) // final print
972  std::cout << "=== rate=" << res.conv_rate << ", T=" << res.elapsed << ", TIT=" << res.elapsed << ", IT=0" << std::endl;
973  return;
974  }
975 
976  if (_verbose>0) // printing
977  {
978  std::cout << "=== MINRESSolver" << std::endl;
979  if (_verbose>1) {
980  this->printHeader(std::cout);
981  this->printOutput(std::cout,0,def0);
982  }
983  }
984 
985  // some local variables
986  real_type def=def0; // the defect/residual norm
987  field_type alpha, // recurrence coefficients as computed in the Lanczos alg making up the matrix T
988  c[2]={0.0, 0.0}, // diagonal entry of Givens rotation
989  s[2]={0.0, 0.0}; // off-diagonal entries of Givens rotation
990  real_type beta;
991 
992  field_type T[3]={0.0, 0.0, 0.0}; // recurrence coefficients (column k of Matrix T)
993 
994  X z(b.size()), // some temporary vectors
995  dummy(b.size());
996 
997  field_type xi[2]={1.0, 0.0};
998 
999  // initialize
1000  z = 0.0; // clear correction
1001 
1002  _prec.apply(z,b); // apply preconditioner z=M^-1*b
1003 
1004  beta = std::sqrt(std::abs(_sp.dot(z,b)));
1005  real_type beta0 = beta;
1006 
1007  X p[3]; // the search directions
1008  X q[3]; // Orthonormal basis vectors (in unpreconditioned case)
1009 
1010  q[0].resize(b.size());
1011  q[1].resize(b.size());
1012  q[2].resize(b.size());
1013  q[0] = 0.0;
1014  q[1] = b;
1015  q[1] /= beta;
1016  q[2] = 0.0;
1017 
1018  p[0].resize(b.size());
1019  p[1].resize(b.size());
1020  p[2].resize(b.size());
1021  p[0] = 0.0;
1022  p[1] = 0.0;
1023  p[2] = 0.0;
1024 
1025 
1026  z /= beta; // this is w_current
1027 
1028  // the loop
1029  int i=1;
1030  for ( ; i<=_maxit; i++)
1031  {
1032  dummy = z; // remember z_old for the computation of the search direction p in the next iteration
1033 
1034  int i1 = i%3,
1035  i0 = (i1+2)%3,
1036  i2 = (i1+1)%3;
1037 
1038  // Symmetrically Preconditioned Lanczos (Greenbaum p.121)
1039  _op.apply(z,q[i2]); // q[i2] = Az
1040  q[i2].axpy(-beta, q[i0]);
1041  alpha = _sp.dot(q[i2],z);
1042  q[i2].axpy(-alpha, q[i1]);
1043 
1044  z=0.0;
1045  _prec.apply(z,q[i2]);
1046 
1047  beta = std::sqrt(std::abs(_sp.dot(q[i2],z)));
1048 
1049  q[i2] /= beta;
1050  z /= beta;
1051 
1052  // QR Factorization of recurrence coefficient matrix
1053  // apply previous Givens rotations to last column of T
1054  T[1] = T[2];
1055  if (i>2)
1056  {
1057  T[0] = s[i%2]*T[1];
1058  T[1] = c[i%2]*T[1];
1059  }
1060  if (i>1)
1061  {
1062  T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1];
1063  T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha;
1064  }
1065  else
1066  T[2] = alpha;
1067 
1068  // recompute c, s -> current Givens rotation \TODO use BLAS-routine drotg instead for greater robustness
1069 // cblas_drotg (a, b, c, s);
1070  c[i%2] = 1.0/std::sqrt(T[2]*T[2] + beta*beta);
1071  s[i%2] = beta*c[i%2];
1072  c[i%2] *= T[2];
1073 
1074  // apply current Givens rotation to T eliminating the last entry...
1075  T[2] = c[i%2]*T[2] + s[i%2]*beta;
1076 
1077  // ...and to xi, the right hand side of the least squares problem min_y||beta*xi-T*y||
1078  xi[i%2] = -s[i%2]*xi[(i+1)%2];
1079  xi[(i+1)%2] *= c[i%2];
1080 
1081  // compute correction direction
1082  p[i2] = dummy;
1083  p[i2].axpy(-T[1],p[i1]);
1084  p[i2].axpy(-T[0],p[i0]);
1085  p[i2] /= T[2];
1086 
1087  // apply correction/update solution
1088  x.axpy(beta0*xi[(i+1)%2], p[i2]);
1089 
1090  // remember beta_old
1091  T[2] = beta;
1092 
1093  // update residual - not necessary if in the preconditioned case we are content with the residual norm of the
1094  // preconditioned system as convergence test
1095 // _op.apply(p[i2],dummy);
1096 // b.axpy(-beta0*xi[(i+1)%2],dummy);
1097 
1098 // convergence test
1099  real_type defnew = std::abs(beta0*xi[i%2]); // the last entry the QR-transformed least squares RHS is the new residual norm
1100 
1101  if (_verbose>1) // print
1102  this->printOutput(std::cout,i,defnew,def);
1103 
1104  def = defnew; // update norm
1105  if (def<def0*_reduction || def<1E-30 || i==_maxit) // convergence check
1106  {
1107  res.converged = true;
1108  break;
1109  }
1110  }
1111 
1112  if (_verbose==1) // printing for non verbose
1113  this->printOutput(std::cout,i,def);
1114 
1115  _prec.post(x); // postprocess preconditioner
1116  res.iterations = i; // fill statistics
1117  res.reduction = def/def0;
1118  res.conv_rate = pow(res.reduction,1.0/i);
1119  res.elapsed = watch.elapsed();
1120 
1121  if (_verbose>0) // final print
1122  {
1123  std::cout << "=== rate=" << res.conv_rate
1124  << ", T=" << res.elapsed
1125  << ", TIT=" << res.elapsed/i
1126  << ", IT=" << i << std::endl;
1127  }
1128 
1129  }
1130 
1136  virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
1137  {
1138  std::swap(_reduction,reduction);
1139  (*this).apply(x,b,res);
1140  std::swap(_reduction,reduction);
1141  }
1142 
1143  private:
1144  SeqScalarProduct<X> ssp;
1145  LinearOperator<X,X>& _op;
1146  Preconditioner<X,X>& _prec;
1147  ScalarProduct<X>& _sp;
1148  double _reduction;
1149  int _maxit;
1150  int _verbose;
1151  };
1152 
1164  template<class X, class Y=X, class F = Y>
1166  {
1167  public:
1169  typedef X domain_type;
1171  typedef Y range_type;
1173  typedef typename X::field_type field_type;
1175  typedef typename FieldTraits<field_type>::real_type real_type;
1177  typedef F basis_type;
1178 
1186  template<class L, class P>
1187  RestartedGMResSolver (L& op, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) :
1188  _A_(op), _M(prec),
1189  ssp(), _sp(ssp), _restart(restart),
1190  _reduction(reduction), _maxit(maxit), _verbose(verbose),
1191  _recalc_defect(recalc_defect)
1192  {
1193  dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1194  "P and L must be the same category!");
1195  dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
1196  "L must be sequential!");
1197  }
1198 
1206  template<class L, class S, class P>
1207  RestartedGMResSolver (L& op, S& sp, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) :
1208  _A_(op), _M(prec),
1209  _sp(sp), _restart(restart),
1210  _reduction(reduction), _maxit(maxit), _verbose(verbose),
1211  _recalc_defect(recalc_defect)
1212  {
1213  dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category),
1214  "P and L must have the same category!");
1215  dune_static_assert(static_cast<int>(P::category) == static_cast<int>(S::category),
1216  "P and S must have the same category!");
1217  }
1218 
1220  virtual void apply (X& x, X& b, InverseOperatorResult& res)
1221  {
1222  apply(x,b,_reduction,res);
1223  }
1224 
1230  virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res)
1231  {
1232  int m = _restart;
1233  real_type norm;
1234  real_type norm_old = 0.0;
1235  real_type norm_0;
1236  real_type beta;
1237  int i, j = 1, k;
1238  std::vector<field_type> s(m+1), cs(m), sn(m);
1239  // helper vector
1240  X w(b);
1241  std::vector< std::vector<field_type> > H(m+1,s);
1242  std::vector<F> v(m+1,b);
1243 
1244  // start timer
1245  Timer watch; // start a timer
1246 
1247  // clear solver statistics
1248  res.clear();
1249  _M.pre(x,b);
1250  if (_recalc_defect)
1251  {
1252  // norm_0 = norm(M^-1 b)
1253  w = 0.0; _M.apply(w,b); // w = M^-1 b
1254  norm_0 = _sp.norm(w);
1255  // r = _M.solve(b - A * x);
1256  w = b;
1257  _A_.applyscaleadd(-1,x, /* => */ w); // w = b - Ax;
1258  v[0] = 0.0; _M.apply(v[0],w); // r = M^-1 w
1259  beta = _sp.norm(v[0]);
1260  }
1261  else
1262  {
1263  // norm_0 = norm(b-Ax)
1264  _A_.applyscaleadd(-1,x, /* => */ b); // b = b - Ax;
1265  norm_0 = _sp.norm(b);
1266  v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b
1267  beta = _sp.norm(v[0]);
1268  }
1269 
1270  // avoid division by zero
1271  if (norm_0 == 0.0)
1272  norm_0 = 1.0;
1273  norm = norm_old = _sp.norm(v[0]);
1274 
1275  // print header
1276  if (_verbose > 0)
1277  {
1278  std::cout << "=== RestartedGMResSolver" << std::endl;
1279  if (_verbose > 1)
1280  {
1281  this->printHeader(std::cout);
1282  this->printOutput(std::cout,0,norm_0);
1283  this->printOutput(std::cout,0,norm, norm_0);
1284  }
1285  }
1286 
1287  // check convergence
1288  if (norm <= reduction * norm_0) {
1289  _M.post(x); // postprocess preconditioner
1290  res.converged = true;
1291  if (_verbose > 0) // final print
1292  print_result(res);
1293  return;
1294  }
1295 
1296  while (j <= _maxit && res.converged != true) {
1297  v[0] *= (1.0 / beta);
1298  for (i=1; i<=m; i++) s[i] = 0.0;
1299  s[0] = beta;
1300 
1301  for (i = 0; i < m && j <= _maxit && res.converged != true; i++, j++) {
1302  w = 0.0;
1303  v[i+1] = 0.0; // use v[i+1] as temporary vector
1304  _A_.apply(v[i], /* => */ v[i+1]);
1305  _M.apply(w, v[i+1]);
1306  for (k = 0; k <= i; k++) {
1307  H[k][i] = _sp.dot(w, v[k]);
1308  // w -= H[k][i] * v[k];
1309  w.axpy(-H[k][i], v[k]);
1310  }
1311  H[i+1][i] = _sp.norm(w);
1312  if (H[i+1][i] == 0.0)
1313  DUNE_THROW(ISTLError,"breakdown in GMRes - |w| "
1314  << w << " == 0.0 after " << j << " iterations");
1315  // v[i+1] = w * (1.0 / H[i+1][i]);
1316  v[i+1] = w; v[i+1] *= (1.0 / H[i+1][i]);
1317 
1318  for (k = 0; k < i; k++)
1319  applyPlaneRotation(H[k][i], H[k+1][i], cs[k], sn[k]);
1320 
1321  generatePlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
1322  applyPlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]);
1323  applyPlaneRotation(s[i], s[i+1], cs[i], sn[i]);
1324 
1325  norm = std::abs(s[i+1]);
1326 
1327  if (_verbose > 1) // print
1328  {
1329  this->printOutput(std::cout,j,norm,norm_old);
1330  }
1331 
1332  norm_old = norm;
1333 
1334  if (norm < reduction * norm_0) {
1335  res.converged = true;
1336  }
1337  }
1338 
1339  if (_recalc_defect)
1340  {
1341  // update x
1342  update(x, i - 1, H, s, v);
1343 
1344  // update residuum
1345  // r = M^-1 (b - A * x);
1346  w = b; _A_.applyscaleadd(-1,x, /* => */ w);
1347  _M.apply(v[0], w);
1348  beta = _sp.norm(v[0]);
1349  norm = beta;
1350  }
1351  else
1352  {
1353  // calc update vector
1354  w = 0;
1355  update(w, i - 1, H, s, v);
1356 
1357  // update x
1358  x += w;
1359 
1360  // r = M^-1 (b - A * x);
1361  // update defect
1362  _A_.applyscaleadd(-1,w, /* => */ b);
1363  // r = M^-1 (b - A * x);
1364  v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b
1365  beta = _sp.norm(v[0]);
1366  norm = beta;
1367 
1368  res.converged = false;
1369  }
1370 
1371  if (_verbose > 1) // print
1372  {
1373  this->printOutput(std::cout,j,norm,norm_old);
1374  }
1375 
1376  norm_old = norm;
1377 
1378  if (norm < reduction * norm_0) {
1379  // fill statistics
1380  res.converged = true;
1381  }
1382 
1383  if (res.converged != true && _verbose > 0)
1384  std::cout << "=== GMRes::restart\n";
1385  }
1386 
1387  _M.post(x); // postprocess preconditioner
1388 
1389  res.iterations = j;
1390  res.reduction = norm / norm_0;
1391  res.conv_rate = pow(res.reduction,1.0/j);
1392  res.elapsed = watch.elapsed();
1393 
1394  if (_verbose>0)
1395  print_result(res);
1396  }
1397  private:
1398 
1399  void
1400  print_result (const InverseOperatorResult & res) const
1401  {
1402  int j = res.iterations>0?res.iterations:1;
1403  std::cout << "=== rate=" << res.conv_rate
1404  << ", T=" << res.elapsed
1405  << ", TIT=" << res.elapsed/j
1406  << ", IT=" << res.iterations
1407  << std::endl;
1408  }
1409 
1410  static void
1411  update(X &x, int k,
1412  std::vector< std::vector<field_type> > & h,
1413  std::vector<field_type> & s, std::vector<F> v)
1414  {
1415  std::vector<field_type> y(s);
1416 
1417  // Backsolve:
1418  for (int i = k; i >= 0; i--) {
1419  y[i] /= h[i][i];
1420  for (int j = i - 1; j >= 0; j--)
1421  y[j] -= h[j][i] * y[i];
1422  }
1423 
1424  for (int j = 0; j <= k; j++)
1425  // x += v[j] * y[j];
1426  x.axpy(y[j],v[j]);
1427  }
1428 
1429  void
1430  generatePlaneRotation(field_type &dx, field_type &dy, field_type &cs, field_type &sn)
1431  {
1432  if (dy == 0.0) {
1433  cs = 1.0;
1434  sn = 0.0;
1435  } else if (std::abs(dy) > std::abs(dx)) {
1436  field_type temp = dx / dy;
1437  sn = 1.0 / std::sqrt( 1.0 + temp*temp );
1438  cs = temp * sn;
1439  } else {
1440  field_type temp = dy / dx;
1441  cs = 1.0 / std::sqrt( 1.0 + temp*temp );
1442  sn = temp * cs;
1443  }
1444  }
1445 
1446 
1447  void
1448  applyPlaneRotation(field_type &dx, field_type &dy, field_type &cs, field_type &sn)
1449  {
1450  field_type temp = cs * dx + sn * dy;
1451  dy = -sn * dx + cs * dy;
1452  dx = temp;
1453  }
1454 
1455  LinearOperator<X,X>& _A_;
1456  Preconditioner<X,X>& _M;
1457  SeqScalarProduct<X> ssp;
1458  ScalarProduct<X>& _sp;
1459  int _restart;
1460  double _reduction;
1461  int _maxit;
1462  int _verbose;
1463  bool _recalc_defect;
1464  };
1465 
1468 } // end namespace
1469 
1470 #endif