Reference documentation for deal.II version 8.4.1
solver_gmres.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2016 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii__solver_gmres_h
17 #define dealii__solver_gmres_h
18 
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/subscriptor.h>
23 #include <deal.II/base/logstream.h>
24 #include <deal.II/lac/householder.h>
25 #include <deal.II/lac/solver.h>
26 #include <deal.II/lac/solver_control.h>
27 #include <deal.II/lac/full_matrix.h>
28 #include <deal.II/lac/lapack_full_matrix.h>
29 #include <deal.II/lac/vector.h>
30 
31 #include <vector>
32 #include <cmath>
33 #include <algorithm>
34 
35 DEAL_II_NAMESPACE_OPEN
36 
39 
40 namespace internal
41 {
45  namespace SolverGMRES
46  {
55  template <typename VectorType>
56  class TmpVectors
57  {
58  public:
63  TmpVectors(const unsigned int max_size,
65 
69  ~TmpVectors();
70 
75  VectorType &operator[] (const unsigned int i) const;
76 
83  VectorType &operator() (const unsigned int i,
84  const VectorType &temp);
85 
86  private:
91 
95  std::vector<VectorType *> data;
96 
101  unsigned int offset;
102  };
103  }
104 }
105 
173 template <class VectorType = Vector<double> >
174 class SolverGMRES : public Solver<VectorType>
175 {
176 public:
181  {
188  explicit
189  AdditionalData (const unsigned int max_n_tmp_vectors = 30,
190  const bool right_preconditioning = false,
191  const bool use_default_residual = true,
192  const bool force_re_orthogonalization = false);
193 
199  AdditionalData (const unsigned int max_n_tmp_vectors,
200  const bool right_preconditioning,
201  const bool use_default_residual,
202  const bool force_re_orthogonalization,
204 
210  unsigned int max_n_tmp_vectors;
211 
220 
225 
233 
244  };
245 
251  const AdditionalData &data=AdditionalData());
252 
258  const AdditionalData &data=AdditionalData());
259 
263  template<typename MatrixType, typename PreconditionerType>
264  void
265  solve (const MatrixType &A,
266  VectorType &x,
267  const VectorType &b,
268  const PreconditionerType &precondition);
269 
276  boost::signals2::connection
277  connect_condition_number_slot(const std_cxx11::function<void (double)> &slot,
278  const bool every_iteration=false);
279 
286  boost::signals2::connection
288  const std_cxx11::function<void (const std::vector<std::complex<double> > &)> &slot,
289  const bool every_iteration=false);
290 
291 
292  DeclException1 (ExcTooFewTmpVectors,
293  int,
294  << "The number of temporary vectors you gave ("
295  << arg1 << ") is too small. It should be at least 10 for "
296  << "any results, and much more for reasonable ones.");
297 
298 protected:
299 
304 
309  boost::signals2::signal<void (double)> condition_number_signal;
310 
315  boost::signals2::signal<void (double)> all_condition_numbers_signal;
316 
321  boost::signals2::signal<void (const std::vector<std::complex<double> > &)> eigenvalues_signal;
322 
327  boost::signals2::signal<void (const std::vector<std::complex<double> > &)> all_eigenvalues_signal;
328 
329 
333  virtual double criterion();
334 
341  int col) const;
342 
352  static double
354  const unsigned int dim,
355  const unsigned int accumulated_iterations,
356  VectorType &vv,
357  Vector<double> &h,
358  bool &re_orthogonalize);
359 
367  static void
369  const FullMatrix<double> &H_orig ,
370  const unsigned int dim,
371  const boost::signals2::signal<void (const std::vector<std::complex<double> > &)> &eigenvalues_signal,
372  const boost::signals2::signal<void(double)> &cond_signal,
373  const bool log_eigenvalues);
374 
379 
384 
385 
386 private:
391 };
392 
414 template <class VectorType = Vector<double> >
415 class SolverFGMRES : public Solver<VectorType>
416 {
417 public:
422  {
426  explicit
427  AdditionalData(const unsigned int max_basis_size = 30,
428  const bool /*use_default_residual*/ = true)
429  :
431  {}
432 
436  unsigned int max_basis_size;
437  };
438 
444  const AdditionalData &data=AdditionalData());
445 
451  const AdditionalData &data=AdditionalData());
452 
456  template<typename MatrixType, typename PreconditionerType>
457  void
458  solve (const MatrixType &A,
459  VectorType &x,
460  const VectorType &b,
461  const PreconditionerType &precondition);
462 
463 private:
464 
469 
474 
479 };
480 
482 /* --------------------- Inline and template functions ------------------- */
483 
484 
485 #ifndef DOXYGEN
486 namespace internal
487 {
488  namespace SolverGMRES
489  {
490  template <class VectorType>
491  inline
493  TmpVectors (const unsigned int max_size,
495  :
496  mem(vmem),
497  data (max_size, 0),
498  offset(0)
499  {}
500 
501 
502  template <class VectorType>
503  inline
505  {
506  for (typename std::vector<VectorType *>::iterator v = data.begin();
507  v != data.end(); ++v)
508  if (*v != 0)
509  mem.free(*v);
510  }
511 
512 
513  template <class VectorType>
514  inline VectorType &
515  TmpVectors<VectorType>::operator[] (const unsigned int i) const
516  {
517  Assert (i+offset<data.size(),
518  ExcIndexRange(i, -offset, data.size()-offset));
519 
520  Assert (data[i-offset] != 0, ExcNotInitialized());
521  return *data[i-offset];
522  }
523 
524 
525  template <class VectorType>
526  inline VectorType &
527  TmpVectors<VectorType>::operator() (const unsigned int i,
528  const VectorType &temp)
529  {
530  Assert (i+offset<data.size(),
531  ExcIndexRange(i,-offset, data.size()-offset));
532  if (data[i-offset] == 0)
533  {
534  data[i-offset] = mem.alloc();
535  data[i-offset]->reinit(temp);
536  }
537  return *data[i-offset];
538  }
539 
540  // A comparator for better printing eigenvalues
541  inline
542  bool complex_less_pred(const std::complex<double> &x,
543  const std::complex<double> &y)
544  {
545  return x.real() < y.real() || (x.real() == y.real() && x.imag() < y.imag());
546  }
547  }
548 }
549 
550 
551 
552 template <class VectorType>
553 inline
555 AdditionalData (const unsigned int max_n_tmp_vectors,
556  const bool right_preconditioning,
557  const bool use_default_residual,
558  const bool force_re_orthogonalization)
559  :
560  max_n_tmp_vectors(max_n_tmp_vectors),
561  right_preconditioning(right_preconditioning),
562  use_default_residual(use_default_residual),
563  force_re_orthogonalization(force_re_orthogonalization),
564  compute_eigenvalues(false)
565 {}
566 
567 
568 
569 template <class VectorType>
570 inline
572 AdditionalData (const unsigned int max_n_tmp_vectors,
573  const bool right_preconditioning,
574  const bool use_default_residual,
575  const bool force_re_orthogonalization,
576  const bool compute_eigenvalues)
577  :
578  max_n_tmp_vectors(max_n_tmp_vectors),
579  right_preconditioning(right_preconditioning),
580  use_default_residual(use_default_residual),
581  force_re_orthogonalization(force_re_orthogonalization),
582  compute_eigenvalues(compute_eigenvalues)
583 {}
584 
585 
586 
587 template <class VectorType>
590  const AdditionalData &data)
591  :
592  Solver<VectorType> (cn,mem),
593  additional_data(data)
594 {}
595 
596 
597 
598 template <class VectorType>
600  const AdditionalData &data) :
601  Solver<VectorType> (cn),
602  additional_data(data)
603 {}
604 
605 
606 
607 template <class VectorType>
608 inline
609 void
611  Vector<double> &b,
612  Vector<double> &ci,
613  Vector<double> &si,
614  int col) const
615 {
616  for (int i=0 ; i<col ; i++)
617  {
618  const double s = si(i);
619  const double c = ci(i);
620  const double dummy = h(i);
621  h(i) = c*dummy + s*h(i+1);
622  h(i+1) = -s*dummy + c*h(i+1);
623  };
624 
625  const double r = 1./std::sqrt(h(col)*h(col) + h(col+1)*h(col+1));
626  si(col) = h(col+1) *r;
627  ci(col) = h(col) *r;
628  h(col) = ci(col)*h(col) + si(col)*h(col+1);
629  b(col+1)= -si(col)*b(col);
630  b(col) *= ci(col);
631 }
632 
633 
634 
635 template <class VectorType>
636 inline
637 double
639 (const internal::SolverGMRES::TmpVectors<VectorType> &orthogonal_vectors,
640  const unsigned int dim,
641  const unsigned int accumulated_iterations,
642  VectorType &vv,
643  Vector<double> &h,
644  bool &re_orthogonalize)
645 {
646  Assert(dim > 0, ExcInternalError());
647  const unsigned int inner_iteration = dim - 1;
648 
649  // need initial norm for detection of re-orthogonalization, see below
650  double norm_vv_start = 0;
651  if (re_orthogonalize == false && inner_iteration % 5 == 4)
652  norm_vv_start = vv.l2_norm();
653 
654  // Orthogonalization
655  h(0) = vv * orthogonal_vectors[0];
656  for (unsigned int i=1 ; i<dim ; ++i)
657  h(i) = vv.add_and_dot(-h(i-1), orthogonal_vectors[i-1], orthogonal_vectors[i]);
658  double norm_vv = std::sqrt(vv.add_and_dot(-h(dim-1), orthogonal_vectors[dim-1], vv));
659 
660  // Re-orthogonalization if loss of orthogonality detected. For the test, use
661  // a strategy discussed in C. T. Kelley, Iterative Methods for Linear and
662  // Nonlinear Equations, SIAM, Philadelphia, 1995: Compare the norm of vv
663  // after orthogonalization with its norm when starting the
664  // orthogonalization. If vv became very small (here: less than the square
665  // root of the machine precision times 10), it is almost in the span of the
666  // previous vectors, which indicates loss of precision.
667  if (re_orthogonalize == false && inner_iteration % 5 == 4)
668  {
669  if (norm_vv > 10. * norm_vv_start *
670  std::sqrt(std::numeric_limits<typename VectorType::value_type>::epsilon()))
671  return norm_vv;
672 
673  else
674  {
675  re_orthogonalize = true;
676  deallog << "Re-orthogonalization enabled at step "
677  << accumulated_iterations << std::endl;
678  }
679  }
680 
681  if (re_orthogonalize == true)
682  {
683  double htmp = vv * orthogonal_vectors[0];
684  h(0) += htmp;
685  for (unsigned int i=1 ; i<dim ; ++i)
686  {
687  htmp = vv.add_and_dot(-htmp, orthogonal_vectors[i-1], orthogonal_vectors[i]);
688  h(i) += htmp;
689  }
690  norm_vv = std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim-1], vv));
691  }
692 
693  return norm_vv;
694 }
695 
696 
697 
698 template<class VectorType>
699 inline void
701 (const FullMatrix<double> &H_orig,
702  const unsigned int dim,
703  const boost::signals2::signal<void (const std::vector<std::complex<double> > &)> &eigenvalues_signal,
704  const boost::signals2::signal<void (double)> &cond_signal,
705  const bool log_eigenvalues)
706 {
707  //Avoid copying the Hessenberg matrix if it isn't needed.
708  if (!eigenvalues_signal.empty() || !cond_signal.empty() || log_eigenvalues )
709  {
710  LAPACKFullMatrix<double> mat(dim,dim);
711  for (unsigned int i=0; i<dim; ++i)
712  for (unsigned int j=0; j<dim; ++j)
713  mat(i,j) = H_orig(i,j);
714  //Avoid computing eigenvalues if they are not needed.
715  if (!eigenvalues_signal.empty() || log_eigenvalues )
716  {
717  //Copy mat so that we can compute svd below. Necessary since
718  //compute_eigenvalues will leave mat in state LAPACKSupport::unusable.
719  LAPACKFullMatrix<double> mat_eig(mat);
720  mat_eig.compute_eigenvalues();
721  std::vector<std::complex<double> > eigenvalues(dim);
722  for (unsigned int i=0; i<mat_eig.n(); ++i)
723  eigenvalues[i] = mat_eig.eigenvalue(i);
724  //Sort eigenvalues for nicer output.
725  std::sort(eigenvalues.begin(), eigenvalues.end(),
726  internal::SolverGMRES::complex_less_pred);
727  eigenvalues_signal(eigenvalues);
728  if (log_eigenvalues)
729  {
730  deallog << "Eigenvalue estimate: ";
731  for (unsigned int i=0; i<mat_eig.n(); ++i)
732  deallog << ' ' << eigenvalues[i];
733  deallog << std::endl;
734  }
735  }
736  //Calculate condition number, avoid calculating the svd if a slot
737  //isn't connected. Need at least a 2-by-2 matrix to do the estimate.
738  if (!cond_signal.empty() && (mat.n()>1))
739  {
740  mat.compute_svd();
741  double condition_number=mat.singular_value(0)/mat.singular_value(mat.n()-1);
742  cond_signal(condition_number);
743  }
744  }
745 }
746 
747 
748 
749 template<class VectorType>
750 template<typename MatrixType, typename PreconditionerType>
751 void
752 SolverGMRES<VectorType>::solve (const MatrixType &A,
753  VectorType &x,
754  const VectorType &b,
755  const PreconditionerType &precondition)
756 {
757  // this code was written a very long time ago by people not associated with
758  // deal.II. we don't make any guarantees to its optimality or that it even
759  // works as expected...
760 
761 //TODO:[?] Check, why there are two different start residuals.
762 //TODO:[GK] Make sure the parameter in the constructor means maximum basis size
763 
764  deallog.push("GMRES");
765  const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors;
766 
767  // Generate an object where basis vectors are stored.
768  internal::SolverGMRES::TmpVectors<VectorType> tmp_vectors (n_tmp_vectors, this->memory);
769 
770  // number of the present iteration; this
771  // number is not reset to zero upon a
772  // restart
773  unsigned int accumulated_iterations = 0;
774 
775  const bool do_eigenvalues=
776  !condition_number_signal.empty()
777  |!all_condition_numbers_signal.empty()
778  |!eigenvalues_signal.empty()
779  |!all_eigenvalues_signal.empty()
780  |additional_data.compute_eigenvalues;
781  // for eigenvalue computation, need to collect the Hessenberg matrix (before
782  // applying Givens rotations)
783  FullMatrix<double> H_orig;
784  if (do_eigenvalues)
785  H_orig.reinit(n_tmp_vectors, n_tmp_vectors-1);
786 
787  // matrix used for the orthogonalization process later
788  H.reinit(n_tmp_vectors, n_tmp_vectors-1);
789 
790  // some additional vectors, also used in the orthogonalization
792  gamma(n_tmp_vectors),
793  ci (n_tmp_vectors-1),
794  si (n_tmp_vectors-1),
795  h (n_tmp_vectors-1);
796 
797 
798  unsigned int dim = 0;
799 
801  double last_res = -std::numeric_limits<double>::max();
802 
803  // switch to determine whether we want a left or a right preconditioner. at
804  // present, left is default, but both ways are implemented
805  const bool left_precondition = !additional_data.right_preconditioning;
806 
807  // Per default the left preconditioned GMRes uses the preconditioned
808  // residual and the right preconditioned GMRes uses the unpreconditioned
809  // residual as stopping criterion.
810  const bool use_default_residual = additional_data.use_default_residual;
811 
812  // define two aliases
813  VectorType &v = tmp_vectors(0, x);
814  VectorType &p = tmp_vectors(n_tmp_vectors-1, x);
815 
816  // Following vectors are needed
817  // when not the default residuals
818  // are used as stopping criterion
819  VectorType *r=0;
820  VectorType *x_=0;
821  ::Vector<double> *gamma_=0;
822  if (!use_default_residual)
823  {
824  r=this->memory.alloc();
825  x_=this->memory.alloc();
826  r->reinit(x);
827  x_->reinit(x);
828 
829  gamma_ = new ::Vector<double> (gamma.size());
830  }
831 
832  bool re_orthogonalize = additional_data.force_re_orthogonalization;
833 
835  // outer iteration: loop until we either reach convergence or the maximum
836  // number of iterations is exceeded. each cycle of this loop amounts to one
837  // restart
838  do
839  {
840  // reset this vector to the right size
841  h.reinit (n_tmp_vectors-1);
842 
843  if (left_precondition)
844  {
845  A.vmult(p,x);
846  p.sadd(-1.,1.,b);
847  precondition.vmult(v,p);
848  }
849  else
850  {
851  A.vmult(v,x);
852  v.sadd(-1.,1.,b);
853  };
854 
855  double rho = v.l2_norm();
856 
857  // check the residual here as well since it may be that we got the exact
858  // (or an almost exact) solution vector at the outset. if we wouldn't
859  // check here, the next scaling operation would produce garbage
860  if (use_default_residual)
861  {
862  last_res = rho;
863  iteration_state = this->iteration_status (accumulated_iterations, rho, x);
864 
865  if (iteration_state != SolverControl::iterate)
866  break;
867  }
868  else
869  {
870  deallog << "default_res=" << rho << std::endl;
871 
872  if (left_precondition)
873  {
874  A.vmult(*r,x);
875  r->sadd(-1.,1.,b);
876  }
877  else
878  precondition.vmult(*r,v);
879 
880  double res = r->l2_norm();
881  last_res = res;
882  iteration_state = this->iteration_status (accumulated_iterations, res, x);
883 
884  if (iteration_state != SolverControl::iterate)
885  break;
886  }
887 
888  gamma(0) = rho;
889 
890  v *= 1./rho;
891 
892  // inner iteration doing at most as many steps as there are temporary
893  // vectors. the number of steps actually been done is propagated outside
894  // through the @p dim variable
895  for (unsigned int inner_iteration=0;
896  ((inner_iteration < n_tmp_vectors-2)
897  &&
898  (iteration_state==SolverControl::iterate));
899  ++inner_iteration)
900  {
901  ++accumulated_iterations;
902  // yet another alias
903  VectorType &vv = tmp_vectors(inner_iteration+1, x);
904 
905  if (left_precondition)
906  {
907  A.vmult(p, tmp_vectors[inner_iteration]);
908  precondition.vmult(vv,p);
909  }
910  else
911  {
912  precondition.vmult(p, tmp_vectors[inner_iteration]);
913  A.vmult(vv,p);
914  }
915 
916  dim = inner_iteration+1;
917 
918  const double s = modified_gram_schmidt(tmp_vectors, dim,
919  accumulated_iterations,
920  vv, h, re_orthogonalize);
921  h(inner_iteration+1) = s;
922 
923  //s=0 is a lucky breakdown, the solver will reach convergence,
924  //but we must not divide by zero here.
925  if (s != 0)
926  vv *= 1./s;
927 
928  // for eigenvalues, get the resulting coefficients from the
929  // orthogonalization process
930  if (do_eigenvalues)
931  for (unsigned int i=0; i<dim+1; ++i)
932  H_orig(i,inner_iteration) = h(i);
933 
934  // Transformation into tridiagonal structure
935  givens_rotation(h,gamma,ci,si,inner_iteration);
936 
937  // append vector on matrix
938  for (unsigned int i=0; i<dim; ++i)
939  H(i,inner_iteration) = h(i);
940 
941  // default residual
942  rho = std::fabs(gamma(dim));
943 
944  if (use_default_residual)
945  {
946  last_res = rho;
947  iteration_state = this->iteration_status (accumulated_iterations, rho, x);
948  }
949  else
950  {
951  deallog << "default_res=" << rho << std::endl;
952 
953  ::Vector<double> h_(dim);
954  *x_=x;
955  *gamma_=gamma;
956  H1.reinit(dim+1,dim);
957 
958  for (unsigned int i=0; i<dim+1; ++i)
959  for (unsigned int j=0; j<dim; ++j)
960  H1(i,j) = H(i,j);
961 
962  H1.backward(h_,*gamma_);
963 
964  if (left_precondition)
965  for (unsigned int i=0 ; i<dim; ++i)
966  x_->add(h_(i), tmp_vectors[i]);
967  else
968  {
969  p = 0.;
970  for (unsigned int i=0; i<dim; ++i)
971  p.add(h_(i), tmp_vectors[i]);
972  precondition.vmult(*r,p);
973  x_->add(1.,*r);
974  };
975  A.vmult(*r,*x_);
976  r->sadd(-1.,1.,b);
977  // Now *r contains the unpreconditioned residual!!
978  if (left_precondition)
979  {
980  const double res=r->l2_norm();
981  last_res = res;
982 
983  iteration_state = this->iteration_status (accumulated_iterations, res, x);
984  }
985  else
986  {
987  precondition.vmult(*x_, *r);
988  const double preconditioned_res=x_->l2_norm();
989  last_res = preconditioned_res;
990 
991  iteration_state = this->iteration_status (accumulated_iterations,
992  preconditioned_res, x);
993  }
994  }
995  };
996  // end of inner iteration. now calculate the solution from the temporary
997  // vectors
998  h.reinit(dim);
999  H1.reinit(dim+1,dim);
1000 
1001  for (unsigned int i=0; i<dim+1; ++i)
1002  for (unsigned int j=0; j<dim; ++j)
1003  H1(i,j) = H(i,j);
1004 
1005  compute_eigs_and_cond(H_orig,dim,all_eigenvalues_signal,
1006  all_condition_numbers_signal,
1007  additional_data.compute_eigenvalues);
1008 
1009  H1.backward(h,gamma);
1010 
1011  if (left_precondition)
1012  for (unsigned int i=0 ; i<dim; ++i)
1013  x.add(h(i), tmp_vectors[i]);
1014  else
1015  {
1016  p = 0.;
1017  for (unsigned int i=0; i<dim; ++i)
1018  p.add(h(i), tmp_vectors[i]);
1019  precondition.vmult(v,p);
1020  x.add(1.,v);
1021  };
1022  // end of outer iteration. restart if no convergence and the number of
1023  // iterations is not exceeded
1024  }
1025  while (iteration_state == SolverControl::iterate);
1026 
1027  compute_eigs_and_cond(H_orig,dim,eigenvalues_signal,condition_number_signal,
1028  false);
1029  if (!use_default_residual)
1030  {
1031  this->memory.free(r);
1032  this->memory.free(x_);
1033 
1034  delete gamma_;
1035  }
1036 
1037  deallog.pop();
1038 
1039  // in case of failure: throw exception
1040  AssertThrow(iteration_state == SolverControl::success,
1041  SolverControl::NoConvergence (accumulated_iterations,
1042  last_res));
1043 }
1044 
1045 
1046 
1047 template<class VectorType>
1048 boost::signals2::connection
1050 (const std_cxx11::function<void(double)> &slot,
1051  const bool every_iteration)
1052 {
1053  if (every_iteration)
1054  {
1055  return all_condition_numbers_signal.connect(slot);
1056  }
1057  else
1058  {
1059  return condition_number_signal.connect(slot);
1060  }
1061 }
1062 
1063 
1064 
1065 template<class VectorType>
1066 boost::signals2::connection
1068 (const std_cxx11::function<void (const std::vector<std::complex<double> > &)> &slot,
1069  const bool every_iteration)
1070 {
1071  if (every_iteration)
1072  {
1073  return all_eigenvalues_signal.connect(slot);
1074  }
1075  else
1076  {
1077  return eigenvalues_signal.connect(slot);
1078  }
1079 }
1080 
1081 
1082 
1083 template<class VectorType>
1084 double
1086 {
1087  // dummy implementation. this function is not needed for the present
1088  // implementation of gmres
1089  Assert (false, ExcInternalError());
1090  return 0;
1091 }
1092 
1093 
1094 //----------------------------------------------------------------------//
1095 
1096 template <class VectorType>
1099  const AdditionalData &data)
1100  :
1101  Solver<VectorType> (cn, mem),
1102  additional_data(data)
1103 {}
1104 
1105 
1106 
1107 template <class VectorType>
1109  const AdditionalData &data)
1110  :
1111  Solver<VectorType> (cn),
1112  additional_data(data)
1113 {}
1114 
1115 
1116 
1117 template<class VectorType>
1118 template<typename MatrixType, typename PreconditionerType>
1119 void
1120 SolverFGMRES<VectorType>::solve (const MatrixType &A,
1121  VectorType &x,
1122  const VectorType &b,
1123  const PreconditionerType &precondition)
1124 {
1125  deallog.push("FGMRES");
1126 
1127  SolverControl::State iteration_state = SolverControl::iterate;
1128 
1129  const unsigned int basis_size = additional_data.max_basis_size;
1130 
1131  // Generate an object where basis vectors are stored.
1132  typename internal::SolverGMRES::TmpVectors<VectorType> v (basis_size, this->memory);
1133  typename internal::SolverGMRES::TmpVectors<VectorType> z (basis_size, this->memory);
1134 
1135  // number of the present iteration; this number is not reset to zero upon a
1136  // restart
1137  unsigned int accumulated_iterations = 0;
1138 
1139  // matrix used for the orthogonalization process later
1140  H.reinit(basis_size+1, basis_size);
1141 
1142  // Vectors for projected system
1143  Vector<double> projected_rhs;
1144  Vector<double> y;
1145 
1146  // Iteration starts here
1147  double res = -std::numeric_limits<double>::max();
1148 
1149  VectorType *aux = this->memory.alloc();
1150  aux->reinit(x);
1151  do
1152  {
1153  A.vmult(*aux, x);
1154  aux->sadd(-1., 1., b);
1155 
1156  double beta = aux->l2_norm();
1157  res = beta;
1158  iteration_state = this->iteration_status(accumulated_iterations, res, x);
1159  if (iteration_state == SolverControl::success)
1160  break;
1161 
1162  H.reinit(basis_size+1, basis_size);
1163  double a = beta;
1164 
1165  for (unsigned int j=0; j<basis_size; ++j)
1166  {
1167  if (a != 0) // treat lucky breakdown
1168  v(j,x).equ(1./a, *aux);
1169  else
1170  v(j,x) = 0.;
1171 
1172 
1173  precondition.vmult(z(j,x), v[j]);
1174  A.vmult(*aux, z[j]);
1175 
1176  // Gram-Schmidt
1177  H(0,j) = *aux * v[0];
1178  for (unsigned int i=1; i<=j; ++i)
1179  H(i,j) = aux->add_and_dot(-H(i-1,j), v[i-1], v[i]);
1180  H(j+1,j) = a = std::sqrt(aux->add_and_dot(-H(j,j), v[j], *aux));
1181 
1182  // Compute projected solution
1183 
1184  if (j>0)
1185  {
1186  H1.reinit(j+1,j);
1187  projected_rhs.reinit(j+1);
1188  y.reinit(j);
1189  projected_rhs(0) = beta;
1190  H1.fill(H);
1191 
1192  // check convergence. note that the vector 'x' we pass to the
1193  // criterion is not the final solution we compute if we
1194  // decide to jump out of the iteration (we update 'x' again
1195  // right after the current loop)
1196  Householder<double> house(H1);
1197  res = house.least_squares(y, projected_rhs);
1198  iteration_state = this->iteration_status(++accumulated_iterations, res, x);
1199  if (iteration_state != SolverControl::iterate)
1200  break;
1201  }
1202  }
1203 
1204  // Update solution vector
1205  for (unsigned int j=0; j<y.size(); ++j)
1206  x.add(y(j), z[j]);
1207  }
1208  while (iteration_state == SolverControl::iterate);
1209 
1210  this->memory.free(aux);
1211 
1212  deallog.pop();
1213  // in case of failure: throw exception
1214  if (iteration_state != SolverControl::success)
1215  AssertThrow(false, SolverControl::NoConvergence (accumulated_iterations,
1216  res));
1217 }
1218 
1219 #endif // DOXYGEN
1220 
1221 DEAL_II_NAMESPACE_CLOSE
1222 
1223 #endif
VectorMemory< VectorType > & mem
Definition: solver_gmres.h:90
SolverFGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
void pop()
Definition: logstream.cc:300
Number add_and_dot(const Number a, const Vector< Number > &V, const Vector< Number > &W)
virtual VectorType * alloc()=0
boost::signals2::signal< void(const std::vector< std::complex< double > > &)> all_eigenvalues_signal
Definition: solver_gmres.h:327
Continue iteration.
FullMatrix< double > H1
Definition: solver_gmres.h:478
AdditionalData additional_data
Definition: solver_gmres.h:303
#define AssertThrow(cond, exc)
Definition: exceptions.h:358
AdditionalData additional_data
Definition: solver_gmres.h:468
FullMatrix< double > H1
Definition: solver_gmres.h:383
virtual void free(const VectorType *const)=0
FullMatrix< double > H
Definition: solver_gmres.h:378
#define DEAL_II_DEPRECATED
Definition: config.h:88
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &precondition)
VectorType & operator()(const unsigned int i, const VectorType &temp)
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 reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:294
VectorType & operator[](const unsigned int i) const
static void compute_eigs_and_cond(const FullMatrix< double > &H_orig, const unsigned int dim, const boost::signals2::signal< void(const std::vector< std::complex< double > > &)> &eigenvalues_signal, const boost::signals2::signal< void(double)> &cond_signal, const bool log_eigenvalues)
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &precondition)
std::size_t size() const
std::vector< VectorType * > data
Definition: solver_gmres.h:95
void givens_rotation(Vector< double > &h, Vector< double > &b, Vector< double > &ci, Vector< double > &si, int col) const
static double modified_gram_schmidt(const internal::SolverGMRES::TmpVectors< VectorType > &orthogonal_vectors, const unsigned int dim, const unsigned int accumulated_iterations, VectorType &vv, Vector< double > &h, bool &re_orthogonalize)
Definition: solver.h:325
boost::signals2::connection connect_condition_number_slot(const std_cxx11::function< void(double)> &slot, const bool every_iteration=false)
virtual double criterion()
void push(const std::string &text)
Definition: logstream.cc:288
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_gmres.h:309
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
boost::signals2::connection connect_eigenvalues_slot(const std_cxx11::function< void(const std::vector< std::complex< double > > &)> &slot, const bool every_iteration=false)
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
FullMatrix< double > H
Definition: solver_gmres.h:473
AdditionalData(const unsigned int max_basis_size=30, const bool=true)
Definition: solver_gmres.h:427
SolverGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
boost::signals2::signal< void(const std::vector< std::complex< double > > &)> eigenvalues_signal
Definition: solver_gmres.h:321
boost::signals2::signal< void(double)> all_condition_numbers_signal
Definition: solver_gmres.h:315
Stop iteration, goal reached.