Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
solver_gmres.h
1 // ---------------------------------------------------------------------
2 // @f$Id: solver_gmres.h 31349 2013-10-20 19:07:06Z maier @f$
3 //
4 // Copyright (C) 1998 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__solver_gmres_h
18 #define __deal2__solver_gmres_h
19 
20 
21 
22 #include <deal.II/base/config.h>
23 #include <deal.II/base/subscriptor.h>
24 #include <deal.II/base/logstream.h>
25 #include <deal.II/lac/householder.h>
26 #include <deal.II/lac/solver.h>
27 #include <deal.II/lac/solver_control.h>
28 #include <deal.II/lac/full_matrix.h>
29 #include <deal.II/lac/vector.h>
30 
31 #include <vector>
32 #include <cmath>
33 
34 DEAL_II_NAMESPACE_OPEN
35 
38 
39 namespace internal
40 {
44  namespace SolverGMRES
45  {
54  template <class VECTOR>
55  class TmpVectors
56  {
57  public:
61  TmpVectors(const unsigned int max_size,
62  VectorMemory<VECTOR> &vmem);
63 
67  ~TmpVectors();
68 
73  VECTOR &operator[] (const unsigned int i) const;
74 
81  VECTOR &operator() (const unsigned int i,
82  const VECTOR &temp);
83 
84  private:
89 
93  std::vector<VECTOR *> data;
94 
99  unsigned int offset;
100  };
101  }
102 }
103 
155 template <class VECTOR = Vector<double> >
156 class SolverGMRES : public Solver<VECTOR>
157 {
158 public:
163  {
170  AdditionalData (const unsigned int max_n_tmp_vectors = 30,
171  const bool right_preconditioning = false,
172  const bool use_default_residual = true,
173  const bool force_re_orthogonalization = false);
174 
180  unsigned int max_n_tmp_vectors;
181 
190 
195 
203  };
204 
210  const AdditionalData &data=AdditionalData());
211 
217  const AdditionalData &data=AdditionalData());
218 
222  template<class MATRIX, class PRECONDITIONER>
223  void
224  solve (const MATRIX &A,
225  VECTOR &x,
226  const VECTOR &b,
227  const PRECONDITIONER &precondition);
228 
229  DeclException1 (ExcTooFewTmpVectors,
230  int,
231  << "The number of temporary vectors you gave ("
232  << arg1 << ") is too small. It should be at least 10 for "
233  << "any results, and much more for reasonable ones.");
234 
235 protected:
240 
244  virtual double criterion();
245 
252  int col) const;
253 
263  static double
265  const unsigned int dim,
266  const unsigned int accumulated_iterations,
267  VECTOR &vv,
268  Vector<double> &h,
269  bool &re_orthogonalize);
270 
275 
280 
281 
282 private:
287 };
288 
308 template <class VECTOR = Vector<double> >
309 class SolverFGMRES : public Solver<VECTOR>
310 {
311 public:
316  {
322  AdditionalData(const unsigned int max_basis_size = 30,
323  const bool /*use_default_residual*/ = true)
324  :
326  {}
327 
331  unsigned int max_basis_size;
332  };
333 
339  const AdditionalData &data=AdditionalData());
340 
346  const AdditionalData &data=AdditionalData());
347 
351  template<class MATRIX, class PRECONDITIONER>
352  void
353  solve (const MATRIX &A,
354  VECTOR &x,
355  const VECTOR &b,
356  const PRECONDITIONER &precondition);
357 
358 private:
371 };
372 
374 /* --------------------- Inline and template functions ------------------- */
375 
376 
377 #ifndef DOXYGEN
378 namespace internal
379 {
380  namespace SolverGMRES
381  {
382  template <class VECTOR>
383  inline
385  TmpVectors (const unsigned int max_size,
386  VectorMemory<VECTOR> &vmem)
387  :
388  mem(vmem),
389  data (max_size, 0),
390  offset(0)
391  {}
392 
393 
394  template <class VECTOR>
395  inline
397  {
398  for (typename std::vector<VECTOR *>::iterator v = data.begin();
399  v != data.end(); ++v)
400  if (*v != 0)
401  mem.free(*v);
402  }
403 
404 
405  template <class VECTOR>
406  inline VECTOR &
407  TmpVectors<VECTOR>::operator[] (const unsigned int i) const
408  {
409  Assert (i+offset<data.size(),
410  ExcIndexRange(i, -offset, data.size()-offset));
411 
412  Assert (data[i-offset] != 0, ExcNotInitialized());
413  return *data[i-offset];
414  }
415 
416 
417  template <class VECTOR>
418  inline VECTOR &
419  TmpVectors<VECTOR>::operator() (const unsigned int i,
420  const VECTOR &temp)
421  {
422  Assert (i+offset<data.size(),
423  ExcIndexRange(i,-offset, data.size()-offset));
424  if (data[i-offset] == 0)
425  {
426  data[i-offset] = mem.alloc();
427  data[i-offset]->reinit(temp);
428  }
429  return *data[i-offset];
430  }
431  }
432 }
433 
434 
435 
436 template <class VECTOR>
437 inline
439 AdditionalData (const unsigned int max_n_tmp_vectors,
440  const bool right_preconditioning,
441  const bool use_default_residual,
442  const bool force_re_orthogonalization)
443  :
444  max_n_tmp_vectors(max_n_tmp_vectors),
445  right_preconditioning(right_preconditioning),
446  use_default_residual(use_default_residual),
447  force_re_orthogonalization(force_re_orthogonalization)
448 {}
449 
450 
451 template <class VECTOR>
454  const AdditionalData &data)
455  :
456  Solver<VECTOR> (cn,mem),
457  additional_data(data)
458 {}
459 
460 
461 
462 template <class VECTOR>
464  const AdditionalData &data) :
465  Solver<VECTOR> (cn),
466  additional_data(data)
467 {}
468 
469 
470 
471 template <class VECTOR>
472 inline
473 void
475  Vector<double> &b,
476  Vector<double> &ci,
477  Vector<double> &si,
478  int col) const
479 {
480  for (int i=0 ; i<col ; i++)
481  {
482  const double s = si(i);
483  const double c = ci(i);
484  const double dummy = h(i);
485  h(i) = c*dummy + s*h(i+1);
486  h(i+1) = -s*dummy + c*h(i+1);
487  };
488 
489  const double r = 1./std::sqrt(h(col)*h(col) + h(col+1)*h(col+1));
490  si(col) = h(col+1) *r;
491  ci(col) = h(col) *r;
492  h(col) = ci(col)*h(col) + si(col)*h(col+1);
493  b(col+1)= -si(col)*b(col);
494  b(col) *= ci(col);
495 }
496 
497 
498 
499 template <class VECTOR>
500 inline
501 double
503  const unsigned int dim,
504  const unsigned int accumulated_iterations,
505  VECTOR &vv,
506  Vector<double> &h,
507  bool &re_orthogonalize)
508 {
509  const unsigned int inner_iteration = dim - 1;
510 
511  // need initial norm for detection of re-orthogonalization, see below
512  double norm_vv_start = 0;
513  if (re_orthogonalize == false && inner_iteration % 5 == 4)
514  norm_vv_start = vv.l2_norm();
515 
516  // Orthogonalization
517  for (unsigned int i=0 ; i<dim ; ++i)
518  {
519  h(i) = vv * orthogonal_vectors[i];
520  vv.add(-h(i), orthogonal_vectors[i]);
521  };
522 
523  // Re-orthogonalization if loss of orthogonality detected. For the test, use
524  // a strategy discussed in C. T. Kelley, Iterative Methods for Linear and
525  // Nonlinear Equations, SIAM, Philadelphia, 1995: Compare the norm of vv
526  // after orthogonalization with its norm when starting the
527  // orthogonalization. If vv became very small (here: less than the square
528  // root of the machine precision times 10), it is almost in the span of the
529  // previous vectors, which indicates loss of precision.
530  if (re_orthogonalize == false && inner_iteration % 5 == 4)
531  {
532  const double norm_vv = vv.l2_norm();
533  if (norm_vv > 10. * norm_vv_start *
534  std::sqrt(std::numeric_limits<typename VECTOR::value_type>::epsilon()))
535  return norm_vv;
536 
537  else
538  {
539  re_orthogonalize = true;
540  deallog << "Re-orthogonalization enabled at step "
541  << accumulated_iterations << std::endl;
542  }
543  }
544 
545  if (re_orthogonalize == true)
546  for (unsigned int i=0 ; i<dim ; ++i)
547  {
548  double htmp = vv * orthogonal_vectors[i];
549  h(i) += htmp;
550  vv.add(-htmp, orthogonal_vectors[i]);
551  }
552 
553  return vv.l2_norm();
554 }
555 
556 
557 
558 template<class VECTOR>
559 template<class MATRIX, class PRECONDITIONER>
560 void
562  VECTOR &x,
563  const VECTOR &b,
564  const PRECONDITIONER &precondition)
565 {
566  // this code was written a very long time ago by people not associated with
567  // deal.II. we don't make any guarantees to its optimality or that it even
568  // works as expected...
569 
570 //TODO:[?] Check, why there are two different start residuals.
571 //TODO:[GK] Make sure the parameter in the constructor means maximum basis size
572 
573  deallog.push("GMRES");
574  const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors;
575 
576  // Generate an object where basis vectors are stored.
577  internal::SolverGMRES::TmpVectors<VECTOR> tmp_vectors (n_tmp_vectors, this->memory);
578 
579  // number of the present iteration; this
580  // number is not reset to zero upon a
581  // restart
582  unsigned int accumulated_iterations = 0;
583 
584  // matrix used for the orthogonalization process later
585  H.reinit(n_tmp_vectors, n_tmp_vectors-1);
586 
587  // some additional vectors, also used in the orthogonalization
589  gamma(n_tmp_vectors),
590  ci (n_tmp_vectors-1),
591  si (n_tmp_vectors-1),
592  h (n_tmp_vectors-1);
593 
594 
595  unsigned int dim = 0;
596 
598 
599  // switch to determine whether we want a left or a right preconditioner. at
600  // present, left is default, but both ways are implemented
601  const bool left_precondition = !additional_data.right_preconditioning;
602 
603  // Per default the left preconditioned GMRes uses the preconditioned
604  // residual and the right preconditioned GMRes uses the unpreconditioned
605  // residual as stopping criterion.
606  const bool use_default_residual = additional_data.use_default_residual;
607 
608  // define two aliases
609  VECTOR &v = tmp_vectors(0, x);
610  VECTOR &p = tmp_vectors(n_tmp_vectors-1, x);
611 
612  // Following vectors are needed
613  // when not the default residuals
614  // are used as stopping criterion
615  VECTOR *r=0;
616  VECTOR *x_=0;
617  ::Vector<double> *gamma_=0;
618  if (!use_default_residual)
619  {
620  r=this->memory.alloc();
621  x_=this->memory.alloc();
622  r->reinit(x);
623  x_->reinit(x);
624 
625  gamma_ = new ::Vector<double> (gamma.size());
626  }
627 
628  bool re_orthogonalize = additional_data.force_re_orthogonalization;
629 
631  // outer iteration: loop until we either reach convergence or the maximum
632  // number of iterations is exceeded. each cycle of this loop amounts to one
633  // restart
634  do
635  {
636  // reset this vector to the right size
637  h.reinit (n_tmp_vectors-1);
638 
639  if (left_precondition)
640  {
641  A.vmult(p,x);
642  p.sadd(-1.,1.,b);
643  precondition.vmult(v,p);
644  }
645  else
646  {
647  A.vmult(v,x);
648  v.sadd(-1.,1.,b);
649  };
650 
651  double rho = v.l2_norm();
652 
653  // check the residual here as well since it may be that we got the exact
654  // (or an almost exact) solution vector at the outset. if we wouldn't
655  // check here, the next scaling operation would produce garbage
656  if (use_default_residual)
657  {
658  iteration_state = this->control().check (
659  accumulated_iterations, rho);
660 
661  if (iteration_state != SolverControl::iterate)
662  break;
663  }
664  else
665  {
666  deallog << "default_res=" << rho << std::endl;
667 
668  if (left_precondition)
669  {
670  A.vmult(*r,x);
671  r->sadd(-1.,1.,b);
672  }
673  else
674  precondition.vmult(*r,v);
675 
676  double res = r->l2_norm();
677  iteration_state = this->control().check (
678  accumulated_iterations, res);
679 
680  if (iteration_state != SolverControl::iterate)
681  {
682  this->memory.free(r);
683  this->memory.free(x_);
684 
685  delete gamma_;
686  break;
687  }
688  }
689 
690  gamma(0) = rho;
691 
692  v *= 1./rho;
693 
694  // inner iteration doing at most as many steps as there are temporary
695  // vectors. the number of steps actually been done is propagated outside
696  // through the @p dim variable
697  for (unsigned int inner_iteration=0;
698  ((inner_iteration < n_tmp_vectors-2)
699  &&
700  (iteration_state==SolverControl::iterate));
701  ++inner_iteration)
702  {
703  ++accumulated_iterations;
704  // yet another alias
705  VECTOR &vv = tmp_vectors(inner_iteration+1, x);
706 
707  if (left_precondition)
708  {
709  A.vmult(p, tmp_vectors[inner_iteration]);
710  precondition.vmult(vv,p);
711  }
712  else
713  {
714  precondition.vmult(p, tmp_vectors[inner_iteration]);
715  A.vmult(vv,p);
716  };
717 
718  dim = inner_iteration+1;
719 
720  const double s = modified_gram_schmidt(tmp_vectors, dim,
721  accumulated_iterations,
722  vv, h, re_orthogonalize);
723  h(inner_iteration+1) = s;
724 
725  //s=0 is a lucky breakdown, the solver will reach convergence,
726  //but we must not divide by zero here.
727  if (numbers::is_finite(1./s))
728  vv *= 1./s;
729 
730  // Transformation into triagonal structure
731  givens_rotation(h,gamma,ci,si,inner_iteration);
732 
733  // append vector on matrix
734  for (unsigned int i=0; i<dim; ++i)
735  H(i,inner_iteration) = h(i);
736 
737  // default residual
738  rho = std::fabs(gamma(dim));
739 
740  if (use_default_residual)
741  iteration_state = this->control().check (
742  accumulated_iterations, rho);
743  else
744  {
745  deallog << "default_res=" << rho << std::endl;
746 
747  ::Vector<double> h_(dim);
748  *x_=x;
749  *gamma_=gamma;
750  H1.reinit(dim+1,dim);
751 
752  for (unsigned int i=0; i<dim+1; ++i)
753  for (unsigned int j=0; j<dim; ++j)
754  H1(i,j) = H(i,j);
755 
756  H1.backward(h_,*gamma_);
757 
758  if (left_precondition)
759  for (unsigned int i=0 ; i<dim; ++i)
760  x_->add(h_(i), tmp_vectors[i]);
761  else
762  {
763  p = 0.;
764  for (unsigned int i=0; i<dim; ++i)
765  p.add(h_(i), tmp_vectors[i]);
766  precondition.vmult(*r,p);
767  x_->add(1.,*r);
768  };
769  A.vmult(*r,*x_);
770  r->sadd(-1.,1.,b);
771  // Now *r contains the unpreconditioned residual!!
772  if (left_precondition)
773  {
774  const double res=r->l2_norm();
775 
776  iteration_state = this->control().check (
777  accumulated_iterations, res);
778  }
779  else
780  {
781  precondition.vmult(*x_, *r);
782  const double preconditioned_res=x_->l2_norm();
783 
784  iteration_state = this->control().check (
785  accumulated_iterations, preconditioned_res);
786  }
787  }
788  };
789  // end of inner iteration. now calculate the solution from the temporary
790  // vectors
791  h.reinit(dim);
792  H1.reinit(dim+1,dim);
793 
794  for (unsigned int i=0; i<dim+1; ++i)
795  for (unsigned int j=0; j<dim; ++j)
796  H1(i,j) = H(i,j);
797 
798  H1.backward(h,gamma);
799 
800  if (left_precondition)
801  for (unsigned int i=0 ; i<dim; ++i)
802  x.add(h(i), tmp_vectors[i]);
803  else
804  {
805  p = 0.;
806  for (unsigned int i=0; i<dim; ++i)
807  p.add(h(i), tmp_vectors[i]);
808  precondition.vmult(v,p);
809  x.add(1.,v);
810  };
811  // end of outer iteration. restart if no convergence and the number of
812  // iterations is not exceeded
813  }
814  while (iteration_state == SolverControl::iterate);
815 
816  if (!use_default_residual)
817  {
818  this->memory.free(r);
819  this->memory.free(x_);
820 
821  delete gamma_;
822  }
823 
824  deallog.pop();
825  // in case of failure: throw exception
826  if (this->control().last_check() != SolverControl::success)
827  AssertThrow(false, SolverControl::NoConvergence (this->control().last_step(),
828  this->control().last_value()));
829  // otherwise exit as normal
830 }
831 
832 
833 
834 template<class VECTOR>
835 double
837 {
838  // dummy implementation. this function is not needed for the present
839  // implementation of gmres
840  Assert (false, ExcInternalError());
841  return 0;
842 }
843 
844 
845 //----------------------------------------------------------------------//
846 
847 template <class VECTOR>
850  const AdditionalData &data)
851  :
852  Solver<VECTOR> (cn, mem),
853  additional_data(data)
854 {}
855 
856 
857 
858 template <class VECTOR>
860  const AdditionalData &data)
861  :
862  Solver<VECTOR> (cn),
863  additional_data(data)
864 {}
865 
866 
867 
868 template<class VECTOR>
869 template<class MATRIX, class PRECONDITIONER>
870 void
872  const MATRIX &A,
873  VECTOR &x,
874  const VECTOR &b,
875  const PRECONDITIONER &precondition)
876 {
877  deallog.push("FGMRES");
878 
880 
881  const unsigned int basis_size = additional_data.max_basis_size;
882 
883  // Generate an object where basis vectors are stored.
884  typename internal::SolverGMRES::TmpVectors<VECTOR> v (basis_size, this->memory);
885  typename internal::SolverGMRES::TmpVectors<VECTOR> z (basis_size, this->memory);
886 
887  // number of the present iteration; this number is not reset to zero upon a
888  // restart
889  unsigned int accumulated_iterations = 0;
890 
891  // matrix used for the orthogonalization process later
892  H.reinit(basis_size+1, basis_size);
893 
894  // Vectors for projected system
895  Vector<double> projected_rhs;
896  Vector<double> y;
897 
898  // Iteration starts here
899 
900  VECTOR *aux = this->memory.alloc();
901  aux->reinit(x);
902  do
903  {
904  A.vmult(*aux, x);
905  aux->sadd(-1., 1., b);
906 
907  double beta = aux->l2_norm();
908  if (this->control().check(accumulated_iterations,beta)
910  break;
911 
912  H.reinit(basis_size+1, basis_size);
913  double a = beta;
914 
915  for (unsigned int j=0; j<basis_size; ++j)
916  {
917  if (numbers::is_finite(1./a)) // treat lucky breakdown
918  v(j,x).equ(1./a, *aux);
919  else
920  v(j,x) = 0.;
921 
922 
923  precondition.vmult(z(j,x), v[j]);
924  A.vmult(*aux, z[j]);
925 
926  // Gram-Schmidt
927  for (unsigned int i=0; i<=j; ++i)
928  {
929  H(i,j) = *aux * v[i];
930  aux->add(-H(i,j), v[i]);
931  }
932  H(j+1,j) = a = aux->l2_norm();
933 
934  // Compute projected solution
935 
936  if (j>0)
937  {
938  H1.reinit(j+1,j);
939  projected_rhs.reinit(j+1);
940  y.reinit(j);
941  projected_rhs(0) = beta;
942  H1.fill(H);
943  Householder<double> house(H1);
944  double res = house.least_squares(y, projected_rhs);
945  iteration_state = this->control().check(++accumulated_iterations, res);
946  if (iteration_state != SolverControl::iterate)
947  break;
948  }
949  }
950  // Update solution vector
951  for (unsigned int j=0; j<y.size(); ++j)
952  x.add(y(j), z[j]);
953 
954  }
955  while (iteration_state == SolverControl::iterate);
956 
957  this->memory.free(aux);
958 
959  deallog.pop();
960  // in case of failure: throw exception
961  if (this->control().last_check() != SolverControl::success)
962  AssertThrow(false, SolverControl::NoConvergence (this->control().last_step(),
963  this->control().last_value()));
964 }
965 
966 #endif // DOXYGEN
967 
968 DEAL_II_NAMESPACE_CLOSE
969 
970 #endif
void pop()
AdditionalData additional_data
Definition: solver_gmres.h:362
void vmult(VECTOR &u, const VECTOR &v) const
Continue iteration.
virtual double criterion()
FullMatrix< double > H1
Definition: solver_gmres.h:279
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
virtual VECTOR * alloc()=0
FullMatrix< double > H
Definition: solver_gmres.h:274
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
SolverFGMRES(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
void givens_rotation(Vector< double > &h, Vector< double > &b, Vector< double > &ci, Vector< double > &si, int col) const
bool is_finite(const double x)
FullMatrix< double > H
Definition: solver_gmres.h:366
AdditionalData(const unsigned int max_basis_size=30, const bool=true)
Definition: solver_gmres.h:322
Stop iteration, goal reached.
FullMatrix< double > H1
Definition: solver_gmres.h:370
#define Assert(cond, exc)
Definition: exceptions.h:299
AdditionalData additional_data
Definition: solver_gmres.h:239
SolverGMRES(SolverControl &cn, VectorMemory< VECTOR > &mem, const AdditionalData &data=AdditionalData())
VECTOR & operator()(const unsigned int i, const VECTOR &temp)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual void reinit(const size_type N, const bool fast=false)
std::vector< VECTOR * > data
Definition: solver_gmres.h:93
TmpVectors(const unsigned int max_size, VectorMemory< VECTOR > &vmem)
VECTOR & operator[](const unsigned int i) const
virtual void free(const VECTOR *const)=0
static double modified_gram_schmidt(const internal::SolverGMRES::TmpVectors< VECTOR > &orthogonal_vectors, const unsigned int dim, const unsigned int accumulated_iterations, VECTOR &vv, Vector< double > &h, bool &re_orthogonalize)
Definition: solver.h:147
void push(const std::string &text)
std::size_t size() const
::ExceptionBase & ExcNotInitialized()
::ExceptionBase & ExcInternalError()
AdditionalData(const unsigned int max_n_tmp_vectors=30, const bool right_preconditioning=false, const bool use_default_residual=true, const bool force_re_orthogonalization=false)
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
VectorMemory< VECTOR > & mem
Definition: solver_gmres.h:88