deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20: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\}}\)
Loading...
Searching...
No Matches
solver_gmres.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_solver_gmres_h
16#define dealii_solver_gmres_h
17
18
19
20#include <deal.II/base/config.h>
21
25
30#include <deal.II/lac/solver.h>
32#include <deal.II/lac/vector.h>
33
34#include <boost/signals2/signal.hpp>
35
36#include <algorithm>
37#include <cmath>
38#include <limits>
39#include <memory>
40#include <vector>
41
43
44// forward declarations
45#ifndef DOXYGEN
46namespace LinearAlgebra
47{
48 namespace distributed
49 {
50 template <typename, typename>
51 class Vector;
52 } // namespace distributed
53} // namespace LinearAlgebra
54#endif
55
61namespace internal
62{
66 namespace SolverGMRESImplementation
67 {
75 template <typename VectorType>
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
126
127
136 template <typename Number>
138 {
139 public:
144 void
147 const unsigned int max_basis_size,
148 const bool force_reorthogonalization);
149
172 template <typename VectorType>
173 double
175 const unsigned int n,
177 const unsigned int accumulated_iterations = 0,
178 const boost::signals2::signal<void(int)> &reorthogonalize_signal =
179 boost::signals2::signal<void(int)>());
180
188 const Vector<double> &
190
195 const FullMatrix<double> &
197
201 std::vector<const Number *> vector_ptrs;
202
203 private:
208
215
220 std::vector<std::pair<double, double>> givens_rotations;
221
226
232
237
247
252
280 double
281 do_givens_rotation(const bool delayed_reorthogonalization,
282 const int col,
283 FullMatrix<double> &matrix,
284 std::vector<std::pair<double, double>> &rotations,
286 };
287 } // namespace SolverGMRESImplementation
288} // namespace internal
289
357template <typename VectorType = Vector<double>>
360{
361public:
366 {
377 explicit AdditionalData(const unsigned int max_basis_size = 30,
378 const bool right_preconditioning = false,
379 const bool use_default_residual = true,
380 const bool force_re_orthogonalization = false,
381 const bool batched_mode = false,
383 orthogonalization_strategy =
385 delayed_classical_gram_schmidt);
386
398 unsigned int max_n_tmp_vectors;
399
407 unsigned int max_basis_size;
408
417
422
430
438
443 };
444
451
457
462
466 template <typename MatrixType, typename PreconditionerType>
470 void solve(const MatrixType &A,
471 VectorType &x,
472 const VectorType &b,
473 const PreconditionerType &preconditioner);
474
481 boost::signals2::connection
482 connect_condition_number_slot(const std::function<void(double)> &slot,
484
491 boost::signals2::connection
492 connect_eigenvalues_slot(
493 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
495
503 boost::signals2::connection
504 connect_hessenberg_slot(
505 const std::function<void(const FullMatrix<double> &)> &slot,
506 const bool every_iteration = true);
507
514 boost::signals2::connection
515 connect_krylov_space_slot(
516 const std::function<
517 void(const internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
518 &slot);
519
520
525 boost::signals2::connection
526 connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
527
528
529 DeclException1(ExcTooFewTmpVectors,
530 int,
531 << "The number of temporary vectors you gave (" << arg1
532 << ") is too small. It should be at least 10 for "
533 << "any results, and much more for reasonable ones.");
534
539 AdditionalData additional_data;
540
545 boost::signals2::signal<void(double)> condition_number_signal;
546
551 boost::signals2::signal<void(double)> all_condition_numbers_signal;
552
557 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
558 eigenvalues_signal;
559
564 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
565 all_eigenvalues_signal;
566
571 boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
572
577 boost::signals2::signal<void(const FullMatrix<double> &)>
578 all_hessenberg_signal;
579
584 boost::signals2::signal<void(
585 const internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
586 krylov_space_signal;
587
592 boost::signals2::signal<void(int)> re_orthogonalize_signal;
593
599 SolverControl &solver_control;
600
604 virtual double
605 criterion();
606
613 static void
614 compute_eigs_and_cond(
615 const FullMatrix<double> &H_orig,
616 const unsigned int n,
618 void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
619 const boost::signals2::signal<void(const FullMatrix<double> &)>
620 &hessenberg_signal,
621 const boost::signals2::signal<void(double)> &cond_signal);
622
627 internal::SolverGMRESImplementation::ArnoldiProcess<
628 typename VectorType::value_type>
629 arnoldi_process;
630};
631
632
633
654template <typename VectorType = Vector<double>>
655DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
657{
658public:
663 {
667 explicit AdditionalData(const unsigned int max_basis_size = 30,
669 orthogonalization_strategy =
671 delayed_classical_gram_schmidt)
672 : max_basis_size(max_basis_size)
673 , orthogonalization_strategy(orthogonalization_strategy)
674 {}
675
679 unsigned int max_basis_size;
680
685 };
686
693
700
704 template <typename MatrixType, typename PreconditionerType>
708 void solve(const MatrixType &A,
709 VectorType &x,
710 const VectorType &b,
711 const PreconditionerType &preconditioner);
712
713private:
717 AdditionalData additional_data;
718
723 internal::SolverGMRESImplementation::ArnoldiProcess<
724 typename VectorType::value_type>
725 arnoldi_process;
726};
727
729/* --------------------- Inline and template functions ------------------- */
730
731
732#ifndef DOXYGEN
733
734template <typename VectorType>
737 const unsigned int max_basis_size,
738 const bool right_preconditioning,
739 const bool use_default_residual,
740 const bool force_re_orthogonalization,
741 const bool batched_mode,
742 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy)
743 : max_n_tmp_vectors(0)
744 , max_basis_size(max_basis_size)
745 , right_preconditioning(right_preconditioning)
746 , use_default_residual(use_default_residual)
747 , force_re_orthogonalization(force_re_orthogonalization)
748 , batched_mode(batched_mode)
749 , orthogonalization_strategy(orthogonalization_strategy)
750{
751 Assert(max_basis_size >= 1,
752 ExcMessage("SolverGMRES needs at least one vector in the "
753 "Arnoldi basis."));
754}
755
756
757
758template <typename VectorType>
762 const AdditionalData &data)
763 : SolverBase<VectorType>(cn, mem)
764 , additional_data(data)
765 , solver_control(cn)
766{}
767
768
769
770template <typename VectorType>
773 const AdditionalData &data)
775 , additional_data(data)
776 , solver_control(cn)
777{}
778
779
780
781namespace internal
782{
783 namespace SolverGMRESImplementation
784 {
785 template <typename VectorType>
786 inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
788 : mem(vmem)
789 , data(max_size)
790 {}
791
792
793
794 template <typename VectorType>
795 inline VectorType &
796 TmpVectors<VectorType>::operator[](const unsigned int i) const
797 {
798 AssertIndexRange(i, data.size());
799
800 Assert(data[i] != nullptr, ExcNotInitialized());
801 return *data[i];
802 }
803
804
805
806 template <typename VectorType>
807 inline VectorType &
808 TmpVectors<VectorType>::operator()(const unsigned int i,
809 const VectorType &temp)
810 {
811 AssertIndexRange(i, data.size());
812 if (data[i] == nullptr)
813 {
814 data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
815 data[i]->reinit(temp, true);
816 }
817 return *data[i];
818 }
819
820
821
822 template <typename VectorType>
823 unsigned int
824 TmpVectors<VectorType>::size() const
825 {
826 return (data.size() > 0 ? data.size() - 1 : 0);
827 }
828
829
830
831 template <typename VectorType, typename Enable = void>
833
834 template <typename VectorType>
837 std::enable_if_t<!internal::is_block_vector<VectorType>>>
838 {
839 static constexpr bool value =
840 std::is_same_v<
842 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
844 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
845 };
846
847
848
849 template <typename VectorType>
852 std::enable_if_t<internal::is_block_vector<VectorType>>>
853 {
854 static constexpr bool value =
855 std::is_same_v<
856 typename VectorType::BlockType,
857 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
859 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
860 };
861
862
863
864 template <typename VectorType,
865 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
866 * = nullptr>
867 unsigned int
868 n_blocks(const VectorType &)
869 {
870 return 1;
871 }
872
873
874
875 template <typename VectorType,
876 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
877 nullptr>
878 unsigned int
879 n_blocks(const VectorType &vector)
880 {
881 return vector.n_blocks();
882 }
883
884
885
886 template <typename VectorType,
887 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
888 * = nullptr>
889 VectorType &
890 block(VectorType &vector, const unsigned int b)
891 {
892 AssertDimension(b, 0);
893 (void)b;
894 return vector;
895 }
896
897
898
899 template <typename VectorType,
900 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
901 * = nullptr>
902 const VectorType &
903 block(const VectorType &vector, const unsigned int b)
904 {
905 AssertDimension(b, 0);
906 (void)b;
907 return vector;
908 }
909
910
911
912 template <typename VectorType,
913 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
914 nullptr>
915 typename VectorType::BlockType &
916 block(VectorType &vector, const unsigned int b)
917 {
918 return vector.block(b);
919 }
920
921
922
923 template <typename VectorType,
924 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
925 nullptr>
926 const typename VectorType::BlockType &
927 block(const VectorType &vector, const unsigned int b)
928 {
929 return vector.block(b);
930 }
931
932
933
934 template <bool delayed_reorthogonalization,
935 typename VectorType,
936 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
937 VectorType> * = nullptr>
938 void
939 Tvmult_add(const unsigned int n,
940 const VectorType &vv,
943 std::vector<const typename VectorType::value_type *> &)
944 {
945 for (unsigned int i = 0; i < n; ++i)
946 {
947 h(i) += vv * orthogonal_vectors[i];
948 if (delayed_reorthogonalization)
949 h(n + i) += orthogonal_vectors[i] * orthogonal_vectors[n - 1];
950 }
951 if (delayed_reorthogonalization)
952 h(n + n) += vv * vv;
953 }
954
955
956
957 // worker method for deal.II's vector types implemented in .cc file
958 template <bool delayed_reorthogonalization, typename Number>
959 void
960 do_Tvmult_add(const unsigned int n_vectors,
961 const std::size_t locally_owned_size,
962 const Number *current_vector,
963 const std::vector<const Number *> &orthogonal_vectors,
964 Vector<double> &b);
965
966
967
968 template <bool delayed_reorthogonalization,
969 typename VectorType,
970 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
971 VectorType> * = nullptr>
972 void
973 Tvmult_add(
974 const unsigned int n,
975 const VectorType &vv,
978 std::vector<const typename VectorType::value_type *> &vector_ptrs)
979 {
980 for (unsigned int b = 0; b < n_blocks(vv); ++b)
981 {
982 vector_ptrs.resize(n);
983 for (unsigned int i = 0; i < n; ++i)
984 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
985
987 block(vv, b).end() -
988 block(vv, b).begin(),
989 block(vv, b).begin(),
990 vector_ptrs,
991 h);
992 }
993
994 Utilities::MPI::sum(h, block(vv, 0).get_mpi_communicator(), h);
995 }
996
997
998
999 template <bool delayed_reorthogonalization,
1000 typename VectorType,
1001 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1002 VectorType> * = nullptr>
1003 double
1004 subtract_and_norm(const unsigned int n,
1006 const Vector<double> &h,
1007 VectorType &vv,
1008 std::vector<const typename VectorType::value_type *> &)
1009 {
1010 Assert(n > 0, ExcInternalError());
1011
1013 const_cast<VectorType &>(orthogonal_vectors[n - 1]);
1014 for (unsigned int i = 0; i < n - 1; ++i)
1015 {
1016 if (delayed_reorthogonalization && i + 2 < n)
1017 last_vector.add(-h(n + i), orthogonal_vectors[i]);
1018 vv.add(-h(i), orthogonal_vectors[i]);
1019 }
1020
1021 if (delayed_reorthogonalization)
1022 {
1023 if (n > 1)
1024 last_vector.sadd(1. / h(n + n - 1),
1025 -h(n + n - 2) / h(n + n - 1),
1026 orthogonal_vectors[n - 2]);
1027
1028 // h(n + n) is lucky breakdown
1029 const double scaling_factor_vv = h(n + n) > 0.0 ?
1030 1. / (h(n + n - 1) * h(n + n)) :
1031 1. / (h(n + n - 1) * h(n + n - 1));
1032 vv.sadd(scaling_factor_vv,
1033 -h(n - 1) * scaling_factor_vv,
1034 last_vector);
1035
1036 // the delayed reorthogonalization computes the norm from other
1037 // quantities
1038 return std::numeric_limits<double>::signaling_NaN();
1039 }
1040 else
1041 return std::sqrt(
1042 vv.add_and_dot(-h(n - 1), orthogonal_vectors[n - 1], vv));
1043 }
1044
1045
1046
1047 // worker method for deal.II's vector types implemented in .cc file
1048 template <bool delayed_reorthogonalization, typename Number>
1049 double
1050 do_subtract_and_norm(const unsigned int n_vectors,
1051 const std::size_t locally_owned_size,
1052 const std::vector<const Number *> &orthogonal_vectors,
1053 const Vector<double> &h,
1054 Number *current_vector);
1055
1056
1057
1058 template <bool delayed_reorthogonalization,
1059 typename VectorType,
1060 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1061 VectorType> * = nullptr>
1062 double
1064 const unsigned int n,
1066 const Vector<double> &h,
1067 VectorType &vv,
1068 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1069 {
1070 double norm_vv_temp = 0.0;
1071
1072 for (unsigned int b = 0; b < n_blocks(vv); ++b)
1073 {
1074 vector_ptrs.resize(n);
1075 for (unsigned int i = 0; i < n; ++i)
1076 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
1077
1079 n,
1080 block(vv, b).end() - block(vv, b).begin(),
1081 vector_ptrs,
1082 h,
1083 block(vv, b).begin());
1084 }
1085
1086 return std::sqrt(
1087 Utilities::MPI::sum(norm_vv_temp, block(vv, 0).get_mpi_communicator()));
1088 }
1089
1090
1091
1092 template <typename VectorType,
1093 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1094 VectorType> * = nullptr>
1095 void
1096 add(VectorType &p,
1097 const unsigned int n,
1098 const Vector<double> &h,
1100 const bool zero_out,
1101 std::vector<const typename VectorType::value_type *> &)
1102 {
1103 if (zero_out)
1104 p.equ(h(0), tmp_vectors[0]);
1105 else
1106 p.add(h(0), tmp_vectors[0]);
1107
1108 for (unsigned int i = 1; i < n; ++i)
1109 p.add(h(i), tmp_vectors[i]);
1110 }
1111
1112
1113
1114 // worker method for deal.II's vector types implemented in .cc file
1115 template <typename Number>
1116 void
1117 do_add(const unsigned int n_vectors,
1118 const std::size_t locally_owned_size,
1119 const std::vector<const Number *> &tmp_vectors,
1120 const Vector<double> &h,
1121 const bool zero_out,
1122 Number *output);
1123
1124
1125
1126 template <typename VectorType,
1127 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1128 VectorType> * = nullptr>
1129 void
1130 add(VectorType &p,
1131 const unsigned int n,
1132 const Vector<double> &h,
1134 const bool zero_out,
1135 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1136 {
1137 for (unsigned int b = 0; b < n_blocks(p); ++b)
1138 {
1139 vector_ptrs.resize(n);
1140 for (unsigned int i = 0; i < n; ++i)
1141 vector_ptrs[i] = block(tmp_vectors[i], b).begin();
1142 do_add(n,
1143 block(p, b).end() - block(p, b).begin(),
1144 vector_ptrs,
1145 h,
1146 zero_out,
1147 block(p, b).begin());
1148 }
1149 }
1150
1151
1152
1153 template <typename Number>
1154 inline void
1155 ArnoldiProcess<Number>::initialize(
1156 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy,
1157 const unsigned int basis_size,
1158 const bool force_reorthogonalization)
1159 {
1160 this->orthogonalization_strategy = orthogonalization_strategy;
1161 this->do_reorthogonalization = force_reorthogonalization;
1162
1163 hessenberg_matrix.reinit(basis_size + 1, basis_size);
1164 triangular_matrix.reinit(basis_size + 1, basis_size, true);
1165
1166 // some additional vectors, also used in the orthogonalization
1167 projected_rhs.reinit(basis_size + 1, true);
1168 givens_rotations.reserve(basis_size);
1169
1170 if (orthogonalization_strategy ==
1173 h.reinit(2 * basis_size + 3);
1174 else
1175 h.reinit(basis_size + 1);
1176 }
1177
1178
1179
1180 template <typename Number>
1181 template <typename VectorType>
1182 inline double
1183 ArnoldiProcess<Number>::orthonormalize_nth_vector(
1184 const unsigned int n,
1186 const unsigned int accumulated_iterations,
1187 const boost::signals2::signal<void(int)> &reorthogonalize_signal)
1188 {
1189 AssertIndexRange(n, hessenberg_matrix.m());
1190 AssertIndexRange(n, orthogonal_vectors.size() + 1);
1191
1193
1194 double residual_estimate = std::numeric_limits<double>::signaling_NaN();
1195 if (n == 0)
1196 {
1197 givens_rotations.clear();
1198 residual_estimate = vv.l2_norm();
1199 if (residual_estimate != 0.)
1201 projected_rhs(0) = residual_estimate;
1202 }
1203 else if (orthogonalization_strategy ==
1206 {
1207 // The algorithm implemented in the following few lines is algorithm
1208 // 4 of Bielich et al. (2022).
1209
1210 // To avoid un-scaled numbers as appearing with the original
1211 // algorithm by Bielich et al., we use a preliminary scaling of the
1212 // last vector. This will be corrected in the delayed step.
1213 const double previous_scaling = n > 0 ? h(n + n - 2) : 1.;
1214
1215 // Reset h to zero
1216 h.reinit(n + n + 1);
1217
1218 // global reduction
1219 Tvmult_add<true>(n, vv, orthogonal_vectors, h, vector_ptrs);
1220
1221 // delayed correction terms
1222 double tmp = 0;
1223 for (unsigned int i = 0; i < n - 1; ++i)
1224 tmp += h(n + i) * h(n + i);
1225 const double alpha_j = h(n + n - 1) > tmp ?
1226 std::sqrt(h(n + n - 1) - tmp) :
1227 std::sqrt(h(n + n - 1));
1228 h(n + n - 1) = alpha_j;
1229
1230 tmp = 0;
1231 for (unsigned int i = 0; i < n - 1; ++i)
1232 tmp += h(i) * h(n + i);
1233 h(n - 1) = (h(n - 1) - tmp) / alpha_j;
1234
1235 // representation of H(j-1)
1236 if (n > 1)
1237 {
1238 for (unsigned int i = 0; i < n - 1; ++i)
1239 hessenberg_matrix(i, n - 2) += h(n + i) * previous_scaling;
1240 hessenberg_matrix(n - 1, n - 2) = alpha_j * previous_scaling;
1241 }
1242 for (unsigned int i = 0; i < n; ++i)
1243 {
1244 double sum = 0;
1245 for (unsigned int j = (i == 0 ? 0 : i - 1); j < n - 1; ++j)
1246 sum += hessenberg_matrix(i, j) * h(n + j);
1247 hessenberg_matrix(i, n - 1) = (h(i) - sum) / alpha_j;
1248 }
1249
1250 // compute norm estimate for approximate convergence criterion
1251 // (value of norm to be corrected in next iteration)
1252 double sum = 0;
1253 for (unsigned int i = 0; i < n - 1; ++i)
1254 sum += h(i) * h(i);
1255 sum += (2. - 1.) * h(n - 1) * h(n - 1);
1256 hessenberg_matrix(n, n - 1) =
1257 std::sqrt(std::abs(h(n + n) - sum)) / alpha_j;
1258
1259 // projection and delayed reorthogonalization. We scale the vector
1260 // vv here by the preliminary norm to avoid working with too large
1261 // values and correct the actual norm in the Hessenberg matrix in
1262 // high precision in the next iteration.
1263 h(n + n) = hessenberg_matrix(n, n - 1);
1264 subtract_and_norm<true>(n, orthogonal_vectors, h, vv, vector_ptrs);
1265
1266 // transform new column of upper Hessenberg matrix into upper
1267 // triangular form by computing the respective factor
1268 residual_estimate = do_givens_rotation(
1269 true, n - 2, triangular_matrix, givens_rotations, projected_rhs);
1270 }
1271 else
1272 {
1273 // need initial norm for detection of re-orthogonalization, see below
1274 double norm_vv = 0.0;
1275 double norm_vv_start = 0;
1276 const bool consider_reorthogonalize =
1277 (do_reorthogonalization == false) && (n % 5 == 0);
1279 norm_vv_start = vv.l2_norm();
1280
1281 // Reset h to zero
1282 h.reinit(n);
1283
1284 // run two loops with index 0: orthogonalize, 1: reorthogonalize
1285 for (unsigned int c = 0; c < 2; ++c)
1286 {
1287 // Orthogonalization
1288 if (orthogonalization_strategy ==
1291 {
1292 double htmp = vv * orthogonal_vectors[0];
1293 h(0) += htmp;
1294 for (unsigned int i = 1; i < n; ++i)
1295 {
1296 htmp = vv.add_and_dot(-htmp,
1297 orthogonal_vectors[i - 1],
1299 h(i) += htmp;
1300 }
1301
1303 vv.add_and_dot(-htmp, orthogonal_vectors[n - 1], vv));
1304 }
1305 else if (orthogonalization_strategy ==
1308 {
1309 Tvmult_add<false>(n, vv, orthogonal_vectors, h, vector_ptrs);
1311 n, orthogonal_vectors, h, vv, vector_ptrs);
1312 }
1313 else
1314 {
1316 }
1317
1318 if (c == 1)
1319 break; // reorthogonalization already performed -> finished
1320
1321 // Re-orthogonalization if loss of orthogonality detected. For the
1322 // test, use a strategy discussed in C. T. Kelley, Iterative
1323 // Methods for Linear and Nonlinear Equations, SIAM, Philadelphia,
1324 // 1995: Compare the norm of vv after orthogonalization with its
1325 // norm when starting the orthogonalization. If vv became very
1326 // small (here: less than the square root of the machine precision
1327 // times 10), it is almost in the span of the previous vectors,
1328 // which indicates loss of precision.
1330 {
1331 if (norm_vv >
1332 10. * norm_vv_start *
1333 std::sqrt(std::numeric_limits<
1334 typename VectorType::value_type>::epsilon()))
1335 break;
1336
1337 else
1338 {
1339 do_reorthogonalization = true;
1340 if (!reorthogonalize_signal.empty())
1342 }
1343 }
1344
1345 if (do_reorthogonalization == false)
1346 break; // no reorthogonalization needed -> finished
1347 }
1348
1349 for (unsigned int i = 0; i < n; ++i)
1350 hessenberg_matrix(i, n - 1) = h(i);
1351 hessenberg_matrix(n, n - 1) = norm_vv;
1352
1353 // norm_vv is a lucky breakdown, the solver will reach convergence,
1354 // but we must not divide by zero here.
1355 if (norm_vv != 0)
1356 vv /= norm_vv;
1357
1358 residual_estimate = do_givens_rotation(
1359 false, n - 1, triangular_matrix, givens_rotations, projected_rhs);
1360 }
1361
1362 return residual_estimate;
1363 }
1364
1365
1366
1367 template <typename Number>
1368 inline double
1369 ArnoldiProcess<Number>::do_givens_rotation(
1370 const bool delayed_reorthogonalization,
1371 const int col,
1372 FullMatrix<double> &matrix,
1373 std::vector<std::pair<double, double>> &rotations,
1375 {
1376 // for the delayed orthogonalization, we can only compute the column of
1377 // the previous iteration (as there will be correction terms added to the
1378 // present column for stability reasons), but we still want to compute
1379 // the residual estimate from the accumulated work; we therefore perform
1380 // givens rotations on two columns simultaneously
1381 if (delayed_reorthogonalization)
1382 {
1383 if (col >= 0)
1384 {
1385 AssertDimension(rotations.size(), static_cast<std::size_t>(col));
1386 matrix(0, col) = hessenberg_matrix(0, col);
1387 }
1388 double H_next = hessenberg_matrix(0, col + 1);
1389 for (int i = 0; i < col; ++i)
1390 {
1391 const double c = rotations[i].first;
1392 const double s = rotations[i].second;
1393 const double Hi = matrix(i, col);
1394 const double Hi1 = hessenberg_matrix(i + 1, col);
1395 H_next = -s * H_next + c * hessenberg_matrix(i + 1, col + 1);
1396 matrix(i, col) = c * Hi + s * Hi1;
1397 matrix(i + 1, col) = -s * Hi + c * Hi1;
1398 }
1399
1400 if (col >= 0)
1401 {
1402 const double H_col1 = hessenberg_matrix(col + 1, col);
1403 const double H_col = matrix(col, col);
1404 const double r = 1. / std::sqrt(H_col * H_col + H_col1 * H_col1);
1405 rotations.emplace_back(H_col * r, H_col1 * r);
1406 matrix(col, col) =
1407 rotations[col].first * H_col + rotations[col].second * H_col1;
1408
1409 rhs(col + 1) = -rotations[col].second * rhs(col);
1410 rhs(col) *= rotations[col].first;
1411
1412 H_next =
1413 -rotations[col].second * H_next +
1414 rotations[col].first * hessenberg_matrix(col + 1, col + 1);
1415 }
1416
1417 const double H_last = hessenberg_matrix(col + 2, col + 1);
1418 const double r = 1. / std::sqrt(H_next * H_next + H_last * H_last);
1419 return std::abs(H_last * r * rhs(col + 1));
1420 }
1421 else
1422 {
1423 AssertDimension(rotations.size(), static_cast<std::size_t>(col));
1424
1425 matrix(0, col) = hessenberg_matrix(0, col);
1426 for (int i = 0; i < col; ++i)
1427 {
1428 const double c = rotations[i].first;
1429 const double s = rotations[i].second;
1430 const double Hi = matrix(i, col);
1431 const double Hi1 = hessenberg_matrix(i + 1, col);
1432 matrix(i, col) = c * Hi + s * Hi1;
1433 matrix(i + 1, col) = -s * Hi + c * Hi1;
1434 }
1435
1436 const double Hi = matrix(col, col);
1437 const double Hi1 = hessenberg_matrix(col + 1, col);
1438 const double r = 1. / std::sqrt(Hi * Hi + Hi1 * Hi1);
1439 rotations.emplace_back(Hi * r, Hi1 * r);
1440 matrix(col, col) =
1441 rotations[col].first * Hi + rotations[col].second * Hi1;
1442
1443 rhs(col + 1) = -rotations[col].second * rhs(col);
1444 rhs(col) *= rotations[col].first;
1445
1446 return std::abs(rhs(col + 1));
1447 }
1448 }
1449
1450
1451
1452 template <typename Number>
1453 inline const Vector<double> &
1454 ArnoldiProcess<Number>::solve_projected_system(
1455 const bool orthogonalization_finished)
1456 {
1459 FullMatrix<double> *matrix = &triangular_matrix;
1460 Vector<double> *rhs = &projected_rhs;
1461 unsigned int n = givens_rotations.size();
1462
1463 // If we solve with the delayed orthogonalization, we still need to
1464 // perform the elimination of the last column before we can solve the
1465 // projected system. We distinguish two cases, one where the
1466 // orthogonalization has finished (i.e., end of inner iteration in
1467 // GMRES) and we can safely overwrite the content of the tridiagonal
1468 // matrix and right hand side, and the case during the inner iterations,
1469 // where we need to create copies of the matrices in the QR
1470 // decomposition as well as the right hand side.
1471 if (orthogonalization_strategy ==
1474 {
1475 n += 1;
1477 {
1478 tmp_triangular_matrix = triangular_matrix;
1479 tmp_rhs = projected_rhs;
1480 std::vector<std::pair<double, double>> tmp_givens_rotations(
1481 givens_rotations);
1482 do_givens_rotation(false,
1483 givens_rotations.size(),
1486 tmp_rhs);
1488 rhs = &tmp_rhs;
1489 }
1490 else
1491 do_givens_rotation(false,
1492 givens_rotations.size(),
1493 triangular_matrix,
1494 givens_rotations,
1495 projected_rhs);
1496 }
1497
1498 // Now solve the triangular system by backward substitution
1499 projected_solution.reinit(n);
1500 for (int i = n - 1; i >= 0; --i)
1501 {
1502 double s = (*rhs)(i);
1503 for (unsigned int j = i + 1; j < n; ++j)
1504 s -= projected_solution(j) * (*matrix)(i, j);
1505 projected_solution(i) = s / (*matrix)(i, i);
1506 AssertIsFinite(projected_solution(i));
1507 }
1508
1509 return projected_solution;
1510 }
1511
1512
1513
1514 template <typename Number>
1515 inline const FullMatrix<double> &
1516 ArnoldiProcess<Number>::get_hessenberg_matrix() const
1517 {
1518 return hessenberg_matrix;
1519 }
1520
1521
1522
1523 // A comparator for better printing eigenvalues
1524 inline bool
1525 complex_less_pred(const std::complex<double> &x,
1526 const std::complex<double> &y)
1527 {
1528 return x.real() < y.real() ||
1529 (x.real() == y.real() && x.imag() < y.imag());
1530 }
1531 } // namespace SolverGMRESImplementation
1532} // namespace internal
1533
1534
1535
1536template <typename VectorType>
1540 const unsigned int n,
1541 const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
1542 &eigenvalues_signal,
1543 const boost::signals2::signal<void(const FullMatrix<double> &)>
1544 &hessenberg_signal,
1545 const boost::signals2::signal<void(double)> &cond_signal)
1546{
1547 // Avoid copying the Hessenberg matrix if it isn't needed.
1548 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1549 !cond_signal.empty()) &&
1550 n > 0)
1551 {
1553 for (unsigned int i = 0; i < n; ++i)
1554 for (unsigned int j = 0; j < n; ++j)
1555 mat(i, j) = H_orig(i, j);
1556 hessenberg_signal(H_orig);
1557 // Avoid computing eigenvalues if they are not needed.
1558 if (!eigenvalues_signal.empty())
1559 {
1560 // Copy mat so that we can compute svd below. Necessary since
1561 // compute_eigenvalues will leave mat in state
1562 // LAPACKSupport::unusable.
1564 mat_eig.compute_eigenvalues();
1565 std::vector<std::complex<double>> eigenvalues(n);
1566 for (unsigned int i = 0; i < mat_eig.n(); ++i)
1567 eigenvalues[i] = mat_eig.eigenvalue(i);
1568 // Sort eigenvalues for nicer output.
1569 std::sort(eigenvalues.begin(),
1570 eigenvalues.end(),
1571 internal::SolverGMRESImplementation::complex_less_pred);
1572 eigenvalues_signal(eigenvalues);
1573 }
1574 // Calculate condition number, avoid calculating the svd if a slot
1575 // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
1576 if (!cond_signal.empty() && (mat.n() > 1))
1577 {
1578 mat.compute_svd();
1579 double condition_number =
1580 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1582 }
1583 }
1584}
1585
1586
1587
1588template <typename VectorType>
1590template <typename MatrixType, typename PreconditionerType>
1594void SolverGMRES<VectorType>::solve(const MatrixType &A,
1595 VectorType &x,
1596 const VectorType &b,
1597 const PreconditionerType &preconditioner)
1598{
1599 std::unique_ptr<LogStream::Prefix> prefix;
1600 if (!additional_data.batched_mode)
1601 prefix = std::make_unique<LogStream::Prefix>("GMRES");
1602
1603 // extra call to std::max to placate static analyzers: coverity rightfully
1604 // complains that data.max_n_tmp_vectors - 2 may overflow
1605 const unsigned int basis_size =
1606 (additional_data.max_basis_size > 0 ?
1607 additional_data.max_basis_size :
1608 std::max(additional_data.max_n_tmp_vectors, 3u) - 2);
1609
1610 // Generate an object where basis vectors are stored.
1612 basis_size + 2, this->memory);
1613
1614 // number of the present iteration; this number is not reset to zero upon a
1615 // restart
1616 unsigned int accumulated_iterations = 0;
1617
1618 const bool do_eigenvalues =
1619 !additional_data.batched_mode &&
1620 (!condition_number_signal.empty() ||
1621 !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1622 !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1623 !all_hessenberg_signal.empty());
1624
1626 double res = std::numeric_limits<double>::lowest();
1627
1628 // switch to determine whether we want a left or a right preconditioner. at
1629 // present, left is default, but both ways are implemented
1630 const bool left_precondition = !additional_data.right_preconditioning;
1631
1632 // Per default the left preconditioned GMRES uses the preconditioned
1633 // residual and the right preconditioned GMRES uses the unpreconditioned
1634 // residual as stopping criterion.
1635 const bool use_default_residual = additional_data.use_default_residual;
1636
1637 // define an alias
1639
1640 // Following vectors are needed when we are not using the default residuals
1641 // as stopping criterion
1644 if (!use_default_residual)
1645 {
1646 r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1647 x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1648 r->reinit(x);
1649 x_->reinit(x);
1650 }
1651
1652 arnoldi_process.initialize(additional_data.orthogonalization_strategy,
1653 basis_size,
1654 additional_data.force_re_orthogonalization);
1655
1657 // outer iteration: loop until we either reach convergence or the maximum
1658 // number of iterations is exceeded. each cycle of this loop amounts to one
1659 // restart
1660 do
1661 {
1662 VectorType &v = basis_vectors(0, x);
1663
1664 // Compute the preconditioned/unpreconditioned residual for left/right
1665 // preconditioning. If 'x' is the zero vector, then we can bypass the
1666 // full computation. But 'x' is only likely to be the zero vector if
1667 // that's what the user provided as the starting guess, so it's only
1668 // worth checking for this in the first iteration. (Calling all_zero()
1669 // costs as much in memory transfer and communication as computing the
1670 // norm of a vector.)
1672 {
1673 if (accumulated_iterations == 0 && x.all_zero())
1674 preconditioner.vmult(v, b);
1675 else
1676 {
1677 A.vmult(p, x);
1678 p.sadd(-1., 1., b);
1679 preconditioner.vmult(v, p);
1680 }
1681 }
1682 else
1683 {
1684 if (accumulated_iterations == 0 && x.all_zero())
1685 v = b;
1686 else
1687 {
1688 A.vmult(v, x);
1689 v.sadd(-1., 1., b);
1690 }
1691 }
1692
1693 const double norm_v = arnoldi_process.orthonormalize_nth_vector(
1694 0, basis_vectors, accumulated_iterations, re_orthogonalize_signal);
1695
1696 // check the residual here as well since it may be that we got the exact
1697 // (or an almost exact) solution vector at the outset. if we wouldn't
1698 // check here, the next scaling operation would produce garbage
1699 if (use_default_residual)
1700 {
1701 res = norm_v;
1702 if (additional_data.batched_mode)
1704 else
1706 this->iteration_status(accumulated_iterations, res, x);
1707
1709 break;
1710 }
1711 else
1712 {
1713 deallog << "default_res=" << norm_v << std::endl;
1714
1716 {
1717 A.vmult(*r, x);
1718 r->sadd(-1., 1., b);
1719 }
1720 else
1721 preconditioner.vmult(*r, v);
1722
1723 res = r->l2_norm();
1724 if (additional_data.batched_mode)
1726 else
1728 this->iteration_status(accumulated_iterations, res, x);
1729
1731 break;
1732 }
1733
1734 // inner iteration doing at most as many steps as the size of the
1735 // Arnoldi basis
1736 unsigned int inner_iteration = 0;
1737 for (; (inner_iteration < basis_size &&
1740 {
1742 // yet another alias
1744
1746 {
1747 A.vmult(p, basis_vectors[inner_iteration]);
1748 preconditioner.vmult(vv, p);
1749 }
1750 else
1751 {
1752 preconditioner.vmult(p, basis_vectors[inner_iteration]);
1753 A.vmult(vv, p);
1754 }
1755
1756 res =
1757 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1,
1760 re_orthogonalize_signal);
1761
1762 if (use_default_residual)
1763 {
1764 if (additional_data.batched_mode)
1766 solver_control.check(accumulated_iterations, res);
1767 else
1769 this->iteration_status(accumulated_iterations, res, x);
1770 }
1771 else
1772 {
1773 if (!additional_data.batched_mode)
1774 deallog << "default_res=" << res << std::endl;
1775
1776 *x_ = x;
1777 const Vector<double> &projected_solution =
1778 arnoldi_process.solve_projected_system(false);
1779
1781 for (unsigned int i = 0; i < inner_iteration + 1; ++i)
1782 x_->add(projected_solution(i), basis_vectors[i]);
1783 else
1784 {
1785 p = 0.;
1786 for (unsigned int i = 0; i < inner_iteration + 1; ++i)
1787 p.add(projected_solution(i), basis_vectors[i]);
1788 preconditioner.vmult(*r, p);
1789 x_->add(1., *r);
1790 };
1791 A.vmult(*r, *x_);
1792 r->sadd(-1., 1., b);
1793
1794 // Now *r contains the unpreconditioned residual!!
1796 {
1797 res = r->l2_norm();
1799 this->iteration_status(accumulated_iterations, res, x);
1800 }
1801 else
1802 {
1803 preconditioner.vmult(*x_, *r);
1804 res = x_->l2_norm();
1805
1806 if (additional_data.batched_mode)
1808 solver_control.check(accumulated_iterations, res);
1809 else
1811 this->iteration_status(accumulated_iterations, res, x);
1812 }
1813 }
1814 }
1815
1816 // end of inner iteration; now update the global solution vector x with
1817 // the solution of the projected system (least-squares solution)
1818 const Vector<double> &projected_solution =
1819 arnoldi_process.solve_projected_system(true);
1820
1821 if (do_eigenvalues)
1822 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
1824 all_eigenvalues_signal,
1825 all_hessenberg_signal,
1826 condition_number_signal);
1827
1829 ::internal::SolverGMRESImplementation::add(
1830 x,
1832 projected_solution,
1834 false,
1835 arnoldi_process.vector_ptrs);
1836 else
1837 {
1838 ::internal::SolverGMRESImplementation::add(
1839 p,
1841 projected_solution,
1843 true,
1844 arnoldi_process.vector_ptrs);
1845 preconditioner.vmult(v, p);
1846 x.add(1., v);
1847 }
1848
1849 // in the last round, print the eigenvalues from the last Arnoldi step
1851 {
1852 if (do_eigenvalues)
1853 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
1855 eigenvalues_signal,
1856 hessenberg_signal,
1857 condition_number_signal);
1858
1859 if (!additional_data.batched_mode && !krylov_space_signal.empty())
1860 krylov_space_signal(basis_vectors);
1861
1862 // end of outer iteration. restart if no convergence and the number of
1863 // iterations is not exceeded
1864 }
1865 }
1867
1868 // in case of failure: throw exception
1871}
1872
1873
1874
1875template <typename VectorType>
1877boost::signals2::connection
1879 const std::function<void(double)> &slot,
1880 const bool every_iteration)
1881{
1882 if (every_iteration)
1883 {
1884 return all_condition_numbers_signal.connect(slot);
1885 }
1886 else
1887 {
1888 return condition_number_signal.connect(slot);
1889 }
1890}
1891
1892
1893
1894template <typename VectorType>
1896boost::signals2::connection SolverGMRES<VectorType>::connect_eigenvalues_slot(
1897 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1898 const bool every_iteration)
1899{
1900 if (every_iteration)
1901 {
1902 return all_eigenvalues_signal.connect(slot);
1903 }
1904 else
1905 {
1906 return eigenvalues_signal.connect(slot);
1907 }
1908}
1909
1910
1911
1912template <typename VectorType>
1914boost::signals2::connection SolverGMRES<VectorType>::connect_hessenberg_slot(
1915 const std::function<void(const FullMatrix<double> &)> &slot,
1916 const bool every_iteration)
1917{
1918 if (every_iteration)
1919 {
1920 return all_hessenberg_signal.connect(slot);
1921 }
1922 else
1923 {
1924 return hessenberg_signal.connect(slot);
1925 }
1926}
1927
1928
1929
1930template <typename VectorType>
1932boost::signals2::connection SolverGMRES<VectorType>::connect_krylov_space_slot(
1933 const std::function<void(
1935{
1936 return krylov_space_signal.connect(slot);
1937}
1938
1939
1940
1941template <typename VectorType>
1943boost::signals2::connection
1945 const std::function<void(int)> &slot)
1946{
1947 return re_orthogonalize_signal.connect(slot);
1948}
1949
1950
1951
1952template <typename VectorType>
1955{
1956 // dummy implementation. this function is not needed for the present
1957 // implementation of gmres
1959 return 0;
1960}
1961
1962
1963
1964//----------------------------------------------------------------------//
1965
1966template <typename VectorType>
1970 const AdditionalData &data)
1971 : SolverBase<VectorType>(cn, mem)
1972 , additional_data(data)
1973{}
1974
1975
1976
1977template <typename VectorType>
1980 const AdditionalData &data)
1982 , additional_data(data)
1983{}
1984
1985
1986
1987template <typename VectorType>
1989template <typename MatrixType, typename PreconditionerType>
1993void SolverFGMRES<VectorType>::solve(const MatrixType &A,
1994 VectorType &x,
1995 const VectorType &b,
1996 const PreconditionerType &preconditioner)
1997{
1998 LogStream::Prefix prefix("FGMRES");
1999
2001
2002 const unsigned int basis_size = additional_data.max_basis_size;
2003
2004 // Generate an object where basis vectors are stored.
2006 basis_size + 1, this->memory);
2008 basis_size, this->memory);
2009
2010 // number of the present iteration; this number is not reset to zero upon a
2011 // restart
2012 unsigned int accumulated_iterations = 0;
2013
2014 // matrix used for the orthogonalization process later
2015 arnoldi_process.initialize(additional_data.orthogonalization_strategy,
2016 basis_size,
2017 false);
2018
2019 // Iteration starts here
2020 double res = std::numeric_limits<double>::lowest();
2021
2022 do
2023 {
2024 // Compute the residual. If 'x' is the zero vector, then we can bypass
2025 // the full computation. But 'x' is only likely to be the zero vector if
2026 // that's what the user provided as the starting guess, so it's only
2027 // worth checking for this in the first iteration. (Calling all_zero()
2028 // costs as much in memory transfer and communication as computing the
2029 // norm of a vector.)
2030 if (accumulated_iterations == 0 && x.all_zero())
2031 v(0, x) = b;
2032 else
2033 {
2034 A.vmult(v(0, x), x);
2035 v[0].sadd(-1., 1., b);
2036 }
2037
2038 res = arnoldi_process.orthonormalize_nth_vector(0, v);
2039 iteration_state = this->iteration_status(accumulated_iterations, res, x);
2041 break;
2042
2043 unsigned int inner_iteration = 0;
2044 for (; (inner_iteration < basis_size &&
2047 {
2048 preconditioner.vmult(z(inner_iteration, x), v[inner_iteration]);
2049 A.vmult(v(inner_iteration + 1, x), z[inner_iteration]);
2050
2051 res =
2052 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1, v);
2053
2054 // check convergence. note that the vector 'x' we pass to the
2055 // criterion is not the final solution we compute if we
2056 // decide to jump out of the iteration (we update 'x' again
2057 // right after the current loop)
2059 this->iteration_status(++accumulated_iterations, res, x);
2060 }
2061
2062 // Solve triangular system with projected quantities and update solution
2063 // vector
2064 const Vector<double> &projected_solution =
2065 arnoldi_process.solve_projected_system(true);
2066 ::internal::SolverGMRESImplementation::add(
2067 x,
2069 projected_solution,
2070 z,
2071 false,
2072 arnoldi_process.vector_ptrs);
2073 }
2075
2076 // in case of failure: throw exception
2078 AssertThrow(false,
2080}
2081
2082#endif // DOXYGEN
2083
2085
2086#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)
SolverFGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverGMRES(const SolverGMRES< VectorType > &)=delete
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::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
static void compute_eigs_and_cond(const FullMatrix< double > &H_orig, const unsigned int n, 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)
SolverGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
boost::signals2::connection connect_hessenberg_slot(const std::function< void(const FullMatrix< double > &)> &slot, const bool every_iteration=true)
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
virtual double criterion()
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< std::complex< double > > &)> &slot, const bool every_iteration=false)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
friend class Tensor
Definition tensor.h:865
const FullMatrix< double > & get_hessenberg_matrix() const
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
double orthonormalize_nth_vector(const unsigned int n, TmpVectors< VectorType > &orthogonal_vectors, const unsigned int accumulated_iterations=0, const boost::signals2::signal< void(int)> &reorthogonalize_signal=boost::signals2::signal< void(int)>())
const Vector< double > & solve_projected_system(const bool orthogonalization_finished)
void initialize(const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy, const unsigned int max_basis_size, const bool force_reorthogonalization)
std::vector< std::pair< double, double > > givens_rotations
double do_givens_rotation(const bool delayed_reorthogonalization, const int col, FullMatrix< double > &matrix, std::vector< std::pair< double, double > > &rotations, Vector< double > &rhs)
std::vector< typename VectorMemory< VectorType >::Pointer > data
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
VectorType & operator[](const unsigned int i) const
VectorType & operator()(const unsigned int i, const VectorType &temp)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:190
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:513
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
LogStream deallog
Definition logstream.cc:38
std::vector< index_type > data
Definition mpi.cc:740
types::global_dof_index locally_owned_size
Definition mpi.cc:827
@ matrix
Contents is actually a matrix.
constexpr char A
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:49
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
void do_Tvmult_add(const unsigned int n_vectors, const std::size_t locally_owned_size, const Number *current_vector, const std::vector< const Number * > &orthogonal_vectors, Vector< double > &h)
double do_subtract_and_norm(const unsigned int n_vectors, const std::size_t locally_owned_size, const std::vector< const Number * > &orthogonal_vectors, const Vector< double > &h, Number *current_vector)
void do_add(const unsigned int n_vectors, const std::size_t locally_owned_size, const std::vector< const Number * > &tmp_vectors, const Vector< double > &h, const bool zero_out, Number *output)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
AdditionalData(const unsigned int max_basis_size=30, const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy=LinearAlgebra::OrthogonalizationStrategy::delayed_classical_gram_schmidt)
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
AdditionalData(const unsigned int max_basis_size=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::delayed_classical_gram_schmidt)
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)