deal.II version GIT relicensing-2013-g7f3fb24d6c 2024-10-21 19:30: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
1296 if (i == 0)
1297 {
1298 Number *dst_ptr = dst.begin();
1299 const Number *src_ptr = src.begin();
1300
1301 preconditioner.vmult(
1302 dst,
1303 src,
1304 [&](const unsigned int start_range, const unsigned int end_range) {
1305 // zero 'dst' before running the vmult operation
1306 if (end_range > start_range)
1307 std::memset(dst.begin() + start_range,
1308 0,
1309 sizeof(Number) * (end_range - start_range));
1310 },
1311 [&](const unsigned int start_range, const unsigned int end_range) {
1312 if (relaxation == 1.0)
1313 return; // nothing to do
1314
1315 const auto src_ptr = src.begin();
1316 const auto dst_ptr = dst.begin();
1317
1319 for (std::size_t i = start_range; i < end_range; ++i)
1320 dst_ptr[i] *= relaxation;
1321 });
1322 }
1323 else
1324 {
1325 tmp.reinit(src, true);
1326
1327 Assert(transposed == false, ExcNotImplemented());
1328
1329 A.vmult(tmp, dst);
1330
1331 preconditioner.vmult(
1332 dst,
1333 tmp,
1334 [&](const unsigned int start_range, const unsigned int end_range) {
1335 const auto src_ptr = src.begin();
1336 const auto tmp_ptr = tmp.begin();
1337
1338 if (relaxation == 1.0)
1339 {
1341 for (std::size_t i = start_range; i < end_range; ++i)
1342 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1343 }
1344 else
1345 {
1346 // note: we scale the residual here to be able to add into
1347 // the dst vector, which contains the solution from the last
1348 // iteration
1350 for (std::size_t i = start_range; i < end_range; ++i)
1351 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1352 }
1353 },
1354 [&](const unsigned int, const unsigned int) {
1355 // nothing to do, since scaling by the relaxation factor
1356 // has been done in the pre operation
1357 });
1358 }
1359 }
1360
1361 // 2) specialized implementation with a preconditioner and a matrix that
1362 // accepts ranges
1363 template <
1364 typename MatrixType,
1365 typename PreconditionerType,
1366 typename VectorType,
1367 std::enable_if_t<
1368 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1369 VectorType> &&
1370 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1371 int> * = nullptr>
1372 void
1373 step_operations(const MatrixType &A,
1374 const PreconditionerType &preconditioner,
1375 VectorType &dst,
1376 const VectorType &src,
1377 const double relaxation,
1378 VectorType &tmp,
1379 VectorType &,
1380 const unsigned int i,
1381 const bool transposed)
1382 {
1383 (void)transposed;
1384 using Number = typename VectorType::value_type;
1385
1386 if (i == 0)
1387 {
1388 Number *dst_ptr = dst.begin();
1389 const Number *src_ptr = src.begin();
1390
1391 preconditioner.vmult(
1392 dst,
1393 src,
1394 [&](const unsigned int start_range, const unsigned int end_range) {
1395 // zero 'dst' before running the vmult operation
1396 if (end_range > start_range)
1397 std::memset(dst.begin() + start_range,
1398 0,
1399 sizeof(Number) * (end_range - start_range));
1400 },
1401 [&](const unsigned int start_range, const unsigned int end_range) {
1402 if (relaxation == 1.0)
1403 return; // nothing to do
1404
1405 const auto src_ptr = src.begin();
1406 const auto dst_ptr = dst.begin();
1407
1409 for (std::size_t i = start_range; i < end_range; ++i)
1410 dst_ptr[i] *= relaxation;
1411 });
1412 }
1413 else
1414 {
1415 tmp.reinit(src, true);
1416
1417 Assert(transposed == false, ExcNotImplemented());
1418
1419 A.vmult(
1420 tmp,
1421 dst,
1422 [&](const unsigned int start_range, const unsigned int end_range) {
1423 // zero 'tmp' before running the vmult
1424 // operation
1425 if (end_range > start_range)
1426 std::memset(tmp.begin() + start_range,
1427 0,
1428 sizeof(Number) * (end_range - start_range));
1429 },
1430 [&](const unsigned int start_range, const unsigned int end_range) {
1431 const auto src_ptr = src.begin();
1432 const auto tmp_ptr = tmp.begin();
1433
1434 if (relaxation == 1.0)
1435 {
1437 for (std::size_t i = start_range; i < end_range; ++i)
1438 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1439 }
1440 else
1441 {
1442 // note: we scale the residual here to be able to add into
1443 // the dst vector, which contains the solution from the last
1444 // iteration
1446 for (std::size_t i = start_range; i < end_range; ++i)
1447 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1448 }
1449 });
1450
1451 preconditioner.vmult(dst, tmp, [](const auto, const auto) {
1452 // note: `dst` vector does not have to be zeroed out
1453 // since we add the result into it
1454 });
1455 }
1456 }
1457
1458 // 3) specialized implementation for inverse-diagonal preconditioner
1459 template <typename MatrixType,
1460 typename VectorType,
1461 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1462 !has_vmult_with_std_functions<
1463 MatrixType,
1464 VectorType,
1466 VectorType> * = nullptr>
1467 void
1468 step_operations(const MatrixType &A,
1469 const ::DiagonalMatrix<VectorType> &preconditioner,
1470 VectorType &dst,
1471 const VectorType &src,
1472 const double relaxation,
1473 VectorType &tmp,
1474 VectorType &,
1475 const unsigned int i,
1476 const bool transposed)
1477 {
1478 using Number = typename VectorType::value_type;
1479
1480 if (i == 0)
1481 {
1482 Number *dst_ptr = dst.begin();
1483 const Number *src_ptr = src.begin();
1484 const Number *diag_ptr = preconditioner.get_vector().begin();
1485
1486 if (relaxation == 1.0)
1487 {
1489 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1490 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1491 }
1492 else
1493 {
1495 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1496 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1497 }
1498 }
1499 else
1500 {
1501 tmp.reinit(src, true);
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 (transposed)
1509 Tvmult(A, tmp, dst);
1510 else
1511 A.vmult(tmp, dst);
1512
1513 if (relaxation == 1.0)
1514 {
1516 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1517 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1518 }
1519 else
1520 {
1522 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1523 dst_ptr[i] +=
1524 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1525 }
1526 }
1527 }
1528
1529 // 4) specialized implementation for inverse-diagonal preconditioner and
1530 // matrix that accepts ranges
1531 template <typename MatrixType,
1532 typename VectorType,
1533 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1534 has_vmult_with_std_functions<
1535 MatrixType,
1536 VectorType,
1538 VectorType> * = nullptr>
1539 void
1540 step_operations(const MatrixType &A,
1541 const ::DiagonalMatrix<VectorType> &preconditioner,
1542 VectorType &dst,
1543 const VectorType &src,
1544 const double relaxation,
1545 VectorType &tmp,
1546 VectorType &,
1547 const unsigned int i,
1548 const bool transposed)
1549 {
1550 (void)transposed;
1551 using Number = typename VectorType::value_type;
1552
1553 if (i == 0)
1554 {
1555 Number *dst_ptr = dst.begin();
1556 const Number *src_ptr = src.begin();
1557 const Number *diag_ptr = preconditioner.get_vector().begin();
1558
1559 if (relaxation == 1.0)
1560 {
1562 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1563 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1564 }
1565 else
1566 {
1568 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1569 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1570 }
1571 }
1572 else
1573 {
1574 tmp.reinit(src, true);
1575
1576 Assert(transposed == false, ExcNotImplemented());
1577
1578 A.vmult(
1579 tmp,
1580 dst,
1581 [&](const unsigned int start_range, const unsigned int end_range) {
1582 // zero 'tmp' before running the vmult operation
1583 if (end_range > start_range)
1584 std::memset(tmp.begin() + start_range,
1585 0,
1586 sizeof(Number) * (end_range - start_range));
1587 },
1588 [&](const unsigned int begin, const unsigned int end) {
1589 const Number *dst_ptr = dst.begin();
1590 const Number *src_ptr = src.begin();
1591 Number *tmp_ptr = tmp.begin();
1592 const Number *diag_ptr = preconditioner.get_vector().begin();
1593
1594 // for efficiency reason, write back to temp_vector that is
1595 // already read (avoid read-for-ownership)
1596 if (relaxation == 1.0)
1597 {
1599 for (std::size_t i = begin; i < end; ++i)
1600 tmp_ptr[i] =
1601 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1602 }
1603 else
1604 {
1606 for (std::size_t i = begin; i < end; ++i)
1607 tmp_ptr[i] = dst_ptr[i] + relaxation *
1608 (src_ptr[i] - tmp_ptr[i]) *
1609 diag_ptr[i];
1610 }
1611 });
1612
1613 tmp.swap(dst);
1614 }
1615 }
1616
1617 } // namespace PreconditionRelaxation
1618} // namespace internal
1619
1620#endif
1621
1622
1623
1650template <typename MatrixType = SparseMatrix<double>>
1652 : public PreconditionRelaxation<
1653 MatrixType,
1654 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1655{
1657 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1659
1660public:
1665
1669 void
1670 initialize(const MatrixType &A,
1671 const AdditionalData &parameters = AdditionalData());
1672};
1673
1674
1720template <typename MatrixType = SparseMatrix<double>>
1722 : public PreconditionRelaxation<
1723 MatrixType,
1724 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1725{
1727 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1729
1730public:
1735
1739 void
1740 initialize(const MatrixType &A,
1741 const AdditionalData &parameters = AdditionalData());
1742};
1743
1744
1745
1772template <typename MatrixType = SparseMatrix<double>>
1774 : public PreconditionRelaxation<
1775 MatrixType,
1776 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1777{
1779 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1781
1782public:
1787
1793 void
1794 initialize(const MatrixType &A,
1795 const AdditionalData &parameters = AdditionalData());
1796};
1797
1798
1828template <typename MatrixType = SparseMatrix<double>>
1830 : public PreconditionRelaxation<
1831 MatrixType,
1832 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1833{
1835 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1837
1838public:
1843
1848 {
1849 public:
1860 AdditionalData(const std::vector<size_type> &permutation,
1861 const std::vector<size_type> &inverse_permutation,
1862 const typename BaseClass::AdditionalData &parameters =
1863 typename BaseClass::AdditionalData());
1864
1868 const std::vector<size_type> &permutation;
1872 const std::vector<size_type> &inverse_permutation;
1877 };
1878
1890 void
1891 initialize(const MatrixType &A,
1892 const std::vector<size_type> &permutation,
1893 const std::vector<size_type> &inverse_permutation,
1894 const typename BaseClass::AdditionalData &parameters =
1895 typename BaseClass::AdditionalData());
1896
1907 void
1908 initialize(const MatrixType &A, const AdditionalData &additional_data);
1909};
1910
1911
1912
2101template <typename MatrixType = SparseMatrix<double>,
2102 typename VectorType = Vector<double>,
2103 typename PreconditionerType = DiagonalMatrix<VectorType>>
2105{
2106public:
2111
2117 : public internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>
2118 {
2120
2125 {
2129 first_kind,
2135 };
2136
2141 const unsigned int degree = 1,
2142 const double smoothing_range = 0.,
2143 const unsigned int eig_cg_n_iterations = 8,
2144 const double eig_cg_residual = 1e-2,
2145 const double max_eigenvalue = 1,
2147 EigenvalueAlgorithm::lanczos,
2149
2162 unsigned int degree;
2163
2168 };
2169
2170
2175
2187 void
2188 initialize(const MatrixType &matrix,
2189 const AdditionalData &additional_data = AdditionalData());
2190
2195 void
2196 vmult(VectorType &dst, const VectorType &src) const;
2197
2202 void
2203 Tvmult(VectorType &dst, const VectorType &src) const;
2204
2208 void
2209 step(VectorType &dst, const VectorType &src) const;
2210
2214 void
2215 Tstep(VectorType &dst, const VectorType &src) const;
2216
2220 void
2222
2227 size_type
2228 m() const;
2229
2234 size_type
2235 n() const;
2236
2238
2252 estimate_eigenvalues(const VectorType &src) const;
2253
2254private:
2259 const MatrixType,
2262
2266 mutable VectorType solution_old;
2267
2271 mutable VectorType temp_vector1;
2272
2276 mutable VectorType temp_vector2;
2277
2283
2287 double theta;
2288
2293 double delta;
2294
2300
2306};
2307
2308
2309
2311/* ---------------------------------- Inline functions ------------------- */
2312
2313#ifndef DOXYGEN
2314
2315
2316namespace internal
2317{
2318 template <typename VectorType>
2319 void
2320 set_initial_guess(VectorType &vector)
2321 {
2322 vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2323 if (vector.locally_owned_elements().is_element(0))
2324 vector(0) = 0.;
2325 }
2326
2327 template <typename Number>
2328 void
2329 set_initial_guess(::Vector<Number> &vector)
2330 {
2331 // Choose a high-frequency mode consisting of numbers between 0 and 1
2332 // that is cheap to compute (cheaper than random numbers) but avoids
2333 // obviously re-occurring numbers in multi-component systems by choosing
2334 // a period of 11
2335 for (unsigned int i = 0; i < vector.size(); ++i)
2336 vector(i) = i % 11;
2337
2338 const Number mean_value = vector.mean_value();
2339 vector.add(-mean_value);
2340 }
2341
2342 template <typename Number>
2343 void
2344 set_initial_guess(
2346 {
2347 for (unsigned int block = 0; block < vector.n_blocks(); ++block)
2348 set_initial_guess(vector.block(block));
2349 }
2350
2351 template <typename Number, typename MemorySpace>
2352 void
2353 set_initial_guess(
2355 {
2356 // Choose a high-frequency mode consisting of numbers between 0 and 1
2357 // that is cheap to compute (cheaper than random numbers) but avoids
2358 // obviously re-occurring numbers in multi-component systems by choosing
2359 // a period of 11.
2360 // Make initial guess robust with respect to number of processors
2361 // by operating on the global index.
2362 types::global_dof_index first_local_range = 0;
2363 if (!vector.locally_owned_elements().is_empty())
2364 first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2365
2366 const auto n_local_elements = vector.locally_owned_size();
2367 Number *values_ptr = vector.get_values();
2368 Kokkos::RangePolicy<typename MemorySpace::kokkos_space::execution_space,
2369 Kokkos::IndexType<types::global_dof_index>>
2370 policy(0, n_local_elements);
2371 Kokkos::parallel_for(
2372 "::PreconditionChebyshev::set_initial_guess",
2373 policy,
2374 KOKKOS_LAMBDA(types::global_dof_index i) {
2375 values_ptr[i] = (i + first_local_range) % 11;
2376 });
2377 const Number mean_value = vector.mean_value();
2378 vector.add(-mean_value);
2379 }
2380
2381 struct EigenvalueTracker
2382 {
2383 public:
2384 void
2385 slot(const std::vector<double> &eigenvalues)
2386 {
2388 }
2389
2390 std::vector<double> values;
2391 };
2392
2393
2394
2395 template <typename MatrixType,
2396 typename VectorType,
2397 typename PreconditionerType>
2398 double
2399 power_iteration(const MatrixType &matrix,
2400 VectorType &eigenvector,
2401 const PreconditionerType &preconditioner,
2402 const unsigned int n_iterations)
2403 {
2404 typename VectorType::value_type eigenvalue_estimate = 0.;
2405 eigenvector /= eigenvector.l2_norm();
2406 VectorType vector1, vector2;
2407 vector1.reinit(eigenvector, true);
2408 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2409 vector2.reinit(eigenvector, true);
2410 for (unsigned int i = 0; i < n_iterations; ++i)
2411 {
2412 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2413 {
2414 matrix.vmult(vector2, eigenvector);
2415 preconditioner.vmult(vector1, vector2);
2416 }
2417 else
2418 matrix.vmult(vector1, eigenvector);
2419
2420 eigenvalue_estimate = eigenvector * vector1;
2421
2422 vector1 /= vector1.l2_norm();
2423 eigenvector.swap(vector1);
2424 }
2425 return std::abs(eigenvalue_estimate);
2426 }
2427
2428
2429
2430 template <typename MatrixType,
2431 typename VectorType,
2432 typename PreconditionerType>
2433 EigenvalueInformation
2434 estimate_eigenvalues(
2435 const EigenvalueAlgorithmAdditionalData<PreconditionerType> &data,
2436 const MatrixType *matrix_ptr,
2437 VectorType &solution_old,
2438 VectorType &temp_vector1,
2439 const unsigned int degree)
2440 {
2441 Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2442
2443 EigenvalueInformation info{};
2444
2445 if (data.eig_cg_n_iterations > 0)
2446 {
2447 Assert(data.eig_cg_n_iterations > 2,
2448 ExcMessage(
2449 "Need to set at least two iterations to find eigenvalues."));
2450
2451 internal::EigenvalueTracker eigenvalue_tracker;
2452
2453 // set an initial guess that contains some high-frequency parts (to the
2454 // extent possible without knowing the discretization and the numbering)
2455 // to trigger high eigenvalues according to the external function
2456 internal::set_initial_guess(temp_vector1);
2457 data.constraints.set_zero(temp_vector1);
2458
2459 if (data.eigenvalue_algorithm == internal::EigenvalueAlgorithm::lanczos)
2460 {
2461 // set a very strict tolerance to force at least two iterations
2462 IterationNumberControl control(data.eig_cg_n_iterations,
2463 1e-10,
2464 false,
2465 false);
2466
2467 ::SolverCG<VectorType> solver(control);
2468 solver.connect_eigenvalues_slot(
2469 [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2470 eigenvalue_tracker.slot(eigenvalues);
2471 });
2472
2473 solver.solve(*matrix_ptr,
2474 solution_old,
2475 temp_vector1,
2476 *data.preconditioner);
2477
2478 info.cg_iterations = control.last_step();
2479 }
2480 else if (data.eigenvalue_algorithm ==
2482 {
2483 (void)degree;
2484
2486 ExcMessage("Cannot estimate the minimal eigenvalue with the "
2487 "power iteration"));
2488
2489 eigenvalue_tracker.values.push_back(
2490 internal::power_iteration(*matrix_ptr,
2491 temp_vector1,
2492 *data.preconditioner,
2493 data.eig_cg_n_iterations));
2494 }
2495 else
2497
2498 // read the eigenvalues from the attached eigenvalue tracker
2499 if (eigenvalue_tracker.values.empty())
2500 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2501 else
2502 {
2503 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2504
2505 // include a safety factor since the CG method will in general not
2506 // be converged
2507 info.max_eigenvalue_estimate =
2508 1.2 * eigenvalue_tracker.values.back();
2509 }
2510 }
2511 else
2512 {
2513 info.max_eigenvalue_estimate = data.max_eigenvalue;
2514 info.min_eigenvalue_estimate =
2515 data.max_eigenvalue / data.smoothing_range;
2516 }
2517
2518 return info;
2519 }
2520} // namespace internal
2521
2522
2523
2525 : n_rows(0)
2526 , n_columns(0)
2527{}
2528
2529template <typename MatrixType>
2530inline void
2531PreconditionIdentity::initialize(const MatrixType &matrix,
2533{
2534 n_rows = matrix.m();
2535 n_columns = matrix.n();
2536}
2537
2538
2539template <typename VectorType>
2540inline void
2541PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
2542{
2543 dst = src;
2544}
2545
2546
2547
2548template <typename VectorType>
2549inline void
2550PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
2551{
2552 dst = src;
2553}
2554
2555template <typename VectorType>
2556inline void
2557PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
2558{
2559 dst += src;
2560}
2561
2562
2563
2564template <typename VectorType>
2565inline void
2566PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
2567{
2568 dst += src;
2569}
2570
2571
2572
2573inline void
2575{}
2576
2577
2578
2581{
2583 return n_rows;
2584}
2585
2588{
2590 return n_columns;
2591}
2592
2593//---------------------------------------------------------------------------
2594
2596 const double relaxation)
2597 : relaxation(relaxation)
2598{}
2599
2600
2602 : relaxation(0)
2603 , n_rows(0)
2604 , n_columns(0)
2605{
2606 AdditionalData add_data;
2607 relaxation = add_data.relaxation;
2608}
2609
2610
2611
2612inline void
2615{
2616 relaxation = parameters.relaxation;
2617}
2618
2619
2620
2621template <typename MatrixType>
2622inline void
2624 const MatrixType &matrix,
2626{
2627 relaxation = parameters.relaxation;
2628 n_rows = matrix.m();
2629 n_columns = matrix.n();
2630}
2631
2632
2633
2634template <typename VectorType>
2635inline void
2636PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2637{
2638 static_assert(
2639 std::is_same_v<size_type, typename VectorType::size_type>,
2640 "PreconditionRichardson and VectorType must have the same size_type.");
2641
2642 dst.equ(relaxation, src);
2643}
2644
2645
2646
2647template <typename VectorType>
2648inline void
2649PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2650{
2651 static_assert(
2652 std::is_same_v<size_type, typename VectorType::size_type>,
2653 "PreconditionRichardson and VectorType must have the same size_type.");
2654
2655 dst.equ(relaxation, src);
2656}
2657
2658template <typename VectorType>
2659inline void
2660PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2661{
2662 static_assert(
2663 std::is_same_v<size_type, typename VectorType::size_type>,
2664 "PreconditionRichardson and VectorType must have the same size_type.");
2665
2666 dst.add(relaxation, src);
2667}
2668
2669
2670
2671template <typename VectorType>
2672inline void
2673PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2674{
2675 static_assert(
2676 std::is_same_v<size_type, typename VectorType::size_type>,
2677 "PreconditionRichardson and VectorType must have the same size_type.");
2678
2679 dst.add(relaxation, src);
2680}
2681
2684{
2686 return n_rows;
2687}
2688
2691{
2693 return n_columns;
2694}
2695
2696//---------------------------------------------------------------------------
2697
2698template <typename MatrixType, typename PreconditionerType>
2699inline void
2701 const MatrixType &rA,
2702 const AdditionalData &parameters)
2703{
2704 A = &rA;
2705 eigenvalues_are_initialized = false;
2706
2707 Assert(parameters.preconditioner, ExcNotInitialized());
2708
2709 this->data = parameters;
2710}
2711
2712
2713template <typename MatrixType, typename PreconditionerType>
2714inline void
2716{
2717 eigenvalues_are_initialized = false;
2718 A = nullptr;
2719 data.relaxation = 1.0;
2720 data.preconditioner = nullptr;
2721}
2722
2723template <typename MatrixType, typename PreconditionerType>
2724inline
2727{
2728 Assert(A != nullptr, ExcNotInitialized());
2729 return A->m();
2730}
2731
2732template <typename MatrixType, typename PreconditionerType>
2733inline
2736{
2737 Assert(A != nullptr, ExcNotInitialized());
2738 return A->n();
2739}
2740
2741template <typename MatrixType, typename PreconditionerType>
2742template <typename VectorType>
2743inline void
2745 VectorType &dst,
2746 const VectorType &src) const
2747{
2748 Assert(this->A != nullptr, ExcNotInitialized());
2749 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2750
2751 if (eigenvalues_are_initialized == false)
2752 estimate_eigenvalues(src);
2753
2754 VectorType tmp1, tmp2;
2755
2756 for (unsigned int i = 0; i < data.n_iterations; ++i)
2757 internal::PreconditionRelaxation::step_operations(*A,
2758 *data.preconditioner,
2759 dst,
2760 src,
2761 data.relaxation,
2762 tmp1,
2763 tmp2,
2764 i,
2765 false);
2766}
2767
2768template <typename MatrixType, typename PreconditionerType>
2769template <typename VectorType>
2770inline void
2772 VectorType &dst,
2773 const VectorType &src) const
2774{
2775 Assert(this->A != nullptr, ExcNotInitialized());
2776 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2777
2778 if (eigenvalues_are_initialized == false)
2779 estimate_eigenvalues(src);
2780
2781 VectorType tmp1, tmp2;
2782
2783 for (unsigned int i = 0; i < data.n_iterations; ++i)
2784 internal::PreconditionRelaxation::step_operations(
2785 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2786}
2787
2788template <typename MatrixType, typename PreconditionerType>
2789template <typename VectorType>
2790inline void
2792 VectorType &dst,
2793 const VectorType &src) const
2794{
2795 Assert(this->A != nullptr, ExcNotInitialized());
2796 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2797
2798 if (eigenvalues_are_initialized == false)
2799 estimate_eigenvalues(src);
2800
2801 VectorType tmp1, tmp2;
2802
2803 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2804 internal::PreconditionRelaxation::step_operations(*A,
2805 *data.preconditioner,
2806 dst,
2807 src,
2808 data.relaxation,
2809 tmp1,
2810 tmp2,
2811 i,
2812 false);
2813}
2814
2815template <typename MatrixType, typename PreconditionerType>
2816template <typename VectorType>
2817inline void
2819 VectorType &dst,
2820 const VectorType &src) const
2821{
2822 Assert(this->A != nullptr, ExcNotInitialized());
2823 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2824
2825 if (eigenvalues_are_initialized == false)
2826 estimate_eigenvalues(src);
2827
2828 VectorType tmp1, tmp2;
2829
2830 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2831 internal::PreconditionRelaxation::step_operations(
2832 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2833}
2834
2835template <typename MatrixType, typename PreconditionerType>
2836template <typename VectorType>
2839 const VectorType &src) const
2840{
2841 Assert(eigenvalues_are_initialized == false, ExcInternalError());
2842
2843 EigenvalueInformation info;
2844
2845 if (data.relaxation == 0.0)
2846 {
2847 VectorType solution_old, temp_vector1;
2848
2849 solution_old.reinit(src);
2850 temp_vector1.reinit(src, true);
2851
2852 info = internal::estimate_eigenvalues<MatrixType>(
2853 data, A, solution_old, temp_vector1, data.n_iterations);
2854
2855 const double alpha =
2856 (data.smoothing_range > 1. ?
2857 info.max_eigenvalue_estimate / data.smoothing_range :
2858 std::min(0.9 * info.max_eigenvalue_estimate,
2859 info.min_eigenvalue_estimate));
2860
2862 ->data.relaxation = 2.0 / (alpha + info.max_eigenvalue_estimate);
2863 }
2864
2866 ->eigenvalues_are_initialized = true;
2867
2868 return info;
2869}
2870
2871template <typename MatrixType, typename PreconditionerType>
2872double
2874{
2875 return data.relaxation;
2876}
2877
2878
2879//---------------------------------------------------------------------------
2880
2881template <typename MatrixType>
2882inline void
2884 const AdditionalData &parameters_in)
2885{
2886 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2887 Assert(
2888 parameters_in.relaxation != 0.0,
2889 ExcMessage(
2890 "Relaxation cannot automatically be determined by PreconditionJacobi."));
2891
2892 AdditionalData parameters;
2893 parameters.relaxation = 1.0;
2894 parameters.n_iterations = parameters_in.n_iterations;
2895 parameters.preconditioner =
2896 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2897
2898 this->BaseClass::initialize(A, parameters);
2899}
2900
2901//---------------------------------------------------------------------------
2902
2903template <typename MatrixType>
2904inline void
2905PreconditionSOR<MatrixType>::initialize(const MatrixType &A,
2906 const AdditionalData &parameters_in)
2907{
2908 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2909 Assert(
2910 parameters_in.relaxation != 0.0,
2911 ExcMessage(
2912 "Relaxation cannot automatically be determined by PreconditionSOR."));
2913
2914 AdditionalData parameters;
2915 parameters.relaxation = 1.0;
2916 parameters.n_iterations = parameters_in.n_iterations;
2917 parameters.preconditioner =
2918 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2919
2920 this->BaseClass::initialize(A, parameters);
2921}
2922
2923//---------------------------------------------------------------------------
2924
2925template <typename MatrixType>
2926inline void
2928 const AdditionalData &parameters_in)
2929{
2930 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2931 Assert(
2932 parameters_in.relaxation != 0.0,
2933 ExcMessage(
2934 "Relaxation cannot automatically be determined by PreconditionSSOR."));
2935
2936 AdditionalData parameters;
2937 parameters.relaxation = 1.0;
2938 parameters.n_iterations = parameters_in.n_iterations;
2939 parameters.preconditioner =
2940 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2941
2942 this->BaseClass::initialize(A, parameters);
2943}
2944
2945
2946
2947//---------------------------------------------------------------------------
2948
2949template <typename MatrixType>
2950inline void
2952 const MatrixType &A,
2953 const std::vector<size_type> &p,
2954 const std::vector<size_type> &ip,
2955 const typename BaseClass::AdditionalData &parameters_in)
2956{
2957 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2958 Assert(
2959 parameters_in.relaxation != 0.0,
2960 ExcMessage(
2961 "Relaxation cannot automatically be determined by PreconditionPSOR."));
2962
2963 typename BaseClass::AdditionalData parameters;
2964 parameters.relaxation = 1.0;
2965 parameters.n_iterations = parameters_in.n_iterations;
2966 parameters.preconditioner =
2967 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2968
2969 this->BaseClass::initialize(A, parameters);
2970}
2971
2972
2973template <typename MatrixType>
2974inline void
2976 const AdditionalData &additional_data)
2977{
2978 initialize(A,
2979 additional_data.permutation,
2980 additional_data.inverse_permutation,
2981 additional_data.parameters);
2982}
2983
2984template <typename MatrixType>
2986 const std::vector<size_type> &permutation,
2987 const std::vector<size_type> &inverse_permutation,
2989 &parameters)
2990 : permutation(permutation)
2991 , inverse_permutation(inverse_permutation)
2992 , parameters(parameters)
2993{}
2994
2995
2996//---------------------------------------------------------------------------
2997
2998
2999template <typename MatrixType, typename VectorType>
3001 const MatrixType &M,
3002 const function_ptr method)
3003 : matrix(M)
3004 , precondition(method)
3005{}
3006
3007
3008
3009template <typename MatrixType, typename VectorType>
3010void
3012 VectorType &dst,
3013 const VectorType &src) const
3014{
3015 (matrix.*precondition)(dst, src);
3016}
3017
3018//---------------------------------------------------------------------------
3019
3020namespace internal
3021{
3022
3023 template <typename PreconditionerType>
3026 const double smoothing_range,
3027 const unsigned int eig_cg_n_iterations,
3028 const double eig_cg_residual,
3029 const double max_eigenvalue,
3030 const EigenvalueAlgorithm eigenvalue_algorithm)
3031 : smoothing_range(smoothing_range)
3032 , eig_cg_n_iterations(eig_cg_n_iterations)
3033 , eig_cg_residual(eig_cg_residual)
3034 , max_eigenvalue(max_eigenvalue)
3035 , eigenvalue_algorithm(eigenvalue_algorithm)
3036 {}
3037
3038
3039
3040 template <typename PreconditionerType>
3041 inline EigenvalueAlgorithmAdditionalData<PreconditionerType> &
3042 EigenvalueAlgorithmAdditionalData<PreconditionerType>::operator=(
3043 const EigenvalueAlgorithmAdditionalData &other_data)
3044 {
3045 smoothing_range = other_data.smoothing_range;
3046 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3047 eig_cg_residual = other_data.eig_cg_residual;
3048 max_eigenvalue = other_data.max_eigenvalue;
3049 preconditioner = other_data.preconditioner;
3050 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3051 constraints.copy_from(other_data.constraints);
3052
3053 return *this;
3054 }
3055} // namespace internal
3056
3057template <typename MatrixType, typename PreconditionerType>
3059 AdditionalData(const double relaxation,
3060 const unsigned int n_iterations,
3061 const double smoothing_range,
3062 const unsigned int eig_cg_n_iterations,
3063 const double eig_cg_residual,
3064 const double max_eigenvalue,
3065 const EigenvalueAlgorithm eigenvalue_algorithm)
3066 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3067 smoothing_range,
3068 eig_cg_n_iterations,
3069 eig_cg_residual,
3070 max_eigenvalue,
3071 eigenvalue_algorithm)
3072 , relaxation(relaxation)
3073 , n_iterations(n_iterations)
3074{}
3075
3076
3077
3078//---------------------------------------------------------------------------
3079
3080namespace internal
3081{
3082 namespace PreconditionChebyshevImplementation
3083 {
3084 // for deal.II vectors, perform updates for Chebyshev preconditioner all
3085 // at once to reduce memory transfer. Here, we select between general
3086 // vectors and deal.II vectors where we expand the loop over the (local)
3087 // size of the vector
3088
3089 // generic part for non-deal.II vectors
3090 template <typename VectorType, typename PreconditionerType>
3091 inline void
3092 vector_updates(const VectorType &rhs,
3093 const PreconditionerType &preconditioner,
3094 const unsigned int iteration_index,
3095 const double factor1,
3096 const double factor2,
3097 VectorType &solution_old,
3098 VectorType &temp_vector1,
3099 VectorType &temp_vector2,
3100 VectorType &solution)
3101 {
3102 if (iteration_index == 0)
3103 {
3104 solution.equ(factor2, rhs);
3105 preconditioner.vmult(solution_old, solution);
3106 }
3107 else if (iteration_index == 1)
3108 {
3109 // compute t = P^{-1} * (b-A*x^{n})
3110 temp_vector1.sadd(-1.0, 1.0, rhs);
3111 preconditioner.vmult(solution_old, temp_vector1);
3112
3113 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
3114 solution_old.sadd(factor2, 1 + factor1, solution);
3115 }
3116 else
3117 {
3118 // compute t = P^{-1} * (b-A*x^{n})
3119 temp_vector1.sadd(-1.0, 1.0, rhs);
3120 preconditioner.vmult(temp_vector2, temp_vector1);
3121
3122 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
3123 solution_old.sadd(-factor1, factor2, temp_vector2);
3124 solution_old.add(1 + factor1, solution);
3125 }
3126
3127 solution.swap(solution_old);
3128 }
3129
3130 // generic part for deal.II vectors
3131 template <
3132 typename Number,
3133 typename PreconditionerType,
3134 std::enable_if_t<
3135 !has_vmult_with_std_functions_for_precondition<
3136 PreconditionerType,
3138 int> * = nullptr>
3139 inline void
3140 vector_updates(
3142 const PreconditionerType &preconditioner,
3143 const unsigned int iteration_index,
3144 const double factor1_,
3145 const double factor2_,
3147 &solution_old,
3149 &temp_vector1,
3151 &temp_vector2,
3153 {
3154 const Number factor1 = factor1_;
3155 const Number factor1_plus_1 = 1. + factor1_;
3156 const Number factor2 = factor2_;
3157
3158 if (iteration_index == 0)
3159 {
3160 const auto solution_old_ptr = solution_old.begin();
3161
3162 // compute t = P^{-1} * (b)
3163 preconditioner.vmult(solution_old, rhs);
3164
3165 // compute x^{n+1} = f_2 * t
3167 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3168 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
3169 }
3170 else if (iteration_index == 1)
3171 {
3172 const auto solution_ptr = solution.begin();
3173 const auto solution_old_ptr = solution_old.begin();
3174
3175 // compute t = P^{-1} * (b-A*x^{n})
3176 temp_vector1.sadd(-1.0, 1.0, rhs);
3177
3178 preconditioner.vmult(solution_old, temp_vector1);
3179
3180 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
3182 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3183 solution_old_ptr[i] =
3184 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
3185 }
3186 else
3187 {
3188 const auto solution_ptr = solution.begin();
3189 const auto solution_old_ptr = solution_old.begin();
3190 const auto temp_vector2_ptr = temp_vector2.begin();
3191
3192 // compute t = P^{-1} * (b-A*x^{n})
3193 temp_vector1.sadd(-1.0, 1.0, rhs);
3194
3195 preconditioner.vmult(temp_vector2, temp_vector1);
3196
3197 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
3199 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3200 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3201 factor1 * solution_old_ptr[i] +
3202 temp_vector2_ptr[i] * factor2;
3203 }
3204
3205 solution.swap(solution_old);
3206 }
3207
3208 template <
3209 typename Number,
3210 typename PreconditionerType,
3211 std::enable_if_t<
3212 has_vmult_with_std_functions_for_precondition<
3213 PreconditionerType,
3215 int> * = nullptr>
3216 inline void
3217 vector_updates(
3219 const PreconditionerType &preconditioner,
3220 const unsigned int iteration_index,
3221 const double factor1_,
3222 const double factor2_,
3224 &solution_old,
3226 &temp_vector1,
3228 &temp_vector2,
3230 {
3231 const Number factor1 = factor1_;
3232 const Number factor1_plus_1 = 1. + factor1_;
3233 const Number factor2 = factor2_;
3234
3235 const auto rhs_ptr = rhs.begin();
3236 const auto temp_vector1_ptr = temp_vector1.begin();
3237 const auto temp_vector2_ptr = temp_vector2.begin();
3238 const auto solution_ptr = solution.begin();
3239 const auto solution_old_ptr = solution_old.begin();
3240
3241 if (iteration_index == 0)
3242 {
3243 preconditioner.vmult(
3244 solution,
3245 rhs,
3246 [&](const auto start_range, const auto end_range) {
3247 if (end_range > start_range)
3248 std::memset(solution.begin() + start_range,
3249 0,
3250 sizeof(Number) * (end_range - start_range));
3251 },
3252 [&](const auto begin, const auto end) {
3254 for (std::size_t i = begin; i < end; ++i)
3255 solution_ptr[i] *= factor2;
3256 });
3257 }
3258 else
3259 {
3260 preconditioner.vmult(
3261 temp_vector2,
3262 temp_vector1,
3263 [&](const auto begin, const auto end) {
3264 if (end > begin)
3265 std::memset(temp_vector2.begin() + begin,
3266 0,
3267 sizeof(Number) * (end - begin));
3268
3270 for (std::size_t i = begin; i < end; ++i)
3271 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3272 },
3273 [&](const auto begin, const auto end) {
3274 if (iteration_index == 1)
3275 {
3277 for (std::size_t i = begin; i < end; ++i)
3278 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
3279 factor2 * temp_vector2_ptr[i];
3280 }
3281 else
3282 {
3284 for (std::size_t i = begin; i < end; ++i)
3285 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3286 factor1 * solution_old_ptr[i] +
3287 factor2 * temp_vector2_ptr[i];
3288 }
3289 });
3290 }
3291
3292 if (iteration_index > 0)
3293 {
3294 solution_old.swap(temp_vector2);
3295 solution_old.swap(solution);
3296 }
3297 }
3298
3299 // worker routine for deal.II vectors. Because of vectorization, we need
3300 // to put the loop into an extra structure because the virtual function of
3301 // VectorUpdatesRange prevents the compiler from applying vectorization.
3302 template <typename Number>
3303 struct VectorUpdater
3304 {
3305 VectorUpdater(const Number *rhs,
3306 const Number *matrix_diagonal_inverse,
3307 const unsigned int iteration_index,
3308 const Number factor1,
3309 const Number factor2,
3310 Number *solution_old,
3311 Number *tmp_vector,
3312 Number *solution)
3313 : rhs(rhs)
3314 , matrix_diagonal_inverse(matrix_diagonal_inverse)
3315 , iteration_index(iteration_index)
3316 , factor1(factor1)
3317 , factor2(factor2)
3318 , solution_old(solution_old)
3319 , tmp_vector(tmp_vector)
3320 , solution(solution)
3321 {}
3322
3323 void
3324 apply_to_subrange(const std::size_t begin, const std::size_t end) const
3325 {
3326 // To circumvent a bug in gcc
3327 // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
3328 // copies of the variables factor1 and factor2 and do not check based on
3329 // factor1.
3330 const Number factor1 = this->factor1;
3331 const Number factor1_plus_1 = 1. + this->factor1;
3332 const Number factor2 = this->factor2;
3333 if (iteration_index == 0)
3334 {
3336 for (std::size_t i = begin; i < end; ++i)
3337 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
3338 }
3339 else if (iteration_index == 1)
3340 {
3341 // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
3343 for (std::size_t i = begin; i < end; ++i)
3344 // for efficiency reason, write back to temp_vector that is
3345 // already read (avoid read-for-ownership)
3346 tmp_vector[i] =
3347 factor1_plus_1 * solution[i] +
3348 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3349 }
3350 else
3351 {
3352 // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
3353 // + f_2 * P^{-1} * (b-A*x^{n})
3355 for (std::size_t i = begin; i < end; ++i)
3356 // for efficiency reason, write back to temp_vector, which is
3357 // already modified during vmult (in best case, the modified
3358 // values are not written back to main memory yet so that
3359 // we do not have to pay additional costs for writing and
3360 // read-for-ownershop)
3361 tmp_vector[i] =
3362 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
3363 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3364 }
3365 }
3366
3367 const Number *rhs;
3368 const Number *matrix_diagonal_inverse;
3369 const unsigned int iteration_index;
3370 const Number factor1;
3371 const Number factor2;
3372 mutable Number *solution_old;
3373 mutable Number *tmp_vector;
3374 mutable Number *solution;
3375 };
3376
3377 template <typename Number>
3378 struct VectorUpdatesRange : public ::parallel::ParallelForInteger
3379 {
3380 VectorUpdatesRange(const VectorUpdater<Number> &updater,
3381 const std::size_t size)
3382 : updater(updater)
3383 {
3385 VectorUpdatesRange::apply_to_subrange(0, size);
3386 else
3387 apply_parallel(
3388 0,
3389 size,
3391 }
3392
3393 ~VectorUpdatesRange() override = default;
3394
3395 virtual void
3396 apply_to_subrange(const std::size_t begin,
3397 const std::size_t end) const override
3398 {
3399 updater.apply_to_subrange(begin, end);
3400 }
3401
3402 const VectorUpdater<Number> &updater;
3403 };
3404
3405 // selection for diagonal matrix around deal.II vector
3406 template <typename Number>
3407 inline void
3408 vector_updates(
3409 const ::Vector<Number> &rhs,
3410 const ::DiagonalMatrix<::Vector<Number>> &jacobi,
3411 const unsigned int iteration_index,
3412 const double factor1,
3413 const double factor2,
3414 ::Vector<Number> &solution_old,
3415 ::Vector<Number> &temp_vector1,
3417 ::Vector<Number> &solution)
3418 {
3419 VectorUpdater<Number> upd(rhs.begin(),
3420 jacobi.get_vector().begin(),
3421 iteration_index,
3422 factor1,
3423 factor2,
3424 solution_old.begin(),
3425 temp_vector1.begin(),
3426 solution.begin());
3427 VectorUpdatesRange<Number>(upd, rhs.size());
3428
3429 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3430 if (iteration_index == 0)
3431 {
3432 // nothing to do here because we can immediately write into the
3433 // solution vector without remembering any of the other vectors
3434 }
3435 else
3436 {
3437 solution.swap(temp_vector1);
3438 solution_old.swap(temp_vector1);
3439 }
3440 }
3441
3442 // selection for diagonal matrix around parallel deal.II vector
3443 template <typename Number>
3444 inline void
3445 vector_updates(
3447 const ::DiagonalMatrix<
3449 const unsigned int iteration_index,
3450 const double factor1,
3451 const double factor2,
3453 &solution_old,
3455 &temp_vector1,
3458 {
3459 VectorUpdater<Number> upd(rhs.begin(),
3460 jacobi.get_vector().begin(),
3461 iteration_index,
3462 factor1,
3463 factor2,
3464 solution_old.begin(),
3465 temp_vector1.begin(),
3466 solution.begin());
3467 VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
3468
3469 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3470 if (iteration_index == 0)
3471 {
3472 // nothing to do here because we can immediately write into the
3473 // solution vector without remembering any of the other vectors
3474 }
3475 else
3476 {
3477 solution.swap(temp_vector1);
3478 solution_old.swap(temp_vector1);
3479 }
3480 }
3481
3482 // We need to have a separate declaration for static const members
3483
3484 // general case and the case that the preconditioner can work on
3485 // ranges (covered by vector_updates())
3486 template <
3487 typename MatrixType,
3488 typename VectorType,
3489 typename PreconditionerType,
3490 std::enable_if_t<
3491 !has_vmult_with_std_functions<MatrixType,
3492 VectorType,
3493 PreconditionerType> &&
3494 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3495 VectorType> &&
3496 has_vmult_with_std_functions_for_precondition<MatrixType,
3497 VectorType>),
3498 int> * = nullptr>
3499 inline void
3500 vmult_and_update(const MatrixType &matrix,
3501 const PreconditionerType &preconditioner,
3502 const VectorType &rhs,
3503 const unsigned int iteration_index,
3504 const double factor1,
3505 const double factor2,
3506 VectorType &solution,
3507 VectorType &solution_old,
3508 VectorType &temp_vector1,
3509 VectorType &temp_vector2)
3510 {
3511 if (iteration_index > 0)
3512 matrix.vmult(temp_vector1, solution);
3513 vector_updates(rhs,
3514 preconditioner,
3515 iteration_index,
3516 factor1,
3517 factor2,
3518 solution_old,
3519 temp_vector1,
3520 temp_vector2,
3521 solution);
3522 }
3523
3524 // case that both the operator and the preconditioner can work on
3525 // subranges
3526 template <
3527 typename MatrixType,
3528 typename VectorType,
3529 typename PreconditionerType,
3530 std::enable_if_t<
3531 !has_vmult_with_std_functions<MatrixType,
3532 VectorType,
3533 PreconditionerType> &&
3534 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3535 VectorType> &&
3536 has_vmult_with_std_functions_for_precondition<MatrixType,
3537 VectorType>),
3538 int> * = nullptr>
3539 inline void
3540 vmult_and_update(const MatrixType &matrix,
3541 const PreconditionerType &preconditioner,
3542 const VectorType &rhs,
3543 const unsigned int iteration_index,
3544 const double factor1_,
3545 const double factor2_,
3546 VectorType &solution,
3547 VectorType &solution_old,
3548 VectorType &temp_vector1,
3549 VectorType &temp_vector2)
3550 {
3551 using Number = typename VectorType::value_type;
3552
3553 const Number factor1 = factor1_;
3554 const Number factor1_plus_1 = 1. + factor1_;
3555 const Number factor2 = factor2_;
3556
3557 if (iteration_index == 0)
3558 {
3559 preconditioner.vmult(
3560 solution,
3561 rhs,
3562 [&](const unsigned int start_range, const unsigned int end_range) {
3563 // zero 'solution' before running the vmult operation
3564 if (end_range > start_range)
3565 std::memset(solution.begin() + start_range,
3566 0,
3567 sizeof(Number) * (end_range - start_range));
3568 },
3569 [&](const unsigned int start_range, const unsigned int end_range) {
3570 const auto solution_ptr = solution.begin();
3571
3573 for (std::size_t i = start_range; i < end_range; ++i)
3574 solution_ptr[i] *= factor2;
3575 });
3576 }
3577 else
3578 {
3579 temp_vector1.reinit(rhs, true);
3580 temp_vector2.reinit(rhs, true);
3581
3582 // 1) compute residual (including operator application)
3583 matrix.vmult(
3584 temp_vector1,
3585 solution,
3586 [&](const unsigned int start_range, const unsigned int end_range) {
3587 // zero 'temp_vector1' before running the vmult
3588 // operation
3589 if (end_range > start_range)
3590 std::memset(temp_vector1.begin() + start_range,
3591 0,
3592 sizeof(Number) * (end_range - start_range));
3593 },
3594 [&](const unsigned int start_range, const unsigned int end_range) {
3595 const auto rhs_ptr = rhs.begin();
3596 const auto tmp_ptr = temp_vector1.begin();
3597
3599 for (std::size_t i = start_range; i < end_range; ++i)
3600 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3601 });
3602
3603 // 2) perform vector updates (including preconditioner application)
3604 preconditioner.vmult(
3605 temp_vector2,
3606 temp_vector1,
3607 [&](const unsigned int start_range, const unsigned int end_range) {
3608 // zero 'temp_vector2' before running the vmult
3609 // operation
3610 if (end_range > start_range)
3611 std::memset(temp_vector2.begin() + start_range,
3612 0,
3613 sizeof(Number) * (end_range - start_range));
3614 },
3615 [&](const unsigned int start_range, const unsigned int end_range) {
3616 const auto solution_ptr = solution.begin();
3617 const auto solution_old_ptr = solution_old.begin();
3618 const auto tmp_ptr = temp_vector2.begin();
3619
3620 if (iteration_index == 1)
3621 {
3623 for (std::size_t i = start_range; i < end_range; ++i)
3624 tmp_ptr[i] =
3625 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3626 }
3627 else
3628 {
3630 for (std::size_t i = start_range; i < end_range; ++i)
3631 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3632 factor1 * solution_old_ptr[i] +
3633 factor2 * tmp_ptr[i];
3634 }
3635 });
3636
3637 solution.swap(temp_vector2);
3638 solution_old.swap(temp_vector2);
3639 }
3640 }
3641
3642 // case that the operator can work on subranges and the preconditioner
3643 // is a diagonal
3644 template <typename MatrixType,
3645 typename VectorType,
3646 typename PreconditionerType,
3647 std::enable_if_t<has_vmult_with_std_functions<MatrixType,
3648 VectorType,
3649 PreconditionerType>,
3650 int> * = nullptr>
3651 inline void
3652 vmult_and_update(const MatrixType &matrix,
3653 const PreconditionerType &preconditioner,
3654 const VectorType &rhs,
3655 const unsigned int iteration_index,
3656 const double factor1,
3657 const double factor2,
3658 VectorType &solution,
3659 VectorType &solution_old,
3660 VectorType &temp_vector1,
3661 VectorType &)
3662 {
3663 using Number = typename VectorType::value_type;
3664 VectorUpdater<Number> updater(rhs.begin(),
3665 preconditioner.get_vector().begin(),
3666 iteration_index,
3667 factor1,
3668 factor2,
3669 solution_old.begin(),
3670 temp_vector1.begin(),
3671 solution.begin());
3672 if (iteration_index > 0)
3673 matrix.vmult(
3674 temp_vector1,
3675 solution,
3676 [&](const unsigned int start_range, const unsigned int end_range) {
3677 // zero 'temp_vector1' before running the vmult
3678 // operation
3679 if (end_range > start_range)
3680 std::memset(temp_vector1.begin() + start_range,
3681 0,
3682 sizeof(Number) * (end_range - start_range));
3683 },
3684 [&](const unsigned int start_range, const unsigned int end_range) {
3685 if (end_range > start_range)
3686 updater.apply_to_subrange(start_range, end_range);
3687 });
3688 else
3689 updater.apply_to_subrange(0U, solution.locally_owned_size());
3690
3691 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3692 if (iteration_index == 0)
3693 {
3694 // nothing to do here because we can immediately write into the
3695 // solution vector without remembering any of the other vectors
3696 }
3697 else
3698 {
3699 solution.swap(temp_vector1);
3700 solution_old.swap(temp_vector1);
3701 }
3702 }
3703
3704 template <typename MatrixType, typename PreconditionerType>
3705 inline void
3706 initialize_preconditioner(
3707 const MatrixType &matrix,
3708 std::shared_ptr<PreconditionerType> &preconditioner)
3709 {
3710 (void)matrix;
3711 (void)preconditioner;
3712 AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
3713 }
3714
3715 template <typename MatrixType, typename VectorType>
3716 inline void
3717 initialize_preconditioner(
3718 const MatrixType &matrix,
3719 std::shared_ptr<::DiagonalMatrix<VectorType>> &preconditioner)
3720 {
3721 if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
3722 {
3723 if (preconditioner.get() == nullptr)
3724 preconditioner =
3725 std::make_shared<::DiagonalMatrix<VectorType>>();
3726
3727 Assert(
3728 preconditioner->m() == 0,
3729 ExcMessage(
3730 "Preconditioner appears to be initialized but not sized correctly"));
3731
3732 // This part only works in serial
3733 if (preconditioner->m() != matrix.m())
3734 {
3735 preconditioner->get_vector().reinit(matrix.m());
3736 for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
3737 preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
3738 }
3739 }
3740 }
3741 } // namespace PreconditionChebyshevImplementation
3742} // namespace internal
3743
3744
3745
3746template <typename MatrixType, typename VectorType, typename PreconditionerType>
3748 AdditionalData::AdditionalData(const unsigned int degree,
3749 const double smoothing_range,
3750 const unsigned int eig_cg_n_iterations,
3751 const double eig_cg_residual,
3752 const double max_eigenvalue,
3753 const EigenvalueAlgorithm eigenvalue_algorithm,
3754 const PolynomialType polynomial_type)
3755 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3756 smoothing_range,
3757 eig_cg_n_iterations,
3758 eig_cg_residual,
3759 max_eigenvalue,
3760 eigenvalue_algorithm)
3761 , degree(degree)
3762 , polynomial_type(polynomial_type)
3763{}
3764
3765
3766
3767template <typename MatrixType, typename VectorType, typename PreconditionerType>
3770 : theta(1.)
3771 , delta(1.)
3772 , eigenvalues_are_initialized(false)
3773{
3774 static_assert(
3775 std::is_same_v<size_type, typename VectorType::size_type>,
3776 "PreconditionChebyshev and VectorType must have the same size_type.");
3777}
3778
3779
3780
3781template <typename MatrixType, typename VectorType, typename PreconditionerType>
3782inline void
3784 const MatrixType &matrix,
3785 const AdditionalData &additional_data)
3786{
3787 matrix_ptr = &matrix;
3788 data = additional_data;
3789 Assert(data.degree > 0,
3790 ExcMessage("The degree of the Chebyshev method must be positive."));
3791 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3792 matrix, data.preconditioner);
3793 eigenvalues_are_initialized = false;
3794}
3795
3796
3797
3798template <typename MatrixType, typename VectorType, typename PreconditionerType>
3799inline void
3801{
3802 eigenvalues_are_initialized = false;
3803 theta = delta = 1.0;
3804 matrix_ptr = nullptr;
3805 {
3806 VectorType empty_vector;
3807 solution_old.reinit(empty_vector);
3808 temp_vector1.reinit(empty_vector);
3809 temp_vector2.reinit(empty_vector);
3810 }
3811 data.preconditioner.reset();
3812}
3813
3814
3815
3816template <typename MatrixType, typename VectorType, typename PreconditionerType>
3817inline typename internal::EigenvalueInformation
3819 estimate_eigenvalues(const VectorType &src) const
3820{
3821 Assert(eigenvalues_are_initialized == false, ExcInternalError());
3822
3823 solution_old.reinit(src);
3824 temp_vector1.reinit(src, true);
3825
3826 auto info = internal::estimate_eigenvalues<MatrixType>(
3827 data, matrix_ptr, solution_old, temp_vector1, data.degree);
3828
3829 const double alpha = (data.smoothing_range > 1. ?
3830 info.max_eigenvalue_estimate / data.smoothing_range :
3831 std::min(0.9 * info.max_eigenvalue_estimate,
3832 info.min_eigenvalue_estimate));
3833
3834 // in case the user set the degree to invalid unsigned int, we have to
3835 // determine the number of necessary iterations from the Chebyshev error
3836 // estimate, given the target tolerance specified by smoothing_range. This
3837 // estimate is based on the error formula given in section 5.1 of
3838 // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
3839 if (data.degree == numbers::invalid_unsigned_int)
3840 {
3841 const double actual_range = info.max_eigenvalue_estimate / alpha;
3842 const double sigma = (1. - std::sqrt(1. / actual_range)) /
3843 (1. + std::sqrt(1. / actual_range));
3844 const double eps = data.smoothing_range;
3845 const_cast<
3847 this)
3848 ->data.degree =
3849 1 + static_cast<unsigned int>(
3850 std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
3851 std::log(1. / sigma));
3852 }
3853
3854 info.degree = data.degree;
3855
3856 const_cast<
3858 ->delta =
3859 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3860 (info.max_eigenvalue_estimate) :
3861 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3862 const_cast<
3864 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3865
3866 // We do not need the second temporary vector in case we have a
3867 // DiagonalMatrix as preconditioner and use deal.II's own vectors
3868 using NumberType = typename VectorType::value_type;
3869 if (std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> ==
3870 false ||
3871 (std::is_same_v<VectorType, ::Vector<NumberType>> == false &&
3872 ((std::is_same_v<
3873 VectorType,
3875 false) ||
3876 (std::is_same_v<VectorType,
3878 Vector<NumberType, MemorySpace::Default>> == false))))
3879 temp_vector2.reinit(src, true);
3880 else
3881 {
3882 VectorType empty_vector;
3883 temp_vector2.reinit(empty_vector);
3884 }
3885
3886 const_cast<
3888 ->eigenvalues_are_initialized = true;
3889
3890 return info;
3891}
3892
3893
3894
3895template <typename MatrixType, typename VectorType, typename PreconditionerType>
3896inline void
3898 VectorType &solution,
3899 const VectorType &rhs) const
3900{
3901 std::lock_guard<std::mutex> lock(mutex);
3902 if (eigenvalues_are_initialized == false)
3903 estimate_eigenvalues(rhs);
3904
3905 internal::PreconditionChebyshevImplementation::vmult_and_update(
3906 *matrix_ptr,
3907 *data.preconditioner,
3908 rhs,
3909 0,
3910 0.,
3911 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3912 (4. / (3. * delta)) :
3913 (1. / theta),
3914 solution,
3915 solution_old,
3916 temp_vector1,
3917 temp_vector2);
3918
3919 // if delta is zero, we do not need to iterate because the updates will be
3920 // zero
3921 if (data.degree < 2 || std::abs(delta) < 1e-40)
3922 return;
3923
3924 double rhok = delta / theta, sigma = theta / delta;
3925 for (unsigned int k = 0; k < data.degree - 1; ++k)
3926 {
3927 double factor1 = 0.0;
3928 double factor2 = 0.0;
3929
3930 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3931 {
3932 factor1 = (2 * k + 1.) / (2 * k + 5.);
3933 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3934 }
3935 else
3936 {
3937 const double rhokp = 1. / (2. * sigma - rhok);
3938 factor1 = rhokp * rhok;
3939 factor2 = 2. * rhokp / delta;
3940 rhok = rhokp;
3941 }
3942
3943 internal::PreconditionChebyshevImplementation::vmult_and_update(
3944 *matrix_ptr,
3945 *data.preconditioner,
3946 rhs,
3947 k + 1,
3948 factor1,
3949 factor2,
3950 solution,
3951 solution_old,
3952 temp_vector1,
3953 temp_vector2);
3954 }
3955}
3956
3957
3958
3959template <typename MatrixType, typename VectorType, typename PreconditionerType>
3960inline void
3962 VectorType &solution,
3963 const VectorType &rhs) const
3964{
3965 std::lock_guard<std::mutex> lock(mutex);
3966 if (eigenvalues_are_initialized == false)
3967 estimate_eigenvalues(rhs);
3968
3969 internal::PreconditionChebyshevImplementation::vector_updates(
3970 rhs,
3971 *data.preconditioner,
3972 0,
3973 0.,
3974 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3975 (4. / (3. * delta)) :
3976 (1. / theta),
3977 solution_old,
3978 temp_vector1,
3979 temp_vector2,
3980 solution);
3981
3982 if (data.degree < 2 || std::abs(delta) < 1e-40)
3983 return;
3984
3985 double rhok = delta / theta, sigma = theta / delta;
3986 for (unsigned int k = 0; k < data.degree - 1; ++k)
3987 {
3988 double factor1 = 0.0;
3989 double factor2 = 0.0;
3990
3991 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3992 {
3993 factor1 = (2 * k + 1.) / (2 * k + 5.);
3994 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3995 }
3996 else
3997 {
3998 const double rhokp = 1. / (2. * sigma - rhok);
3999 factor1 = rhokp * rhok;
4000 factor2 = 2. * rhokp / delta;
4001 rhok = rhokp;
4002 }
4003
4004 matrix_ptr->Tvmult(temp_vector1, solution);
4005 internal::PreconditionChebyshevImplementation::vector_updates(
4006 rhs,
4007 *data.preconditioner,
4008 k + 1,
4009 factor1,
4010 factor2,
4011 solution_old,
4012 temp_vector1,
4013 temp_vector2,
4014 solution);
4015 }
4016}
4017
4018
4019
4020template <typename MatrixType, typename VectorType, typename PreconditionerType>
4021inline void
4023 VectorType &solution,
4024 const VectorType &rhs) const
4025{
4026 std::lock_guard<std::mutex> lock(mutex);
4027 if (eigenvalues_are_initialized == false)
4028 estimate_eigenvalues(rhs);
4029
4030 internal::PreconditionChebyshevImplementation::vmult_and_update(
4031 *matrix_ptr,
4032 *data.preconditioner,
4033 rhs,
4034 1,
4035 0.,
4036 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4037 (4. / (3. * delta)) :
4038 (1. / theta),
4039 solution,
4040 solution_old,
4041 temp_vector1,
4042 temp_vector2);
4043
4044 if (data.degree < 2 || std::abs(delta) < 1e-40)
4045 return;
4046
4047 double rhok = delta / theta, sigma = theta / delta;
4048 for (unsigned int k = 0; k < data.degree - 1; ++k)
4049 {
4050 double factor1 = 0.0;
4051 double factor2 = 0.0;
4052
4053 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4054 {
4055 factor1 = (2 * k + 1.) / (2 * k + 5.);
4056 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4057 }
4058 else
4059 {
4060 const double rhokp = 1. / (2. * sigma - rhok);
4061 factor1 = rhokp * rhok;
4062 factor2 = 2. * rhokp / delta;
4063 rhok = rhokp;
4064 }
4065
4066 internal::PreconditionChebyshevImplementation::vmult_and_update(
4067 *matrix_ptr,
4068 *data.preconditioner,
4069 rhs,
4070 k + 2,
4071 factor1,
4072 factor2,
4073 solution,
4074 solution_old,
4075 temp_vector1,
4076 temp_vector2);
4077 }
4078}
4079
4080
4081
4082template <typename MatrixType, typename VectorType, typename PreconditionerType>
4083inline void
4085 VectorType &solution,
4086 const VectorType &rhs) const
4087{
4088 std::lock_guard<std::mutex> lock(mutex);
4089 if (eigenvalues_are_initialized == false)
4090 estimate_eigenvalues(rhs);
4091
4092 matrix_ptr->Tvmult(temp_vector1, solution);
4093 internal::PreconditionChebyshevImplementation::vector_updates(
4094 rhs,
4095 *data.preconditioner,
4096 1,
4097 0.,
4098 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4099 (4. / (3. * delta)) :
4100 (1. / theta),
4101 solution_old,
4102 temp_vector1,
4103 temp_vector2,
4104 solution);
4105
4106 if (data.degree < 2 || std::abs(delta) < 1e-40)
4107 return;
4108
4109 double rhok = delta / theta, sigma = theta / delta;
4110 for (unsigned int k = 0; k < data.degree - 1; ++k)
4111 {
4112 double factor1 = 0.0;
4113 double factor2 = 0.0;
4114
4115 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4116 {
4117 factor1 = (2 * k + 1.) / (2 * k + 5.);
4118 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4119 }
4120 else
4121 {
4122 const double rhokp = 1. / (2. * sigma - rhok);
4123 factor1 = rhokp * rhok;
4124 factor2 = 2. * rhokp / delta;
4125 rhok = rhokp;
4126 }
4127
4128 matrix_ptr->Tvmult(temp_vector1, solution);
4129 internal::PreconditionChebyshevImplementation::vector_updates(
4130 rhs,
4131 *data.preconditioner,
4132 k + 2,
4133 factor1,
4134 factor2,
4135 solution_old,
4136 temp_vector1,
4137 temp_vector2,
4138 solution);
4139 }
4140}
4141
4142
4143
4144template <typename MatrixType, typename VectorType, typename PreconditionerType>
4145inline typename PreconditionChebyshev<MatrixType,
4146 VectorType,
4147 PreconditionerType>::size_type
4149{
4150 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4151 return matrix_ptr->m();
4152}
4153
4154
4155
4156template <typename MatrixType, typename VectorType, typename PreconditionerType>
4157inline typename PreconditionChebyshev<MatrixType,
4158 VectorType,
4159 PreconditionerType>::size_type
4161{
4162 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4163 return matrix_ptr->n();
4164}
4165
4166#endif // DOXYGEN
4167
4169
4170#endif
unsigned int n_blocks() const
BlockType & block(const unsigned int i)
bool is_empty() const
Definition index_set.h:1926
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1982
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()
@ 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)
static const unsigned int invalid_unsigned_int
Definition types.h:220
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:81
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)