Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10: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 <algorithm>
35#include <cmath>
36#include <limits>
37#include <memory>
38#include <vector>
39
41
42// forward declarations
43#ifndef DOXYGEN
44namespace LinearAlgebra
45{
46 namespace distributed
47 {
48 template <typename, typename>
49 class Vector;
50 } // namespace distributed
51} // namespace LinearAlgebra
52#endif
53
59namespace internal
60{
64 namespace SolverGMRESImplementation
65 {
73 template <typename VectorType>
75 {
76 public:
81 TmpVectors(const unsigned int max_size, VectorMemory<VectorType> &vmem);
82
86 ~TmpVectors() = default;
87
92 VectorType &
93 operator[](const unsigned int i) const;
94
101 VectorType &
102 operator()(const unsigned int i, const VectorType &temp);
103
108 unsigned int
109 size() const;
110
111
112 private:
117
121 std::vector<typename VectorMemory<VectorType>::Pointer> data;
122 };
123
124
125
135 {
136 public:
141 void
144 const unsigned int max_basis_size,
145 const bool force_reorthogonalization);
146
169 template <typename VectorType>
170 double
172 const unsigned int n,
173 TmpVectors<VectorType> &orthogonal_vectors,
174 const unsigned int accumulated_iterations = 0,
175 const boost::signals2::signal<void(int)> &reorthogonalize_signal =
176 boost::signals2::signal<void(int)>());
177
185 const Vector<double> &
186 solve_projected_system(const bool orthogonalization_finished);
187
192 const FullMatrix<double> &
194
195 private:
200
207
212 std::vector<std::pair<double, double>> givens_rotations;
213
218
224
229
239
244
272 double
273 do_givens_rotation(const bool delayed_reorthogonalization,
274 const int col,
275 FullMatrix<double> &matrix,
276 std::vector<std::pair<double, double>> &rotations,
277 Vector<double> &rhs);
278 };
279 } // namespace SolverGMRESImplementation
280} // namespace internal
281
349template <typename VectorType = Vector<double>>
350class SolverGMRES : public SolverBase<VectorType>
351{
352public:
357 {
368 explicit AdditionalData(const unsigned int max_basis_size = 30,
369 const bool right_preconditioning = false,
370 const bool use_default_residual = true,
371 const bool force_re_orthogonalization = false,
372 const bool batched_mode = false,
376 delayed_classical_gram_schmidt);
377
389 unsigned int max_n_tmp_vectors;
390
398 unsigned int max_basis_size;
399
408
413
421
429
434 };
435
441 const AdditionalData &data = AdditionalData());
442
448
453
457 template <typename MatrixType, typename PreconditionerType>
458 void
459 solve(const MatrixType &A,
460 VectorType &x,
461 const VectorType &b,
462 const PreconditionerType &preconditioner);
463
470 boost::signals2::connection
471 connect_condition_number_slot(const std::function<void(double)> &slot,
472 const bool every_iteration = false);
473
480 boost::signals2::connection
482 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
483 const bool every_iteration = false);
484
492 boost::signals2::connection
494 const std::function<void(const FullMatrix<double> &)> &slot,
495 const bool every_iteration = true);
496
503 boost::signals2::connection
505 const std::function<
507 &slot);
508
509
514 boost::signals2::connection
515 connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
516
517
519 int,
520 << "The number of temporary vectors you gave (" << arg1
521 << ") is too small. It should be at least 10 for "
522 << "any results, and much more for reasonable ones.");
523
524protected:
529
534 boost::signals2::signal<void(double)> condition_number_signal;
535
540 boost::signals2::signal<void(double)> all_condition_numbers_signal;
541
546 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
548
553 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
555
560 boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
561
566 boost::signals2::signal<void(const FullMatrix<double> &)>
568
573 boost::signals2::signal<void(
576
581 boost::signals2::signal<void(int)> re_orthogonalize_signal;
582
589
593 virtual double
595
602 static void
604 const FullMatrix<double> &H_orig,
605 const unsigned int n,
606 const boost::signals2::signal<
607 void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
608 const boost::signals2::signal<void(const FullMatrix<double> &)>
610 const boost::signals2::signal<void(double)> &cond_signal);
611
617};
618
619
620
641template <typename VectorType = Vector<double>>
642class SolverFGMRES : public SolverBase<VectorType>
643{
644public:
672
678 const AdditionalData &data = AdditionalData());
679
685 const AdditionalData &data = AdditionalData());
686
690 template <typename MatrixType, typename PreconditionerType>
691 void
692 solve(const MatrixType &A,
693 VectorType &x,
694 const VectorType &b,
695 const PreconditionerType &preconditioner);
696
697private:
702
708};
709
711/* --------------------- Inline and template functions ------------------- */
712
713
714#ifndef DOXYGEN
715
716template <typename VectorType>
718 const unsigned int max_basis_size,
719 const bool right_preconditioning,
720 const bool use_default_residual,
721 const bool force_re_orthogonalization,
722 const bool batched_mode,
723 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy)
724 : max_n_tmp_vectors(0)
725 , max_basis_size(max_basis_size)
726 , right_preconditioning(right_preconditioning)
727 , use_default_residual(use_default_residual)
728 , force_re_orthogonalization(force_re_orthogonalization)
729 , batched_mode(batched_mode)
730 , orthogonalization_strategy(orthogonalization_strategy)
731{
732 Assert(max_basis_size >= 1,
733 ExcMessage("SolverGMRES needs at least one vector in the "
734 "Arnoldi basis."));
735}
736
737
738
739template <typename VectorType>
742 const AdditionalData &data)
743 : SolverBase<VectorType>(cn, mem)
744 , additional_data(data)
745 , solver_control(cn)
746{}
747
748
749
750template <typename VectorType>
752 const AdditionalData &data)
753 : SolverBase<VectorType>(cn)
754 , additional_data(data)
755 , solver_control(cn)
756{}
757
758
759
760namespace internal
761{
762 namespace SolverGMRESImplementation
763 {
764 template <typename VectorType>
765 inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
767 : mem(vmem)
768 , data(max_size)
769 {}
770
771
772
773 template <typename VectorType>
774 inline VectorType &
775 TmpVectors<VectorType>::operator[](const unsigned int i) const
776 {
777 AssertIndexRange(i, data.size());
778
779 Assert(data[i] != nullptr, ExcNotInitialized());
780 return *data[i];
781 }
782
783
784
785 template <typename VectorType>
786 inline VectorType &
787 TmpVectors<VectorType>::operator()(const unsigned int i,
788 const VectorType &temp)
789 {
790 AssertIndexRange(i, data.size());
791 if (data[i] == nullptr)
792 {
793 data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
794 data[i]->reinit(temp, true);
795 }
796 return *data[i];
797 }
798
799
800
801 template <typename VectorType>
802 unsigned int
803 TmpVectors<VectorType>::size() const
804 {
805 return (data.size() > 0 ? data.size() - 1 : 0);
806 }
807
808
809
810 template <typename VectorType, typename Enable = void>
811 struct is_dealii_compatible_distributed_vector;
812
813 template <typename VectorType>
814 struct is_dealii_compatible_distributed_vector<
815 VectorType,
816 std::enable_if_t<!internal::is_block_vector<VectorType>>>
817 {
818 static constexpr bool value = std::is_same_v<
819 VectorType,
820 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
822 };
823
824
825
826 template <typename VectorType>
827 struct is_dealii_compatible_distributed_vector<
828 VectorType,
829 std::enable_if_t<internal::is_block_vector<VectorType>>>
830 {
831 static constexpr bool value = std::is_same_v<
832 typename VectorType::BlockType,
833 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
835 };
836
837
838
839 template <typename VectorType,
840 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
841 * = nullptr>
842 unsigned int
843 n_blocks(const VectorType &)
844 {
845 return 1;
846 }
847
848
849
850 template <typename VectorType,
851 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
852 nullptr>
853 unsigned int
854 n_blocks(const VectorType &vector)
855 {
856 return vector.n_blocks();
857 }
858
859
860
861 template <typename VectorType,
862 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
863 * = nullptr>
864 VectorType &
865 block(VectorType &vector, const unsigned int b)
866 {
867 AssertDimension(b, 0);
868 (void)b;
869 return vector;
870 }
871
872
873
874 template <typename VectorType,
875 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
876 * = nullptr>
877 const VectorType &
878 block(const VectorType &vector, const unsigned int b)
879 {
880 AssertDimension(b, 0);
881 (void)b;
882 return vector;
883 }
884
885
886
887 template <typename VectorType,
888 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
889 nullptr>
890 typename VectorType::BlockType &
891 block(VectorType &vector, const unsigned int b)
892 {
893 return vector.block(b);
894 }
895
896
897
898 template <typename VectorType,
899 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
900 nullptr>
901 const typename VectorType::BlockType &
902 block(const VectorType &vector, const unsigned int b)
903 {
904 return vector.block(b);
905 }
906
907
908
909 template <bool delayed_reorthogonalization,
910 typename VectorType,
911 std::enable_if_t<
912 !is_dealii_compatible_distributed_vector<VectorType>::value,
913 VectorType> * = nullptr>
914 void
915 Tvmult_add(const unsigned int n,
916 const VectorType &vv,
917 const TmpVectors<VectorType> &orthogonal_vectors,
919 {
920 for (unsigned int i = 0; i < n; ++i)
921 {
922 h(i) += vv * orthogonal_vectors[i];
923 if (delayed_reorthogonalization)
924 h(n + i) += orthogonal_vectors[i] * orthogonal_vectors[n - 1];
925 }
926 if (delayed_reorthogonalization)
927 h(n + n) += vv * vv;
928 }
929
930
931
932 template <bool delayed_reorthogonalization,
933 typename VectorType,
934 std::enable_if_t<
935 is_dealii_compatible_distributed_vector<VectorType>::value,
936 VectorType> * = nullptr>
937 void
938 Tvmult_add(const unsigned int n,
939 const VectorType &vv,
940 const TmpVectors<VectorType> &orthogonal_vectors,
942 {
943 for (unsigned int b = 0; b < n_blocks(vv); ++b)
944 {
945 unsigned int j = 0;
946
947 if (n <= 128)
948 {
949 // optimized path
950 static constexpr unsigned int n_lanes =
952
954 for (unsigned int i = 0; i < n; ++i)
955 hs[i] = 0.0;
957 correct[delayed_reorthogonalization ? 129 : 1];
958 if (delayed_reorthogonalization)
959 for (unsigned int i = 0; i < n + 1; ++i)
960 correct[i] = 0.0;
961
962 unsigned int c = 0;
963
964 constexpr unsigned int inner_batch_size =
965 delayed_reorthogonalization ? 6 : 12;
966
967 for (; c < block(vv, b).locally_owned_size() / n_lanes /
968 inner_batch_size;
969 ++c, j += n_lanes * inner_batch_size)
970 {
971 VectorizedArray<double> vvec[inner_batch_size];
972 for (unsigned int k = 0; k < inner_batch_size; ++k)
973 vvec[k].load(block(vv, b).begin() + j + k * n_lanes);
974 VectorizedArray<double> last_vector[inner_batch_size];
975 for (unsigned int k = 0; k < inner_batch_size; ++k)
976 last_vector[k].load(
977 block(orthogonal_vectors[n - 1], b).begin() + j +
978 k * n_lanes);
979
980 {
981 VectorizedArray<double> local_sum_0 =
982 last_vector[0] * vvec[0];
983 VectorizedArray<double> local_sum_1 =
984 last_vector[0] * last_vector[0];
985 VectorizedArray<double> local_sum_2 = vvec[0] * vvec[0];
986 for (unsigned int k = 1; k < inner_batch_size; ++k)
987 {
988 local_sum_0 += last_vector[k] * vvec[k];
989 if (delayed_reorthogonalization)
990 {
991 local_sum_1 += last_vector[k] * last_vector[k];
992 local_sum_2 += vvec[k] * vvec[k];
993 }
994 }
995 hs[n - 1] += local_sum_0;
996 if (delayed_reorthogonalization)
997 {
998 correct[n - 1] += local_sum_1;
999 correct[n] += local_sum_2;
1000 }
1001 }
1002
1003 for (unsigned int i = 0; i < n - 1; ++i)
1004 {
1005 // break the dependency chain into the field hs[i] for
1006 // small sizes i by first accumulating 4 or 8 results
1007 // into a local variable
1009 temp.load(block(orthogonal_vectors[i], b).begin() + j);
1010 VectorizedArray<double> local_sum_0 = temp * vvec[0];
1011 VectorizedArray<double> local_sum_1 =
1012 delayed_reorthogonalization ? temp * last_vector[0] :
1013 0.;
1014 for (unsigned int k = 1; k < inner_batch_size; ++k)
1015 {
1016 temp.load(block(orthogonal_vectors[i], b).begin() +
1017 j + k * n_lanes);
1018 local_sum_0 += temp * vvec[k];
1019 if (delayed_reorthogonalization)
1020 local_sum_1 += temp * last_vector[k];
1021 }
1022 hs[i] += local_sum_0;
1023 if (delayed_reorthogonalization)
1024 correct[i] += local_sum_1;
1025 }
1026 }
1027
1028 c *= inner_batch_size;
1029 for (; c < block(vv, b).locally_owned_size() / n_lanes;
1030 ++c, j += n_lanes)
1031 {
1032 VectorizedArray<double> vvec, last_vector;
1033 vvec.load(block(vv, b).begin() + j);
1034 last_vector.load(block(orthogonal_vectors[n - 1], b).begin() +
1035 j);
1036 hs[n - 1] += last_vector * vvec;
1037 if (delayed_reorthogonalization)
1038 {
1039 correct[n - 1] += last_vector * last_vector;
1040 correct[n] += vvec * vvec;
1041 }
1042
1043 for (unsigned int i = 0; i < n - 1; ++i)
1044 {
1046 temp.load(block(orthogonal_vectors[i], b).begin() + j);
1047 hs[i] += temp * vvec;
1048 if (delayed_reorthogonalization)
1049 correct[i] += temp * last_vector;
1050 }
1051 }
1052
1053 for (unsigned int i = 0; i < n; ++i)
1054 {
1055 h(i) += hs[i].sum();
1056 if (delayed_reorthogonalization)
1057 h(i + n) += correct[i].sum();
1058 }
1059 if (delayed_reorthogonalization)
1060 h(n + n) += correct[n].sum();
1061 }
1062
1063 // remainder loop of optimized path or non-optimized path (if
1064 // n>128)
1065 for (; j < block(vv, b).locally_owned_size(); ++j)
1066 {
1067 const double vvec = block(vv, b).local_element(j);
1068 const double last_vector =
1069 block(orthogonal_vectors[n - 1], b).local_element(j);
1070 h(n - 1) += last_vector * vvec;
1071 if (delayed_reorthogonalization)
1072 {
1073 h(n + n - 1) += last_vector * last_vector;
1074 h(n + n) += vvec * vvec;
1075 }
1076 for (unsigned int i = 0; i < n - 1; ++i)
1077 {
1078 const double temp =
1079 block(orthogonal_vectors[i], b).local_element(j);
1080 h(i) += temp * vvec;
1081 if (delayed_reorthogonalization)
1082 h(n + i) += temp * last_vector;
1083 }
1084 }
1085 }
1086
1087 Utilities::MPI::sum(h, block(vv, 0).get_mpi_communicator(), h);
1088 }
1089
1090
1091
1092 template <bool delayed_reorthogonalization,
1093 typename VectorType,
1094 std::enable_if_t<
1095 !is_dealii_compatible_distributed_vector<VectorType>::value,
1096 VectorType> * = nullptr>
1097 double
1098 subtract_and_norm(const unsigned int n,
1099 const TmpVectors<VectorType> &orthogonal_vectors,
1100 const Vector<double> &h,
1101 VectorType &vv)
1102 {
1103 Assert(n > 0, ExcInternalError());
1104
1105 VectorType &last_vector =
1106 const_cast<VectorType &>(orthogonal_vectors[n - 1]);
1107 for (unsigned int i = 0; i < n - 1; ++i)
1108 {
1109 if (delayed_reorthogonalization && i + 2 < n)
1110 last_vector.add(-h(n + i), orthogonal_vectors[i]);
1111 vv.add(-h(i), orthogonal_vectors[i]);
1112 }
1113
1114 if (delayed_reorthogonalization)
1115 {
1116 if (n > 1)
1117 last_vector.sadd(1. / h(n + n - 1),
1118 -h(n + n - 2) / h(n + n - 1),
1119 orthogonal_vectors[n - 2]);
1120
1121 // h(n + n) is lucky breakdown
1122 const double scaling_factor_vv = h(n + n) > 0.0 ?
1123 1. / (h(n + n - 1) * h(n + n)) :
1124 1. / (h(n + n - 1) * h(n + n - 1));
1125 vv.sadd(scaling_factor_vv,
1126 -h(n - 1) * scaling_factor_vv,
1127 last_vector);
1128
1129 // the delayed reorthogonalization computes the norm from other
1130 // quantities
1131 return std::numeric_limits<double>::signaling_NaN();
1132 }
1133 else
1134 return std::sqrt(
1135 vv.add_and_dot(-h(n - 1), orthogonal_vectors[n - 1], vv));
1136 }
1137
1138
1139
1140 template <bool delayed_reorthogonalization,
1141 typename VectorType,
1142 std::enable_if_t<
1143 is_dealii_compatible_distributed_vector<VectorType>::value,
1144 VectorType> * = nullptr>
1145 double
1146 subtract_and_norm(const unsigned int n,
1147 const TmpVectors<VectorType> &orthogonal_vectors,
1148 const Vector<double> &h,
1149 VectorType &vv)
1150 {
1151 static constexpr unsigned int n_lanes = VectorizedArray<double>::size();
1152
1153 double norm_vv_temp = 0.0;
1154 VectorType &last_vector =
1155 const_cast<VectorType &>(orthogonal_vectors[n - 1]);
1156 const double inverse_norm_previous =
1157 delayed_reorthogonalization ? 1. / h(n + n - 1) : 0.;
1158 const double scaling_factor_vv =
1159 delayed_reorthogonalization ?
1160 (h(n + n) > 0.0 ? inverse_norm_previous / h(n + n) :
1161 inverse_norm_previous / h(n + n - 1)) :
1162 0.;
1163
1164 for (unsigned int b = 0; b < n_blocks(vv); ++b)
1165 {
1166 VectorizedArray<double> norm_vv_temp_vectorized = 0.0;
1167
1168 constexpr unsigned int inner_batch_size =
1169 delayed_reorthogonalization ? 6 : 12;
1170
1171 unsigned int j = 0;
1172 unsigned int c = 0;
1173 for (; c <
1174 block(vv, b).locally_owned_size() / n_lanes / inner_batch_size;
1175 ++c, j += n_lanes * inner_batch_size)
1176 {
1177 VectorizedArray<double> temp[inner_batch_size];
1178 VectorizedArray<double> last_vec[inner_batch_size];
1179
1180 const double last_factor = h(n - 1);
1181 for (unsigned int k = 0; k < inner_batch_size; ++k)
1182 {
1183 temp[k].load(block(vv, b).begin() + j + k * n_lanes);
1184 last_vec[k].load(block(last_vector, b).begin() + j +
1185 k * n_lanes);
1186 if (!delayed_reorthogonalization)
1187 temp[k] -= last_factor * last_vec[k];
1188 }
1189
1190 for (unsigned int i = 0; i < n - 1; ++i)
1191 {
1192 const double factor = h(i);
1193 const double correction_factor =
1194 (delayed_reorthogonalization ? h(n + i) : 0.0);
1195 for (unsigned int k = 0; k < inner_batch_size; ++k)
1196 {
1198 vec.load(block(orthogonal_vectors[i], b).begin() + j +
1199 k * n_lanes);
1200 temp[k] -= factor * vec;
1201 if (delayed_reorthogonalization)
1202 last_vec[k] -= correction_factor * vec;
1203 }
1204 }
1205
1206 if (delayed_reorthogonalization)
1207 for (unsigned int k = 0; k < inner_batch_size; ++k)
1208 {
1209 last_vec[k] = last_vec[k] * inverse_norm_previous;
1210 last_vec[k].store(block(last_vector, b).begin() + j +
1211 k * n_lanes);
1212 temp[k] -= last_factor * last_vec[k];
1213 temp[k] = temp[k] * scaling_factor_vv;
1214 temp[k].store(block(vv, b).begin() + j + k * n_lanes);
1215 }
1216 else
1217 for (unsigned int k = 0; k < inner_batch_size; ++k)
1218 {
1219 temp[k].store(block(vv, b).begin() + j + k * n_lanes);
1220 norm_vv_temp_vectorized += temp[k] * temp[k];
1221 }
1222 }
1223
1224 c *= inner_batch_size;
1225 for (; c < block(vv, b).locally_owned_size() / n_lanes;
1226 ++c, j += n_lanes)
1227 {
1228 VectorizedArray<double> temp, last_vec;
1229 temp.load(block(vv, b).begin() + j);
1230 last_vec.load(block(last_vector, b).begin() + j);
1231 if (!delayed_reorthogonalization)
1232 temp -= h(n - 1) * last_vec;
1233
1234 for (unsigned int i = 0; i < n - 1; ++i)
1235 {
1237 vec.load(block(orthogonal_vectors[i], b).begin() + j);
1238 temp -= h(i) * vec;
1239 if (delayed_reorthogonalization)
1240 last_vec -= h(n + i) * vec;
1241 }
1242
1243 if (delayed_reorthogonalization)
1244 {
1245 last_vec = last_vec * inverse_norm_previous;
1246 last_vec.store(block(last_vector, b).begin() + j);
1247 temp -= h(n - 1) * last_vec;
1248 temp = temp * scaling_factor_vv;
1249 temp.store(block(vv, b).begin() + j);
1250 }
1251 else
1252 {
1253 temp.store(block(vv, b).begin() + j);
1254 norm_vv_temp_vectorized += temp * temp;
1255 }
1256 }
1257
1258 if (!delayed_reorthogonalization)
1259 norm_vv_temp += norm_vv_temp_vectorized.sum();
1260
1261 for (; j < block(vv, b).locally_owned_size(); ++j)
1262 {
1263 double temp = block(vv, b).local_element(j);
1264 double last_vec = block(last_vector, b).local_element(j);
1265 if (delayed_reorthogonalization)
1266 {
1267 for (unsigned int i = 0; i < n - 1; ++i)
1268 {
1269 const double vec =
1270 block(orthogonal_vectors[i], b).local_element(j);
1271 temp -= h(i) * vec;
1272 last_vec -= h(n + i) * vec;
1273 }
1274 last_vec *= inverse_norm_previous;
1275 block(last_vector, b).local_element(j) = last_vec;
1276 temp -= h(n - 1) * last_vec;
1277 temp *= scaling_factor_vv;
1278 }
1279 else
1280 {
1281 temp -= h(n - 1) * last_vec;
1282 for (unsigned int i = 0; i < n - 1; ++i)
1283 temp -=
1284 h(i) * block(orthogonal_vectors[i], b).local_element(j);
1285 norm_vv_temp += temp * temp;
1286 }
1287 block(vv, b).local_element(j) = temp;
1288 }
1289 }
1290
1291 return std::sqrt(
1292 Utilities::MPI::sum(norm_vv_temp, block(vv, 0).get_mpi_communicator()));
1293 }
1294
1295
1296
1297 template <typename VectorType,
1298 std::enable_if_t<
1299 !is_dealii_compatible_distributed_vector<VectorType>::value,
1300 VectorType> * = nullptr>
1301 void
1302 add(VectorType &p,
1303 const unsigned int n,
1304 const Vector<double> &h,
1305 const TmpVectors<VectorType> &tmp_vectors,
1306 const bool zero_out)
1307 {
1308 if (zero_out)
1309 p.equ(h(0), tmp_vectors[0]);
1310 else
1311 p.add(h(0), tmp_vectors[0]);
1312
1313 for (unsigned int i = 1; i < n; ++i)
1314 p.add(h(i), tmp_vectors[i]);
1315 }
1316
1317
1318
1319 template <typename VectorType,
1320 std::enable_if_t<
1321 is_dealii_compatible_distributed_vector<VectorType>::value,
1322 VectorType> * = nullptr>
1323 void
1324 add(VectorType &p,
1325 const unsigned int n,
1326 const Vector<double> &h,
1327 const TmpVectors<VectorType> &tmp_vectors,
1328 const bool zero_out)
1329 {
1330 for (unsigned int b = 0; b < n_blocks(p); ++b)
1331 for (unsigned int j = 0; j < block(p, b).locally_owned_size(); ++j)
1332 {
1333 double temp = zero_out ? 0 : block(p, b).local_element(j);
1334 for (unsigned int i = 0; i < n; ++i)
1335 temp += block(tmp_vectors[i], b).local_element(j) * h(i);
1336 block(p, b).local_element(j) = temp;
1337 }
1338 }
1339
1340
1341
1342 inline void
1343 ArnoldiProcess::initialize(
1344 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy,
1345 const unsigned int basis_size,
1346 const bool force_reorthogonalization)
1347 {
1348 this->orthogonalization_strategy = orthogonalization_strategy;
1349 this->do_reorthogonalization = force_reorthogonalization;
1350
1351 hessenberg_matrix.reinit(basis_size + 1, basis_size);
1352 triangular_matrix.reinit(basis_size + 1, basis_size, true);
1353
1354 // some additional vectors, also used in the orthogonalization
1355 projected_rhs.reinit(basis_size + 1, true);
1356 givens_rotations.reserve(basis_size);
1357
1358 if (orthogonalization_strategy ==
1361 h.reinit(2 * basis_size + 3);
1362 else
1363 h.reinit(basis_size + 1);
1364 }
1365
1366
1367
1368 template <typename VectorType>
1369 inline double
1370 ArnoldiProcess::orthonormalize_nth_vector(
1371 const unsigned int n,
1372 TmpVectors<VectorType> &orthogonal_vectors,
1373 const unsigned int accumulated_iterations,
1374 const boost::signals2::signal<void(int)> &reorthogonalize_signal)
1375 {
1376 AssertIndexRange(n, hessenberg_matrix.m());
1377 AssertIndexRange(n, orthogonal_vectors.size() + 1);
1378
1379 VectorType &vv = orthogonal_vectors[n];
1380
1381 double residual_estimate = std::numeric_limits<double>::signaling_NaN();
1382 if (n == 0)
1383 {
1384 givens_rotations.clear();
1385 residual_estimate = vv.l2_norm();
1386 if (residual_estimate != 0.)
1387 vv /= residual_estimate;
1388 projected_rhs(0) = residual_estimate;
1389 }
1390 else if (orthogonalization_strategy ==
1393 {
1394 // The algorithm implemented in the following few lines is algorithm
1395 // 4 of Bielich et al. (2022).
1396
1397 // To avoid un-scaled numbers as appearing with the original
1398 // algorithm of Bielich et al., we use a preliminary scaling of the
1399 // last vector. This will be corrected in the delayed step.
1400 const double previous_scaling = n > 0 ? h(n + n - 2) : 1.;
1401
1402 // Reset h to zero
1403 h.reinit(n + n + 1);
1404
1405 // global reduction
1406 Tvmult_add<true>(n, vv, orthogonal_vectors, h);
1407
1408 // delayed correction terms
1409 double tmp = 0;
1410 for (unsigned int i = 0; i < n - 1; ++i)
1411 tmp += h(n + i) * h(n + i);
1412 const double alpha_j = h(n + n - 1) > tmp ?
1413 std::sqrt(h(n + n - 1) - tmp) :
1414 std::sqrt(h(n + n - 1));
1415 h(n + n - 1) = alpha_j;
1416
1417 tmp = 0;
1418 for (unsigned int i = 0; i < n - 1; ++i)
1419 tmp += h(i) * h(n + i);
1420 h(n - 1) = (h(n - 1) - tmp) / alpha_j;
1421
1422 // representation of H(j-1)
1423 if (n > 1)
1424 {
1425 for (unsigned int i = 0; i < n - 1; ++i)
1426 hessenberg_matrix(i, n - 2) += h(n + i) * previous_scaling;
1427 hessenberg_matrix(n - 1, n - 2) = alpha_j * previous_scaling;
1428 }
1429 for (unsigned int i = 0; i < n; ++i)
1430 {
1431 double sum = 0;
1432 for (unsigned int j = (i == 0 ? 0 : i - 1); j < n - 1; ++j)
1433 sum += hessenberg_matrix(i, j) * h(n + j);
1434 hessenberg_matrix(i, n - 1) = (h(i) - sum) / alpha_j;
1435 }
1436
1437 // Compute estimate norm for approximate convergence criterion (to
1438 // be corrected in next iteration)
1439 double sum = 0;
1440 for (unsigned int i = 0; i < n - 1; ++i)
1441 sum += h(i) * h(i);
1442 sum += (2. - 1.) * h(n - 1) * h(n - 1);
1443 hessenberg_matrix(n, n - 1) =
1444 std::sqrt(std::abs(h(n + n) - sum)) / alpha_j;
1445
1446 // projection and delayed reorthogonalization. We scale the vector
1447 // vv here by the preliminary norm to avoid working with too large
1448 // values and correct to the actual norm in high precision in the
1449 // next iteration.
1450 h(n + n) = hessenberg_matrix(n, n - 1);
1451 subtract_and_norm<true>(n, orthogonal_vectors, h, vv);
1452
1453 // transform new column of upper Hessenberg matrix into upper
1454 // triangular form by computing the respective factor
1455 residual_estimate = do_givens_rotation(
1456 true, n - 2, triangular_matrix, givens_rotations, projected_rhs);
1457 }
1458 else
1459 {
1460 // need initial norm for detection of re-orthogonalization, see below
1461 double norm_vv = 0.0;
1462 double norm_vv_start = 0;
1463 const bool consider_reorthogonalize =
1464 (do_reorthogonalization == false) && (n % 5 == 0);
1465 if (consider_reorthogonalize)
1466 norm_vv_start = vv.l2_norm();
1467
1468 // Reset h to zero
1469 h.reinit(n);
1470
1471 // run two loops with index 0: orthogonalize, 1: reorthogonalize
1472 for (unsigned int c = 0; c < 2; ++c)
1473 {
1474 // Orthogonalization
1475 if (orthogonalization_strategy ==
1478 {
1479 double htmp = vv * orthogonal_vectors[0];
1480 h(0) += htmp;
1481 for (unsigned int i = 1; i < n; ++i)
1482 {
1483 htmp = vv.add_and_dot(-htmp,
1484 orthogonal_vectors[i - 1],
1485 orthogonal_vectors[i]);
1486 h(i) += htmp;
1487 }
1488
1489 norm_vv = std::sqrt(
1490 vv.add_and_dot(-htmp, orthogonal_vectors[n - 1], vv));
1491 }
1492 else if (orthogonalization_strategy ==
1495 {
1496 Tvmult_add<false>(n, vv, orthogonal_vectors, h);
1497 norm_vv =
1498 subtract_and_norm<false>(n, orthogonal_vectors, h, vv);
1499 }
1500 else
1501 {
1503 }
1504
1505 if (c == 1)
1506 break; // reorthogonalization already performed -> finished
1507
1508 // Re-orthogonalization if loss of orthogonality detected. For the
1509 // test, use a strategy discussed in C. T. Kelley, Iterative
1510 // Methods for Linear and Nonlinear Equations, SIAM, Philadelphia,
1511 // 1995: Compare the norm of vv after orthogonalization with its
1512 // norm when starting the orthogonalization. If vv became very
1513 // small (here: less than the square root of the machine precision
1514 // times 10), it is almost in the span of the previous vectors,
1515 // which indicates loss of precision.
1516 if (consider_reorthogonalize)
1517 {
1518 if (norm_vv >
1519 10. * norm_vv_start *
1520 std::sqrt(std::numeric_limits<
1521 typename VectorType::value_type>::epsilon()))
1522 break;
1523
1524 else
1525 {
1526 do_reorthogonalization = true;
1527 if (!reorthogonalize_signal.empty())
1528 reorthogonalize_signal(accumulated_iterations);
1529 }
1530 }
1531
1532 if (do_reorthogonalization == false)
1533 break; // no reorthogonalization needed -> finished
1534 }
1535
1536 for (unsigned int i = 0; i < n; ++i)
1537 hessenberg_matrix(i, n - 1) = h(i);
1538 hessenberg_matrix(n, n - 1) = norm_vv;
1539
1540 // norm_vv is a lucky breakdown, the solver will reach convergence,
1541 // but we must not divide by zero here.
1542 if (norm_vv != 0)
1543 vv /= norm_vv;
1544
1545 residual_estimate = do_givens_rotation(
1546 false, n - 1, triangular_matrix, givens_rotations, projected_rhs);
1547 }
1548
1549 return residual_estimate;
1550 }
1551
1552
1553
1554 inline double
1555 ArnoldiProcess::do_givens_rotation(
1556 const bool delayed_reorthogonalization,
1557 const int col,
1558 FullMatrix<double> &matrix,
1559 std::vector<std::pair<double, double>> &rotations,
1560 Vector<double> &rhs)
1561 {
1562 // for the delayed orthogonalization, we can only compute the column of
1563 // the previous column (as there will be correction terms added to the
1564 // present column for stability reasons), but we still want to compute
1565 // the residual estimate from the accumulated work; we therefore perform
1566 // givens rotations on two columns simultaneously
1567 if (delayed_reorthogonalization)
1568 {
1569 if (col >= 0)
1570 {
1571 AssertDimension(rotations.size(), static_cast<std::size_t>(col));
1572 matrix(0, col) = hessenberg_matrix(0, col);
1573 }
1574 double H_next = hessenberg_matrix(0, col + 1);
1575 for (int i = 0; i < col; ++i)
1576 {
1577 const double c = rotations[i].first;
1578 const double s = rotations[i].second;
1579 const double Hi = matrix(i, col);
1580 const double Hi1 = hessenberg_matrix(i + 1, col);
1581 H_next = -s * H_next + c * hessenberg_matrix(i + 1, col + 1);
1582 matrix(i, col) = c * Hi + s * Hi1;
1583 matrix(i + 1, col) = -s * Hi + c * Hi1;
1584 }
1585
1586 if (col >= 0)
1587 {
1588 const double H_col1 = hessenberg_matrix(col + 1, col);
1589 const double H_col = matrix(col, col);
1590 const double r = 1. / std::sqrt(H_col * H_col + H_col1 * H_col1);
1591 rotations.emplace_back(H_col * r, H_col1 * r);
1592 matrix(col, col) =
1593 rotations[col].first * H_col + rotations[col].second * H_col1;
1594
1595 rhs(col + 1) = -rotations[col].second * rhs(col);
1596 rhs(col) *= rotations[col].first;
1597
1598 H_next =
1599 -rotations[col].second * H_next +
1600 rotations[col].first * hessenberg_matrix(col + 1, col + 1);
1601 }
1602
1603 const double H_last = hessenberg_matrix(col + 2, col + 1);
1604 const double r = 1. / std::sqrt(H_next * H_next + H_last * H_last);
1605 return std::abs(H_last * r * rhs(col + 1));
1606 }
1607 else
1608 {
1609 AssertDimension(rotations.size(), static_cast<std::size_t>(col));
1610
1611 matrix(0, col) = hessenberg_matrix(0, col);
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 matrix(i, col) = c * Hi + s * Hi1;
1619 matrix(i + 1, col) = -s * Hi + c * Hi1;
1620 }
1621
1622 const double Hi = matrix(col, col);
1623 const double Hi1 = hessenberg_matrix(col + 1, col);
1624 const double r = 1. / std::sqrt(Hi * Hi + Hi1 * Hi1);
1625 rotations.emplace_back(Hi * r, Hi1 * r);
1626 matrix(col, col) =
1627 rotations[col].first * Hi + rotations[col].second * Hi1;
1628
1629 rhs(col + 1) = -rotations[col].second * rhs(col);
1630 rhs(col) *= rotations[col].first;
1631
1632 return std::abs(rhs(col + 1));
1633 }
1634 }
1635
1636
1637
1638 inline const Vector<double> &
1639 ArnoldiProcess::solve_projected_system(
1640 const bool orthogonalization_finished)
1641 {
1642 FullMatrix<double> tmp_triangular_matrix;
1643 Vector<double> tmp_rhs;
1644 FullMatrix<double> *matrix = &triangular_matrix;
1645 Vector<double> *rhs = &projected_rhs;
1646 unsigned int n = givens_rotations.size();
1647
1648 // If we solve with the delayed orthogonalization, we still need to
1649 // perform the elimination of the last column. We distinguish two cases,
1650 // one where the orthogonalization has finished (i.e., end of inner
1651 // iteration in GMRES) and we can safely overwrite the content of the
1652 // tridiagonal matrix and right hand side, and the case during the inner
1653 // iterations where we need to create copies of the matrices in the QR
1654 // decomposition as well as the right hand side.
1655 if (orthogonalization_strategy ==
1658 {
1659 n += 1;
1660 if (!orthogonalization_finished)
1661 {
1662 tmp_triangular_matrix = triangular_matrix;
1663 tmp_rhs = projected_rhs;
1664 std::vector<std::pair<double, double>> tmp_givens_rotations(
1665 givens_rotations);
1666 do_givens_rotation(false,
1667 givens_rotations.size(),
1668 tmp_triangular_matrix,
1669 tmp_givens_rotations,
1670 tmp_rhs);
1671 matrix = &tmp_triangular_matrix;
1672 rhs = &tmp_rhs;
1673 }
1674 else
1675 do_givens_rotation(false,
1676 givens_rotations.size(),
1677 triangular_matrix,
1678 givens_rotations,
1679 projected_rhs);
1680 }
1681
1682 // Now solve the triangular system by backward substitution
1683 projected_solution.reinit(n);
1684 for (int i = n - 1; i >= 0; --i)
1685 {
1686 double s = (*rhs)(i);
1687 for (unsigned int j = i + 1; j < n; ++j)
1688 s -= projected_solution(j) * (*matrix)(i, j);
1689 projected_solution(i) = s / (*matrix)(i, i);
1690 AssertIsFinite(projected_solution(i));
1691 }
1692
1693 return projected_solution;
1694 }
1695
1696
1697
1698 inline const FullMatrix<double> &
1699 ArnoldiProcess::get_hessenberg_matrix() const
1700 {
1701 return hessenberg_matrix;
1702 }
1703
1704
1705
1706 // A comparator for better printing eigenvalues
1707 inline bool
1708 complex_less_pred(const std::complex<double> &x,
1709 const std::complex<double> &y)
1710 {
1711 return x.real() < y.real() ||
1712 (x.real() == y.real() && x.imag() < y.imag());
1713 }
1714 } // namespace SolverGMRESImplementation
1715} // namespace internal
1716
1717
1718
1719template <typename VectorType>
1720inline void
1722 const FullMatrix<double> &H_orig,
1723 const unsigned int n,
1724 const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
1725 &eigenvalues_signal,
1726 const boost::signals2::signal<void(const FullMatrix<double> &)>
1727 &hessenberg_signal,
1728 const boost::signals2::signal<void(double)> &cond_signal)
1729{
1730 // Avoid copying the Hessenberg matrix if it isn't needed.
1731 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1732 !cond_signal.empty()) &&
1733 n > 0)
1734 {
1735 LAPACKFullMatrix<double> mat(n, n);
1736 for (unsigned int i = 0; i < n; ++i)
1737 for (unsigned int j = 0; j < n; ++j)
1738 mat(i, j) = H_orig(i, j);
1739 hessenberg_signal(H_orig);
1740 // Avoid computing eigenvalues if they are not needed.
1741 if (!eigenvalues_signal.empty())
1742 {
1743 // Copy mat so that we can compute svd below. Necessary since
1744 // compute_eigenvalues will leave mat in state
1745 // LAPACKSupport::unusable.
1746 LAPACKFullMatrix<double> mat_eig(mat);
1747 mat_eig.compute_eigenvalues();
1748 std::vector<std::complex<double>> eigenvalues(n);
1749 for (unsigned int i = 0; i < mat_eig.n(); ++i)
1750 eigenvalues[i] = mat_eig.eigenvalue(i);
1751 // Sort eigenvalues for nicer output.
1752 std::sort(eigenvalues.begin(),
1753 eigenvalues.end(),
1754 internal::SolverGMRESImplementation::complex_less_pred);
1755 eigenvalues_signal(eigenvalues);
1756 }
1757 // Calculate condition number, avoid calculating the svd if a slot
1758 // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
1759 if (!cond_signal.empty() && (mat.n() > 1))
1760 {
1761 mat.compute_svd();
1762 double condition_number =
1763 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1764 cond_signal(condition_number);
1765 }
1766 }
1767}
1768
1769
1770
1771template <typename VectorType>
1772template <typename MatrixType, typename PreconditionerType>
1773void
1774SolverGMRES<VectorType>::solve(const MatrixType &A,
1775 VectorType &x,
1776 const VectorType &b,
1777 const PreconditionerType &preconditioner)
1778{
1779 std::unique_ptr<LogStream::Prefix> prefix;
1780 if (!additional_data.batched_mode)
1781 prefix = std::make_unique<LogStream::Prefix>("GMRES");
1782
1783 // extra call to std::max to placate static analyzers: coverity rightfully
1784 // complains that data.max_n_tmp_vectors - 2 may overflow
1785 const unsigned int basis_size =
1786 (additional_data.max_basis_size > 0 ?
1787 additional_data.max_basis_size :
1788 std::max(additional_data.max_n_tmp_vectors, 3u) - 2);
1789
1790 // Generate an object where basis vectors are stored.
1792 basis_size + 2, this->memory);
1793
1794 // number of the present iteration; this number is not reset to zero upon a
1795 // restart
1796 unsigned int accumulated_iterations = 0;
1797
1798 const bool do_eigenvalues =
1799 !additional_data.batched_mode &&
1800 (!condition_number_signal.empty() ||
1801 !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1802 !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1803 !all_hessenberg_signal.empty());
1804
1806 double res = std::numeric_limits<double>::lowest();
1807
1808 // switch to determine whether we want a left or a right preconditioner. at
1809 // present, left is default, but both ways are implemented
1810 const bool left_precondition = !additional_data.right_preconditioning;
1811
1812 // Per default the left preconditioned GMRES uses the preconditioned
1813 // residual and the right preconditioned GMRES uses the unpreconditioned
1814 // residual as stopping criterion.
1815 const bool use_default_residual = additional_data.use_default_residual;
1816
1817 // define an alias
1818 VectorType &p = basis_vectors(basis_size + 1, x);
1819
1820 // Following vectors are needed when we are not using the default residuals
1821 // as stopping criterion
1824 if (!use_default_residual)
1825 {
1826 r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1827 x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1828 r->reinit(x);
1829 x_->reinit(x);
1830 }
1831
1832 arnoldi_process.initialize(additional_data.orthogonalization_strategy,
1833 basis_size,
1834 additional_data.force_re_orthogonalization);
1835
1837 // outer iteration: loop until we either reach convergence or the maximum
1838 // number of iterations is exceeded. each cycle of this loop amounts to one
1839 // restart
1840 do
1841 {
1842 VectorType &v = basis_vectors(0, x);
1843
1844 // Compute the preconditioned/unpreconditioned residual for left/right
1845 // preconditioning. If 'x' is the zero vector, then we can bypass the
1846 // full computation. But 'x' is only likely to be the zero vector if
1847 // that's what the user provided as the starting guess, so it's only
1848 // worth checking for this in the first iteration. (Calling all_zero()
1849 // costs as much in memory transfer and communication as computing the
1850 // norm of a vector.)
1851 if (left_precondition)
1852 {
1853 if (accumulated_iterations == 0 && x.all_zero())
1854 preconditioner.vmult(v, b);
1855 else
1856 {
1857 A.vmult(p, x);
1858 p.sadd(-1., 1., b);
1859 preconditioner.vmult(v, p);
1860 }
1861 }
1862 else
1863 {
1864 if (accumulated_iterations == 0 && x.all_zero())
1865 v = b;
1866 else
1867 {
1868 A.vmult(v, x);
1869 v.sadd(-1., 1., b);
1870 }
1871 }
1872
1873 const double norm_v =
1874 arnoldi_process.orthonormalize_nth_vector(0,
1875 basis_vectors,
1876 accumulated_iterations);
1877
1878 // check the residual here as well since it may be that we got the exact
1879 // (or an almost exact) solution vector at the outset. if we wouldn't
1880 // check here, the next scaling operation would produce garbage
1881 if (use_default_residual)
1882 {
1883 res = norm_v;
1884 if (additional_data.batched_mode)
1885 iteration_state = solver_control.check(accumulated_iterations, res);
1886 else
1887 iteration_state =
1888 this->iteration_status(accumulated_iterations, res, x);
1889
1890 if (iteration_state != SolverControl::iterate)
1891 break;
1892 }
1893 else
1894 {
1895 deallog << "default_res=" << norm_v << std::endl;
1896
1897 if (left_precondition)
1898 {
1899 A.vmult(*r, x);
1900 r->sadd(-1., 1., b);
1901 }
1902 else
1903 preconditioner.vmult(*r, v);
1904
1905 res = r->l2_norm();
1906 if (additional_data.batched_mode)
1907 iteration_state = solver_control.check(accumulated_iterations, res);
1908 else
1909 iteration_state =
1910 this->iteration_status(accumulated_iterations, res, x);
1911
1912 if (iteration_state != SolverControl::iterate)
1913 break;
1914 }
1915
1916 // inner iteration doing at most as many steps as the size of the
1917 // Arnoldi basis
1918 unsigned int inner_iteration = 0;
1919 for (; (inner_iteration < basis_size &&
1920 iteration_state == SolverControl::iterate);
1921 ++inner_iteration)
1922 {
1923 ++accumulated_iterations;
1924 // yet another alias
1925 VectorType &vv = basis_vectors(inner_iteration + 1, x);
1926
1927 if (left_precondition)
1928 {
1929 A.vmult(p, basis_vectors[inner_iteration]);
1930 preconditioner.vmult(vv, p);
1931 }
1932 else
1933 {
1934 preconditioner.vmult(p, basis_vectors[inner_iteration]);
1935 A.vmult(vv, p);
1936 }
1937
1938 res =
1939 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1,
1940 basis_vectors,
1941 accumulated_iterations,
1942 re_orthogonalize_signal);
1943
1944 if (use_default_residual)
1945 {
1946 if (additional_data.batched_mode)
1947 iteration_state =
1948 solver_control.check(accumulated_iterations, res);
1949 else
1950 iteration_state =
1951 this->iteration_status(accumulated_iterations, res, x);
1952 }
1953 else
1954 {
1955 if (!additional_data.batched_mode)
1956 deallog << "default_res=" << res << std::endl;
1957
1958 *x_ = x;
1959 const Vector<double> &projected_solution =
1960 arnoldi_process.solve_projected_system(false);
1961
1962 if (left_precondition)
1963 for (unsigned int i = 0; i < inner_iteration + 1; ++i)
1964 x_->add(projected_solution(i), basis_vectors[i]);
1965 else
1966 {
1967 p = 0.;
1968 for (unsigned int i = 0; i < inner_iteration + 1; ++i)
1969 p.add(projected_solution(i), basis_vectors[i]);
1970 preconditioner.vmult(*r, p);
1971 x_->add(1., *r);
1972 };
1973 A.vmult(*r, *x_);
1974 r->sadd(-1., 1., b);
1975
1976 // Now *r contains the unpreconditioned residual!!
1977 if (left_precondition)
1978 {
1979 res = r->l2_norm();
1980 iteration_state =
1981 this->iteration_status(accumulated_iterations, res, x);
1982 }
1983 else
1984 {
1985 preconditioner.vmult(*x_, *r);
1986 res = x_->l2_norm();
1987
1988 if (additional_data.batched_mode)
1989 iteration_state =
1990 solver_control.check(accumulated_iterations, res);
1991 else
1992 iteration_state =
1993 this->iteration_status(accumulated_iterations, res, x);
1994 }
1995 }
1996 }
1997
1998 // end of inner iteration; now update the global solution vector x with
1999 // the solution of the projected system (least-squares solution)
2000 const Vector<double> &projected_solution =
2001 arnoldi_process.solve_projected_system(true);
2002
2003 if (do_eigenvalues)
2004 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
2005 inner_iteration,
2006 all_eigenvalues_signal,
2007 all_hessenberg_signal,
2008 condition_number_signal);
2009
2010 if (left_precondition)
2011 ::internal::SolverGMRESImplementation::add(
2012 x, inner_iteration, projected_solution, basis_vectors, false);
2013 else
2014 {
2015 ::internal::SolverGMRESImplementation::add(
2016 p, inner_iteration, projected_solution, basis_vectors, true);
2017 preconditioner.vmult(v, p);
2018 x.add(1., v);
2019 }
2020
2021 // in the last round, print the eigenvalues from the last Arnoldi step
2022 if (iteration_state != SolverControl::iterate)
2023 {
2024 if (do_eigenvalues)
2025 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
2026 inner_iteration,
2027 eigenvalues_signal,
2028 hessenberg_signal,
2029 condition_number_signal);
2030
2031 if (!additional_data.batched_mode && !krylov_space_signal.empty())
2032 krylov_space_signal(basis_vectors);
2033
2034 // end of outer iteration. restart if no convergence and the number of
2035 // iterations is not exceeded
2036 }
2037 }
2038 while (iteration_state == SolverControl::iterate);
2039
2040 // in case of failure: throw exception
2041 AssertThrow(iteration_state == SolverControl::success,
2042 SolverControl::NoConvergence(accumulated_iterations, res));
2043}
2044
2045
2046
2047template <typename VectorType>
2048boost::signals2::connection
2050 const std::function<void(double)> &slot,
2051 const bool every_iteration)
2052{
2053 if (every_iteration)
2054 {
2055 return all_condition_numbers_signal.connect(slot);
2056 }
2057 else
2058 {
2059 return condition_number_signal.connect(slot);
2060 }
2061}
2062
2063
2064
2065template <typename VectorType>
2066boost::signals2::connection
2068 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
2069 const bool every_iteration)
2070{
2071 if (every_iteration)
2072 {
2073 return all_eigenvalues_signal.connect(slot);
2074 }
2075 else
2076 {
2077 return eigenvalues_signal.connect(slot);
2078 }
2079}
2080
2081
2082
2083template <typename VectorType>
2084boost::signals2::connection
2086 const std::function<void(const FullMatrix<double> &)> &slot,
2087 const bool every_iteration)
2088{
2089 if (every_iteration)
2090 {
2091 return all_hessenberg_signal.connect(slot);
2092 }
2093 else
2094 {
2095 return hessenberg_signal.connect(slot);
2096 }
2097}
2098
2099
2100
2101template <typename VectorType>
2102boost::signals2::connection
2104 const std::function<void(
2106{
2107 return krylov_space_signal.connect(slot);
2108}
2109
2110
2111
2112template <typename VectorType>
2113boost::signals2::connection
2115 const std::function<void(int)> &slot)
2116{
2117 return re_orthogonalize_signal.connect(slot);
2118}
2119
2120
2121
2122template <typename VectorType>
2123double
2125{
2126 // dummy implementation. this function is not needed for the present
2127 // implementation of gmres
2129 return 0;
2130}
2131
2132
2133//----------------------------------------------------------------------//
2134
2135template <typename VectorType>
2138 const AdditionalData &data)
2139 : SolverBase<VectorType>(cn, mem)
2140 , additional_data(data)
2141{}
2142
2143
2144
2145template <typename VectorType>
2147 const AdditionalData &data)
2148 : SolverBase<VectorType>(cn)
2149 , additional_data(data)
2150{}
2151
2152
2153
2154template <typename VectorType>
2155template <typename MatrixType, typename PreconditionerType>
2156void
2157SolverFGMRES<VectorType>::solve(const MatrixType &A,
2158 VectorType &x,
2159 const VectorType &b,
2160 const PreconditionerType &preconditioner)
2161{
2162 LogStream::Prefix prefix("FGMRES");
2163
2165
2166 const unsigned int basis_size = additional_data.max_basis_size;
2167
2168 // Generate an object where basis vectors are stored.
2170 basis_size + 1, this->memory);
2172 basis_size, this->memory);
2173
2174 // number of the present iteration; this number is not reset to zero upon a
2175 // restart
2176 unsigned int accumulated_iterations = 0;
2177
2178 // matrix used for the orthogonalization process later
2179 arnoldi_process.initialize(additional_data.orthogonalization_strategy,
2180 basis_size,
2181 false);
2182
2183 // Iteration starts here
2184 double res = std::numeric_limits<double>::lowest();
2185
2186 do
2187 {
2188 // Compute the residual. If 'x' is the zero vector, then we can bypass
2189 // the full computation. But 'x' is only likely to be the zero vector if
2190 // that's what the user provided as the starting guess, so it's only
2191 // worth checking for this in the first iteration. (Calling all_zero()
2192 // costs as much in memory transfer and communication as computing the
2193 // norm of a vector.)
2194 if (accumulated_iterations == 0 && x.all_zero())
2195 v(0, x) = b;
2196 else
2197 {
2198 A.vmult(v(0, x), x);
2199 v[0].sadd(-1., 1., b);
2200 }
2201
2202 res = arnoldi_process.orthonormalize_nth_vector(0, v);
2203 iteration_state = this->iteration_status(accumulated_iterations, res, x);
2204 if (iteration_state == SolverControl::success)
2205 break;
2206
2207 unsigned int inner_iteration = 0;
2208 for (; (inner_iteration < basis_size &&
2209 iteration_state == SolverControl::iterate);
2210 ++inner_iteration)
2211 {
2212 preconditioner.vmult(z(inner_iteration, x), v[inner_iteration]);
2213 A.vmult(v(inner_iteration + 1, x), z[inner_iteration]);
2214
2215 res =
2216 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1, v);
2217
2218 // check convergence. note that the vector 'x' we pass to the
2219 // criterion is not the final solution we compute if we
2220 // decide to jump out of the iteration (we update 'x' again
2221 // right after the current loop)
2222 iteration_state =
2223 this->iteration_status(++accumulated_iterations, res, x);
2224 }
2225
2226 // Solve triangular system with projected quantities and update solution
2227 // vector
2228 const Vector<double> &projected_solution =
2229 arnoldi_process.solve_projected_system(true);
2230 ::internal::SolverGMRESImplementation::add(
2231 x, inner_iteration, projected_solution, z, false);
2232 }
2233 while (iteration_state == SolverControl::iterate);
2234
2235 // in case of failure: throw exception
2236 if (iteration_state != SolverControl::success)
2237 AssertThrow(false,
2238 SolverControl::NoConvergence(accumulated_iterations, res));
2239}
2240
2241#endif // DOXYGEN
2242
2244
2245#endif
virtual State check(const unsigned int step, const double check_value)
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
internal::SolverGMRESImplementation::ArnoldiProcess arnoldi_process
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())
AdditionalData additional_data
boost::signals2::signal< void(const FullMatrix< double > &)> hessenberg_signal
AdditionalData additional_data
boost::signals2::signal< void(double)> condition_number_signal
SolverGMRES(const SolverGMRES< VectorType > &)=delete
boost::signals2::signal< void(const std::vector< std::complex< double > > &)> eigenvalues_signal
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
boost::signals2::connection connect_krylov_space_slot(const std::function< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> &slot)
boost::signals2::signal< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> krylov_space_signal
boost::signals2::signal< void(const FullMatrix< double > &)> all_hessenberg_signal
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
SolverControl & solver_control
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)
internal::SolverGMRESImplementation::ArnoldiProcess arnoldi_process
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)
boost::signals2::signal< void(const std::vector< std::complex< double > > &)> all_eigenvalues_signal
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::signal< void(double)> all_condition_numbers_signal
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< std::complex< double > > &)> &slot, const bool every_iteration=false)
boost::signals2::signal< void(int)> re_orthogonalize_signal
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
Number sum() const
void store(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
void initialize(const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy, const unsigned int max_basis_size, const bool force_reorthogonalization)
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)>())
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
std::vector< std::pair< double, double > > givens_rotations
const Vector< double > & solve_projected_system(const bool orthogonalization_finished)
double do_givens_rotation(const bool delayed_reorthogonalization, const int col, FullMatrix< double > &matrix, std::vector< std::pair< double, double > > &rotations, Vector< double > &rhs)
const FullMatrix< double > & get_hessenberg_matrix() const
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcTooFewTmpVectors(int arg1)
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:516
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
LogStream deallog
Definition logstream.cc:36
spacedim const Point< spacedim > & p
Definition grid_tools.h:981
@ matrix
Contents is actually a matrix.
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)
T sum(const T &t, const MPI_Comm mpi_communicator)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const Iterator & begin
Definition parallel.h:609
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::modified_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)