Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
solver_gmres.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 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.md at
12 // the top level directory of deal.II.
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 
23 #include <deal.II/base/logstream.h>
24 #include <deal.II/base/std_cxx14/memory.h>
25 #include <deal.II/base/subscriptor.h>
26 
27 #include <deal.II/lac/full_matrix.h>
28 #include <deal.II/lac/householder.h>
29 #include <deal.II/lac/lapack_full_matrix.h>
30 #include <deal.II/lac/solver.h>
31 #include <deal.II/lac/solver_control.h>
32 #include <deal.II/lac/vector.h>
33 
34 #include <algorithm>
35 #include <cmath>
36 #include <vector>
37 
38 DEAL_II_NAMESPACE_OPEN
39 
42 
43 namespace internal
44 {
48  namespace SolverGMRESImplementation
49  {
58  template <typename VectorType>
59  class TmpVectors
60  {
61  public:
66  TmpVectors(const unsigned int max_size, VectorMemory<VectorType> &vmem);
67 
71  ~TmpVectors() = default;
72 
77  VectorType &operator[](const unsigned int i) const;
78 
85  VectorType &
86  operator()(const unsigned int i, const VectorType &temp);
87 
92  unsigned int
93  size() const;
94 
95 
96  private:
101 
105  std::vector<typename VectorMemory<VectorType>::Pointer> data;
106  };
107  } // namespace SolverGMRESImplementation
108 } // namespace internal
109 
177 template <class VectorType = Vector<double>>
178 class SolverGMRES : public SolverBase<VectorType>
179 {
180 public:
185  {
192  explicit AdditionalData(const unsigned int max_n_tmp_vectors = 30,
193  const bool right_preconditioning = false,
194  const bool use_default_residual = true,
195  const bool force_re_orthogonalization = false);
196 
203  unsigned int max_n_tmp_vectors;
204 
213 
218 
226  };
227 
233  const AdditionalData & data = AdditionalData());
234 
240 
244  SolverGMRES(const SolverGMRES<VectorType> &) = delete;
245 
249  template <typename MatrixType, typename PreconditionerType>
250  void
251  solve(const MatrixType & A,
252  VectorType & x,
253  const VectorType & b,
254  const PreconditionerType &preconditioner);
255 
262  boost::signals2::connection
263  connect_condition_number_slot(const std::function<void(double)> &slot,
264  const bool every_iteration = false);
265 
272  boost::signals2::connection
274  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
275  const bool every_iteration = false);
276 
284  boost::signals2::connection
286  const std::function<void(const FullMatrix<double> &)> &slot,
287  const bool every_iteration = true);
288 
295  boost::signals2::connection
297  const std::function<
299  &slot);
300 
301 
306  boost::signals2::connection
307  connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
308 
309 
311  int,
312  << "The number of temporary vectors you gave (" << arg1
313  << ") is too small. It should be at least 10 for "
314  << "any results, and much more for reasonable ones.");
315 
316 protected:
321 
326  boost::signals2::signal<void(double)> condition_number_signal;
327 
332  boost::signals2::signal<void(double)> all_condition_numbers_signal;
333 
338  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
340 
345  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
347 
352  boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
353 
358  boost::signals2::signal<void(const FullMatrix<double> &)>
360 
365  boost::signals2::signal<void(
368 
373  boost::signals2::signal<void(int)> re_orthogonalize_signal;
374 
378  virtual double
379  criterion();
380 
385  void
387  Vector<double> &b,
388  Vector<double> &ci,
389  Vector<double> &si,
390  int col) const;
391 
402  static double
405  & orthogonal_vectors,
406  const unsigned int dim,
407  const unsigned int accumulated_iterations,
408  VectorType & vv,
409  Vector<double> & h,
410  bool & re_orthogonalize,
411  const boost::signals2::signal<void(int)> &re_orthogonalize_signal =
412  boost::signals2::signal<void(int)>());
413 
420  static void
422  const FullMatrix<double> &H_orig,
423  const unsigned int dim,
424  const boost::signals2::signal<
425  void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
426  const boost::signals2::signal<void(const FullMatrix<double> &)>
428  const boost::signals2::signal<void(double)> &cond_signal);
429 
434 
439 };
440 
459 template <class VectorType = Vector<double>>
460 class SolverFGMRES : public SolverBase<VectorType>
461 {
462 public:
467  {
471  explicit AdditionalData(const unsigned int max_basis_size = 30)
473  {}
474 
479  DEAL_II_DEPRECATED
480  AdditionalData(const unsigned int max_basis_size,
481  const bool use_default_residual)
483  {
484  (void)use_default_residual;
485  }
486 
490  unsigned int max_basis_size;
491  };
492 
498  const AdditionalData & data = AdditionalData());
499 
505  const AdditionalData &data = AdditionalData());
506 
510  template <typename MatrixType, typename PreconditionerType>
511  void
512  solve(const MatrixType & A,
513  VectorType & x,
514  const VectorType & b,
515  const PreconditionerType &preconditioner);
516 
517 private:
522 
527 
532 };
533 
535 /* --------------------- Inline and template functions ------------------- */
536 
537 
538 #ifndef DOXYGEN
539 namespace internal
540 {
541  namespace SolverGMRESImplementation
542  {
543  template <class VectorType>
544  inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
546  : mem(vmem)
547  , data(max_size)
548  {}
549 
550 
551 
552  template <class VectorType>
553  inline VectorType &TmpVectors<VectorType>::
554  operator[](const unsigned int i) const
555  {
556  Assert(i < data.size(), ExcIndexRange(i, 0, data.size()));
557 
558  Assert(data[i] != nullptr, ExcNotInitialized());
559  return *data[i];
560  }
561 
562 
563 
564  template <class VectorType>
565  inline VectorType &
566  TmpVectors<VectorType>::operator()(const unsigned int i,
567  const VectorType & temp)
568  {
569  Assert(i < data.size(), ExcIndexRange(i, 0, data.size()));
570  if (data[i] == nullptr)
571  {
572  data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
573  data[i]->reinit(temp);
574  }
575  return *data[i];
576  }
577 
578 
579 
580  template <class VectorType>
581  unsigned int
583  {
584  return (data.size() > 0 ? data.size() - 1 : 0);
585  }
586 
587 
588 
589  // A comparator for better printing eigenvalues
590  inline bool
591  complex_less_pred(const std::complex<double> &x,
592  const std::complex<double> &y)
593  {
594  return x.real() < y.real() ||
595  (x.real() == y.real() && x.imag() < y.imag());
596  }
597  } // namespace SolverGMRESImplementation
598 } // namespace internal
599 
600 
601 
602 template <class VectorType>
604  const unsigned int max_n_tmp_vectors,
605  const bool right_preconditioning,
606  const bool use_default_residual,
607  const bool force_re_orthogonalization)
608  : max_n_tmp_vectors(max_n_tmp_vectors)
609  , right_preconditioning(right_preconditioning)
610  , use_default_residual(use_default_residual)
611  , force_re_orthogonalization(force_re_orthogonalization)
612 {
613  Assert(3 <= max_n_tmp_vectors,
614  ExcMessage("SolverGMRES needs at least three "
615  "temporary vectors."));
616 }
617 
618 
619 
620 template <class VectorType>
623  const AdditionalData & data)
624  : SolverBase<VectorType>(cn, mem)
625  , additional_data(data)
626 {}
627 
628 
629 
630 template <class VectorType>
632  const AdditionalData &data)
633  : SolverBase<VectorType>(cn)
634  , additional_data(data)
635 {}
636 
637 
638 
639 template <class VectorType>
640 inline void
642  Vector<double> &b,
643  Vector<double> &ci,
644  Vector<double> &si,
645  int col) const
646 {
647  for (int i = 0; i < col; i++)
648  {
649  const double s = si(i);
650  const double c = ci(i);
651  const double dummy = h(i);
652  h(i) = c * dummy + s * h(i + 1);
653  h(i + 1) = -s * dummy + c * h(i + 1);
654  };
655 
656  const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
657  si(col) = h(col + 1) * r;
658  ci(col) = h(col) * r;
659  h(col) = ci(col) * h(col) + si(col) * h(col + 1);
660  b(col + 1) = -si(col) * b(col);
661  b(col) *= ci(col);
662 }
663 
664 
665 
666 template <class VectorType>
667 inline double
670  & orthogonal_vectors,
671  const unsigned int dim,
672  const unsigned int accumulated_iterations,
673  VectorType & vv,
674  Vector<double> & h,
675  bool & reorthogonalize,
676  const boost::signals2::signal<void(int)> &reorthogonalize_signal)
677 {
678  Assert(dim > 0, ExcInternalError());
679  const unsigned int inner_iteration = dim - 1;
680 
681  // need initial norm for detection of re-orthogonalization, see below
682  double norm_vv_start = 0;
683  const bool consider_reorthogonalize =
684  (reorthogonalize == false) && (inner_iteration % 5 == 4);
685  if (consider_reorthogonalize)
686  norm_vv_start = vv.l2_norm();
687 
688  // Orthogonalization
689  h(0) = vv * orthogonal_vectors[0];
690  for (unsigned int i = 1; i < dim; ++i)
691  h(i) = vv.add_and_dot(-h(i - 1),
692  orthogonal_vectors[i - 1],
693  orthogonal_vectors[i]);
694  double norm_vv =
695  std::sqrt(vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv));
696 
697  // Re-orthogonalization if loss of orthogonality detected. For the test, use
698  // a strategy discussed in C. T. Kelley, Iterative Methods for Linear and
699  // Nonlinear Equations, SIAM, Philadelphia, 1995: Compare the norm of vv
700  // after orthogonalization with its norm when starting the
701  // orthogonalization. If vv became very small (here: less than the square
702  // root of the machine precision times 10), it is almost in the span of the
703  // previous vectors, which indicates loss of precision.
704  if (consider_reorthogonalize)
705  {
706  if (norm_vv >
707  10. * norm_vv_start *
708  std::sqrt(
709  std::numeric_limits<typename VectorType::value_type>::epsilon()))
710  return norm_vv;
711 
712  else
713  {
714  reorthogonalize = true;
715  if (!reorthogonalize_signal.empty())
716  reorthogonalize_signal(accumulated_iterations);
717  }
718  }
719 
720  if (reorthogonalize == true)
721  {
722  double htmp = vv * orthogonal_vectors[0];
723  h(0) += htmp;
724  for (unsigned int i = 1; i < dim; ++i)
725  {
726  htmp = vv.add_and_dot(-htmp,
727  orthogonal_vectors[i - 1],
728  orthogonal_vectors[i]);
729  h(i) += htmp;
730  }
731  norm_vv =
732  std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
733  }
734 
735  return norm_vv;
736 }
737 
738 
739 
740 template <class VectorType>
741 inline void
743  const FullMatrix<double> &H_orig,
744  const unsigned int dim,
745  const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
746  &eigenvalues_signal,
747  const boost::signals2::signal<void(const FullMatrix<double> &)>
748  & hessenberg_signal,
749  const boost::signals2::signal<void(double)> &cond_signal)
750 {
751  // Avoid copying the Hessenberg matrix if it isn't needed.
752  if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
753  !cond_signal.empty()) &&
754  dim > 0)
755  {
756  LAPACKFullMatrix<double> mat(dim, dim);
757  for (unsigned int i = 0; i < dim; ++i)
758  for (unsigned int j = 0; j < dim; ++j)
759  mat(i, j) = H_orig(i, j);
760  hessenberg_signal(H_orig);
761  // Avoid computing eigenvalues if they are not needed.
762  if (!eigenvalues_signal.empty())
763  {
764  // Copy mat so that we can compute svd below. Necessary since
765  // compute_eigenvalues will leave mat in state
766  // LAPACKSupport::unusable.
767  LAPACKFullMatrix<double> mat_eig(mat);
768  mat_eig.compute_eigenvalues();
769  std::vector<std::complex<double>> eigenvalues(dim);
770  for (unsigned int i = 0; i < mat_eig.n(); ++i)
771  eigenvalues[i] = mat_eig.eigenvalue(i);
772  // Sort eigenvalues for nicer output.
773  std::sort(eigenvalues.begin(),
774  eigenvalues.end(),
775  internal::SolverGMRESImplementation::complex_less_pred);
776  eigenvalues_signal(eigenvalues);
777  }
778  // Calculate condition number, avoid calculating the svd if a slot
779  // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
780  if (!cond_signal.empty() && (mat.n() > 1))
781  {
782  mat.compute_svd();
783  double condition_number =
784  mat.singular_value(0) / mat.singular_value(mat.n() - 1);
785  cond_signal(condition_number);
786  }
787  }
788 }
789 
790 
791 
792 template <class VectorType>
793 template <typename MatrixType, typename PreconditionerType>
794 void
795 SolverGMRES<VectorType>::solve(const MatrixType & A,
796  VectorType & x,
797  const VectorType & b,
798  const PreconditionerType &preconditioner)
799 {
800  // TODO:[?] Check, why there are two different start residuals.
801  // TODO:[GK] Make sure the parameter in the constructor means maximum basis
802  // size
803 
804  LogStream::Prefix prefix("GMRES");
805 
806  // extra call to std::max to placate static analyzers: coverity rightfully
807  // complains that data.max_n_tmp_vectors - 2 may overflow
808  const unsigned int n_tmp_vectors =
809  std::max(additional_data.max_n_tmp_vectors, 3u);
810 
811  // Generate an object where basis vectors are stored.
813  n_tmp_vectors, this->memory);
814 
815  // number of the present iteration; this
816  // number is not reset to zero upon a
817  // restart
818  unsigned int accumulated_iterations = 0;
819 
820  const bool do_eigenvalues =
821  !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
822  !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty() ||
823  !hessenberg_signal.empty() || !all_hessenberg_signal.empty();
824  // for eigenvalue computation, need to collect the Hessenberg matrix (before
825  // applying Givens rotations)
826  FullMatrix<double> H_orig;
827  if (do_eigenvalues)
828  H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
829 
830  // matrix used for the orthogonalization process later
831  H.reinit(n_tmp_vectors, n_tmp_vectors - 1);
832 
833  // some additional vectors, also used in the orthogonalization
834  ::Vector<double> gamma(n_tmp_vectors), ci(n_tmp_vectors - 1),
835  si(n_tmp_vectors - 1), h(n_tmp_vectors - 1);
836 
837 
838  unsigned int dim = 0;
839 
841  double last_res = -std::numeric_limits<double>::max();
842 
843  // switch to determine whether we want a left or a right preconditioner. at
844  // present, left is default, but both ways are implemented
845  const bool left_precondition = !additional_data.right_preconditioning;
846 
847  // Per default the left preconditioned GMRes uses the preconditioned
848  // residual and the right preconditioned GMRes uses the unpreconditioned
849  // residual as stopping criterion.
850  const bool use_default_residual = additional_data.use_default_residual;
851 
852  // define two aliases
853  VectorType &v = tmp_vectors(0, x);
854  VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
855 
856  // Following vectors are needed when we are not using the default residuals
857  // as stopping criterion
860  std::unique_ptr<::Vector<double>> gamma_;
861  if (!use_default_residual)
862  {
863  r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
864  x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
865  r->reinit(x);
866  x_->reinit(x);
867 
868  gamma_ = std_cxx14::make_unique<::Vector<double>>(gamma.size());
869  }
870 
871  bool re_orthogonalize = additional_data.force_re_orthogonalization;
872 
874  // outer iteration: loop until we either reach convergence or the maximum
875  // number of iterations is exceeded. each cycle of this loop amounts to one
876  // restart
877  do
878  {
879  // reset this vector to the right size
880  h.reinit(n_tmp_vectors - 1);
881 
882  if (left_precondition)
883  {
884  A.vmult(p, x);
885  p.sadd(-1., 1., b);
886  preconditioner.vmult(v, p);
887  }
888  else
889  {
890  A.vmult(v, x);
891  v.sadd(-1., 1., b);
892  };
893 
894  double rho = v.l2_norm();
895 
896  // check the residual here as well since it may be that we got the exact
897  // (or an almost exact) solution vector at the outset. if we wouldn't
898  // check here, the next scaling operation would produce garbage
899  if (use_default_residual)
900  {
901  last_res = rho;
902  iteration_state =
903  this->iteration_status(accumulated_iterations, rho, x);
904 
905  if (iteration_state != SolverControl::iterate)
906  break;
907  }
908  else
909  {
910  deallog << "default_res=" << rho << std::endl;
911 
912  if (left_precondition)
913  {
914  A.vmult(*r, x);
915  r->sadd(-1., 1., b);
916  }
917  else
918  preconditioner.vmult(*r, v);
919 
920  double res = r->l2_norm();
921  last_res = res;
922  iteration_state =
923  this->iteration_status(accumulated_iterations, res, x);
924 
925  if (iteration_state != SolverControl::iterate)
926  break;
927  }
928 
929  gamma(0) = rho;
930 
931  v *= 1. / rho;
932 
933  // inner iteration doing at most as many steps as there are temporary
934  // vectors. the number of steps actually been done is propagated outside
935  // through the @p dim variable
936  for (unsigned int inner_iteration = 0;
937  ((inner_iteration < n_tmp_vectors - 2) &&
938  (iteration_state == SolverControl::iterate));
939  ++inner_iteration)
940  {
941  ++accumulated_iterations;
942  // yet another alias
943  VectorType &vv = tmp_vectors(inner_iteration + 1, x);
944 
945  if (left_precondition)
946  {
947  A.vmult(p, tmp_vectors[inner_iteration]);
948  preconditioner.vmult(vv, p);
949  }
950  else
951  {
952  preconditioner.vmult(p, tmp_vectors[inner_iteration]);
953  A.vmult(vv, p);
954  }
955 
956  dim = inner_iteration + 1;
957 
958  const double s = modified_gram_schmidt(tmp_vectors,
959  dim,
960  accumulated_iterations,
961  vv,
962  h,
963  re_orthogonalize,
964  re_orthogonalize_signal);
965  h(inner_iteration + 1) = s;
966 
967  // s=0 is a lucky breakdown, the solver will reach convergence,
968  // but we must not divide by zero here.
969  if (s != 0)
970  vv *= 1. / s;
971 
972  // for eigenvalues, get the resulting coefficients from the
973  // orthogonalization process
974  if (do_eigenvalues)
975  for (unsigned int i = 0; i < dim + 1; ++i)
976  H_orig(i, inner_iteration) = h(i);
977 
978  // Transformation into tridiagonal structure
979  givens_rotation(h, gamma, ci, si, inner_iteration);
980 
981  // append vector on matrix
982  for (unsigned int i = 0; i < dim; ++i)
983  H(i, inner_iteration) = h(i);
984 
985  // default residual
986  rho = std::fabs(gamma(dim));
987 
988  if (use_default_residual)
989  {
990  last_res = rho;
991  iteration_state =
992  this->iteration_status(accumulated_iterations, rho, x);
993  }
994  else
995  {
996  deallog << "default_res=" << rho << std::endl;
997 
998  ::Vector<double> h_(dim);
999  *x_ = x;
1000  *gamma_ = gamma;
1001  H1.reinit(dim + 1, dim);
1002 
1003  for (unsigned int i = 0; i < dim + 1; ++i)
1004  for (unsigned int j = 0; j < dim; ++j)
1005  H1(i, j) = H(i, j);
1006 
1007  H1.backward(h_, *gamma_);
1008 
1009  if (left_precondition)
1010  for (unsigned int i = 0; i < dim; ++i)
1011  x_->add(h_(i), tmp_vectors[i]);
1012  else
1013  {
1014  p = 0.;
1015  for (unsigned int i = 0; i < dim; ++i)
1016  p.add(h_(i), tmp_vectors[i]);
1017  preconditioner.vmult(*r, p);
1018  x_->add(1., *r);
1019  };
1020  A.vmult(*r, *x_);
1021  r->sadd(-1., 1., b);
1022  // Now *r contains the unpreconditioned residual!!
1023  if (left_precondition)
1024  {
1025  const double res = r->l2_norm();
1026  last_res = res;
1027 
1028  iteration_state =
1029  this->iteration_status(accumulated_iterations, res, x);
1030  }
1031  else
1032  {
1033  preconditioner.vmult(*x_, *r);
1034  const double preconditioned_res = x_->l2_norm();
1035  last_res = preconditioned_res;
1036 
1037  iteration_state =
1038  this->iteration_status(accumulated_iterations,
1039  preconditioned_res,
1040  x);
1041  }
1042  }
1043  };
1044  // end of inner iteration. now calculate the solution from the temporary
1045  // vectors
1046  h.reinit(dim);
1047  H1.reinit(dim + 1, dim);
1048 
1049  for (unsigned int i = 0; i < dim + 1; ++i)
1050  for (unsigned int j = 0; j < dim; ++j)
1051  H1(i, j) = H(i, j);
1052 
1053  compute_eigs_and_cond(H_orig,
1054  dim,
1055  all_eigenvalues_signal,
1056  all_hessenberg_signal,
1057  condition_number_signal);
1058 
1059  H1.backward(h, gamma);
1060 
1061  if (left_precondition)
1062  for (unsigned int i = 0; i < dim; ++i)
1063  x.add(h(i), tmp_vectors[i]);
1064  else
1065  {
1066  p = 0.;
1067  for (unsigned int i = 0; i < dim; ++i)
1068  p.add(h(i), tmp_vectors[i]);
1069  preconditioner.vmult(v, p);
1070  x.add(1., v);
1071  };
1072  // end of outer iteration. restart if no convergence and the number of
1073  // iterations is not exceeded
1074  }
1075  while (iteration_state == SolverControl::iterate);
1076 
1077  compute_eigs_and_cond(H_orig,
1078  dim,
1079  eigenvalues_signal,
1080  hessenberg_signal,
1081  condition_number_signal);
1082 
1083  if (!krylov_space_signal.empty())
1084  krylov_space_signal(tmp_vectors);
1085 
1086  // in case of failure: throw exception
1087  AssertThrow(iteration_state == SolverControl::success,
1088  SolverControl::NoConvergence(accumulated_iterations, last_res));
1089 }
1090 
1091 
1092 
1093 template <class VectorType>
1094 boost::signals2::connection
1096  const std::function<void(double)> &slot,
1097  const bool every_iteration)
1098 {
1099  if (every_iteration)
1100  {
1101  return all_condition_numbers_signal.connect(slot);
1102  }
1103  else
1104  {
1105  return condition_number_signal.connect(slot);
1106  }
1107 }
1108 
1109 
1110 
1111 template <class VectorType>
1112 boost::signals2::connection
1114  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1115  const bool every_iteration)
1116 {
1117  if (every_iteration)
1118  {
1119  return all_eigenvalues_signal.connect(slot);
1120  }
1121  else
1122  {
1123  return eigenvalues_signal.connect(slot);
1124  }
1125 }
1126 
1127 
1128 
1129 template <class VectorType>
1130 boost::signals2::connection
1132  const std::function<void(const FullMatrix<double> &)> &slot,
1133  const bool every_iteration)
1134 {
1135  if (every_iteration)
1136  {
1137  return all_hessenberg_signal.connect(slot);
1138  }
1139  else
1140  {
1141  return hessenberg_signal.connect(slot);
1142  }
1143 }
1144 
1145 
1146 
1147 template <class VectorType>
1148 boost::signals2::connection
1150  const std::function<void(
1152 {
1153  return krylov_space_signal.connect(slot);
1154 }
1155 
1156 
1157 
1158 template <class VectorType>
1159 boost::signals2::connection
1161  const std::function<void(int)> &slot)
1162 {
1163  return re_orthogonalize_signal.connect(slot);
1164 }
1165 
1166 
1167 
1168 template <class VectorType>
1169 double
1171 {
1172  // dummy implementation. this function is not needed for the present
1173  // implementation of gmres
1174  Assert(false, ExcInternalError());
1175  return 0;
1176 }
1177 
1178 
1179 //----------------------------------------------------------------------//
1180 
1181 template <class VectorType>
1184  const AdditionalData & data)
1185  : SolverBase<VectorType>(cn, mem)
1186  , additional_data(data)
1187 {}
1188 
1189 
1190 
1191 template <class VectorType>
1193  const AdditionalData &data)
1194  : SolverBase<VectorType>(cn)
1195  , additional_data(data)
1196 {}
1197 
1198 
1199 
1200 template <class VectorType>
1201 template <typename MatrixType, typename PreconditionerType>
1202 void
1203 SolverFGMRES<VectorType>::solve(const MatrixType & A,
1204  VectorType & x,
1205  const VectorType & b,
1206  const PreconditionerType &preconditioner)
1207 {
1208  LogStream::Prefix prefix("FGMRES");
1209 
1210  SolverControl::State iteration_state = SolverControl::iterate;
1211 
1212  const unsigned int basis_size = additional_data.max_basis_size;
1213 
1214  // Generate an object where basis vectors are stored.
1216  basis_size, this->memory);
1218  basis_size, this->memory);
1219 
1220  // number of the present iteration; this number is not reset to zero upon a
1221  // restart
1222  unsigned int accumulated_iterations = 0;
1223 
1224  // matrix used for the orthogonalization process later
1225  H.reinit(basis_size + 1, basis_size);
1226 
1227  // Vectors for projected system
1228  Vector<double> projected_rhs;
1229  Vector<double> y;
1230 
1231  // Iteration starts here
1232  double res = -std::numeric_limits<double>::max();
1233 
1234  typename VectorMemory<VectorType>::Pointer aux(this->memory);
1235  aux->reinit(x);
1236  do
1237  {
1238  A.vmult(*aux, x);
1239  aux->sadd(-1., 1., b);
1240 
1241  double beta = aux->l2_norm();
1242  res = beta;
1243  iteration_state = this->iteration_status(accumulated_iterations, res, x);
1244  if (iteration_state == SolverControl::success)
1245  break;
1246 
1247  H.reinit(basis_size + 1, basis_size);
1248  double a = beta;
1249 
1250  for (unsigned int j = 0; j < basis_size; ++j)
1251  {
1252  if (a != 0) // treat lucky breakdown
1253  v(j, x).equ(1. / a, *aux);
1254  else
1255  v(j, x) = 0.;
1256 
1257 
1258  preconditioner.vmult(z(j, x), v[j]);
1259  A.vmult(*aux, z[j]);
1260 
1261  // Gram-Schmidt
1262  H(0, j) = *aux * v[0];
1263  for (unsigned int i = 1; i <= j; ++i)
1264  H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1265  H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1266 
1267  // Compute projected solution
1268 
1269  if (j > 0)
1270  {
1271  H1.reinit(j + 1, j);
1272  projected_rhs.reinit(j + 1);
1273  y.reinit(j);
1274  projected_rhs(0) = beta;
1275  H1.fill(H);
1276 
1277  // check convergence. note that the vector 'x' we pass to the
1278  // criterion is not the final solution we compute if we
1279  // decide to jump out of the iteration (we update 'x' again
1280  // right after the current loop)
1281  Householder<double> house(H1);
1282  res = house.least_squares(y, projected_rhs);
1283  iteration_state =
1284  this->iteration_status(++accumulated_iterations, res, x);
1285  if (iteration_state != SolverControl::iterate)
1286  break;
1287  }
1288  }
1289 
1290  // Update solution vector
1291  for (unsigned int j = 0; j < y.size(); ++j)
1292  x.add(y(j), z[j]);
1293  }
1294  while (iteration_state == SolverControl::iterate);
1295 
1296  // in case of failure: throw exception
1297  if (iteration_state != SolverControl::success)
1298  AssertThrow(false,
1299  SolverControl::NoConvergence(accumulated_iterations, res));
1300 }
1301 
1302 #endif // DOXYGEN
1303 
1304 DEAL_II_NAMESPACE_CLOSE
1305 
1306 #endif
SolverFGMRES(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:339
Number add_and_dot(const Number a, const Vector< Number > &V, const Vector< Number > &W)
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
Continue iteration.
boost::signals2::signal< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> krylov_space_signal
Definition: solver_gmres.h:367
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
FullMatrix< double > H1
Definition: solver_gmres.h:531
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_gmres.h:326
AdditionalData additional_data
Definition: solver_gmres.h:320
static ::ExceptionBase & ExcNotInitialized()
std::vector< typename VectorMemory< VectorType >::Pointer > data
Definition: solver_gmres.h:105
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
AdditionalData additional_data
Definition: solver_gmres.h:521
FullMatrix< double > H1
Definition: solver_gmres.h:438
VectorType & operator()(const unsigned int i, const VectorType &temp)
AdditionalData(const unsigned int max_basis_size=30)
Definition: solver_gmres.h:471
boost::signals2::signal< void(int)> re_orthogonalize_signal
Definition: solver_gmres.h:373
FullMatrix< double > H
Definition: solver_gmres.h:433
static ::ExceptionBase & ExcTooFewTmpVectors(int arg1)
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
static ::ExceptionBase & ExcMessage(std::string arg1)
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)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
Stop iteration, goal reached.
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1407
boost::signals2::connection connect_krylov_space_slot(const std::function< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> &slot)
static double modified_gram_schmidt(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &orthogonal_vectors, const unsigned int dim, const unsigned int accumulated_iterations, VectorType &vv, Vector< double > &h, bool &re_orthogonalize, const boost::signals2::signal< void(int)> &re_orthogonalize_signal=boost::signals2::signal< void(int)>())
boost::signals2::signal< void(const FullMatrix< double > &)> hessenberg_signal
Definition: solver_gmres.h:352
AdditionalData(const unsigned int max_basis_size, const bool use_default_residual)
Definition: solver_gmres.h:480
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
boost::signals2::signal< void(const std::vector< std::complex< double >> &)> all_eigenvalues_signal
Definition: solver_gmres.h:346
virtual double criterion()
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
size_type size() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
boost::signals2::signal< void(const FullMatrix< double > &)> all_hessenberg_signal
Definition: solver_gmres.h:359
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
VectorType & operator[](const unsigned int i) const
boost::signals2::connection connect_hessenberg_slot(const std::function< void(const FullMatrix< double > &)> &slot, const bool every_iteration=true)
FullMatrix< double > H
Definition: solver_gmres.h:526
SolverGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
void givens_rotation(Vector< double > &h, Vector< double > &b, Vector< double > &ci, Vector< double > &si, int col) const
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< std::complex< double >> &)> &slot, const bool every_iteration=false)
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(const FullMatrix< double > &)> &hessenberg_signal, const boost::signals2::signal< void(double)> &cond_signal)
boost::signals2::signal< void(double)> all_condition_numbers_signal
Definition: solver_gmres.h:332
static ::ExceptionBase & ExcInternalError()