Reference documentation for deal.II version GIT 2f5445400b 2023-02-05 22: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  {
213  };
214 
222  explicit AdditionalData(
223  const unsigned int max_n_tmp_vectors = 30,
224  const bool right_preconditioning = false,
225  const bool use_default_residual = true,
226  const bool force_re_orthogonalization = false,
227  const bool batched_mode = false,
230 
237  unsigned int max_n_tmp_vectors;
238 
247 
252 
260 
268 
273  };
274 
280  const AdditionalData & data = AdditionalData());
281 
287 
292 
296  template <typename MatrixType, typename PreconditionerType>
297  void
298  solve(const MatrixType & A,
299  VectorType & x,
300  const VectorType & b,
301  const PreconditionerType &preconditioner);
302 
309  boost::signals2::connection
310  connect_condition_number_slot(const std::function<void(double)> &slot,
311  const bool every_iteration = false);
312 
319  boost::signals2::connection
321  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
322  const bool every_iteration = false);
323 
331  boost::signals2::connection
333  const std::function<void(const FullMatrix<double> &)> &slot,
334  const bool every_iteration = true);
335 
342  boost::signals2::connection
344  const std::function<
346  &slot);
347 
348 
353  boost::signals2::connection
354  connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
355 
356 
358  int,
359  << "The number of temporary vectors you gave (" << arg1
360  << ") is too small. It should be at least 10 for "
361  << "any results, and much more for reasonable ones.");
362 
363 protected:
368 
373  boost::signals2::signal<void(double)> condition_number_signal;
374 
379  boost::signals2::signal<void(double)> all_condition_numbers_signal;
380 
385  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
387 
392  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
394 
399  boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
400 
405  boost::signals2::signal<void(const FullMatrix<double> &)>
407 
412  boost::signals2::signal<void(
415 
420  boost::signals2::signal<void(int)> re_orthogonalize_signal;
421 
428 
432  virtual double
434 
439  void
441  Vector<double> &b,
444  int col) const;
445 
452  static void
454  const FullMatrix<double> &H_orig,
455  const unsigned int dim,
456  const boost::signals2::signal<
457  void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
458  const boost::signals2::signal<void(const FullMatrix<double> &)>
460  const boost::signals2::signal<void(double)> &cond_signal);
461 
466 
471 
476 
481 
486 };
487 
488 
489 
510 template <class VectorType = Vector<double>>
511 class SolverFGMRES : public SolverBase<VectorType>
512 {
513 public:
518  {
522  explicit AdditionalData(const unsigned int max_basis_size = 30)
524  {}
525 
529  unsigned int max_basis_size;
530  };
531 
537  const AdditionalData & data = AdditionalData());
538 
544  const AdditionalData &data = AdditionalData());
545 
549  template <typename MatrixType, typename PreconditionerType>
550  void
551  solve(const MatrixType & A,
552  VectorType & x,
553  const VectorType & b,
554  const PreconditionerType &preconditioner);
555 
556 private:
561 
566 
571 };
572 
574 /* --------------------- Inline and template functions ------------------- */
575 
576 
577 #ifndef DOXYGEN
578 namespace internal
579 {
580  namespace SolverGMRESImplementation
581  {
582  template <class VectorType>
583  inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
585  : mem(vmem)
586  , data(max_size)
587  {}
588 
589 
590 
591  template <class VectorType>
592  inline VectorType &
593  TmpVectors<VectorType>::operator[](const unsigned int i) const
594  {
595  AssertIndexRange(i, data.size());
596 
597  Assert(data[i] != nullptr, ExcNotInitialized());
598  return *data[i];
599  }
600 
601 
602 
603  template <class VectorType>
604  inline VectorType &
605  TmpVectors<VectorType>::operator()(const unsigned int i,
606  const VectorType & temp)
607  {
608  AssertIndexRange(i, data.size());
609  if (data[i] == nullptr)
610  {
611  data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
612  data[i]->reinit(temp, true);
613  }
614  return *data[i];
615  }
616 
617 
618 
619  template <class VectorType>
620  unsigned int
622  {
623  return (data.size() > 0 ? data.size() - 1 : 0);
624  }
625 
626 
627 
628  // A comparator for better printing eigenvalues
629  inline bool
630  complex_less_pred(const std::complex<double> &x,
631  const std::complex<double> &y)
632  {
633  return x.real() < y.real() ||
634  (x.real() == y.real() && x.imag() < y.imag());
635  }
636 
637  // A function to solve the (upper) triangular system after Givens
638  // rotations on a matrix that has possibly unused rows and columns
639  inline void
640  solve_triangular(const unsigned int dim,
641  const FullMatrix<double> &H,
642  const Vector<double> & rhs,
643  Vector<double> & solution)
644  {
645  for (int i = dim - 1; i >= 0; --i)
646  {
647  double s = rhs(i);
648  for (unsigned int j = i + 1; j < dim; ++j)
649  s -= solution(j) * H(i, j);
650  solution(i) = s / H(i, i);
651  AssertIsFinite(solution(i));
652  }
653  }
654  } // namespace SolverGMRESImplementation
655 } // namespace internal
656 
657 
658 
659 template <class VectorType>
661  const unsigned int max_n_tmp_vectors,
662  const bool right_preconditioning,
663  const bool use_default_residual,
664  const bool force_re_orthogonalization,
665  const bool batched_mode,
666  const OrthogonalizationStrategy orthogonalization_strategy)
667  : max_n_tmp_vectors(max_n_tmp_vectors)
668  , right_preconditioning(right_preconditioning)
669  , use_default_residual(use_default_residual)
670  , force_re_orthogonalization(force_re_orthogonalization)
671  , batched_mode(batched_mode)
672  , orthogonalization_strategy(orthogonalization_strategy)
673 {
674  Assert(3 <= max_n_tmp_vectors,
675  ExcMessage("SolverGMRES needs at least three "
676  "temporary vectors."));
677 }
678 
679 
680 
681 template <class VectorType>
684  const AdditionalData & data)
685  : SolverBase<VectorType>(cn, mem)
686  , additional_data(data)
687  , solver_control(cn)
688 {}
689 
690 
691 
692 template <class VectorType>
694  const AdditionalData &data)
695  : SolverBase<VectorType>(cn)
696  , additional_data(data)
697  , solver_control(cn)
698 {}
699 
700 
701 
702 template <class VectorType>
703 inline void
705  Vector<double> &b,
706  Vector<double> &ci,
707  Vector<double> &si,
708  int col) const
709 {
710  for (int i = 0; i < col; ++i)
711  {
712  const double s = si(i);
713  const double c = ci(i);
714  const double dummy = h(i);
715  h(i) = c * dummy + s * h(i + 1);
716  h(i + 1) = -s * dummy + c * h(i + 1);
717  };
718 
719  const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
720  si(col) = h(col + 1) * r;
721  ci(col) = h(col) * r;
722  h(col) = ci(col) * h(col) + si(col) * h(col + 1);
723  b(col + 1) = -si(col) * b(col);
724  b(col) *= ci(col);
725 }
726 
727 
728 
729 namespace internal
730 {
731  namespace SolverGMRESImplementation
732  {
733  template <class VectorType>
734  void
735  Tvmult_add(const unsigned int dim,
736  const VectorType & vv,
738  & orthogonal_vectors,
739  Vector<double> &h)
740  {
741  for (unsigned int i = 0; i < dim; ++i)
742  h[i] += vv * orthogonal_vectors[i];
743  }
744 
745 
746 
747  template <class Number>
748  void
749  Tvmult_add(
750  const unsigned int dim,
754  & orthogonal_vectors,
755  Vector<double> &h)
756  {
757  unsigned int j = 0;
758 
759  if (dim <= 128)
760  {
761  // optimized path
762  static constexpr unsigned int n_lanes =
764 
765  VectorizedArray<double> hs[128];
766  for (unsigned int d = 0; d < dim; ++d)
767  hs[d] = 0.0;
768 
769  unsigned int c = 0;
770 
771  for (; c < vv.locally_owned_size() / n_lanes / 4;
772  ++c, j += n_lanes * 4)
773  for (unsigned int i = 0; i < dim; ++i)
774  {
775  VectorizedArray<double> vvec[4];
776  for (unsigned int k = 0; k < 4; ++k)
777  vvec[k].load(vv.begin() + j + k * n_lanes);
778 
779  for (unsigned int k = 0; k < 4; ++k)
780  {
782  temp.load(orthogonal_vectors[i].begin() + j + k * n_lanes);
783  hs[i] += temp * vvec[k];
784  }
785  }
786 
787  c *= 4;
788  for (; c < vv.locally_owned_size() / n_lanes; ++c, j += n_lanes)
789  for (unsigned int i = 0; i < dim; ++i)
790  {
791  VectorizedArray<double> vvec, temp;
792  vvec.load(vv.begin() + j);
793  temp.load(orthogonal_vectors[i].begin() + j);
794  hs[i] += temp * vvec;
795  }
796 
797  for (unsigned int i = 0; i < dim; ++i)
798  for (unsigned int v = 0; v < n_lanes; ++v)
799  h(i) += hs[i][v];
800  }
801 
802  // remainder loop of optimized path or non-optimized path (if
803  // dim>128)
804  for (; j < vv.locally_owned_size(); ++j)
805  for (unsigned int i = 0; i < dim; ++i)
806  h(i) += orthogonal_vectors[i].local_element(j) * vv.local_element(j);
807 
808  Utilities::MPI::sum(h, MPI_COMM_WORLD, h);
809  }
810 
811 
812 
813  template <class VectorType>
814  double
815  substract_and_norm(
816  const unsigned int dim,
818  & orthogonal_vectors,
819  const Vector<double> &h,
820  VectorType & vv)
821  {
822  Assert(dim > 0, ExcInternalError());
823 
824  for (unsigned int i = 0; i < dim; ++i)
825  vv.add(-h(i), orthogonal_vectors[i]);
826 
827  return std::sqrt(vv.add_and_dot(-h(dim), orthogonal_vectors[dim], vv));
828  }
829 
830 
831 
832  template <class Number>
833  double
834  substract_and_norm(
835  const unsigned int dim,
838  & orthogonal_vectors,
839  const Vector<double> &h,
841  {
842  static constexpr unsigned int n_lanes = VectorizedArray<double>::size();
843 
844  double norm_vv_temp = 0.0;
845  VectorizedArray<double> norm_vv_temp_vectorized = 0.0;
846 
847  unsigned int j = 0;
848  unsigned int c = 0;
849  for (; c < vv.locally_owned_size() / n_lanes / 4; ++c, j += n_lanes * 4)
850  {
851  VectorizedArray<double> temp[4];
852 
853  for (unsigned int k = 0; k < 4; ++k)
854  temp[k].load(vv.begin() + j + k * n_lanes);
855 
856  for (unsigned int i = 0; i < dim; ++i)
857  {
858  const double factor = h(i);
859  for (unsigned int k = 0; k < 4; ++k)
860  {
862  vec.load(orthogonal_vectors[i].begin() + j + k * n_lanes);
863  temp[k] -= factor * vec;
864  }
865  }
866 
867  for (unsigned int k = 0; k < 4; ++k)
868  temp[k].store(vv.begin() + j + k * n_lanes);
869 
870  norm_vv_temp_vectorized += (temp[0] * temp[0] + temp[1] * temp[1]) +
871  (temp[2] * temp[2] + temp[3] * temp[3]);
872  }
873 
874  c *= 4;
875  for (; c < vv.locally_owned_size() / n_lanes; ++c, j += n_lanes)
876  {
878  temp.load(vv.begin() + j);
879 
880  for (unsigned int i = 0; i < dim; ++i)
881  {
883  vec.load(orthogonal_vectors[i].begin() + j);
884  temp -= h(i) * vec;
885  }
886 
887  temp.store(vv.begin() + j);
888 
889  norm_vv_temp_vectorized += temp * temp;
890  }
891 
892  for (unsigned int v = 0; v < n_lanes; ++v)
893  norm_vv_temp += norm_vv_temp_vectorized[v];
894 
895  for (; j < vv.locally_owned_size(); ++j)
896  {
897  double temp = vv.local_element(j);
898  for (unsigned int i = 0; i < dim; ++i)
899  temp -= h(i) * orthogonal_vectors[i].local_element(j);
900  vv.local_element(j) = temp;
901 
902  norm_vv_temp += temp * temp;
903  }
904 
905  return std::sqrt(Utilities::MPI::sum(norm_vv_temp, MPI_COMM_WORLD));
906  }
907 
908 
909  template <class VectorType>
910  double
911  sadd_and_norm(VectorType & v,
912  const double factor_a,
913  const VectorType &b,
914  const double factor_b)
915  {
916  v.sadd(factor_a, factor_b, b);
917  return v.l2_norm();
918  }
919 
920 
921  template <class Number>
922  double
923  sadd_and_norm(
925  const double factor_a,
927  const double factor_b)
928  {
929  double norm = 0;
930 
931  for (unsigned int j = 0; j < v.locally_owned_size(); ++j)
932  {
933  const double temp =
934  v.local_element(j) * factor_a + b.local_element(j) * factor_b;
935 
936  v.local_element(j) = temp;
937 
938  norm += temp * temp;
939  }
940 
941  return std::sqrt(Utilities::MPI::sum(norm, MPI_COMM_WORLD));
942  }
943 
944 
945 
946  template <class VectorType>
947  void
948  add(VectorType & p,
949  const unsigned int dim,
950  const Vector<double> &h,
952  & tmp_vectors,
953  const bool zero_out)
954  {
955  if (zero_out)
956  p.equ(h(0), tmp_vectors[0]);
957  else
958  p.add(h(0), tmp_vectors[0]);
959 
960  for (unsigned int i = 1; i < dim; ++i)
961  p.add(h(i), tmp_vectors[i]);
962  }
963 
964 
965 
966  template <class Number>
967  void
969  const unsigned int dim,
970  const Vector<double> & h,
973  & tmp_vectors,
974  const bool zero_out)
975  {
976  for (unsigned int j = 0; j < p.locally_owned_size(); ++j)
977  {
978  double temp = zero_out ? 0 : p.local_element(j);
979  for (unsigned int i = 0; i < dim; ++i)
980  temp += tmp_vectors[i].local_element(j) * h(i);
981  p.local_element(j) = temp;
982  }
983  }
984 
985 
986 
998  template <class VectorType>
999  inline double
1000  iterated_gram_schmidt(
1002  OrthogonalizationStrategy orthogonalization_strategy,
1004  & orthogonal_vectors,
1005  const unsigned int dim,
1006  const unsigned int accumulated_iterations,
1007  VectorType & vv,
1008  Vector<double> & h,
1009  bool & reorthogonalize,
1010  const boost::signals2::signal<void(int)> &reorthogonalize_signal =
1011  boost::signals2::signal<void(int)>())
1012  {
1013  Assert(dim > 0, ExcInternalError());
1014  const unsigned int inner_iteration = dim - 1;
1015 
1016  // need initial norm for detection of re-orthogonalization, see below
1017  double norm_vv_start = 0;
1018  const bool consider_reorthogonalize =
1019  (reorthogonalize == false) && (inner_iteration % 5 == 4);
1020  if (consider_reorthogonalize)
1021  norm_vv_start = vv.l2_norm();
1022 
1023  for (unsigned int i = 0; i < dim; ++i)
1024  h[i] = 0;
1025 
1026  for (unsigned int c = 0; c < 2;
1027  ++c) // 0: orthogonalize, 1: reorthogonalize
1028  {
1029  // Orthogonalization
1030  double norm_vv = 0.0;
1031 
1032  if (orthogonalization_strategy ==
1034  OrthogonalizationStrategy::modified_gram_schmidt)
1035  {
1036  double htmp = vv * orthogonal_vectors[0];
1037  h(0) += htmp;
1038  for (unsigned int i = 1; i < dim; ++i)
1039  {
1040  htmp = vv.add_and_dot(-htmp,
1041  orthogonal_vectors[i - 1],
1042  orthogonal_vectors[i]);
1043  h(i) += htmp;
1044  }
1045 
1046  norm_vv = std::sqrt(
1047  vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
1048  }
1049  else if (orthogonalization_strategy ==
1051  OrthogonalizationStrategy::classical_gram_schmidt)
1052  {
1053  Tvmult_add(dim, vv, orthogonal_vectors, h);
1054  norm_vv = substract_and_norm(dim, orthogonal_vectors, h, vv);
1055  }
1056  else
1057  {
1058  AssertThrow(false, ExcNotImplemented());
1059  }
1060 
1061  if (c == 1)
1062  return norm_vv; // reorthogonalization already performed -> finished
1063 
1064  // Re-orthogonalization if loss of orthogonality detected. For the
1065  // test, use a strategy discussed in C. T. Kelley, Iterative Methods
1066  // for Linear and Nonlinear Equations, SIAM, Philadelphia, 1995:
1067  // Compare the norm of vv after orthogonalization with its norm when
1068  // starting the orthogonalization. If vv became very small (here: less
1069  // than the square root of the machine precision times 10), it is
1070  // almost in the span of the previous vectors, which indicates loss of
1071  // precision.
1072  if (consider_reorthogonalize)
1073  {
1074  if (norm_vv >
1075  10. * norm_vv_start *
1076  std::sqrt(std::numeric_limits<
1077  typename VectorType::value_type>::epsilon()))
1078  return norm_vv;
1079 
1080  else
1081  {
1082  reorthogonalize = true;
1083  if (!reorthogonalize_signal.empty())
1084  reorthogonalize_signal(accumulated_iterations);
1085  }
1086  }
1087 
1088  if (reorthogonalize == false)
1089  return norm_vv; // no reorthogonalization needed -> finished
1090  }
1091 
1092  AssertThrow(false, ExcInternalError());
1093 
1094  return 0.0;
1095  }
1096  } // namespace SolverGMRESImplementation
1097 } // namespace internal
1098 
1099 
1100 
1101 template <class VectorType>
1102 inline void
1104  const FullMatrix<double> &H_orig,
1105  const unsigned int dim,
1106  const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
1107  &eigenvalues_signal,
1108  const boost::signals2::signal<void(const FullMatrix<double> &)>
1109  & hessenberg_signal,
1110  const boost::signals2::signal<void(double)> &cond_signal)
1111 {
1112  // Avoid copying the Hessenberg matrix if it isn't needed.
1113  if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1114  !cond_signal.empty()) &&
1115  dim > 0)
1116  {
1117  LAPACKFullMatrix<double> mat(dim, dim);
1118  for (unsigned int i = 0; i < dim; ++i)
1119  for (unsigned int j = 0; j < dim; ++j)
1120  mat(i, j) = H_orig(i, j);
1121  hessenberg_signal(H_orig);
1122  // Avoid computing eigenvalues if they are not needed.
1123  if (!eigenvalues_signal.empty())
1124  {
1125  // Copy mat so that we can compute svd below. Necessary since
1126  // compute_eigenvalues will leave mat in state
1127  // LAPACKSupport::unusable.
1128  LAPACKFullMatrix<double> mat_eig(mat);
1129  mat_eig.compute_eigenvalues();
1130  std::vector<std::complex<double>> eigenvalues(dim);
1131  for (unsigned int i = 0; i < mat_eig.n(); ++i)
1132  eigenvalues[i] = mat_eig.eigenvalue(i);
1133  // Sort eigenvalues for nicer output.
1134  std::sort(eigenvalues.begin(),
1135  eigenvalues.end(),
1136  internal::SolverGMRESImplementation::complex_less_pred);
1137  eigenvalues_signal(eigenvalues);
1138  }
1139  // Calculate condition number, avoid calculating the svd if a slot
1140  // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
1141  if (!cond_signal.empty() && (mat.n() > 1))
1142  {
1143  mat.compute_svd();
1144  double condition_number =
1145  mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1146  cond_signal(condition_number);
1147  }
1148  }
1149 }
1150 
1151 
1152 
1153 template <class VectorType>
1154 template <typename MatrixType, typename PreconditionerType>
1155 void
1156 SolverGMRES<VectorType>::solve(const MatrixType & A,
1157  VectorType & x,
1158  const VectorType & b,
1159  const PreconditionerType &preconditioner)
1160 {
1161  // TODO:[?] Check, why there are two different start residuals.
1162  // TODO:[GK] Make sure the parameter in the constructor means maximum basis
1163  // size
1164 
1165  std::unique_ptr<LogStream::Prefix> prefix;
1166  if (!additional_data.batched_mode)
1167  prefix = std::make_unique<LogStream::Prefix>("GMRES");
1168 
1169  // extra call to std::max to placate static analyzers: coverity rightfully
1170  // complains that data.max_n_tmp_vectors - 2 may overflow
1171  const unsigned int n_tmp_vectors =
1172  std::max(additional_data.max_n_tmp_vectors, 3u);
1173 
1174  // Generate an object where basis vectors are stored.
1176  n_tmp_vectors, this->memory);
1177 
1178  // number of the present iteration; this
1179  // number is not reset to zero upon a
1180  // restart
1181  unsigned int accumulated_iterations = 0;
1182 
1183  const bool do_eigenvalues =
1184  !additional_data.batched_mode &&
1185  (!condition_number_signal.empty() ||
1186  !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1187  !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1188  !all_hessenberg_signal.empty());
1189  // for eigenvalue computation, need to collect the Hessenberg matrix (before
1190  // applying Givens rotations)
1191  FullMatrix<double> H_orig;
1192  if (do_eigenvalues)
1193  H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
1194 
1195  // matrix used for the orthogonalization process later
1196  H.reinit(n_tmp_vectors, n_tmp_vectors - 1, /* omit_initialization */ true);
1197 
1198  // some additional vectors, also used in the orthogonalization
1199  gamma.reinit(n_tmp_vectors);
1200  ci.reinit(n_tmp_vectors - 1);
1201  si.reinit(n_tmp_vectors - 1);
1202  h.reinit(n_tmp_vectors - 1);
1203 
1204  unsigned int dim = 0;
1205 
1206  SolverControl::State iteration_state = SolverControl::iterate;
1207  double last_res = std::numeric_limits<double>::lowest();
1208 
1209  // switch to determine whether we want a left or a right preconditioner. at
1210  // present, left is default, but both ways are implemented
1211  const bool left_precondition = !additional_data.right_preconditioning;
1212 
1213  // Per default the left preconditioned GMRes uses the preconditioned
1214  // residual and the right preconditioned GMRes uses the unpreconditioned
1215  // residual as stopping criterion.
1216  const bool use_default_residual = additional_data.use_default_residual;
1217 
1218  // define two aliases
1219  VectorType &v = tmp_vectors(0, x);
1220  VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
1221 
1222  // Following vectors are needed when we are not using the default residuals
1223  // as stopping criterion
1226  std::unique_ptr<::Vector<double>> gamma_;
1227  if (!use_default_residual)
1228  {
1229  r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1230  x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1231  r->reinit(x);
1232  x_->reinit(x);
1233 
1234  gamma_ = std::make_unique<::Vector<double>>(gamma.size());
1235  }
1236 
1237  bool re_orthogonalize = additional_data.force_re_orthogonalization;
1238 
1240  // outer iteration: loop until we either reach convergence or the maximum
1241  // number of iterations is exceeded. each cycle of this loop amounts to one
1242  // restart
1243  do
1244  {
1245  // reset this vector to the right size
1246  h.reinit(n_tmp_vectors - 1);
1247 
1248  double rho = 0.0;
1249 
1250  if (left_precondition)
1251  {
1252  A.vmult(p, x);
1253  p.sadd(-1., 1., b);
1254  preconditioner.vmult(v, p);
1255  rho = v.l2_norm();
1256  }
1257  else
1258  {
1259  A.vmult(v, x);
1260  rho = ::internal::SolverGMRESImplementation::sadd_and_norm(v,
1261  -1,
1262  b,
1263  1.0);
1264  }
1265 
1266  // check the residual here as well since it may be that we got the exact
1267  // (or an almost exact) solution vector at the outset. if we wouldn't
1268  // check here, the next scaling operation would produce garbage
1269  if (use_default_residual)
1270  {
1271  last_res = rho;
1272  if (additional_data.batched_mode)
1273  iteration_state = solver_control.check(accumulated_iterations, rho);
1274  else
1275  iteration_state =
1276  this->iteration_status(accumulated_iterations, rho, x);
1277 
1278  if (iteration_state != SolverControl::iterate)
1279  break;
1280  }
1281  else
1282  {
1283  deallog << "default_res=" << rho << std::endl;
1284 
1285  if (left_precondition)
1286  {
1287  A.vmult(*r, x);
1288  r->sadd(-1., 1., b);
1289  }
1290  else
1291  preconditioner.vmult(*r, v);
1292 
1293  double res = r->l2_norm();
1294  last_res = res;
1295  if (additional_data.batched_mode)
1296  iteration_state = solver_control.check(accumulated_iterations, rho);
1297  else
1298  iteration_state =
1299  this->iteration_status(accumulated_iterations, res, x);
1300 
1301  if (iteration_state != SolverControl::iterate)
1302  break;
1303  }
1304 
1305  gamma(0) = rho;
1306 
1307  v *= 1. / rho;
1308 
1309  // inner iteration doing at most as many steps as there are temporary
1310  // vectors. the number of steps actually been done is propagated outside
1311  // through the @p dim variable
1312  for (unsigned int inner_iteration = 0;
1313  ((inner_iteration < n_tmp_vectors - 2) &&
1314  (iteration_state == SolverControl::iterate));
1315  ++inner_iteration)
1316  {
1317  ++accumulated_iterations;
1318  // yet another alias
1319  VectorType &vv = tmp_vectors(inner_iteration + 1, x);
1320 
1321  if (left_precondition)
1322  {
1323  A.vmult(p, tmp_vectors[inner_iteration]);
1324  preconditioner.vmult(vv, p);
1325  }
1326  else
1327  {
1328  preconditioner.vmult(p, tmp_vectors[inner_iteration]);
1329  A.vmult(vv, p);
1330  }
1331 
1332  dim = inner_iteration + 1;
1333 
1334  const double s =
1335  internal::SolverGMRESImplementation::iterated_gram_schmidt(
1336  additional_data.orthogonalization_strategy,
1337  tmp_vectors,
1338  dim,
1339  accumulated_iterations,
1340  vv,
1341  h,
1342  re_orthogonalize,
1343  re_orthogonalize_signal);
1344  h(inner_iteration + 1) = s;
1345 
1346  // s=0 is a lucky breakdown, the solver will reach convergence,
1347  // but we must not divide by zero here.
1348  if (s != 0)
1349  vv *= 1. / s;
1350 
1351  // for eigenvalues, get the resulting coefficients from the
1352  // orthogonalization process
1353  if (do_eigenvalues)
1354  for (unsigned int i = 0; i < dim + 1; ++i)
1355  H_orig(i, inner_iteration) = h(i);
1356 
1357  // Transformation into tridiagonal structure
1358  givens_rotation(h, gamma, ci, si, inner_iteration);
1359 
1360  // append vector on matrix
1361  for (unsigned int i = 0; i < dim; ++i)
1362  H(i, inner_iteration) = h(i);
1363 
1364  // default residual
1365  rho = std::fabs(gamma(dim));
1366 
1367  if (use_default_residual)
1368  {
1369  last_res = rho;
1370  if (additional_data.batched_mode)
1371  iteration_state =
1372  solver_control.check(accumulated_iterations, rho);
1373  else
1374  iteration_state =
1375  this->iteration_status(accumulated_iterations, rho, x);
1376  }
1377  else
1378  {
1379  if (!additional_data.batched_mode)
1380  deallog << "default_res=" << rho << std::endl;
1381 
1382  *x_ = x;
1383  *gamma_ = gamma;
1384  internal::SolverGMRESImplementation::solve_triangular(dim,
1385  H,
1386  *gamma_,
1387  h);
1388 
1389  if (left_precondition)
1390  for (unsigned int i = 0; i < dim; ++i)
1391  x_->add(h(i), tmp_vectors[i]);
1392  else
1393  {
1394  p = 0.;
1395  for (unsigned int i = 0; i < dim; ++i)
1396  p.add(h(i), tmp_vectors[i]);
1397  preconditioner.vmult(*r, p);
1398  x_->add(1., *r);
1399  };
1400  A.vmult(*r, *x_);
1401  r->sadd(-1., 1., b);
1402  // Now *r contains the unpreconditioned residual!!
1403  if (left_precondition)
1404  {
1405  const double res = r->l2_norm();
1406  last_res = res;
1407 
1408  iteration_state =
1409  this->iteration_status(accumulated_iterations, res, x);
1410  }
1411  else
1412  {
1413  preconditioner.vmult(*x_, *r);
1414  const double preconditioned_res = x_->l2_norm();
1415  last_res = preconditioned_res;
1416 
1417  if (additional_data.batched_mode)
1418  iteration_state =
1419  solver_control.check(accumulated_iterations, rho);
1420  else
1421  iteration_state =
1422  this->iteration_status(accumulated_iterations,
1423  preconditioned_res,
1424  x);
1425  }
1426  }
1427  }
1428 
1429  // end of inner iteration. now calculate the solution from the temporary
1430  // vectors
1431  internal::SolverGMRESImplementation::solve_triangular(dim, H, gamma, h);
1432 
1433  if (do_eigenvalues)
1434  compute_eigs_and_cond(H_orig,
1435  dim,
1436  all_eigenvalues_signal,
1437  all_hessenberg_signal,
1438  condition_number_signal);
1439 
1440  if (left_precondition)
1441  ::internal::SolverGMRESImplementation::add(
1442  x, dim, h, tmp_vectors, false);
1443  else
1444  {
1445  ::internal::SolverGMRESImplementation::add(
1446  p, dim, h, tmp_vectors, true);
1447  preconditioner.vmult(v, p);
1448  x.add(1., v);
1449  };
1450  // end of outer iteration. restart if no convergence and the number of
1451  // iterations is not exceeded
1452  }
1453  while (iteration_state == SolverControl::iterate);
1454 
1455  if (do_eigenvalues)
1456  compute_eigs_and_cond(H_orig,
1457  dim,
1458  eigenvalues_signal,
1459  hessenberg_signal,
1460  condition_number_signal);
1461 
1462  if (!additional_data.batched_mode && !krylov_space_signal.empty())
1463  krylov_space_signal(tmp_vectors);
1464 
1465  // in case of failure: throw exception
1466  AssertThrow(iteration_state == SolverControl::success,
1467  SolverControl::NoConvergence(accumulated_iterations, last_res));
1468 }
1469 
1470 
1471 
1472 template <class VectorType>
1473 boost::signals2::connection
1475  const std::function<void(double)> &slot,
1476  const bool every_iteration)
1477 {
1478  if (every_iteration)
1479  {
1480  return all_condition_numbers_signal.connect(slot);
1481  }
1482  else
1483  {
1484  return condition_number_signal.connect(slot);
1485  }
1486 }
1487 
1488 
1489 
1490 template <class VectorType>
1491 boost::signals2::connection
1493  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1494  const bool every_iteration)
1495 {
1496  if (every_iteration)
1497  {
1498  return all_eigenvalues_signal.connect(slot);
1499  }
1500  else
1501  {
1502  return eigenvalues_signal.connect(slot);
1503  }
1504 }
1505 
1506 
1507 
1508 template <class VectorType>
1509 boost::signals2::connection
1511  const std::function<void(const FullMatrix<double> &)> &slot,
1512  const bool every_iteration)
1513 {
1514  if (every_iteration)
1515  {
1516  return all_hessenberg_signal.connect(slot);
1517  }
1518  else
1519  {
1520  return hessenberg_signal.connect(slot);
1521  }
1522 }
1523 
1524 
1525 
1526 template <class VectorType>
1527 boost::signals2::connection
1529  const std::function<void(
1531 {
1532  return krylov_space_signal.connect(slot);
1533 }
1534 
1535 
1536 
1537 template <class VectorType>
1538 boost::signals2::connection
1540  const std::function<void(int)> &slot)
1541 {
1542  return re_orthogonalize_signal.connect(slot);
1543 }
1544 
1545 
1546 
1547 template <class VectorType>
1548 double
1550 {
1551  // dummy implementation. this function is not needed for the present
1552  // implementation of gmres
1553  Assert(false, ExcInternalError());
1554  return 0;
1555 }
1556 
1557 
1558 //----------------------------------------------------------------------//
1559 
1560 template <class VectorType>
1563  const AdditionalData & data)
1564  : SolverBase<VectorType>(cn, mem)
1565  , additional_data(data)
1566 {}
1567 
1568 
1569 
1570 template <class VectorType>
1572  const AdditionalData &data)
1573  : SolverBase<VectorType>(cn)
1574  , additional_data(data)
1575 {}
1576 
1577 
1578 
1579 template <class VectorType>
1580 template <typename MatrixType, typename PreconditionerType>
1581 void
1582 SolverFGMRES<VectorType>::solve(const MatrixType & A,
1583  VectorType & x,
1584  const VectorType & b,
1585  const PreconditionerType &preconditioner)
1586 {
1587  LogStream::Prefix prefix("FGMRES");
1588 
1589  SolverControl::State iteration_state = SolverControl::iterate;
1590 
1591  const unsigned int basis_size = additional_data.max_basis_size;
1592 
1593  // Generate an object where basis vectors are stored.
1595  basis_size, this->memory);
1597  basis_size, this->memory);
1598 
1599  // number of the present iteration; this number is not reset to zero upon a
1600  // restart
1601  unsigned int accumulated_iterations = 0;
1602 
1603  // matrix used for the orthogonalization process later
1604  H.reinit(basis_size + 1, basis_size);
1605 
1606  // Vectors for projected system
1607  Vector<double> projected_rhs;
1608  Vector<double> y;
1609 
1610  // Iteration starts here
1611  double res = std::numeric_limits<double>::lowest();
1612 
1613  typename VectorMemory<VectorType>::Pointer aux(this->memory);
1614  aux->reinit(x);
1615  do
1616  {
1617  A.vmult(*aux, x);
1618  aux->sadd(-1., 1., b);
1619 
1620  double beta = aux->l2_norm();
1621  res = beta;
1622  iteration_state = this->iteration_status(accumulated_iterations, res, x);
1623  if (iteration_state == SolverControl::success)
1624  break;
1625 
1626  H.reinit(basis_size + 1, basis_size);
1627  double a = beta;
1628 
1629  for (unsigned int j = 0; j < basis_size; ++j)
1630  {
1631  if (a != 0) // treat lucky breakdown
1632  v(j, x).equ(1. / a, *aux);
1633  else
1634  v(j, x) = 0.;
1635 
1636 
1637  preconditioner.vmult(z(j, x), v[j]);
1638  A.vmult(*aux, z[j]);
1639 
1640  // Gram-Schmidt
1641  H(0, j) = *aux * v[0];
1642  for (unsigned int i = 1; i <= j; ++i)
1643  H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1644  H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1645 
1646  // Compute projected solution
1647 
1648  if (j > 0)
1649  {
1650  H1.reinit(j + 1, j);
1651  projected_rhs.reinit(j + 1);
1652  y.reinit(j);
1653  projected_rhs(0) = beta;
1654  H1.fill(H);
1655 
1656  // check convergence. note that the vector 'x' we pass to the
1657  // criterion is not the final solution we compute if we
1658  // decide to jump out of the iteration (we update 'x' again
1659  // right after the current loop)
1660  Householder<double> house(H1);
1661  res = house.least_squares(y, projected_rhs);
1662  iteration_state =
1663  this->iteration_status(++accumulated_iterations, res, x);
1664  if (iteration_state != SolverControl::iterate)
1665  break;
1666  }
1667  }
1668 
1669  // Update solution vector
1670  for (unsigned int j = 0; j < y.size(); ++j)
1671  x.add(y(j), z[j]);
1672  }
1673  while (iteration_state == SolverControl::iterate);
1674 
1675  // in case of failure: throw exception
1676  if (iteration_state != SolverControl::success)
1677  AssertThrow(false,
1678  SolverControl::NoConvergence(accumulated_iterations, res));
1679 }
1680 
1681 #endif // DOXYGEN
1682 
1684 
1685 #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:565
SolverFGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
Definition: solver_gmres.h:560
FullMatrix< double > H1
Definition: solver_gmres.h:570
boost::signals2::signal< void(const FullMatrix< double > &)> hessenberg_signal
Definition: solver_gmres.h:399
Vector< double > si
Definition: solver_gmres.h:480
AdditionalData additional_data
Definition: solver_gmres.h:367
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_gmres.h:373
SolverGMRES(const SolverGMRES< VectorType > &)=delete
boost::signals2::signal< void(const std::vector< std::complex< double >> &)> eigenvalues_signal
Definition: solver_gmres.h:386
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:414
boost::signals2::signal< void(const FullMatrix< double > &)> all_hessenberg_signal
Definition: solver_gmres.h:406
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:475
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
SolverControl & solver_control
Definition: solver_gmres.h:427
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:470
boost::signals2::signal< void(const std::vector< std::complex< double >> &)> all_eigenvalues_signal
Definition: solver_gmres.h:393
Vector< double > h
Definition: solver_gmres.h:485
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
FullMatrix< double > H
Definition: solver_gmres.h:465
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:379
boost::signals2::signal< void(int)> re_orthogonalize_signal
Definition: solver_gmres.h:420
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:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1583
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1847
static ::ExceptionBase & ExcTooFewTmpVectors(int arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1821
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:510
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1672
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:522
OrthogonalizationStrategy orthogonalization_strategy
Definition: solver_gmres.h:272
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)