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