deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+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: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 1999 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#ifndef dealii_precondition_h
14#define dealii_precondition_h
15
16#include <deal.II/base/config.h>
17
18#include <deal.II/base/memory_space.h>
19#include <deal.II/base/mutex.h>
20#include <deal.II/base/observer_pointer.h>
21#include <deal.II/base/parallel.h>
22#include <deal.II/base/template_constraints.h>
23
24#include <deal.II/lac/affine_constraints.h>
25#include <deal.II/lac/block_vector_base.h>
26#include <deal.II/lac/diagonal_matrix.h>
27#include <deal.II/lac/identity_matrix.h>
28#include <deal.II/lac/solver_cg.h>
29#include <deal.II/lac/vector_memory.h>
30
31#include <Kokkos_Core.hpp>
32
33#include <limits>
34
35DEAL_II_NAMESPACE_OPEN
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, typename>
50 class BlockVector;
51 } // namespace distributed
52} // namespace LinearAlgebra
53#endif
54
55
61namespace internal
62{
67 enum class EigenvalueAlgorithm
68 {
74 lanczos,
84 power_iteration
85 };
86
91 struct EigenvalueInformation
92 {
96 double min_eigenvalue_estimate;
100 double max_eigenvalue_estimate;
104 unsigned int cg_iterations;
109 unsigned int degree;
113 EigenvalueInformation()
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>
205} // namespace internal
206
207
227class PreconditionIdentity : public EnableObserverPointer
228{
229public:
233 using size_type = types::global_dof_index;
234
239 struct AdditionalData
240 {
244 AdditionalData() = default;
245 };
246
250 PreconditionIdentity();
251
256 template <typename MatrixType>
257 void
258 initialize(const MatrixType &matrix,
259 const AdditionalData &additional_data = AdditionalData());
260
264 template <typename VectorType>
265 void
266 vmult(VectorType &, const VectorType &) const;
267
272 template <typename VectorType>
273 void
274 Tvmult(VectorType &, const VectorType &) const;
275
279 template <typename VectorType>
280 void
281 vmult_add(VectorType &, const VectorType &) const;
282
287 template <typename VectorType>
288 void
289 Tvmult_add(VectorType &, const VectorType &) const;
290
295 void
296 clear();
297
305 size_type
306 m() const;
307
315 size_type
316 n() const;
317
318private:
322 size_type n_rows;
323
327 size_type n_columns;
328};
329
330
331
342class PreconditionRichardson : public EnableObserverPointer
343{
344public:
348 using size_type = types::global_dof_index;
349
353 class AdditionalData
354 {
355 public:
360 AdditionalData(const double relaxation = 1.);
361
365 double relaxation;
366 };
367
372 PreconditionRichardson();
373
377 void
378 initialize(const AdditionalData &parameters);
379
385 template <typename MatrixType>
386 void
387 initialize(const MatrixType &matrix, const AdditionalData &parameters);
388
392 template <typename VectorType>
393 void
394 vmult(VectorType &, const VectorType &) const;
395
400 template <typename VectorType>
401 void
402 Tvmult(VectorType &, const VectorType &) const;
406 template <typename VectorType>
407 void
408 vmult_add(VectorType &, const VectorType &) const;
409
414 template <typename VectorType>
415 void
416 Tvmult_add(VectorType &, const VectorType &) const;
417
422 void
423 clear()
424 {}
425
433 size_type
434 m() const;
435
443 size_type
444 n() const;
445
446private:
450 double relaxation;
451
455 size_type n_rows;
456
460 size_type n_columns;
461};
462
463
464
504template <typename MatrixType = SparseMatrix<double>,
505 typename VectorType = Vector<double>>
506class PreconditionUseMatrix : public EnableObserverPointer
507{
508public:
512 using function_ptr = void (MatrixType::*)(VectorType &,
513 const VectorType &) const;
514
520 PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
521
526 void
527 vmult(VectorType &dst, const VectorType &src) const;
528
529private:
533 const MatrixType &matrix;
534
538 const function_ptr precondition;
539};
540
541
542
572template <typename MatrixType = SparseMatrix<double>,
573 typename PreconditionerType = IdentityMatrix>
574class PreconditionRelaxation : public EnableObserverPointer
575{
576public:
580 using size_type = types::global_dof_index;
581
585 class AdditionalData
586 : public internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>
587 {
588 public:
589 using EigenvalueAlgorithm = internal::EigenvalueAlgorithm;
590
594 AdditionalData(const double relaxation = 1.,
595 const unsigned int n_iterations = 1,
596 const double smoothing_range = 0.,
597 const unsigned int eig_cg_n_iterations = 8,
598 const double eig_cg_residual = 1e-2,
599 const double max_eigenvalue = 1,
600 const EigenvalueAlgorithm eigenvalue_algorithm =
601 EigenvalueAlgorithm::lanczos,
602 const double safety_factor = 1.2);
603
607 double relaxation;
608
613 unsigned int n_iterations;
614 };
615
621 void
622 initialize(const MatrixType &A,
623 const AdditionalData &parameters = AdditionalData());
624
628 void
629 clear();
630
635 size_type
636 m() const;
637
642 size_type
643 n() const;
644
648 template <typename VectorType>
649 void
650 vmult(VectorType &, const VectorType &) const;
651
656 template <typename VectorType>
657 void
658 Tvmult(VectorType &, const VectorType &) const;
659
663 template <typename VectorType>
664 void
665 step(VectorType &x, const VectorType &rhs) const;
666
670 template <typename VectorType>
671 void
672 Tstep(VectorType &x, const VectorType &rhs) const;
673
674 using EigenvalueInformation = internal::EigenvalueInformation;
675
679 template <typename VectorType>
680 EigenvalueInformation
681 estimate_eigenvalues(const VectorType &src) const;
682
688 double
689 get_relaxation() const;
690
691protected:
695 ObserverPointer<const MatrixType, PreconditionRelaxation<MatrixType>> A;
696
701 AdditionalData data;
702
706 std::shared_ptr<PreconditionerType> preconditioner;
707
712 bool eigenvalues_are_initialized;
713};
714
715
716
717#ifndef DOXYGEN
718
719namespace internal
720{
721 // a helper type-trait that leverage SFINAE to figure out if MatrixType has
722 // ... MatrixType::vmult(VectorType &, const VectorType&,
723 // std::function<...>, std::function<...>) const
724 template <typename MatrixType, typename VectorType>
725 using vmult_functions_t = decltype(std::declval<const MatrixType>().vmult(
726 std::declval<VectorType &>(),
727 std::declval<const VectorType &>(),
728 std::declval<
729 const std::function<void(const unsigned int, const unsigned int)> &>(),
730 std::declval<
731 const std::function<void(const unsigned int, const unsigned int)> &>()));
732
733 template <typename MatrixType,
734 typename VectorType,
735 typename PreconditionerType>
736 constexpr bool has_vmult_with_std_functions =
737 is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
738 std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> &&
739 (std::is_same_v<VectorType,
740 ::Vector<typename VectorType::value_type>> ||
741 std::is_same_v<
742 VectorType,
743 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
744 MemorySpace::Host>>);
745
746
747 template <typename MatrixType, typename VectorType>
748 constexpr bool has_vmult_with_std_functions_for_precondition =
749 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
750
751 namespace PreconditionRelaxation
752 {
753 template <typename T, typename VectorType>
754 using Tvmult_t = decltype(std::declval<const T>().Tvmult(
755 std::declval<VectorType &>(),
756 std::declval<const VectorType &>()));
757
758 template <typename T, typename VectorType>
759 constexpr bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
760
761 template <typename T, typename VectorType>
762 using step_t = decltype(std::declval<const T>().step(
763 std::declval<VectorType &>(),
764 std::declval<const VectorType &>()));
765
766 template <typename T, typename VectorType>
767 constexpr bool has_step = is_supported_operation<step_t, T, VectorType>;
768
769 template <typename T, typename VectorType>
770 using step_omega_t =
771 decltype(std::declval<const T>().step(std::declval<VectorType &>(),
772 std::declval<const VectorType &>(),
773 std::declval<const double>()));
774
775 template <typename T, typename VectorType>
776 constexpr bool has_step_omega =
777 is_supported_operation<step_omega_t, T, VectorType>;
778
779 template <typename T, typename VectorType>
780 using Tstep_t = decltype(std::declval<const T>().Tstep(
781 std::declval<VectorType &>(),
782 std::declval<const VectorType &>()));
783
784 template <typename T, typename VectorType>
785 constexpr bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
786
787 template <typename T, typename VectorType>
788 using Tstep_omega_t =
789 decltype(std::declval<const T>().Tstep(std::declval<VectorType &>(),
790 std::declval<const VectorType &>(),
791 std::declval<const double>()));
792
793 template <typename T, typename VectorType>
794 constexpr bool has_Tstep_omega =
795 is_supported_operation<Tstep_omega_t, T, VectorType>;
796
797 template <typename T, typename VectorType>
798 using jacobi_step_t = decltype(std::declval<const T>().Jacobi_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_jacobi_step =
805 is_supported_operation<jacobi_step_t, T, VectorType>;
806
807 template <typename T, typename VectorType>
808 using SOR_step_t = decltype(std::declval<const T>().SOR_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_SOR_step =
815 is_supported_operation<SOR_step_t, T, VectorType>;
816
817 template <typename T, typename VectorType>
818 using SSOR_step_t = decltype(std::declval<const T>().SSOR_step(
819 std::declval<VectorType &>(),
820 std::declval<const VectorType &>(),
821 std::declval<const double>()));
822
823 template <typename T, typename VectorType>
824 constexpr bool has_SSOR_step =
825 is_supported_operation<SSOR_step_t, T, VectorType>;
826
827 template <typename MatrixType>
828 class PreconditionJacobiImpl
829 {
830 public:
831 PreconditionJacobiImpl(const MatrixType &A, const double relaxation)
832 : A(&A)
833 , relaxation(relaxation)
834 {}
835
836 template <typename VectorType>
837 void
838 vmult(VectorType &dst, const VectorType &src) const
839 {
840 this->A->precondition_Jacobi(dst, src, this->relaxation);
841 }
842
843 template <typename VectorType>
844 void
845 Tvmult(VectorType &dst, const VectorType &src) const
846 {
847 // call vmult, since preconditioner is symmetrical
848 this->vmult(dst, src);
849 }
850
851 template <typename VectorType,
852 std::enable_if_t<has_jacobi_step<MatrixType, VectorType>,
853 MatrixType> * = nullptr>
854 void
855 step(VectorType &dst, const VectorType &src) const
856 {
857 this->A->Jacobi_step(dst, src, this->relaxation);
858 }
859
860 template <typename VectorType,
861 std::enable_if_t<!has_jacobi_step<MatrixType, VectorType>,
862 MatrixType> * = nullptr>
863 void
864 step(VectorType &, const VectorType &) const
865 {
866 AssertThrow(false,
867 ExcMessage(
868 "Matrix A does not provide a Jacobi_step() function!"));
869 }
870
871 template <typename VectorType>
872 void
873 Tstep(VectorType &dst, const VectorType &src) const
874 {
875 // call step, since preconditioner is symmetrical
876 this->step(dst, src);
877 }
878
879 private:
880 const ObserverPointer<const MatrixType> A;
881 const double relaxation;
882 };
883
884 template <typename MatrixType>
885 class PreconditionSORImpl
886 {
887 public:
888 PreconditionSORImpl(const MatrixType &A, const double relaxation)
889 : A(&A)
890 , relaxation(relaxation)
891 {}
892
893 template <typename VectorType>
894 void
895 vmult(VectorType &dst, const VectorType &src) const
896 {
897 this->A->precondition_SOR(dst, src, this->relaxation);
898 }
899
900 template <typename VectorType>
901 void
902 Tvmult(VectorType &dst, const VectorType &src) const
903 {
904 this->A->precondition_TSOR(dst, src, this->relaxation);
905 }
906
907 template <typename VectorType,
908 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
909 MatrixType> * = nullptr>
910 void
911 step(VectorType &dst, const VectorType &src) const
912 {
913 this->A->SOR_step(dst, src, this->relaxation);
914 }
915
916 template <typename VectorType,
917 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
918 MatrixType> * = nullptr>
919 void
920 step(VectorType &, const VectorType &) const
921 {
922 AssertThrow(false,
923 ExcMessage(
924 "Matrix A does not provide a SOR_step() function!"));
925 }
926
927 template <typename VectorType,
928 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
929 MatrixType> * = nullptr>
930 void
931 Tstep(VectorType &dst, const VectorType &src) const
932 {
933 this->A->TSOR_step(dst, src, this->relaxation);
934 }
935
936 template <typename VectorType,
937 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
938 MatrixType> * = nullptr>
939 void
940 Tstep(VectorType &, const VectorType &) const
941 {
942 AssertThrow(false,
943 ExcMessage(
944 "Matrix A does not provide a TSOR_step() function!"));
945 }
946
947 private:
948 const ObserverPointer<const MatrixType> A;
949 const double relaxation;
950 };
951
952 template <typename MatrixType>
953 class PreconditionSSORImpl
954 {
955 public:
956 using size_type = typename MatrixType::size_type;
957
958 PreconditionSSORImpl(const MatrixType &A, const double relaxation)
959 : A(&A)
960 , relaxation(relaxation)
961 {
962 // in case we have a SparseMatrix class, we can extract information
963 // about the diagonal.
964 const SparseMatrix<typename MatrixType::value_type> *mat =
965 dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
966 &*this->A);
967
968 // calculate the positions first after the diagonal.
969 if (mat != nullptr)
970 {
971 const size_type n = this->A->n();
972 pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
973 for (size_type row = 0; row < n; ++row)
974 {
975 // find the first element in this line which is on the right of
976 // the diagonal. we need to precondition with the elements on
977 // the left only. note: the first entry in each line denotes the
978 // diagonal element, which we need not check.
979 typename SparseMatrix<
980 typename MatrixType::value_type>::const_iterator it =
981 mat->begin(row) + 1;
982 for (; it < mat->end(row); ++it)
983 if (it->column() > row)
984 break;
985 pos_right_of_diagonal[row] = it - mat->begin();
986 }
987 }
988 }
989
990 template <typename VectorType>
991 void
992 vmult(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 void
1002 Tvmult(VectorType &dst, const VectorType &src) const
1003 {
1004 this->A->precondition_SSOR(dst,
1005 src,
1006 this->relaxation,
1007 pos_right_of_diagonal);
1008 }
1009
1010 template <typename VectorType,
1011 std::enable_if_t<has_SSOR_step<MatrixType, VectorType>,
1012 MatrixType> * = nullptr>
1013 void
1014 step(VectorType &dst, const VectorType &src) const
1015 {
1016 this->A->SSOR_step(dst, src, this->relaxation);
1017 }
1018
1019 template <typename VectorType,
1020 std::enable_if_t<!has_SSOR_step<MatrixType, VectorType>,
1021 MatrixType> * = nullptr>
1022 void
1023 step(VectorType &, const VectorType &) const
1024 {
1025 AssertThrow(false,
1026 ExcMessage(
1027 "Matrix A does not provide a SSOR_step() function!"));
1028 }
1029
1030 template <typename VectorType>
1031 void
1032 Tstep(VectorType &dst, const VectorType &src) const
1033 {
1034 // call step, since preconditioner is symmetrical
1035 this->step(dst, src);
1036 }
1037
1038 private:
1039 const ObserverPointer<const MatrixType> A;
1040 const double relaxation;
1041
1046 std::vector<std::size_t> pos_right_of_diagonal;
1047 };
1048
1049 template <typename MatrixType>
1050 class PreconditionPSORImpl
1051 {
1052 public:
1053 using size_type = typename MatrixType::size_type;
1054
1055 PreconditionPSORImpl(const MatrixType &A,
1056 const double relaxation,
1057 const std::vector<size_type> &permutation,
1058 const std::vector<size_type> &inverse_permutation)
1059 : A(&A)
1060 , relaxation(relaxation)
1061 , permutation(permutation)
1062 , inverse_permutation(inverse_permutation)
1063 {}
1064
1065 template <typename VectorType>
1066 void
1067 vmult(VectorType &dst, const VectorType &src) const
1068 {
1069 dst = src;
1070 this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
1071 }
1072
1073 template <typename VectorType>
1074 void
1075 Tvmult(VectorType &dst, const VectorType &src) const
1076 {
1077 dst = src;
1078 this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
1079 }
1080
1081 private:
1082 const ObserverPointer<const MatrixType> A;
1083 const double relaxation;
1084
1085 const std::vector<size_type> &permutation;
1086 const std::vector<size_type> &inverse_permutation;
1087 };
1088
1089 template <typename MatrixType,
1090 typename PreconditionerType,
1091 typename VectorType,
1092 std::enable_if_t<has_step_omega<PreconditionerType, VectorType>,
1093 PreconditionerType> * = nullptr>
1094 void
1095 step(const MatrixType &,
1096 const PreconditionerType &preconditioner,
1097 VectorType &dst,
1098 const VectorType &src,
1099 const double relaxation,
1100 VectorType &,
1101 VectorType &)
1102 {
1103 preconditioner.step(dst, src, relaxation);
1104 }
1105
1106 template <
1107 typename MatrixType,
1108 typename PreconditionerType,
1109 typename VectorType,
1110 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1111 has_step<PreconditionerType, VectorType>,
1112 PreconditionerType> * = nullptr>
1113 void
1114 step(const MatrixType &,
1115 const PreconditionerType &preconditioner,
1116 VectorType &dst,
1117 const VectorType &src,
1118 const double relaxation,
1119 VectorType &,
1120 VectorType &)
1121 {
1122 Assert(relaxation == 1.0, ExcInternalError());
1123
1124 (void)relaxation;
1125
1126 preconditioner.step(dst, src);
1127 }
1128
1129 template <
1130 typename MatrixType,
1131 typename PreconditionerType,
1132 typename VectorType,
1133 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1134 !has_step<PreconditionerType, VectorType>,
1135 PreconditionerType> * = nullptr>
1136 void
1137 step(const MatrixType &A,
1138 const PreconditionerType &preconditioner,
1139 VectorType &dst,
1140 const VectorType &src,
1141 const double relaxation,
1142 VectorType &residual,
1143 VectorType &tmp)
1144 {
1145 residual.reinit(dst, true);
1146 tmp.reinit(dst, true);
1147
1148 A.vmult(residual, dst);
1149 residual.sadd(-1.0, 1.0, src);
1150
1151 preconditioner.vmult(tmp, residual);
1152 dst.add(relaxation, tmp);
1153 }
1154
1155 template <typename MatrixType,
1156 typename PreconditionerType,
1157 typename VectorType,
1158 std::enable_if_t<has_Tstep_omega<PreconditionerType, VectorType>,
1159 PreconditionerType> * = nullptr>
1160 void
1161 Tstep(const MatrixType &,
1162 const PreconditionerType &preconditioner,
1163 VectorType &dst,
1164 const VectorType &src,
1165 const double relaxation,
1166 VectorType &,
1167 VectorType &)
1168 {
1169 preconditioner.Tstep(dst, src, relaxation);
1170 }
1171
1172 template <
1173 typename MatrixType,
1174 typename PreconditionerType,
1175 typename VectorType,
1176 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1177 has_Tstep<PreconditionerType, VectorType>,
1178 PreconditionerType> * = nullptr>
1179 void
1180 Tstep(const MatrixType &,
1181 const PreconditionerType &preconditioner,
1182 VectorType &dst,
1183 const VectorType &src,
1184 const double relaxation,
1185 VectorType &,
1186 VectorType &)
1187 {
1188 Assert(relaxation == 1.0, ExcInternalError());
1189
1190 (void)relaxation;
1191
1192 preconditioner.Tstep(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 &A, VectorType &dst, const VectorType &src)
1201 {
1202 A.Tvmult(dst, src);
1203 }
1204
1205 template <typename MatrixType,
1206 typename VectorType,
1207 std::enable_if_t<!has_Tvmult<MatrixType, VectorType>, MatrixType>
1208 * = nullptr>
1209 void
1210 Tvmult(const MatrixType &, VectorType &, const VectorType &)
1211 {
1212 AssertThrow(false,
1213 ExcMessage("Matrix A does not provide a Tvmult() function!"));
1214 }
1215
1216 template <
1217 typename MatrixType,
1218 typename PreconditionerType,
1219 typename VectorType,
1220 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1221 !has_Tstep<PreconditionerType, VectorType>,
1222 PreconditionerType> * = nullptr>
1223 void
1224 Tstep(const MatrixType &A,
1225 const PreconditionerType &preconditioner,
1226 VectorType &dst,
1227 const VectorType &src,
1228 const double relaxation,
1229 VectorType &residual,
1230 VectorType &tmp)
1231 {
1232 residual.reinit(dst, true);
1233 tmp.reinit(dst, true);
1234
1235 Tvmult(A, residual, dst);
1236 residual.sadd(-1.0, 1.0, src);
1237
1238 Tvmult(preconditioner, tmp, residual);
1239 dst.add(relaxation, tmp);
1240 }
1241
1242 // 0) general implementation
1243 template <typename MatrixType,
1244 typename PreconditionerType,
1245 typename VectorType,
1246 std::enable_if_t<!has_vmult_with_std_functions_for_precondition<
1247 PreconditionerType,
1248 VectorType>,
1249 int> * = nullptr>
1250 void
1251 step_operations(const MatrixType &A,
1252 const PreconditionerType &preconditioner,
1253 VectorType &dst,
1254 const VectorType &src,
1255 const double relaxation,
1256 VectorType &tmp1,
1257 VectorType &tmp2,
1258 const unsigned int i,
1259 const bool transposed)
1260 {
1261 if (i == 0)
1262 {
1263 if (transposed)
1264 Tvmult(preconditioner, dst, src);
1265 else
1266 preconditioner.vmult(dst, src);
1267
1268 if (relaxation != 1.0)
1269 dst *= relaxation;
1270 }
1271 else
1272 {
1273 if (transposed)
1274 Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1275 else
1276 step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1277 }
1278 }
1279
1280 // 1) specialized implementation with a preconditioner that accepts
1281 // ranges
1282 template <
1283 typename MatrixType,
1284 typename PreconditionerType,
1285 typename VectorType,
1286 std::enable_if_t<
1287 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1288 VectorType> &&
1289 !has_vmult_with_std_functions_for_precondition<MatrixType,
1290 VectorType>,
1291 int> * = nullptr>
1292 void
1293 step_operations(const MatrixType &A,
1294 const PreconditionerType &preconditioner,
1295 VectorType &dst,
1296 const VectorType &src,
1297 const double relaxation,
1298 VectorType &tmp,
1299 VectorType &,
1300 const unsigned int i,
1301 const bool transposed)
1302 {
1303 (void)transposed;
1304 using Number = typename VectorType::value_type;
1305 Number *dst_ptr = dst.begin();
1306 const Number *src_ptr = src.begin();
1307
1308 if (i == 0)
1309 {
1310 preconditioner.vmult(
1311 dst,
1312 src,
1313 [&](const unsigned int start_range, const unsigned int end_range) {
1314 // zero 'dst' before running the vmult operation
1315 if (end_range > start_range)
1316 std::memset(dst_ptr + start_range,
1317 0,
1318 sizeof(Number) * (end_range - start_range));
1319 },
1320 [&](const unsigned int start_range, const unsigned int end_range) {
1321 if (relaxation == 1.0)
1322 return; // nothing to do
1323
1324 DEAL_II_OPENMP_SIMD_PRAGMA
1325 for (std::size_t i = start_range; i < end_range; ++i)
1326 dst_ptr[i] *= relaxation;
1327 });
1328 }
1329 else
1330 {
1331 tmp.reinit(src, true);
1332
1333 Assert(transposed == false, ExcNotImplemented());
1334
1335 A.vmult(tmp, dst);
1336
1337 preconditioner.vmult(
1338 dst,
1339 tmp,
1340 [&](const unsigned int start_range, const unsigned int end_range) {
1341 const auto tmp_ptr = tmp.begin();
1342
1343 if (relaxation == 1.0)
1344 {
1345 DEAL_II_OPENMP_SIMD_PRAGMA
1346 for (std::size_t i = start_range; i < end_range; ++i)
1347 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1348 }
1349 else
1350 {
1351 // note: we scale the residual here to be able to add into
1352 // the dst vector, which contains the solution from the last
1353 // iteration
1354 DEAL_II_OPENMP_SIMD_PRAGMA
1355 for (std::size_t i = start_range; i < end_range; ++i)
1356 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1357 }
1358 },
1359 [&](const unsigned int, const unsigned int) {
1360 // nothing to do, since scaling by the relaxation factor
1361 // has been done in the pre operation
1362 });
1363 }
1364 }
1365
1366 // 2) specialized implementation with a preconditioner and a matrix that
1367 // accepts ranges
1368 template <
1369 typename MatrixType,
1370 typename PreconditionerType,
1371 typename VectorType,
1372 std::enable_if_t<
1373 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1374 VectorType> &&
1375 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1376 int> * = nullptr>
1377 void
1378 step_operations(const MatrixType &A,
1379 const PreconditionerType &preconditioner,
1380 VectorType &dst,
1381 const VectorType &src,
1382 const double relaxation,
1383 VectorType &tmp,
1384 VectorType &,
1385 const unsigned int i,
1386 const bool transposed)
1387 {
1388 (void)transposed;
1389 using Number = typename VectorType::value_type;
1390
1391 Number *dst_ptr = dst.begin();
1392 const Number *src_ptr = src.begin();
1393
1394 if (i == 0)
1395 {
1396 preconditioner.vmult(
1397 dst,
1398 src,
1399 [&](const unsigned int start_range, const unsigned int end_range) {
1400 // zero 'dst' before running the vmult operation
1401 if (end_range > start_range)
1402 std::memset(dst_ptr + start_range,
1403 0,
1404 sizeof(Number) * (end_range - start_range));
1405 },
1406 [&](const unsigned int start_range, const unsigned int end_range) {
1407 if (relaxation == 1.0)
1408 return; // nothing to do
1409
1410 DEAL_II_OPENMP_SIMD_PRAGMA
1411 for (std::size_t i = start_range; i < end_range; ++i)
1412 dst_ptr[i] *= relaxation;
1413 });
1414 }
1415 else
1416 {
1417 tmp.reinit(src, true);
1418 const auto tmp_ptr = tmp.begin();
1419
1420 Assert(transposed == false, ExcNotImplemented());
1421
1422 A.vmult(
1423 tmp,
1424 dst,
1425 [&](const unsigned int start_range, const unsigned int end_range) {
1426 // zero 'tmp' before running the vmult
1427 // operation
1428 if (end_range > start_range)
1429 std::memset(tmp_ptr + start_range,
1430 0,
1431 sizeof(Number) * (end_range - start_range));
1432 },
1433 [&](const unsigned int start_range, const unsigned int end_range) {
1434 if (relaxation == 1.0)
1435 {
1436 DEAL_II_OPENMP_SIMD_PRAGMA
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
1445 DEAL_II_OPENMP_SIMD_PRAGMA
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 <
1460 typename MatrixType,
1461 typename VectorType,
1462 std::enable_if_t<
1463 !IsBlockVector<VectorType>::value &&
1464 !std::is_same_v<
1465 VectorType,
1466 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
1467 MemorySpace::Default>> &&
1468 !has_vmult_with_std_functions<MatrixType,
1469 VectorType,
1470 ::DiagonalMatrix<VectorType>>,
1471 VectorType> * = nullptr>
1472 void
1473 step_operations(const MatrixType &A,
1474 const ::DiagonalMatrix<VectorType> &preconditioner,
1475 VectorType &dst,
1476 const VectorType &src,
1477 const double relaxation,
1478 VectorType &tmp,
1479 VectorType &,
1480 const unsigned int i,
1481 const bool transposed)
1482 {
1483 using Number = typename VectorType::value_type;
1484
1485 if (i == 0)
1486 {
1487 Number *dst_ptr = dst.begin();
1488 const Number *src_ptr = src.begin();
1489 const Number *diag_ptr = preconditioner.get_vector().begin();
1490
1491 if (relaxation == 1.0)
1492 {
1493 DEAL_II_OPENMP_SIMD_PRAGMA
1494 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1495 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1496 }
1497 else
1498 {
1499 DEAL_II_OPENMP_SIMD_PRAGMA
1500 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1501 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1502 }
1503 }
1504 else
1505 {
1506 tmp.reinit(src, true);
1507
1508 if (transposed)
1509 Tvmult(A, tmp, dst);
1510 else
1511 A.vmult(tmp, dst);
1512
1513 Number *dst_ptr = dst.begin();
1514 const Number *src_ptr = src.begin();
1515 const Number *tmp_ptr = tmp.begin();
1516 const Number *diag_ptr = preconditioner.get_vector().begin();
1517
1518 if (relaxation == 1.0)
1519 {
1520 DEAL_II_OPENMP_SIMD_PRAGMA
1521 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1522 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1523 }
1524 else
1525 {
1526 DEAL_II_OPENMP_SIMD_PRAGMA
1527 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1528 dst_ptr[i] +=
1529 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1530 }
1531 }
1532 }
1533
1534 // 4) specialized implementation for inverse-diagonal preconditioner and
1535 // matrix that accepts ranges
1536 template <typename MatrixType,
1537 typename VectorType,
1538 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1539 has_vmult_with_std_functions<
1540 MatrixType,
1541 VectorType,
1542 ::DiagonalMatrix<VectorType>>,
1543 VectorType> * = nullptr>
1544 void
1545 step_operations(const MatrixType &A,
1546 const ::DiagonalMatrix<VectorType> &preconditioner,
1547 VectorType &dst,
1548 const VectorType &src,
1549 const double relaxation,
1550 VectorType &tmp,
1551 VectorType &,
1552 const unsigned int i,
1553 const bool transposed)
1554 {
1555 (void)transposed;
1556 using Number = typename VectorType::value_type;
1557
1558 if (i == 0)
1559 {
1560 Number *dst_ptr = dst.begin();
1561 const Number *src_ptr = src.begin();
1562 const Number *diag_ptr = preconditioner.get_vector().begin();
1563
1564 if (relaxation == 1.0)
1565 {
1566 DEAL_II_OPENMP_SIMD_PRAGMA
1567 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1568 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1569 }
1570 else
1571 {
1572 DEAL_II_OPENMP_SIMD_PRAGMA
1573 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1574 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1575 }
1576 }
1577 else
1578 {
1579 tmp.reinit(src, true);
1580
1581 Assert(transposed == false, ExcNotImplemented());
1582
1583 A.vmult(
1584 tmp,
1585 dst,
1586 [&](const unsigned int start_range, const unsigned int end_range) {
1587 // zero 'tmp' before running the vmult operation
1588 if (end_range > start_range)
1589 std::memset(tmp.begin() + start_range,
1590 0,
1591 sizeof(Number) * (end_range - start_range));
1592 },
1593 [&](const unsigned int begin, const unsigned int end) {
1594 const Number *dst_ptr = dst.begin();
1595 const Number *src_ptr = src.begin();
1596 Number *tmp_ptr = tmp.begin();
1597 const Number *diag_ptr = preconditioner.get_vector().begin();
1598
1599 // for efficiency reason, write back to temp_vector that is
1600 // already read (avoid read-for-ownership)
1601 if (relaxation == 1.0)
1602 {
1603 DEAL_II_OPENMP_SIMD_PRAGMA
1604 for (std::size_t i = begin; i < end; ++i)
1605 tmp_ptr[i] =
1606 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1607 }
1608 else
1609 {
1610 DEAL_II_OPENMP_SIMD_PRAGMA
1611 for (std::size_t i = begin; i < end; ++i)
1612 tmp_ptr[i] = dst_ptr[i] + relaxation *
1613 (src_ptr[i] - tmp_ptr[i]) *
1614 diag_ptr[i];
1615 }
1616 });
1617
1618 tmp.swap(dst);
1619 }
1620 }
1621
1622 } // namespace PreconditionRelaxation
1623} // namespace internal
1624
1625#endif
1626
1627
1628
1655template <typename MatrixType = SparseMatrix<double>>
1656class PreconditionJacobi
1657 : public PreconditionRelaxation<
1658 MatrixType,
1659 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1660{
1661 using PreconditionerType =
1662 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1663 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1664
1665public:
1669 using AdditionalData = typename BaseClass::AdditionalData;
1670
1674 void
1675 initialize(const MatrixType &A,
1676 const AdditionalData &parameters = AdditionalData());
1677};
1678
1679
1725template <typename MatrixType = SparseMatrix<double>>
1726class PreconditionSOR
1727 : public PreconditionRelaxation<
1728 MatrixType,
1729 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1730{
1731 using PreconditionerType =
1732 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1733 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1734
1735public:
1739 using AdditionalData = typename BaseClass::AdditionalData;
1740
1744 void
1745 initialize(const MatrixType &A,
1746 const AdditionalData &parameters = AdditionalData());
1747};
1748
1749
1750
1777template <typename MatrixType = SparseMatrix<double>>
1778class PreconditionSSOR
1779 : public PreconditionRelaxation<
1780 MatrixType,
1781 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1782{
1783 using PreconditionerType =
1784 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1785 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1786
1787public:
1791 using AdditionalData = typename BaseClass::AdditionalData;
1792
1798 void
1799 initialize(const MatrixType &A,
1800 const AdditionalData &parameters = AdditionalData());
1801};
1802
1803
1833template <typename MatrixType = SparseMatrix<double>>
1834class PreconditionPSOR
1835 : public PreconditionRelaxation<
1836 MatrixType,
1837 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1838{
1839 using PreconditionerType =
1840 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1841 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1842
1843public:
1847 using size_type = typename BaseClass::size_type;
1848
1852 class AdditionalData
1853 {
1854 public:
1865 AdditionalData(const std::vector<size_type> &permutation,
1866 const std::vector<size_type> &inverse_permutation,
1867 const typename BaseClass::AdditionalData &parameters =
1868 typename BaseClass::AdditionalData());
1869
1873 const std::vector<size_type> &permutation;
1877 const std::vector<size_type> &inverse_permutation;
1881 typename BaseClass::AdditionalData parameters;
1882 };
1883
1895 void
1896 initialize(const MatrixType &A,
1897 const std::vector<size_type> &permutation,
1898 const std::vector<size_type> &inverse_permutation,
1899 const typename BaseClass::AdditionalData &parameters =
1900 typename BaseClass::AdditionalData());
1901
1912 void
1913 initialize(const MatrixType &A, const AdditionalData &additional_data);
1914};
1915
1916
1917
2106template <typename MatrixType = SparseMatrix<double>,
2107 typename VectorType = Vector<double>,
2108 typename PreconditionerType = DiagonalMatrix<VectorType>>
2109class PreconditionChebyshev : public EnableObserverPointer
2110{
2111public:
2115 using size_type = types::global_dof_index;
2116
2121 struct AdditionalData
2122 : public internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>
2123 {
2124 using EigenvalueAlgorithm = internal::EigenvalueAlgorithm;
2125
2129 enum class PolynomialType
2130 {
2134 first_kind,
2139 fourth_kind
2140 };
2141
2145 AdditionalData(
2146 const unsigned int degree = 1,
2147 const double smoothing_range = 0.,
2148 const unsigned int eig_cg_n_iterations = 8,
2149 const double eig_cg_residual = 1e-2,
2150 const double max_eigenvalue = 1,
2151 const EigenvalueAlgorithm eigenvalue_algorithm =
2152 EigenvalueAlgorithm::lanczos,
2153 const PolynomialType polynomial_type = PolynomialType::first_kind,
2154 const double safety_factor = 1.2);
2155
2168 unsigned int degree;
2169
2173 PolynomialType polynomial_type;
2174 };
2175
2176
2180 PreconditionChebyshev();
2181
2193 void
2194 initialize(const MatrixType &matrix,
2195 const AdditionalData &additional_data = AdditionalData());
2196
2201 void
2202 vmult(VectorType &dst, const VectorType &src) const;
2203
2211 double
2212 vmult_with_last_residual_norm(VectorType &dst, const VectorType &src) const;
2213
2218 void
2219 Tvmult(VectorType &dst, const VectorType &src) const;
2220
2224 void
2225 step(VectorType &dst, const VectorType &src) const;
2226
2234 double
2235 step_with_last_residual_norm(VectorType &dst, const VectorType &src) const;
2236
2240 void
2241 Tstep(VectorType &dst, const VectorType &src) const;
2242
2246 void
2247 clear();
2248
2253 size_type
2254 m() const;
2255
2260 size_type
2261 n() const;
2262
2263 using EigenvalueInformation = internal::EigenvalueInformation;
2264
2272 void
2273 set_degree(const unsigned int degree);
2274
2287 EigenvalueInformation
2288 estimate_eigenvalues(const VectorType &src) const;
2289
2290private:
2295 template <bool use_transpose>
2296 double
2297 apply_internal(const bool zero_out_dst,
2298 const VectorType &src,
2299 VectorType &dst,
2300 const bool compute_final_norm = false) const;
2301
2305 ObserverPointer<
2306 const MatrixType,
2307 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>>
2308 matrix_ptr;
2309
2313 mutable VectorType solution_old;
2314
2318 mutable VectorType temp_vector1;
2319
2323 mutable VectorType temp_vector2;
2324
2329 AdditionalData data;
2330
2334 double theta;
2335
2340 double delta;
2341
2346 bool eigenvalues_are_initialized;
2347
2352 mutable Threads::Mutex mutex;
2353};
2354
2355
2356
2358/* ---------------------------------- Inline functions ------------------- */
2359
2360#ifndef DOXYGEN
2361
2362
2363namespace internal
2364{
2365 template <typename VectorType>
2366 void
2367 set_initial_guess(VectorType &vector)
2368 {
2369 vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2370 if (vector.locally_owned_elements().is_element(0))
2371 vector(0) = 0.;
2372 }
2373
2374
2375
2376 template <typename Number>
2377 void
2378 set_initial_guess(::Vector<Number> &vector)
2379 {
2380 // Choose a high-frequency mode consisting of numbers between 0 and 10
2381 // that is cheap to compute (cheaper than random numbers) but avoids
2382 // obviously re-occurring numbers in multi-component systems by choosing
2383 // a period of 11
2384 for (unsigned int i = 0; i < vector.size(); ++i)
2385 vector(i) = i % 11;
2386
2387 // Then make sure the vector has mean value zero:
2388 const Number mean_value = vector.mean_value();
2389 vector.add(-mean_value);
2390 }
2391
2392
2393
2394 template <typename Number, typename MemorySpace>
2395 void
2396 set_initial_guess(
2397 ::LinearAlgebra::distributed::BlockVector<Number, MemorySpace>
2398 &vector)
2399 {
2400 for (unsigned int block = 0; block < vector.n_blocks(); ++block)
2401 set_initial_guess(vector.block(block));
2402 }
2403
2404
2405
2406 template <typename Number, typename MemorySpace>
2407 void
2408 set_initial_guess(
2409 ::LinearAlgebra::distributed::Vector<Number, MemorySpace> &vector)
2410 {
2411 // Choose a high-frequency mode consisting of numbers between 0 and 10
2412 // that is cheap to compute (cheaper than random numbers) but avoids
2413 // obviously re-occurring numbers in multi-component systems by choosing
2414 // a period of 11.
2415 // Make initial guess robust with respect to number of processors
2416 // by operating on the global index.
2417 types::global_dof_index first_local_range = 0;
2418 if (!vector.locally_owned_elements().is_empty())
2419 first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2420
2421 const auto n_local_elements = vector.locally_owned_size();
2422 Number *values_ptr = vector.get_values();
2423 typename MemorySpace::kokkos_space::execution_space exec;
2424 Kokkos::RangePolicy<typename MemorySpace::kokkos_space::execution_space,
2425 Kokkos::IndexType<types::global_dof_index>>
2426 policy(exec, 0, n_local_elements);
2427
2428 Kokkos::parallel_for(
2429 "::PreconditionChebyshev::set_initial_guess",
2430 policy,
2431 KOKKOS_LAMBDA(types::global_dof_index i) {
2432 values_ptr[i] = (i + first_local_range) % 11;
2433 });
2434 exec.fence();
2435
2436 // Then make sure the vector has mean value zero:
2437 const Number mean_value = vector.mean_value();
2438 vector.add(-mean_value);
2439 }
2440
2441
2442
2443 struct EigenvalueTracker
2444 {
2445 public:
2446 void
2447 slot(const std::vector<double> &eigenvalues)
2448 {
2449 values = eigenvalues;
2450 }
2451
2452 std::vector<double> values;
2453 };
2454
2455
2456
2457 template <typename MatrixType,
2458 typename VectorType,
2459 typename PreconditionerType>
2460 double
2461 power_iteration(const MatrixType &matrix,
2462 VectorType &eigenvector,
2463 const PreconditionerType &preconditioner,
2464 const unsigned int n_iterations)
2465 {
2466 typename VectorType::value_type eigenvalue_estimate = 0.;
2467 eigenvector /= eigenvector.l2_norm();
2468 VectorType vector1, vector2;
2469 vector1.reinit(eigenvector, true);
2470 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2471 vector2.reinit(eigenvector, true);
2472 for (unsigned int i = 0; i < n_iterations; ++i)
2473 {
2474 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2475 {
2476 matrix.vmult(vector2, eigenvector);
2477 preconditioner.vmult(vector1, vector2);
2478 }
2479 else
2480 matrix.vmult(vector1, eigenvector);
2481
2482 eigenvalue_estimate = eigenvector * vector1;
2483
2484 vector1 /= vector1.l2_norm();
2485 eigenvector.swap(vector1);
2486 }
2487 return std::abs(eigenvalue_estimate);
2488 }
2489
2490
2491
2492 template <typename MatrixType,
2493 typename VectorType,
2494 typename PreconditionerType>
2495 EigenvalueInformation
2496 estimate_eigenvalues(
2497 const EigenvalueAlgorithmAdditionalData<PreconditionerType> &data,
2498 const MatrixType *matrix_ptr,
2499 VectorType &solution_old,
2500 VectorType &temp_vector1,
2501 const unsigned int degree,
2502 const double safety_factor)
2503 {
2504 Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2505
2506 Assert(
2507 safety_factor >= 1.,
2508 ExcMessage(
2509 "The safety factor must be at least 1.0 to ensure the maximum eigenvalue is "
2510 "not underestimated. Please set safety_factor >= 1.0 in your AdditionalData."));
2511
2512 EigenvalueInformation info{};
2513
2514 if (data.eig_cg_n_iterations > 0)
2515 {
2516 Assert(data.eig_cg_n_iterations > 2,
2517 ExcMessage(
2518 "Need to set at least two iterations to find eigenvalues."));
2519
2520 internal::EigenvalueTracker eigenvalue_tracker;
2521
2522 // set an initial guess that contains some high-frequency parts (to the
2523 // extent possible without knowing the discretization and the numbering)
2524 // to trigger high eigenvalues according to the external function
2525 internal::set_initial_guess(temp_vector1);
2526 data.constraints.set_zero(temp_vector1);
2527
2528 // We want to run the power iteration starting with temp_vector1, but
2529 // for that it needs to be a nonzero vector. However, there are
2530 // situations where it only contains zeroes (e.g., if all dofs are
2531 // constrained). If that is the case, there is nothing we can do
2532 // because in both of the two methods below, we either solve a linear
2533 // system with temp_vector1 as the right hand side, or use it as a
2534 // starting guess for a power iteration, but in neither case a
2535 // zero vector yields anything useful. So skip the following code and
2536 // instead use the fallback below if temp_vector1 is the zero vector.
2537 if (temp_vector1.all_zero() == false)
2538 {
2539 if (data.eigenvalue_algorithm ==
2540 internal::EigenvalueAlgorithm::lanczos)
2541 {
2542 // set a very strict tolerance to force at least two iterations
2543 IterationNumberControl control(data.eig_cg_n_iterations,
2544 1e-10,
2545 false,
2546 false);
2547
2548 ::SolverCG<VectorType> solver(control);
2549 solver.connect_eigenvalues_slot(
2550 [&eigenvalue_tracker](
2551 const std::vector<double> &eigenvalues) {
2552 eigenvalue_tracker.slot(eigenvalues);
2553 });
2554
2555 solver.solve(*matrix_ptr,
2556 solution_old,
2557 temp_vector1,
2558 *data.preconditioner);
2559
2560 info.cg_iterations = control.last_step();
2561 }
2562 else if (data.eigenvalue_algorithm ==
2563 internal::EigenvalueAlgorithm::power_iteration)
2564 {
2565 (void)degree;
2566
2567 Assert(degree != numbers::invalid_unsigned_int,
2568 ExcMessage(
2569 "Cannot estimate the minimal eigenvalue with the "
2570 "power iteration"));
2571
2572 eigenvalue_tracker.values.push_back(
2573 internal::power_iteration(*matrix_ptr,
2574 temp_vector1,
2575 *data.preconditioner,
2576 data.eig_cg_n_iterations));
2577 }
2578 else
2579 DEAL_II_NOT_IMPLEMENTED();
2580 }
2581
2582 // read the eigenvalues from the attached eigenvalue tracker
2583 if (eigenvalue_tracker.values.empty())
2584 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2585 else
2586 {
2587 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2588
2589 // include a safety factor since the CG method will in general not
2590 // be converged
2591 info.max_eigenvalue_estimate =
2592 safety_factor * eigenvalue_tracker.values.back();
2593 }
2594 }
2595 else
2596 {
2597 info.max_eigenvalue_estimate = data.max_eigenvalue;
2598 info.min_eigenvalue_estimate =
2599 data.max_eigenvalue / data.smoothing_range;
2600 }
2601
2602 return info;
2603 }
2604} // namespace internal
2605
2606
2607
2608inline PreconditionIdentity::PreconditionIdentity()
2609 : n_rows(0)
2610 , n_columns(0)
2611{}
2612
2613template <typename MatrixType>
2614inline void
2615PreconditionIdentity::initialize(const MatrixType &matrix,
2616 const PreconditionIdentity::AdditionalData &)
2617{
2618 n_rows = matrix.m();
2619 n_columns = matrix.n();
2620}
2621
2622
2623template <typename VectorType>
2624inline void
2625PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
2626{
2627 dst = src;
2628}
2629
2630
2631
2632template <typename VectorType>
2633inline void
2634PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
2635{
2636 dst = src;
2637}
2638
2639template <typename VectorType>
2640inline void
2641PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
2642{
2643 dst += src;
2644}
2645
2646
2647
2648template <typename VectorType>
2649inline void
2650PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
2651{
2652 dst += src;
2653}
2654
2655
2656
2657inline void
2658PreconditionIdentity::clear()
2659{}
2660
2661
2662
2663inline PreconditionIdentity::size_type
2664PreconditionIdentity::m() const
2665{
2666 Assert(n_rows != 0, ExcNotInitialized());
2667 return n_rows;
2668}
2669
2670inline PreconditionIdentity::size_type
2671PreconditionIdentity::n() const
2672{
2673 Assert(n_columns != 0, ExcNotInitialized());
2674 return n_columns;
2675}
2676
2677//---------------------------------------------------------------------------
2678
2679inline PreconditionRichardson::AdditionalData::AdditionalData(
2680 const double relaxation)
2681 : relaxation(relaxation)
2682{}
2683
2684
2685inline PreconditionRichardson::PreconditionRichardson()
2686 : relaxation(0)
2687 , n_rows(0)
2688 , n_columns(0)
2689{
2690 AdditionalData add_data;
2691 relaxation = add_data.relaxation;
2692}
2693
2694
2695
2696inline void
2697PreconditionRichardson::initialize(
2698 const PreconditionRichardson::AdditionalData &parameters)
2699{
2700 relaxation = parameters.relaxation;
2701}
2702
2703
2704
2705template <typename MatrixType>
2706inline void
2707PreconditionRichardson::initialize(
2708 const MatrixType &matrix,
2709 const PreconditionRichardson::AdditionalData &parameters)
2710{
2711 relaxation = parameters.relaxation;
2712 n_rows = matrix.m();
2713 n_columns = matrix.n();
2714}
2715
2716
2717
2718template <typename VectorType>
2719inline void
2720PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2721{
2722 static_assert(
2723 std::is_same_v<size_type, typename VectorType::size_type>,
2724 "PreconditionRichardson and VectorType must have the same size_type.");
2725
2726 dst.equ(relaxation, src);
2727}
2728
2729
2730
2731template <typename VectorType>
2732inline void
2733PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2734{
2735 static_assert(
2736 std::is_same_v<size_type, typename VectorType::size_type>,
2737 "PreconditionRichardson and VectorType must have the same size_type.");
2738
2739 dst.equ(relaxation, src);
2740}
2741
2742template <typename VectorType>
2743inline void
2744PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2745{
2746 static_assert(
2747 std::is_same_v<size_type, typename VectorType::size_type>,
2748 "PreconditionRichardson and VectorType must have the same size_type.");
2749
2750 dst.add(relaxation, src);
2751}
2752
2753
2754
2755template <typename VectorType>
2756inline void
2757PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2758{
2759 static_assert(
2760 std::is_same_v<size_type, typename VectorType::size_type>,
2761 "PreconditionRichardson and VectorType must have the same size_type.");
2762
2763 dst.add(relaxation, src);
2764}
2765
2766inline PreconditionRichardson::size_type
2767PreconditionRichardson::m() const
2768{
2769 Assert(n_rows != 0, ExcNotInitialized());
2770 return n_rows;
2771}
2772
2773inline PreconditionRichardson::size_type
2774PreconditionRichardson::n() const
2775{
2776 Assert(n_columns != 0, ExcNotInitialized());
2777 return n_columns;
2778}
2779
2780//---------------------------------------------------------------------------
2781
2782template <typename MatrixType, typename PreconditionerType>
2783inline void
2784PreconditionRelaxation<MatrixType, PreconditionerType>::initialize(
2785 const MatrixType &rA,
2786 const AdditionalData &parameters)
2787{
2788 A = &rA;
2789 eigenvalues_are_initialized = false;
2790
2791 Assert(parameters.preconditioner, ExcNotInitialized());
2792
2793 this->data = parameters;
2794}
2795
2796
2797template <typename MatrixType, typename PreconditionerType>
2798inline void
2799PreconditionRelaxation<MatrixType, PreconditionerType>::clear()
2800{
2801 eigenvalues_are_initialized = false;
2802 A = nullptr;
2803 data.relaxation = 1.0;
2804 data.preconditioner = nullptr;
2805}
2806
2807template <typename MatrixType, typename PreconditionerType>
2808inline
2809 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2810 PreconditionRelaxation<MatrixType, PreconditionerType>::m() const
2811{
2812 Assert(A != nullptr, ExcNotInitialized());
2813 return A->m();
2814}
2815
2816template <typename MatrixType, typename PreconditionerType>
2817inline
2818 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2819 PreconditionRelaxation<MatrixType, PreconditionerType>::n() const
2820{
2821 Assert(A != nullptr, ExcNotInitialized());
2822 return A->n();
2823}
2824
2825template <typename MatrixType, typename PreconditionerType>
2826template <typename VectorType>
2827inline void
2828PreconditionRelaxation<MatrixType, PreconditionerType>::vmult(
2829 VectorType &dst,
2830 const VectorType &src) const
2831{
2832 Assert(this->A != nullptr, ExcNotInitialized());
2833 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2834
2835 if (eigenvalues_are_initialized == false)
2836 estimate_eigenvalues(src);
2837
2838 VectorType tmp1, tmp2;
2839
2840 for (unsigned int i = 0; i < data.n_iterations; ++i)
2841 internal::PreconditionRelaxation::step_operations(*A,
2842 *data.preconditioner,
2843 dst,
2844 src,
2845 data.relaxation,
2846 tmp1,
2847 tmp2,
2848 i,
2849 false);
2850}
2851
2852template <typename MatrixType, typename PreconditionerType>
2853template <typename VectorType>
2854inline void
2855PreconditionRelaxation<MatrixType, PreconditionerType>::Tvmult(
2856 VectorType &dst,
2857 const VectorType &src) const
2858{
2859 Assert(this->A != nullptr, ExcNotInitialized());
2860 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2861
2862 if (eigenvalues_are_initialized == false)
2863 estimate_eigenvalues(src);
2864
2865 VectorType tmp1, tmp2;
2866
2867 for (unsigned int i = 0; i < data.n_iterations; ++i)
2868 internal::PreconditionRelaxation::step_operations(
2869 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2870}
2871
2872template <typename MatrixType, typename PreconditionerType>
2873template <typename VectorType>
2874inline void
2875PreconditionRelaxation<MatrixType, PreconditionerType>::step(
2876 VectorType &dst,
2877 const VectorType &src) const
2878{
2879 Assert(this->A != nullptr, ExcNotInitialized());
2880 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2881
2882 if (eigenvalues_are_initialized == false)
2883 estimate_eigenvalues(src);
2884
2885 VectorType tmp1, tmp2;
2886
2887 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2888 internal::PreconditionRelaxation::step_operations(*A,
2889 *data.preconditioner,
2890 dst,
2891 src,
2892 data.relaxation,
2893 tmp1,
2894 tmp2,
2895 i,
2896 false);
2897}
2898
2899template <typename MatrixType, typename PreconditionerType>
2900template <typename VectorType>
2901inline void
2902PreconditionRelaxation<MatrixType, PreconditionerType>::Tstep(
2903 VectorType &dst,
2904 const VectorType &src) const
2905{
2906 Assert(this->A != nullptr, ExcNotInitialized());
2907 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2908
2909 if (eigenvalues_are_initialized == false)
2910 estimate_eigenvalues(src);
2911
2912 VectorType tmp1, tmp2;
2913
2914 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2915 internal::PreconditionRelaxation::step_operations(
2916 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2917}
2918
2919template <typename MatrixType, typename PreconditionerType>
2920template <typename VectorType>
2921inline internal::EigenvalueInformation
2922PreconditionRelaxation<MatrixType, PreconditionerType>::estimate_eigenvalues(
2923 const VectorType &src) const
2924{
2925 Assert(eigenvalues_are_initialized == false, ExcInternalError());
2926
2927 EigenvalueInformation info;
2928
2929 if (data.relaxation == 0.0)
2930 {
2931 VectorType solution_old, temp_vector1;
2932
2933 solution_old.reinit(src);
2934 temp_vector1.reinit(src, true);
2935
2936 info = internal::estimate_eigenvalues<MatrixType>(data,
2937 A,
2938 solution_old,
2939 temp_vector1,
2940 data.n_iterations,
2941 data.safety_factor);
2942
2943 const double alpha =
2944 (data.smoothing_range > 1. ?
2945 info.max_eigenvalue_estimate / data.smoothing_range :
2946 std::min(0.9 * info.max_eigenvalue_estimate,
2947 info.min_eigenvalue_estimate));
2948
2949 const_cast<PreconditionRelaxation<MatrixType, PreconditionerType> *>(this)
2950 ->data.relaxation = 2.0 / (alpha + info.max_eigenvalue_estimate);
2951 }
2952
2953 const_cast<PreconditionRelaxation<MatrixType, PreconditionerType> *>(this)
2954 ->eigenvalues_are_initialized = true;
2955
2956 return info;
2957}
2958
2959template <typename MatrixType, typename PreconditionerType>
2960double
2961PreconditionRelaxation<MatrixType, PreconditionerType>::get_relaxation() const
2962{
2963 return data.relaxation;
2964}
2965
2966
2967//---------------------------------------------------------------------------
2968
2969template <typename MatrixType>
2970inline void
2971PreconditionJacobi<MatrixType>::initialize(const MatrixType &A,
2972 const AdditionalData &parameters_in)
2973{
2974 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2975 Assert(
2976 parameters_in.relaxation != 0.0,
2977 ExcMessage(
2978 "Relaxation cannot automatically be determined by PreconditionJacobi."));
2979
2980 AdditionalData parameters;
2981 parameters.relaxation = 1.0;
2982 parameters.n_iterations = parameters_in.n_iterations;
2983 parameters.preconditioner =
2984 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2985
2986 this->BaseClass::initialize(A, parameters);
2987}
2988
2989//---------------------------------------------------------------------------
2990
2991template <typename MatrixType>
2992inline void
2993PreconditionSOR<MatrixType>::initialize(const MatrixType &A,
2994 const AdditionalData &parameters_in)
2995{
2996 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2997 Assert(
2998 parameters_in.relaxation != 0.0,
2999 ExcMessage(
3000 "Relaxation cannot automatically be determined by PreconditionSOR."));
3001
3002 AdditionalData parameters;
3003 parameters.relaxation = 1.0;
3004 parameters.n_iterations = parameters_in.n_iterations;
3005 parameters.preconditioner =
3006 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
3007
3008 this->BaseClass::initialize(A, parameters);
3009}
3010
3011//---------------------------------------------------------------------------
3012
3013template <typename MatrixType>
3014inline void
3015PreconditionSSOR<MatrixType>::initialize(const MatrixType &A,
3016 const AdditionalData &parameters_in)
3017{
3018 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
3019 Assert(
3020 parameters_in.relaxation != 0.0,
3021 ExcMessage(
3022 "Relaxation cannot automatically be determined by PreconditionSSOR."));
3023
3024 AdditionalData parameters;
3025 parameters.relaxation = 1.0;
3026 parameters.n_iterations = parameters_in.n_iterations;
3027 parameters.preconditioner =
3028 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
3029
3030 this->BaseClass::initialize(A, parameters);
3031}
3032
3033
3034
3035//---------------------------------------------------------------------------
3036
3037template <typename MatrixType>
3038inline void
3039PreconditionPSOR<MatrixType>::initialize(
3040 const MatrixType &A,
3041 const std::vector<size_type> &p,
3042 const std::vector<size_type> &ip,
3043 const typename BaseClass::AdditionalData &parameters_in)
3044{
3045 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
3046 Assert(
3047 parameters_in.relaxation != 0.0,
3048 ExcMessage(
3049 "Relaxation cannot automatically be determined by PreconditionPSOR."));
3050
3051 typename BaseClass::AdditionalData parameters;
3052 parameters.relaxation = 1.0;
3053 parameters.n_iterations = parameters_in.n_iterations;
3054 parameters.preconditioner =
3055 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
3056
3057 this->BaseClass::initialize(A, parameters);
3058}
3059
3060
3061template <typename MatrixType>
3062inline void
3063PreconditionPSOR<MatrixType>::initialize(const MatrixType &A,
3064 const AdditionalData &additional_data)
3065{
3066 initialize(A,
3067 additional_data.permutation,
3068 additional_data.inverse_permutation,
3069 additional_data.parameters);
3070}
3071
3072template <typename MatrixType>
3073PreconditionPSOR<MatrixType>::AdditionalData::AdditionalData(
3074 const std::vector<size_type> &permutation,
3075 const std::vector<size_type> &inverse_permutation,
3076 const typename PreconditionPSOR<MatrixType>::BaseClass::AdditionalData
3077 &parameters)
3078 : permutation(permutation)
3079 , inverse_permutation(inverse_permutation)
3080 , parameters(parameters)
3081{}
3082
3083
3084//---------------------------------------------------------------------------
3085
3086
3087template <typename MatrixType, typename VectorType>
3088PreconditionUseMatrix<MatrixType, VectorType>::PreconditionUseMatrix(
3089 const MatrixType &M,
3090 const function_ptr method)
3091 : matrix(M)
3092 , precondition(method)
3093{}
3094
3095
3096
3097template <typename MatrixType, typename VectorType>
3098void
3099PreconditionUseMatrix<MatrixType, VectorType>::vmult(
3100 VectorType &dst,
3101 const VectorType &src) const
3102{
3103 (matrix.*precondition)(dst, src);
3104}
3105
3106//---------------------------------------------------------------------------
3107
3108namespace internal
3109{
3110
3111 template <typename PreconditionerType>
3112 inline EigenvalueAlgorithmAdditionalData<PreconditionerType>::
3113 EigenvalueAlgorithmAdditionalData(
3114 const double smoothing_range,
3115 const unsigned int eig_cg_n_iterations,
3116 const double eig_cg_residual,
3117 const double max_eigenvalue,
3118 const EigenvalueAlgorithm eigenvalue_algorithm,
3119 const double safety_factor)
3120 : smoothing_range(smoothing_range)
3121 , eig_cg_n_iterations(eig_cg_n_iterations)
3122 , eig_cg_residual(eig_cg_residual)
3123 , max_eigenvalue(max_eigenvalue)
3124 , eigenvalue_algorithm(eigenvalue_algorithm)
3125 , safety_factor(safety_factor)
3126 {}
3127
3128
3129
3130 template <typename PreconditionerType>
3131 inline EigenvalueAlgorithmAdditionalData<PreconditionerType> &
3132 EigenvalueAlgorithmAdditionalData<PreconditionerType>::operator=(
3133 const EigenvalueAlgorithmAdditionalData &other_data)
3134 {
3135 smoothing_range = other_data.smoothing_range;
3136 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3137 eig_cg_residual = other_data.eig_cg_residual;
3138 max_eigenvalue = other_data.max_eigenvalue;
3139 preconditioner = other_data.preconditioner;
3140 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3141 safety_factor = other_data.safety_factor;
3142 constraints.copy_from(other_data.constraints);
3143
3144 return *this;
3145 }
3146} // namespace internal
3147
3148template <typename MatrixType, typename PreconditionerType>
3149inline PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData::
3150 AdditionalData(const double relaxation,
3151 const unsigned int n_iterations,
3152 const double smoothing_range,
3153 const unsigned int eig_cg_n_iterations,
3154 const double eig_cg_residual,
3155 const double max_eigenvalue,
3156 const EigenvalueAlgorithm eigenvalue_algorithm,
3157 const double safety_factor)
3158 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3159 smoothing_range,
3160 eig_cg_n_iterations,
3161 eig_cg_residual,
3162 max_eigenvalue,
3163 eigenvalue_algorithm,
3164 safety_factor)
3165 , relaxation(relaxation)
3166 , n_iterations(n_iterations)
3167{}
3168
3169
3170
3171//---------------------------------------------------------------------------
3172
3173namespace internal
3174{
3175 namespace PreconditionChebyshevImplementation
3176 {
3177 // for deal.II vectors, perform updates for Chebyshev preconditioner all
3178 // at once to reduce memory transfer. Here, we select between general
3179 // vectors and deal.II vectors where we expand the loop over the (local)
3180 // size of the vector
3181
3182 // generic part for non-deal.II vectors
3183 template <typename VectorType, typename PreconditionerType>
3184 inline double
3185 vector_updates(const VectorType &rhs,
3186 const PreconditionerType &preconditioner,
3187 const unsigned int iteration_index,
3188 const double factor1,
3189 const double factor2,
3190 VectorType &solution_old,
3191 VectorType &temp_vector1,
3192 VectorType &temp_vector2,
3193 VectorType &solution,
3194 const bool compute_residual_norm)
3195 {
3196 double residual_norm = 0;
3197 if (iteration_index == 0)
3198 {
3199 if (compute_residual_norm)
3200 residual_norm = rhs.l2_norm();
3201 solution.equ(factor2, rhs);
3202 preconditioner.vmult(solution_old, solution);
3203 }
3204 else if (iteration_index == 1)
3205 {
3206 // compute t = P^{-1} * (A*x^{n} - b)
3207 if (compute_residual_norm)
3208 residual_norm =
3209 std::sqrt(temp_vector1.add_and_dot(-1.0, rhs, temp_vector1));
3210 else
3211 temp_vector1.add(-1.0, rhs);
3212 preconditioner.vmult(solution_old, temp_vector1);
3213
3214 // compute x^{n+1} = x^{n} + f_1 * x^{n} - f_2 * t
3215 solution_old.sadd(-factor2, 1 + factor1, solution);
3216 }
3217 else
3218 {
3219 // compute t = P^{-1} * (A*x^{n} - b)
3220 if (compute_residual_norm)
3221 residual_norm =
3222 std::sqrt(temp_vector1.add_and_dot(-1.0, rhs, temp_vector1));
3223 else
3224 temp_vector1.add(-1.0, rhs);
3225
3226 preconditioner.vmult(temp_vector2, temp_vector1);
3227
3228 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) - f_2 * t
3229 solution_old.sadd(-factor1, -factor2, temp_vector2);
3230 solution_old.add(1 + factor1, solution);
3231 }
3232
3233 solution.swap(solution_old);
3234 return residual_norm;
3235 }
3236
3237 // generic part for deal.II vectors
3238 template <
3239 typename Number,
3240 typename PreconditionerType,
3241 std::enable_if_t<
3242 !has_vmult_with_std_functions_for_precondition<
3243 PreconditionerType,
3244 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>>,
3245 int> * = nullptr>
3246 inline double
3247 vector_updates(
3248 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3249 const PreconditionerType &preconditioner,
3250 const unsigned int iteration_index,
3251 const double factor1_,
3252 const double factor2_,
3253 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3254 &solution_old,
3255 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3256 &temp_vector1,
3257 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3258 &temp_vector2,
3259 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution,
3260 const bool compute_residual_norm)
3261 {
3262 const Number factor1 = factor1_;
3263 const Number factor1_plus_1 = 1. + factor1_;
3264 const Number factor2 = factor2_;
3265 double residual_norm = 0;
3266
3267 if (iteration_index == 0)
3268 {
3269 if (compute_residual_norm)
3270 residual_norm = rhs.l2_norm();
3271
3272 // compute t = P^{-1} * (b)
3273 preconditioner.vmult(solution_old, rhs);
3274
3275 // compute x^{n+1} = f_2 * t
3276 const auto solution_old_ptr = solution_old.begin();
3277 DEAL_II_OPENMP_SIMD_PRAGMA
3278 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3279 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
3280 }
3281 else if (iteration_index == 1)
3282 {
3283 // compute t = P^{-1} * (A*x^{n} - b)
3284 if (compute_residual_norm)
3285 residual_norm =
3286 std::sqrt(temp_vector1.add_and_dot(-1.0, rhs, temp_vector1));
3287 else
3288 temp_vector1.add(-1.0, rhs);
3289
3290 preconditioner.vmult(solution_old, temp_vector1);
3291
3292 // compute x^{n+1} = x^{n} + f_1 * x^{n} - f_2 * t
3293 const auto solution_ptr = solution.begin();
3294 const auto solution_old_ptr = solution_old.begin();
3295 DEAL_II_OPENMP_SIMD_PRAGMA
3296 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3297 solution_old_ptr[i] =
3298 factor1_plus_1 * solution_ptr[i] - solution_old_ptr[i] * factor2;
3299 }
3300 else
3301 {
3302 // compute t = P^{-1} * (A*x^{n} - b)
3303 if (compute_residual_norm)
3304 residual_norm =
3305 std::sqrt(temp_vector1.add_and_dot(-1.0, rhs, temp_vector1));
3306 else
3307 temp_vector1.add(-1.0, rhs);
3308
3309 preconditioner.vmult(temp_vector2, temp_vector1);
3310
3311 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) - f_2 * t
3312 const auto solution_ptr = solution.begin();
3313 const auto solution_old_ptr = solution_old.begin();
3314 const auto temp_vector2_ptr = temp_vector2.begin();
3315 DEAL_II_OPENMP_SIMD_PRAGMA
3316 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3317 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3318 factor1 * solution_old_ptr[i] -
3319 temp_vector2_ptr[i] * factor2;
3320 }
3321
3322 solution.swap(solution_old);
3323
3324 return residual_norm;
3325 }
3326
3327 template <
3328 typename Number,
3329 typename PreconditionerType,
3330 std::enable_if_t<
3331 has_vmult_with_std_functions_for_precondition<
3332 PreconditionerType,
3333 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>>,
3334 int> * = nullptr>
3335 inline double
3336 vector_updates(
3337 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3338 const PreconditionerType &preconditioner,
3339 const unsigned int iteration_index,
3340 const double factor1_,
3341 const double factor2_,
3342 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3343 &solution_old,
3344 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3345 &temp_vector1,
3346 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3347 &temp_vector2,
3348 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution,
3349 const bool compute_residual_norm)
3350 {
3351 const Number factor1 = factor1_;
3352 const Number factor1_plus_1 = 1. + factor1_;
3353 const Number factor2 = factor2_;
3354 double residual_norm = 0;
3355
3356 const auto rhs_ptr = rhs.begin();
3357 const auto temp_vector1_ptr = temp_vector1.begin();
3358 const auto temp_vector2_ptr = temp_vector2.begin();
3359 const auto solution_ptr = solution.begin();
3360 const auto solution_old_ptr = solution_old.begin();
3361
3362 if (iteration_index == 0)
3363 {
3364 if (compute_residual_norm)
3365 residual_norm = rhs.l2_norm();
3366 preconditioner.vmult(
3367 temp_vector2,
3368 rhs,
3369 [&](const auto start_range, const auto end_range) {
3370 if (end_range > start_range)
3371 std::memset(temp_vector2_ptr + start_range,
3372 0,
3373 sizeof(Number) * (end_range - start_range));
3374 },
3375 [&](const auto begin, const auto end) {
3376 DEAL_II_OPENMP_SIMD_PRAGMA
3377 for (std::size_t i = begin; i < end; ++i)
3378 temp_vector2_ptr[i] *= factor2;
3379 });
3380 }
3381 else
3382 {
3383 VectorizedArray<Number> sum = 0;
3384 const auto post_update = [&](const unsigned int begin,
3385 const unsigned int end) {
3386 if (iteration_index == 1)
3387 {
3388 DEAL_II_OPENMP_SIMD_PRAGMA
3389 for (std::size_t i = begin; i < end; ++i)
3390 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
3391 factor2 * temp_vector2_ptr[i];
3392 }
3393 else
3394 {
3395 DEAL_II_OPENMP_SIMD_PRAGMA
3396 for (std::size_t i = begin; i < end; ++i)
3397 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3398 factor1 * solution_old_ptr[i] +
3399 factor2 * temp_vector2_ptr[i];
3400 }
3401 };
3402 if (compute_residual_norm)
3403 {
3404 preconditioner.vmult(
3405 temp_vector2,
3406 temp_vector1,
3407 [&](const auto begin, const auto end) {
3408 if (end > begin)
3409 std::memset(temp_vector2_ptr + begin,
3410 0,
3411 sizeof(Number) * (end - begin));
3412
3413 VectorizedArray<Number> local_sum = 0;
3414 constexpr unsigned int n_lanes =
3415 VectorizedArray<Number>::size();
3416 const unsigned int end_regular = end / n_lanes * n_lanes;
3417 for (unsigned int i = begin; i < end_regular; i += n_lanes)
3418 {
3419 VectorizedArray<Number> rhs_i, tmp_i;
3420 rhs_i.load(rhs_ptr + i);
3421 tmp_i.load(temp_vector1_ptr + i);
3422 const VectorizedArray<Number> residual_i = rhs_i - tmp_i;
3423 residual_i.store(temp_vector1_ptr + i);
3424 local_sum += residual_i * residual_i;
3425 }
3426 for (unsigned int i = end_regular; i < end; ++i)
3427 {
3428 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3429 local_sum[i - end_regular] +=
3430 temp_vector1_ptr[i] * temp_vector1_ptr[i];
3431 }
3432 sum += local_sum;
3433 },
3434 post_update);
3435 residual_norm = std::sqrt(
3436 Utilities::MPI::sum(sum.sum(),
3437 solution_old.get_mpi_communicator()));
3438 }
3439 else
3440 preconditioner.vmult(
3441 temp_vector2,
3442 temp_vector1,
3443 [&](const auto begin, const auto end) {
3444 if (end > begin)
3445 std::memset(temp_vector2_ptr + begin,
3446 0,
3447 sizeof(Number) * (end - begin));
3448
3449 for (std::size_t i = begin; i < end; ++i)
3450 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3451 },
3452 post_update);
3453 }
3454
3455 solution_old.swap(temp_vector2);
3456 solution_old.swap(solution);
3457
3458 return residual_norm;
3459 }
3460
3461 // vector updates for device vectors and DiagonalMatrix as preconditioner
3462 template <typename Number>
3463 inline double
3464 vector_updates(
3465 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>
3466 &rhs,
3467 const ::DiagonalMatrix<
3468 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>>
3469 &preconditioner,
3470 const unsigned int iteration_index,
3471 const double factor1_,
3472 const double factor2_,
3473 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>
3474 &solution_old,
3475 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>
3476 &temp_vector1,
3477 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default> &,
3478 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>
3479 &solution,
3480 const bool compute_residual_norm)
3481 {
3482 const Number factor1 = factor1_;
3483 const Number factor1_plus_1 = 1. + factor1_;
3484 const Number factor2 = factor2_;
3485 double residual_norm = 0;
3486
3487 auto exec = typename ::MemorySpace::Default::kokkos_space::
3488 execution_space{};
3489 Number *sol_data = solution.begin();
3490 const Number *rhs_data = rhs.begin();
3491 const Number *prec_data = preconditioner.get_vector().begin();
3492 Number *sol_old_data = solution_old.begin();
3493 Number *tmp_data = temp_vector1.begin();
3494 using size_type = typename LinearAlgebra::distributed::
3495 Vector<Number, MemorySpace::Default>::size_type;
3496 const size_type size = solution.locally_owned_size();
3497 if (iteration_index == 0)
3498 {
3499 Kokkos::parallel_for(
3500 "::ChebyshevIt0",
3501 Kokkos::RangePolicy<
3502 ::MemorySpace::Default::kokkos_space::execution_space>(
3503 exec, 0, size),
3504 KOKKOS_LAMBDA(size_type i) {
3505 sol_old_data[i] = prec_data[i] * factor2 * rhs_data[i];
3506 });
3507 if (compute_residual_norm)
3508 residual_norm = rhs.l2_norm();
3509 }
3510 else if (iteration_index == 1)
3511 {
3512 Kokkos::parallel_for(
3513 "::ChebyshevIt1",
3514 Kokkos::RangePolicy<
3515 ::MemorySpace::Default::kokkos_space::execution_space>(
3516 exec, 0, size),
3517 KOKKOS_LAMBDA(size_type i) {
3518 const Number residual = rhs_data[i] - tmp_data[i];
3519 sol_old_data[i] = factor2 * prec_data[i] * residual +
3520 factor1_plus_1 * sol_data[i];
3521 });
3522 if (compute_residual_norm)
3523 residual_norm =
3524 std::sqrt(temp_vector1.add_and_dot(-1.0, rhs, temp_vector1));
3525 }
3526 else
3527 {
3528 Kokkos::parallel_for(
3529 "::ChebyshevItk",
3530 Kokkos::RangePolicy<
3531 ::MemorySpace::Default::kokkos_space::execution_space>(
3532 exec, 0, size),
3533 KOKKOS_LAMBDA(size_type i) {
3534 const Number residual = rhs_data[i] - tmp_data[i];
3535 sol_old_data[i] = factor2 * prec_data[i] * residual -
3536 factor1 * sol_old_data[i] +
3537 factor1_plus_1 * sol_data[i];
3538 });
3539 if (compute_residual_norm)
3540 residual_norm =
3541 std::sqrt(temp_vector1.add_and_dot(-1.0, rhs, temp_vector1));
3542 }
3543 exec.fence();
3544
3545 solution.swap(solution_old);
3546 return residual_norm;
3547 }
3548
3549 // worker routine for deal.II vectors. Because of vectorization, we need
3550 // to put the loop into an extra structure because the virtual function of
3551 // VectorUpdatesRange prevents the compiler from applying vectorization.
3552 template <typename Number>
3553 struct VectorUpdater
3554 {
3555 VectorUpdater(const Number *rhs,
3556 const Number *matrix_diagonal_inverse,
3557 const unsigned int iteration_index,
3558 const Number factor1,
3559 const Number factor2,
3560 Number *solution_old,
3561 Number *tmp_vector,
3562 Number *solution,
3563 const bool compute_residual_norm)
3564 : rhs(rhs)
3565 , matrix_diagonal_inverse(matrix_diagonal_inverse)
3566 , iteration_index(iteration_index)
3567 , factor1(factor1)
3568 , factor2(factor2)
3569 , solution_old(solution_old)
3570 , tmp_vector(tmp_vector)
3571 , solution(solution)
3572 , compute_residual_norm(compute_residual_norm)
3573 , sum_accumulator(0.0)
3574 {}
3575
3576 void
3577 apply_to_subrange(const std::size_t begin, const std::size_t end) const
3578 {
3579 // Create local copies to help the aliasing detection.
3580 const Number factor1 = this->factor1;
3581 const Number factor1_plus_1 = 1. + this->factor1;
3582 const Number factor2 = this->factor2;
3583 if (compute_residual_norm)
3584 {
3585 VectorizedArray<Number> local_sum = 0;
3586 constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
3587 const unsigned int end_regular = end / n_lanes * n_lanes;
3588 if (iteration_index == 0)
3589 {
3590 for (unsigned int i = begin; i < end_regular; i += n_lanes)
3591 {
3592 VectorizedArray<Number> rhs_i, diag_i;
3593 rhs_i.load(rhs + i);
3594 diag_i.load(matrix_diagonal_inverse + i);
3595 const VectorizedArray<Number> res_i =
3596 factor2 * diag_i * rhs_i;
3597 res_i.store(tmp_vector + i);
3598 local_sum += rhs_i * rhs_i;
3599 }
3600 for (unsigned int i = end_regular; i < end; ++i)
3601 {
3602 tmp_vector[i] =
3603 factor2 * matrix_diagonal_inverse[i] * rhs[i];
3604 local_sum[i - end_regular] += rhs[i] * rhs[i];
3605 }
3606 }
3607 else if (iteration_index == 1)
3608 {
3609 // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
3610 for (unsigned int i = begin; i < end_regular; i += n_lanes)
3611 {
3612 VectorizedArray<Number> rhs_i, sol_i, tmp_i, diag_i;
3613 rhs_i.load(rhs + i);
3614 diag_i.load(matrix_diagonal_inverse + i);
3615 sol_i.load(solution + i);
3616 tmp_i.load(tmp_vector + i);
3617 const VectorizedArray<Number> residual_i = rhs_i - tmp_i;
3618 local_sum += residual_i * residual_i;
3619 const VectorizedArray<Number> res_i =
3620 factor1_plus_1 * sol_i + factor2 * diag_i * residual_i;
3621 res_i.store(tmp_vector + i);
3622 }
3623 for (unsigned int i = end_regular; i < end; ++i)
3624 {
3625 const Number residual_i = rhs[i] - tmp_vector[i];
3626 local_sum[i - end_regular] += residual_i * residual_i;
3627 tmp_vector[i] =
3628 factor1_plus_1 * solution[i] +
3629 factor2 * matrix_diagonal_inverse[i] * residual_i;
3630 }
3631 }
3632 else
3633 {
3634 // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
3635 // + f_2 * P^{-1} * (b-A*x^{n})
3636 for (unsigned int i = begin; i < end_regular; i += n_lanes)
3637 {
3638 VectorizedArray<Number> rhs_i, sol_i, sol_old_i, tmp_i,
3639 diag_i;
3640 rhs_i.load(rhs + i);
3641 diag_i.load(matrix_diagonal_inverse + i);
3642 sol_i.load(solution + i);
3643 sol_old_i.load(solution_old + i);
3644 tmp_i.load(tmp_vector + i);
3645 const VectorizedArray<Number> residual_i = rhs_i - tmp_i;
3646 local_sum += residual_i * residual_i;
3647 const VectorizedArray<Number> res_i =
3648 factor1_plus_1 * sol_i - factor1 * sol_old_i +
3649 factor2 * diag_i * residual_i;
3650 res_i.store(tmp_vector + i);
3651 }
3652 for (unsigned int i = end_regular; i < end; ++i)
3653 {
3654 const Number residual_i = rhs[i] - tmp_vector[i];
3655 local_sum[i - end_regular] += residual_i * residual_i;
3656 tmp_vector[i] =
3657 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
3658 factor2 * matrix_diagonal_inverse[i] * residual_i;
3659 }
3660 }
3661 sum_accumulator += local_sum;
3662 }
3663 else
3664 {
3665 if (iteration_index == 0)
3666 {
3667 DEAL_II_OPENMP_SIMD_PRAGMA
3668 for (std::size_t i = begin; i < end; ++i)
3669 tmp_vector[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
3670 }
3671 else if (iteration_index == 1)
3672 {
3673 // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
3674 DEAL_II_OPENMP_SIMD_PRAGMA
3675 for (std::size_t i = begin; i < end; ++i)
3676 // for efficiency reason, write back to temp_vector that is
3677 // already read (avoid read-for-ownership)
3678 tmp_vector[i] = factor1_plus_1 * solution[i] +
3679 factor2 * matrix_diagonal_inverse[i] *
3680 (rhs[i] - tmp_vector[i]);
3681 }
3682 else
3683 {
3684 // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
3685 // + f_2 * P^{-1} * (b-A*x^{n})
3686 DEAL_II_OPENMP_SIMD_PRAGMA
3687 for (std::size_t i = begin; i < end; ++i)
3688 // for efficiency reason, write back to temp_vector, which is
3689 // already modified during vmult (in best case, the modified
3690 // values are not written back to main memory yet so that
3691 // we do not have to pay additional costs for writing and
3692 // read-for-ownership)
3693 tmp_vector[i] = factor1_plus_1 * solution[i] -
3694 factor1 * solution_old[i] +
3695 factor2 * matrix_diagonal_inverse[i] *
3696 (rhs[i] - tmp_vector[i]);
3697 }
3698 }
3699 }
3700
3701 const Number *rhs;
3702 const Number *matrix_diagonal_inverse;
3703 const unsigned int iteration_index;
3704 const Number factor1;
3705 const Number factor2;
3706 mutable Number *solution_old;
3707 mutable Number *tmp_vector;
3708 mutable Number *solution;
3709 bool compute_residual_norm;
3710 mutable VectorizedArray<Number> sum_accumulator;
3711 };
3712
3713 template <typename Number>
3714 struct VectorUpdatesRange : public ::parallel::ParallelForInteger
3715 {
3716 VectorUpdatesRange(const VectorUpdater<Number> &updater,
3717 const std::size_t size)
3718 : updater(updater)
3719 {
3720 Assert(updater.compute_residual_norm == false ||
3721 MultithreadInfo::n_threads() == 1,
3722 ExcNotImplemented(
3723 "Currently, we cannot compute residual in parallel"));
3724 if (size <
3725 internal::VectorImplementation::minimum_parallel_grain_size ||
3726 MultithreadInfo::n_threads() == 1)
3727 VectorUpdatesRange::apply_to_subrange(0, size);
3728 else
3729 apply_parallel(
3730 0,
3731 size,
3732 internal::VectorImplementation::minimum_parallel_grain_size);
3733 }
3734
3735 ~VectorUpdatesRange() override = default;
3736
3737 virtual void
3738 apply_to_subrange(const std::size_t begin,
3739 const std::size_t end) const override
3740 {
3741 updater.apply_to_subrange(begin, end);
3742 }
3743
3744 const VectorUpdater<Number> &updater;
3745 };
3746
3747 // selection for diagonal matrix around deal.II vector
3748 template <typename Number>
3749 inline double
3750 vector_updates(
3751 const ::Vector<Number> &rhs,
3752 const ::DiagonalMatrix<::Vector<Number>> &jacobi,
3753 const unsigned int iteration_index,
3754 const double factor1,
3755 const double factor2,
3756 ::Vector<Number> &solution_old,
3757 ::Vector<Number> &temp_vector1,
3758 ::Vector<Number> &,
3759 ::Vector<Number> &solution,
3760 const bool compute_residual_norm)
3761 {
3762 VectorUpdater<Number> upd(rhs.begin(),
3763 jacobi.get_vector().begin(),
3764 iteration_index,
3765 factor1,
3766 factor2,
3767 solution_old.begin(),
3768 temp_vector1.begin(),
3769 solution.begin(),
3770 compute_residual_norm);
3771
3772 VectorUpdatesRange<Number>(upd, rhs.size());
3773
3774 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3775 solution.swap(temp_vector1);
3776 solution_old.swap(temp_vector1);
3777
3778 if (compute_residual_norm)
3779 return std::sqrt(
3780 Utilities::MPI::sum(upd.sum_accumulator.sum(),
3781 solution_old.get_mpi_communicator()));
3782 else
3783 return 0;
3784 }
3785
3786 // selection for diagonal matrix around parallel deal.II vector
3787 template <typename Number>
3788 inline double
3789 vector_updates(
3790 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3791 const ::DiagonalMatrix<
3792 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>> &jacobi,
3793 const unsigned int iteration_index,
3794 const double factor1,
3795 const double factor2,
3796 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3797 &solution_old,
3798 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3799 &temp_vector1,
3800 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &,
3801 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution,
3802 const bool compute_residual_norm)
3803 {
3804 VectorUpdater<Number> upd(rhs.begin(),
3805 jacobi.get_vector().begin(),
3806 iteration_index,
3807 factor1,
3808 factor2,
3809 solution_old.begin(),
3810 temp_vector1.begin(),
3811 solution.begin(),
3812 compute_residual_norm);
3813
3814 VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
3815
3816 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3817 solution.swap(temp_vector1);
3818 solution_old.swap(temp_vector1);
3819
3820 if (compute_residual_norm)
3821 return std::sqrt(
3822 Utilities::MPI::sum(upd.sum_accumulator.sum(),
3823 solution_old.get_mpi_communicator()));
3824 else
3825 return 0;
3826 }
3827
3828 // We need to have a separate declaration for static const members
3829
3830 // general case and the case that the preconditioner can work on
3831 // ranges (covered by vector_updates())
3832 template <
3833 typename MatrixType,
3834 typename VectorType,
3835 typename PreconditionerType,
3836 std::enable_if_t<
3837 !has_vmult_with_std_functions<MatrixType,
3838 VectorType,
3839 PreconditionerType> &&
3840 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3841 VectorType> &&
3842 has_vmult_with_std_functions_for_precondition<MatrixType,
3843 VectorType>),
3844 int> * = nullptr>
3845 inline double
3846 vmult_and_update(const MatrixType &matrix,
3847 const PreconditionerType &preconditioner,
3848 const VectorType &rhs,
3849 const unsigned int iteration_index,
3850 const double factor1,
3851 const double factor2,
3852 VectorType &solution,
3853 VectorType &solution_old,
3854 VectorType &temp_vector1,
3855 VectorType &temp_vector2,
3856 const bool compute_residual_norm)
3857 {
3858 if (iteration_index > 0)
3859 matrix.vmult(temp_vector1, solution);
3860 return vector_updates(rhs,
3861 preconditioner,
3862 iteration_index,
3863 factor1,
3864 factor2,
3865 solution_old,
3866 temp_vector1,
3867 temp_vector2,
3868 solution,
3869 compute_residual_norm);
3870 }
3871
3872 // case that both the operator and the preconditioner can work on
3873 // subranges
3874 template <
3875 typename MatrixType,
3876 typename VectorType,
3877 typename PreconditionerType,
3878 std::enable_if_t<
3879 !has_vmult_with_std_functions<MatrixType,
3880 VectorType,
3881 PreconditionerType> &&
3882 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3883 VectorType> &&
3884 has_vmult_with_std_functions_for_precondition<MatrixType,
3885 VectorType>),
3886 int> * = nullptr>
3887 inline double
3888 vmult_and_update(const MatrixType &matrix,
3889 const PreconditionerType &preconditioner,
3890 const VectorType &rhs,
3891 const unsigned int iteration_index,
3892 const double factor1_,
3893 const double factor2_,
3894 VectorType &solution,
3895 VectorType &solution_old,
3896 VectorType &temp_vector1,
3897 VectorType &temp_vector2,
3898 const bool compute_residual_norm)
3899 {
3900 using Number = typename VectorType::value_type;
3901
3902 const Number factor1 = factor1_;
3903 const Number factor1_plus_1 = 1. + factor1_;
3904 const Number factor2 = factor2_;
3905 double residual_norm = 0;
3906
3907 if (iteration_index == 0)
3908 {
3909 if (compute_residual_norm)
3910 residual_norm = rhs.l2_norm();
3911
3912 preconditioner.vmult(
3913 temp_vector2,
3914 rhs,
3915 [&](const unsigned int start_range, const unsigned int end_range) {
3916 // zero 'solution' before running the vmult operation
3917 if (end_range > start_range)
3918 std::memset(temp_vector2.begin() + start_range,
3919 0,
3920 sizeof(Number) * (end_range - start_range));
3921 },
3922 [&](const unsigned int start_range, const unsigned int end_range) {
3923 const auto tmp_ptr = temp_vector2.begin();
3924
3925 DEAL_II_OPENMP_SIMD_PRAGMA
3926 for (std::size_t i = start_range; i < end_range; ++i)
3927 tmp_ptr[i] *= factor2;
3928 });
3929 }
3930 else
3931 {
3932 temp_vector1.reinit(rhs, true);
3933 temp_vector2.reinit(rhs, true);
3934
3935 VectorizedArray<Number> sum = 0;
3936 const auto post_operation = [&](const unsigned int start_range,
3937 const unsigned int end_range) {
3938 const auto rhs_ptr = rhs.begin();
3939 const auto tmp_ptr = temp_vector1.begin();
3940
3941 DEAL_II_OPENMP_SIMD_PRAGMA
3942 for (std::size_t i = start_range; i < end_range; ++i)
3943 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3944 };
3945 // 1) compute residual (including operator application)
3946 if (compute_residual_norm)
3947 matrix.vmult(
3948 temp_vector1,
3949 solution,
3950 [&](const unsigned int begin, const unsigned int end) {
3951 const auto rhs_ptr = rhs.begin();
3952 const auto tmp_ptr = temp_vector1.begin();
3953 VectorizedArray<Number> local_sum = 0;
3954 constexpr unsigned int n_lanes =
3955 VectorizedArray<Number>::size();
3956 const unsigned int end_regular = end / n_lanes * n_lanes;
3957 for (unsigned int i = begin; i < end_regular; i += n_lanes)
3958 {
3959 VectorizedArray<Number> rhs_i, tmp_i;
3960 rhs_i.load(rhs_ptr + i);
3961 tmp_i.load(tmp_ptr + i);
3962 const VectorizedArray<Number> residual_i = rhs_i - tmp_i;
3963 residual_i.store(tmp_ptr + i);
3964 local_sum += residual_i * residual_i;
3965 }
3966 for (unsigned int i = end_regular; i < end; ++i)
3967 {
3968 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3969 local_sum[i - end_regular] += tmp_ptr[i] * tmp_ptr[i];
3970 }
3971 sum += local_sum;
3972 },
3973 post_operation
3974
3975 );
3976 else
3977 matrix.vmult(
3978 temp_vector1,
3979 solution,
3980 [&](const unsigned int start_range,
3981 const unsigned int end_range) {
3982 const auto rhs_ptr = rhs.begin();
3983 const auto tmp_ptr = temp_vector1.begin();
3984
3985 DEAL_II_OPENMP_SIMD_PRAGMA
3986 for (std::size_t i = start_range; i < end_range; ++i)
3987 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3988 },
3989 post_operation
3990
3991 );
3992
3993 // 2) perform vector updates (including preconditioner application)
3994 preconditioner.vmult(
3995 temp_vector2,
3996 temp_vector1,
3997 [&](const unsigned int start_range, const unsigned int end_range) {
3998 // zero 'temp_vector2' before running the vmult
3999 // operation
4000 if (end_range > start_range)
4001 std::memset(temp_vector2.begin() + start_range,
4002 0,
4003 sizeof(Number) * (end_range - start_range));
4004 },
4005 [&](const unsigned int start_range, const unsigned int end_range) {
4006 const auto solution_ptr = solution.begin();
4007 const auto solution_old_ptr = solution_old.begin();
4008 const auto tmp_ptr = temp_vector2.begin();
4009
4010 if (iteration_index == 1)
4011 {
4012 DEAL_II_OPENMP_SIMD_PRAGMA
4013 for (std::size_t i = start_range; i < end_range; ++i)
4014 tmp_ptr[i] =
4015 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
4016 }
4017 else
4018 {
4019 DEAL_II_OPENMP_SIMD_PRAGMA
4020 for (std::size_t i = start_range; i < end_range; ++i)
4021 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
4022 factor1 * solution_old_ptr[i] +
4023 factor2 * tmp_ptr[i];
4024 }
4025 });
4026
4027 if (compute_residual_norm)
4028 residual_norm = std::sqrt(
4029 Utilities::MPI::sum(sum.sum(), solution.get_mpi_communicator()));
4030 }
4031 solution.swap(temp_vector2);
4032 solution_old.swap(temp_vector2);
4033 return residual_norm;
4034 }
4035
4036 // case that the operator can work on subranges and the preconditioner
4037 // is a diagonal
4038 template <typename MatrixType,
4039 typename VectorType,
4040 typename PreconditionerType,
4041 std::enable_if_t<has_vmult_with_std_functions<MatrixType,
4042 VectorType,
4043 PreconditionerType>,
4044 int> * = nullptr>
4045 inline double
4046 vmult_and_update(const MatrixType &matrix,
4047 const PreconditionerType &preconditioner,
4048 const VectorType &rhs,
4049 const unsigned int iteration_index,
4050 const double factor1,
4051 const double factor2,
4052 VectorType &solution,
4053 VectorType &solution_old,
4054 VectorType &temp_vector1,
4055 VectorType &,
4056 const bool compute_residual_norm)
4057 {
4058 using Number = typename VectorType::value_type;
4059 VectorUpdater<Number> updater(rhs.begin(),
4060 preconditioner.get_vector().begin(),
4061 iteration_index,
4062 factor1,
4063 factor2,
4064 solution_old.begin(),
4065 temp_vector1.begin(),
4066 solution.begin(),
4067 compute_residual_norm);
4068 if (iteration_index > 0)
4069 matrix.vmult(
4070 temp_vector1,
4071 solution,
4072 [&](const unsigned int start_range, const unsigned int end_range) {
4073 // zero 'temp_vector1' before running the vmult
4074 // operation
4075 if (end_range > start_range)
4076 std::memset(temp_vector1.begin() + start_range,
4077 0,
4078 sizeof(Number) * (end_range - start_range));
4079 },
4080 [&](const unsigned int start_range, const unsigned int end_range) {
4081 if (end_range > start_range)
4082 updater.apply_to_subrange(start_range, end_range);
4083 });
4084 else
4085 updater.apply_to_subrange(0U, solution.locally_owned_size());
4086
4087 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
4088 solution.swap(temp_vector1);
4089 solution_old.swap(temp_vector1);
4090
4091 if (compute_residual_norm)
4092 return std::sqrt(Utilities::MPI::sum(updater.sum_accumulator.sum(),
4093 solution.get_mpi_communicator()));
4094 else
4095 return 0;
4096 }
4097
4098 template <typename MatrixType, typename PreconditionerType>
4099 inline void
4100 initialize_preconditioner(
4101 const MatrixType & /*matrix*/,
4102 std::shared_ptr<PreconditionerType> &preconditioner)
4103 {
4104 (void)preconditioner;
4105 AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
4106 }
4107
4108 template <typename MatrixType, typename VectorType>
4109 inline void
4110 initialize_preconditioner(
4111 const MatrixType &matrix,
4112 std::shared_ptr<::DiagonalMatrix<VectorType>> &preconditioner)
4113 {
4114 if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
4115 {
4116 if (preconditioner.get() == nullptr)
4117 preconditioner =
4118 std::make_shared<::DiagonalMatrix<VectorType>>();
4119
4120 Assert(
4121 preconditioner->m() == 0,
4122 ExcMessage(
4123 "Preconditioner appears to be initialized but not sized correctly"));
4124
4125 // This part only works in serial
4126 if (preconditioner->m() != matrix.m())
4127 {
4128 preconditioner->get_vector().reinit(matrix.m());
4129 for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
4130 preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
4131 }
4132 }
4133 }
4134 } // namespace PreconditionChebyshevImplementation
4135} // namespace internal
4136
4137
4138
4139template <typename MatrixType, typename VectorType, typename PreconditionerType>
4140inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
4141 AdditionalData::AdditionalData(const unsigned int degree,
4142 const double smoothing_range,
4143 const unsigned int eig_cg_n_iterations,
4144 const double eig_cg_residual,
4145 const double max_eigenvalue,
4146 const EigenvalueAlgorithm eigenvalue_algorithm,
4147 const PolynomialType polynomial_type,
4148 const double safety_factor)
4149 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
4150 smoothing_range,
4151 eig_cg_n_iterations,
4152 eig_cg_residual,
4153 max_eigenvalue,
4154 eigenvalue_algorithm,
4155 safety_factor)
4156 , degree(degree)
4157 , polynomial_type(polynomial_type)
4158{}
4159
4160
4161
4162template <typename MatrixType, typename VectorType, typename PreconditionerType>
4163inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
4164 PreconditionChebyshev()
4165 : theta(1.)
4166 , delta(1.)
4167 , eigenvalues_are_initialized(false)
4168{
4169 static_assert(
4170 std::is_same_v<size_type, typename VectorType::size_type>,
4171 "PreconditionChebyshev and VectorType must have the same size_type.");
4172}
4173
4174
4175
4176template <typename MatrixType, typename VectorType, typename PreconditionerType>
4177inline void
4178PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::initialize(
4179 const MatrixType &matrix,
4180 const AdditionalData &additional_data)
4181{
4182 matrix_ptr = &matrix;
4183 data = additional_data;
4184 Assert(data.degree > 0,
4185 ExcMessage("The degree of the Chebyshev method must be positive."));
4186 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
4187 matrix, data.preconditioner);
4188 eigenvalues_are_initialized = false;
4189}
4190
4191
4192
4193template <typename MatrixType, typename VectorType, typename PreconditionerType>
4194inline void
4195PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::clear()
4196{
4197 eigenvalues_are_initialized = false;
4198 theta = delta = 1.0;
4199 matrix_ptr = nullptr;
4200 {
4201 VectorType empty_vector;
4202 solution_old.reinit(empty_vector);
4203 temp_vector1.reinit(empty_vector);
4204 temp_vector2.reinit(empty_vector);
4205 }
4206 data.preconditioner.reset();
4207}
4208
4209
4210
4211template <typename MatrixType, typename VectorType, typename PreconditionerType>
4212inline typename internal::EigenvalueInformation
4213PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
4214 estimate_eigenvalues(const VectorType &src) const
4215{
4216 Assert(eigenvalues_are_initialized == false, ExcInternalError());
4217
4218 solution_old.reinit(src);
4219 temp_vector1.reinit(src, true);
4220
4221 auto info = internal::estimate_eigenvalues<MatrixType>(data,
4222 matrix_ptr,
4223 solution_old,
4224 temp_vector1,
4225 data.degree,
4226 data.safety_factor);
4227
4228 const double alpha = (data.smoothing_range > 1. ?
4229 info.max_eigenvalue_estimate / data.smoothing_range :
4230 std::min(0.9 * info.max_eigenvalue_estimate,
4231 info.min_eigenvalue_estimate));
4232
4233 // in case the user set the degree to invalid unsigned int, we have to
4234 // determine the number of necessary iterations from the Chebyshev error
4235 // estimate, given the target tolerance specified by smoothing_range. This
4236 // estimate is based on the error formula given in section 5.1 of
4237 // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
4238 if (data.degree == numbers::invalid_unsigned_int)
4239 {
4240 // In solver mode, smoothing_range is interpreted as the relative
4241 // target tolerance and must be strictly less than one. Otherwise the
4242 // Chebyshev error formula below evaluates the square root of a
4243 // negative number, producing a NaN that silently gets cast to an
4244 // unsigned int with implementation-defined (typically zero) result.
4245 // This would turn the "solver" into a single damped Jacobi step,
4246 // which is almost certainly not what the user intended.
4247 Assert(data.smoothing_range < 1.,
4248 ExcMessage(
4249 "When using PreconditionChebyshev as a solver (i.e., when "
4250 "AdditionalData::degree == numbers::invalid_unsigned_int), "
4251 "AdditionalData::smoothing_range must be strictly less than "
4252 "one because it is interpreted as the relative target "
4253 "tolerance of the Chebyshev iteration. A value >= 1 is only "
4254 "meaningful when PreconditionChebyshev is used as a smoother "
4255 "with a fixed polynomial degree."));
4256
4257 const double actual_range = info.max_eigenvalue_estimate / alpha;
4258 const double sigma = (1. - std::sqrt(1. / actual_range)) /
4259 (1. + std::sqrt(1. / actual_range));
4260 const double eps = data.smoothing_range;
4261 const_cast<
4262 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(
4263 this)
4264 ->data.degree =
4265 1 + static_cast<unsigned int>(
4266 std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
4267 std::log(1. / sigma));
4268 }
4269
4270 info.degree = data.degree;
4271
4272 const_cast<
4273 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
4274 ->delta =
4275 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4276 (info.max_eigenvalue_estimate) :
4277 ((info.max_eigenvalue_estimate - alpha) * 0.5);
4278 const_cast<
4279 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
4280 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
4281
4282 // We do not need the second temporary vector in case we have a
4283 // DiagonalMatrix as preconditioner and use deal.II's own vectors
4284 using NumberType = typename VectorType::value_type;
4285 if (std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> ==
4286 false ||
4287 (std::is_same_v<VectorType, ::Vector<NumberType>> == false &&
4288 ((std::is_same_v<
4289 VectorType,
4290 LinearAlgebra::distributed::Vector<NumberType, MemorySpace::Host>> ==
4291 false) ||
4292 (std::is_same_v<VectorType,
4293 LinearAlgebra::distributed::
4294 Vector<NumberType, MemorySpace::Default>> == false))))
4295 temp_vector2.reinit(src, true);
4296 else
4297 {
4298 VectorType empty_vector;
4299 temp_vector2.reinit(empty_vector);
4300 }
4301
4302 const_cast<
4303 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
4304 ->eigenvalues_are_initialized = true;
4305
4306 return info;
4307}
4308
4309
4310
4311template <typename MatrixType, typename VectorType, typename PreconditionerType>
4312template <bool do_transpose>
4313inline double
4314PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
4315 apply_internal(const bool zero_out_dst,
4316 const VectorType &rhs,
4317 VectorType &solution,
4318 const bool compute_final_norm) const
4319{
4320 std::scoped_lock lock(mutex);
4321 if (eigenvalues_are_initialized == false)
4322 estimate_eigenvalues(rhs);
4323
4324 double residual_norm = 0;
4325 if constexpr (do_transpose)
4326 {
4327 matrix_ptr->Tvmult(temp_vector1, solution);
4328 residual_norm =
4329 internal::PreconditionChebyshevImplementation::vector_updates(
4330 rhs,
4331 *data.preconditioner,
4332 zero_out_dst ? 0 : 1,
4333 0.,
4334 (data.polynomial_type ==
4335 AdditionalData::PolynomialType::fourth_kind) ?
4336 (4. / (3. * delta)) :
4337 (1. / theta),
4338 solution_old,
4339 temp_vector1,
4340 temp_vector2,
4341 solution,
4342 compute_final_norm && data.degree < 2);
4343 }
4344 else
4345 residual_norm =
4346 internal::PreconditionChebyshevImplementation::vmult_and_update(
4347 *matrix_ptr,
4348 *data.preconditioner,
4349 rhs,
4350 zero_out_dst ? 0 : 1,
4351 0.,
4352 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4353 (4. / (3. * delta)) :
4354 (1. / theta),
4355 solution,
4356 solution_old,
4357 temp_vector1,
4358 temp_vector2,
4359 compute_final_norm && data.degree < 2);
4360
4361 // if delta is zero, we do not need to iterate because the updates will be
4362 // zero
4363 if (data.degree < 2 || std::abs(delta) < 1e-40)
4364 return residual_norm;
4365
4366 double rhok = delta / theta, sigma = theta / delta;
4367 for (unsigned int k = 0; k < data.degree - 1; ++k)
4368 {
4369 double factor1 = 0.0;
4370 double factor2 = 0.0;
4371
4372 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4373 {
4374 factor1 = (2 * k + 1.) / (2 * k + 5.);
4375 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4376 }
4377 else
4378 {
4379 const double rhokp = 1. / (2. * sigma - rhok);
4380 factor1 = rhokp * rhok;
4381 factor2 = 2. * rhokp / delta;
4382 rhok = rhokp;
4383 }
4384
4385 if constexpr (do_transpose)
4386 {
4387 matrix_ptr->Tvmult(temp_vector1, solution);
4388 residual_norm =
4389 internal::PreconditionChebyshevImplementation::vector_updates(
4390 rhs,
4391 *data.preconditioner,
4392 k + (zero_out_dst ? 1 : 2),
4393 factor1,
4394 factor2,
4395 solution_old,
4396 temp_vector1,
4397 temp_vector2,
4398 solution,
4399 compute_final_norm && k == data.degree - 2);
4400 }
4401 else
4402 {
4403 residual_norm =
4404 internal::PreconditionChebyshevImplementation::vmult_and_update(
4405 *matrix_ptr,
4406 *data.preconditioner,
4407 rhs,
4408 k + (zero_out_dst ? 1 : 2),
4409 factor1,
4410 factor2,
4411 solution,
4412 solution_old,
4413 temp_vector1,
4414 temp_vector2,
4415 compute_final_norm && k == data.degree - 2);
4416 }
4417 }
4418
4419 return residual_norm;
4420}
4421
4422
4423
4424template <typename MatrixType, typename VectorType, typename PreconditionerType>
4425inline void
4426PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::vmult(
4427 VectorType &solution,
4428 const VectorType &rhs) const
4429{
4430 apply_internal<false>(/* zero_out_dst */ true, rhs, solution);
4431}
4432
4433
4434
4435template <typename MatrixType, typename VectorType, typename PreconditionerType>
4436inline void
4437PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tvmult(
4438 VectorType &solution,
4439 const VectorType &rhs) const
4440{
4441 apply_internal<true>(/* zero_out_dst */ true, rhs, solution);
4442}
4443
4444
4445
4446template <typename MatrixType, typename VectorType, typename PreconditionerType>
4447inline void
4448PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::step(
4449 VectorType &solution,
4450 const VectorType &rhs) const
4451{
4452 apply_internal<false>(/* zero_out_dst */ false, rhs, solution);
4453}
4454
4455
4456
4457template <typename MatrixType, typename VectorType, typename PreconditionerType>
4458inline void
4459PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tstep(
4460 VectorType &solution,
4461 const VectorType &rhs) const
4462{
4463 apply_internal<true>(/* zero_out_dst */ false, rhs, solution);
4464}
4465
4466
4467
4468template <typename MatrixType, typename VectorType, typename PreconditionerType>
4469inline double
4470PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
4471 vmult_with_last_residual_norm(VectorType &solution,
4472 const VectorType &rhs) const
4473{
4474 return apply_internal<false>(/* zero_out_dst */ true,
4475 rhs,
4476 solution,
4477 /* compute_norm */ true);
4478}
4479
4480
4481
4482template <typename MatrixType, typename VectorType, typename PreconditionerType>
4483inline double
4484PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
4485 step_with_last_residual_norm(VectorType &solution,
4486 const VectorType &rhs) const
4487{
4488 return apply_internal<false>(/* zero_out_dst */ false,
4489 rhs,
4490 solution,
4491 /* compute_norm */ true);
4492}
4493
4494
4495
4496template <typename MatrixType, typename VectorType, typename PreconditionerType>
4497inline typename PreconditionChebyshev<MatrixType,
4498 VectorType,
4499 PreconditionerType>::size_type
4500PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::m() const
4501{
4502 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4503 return matrix_ptr->m();
4504}
4505
4506
4507
4508template <typename MatrixType, typename VectorType, typename PreconditionerType>
4509inline typename PreconditionChebyshev<MatrixType,
4510 VectorType,
4511 PreconditionerType>::size_type
4512PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::n() const
4513{
4514 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4515 return matrix_ptr->n();
4516}
4517
4518
4519
4520template <typename MatrixType, typename VectorType, typename PreconditionerType>
4521inline void
4522PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::set_degree(
4523 const unsigned int degree)
4524{
4525 Assert(data.degree > 0,
4526 ExcMessage("The degree of the Chebyshev method must be positive."));
4527 data.degree = degree;
4528}
4529
4530#endif // DOXYGEN
4531
4532DEAL_II_NAMESPACE_CLOSE
4533
4534#endif
end
*  iterator end()
Definition ALE_Finite_Strain_Plasticity.h:10398
begin
*  *  iterator begin()
Definition ALE_Finite_Strain_Plasticity.h:10397
const_iterator
*  *  const_iterator()=default
Definition petsc_matrix_base.h:183
affine_constraints.h
block_vector_base.h
BlockVectorBase::n_blocks
unsigned int n_blocks() const
BlockVectorBase::block
BlockType & block(const unsigned int i)
IndexSet::is_empty
bool is_empty() const
Definition index_set.h:1909
IndexSet::nth_index_in_set
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1958
LinearAlgebra::distributed::Vector::locally_owned_size
size_type locally_owned_size() const
LinearAlgebra::distributed::Vector::swap
void swap(Vector< Number, MemorySpace > &v) noexcept
LinearAlgebra::distributed::Vector::locally_owned_elements
::IndexSet locally_owned_elements() const
LinearAlgebra::distributed::Vector::add_and_dot
Number add_and_dot(const Number a, const Vector< Number, MemorySpace > &V, const Vector< Number, MemorySpace > &W)
MultithreadInfo::n_threads
static unsigned int n_threads()
Definition multithread_info.cc:127
PreconditionChebyshev::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::m
size_type m() const
PreconditionChebyshev::n
size_type n() const
PreconditionChebyshev::step
void step(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::step_with_last_residual_norm
double step_with_last_residual_norm(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::estimate_eigenvalues
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
PreconditionChebyshev::Tstep
void Tstep(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::apply_internal
double apply_internal(const bool zero_out_dst, const VectorType &src, VectorType &dst, const bool compute_final_norm=false) const
PreconditionChebyshev::vmult_with_last_residual_norm
double vmult_with_last_residual_norm(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::set_degree
void set_degree(const unsigned int degree)
PreconditionChebyshev::vmult
void vmult(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::initialize
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
PreconditionChebyshev::matrix_ptr
ObserverPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
Definition precondition.h:2308
PreconditionIdentity::vmult_add
void vmult_add(VectorType &, const VectorType &) const
PreconditionIdentity::vmult
void vmult(VectorType &, const VectorType &) const
PreconditionIdentity::m
size_type m() const
PreconditionIdentity::initialize
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
PreconditionIdentity::n
size_type n() const
PreconditionIdentity::Tvmult
void Tvmult(VectorType &, const VectorType &) const
PreconditionIdentity::Tvmult_add
void Tvmult_add(VectorType &, const VectorType &) const
PreconditionJacobi::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1669
PreconditionJacobi::PreconditionerType
internal::PreconditionRelaxation::PreconditionJacobiImpl< MatrixType > PreconditionerType
Definition precondition.h:1662
PreconditionJacobi::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionPSOR::AdditionalData::AdditionalData
AdditionalData(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
PreconditionPSOR::AdditionalData::parameters
BaseClass::AdditionalData parameters
Definition precondition.h:1881
PreconditionPSOR::AdditionalData::inverse_permutation
const std::vector< size_type > & inverse_permutation
Definition precondition.h:1877
PreconditionPSOR::AdditionalData::permutation
const std::vector< size_type > & permutation
Definition precondition.h:1873
PreconditionPSOR::size_type
typename BaseClass::size_type size_type
Definition precondition.h:1847
PreconditionPSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &additional_data)
PreconditionPSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionPSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1840
PreconditionPSOR::initialize
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())
PreconditionRelaxation::AdditionalData::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, const double safety_factor=1.2)
PreconditionRelaxation::get_relaxation
double get_relaxation() const
PreconditionRelaxation::Tvmult
void Tvmult(VectorType &, const VectorType &) const
PreconditionRelaxation::preconditioner
std::shared_ptr< PreconditionerType > preconditioner
Definition precondition.h:706
PreconditionRelaxation::estimate_eigenvalues
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
PreconditionRelaxation::n
size_type n() const
PreconditionRelaxation::m
size_type m() const
PreconditionRelaxation::step
void step(VectorType &x, const VectorType &rhs) const
PreconditionRelaxation::A
ObserverPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
Definition precondition.h:695
PreconditionRelaxation::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionRelaxation::Tstep
void Tstep(VectorType &x, const VectorType &rhs) const
PreconditionRelaxation::vmult
void vmult(VectorType &, const VectorType &) const
PreconditionRelaxation::size_type
types::global_dof_index size_type
Definition precondition.h:580
PreconditionRichardson::AdditionalData::AdditionalData
AdditionalData(const double relaxation=1.)
PreconditionRichardson::size_type
types::global_dof_index size_type
Definition precondition.h:348
PreconditionRichardson::vmult_add
void vmult_add(VectorType &, const VectorType &) const
PreconditionRichardson::n
size_type n() const
PreconditionRichardson::initialize
void initialize(const AdditionalData &parameters)
PreconditionRichardson::initialize
void initialize(const MatrixType &matrix, const AdditionalData &parameters)
PreconditionRichardson::m
size_type m() const
PreconditionRichardson::vmult
void vmult(VectorType &, const VectorType &) const
PreconditionRichardson::Tvmult
void Tvmult(VectorType &, const VectorType &) const
PreconditionRichardson::Tvmult_add
void Tvmult_add(VectorType &, const VectorType &) const
PreconditionSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1732
PreconditionSOR::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1739
PreconditionSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionSSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionSSOR::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1791
PreconditionSSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionSSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1784
PreconditionUseMatrix::function_ptr
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition precondition.h:513
PreconditionUseMatrix::precondition
const function_ptr precondition
Definition precondition.h:538
PreconditionUseMatrix::matrix
const MatrixType & matrix
Definition precondition.h:533
PreconditionUseMatrix::vmult
void vmult(VectorType &dst, const VectorType &src) const
PreconditionUseMatrix::PreconditionUseMatrix
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
SparseMatrix::end
const_iterator end() const
SparseMatrix::begin
const_iterator begin() const
Vector::add
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
Vector::mean_value
Number mean_value() const
Vector::get_mpi_communicator
MPI_Comm get_mpi_communicator() const
Vector::size
virtual size_type size() const override
Vector::swap
virtual void swap(Vector< Number > &v) noexcept
Vector::begin
iterator begin()
VectorizedArray::store
void store(OtherNumber *ptr) const
Definition vectorization.h:592
VectorizedArray::load
void load(const OtherNumber *ptr)
Definition vectorization.h:579
DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:214
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
diagonal_matrix.h
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exception_macros.h:705
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
identity_matrix.h
vector_memory.h
memory_space.h
data
std::vector< index_type > data
Definition mpi.cc:734
size
std::size_t size
Definition mpi.cc:733
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:56
LAPACKSupport::A
constexpr char A
Definition lapack_support.h:151
LinearAlgebra::TpetraWrappers::TpetraTypes::VectorType
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Definition trilinos_tpetra_types.h:113
LinearAlgebra::TpetraWrappers::TpetraTypes::MatrixType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
Definition trilinos_tpetra_types.h:139
PlasticityLab::Number
*  *  *  RotationFunction< dim, Number >::RotationFunction Number(dim)
Definition ALE_Finite_Strain_Plasticity.h:8733
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm mpi_communicator)
internal::VectorImplementation::minimum_parallel_grain_size
unsigned int minimum_parallel_grain_size
Definition parallel.cc:48
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition matrix_block.h:595
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
std
STL namespace.
std::log
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6485
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6596
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6505
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6563
types::global_dof_index
unsigned int global_dof_index
Definition types.h:92
observer_pointer.h
solver_cg.h
PreconditionChebyshev::AdditionalData::AdditionalData
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, const double safety_factor=1.2)
internal::EigenvalueAlgorithmAdditionalData::preconditioner
std::shared_ptr< PreconditionerType > preconditioner
Definition precondition.h:203
internal::EigenvalueAlgorithmAdditionalData::EigenvalueAlgorithmAdditionalData
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, const double safety_factor=1.2)
internal::EigenvalueAlgorithmAdditionalData::operator=
EigenvalueAlgorithmAdditionalData< PreconditionerType > & operator=(const EigenvalueAlgorithmAdditionalData< PreconditionerType > &other_data)
internal::EigenvalueAlgorithmAdditionalData::constraints
::AffineConstraints< double > constraints
Definition precondition.h:185
eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
SymmetricTensorEigenvectorMethod::jacobi
@ jacobi
template_constraints.h
Generated by doxygen 1.9.8