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