Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3067-g24f3489fdc 2025-04-14 21:10: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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
precondition.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_precondition_h
16#define dealii_precondition_h
17
18#include <deal.II/base/config.h>
19
21#include <deal.II/base/mutex.h>
25
32
33#include <Kokkos_Core.hpp>
34
35#include <limits>
36
38
39// forward declarations
40#ifndef DOXYGEN
41template <typename number>
42class Vector;
43template <typename number>
44class SparseMatrix;
45namespace LinearAlgebra
46{
47 namespace distributed
48 {
49 template <typename, typename>
50 class Vector;
51 template <typename, typename>
52 class BlockVector;
53 } // namespace distributed
54} // namespace LinearAlgebra
55#endif
56
57
63namespace internal
64{
70 {
76 lanczos,
87 };
88
94 {
106 unsigned int cg_iterations;
111 unsigned int degree;
116 : min_eigenvalue_estimate{std::numeric_limits<double>::max()}
117 , max_eigenvalue_estimate{std::numeric_limits<double>::lowest()}
118 , cg_iterations{0}
119 , degree{0}
120 {}
121 };
122
128 template <typename PreconditionerType>
198} // namespace internal
199
200
221{
222public:
227
233 {
237 AdditionalData() = default;
238 };
239
244
249 template <typename MatrixType>
250 void
251 initialize(const MatrixType &matrix,
252 const AdditionalData &additional_data = AdditionalData());
253
257 template <typename VectorType>
258 void
259 vmult(VectorType &, const VectorType &) const;
260
265 template <typename VectorType>
266 void
267 Tvmult(VectorType &, const VectorType &) const;
268
272 template <typename VectorType>
273 void
274 vmult_add(VectorType &, const VectorType &) const;
275
280 template <typename VectorType>
281 void
282 Tvmult_add(VectorType &, const VectorType &) const;
283
288 void
290
299 m() const;
300
309 n() const;
310
311private:
316
321};
322
323
324
336{
337public:
342
347 {
348 public:
353 AdditionalData(const double relaxation = 1.);
354
359 };
360
366
370 void
371 initialize(const AdditionalData &parameters);
372
378 template <typename MatrixType>
379 void
380 initialize(const MatrixType &matrix, const AdditionalData &parameters);
381
385 template <typename VectorType>
386 void
387 vmult(VectorType &, const VectorType &) const;
388
393 template <typename VectorType>
394 void
395 Tvmult(VectorType &, const VectorType &) const;
399 template <typename VectorType>
400 void
401 vmult_add(VectorType &, const VectorType &) const;
402
407 template <typename VectorType>
408 void
409 Tvmult_add(VectorType &, const VectorType &) const;
410
415 void
417 {}
418
427 m() const;
428
437 n() const;
438
439private:
444
449
454};
455
456
457
497template <typename MatrixType = SparseMatrix<double>,
498 typename VectorType = Vector<double>>
500{
501public:
505 using function_ptr = void (MatrixType::*)(VectorType &,
506 const VectorType &) const;
507
513 PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
514
519 void
520 vmult(VectorType &dst, const VectorType &src) const;
521
522private:
526 const MatrixType &matrix;
527
532};
533
534
535
565template <typename MatrixType = SparseMatrix<double>,
566 typename PreconditionerType = IdentityMatrix>
568{
569public:
574
579 : public internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>
580 {
581 public:
583
587 AdditionalData(const double relaxation = 1.,
588 const unsigned int n_iterations = 1,
589 const double smoothing_range = 0.,
590 const unsigned int eig_cg_n_iterations = 8,
591 const double eig_cg_residual = 1e-2,
592 const double max_eigenvalue = 1,
594 EigenvalueAlgorithm::lanczos);
595
600
605 unsigned int n_iterations;
606 };
607
613 void
614 initialize(const MatrixType &A,
615 const AdditionalData &parameters = AdditionalData());
616
620 void
622
628 m() const;
629
635 n() const;
636
640 template <typename VectorType>
641 void
642 vmult(VectorType &, const VectorType &) const;
643
648 template <typename VectorType>
649 void
650 Tvmult(VectorType &, const VectorType &) const;
651
655 template <typename VectorType>
656 void
657 step(VectorType &x, const VectorType &rhs) const;
658
662 template <typename VectorType>
663 void
664 Tstep(VectorType &x, const VectorType &rhs) const;
665
667
671 template <typename VectorType>
673 estimate_eigenvalues(const VectorType &src) const;
674
680 double
682
683protected:
688
694
698 std::shared_ptr<PreconditionerType> preconditioner;
699
705};
706
707
708
709#ifndef DOXYGEN
710
711namespace internal
712{
713 // a helper type-trait that leverage SFINAE to figure out if MatrixType has
714 // ... MatrixType::vmult(VectorType &, const VectorType&,
715 // std::function<...>, std::function<...>) const
716 template <typename MatrixType, typename VectorType>
717 using vmult_functions_t = decltype(std::declval<const MatrixType>().vmult(
718 std::declval<VectorType &>(),
719 std::declval<const VectorType &>(),
720 std::declval<
721 const std::function<void(const unsigned int, const unsigned int)> &>(),
722 std::declval<
723 const std::function<void(const unsigned int, const unsigned int)> &>()));
724
725 template <typename MatrixType,
726 typename VectorType,
727 typename PreconditionerType>
728 constexpr bool has_vmult_with_std_functions =
729 is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
730 std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> &&
731 (std::is_same_v<VectorType,
733 std::is_same_v<
734 VectorType,
735 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
737
738
739 template <typename MatrixType, typename VectorType>
740 constexpr bool has_vmult_with_std_functions_for_precondition =
741 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
742
743 namespace PreconditionRelaxation
744 {
745 template <typename T, typename VectorType>
746 using Tvmult_t = decltype(std::declval<const T>().Tvmult(
747 std::declval<VectorType &>(),
748 std::declval<const VectorType &>()));
749
750 template <typename T, typename VectorType>
751 constexpr bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
752
753 template <typename T, typename VectorType>
754 using step_t = decltype(std::declval<const T>().step(
755 std::declval<VectorType &>(),
756 std::declval<const VectorType &>()));
757
758 template <typename T, typename VectorType>
759 constexpr bool has_step = is_supported_operation<step_t, T, VectorType>;
760
761 template <typename T, typename VectorType>
762 using step_omega_t =
763 decltype(std::declval<const T>().step(std::declval<VectorType &>(),
764 std::declval<const VectorType &>(),
765 std::declval<const double>()));
766
767 template <typename T, typename VectorType>
768 constexpr bool has_step_omega =
769 is_supported_operation<step_omega_t, T, VectorType>;
770
771 template <typename T, typename VectorType>
772 using Tstep_t = decltype(std::declval<const T>().Tstep(
773 std::declval<VectorType &>(),
774 std::declval<const VectorType &>()));
775
776 template <typename T, typename VectorType>
777 constexpr bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
778
779 template <typename T, typename VectorType>
780 using Tstep_omega_t =
781 decltype(std::declval<const T>().Tstep(std::declval<VectorType &>(),
782 std::declval<const VectorType &>(),
783 std::declval<const double>()));
784
785 template <typename T, typename VectorType>
786 constexpr bool has_Tstep_omega =
787 is_supported_operation<Tstep_omega_t, T, VectorType>;
788
789 template <typename T, typename VectorType>
790 using jacobi_step_t = decltype(std::declval<const T>().Jacobi_step(
791 std::declval<VectorType &>(),
792 std::declval<const VectorType &>(),
793 std::declval<const double>()));
794
795 template <typename T, typename VectorType>
796 constexpr bool has_jacobi_step =
797 is_supported_operation<jacobi_step_t, T, VectorType>;
798
799 template <typename T, typename VectorType>
800 using SOR_step_t = decltype(std::declval<const T>().SOR_step(
801 std::declval<VectorType &>(),
802 std::declval<const VectorType &>(),
803 std::declval<const double>()));
804
805 template <typename T, typename VectorType>
806 constexpr bool has_SOR_step =
807 is_supported_operation<SOR_step_t, T, VectorType>;
808
809 template <typename T, typename VectorType>
810 using SSOR_step_t = decltype(std::declval<const T>().SSOR_step(
811 std::declval<VectorType &>(),
812 std::declval<const VectorType &>(),
813 std::declval<const double>()));
814
815 template <typename T, typename VectorType>
816 constexpr bool has_SSOR_step =
817 is_supported_operation<SSOR_step_t, T, VectorType>;
818
819 template <typename MatrixType>
820 class PreconditionJacobiImpl
821 {
822 public:
823 PreconditionJacobiImpl(const MatrixType &A, const double relaxation)
824 : A(&A)
825 , relaxation(relaxation)
826 {}
827
828 template <typename VectorType>
829 void
830 vmult(VectorType &dst, const VectorType &src) const
831 {
832 this->A->precondition_Jacobi(dst, src, this->relaxation);
833 }
834
835 template <typename VectorType>
836 void
837 Tvmult(VectorType &dst, const VectorType &src) const
838 {
839 // call vmult, since preconditioner is symmetrical
840 this->vmult(dst, src);
841 }
842
843 template <typename VectorType,
844 std::enable_if_t<has_jacobi_step<MatrixType, VectorType>,
845 MatrixType> * = nullptr>
846 void
847 step(VectorType &dst, const VectorType &src) const
848 {
849 this->A->Jacobi_step(dst, src, this->relaxation);
850 }
851
852 template <typename VectorType,
853 std::enable_if_t<!has_jacobi_step<MatrixType, VectorType>,
854 MatrixType> * = nullptr>
855 void
856 step(VectorType &, const VectorType &) const
857 {
858 AssertThrow(false,
860 "Matrix A does not provide a Jacobi_step() function!"));
861 }
862
863 template <typename VectorType>
864 void
865 Tstep(VectorType &dst, const VectorType &src) const
866 {
867 // call step, since preconditioner is symmetrical
868 this->step(dst, src);
869 }
870
871 private:
873 const double relaxation;
874 };
875
876 template <typename MatrixType>
877 class PreconditionSORImpl
878 {
879 public:
880 PreconditionSORImpl(const MatrixType &A, const double relaxation)
881 : A(&A)
882 , relaxation(relaxation)
883 {}
884
885 template <typename VectorType>
886 void
887 vmult(VectorType &dst, const VectorType &src) const
888 {
889 this->A->precondition_SOR(dst, src, this->relaxation);
890 }
891
892 template <typename VectorType>
893 void
894 Tvmult(VectorType &dst, const VectorType &src) const
895 {
896 this->A->precondition_TSOR(dst, src, this->relaxation);
897 }
898
899 template <typename VectorType,
900 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
901 MatrixType> * = nullptr>
902 void
903 step(VectorType &dst, const VectorType &src) const
904 {
905 this->A->SOR_step(dst, src, this->relaxation);
906 }
907
908 template <typename VectorType,
909 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
910 MatrixType> * = nullptr>
911 void
912 step(VectorType &, const VectorType &) const
913 {
914 AssertThrow(false,
916 "Matrix A does not provide a SOR_step() function!"));
917 }
918
919 template <typename VectorType,
920 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
921 MatrixType> * = nullptr>
922 void
923 Tstep(VectorType &dst, const VectorType &src) const
924 {
925 this->A->TSOR_step(dst, src, this->relaxation);
926 }
927
928 template <typename VectorType,
929 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
930 MatrixType> * = nullptr>
931 void
932 Tstep(VectorType &, const VectorType &) const
933 {
934 AssertThrow(false,
936 "Matrix A does not provide a TSOR_step() function!"));
937 }
938
939 private:
941 const double relaxation;
942 };
943
944 template <typename MatrixType>
945 class PreconditionSSORImpl
946 {
947 public:
948 using size_type = typename MatrixType::size_type;
949
950 PreconditionSSORImpl(const MatrixType &A, const double relaxation)
951 : A(&A)
952 , relaxation(relaxation)
953 {
954 // in case we have a SparseMatrix class, we can extract information
955 // about the diagonal.
958 &*this->A);
959
960 // calculate the positions first after the diagonal.
961 if (mat != nullptr)
962 {
963 const size_type n = this->A->n();
964 pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
965 for (size_type row = 0; row < n; ++row)
966 {
967 // find the first element in this line which is on the right of
968 // the diagonal. we need to precondition with the elements on
969 // the left only. note: the first entry in each line denotes the
970 // diagonal element, which we need not check.
971 typename SparseMatrix<
972 typename MatrixType::value_type>::const_iterator it =
973 mat->begin(row) + 1;
974 for (; it < mat->end(row); ++it)
975 if (it->column() > row)
976 break;
977 pos_right_of_diagonal[row] = it - mat->begin();
978 }
979 }
980 }
981
982 template <typename VectorType>
983 void
984 vmult(VectorType &dst, const VectorType &src) const
985 {
986 this->A->precondition_SSOR(dst,
987 src,
988 this->relaxation,
989 pos_right_of_diagonal);
990 }
991
992 template <typename VectorType>
993 void
994 Tvmult(VectorType &dst, const VectorType &src) const
995 {
996 this->A->precondition_SSOR(dst,
997 src,
998 this->relaxation,
999 pos_right_of_diagonal);
1000 }
1001
1002 template <typename VectorType,
1003 std::enable_if_t<has_SSOR_step<MatrixType, VectorType>,
1004 MatrixType> * = nullptr>
1005 void
1006 step(VectorType &dst, const VectorType &src) const
1007 {
1008 this->A->SSOR_step(dst, src, this->relaxation);
1009 }
1010
1011 template <typename VectorType,
1012 std::enable_if_t<!has_SSOR_step<MatrixType, VectorType>,
1013 MatrixType> * = nullptr>
1014 void
1015 step(VectorType &, const VectorType &) const
1016 {
1017 AssertThrow(false,
1018 ExcMessage(
1019 "Matrix A does not provide a SSOR_step() function!"));
1020 }
1021
1022 template <typename VectorType>
1023 void
1024 Tstep(VectorType &dst, const VectorType &src) const
1025 {
1026 // call step, since preconditioner is symmetrical
1027 this->step(dst, src);
1028 }
1029
1030 private:
1032 const double relaxation;
1033
1038 std::vector<std::size_t> pos_right_of_diagonal;
1039 };
1040
1041 template <typename MatrixType>
1042 class PreconditionPSORImpl
1043 {
1044 public:
1045 using size_type = typename MatrixType::size_type;
1046
1047 PreconditionPSORImpl(const MatrixType &A,
1048 const double relaxation,
1049 const std::vector<size_type> &permutation,
1050 const std::vector<size_type> &inverse_permutation)
1051 : A(&A)
1052 , relaxation(relaxation)
1053 , permutation(permutation)
1054 , inverse_permutation(inverse_permutation)
1055 {}
1056
1057 template <typename VectorType>
1058 void
1059 vmult(VectorType &dst, const VectorType &src) const
1060 {
1061 dst = src;
1062 this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
1063 }
1064
1065 template <typename VectorType>
1066 void
1067 Tvmult(VectorType &dst, const VectorType &src) const
1068 {
1069 dst = src;
1070 this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
1071 }
1072
1073 private:
1075 const double relaxation;
1076
1077 const std::vector<size_type> &permutation;
1078 const std::vector<size_type> &inverse_permutation;
1079 };
1080
1081 template <typename MatrixType,
1082 typename PreconditionerType,
1083 typename VectorType,
1084 std::enable_if_t<has_step_omega<PreconditionerType, VectorType>,
1085 PreconditionerType> * = nullptr>
1086 void
1087 step(const MatrixType &,
1088 const PreconditionerType &preconditioner,
1089 VectorType &dst,
1090 const VectorType &src,
1091 const double relaxation,
1092 VectorType &,
1093 VectorType &)
1094 {
1095 preconditioner.step(dst, src, relaxation);
1096 }
1097
1098 template <
1099 typename MatrixType,
1100 typename PreconditionerType,
1101 typename VectorType,
1102 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1103 has_step<PreconditionerType, VectorType>,
1104 PreconditionerType> * = nullptr>
1105 void
1106 step(const MatrixType &,
1107 const PreconditionerType &preconditioner,
1108 VectorType &dst,
1109 const VectorType &src,
1110 const double relaxation,
1111 VectorType &,
1112 VectorType &)
1113 {
1114 Assert(relaxation == 1.0, ExcInternalError());
1115
1116 (void)relaxation;
1117
1118 preconditioner.step(dst, src);
1119 }
1120
1121 template <
1122 typename MatrixType,
1123 typename PreconditionerType,
1124 typename VectorType,
1125 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1126 !has_step<PreconditionerType, VectorType>,
1127 PreconditionerType> * = nullptr>
1128 void
1129 step(const MatrixType &A,
1130 const PreconditionerType &preconditioner,
1131 VectorType &dst,
1132 const VectorType &src,
1133 const double relaxation,
1134 VectorType &residual,
1135 VectorType &tmp)
1136 {
1137 residual.reinit(dst, true);
1138 tmp.reinit(dst, true);
1139
1140 A.vmult(residual, dst);
1141 residual.sadd(-1.0, 1.0, src);
1142
1143 preconditioner.vmult(tmp, residual);
1144 dst.add(relaxation, tmp);
1145 }
1146
1147 template <typename MatrixType,
1148 typename PreconditionerType,
1149 typename VectorType,
1150 std::enable_if_t<has_Tstep_omega<PreconditionerType, VectorType>,
1151 PreconditionerType> * = nullptr>
1152 void
1153 Tstep(const MatrixType &,
1154 const PreconditionerType &preconditioner,
1155 VectorType &dst,
1156 const VectorType &src,
1157 const double relaxation,
1158 VectorType &,
1159 VectorType &)
1160 {
1161 preconditioner.Tstep(dst, src, relaxation);
1162 }
1163
1164 template <
1165 typename MatrixType,
1166 typename PreconditionerType,
1167 typename VectorType,
1168 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1169 has_Tstep<PreconditionerType, VectorType>,
1170 PreconditionerType> * = nullptr>
1171 void
1172 Tstep(const MatrixType &,
1173 const PreconditionerType &preconditioner,
1174 VectorType &dst,
1175 const VectorType &src,
1176 const double relaxation,
1177 VectorType &,
1178 VectorType &)
1179 {
1180 Assert(relaxation == 1.0, ExcInternalError());
1181
1182 (void)relaxation;
1183
1184 preconditioner.Tstep(dst, src);
1185 }
1186
1187 template <typename MatrixType,
1188 typename VectorType,
1189 std::enable_if_t<has_Tvmult<MatrixType, VectorType>, MatrixType>
1190 * = nullptr>
1191 void
1192 Tvmult(const MatrixType &A, VectorType &dst, const VectorType &src)
1193 {
1194 A.Tvmult(dst, src);
1195 }
1196
1197 template <typename MatrixType,
1198 typename VectorType,
1199 std::enable_if_t<!has_Tvmult<MatrixType, VectorType>, MatrixType>
1200 * = nullptr>
1201 void
1202 Tvmult(const MatrixType &, VectorType &, const VectorType &)
1203 {
1204 AssertThrow(false,
1205 ExcMessage("Matrix A does not provide a Tvmult() function!"));
1206 }
1207
1208 template <
1209 typename MatrixType,
1210 typename PreconditionerType,
1211 typename VectorType,
1212 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1213 !has_Tstep<PreconditionerType, VectorType>,
1214 PreconditionerType> * = nullptr>
1215 void
1216 Tstep(const MatrixType &A,
1217 const PreconditionerType &preconditioner,
1218 VectorType &dst,
1219 const VectorType &src,
1220 const double relaxation,
1221 VectorType &residual,
1222 VectorType &tmp)
1223 {
1224 residual.reinit(dst, true);
1225 tmp.reinit(dst, true);
1226
1227 Tvmult(A, residual, dst);
1228 residual.sadd(-1.0, 1.0, src);
1229
1230 Tvmult(preconditioner, tmp, residual);
1231 dst.add(relaxation, tmp);
1232 }
1233
1234 // 0) general implementation
1235 template <typename MatrixType,
1236 typename PreconditionerType,
1237 typename VectorType,
1238 std::enable_if_t<!has_vmult_with_std_functions_for_precondition<
1239 PreconditionerType,
1240 VectorType>,
1241 int> * = nullptr>
1242 void
1243 step_operations(const MatrixType &A,
1244 const PreconditionerType &preconditioner,
1245 VectorType &dst,
1246 const VectorType &src,
1247 const double relaxation,
1248 VectorType &tmp1,
1249 VectorType &tmp2,
1250 const unsigned int i,
1251 const bool transposed)
1252 {
1253 if (i == 0)
1254 {
1255 if (transposed)
1256 Tvmult(preconditioner, dst, src);
1257 else
1258 preconditioner.vmult(dst, src);
1259
1260 if (relaxation != 1.0)
1261 dst *= relaxation;
1262 }
1263 else
1264 {
1265 if (transposed)
1266 Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1267 else
1268 step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1269 }
1270 }
1271
1272 // 1) specialized implementation with a preconditioner that accepts
1273 // ranges
1274 template <
1275 typename MatrixType,
1276 typename PreconditionerType,
1277 typename VectorType,
1278 std::enable_if_t<
1279 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1280 VectorType> &&
1281 !has_vmult_with_std_functions_for_precondition<MatrixType,
1282 VectorType>,
1283 int> * = nullptr>
1284 void
1285 step_operations(const MatrixType &A,
1286 const PreconditionerType &preconditioner,
1287 VectorType &dst,
1288 const VectorType &src,
1289 const double relaxation,
1290 VectorType &tmp,
1291 VectorType &,
1292 const unsigned int i,
1293 const bool transposed)
1294 {
1295 (void)transposed;
1296 using Number = typename VectorType::value_type;
1297 Number *dst_ptr = dst.begin();
1298 const Number *src_ptr = src.begin();
1299
1300 if (i == 0)
1301 {
1302 preconditioner.vmult(
1303 dst,
1304 src,
1305 [&](const unsigned int start_range, const unsigned int end_range) {
1306 // zero 'dst' before running the vmult operation
1307 if (end_range > start_range)
1308 std::memset(dst_ptr + start_range,
1309 0,
1310 sizeof(Number) * (end_range - start_range));
1311 },
1312 [&](const unsigned int start_range, const unsigned int end_range) {
1313 if (relaxation == 1.0)
1314 return; // nothing to do
1315
1317 for (std::size_t i = start_range; i < end_range; ++i)
1318 dst_ptr[i] *= relaxation;
1319 });
1320 }
1321 else
1322 {
1323 tmp.reinit(src, true);
1324
1325 Assert(transposed == false, ExcNotImplemented());
1326
1327 A.vmult(tmp, dst);
1328
1329 preconditioner.vmult(
1330 dst,
1331 tmp,
1332 [&](const unsigned int start_range, const unsigned int end_range) {
1333 const auto tmp_ptr = tmp.begin();
1334
1335 if (relaxation == 1.0)
1336 {
1338 for (std::size_t i = start_range; i < end_range; ++i)
1339 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1340 }
1341 else
1342 {
1343 // note: we scale the residual here to be able to add into
1344 // the dst vector, which contains the solution from the last
1345 // iteration
1347 for (std::size_t i = start_range; i < end_range; ++i)
1348 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1349 }
1350 },
1351 [&](const unsigned int, const unsigned int) {
1352 // nothing to do, since scaling by the relaxation factor
1353 // has been done in the pre operation
1354 });
1355 }
1356 }
1357
1358 // 2) specialized implementation with a preconditioner and a matrix that
1359 // accepts ranges
1360 template <
1361 typename MatrixType,
1362 typename PreconditionerType,
1363 typename VectorType,
1364 std::enable_if_t<
1365 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1366 VectorType> &&
1367 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1368 int> * = nullptr>
1369 void
1370 step_operations(const MatrixType &A,
1371 const PreconditionerType &preconditioner,
1372 VectorType &dst,
1373 const VectorType &src,
1374 const double relaxation,
1375 VectorType &tmp,
1376 VectorType &,
1377 const unsigned int i,
1378 const bool transposed)
1379 {
1380 (void)transposed;
1381 using Number = typename VectorType::value_type;
1382
1383 Number *dst_ptr = dst.begin();
1384 const Number *src_ptr = src.begin();
1385
1386 if (i == 0)
1387 {
1388 preconditioner.vmult(
1389 dst,
1390 src,
1391 [&](const unsigned int start_range, const unsigned int end_range) {
1392 // zero 'dst' before running the vmult operation
1393 if (end_range > start_range)
1394 std::memset(dst_ptr + start_range,
1395 0,
1396 sizeof(Number) * (end_range - start_range));
1397 },
1398 [&](const unsigned int start_range, const unsigned int end_range) {
1399 if (relaxation == 1.0)
1400 return; // nothing to do
1401
1403 for (std::size_t i = start_range; i < end_range; ++i)
1404 dst_ptr[i] *= relaxation;
1405 });
1406 }
1407 else
1408 {
1409 tmp.reinit(src, true);
1410 const auto tmp_ptr = tmp.begin();
1411
1412 Assert(transposed == false, ExcNotImplemented());
1413
1414 A.vmult(
1415 tmp,
1416 dst,
1417 [&](const unsigned int start_range, const unsigned int end_range) {
1418 // zero 'tmp' before running the vmult
1419 // operation
1420 if (end_range > start_range)
1421 std::memset(tmp_ptr + start_range,
1422 0,
1423 sizeof(Number) * (end_range - start_range));
1424 },
1425 [&](const unsigned int start_range, const unsigned int end_range) {
1426 if (relaxation == 1.0)
1427 {
1429 for (std::size_t i = start_range; i < end_range; ++i)
1430 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1431 }
1432 else
1433 {
1434 // note: we scale the residual here to be able to add into
1435 // the dst vector, which contains the solution from the last
1436 // iteration
1438 for (std::size_t i = start_range; i < end_range; ++i)
1439 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1440 }
1441 });
1442
1443 preconditioner.vmult(dst, tmp, [](const auto, const auto) {
1444 // note: `dst` vector does not have to be zeroed out
1445 // since we add the result into it
1446 });
1447 }
1448 }
1449
1450 // 3) specialized implementation for inverse-diagonal preconditioner
1451 template <
1452 typename MatrixType,
1453 typename VectorType,
1454 std::enable_if_t<
1456 !std::is_same_v<
1457 VectorType,
1458 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
1460 !has_vmult_with_std_functions<MatrixType,
1461 VectorType,
1463 VectorType> * = nullptr>
1464 void
1465 step_operations(const MatrixType &A,
1466 const ::DiagonalMatrix<VectorType> &preconditioner,
1467 VectorType &dst,
1468 const VectorType &src,
1469 const double relaxation,
1470 VectorType &tmp,
1471 VectorType &,
1472 const unsigned int i,
1473 const bool transposed)
1474 {
1475 using Number = typename VectorType::value_type;
1476
1477 if (i == 0)
1478 {
1479 Number *dst_ptr = dst.begin();
1480 const Number *src_ptr = src.begin();
1481 const Number *diag_ptr = preconditioner.get_vector().begin();
1482
1483 if (relaxation == 1.0)
1484 {
1486 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1487 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1488 }
1489 else
1490 {
1492 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1493 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1494 }
1495 }
1496 else
1497 {
1498 tmp.reinit(src, true);
1499
1500 if (transposed)
1501 Tvmult(A, tmp, dst);
1502 else
1503 A.vmult(tmp, dst);
1504
1505 Number *dst_ptr = dst.begin();
1506 const Number *src_ptr = src.begin();
1507 const Number *tmp_ptr = tmp.begin();
1508 const Number *diag_ptr = preconditioner.get_vector().begin();
1509
1510 if (relaxation == 1.0)
1511 {
1513 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1514 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1515 }
1516 else
1517 {
1519 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1520 dst_ptr[i] +=
1521 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1522 }
1523 }
1524 }
1525
1526 // 4) specialized implementation for inverse-diagonal preconditioner and
1527 // matrix that accepts ranges
1528 template <typename MatrixType,
1529 typename VectorType,
1530 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1531 has_vmult_with_std_functions<
1532 MatrixType,
1533 VectorType,
1535 VectorType> * = nullptr>
1536 void
1537 step_operations(const MatrixType &A,
1538 const ::DiagonalMatrix<VectorType> &preconditioner,
1539 VectorType &dst,
1540 const VectorType &src,
1541 const double relaxation,
1542 VectorType &tmp,
1543 VectorType &,
1544 const unsigned int i,
1545 const bool transposed)
1546 {
1547 (void)transposed;
1548 using Number = typename VectorType::value_type;
1549
1550 if (i == 0)
1551 {
1552 Number *dst_ptr = dst.begin();
1553 const Number *src_ptr = src.begin();
1554 const Number *diag_ptr = preconditioner.get_vector().begin();
1555
1556 if (relaxation == 1.0)
1557 {
1559 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1560 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1561 }
1562 else
1563 {
1565 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1566 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1567 }
1568 }
1569 else
1570 {
1571 tmp.reinit(src, true);
1572
1573 Assert(transposed == false, ExcNotImplemented());
1574
1575 A.vmult(
1576 tmp,
1577 dst,
1578 [&](const unsigned int start_range, const unsigned int end_range) {
1579 // zero 'tmp' before running the vmult operation
1580 if (end_range > start_range)
1581 std::memset(tmp.begin() + start_range,
1582 0,
1583 sizeof(Number) * (end_range - start_range));
1584 },
1585 [&](const unsigned int begin, const unsigned int end) {
1586 const Number *dst_ptr = dst.begin();
1587 const Number *src_ptr = src.begin();
1588 Number *tmp_ptr = tmp.begin();
1589 const Number *diag_ptr = preconditioner.get_vector().begin();
1590
1591 // for efficiency reason, write back to temp_vector that is
1592 // already read (avoid read-for-ownership)
1593 if (relaxation == 1.0)
1594 {
1596 for (std::size_t i = begin; i < end; ++i)
1597 tmp_ptr[i] =
1598 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1599 }
1600 else
1601 {
1603 for (std::size_t i = begin; i < end; ++i)
1604 tmp_ptr[i] = dst_ptr[i] + relaxation *
1605 (src_ptr[i] - tmp_ptr[i]) *
1606 diag_ptr[i];
1607 }
1608 });
1609
1610 tmp.swap(dst);
1611 }
1612 }
1613
1614 } // namespace PreconditionRelaxation
1615} // namespace internal
1616
1617#endif
1618
1619
1620
1647template <typename MatrixType = SparseMatrix<double>>
1649 : public PreconditionRelaxation<
1650 MatrixType,
1651 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1652{
1654 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1656
1657public:
1662
1666 void
1667 initialize(const MatrixType &A,
1668 const AdditionalData &parameters = AdditionalData());
1669};
1670
1671
1717template <typename MatrixType = SparseMatrix<double>>
1719 : public PreconditionRelaxation<
1720 MatrixType,
1721 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1722{
1724 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1726
1727public:
1732
1736 void
1737 initialize(const MatrixType &A,
1738 const AdditionalData &parameters = AdditionalData());
1739};
1740
1741
1742
1769template <typename MatrixType = SparseMatrix<double>>
1771 : public PreconditionRelaxation<
1772 MatrixType,
1773 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1774{
1776 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1778
1779public:
1784
1790 void
1791 initialize(const MatrixType &A,
1792 const AdditionalData &parameters = AdditionalData());
1793};
1794
1795
1825template <typename MatrixType = SparseMatrix<double>>
1827 : public PreconditionRelaxation<
1828 MatrixType,
1829 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1830{
1832 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1834
1835public:
1840
1845 {
1846 public:
1857 AdditionalData(const std::vector<size_type> &permutation,
1858 const std::vector<size_type> &inverse_permutation,
1859 const typename BaseClass::AdditionalData &parameters =
1860 typename BaseClass::AdditionalData());
1861
1865 const std::vector<size_type> &permutation;
1869 const std::vector<size_type> &inverse_permutation;
1874 };
1875
1887 void
1888 initialize(const MatrixType &A,
1889 const std::vector<size_type> &permutation,
1890 const std::vector<size_type> &inverse_permutation,
1891 const typename BaseClass::AdditionalData &parameters =
1892 typename BaseClass::AdditionalData());
1893
1904 void
1905 initialize(const MatrixType &A, const AdditionalData &additional_data);
1906};
1907
1908
1909
2098template <typename MatrixType = SparseMatrix<double>,
2099 typename VectorType = Vector<double>,
2100 typename PreconditionerType = DiagonalMatrix<VectorType>>
2102{
2103public:
2108
2114 : public internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>
2115 {
2117
2122 {
2126 first_kind,
2132 };
2133
2138 const unsigned int degree = 1,
2139 const double smoothing_range = 0.,
2140 const unsigned int eig_cg_n_iterations = 8,
2141 const double eig_cg_residual = 1e-2,
2142 const double max_eigenvalue = 1,
2144 EigenvalueAlgorithm::lanczos,
2146
2159 unsigned int degree;
2160
2165 };
2166
2167
2172
2184 void
2185 initialize(const MatrixType &matrix,
2186 const AdditionalData &additional_data = AdditionalData());
2187
2192 void
2193 vmult(VectorType &dst, const VectorType &src) const;
2194
2199 void
2200 Tvmult(VectorType &dst, const VectorType &src) const;
2201
2205 void
2206 step(VectorType &dst, const VectorType &src) const;
2207
2211 void
2212 Tstep(VectorType &dst, const VectorType &src) const;
2213
2217 void
2219
2224 size_type
2225 m() const;
2226
2231 size_type
2232 n() const;
2233
2235
2249 estimate_eigenvalues(const VectorType &src) const;
2250
2251private:
2256 const MatrixType,
2259
2263 mutable VectorType solution_old;
2264
2268 mutable VectorType temp_vector1;
2269
2273 mutable VectorType temp_vector2;
2274
2280
2284 double theta;
2285
2290 double delta;
2291
2297
2303};
2304
2305
2306
2308/* ---------------------------------- Inline functions ------------------- */
2309
2310#ifndef DOXYGEN
2311
2312
2313namespace internal
2314{
2315 template <typename VectorType>
2316 void
2317 set_initial_guess(VectorType &vector)
2318 {
2319 vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2320 if (vector.locally_owned_elements().is_element(0))
2321 vector(0) = 0.;
2322 }
2323
2324 template <typename Number>
2325 void
2326 set_initial_guess(::Vector<Number> &vector)
2327 {
2328 // Choose a high-frequency mode consisting of numbers between 0 and 1
2329 // that is cheap to compute (cheaper than random numbers) but avoids
2330 // obviously re-occurring numbers in multi-component systems by choosing
2331 // a period of 11
2332 for (unsigned int i = 0; i < vector.size(); ++i)
2333 vector(i) = i % 11;
2334
2335 const Number mean_value = vector.mean_value();
2336 vector.add(-mean_value);
2337 }
2338
2339 template <typename Number, typename MemorySpace>
2340 void
2341 set_initial_guess(
2343 &vector)
2344 {
2345 for (unsigned int block = 0; block < vector.n_blocks(); ++block)
2346 set_initial_guess(vector.block(block));
2347 }
2348
2349 template <typename Number, typename MemorySpace>
2350 void
2351 set_initial_guess(
2353 {
2354 // Choose a high-frequency mode consisting of numbers between 0 and 1
2355 // that is cheap to compute (cheaper than random numbers) but avoids
2356 // obviously re-occurring numbers in multi-component systems by choosing
2357 // a period of 11.
2358 // Make initial guess robust with respect to number of processors
2359 // by operating on the global index.
2360 types::global_dof_index first_local_range = 0;
2361 if (!vector.locally_owned_elements().is_empty())
2362 first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2363
2364 const auto n_local_elements = vector.locally_owned_size();
2365 Number *values_ptr = vector.get_values();
2366 Kokkos::RangePolicy<typename MemorySpace::kokkos_space::execution_space,
2367 Kokkos::IndexType<types::global_dof_index>>
2368 policy(0, n_local_elements);
2369 Kokkos::parallel_for(
2370 "::PreconditionChebyshev::set_initial_guess",
2371 policy,
2372 KOKKOS_LAMBDA(types::global_dof_index i) {
2373 values_ptr[i] = (i + first_local_range) % 11;
2374 });
2375 const Number mean_value = vector.mean_value();
2376 vector.add(-mean_value);
2377 }
2378
2379 struct EigenvalueTracker
2380 {
2381 public:
2382 void
2383 slot(const std::vector<double> &eigenvalues)
2384 {
2386 }
2387
2388 std::vector<double> values;
2389 };
2390
2391
2392
2393 template <typename MatrixType,
2394 typename VectorType,
2395 typename PreconditionerType>
2396 double
2397 power_iteration(const MatrixType &matrix,
2398 VectorType &eigenvector,
2399 const PreconditionerType &preconditioner,
2400 const unsigned int n_iterations)
2401 {
2402 typename VectorType::value_type eigenvalue_estimate = 0.;
2403 eigenvector /= eigenvector.l2_norm();
2404 VectorType vector1, vector2;
2405 vector1.reinit(eigenvector, true);
2406 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2407 vector2.reinit(eigenvector, true);
2408 for (unsigned int i = 0; i < n_iterations; ++i)
2409 {
2410 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2411 {
2412 matrix.vmult(vector2, eigenvector);
2413 preconditioner.vmult(vector1, vector2);
2414 }
2415 else
2416 matrix.vmult(vector1, eigenvector);
2417
2418 eigenvalue_estimate = eigenvector * vector1;
2419
2420 vector1 /= vector1.l2_norm();
2421 eigenvector.swap(vector1);
2422 }
2423 return std::abs(eigenvalue_estimate);
2424 }
2425
2426
2427
2428 template <typename MatrixType,
2429 typename VectorType,
2430 typename PreconditionerType>
2431 EigenvalueInformation
2432 estimate_eigenvalues(
2433 const EigenvalueAlgorithmAdditionalData<PreconditionerType> &data,
2434 const MatrixType *matrix_ptr,
2435 VectorType &solution_old,
2436 VectorType &temp_vector1,
2437 const unsigned int degree)
2438 {
2439 Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2440
2441 EigenvalueInformation info{};
2442
2443 if (data.eig_cg_n_iterations > 0)
2444 {
2445 Assert(data.eig_cg_n_iterations > 2,
2446 ExcMessage(
2447 "Need to set at least two iterations to find eigenvalues."));
2448
2449 internal::EigenvalueTracker eigenvalue_tracker;
2450
2451 // set an initial guess that contains some high-frequency parts (to the
2452 // extent possible without knowing the discretization and the numbering)
2453 // to trigger high eigenvalues according to the external function
2454 internal::set_initial_guess(temp_vector1);
2455 data.constraints.set_zero(temp_vector1);
2456
2457 if (data.eigenvalue_algorithm == internal::EigenvalueAlgorithm::lanczos)
2458 {
2459 // set a very strict tolerance to force at least two iterations
2460 IterationNumberControl control(data.eig_cg_n_iterations,
2461 1e-10,
2462 false,
2463 false);
2464
2465 ::SolverCG<VectorType> solver(control);
2466 solver.connect_eigenvalues_slot(
2467 [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2468 eigenvalue_tracker.slot(eigenvalues);
2469 });
2470
2471 solver.solve(*matrix_ptr,
2472 solution_old,
2473 temp_vector1,
2474 *data.preconditioner);
2475
2476 info.cg_iterations = control.last_step();
2477 }
2478 else if (data.eigenvalue_algorithm ==
2480 {
2481 (void)degree;
2482
2484 ExcMessage("Cannot estimate the minimal eigenvalue with the "
2485 "power iteration"));
2486
2487 eigenvalue_tracker.values.push_back(
2488 internal::power_iteration(*matrix_ptr,
2489 temp_vector1,
2490 *data.preconditioner,
2491 data.eig_cg_n_iterations));
2492 }
2493 else
2495
2496 // read the eigenvalues from the attached eigenvalue tracker
2497 if (eigenvalue_tracker.values.empty())
2498 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2499 else
2500 {
2501 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2502
2503 // include a safety factor since the CG method will in general not
2504 // be converged
2505 info.max_eigenvalue_estimate =
2506 1.2 * eigenvalue_tracker.values.back();
2507 }
2508 }
2509 else
2510 {
2511 info.max_eigenvalue_estimate = data.max_eigenvalue;
2512 info.min_eigenvalue_estimate =
2513 data.max_eigenvalue / data.smoothing_range;
2514 }
2515
2516 return info;
2517 }
2518} // namespace internal
2519
2520
2521
2523 : n_rows(0)
2524 , n_columns(0)
2525{}
2526
2527template <typename MatrixType>
2528inline void
2529PreconditionIdentity::initialize(const MatrixType &matrix,
2531{
2532 n_rows = matrix.m();
2533 n_columns = matrix.n();
2534}
2535
2536
2537template <typename VectorType>
2538inline void
2539PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
2540{
2541 dst = src;
2542}
2543
2544
2545
2546template <typename VectorType>
2547inline void
2548PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
2549{
2550 dst = src;
2551}
2552
2553template <typename VectorType>
2554inline void
2555PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
2556{
2557 dst += src;
2558}
2559
2560
2561
2562template <typename VectorType>
2563inline void
2564PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
2565{
2566 dst += src;
2567}
2568
2569
2570
2571inline void
2573{}
2574
2575
2576
2579{
2581 return n_rows;
2582}
2583
2586{
2588 return n_columns;
2589}
2590
2591//---------------------------------------------------------------------------
2592
2594 const double relaxation)
2595 : relaxation(relaxation)
2596{}
2597
2598
2600 : relaxation(0)
2601 , n_rows(0)
2602 , n_columns(0)
2603{
2604 AdditionalData add_data;
2605 relaxation = add_data.relaxation;
2606}
2607
2608
2609
2610inline void
2613{
2614 relaxation = parameters.relaxation;
2615}
2616
2617
2618
2619template <typename MatrixType>
2620inline void
2622 const MatrixType &matrix,
2624{
2625 relaxation = parameters.relaxation;
2626 n_rows = matrix.m();
2627 n_columns = matrix.n();
2628}
2629
2630
2631
2632template <typename VectorType>
2633inline void
2634PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2635{
2636 static_assert(
2637 std::is_same_v<size_type, typename VectorType::size_type>,
2638 "PreconditionRichardson and VectorType must have the same size_type.");
2639
2640 dst.equ(relaxation, src);
2641}
2642
2643
2644
2645template <typename VectorType>
2646inline void
2647PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2648{
2649 static_assert(
2650 std::is_same_v<size_type, typename VectorType::size_type>,
2651 "PreconditionRichardson and VectorType must have the same size_type.");
2652
2653 dst.equ(relaxation, src);
2654}
2655
2656template <typename VectorType>
2657inline void
2658PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2659{
2660 static_assert(
2661 std::is_same_v<size_type, typename VectorType::size_type>,
2662 "PreconditionRichardson and VectorType must have the same size_type.");
2663
2664 dst.add(relaxation, src);
2665}
2666
2667
2668
2669template <typename VectorType>
2670inline void
2671PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2672{
2673 static_assert(
2674 std::is_same_v<size_type, typename VectorType::size_type>,
2675 "PreconditionRichardson and VectorType must have the same size_type.");
2676
2677 dst.add(relaxation, src);
2678}
2679
2682{
2684 return n_rows;
2685}
2686
2689{
2691 return n_columns;
2692}
2693
2694//---------------------------------------------------------------------------
2695
2696template <typename MatrixType, typename PreconditionerType>
2697inline void
2699 const MatrixType &rA,
2700 const AdditionalData &parameters)
2701{
2702 A = &rA;
2703 eigenvalues_are_initialized = false;
2704
2705 Assert(parameters.preconditioner, ExcNotInitialized());
2706
2707 this->data = parameters;
2708}
2709
2710
2711template <typename MatrixType, typename PreconditionerType>
2712inline void
2714{
2715 eigenvalues_are_initialized = false;
2716 A = nullptr;
2717 data.relaxation = 1.0;
2718 data.preconditioner = nullptr;
2719}
2720
2721template <typename MatrixType, typename PreconditionerType>
2722inline
2725{
2726 Assert(A != nullptr, ExcNotInitialized());
2727 return A->m();
2728}
2729
2730template <typename MatrixType, typename PreconditionerType>
2731inline
2734{
2735 Assert(A != nullptr, ExcNotInitialized());
2736 return A->n();
2737}
2738
2739template <typename MatrixType, typename PreconditionerType>
2740template <typename VectorType>
2741inline void
2743 VectorType &dst,
2744 const VectorType &src) const
2745{
2746 Assert(this->A != nullptr, ExcNotInitialized());
2747 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2748
2749 if (eigenvalues_are_initialized == false)
2750 estimate_eigenvalues(src);
2751
2752 VectorType tmp1, tmp2;
2753
2754 for (unsigned int i = 0; i < data.n_iterations; ++i)
2755 internal::PreconditionRelaxation::step_operations(*A,
2756 *data.preconditioner,
2757 dst,
2758 src,
2759 data.relaxation,
2760 tmp1,
2761 tmp2,
2762 i,
2763 false);
2764}
2765
2766template <typename MatrixType, typename PreconditionerType>
2767template <typename VectorType>
2768inline void
2770 VectorType &dst,
2771 const VectorType &src) const
2772{
2773 Assert(this->A != nullptr, ExcNotInitialized());
2774 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2775
2776 if (eigenvalues_are_initialized == false)
2777 estimate_eigenvalues(src);
2778
2779 VectorType tmp1, tmp2;
2780
2781 for (unsigned int i = 0; i < data.n_iterations; ++i)
2782 internal::PreconditionRelaxation::step_operations(
2783 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2784}
2785
2786template <typename MatrixType, typename PreconditionerType>
2787template <typename VectorType>
2788inline void
2790 VectorType &dst,
2791 const VectorType &src) const
2792{
2793 Assert(this->A != nullptr, ExcNotInitialized());
2794 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2795
2796 if (eigenvalues_are_initialized == false)
2797 estimate_eigenvalues(src);
2798
2799 VectorType tmp1, tmp2;
2800
2801 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2802 internal::PreconditionRelaxation::step_operations(*A,
2803 *data.preconditioner,
2804 dst,
2805 src,
2806 data.relaxation,
2807 tmp1,
2808 tmp2,
2809 i,
2810 false);
2811}
2812
2813template <typename MatrixType, typename PreconditionerType>
2814template <typename VectorType>
2815inline void
2817 VectorType &dst,
2818 const VectorType &src) const
2819{
2820 Assert(this->A != nullptr, ExcNotInitialized());
2821 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2822
2823 if (eigenvalues_are_initialized == false)
2824 estimate_eigenvalues(src);
2825
2826 VectorType tmp1, tmp2;
2827
2828 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2829 internal::PreconditionRelaxation::step_operations(
2830 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2831}
2832
2833template <typename MatrixType, typename PreconditionerType>
2834template <typename VectorType>
2837 const VectorType &src) const
2838{
2839 Assert(eigenvalues_are_initialized == false, ExcInternalError());
2840
2841 EigenvalueInformation info;
2842
2843 if (data.relaxation == 0.0)
2844 {
2845 VectorType solution_old, temp_vector1;
2846
2847 solution_old.reinit(src);
2848 temp_vector1.reinit(src, true);
2849
2850 info = internal::estimate_eigenvalues<MatrixType>(
2851 data, A, solution_old, temp_vector1, data.n_iterations);
2852
2853 const double alpha =
2854 (data.smoothing_range > 1. ?
2855 info.max_eigenvalue_estimate / data.smoothing_range :
2856 std::min(0.9 * info.max_eigenvalue_estimate,
2857 info.min_eigenvalue_estimate));
2858
2860 ->data.relaxation = 2.0 / (alpha + info.max_eigenvalue_estimate);
2861 }
2862
2864 ->eigenvalues_are_initialized = true;
2865
2866 return info;
2867}
2868
2869template <typename MatrixType, typename PreconditionerType>
2870double
2872{
2873 return data.relaxation;
2874}
2875
2876
2877//---------------------------------------------------------------------------
2878
2879template <typename MatrixType>
2880inline void
2882 const AdditionalData &parameters_in)
2883{
2884 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2885 Assert(
2886 parameters_in.relaxation != 0.0,
2887 ExcMessage(
2888 "Relaxation cannot automatically be determined by PreconditionJacobi."));
2889
2890 AdditionalData parameters;
2891 parameters.relaxation = 1.0;
2892 parameters.n_iterations = parameters_in.n_iterations;
2893 parameters.preconditioner =
2894 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2895
2896 this->BaseClass::initialize(A, parameters);
2897}
2898
2899//---------------------------------------------------------------------------
2900
2901template <typename MatrixType>
2902inline void
2903PreconditionSOR<MatrixType>::initialize(const MatrixType &A,
2904 const AdditionalData &parameters_in)
2905{
2906 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2907 Assert(
2908 parameters_in.relaxation != 0.0,
2909 ExcMessage(
2910 "Relaxation cannot automatically be determined by PreconditionSOR."));
2911
2912 AdditionalData parameters;
2913 parameters.relaxation = 1.0;
2914 parameters.n_iterations = parameters_in.n_iterations;
2915 parameters.preconditioner =
2916 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2917
2918 this->BaseClass::initialize(A, parameters);
2919}
2920
2921//---------------------------------------------------------------------------
2922
2923template <typename MatrixType>
2924inline void
2926 const AdditionalData &parameters_in)
2927{
2928 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2929 Assert(
2930 parameters_in.relaxation != 0.0,
2931 ExcMessage(
2932 "Relaxation cannot automatically be determined by PreconditionSSOR."));
2933
2934 AdditionalData parameters;
2935 parameters.relaxation = 1.0;
2936 parameters.n_iterations = parameters_in.n_iterations;
2937 parameters.preconditioner =
2938 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2939
2940 this->BaseClass::initialize(A, parameters);
2941}
2942
2943
2944
2945//---------------------------------------------------------------------------
2946
2947template <typename MatrixType>
2948inline void
2950 const MatrixType &A,
2951 const std::vector<size_type> &p,
2952 const std::vector<size_type> &ip,
2953 const typename BaseClass::AdditionalData &parameters_in)
2954{
2955 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2956 Assert(
2957 parameters_in.relaxation != 0.0,
2958 ExcMessage(
2959 "Relaxation cannot automatically be determined by PreconditionPSOR."));
2960
2961 typename BaseClass::AdditionalData parameters;
2962 parameters.relaxation = 1.0;
2963 parameters.n_iterations = parameters_in.n_iterations;
2964 parameters.preconditioner =
2965 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2966
2967 this->BaseClass::initialize(A, parameters);
2968}
2969
2970
2971template <typename MatrixType>
2972inline void
2974 const AdditionalData &additional_data)
2975{
2976 initialize(A,
2977 additional_data.permutation,
2978 additional_data.inverse_permutation,
2979 additional_data.parameters);
2980}
2981
2982template <typename MatrixType>
2984 const std::vector<size_type> &permutation,
2985 const std::vector<size_type> &inverse_permutation,
2987 &parameters)
2988 : permutation(permutation)
2989 , inverse_permutation(inverse_permutation)
2990 , parameters(parameters)
2991{}
2992
2993
2994//---------------------------------------------------------------------------
2995
2996
2997template <typename MatrixType, typename VectorType>
2999 const MatrixType &M,
3000 const function_ptr method)
3001 : matrix(M)
3002 , precondition(method)
3003{}
3004
3005
3006
3007template <typename MatrixType, typename VectorType>
3008void
3010 VectorType &dst,
3011 const VectorType &src) const
3012{
3013 (matrix.*precondition)(dst, src);
3014}
3015
3016//---------------------------------------------------------------------------
3017
3018namespace internal
3019{
3020
3021 template <typename PreconditionerType>
3024 const double smoothing_range,
3025 const unsigned int eig_cg_n_iterations,
3026 const double eig_cg_residual,
3027 const double max_eigenvalue,
3028 const EigenvalueAlgorithm eigenvalue_algorithm)
3029 : smoothing_range(smoothing_range)
3030 , eig_cg_n_iterations(eig_cg_n_iterations)
3031 , eig_cg_residual(eig_cg_residual)
3032 , max_eigenvalue(max_eigenvalue)
3033 , eigenvalue_algorithm(eigenvalue_algorithm)
3034 {}
3035
3036
3037
3038 template <typename PreconditionerType>
3039 inline EigenvalueAlgorithmAdditionalData<PreconditionerType> &
3040 EigenvalueAlgorithmAdditionalData<PreconditionerType>::operator=(
3041 const EigenvalueAlgorithmAdditionalData &other_data)
3042 {
3043 smoothing_range = other_data.smoothing_range;
3044 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3045 eig_cg_residual = other_data.eig_cg_residual;
3046 max_eigenvalue = other_data.max_eigenvalue;
3047 preconditioner = other_data.preconditioner;
3048 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3049 constraints.copy_from(other_data.constraints);
3050
3051 return *this;
3052 }
3053} // namespace internal
3054
3055template <typename MatrixType, typename PreconditionerType>
3057 AdditionalData(const double relaxation,
3058 const unsigned int n_iterations,
3059 const double smoothing_range,
3060 const unsigned int eig_cg_n_iterations,
3061 const double eig_cg_residual,
3062 const double max_eigenvalue,
3063 const EigenvalueAlgorithm eigenvalue_algorithm)
3064 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3065 smoothing_range,
3066 eig_cg_n_iterations,
3067 eig_cg_residual,
3068 max_eigenvalue,
3069 eigenvalue_algorithm)
3070 , relaxation(relaxation)
3071 , n_iterations(n_iterations)
3072{}
3073
3074
3075
3076//---------------------------------------------------------------------------
3077
3078namespace internal
3079{
3080 namespace PreconditionChebyshevImplementation
3081 {
3082 // for deal.II vectors, perform updates for Chebyshev preconditioner all
3083 // at once to reduce memory transfer. Here, we select between general
3084 // vectors and deal.II vectors where we expand the loop over the (local)
3085 // size of the vector
3086
3087 // generic part for non-deal.II vectors
3088 template <typename VectorType, typename PreconditionerType>
3089 inline void
3090 vector_updates(const VectorType &rhs,
3091 const PreconditionerType &preconditioner,
3092 const unsigned int iteration_index,
3093 const double factor1,
3094 const double factor2,
3095 VectorType &solution_old,
3096 VectorType &temp_vector1,
3097 VectorType &temp_vector2,
3098 VectorType &solution)
3099 {
3100 if (iteration_index == 0)
3101 {
3102 solution.equ(factor2, rhs);
3103 preconditioner.vmult(solution_old, solution);
3104 }
3105 else if (iteration_index == 1)
3106 {
3107 // compute t = P^{-1} * (b-A*x^{n})
3108 temp_vector1.sadd(-1.0, 1.0, rhs);
3109 preconditioner.vmult(solution_old, temp_vector1);
3110
3111 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
3112 solution_old.sadd(factor2, 1 + factor1, solution);
3113 }
3114 else
3115 {
3116 // compute t = P^{-1} * (b-A*x^{n})
3117 temp_vector1.sadd(-1.0, 1.0, rhs);
3118 preconditioner.vmult(temp_vector2, temp_vector1);
3119
3120 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
3121 solution_old.sadd(-factor1, factor2, temp_vector2);
3122 solution_old.add(1 + factor1, solution);
3123 }
3124
3125 solution.swap(solution_old);
3126 }
3127
3128 // generic part for deal.II vectors
3129 template <
3130 typename Number,
3131 typename PreconditionerType,
3132 std::enable_if_t<
3133 !has_vmult_with_std_functions_for_precondition<
3134 PreconditionerType,
3136 int> * = nullptr>
3137 inline void
3138 vector_updates(
3140 const PreconditionerType &preconditioner,
3141 const unsigned int iteration_index,
3142 const double factor1_,
3143 const double factor2_,
3145 &solution_old,
3147 &temp_vector1,
3149 &temp_vector2,
3151 {
3152 const Number factor1 = factor1_;
3153 const Number factor1_plus_1 = 1. + factor1_;
3154 const Number factor2 = factor2_;
3155
3156 if (iteration_index == 0)
3157 {
3158 // compute t = P^{-1} * (b)
3159 preconditioner.vmult(solution_old, rhs);
3160
3161 // compute x^{n+1} = f_2 * t
3162 const auto solution_old_ptr = solution_old.begin();
3164 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3165 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
3166 }
3167 else if (iteration_index == 1)
3168 {
3169 // compute t = P^{-1} * (b-A*x^{n})
3170 temp_vector1.sadd(-1.0, 1.0, rhs);
3171
3172 preconditioner.vmult(solution_old, temp_vector1);
3173
3174 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
3175 const auto solution_ptr = solution.begin();
3176 const auto solution_old_ptr = solution_old.begin();
3178 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3179 solution_old_ptr[i] =
3180 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
3181 }
3182 else
3183 {
3184 // compute t = P^{-1} * (b-A*x^{n})
3185 temp_vector1.sadd(-1.0, 1.0, rhs);
3186
3187 preconditioner.vmult(temp_vector2, temp_vector1);
3188
3189 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
3190 const auto solution_ptr = solution.begin();
3191 const auto solution_old_ptr = solution_old.begin();
3192 const auto temp_vector2_ptr = temp_vector2.begin();
3194 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3195 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3196 factor1 * solution_old_ptr[i] +
3197 temp_vector2_ptr[i] * factor2;
3198 }
3199
3200 solution.swap(solution_old);
3201 }
3202
3203 template <
3204 typename Number,
3205 typename PreconditionerType,
3206 std::enable_if_t<
3207 has_vmult_with_std_functions_for_precondition<
3208 PreconditionerType,
3210 int> * = nullptr>
3211 inline void
3212 vector_updates(
3214 const PreconditionerType &preconditioner,
3215 const unsigned int iteration_index,
3216 const double factor1_,
3217 const double factor2_,
3219 &solution_old,
3221 &temp_vector1,
3223 &temp_vector2,
3225 {
3226 const Number factor1 = factor1_;
3227 const Number factor1_plus_1 = 1. + factor1_;
3228 const Number factor2 = factor2_;
3229
3230 const auto rhs_ptr = rhs.begin();
3231 const auto temp_vector1_ptr = temp_vector1.begin();
3232 const auto temp_vector2_ptr = temp_vector2.begin();
3233 const auto solution_ptr = solution.begin();
3234 const auto solution_old_ptr = solution_old.begin();
3235
3236 if (iteration_index == 0)
3237 {
3238 preconditioner.vmult(
3239 solution,
3240 rhs,
3241 [&](const auto start_range, const auto end_range) {
3242 if (end_range > start_range)
3243 std::memset(solution_ptr + start_range,
3244 0,
3245 sizeof(Number) * (end_range - start_range));
3246 },
3247 [&](const auto begin, const auto end) {
3249 for (std::size_t i = begin; i < end; ++i)
3250 solution_ptr[i] *= factor2;
3251 });
3252 }
3253 else
3254 {
3255 preconditioner.vmult(
3256 temp_vector2,
3257 temp_vector1,
3258 [&](const auto begin, const auto end) {
3259 if (end > begin)
3260 std::memset(temp_vector2_ptr + begin,
3261 0,
3262 sizeof(Number) * (end - begin));
3263
3265 for (std::size_t i = begin; i < end; ++i)
3266 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3267 },
3268 [&](const auto begin, const auto end) {
3269 if (iteration_index == 1)
3270 {
3272 for (std::size_t i = begin; i < end; ++i)
3273 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
3274 factor2 * temp_vector2_ptr[i];
3275 }
3276 else
3277 {
3279 for (std::size_t i = begin; i < end; ++i)
3280 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3281 factor1 * solution_old_ptr[i] +
3282 factor2 * temp_vector2_ptr[i];
3283 }
3284 });
3285 }
3286
3287 if (iteration_index > 0)
3288 {
3289 solution_old.swap(temp_vector2);
3290 solution_old.swap(solution);
3291 }
3292 }
3293
3294 // worker routine for deal.II vectors. Because of vectorization, we need
3295 // to put the loop into an extra structure because the virtual function of
3296 // VectorUpdatesRange prevents the compiler from applying vectorization.
3297 template <typename Number>
3298 struct VectorUpdater
3299 {
3300 VectorUpdater(const Number *rhs,
3301 const Number *matrix_diagonal_inverse,
3302 const unsigned int iteration_index,
3303 const Number factor1,
3304 const Number factor2,
3305 Number *solution_old,
3306 Number *tmp_vector,
3307 Number *solution)
3308 : rhs(rhs)
3309 , matrix_diagonal_inverse(matrix_diagonal_inverse)
3310 , iteration_index(iteration_index)
3311 , factor1(factor1)
3312 , factor2(factor2)
3313 , solution_old(solution_old)
3314 , tmp_vector(tmp_vector)
3315 , solution(solution)
3316 {}
3317
3318 void
3319 apply_to_subrange(const std::size_t begin, const std::size_t end) const
3320 {
3321 // To circumvent a bug in gcc
3322 // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
3323 // copies of the variables factor1 and factor2 and do not check based on
3324 // factor1.
3325 const Number factor1 = this->factor1;
3326 const Number factor1_plus_1 = 1. + this->factor1;
3327 const Number factor2 = this->factor2;
3328 if (iteration_index == 0)
3329 {
3331 for (std::size_t i = begin; i < end; ++i)
3332 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
3333 }
3334 else if (iteration_index == 1)
3335 {
3336 // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
3338 for (std::size_t i = begin; i < end; ++i)
3339 // for efficiency reason, write back to temp_vector that is
3340 // already read (avoid read-for-ownership)
3341 tmp_vector[i] =
3342 factor1_plus_1 * solution[i] +
3343 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3344 }
3345 else
3346 {
3347 // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
3348 // + f_2 * P^{-1} * (b-A*x^{n})
3350 for (std::size_t i = begin; i < end; ++i)
3351 // for efficiency reason, write back to temp_vector, which is
3352 // already modified during vmult (in best case, the modified
3353 // values are not written back to main memory yet so that
3354 // we do not have to pay additional costs for writing and
3355 // read-for-ownership)
3356 tmp_vector[i] =
3357 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
3358 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3359 }
3360 }
3361
3362 const Number *rhs;
3363 const Number *matrix_diagonal_inverse;
3364 const unsigned int iteration_index;
3365 const Number factor1;
3366 const Number factor2;
3367 mutable Number *solution_old;
3368 mutable Number *tmp_vector;
3369 mutable Number *solution;
3370 };
3371
3372 template <typename Number>
3373 struct VectorUpdatesRange : public ::parallel::ParallelForInteger
3374 {
3375 VectorUpdatesRange(const VectorUpdater<Number> &updater,
3376 const std::size_t size)
3377 : updater(updater)
3378 {
3380 VectorUpdatesRange::apply_to_subrange(0, size);
3381 else
3382 apply_parallel(
3383 0,
3384 size,
3386 }
3387
3388 ~VectorUpdatesRange() override = default;
3389
3390 virtual void
3391 apply_to_subrange(const std::size_t begin,
3392 const std::size_t end) const override
3393 {
3394 updater.apply_to_subrange(begin, end);
3395 }
3396
3397 const VectorUpdater<Number> &updater;
3398 };
3399
3400 // selection for diagonal matrix around deal.II vector
3401 template <typename Number>
3402 inline void
3403 vector_updates(
3404 const ::Vector<Number> &rhs,
3405 const ::DiagonalMatrix<::Vector<Number>> &jacobi,
3406 const unsigned int iteration_index,
3407 const double factor1,
3408 const double factor2,
3409 ::Vector<Number> &solution_old,
3410 ::Vector<Number> &temp_vector1,
3412 ::Vector<Number> &solution)
3413 {
3414 VectorUpdater<Number> upd(rhs.begin(),
3415 jacobi.get_vector().begin(),
3416 iteration_index,
3417 factor1,
3418 factor2,
3419 solution_old.begin(),
3420 temp_vector1.begin(),
3421 solution.begin());
3422 VectorUpdatesRange<Number>(upd, rhs.size());
3423
3424 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3425 if (iteration_index == 0)
3426 {
3427 // nothing to do here because we can immediately write into the
3428 // solution vector without remembering any of the other vectors
3429 }
3430 else
3431 {
3432 solution.swap(temp_vector1);
3433 solution_old.swap(temp_vector1);
3434 }
3435 }
3436
3437 // selection for diagonal matrix around parallel deal.II vector
3438 template <typename Number>
3439 inline void
3440 vector_updates(
3442 const ::DiagonalMatrix<
3444 const unsigned int iteration_index,
3445 const double factor1,
3446 const double factor2,
3448 &solution_old,
3450 &temp_vector1,
3453 {
3454 VectorUpdater<Number> upd(rhs.begin(),
3455 jacobi.get_vector().begin(),
3456 iteration_index,
3457 factor1,
3458 factor2,
3459 solution_old.begin(),
3460 temp_vector1.begin(),
3461 solution.begin());
3462 VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
3463
3464 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3465 if (iteration_index == 0)
3466 {
3467 // nothing to do here because we can immediately write into the
3468 // solution vector without remembering any of the other vectors
3469 }
3470 else
3471 {
3472 solution.swap(temp_vector1);
3473 solution_old.swap(temp_vector1);
3474 }
3475 }
3476
3477 // We need to have a separate declaration for static const members
3478
3479 // general case and the case that the preconditioner can work on
3480 // ranges (covered by vector_updates())
3481 template <
3482 typename MatrixType,
3483 typename VectorType,
3484 typename PreconditionerType,
3485 std::enable_if_t<
3486 !has_vmult_with_std_functions<MatrixType,
3487 VectorType,
3488 PreconditionerType> &&
3489 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3490 VectorType> &&
3491 has_vmult_with_std_functions_for_precondition<MatrixType,
3492 VectorType>),
3493 int> * = nullptr>
3494 inline void
3495 vmult_and_update(const MatrixType &matrix,
3496 const PreconditionerType &preconditioner,
3497 const VectorType &rhs,
3498 const unsigned int iteration_index,
3499 const double factor1,
3500 const double factor2,
3501 VectorType &solution,
3502 VectorType &solution_old,
3503 VectorType &temp_vector1,
3504 VectorType &temp_vector2)
3505 {
3506 if (iteration_index > 0)
3507 matrix.vmult(temp_vector1, solution);
3508 vector_updates(rhs,
3509 preconditioner,
3510 iteration_index,
3511 factor1,
3512 factor2,
3513 solution_old,
3514 temp_vector1,
3515 temp_vector2,
3516 solution);
3517 }
3518
3519 // case that both the operator and the preconditioner can work on
3520 // subranges
3521 template <
3522 typename MatrixType,
3523 typename VectorType,
3524 typename PreconditionerType,
3525 std::enable_if_t<
3526 !has_vmult_with_std_functions<MatrixType,
3527 VectorType,
3528 PreconditionerType> &&
3529 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3530 VectorType> &&
3531 has_vmult_with_std_functions_for_precondition<MatrixType,
3532 VectorType>),
3533 int> * = nullptr>
3534 inline void
3535 vmult_and_update(const MatrixType &matrix,
3536 const PreconditionerType &preconditioner,
3537 const VectorType &rhs,
3538 const unsigned int iteration_index,
3539 const double factor1_,
3540 const double factor2_,
3541 VectorType &solution,
3542 VectorType &solution_old,
3543 VectorType &temp_vector1,
3544 VectorType &temp_vector2)
3545 {
3546 using Number = typename VectorType::value_type;
3547
3548 const Number factor1 = factor1_;
3549 const Number factor1_plus_1 = 1. + factor1_;
3550 const Number factor2 = factor2_;
3551
3552 if (iteration_index == 0)
3553 {
3554 preconditioner.vmult(
3555 solution,
3556 rhs,
3557 [&](const unsigned int start_range, const unsigned int end_range) {
3558 // zero 'solution' before running the vmult operation
3559 if (end_range > start_range)
3560 std::memset(solution.begin() + start_range,
3561 0,
3562 sizeof(Number) * (end_range - start_range));
3563 },
3564 [&](const unsigned int start_range, const unsigned int end_range) {
3565 const auto solution_ptr = solution.begin();
3566
3568 for (std::size_t i = start_range; i < end_range; ++i)
3569 solution_ptr[i] *= factor2;
3570 });
3571 }
3572 else
3573 {
3574 temp_vector1.reinit(rhs, true);
3575 temp_vector2.reinit(rhs, true);
3576
3577 // 1) compute residual (including operator application)
3578 matrix.vmult(
3579 temp_vector1,
3580 solution,
3581 [&](const unsigned int start_range, const unsigned int end_range) {
3582 // zero 'temp_vector1' before running the vmult
3583 // operation
3584 if (end_range > start_range)
3585 std::memset(temp_vector1.begin() + start_range,
3586 0,
3587 sizeof(Number) * (end_range - start_range));
3588 },
3589 [&](const unsigned int start_range, const unsigned int end_range) {
3590 const auto rhs_ptr = rhs.begin();
3591 const auto tmp_ptr = temp_vector1.begin();
3592
3594 for (std::size_t i = start_range; i < end_range; ++i)
3595 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3596 });
3597
3598 // 2) perform vector updates (including preconditioner application)
3599 preconditioner.vmult(
3600 temp_vector2,
3601 temp_vector1,
3602 [&](const unsigned int start_range, const unsigned int end_range) {
3603 // zero 'temp_vector2' before running the vmult
3604 // operation
3605 if (end_range > start_range)
3606 std::memset(temp_vector2.begin() + start_range,
3607 0,
3608 sizeof(Number) * (end_range - start_range));
3609 },
3610 [&](const unsigned int start_range, const unsigned int end_range) {
3611 const auto solution_ptr = solution.begin();
3612 const auto solution_old_ptr = solution_old.begin();
3613 const auto tmp_ptr = temp_vector2.begin();
3614
3615 if (iteration_index == 1)
3616 {
3618 for (std::size_t i = start_range; i < end_range; ++i)
3619 tmp_ptr[i] =
3620 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3621 }
3622 else
3623 {
3625 for (std::size_t i = start_range; i < end_range; ++i)
3626 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3627 factor1 * solution_old_ptr[i] +
3628 factor2 * tmp_ptr[i];
3629 }
3630 });
3631
3632 solution.swap(temp_vector2);
3633 solution_old.swap(temp_vector2);
3634 }
3635 }
3636
3637 // case that the operator can work on subranges and the preconditioner
3638 // is a diagonal
3639 template <typename MatrixType,
3640 typename VectorType,
3641 typename PreconditionerType,
3642 std::enable_if_t<has_vmult_with_std_functions<MatrixType,
3643 VectorType,
3644 PreconditionerType>,
3645 int> * = nullptr>
3646 inline void
3647 vmult_and_update(const MatrixType &matrix,
3648 const PreconditionerType &preconditioner,
3649 const VectorType &rhs,
3650 const unsigned int iteration_index,
3651 const double factor1,
3652 const double factor2,
3653 VectorType &solution,
3654 VectorType &solution_old,
3655 VectorType &temp_vector1,
3656 VectorType &)
3657 {
3658 using Number = typename VectorType::value_type;
3659 VectorUpdater<Number> updater(rhs.begin(),
3660 preconditioner.get_vector().begin(),
3661 iteration_index,
3662 factor1,
3663 factor2,
3664 solution_old.begin(),
3665 temp_vector1.begin(),
3666 solution.begin());
3667 if (iteration_index > 0)
3668 matrix.vmult(
3669 temp_vector1,
3670 solution,
3671 [&](const unsigned int start_range, const unsigned int end_range) {
3672 // zero 'temp_vector1' before running the vmult
3673 // operation
3674 if (end_range > start_range)
3675 std::memset(temp_vector1.begin() + start_range,
3676 0,
3677 sizeof(Number) * (end_range - start_range));
3678 },
3679 [&](const unsigned int start_range, const unsigned int end_range) {
3680 if (end_range > start_range)
3681 updater.apply_to_subrange(start_range, end_range);
3682 });
3683 else
3684 updater.apply_to_subrange(0U, solution.locally_owned_size());
3685
3686 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3687 if (iteration_index == 0)
3688 {
3689 // nothing to do here because we can immediately write into the
3690 // solution vector without remembering any of the other vectors
3691 }
3692 else
3693 {
3694 solution.swap(temp_vector1);
3695 solution_old.swap(temp_vector1);
3696 }
3697 }
3698
3699 template <typename MatrixType, typename PreconditionerType>
3700 inline void
3701 initialize_preconditioner(
3702 const MatrixType &matrix,
3703 std::shared_ptr<PreconditionerType> &preconditioner)
3704 {
3705 (void)matrix;
3706 (void)preconditioner;
3707 AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
3708 }
3709
3710 template <typename MatrixType, typename VectorType>
3711 inline void
3712 initialize_preconditioner(
3713 const MatrixType &matrix,
3714 std::shared_ptr<::DiagonalMatrix<VectorType>> &preconditioner)
3715 {
3716 if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
3717 {
3718 if (preconditioner.get() == nullptr)
3719 preconditioner =
3720 std::make_shared<::DiagonalMatrix<VectorType>>();
3721
3722 Assert(
3723 preconditioner->m() == 0,
3724 ExcMessage(
3725 "Preconditioner appears to be initialized but not sized correctly"));
3726
3727 // This part only works in serial
3728 if (preconditioner->m() != matrix.m())
3729 {
3730 preconditioner->get_vector().reinit(matrix.m());
3731 for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
3732 preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
3733 }
3734 }
3735 }
3736 } // namespace PreconditionChebyshevImplementation
3737} // namespace internal
3738
3739
3740
3741template <typename MatrixType, typename VectorType, typename PreconditionerType>
3743 AdditionalData::AdditionalData(const unsigned int degree,
3744 const double smoothing_range,
3745 const unsigned int eig_cg_n_iterations,
3746 const double eig_cg_residual,
3747 const double max_eigenvalue,
3748 const EigenvalueAlgorithm eigenvalue_algorithm,
3749 const PolynomialType polynomial_type)
3750 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3751 smoothing_range,
3752 eig_cg_n_iterations,
3753 eig_cg_residual,
3754 max_eigenvalue,
3755 eigenvalue_algorithm)
3756 , degree(degree)
3757 , polynomial_type(polynomial_type)
3758{}
3759
3760
3761
3762template <typename MatrixType, typename VectorType, typename PreconditionerType>
3765 : theta(1.)
3766 , delta(1.)
3767 , eigenvalues_are_initialized(false)
3768{
3769 static_assert(
3770 std::is_same_v<size_type, typename VectorType::size_type>,
3771 "PreconditionChebyshev and VectorType must have the same size_type.");
3772}
3773
3774
3775
3776template <typename MatrixType, typename VectorType, typename PreconditionerType>
3777inline void
3779 const MatrixType &matrix,
3780 const AdditionalData &additional_data)
3781{
3782 matrix_ptr = &matrix;
3783 data = additional_data;
3784 Assert(data.degree > 0,
3785 ExcMessage("The degree of the Chebyshev method must be positive."));
3786 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3787 matrix, data.preconditioner);
3788 eigenvalues_are_initialized = false;
3789}
3790
3791
3792
3793template <typename MatrixType, typename VectorType, typename PreconditionerType>
3794inline void
3796{
3797 eigenvalues_are_initialized = false;
3798 theta = delta = 1.0;
3799 matrix_ptr = nullptr;
3800 {
3801 VectorType empty_vector;
3802 solution_old.reinit(empty_vector);
3803 temp_vector1.reinit(empty_vector);
3804 temp_vector2.reinit(empty_vector);
3805 }
3806 data.preconditioner.reset();
3807}
3808
3809
3810
3811template <typename MatrixType, typename VectorType, typename PreconditionerType>
3812inline typename internal::EigenvalueInformation
3814 estimate_eigenvalues(const VectorType &src) const
3815{
3816 Assert(eigenvalues_are_initialized == false, ExcInternalError());
3817
3818 solution_old.reinit(src);
3819 temp_vector1.reinit(src, true);
3820
3821 auto info = internal::estimate_eigenvalues<MatrixType>(
3822 data, matrix_ptr, solution_old, temp_vector1, data.degree);
3823
3824 const double alpha = (data.smoothing_range > 1. ?
3825 info.max_eigenvalue_estimate / data.smoothing_range :
3826 std::min(0.9 * info.max_eigenvalue_estimate,
3827 info.min_eigenvalue_estimate));
3828
3829 // in case the user set the degree to invalid unsigned int, we have to
3830 // determine the number of necessary iterations from the Chebyshev error
3831 // estimate, given the target tolerance specified by smoothing_range. This
3832 // estimate is based on the error formula given in section 5.1 of
3833 // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
3834 if (data.degree == numbers::invalid_unsigned_int)
3835 {
3836 const double actual_range = info.max_eigenvalue_estimate / alpha;
3837 const double sigma = (1. - std::sqrt(1. / actual_range)) /
3838 (1. + std::sqrt(1. / actual_range));
3839 const double eps = data.smoothing_range;
3840 const_cast<
3842 this)
3843 ->data.degree =
3844 1 + static_cast<unsigned int>(
3845 std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
3846 std::log(1. / sigma));
3847 }
3848
3849 info.degree = data.degree;
3850
3851 const_cast<
3853 ->delta =
3854 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3855 (info.max_eigenvalue_estimate) :
3856 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3857 const_cast<
3859 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3860
3861 // We do not need the second temporary vector in case we have a
3862 // DiagonalMatrix as preconditioner and use deal.II's own vectors
3863 using NumberType = typename VectorType::value_type;
3864 if (std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> ==
3865 false ||
3866 (std::is_same_v<VectorType, ::Vector<NumberType>> == false &&
3867 ((std::is_same_v<
3868 VectorType,
3870 false) ||
3871 (std::is_same_v<VectorType,
3873 Vector<NumberType, MemorySpace::Default>> == false))))
3874 temp_vector2.reinit(src, true);
3875 else
3876 {
3877 VectorType empty_vector;
3878 temp_vector2.reinit(empty_vector);
3879 }
3880
3881 const_cast<
3883 ->eigenvalues_are_initialized = true;
3884
3885 return info;
3886}
3887
3888
3889
3890template <typename MatrixType, typename VectorType, typename PreconditionerType>
3891inline void
3893 VectorType &solution,
3894 const VectorType &rhs) const
3895{
3896 std::lock_guard<std::mutex> lock(mutex);
3897 if (eigenvalues_are_initialized == false)
3898 estimate_eigenvalues(rhs);
3899
3900 internal::PreconditionChebyshevImplementation::vmult_and_update(
3901 *matrix_ptr,
3902 *data.preconditioner,
3903 rhs,
3904 0,
3905 0.,
3906 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3907 (4. / (3. * delta)) :
3908 (1. / theta),
3909 solution,
3910 solution_old,
3911 temp_vector1,
3912 temp_vector2);
3913
3914 // if delta is zero, we do not need to iterate because the updates will be
3915 // zero
3916 if (data.degree < 2 || std::abs(delta) < 1e-40)
3917 return;
3918
3919 double rhok = delta / theta, sigma = theta / delta;
3920 for (unsigned int k = 0; k < data.degree - 1; ++k)
3921 {
3922 double factor1 = 0.0;
3923 double factor2 = 0.0;
3924
3925 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3926 {
3927 factor1 = (2 * k + 1.) / (2 * k + 5.);
3928 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3929 }
3930 else
3931 {
3932 const double rhokp = 1. / (2. * sigma - rhok);
3933 factor1 = rhokp * rhok;
3934 factor2 = 2. * rhokp / delta;
3935 rhok = rhokp;
3936 }
3937
3938 internal::PreconditionChebyshevImplementation::vmult_and_update(
3939 *matrix_ptr,
3940 *data.preconditioner,
3941 rhs,
3942 k + 1,
3943 factor1,
3944 factor2,
3945 solution,
3946 solution_old,
3947 temp_vector1,
3948 temp_vector2);
3949 }
3950}
3951
3952
3953
3954template <typename MatrixType, typename VectorType, typename PreconditionerType>
3955inline void
3957 VectorType &solution,
3958 const VectorType &rhs) const
3959{
3960 std::lock_guard<std::mutex> lock(mutex);
3961 if (eigenvalues_are_initialized == false)
3962 estimate_eigenvalues(rhs);
3963
3964 internal::PreconditionChebyshevImplementation::vector_updates(
3965 rhs,
3966 *data.preconditioner,
3967 0,
3968 0.,
3969 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3970 (4. / (3. * delta)) :
3971 (1. / theta),
3972 solution_old,
3973 temp_vector1,
3974 temp_vector2,
3975 solution);
3976
3977 if (data.degree < 2 || std::abs(delta) < 1e-40)
3978 return;
3979
3980 double rhok = delta / theta, sigma = theta / delta;
3981 for (unsigned int k = 0; k < data.degree - 1; ++k)
3982 {
3983 double factor1 = 0.0;
3984 double factor2 = 0.0;
3985
3986 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3987 {
3988 factor1 = (2 * k + 1.) / (2 * k + 5.);
3989 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3990 }
3991 else
3992 {
3993 const double rhokp = 1. / (2. * sigma - rhok);
3994 factor1 = rhokp * rhok;
3995 factor2 = 2. * rhokp / delta;
3996 rhok = rhokp;
3997 }
3998
3999 matrix_ptr->Tvmult(temp_vector1, solution);
4000 internal::PreconditionChebyshevImplementation::vector_updates(
4001 rhs,
4002 *data.preconditioner,
4003 k + 1,
4004 factor1,
4005 factor2,
4006 solution_old,
4007 temp_vector1,
4008 temp_vector2,
4009 solution);
4010 }
4011}
4012
4013
4014
4015template <typename MatrixType, typename VectorType, typename PreconditionerType>
4016inline void
4018 VectorType &solution,
4019 const VectorType &rhs) const
4020{
4021 std::lock_guard<std::mutex> lock(mutex);
4022 if (eigenvalues_are_initialized == false)
4023 estimate_eigenvalues(rhs);
4024
4025 internal::PreconditionChebyshevImplementation::vmult_and_update(
4026 *matrix_ptr,
4027 *data.preconditioner,
4028 rhs,
4029 1,
4030 0.,
4031 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4032 (4. / (3. * delta)) :
4033 (1. / theta),
4034 solution,
4035 solution_old,
4036 temp_vector1,
4037 temp_vector2);
4038
4039 if (data.degree < 2 || std::abs(delta) < 1e-40)
4040 return;
4041
4042 double rhok = delta / theta, sigma = theta / delta;
4043 for (unsigned int k = 0; k < data.degree - 1; ++k)
4044 {
4045 double factor1 = 0.0;
4046 double factor2 = 0.0;
4047
4048 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4049 {
4050 factor1 = (2 * k + 1.) / (2 * k + 5.);
4051 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4052 }
4053 else
4054 {
4055 const double rhokp = 1. / (2. * sigma - rhok);
4056 factor1 = rhokp * rhok;
4057 factor2 = 2. * rhokp / delta;
4058 rhok = rhokp;
4059 }
4060
4061 internal::PreconditionChebyshevImplementation::vmult_and_update(
4062 *matrix_ptr,
4063 *data.preconditioner,
4064 rhs,
4065 k + 2,
4066 factor1,
4067 factor2,
4068 solution,
4069 solution_old,
4070 temp_vector1,
4071 temp_vector2);
4072 }
4073}
4074
4075
4076
4077template <typename MatrixType, typename VectorType, typename PreconditionerType>
4078inline void
4080 VectorType &solution,
4081 const VectorType &rhs) const
4082{
4083 std::lock_guard<std::mutex> lock(mutex);
4084 if (eigenvalues_are_initialized == false)
4085 estimate_eigenvalues(rhs);
4086
4087 matrix_ptr->Tvmult(temp_vector1, solution);
4088 internal::PreconditionChebyshevImplementation::vector_updates(
4089 rhs,
4090 *data.preconditioner,
4091 1,
4092 0.,
4093 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4094 (4. / (3. * delta)) :
4095 (1. / theta),
4096 solution_old,
4097 temp_vector1,
4098 temp_vector2,
4099 solution);
4100
4101 if (data.degree < 2 || std::abs(delta) < 1e-40)
4102 return;
4103
4104 double rhok = delta / theta, sigma = theta / delta;
4105 for (unsigned int k = 0; k < data.degree - 1; ++k)
4106 {
4107 double factor1 = 0.0;
4108 double factor2 = 0.0;
4109
4110 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4111 {
4112 factor1 = (2 * k + 1.) / (2 * k + 5.);
4113 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4114 }
4115 else
4116 {
4117 const double rhokp = 1. / (2. * sigma - rhok);
4118 factor1 = rhokp * rhok;
4119 factor2 = 2. * rhokp / delta;
4120 rhok = rhokp;
4121 }
4122
4123 matrix_ptr->Tvmult(temp_vector1, solution);
4124 internal::PreconditionChebyshevImplementation::vector_updates(
4125 rhs,
4126 *data.preconditioner,
4127 k + 2,
4128 factor1,
4129 factor2,
4130 solution_old,
4131 temp_vector1,
4132 temp_vector2,
4133 solution);
4134 }
4135}
4136
4137
4138
4139template <typename MatrixType, typename VectorType, typename PreconditionerType>
4140inline typename PreconditionChebyshev<MatrixType,
4141 VectorType,
4142 PreconditionerType>::size_type
4144{
4145 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4146 return matrix_ptr->m();
4147}
4148
4149
4150
4151template <typename MatrixType, typename VectorType, typename PreconditionerType>
4152inline typename PreconditionChebyshev<MatrixType,
4153 VectorType,
4154 PreconditionerType>::size_type
4156{
4157 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4158 return matrix_ptr->n();
4159}
4160
4161#endif // DOXYGEN
4162
4164
4165#endif
unsigned int n_blocks() const
BlockType & block(const unsigned int i)
bool is_empty() const
Definition index_set.h:1937
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1994
void sadd(const Number s, const Number a, const Vector< Number, MemorySpace > &V)
size_type locally_owned_size() const
void swap(Vector< Number, MemorySpace > &v) noexcept
::IndexSet locally_owned_elements() const
void Tvmult(VectorType &dst, const VectorType &src) const
size_type m() const
size_type n() const
void step(VectorType &dst, const VectorType &src) const
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
void Tstep(VectorType &dst, const VectorType &src) const
void vmult(VectorType &dst, const VectorType &src) const
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
ObserverPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
void vmult_add(VectorType &, const VectorType &) const
void vmult(VectorType &, const VectorType &) const
size_type m() const
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
size_type n() const
void Tvmult(VectorType &, const VectorType &) const
void Tvmult_add(VectorType &, const VectorType &) const
typename BaseClass::AdditionalData AdditionalData
internal::PreconditionRelaxation::PreconditionJacobiImpl< MatrixType > PreconditionerType
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
AdditionalData(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
BaseClass::AdditionalData parameters
const std::vector< size_type > & inverse_permutation
const std::vector< size_type > & permutation
typename BaseClass::size_type size_type
void initialize(const MatrixType &A, const AdditionalData &additional_data)
internal::PreconditionRelaxation::PreconditionPSORImpl< MatrixType > PreconditionerType
void initialize(const MatrixType &A, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
AdditionalData(const double relaxation=1., const unsigned int n_iterations=1, const double smoothing_range=0., const unsigned int eig_cg_n_iterations=8, const double eig_cg_residual=1e-2, const double max_eigenvalue=1, const EigenvalueAlgorithm eigenvalue_algorithm=EigenvalueAlgorithm::lanczos)
double get_relaxation() const
void Tvmult(VectorType &, const VectorType &) const
std::shared_ptr< PreconditionerType > preconditioner
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
size_type n() const
size_type m() const
void step(VectorType &x, const VectorType &rhs) const
ObserverPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
void Tstep(VectorType &x, const VectorType &rhs) const
void vmult(VectorType &, const VectorType &) const
types::global_dof_index size_type
AdditionalData(const double relaxation=1.)
types::global_dof_index size_type
void vmult_add(VectorType &, const VectorType &) const
size_type n() const
void initialize(const AdditionalData &parameters)
void initialize(const MatrixType &matrix, const AdditionalData &parameters)
size_type m() const
void vmult(VectorType &, const VectorType &) const
void Tvmult(VectorType &, const VectorType &) const
void Tvmult_add(VectorType &, const VectorType &) const
internal::PreconditionRelaxation::PreconditionSORImpl< MatrixType > PreconditionerType
typename BaseClass::AdditionalData AdditionalData
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
typename BaseClass::AdditionalData AdditionalData
internal::PreconditionRelaxation::PreconditionSSORImpl< MatrixType > PreconditionerType
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
const function_ptr precondition
const MatrixType & matrix
void vmult(VectorType &dst, const VectorType &src) const
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
const_iterator end() const
const_iterator begin() const
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
Number mean_value() const
virtual size_type size() const override
virtual void swap(Vector< Number > &v) noexcept
iterator begin()
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:208
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
Definition mpi.cc:746
std::size_t size
Definition mpi.cc:745
@ matrix
Contents is actually a matrix.
constexpr char A
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
unsigned int minimum_parallel_grain_size
Definition parallel.cc:50
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
STL namespace.
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:94
AdditionalData(const unsigned int degree=1, const double smoothing_range=0., const unsigned int eig_cg_n_iterations=8, const double eig_cg_residual=1e-2, const double max_eigenvalue=1, const EigenvalueAlgorithm eigenvalue_algorithm=EigenvalueAlgorithm::lanczos, const PolynomialType polynomial_type=PolynomialType::first_kind)
std::shared_ptr< PreconditionerType > preconditioner
EigenvalueAlgorithmAdditionalData< PreconditionerType > & operator=(const EigenvalueAlgorithmAdditionalData< PreconditionerType > &other_data)
::AffineConstraints< double > constraints
EigenvalueAlgorithmAdditionalData(const double smoothing_range, const unsigned int eig_cg_n_iterations, const double eig_cg_residual, const double max_eigenvalue, const EigenvalueAlgorithm eigenvalue_algorithm)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)