Reference documentation for deal.II version GIT ca9e7e8105 2023-03-24 14:15:03+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 
31 #include <deal.II/lac/solver.h>
33 #include <deal.II/lac/vector.h>
34 
35 #include <algorithm>
36 #include <cmath>
37 #include <limits>
38 #include <memory>
39 #include <vector>
40 
42 
43 // forward declarations
44 #ifndef DOXYGEN
45 namespace LinearAlgebra
46 {
47  namespace distributed
48  {
49  template <typename, typename>
50  class Vector;
51  } // namespace distributed
52 } // namespace LinearAlgebra
53 #endif
54 
60 namespace internal
61 {
65  namespace SolverGMRESImplementation
66  {
75  template <typename VectorType>
76  class TmpVectors
77  {
78  public:
83  TmpVectors(const unsigned int max_size, VectorMemory<VectorType> &vmem);
84 
88  ~TmpVectors() = default;
89 
94  VectorType &
95  operator[](const unsigned int i) const;
96 
103  VectorType &
104  operator()(const unsigned int i, const VectorType &temp);
105 
110  unsigned int
111  size() const;
112 
113 
114  private:
119 
123  std::vector<typename VectorMemory<VectorType>::Pointer> data;
124  };
125  } // namespace SolverGMRESImplementation
126 } // namespace internal
127 
192 template <class VectorType = Vector<double>>
193 class SolverGMRES : public SolverBase<VectorType>
194 {
195 public:
200  {
202  {
214  };
215 
223  explicit AdditionalData(
224  const unsigned int max_n_tmp_vectors = 30,
225  const bool right_preconditioning = false,
226  const bool use_default_residual = true,
227  const bool force_re_orthogonalization = false,
228  const bool batched_mode = false,
231 
238  unsigned int max_n_tmp_vectors;
239 
248 
253 
261 
269 
274  };
275 
281  const AdditionalData & data = AdditionalData());
282 
288 
293 
297  template <typename MatrixType, typename PreconditionerType>
298  void
299  solve(const MatrixType & A,
300  VectorType & x,
301  const VectorType & b,
302  const PreconditionerType &preconditioner);
303 
310  boost::signals2::connection
311  connect_condition_number_slot(const std::function<void(double)> &slot,
312  const bool every_iteration = false);
313 
320  boost::signals2::connection
322  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
323  const bool every_iteration = false);
324 
332  boost::signals2::connection
334  const std::function<void(const FullMatrix<double> &)> &slot,
335  const bool every_iteration = true);
336 
343  boost::signals2::connection
345  const std::function<
347  &slot);
348 
349 
354  boost::signals2::connection
355  connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
356 
357 
359  int,
360  << "The number of temporary vectors you gave (" << arg1
361  << ") is too small. It should be at least 10 for "
362  << "any results, and much more for reasonable ones.");
363 
364 protected:
369 
374  boost::signals2::signal<void(double)> condition_number_signal;
375 
380  boost::signals2::signal<void(double)> all_condition_numbers_signal;
381 
386  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
388 
393  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
395 
400  boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
401 
406  boost::signals2::signal<void(const FullMatrix<double> &)>
408 
413  boost::signals2::signal<void(
416 
421  boost::signals2::signal<void(int)> re_orthogonalize_signal;
422 
429 
433  virtual double
435 
440  void
442  Vector<double> &b,
445  int col) const;
446 
453  static void
455  const FullMatrix<double> &H_orig,
456  const unsigned int dim,
457  const boost::signals2::signal<
458  void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
459  const boost::signals2::signal<void(const FullMatrix<double> &)>
461  const boost::signals2::signal<void(double)> &cond_signal);
462 
467 
472 
477 
482 
487 };
488 
489 
490 
511 template <class VectorType = Vector<double>>
512 class SolverFGMRES : public SolverBase<VectorType>
513 {
514 public:
519  {
523  explicit AdditionalData(const unsigned int max_basis_size = 30)
525  {}
526 
530  unsigned int max_basis_size;
531  };
532 
538  const AdditionalData & data = AdditionalData());
539 
545  const AdditionalData &data = AdditionalData());
546 
550  template <typename MatrixType, typename PreconditionerType>
551  void
552  solve(const MatrixType & A,
553  VectorType & x,
554  const VectorType & b,
555  const PreconditionerType &preconditioner);
556 
557 private:
562 
567 
572 };
573 
575 /* --------------------- Inline and template functions ------------------- */
576 
577 
578 #ifndef DOXYGEN
579 namespace internal
580 {
581  namespace SolverGMRESImplementation
582  {
583  template <class VectorType>
584  inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
586  : mem(vmem)
587  , data(max_size)
588  {}
589 
590 
591 
592  template <class VectorType>
593  inline VectorType &
594  TmpVectors<VectorType>::operator[](const unsigned int i) const
595  {
596  AssertIndexRange(i, data.size());
597 
598  Assert(data[i] != nullptr, ExcNotInitialized());
599  return *data[i];
600  }
601 
602 
603 
604  template <class VectorType>
605  inline VectorType &
606  TmpVectors<VectorType>::operator()(const unsigned int i,
607  const VectorType & temp)
608  {
609  AssertIndexRange(i, data.size());
610  if (data[i] == nullptr)
611  {
612  data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
613  data[i]->reinit(temp, true);
614  }
615  return *data[i];
616  }
617 
618 
619 
620  template <class VectorType>
621  unsigned int
623  {
624  return (data.size() > 0 ? data.size() - 1 : 0);
625  }
626 
627 
628 
629  // A comparator for better printing eigenvalues
630  inline bool
631  complex_less_pred(const std::complex<double> &x,
632  const std::complex<double> &y)
633  {
634  return x.real() < y.real() ||
635  (x.real() == y.real() && x.imag() < y.imag());
636  }
637 
638  // A function to solve the (upper) triangular system after Givens
639  // rotations on a matrix that has possibly unused rows and columns
640  inline void
641  solve_triangular(const unsigned int dim,
642  const FullMatrix<double> &H,
643  const Vector<double> & rhs,
644  Vector<double> & solution)
645  {
646  for (int i = dim - 1; i >= 0; --i)
647  {
648  double s = rhs(i);
649  for (unsigned int j = i + 1; j < dim; ++j)
650  s -= solution(j) * H(i, j);
651  solution(i) = s / H(i, i);
652  AssertIsFinite(solution(i));
653  }
654  }
655  } // namespace SolverGMRESImplementation
656 } // namespace internal
657 
658 
659 
660 template <class VectorType>
662  const unsigned int max_n_tmp_vectors,
663  const bool right_preconditioning,
664  const bool use_default_residual,
665  const bool force_re_orthogonalization,
666  const bool batched_mode,
667  const OrthogonalizationStrategy orthogonalization_strategy)
668  : max_n_tmp_vectors(max_n_tmp_vectors)
669  , right_preconditioning(right_preconditioning)
670  , use_default_residual(use_default_residual)
671  , force_re_orthogonalization(force_re_orthogonalization)
672  , batched_mode(batched_mode)
673  , orthogonalization_strategy(orthogonalization_strategy)
674 {
675  Assert(3 <= max_n_tmp_vectors,
676  ExcMessage("SolverGMRES needs at least three "
677  "temporary vectors."));
678 }
679 
680 
681 
682 template <class VectorType>
685  const AdditionalData & data)
686  : SolverBase<VectorType>(cn, mem)
687  , additional_data(data)
688  , solver_control(cn)
689 {}
690 
691 
692 
693 template <class VectorType>
695  const AdditionalData &data)
696  : SolverBase<VectorType>(cn)
697  , additional_data(data)
698  , solver_control(cn)
699 {}
700 
701 
702 
703 template <class VectorType>
704 inline void
706  Vector<double> &b,
707  Vector<double> &ci,
708  Vector<double> &si,
709  int col) const
710 {
711  for (int i = 0; i < col; ++i)
712  {
713  const double s = si(i);
714  const double c = ci(i);
715  const double dummy = h(i);
716  h(i) = c * dummy + s * h(i + 1);
717  h(i + 1) = -s * dummy + c * h(i + 1);
718  };
719 
720  const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
721  si(col) = h(col + 1) * r;
722  ci(col) = h(col) * r;
723  h(col) = ci(col) * h(col) + si(col) * h(col + 1);
724  b(col + 1) = -si(col) * b(col);
725  b(col) *= ci(col);
726 }
727 
728 
729 
730 namespace internal
731 {
732  namespace SolverGMRESImplementation
733  {
734  template <typename VectorType, typename Enable = void>
735  struct is_dealii_compatible_distributed_vector;
736 
737  template <typename VectorType>
738  struct is_dealii_compatible_distributed_vector<
739  VectorType,
740  typename std::enable_if<!internal::is_block_vector<VectorType>>::type>
741  {
742  static constexpr bool value = std::is_same<
743  VectorType,
744  LinearAlgebra::distributed::Vector<typename VectorType::value_type,
745  MemorySpace::Host>>::value;
746  };
747 
748 
749 
750  template <typename VectorType>
751  struct is_dealii_compatible_distributed_vector<
752  VectorType,
753  typename std::enable_if<internal::is_block_vector<VectorType>>::type>
754  {
755  static constexpr bool value = std::is_same<
756  typename VectorType::BlockType,
757  LinearAlgebra::distributed::Vector<typename VectorType::value_type,
758  MemorySpace::Host>>::value;
759  };
760 
761 
762 
763  template <typename VectorType,
764  std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
765  * = nullptr>
766  unsigned int
767  n_blocks(const VectorType &)
768  {
769  return 1;
770  }
771 
772 
773 
774  template <typename VectorType,
775  std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
776  nullptr>
777  unsigned int
778  n_blocks(const VectorType &vector)
779  {
780  return vector.n_blocks();
781  }
782 
783 
784 
785  template <typename VectorType,
786  std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
787  * = nullptr>
788  VectorType &
789  block(VectorType &vector, const unsigned int b)
790  {
791  AssertDimension(b, 0);
792  (void)b;
793  return vector;
794  }
795 
796 
797 
798  template <typename VectorType,
799  std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
800  * = nullptr>
801  const VectorType &
802  block(const VectorType &vector, const unsigned int b)
803  {
804  AssertDimension(b, 0);
805  (void)b;
806  return vector;
807  }
808 
809 
810 
811  template <typename VectorType,
812  std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
813  nullptr>
814  typename VectorType::BlockType &
815  block(VectorType &vector, const unsigned int b)
816  {
817  return vector.block(b);
818  }
819 
820 
821 
822  template <typename VectorType,
823  std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
824  nullptr>
825  const typename VectorType::BlockType &
826  block(const VectorType &vector, const unsigned int b)
827  {
828  return vector.block(b);
829  }
830 
831 
832 
833  template <class VectorType,
834  std::enable_if_t<
835  !is_dealii_compatible_distributed_vector<VectorType>::value,
836  VectorType> * = nullptr>
837  void
838  Tvmult_add(const unsigned int dim,
839  const VectorType & vv,
841  & orthogonal_vectors,
842  Vector<double> &h)
843  {
844  for (unsigned int i = 0; i < dim; ++i)
845  h[i] += vv * orthogonal_vectors[i];
846  }
847 
848 
849 
850  template <class VectorType,
851  std::enable_if_t<
852  is_dealii_compatible_distributed_vector<VectorType>::value,
853  VectorType> * = nullptr>
854  void
855  Tvmult_add(const unsigned int dim,
856  const VectorType & vv,
858  & orthogonal_vectors,
859  Vector<double> &h)
860  {
861  for (unsigned int b = 0; b < n_blocks(vv); ++b)
862  {
863  unsigned int j = 0;
864 
865  if (dim <= 128)
866  {
867  // optimized path
868  static constexpr unsigned int n_lanes =
870 
871  VectorizedArray<double> hs[128];
872  for (unsigned int d = 0; d < dim; ++d)
873  hs[d] = 0.0;
874 
875  unsigned int c = 0;
876 
877  for (; c < block(vv, b).locally_owned_size() / n_lanes / 4;
878  ++c, j += n_lanes * 4)
879  for (unsigned int i = 0; i < dim; ++i)
880  {
881  VectorizedArray<double> vvec[4];
882  for (unsigned int k = 0; k < 4; ++k)
883  vvec[k].load(block(vv, b).begin() + j + k * n_lanes);
884 
885  for (unsigned int k = 0; k < 4; ++k)
886  {
888  temp.load(block(orthogonal_vectors[i], b).begin() + j +
889  k * n_lanes);
890  hs[i] += temp * vvec[k];
891  }
892  }
893 
894  c *= 4;
895  for (; c < block(vv, b).locally_owned_size() / n_lanes;
896  ++c, j += n_lanes)
897  for (unsigned int i = 0; i < dim; ++i)
898  {
899  VectorizedArray<double> vvec, temp;
900  vvec.load(block(vv, b).begin() + j);
901  temp.load(block(orthogonal_vectors[i], b).begin() + j);
902  hs[i] += temp * vvec;
903  }
904 
905  for (unsigned int i = 0; i < dim; ++i)
906  for (unsigned int v = 0; v < n_lanes; ++v)
907  h(i) += hs[i][v];
908  }
909 
910  // remainder loop of optimized path or non-optimized path (if
911  // dim>128)
912  for (; j < block(vv, b).locally_owned_size(); ++j)
913  for (unsigned int i = 0; i < dim; ++i)
914  h(i) += block(orthogonal_vectors[i], b).local_element(j) *
915  block(vv, b).local_element(j);
916  }
917 
918  Utilities::MPI::sum(h, MPI_COMM_WORLD, h);
919  }
920 
921 
922 
923  template <class VectorType,
924  std::enable_if_t<
925  !is_dealii_compatible_distributed_vector<VectorType>::value,
926  VectorType> * = nullptr>
927  double
928  substract_and_norm(
929  const unsigned int dim,
931  & orthogonal_vectors,
932  const Vector<double> &h,
933  VectorType & vv)
934  {
935  Assert(dim > 0, ExcInternalError());
936 
937  for (unsigned int i = 0; i < dim; ++i)
938  vv.add(-h(i), orthogonal_vectors[i]);
939 
940  return std::sqrt(vv.add_and_dot(-h(dim), orthogonal_vectors[dim], vv));
941  }
942 
943 
944 
945  template <class VectorType,
946  std::enable_if_t<
947  is_dealii_compatible_distributed_vector<VectorType>::value,
948  VectorType> * = nullptr>
949  double
950  substract_and_norm(
951  const unsigned int dim,
953  & orthogonal_vectors,
954  const Vector<double> &h,
955  VectorType & vv)
956  {
957  static constexpr unsigned int n_lanes = VectorizedArray<double>::size();
958 
959  double norm_vv_temp = 0.0;
960 
961  for (unsigned int b = 0; b < n_blocks(vv); ++b)
962  {
963  VectorizedArray<double> norm_vv_temp_vectorized = 0.0;
964 
965  unsigned int j = 0;
966  unsigned int c = 0;
967  for (; c < block(vv, b).locally_owned_size() / n_lanes / 4;
968  ++c, j += n_lanes * 4)
969  {
970  VectorizedArray<double> temp[4];
971 
972  for (unsigned int k = 0; k < 4; ++k)
973  temp[k].load(block(vv, b).begin() + j + k * n_lanes);
974 
975  for (unsigned int i = 0; i < dim; ++i)
976  {
977  const double factor = h(i);
978  for (unsigned int k = 0; k < 4; ++k)
979  {
981  vec.load(block(orthogonal_vectors[i], b).begin() + j +
982  k * n_lanes);
983  temp[k] -= factor * vec;
984  }
985  }
986 
987  for (unsigned int k = 0; k < 4; ++k)
988  temp[k].store(block(vv, b).begin() + j + k * n_lanes);
989 
990  norm_vv_temp_vectorized +=
991  (temp[0] * temp[0] + temp[1] * temp[1]) +
992  (temp[2] * temp[2] + temp[3] * temp[3]);
993  }
994 
995  c *= 4;
996  for (; c < block(vv, b).locally_owned_size() / n_lanes;
997  ++c, j += n_lanes)
998  {
1000  temp.load(block(vv, b).begin() + j);
1001 
1002  for (unsigned int i = 0; i < dim; ++i)
1003  {
1005  vec.load(block(orthogonal_vectors[i], b).begin() + j);
1006  temp -= h(i) * vec;
1007  }
1008 
1009  temp.store(block(vv, b).begin() + j);
1010 
1011  norm_vv_temp_vectorized += temp * temp;
1012  }
1013 
1014  for (unsigned int v = 0; v < n_lanes; ++v)
1015  norm_vv_temp += norm_vv_temp_vectorized[v];
1016 
1017  for (; j < block(vv, b).locally_owned_size(); ++j)
1018  {
1019  double temp = block(vv, b).local_element(j);
1020  for (unsigned int i = 0; i < dim; ++i)
1021  temp -= h(i) * block(orthogonal_vectors[i], b).local_element(j);
1022  block(vv, b).local_element(j) = temp;
1023 
1024  norm_vv_temp += temp * temp;
1025  }
1026  }
1027 
1028  return std::sqrt(Utilities::MPI::sum(norm_vv_temp, MPI_COMM_WORLD));
1029  }
1030 
1031 
1032  template <class VectorType,
1033  std::enable_if_t<
1034  !is_dealii_compatible_distributed_vector<VectorType>::value,
1035  VectorType> * = nullptr>
1036  double
1037  sadd_and_norm(VectorType & v,
1038  const double factor_a,
1039  const VectorType &b,
1040  const double factor_b)
1041  {
1042  v.sadd(factor_a, factor_b, b);
1043  return v.l2_norm();
1044  }
1045 
1046 
1047  template <class VectorType,
1048  std::enable_if_t<
1049  is_dealii_compatible_distributed_vector<VectorType>::value,
1050  VectorType> * = nullptr>
1051  double
1052  sadd_and_norm(VectorType & v,
1053  const double factor_a,
1054  const VectorType &w,
1055  const double factor_b)
1056  {
1057  double norm = 0;
1058 
1059  for (unsigned int b = 0; b < n_blocks(v); ++b)
1060  for (unsigned int j = 0; j < block(v, b).locally_owned_size(); ++j)
1061  {
1062  const double temp = block(v, b).local_element(j) * factor_a +
1063  block(w, b).local_element(j) * factor_b;
1064 
1065  block(v, b).local_element(j) = temp;
1066 
1067  norm += temp * temp;
1068  }
1069 
1070  return std::sqrt(Utilities::MPI::sum(norm, MPI_COMM_WORLD));
1071  }
1072 
1073 
1074 
1075  template <class VectorType,
1076  std::enable_if_t<
1077  !is_dealii_compatible_distributed_vector<VectorType>::value,
1078  VectorType> * = nullptr>
1079  void
1080  add(VectorType & p,
1081  const unsigned int dim,
1082  const Vector<double> &h,
1084  & tmp_vectors,
1085  const bool zero_out)
1086  {
1087  if (zero_out)
1088  p.equ(h(0), tmp_vectors[0]);
1089  else
1090  p.add(h(0), tmp_vectors[0]);
1091 
1092  for (unsigned int i = 1; i < dim; ++i)
1093  p.add(h(i), tmp_vectors[i]);
1094  }
1095 
1096 
1097 
1098  template <class VectorType,
1099  std::enable_if_t<
1100  is_dealii_compatible_distributed_vector<VectorType>::value,
1101  VectorType> * = nullptr>
1102  void
1103  add(VectorType & p,
1104  const unsigned int dim,
1105  const Vector<double> &h,
1107  & tmp_vectors,
1108  const bool zero_out)
1109  {
1110  for (unsigned int b = 0; b < n_blocks(p); ++b)
1111  for (unsigned int j = 0; j < block(p, b).locally_owned_size(); ++j)
1112  {
1113  double temp = zero_out ? 0 : block(p, b).local_element(j);
1114  for (unsigned int i = 0; i < dim; ++i)
1115  temp += block(tmp_vectors[i], b).local_element(j) * h(i);
1116  block(p, b).local_element(j) = temp;
1117  }
1118  }
1119 
1120 
1121 
1133  template <class VectorType>
1134  inline double
1135  iterated_gram_schmidt(
1137  OrthogonalizationStrategy orthogonalization_strategy,
1139  & orthogonal_vectors,
1140  const unsigned int dim,
1141  const unsigned int accumulated_iterations,
1142  VectorType & vv,
1143  Vector<double> & h,
1144  bool & reorthogonalize,
1145  const boost::signals2::signal<void(int)> &reorthogonalize_signal =
1146  boost::signals2::signal<void(int)>())
1147  {
1148  Assert(dim > 0, ExcInternalError());
1149  const unsigned int inner_iteration = dim - 1;
1150 
1151  // need initial norm for detection of re-orthogonalization, see below
1152  double norm_vv_start = 0;
1153  const bool consider_reorthogonalize =
1154  (reorthogonalize == false) && (inner_iteration % 5 == 4);
1155  if (consider_reorthogonalize)
1156  norm_vv_start = vv.l2_norm();
1157 
1158  for (unsigned int i = 0; i < dim; ++i)
1159  h[i] = 0;
1160 
1161  for (unsigned int c = 0; c < 2;
1162  ++c) // 0: orthogonalize, 1: reorthogonalize
1163  {
1164  // Orthogonalization
1165  double norm_vv = 0.0;
1166 
1167  if (orthogonalization_strategy ==
1169  OrthogonalizationStrategy::modified_gram_schmidt)
1170  {
1171  double htmp = vv * orthogonal_vectors[0];
1172  h(0) += htmp;
1173  for (unsigned int i = 1; i < dim; ++i)
1174  {
1175  htmp = vv.add_and_dot(-htmp,
1176  orthogonal_vectors[i - 1],
1177  orthogonal_vectors[i]);
1178  h(i) += htmp;
1179  }
1180 
1181  norm_vv = std::sqrt(
1182  vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
1183  }
1184  else if (orthogonalization_strategy ==
1186  OrthogonalizationStrategy::classical_gram_schmidt)
1187  {
1188  Tvmult_add(dim, vv, orthogonal_vectors, h);
1189  norm_vv = substract_and_norm(dim, orthogonal_vectors, h, vv);
1190  }
1191  else
1192  {
1193  AssertThrow(false, ExcNotImplemented());
1194  }
1195 
1196  if (c == 1)
1197  return norm_vv; // reorthogonalization already performed -> finished
1198 
1199  // Re-orthogonalization if loss of orthogonality detected. For the
1200  // test, use a strategy discussed in C. T. Kelley, Iterative Methods
1201  // for Linear and Nonlinear Equations, SIAM, Philadelphia, 1995:
1202  // Compare the norm of vv after orthogonalization with its norm when
1203  // starting the orthogonalization. If vv became very small (here: less
1204  // than the square root of the machine precision times 10), it is
1205  // almost in the span of the previous vectors, which indicates loss of
1206  // precision.
1207  if (consider_reorthogonalize)
1208  {
1209  if (norm_vv >
1210  10. * norm_vv_start *
1211  std::sqrt(std::numeric_limits<
1212  typename VectorType::value_type>::epsilon()))
1213  return norm_vv;
1214 
1215  else
1216  {
1217  reorthogonalize = true;
1218  if (!reorthogonalize_signal.empty())
1219  reorthogonalize_signal(accumulated_iterations);
1220  }
1221  }
1222 
1223  if (reorthogonalize == false)
1224  return norm_vv; // no reorthogonalization needed -> finished
1225  }
1226 
1227  AssertThrow(false, ExcInternalError());
1228 
1229  return 0.0;
1230  }
1231  } // namespace SolverGMRESImplementation
1232 } // namespace internal
1233 
1234 
1235 
1236 template <class VectorType>
1237 inline void
1239  const FullMatrix<double> &H_orig,
1240  const unsigned int dim,
1241  const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
1242  &eigenvalues_signal,
1243  const boost::signals2::signal<void(const FullMatrix<double> &)>
1244  & hessenberg_signal,
1245  const boost::signals2::signal<void(double)> &cond_signal)
1246 {
1247  // Avoid copying the Hessenberg matrix if it isn't needed.
1248  if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1249  !cond_signal.empty()) &&
1250  dim > 0)
1251  {
1252  LAPACKFullMatrix<double> mat(dim, dim);
1253  for (unsigned int i = 0; i < dim; ++i)
1254  for (unsigned int j = 0; j < dim; ++j)
1255  mat(i, j) = H_orig(i, j);
1256  hessenberg_signal(H_orig);
1257  // Avoid computing eigenvalues if they are not needed.
1258  if (!eigenvalues_signal.empty())
1259  {
1260  // Copy mat so that we can compute svd below. Necessary since
1261  // compute_eigenvalues will leave mat in state
1262  // LAPACKSupport::unusable.
1263  LAPACKFullMatrix<double> mat_eig(mat);
1264  mat_eig.compute_eigenvalues();
1265  std::vector<std::complex<double>> eigenvalues(dim);
1266  for (unsigned int i = 0; i < mat_eig.n(); ++i)
1267  eigenvalues[i] = mat_eig.eigenvalue(i);
1268  // Sort eigenvalues for nicer output.
1269  std::sort(eigenvalues.begin(),
1270  eigenvalues.end(),
1271  internal::SolverGMRESImplementation::complex_less_pred);
1272  eigenvalues_signal(eigenvalues);
1273  }
1274  // Calculate condition number, avoid calculating the svd if a slot
1275  // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
1276  if (!cond_signal.empty() && (mat.n() > 1))
1277  {
1278  mat.compute_svd();
1279  double condition_number =
1280  mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1281  cond_signal(condition_number);
1282  }
1283  }
1284 }
1285 
1286 
1287 
1288 template <class VectorType>
1289 template <typename MatrixType, typename PreconditionerType>
1290 void
1291 SolverGMRES<VectorType>::solve(const MatrixType & A,
1292  VectorType & x,
1293  const VectorType & b,
1294  const PreconditionerType &preconditioner)
1295 {
1296  // TODO:[?] Check, why there are two different start residuals.
1297  // TODO:[GK] Make sure the parameter in the constructor means maximum basis
1298  // size
1299 
1300  std::unique_ptr<LogStream::Prefix> prefix;
1301  if (!additional_data.batched_mode)
1302  prefix = std::make_unique<LogStream::Prefix>("GMRES");
1303 
1304  // extra call to std::max to placate static analyzers: coverity rightfully
1305  // complains that data.max_n_tmp_vectors - 2 may overflow
1306  const unsigned int n_tmp_vectors =
1307  std::max(additional_data.max_n_tmp_vectors, 3u);
1308 
1309  // Generate an object where basis vectors are stored.
1311  n_tmp_vectors, this->memory);
1312 
1313  // number of the present iteration; this
1314  // number is not reset to zero upon a
1315  // restart
1316  unsigned int accumulated_iterations = 0;
1317 
1318  const bool do_eigenvalues =
1319  !additional_data.batched_mode &&
1320  (!condition_number_signal.empty() ||
1321  !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1322  !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1323  !all_hessenberg_signal.empty());
1324  // for eigenvalue computation, need to collect the Hessenberg matrix (before
1325  // applying Givens rotations)
1326  FullMatrix<double> H_orig;
1327  if (do_eigenvalues)
1328  H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
1329 
1330  // matrix used for the orthogonalization process later
1331  H.reinit(n_tmp_vectors, n_tmp_vectors - 1, /* omit_initialization */ true);
1332 
1333  // some additional vectors, also used in the orthogonalization
1334  gamma.reinit(n_tmp_vectors);
1335  ci.reinit(n_tmp_vectors - 1);
1336  si.reinit(n_tmp_vectors - 1);
1337  h.reinit(n_tmp_vectors - 1);
1338 
1339  unsigned int dim = 0;
1340 
1341  SolverControl::State iteration_state = SolverControl::iterate;
1342  double last_res = std::numeric_limits<double>::lowest();
1343 
1344  // switch to determine whether we want a left or a right preconditioner. at
1345  // present, left is default, but both ways are implemented
1346  const bool left_precondition = !additional_data.right_preconditioning;
1347 
1348  // Per default the left preconditioned GMRes uses the preconditioned
1349  // residual and the right preconditioned GMRes uses the unpreconditioned
1350  // residual as stopping criterion.
1351  const bool use_default_residual = additional_data.use_default_residual;
1352 
1353  // define two aliases
1354  VectorType &v = tmp_vectors(0, x);
1355  VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
1356 
1357  // Following vectors are needed when we are not using the default residuals
1358  // as stopping criterion
1361  std::unique_ptr<::Vector<double>> gamma_;
1362  if (!use_default_residual)
1363  {
1364  r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1365  x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1366  r->reinit(x);
1367  x_->reinit(x);
1368 
1369  gamma_ = std::make_unique<::Vector<double>>(gamma.size());
1370  }
1371 
1372  bool re_orthogonalize = additional_data.force_re_orthogonalization;
1373 
1375  // outer iteration: loop until we either reach convergence or the maximum
1376  // number of iterations is exceeded. each cycle of this loop amounts to one
1377  // restart
1378  do
1379  {
1380  // reset this vector to the right size
1381  h.reinit(n_tmp_vectors - 1);
1382 
1383  double rho = 0.0;
1384 
1385  if (left_precondition)
1386  {
1387  A.vmult(p, x);
1388  p.sadd(-1., 1., b);
1389  preconditioner.vmult(v, p);
1390  rho = v.l2_norm();
1391  }
1392  else
1393  {
1394  A.vmult(v, x);
1395  rho = ::internal::SolverGMRESImplementation::sadd_and_norm(v,
1396  -1,
1397  b,
1398  1.0);
1399  }
1400 
1401  // check the residual here as well since it may be that we got the exact
1402  // (or an almost exact) solution vector at the outset. if we wouldn't
1403  // check here, the next scaling operation would produce garbage
1404  if (use_default_residual)
1405  {
1406  last_res = rho;
1407  if (additional_data.batched_mode)
1408  iteration_state = solver_control.check(accumulated_iterations, rho);
1409  else
1410  iteration_state =
1411  this->iteration_status(accumulated_iterations, rho, x);
1412 
1413  if (iteration_state != SolverControl::iterate)
1414  break;
1415  }
1416  else
1417  {
1418  deallog << "default_res=" << rho << std::endl;
1419 
1420  if (left_precondition)
1421  {
1422  A.vmult(*r, x);
1423  r->sadd(-1., 1., b);
1424  }
1425  else
1426  preconditioner.vmult(*r, v);
1427 
1428  double res = r->l2_norm();
1429  last_res = res;
1430  if (additional_data.batched_mode)
1431  iteration_state = solver_control.check(accumulated_iterations, rho);
1432  else
1433  iteration_state =
1434  this->iteration_status(accumulated_iterations, res, x);
1435 
1436  if (iteration_state != SolverControl::iterate)
1437  break;
1438  }
1439 
1440  gamma(0) = rho;
1441 
1442  v *= 1. / rho;
1443 
1444  // inner iteration doing at most as many steps as there are temporary
1445  // vectors. the number of steps actually been done is propagated outside
1446  // through the @p dim variable
1447  for (unsigned int inner_iteration = 0;
1448  ((inner_iteration < n_tmp_vectors - 2) &&
1449  (iteration_state == SolverControl::iterate));
1450  ++inner_iteration)
1451  {
1452  ++accumulated_iterations;
1453  // yet another alias
1454  VectorType &vv = tmp_vectors(inner_iteration + 1, x);
1455 
1456  if (left_precondition)
1457  {
1458  A.vmult(p, tmp_vectors[inner_iteration]);
1459  preconditioner.vmult(vv, p);
1460  }
1461  else
1462  {
1463  preconditioner.vmult(p, tmp_vectors[inner_iteration]);
1464  A.vmult(vv, p);
1465  }
1466 
1467  dim = inner_iteration + 1;
1468 
1469  const double s =
1470  internal::SolverGMRESImplementation::iterated_gram_schmidt(
1471  additional_data.orthogonalization_strategy,
1472  tmp_vectors,
1473  dim,
1474  accumulated_iterations,
1475  vv,
1476  h,
1477  re_orthogonalize,
1478  re_orthogonalize_signal);
1479  h(inner_iteration + 1) = s;
1480 
1481  // s=0 is a lucky breakdown, the solver will reach convergence,
1482  // but we must not divide by zero here.
1483  if (s != 0)
1484  vv *= 1. / s;
1485 
1486  // for eigenvalues, get the resulting coefficients from the
1487  // orthogonalization process
1488  if (do_eigenvalues)
1489  for (unsigned int i = 0; i < dim + 1; ++i)
1490  H_orig(i, inner_iteration) = h(i);
1491 
1492  // Transformation into tridiagonal structure
1493  givens_rotation(h, gamma, ci, si, inner_iteration);
1494 
1495  // append vector on matrix
1496  for (unsigned int i = 0; i < dim; ++i)
1497  H(i, inner_iteration) = h(i);
1498 
1499  // default residual
1500  rho = std::fabs(gamma(dim));
1501 
1502  if (use_default_residual)
1503  {
1504  last_res = rho;
1505  if (additional_data.batched_mode)
1506  iteration_state =
1507  solver_control.check(accumulated_iterations, rho);
1508  else
1509  iteration_state =
1510  this->iteration_status(accumulated_iterations, rho, x);
1511  }
1512  else
1513  {
1514  if (!additional_data.batched_mode)
1515  deallog << "default_res=" << rho << std::endl;
1516 
1517  *x_ = x;
1518  *gamma_ = gamma;
1519  internal::SolverGMRESImplementation::solve_triangular(dim,
1520  H,
1521  *gamma_,
1522  h);
1523 
1524  if (left_precondition)
1525  for (unsigned int i = 0; i < dim; ++i)
1526  x_->add(h(i), tmp_vectors[i]);
1527  else
1528  {
1529  p = 0.;
1530  for (unsigned int i = 0; i < dim; ++i)
1531  p.add(h(i), tmp_vectors[i]);
1532  preconditioner.vmult(*r, p);
1533  x_->add(1., *r);
1534  };
1535  A.vmult(*r, *x_);
1536  r->sadd(-1., 1., b);
1537  // Now *r contains the unpreconditioned residual!!
1538  if (left_precondition)
1539  {
1540  const double res = r->l2_norm();
1541  last_res = res;
1542 
1543  iteration_state =
1544  this->iteration_status(accumulated_iterations, res, x);
1545  }
1546  else
1547  {
1548  preconditioner.vmult(*x_, *r);
1549  const double preconditioned_res = x_->l2_norm();
1550  last_res = preconditioned_res;
1551 
1552  if (additional_data.batched_mode)
1553  iteration_state =
1554  solver_control.check(accumulated_iterations, rho);
1555  else
1556  iteration_state =
1557  this->iteration_status(accumulated_iterations,
1558  preconditioned_res,
1559  x);
1560  }
1561  }
1562  }
1563 
1564  // end of inner iteration. now calculate the solution from the temporary
1565  // vectors
1566  internal::SolverGMRESImplementation::solve_triangular(dim, H, gamma, h);
1567 
1568  if (do_eigenvalues)
1569  compute_eigs_and_cond(H_orig,
1570  dim,
1571  all_eigenvalues_signal,
1572  all_hessenberg_signal,
1573  condition_number_signal);
1574 
1575  if (left_precondition)
1576  ::internal::SolverGMRESImplementation::add(
1577  x, dim, h, tmp_vectors, false);
1578  else
1579  {
1580  ::internal::SolverGMRESImplementation::add(
1581  p, dim, h, tmp_vectors, true);
1582  preconditioner.vmult(v, p);
1583  x.add(1., v);
1584  };
1585  // end of outer iteration. restart if no convergence and the number of
1586  // iterations is not exceeded
1587  }
1588  while (iteration_state == SolverControl::iterate);
1589 
1590  if (do_eigenvalues)
1591  compute_eigs_and_cond(H_orig,
1592  dim,
1593  eigenvalues_signal,
1594  hessenberg_signal,
1595  condition_number_signal);
1596 
1597  if (!additional_data.batched_mode && !krylov_space_signal.empty())
1598  krylov_space_signal(tmp_vectors);
1599 
1600  // in case of failure: throw exception
1601  AssertThrow(iteration_state == SolverControl::success,
1602  SolverControl::NoConvergence(accumulated_iterations, last_res));
1603 }
1604 
1605 
1606 
1607 template <class VectorType>
1608 boost::signals2::connection
1610  const std::function<void(double)> &slot,
1611  const bool every_iteration)
1612 {
1613  if (every_iteration)
1614  {
1615  return all_condition_numbers_signal.connect(slot);
1616  }
1617  else
1618  {
1619  return condition_number_signal.connect(slot);
1620  }
1621 }
1622 
1623 
1624 
1625 template <class VectorType>
1626 boost::signals2::connection
1628  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1629  const bool every_iteration)
1630 {
1631  if (every_iteration)
1632  {
1633  return all_eigenvalues_signal.connect(slot);
1634  }
1635  else
1636  {
1637  return eigenvalues_signal.connect(slot);
1638  }
1639 }
1640 
1641 
1642 
1643 template <class VectorType>
1644 boost::signals2::connection
1646  const std::function<void(const FullMatrix<double> &)> &slot,
1647  const bool every_iteration)
1648 {
1649  if (every_iteration)
1650  {
1651  return all_hessenberg_signal.connect(slot);
1652  }
1653  else
1654  {
1655  return hessenberg_signal.connect(slot);
1656  }
1657 }
1658 
1659 
1660 
1661 template <class VectorType>
1662 boost::signals2::connection
1664  const std::function<void(
1666 {
1667  return krylov_space_signal.connect(slot);
1668 }
1669 
1670 
1671 
1672 template <class VectorType>
1673 boost::signals2::connection
1675  const std::function<void(int)> &slot)
1676 {
1677  return re_orthogonalize_signal.connect(slot);
1678 }
1679 
1680 
1681 
1682 template <class VectorType>
1683 double
1685 {
1686  // dummy implementation. this function is not needed for the present
1687  // implementation of gmres
1688  Assert(false, ExcInternalError());
1689  return 0;
1690 }
1691 
1692 
1693 //----------------------------------------------------------------------//
1694 
1695 template <class VectorType>
1698  const AdditionalData & data)
1699  : SolverBase<VectorType>(cn, mem)
1700  , additional_data(data)
1701 {}
1702 
1703 
1704 
1705 template <class VectorType>
1707  const AdditionalData &data)
1708  : SolverBase<VectorType>(cn)
1709  , additional_data(data)
1710 {}
1711 
1712 
1713 
1714 template <class VectorType>
1715 template <typename MatrixType, typename PreconditionerType>
1716 void
1717 SolverFGMRES<VectorType>::solve(const MatrixType & A,
1718  VectorType & x,
1719  const VectorType & b,
1720  const PreconditionerType &preconditioner)
1721 {
1722  LogStream::Prefix prefix("FGMRES");
1723 
1724  SolverControl::State iteration_state = SolverControl::iterate;
1725 
1726  const unsigned int basis_size = additional_data.max_basis_size;
1727 
1728  // Generate an object where basis vectors are stored.
1730  basis_size, this->memory);
1732  basis_size, this->memory);
1733 
1734  // number of the present iteration; this number is not reset to zero upon a
1735  // restart
1736  unsigned int accumulated_iterations = 0;
1737 
1738  // matrix used for the orthogonalization process later
1739  H.reinit(basis_size + 1, basis_size);
1740 
1741  // Vectors for projected system
1742  Vector<double> projected_rhs;
1743  Vector<double> y;
1744 
1745  // Iteration starts here
1746  double res = std::numeric_limits<double>::lowest();
1747 
1748  typename VectorMemory<VectorType>::Pointer aux(this->memory);
1749  aux->reinit(x);
1750  do
1751  {
1752  A.vmult(*aux, x);
1753  aux->sadd(-1., 1., b);
1754 
1755  double beta = aux->l2_norm();
1756  res = beta;
1757  iteration_state = this->iteration_status(accumulated_iterations, res, x);
1758  if (iteration_state == SolverControl::success)
1759  break;
1760 
1761  H.reinit(basis_size + 1, basis_size);
1762  double a = beta;
1763 
1764  for (unsigned int j = 0; j < basis_size; ++j)
1765  {
1766  if (a != 0) // treat lucky breakdown
1767  v(j, x).equ(1. / a, *aux);
1768  else
1769  v(j, x) = 0.;
1770 
1771 
1772  preconditioner.vmult(z(j, x), v[j]);
1773  A.vmult(*aux, z[j]);
1774 
1775  // Gram-Schmidt
1776  H(0, j) = *aux * v[0];
1777  for (unsigned int i = 1; i <= j; ++i)
1778  H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1779  H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1780 
1781  // Compute projected solution
1782 
1783  if (j > 0)
1784  {
1785  H1.reinit(j + 1, j);
1786  projected_rhs.reinit(j + 1);
1787  y.reinit(j);
1788  projected_rhs(0) = beta;
1789  H1.fill(H);
1790 
1791  // check convergence. note that the vector 'x' we pass to the
1792  // criterion is not the final solution we compute if we
1793  // decide to jump out of the iteration (we update 'x' again
1794  // right after the current loop)
1795  Householder<double> house(H1);
1796  res = house.least_squares(y, projected_rhs);
1797  iteration_state =
1798  this->iteration_status(++accumulated_iterations, res, x);
1799  if (iteration_state != SolverControl::iterate)
1800  break;
1801  }
1802  }
1803 
1804  // Update solution vector
1805  for (unsigned int j = 0; j < y.size(); ++j)
1806  x.add(y(j), z[j]);
1807  }
1808  while (iteration_state == SolverControl::iterate);
1809 
1810  // in case of failure: throw exception
1811  if (iteration_state != SolverControl::success)
1812  AssertThrow(false,
1813  SolverControl::NoConvergence(accumulated_iterations, res));
1814 }
1815 
1816 #endif // DOXYGEN
1817 
1819 
1820 #endif
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:566
SolverFGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
Definition: solver_gmres.h:561
FullMatrix< double > H1
Definition: solver_gmres.h:571
boost::signals2::signal< void(const FullMatrix< double > &)> hessenberg_signal
Definition: solver_gmres.h:400
Vector< double > si
Definition: solver_gmres.h:481
AdditionalData additional_data
Definition: solver_gmres.h:368
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_gmres.h:374
SolverGMRES(const SolverGMRES< VectorType > &)=delete
boost::signals2::signal< void(const std::vector< std::complex< double >> &)> eigenvalues_signal
Definition: solver_gmres.h:387
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:415
boost::signals2::signal< void(const FullMatrix< double > &)> all_hessenberg_signal
Definition: solver_gmres.h:407
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:476
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
SolverControl & solver_control
Definition: solver_gmres.h:428
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:471
boost::signals2::signal< void(const std::vector< std::complex< double >> &)> all_eigenvalues_signal
Definition: solver_gmres.h:394
Vector< double > h
Definition: solver_gmres.h:486
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
FullMatrix< double > H
Definition: solver_gmres.h:466
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:380
boost::signals2::signal< void(int)> re_orthogonalize_signal
Definition: solver_gmres.h:421
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:123
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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1854
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcTooFewTmpVectors(int arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:510
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
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
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
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)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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:523
OrthogonalizationStrategy orthogonalization_strategy
Definition: solver_gmres.h:273
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)