Reference documentation for deal.II version GIT relicensing-43-g3fc2fae7c3 2024-03-03 21:10: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
25
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 {
75 template <typename VectorType>
77 {
78 public:
83 TmpVectors(const unsigned int max_size, VectorMemory<VectorType> &vmem);
84
88 ~TmpVectors() = default;
89
94 VectorType &
95 operator[](const unsigned int i) const;
96
103 VectorType &
104 operator()(const unsigned int i, const VectorType &temp);
105
110 unsigned int
111 size() const;
112
113
114 private:
119
123 std::vector<typename VectorMemory<VectorType>::Pointer> data;
124 };
125 } // namespace SolverGMRESImplementation
126} // namespace internal
127
192template <typename VectorType = Vector<double>>
193class SolverGMRES : public SolverBase<VectorType>
194{
195public:
261
267 const AdditionalData &data = AdditionalData());
268
274
279
283 template <typename MatrixType, typename PreconditionerType>
284 void
285 solve(const MatrixType &A,
286 VectorType &x,
287 const VectorType &b,
288 const PreconditionerType &preconditioner);
289
296 boost::signals2::connection
297 connect_condition_number_slot(const std::function<void(double)> &slot,
298 const bool every_iteration = false);
299
306 boost::signals2::connection
308 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
309 const bool every_iteration = false);
310
318 boost::signals2::connection
320 const std::function<void(const FullMatrix<double> &)> &slot,
321 const bool every_iteration = true);
322
329 boost::signals2::connection
331 const std::function<
333 &slot);
334
335
340 boost::signals2::connection
341 connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
342
343
345 int,
346 << "The number of temporary vectors you gave (" << arg1
347 << ") is too small. It should be at least 10 for "
348 << "any results, and much more for reasonable ones.");
349
350protected:
355
360 boost::signals2::signal<void(double)> condition_number_signal;
361
366 boost::signals2::signal<void(double)> all_condition_numbers_signal;
367
372 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
374
379 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
381
386 boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
387
392 boost::signals2::signal<void(const FullMatrix<double> &)>
394
399 boost::signals2::signal<void(
402
407 boost::signals2::signal<void(int)> re_orthogonalize_signal;
408
415
419 virtual double
421
426 void
431 int col) const;
432
439 static void
441 const FullMatrix<double> &H_orig,
442 const unsigned int dim,
443 const boost::signals2::signal<
444 void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
445 const boost::signals2::signal<void(const FullMatrix<double> &)>
447 const boost::signals2::signal<void(double)> &cond_signal);
448
453
458
463
468
473};
474
475
476
497template <typename VectorType = Vector<double>>
498class SolverFGMRES : public SolverBase<VectorType>
499{
500public:
528
534 const AdditionalData &data = AdditionalData());
535
541 const AdditionalData &data = AdditionalData());
542
546 template <typename MatrixType, typename PreconditionerType>
547 void
548 solve(const MatrixType &A,
549 VectorType &x,
550 const VectorType &b,
551 const PreconditionerType &preconditioner);
552
553private:
558
563
568};
569
571/* --------------------- Inline and template functions ------------------- */
572
573
574#ifndef DOXYGEN
575namespace internal
576{
577 namespace SolverGMRESImplementation
578 {
579 template <typename VectorType>
580 inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
582 : mem(vmem)
583 , data(max_size)
584 {}
585
586
587
588 template <typename VectorType>
589 inline VectorType &
590 TmpVectors<VectorType>::operator[](const unsigned int i) const
591 {
592 AssertIndexRange(i, data.size());
593
594 Assert(data[i] != nullptr, ExcNotInitialized());
595 return *data[i];
596 }
597
598
599
600 template <typename VectorType>
601 inline VectorType &
602 TmpVectors<VectorType>::operator()(const unsigned int i,
603 const VectorType &temp)
604 {
605 AssertIndexRange(i, data.size());
606 if (data[i] == nullptr)
607 {
608 data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
609 data[i]->reinit(temp, true);
610 }
611 return *data[i];
612 }
613
614
615
616 template <typename VectorType>
617 unsigned int
619 {
620 return (data.size() > 0 ? data.size() - 1 : 0);
621 }
622
623
624
625 // A comparator for better printing eigenvalues
626 inline bool
627 complex_less_pred(const std::complex<double> &x,
628 const std::complex<double> &y)
629 {
630 return x.real() < y.real() ||
631 (x.real() == y.real() && x.imag() < y.imag());
632 }
633
634 // A function to solve the (upper) triangular system after Givens
635 // rotations on a matrix that has possibly unused rows and columns
636 inline void
637 solve_triangular(const unsigned int dim,
638 const FullMatrix<double> &H,
639 const Vector<double> &rhs,
640 Vector<double> &solution)
641 {
642 for (int i = dim - 1; i >= 0; --i)
643 {
644 double s = rhs(i);
645 for (unsigned int j = i + 1; j < dim; ++j)
646 s -= solution(j) * H(i, j);
647 solution(i) = s / H(i, i);
648 AssertIsFinite(solution(i));
649 }
650 }
651 } // namespace SolverGMRESImplementation
652} // namespace internal
653
654
655
656template <typename VectorType>
658 const unsigned int max_n_tmp_vectors,
659 const bool right_preconditioning,
660 const bool use_default_residual,
661 const bool force_re_orthogonalization,
662 const bool batched_mode,
663 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy)
664 : max_n_tmp_vectors(max_n_tmp_vectors)
665 , right_preconditioning(right_preconditioning)
666 , use_default_residual(use_default_residual)
667 , force_re_orthogonalization(force_re_orthogonalization)
668 , batched_mode(batched_mode)
669 , orthogonalization_strategy(orthogonalization_strategy)
670{
671 Assert(3 <= max_n_tmp_vectors,
672 ExcMessage("SolverGMRES needs at least three "
673 "temporary vectors."));
674}
675
676
677
678template <typename VectorType>
681 const AdditionalData &data)
682 : SolverBase<VectorType>(cn, mem)
683 , additional_data(data)
684 , solver_control(cn)
685{}
686
687
688
689template <typename VectorType>
691 const AdditionalData &data)
692 : SolverBase<VectorType>(cn)
693 , additional_data(data)
694 , solver_control(cn)
695{}
696
697
698
699template <typename VectorType>
700inline void
703 Vector<double> &ci,
704 Vector<double> &si,
705 int col) const
706{
707 for (int i = 0; i < col; ++i)
708 {
709 const double s = si(i);
710 const double c = ci(i);
711 const double dummy = h(i);
712 h(i) = c * dummy + s * h(i + 1);
713 h(i + 1) = -s * dummy + c * h(i + 1);
714 };
715
716 const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
717 si(col) = h(col + 1) * r;
718 ci(col) = h(col) * r;
719 h(col) = ci(col) * h(col) + si(col) * h(col + 1);
720 b(col + 1) = -si(col) * b(col);
721 b(col) *= ci(col);
722}
723
724
725
726namespace internal
727{
728 namespace SolverGMRESImplementation
729 {
730 template <typename VectorType, typename Enable = void>
731 struct is_dealii_compatible_distributed_vector;
732
733 template <typename VectorType>
734 struct is_dealii_compatible_distributed_vector<
735 VectorType,
736 std::enable_if_t<!internal::is_block_vector<VectorType>>>
737 {
738 static constexpr bool value = std::is_same_v<
739 VectorType,
740 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
742 };
743
744
745
746 template <typename VectorType>
747 struct is_dealii_compatible_distributed_vector<
748 VectorType,
749 std::enable_if_t<internal::is_block_vector<VectorType>>>
750 {
751 static constexpr bool value = std::is_same_v<
752 typename VectorType::BlockType,
753 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
755 };
756
757
758
759 template <typename VectorType,
760 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
761 * = nullptr>
762 unsigned int
763 n_blocks(const VectorType &)
764 {
765 return 1;
766 }
767
768
769
770 template <typename VectorType,
771 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
772 nullptr>
773 unsigned int
774 n_blocks(const VectorType &vector)
775 {
776 return vector.n_blocks();
777 }
778
779
780
781 template <typename VectorType,
782 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
783 * = nullptr>
784 VectorType &
785 block(VectorType &vector, const unsigned int b)
786 {
787 AssertDimension(b, 0);
788 (void)b;
789 return vector;
790 }
791
792
793
794 template <typename VectorType,
795 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
796 * = nullptr>
797 const VectorType &
798 block(const VectorType &vector, const unsigned int b)
799 {
800 AssertDimension(b, 0);
801 (void)b;
802 return vector;
803 }
804
805
806
807 template <typename VectorType,
808 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
809 nullptr>
810 typename VectorType::BlockType &
811 block(VectorType &vector, const unsigned int b)
812 {
813 return vector.block(b);
814 }
815
816
817
818 template <typename VectorType,
819 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
820 nullptr>
821 const typename VectorType::BlockType &
822 block(const VectorType &vector, const unsigned int b)
823 {
824 return vector.block(b);
825 }
826
827
828
829 template <typename VectorType,
830 std::enable_if_t<
831 !is_dealii_compatible_distributed_vector<VectorType>::value,
832 VectorType> * = nullptr>
833 void
834 Tvmult_add(const unsigned int dim,
835 const VectorType &vv,
837 &orthogonal_vectors,
839 {
840 for (unsigned int i = 0; i < dim; ++i)
841 h[i] += vv * orthogonal_vectors[i];
842 }
843
844
845
846 template <typename VectorType,
847 std::enable_if_t<
848 is_dealii_compatible_distributed_vector<VectorType>::value,
849 VectorType> * = nullptr>
850 void
851 Tvmult_add(const unsigned int dim,
852 const VectorType &vv,
854 &orthogonal_vectors,
856 {
857 for (unsigned int b = 0; b < n_blocks(vv); ++b)
858 {
859 unsigned int j = 0;
860
861 if (dim <= 128)
862 {
863 // optimized path
864 static constexpr unsigned int n_lanes =
866
868 for (unsigned int d = 0; d < dim; ++d)
869 hs[d] = 0.0;
870
871 unsigned int c = 0;
872
873 for (; c < block(vv, b).locally_owned_size() / n_lanes / 4;
874 ++c, j += n_lanes * 4)
875 for (unsigned int i = 0; i < dim; ++i)
876 {
878 for (unsigned int k = 0; k < 4; ++k)
879 vvec[k].load(block(vv, b).begin() + j + k * n_lanes);
880
881 for (unsigned int k = 0; k < 4; ++k)
882 {
884 temp.load(block(orthogonal_vectors[i], b).begin() + j +
885 k * n_lanes);
886 hs[i] += temp * vvec[k];
887 }
888 }
889
890 c *= 4;
891 for (; c < block(vv, b).locally_owned_size() / n_lanes;
892 ++c, j += n_lanes)
893 for (unsigned int i = 0; i < dim; ++i)
894 {
895 VectorizedArray<double> vvec, temp;
896 vvec.load(block(vv, b).begin() + j);
897 temp.load(block(orthogonal_vectors[i], b).begin() + j);
898 hs[i] += temp * vvec;
899 }
900
901 for (unsigned int i = 0; i < dim; ++i)
902 for (unsigned int v = 0; v < n_lanes; ++v)
903 h(i) += hs[i][v];
904 }
905
906 // remainder loop of optimized path or non-optimized path (if
907 // dim>128)
908 for (; j < block(vv, b).locally_owned_size(); ++j)
909 for (unsigned int i = 0; i < dim; ++i)
910 h(i) += block(orthogonal_vectors[i], b).local_element(j) *
911 block(vv, b).local_element(j);
912 }
913
914 Utilities::MPI::sum(h, block(vv, 0).get_mpi_communicator(), h);
915 }
916
917
918
919 template <typename VectorType,
920 std::enable_if_t<
921 !is_dealii_compatible_distributed_vector<VectorType>::value,
922 VectorType> * = nullptr>
923 double
924 subtract_and_norm(
925 const unsigned int dim,
927 &orthogonal_vectors,
928 const Vector<double> &h,
929 VectorType &vv)
930 {
931 Assert(dim > 0, ExcInternalError());
932
933 for (unsigned int i = 0; i < dim; ++i)
934 vv.add(-h(i), orthogonal_vectors[i]);
935
936 return std::sqrt(vv.add_and_dot(-h(dim), orthogonal_vectors[dim], vv));
937 }
938
939
940
941 template <typename VectorType,
942 std::enable_if_t<
943 is_dealii_compatible_distributed_vector<VectorType>::value,
944 VectorType> * = nullptr>
945 double
946 subtract_and_norm(
947 const unsigned int dim,
949 &orthogonal_vectors,
950 const Vector<double> &h,
951 VectorType &vv)
952 {
953 static constexpr unsigned int n_lanes = VectorizedArray<double>::size();
954
955 double norm_vv_temp = 0.0;
956
957 for (unsigned int b = 0; b < n_blocks(vv); ++b)
958 {
959 VectorizedArray<double> norm_vv_temp_vectorized = 0.0;
960
961 unsigned int j = 0;
962 unsigned int c = 0;
963 for (; c < block(vv, b).locally_owned_size() / n_lanes / 4;
964 ++c, j += n_lanes * 4)
965 {
967
968 for (unsigned int k = 0; k < 4; ++k)
969 temp[k].load(block(vv, b).begin() + j + k * n_lanes);
970
971 for (unsigned int i = 0; i < dim; ++i)
972 {
973 const double factor = h(i);
974 for (unsigned int k = 0; k < 4; ++k)
975 {
977 vec.load(block(orthogonal_vectors[i], b).begin() + j +
978 k * n_lanes);
979 temp[k] -= factor * vec;
980 }
981 }
982
983 for (unsigned int k = 0; k < 4; ++k)
984 temp[k].store(block(vv, b).begin() + j + k * n_lanes);
985
986 norm_vv_temp_vectorized +=
987 (temp[0] * temp[0] + temp[1] * temp[1]) +
988 (temp[2] * temp[2] + temp[3] * temp[3]);
989 }
990
991 c *= 4;
992 for (; c < block(vv, b).locally_owned_size() / n_lanes;
993 ++c, j += n_lanes)
994 {
996 temp.load(block(vv, b).begin() + j);
997
998 for (unsigned int i = 0; i < dim; ++i)
999 {
1001 vec.load(block(orthogonal_vectors[i], b).begin() + j);
1002 temp -= h(i) * vec;
1003 }
1004
1005 temp.store(block(vv, b).begin() + j);
1006
1007 norm_vv_temp_vectorized += temp * temp;
1008 }
1009
1010 for (unsigned int v = 0; v < n_lanes; ++v)
1011 norm_vv_temp += norm_vv_temp_vectorized[v];
1012
1013 for (; j < block(vv, b).locally_owned_size(); ++j)
1014 {
1015 double temp = block(vv, b).local_element(j);
1016 for (unsigned int i = 0; i < dim; ++i)
1017 temp -= h(i) * block(orthogonal_vectors[i], b).local_element(j);
1018 block(vv, b).local_element(j) = temp;
1019
1020 norm_vv_temp += temp * temp;
1021 }
1022 }
1023
1024 return std::sqrt(
1025 Utilities::MPI::sum(norm_vv_temp, block(vv, 0).get_mpi_communicator()));
1026 }
1027
1028
1029 template <typename VectorType,
1030 std::enable_if_t<
1031 !is_dealii_compatible_distributed_vector<VectorType>::value,
1032 VectorType> * = nullptr>
1033 double
1034 sadd_and_norm(VectorType &v,
1035 const double factor_a,
1036 const VectorType &b,
1037 const double factor_b)
1038 {
1039 v.sadd(factor_a, factor_b, b);
1040 return v.l2_norm();
1041 }
1042
1043
1044 template <typename VectorType,
1045 std::enable_if_t<
1046 is_dealii_compatible_distributed_vector<VectorType>::value,
1047 VectorType> * = nullptr>
1048 double
1049 sadd_and_norm(VectorType &v,
1050 const double factor_a,
1051 const VectorType &w,
1052 const double factor_b)
1053 {
1054 double norm = 0;
1055
1056 for (unsigned int b = 0; b < n_blocks(v); ++b)
1057 for (unsigned int j = 0; j < block(v, b).locally_owned_size(); ++j)
1058 {
1059 const double temp = block(v, b).local_element(j) * factor_a +
1060 block(w, b).local_element(j) * factor_b;
1061
1062 block(v, b).local_element(j) = temp;
1063
1064 norm += temp * temp;
1065 }
1066
1067 return std::sqrt(
1068 Utilities::MPI::sum(norm, block(v, 0).get_mpi_communicator()));
1069 }
1070
1071
1072
1073 template <typename VectorType,
1074 std::enable_if_t<
1075 !is_dealii_compatible_distributed_vector<VectorType>::value,
1076 VectorType> * = nullptr>
1077 void
1078 add(VectorType &p,
1079 const unsigned int dim,
1080 const Vector<double> &h,
1082 &tmp_vectors,
1083 const bool zero_out)
1084 {
1085 if (zero_out)
1086 p.equ(h(0), tmp_vectors[0]);
1087 else
1088 p.add(h(0), tmp_vectors[0]);
1089
1090 for (unsigned int i = 1; i < dim; ++i)
1091 p.add(h(i), tmp_vectors[i]);
1092 }
1093
1094
1095
1096 template <typename VectorType,
1097 std::enable_if_t<
1098 is_dealii_compatible_distributed_vector<VectorType>::value,
1099 VectorType> * = nullptr>
1100 void
1101 add(VectorType &p,
1102 const unsigned int dim,
1103 const Vector<double> &h,
1105 &tmp_vectors,
1106 const bool zero_out)
1107 {
1108 for (unsigned int b = 0; b < n_blocks(p); ++b)
1109 for (unsigned int j = 0; j < block(p, b).locally_owned_size(); ++j)
1110 {
1111 double temp = zero_out ? 0 : block(p, b).local_element(j);
1112 for (unsigned int i = 0; i < dim; ++i)
1113 temp += block(tmp_vectors[i], b).local_element(j) * h(i);
1114 block(p, b).local_element(j) = temp;
1115 }
1116 }
1117
1118
1119
1131 template <typename VectorType>
1132 inline double
1133 iterated_gram_schmidt(
1134 const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy,
1136 &orthogonal_vectors,
1137 const unsigned int dim,
1138 const unsigned int accumulated_iterations,
1139 VectorType &vv,
1140 Vector<double> &h,
1141 bool &reorthogonalize,
1142 const boost::signals2::signal<void(int)> &reorthogonalize_signal =
1143 boost::signals2::signal<void(int)>())
1144 {
1145 Assert(dim > 0, ExcInternalError());
1146 const unsigned int inner_iteration = dim - 1;
1147
1148 // need initial norm for detection of re-orthogonalization, see below
1149 double norm_vv_start = 0;
1150 const bool consider_reorthogonalize =
1151 (reorthogonalize == false) && (inner_iteration % 5 == 4);
1152 if (consider_reorthogonalize)
1153 norm_vv_start = vv.l2_norm();
1154
1155 for (unsigned int i = 0; i < dim; ++i)
1156 h[i] = 0;
1157
1158 for (unsigned int c = 0; c < 2;
1159 ++c) // 0: orthogonalize, 1: reorthogonalize
1160 {
1161 // Orthogonalization
1162 double norm_vv = 0.0;
1163
1164 if (orthogonalization_strategy ==
1166 {
1167 double htmp = vv * orthogonal_vectors[0];
1168 h(0) += htmp;
1169 for (unsigned int i = 1; i < dim; ++i)
1170 {
1171 htmp = vv.add_and_dot(-htmp,
1172 orthogonal_vectors[i - 1],
1173 orthogonal_vectors[i]);
1174 h(i) += htmp;
1175 }
1176
1177 norm_vv = std::sqrt(
1178 vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
1179 }
1180 else if (orthogonalization_strategy ==
1183 {
1184 Tvmult_add(dim, vv, orthogonal_vectors, h);
1185 norm_vv = subtract_and_norm(dim, orthogonal_vectors, h, vv);
1186 }
1187 else
1188 {
1190 }
1191
1192 if (c == 1)
1193 return norm_vv; // reorthogonalization already performed -> finished
1194
1195 // Re-orthogonalization if loss of orthogonality detected. For the
1196 // test, use a strategy discussed in C. T. Kelley, Iterative Methods
1197 // for Linear and Nonlinear Equations, SIAM, Philadelphia, 1995:
1198 // Compare the norm of vv after orthogonalization with its norm when
1199 // starting the orthogonalization. If vv became very small (here: less
1200 // than the square root of the machine precision times 10), it is
1201 // almost in the span of the previous vectors, which indicates loss of
1202 // precision.
1203 if (consider_reorthogonalize)
1204 {
1205 if (norm_vv >
1206 10. * norm_vv_start *
1207 std::sqrt(std::numeric_limits<
1208 typename VectorType::value_type>::epsilon()))
1209 return norm_vv;
1210
1211 else
1212 {
1213 reorthogonalize = true;
1214 if (!reorthogonalize_signal.empty())
1215 reorthogonalize_signal(accumulated_iterations);
1216 }
1217 }
1218
1219 if (reorthogonalize == false)
1220 return norm_vv; // no reorthogonalization needed -> finished
1221 }
1222
1223 AssertThrow(false, ExcInternalError());
1224
1225 return 0.0;
1226 }
1227 } // namespace SolverGMRESImplementation
1228} // namespace internal
1229
1230
1231
1232template <typename VectorType>
1233inline void
1235 const FullMatrix<double> &H_orig,
1236 const unsigned int dim,
1237 const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
1238 &eigenvalues_signal,
1239 const boost::signals2::signal<void(const FullMatrix<double> &)>
1240 &hessenberg_signal,
1241 const boost::signals2::signal<void(double)> &cond_signal)
1242{
1243 // Avoid copying the Hessenberg matrix if it isn't needed.
1244 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1245 !cond_signal.empty()) &&
1246 dim > 0)
1247 {
1248 LAPACKFullMatrix<double> mat(dim, dim);
1249 for (unsigned int i = 0; i < dim; ++i)
1250 for (unsigned int j = 0; j < dim; ++j)
1251 mat(i, j) = H_orig(i, j);
1252 hessenberg_signal(H_orig);
1253 // Avoid computing eigenvalues if they are not needed.
1254 if (!eigenvalues_signal.empty())
1255 {
1256 // Copy mat so that we can compute svd below. Necessary since
1257 // compute_eigenvalues will leave mat in state
1258 // LAPACKSupport::unusable.
1259 LAPACKFullMatrix<double> mat_eig(mat);
1260 mat_eig.compute_eigenvalues();
1261 std::vector<std::complex<double>> eigenvalues(dim);
1262 for (unsigned int i = 0; i < mat_eig.n(); ++i)
1263 eigenvalues[i] = mat_eig.eigenvalue(i);
1264 // Sort eigenvalues for nicer output.
1265 std::sort(eigenvalues.begin(),
1266 eigenvalues.end(),
1267 internal::SolverGMRESImplementation::complex_less_pred);
1268 eigenvalues_signal(eigenvalues);
1269 }
1270 // Calculate condition number, avoid calculating the svd if a slot
1271 // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
1272 if (!cond_signal.empty() && (mat.n() > 1))
1273 {
1274 mat.compute_svd();
1275 double condition_number =
1276 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1277 cond_signal(condition_number);
1278 }
1279 }
1280}
1281
1282
1283
1284template <typename VectorType>
1285template <typename MatrixType, typename PreconditionerType>
1286void
1287SolverGMRES<VectorType>::solve(const MatrixType &A,
1288 VectorType &x,
1289 const VectorType &b,
1290 const PreconditionerType &preconditioner)
1291{
1292 // TODO:[?] Check, why there are two different start residuals.
1293 // TODO:[GK] Make sure the parameter in the constructor means maximum basis
1294 // size
1295
1296 std::unique_ptr<LogStream::Prefix> prefix;
1297 if (!additional_data.batched_mode)
1298 prefix = std::make_unique<LogStream::Prefix>("GMRES");
1299
1300 // extra call to std::max to placate static analyzers: coverity rightfully
1301 // complains that data.max_n_tmp_vectors - 2 may overflow
1302 const unsigned int n_tmp_vectors =
1303 std::max(additional_data.max_n_tmp_vectors, 3u);
1304
1305 // Generate an object where basis vectors are stored.
1307 n_tmp_vectors, this->memory);
1308
1309 // number of the present iteration; this
1310 // number is not reset to zero upon a
1311 // restart
1312 unsigned int accumulated_iterations = 0;
1313
1314 const bool do_eigenvalues =
1315 !additional_data.batched_mode &&
1316 (!condition_number_signal.empty() ||
1317 !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1318 !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1319 !all_hessenberg_signal.empty());
1320 // for eigenvalue computation, need to collect the Hessenberg matrix (before
1321 // applying Givens rotations)
1322 FullMatrix<double> H_orig;
1323 if (do_eigenvalues)
1324 H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
1325
1326 // matrix used for the orthogonalization process later
1327 H.reinit(n_tmp_vectors, n_tmp_vectors - 1, /* omit_initialization */ true);
1328
1329 // some additional vectors, also used in the orthogonalization
1330 gamma.reinit(n_tmp_vectors);
1331 ci.reinit(n_tmp_vectors - 1);
1332 si.reinit(n_tmp_vectors - 1);
1333 h.reinit(n_tmp_vectors - 1);
1334
1335 unsigned int dim = 0;
1336
1338 double last_res = std::numeric_limits<double>::lowest();
1339
1340 // switch to determine whether we want a left or a right preconditioner. at
1341 // present, left is default, but both ways are implemented
1342 const bool left_precondition = !additional_data.right_preconditioning;
1343
1344 // Per default the left preconditioned GMRes uses the preconditioned
1345 // residual and the right preconditioned GMRes uses the unpreconditioned
1346 // residual as stopping criterion.
1347 const bool use_default_residual = additional_data.use_default_residual;
1348
1349 // define two aliases
1350 VectorType &v = tmp_vectors(0, x);
1351 VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
1352
1353 // Following vectors are needed when we are not using the default residuals
1354 // as stopping criterion
1357 std::unique_ptr<::Vector<double>> gamma_;
1358 if (!use_default_residual)
1359 {
1360 r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1361 x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
1362 r->reinit(x);
1363 x_->reinit(x);
1364
1365 gamma_ = std::make_unique<::Vector<double>>(gamma.size());
1366 }
1367
1368 bool re_orthogonalize = additional_data.force_re_orthogonalization;
1369
1371 // outer iteration: loop until we either reach convergence or the maximum
1372 // number of iterations is exceeded. each cycle of this loop amounts to one
1373 // restart
1374 do
1375 {
1376 // reset this vector to the right size
1377 h.reinit(n_tmp_vectors - 1);
1378
1379 double rho = 0.0;
1380
1381 if (left_precondition)
1382 {
1383 A.vmult(p, x);
1384 p.sadd(-1., 1., b);
1385 preconditioner.vmult(v, p);
1386 rho = v.l2_norm();
1387 }
1388 else
1389 {
1390 A.vmult(v, x);
1391 rho = ::internal::SolverGMRESImplementation::sadd_and_norm(v,
1392 -1,
1393 b,
1394 1.0);
1395 }
1396
1397 // check the residual here as well since it may be that we got the exact
1398 // (or an almost exact) solution vector at the outset. if we wouldn't
1399 // check here, the next scaling operation would produce garbage
1400 if (use_default_residual)
1401 {
1402 last_res = rho;
1403 if (additional_data.batched_mode)
1404 iteration_state = solver_control.check(accumulated_iterations, rho);
1405 else
1406 iteration_state =
1407 this->iteration_status(accumulated_iterations, rho, x);
1408
1409 if (iteration_state != SolverControl::iterate)
1410 break;
1411 }
1412 else
1413 {
1414 deallog << "default_res=" << rho << std::endl;
1415
1416 if (left_precondition)
1417 {
1418 A.vmult(*r, x);
1419 r->sadd(-1., 1., b);
1420 }
1421 else
1422 preconditioner.vmult(*r, v);
1423
1424 double res = r->l2_norm();
1425 last_res = res;
1426 if (additional_data.batched_mode)
1427 iteration_state = solver_control.check(accumulated_iterations, rho);
1428 else
1429 iteration_state =
1430 this->iteration_status(accumulated_iterations, res, x);
1431
1432 if (iteration_state != SolverControl::iterate)
1433 break;
1434 }
1435
1436 gamma(0) = rho;
1437
1438 v *= 1. / rho;
1439
1440 // inner iteration doing at most as many steps as there are temporary
1441 // vectors. the number of steps actually been done is propagated outside
1442 // through the @p dim variable
1443 for (unsigned int inner_iteration = 0;
1444 ((inner_iteration < n_tmp_vectors - 2) &&
1445 (iteration_state == SolverControl::iterate));
1446 ++inner_iteration)
1447 {
1448 ++accumulated_iterations;
1449 // yet another alias
1450 VectorType &vv = tmp_vectors(inner_iteration + 1, x);
1451
1452 if (left_precondition)
1453 {
1454 A.vmult(p, tmp_vectors[inner_iteration]);
1455 preconditioner.vmult(vv, p);
1456 }
1457 else
1458 {
1459 preconditioner.vmult(p, tmp_vectors[inner_iteration]);
1460 A.vmult(vv, p);
1461 }
1462
1463 dim = inner_iteration + 1;
1464
1465 const double s =
1466 internal::SolverGMRESImplementation::iterated_gram_schmidt(
1467 additional_data.orthogonalization_strategy,
1468 tmp_vectors,
1469 dim,
1470 accumulated_iterations,
1471 vv,
1472 h,
1473 re_orthogonalize,
1474 re_orthogonalize_signal);
1475 h(inner_iteration + 1) = s;
1476
1477 // s=0 is a lucky breakdown, the solver will reach convergence,
1478 // but we must not divide by zero here.
1479 if (s != 0)
1480 vv *= 1. / s;
1481
1482 // for eigenvalues, get the resulting coefficients from the
1483 // orthogonalization process
1484 if (do_eigenvalues)
1485 for (unsigned int i = 0; i < dim + 1; ++i)
1486 H_orig(i, inner_iteration) = h(i);
1487
1488 // Transformation into tridiagonal structure
1489 givens_rotation(h, gamma, ci, si, inner_iteration);
1490
1491 // append vector on matrix
1492 for (unsigned int i = 0; i < dim; ++i)
1493 H(i, inner_iteration) = h(i);
1494
1495 // default residual
1496 rho = std::fabs(gamma(dim));
1497
1498 if (use_default_residual)
1499 {
1500 last_res = rho;
1501 if (additional_data.batched_mode)
1502 iteration_state =
1503 solver_control.check(accumulated_iterations, rho);
1504 else
1505 iteration_state =
1506 this->iteration_status(accumulated_iterations, rho, x);
1507 }
1508 else
1509 {
1510 if (!additional_data.batched_mode)
1511 deallog << "default_res=" << rho << std::endl;
1512
1513 *x_ = x;
1514 *gamma_ = gamma;
1515 internal::SolverGMRESImplementation::solve_triangular(dim,
1516 H,
1517 *gamma_,
1518 h);
1519
1520 if (left_precondition)
1521 for (unsigned int i = 0; i < dim; ++i)
1522 x_->add(h(i), tmp_vectors[i]);
1523 else
1524 {
1525 p = 0.;
1526 for (unsigned int i = 0; i < dim; ++i)
1527 p.add(h(i), tmp_vectors[i]);
1528 preconditioner.vmult(*r, p);
1529 x_->add(1., *r);
1530 };
1531 A.vmult(*r, *x_);
1532 r->sadd(-1., 1., b);
1533 // Now *r contains the unpreconditioned residual!!
1534 if (left_precondition)
1535 {
1536 const double res = r->l2_norm();
1537 last_res = res;
1538
1539 iteration_state =
1540 this->iteration_status(accumulated_iterations, res, x);
1541 }
1542 else
1543 {
1544 preconditioner.vmult(*x_, *r);
1545 const double preconditioned_res = x_->l2_norm();
1546 last_res = preconditioned_res;
1547
1548 if (additional_data.batched_mode)
1549 iteration_state =
1550 solver_control.check(accumulated_iterations, rho);
1551 else
1552 iteration_state =
1553 this->iteration_status(accumulated_iterations,
1554 preconditioned_res,
1555 x);
1556 }
1557 }
1558 }
1559
1560 // end of inner iteration. now calculate the solution from the temporary
1561 // vectors
1562 internal::SolverGMRESImplementation::solve_triangular(dim, H, gamma, h);
1563
1564 if (do_eigenvalues)
1565 compute_eigs_and_cond(H_orig,
1566 dim,
1567 all_eigenvalues_signal,
1568 all_hessenberg_signal,
1569 condition_number_signal);
1570
1571 if (left_precondition)
1572 ::internal::SolverGMRESImplementation::add(
1573 x, dim, h, tmp_vectors, false);
1574 else
1575 {
1576 ::internal::SolverGMRESImplementation::add(
1577 p, dim, h, tmp_vectors, true);
1578 preconditioner.vmult(v, p);
1579 x.add(1., v);
1580 };
1581 // end of outer iteration. restart if no convergence and the number of
1582 // iterations is not exceeded
1583 }
1584 while (iteration_state == SolverControl::iterate);
1585
1586 if (do_eigenvalues)
1587 compute_eigs_and_cond(H_orig,
1588 dim,
1589 eigenvalues_signal,
1590 hessenberg_signal,
1591 condition_number_signal);
1592
1593 if (!additional_data.batched_mode && !krylov_space_signal.empty())
1594 krylov_space_signal(tmp_vectors);
1595
1596 // in case of failure: throw exception
1597 AssertThrow(iteration_state == SolverControl::success,
1598 SolverControl::NoConvergence(accumulated_iterations, last_res));
1599}
1600
1601
1602
1603template <typename VectorType>
1604boost::signals2::connection
1606 const std::function<void(double)> &slot,
1607 const bool every_iteration)
1608{
1609 if (every_iteration)
1610 {
1611 return all_condition_numbers_signal.connect(slot);
1612 }
1613 else
1614 {
1615 return condition_number_signal.connect(slot);
1616 }
1617}
1618
1619
1620
1621template <typename VectorType>
1622boost::signals2::connection
1624 const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1625 const bool every_iteration)
1626{
1627 if (every_iteration)
1628 {
1629 return all_eigenvalues_signal.connect(slot);
1630 }
1631 else
1632 {
1633 return eigenvalues_signal.connect(slot);
1634 }
1635}
1636
1637
1638
1639template <typename VectorType>
1640boost::signals2::connection
1642 const std::function<void(const FullMatrix<double> &)> &slot,
1643 const bool every_iteration)
1644{
1645 if (every_iteration)
1646 {
1647 return all_hessenberg_signal.connect(slot);
1648 }
1649 else
1650 {
1651 return hessenberg_signal.connect(slot);
1652 }
1653}
1654
1655
1656
1657template <typename VectorType>
1658boost::signals2::connection
1660 const std::function<void(
1662{
1663 return krylov_space_signal.connect(slot);
1664}
1665
1666
1667
1668template <typename VectorType>
1669boost::signals2::connection
1671 const std::function<void(int)> &slot)
1672{
1673 return re_orthogonalize_signal.connect(slot);
1674}
1675
1676
1677
1678template <typename VectorType>
1679double
1681{
1682 // dummy implementation. this function is not needed for the present
1683 // implementation of gmres
1685 return 0;
1686}
1687
1688
1689//----------------------------------------------------------------------//
1690
1691template <typename VectorType>
1694 const AdditionalData &data)
1695 : SolverBase<VectorType>(cn, mem)
1696 , additional_data(data)
1697{}
1698
1699
1700
1701template <typename VectorType>
1703 const AdditionalData &data)
1704 : SolverBase<VectorType>(cn)
1705 , additional_data(data)
1706{}
1707
1708
1709
1710template <typename VectorType>
1711template <typename MatrixType, typename PreconditionerType>
1712void
1713SolverFGMRES<VectorType>::solve(const MatrixType &A,
1714 VectorType &x,
1715 const VectorType &b,
1716 const PreconditionerType &preconditioner)
1717{
1718 LogStream::Prefix prefix("FGMRES");
1719
1721
1722 const unsigned int basis_size = additional_data.max_basis_size;
1723
1724 // Generate an object where basis vectors are stored.
1726 basis_size, this->memory);
1728 basis_size, this->memory);
1729
1730 // number of the present iteration; this number is not reset to zero upon a
1731 // restart
1732 unsigned int accumulated_iterations = 0;
1733
1734 // matrix used for the orthogonalization process later
1735 H.reinit(basis_size + 1, basis_size);
1736
1737 Vector<double> h(basis_size + 1);
1738
1739 // Vectors for projected system
1740 Vector<double> projected_rhs;
1742
1743 // Iteration starts here
1744 double res = std::numeric_limits<double>::lowest();
1745
1746 typename VectorMemory<VectorType>::Pointer aux(this->memory);
1747 aux->reinit(x);
1748 do
1749 {
1750 A.vmult(*aux, x);
1751 aux->sadd(-1., 1., b);
1752
1753 double beta = aux->l2_norm();
1754 res = beta;
1755 iteration_state = this->iteration_status(accumulated_iterations, res, x);
1756 if (iteration_state == SolverControl::success)
1757 break;
1758
1759 H.reinit(basis_size + 1, basis_size);
1760 double a = beta;
1761
1762 for (unsigned int j = 0; j < basis_size; ++j)
1763 {
1764 if (a != 0) // treat lucky breakdown
1765 v(j, x).equ(1. / a, *aux);
1766 else
1767 v(j, x) = 0.;
1768
1769
1770 preconditioner.vmult(z(j, x), v[j]);
1771 A.vmult(*aux, z[j]);
1772
1773 // Gram-Schmidt
1774 bool re_orthogonalize = false;
1775 const double s =
1776 internal::SolverGMRESImplementation::iterated_gram_schmidt<
1777 VectorType>(additional_data.orthogonalization_strategy,
1778 v,
1779 j + 1,
1780 0,
1781 *aux,
1782 h,
1783 re_orthogonalize);
1784 for (unsigned int i = 0; i <= j; ++i)
1785 H(i, j) = h(i);
1786 H(j + 1, j) = a = s;
1787
1788 // Compute projected solution
1789
1790 if (j > 0)
1791 {
1792 H1.reinit(j + 1, j);
1793 projected_rhs.reinit(j + 1);
1794 y.reinit(j);
1795 projected_rhs(0) = beta;
1796 H1.fill(H);
1797
1798 // check convergence. note that the vector 'x' we pass to the
1799 // criterion is not the final solution we compute if we
1800 // decide to jump out of the iteration (we update 'x' again
1801 // right after the current loop)
1802 Householder<double> house(H1);
1803 res = house.least_squares(y, projected_rhs);
1804 iteration_state =
1805 this->iteration_status(++accumulated_iterations, res, x);
1806 if (iteration_state != SolverControl::iterate)
1807 break;
1808 }
1809 }
1810
1811 // Update solution vector
1812 for (unsigned int j = 0; j < y.size(); ++j)
1813 x.add(y(j), z[j]);
1814 }
1815 while (iteration_state == SolverControl::iterate);
1816
1817 // in case of failure: throw exception
1818 if (iteration_state != SolverControl::success)
1819 AssertThrow(false,
1820 SolverControl::NoConvergence(accumulated_iterations, res));
1821}
1822
1823#endif // DOXYGEN
1824
1826
1827#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)
FullMatrix< double > H
SolverFGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
FullMatrix< double > H1
boost::signals2::signal< void(const FullMatrix< double > &)> hessenberg_signal
Vector< double > si
AdditionalData additional_data
boost::signals2::signal< void(double)> condition_number_signal
SolverGMRES(const SolverGMRES< VectorType > &)=delete
boost::signals2::signal< void(const std::vector< std::complex< double > > &)> eigenvalues_signal
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
boost::signals2::connection connect_krylov_space_slot(const std::function< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> &slot)
boost::signals2::signal< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> krylov_space_signal
boost::signals2::signal< void(const FullMatrix< double > &)> all_hessenberg_signal
Vector< double > ci
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
SolverControl & solver_control
SolverGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
void givens_rotation(Vector< double > &h, Vector< double > &b, Vector< double > &ci, Vector< double > &si, int col) const
boost::signals2::connection connect_hessenberg_slot(const std::function< void(const FullMatrix< double > &)> &slot, const bool every_iteration=true)
Vector< double > gamma
boost::signals2::signal< void(const std::vector< std::complex< double > > &)> all_eigenvalues_signal
Vector< double > h
static void compute_eigs_and_cond(const FullMatrix< double > &H_orig, const unsigned int dim, 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, const AdditionalData &data=AdditionalData())
FullMatrix< double > H
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
virtual double criterion()
boost::signals2::signal< void(double)> all_condition_numbers_signal
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< std::complex< double > > &)> &slot, const bool every_iteration=false)
boost::signals2::signal< void(int)> re_orthogonalize_signal
Number add_and_dot(const Number a, const Vector< Number > &V, const Vector< Number > &W)
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
void store(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
std::vector< typename VectorMemory< VectorType >::Pointer > data
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
VectorType & operator[](const unsigned int i) const
VectorType & operator()(const unsigned int i, const VectorType &temp)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcTooFewTmpVectors(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
LogStream deallog
Definition logstream.cc:36
spacedim const Point< spacedim > & p
Definition grid_tools.h:981
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
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)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
T sum(const T &t, const MPI_Comm mpi_communicator)
long double gamma(const unsigned int n)
const Iterator & begin
Definition parallel.h:609
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
AdditionalData(const unsigned int max_basis_size=30, const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy=LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt)
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
AdditionalData(const unsigned int max_n_tmp_vectors=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::modified_gram_schmidt)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)