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