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