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