deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
solver_gmres.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 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
134 template <typename Number>
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
199 std::vector<const Number *> vector_ptrs;
200
201 private:
206
213
218 std::vector<std::pair<double, double>> givens_rotations;
219
224
230
235
245
250
278 double
279 do_givens_rotation(const bool delayed_reorthogonalization,
280 const int col,
281 FullMatrix<double> &matrix,
282 std::vector<std::pair<double, double>> &rotations,
283 Vector<double> &rhs);
284 };
285 } // namespace SolverGMRESImplementation
286} // namespace internal
287
355template <typename VectorType = Vector<double>>
357class SolverGMRES : public SolverBase<VectorType>
358{
359public:
364 {
375 explicit AdditionalData(const unsigned int max_basis_size = 30,
376 const bool right_preconditioning = false,
377 const bool use_default_residual = true,
378 const bool force_re_orthogonalization = false,
379 const bool batched_mode = false,
381 orthogonalization_strategy =
383 delayed_classical_gram_schmidt);
384
396 unsigned int max_n_tmp_vectors;
397
405 unsigned int max_basis_size;
406
415
420
428
436
441 };
442
448 const AdditionalData &data = AdditionalData());
449
455
460
464 template <typename MatrixType, typename PreconditionerType>
468 void solve(const MatrixType &A,
469 VectorType &x,
470 const VectorType &b,
471 const PreconditionerType &preconditioner);
472
479 boost::signals2::connection
480 connect_condition_number_slot(const std::function<void(double)> &slot,
481 const bool every_iteration = false);
482
489 boost::signals2::connection
490 connect_eigenvalues_slot(
491 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
492 const bool every_iteration = false);
493
501 boost::signals2::connection
502 connect_hessenberg_slot(
503 const std::function<void(const FullMatrix<double> &)> &slot,
504 const bool every_iteration = true);
505
512 boost::signals2::connection
513 connect_krylov_space_slot(
514 const std::function<
515 void(const internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
516 &slot);
517
518
523 boost::signals2::connection
524 connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
525
526
527 DeclException1(ExcTooFewTmpVectors,
528 int,
529 << "The number of temporary vectors you gave (" << arg1
530 << ") is too small. It should be at least 10 for "
531 << "any results, and much more for reasonable ones.");
532
533protected:
537 AdditionalData additional_data;
538
543 boost::signals2::signal<void(double)> condition_number_signal;
544
549 boost::signals2::signal<void(double)> all_condition_numbers_signal;
550
555 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
556 eigenvalues_signal;
557
562 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
563 all_eigenvalues_signal;
564
569 boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
570
575 boost::signals2::signal<void(const FullMatrix<double> &)>
576 all_hessenberg_signal;
577
582 boost::signals2::signal<void(
583 const internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
584 krylov_space_signal;
585
590 boost::signals2::signal<void(int)> re_orthogonalize_signal;
591
597 SolverControl &solver_control;
598
602 virtual double
603 criterion();
604
611 static void
612 compute_eigs_and_cond(
613 const FullMatrix<double> &H_orig,
614 const unsigned int n,
615 const boost::signals2::signal<
616 void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
617 const boost::signals2::signal<void(const FullMatrix<double> &)>
618 &hessenberg_signal,
619 const boost::signals2::signal<void(double)> &cond_signal);
620
625 internal::SolverGMRESImplementation::ArnoldiProcess<
626 typename VectorType::value_type>
627 arnoldi_process;
628};
629
630
631
652template <typename VectorType = Vector<double>>
653DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
654class SolverFGMRES : public SolverBase<VectorType>
655{
656public:
661 {
665 explicit AdditionalData(const unsigned int max_basis_size = 30,
667 orthogonalization_strategy =
669 delayed_classical_gram_schmidt)
670 : max_basis_size(max_basis_size)
671 , orthogonalization_strategy(orthogonalization_strategy)
672 {}
673
677 unsigned int max_basis_size;
678
683 };
684
690 const AdditionalData &data = AdditionalData());
691
697 const AdditionalData &data = AdditionalData());
698
702 template <typename MatrixType, typename PreconditionerType>
706 void solve(const MatrixType &A,
707 VectorType &x,
708 const VectorType &b,
709 const PreconditionerType &preconditioner);
710
711private:
715 AdditionalData additional_data;
716
721 internal::SolverGMRESImplementation::ArnoldiProcess<
722 typename VectorType::value_type>
723 arnoldi_process;
724};
725
727/* --------------------- Inline and template functions ------------------- */
728
729
730#ifndef DOXYGEN
731
732template <typename VectorType>
735 const unsigned int max_basis_size,
736 const bool right_preconditioning,
737 const bool use_default_residual,
738 const bool force_re_orthogonalization,
739 const bool batched_mode,
740 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy)
741 : max_n_tmp_vectors(0)
742 , max_basis_size(max_basis_size)
743 , right_preconditioning(right_preconditioning)
744 , use_default_residual(use_default_residual)
745 , force_re_orthogonalization(force_re_orthogonalization)
746 , batched_mode(batched_mode)
747 , orthogonalization_strategy(orthogonalization_strategy)
748{
749 Assert(max_basis_size >= 1,
750 ExcMessage("SolverGMRES needs at least one vector in the "
751 "Arnoldi basis."));
752}
753
754
755
756template <typename VectorType>
760 const AdditionalData &data)
761 : SolverBase<VectorType>(cn, mem)
762 , additional_data(data)
763 , solver_control(cn)
764{}
765
766
767
768template <typename VectorType>
771 const AdditionalData &data)
773 , additional_data(data)
774 , solver_control(cn)
775{}
776
777
778
779namespace internal
780{
781 namespace SolverGMRESImplementation
782 {
783 template <typename VectorType>
784 inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
786 : mem(vmem)
787 , data(max_size)
788 {}
789
790
791
792 template <typename VectorType>
793 inline VectorType &
794 TmpVectors<VectorType>::operator[](const unsigned int i) const
795 {
796 AssertIndexRange(i, data.size());
797
798 Assert(data[i] != nullptr, ExcNotInitialized());
799 return *data[i];
800 }
801
802
803
804 template <typename VectorType>
805 inline VectorType &
806 TmpVectors<VectorType>::operator()(const unsigned int i,
807 const VectorType &temp)
808 {
809 AssertIndexRange(i, data.size());
810 if (data[i] == nullptr)
811 {
812 data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
813 data[i]->reinit(temp, true);
814 }
815 return *data[i];
816 }
817
818
819
820 template <typename VectorType>
821 unsigned int
822 TmpVectors<VectorType>::size() const
823 {
824 return (data.size() > 0 ? data.size() - 1 : 0);
825 }
826
827
828
829 template <typename VectorType, typename Enable = void>
830 struct is_dealii_compatible_vector;
831
832 template <typename VectorType>
833 struct is_dealii_compatible_vector<
835 std::enable_if_t<!internal::is_block_vector<VectorType>>>
836 {
837 static constexpr bool value =
838 std::is_same_v<
840 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
842 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
843 };
844
845
846
847 template <typename VectorType>
848 struct is_dealii_compatible_vector<
850 std::enable_if_t<internal::is_block_vector<VectorType>>>
851 {
852 static constexpr bool value =
853 std::is_same_v<
854 typename VectorType::BlockType,
855 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
857 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
858 };
859
860
861
862 template <typename VectorType,
863 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
864 * = nullptr>
865 unsigned int
866 n_blocks(const VectorType &)
867 {
868 return 1;
869 }
870
871
872
873 template <typename VectorType,
874 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
875 nullptr>
876 unsigned int
877 n_blocks(const VectorType &vector)
878 {
879 return vector.n_blocks();
880 }
881
882
883
884 template <typename VectorType,
885 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
886 * = nullptr>
887 VectorType &
888 block(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 const VectorType &
901 block(const VectorType &vector, const unsigned int b)
902 {
903 AssertDimension(b, 0);
904 (void)b;
905 return vector;
906 }
907
908
909
910 template <typename VectorType,
911 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
912 nullptr>
913 typename VectorType::BlockType &
914 block(VectorType &vector, const unsigned int b)
915 {
916 return vector.block(b);
917 }
918
919
920
921 template <typename VectorType,
922 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
923 nullptr>
924 const typename VectorType::BlockType &
925 block(const VectorType &vector, const unsigned int b)
926 {
927 return vector.block(b);
928 }
929
930
931
932 template <bool delayed_reorthogonalization,
933 typename VectorType,
934 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
935 VectorType> * = nullptr>
936 void
937 Tvmult_add(const unsigned int n,
938 const VectorType &vv,
939 const TmpVectors<VectorType> &orthogonal_vectors,
941 std::vector<const typename VectorType::value_type *> &)
942 {
943 for (unsigned int i = 0; i < n; ++i)
944 {
945 h(i) += vv * orthogonal_vectors[i];
946 if (delayed_reorthogonalization)
947 h(n + i) += orthogonal_vectors[i] * orthogonal_vectors[n - 1];
948 }
949 if (delayed_reorthogonalization)
950 h(n + n) += vv * vv;
951 }
952
953
954
955 // worker method for deal.II's vector types implemented in .cc file
956 template <bool delayed_reorthogonalization, typename Number>
957 void
958 do_Tvmult_add(const unsigned int n_vectors,
959 const std::size_t locally_owned_size,
960 const Number *current_vector,
961 const std::vector<const Number *> &orthogonal_vectors,
962 Vector<double> &b);
963
964
965
966 template <bool delayed_reorthogonalization,
967 typename VectorType,
968 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
969 VectorType> * = nullptr>
970 void
971 Tvmult_add(
972 const unsigned int n,
973 const VectorType &vv,
974 const TmpVectors<VectorType> &orthogonal_vectors,
976 std::vector<const typename VectorType::value_type *> &vector_ptrs)
977 {
978 for (unsigned int b = 0; b < n_blocks(vv); ++b)
979 {
980 vector_ptrs.resize(n);
981 for (unsigned int i = 0; i < n; ++i)
982 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
983
984 do_Tvmult_add<delayed_reorthogonalization>(n,
985 block(vv, b).end() -
986 block(vv, b).begin(),
987 block(vv, b).begin(),
988 vector_ptrs,
989 h);
990 }
991
992 Utilities::MPI::sum(h, block(vv, 0).get_mpi_communicator(), h);
993 }
994
995
996
997 template <bool delayed_reorthogonalization,
998 typename VectorType,
999 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1000 VectorType> * = nullptr>
1001 double
1002 subtract_and_norm(const unsigned int n,
1003 const TmpVectors<VectorType> &orthogonal_vectors,
1004 const Vector<double> &h,
1005 VectorType &vv,
1006 std::vector<const typename VectorType::value_type *> &)
1007 {
1008 Assert(n > 0, ExcInternalError());
1009
1010 VectorType &last_vector =
1011 const_cast<VectorType &>(orthogonal_vectors[n - 1]);
1012 for (unsigned int i = 0; i < n - 1; ++i)
1013 {
1014 if (delayed_reorthogonalization && i + 2 < n)
1015 last_vector.add(-h(n + i), orthogonal_vectors[i]);
1016 vv.add(-h(i), orthogonal_vectors[i]);
1017 }
1018
1019 if (delayed_reorthogonalization)
1020 {
1021 if (n > 1)
1022 last_vector.sadd(1. / h(n + n - 1),
1023 -h(n + n - 2) / h(n + n - 1),
1024 orthogonal_vectors[n - 2]);
1025
1026 // h(n + n) is lucky breakdown
1027 const double scaling_factor_vv = h(n + n) > 0.0 ?
1028 1. / (h(n + n - 1) * h(n + n)) :
1029 1. / (h(n + n - 1) * h(n + n - 1));
1030 vv.sadd(scaling_factor_vv,
1031 -h(n - 1) * scaling_factor_vv,
1032 last_vector);
1033
1034 // the delayed reorthogonalization computes the norm from other
1035 // quantities
1036 return std::numeric_limits<double>::signaling_NaN();
1037 }
1038 else
1039 return std::sqrt(
1040 vv.add_and_dot(-h(n - 1), orthogonal_vectors[n - 1], vv));
1041 }
1042
1043
1044
1045 // worker method for deal.II's vector types implemented in .cc file
1046 template <bool delayed_reorthogonalization, typename Number>
1047 double
1048 do_subtract_and_norm(const unsigned int n_vectors,
1049 const std::size_t locally_owned_size,
1050 const std::vector<const Number *> &orthogonal_vectors,
1051 const Vector<double> &h,
1052 Number *current_vector);
1053
1054
1055
1056 template <bool delayed_reorthogonalization,
1057 typename VectorType,
1058 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1059 VectorType> * = nullptr>
1060 double
1061 subtract_and_norm(
1062 const unsigned int n,
1063 const TmpVectors<VectorType> &orthogonal_vectors,
1064 const Vector<double> &h,
1065 VectorType &vv,
1066 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1067 {
1068 double norm_vv_temp = 0.0;
1069
1070 for (unsigned int b = 0; b < n_blocks(vv); ++b)
1071 {
1072 vector_ptrs.resize(n);
1073 for (unsigned int i = 0; i < n; ++i)
1074 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
1075
1076 norm_vv_temp += do_subtract_and_norm<delayed_reorthogonalization>(
1077 n,
1078 block(vv, b).end() - block(vv, b).begin(),
1079 vector_ptrs,
1080 h,
1081 block(vv, b).begin());
1082 }
1083
1084 return std::sqrt(
1085 Utilities::MPI::sum(norm_vv_temp, block(vv, 0).get_mpi_communicator()));
1086 }
1087
1088
1089
1090 template <typename VectorType,
1091 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1092 VectorType> * = nullptr>
1093 void
1094 add(VectorType &p,
1095 const unsigned int n,
1096 const Vector<double> &h,
1097 const TmpVectors<VectorType> &tmp_vectors,
1098 const bool zero_out,
1099 std::vector<const typename VectorType::value_type *> &)
1100 {
1101 if (zero_out)
1102 p.equ(h(0), tmp_vectors[0]);
1103 else
1104 p.add(h(0), tmp_vectors[0]);
1105
1106 for (unsigned int i = 1; i < n; ++i)
1107 p.add(h(i), tmp_vectors[i]);
1108 }
1109
1110
1111
1112 // worker method for deal.II's vector types implemented in .cc file
1113 template <typename Number>
1114 void
1115 do_add(const unsigned int n_vectors,
1116 const std::size_t locally_owned_size,
1117 const std::vector<const Number *> &tmp_vectors,
1118 const Vector<double> &h,
1119 const bool zero_out,
1120 Number *output);
1121
1122
1123
1124 template <typename VectorType,
1125 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1126 VectorType> * = nullptr>
1127 void
1128 add(VectorType &p,
1129 const unsigned int n,
1130 const Vector<double> &h,
1131 const TmpVectors<VectorType> &tmp_vectors,
1132 const bool zero_out,
1133 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1134 {
1135 for (unsigned int b = 0; b < n_blocks(p); ++b)
1136 {
1137 vector_ptrs.resize(n);
1138 for (unsigned int i = 0; i < n; ++i)
1139 vector_ptrs[i] = block(tmp_vectors[i], b).begin();
1140 do_add(n,
1141 block(p, b).end() - block(p, b).begin(),
1142 vector_ptrs,
1143 h,
1144 zero_out,
1145 block(p, b).begin());
1146 }
1147 }
1148
1149
1150
1151 template <typename Number>
1152 inline void
1153 ArnoldiProcess<Number>::initialize(
1154 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy,
1155 const unsigned int basis_size,
1156 const bool force_reorthogonalization)
1157 {
1158 this->orthogonalization_strategy = orthogonalization_strategy;
1159 this->do_reorthogonalization = force_reorthogonalization;
1160
1161 hessenberg_matrix.reinit(basis_size + 1, basis_size);
1162 triangular_matrix.reinit(basis_size + 1, basis_size, true);
1163
1164 // some additional vectors, also used in the orthogonalization
1165 projected_rhs.reinit(basis_size + 1, true);
1166 givens_rotations.reserve(basis_size);
1167
1168 if (orthogonalization_strategy ==
1171 h.reinit(2 * basis_size + 3);
1172 else
1173 h.reinit(basis_size + 1);
1174 }
1175
1176
1177
1178 template <typename Number>
1179 template <typename VectorType>
1180 inline double
1181 ArnoldiProcess<Number>::orthonormalize_nth_vector(
1182 const unsigned int n,
1183 TmpVectors<VectorType> &orthogonal_vectors,
1184 const unsigned int accumulated_iterations,
1185 const boost::signals2::signal<void(int)> &reorthogonalize_signal)
1186 {
1187 AssertIndexRange(n, hessenberg_matrix.m());
1188 AssertIndexRange(n, orthogonal_vectors.size() + 1);
1189
1190 VectorType &vv = orthogonal_vectors[n];
1191
1192 double residual_estimate = std::numeric_limits<double>::signaling_NaN();
1193 if (n == 0)
1194 {
1195 givens_rotations.clear();
1196 residual_estimate = vv.l2_norm();
1197 if (residual_estimate != 0.)
1198 vv /= residual_estimate;
1199 projected_rhs(0) = residual_estimate;
1200 }
1201 else if (orthogonalization_strategy ==
1204 {
1205 // The algorithm implemented in the following few lines is algorithm
1206 // 4 of Bielich et al. (2022).
1207
1208 // To avoid un-scaled numbers as appearing with the original
1209 // algorithm by Bielich et al., we use a preliminary scaling of the
1210 // last vector. This will be corrected in the delayed step.
1211 const double previous_scaling = n > 0 ? h(n + n - 2) : 1.;
1212
1213 // Reset h to zero
1214 h.reinit(n + n + 1);
1215
1216 // global reduction
1217 Tvmult_add<true>(n, vv, orthogonal_vectors, h, vector_ptrs);
1218
1219 // delayed correction terms
1220 double tmp = 0;
1221 for (unsigned int i = 0; i < n - 1; ++i)
1222 tmp += h(n + i) * h(n + i);
1223 const double alpha_j = h(n + n - 1) > tmp ?
1224 std::sqrt(h(n + n - 1) - tmp) :
1225 std::sqrt(h(n + n - 1));
1226 h(n + n - 1) = alpha_j;
1227
1228 tmp = 0;
1229 for (unsigned int i = 0; i < n - 1; ++i)
1230 tmp += h(i) * h(n + i);
1231 h(n - 1) = (h(n - 1) - tmp) / alpha_j;
1232
1233 // representation of H(j-1)
1234 if (n > 1)
1235 {
1236 for (unsigned int i = 0; i < n - 1; ++i)
1237 hessenberg_matrix(i, n - 2) += h(n + i) * previous_scaling;
1238 hessenberg_matrix(n - 1, n - 2) = alpha_j * previous_scaling;
1239 }
1240 for (unsigned int i = 0; i < n; ++i)
1241 {
1242 double sum = 0;
1243 for (unsigned int j = (i == 0 ? 0 : i - 1); j < n - 1; ++j)
1244 sum += hessenberg_matrix(i, j) * h(n + j);
1245 hessenberg_matrix(i, n - 1) = (h(i) - sum) / alpha_j;
1246 }
1247
1248 // compute norm estimate for approximate convergence criterion
1249 // (value of norm to be corrected in next iteration)
1250 double sum = 0;
1251 for (unsigned int i = 0; i < n - 1; ++i)
1252 sum += h(i) * h(i);
1253 sum += (2. - 1.) * h(n - 1) * h(n - 1);
1254 hessenberg_matrix(n, n - 1) =
1255 std::sqrt(std::abs(h(n + n) - sum)) / alpha_j;
1256
1257 // projection and delayed reorthogonalization. We scale the vector
1258 // vv here by the preliminary norm to avoid working with too large
1259 // values and correct the actual norm in the Hessenberg matrix in
1260 // high precision in the next iteration.
1261 h(n + n) = hessenberg_matrix(n, n - 1);
1262 subtract_and_norm<true>(n, orthogonal_vectors, h, vv, vector_ptrs);
1263
1264 // transform new column of upper Hessenberg matrix into upper
1265 // triangular form by computing the respective factor
1266 residual_estimate = do_givens_rotation(
1267 true, n - 2, triangular_matrix, givens_rotations, projected_rhs);
1268 }
1269 else
1270 {
1271 // need initial norm for detection of re-orthogonalization, see below
1272 double norm_vv = 0.0;
1273 double norm_vv_start = 0;
1274 const bool consider_reorthogonalize =
1275 (do_reorthogonalization == false) && (n % 5 == 0);
1276 if (consider_reorthogonalize)
1277 norm_vv_start = vv.l2_norm();
1278
1279 // Reset h to zero
1280 h.reinit(n);
1281
1282 // run two loops with index 0: orthogonalize, 1: reorthogonalize
1283 for (unsigned int c = 0; c < 2; ++c)
1284 {
1285 // Orthogonalization
1286 if (orthogonalization_strategy ==
1289 {
1290 double htmp = vv * orthogonal_vectors[0];
1291 h(0) += htmp;
1292 for (unsigned int i = 1; i < n; ++i)
1293 {
1294 htmp = vv.add_and_dot(-htmp,
1295 orthogonal_vectors[i - 1],
1296 orthogonal_vectors[i]);
1297 h(i) += htmp;
1298 }
1299
1300 norm_vv = std::sqrt(
1301 vv.add_and_dot(-htmp, orthogonal_vectors[n - 1], vv));
1302 }
1303 else if (orthogonalization_strategy ==
1306 {
1307 Tvmult_add<false>(n, vv, orthogonal_vectors, h, vector_ptrs);
1308 norm_vv = subtract_and_norm<false>(
1309 n, orthogonal_vectors, h, vv, vector_ptrs);
1310 }
1311 else
1312 {
1314 }
1315
1316 if (c == 1)
1317 break; // reorthogonalization already performed -> finished
1318
1319 // Re-orthogonalization if loss of orthogonality detected. For the
1320 // test, use a strategy discussed in C. T. Kelley, Iterative
1321 // Methods for Linear and Nonlinear Equations, SIAM, Philadelphia,
1322 // 1995: Compare the norm of vv after orthogonalization with its
1323 // norm when starting the orthogonalization. If vv became very
1324 // small (here: less than the square root of the machine precision
1325 // times 10), it is almost in the span of the previous vectors,
1326 // which indicates loss of precision.
1327 if (consider_reorthogonalize)
1328 {
1329 if (norm_vv >
1330 10. * norm_vv_start *
1331 std::sqrt(std::numeric_limits<
1332 typename VectorType::value_type>::epsilon()))
1333 break;
1334
1335 else
1336 {
1337 do_reorthogonalization = true;
1338 if (!reorthogonalize_signal.empty())
1339 reorthogonalize_signal(accumulated_iterations);
1340 }
1341 }
1342
1343 if (do_reorthogonalization == false)
1344 break; // no reorthogonalization needed -> finished
1345 }
1346
1347 for (unsigned int i = 0; i < n; ++i)
1348 hessenberg_matrix(i, n - 1) = h(i);
1349 hessenberg_matrix(n, n - 1) = norm_vv;
1350
1351 // norm_vv is a lucky breakdown, the solver will reach convergence,
1352 // but we must not divide by zero here.
1353 if (norm_vv != 0)
1354 vv /= norm_vv;
1355
1356 residual_estimate = do_givens_rotation(
1357 false, n - 1, triangular_matrix, givens_rotations, projected_rhs);
1358 }
1359
1360 return residual_estimate;
1361 }
1362
1363
1364
1365 template <typename Number>
1366 inline double
1367 ArnoldiProcess<Number>::do_givens_rotation(
1368 const bool delayed_reorthogonalization,
1369 const int col,
1370 FullMatrix<double> &matrix,
1371 std::vector<std::pair<double, double>> &rotations,
1372 Vector<double> &rhs)
1373 {
1374 // for the delayed orthogonalization, we can only compute the column of
1375 // the previous iteration (as there will be correction terms added to the
1376 // present column for stability reasons), but we still want to compute
1377 // the residual estimate from the accumulated work; we therefore perform
1378 // givens rotations on two columns simultaneously
1379 if (delayed_reorthogonalization)
1380 {
1381 if (col >= 0)
1382 {
1383 AssertDimension(rotations.size(), static_cast<std::size_t>(col));
1384 matrix(0, col) = hessenberg_matrix(0, col);
1385 }
1386 double H_next = hessenberg_matrix(0, col + 1);
1387 for (int i = 0; i < col; ++i)
1388 {
1389 const double c = rotations[i].first;
1390 const double s = rotations[i].second;
1391 const double Hi = matrix(i, col);
1392 const double Hi1 = hessenberg_matrix(i + 1, col);
1393 H_next = -s * H_next + c * hessenberg_matrix(i + 1, col + 1);
1394 matrix(i, col) = c * Hi + s * Hi1;
1395 matrix(i + 1, col) = -s * Hi + c * Hi1;
1396 }
1397
1398 if (col >= 0)
1399 {
1400 const double H_col1 = hessenberg_matrix(col + 1, col);
1401 const double H_col = matrix(col, col);
1402 const double r = 1. / std::sqrt(H_col * H_col + H_col1 * H_col1);
1403 rotations.emplace_back(H_col * r, H_col1 * r);
1404 matrix(col, col) =
1405 rotations[col].first * H_col + rotations[col].second * H_col1;
1406
1407 rhs(col + 1) = -rotations[col].second * rhs(col);
1408 rhs(col) *= rotations[col].first;
1409
1410 H_next =
1411 -rotations[col].second * H_next +
1412 rotations[col].first * hessenberg_matrix(col + 1, col + 1);
1413 }
1414
1415 const double H_last = hessenberg_matrix(col + 2, col + 1);
1416 const double r = 1. / std::sqrt(H_next * H_next + H_last * H_last);
1417 return std::abs(H_last * r * rhs(col + 1));
1418 }
1419 else
1420 {
1421 AssertDimension(rotations.size(), static_cast<std::size_t>(col));
1422
1423 matrix(0, col) = hessenberg_matrix(0, col);
1424 for (int i = 0; i < col; ++i)
1425 {
1426 const double c = rotations[i].first;
1427 const double s = rotations[i].second;
1428 const double Hi = matrix(i, col);
1429 const double Hi1 = hessenberg_matrix(i + 1, col);
1430 matrix(i, col) = c * Hi + s * Hi1;
1431 matrix(i + 1, col) = -s * Hi + c * Hi1;
1432 }
1433
1434 const double Hi = matrix(col, col);
1435 const double Hi1 = hessenberg_matrix(col + 1, col);
1436 const double r = 1. / std::sqrt(Hi * Hi + Hi1 * Hi1);
1437 rotations.emplace_back(Hi * r, Hi1 * r);
1438 matrix(col, col) =
1439 rotations[col].first * Hi + rotations[col].second * Hi1;
1440
1441 rhs(col + 1) = -rotations[col].second * rhs(col);
1442 rhs(col) *= rotations[col].first;
1443
1444 return std::abs(rhs(col + 1));
1445 }
1446 }
1447
1448
1449
1450 template <typename Number>
1451 inline const Vector<double> &
1452 ArnoldiProcess<Number>::solve_projected_system(
1453 const bool orthogonalization_finished)
1454 {
1455 FullMatrix<double> tmp_triangular_matrix;
1456 Vector<double> tmp_rhs;
1457 FullMatrix<double> *matrix = &triangular_matrix;
1458 Vector<double> *rhs = &projected_rhs;
1459 unsigned int n = givens_rotations.size();
1460
1461 // If we solve with the delayed orthogonalization, we still need to
1462 // perform the elimination of the last column before we can solve the
1463 // projected system. We distinguish two cases, one where the
1464 // orthogonalization has finished (i.e., end of inner iteration in
1465 // GMRES) and we can safely overwrite the content of the tridiagonal
1466 // matrix and right hand side, and the case during the inner iterations,
1467 // where we need to create copies of the matrices in the QR
1468 // decomposition as well as the right hand side.
1469 if (orthogonalization_strategy ==
1472 {
1473 n += 1;
1474 if (!orthogonalization_finished)
1475 {
1476 tmp_triangular_matrix = triangular_matrix;
1477 tmp_rhs = projected_rhs;
1478 std::vector<std::pair<double, double>> tmp_givens_rotations(
1479 givens_rotations);
1480 do_givens_rotation(false,
1481 givens_rotations.size(),
1482 tmp_triangular_matrix,
1483 tmp_givens_rotations,
1484 tmp_rhs);
1485 matrix = &tmp_triangular_matrix;
1486 rhs = &tmp_rhs;
1487 }
1488 else
1489 do_givens_rotation(false,
1490 givens_rotations.size(),
1491 triangular_matrix,
1492 givens_rotations,
1493 projected_rhs);
1494 }
1495
1496 // Now solve the triangular system by backward substitution
1497 projected_solution.reinit(n);
1498 for (int i = n - 1; i >= 0; --i)
1499 {
1500 double s = (*rhs)(i);
1501 for (unsigned int j = i + 1; j < n; ++j)
1502 s -= projected_solution(j) * (*matrix)(i, j);
1503 projected_solution(i) = s / (*matrix)(i, i);
1504 AssertIsFinite(projected_solution(i));
1505 }
1506
1507 return projected_solution;
1508 }
1509
1510
1511
1512 template <typename Number>
1513 inline const FullMatrix<double> &
1514 ArnoldiProcess<Number>::get_hessenberg_matrix() const
1515 {
1516 return hessenberg_matrix;
1517 }
1518
1519
1520
1521 // A comparator for better printing eigenvalues
1522 inline bool
1523 complex_less_pred(const std::complex<double> &x,
1524 const std::complex<double> &y)
1525 {
1526 return x.real() < y.real() ||
1527 (x.real() == y.real() && x.imag() < y.imag());
1528 }
1529 } // namespace SolverGMRESImplementation
1530} // namespace internal
1531
1532
1533
1534template <typename VectorType>
1537 const FullMatrix<double> &H_orig,
1538 const unsigned int n,
1539 const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
1540 &eigenvalues_signal,
1541 const boost::signals2::signal<void(const FullMatrix<double> &)>
1542 &hessenberg_signal,
1543 const boost::signals2::signal<void(double)> &cond_signal)
1544{
1545 // Avoid copying the Hessenberg matrix if it isn't needed.
1546 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1547 !cond_signal.empty()) &&
1548 n > 0)
1549 {
1550 LAPACKFullMatrix<double> mat(n, n);
1551 for (unsigned int i = 0; i < n; ++i)
1552 for (unsigned int j = 0; j < n; ++j)
1553 mat(i, j) = H_orig(i, j);
1554 hessenberg_signal(H_orig);
1555 // Avoid computing eigenvalues if they are not needed.
1556 if (!eigenvalues_signal.empty())
1557 {
1558 // Copy mat so that we can compute svd below. Necessary since
1559 // compute_eigenvalues will leave mat in state
1560 // LAPACKSupport::unusable.
1561 LAPACKFullMatrix<double> mat_eig(mat);
1562 mat_eig.compute_eigenvalues();
1563 std::vector<std::complex<double>> eigenvalues(n);
1564 for (unsigned int i = 0; i < mat_eig.n(); ++i)
1565 eigenvalues[i] = mat_eig.eigenvalue(i);
1566 // Sort eigenvalues for nicer output.
1567 std::sort(eigenvalues.begin(),
1568 eigenvalues.end(),
1569 internal::SolverGMRESImplementation::complex_less_pred);
1570 eigenvalues_signal(eigenvalues);
1571 }
1572 // Calculate condition number, avoid calculating the svd if a slot
1573 // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
1574 if (!cond_signal.empty() && (mat.n() > 1))
1575 {
1576 mat.compute_svd();
1577 double condition_number =
1578 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1579 cond_signal(condition_number);
1580 }
1581 }
1582}
1583
1584
1585
1586template <typename VectorType>
1588template <typename MatrixType, typename PreconditionerType>
1592void SolverGMRES<VectorType>::solve(const MatrixType &A,
1593 VectorType &x,
1594 const VectorType &b,
1595 const PreconditionerType &preconditioner)
1596{
1597 std::unique_ptr<LogStream::Prefix> prefix;
1598 if (!additional_data.batched_mode)
1599 prefix = std::make_unique<LogStream::Prefix>("GMRES");
1600
1601 // extra call to std::max to placate static analyzers: coverity rightfully
1602 // complains that data.max_n_tmp_vectors - 2 may overflow
1603 const unsigned int basis_size =
1604 (additional_data.max_basis_size > 0 ?
1605 additional_data.max_basis_size :
1606 std::max(additional_data.max_n_tmp_vectors, 3u) - 2);
1607
1608 // Generate an object where basis vectors are stored.
1610 basis_size + 2, this->memory);
1611
1612 // number of the present iteration; this number is not reset to zero upon a
1613 // restart
1614 unsigned int accumulated_iterations = 0;
1615
1616 const bool do_eigenvalues =
1617 !additional_data.batched_mode &&
1618 (!condition_number_signal.empty() ||
1619 !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1620 !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1621 !all_hessenberg_signal.empty());
1622
1624 double res = std::numeric_limits<double>::lowest();
1625
1626 // switch to determine whether we want a left or a right preconditioner. at
1627 // present, left is default, but both ways are implemented
1628 const bool left_precondition = !additional_data.right_preconditioning;
1629
1630 // Per default the left preconditioned GMRES uses the preconditioned
1631 // residual and the right preconditioned GMRES uses the unpreconditioned
1632 // residual as stopping criterion.
1633 const bool use_default_residual = additional_data.use_default_residual;
1634
1635 // define an alias
1636 VectorType &p = basis_vectors(basis_size + 1, x);
1637
1638 // Following vectors are needed when we are not using the default residuals
1639 // as stopping criterion
1642 if (!use_default_residual)
1643 {
1644 r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1645 x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1646 r->reinit(x);
1647 x_->reinit(x);
1648 }
1649
1650 arnoldi_process.initialize(additional_data.orthogonalization_strategy,
1651 basis_size,
1652 additional_data.force_re_orthogonalization);
1653
1655 // outer iteration: loop until we either reach convergence or the maximum
1656 // number of iterations is exceeded. each cycle of this loop amounts to one
1657 // restart
1658 do
1659 {
1660 VectorType &v = basis_vectors(0, x);
1661
1662 // Compute the preconditioned/unpreconditioned residual for left/right
1663 // preconditioning. If 'x' is the zero vector, then we can bypass the
1664 // full computation. But 'x' is only likely to be the zero vector if
1665 // that's what the user provided as the starting guess, so it's only
1666 // worth checking for this in the first iteration. (Calling all_zero()
1667 // costs as much in memory transfer and communication as computing the
1668 // norm of a vector.)
1669 if (left_precondition)
1670 {
1671 if (accumulated_iterations == 0 && x.all_zero())
1672 preconditioner.vmult(v, b);
1673 else
1674 {
1675 A.vmult(p, x);
1676 p.sadd(-1., 1., b);
1677 preconditioner.vmult(v, p);
1678 }
1679 }
1680 else
1681 {
1682 if (accumulated_iterations == 0 && x.all_zero())
1683 v = b;
1684 else
1685 {
1686 A.vmult(v, x);
1687 v.sadd(-1., 1., b);
1688 }
1689 }
1690
1691 const double norm_v = arnoldi_process.orthonormalize_nth_vector(
1692 0, basis_vectors, accumulated_iterations, re_orthogonalize_signal);
1693
1694 // check the residual here as well since it may be that we got the exact
1695 // (or an almost exact) solution vector at the outset. if we wouldn't
1696 // check here, the next scaling operation would produce garbage
1697 if (use_default_residual)
1698 {
1699 res = norm_v;
1700 if (additional_data.batched_mode)
1701 iteration_state = solver_control.check(accumulated_iterations, res);
1702 else
1703 iteration_state =
1704 this->iteration_status(accumulated_iterations, res, x);
1705
1706 if (iteration_state != SolverControl::iterate)
1707 break;
1708 }
1709 else
1710 {
1711 deallog << "default_res=" << norm_v << std::endl;
1712
1713 if (left_precondition)
1714 {
1715 A.vmult(*r, x);
1716 r->sadd(-1., 1., b);
1717 }
1718 else
1719 preconditioner.vmult(*r, v);
1720
1721 res = r->l2_norm();
1722 if (additional_data.batched_mode)
1723 iteration_state = solver_control.check(accumulated_iterations, res);
1724 else
1725 iteration_state =
1726 this->iteration_status(accumulated_iterations, res, x);
1727
1728 if (iteration_state != SolverControl::iterate)
1729 break;
1730 }
1731
1732 // inner iteration doing at most as many steps as the size of the
1733 // Arnoldi basis
1734 unsigned int inner_iteration = 0;
1735 for (; (inner_iteration < basis_size &&
1736 iteration_state == SolverControl::iterate);
1737 ++inner_iteration)
1738 {
1739 ++accumulated_iterations;
1740 // yet another alias
1741 VectorType &vv = basis_vectors(inner_iteration + 1, x);
1742
1743 if (left_precondition)
1744 {
1745 A.vmult(p, basis_vectors[inner_iteration]);
1746 preconditioner.vmult(vv, p);
1747 }
1748 else
1749 {
1750 preconditioner.vmult(p, basis_vectors[inner_iteration]);
1751 A.vmult(vv, p);
1752 }
1753
1754 res =
1755 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1,
1756 basis_vectors,
1757 accumulated_iterations,
1758 re_orthogonalize_signal);
1759
1760 if (use_default_residual)
1761 {
1762 if (additional_data.batched_mode)
1763 iteration_state =
1764 solver_control.check(accumulated_iterations, res);
1765 else
1766 iteration_state =
1767 this->iteration_status(accumulated_iterations, res, x);
1768 }
1769 else
1770 {
1771 if (!additional_data.batched_mode)
1772 deallog << "default_res=" << res << std::endl;
1773
1774 *x_ = x;
1775 const Vector<double> &projected_solution =
1776 arnoldi_process.solve_projected_system(false);
1777
1778 if (left_precondition)
1779 for (unsigned int i = 0; i < inner_iteration + 1; ++i)
1780 x_->add(projected_solution(i), basis_vectors[i]);
1781 else
1782 {
1783 p = 0.;
1784 for (unsigned int i = 0; i < inner_iteration + 1; ++i)
1785 p.add(projected_solution(i), basis_vectors[i]);
1786 preconditioner.vmult(*r, p);
1787 x_->add(1., *r);
1788 };
1789 A.vmult(*r, *x_);
1790 r->sadd(-1., 1., b);
1791
1792 // Now *r contains the unpreconditioned residual!!
1793 if (left_precondition)
1794 {
1795 res = r->l2_norm();
1796 iteration_state =
1797 this->iteration_status(accumulated_iterations, res, x);
1798 }
1799 else
1800 {
1801 preconditioner.vmult(*x_, *r);
1802 res = x_->l2_norm();
1803
1804 if (additional_data.batched_mode)
1805 iteration_state =
1806 solver_control.check(accumulated_iterations, res);
1807 else
1808 iteration_state =
1809 this->iteration_status(accumulated_iterations, res, x);
1810 }
1811 }
1812 }
1813
1814 // end of inner iteration; now update the global solution vector x with
1815 // the solution of the projected system (least-squares solution)
1816 const Vector<double> &projected_solution =
1817 arnoldi_process.solve_projected_system(true);
1818
1819 if (do_eigenvalues)
1820 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
1821 inner_iteration,
1822 all_eigenvalues_signal,
1823 all_hessenberg_signal,
1824 condition_number_signal);
1825
1826 if (left_precondition)
1827 ::internal::SolverGMRESImplementation::add(
1828 x,
1829 inner_iteration,
1830 projected_solution,
1831 basis_vectors,
1832 false,
1833 arnoldi_process.vector_ptrs);
1834 else
1835 {
1836 ::internal::SolverGMRESImplementation::add(
1837 p,
1838 inner_iteration,
1839 projected_solution,
1840 basis_vectors,
1841 true,
1842 arnoldi_process.vector_ptrs);
1843 preconditioner.vmult(v, p);
1844 x.add(1., v);
1845 }
1846
1847 // in the last round, print the eigenvalues from the last Arnoldi step
1848 if (iteration_state != SolverControl::iterate)
1849 {
1850 if (do_eigenvalues)
1851 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
1852 inner_iteration,
1853 eigenvalues_signal,
1854 hessenberg_signal,
1855 condition_number_signal);
1856
1857 if (!additional_data.batched_mode && !krylov_space_signal.empty())
1858 krylov_space_signal(basis_vectors);
1859
1860 // end of outer iteration. restart if no convergence and the number of
1861 // iterations is not exceeded
1862 }
1863 }
1864 while (iteration_state == SolverControl::iterate);
1865
1866 // in case of failure: throw exception
1867 AssertThrow(iteration_state == SolverControl::success,
1868 SolverControl::NoConvergence(accumulated_iterations, res));
1869}
1870
1871
1872
1873template <typename VectorType>
1875boost::signals2::connection
1877 const std::function<void(double)> &slot,
1878 const bool every_iteration)
1879{
1880 if (every_iteration)
1881 {
1882 return all_condition_numbers_signal.connect(slot);
1883 }
1884 else
1885 {
1886 return condition_number_signal.connect(slot);
1887 }
1888}
1889
1890
1891
1892template <typename VectorType>
1894boost::signals2::connection SolverGMRES<VectorType>::connect_eigenvalues_slot(
1895 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1896 const bool every_iteration)
1897{
1898 if (every_iteration)
1899 {
1900 return all_eigenvalues_signal.connect(slot);
1901 }
1902 else
1903 {
1904 return eigenvalues_signal.connect(slot);
1905 }
1906}
1907
1908
1909
1910template <typename VectorType>
1912boost::signals2::connection SolverGMRES<VectorType>::connect_hessenberg_slot(
1913 const std::function<void(const FullMatrix<double> &)> &slot,
1914 const bool every_iteration)
1915{
1916 if (every_iteration)
1917 {
1918 return all_hessenberg_signal.connect(slot);
1919 }
1920 else
1921 {
1922 return hessenberg_signal.connect(slot);
1923 }
1924}
1925
1926
1927
1928template <typename VectorType>
1930boost::signals2::connection SolverGMRES<VectorType>::connect_krylov_space_slot(
1931 const std::function<void(
1933{
1934 return krylov_space_signal.connect(slot);
1935}
1936
1937
1938
1939template <typename VectorType>
1941boost::signals2::connection
1943 const std::function<void(int)> &slot)
1944{
1945 return re_orthogonalize_signal.connect(slot);
1946}
1947
1948
1949
1950template <typename VectorType>
1953{
1954 // dummy implementation. this function is not needed for the present
1955 // implementation of gmres
1957 return 0;
1958}
1959
1960
1961
1962//----------------------------------------------------------------------//
1963
1964template <typename VectorType>
1968 const AdditionalData &data)
1969 : SolverBase<VectorType>(cn, mem)
1970 , additional_data(data)
1971{}
1972
1973
1974
1975template <typename VectorType>
1978 const AdditionalData &data)
1979 : SolverBase<VectorType>(cn)
1980 , additional_data(data)
1981{}
1982
1983
1984
1985template <typename VectorType>
1987template <typename MatrixType, typename PreconditionerType>
1991void SolverFGMRES<VectorType>::solve(const MatrixType &A,
1992 VectorType &x,
1993 const VectorType &b,
1994 const PreconditionerType &preconditioner)
1995{
1996 LogStream::Prefix prefix("FGMRES");
1997
1999
2000 const unsigned int basis_size = additional_data.max_basis_size;
2001
2002 // Generate an object where basis vectors are stored.
2004 basis_size + 1, this->memory);
2006 basis_size, this->memory);
2007
2008 // number of the present iteration; this number is not reset to zero upon a
2009 // restart
2010 unsigned int accumulated_iterations = 0;
2011
2012 // matrix used for the orthogonalization process later
2013 arnoldi_process.initialize(additional_data.orthogonalization_strategy,
2014 basis_size,
2015 false);
2016
2017 // Iteration starts here
2018 double res = std::numeric_limits<double>::lowest();
2019
2020 do
2021 {
2022 // Compute the residual. If 'x' is the zero vector, then we can bypass
2023 // the full computation. But 'x' is only likely to be the zero vector if
2024 // that's what the user provided as the starting guess, so it's only
2025 // worth checking for this in the first iteration. (Calling all_zero()
2026 // costs as much in memory transfer and communication as computing the
2027 // norm of a vector.)
2028 if (accumulated_iterations == 0 && x.all_zero())
2029 v(0, x) = b;
2030 else
2031 {
2032 A.vmult(v(0, x), x);
2033 v[0].sadd(-1., 1., b);
2034 }
2035
2036 res = arnoldi_process.orthonormalize_nth_vector(0, v);
2037 iteration_state = this->iteration_status(accumulated_iterations, res, x);
2038 if (iteration_state == SolverControl::success)
2039 break;
2040
2041 unsigned int inner_iteration = 0;
2042 for (; (inner_iteration < basis_size &&
2043 iteration_state == SolverControl::iterate);
2044 ++inner_iteration)
2045 {
2046 preconditioner.vmult(z(inner_iteration, x), v[inner_iteration]);
2047 A.vmult(v(inner_iteration + 1, x), z[inner_iteration]);
2048
2049 res =
2050 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1, v);
2051
2052 // check convergence. note that the vector 'x' we pass to the
2053 // criterion is not the final solution we compute if we
2054 // decide to jump out of the iteration (we update 'x' again
2055 // right after the current loop)
2056 iteration_state =
2057 this->iteration_status(++accumulated_iterations, res, x);
2058 }
2059
2060 // Solve triangular system with projected quantities and update solution
2061 // vector
2062 const Vector<double> &projected_solution =
2063 arnoldi_process.solve_projected_system(true);
2064 ::internal::SolverGMRESImplementation::add(
2065 x,
2066 inner_iteration,
2067 projected_solution,
2068 z,
2069 false,
2070 arnoldi_process.vector_ptrs);
2071 }
2072 while (iteration_state == SolverControl::iterate);
2073
2074 // in case of failure: throw exception
2075 if (iteration_state != SolverControl::success)
2076 AssertThrow(false,
2077 SolverControl::NoConvergence(accumulated_iterations, res));
2078}
2079
2080#endif // DOXYGEN
2081
2083
2084#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)
const FullMatrix< double > & get_hessenberg_matrix() const
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
double orthonormalize_nth_vector(const unsigned int n, TmpVectors< VectorType > &orthogonal_vectors, const unsigned int accumulated_iterations=0, const boost::signals2::signal< void(int)> &reorthogonalize_signal=boost::signals2::signal< void(int)>())
const Vector< double > & solve_projected_system(const bool orthogonalization_finished)
void initialize(const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy, const unsigned int max_basis_size, const bool force_reorthogonalization)
std::vector< std::pair< double, double > > givens_rotations
double do_givens_rotation(const bool delayed_reorthogonalization, const int col, FullMatrix< double > &matrix, std::vector< std::pair< double, double > > &rotations, Vector< double > &rhs)
std::vector< typename VectorMemory< VectorType >::Pointer > data
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
VectorType & operator[](const unsigned int i) const
VectorType & operator()(const unsigned int i, const VectorType &temp)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:175
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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:511
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 * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
void do_Tvmult_add(const unsigned int n_vectors, const std::size_t locally_owned_size, const Number *current_vector, const std::vector< const Number * > &orthogonal_vectors, Vector< double > &h)
double do_subtract_and_norm(const unsigned int n_vectors, const std::size_t locally_owned_size, const std::vector< const Number * > &orthogonal_vectors, const Vector< double > &h, Number *current_vector)
void do_add(const unsigned int n_vectors, const std::size_t locally_owned_size, const std::vector< const Number * > &tmp_vectors, const Vector< double > &h, const bool zero_out, Number *output)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
AdditionalData(const unsigned int max_basis_size=30, const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy=LinearAlgebra::OrthogonalizationStrategy::delayed_classical_gram_schmidt)
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
AdditionalData(const unsigned int max_basis_size=30, const bool right_preconditioning=false, const bool use_default_residual=true, const bool force_re_orthogonalization=false, const bool batched_mode=false, const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy=LinearAlgebra::OrthogonalizationStrategy::delayed_classical_gram_schmidt)
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)