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