deal.II version GIT relicensing-3689-g8ed7910be7 2025-07-08 08:50:00+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 - 2025 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 <utility>
41#include <vector>
42
44
45// forward declarations
46#ifndef DOXYGEN
47namespace LinearAlgebra
48{
49 namespace distributed
50 {
51 template <typename, typename>
52 class Vector;
53 } // namespace distributed
54} // namespace LinearAlgebra
55#endif
56
62namespace internal
63{
67 namespace SolverGMRESImplementation
68 {
76 template <typename VectorType>
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
127
128
138 template <typename Number>
140 {
141 public:
146 void
149 const unsigned int max_basis_size,
150 const bool force_reorthogonalization);
151
174 template <typename VectorType>
175 double
177 const unsigned int n,
178 TmpVectors<VectorType> &orthogonal_vectors,
179 const unsigned int accumulated_iterations = 0,
180 const boost::signals2::signal<void(int)> &reorthogonalize_signal =
181 boost::signals2::signal<void(int)>());
182
190 const Vector<double> &
191 solve_projected_system(const bool orthogonalization_finished);
192
197 const FullMatrix<double> &
199
203 std::vector<const Number *> vector_ptrs;
204
205 private:
210
217
222 std::vector<std::pair<double, double>> givens_rotations;
223
228
234
239
249
254
282 double
283 do_givens_rotation(const bool delayed_reorthogonalization,
284 const int col,
285 FullMatrix<double> &matrix,
286 std::vector<std::pair<double, double>> &rotations,
287 Vector<double> &rhs);
288 };
289 } // namespace SolverGMRESImplementation
290} // namespace internal
291
292
293
413template <typename VectorType = Vector<double>>
415class SolverGMRES : public SolverBase<VectorType>
416{
417public:
422 {
433 explicit AdditionalData(const unsigned int max_basis_size = 30,
434 const bool right_preconditioning = false,
435 const bool use_default_residual = true,
436 const bool force_re_orthogonalization = false,
437 const bool batched_mode = false,
439 orthogonalization_strategy =
441 delayed_classical_gram_schmidt);
442
454 unsigned int max_n_tmp_vectors;
455
463 unsigned int max_basis_size;
464
473
478
486
494
499 };
500
507
513
518
522 template <typename MatrixType, typename PreconditionerType>
526 void solve(const MatrixType &A,
527 VectorType &x,
528 const VectorType &b,
529 const PreconditionerType &preconditioner);
530
537 boost::signals2::connection
538 connect_condition_number_slot(const std::function<void(double)> &slot,
539 const bool every_iteration = false);
540
547 boost::signals2::connection
548 connect_eigenvalues_slot(
549 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
550 const bool every_iteration = false);
551
559 boost::signals2::connection
560 connect_hessenberg_slot(
561 const std::function<void(const FullMatrix<double> &)> &slot,
562 const bool every_iteration = true);
563
570 boost::signals2::connection
571 connect_krylov_space_slot(
572 const std::function<
573 void(const internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
574 &slot);
575
576
581 boost::signals2::connection
582 connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
583
584
585 DeclException1(ExcTooFewTmpVectors,
586 int,
587 << "The number of temporary vectors you gave (" << arg1
588 << ") is too small. It should be at least 10 for "
589 << "any results, and much more for reasonable ones.");
590
591protected:
595 AdditionalData additional_data;
596
601 boost::signals2::signal<void(double)> condition_number_signal;
602
607 boost::signals2::signal<void(double)> all_condition_numbers_signal;
608
613 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
614 eigenvalues_signal;
615
620 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
621 all_eigenvalues_signal;
622
627 boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
628
633 boost::signals2::signal<void(const FullMatrix<double> &)>
634 all_hessenberg_signal;
635
640 boost::signals2::signal<void(
641 const internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
642 krylov_space_signal;
643
648 boost::signals2::signal<void(int)> re_orthogonalize_signal;
649
655 SolverControl &solver_control;
656
660 virtual double
661 criterion();
662
669 static void
670 compute_eigs_and_cond(
671 const FullMatrix<double> &H_orig,
672 const unsigned int n,
673 const boost::signals2::signal<
674 void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
675 const boost::signals2::signal<void(const FullMatrix<double> &)>
676 &hessenberg_signal,
677 const boost::signals2::signal<void(double)> &cond_signal);
678
683 internal::SolverGMRESImplementation::ArnoldiProcess<
684 typename VectorType::value_type>
685 arnoldi_process;
686};
687
688
752template <typename VectorType = Vector<double>>
753DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
754class SolverMPGMRES : public SolverBase<VectorType>
755{
756public:
761 {
765 explicit AdditionalData(const unsigned int max_basis_size = 30,
767 orthogonalization_strategy =
769 delayed_classical_gram_schmidt,
770 const bool use_truncated_mpgmres_strategy = true)
771 : max_basis_size(max_basis_size)
772 , orthogonalization_strategy(orthogonalization_strategy)
773 , use_truncated_mpgmres_strategy(use_truncated_mpgmres_strategy)
774 {}
775
779 unsigned int max_basis_size;
780
785
796 };
797
804
811
815 template <typename MatrixType, typename... PreconditionerTypes>
819 void solve(const MatrixType &A,
820 VectorType &x,
821 const VectorType &b,
822 const PreconditionerTypes &...preconditioners);
823
824protected:
832 {
833 fgmres,
834 truncated_mpgmres,
835 full_mpgmres,
836 };
837
841 template <typename MatrixType, typename... PreconditionerTypes>
842 void
843 solve_internal(const MatrixType &A,
844 VectorType &x,
845 const VectorType &b,
846 const IndexingStrategy &indexing_strategy,
847 const PreconditionerTypes &...preconditioners);
848
849private:
854
860 typename VectorType::value_type>
862};
863
864
865
891template <typename VectorType = Vector<double>>
893class SolverFGMRES : public SolverMPGMRES<VectorType>
894{
895public:
900 {
904 explicit AdditionalData( //
905 const unsigned int max_basis_size = 30,
907 orthogonalization_strategy = LinearAlgebra::OrthogonalizationStrategy::
908 delayed_classical_gram_schmidt)
909 : max_basis_size(max_basis_size)
910 , orthogonalization_strategy(orthogonalization_strategy)
911 {}
912
916 unsigned int max_basis_size;
917
922 };
923
924
931
938
942 template <typename MatrixType, typename... PreconditionerTypes>
946 void solve(const MatrixType &A,
947 VectorType &x,
948 const VectorType &b,
949 const PreconditionerTypes &...preconditioners);
950};
951
955/* --------------------- Inline and template functions ------------------- */
956
957#ifndef DOXYGEN
958
959template <typename VectorType>
962 const unsigned int max_basis_size,
963 const bool right_preconditioning,
964 const bool use_default_residual,
965 const bool force_re_orthogonalization,
966 const bool batched_mode,
967 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy)
968 : max_n_tmp_vectors(0)
969 , max_basis_size(max_basis_size)
970 , right_preconditioning(right_preconditioning)
971 , use_default_residual(use_default_residual)
972 , force_re_orthogonalization(force_re_orthogonalization)
973 , batched_mode(batched_mode)
974 , orthogonalization_strategy(orthogonalization_strategy)
975{
976 Assert(max_basis_size >= 1,
977 ExcMessage("SolverGMRES needs at least one vector in the "
978 "Arnoldi basis."));
979}
980
981
982
983template <typename VectorType>
987 const AdditionalData &data)
988 : SolverBase<VectorType>(cn, mem)
989 , additional_data(data)
990 , solver_control(cn)
991{}
992
993
994
995template <typename VectorType>
998 const AdditionalData &data)
1000 , additional_data(data)
1001 , solver_control(cn)
1002{}
1003
1004
1005
1006namespace internal
1007{
1008 namespace SolverGMRESImplementation
1009 {
1010 template <typename VectorType>
1011 inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
1013 : mem(vmem)
1014 , data(max_size)
1015 {}
1016
1017
1018
1019 template <typename VectorType>
1020 inline VectorType &
1021 TmpVectors<VectorType>::operator[](const unsigned int i) const
1022 {
1023 AssertIndexRange(i, data.size());
1024
1025 Assert(data[i] != nullptr, ExcNotInitialized());
1026 return *data[i];
1027 }
1028
1029
1030
1031 template <typename VectorType>
1032 inline VectorType &
1033 TmpVectors<VectorType>::operator()(const unsigned int i,
1034 const VectorType &temp)
1035 {
1036 AssertIndexRange(i, data.size());
1037 if (data[i] == nullptr)
1038 {
1039 data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
1040 data[i]->reinit(temp, true);
1041 }
1042 return *data[i];
1043 }
1044
1045
1046
1047 template <typename VectorType>
1048 unsigned int
1049 TmpVectors<VectorType>::size() const
1050 {
1051 return (data.size() > 0 ? data.size() - 1 : 0);
1052 }
1053
1054
1055
1056 template <typename VectorType, typename Enable = void>
1057 struct is_dealii_compatible_vector;
1058
1059 template <typename VectorType>
1060 struct is_dealii_compatible_vector<
1061 VectorType,
1062 std::enable_if_t<!internal::is_block_vector<VectorType>>>
1063 {
1064 static constexpr bool value =
1065 std::is_same_v<
1066 VectorType,
1067 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
1069 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
1070 };
1071
1072
1073
1074 template <typename VectorType>
1075 struct is_dealii_compatible_vector<
1076 VectorType,
1077 std::enable_if_t<internal::is_block_vector<VectorType>>>
1078 {
1079 static constexpr bool value =
1080 std::is_same_v<
1081 typename VectorType::BlockType,
1082 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
1084 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
1085 };
1086
1087
1088
1089 template <typename VectorType,
1090 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
1091 * = nullptr>
1092 unsigned int
1093 n_blocks(const VectorType &)
1094 {
1095 return 1;
1096 }
1097
1098
1099
1100 template <typename VectorType,
1101 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
1102 nullptr>
1103 unsigned int
1104 n_blocks(const VectorType &vector)
1105 {
1106 return vector.n_blocks();
1107 }
1108
1109
1110
1111 template <typename VectorType,
1112 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
1113 * = nullptr>
1114 VectorType &
1115 block(VectorType &vector, const unsigned int b)
1116 {
1117 AssertDimension(b, 0);
1118 return vector;
1119 }
1120
1121
1122
1123 template <typename VectorType,
1124 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
1125 * = nullptr>
1126 const VectorType &
1127 block(const VectorType &vector, const unsigned int b)
1128 {
1129 AssertDimension(b, 0);
1130 return vector;
1131 }
1132
1133
1134
1135 template <typename VectorType,
1136 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
1137 nullptr>
1138 typename VectorType::BlockType &
1139 block(VectorType &vector, const unsigned int b)
1140 {
1141 return vector.block(b);
1142 }
1143
1144
1145
1146 template <typename VectorType,
1147 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
1148 nullptr>
1149 const typename VectorType::BlockType &
1150 block(const VectorType &vector, const unsigned int b)
1151 {
1152 return vector.block(b);
1153 }
1154
1155
1156
1157 template <bool delayed_reorthogonalization,
1158 typename VectorType,
1159 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1160 VectorType> * = nullptr>
1161 void
1162 Tvmult_add(const unsigned int n,
1163 const VectorType &vv,
1164 const TmpVectors<VectorType> &orthogonal_vectors,
1165 Vector<double> &h,
1166 std::vector<const typename VectorType::value_type *> &)
1167 {
1168 for (unsigned int i = 0; i < n; ++i)
1169 {
1170 h(i) += vv * orthogonal_vectors[i];
1171 if (delayed_reorthogonalization)
1172 h(n + i) += orthogonal_vectors[i] * orthogonal_vectors[n - 1];
1173 }
1174 if (delayed_reorthogonalization)
1175 h(n + n) += vv * vv;
1176 }
1177
1178
1179
1180 // worker method for deal.II's vector types implemented in .cc file
1181 template <bool delayed_reorthogonalization, typename Number>
1182 void
1183 do_Tvmult_add(const unsigned int n_vectors,
1184 const std::size_t locally_owned_size,
1185 const Number *current_vector,
1186 const std::vector<const Number *> &orthogonal_vectors,
1187 Vector<double> &b);
1188
1189
1190
1191 template <bool delayed_reorthogonalization,
1192 typename VectorType,
1193 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1194 VectorType> * = nullptr>
1195 void
1196 Tvmult_add(
1197 const unsigned int n,
1198 const VectorType &vv,
1199 const TmpVectors<VectorType> &orthogonal_vectors,
1200 Vector<double> &h,
1201 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1202 {
1203 for (unsigned int b = 0; b < n_blocks(vv); ++b)
1204 {
1205 vector_ptrs.resize(n);
1206 for (unsigned int i = 0; i < n; ++i)
1207 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
1208
1209 do_Tvmult_add<delayed_reorthogonalization>(n,
1210 block(vv, b).end() -
1211 block(vv, b).begin(),
1212 block(vv, b).begin(),
1213 vector_ptrs,
1214 h);
1215 }
1216
1217 Utilities::MPI::sum(h, block(vv, 0).get_mpi_communicator(), h);
1218 }
1219
1220
1221
1222 template <bool delayed_reorthogonalization,
1223 typename VectorType,
1224 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1225 VectorType> * = nullptr>
1226 double
1227 subtract_and_norm(const unsigned int n,
1228 const TmpVectors<VectorType> &orthogonal_vectors,
1229 const Vector<double> &h,
1230 VectorType &vv,
1231 std::vector<const typename VectorType::value_type *> &)
1232 {
1233 Assert(n > 0, ExcInternalError());
1234
1235 VectorType &last_vector =
1236 const_cast<VectorType &>(orthogonal_vectors[n - 1]);
1237 for (unsigned int i = 0; i < n - 1; ++i)
1238 {
1239 if (delayed_reorthogonalization && i + 2 < n)
1240 last_vector.add(-h(n + i), orthogonal_vectors[i]);
1241 vv.add(-h(i), orthogonal_vectors[i]);
1242 }
1243
1244 if (delayed_reorthogonalization)
1245 {
1246 if (n > 1)
1247 last_vector.sadd(1. / h(n + n - 1),
1248 -h(n + n - 2) / h(n + n - 1),
1249 orthogonal_vectors[n - 2]);
1250
1251 // h(n + n) is lucky breakdown
1252 const double scaling_factor_vv = h(n + n) > 0.0 ?
1253 1. / (h(n + n - 1) * h(n + n)) :
1254 1. / (h(n + n - 1) * h(n + n - 1));
1255 vv.sadd(scaling_factor_vv,
1256 -h(n - 1) * scaling_factor_vv,
1257 last_vector);
1258
1259 // the delayed reorthogonalization computes the norm from other
1260 // quantities
1261 return std::numeric_limits<double>::signaling_NaN();
1262 }
1263 else
1264 return std::sqrt(
1265 vv.add_and_dot(-h(n - 1), orthogonal_vectors[n - 1], vv));
1266 }
1267
1268
1269
1270 // worker method for deal.II's vector types implemented in .cc file
1271 template <bool delayed_reorthogonalization, typename Number>
1272 double
1273 do_subtract_and_norm(const unsigned int n_vectors,
1274 const std::size_t locally_owned_size,
1275 const std::vector<const Number *> &orthogonal_vectors,
1276 const Vector<double> &h,
1277 Number *current_vector);
1278
1279
1280
1281 template <bool delayed_reorthogonalization,
1282 typename VectorType,
1283 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1284 VectorType> * = nullptr>
1285 double
1286 subtract_and_norm(
1287 const unsigned int n,
1288 const TmpVectors<VectorType> &orthogonal_vectors,
1289 const Vector<double> &h,
1290 VectorType &vv,
1291 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1292 {
1293 double norm_vv_temp = 0.0;
1294
1295 for (unsigned int b = 0; b < n_blocks(vv); ++b)
1296 {
1297 vector_ptrs.resize(n);
1298 for (unsigned int i = 0; i < n; ++i)
1299 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
1300
1301 norm_vv_temp += do_subtract_and_norm<delayed_reorthogonalization>(
1302 n,
1303 block(vv, b).end() - block(vv, b).begin(),
1304 vector_ptrs,
1305 h,
1306 block(vv, b).begin());
1307 }
1308
1309 return std::sqrt(
1310 Utilities::MPI::sum(norm_vv_temp, block(vv, 0).get_mpi_communicator()));
1311 }
1312
1313
1314
1315 template <typename VectorType,
1316 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1317 VectorType> * = nullptr>
1318 void
1319 add(VectorType &p,
1320 const unsigned int n,
1321 const Vector<double> &h,
1322 const TmpVectors<VectorType> &tmp_vectors,
1323 const bool zero_out,
1324 std::vector<const typename VectorType::value_type *> &)
1325 {
1326 if (zero_out)
1327 p.equ(h(0), tmp_vectors[0]);
1328 else
1329 p.add(h(0), tmp_vectors[0]);
1330
1331 for (unsigned int i = 1; i < n; ++i)
1332 p.add(h(i), tmp_vectors[i]);
1333 }
1334
1335
1336
1337 // worker method for deal.II's vector types implemented in .cc file
1338 template <typename Number>
1339 void
1340 do_add(const unsigned int n_vectors,
1341 const std::size_t locally_owned_size,
1342 const std::vector<const Number *> &tmp_vectors,
1343 const Vector<double> &h,
1344 const bool zero_out,
1345 Number *output);
1346
1347
1348
1349 template <typename VectorType,
1350 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1351 VectorType> * = nullptr>
1352 void
1353 add(VectorType &p,
1354 const unsigned int n,
1355 const Vector<double> &h,
1356 const TmpVectors<VectorType> &tmp_vectors,
1357 const bool zero_out,
1358 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1359 {
1360 for (unsigned int b = 0; b < n_blocks(p); ++b)
1361 {
1362 vector_ptrs.resize(n);
1363 for (unsigned int i = 0; i < n; ++i)
1364 vector_ptrs[i] = block(tmp_vectors[i], b).begin();
1365 do_add(n,
1366 block(p, b).end() - block(p, b).begin(),
1367 vector_ptrs,
1368 h,
1369 zero_out,
1370 block(p, b).begin());
1371 }
1372 }
1373
1374
1375
1376 template <typename Number>
1377 inline void
1378 ArnoldiProcess<Number>::initialize(
1379 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy,
1380 const unsigned int basis_size,
1381 const bool force_reorthogonalization)
1382 {
1383 this->orthogonalization_strategy = orthogonalization_strategy;
1384 this->do_reorthogonalization = force_reorthogonalization;
1385
1386 hessenberg_matrix.reinit(basis_size + 1, basis_size);
1387 triangular_matrix.reinit(basis_size + 1, basis_size, true);
1388
1389 // some additional vectors, also used in the orthogonalization
1390 projected_rhs.reinit(basis_size + 1, true);
1391 givens_rotations.reserve(basis_size);
1392
1393 if (orthogonalization_strategy ==
1396 h.reinit(2 * basis_size + 3);
1397 else
1398 h.reinit(basis_size + 1);
1399 }
1400
1401
1402
1403 template <typename Number>
1404 template <typename VectorType>
1405 inline double
1406 ArnoldiProcess<Number>::orthonormalize_nth_vector(
1407 const unsigned int n,
1408 TmpVectors<VectorType> &orthogonal_vectors,
1409 const unsigned int accumulated_iterations,
1410 const boost::signals2::signal<void(int)> &reorthogonalize_signal)
1411 {
1412 AssertIndexRange(n, hessenberg_matrix.m());
1413 AssertIndexRange(n, orthogonal_vectors.size() + 1);
1414
1415 VectorType &vv = orthogonal_vectors[n];
1416
1417 double residual_estimate = std::numeric_limits<double>::signaling_NaN();
1418 if (n == 0)
1419 {
1420 givens_rotations.clear();
1421 residual_estimate = vv.l2_norm();
1422 if (residual_estimate != 0.)
1423 vv /= residual_estimate;
1424 projected_rhs(0) = residual_estimate;
1425 }
1426 else if (orthogonalization_strategy ==
1429 {
1430 // The algorithm implemented in the following few lines is algorithm
1431 // 4 of Bielich et al. (2022).
1432
1433 // To avoid un-scaled numbers as appearing with the original
1434 // algorithm by Bielich et al., we use a preliminary scaling of the
1435 // last vector. This will be corrected in the delayed step.
1436 const double previous_scaling = n > 0 ? h(n + n - 2) : 1.;
1437
1438 // Reset h to zero
1439 h.reinit(n + n + 1);
1440
1441 // global reduction
1442 Tvmult_add<true>(n, vv, orthogonal_vectors, h, vector_ptrs);
1443
1444 // delayed correction terms
1445 double tmp = 0;
1446 for (unsigned int i = 0; i < n - 1; ++i)
1447 tmp += h(n + i) * h(n + i);
1448 const double alpha_j = h(n + n - 1) > tmp ?
1449 std::sqrt(h(n + n - 1) - tmp) :
1450 std::sqrt(h(n + n - 1));
1451 h(n + n - 1) = alpha_j;
1452
1453 tmp = 0;
1454 for (unsigned int i = 0; i < n - 1; ++i)
1455 tmp += h(i) * h(n + i);
1456 h(n - 1) = (h(n - 1) - tmp) / alpha_j;
1457
1458 // representation of H(j-1)
1459 if (n > 1)
1460 {
1461 for (unsigned int i = 0; i < n - 1; ++i)
1462 hessenberg_matrix(i, n - 2) += h(n + i) * previous_scaling;
1463 hessenberg_matrix(n - 1, n - 2) = alpha_j * previous_scaling;
1464 }
1465 for (unsigned int i = 0; i < n; ++i)
1466 {
1467 double sum = 0;
1468 for (unsigned int j = (i == 0 ? 0 : i - 1); j < n - 1; ++j)
1469 sum += hessenberg_matrix(i, j) * h(n + j);
1470 hessenberg_matrix(i, n - 1) = (h(i) - sum) / alpha_j;
1471 }
1472
1473 // compute norm estimate for approximate convergence criterion
1474 // (value of norm to be corrected in next iteration)
1475 double sum = 0;
1476 for (unsigned int i = 0; i < n - 1; ++i)
1477 sum += h(i) * h(i);
1478 sum += (2. - 1.) * h(n - 1) * h(n - 1);
1479 hessenberg_matrix(n, n - 1) =
1480 std::sqrt(std::abs(h(n + n) - sum)) / alpha_j;
1481
1482 // projection and delayed reorthogonalization. We scale the vector
1483 // vv here by the preliminary norm to avoid working with too large
1484 // values and correct the actual norm in the Hessenberg matrix in
1485 // high precision in the next iteration.
1486 h(n + n) = hessenberg_matrix(n, n - 1);
1487 subtract_and_norm<true>(n, orthogonal_vectors, h, vv, vector_ptrs);
1488
1489 // transform new column of upper Hessenberg matrix into upper
1490 // triangular form by computing the respective factor
1491 residual_estimate = do_givens_rotation(
1492 true, n - 2, triangular_matrix, givens_rotations, projected_rhs);
1493 }
1494 else
1495 {
1496 // need initial norm for detection of re-orthogonalization, see below
1497 double norm_vv = 0.0;
1498 double norm_vv_start = 0;
1499 const bool consider_reorthogonalize =
1500 (do_reorthogonalization == false) && (n % 5 == 0);
1501 if (consider_reorthogonalize)
1502 norm_vv_start = vv.l2_norm();
1503
1504 // Reset h to zero
1505 h.reinit(n);
1506
1507 // run two loops with index 0: orthogonalize, 1: reorthogonalize
1508 for (unsigned int c = 0; c < 2; ++c)
1509 {
1510 // Orthogonalization
1511 if (orthogonalization_strategy ==
1514 {
1515 double htmp = vv * orthogonal_vectors[0];
1516 h(0) += htmp;
1517 for (unsigned int i = 1; i < n; ++i)
1518 {
1519 htmp = vv.add_and_dot(-htmp,
1520 orthogonal_vectors[i - 1],
1521 orthogonal_vectors[i]);
1522 h(i) += htmp;
1523 }
1524
1525 norm_vv = std::sqrt(
1526 vv.add_and_dot(-htmp, orthogonal_vectors[n - 1], vv));
1527 }
1528 else if (orthogonalization_strategy ==
1531 {
1532 Tvmult_add<false>(n, vv, orthogonal_vectors, h, vector_ptrs);
1533 norm_vv = subtract_and_norm<false>(
1534 n, orthogonal_vectors, h, vv, vector_ptrs);
1535 }
1536 else
1537 {
1539 }
1540
1541 if (c == 1)
1542 break; // reorthogonalization already performed -> finished
1543
1544 // Re-orthogonalization if loss of orthogonality detected. For the
1545 // test, use a strategy discussed in C. T. Kelley, Iterative
1546 // Methods for Linear and Nonlinear Equations, SIAM, Philadelphia,
1547 // 1995: Compare the norm of vv after orthogonalization with its
1548 // norm when starting the orthogonalization. If vv became very
1549 // small (here: less than the square root of the machine precision
1550 // times 10), it is almost in the span of the previous vectors,
1551 // which indicates loss of precision.
1552 if (consider_reorthogonalize)
1553 {
1554 if (norm_vv >
1555 10. * norm_vv_start *
1556 std::sqrt(std::numeric_limits<
1557 typename VectorType::value_type>::epsilon()))
1558 break;
1559
1560 else
1561 {
1562 do_reorthogonalization = true;
1563 if (!reorthogonalize_signal.empty())
1564 reorthogonalize_signal(accumulated_iterations);
1565 }
1566 }
1567
1568 if (do_reorthogonalization == false)
1569 break; // no reorthogonalization needed -> finished
1570 }
1571
1572 for (unsigned int i = 0; i < n; ++i)
1573 hessenberg_matrix(i, n - 1) = h(i);
1574 hessenberg_matrix(n, n - 1) = norm_vv;
1575
1576 // norm_vv is a lucky breakdown, the solver will reach convergence,
1577 // but we must not divide by zero here.
1578 if (norm_vv != 0)
1579 vv /= norm_vv;
1580
1581 residual_estimate = do_givens_rotation(
1582 false, n - 1, triangular_matrix, givens_rotations, projected_rhs);
1583 }
1584
1585 return residual_estimate;
1586 }
1587
1588
1589
1590 template <typename Number>
1591 inline double
1592 ArnoldiProcess<Number>::do_givens_rotation(
1593 const bool delayed_reorthogonalization,
1594 const int col,
1595 FullMatrix<double> &matrix,
1596 std::vector<std::pair<double, double>> &rotations,
1597 Vector<double> &rhs)
1598 {
1599 // for the delayed orthogonalization, we can only compute the column of
1600 // the previous iteration (as there will be correction terms added to the
1601 // present column for stability reasons), but we still want to compute
1602 // the residual estimate from the accumulated work; we therefore perform
1603 // givens rotations on two columns simultaneously
1604 if (delayed_reorthogonalization)
1605 {
1606 if (col >= 0)
1607 {
1608 AssertDimension(rotations.size(), static_cast<std::size_t>(col));
1609 matrix(0, col) = hessenberg_matrix(0, col);
1610 }
1611 double H_next = hessenberg_matrix(0, col + 1);
1612 for (int i = 0; i < col; ++i)
1613 {
1614 const double c = rotations[i].first;
1615 const double s = rotations[i].second;
1616 const double Hi = matrix(i, col);
1617 const double Hi1 = hessenberg_matrix(i + 1, col);
1618 H_next = -s * H_next + c * hessenberg_matrix(i + 1, col + 1);
1619 matrix(i, col) = c * Hi + s * Hi1;
1620 matrix(i + 1, col) = -s * Hi + c * Hi1;
1621 }
1622
1623 if (col >= 0)
1624 {
1625 const double H_col1 = hessenberg_matrix(col + 1, col);
1626 const double H_col = matrix(col, col);
1627 const double r = 1. / std::sqrt(H_col * H_col + H_col1 * H_col1);
1628 rotations.emplace_back(H_col * r, H_col1 * r);
1629 matrix(col, col) =
1630 rotations[col].first * H_col + rotations[col].second * H_col1;
1631
1632 rhs(col + 1) = -rotations[col].second * rhs(col);
1633 rhs(col) *= rotations[col].first;
1634
1635 H_next =
1636 -rotations[col].second * H_next +
1637 rotations[col].first * hessenberg_matrix(col + 1, col + 1);
1638 }
1639
1640 const double H_last = hessenberg_matrix(col + 2, col + 1);
1641 const double r = 1. / std::sqrt(H_next * H_next + H_last * H_last);
1642 return std::abs(H_last * r * rhs(col + 1));
1643 }
1644 else
1645 {
1646 AssertDimension(rotations.size(), static_cast<std::size_t>(col));
1647
1648 matrix(0, col) = hessenberg_matrix(0, col);
1649 for (int i = 0; i < col; ++i)
1650 {
1651 const double c = rotations[i].first;
1652 const double s = rotations[i].second;
1653 const double Hi = matrix(i, col);
1654 const double Hi1 = hessenberg_matrix(i + 1, col);
1655 matrix(i, col) = c * Hi + s * Hi1;
1656 matrix(i + 1, col) = -s * Hi + c * Hi1;
1657 }
1658
1659 const double Hi = matrix(col, col);
1660 const double Hi1 = hessenberg_matrix(col + 1, col);
1661 const double r = 1. / std::sqrt(Hi * Hi + Hi1 * Hi1);
1662 rotations.emplace_back(Hi * r, Hi1 * r);
1663 matrix(col, col) =
1664 rotations[col].first * Hi + rotations[col].second * Hi1;
1665
1666 rhs(col + 1) = -rotations[col].second * rhs(col);
1667 rhs(col) *= rotations[col].first;
1668
1669 return std::abs(rhs(col + 1));
1670 }
1671 }
1672
1673
1674
1675 template <typename Number>
1676 inline const Vector<double> &
1677 ArnoldiProcess<Number>::solve_projected_system(
1678 const bool orthogonalization_finished)
1679 {
1680 FullMatrix<double> tmp_triangular_matrix;
1681 Vector<double> tmp_rhs;
1682 FullMatrix<double> *matrix = &triangular_matrix;
1683 Vector<double> *rhs = &projected_rhs;
1684 unsigned int n = givens_rotations.size();
1685
1686 // If we solve with the delayed orthogonalization, we still need to
1687 // perform the elimination of the last column before we can solve the
1688 // projected system. We distinguish two cases, one where the
1689 // orthogonalization has finished (i.e., end of inner iteration in
1690 // GMRES) and we can safely overwrite the content of the tridiagonal
1691 // matrix and right hand side, and the case during the inner iterations,
1692 // where we need to create copies of the matrices in the QR
1693 // decomposition as well as the right hand side.
1694 if (orthogonalization_strategy ==
1697 {
1698 n += 1;
1699 if (!orthogonalization_finished)
1700 {
1701 tmp_triangular_matrix = triangular_matrix;
1702 tmp_rhs = projected_rhs;
1703 std::vector<std::pair<double, double>> tmp_givens_rotations(
1704 givens_rotations);
1705 do_givens_rotation(false,
1706 givens_rotations.size(),
1707 tmp_triangular_matrix,
1708 tmp_givens_rotations,
1709 tmp_rhs);
1710 matrix = &tmp_triangular_matrix;
1711 rhs = &tmp_rhs;
1712 }
1713 else
1714 do_givens_rotation(false,
1715 givens_rotations.size(),
1716 triangular_matrix,
1717 givens_rotations,
1718 projected_rhs);
1719 }
1720
1721 // Now solve the triangular system by backward substitution
1722 projected_solution.reinit(n);
1723 for (int i = n - 1; i >= 0; --i)
1724 {
1725 double s = (*rhs)(i);
1726 for (unsigned int j = i + 1; j < n; ++j)
1727 s -= projected_solution(j) * (*matrix)(i, j);
1728 projected_solution(i) = s / (*matrix)(i, i);
1729 AssertIsFinite(projected_solution(i));
1730 }
1731
1732 return projected_solution;
1733 }
1734
1735
1736
1737 template <typename Number>
1738 inline const FullMatrix<double> &
1739 ArnoldiProcess<Number>::get_hessenberg_matrix() const
1740 {
1741 return hessenberg_matrix;
1742 }
1743
1744
1745
1746 // A comparator for better printing eigenvalues
1747 inline bool
1748 complex_less_pred(const std::complex<double> &x,
1749 const std::complex<double> &y)
1750 {
1751 return x.real() < y.real() ||
1752 (x.real() == y.real() && x.imag() < y.imag());
1753 }
1754 } // namespace SolverGMRESImplementation
1755} // namespace internal
1756
1757
1758
1759template <typename VectorType>
1762 const FullMatrix<double> &H_orig,
1763 const unsigned int n,
1764 const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
1765 &eigenvalues_signal,
1766 const boost::signals2::signal<void(const FullMatrix<double> &)>
1767 &hessenberg_signal,
1768 const boost::signals2::signal<void(double)> &cond_signal)
1769{
1770 // Avoid copying the Hessenberg matrix if it isn't needed.
1771 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1772 !cond_signal.empty()) &&
1773 n > 0)
1774 {
1775 LAPACKFullMatrix<double> mat(n, n);
1776 for (unsigned int i = 0; i < n; ++i)
1777 for (unsigned int j = 0; j < n; ++j)
1778 mat(i, j) = H_orig(i, j);
1779 hessenberg_signal(H_orig);
1780 // Avoid computing eigenvalues if they are not needed.
1781 if (!eigenvalues_signal.empty())
1782 {
1783 // Copy mat so that we can compute svd below. Necessary since
1784 // compute_eigenvalues will leave mat in state
1785 // LAPACKSupport::unusable.
1786 LAPACKFullMatrix<double> mat_eig(mat);
1787 mat_eig.compute_eigenvalues();
1788 std::vector<std::complex<double>> eigenvalues(n);
1789 for (unsigned int i = 0; i < mat_eig.n(); ++i)
1790 eigenvalues[i] = mat_eig.eigenvalue(i);
1791 // Sort eigenvalues for nicer output.
1792 std::sort(eigenvalues.begin(),
1793 eigenvalues.end(),
1794 internal::SolverGMRESImplementation::complex_less_pred);
1795 eigenvalues_signal(eigenvalues);
1796 }
1797 // Calculate condition number, avoid calculating the svd if a slot
1798 // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
1799 if (!cond_signal.empty() && (mat.n() > 1))
1800 {
1801 mat.compute_svd();
1802 double condition_number =
1803 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1804 cond_signal(condition_number);
1805 }
1806 }
1807}
1808
1809
1810
1811template <typename VectorType>
1813template <typename MatrixType, typename PreconditionerType>
1817void SolverGMRES<VectorType>::solve(const MatrixType &A,
1818 VectorType &x,
1819 const VectorType &b,
1820 const PreconditionerType &preconditioner)
1821{
1822 std::unique_ptr<LogStream::Prefix> prefix;
1823 if (!additional_data.batched_mode)
1824 prefix = std::make_unique<LogStream::Prefix>("GMRES");
1825
1826 // extra call to std::max to placate static analyzers: coverity rightfully
1827 // complains that data.max_n_tmp_vectors - 2 may overflow
1828 const unsigned int basis_size =
1829 (additional_data.max_basis_size > 0 ?
1830 additional_data.max_basis_size :
1831 std::max(additional_data.max_n_tmp_vectors, 3u) - 2);
1832
1833 // Generate an object where basis vectors are stored.
1835 basis_size + 2, this->memory);
1836
1837 // number of the present iteration; this number is not reset to zero upon a
1838 // restart
1839 unsigned int accumulated_iterations = 0;
1840
1841 const bool do_eigenvalues =
1842 !additional_data.batched_mode &&
1843 (!condition_number_signal.empty() ||
1844 !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1845 !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1846 !all_hessenberg_signal.empty());
1847
1849 double res = std::numeric_limits<double>::lowest();
1850
1851 // switch to determine whether we want a left or a right preconditioner. at
1852 // present, left is default, but both ways are implemented
1853 const bool left_precondition = !additional_data.right_preconditioning;
1854
1855 // Per default the left preconditioned GMRES uses the preconditioned
1856 // residual and the right preconditioned GMRES uses the unpreconditioned
1857 // residual as stopping criterion.
1858 const bool use_default_residual = additional_data.use_default_residual;
1859
1860 // define an alias
1861 VectorType &p = basis_vectors(basis_size + 1, x);
1862
1863 // Following vectors are needed when we are not using the default residuals
1864 // as stopping criterion
1867 if (!use_default_residual)
1868 {
1869 r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1870 x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1871 r->reinit(x);
1872 x_->reinit(x);
1873 }
1874
1875 arnoldi_process.initialize(additional_data.orthogonalization_strategy,
1876 basis_size,
1877 additional_data.force_re_orthogonalization);
1878
1880 // outer iteration: loop until we either reach convergence or the maximum
1881 // number of iterations is exceeded. each cycle of this loop amounts to one
1882 // restart
1883 do
1884 {
1885 VectorType &v = basis_vectors(0, x);
1886
1887 // Compute the preconditioned/unpreconditioned residual for left/right
1888 // preconditioning. If 'x' is the zero vector, then we can bypass the
1889 // full computation. But 'x' is only likely to be the zero vector if
1890 // that's what the user provided as the starting guess, so it's only
1891 // worth checking for this in the first iteration. (Calling all_zero()
1892 // costs as much in memory transfer and communication as computing the
1893 // norm of a vector.)
1894 if (left_precondition)
1895 {
1896 if (accumulated_iterations == 0 && x.all_zero())
1897 preconditioner.vmult(v, b);
1898 else
1899 {
1900 A.vmult(p, x);
1901 p.sadd(-1., 1., b);
1902 preconditioner.vmult(v, p);
1903 }
1904 }
1905 else
1906 {
1907 if (accumulated_iterations == 0 && x.all_zero())
1908 v = b;
1909 else
1910 {
1911 A.vmult(v, x);
1912 v.sadd(-1., 1., b);
1913 }
1914 }
1915
1916 const double norm_v = arnoldi_process.orthonormalize_nth_vector(
1917 0, basis_vectors, accumulated_iterations, re_orthogonalize_signal);
1918
1919 // check the residual here as well since it may be that we got the exact
1920 // (or an almost exact) solution vector at the outset. if we wouldn't
1921 // check here, the next scaling operation would produce garbage
1922 if (use_default_residual)
1923 {
1924 res = norm_v;
1925 if (additional_data.batched_mode)
1926 iteration_state = solver_control.check(accumulated_iterations, res);
1927 else
1928 iteration_state =
1929 this->iteration_status(accumulated_iterations, res, x);
1930
1931 if (iteration_state != SolverControl::iterate)
1932 break;
1933 }
1934 else
1935 {
1936 deallog << "default_res=" << norm_v << std::endl;
1937
1938 if (left_precondition)
1939 {
1940 A.vmult(*r, x);
1941 r->sadd(-1., 1., b);
1942 }
1943 else
1944 preconditioner.vmult(*r, v);
1945
1946 res = r->l2_norm();
1947 if (additional_data.batched_mode)
1948 iteration_state = solver_control.check(accumulated_iterations, res);
1949 else
1950 iteration_state =
1951 this->iteration_status(accumulated_iterations, res, x);
1952
1953 if (iteration_state != SolverControl::iterate)
1954 break;
1955 }
1956
1957 // inner iteration doing at most as many steps as the size of the
1958 // Arnoldi basis
1959 unsigned int inner_iteration = 0;
1960 for (; (inner_iteration < basis_size &&
1961 iteration_state == SolverControl::iterate);
1962 ++inner_iteration)
1963 {
1964 ++accumulated_iterations;
1965 // yet another alias
1966 VectorType &vv = basis_vectors(inner_iteration + 1, x);
1967
1968 if (left_precondition)
1969 {
1970 A.vmult(p, basis_vectors[inner_iteration]);
1971 preconditioner.vmult(vv, p);
1972 }
1973 else
1974 {
1975 preconditioner.vmult(p, basis_vectors[inner_iteration]);
1976 A.vmult(vv, p);
1977 }
1978
1979 res =
1980 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1,
1981 basis_vectors,
1982 accumulated_iterations,
1983 re_orthogonalize_signal);
1984
1985 if (use_default_residual)
1986 {
1987 if (additional_data.batched_mode)
1988 iteration_state =
1989 solver_control.check(accumulated_iterations, res);
1990 else
1991 iteration_state =
1992 this->iteration_status(accumulated_iterations, res, x);
1993 }
1994 else
1995 {
1996 if (!additional_data.batched_mode)
1997 deallog << "default_res=" << res << std::endl;
1998
1999 *x_ = x;
2000 const Vector<double> &projected_solution =
2001 arnoldi_process.solve_projected_system(false);
2002
2003 if (left_precondition)
2004 for (unsigned int i = 0; i < inner_iteration + 1; ++i)
2005 x_->add(projected_solution(i), basis_vectors[i]);
2006 else
2007 {
2008 p = 0.;
2009 for (unsigned int i = 0; i < inner_iteration + 1; ++i)
2010 p.add(projected_solution(i), basis_vectors[i]);
2011 preconditioner.vmult(*r, p);
2012 x_->add(1., *r);
2013 };
2014 A.vmult(*r, *x_);
2015 r->sadd(-1., 1., b);
2016
2017 // Now *r contains the unpreconditioned residual!!
2018 if (left_precondition)
2019 {
2020 res = r->l2_norm();
2021 iteration_state =
2022 this->iteration_status(accumulated_iterations, res, x);
2023 }
2024 else
2025 {
2026 preconditioner.vmult(*x_, *r);
2027 res = x_->l2_norm();
2028
2029 if (additional_data.batched_mode)
2030 iteration_state =
2031 solver_control.check(accumulated_iterations, res);
2032 else
2033 iteration_state =
2034 this->iteration_status(accumulated_iterations, res, x);
2035 }
2036 }
2037 }
2038
2039 // end of inner iteration; now update the global solution vector x with
2040 // the solution of the projected system (least-squares solution)
2041 const Vector<double> &projected_solution =
2042 arnoldi_process.solve_projected_system(true);
2043
2044 if (do_eigenvalues)
2045 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
2046 inner_iteration,
2047 all_eigenvalues_signal,
2048 all_hessenberg_signal,
2049 condition_number_signal);
2050
2051 if (left_precondition)
2052 ::internal::SolverGMRESImplementation::add(
2053 x,
2054 inner_iteration,
2055 projected_solution,
2056 basis_vectors,
2057 false,
2058 arnoldi_process.vector_ptrs);
2059 else
2060 {
2061 ::internal::SolverGMRESImplementation::add(
2062 p,
2063 inner_iteration,
2064 projected_solution,
2065 basis_vectors,
2066 true,
2067 arnoldi_process.vector_ptrs);
2068 preconditioner.vmult(v, p);
2069 x.add(1., v);
2070 }
2071
2072 // in the last round, print the eigenvalues from the last Arnoldi step
2073 if (iteration_state != SolverControl::iterate)
2074 {
2075 if (do_eigenvalues)
2076 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
2077 inner_iteration,
2078 eigenvalues_signal,
2079 hessenberg_signal,
2080 condition_number_signal);
2081
2082 if (!additional_data.batched_mode && !krylov_space_signal.empty())
2083 krylov_space_signal(basis_vectors);
2084
2085 // end of outer iteration. restart if no convergence and the number of
2086 // iterations is not exceeded
2087 }
2088 }
2089 while (iteration_state == SolverControl::iterate);
2090
2091 // in case of failure: throw exception
2092 AssertThrow(iteration_state == SolverControl::success,
2093 SolverControl::NoConvergence(accumulated_iterations, res));
2094}
2095
2096
2097
2098template <typename VectorType>
2100boost::signals2::connection
2102 const std::function<void(double)> &slot,
2103 const bool every_iteration)
2104{
2105 if (every_iteration)
2106 {
2107 return all_condition_numbers_signal.connect(slot);
2108 }
2109 else
2110 {
2111 return condition_number_signal.connect(slot);
2112 }
2113}
2114
2115
2116
2117template <typename VectorType>
2119boost::signals2::connection SolverGMRES<VectorType>::connect_eigenvalues_slot(
2120 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
2121 const bool every_iteration)
2122{
2123 if (every_iteration)
2124 {
2125 return all_eigenvalues_signal.connect(slot);
2126 }
2127 else
2128 {
2129 return eigenvalues_signal.connect(slot);
2130 }
2131}
2132
2133
2134
2135template <typename VectorType>
2137boost::signals2::connection SolverGMRES<VectorType>::connect_hessenberg_slot(
2138 const std::function<void(const FullMatrix<double> &)> &slot,
2139 const bool every_iteration)
2140{
2141 if (every_iteration)
2142 {
2143 return all_hessenberg_signal.connect(slot);
2144 }
2145 else
2146 {
2147 return hessenberg_signal.connect(slot);
2148 }
2149}
2150
2151
2152
2153template <typename VectorType>
2155boost::signals2::connection SolverGMRES<VectorType>::connect_krylov_space_slot(
2156 const std::function<void(
2158{
2159 return krylov_space_signal.connect(slot);
2160}
2161
2162
2163
2164template <typename VectorType>
2166boost::signals2::connection
2168 const std::function<void(int)> &slot)
2169{
2170 return re_orthogonalize_signal.connect(slot);
2171}
2172
2173
2174
2175template <typename VectorType>
2178{
2179 // dummy implementation. this function is not needed for the present
2180 // implementation of gmres
2182 return 0;
2183}
2184
2185
2186
2187//----------------------------------------------------------------------//
2188
2189
2190
2191template <typename VectorType>
2195 const AdditionalData &data)
2196 : SolverBase<VectorType>(cn, mem)
2197 , additional_data(data)
2198{}
2199
2200
2201
2202template <typename VectorType>
2205 const AdditionalData &data)
2206 : SolverBase<VectorType>(cn)
2207 , additional_data(data)
2208{}
2209
2210
2211
2212template <typename VectorType>
2214template <typename MatrixType, typename... PreconditionerTypes>
2218void SolverMPGMRES<VectorType>::solve(
2219 const MatrixType &A,
2220 VectorType &x,
2221 const VectorType &b,
2222 const PreconditionerTypes &...preconditioners)
2223{
2224 LogStream::Prefix prefix("MPGMRES");
2225
2226 if (additional_data.use_truncated_mpgmres_strategy)
2228 A, x, b, IndexingStrategy::truncated_mpgmres, preconditioners...);
2229 else
2231 A, x, b, IndexingStrategy::full_mpgmres, preconditioners...);
2232}
2233
2234
2235
2236template <typename VectorType>
2238template <typename MatrixType, typename... PreconditionerTypes>
2240 const MatrixType &A,
2241 VectorType &x,
2242 const VectorType &b,
2243 const IndexingStrategy &indexing_strategy,
2244 const PreconditionerTypes &...preconditioners)
2245{
2246 constexpr std::size_t n_preconditioners = sizeof...(PreconditionerTypes);
2247
2248 // A lambda for applying the nth preconditioner to a vector src storing
2249 // the result in dst:
2250
2251 const auto apply_nth_preconditioner = [&](unsigned int n,
2252 auto &dst,
2253 const auto &src) {
2254 // We cycle through all preconditioners and call the nth one:
2255 std::size_t i = 0;
2256
2257 [[maybe_unused]] bool preconditioner_called = false;
2258
2259 const auto call_matching_preconditioner = [&](const auto &preconditioner) {
2260 if (i++ == n)
2261 {
2262 Assert(!preconditioner_called, ::ExcInternalError());
2263 preconditioner_called = true;
2264 preconditioner.vmult(dst, src);
2265 }
2266 };
2267
2268 // https://en.cppreference.com/w/cpp/language/fold
2269 (call_matching_preconditioner(preconditioners), ...);
2270 Assert(preconditioner_called, ::ExcInternalError());
2271 };
2272
2273 std::size_t current_index = 0;
2274
2275 // A lambda that cycles through all preconditioners in sequence while
2276 // applying exactly one preconditioner with each function invocation to
2277 // the vector src and storing the result in dst:
2278
2279 const auto preconditioner_vmult = [&](auto &dst, const auto &src) {
2280 // We have no preconditioner that we could apply
2281 if (n_preconditioners == 0)
2282 dst = src;
2283 else
2284 {
2285 apply_nth_preconditioner(current_index, dst, src);
2286 current_index = (current_index + 1) % n_preconditioners;
2287 }
2288 };
2289
2290 // Return the correct index for constructing the next vector in the
2291 // Krylov space sequence according to the chosen indexing strategy
2292
2293 const auto previous_vector_index =
2294 [n_preconditioners, indexing_strategy](unsigned int i) -> unsigned int {
2295 // In the special case of no preconditioners we simply fall back to the
2296 // FGMRES indexing strategy.
2297 if (n_preconditioners == 0)
2298 {
2299 return i;
2300 }
2301
2302 switch (indexing_strategy)
2303 {
2304 case IndexingStrategy::fgmres:
2305 // 0, 1, 2, 3, ...
2306 return i;
2307 case IndexingStrategy::full_mpgmres:
2308 // 0, 0, ..., 1, 1, ..., 2, 2, ..., 3, 3, ...
2309 return i / n_preconditioners;
2310 case IndexingStrategy::truncated_mpgmres:
2311 // 0, 0, ..., 1, 2, 3, ...
2312 return (1 + i >= n_preconditioners) ? (1 + i - n_preconditioners) : 0;
2313 default:
2315 return 0;
2316 }
2317 };
2318
2320
2321 const unsigned int basis_size = additional_data.max_basis_size;
2322
2323 // Generate an object where basis vectors are stored.
2325 basis_size + 1, this->memory);
2327 basis_size, this->memory);
2328
2329 // number of the present iteration; this number is not reset to zero upon a
2330 // restart
2331 unsigned int accumulated_iterations = 0;
2332
2333 // matrix used for the orthogonalization process later
2334 arnoldi_process.initialize(additional_data.orthogonalization_strategy,
2335 basis_size,
2336 false);
2337
2338 // Iteration starts here
2339 double res = std::numeric_limits<double>::lowest();
2340
2341 do
2342 {
2343 // Compute the residual. If 'x' is the zero vector, then we can bypass
2344 // the full computation. But 'x' is only likely to be the zero vector if
2345 // that's what the user provided as the starting guess, so it's only
2346 // worth checking for this in the first iteration. (Calling all_zero()
2347 // costs as much in memory transfer and communication as computing the
2348 // norm of a vector.)
2349 if (accumulated_iterations == 0 && x.all_zero())
2350 v(0, x) = b;
2351 else
2352 {
2353 A.vmult(v(0, x), x);
2354 v[0].sadd(-1., 1., b);
2355 }
2356
2357 res = arnoldi_process.orthonormalize_nth_vector(0, v);
2358 iteration_state = this->iteration_status(accumulated_iterations, res, x);
2359 if (iteration_state == SolverControl::success)
2360 break;
2361
2362 unsigned int inner_iteration = 0;
2363 for (; (inner_iteration < basis_size &&
2364 iteration_state == SolverControl::iterate);
2365 ++inner_iteration)
2366 {
2367 preconditioner_vmult(z(inner_iteration, x),
2368 v[previous_vector_index(inner_iteration)]);
2369 A.vmult(v(inner_iteration + 1, x), z[inner_iteration]);
2370
2371 res =
2372 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1, v);
2373
2374 // check convergence. note that the vector 'x' we pass to the
2375 // criterion is not the final solution we compute if we
2376 // decide to jump out of the iteration (we update 'x' again
2377 // right after the current loop)
2378 iteration_state =
2379 this->iteration_status(++accumulated_iterations, res, x);
2380 }
2381
2382 // Solve triangular system with projected quantities and update solution
2383 // vector
2384 const Vector<double> &projected_solution =
2385 arnoldi_process.solve_projected_system(true);
2386 ::internal::SolverGMRESImplementation::add(
2387 x,
2388 inner_iteration,
2389 projected_solution,
2390 z,
2391 false,
2392 arnoldi_process.vector_ptrs);
2393 }
2394 while (iteration_state == SolverControl::iterate);
2395
2396 // in case of failure: throw exception
2397 if (iteration_state != SolverControl::success)
2398 AssertThrow(false,
2399 SolverControl::NoConvergence(accumulated_iterations, res));
2400}
2401
2402
2403
2404//----------------------------------------------------------------------//
2405
2406
2407
2408template <typename VectorType>
2412 const AdditionalData &data)
2414 cn,
2415 mem,
2416 typename SolverMPGMRES<VectorType>::AdditionalData{
2417 data.max_basis_size,
2418 data.orthogonalization_strategy,
2419 true})
2420{}
2421
2422
2423
2424template <typename VectorType>
2427 const AdditionalData &data)
2429 cn,
2430 typename SolverMPGMRES<VectorType>::AdditionalData{
2431 data.max_basis_size,
2432 data.orthogonalization_strategy,
2433 true})
2434{}
2435
2436
2437
2438template <typename VectorType>
2440template <typename MatrixType, typename... PreconditionerTypes>
2444void SolverFGMRES<VectorType>::solve(
2445 const MatrixType &A,
2446 VectorType &x,
2447 const VectorType &b,
2448 const PreconditionerTypes &...preconditioners)
2449{
2450 LogStream::Prefix prefix("FGMRES");
2452 A, x, b, SolverFGMRES::IndexingStrategy::fgmres, preconditioners...);
2453}
2454
2455
2456#endif // DOXYGEN
2457
2459
2460#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())
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)
SolverMPGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
internal::SolverGMRESImplementation::ArnoldiProcess< typename VectorType::value_type > arnoldi_process
SolverMPGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
void solve_internal(const MatrixType &A, VectorType &x, const VectorType &b, const IndexingStrategy &indexing_strategy, const PreconditionerTypes &...preconditioners)
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
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:35
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:243
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
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)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
LogStream deallog
Definition logstream.cc:38
std::vector< index_type > data
Definition mpi.cc:750
types::global_dof_index locally_owned_size
Definition mpi.cc:837
@ matrix
Contents is actually a matrix.
constexpr char A
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
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
AdditionalData(const unsigned int max_basis_size=30, const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy=LinearAlgebra::OrthogonalizationStrategy::delayed_classical_gram_schmidt, const bool use_truncated_mpgmres_strategy=true)
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)