Reference documentation for deal.II version GIT 05ffa62ef0 2022-11-26 15:30: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>
26 
30 #include <deal.II/lac/solver.h>
32 #include <deal.II/lac/vector.h>
33 
34 #include <algorithm>
35 #include <cmath>
36 #include <limits>
37 #include <memory>
38 #include <vector>
39 
41 
42 // forward declarations
43 #ifndef DOXYGEN
44 namespace LinearAlgebra
45 {
46  namespace distributed
47  {
48  template <typename, typename>
49  class Vector;
50  } // namespace distributed
51 } // namespace LinearAlgebra
52 #endif
53 
59 namespace internal
60 {
64  namespace SolverGMRESImplementation
65  {
74  template <typename VectorType>
75  class TmpVectors
76  {
77  public:
82  TmpVectors(const unsigned int max_size, VectorMemory<VectorType> &vmem);
83 
87  ~TmpVectors() = default;
88 
93  VectorType &
94  operator[](const unsigned int i) const;
95 
102  VectorType &
103  operator()(const unsigned int i, const VectorType &temp);
104 
109  unsigned int
110  size() const;
111 
112 
113  private:
118 
122  std::vector<typename VectorMemory<VectorType>::Pointer> data;
123  };
124  } // namespace SolverGMRESImplementation
125 } // namespace internal
126 
191 template <class VectorType = Vector<double>>
192 class SolverGMRES : public SolverBase<VectorType>
193 {
194 public:
199  {
201  {
210  };
211 
219  explicit AdditionalData(
220  const unsigned int max_n_tmp_vectors = 30,
221  const bool right_preconditioning = false,
222  const bool use_default_residual = true,
223  const bool force_re_orthogonalization = false,
224  const bool batched_mode = false,
227 
234  unsigned int max_n_tmp_vectors;
235 
244 
249 
257 
265 
270  };
271 
277  const AdditionalData & data = AdditionalData());
278 
284 
289 
293  template <typename MatrixType, typename PreconditionerType>
294  void
295  solve(const MatrixType & A,
296  VectorType & x,
297  const VectorType & b,
298  const PreconditionerType &preconditioner);
299 
306  boost::signals2::connection
307  connect_condition_number_slot(const std::function<void(double)> &slot,
308  const bool every_iteration = false);
309 
316  boost::signals2::connection
318  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
319  const bool every_iteration = false);
320 
328  boost::signals2::connection
330  const std::function<void(const FullMatrix<double> &)> &slot,
331  const bool every_iteration = true);
332 
339  boost::signals2::connection
341  const std::function<
343  &slot);
344 
345 
350  boost::signals2::connection
351  connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
352 
353 
355  int,
356  << "The number of temporary vectors you gave (" << arg1
357  << ") is too small. It should be at least 10 for "
358  << "any results, and much more for reasonable ones.");
359 
360 protected:
365 
370  boost::signals2::signal<void(double)> condition_number_signal;
371 
376  boost::signals2::signal<void(double)> all_condition_numbers_signal;
377 
382  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
384 
389  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
391 
396  boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
397 
402  boost::signals2::signal<void(const FullMatrix<double> &)>
404 
409  boost::signals2::signal<void(
412 
417  boost::signals2::signal<void(int)> re_orthogonalize_signal;
418 
425 
429  virtual double
431 
436  void
438  Vector<double> &b,
441  int col) const;
442 
449  static void
451  const FullMatrix<double> &H_orig,
452  const unsigned int dim,
453  const boost::signals2::signal<
454  void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
455  const boost::signals2::signal<void(const FullMatrix<double> &)>
457  const boost::signals2::signal<void(double)> &cond_signal);
458 
463 
468 
473 
478 
483 };
484 
485 
486 
507 template <class VectorType = Vector<double>>
508 class SolverFGMRES : public SolverBase<VectorType>
509 {
510 public:
515  {
519  explicit AdditionalData(const unsigned int max_basis_size = 30)
521  {}
522 
526  unsigned int max_basis_size;
527  };
528 
534  const AdditionalData & data = AdditionalData());
535 
541  const AdditionalData &data = AdditionalData());
542 
546  template <typename MatrixType, typename PreconditionerType>
547  void
548  solve(const MatrixType & A,
549  VectorType & x,
550  const VectorType & b,
551  const PreconditionerType &preconditioner);
552 
553 private:
558 
563 
568 };
569 
571 /* --------------------- Inline and template functions ------------------- */
572 
573 
574 #ifndef DOXYGEN
575 namespace internal
576 {
577  namespace SolverGMRESImplementation
578  {
579  template <class VectorType>
580  inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
582  : mem(vmem)
583  , data(max_size)
584  {}
585 
586 
587 
588  template <class VectorType>
589  inline VectorType &
590  TmpVectors<VectorType>::operator[](const unsigned int i) const
591  {
592  AssertIndexRange(i, data.size());
593 
594  Assert(data[i] != nullptr, ExcNotInitialized());
595  return *data[i];
596  }
597 
598 
599 
600  template <class VectorType>
601  inline VectorType &
602  TmpVectors<VectorType>::operator()(const unsigned int i,
603  const VectorType & temp)
604  {
605  AssertIndexRange(i, data.size());
606  if (data[i] == nullptr)
607  {
608  data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
609  data[i]->reinit(temp, true);
610  }
611  return *data[i];
612  }
613 
614 
615 
616  template <class VectorType>
617  unsigned int
619  {
620  return (data.size() > 0 ? data.size() - 1 : 0);
621  }
622 
623 
624 
625  // A comparator for better printing eigenvalues
626  inline bool
627  complex_less_pred(const std::complex<double> &x,
628  const std::complex<double> &y)
629  {
630  return x.real() < y.real() ||
631  (x.real() == y.real() && x.imag() < y.imag());
632  }
633 
634  // A function to solve the (upper) triangular system after Givens
635  // rotations on a matrix that has possibly unused rows and columns
636  inline void
637  solve_triangular(const unsigned int dim,
638  const FullMatrix<double> &H,
639  const Vector<double> & rhs,
640  Vector<double> & solution)
641  {
642  for (int i = dim - 1; i >= 0; --i)
643  {
644  double s = rhs(i);
645  for (unsigned int j = i + 1; j < dim; ++j)
646  s -= solution(j) * H(i, j);
647  solution(i) = s / H(i, i);
648  AssertIsFinite(solution(i));
649  }
650  }
651  } // namespace SolverGMRESImplementation
652 } // namespace internal
653 
654 
655 
656 template <class VectorType>
658  const unsigned int max_n_tmp_vectors,
659  const bool right_preconditioning,
660  const bool use_default_residual,
661  const bool force_re_orthogonalization,
662  const bool batched_mode,
663  const OrthogonalizationStrategy orthogonalization_strategy)
664  : max_n_tmp_vectors(max_n_tmp_vectors)
665  , right_preconditioning(right_preconditioning)
666  , use_default_residual(use_default_residual)
667  , force_re_orthogonalization(force_re_orthogonalization)
668  , batched_mode(batched_mode)
669  , orthogonalization_strategy(orthogonalization_strategy)
670 {
671  Assert(3 <= max_n_tmp_vectors,
672  ExcMessage("SolverGMRES needs at least three "
673  "temporary vectors."));
674 }
675 
676 
677 
678 template <class VectorType>
681  const AdditionalData & data)
682  : SolverBase<VectorType>(cn, mem)
683  , additional_data(data)
684  , solver_control(cn)
685 {}
686 
687 
688 
689 template <class VectorType>
691  const AdditionalData &data)
692  : SolverBase<VectorType>(cn)
693  , additional_data(data)
694  , solver_control(cn)
695 {}
696 
697 
698 
699 template <class VectorType>
700 inline void
702  Vector<double> &b,
703  Vector<double> &ci,
704  Vector<double> &si,
705  int col) const
706 {
707  for (int i = 0; i < col; ++i)
708  {
709  const double s = si(i);
710  const double c = ci(i);
711  const double dummy = h(i);
712  h(i) = c * dummy + s * h(i + 1);
713  h(i + 1) = -s * dummy + c * h(i + 1);
714  };
715 
716  const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
717  si(col) = h(col + 1) * r;
718  ci(col) = h(col) * r;
719  h(col) = ci(col) * h(col) + si(col) * h(col + 1);
720  b(col + 1) = -si(col) * b(col);
721  b(col) *= ci(col);
722 }
723 
724 
725 
726 namespace internal
727 {
728  namespace SolverGMRESImplementation
729  {
730  template <class VectorType>
731  void
732  Tvmult_add(const unsigned int dim,
733  const VectorType & vv,
735  & orthogonal_vectors,
736  Vector<double> &h)
737  {
738  for (unsigned int i = 0; i < dim; ++i)
739  h[i] += vv * orthogonal_vectors[i];
740  }
741 
742 
743 
744  template <class Number>
745  void
746  Tvmult_add(
747  const unsigned int dim,
751  & orthogonal_vectors,
752  Vector<double> &h)
753  {
754  unsigned int j = 0;
755 
756  if (dim <= 128)
757  {
758  // optimized path
759  static constexpr unsigned int n_lanes =
761 
762  VectorizedArray<double> hs[128];
763  for (unsigned int d = 0; d < dim; ++d)
764  hs[d] = 0.0;
765 
766  unsigned int c = 0;
767 
768  for (; c < vv.locally_owned_size() / n_lanes / 4;
769  ++c, j += n_lanes * 4)
770  for (unsigned int i = 0; i < dim; ++i)
771  {
772  VectorizedArray<double> vvec[4];
773  for (unsigned int k = 0; k < 4; ++k)
774  vvec[k].load(vv.begin() + j + k * n_lanes);
775 
776  for (unsigned int k = 0; k < 4; ++k)
777  {
779  temp.load(orthogonal_vectors[i].begin() + j + k * n_lanes);
780  hs[i] += temp * vvec[k];
781  }
782  }
783 
784  c *= 4;
785  for (; c < vv.locally_owned_size() / n_lanes; ++c, j += n_lanes)
786  for (unsigned int i = 0; i < dim; ++i)
787  {
788  VectorizedArray<double> vvec, temp;
789  vvec.load(vv.begin() + j);
790  temp.load(orthogonal_vectors[i].begin() + j);
791  hs[i] += temp * vvec;
792  }
793 
794  for (unsigned int i = 0; i < dim; ++i)
795  for (unsigned int v = 0; v < n_lanes; ++v)
796  h(i) += hs[i][v];
797  }
798 
799  // remainder loop of optimized path or non-optimized path (if
800  // dim>128)
801  for (; j < vv.locally_owned_size(); ++j)
802  for (unsigned int i = 0; i < dim; ++i)
803  h(i) += orthogonal_vectors[i].local_element(j) * vv.local_element(j);
804 
805  Utilities::MPI::sum(h, MPI_COMM_WORLD, h);
806  }
807 
808 
809 
810  template <class VectorType>
811  double
812  substract_and_norm(
813  const unsigned int dim,
815  & orthogonal_vectors,
816  const Vector<double> &h,
817  VectorType & vv)
818  {
819  Assert(dim > 0, ExcInternalError());
820 
821  for (unsigned int i = 0; i < dim - 1; ++i)
822  vv.add(-h(i), orthogonal_vectors[i]);
823 
824  return vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv);
825  }
826 
827 
828 
829  template <class Number>
830  double
831  substract_and_norm(
832  const unsigned int dim,
835  & orthogonal_vectors,
836  const Vector<double> &h,
838  {
839  static constexpr unsigned int n_lanes = VectorizedArray<double>::size();
840 
841  double norm_vv_temp = 0.0;
842  VectorizedArray<double> norm_vv_temp_vectorized = 0.0;
843 
844  unsigned int j = 0;
845  unsigned int c = 0;
846  for (; c < vv.locally_owned_size() / n_lanes / 4; ++c, j += n_lanes * 4)
847  {
848  VectorizedArray<double> temp[4];
849 
850  for (unsigned int k = 0; k < 4; ++k)
851  temp[k].load(vv.begin() + j + k * n_lanes);
852 
853  for (unsigned int i = 0; i < dim; ++i)
854  {
855  const double factor = h(i);
856  for (unsigned int k = 0; k < 4; ++k)
857  {
859  vec.load(orthogonal_vectors[i].begin() + j + k * n_lanes);
860  temp[k] -= factor * vec;
861  }
862  }
863 
864  for (unsigned int k = 0; k < 4; ++k)
865  temp[k].store(vv.begin() + j + k * n_lanes);
866 
867  norm_vv_temp_vectorized += (temp[0] * temp[0] + temp[1] * temp[1]) +
868  (temp[2] * temp[2] + temp[3] * temp[3]);
869  }
870 
871  c *= 4;
872  for (; c < vv.locally_owned_size() / n_lanes; ++c, j += n_lanes)
873  {
875  temp.load(vv.begin() + j);
876 
877  for (unsigned int i = 0; i < dim; ++i)
878  {
880  vec.load(orthogonal_vectors[i].begin() + j);
881  temp -= h(i) * vec;
882  }
883 
884  temp.store(vv.begin() + j);
885 
886  norm_vv_temp_vectorized += temp * temp;
887  }
888 
889  for (unsigned int v = 0; v < n_lanes; ++v)
890  norm_vv_temp += norm_vv_temp_vectorized[v];
891 
892  for (; j < vv.locally_owned_size(); ++j)
893  {
894  double temp = vv(j);
895  for (unsigned int i = 0; i < dim; ++i)
896  temp -= h(i) * orthogonal_vectors[i](j);
897  vv(j) = temp;
898 
899  norm_vv_temp += temp * temp;
900  }
901 
902  return std::sqrt(Utilities::MPI::sum(norm_vv_temp, MPI_COMM_WORLD));
903  }
904 
905 
906  template <class VectorType>
907  double
908  sadd_and_norm(VectorType & v,
909  const double factor_a,
910  const VectorType &b,
911  const double factor_b)
912  {
913  v.sadd(factor_a, factor_b, b);
914  return v.l2_norm();
915  }
916 
917 
918  template <class Number>
919  double
920  sadd_and_norm(
922  const double factor_a,
924  const double factor_b)
925  {
926  double norm = 0;
927 
928  for (unsigned int j = 0; j < v.locally_owned_size(); ++j)
929  {
930  const double temp =
931  v.local_element(j) * factor_a + b.local_element(j) * factor_b;
932 
933  v.local_element(j) = temp;
934 
935  norm += temp * temp;
936  }
937 
938  return std::sqrt(Utilities::MPI::sum(norm, MPI_COMM_WORLD));
939  }
940 
941 
942 
943  template <class VectorType>
944  void
945  add(VectorType & p,
946  const unsigned int dim,
947  const Vector<double> &h,
949  & tmp_vectors,
950  const bool zero_out)
951  {
952  if (zero_out)
953  p.equ(h(0), tmp_vectors[0]);
954  else
955  p.add(h(0), tmp_vectors[0]);
956 
957  for (unsigned int i = 1; i < dim; ++i)
958  p.add(h(i), tmp_vectors[i]);
959  }
960 
961 
962 
963  template <class Number>
964  void
966  const unsigned int dim,
967  const Vector<double> & h,
970  & tmp_vectors,
971  const bool zero_out)
972  {
973  for (unsigned int j = 0; j < p.locally_owned_size(); ++j)
974  {
975  double temp = zero_out ? 0 : p.local_element(j);
976  for (unsigned int i = 0; i < dim; ++i)
977  temp += tmp_vectors[i].local_element(j) * h(i);
978  p.local_element(j) = temp;
979  }
980  }
981 
982 
983 
995  template <class VectorType>
996  inline double
997  iterated_gram_schmidt(
999  OrthogonalizationStrategy orthogonalization_strategy,
1001  & orthogonal_vectors,
1002  const unsigned int dim,
1003  const unsigned int accumulated_iterations,
1004  VectorType & vv,
1005  Vector<double> & h,
1006  bool & reorthogonalize,
1007  const boost::signals2::signal<void(int)> &reorthogonalize_signal =
1008  boost::signals2::signal<void(int)>())
1009  {
1010  Assert(dim > 0, ExcInternalError());
1011  const unsigned int inner_iteration = dim - 1;
1012 
1013  // need initial norm for detection of re-orthogonalization, see below
1014  double norm_vv_start = 0;
1015  const bool consider_reorthogonalize =
1016  (reorthogonalize == false) && (inner_iteration % 5 == 4);
1017  if (consider_reorthogonalize)
1018  norm_vv_start = vv.l2_norm();
1019 
1020  for (unsigned int i = 0; i < dim; ++i)
1021  h[i] = 0;
1022 
1023  for (unsigned int c = 0; c < 2;
1024  ++c) // 0: orthogonalize, 1: reorthogonalize
1025  {
1026  // Orthogonalization
1027  double norm_vv = 0.0;
1028 
1029  if (orthogonalization_strategy ==
1031  OrthogonalizationStrategy::modified_gram_schmidt)
1032  {
1033  double htmp = vv * orthogonal_vectors[0];
1034  h(0) += htmp;
1035  for (unsigned int i = 1; i < dim; ++i)
1036  {
1037  htmp = vv.add_and_dot(-htmp,
1038  orthogonal_vectors[i - 1],
1039  orthogonal_vectors[i]);
1040  h(i) += htmp;
1041  }
1042 
1043  norm_vv = std::sqrt(
1044  vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
1045  }
1046  else if (orthogonalization_strategy ==
1048  OrthogonalizationStrategy::classical_gram_schmidt)
1049  {
1050  Tvmult_add(dim, vv, orthogonal_vectors, h);
1051  norm_vv = substract_and_norm(dim, orthogonal_vectors, h, vv);
1052  }
1053  else
1054  {
1055  AssertThrow(false, ExcNotImplemented());
1056  }
1057 
1058  if (c == 1)
1059  return norm_vv; // reorthogonalization already performed -> finished
1060 
1061  // Re-orthogonalization if loss of orthogonality detected. For the
1062  // test, use a strategy discussed in C. T. Kelley, Iterative Methods
1063  // for Linear and Nonlinear Equations, SIAM, Philadelphia, 1995:
1064  // Compare the norm of vv after orthogonalization with its norm when
1065  // starting the orthogonalization. If vv became very small (here: less
1066  // than the square root of the machine precision times 10), it is
1067  // almost in the span of the previous vectors, which indicates loss of
1068  // precision.
1069  if (consider_reorthogonalize)
1070  {
1071  if (norm_vv >
1072  10. * norm_vv_start *
1073  std::sqrt(std::numeric_limits<
1074  typename VectorType::value_type>::epsilon()))
1075  return norm_vv;
1076 
1077  else
1078  {
1079  reorthogonalize = true;
1080  if (!reorthogonalize_signal.empty())
1081  reorthogonalize_signal(accumulated_iterations);
1082  }
1083  }
1084 
1085  if (reorthogonalize == false)
1086  return norm_vv; // no reorthogonalization needed -> finished
1087  }
1088 
1089  AssertThrow(false, ExcInternalError());
1090 
1091  return 0.0;
1092  }
1093  } // namespace SolverGMRESImplementation
1094 } // namespace internal
1095 
1096 
1097 
1098 template <class VectorType>
1099 inline void
1101  const FullMatrix<double> &H_orig,
1102  const unsigned int dim,
1103  const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
1104  &eigenvalues_signal,
1105  const boost::signals2::signal<void(const FullMatrix<double> &)>
1106  & hessenberg_signal,
1107  const boost::signals2::signal<void(double)> &cond_signal)
1108 {
1109  // Avoid copying the Hessenberg matrix if it isn't needed.
1110  if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1111  !cond_signal.empty()) &&
1112  dim > 0)
1113  {
1114  LAPACKFullMatrix<double> mat(dim, dim);
1115  for (unsigned int i = 0; i < dim; ++i)
1116  for (unsigned int j = 0; j < dim; ++j)
1117  mat(i, j) = H_orig(i, j);
1118  hessenberg_signal(H_orig);
1119  // Avoid computing eigenvalues if they are not needed.
1120  if (!eigenvalues_signal.empty())
1121  {
1122  // Copy mat so that we can compute svd below. Necessary since
1123  // compute_eigenvalues will leave mat in state
1124  // LAPACKSupport::unusable.
1125  LAPACKFullMatrix<double> mat_eig(mat);
1126  mat_eig.compute_eigenvalues();
1127  std::vector<std::complex<double>> eigenvalues(dim);
1128  for (unsigned int i = 0; i < mat_eig.n(); ++i)
1129  eigenvalues[i] = mat_eig.eigenvalue(i);
1130  // Sort eigenvalues for nicer output.
1131  std::sort(eigenvalues.begin(),
1132  eigenvalues.end(),
1133  internal::SolverGMRESImplementation::complex_less_pred);
1134  eigenvalues_signal(eigenvalues);
1135  }
1136  // Calculate condition number, avoid calculating the svd if a slot
1137  // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
1138  if (!cond_signal.empty() && (mat.n() > 1))
1139  {
1140  mat.compute_svd();
1141  double condition_number =
1142  mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1143  cond_signal(condition_number);
1144  }
1145  }
1146 }
1147 
1148 
1149 
1150 template <class VectorType>
1151 template <typename MatrixType, typename PreconditionerType>
1152 void
1153 SolverGMRES<VectorType>::solve(const MatrixType & A,
1154  VectorType & x,
1155  const VectorType & b,
1156  const PreconditionerType &preconditioner)
1157 {
1158  // TODO:[?] Check, why there are two different start residuals.
1159  // TODO:[GK] Make sure the parameter in the constructor means maximum basis
1160  // size
1161 
1162  std::unique_ptr<LogStream::Prefix> prefix;
1163  if (!additional_data.batched_mode)
1164  prefix = std::make_unique<LogStream::Prefix>("GMRES");
1165 
1166  // extra call to std::max to placate static analyzers: coverity rightfully
1167  // complains that data.max_n_tmp_vectors - 2 may overflow
1168  const unsigned int n_tmp_vectors =
1169  std::max(additional_data.max_n_tmp_vectors, 3u);
1170 
1171  // Generate an object where basis vectors are stored.
1173  n_tmp_vectors, this->memory);
1174 
1175  // number of the present iteration; this
1176  // number is not reset to zero upon a
1177  // restart
1178  unsigned int accumulated_iterations = 0;
1179 
1180  const bool do_eigenvalues =
1181  !additional_data.batched_mode &&
1182  (!condition_number_signal.empty() ||
1183  !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1184  !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1185  !all_hessenberg_signal.empty());
1186  // for eigenvalue computation, need to collect the Hessenberg matrix (before
1187  // applying Givens rotations)
1188  FullMatrix<double> H_orig;
1189  if (do_eigenvalues)
1190  H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
1191 
1192  // matrix used for the orthogonalization process later
1193  H.reinit(n_tmp_vectors, n_tmp_vectors - 1, /* omit_initialization */ true);
1194 
1195  // some additional vectors, also used in the orthogonalization
1196  gamma.reinit(n_tmp_vectors);
1197  ci.reinit(n_tmp_vectors - 1);
1198  si.reinit(n_tmp_vectors - 1);
1199  h.reinit(n_tmp_vectors - 1);
1200 
1201  unsigned int dim = 0;
1202 
1203  SolverControl::State iteration_state = SolverControl::iterate;
1204  double last_res = std::numeric_limits<double>::lowest();
1205 
1206  // switch to determine whether we want a left or a right preconditioner. at
1207  // present, left is default, but both ways are implemented
1208  const bool left_precondition = !additional_data.right_preconditioning;
1209 
1210  // Per default the left preconditioned GMRes uses the preconditioned
1211  // residual and the right preconditioned GMRes uses the unpreconditioned
1212  // residual as stopping criterion.
1213  const bool use_default_residual = additional_data.use_default_residual;
1214 
1215  // define two aliases
1216  VectorType &v = tmp_vectors(0, x);
1217  VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
1218 
1219  // Following vectors are needed when we are not using the default residuals
1220  // as stopping criterion
1223  std::unique_ptr<::Vector<double>> gamma_;
1224  if (!use_default_residual)
1225  {
1226  r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1227  x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1228  r->reinit(x);
1229  x_->reinit(x);
1230 
1231  gamma_ = std::make_unique<::Vector<double>>(gamma.size());
1232  }
1233 
1234  bool re_orthogonalize = additional_data.force_re_orthogonalization;
1235 
1237  // outer iteration: loop until we either reach convergence or the maximum
1238  // number of iterations is exceeded. each cycle of this loop amounts to one
1239  // restart
1240  do
1241  {
1242  // reset this vector to the right size
1243  h.reinit(n_tmp_vectors - 1);
1244 
1245  double rho = 0.0;
1246 
1247  if (left_precondition)
1248  {
1249  A.vmult(p, x);
1250  p.sadd(-1., 1., b);
1251  preconditioner.vmult(v, p);
1252  rho = v.l2_norm();
1253  }
1254  else
1255  {
1256  A.vmult(v, x);
1257  rho = ::internal::SolverGMRESImplementation::sadd_and_norm(v,
1258  -1,
1259  b,
1260  1.0);
1261  }
1262 
1263  // check the residual here as well since it may be that we got the exact
1264  // (or an almost exact) solution vector at the outset. if we wouldn't
1265  // check here, the next scaling operation would produce garbage
1266  if (use_default_residual)
1267  {
1268  last_res = rho;
1269  if (additional_data.batched_mode)
1270  iteration_state = solver_control.check(accumulated_iterations, rho);
1271  else
1272  iteration_state =
1273  this->iteration_status(accumulated_iterations, rho, x);
1274 
1275  if (iteration_state != SolverControl::iterate)
1276  break;
1277  }
1278  else
1279  {
1280  deallog << "default_res=" << rho << std::endl;
1281 
1282  if (left_precondition)
1283  {
1284  A.vmult(*r, x);
1285  r->sadd(-1., 1., b);
1286  }
1287  else
1288  preconditioner.vmult(*r, v);
1289 
1290  double res = r->l2_norm();
1291  last_res = res;
1292  if (additional_data.batched_mode)
1293  iteration_state = solver_control.check(accumulated_iterations, rho);
1294  else
1295  iteration_state =
1296  this->iteration_status(accumulated_iterations, res, x);
1297 
1298  if (iteration_state != SolverControl::iterate)
1299  break;
1300  }
1301 
1302  gamma(0) = rho;
1303 
1304  v *= 1. / rho;
1305 
1306  // inner iteration doing at most as many steps as there are temporary
1307  // vectors. the number of steps actually been done is propagated outside
1308  // through the @p dim variable
1309  for (unsigned int inner_iteration = 0;
1310  ((inner_iteration < n_tmp_vectors - 2) &&
1311  (iteration_state == SolverControl::iterate));
1312  ++inner_iteration)
1313  {
1314  ++accumulated_iterations;
1315  // yet another alias
1316  VectorType &vv = tmp_vectors(inner_iteration + 1, x);
1317 
1318  if (left_precondition)
1319  {
1320  A.vmult(p, tmp_vectors[inner_iteration]);
1321  preconditioner.vmult(vv, p);
1322  }
1323  else
1324  {
1325  preconditioner.vmult(p, tmp_vectors[inner_iteration]);
1326  A.vmult(vv, p);
1327  }
1328 
1329  dim = inner_iteration + 1;
1330 
1331  const double s =
1332  internal::SolverGMRESImplementation::iterated_gram_schmidt(
1333  additional_data.orthogonalization_strategy,
1334  tmp_vectors,
1335  dim,
1336  accumulated_iterations,
1337  vv,
1338  h,
1339  re_orthogonalize,
1340  re_orthogonalize_signal);
1341  h(inner_iteration + 1) = s;
1342 
1343  // s=0 is a lucky breakdown, the solver will reach convergence,
1344  // but we must not divide by zero here.
1345  if (s != 0)
1346  vv *= 1. / s;
1347 
1348  // for eigenvalues, get the resulting coefficients from the
1349  // orthogonalization process
1350  if (do_eigenvalues)
1351  for (unsigned int i = 0; i < dim + 1; ++i)
1352  H_orig(i, inner_iteration) = h(i);
1353 
1354  // Transformation into tridiagonal structure
1355  givens_rotation(h, gamma, ci, si, inner_iteration);
1356 
1357  // append vector on matrix
1358  for (unsigned int i = 0; i < dim; ++i)
1359  H(i, inner_iteration) = h(i);
1360 
1361  // default residual
1362  rho = std::fabs(gamma(dim));
1363 
1364  if (use_default_residual)
1365  {
1366  last_res = rho;
1367  if (additional_data.batched_mode)
1368  iteration_state =
1369  solver_control.check(accumulated_iterations, rho);
1370  else
1371  iteration_state =
1372  this->iteration_status(accumulated_iterations, rho, x);
1373  }
1374  else
1375  {
1376  if (!additional_data.batched_mode)
1377  deallog << "default_res=" << rho << std::endl;
1378 
1379  *x_ = x;
1380  *gamma_ = gamma;
1381  internal::SolverGMRESImplementation::solve_triangular(dim,
1382  H,
1383  *gamma_,
1384  h);
1385 
1386  if (left_precondition)
1387  for (unsigned int i = 0; i < dim; ++i)
1388  x_->add(h(i), tmp_vectors[i]);
1389  else
1390  {
1391  p = 0.;
1392  for (unsigned int i = 0; i < dim; ++i)
1393  p.add(h(i), tmp_vectors[i]);
1394  preconditioner.vmult(*r, p);
1395  x_->add(1., *r);
1396  };
1397  A.vmult(*r, *x_);
1398  r->sadd(-1., 1., b);
1399  // Now *r contains the unpreconditioned residual!!
1400  if (left_precondition)
1401  {
1402  const double res = r->l2_norm();
1403  last_res = res;
1404 
1405  iteration_state =
1406  this->iteration_status(accumulated_iterations, res, x);
1407  }
1408  else
1409  {
1410  preconditioner.vmult(*x_, *r);
1411  const double preconditioned_res = x_->l2_norm();
1412  last_res = preconditioned_res;
1413 
1414  if (additional_data.batched_mode)
1415  iteration_state =
1416  solver_control.check(accumulated_iterations, rho);
1417  else
1418  iteration_state =
1419  this->iteration_status(accumulated_iterations,
1420  preconditioned_res,
1421  x);
1422  }
1423  }
1424  }
1425 
1426  // end of inner iteration. now calculate the solution from the temporary
1427  // vectors
1428  internal::SolverGMRESImplementation::solve_triangular(dim, H, gamma, h);
1429 
1430  if (do_eigenvalues)
1431  compute_eigs_and_cond(H_orig,
1432  dim,
1433  all_eigenvalues_signal,
1434  all_hessenberg_signal,
1435  condition_number_signal);
1436 
1437  if (left_precondition)
1438  ::internal::SolverGMRESImplementation::add(
1439  x, dim, h, tmp_vectors, false);
1440  else
1441  {
1442  ::internal::SolverGMRESImplementation::add(
1443  p, dim, h, tmp_vectors, true);
1444  preconditioner.vmult(v, p);
1445  x.add(1., v);
1446  };
1447  // end of outer iteration. restart if no convergence and the number of
1448  // iterations is not exceeded
1449  }
1450  while (iteration_state == SolverControl::iterate);
1451 
1452  if (do_eigenvalues)
1453  compute_eigs_and_cond(H_orig,
1454  dim,
1455  eigenvalues_signal,
1456  hessenberg_signal,
1457  condition_number_signal);
1458 
1459  if (!additional_data.batched_mode && !krylov_space_signal.empty())
1460  krylov_space_signal(tmp_vectors);
1461 
1462  // in case of failure: throw exception
1463  AssertThrow(iteration_state == SolverControl::success,
1464  SolverControl::NoConvergence(accumulated_iterations, last_res));
1465 }
1466 
1467 
1468 
1469 template <class VectorType>
1470 boost::signals2::connection
1472  const std::function<void(double)> &slot,
1473  const bool every_iteration)
1474 {
1475  if (every_iteration)
1476  {
1477  return all_condition_numbers_signal.connect(slot);
1478  }
1479  else
1480  {
1481  return condition_number_signal.connect(slot);
1482  }
1483 }
1484 
1485 
1486 
1487 template <class VectorType>
1488 boost::signals2::connection
1490  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1491  const bool every_iteration)
1492 {
1493  if (every_iteration)
1494  {
1495  return all_eigenvalues_signal.connect(slot);
1496  }
1497  else
1498  {
1499  return eigenvalues_signal.connect(slot);
1500  }
1501 }
1502 
1503 
1504 
1505 template <class VectorType>
1506 boost::signals2::connection
1508  const std::function<void(const FullMatrix<double> &)> &slot,
1509  const bool every_iteration)
1510 {
1511  if (every_iteration)
1512  {
1513  return all_hessenberg_signal.connect(slot);
1514  }
1515  else
1516  {
1517  return hessenberg_signal.connect(slot);
1518  }
1519 }
1520 
1521 
1522 
1523 template <class VectorType>
1524 boost::signals2::connection
1526  const std::function<void(
1528 {
1529  return krylov_space_signal.connect(slot);
1530 }
1531 
1532 
1533 
1534 template <class VectorType>
1535 boost::signals2::connection
1537  const std::function<void(int)> &slot)
1538 {
1539  return re_orthogonalize_signal.connect(slot);
1540 }
1541 
1542 
1543 
1544 template <class VectorType>
1545 double
1547 {
1548  // dummy implementation. this function is not needed for the present
1549  // implementation of gmres
1550  Assert(false, ExcInternalError());
1551  return 0;
1552 }
1553 
1554 
1555 //----------------------------------------------------------------------//
1556 
1557 template <class VectorType>
1560  const AdditionalData & data)
1561  : SolverBase<VectorType>(cn, mem)
1562  , additional_data(data)
1563 {}
1564 
1565 
1566 
1567 template <class VectorType>
1569  const AdditionalData &data)
1570  : SolverBase<VectorType>(cn)
1571  , additional_data(data)
1572 {}
1573 
1574 
1575 
1576 template <class VectorType>
1577 template <typename MatrixType, typename PreconditionerType>
1578 void
1579 SolverFGMRES<VectorType>::solve(const MatrixType & A,
1580  VectorType & x,
1581  const VectorType & b,
1582  const PreconditionerType &preconditioner)
1583 {
1584  LogStream::Prefix prefix("FGMRES");
1585 
1586  SolverControl::State iteration_state = SolverControl::iterate;
1587 
1588  const unsigned int basis_size = additional_data.max_basis_size;
1589 
1590  // Generate an object where basis vectors are stored.
1592  basis_size, this->memory);
1594  basis_size, this->memory);
1595 
1596  // number of the present iteration; this number is not reset to zero upon a
1597  // restart
1598  unsigned int accumulated_iterations = 0;
1599 
1600  // matrix used for the orthogonalization process later
1601  H.reinit(basis_size + 1, basis_size);
1602 
1603  // Vectors for projected system
1604  Vector<double> projected_rhs;
1605  Vector<double> y;
1606 
1607  // Iteration starts here
1608  double res = std::numeric_limits<double>::lowest();
1609 
1610  typename VectorMemory<VectorType>::Pointer aux(this->memory);
1611  aux->reinit(x);
1612  do
1613  {
1614  A.vmult(*aux, x);
1615  aux->sadd(-1., 1., b);
1616 
1617  double beta = aux->l2_norm();
1618  res = beta;
1619  iteration_state = this->iteration_status(accumulated_iterations, res, x);
1620  if (iteration_state == SolverControl::success)
1621  break;
1622 
1623  H.reinit(basis_size + 1, basis_size);
1624  double a = beta;
1625 
1626  for (unsigned int j = 0; j < basis_size; ++j)
1627  {
1628  if (a != 0) // treat lucky breakdown
1629  v(j, x).equ(1. / a, *aux);
1630  else
1631  v(j, x) = 0.;
1632 
1633 
1634  preconditioner.vmult(z(j, x), v[j]);
1635  A.vmult(*aux, z[j]);
1636 
1637  // Gram-Schmidt
1638  H(0, j) = *aux * v[0];
1639  for (unsigned int i = 1; i <= j; ++i)
1640  H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1641  H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1642 
1643  // Compute projected solution
1644 
1645  if (j > 0)
1646  {
1647  H1.reinit(j + 1, j);
1648  projected_rhs.reinit(j + 1);
1649  y.reinit(j);
1650  projected_rhs(0) = beta;
1651  H1.fill(H);
1652 
1653  // check convergence. note that the vector 'x' we pass to the
1654  // criterion is not the final solution we compute if we
1655  // decide to jump out of the iteration (we update 'x' again
1656  // right after the current loop)
1657  Householder<double> house(H1);
1658  res = house.least_squares(y, projected_rhs);
1659  iteration_state =
1660  this->iteration_status(++accumulated_iterations, res, x);
1661  if (iteration_state != SolverControl::iterate)
1662  break;
1663  }
1664  }
1665 
1666  // Update solution vector
1667  for (unsigned int j = 0; j < y.size(); ++j)
1668  x.add(y(j), z[j]);
1669  }
1670  while (iteration_state == SolverControl::iterate);
1671 
1672  // in case of failure: throw exception
1673  if (iteration_state != SolverControl::success)
1674  AssertThrow(false,
1675  SolverControl::NoConvergence(accumulated_iterations, res));
1676 }
1677 
1678 #endif // DOXYGEN
1679 
1681 
1682 #endif
Number local_element(const size_type local_index) const
size_type locally_owned_size() const
virtual State check(const unsigned int step, const double check_value)
@ 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:562
SolverFGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
Definition: solver_gmres.h:557
FullMatrix< double > H1
Definition: solver_gmres.h:567
boost::signals2::signal< void(const FullMatrix< double > &)> hessenberg_signal
Definition: solver_gmres.h:396
Vector< double > si
Definition: solver_gmres.h:477
AdditionalData additional_data
Definition: solver_gmres.h:364
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_gmres.h:370
SolverGMRES(const SolverGMRES< VectorType > &)=delete
boost::signals2::signal< void(const std::vector< std::complex< double >> &)> eigenvalues_signal
Definition: solver_gmres.h:383
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:411
boost::signals2::signal< void(const FullMatrix< double > &)> all_hessenberg_signal
Definition: solver_gmres.h:403
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)
Vector< double > ci
Definition: solver_gmres.h:472
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
SolverControl & solver_control
Definition: solver_gmres.h:424
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)
Vector< double > gamma
Definition: solver_gmres.h:467
boost::signals2::signal< void(const std::vector< std::complex< double >> &)> all_eigenvalues_signal
Definition: solver_gmres.h:390
Vector< double > h
Definition: solver_gmres.h:482
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
FullMatrix< double > H
Definition: solver_gmres.h:462
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:376
boost::signals2::signal< void(int)> re_orthogonalize_signal
Definition: solver_gmres.h:417
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
Definition: vector.h:109
size_type size() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
void store(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
std::vector< typename VectorMemory< VectorType >::Pointer > data
Definition: solver_gmres.h:122
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:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1786
static ::ExceptionBase & ExcTooFewTmpVectors(int arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
LogStream deallog
Definition: logstream.cc:37
Expression fabs(const Expression &x)
static const char A
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
T sum(const T &t, const MPI_Comm &mpi_communicator)
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:519
OrthogonalizationStrategy orthogonalization_strategy
Definition: solver_gmres.h:269
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, const bool batched_mode=false, const OrthogonalizationStrategy orthogonalization_strategy=OrthogonalizationStrategy::modified_gram_schmidt)