Reference documentation for deal.II version GIT 6bdf9a4b6c 2022-08-12 19:20: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 
45 namespace internal
46 {
50  namespace SolverGMRESImplementation
51  {
60  template <typename VectorType>
61  class TmpVectors
62  {
63  public:
68  TmpVectors(const unsigned int max_size, VectorMemory<VectorType> &vmem);
69 
73  ~TmpVectors() = default;
74 
79  VectorType &
80  operator[](const unsigned int i) const;
81 
88  VectorType &
89  operator()(const unsigned int i, const VectorType &temp);
90 
95  unsigned int
96  size() const;
97 
98 
99  private:
104 
108  std::vector<typename VectorMemory<VectorType>::Pointer> data;
109  };
110  } // namespace SolverGMRESImplementation
111 } // namespace internal
112 
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 
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
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 
461 template <class VectorType = Vector<double>>
462 class SolverFGMRES : public SolverBase<VectorType>
463 {
464 public:
469  {
473  explicit AdditionalData(const unsigned int max_basis_size = 30)
475  {}
476 
480  unsigned int max_basis_size;
481  };
482 
488  const AdditionalData & data = AdditionalData());
489 
495  const AdditionalData &data = AdditionalData());
496 
500  template <typename MatrixType, typename PreconditionerType>
501  void
502  solve(const MatrixType & A,
503  VectorType & x,
504  const VectorType & b,
505  const PreconditionerType &preconditioner);
506 
507 private:
512 
517 
522 };
523 
525 /* --------------------- Inline and template functions ------------------- */
526 
527 
528 #ifndef DOXYGEN
529 namespace internal
530 {
531  namespace SolverGMRESImplementation
532  {
533  template <class VectorType>
534  inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
536  : mem(vmem)
537  , data(max_size)
538  {}
539 
540 
541 
542  template <class VectorType>
543  inline VectorType &
544  TmpVectors<VectorType>::operator[](const unsigned int i) const
545  {
546  AssertIndexRange(i, data.size());
547 
548  Assert(data[i] != nullptr, ExcNotInitialized());
549  return *data[i];
550  }
551 
552 
553 
554  template <class VectorType>
555  inline VectorType &
556  TmpVectors<VectorType>::operator()(const unsigned int i,
557  const VectorType & temp)
558  {
559  AssertIndexRange(i, data.size());
560  if (data[i] == nullptr)
561  {
562  data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
563  data[i]->reinit(temp);
564  }
565  return *data[i];
566  }
567 
568 
569 
570  template <class VectorType>
571  unsigned int
573  {
574  return (data.size() > 0 ? data.size() - 1 : 0);
575  }
576 
577 
578 
579  // A comparator for better printing eigenvalues
580  inline bool
581  complex_less_pred(const std::complex<double> &x,
582  const std::complex<double> &y)
583  {
584  return x.real() < y.real() ||
585  (x.real() == y.real() && x.imag() < y.imag());
586  }
587  } // namespace SolverGMRESImplementation
588 } // namespace internal
589 
590 
591 
592 template <class VectorType>
594  const unsigned int max_n_tmp_vectors,
595  const bool right_preconditioning,
596  const bool use_default_residual,
597  const bool force_re_orthogonalization)
598  : max_n_tmp_vectors(max_n_tmp_vectors)
599  , right_preconditioning(right_preconditioning)
600  , use_default_residual(use_default_residual)
601  , force_re_orthogonalization(force_re_orthogonalization)
602 {
603  Assert(3 <= max_n_tmp_vectors,
604  ExcMessage("SolverGMRES needs at least three "
605  "temporary vectors."));
606 }
607 
608 
609 
610 template <class VectorType>
613  const AdditionalData & data)
614  : SolverBase<VectorType>(cn, mem)
615  , additional_data(data)
616 {}
617 
618 
619 
620 template <class VectorType>
622  const AdditionalData &data)
623  : SolverBase<VectorType>(cn)
624  , additional_data(data)
625 {}
626 
627 
628 
629 template <class VectorType>
630 inline void
632  Vector<double> &b,
633  Vector<double> &ci,
634  Vector<double> &si,
635  int col) const
636 {
637  for (int i = 0; i < col; ++i)
638  {
639  const double s = si(i);
640  const double c = ci(i);
641  const double dummy = h(i);
642  h(i) = c * dummy + s * h(i + 1);
643  h(i + 1) = -s * dummy + c * h(i + 1);
644  };
645 
646  const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
647  si(col) = h(col + 1) * r;
648  ci(col) = h(col) * r;
649  h(col) = ci(col) * h(col) + si(col) * h(col + 1);
650  b(col + 1) = -si(col) * b(col);
651  b(col) *= ci(col);
652 }
653 
654 
655 
656 template <class VectorType>
657 inline double
660  & orthogonal_vectors,
661  const unsigned int dim,
662  const unsigned int accumulated_iterations,
663  VectorType & vv,
664  Vector<double> & h,
665  bool & reorthogonalize,
666  const boost::signals2::signal<void(int)> &reorthogonalize_signal)
667 {
668  Assert(dim > 0, ExcInternalError());
669  const unsigned int inner_iteration = dim - 1;
670 
671  // need initial norm for detection of re-orthogonalization, see below
672  double norm_vv_start = 0;
673  const bool consider_reorthogonalize =
674  (reorthogonalize == false) && (inner_iteration % 5 == 4);
675  if (consider_reorthogonalize)
676  norm_vv_start = vv.l2_norm();
677 
678  // Orthogonalization
679  h(0) = vv * orthogonal_vectors[0];
680  for (unsigned int i = 1; i < dim; ++i)
681  h(i) = vv.add_and_dot(-h(i - 1),
682  orthogonal_vectors[i - 1],
683  orthogonal_vectors[i]);
684  double norm_vv =
685  std::sqrt(vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv));
686 
687  // Re-orthogonalization if loss of orthogonality detected. For the test, use
688  // a strategy discussed in C. T. Kelley, Iterative Methods for Linear and
689  // Nonlinear Equations, SIAM, Philadelphia, 1995: Compare the norm of vv
690  // after orthogonalization with its norm when starting the
691  // orthogonalization. If vv became very small (here: less than the square
692  // root of the machine precision times 10), it is almost in the span of the
693  // previous vectors, which indicates loss of precision.
694  if (consider_reorthogonalize)
695  {
696  if (norm_vv >
697  10. * norm_vv_start *
698  std::sqrt(
700  return norm_vv;
701 
702  else
703  {
704  reorthogonalize = true;
705  if (!reorthogonalize_signal.empty())
706  reorthogonalize_signal(accumulated_iterations);
707  }
708  }
709 
710  if (reorthogonalize == true)
711  {
712  double htmp = vv * orthogonal_vectors[0];
713  h(0) += htmp;
714  for (unsigned int i = 1; i < dim; ++i)
715  {
716  htmp = vv.add_and_dot(-htmp,
717  orthogonal_vectors[i - 1],
718  orthogonal_vectors[i]);
719  h(i) += htmp;
720  }
721  norm_vv =
722  std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
723  }
724 
725  return norm_vv;
726 }
727 
728 
729 
730 template <class VectorType>
731 inline void
733  const FullMatrix<double> &H_orig,
734  const unsigned int dim,
735  const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
736  &eigenvalues_signal,
737  const boost::signals2::signal<void(const FullMatrix<double> &)>
738  & hessenberg_signal,
739  const boost::signals2::signal<void(double)> &cond_signal)
740 {
741  // Avoid copying the Hessenberg matrix if it isn't needed.
742  if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
743  !cond_signal.empty()) &&
744  dim > 0)
745  {
746  LAPACKFullMatrix<double> mat(dim, dim);
747  for (unsigned int i = 0; i < dim; ++i)
748  for (unsigned int j = 0; j < dim; ++j)
749  mat(i, j) = H_orig(i, j);
750  hessenberg_signal(H_orig);
751  // Avoid computing eigenvalues if they are not needed.
752  if (!eigenvalues_signal.empty())
753  {
754  // Copy mat so that we can compute svd below. Necessary since
755  // compute_eigenvalues will leave mat in state
756  // LAPACKSupport::unusable.
757  LAPACKFullMatrix<double> mat_eig(mat);
758  mat_eig.compute_eigenvalues();
759  std::vector<std::complex<double>> eigenvalues(dim);
760  for (unsigned int i = 0; i < mat_eig.n(); ++i)
761  eigenvalues[i] = mat_eig.eigenvalue(i);
762  // Sort eigenvalues for nicer output.
763  std::sort(eigenvalues.begin(),
764  eigenvalues.end(),
765  internal::SolverGMRESImplementation::complex_less_pred);
766  eigenvalues_signal(eigenvalues);
767  }
768  // Calculate condition number, avoid calculating the svd if a slot
769  // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
770  if (!cond_signal.empty() && (mat.n() > 1))
771  {
772  mat.compute_svd();
773  double condition_number =
774  mat.singular_value(0) / mat.singular_value(mat.n() - 1);
775  cond_signal(condition_number);
776  }
777  }
778 }
779 
780 
781 
782 template <class VectorType>
783 template <typename MatrixType, typename PreconditionerType>
784 void
785 SolverGMRES<VectorType>::solve(const MatrixType & A,
786  VectorType & x,
787  const VectorType & b,
788  const PreconditionerType &preconditioner)
789 {
790  // TODO:[?] Check, why there are two different start residuals.
791  // TODO:[GK] Make sure the parameter in the constructor means maximum basis
792  // size
793 
794  LogStream::Prefix prefix("GMRES");
795 
796  // extra call to std::max to placate static analyzers: coverity rightfully
797  // complains that data.max_n_tmp_vectors - 2 may overflow
798  const unsigned int n_tmp_vectors =
799  std::max(additional_data.max_n_tmp_vectors, 3u);
800 
801  // Generate an object where basis vectors are stored.
803  n_tmp_vectors, this->memory);
804 
805  // number of the present iteration; this
806  // number is not reset to zero upon a
807  // restart
808  unsigned int accumulated_iterations = 0;
809 
810  const bool do_eigenvalues =
811  !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
812  !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty() ||
813  !hessenberg_signal.empty() || !all_hessenberg_signal.empty();
814  // for eigenvalue computation, need to collect the Hessenberg matrix (before
815  // applying Givens rotations)
816  FullMatrix<double> H_orig;
817  if (do_eigenvalues)
818  H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
819 
820  // matrix used for the orthogonalization process later
821  H.reinit(n_tmp_vectors, n_tmp_vectors - 1);
822 
823  // some additional vectors, also used in the orthogonalization
824  ::Vector<double> gamma(n_tmp_vectors), ci(n_tmp_vectors - 1),
825  si(n_tmp_vectors - 1), h(n_tmp_vectors - 1);
826 
827 
828  unsigned int dim = 0;
829 
831  double last_res = std::numeric_limits<double>::lowest();
832 
833  // switch to determine whether we want a left or a right preconditioner. at
834  // present, left is default, but both ways are implemented
835  const bool left_precondition = !additional_data.right_preconditioning;
836 
837  // Per default the left preconditioned GMRes uses the preconditioned
838  // residual and the right preconditioned GMRes uses the unpreconditioned
839  // residual as stopping criterion.
840  const bool use_default_residual = additional_data.use_default_residual;
841 
842  // define two aliases
843  VectorType &v = tmp_vectors(0, x);
844  VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
845 
846  // Following vectors are needed when we are not using the default residuals
847  // as stopping criterion
850  std::unique_ptr<::Vector<double>> gamma_;
851  if (!use_default_residual)
852  {
853  r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
854  x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
855  r->reinit(x);
856  x_->reinit(x);
857 
858  gamma_ = std::make_unique<::Vector<double>>(gamma.size());
859  }
860 
861  bool re_orthogonalize = additional_data.force_re_orthogonalization;
862 
864  // outer iteration: loop until we either reach convergence or the maximum
865  // number of iterations is exceeded. each cycle of this loop amounts to one
866  // restart
867  do
868  {
869  // reset this vector to the right size
870  h.reinit(n_tmp_vectors - 1);
871 
872  if (left_precondition)
873  {
874  A.vmult(p, x);
875  p.sadd(-1., 1., b);
876  preconditioner.vmult(v, p);
877  }
878  else
879  {
880  A.vmult(v, x);
881  v.sadd(-1., 1., b);
882  };
883 
884  double rho = v.l2_norm();
885 
886  // check the residual here as well since it may be that we got the exact
887  // (or an almost exact) solution vector at the outset. if we wouldn't
888  // check here, the next scaling operation would produce garbage
889  if (use_default_residual)
890  {
891  last_res = rho;
892  iteration_state =
893  this->iteration_status(accumulated_iterations, rho, x);
894 
895  if (iteration_state != SolverControl::iterate)
896  break;
897  }
898  else
899  {
900  deallog << "default_res=" << rho << std::endl;
901 
902  if (left_precondition)
903  {
904  A.vmult(*r, x);
905  r->sadd(-1., 1., b);
906  }
907  else
908  preconditioner.vmult(*r, v);
909 
910  double res = r->l2_norm();
911  last_res = res;
912  iteration_state =
913  this->iteration_status(accumulated_iterations, res, x);
914 
915  if (iteration_state != SolverControl::iterate)
916  break;
917  }
918 
919  gamma(0) = rho;
920 
921  v *= 1. / rho;
922 
923  // inner iteration doing at most as many steps as there are temporary
924  // vectors. the number of steps actually been done is propagated outside
925  // through the @p dim variable
926  for (unsigned int inner_iteration = 0;
927  ((inner_iteration < n_tmp_vectors - 2) &&
928  (iteration_state == SolverControl::iterate));
929  ++inner_iteration)
930  {
931  ++accumulated_iterations;
932  // yet another alias
933  VectorType &vv = tmp_vectors(inner_iteration + 1, x);
934 
935  if (left_precondition)
936  {
937  A.vmult(p, tmp_vectors[inner_iteration]);
938  preconditioner.vmult(vv, p);
939  }
940  else
941  {
942  preconditioner.vmult(p, tmp_vectors[inner_iteration]);
943  A.vmult(vv, p);
944  }
945 
946  dim = inner_iteration + 1;
947 
948  const double s = modified_gram_schmidt(tmp_vectors,
949  dim,
950  accumulated_iterations,
951  vv,
952  h,
953  re_orthogonalize,
954  re_orthogonalize_signal);
955  h(inner_iteration + 1) = s;
956 
957  // s=0 is a lucky breakdown, the solver will reach convergence,
958  // but we must not divide by zero here.
959  if (s != 0)
960  vv *= 1. / s;
961 
962  // for eigenvalues, get the resulting coefficients from the
963  // orthogonalization process
964  if (do_eigenvalues)
965  for (unsigned int i = 0; i < dim + 1; ++i)
966  H_orig(i, inner_iteration) = h(i);
967 
968  // Transformation into tridiagonal structure
969  givens_rotation(h, gamma, ci, si, inner_iteration);
970 
971  // append vector on matrix
972  for (unsigned int i = 0; i < dim; ++i)
973  H(i, inner_iteration) = h(i);
974 
975  // default residual
976  rho = std::fabs(gamma(dim));
977 
978  if (use_default_residual)
979  {
980  last_res = rho;
981  iteration_state =
982  this->iteration_status(accumulated_iterations, rho, x);
983  }
984  else
985  {
986  deallog << "default_res=" << rho << std::endl;
987 
988  ::Vector<double> h_(dim);
989  *x_ = x;
990  *gamma_ = gamma;
991  H1.reinit(dim + 1, dim);
992 
993  for (unsigned int i = 0; i < dim + 1; ++i)
994  for (unsigned int j = 0; j < dim; ++j)
995  H1(i, j) = H(i, j);
996 
997  H1.backward(h_, *gamma_);
998 
999  if (left_precondition)
1000  for (unsigned int i = 0; i < dim; ++i)
1001  x_->add(h_(i), tmp_vectors[i]);
1002  else
1003  {
1004  p = 0.;
1005  for (unsigned int i = 0; i < dim; ++i)
1006  p.add(h_(i), tmp_vectors[i]);
1007  preconditioner.vmult(*r, p);
1008  x_->add(1., *r);
1009  };
1010  A.vmult(*r, *x_);
1011  r->sadd(-1., 1., b);
1012  // Now *r contains the unpreconditioned residual!!
1013  if (left_precondition)
1014  {
1015  const double res = r->l2_norm();
1016  last_res = res;
1017 
1018  iteration_state =
1019  this->iteration_status(accumulated_iterations, res, x);
1020  }
1021  else
1022  {
1023  preconditioner.vmult(*x_, *r);
1024  const double preconditioned_res = x_->l2_norm();
1025  last_res = preconditioned_res;
1026 
1027  iteration_state =
1028  this->iteration_status(accumulated_iterations,
1029  preconditioned_res,
1030  x);
1031  }
1032  }
1033  };
1034  // end of inner iteration. now calculate the solution from the temporary
1035  // vectors
1036  h.reinit(dim);
1037  H1.reinit(dim + 1, dim);
1038 
1039  for (unsigned int i = 0; i < dim + 1; ++i)
1040  for (unsigned int j = 0; j < dim; ++j)
1041  H1(i, j) = H(i, j);
1042 
1043  compute_eigs_and_cond(H_orig,
1044  dim,
1045  all_eigenvalues_signal,
1046  all_hessenberg_signal,
1047  condition_number_signal);
1048 
1049  H1.backward(h, gamma);
1050 
1051  if (left_precondition)
1052  for (unsigned int i = 0; i < dim; ++i)
1053  x.add(h(i), tmp_vectors[i]);
1054  else
1055  {
1056  p = 0.;
1057  for (unsigned int i = 0; i < dim; ++i)
1058  p.add(h(i), tmp_vectors[i]);
1059  preconditioner.vmult(v, p);
1060  x.add(1., v);
1061  };
1062  // end of outer iteration. restart if no convergence and the number of
1063  // iterations is not exceeded
1064  }
1065  while (iteration_state == SolverControl::iterate);
1066 
1067  compute_eigs_and_cond(H_orig,
1068  dim,
1069  eigenvalues_signal,
1070  hessenberg_signal,
1071  condition_number_signal);
1072 
1073  if (!krylov_space_signal.empty())
1074  krylov_space_signal(tmp_vectors);
1075 
1076  // in case of failure: throw exception
1077  AssertThrow(iteration_state == SolverControl::success,
1078  SolverControl::NoConvergence(accumulated_iterations, last_res));
1079 }
1080 
1081 
1082 
1083 template <class VectorType>
1084 boost::signals2::connection
1086  const std::function<void(double)> &slot,
1087  const bool every_iteration)
1088 {
1089  if (every_iteration)
1090  {
1091  return all_condition_numbers_signal.connect(slot);
1092  }
1093  else
1094  {
1095  return condition_number_signal.connect(slot);
1096  }
1097 }
1098 
1099 
1100 
1101 template <class VectorType>
1102 boost::signals2::connection
1104  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1105  const bool every_iteration)
1106 {
1107  if (every_iteration)
1108  {
1109  return all_eigenvalues_signal.connect(slot);
1110  }
1111  else
1112  {
1113  return eigenvalues_signal.connect(slot);
1114  }
1115 }
1116 
1117 
1118 
1119 template <class VectorType>
1120 boost::signals2::connection
1122  const std::function<void(const FullMatrix<double> &)> &slot,
1123  const bool every_iteration)
1124 {
1125  if (every_iteration)
1126  {
1127  return all_hessenberg_signal.connect(slot);
1128  }
1129  else
1130  {
1131  return hessenberg_signal.connect(slot);
1132  }
1133 }
1134 
1135 
1136 
1137 template <class VectorType>
1138 boost::signals2::connection
1140  const std::function<void(
1142 {
1143  return krylov_space_signal.connect(slot);
1144 }
1145 
1146 
1147 
1148 template <class VectorType>
1149 boost::signals2::connection
1151  const std::function<void(int)> &slot)
1152 {
1153  return re_orthogonalize_signal.connect(slot);
1154 }
1155 
1156 
1157 
1158 template <class VectorType>
1159 double
1161 {
1162  // dummy implementation. this function is not needed for the present
1163  // implementation of gmres
1164  Assert(false, ExcInternalError());
1165  return 0;
1166 }
1167 
1168 
1169 //----------------------------------------------------------------------//
1170 
1171 template <class VectorType>
1174  const AdditionalData & data)
1175  : SolverBase<VectorType>(cn, mem)
1176  , additional_data(data)
1177 {}
1178 
1179 
1180 
1181 template <class VectorType>
1183  const AdditionalData &data)
1184  : SolverBase<VectorType>(cn)
1185  , additional_data(data)
1186 {}
1187 
1188 
1189 
1190 template <class VectorType>
1191 template <typename MatrixType, typename PreconditionerType>
1192 void
1193 SolverFGMRES<VectorType>::solve(const MatrixType & A,
1194  VectorType & x,
1195  const VectorType & b,
1196  const PreconditionerType &preconditioner)
1197 {
1198  LogStream::Prefix prefix("FGMRES");
1199 
1200  SolverControl::State iteration_state = SolverControl::iterate;
1201 
1202  const unsigned int basis_size = additional_data.max_basis_size;
1203 
1204  // Generate an object where basis vectors are stored.
1206  basis_size, this->memory);
1208  basis_size, this->memory);
1209 
1210  // number of the present iteration; this number is not reset to zero upon a
1211  // restart
1212  unsigned int accumulated_iterations = 0;
1213 
1214  // matrix used for the orthogonalization process later
1215  H.reinit(basis_size + 1, basis_size);
1216 
1217  // Vectors for projected system
1218  Vector<double> projected_rhs;
1219  Vector<double> y;
1220 
1221  // Iteration starts here
1222  double res = std::numeric_limits<double>::lowest();
1223 
1224  typename VectorMemory<VectorType>::Pointer aux(this->memory);
1225  aux->reinit(x);
1226  do
1227  {
1228  A.vmult(*aux, x);
1229  aux->sadd(-1., 1., b);
1230 
1231  double beta = aux->l2_norm();
1232  res = beta;
1233  iteration_state = this->iteration_status(accumulated_iterations, res, x);
1234  if (iteration_state == SolverControl::success)
1235  break;
1236 
1237  H.reinit(basis_size + 1, basis_size);
1238  double a = beta;
1239 
1240  for (unsigned int j = 0; j < basis_size; ++j)
1241  {
1242  if (a != 0) // treat lucky breakdown
1243  v(j, x).equ(1. / a, *aux);
1244  else
1245  v(j, x) = 0.;
1246 
1247 
1248  preconditioner.vmult(z(j, x), v[j]);
1249  A.vmult(*aux, z[j]);
1250 
1251  // Gram-Schmidt
1252  H(0, j) = *aux * v[0];
1253  for (unsigned int i = 1; i <= j; ++i)
1254  H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1255  H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1256 
1257  // Compute projected solution
1258 
1259  if (j > 0)
1260  {
1261  H1.reinit(j + 1, j);
1262  projected_rhs.reinit(j + 1);
1263  y.reinit(j);
1264  projected_rhs(0) = beta;
1265  H1.fill(H);
1266 
1267  // check convergence. note that the vector 'x' we pass to the
1268  // criterion is not the final solution we compute if we
1269  // decide to jump out of the iteration (we update 'x' again
1270  // right after the current loop)
1271  Householder<double> house(H1);
1272  res = house.least_squares(y, projected_rhs);
1273  iteration_state =
1274  this->iteration_status(++accumulated_iterations, res, x);
1275  if (iteration_state != SolverControl::iterate)
1276  break;
1277  }
1278  }
1279 
1280  // Update solution vector
1281  for (unsigned int j = 0; j < y.size(); ++j)
1282  x.add(y(j), z[j]);
1283  }
1284  while (iteration_state == SolverControl::iterate);
1285 
1286  // in case of failure: throw exception
1287  if (iteration_state != SolverControl::success)
1288  AssertThrow(false,
1289  SolverControl::NoConvergence(accumulated_iterations, res));
1290 }
1291 
1292 #endif // DOXYGEN
1293 
1295 
1296 #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:516
SolverFGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
Definition: solver_gmres.h:511
FullMatrix< double > H1
Definition: solver_gmres.h:521
boost::signals2::signal< void(const FullMatrix< double > &)> hessenberg_signal
Definition: solver_gmres.h:352
AdditionalData additional_data
Definition: solver_gmres.h:320
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_gmres.h:326
SolverGMRES(const SolverGMRES< VectorType > &)=delete
boost::signals2::signal< void(const std::vector< std::complex< double >> &)> eigenvalues_signal
Definition: solver_gmres.h:339
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:367
boost::signals2::signal< void(const FullMatrix< double > &)> all_hessenberg_signal
Definition: solver_gmres.h:359
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:438
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:346
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:433
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:332
boost::signals2::signal< void(int)> re_orthogonalize_signal
Definition: solver_gmres.h:373
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
Definition: vector.h:110
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)
std::vector< typename VectorMemory< VectorType >::Pointer > data
Definition: solver_gmres.h:108
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
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:473
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)