deal.II version 9.8.0
\(\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
60
61namespace internal
62{
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
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
2208 void
2209 Tvmult(VectorType &dst, const VectorType &src) const;
2210
2214 void
2215 step(VectorType &dst, const VectorType &src) const;
2216
2220 void
2221 Tstep(VectorType &dst, const VectorType &src) const;
2222
2226 void
2227 clear();
2228
2233 size_type
2234 m() const;
2235
2240 size_type
2241 n() const;
2242
2243 using EigenvalueInformation = internal::EigenvalueInformation;
2244
2257 EigenvalueInformation
2258 estimate_eigenvalues(const VectorType &src) const;
2259
2260private:
2265 template <bool use_transpose>
2266 void
2267 apply_internal(const bool zero_out_dst,
2268 const VectorType &src,
2269 VectorType &dst) const;
2270
2274 ObserverPointer<
2275 const MatrixType,
2276 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>>
2277 matrix_ptr;
2278
2282 mutable VectorType solution_old;
2283
2287 mutable VectorType temp_vector1;
2288
2292 mutable VectorType temp_vector2;
2293
2298 AdditionalData data;
2299
2303 double theta;
2304
2309 double delta;
2310
2315 bool eigenvalues_are_initialized;
2316
2321 mutable Threads::Mutex mutex;
2322};
2323
2324
2325
2327/* ---------------------------------- Inline functions ------------------- */
2328
2329#ifndef DOXYGEN
2330
2331
2332namespace internal
2333{
2334 template <typename VectorType>
2335 void
2336 set_initial_guess(VectorType &vector)
2337 {
2338 vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2339 if (vector.locally_owned_elements().is_element(0))
2340 vector(0) = 0.;
2341 }
2342
2343
2344
2345 template <typename Number>
2346 void
2347 set_initial_guess(::Vector<Number> &vector)
2348 {
2349 // Choose a high-frequency mode consisting of numbers between 0 and 10
2350 // that is cheap to compute (cheaper than random numbers) but avoids
2351 // obviously re-occurring numbers in multi-component systems by choosing
2352 // a period of 11
2353 for (unsigned int i = 0; i < vector.size(); ++i)
2354 vector(i) = i % 11;
2355
2356 // Then make sure the vector has mean value zero:
2357 const Number mean_value = vector.mean_value();
2358 vector.add(-mean_value);
2359 }
2360
2361
2362
2363 template <typename Number, typename MemorySpace>
2364 void
2365 set_initial_guess(
2366 ::LinearAlgebra::distributed::BlockVector<Number, MemorySpace>
2367 &vector)
2368 {
2369 for (unsigned int block = 0; block < vector.n_blocks(); ++block)
2370 set_initial_guess(vector.block(block));
2371 }
2372
2373
2374
2375 template <typename Number, typename MemorySpace>
2376 void
2377 set_initial_guess(
2378 ::LinearAlgebra::distributed::Vector<Number, MemorySpace> &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 // Make initial guess robust with respect to number of processors
2385 // by operating on the global index.
2386 types::global_dof_index first_local_range = 0;
2387 if (!vector.locally_owned_elements().is_empty())
2388 first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2389
2390 const auto n_local_elements = vector.locally_owned_size();
2391 Number *values_ptr = vector.get_values();
2392 typename MemorySpace::kokkos_space::execution_space exec;
2393 Kokkos::RangePolicy<typename MemorySpace::kokkos_space::execution_space,
2394 Kokkos::IndexType<types::global_dof_index>>
2395 policy(exec, 0, n_local_elements);
2396
2397 Kokkos::parallel_for(
2398 "::PreconditionChebyshev::set_initial_guess",
2399 policy,
2400 KOKKOS_LAMBDA(types::global_dof_index i) {
2401 values_ptr[i] = (i + first_local_range) % 11;
2402 });
2403 exec.fence();
2404
2405 // Then make sure the vector has mean value zero:
2406 const Number mean_value = vector.mean_value();
2407 vector.add(-mean_value);
2408 }
2409
2410
2411
2412 struct EigenvalueTracker
2413 {
2414 public:
2415 void
2416 slot(const std::vector<double> &eigenvalues)
2417 {
2418 values = eigenvalues;
2419 }
2420
2421 std::vector<double> values;
2422 };
2423
2424
2425
2426 template <typename MatrixType,
2427 typename VectorType,
2428 typename PreconditionerType>
2429 double
2430 power_iteration(const MatrixType &matrix,
2431 VectorType &eigenvector,
2432 const PreconditionerType &preconditioner,
2433 const unsigned int n_iterations)
2434 {
2435 typename VectorType::value_type eigenvalue_estimate = 0.;
2436 eigenvector /= eigenvector.l2_norm();
2437 VectorType vector1, vector2;
2438 vector1.reinit(eigenvector, true);
2439 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2440 vector2.reinit(eigenvector, true);
2441 for (unsigned int i = 0; i < n_iterations; ++i)
2442 {
2443 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2444 {
2445 matrix.vmult(vector2, eigenvector);
2446 preconditioner.vmult(vector1, vector2);
2447 }
2448 else
2449 matrix.vmult(vector1, eigenvector);
2450
2451 eigenvalue_estimate = eigenvector * vector1;
2452
2453 vector1 /= vector1.l2_norm();
2454 eigenvector.swap(vector1);
2455 }
2456 return std::abs(eigenvalue_estimate);
2457 }
2458
2459
2460
2461 template <typename MatrixType,
2462 typename VectorType,
2463 typename PreconditionerType>
2464 EigenvalueInformation
2465 estimate_eigenvalues(
2466 const EigenvalueAlgorithmAdditionalData<PreconditionerType> &data,
2467 const MatrixType *matrix_ptr,
2468 VectorType &solution_old,
2469 VectorType &temp_vector1,
2470 const unsigned int degree,
2471 const double safety_factor)
2472 {
2473 Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2474
2475 Assert(
2476 safety_factor >= 1.,
2477 ExcMessage(
2478 "The safety factor must be at least 1.0 to ensure the maximum eigenvalue is "
2479 "not underestimated. Please set safety_factor >= 1.0 in your AdditionalData."));
2480
2481 EigenvalueInformation info{};
2482
2483 if (data.eig_cg_n_iterations > 0)
2484 {
2485 Assert(data.eig_cg_n_iterations > 2,
2486 ExcMessage(
2487 "Need to set at least two iterations to find eigenvalues."));
2488
2489 internal::EigenvalueTracker eigenvalue_tracker;
2490
2491 // set an initial guess that contains some high-frequency parts (to the
2492 // extent possible without knowing the discretization and the numbering)
2493 // to trigger high eigenvalues according to the external function
2494 internal::set_initial_guess(temp_vector1);
2495 data.constraints.set_zero(temp_vector1);
2496
2497 // We want to run the power iteration starting with temp_vector1, but
2498 // for that it needs to be a nonzero vector. However, there are
2499 // situations where it only contains zeroes (e.g., if all dofs are
2500 // constrained). If that is the case, there is nothing we can do
2501 // because in both of the two methods below, we either solve a linear
2502 // system with temp_vector1 as the right hand side, or use it as a
2503 // starting guess for a power iteration, but in neither case a
2504 // zero vector yields anything useful. So skip the following code and
2505 // instead use the fallback below if temp_vector1 is the zero vector.
2506 if (temp_vector1.all_zero() == false)
2507 {
2508 if (data.eigenvalue_algorithm ==
2509 internal::EigenvalueAlgorithm::lanczos)
2510 {
2511 // set a very strict tolerance to force at least two iterations
2512 IterationNumberControl control(data.eig_cg_n_iterations,
2513 1e-10,
2514 false,
2515 false);
2516
2517 ::SolverCG<VectorType> solver(control);
2518 solver.connect_eigenvalues_slot(
2519 [&eigenvalue_tracker](
2520 const std::vector<double> &eigenvalues) {
2521 eigenvalue_tracker.slot(eigenvalues);
2522 });
2523
2524 solver.solve(*matrix_ptr,
2525 solution_old,
2526 temp_vector1,
2527 *data.preconditioner);
2528
2529 info.cg_iterations = control.last_step();
2530 }
2531 else if (data.eigenvalue_algorithm ==
2532 internal::EigenvalueAlgorithm::power_iteration)
2533 {
2534 (void)degree;
2535
2536 Assert(degree != numbers::invalid_unsigned_int,
2537 ExcMessage(
2538 "Cannot estimate the minimal eigenvalue with the "
2539 "power iteration"));
2540
2541 eigenvalue_tracker.values.push_back(
2542 internal::power_iteration(*matrix_ptr,
2543 temp_vector1,
2544 *data.preconditioner,
2545 data.eig_cg_n_iterations));
2546 }
2547 else
2548 DEAL_II_NOT_IMPLEMENTED();
2549 }
2550
2551 // read the eigenvalues from the attached eigenvalue tracker
2552 if (eigenvalue_tracker.values.empty())
2553 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2554 else
2555 {
2556 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2557
2558 // include a safety factor since the CG method will in general not
2559 // be converged
2560 info.max_eigenvalue_estimate =
2561 safety_factor * eigenvalue_tracker.values.back();
2562 }
2563 }
2564 else
2565 {
2566 info.max_eigenvalue_estimate = data.max_eigenvalue;
2567 info.min_eigenvalue_estimate =
2568 data.max_eigenvalue / data.smoothing_range;
2569 }
2570
2571 return info;
2572 }
2573} // namespace internal
2574
2575
2576
2577inline PreconditionIdentity::PreconditionIdentity()
2578 : n_rows(0)
2579 , n_columns(0)
2580{}
2581
2582template <typename MatrixType>
2583inline void
2584PreconditionIdentity::initialize(const MatrixType &matrix,
2585 const PreconditionIdentity::AdditionalData &)
2586{
2587 n_rows = matrix.m();
2588 n_columns = matrix.n();
2589}
2590
2591
2592template <typename VectorType>
2593inline void
2594PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
2595{
2596 dst = src;
2597}
2598
2599
2600
2601template <typename VectorType>
2602inline void
2603PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
2604{
2605 dst = src;
2606}
2607
2608template <typename VectorType>
2609inline void
2610PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
2611{
2612 dst += src;
2613}
2614
2615
2616
2617template <typename VectorType>
2618inline void
2619PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
2620{
2621 dst += src;
2622}
2623
2624
2625
2626inline void
2627PreconditionIdentity::clear()
2628{}
2629
2630
2631
2632inline PreconditionIdentity::size_type
2633PreconditionIdentity::m() const
2634{
2635 Assert(n_rows != 0, ExcNotInitialized());
2636 return n_rows;
2637}
2638
2639inline PreconditionIdentity::size_type
2640PreconditionIdentity::n() const
2641{
2642 Assert(n_columns != 0, ExcNotInitialized());
2643 return n_columns;
2644}
2645
2646//---------------------------------------------------------------------------
2647
2648inline PreconditionRichardson::AdditionalData::AdditionalData(
2649 const double relaxation)
2650 : relaxation(relaxation)
2651{}
2652
2653
2654inline PreconditionRichardson::PreconditionRichardson()
2655 : relaxation(0)
2656 , n_rows(0)
2657 , n_columns(0)
2658{
2659 AdditionalData add_data;
2660 relaxation = add_data.relaxation;
2661}
2662
2663
2664
2665inline void
2666PreconditionRichardson::initialize(
2667 const PreconditionRichardson::AdditionalData &parameters)
2668{
2669 relaxation = parameters.relaxation;
2670}
2671
2672
2673
2674template <typename MatrixType>
2675inline void
2676PreconditionRichardson::initialize(
2677 const MatrixType &matrix,
2678 const PreconditionRichardson::AdditionalData &parameters)
2679{
2680 relaxation = parameters.relaxation;
2681 n_rows = matrix.m();
2682 n_columns = matrix.n();
2683}
2684
2685
2686
2687template <typename VectorType>
2688inline void
2689PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2690{
2691 static_assert(
2692 std::is_same_v<size_type, typename VectorType::size_type>,
2693 "PreconditionRichardson and VectorType must have the same size_type.");
2694
2695 dst.equ(relaxation, src);
2696}
2697
2698
2699
2700template <typename VectorType>
2701inline void
2702PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2703{
2704 static_assert(
2705 std::is_same_v<size_type, typename VectorType::size_type>,
2706 "PreconditionRichardson and VectorType must have the same size_type.");
2707
2708 dst.equ(relaxation, src);
2709}
2710
2711template <typename VectorType>
2712inline void
2713PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2714{
2715 static_assert(
2716 std::is_same_v<size_type, typename VectorType::size_type>,
2717 "PreconditionRichardson and VectorType must have the same size_type.");
2718
2719 dst.add(relaxation, src);
2720}
2721
2722
2723
2724template <typename VectorType>
2725inline void
2726PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2727{
2728 static_assert(
2729 std::is_same_v<size_type, typename VectorType::size_type>,
2730 "PreconditionRichardson and VectorType must have the same size_type.");
2731
2732 dst.add(relaxation, src);
2733}
2734
2735inline PreconditionRichardson::size_type
2736PreconditionRichardson::m() const
2737{
2738 Assert(n_rows != 0, ExcNotInitialized());
2739 return n_rows;
2740}
2741
2742inline PreconditionRichardson::size_type
2743PreconditionRichardson::n() const
2744{
2745 Assert(n_columns != 0, ExcNotInitialized());
2746 return n_columns;
2747}
2748
2749//---------------------------------------------------------------------------
2750
2751template <typename MatrixType, typename PreconditionerType>
2752inline void
2753PreconditionRelaxation<MatrixType, PreconditionerType>::initialize(
2754 const MatrixType &rA,
2755 const AdditionalData &parameters)
2756{
2757 A = &rA;
2758 eigenvalues_are_initialized = false;
2759
2760 Assert(parameters.preconditioner, ExcNotInitialized());
2761
2762 this->data = parameters;
2763}
2764
2765
2766template <typename MatrixType, typename PreconditionerType>
2767inline void
2768PreconditionRelaxation<MatrixType, PreconditionerType>::clear()
2769{
2770 eigenvalues_are_initialized = false;
2771 A = nullptr;
2772 data.relaxation = 1.0;
2773 data.preconditioner = nullptr;
2774}
2775
2776template <typename MatrixType, typename PreconditionerType>
2777inline
2778 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2779 PreconditionRelaxation<MatrixType, PreconditionerType>::m() const
2780{
2781 Assert(A != nullptr, ExcNotInitialized());
2782 return A->m();
2783}
2784
2785template <typename MatrixType, typename PreconditionerType>
2786inline
2787 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2788 PreconditionRelaxation<MatrixType, PreconditionerType>::n() const
2789{
2790 Assert(A != nullptr, ExcNotInitialized());
2791 return A->n();
2792}
2793
2794template <typename MatrixType, typename PreconditionerType>
2795template <typename VectorType>
2796inline void
2797PreconditionRelaxation<MatrixType, PreconditionerType>::vmult(
2798 VectorType &dst,
2799 const VectorType &src) const
2800{
2801 Assert(this->A != nullptr, ExcNotInitialized());
2802 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2803
2804 if (eigenvalues_are_initialized == false)
2805 estimate_eigenvalues(src);
2806
2807 VectorType tmp1, tmp2;
2808
2809 for (unsigned int i = 0; i < data.n_iterations; ++i)
2810 internal::PreconditionRelaxation::step_operations(*A,
2811 *data.preconditioner,
2812 dst,
2813 src,
2814 data.relaxation,
2815 tmp1,
2816 tmp2,
2817 i,
2818 false);
2819}
2820
2821template <typename MatrixType, typename PreconditionerType>
2822template <typename VectorType>
2823inline void
2824PreconditionRelaxation<MatrixType, PreconditionerType>::Tvmult(
2825 VectorType &dst,
2826 const VectorType &src) const
2827{
2828 Assert(this->A != nullptr, ExcNotInitialized());
2829 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2830
2831 if (eigenvalues_are_initialized == false)
2832 estimate_eigenvalues(src);
2833
2834 VectorType tmp1, tmp2;
2835
2836 for (unsigned int i = 0; i < data.n_iterations; ++i)
2837 internal::PreconditionRelaxation::step_operations(
2838 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2839}
2840
2841template <typename MatrixType, typename PreconditionerType>
2842template <typename VectorType>
2843inline void
2844PreconditionRelaxation<MatrixType, PreconditionerType>::step(
2845 VectorType &dst,
2846 const VectorType &src) const
2847{
2848 Assert(this->A != nullptr, ExcNotInitialized());
2849 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2850
2851 if (eigenvalues_are_initialized == false)
2852 estimate_eigenvalues(src);
2853
2854 VectorType tmp1, tmp2;
2855
2856 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2857 internal::PreconditionRelaxation::step_operations(*A,
2858 *data.preconditioner,
2859 dst,
2860 src,
2861 data.relaxation,
2862 tmp1,
2863 tmp2,
2864 i,
2865 false);
2866}
2867
2868template <typename MatrixType, typename PreconditionerType>
2869template <typename VectorType>
2870inline void
2871PreconditionRelaxation<MatrixType, PreconditionerType>::Tstep(
2872 VectorType &dst,
2873 const VectorType &src) const
2874{
2875 Assert(this->A != nullptr, ExcNotInitialized());
2876 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2877
2878 if (eigenvalues_are_initialized == false)
2879 estimate_eigenvalues(src);
2880
2881 VectorType tmp1, tmp2;
2882
2883 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2884 internal::PreconditionRelaxation::step_operations(
2885 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2886}
2887
2888template <typename MatrixType, typename PreconditionerType>
2889template <typename VectorType>
2890inline internal::EigenvalueInformation
2891PreconditionRelaxation<MatrixType, PreconditionerType>::estimate_eigenvalues(
2892 const VectorType &src) const
2893{
2894 Assert(eigenvalues_are_initialized == false, ExcInternalError());
2895
2896 EigenvalueInformation info;
2897
2898 if (data.relaxation == 0.0)
2899 {
2900 VectorType solution_old, temp_vector1;
2901
2902 solution_old.reinit(src);
2903 temp_vector1.reinit(src, true);
2904
2905 info = internal::estimate_eigenvalues<MatrixType>(data,
2906 A,
2907 solution_old,
2908 temp_vector1,
2909 data.n_iterations,
2910 data.safety_factor);
2911
2912 const double alpha =
2913 (data.smoothing_range > 1. ?
2914 info.max_eigenvalue_estimate / data.smoothing_range :
2915 std::min(0.9 * info.max_eigenvalue_estimate,
2916 info.min_eigenvalue_estimate));
2917
2918 const_cast<PreconditionRelaxation<MatrixType, PreconditionerType> *>(this)
2919 ->data.relaxation = 2.0 / (alpha + info.max_eigenvalue_estimate);
2920 }
2921
2922 const_cast<PreconditionRelaxation<MatrixType, PreconditionerType> *>(this)
2923 ->eigenvalues_are_initialized = true;
2924
2925 return info;
2926}
2927
2928template <typename MatrixType, typename PreconditionerType>
2929double
2930PreconditionRelaxation<MatrixType, PreconditionerType>::get_relaxation() const
2931{
2932 return data.relaxation;
2933}
2934
2935
2936//---------------------------------------------------------------------------
2937
2938template <typename MatrixType>
2939inline void
2940PreconditionJacobi<MatrixType>::initialize(const MatrixType &A,
2941 const AdditionalData &parameters_in)
2942{
2943 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2944 Assert(
2945 parameters_in.relaxation != 0.0,
2946 ExcMessage(
2947 "Relaxation cannot automatically be determined by PreconditionJacobi."));
2948
2949 AdditionalData parameters;
2950 parameters.relaxation = 1.0;
2951 parameters.n_iterations = parameters_in.n_iterations;
2952 parameters.preconditioner =
2953 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2954
2955 this->BaseClass::initialize(A, parameters);
2956}
2957
2958//---------------------------------------------------------------------------
2959
2960template <typename MatrixType>
2961inline void
2962PreconditionSOR<MatrixType>::initialize(const MatrixType &A,
2963 const AdditionalData &parameters_in)
2964{
2965 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2966 Assert(
2967 parameters_in.relaxation != 0.0,
2968 ExcMessage(
2969 "Relaxation cannot automatically be determined by PreconditionSOR."));
2970
2971 AdditionalData parameters;
2972 parameters.relaxation = 1.0;
2973 parameters.n_iterations = parameters_in.n_iterations;
2974 parameters.preconditioner =
2975 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2976
2977 this->BaseClass::initialize(A, parameters);
2978}
2979
2980//---------------------------------------------------------------------------
2981
2982template <typename MatrixType>
2983inline void
2984PreconditionSSOR<MatrixType>::initialize(const MatrixType &A,
2985 const AdditionalData &parameters_in)
2986{
2987 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2988 Assert(
2989 parameters_in.relaxation != 0.0,
2990 ExcMessage(
2991 "Relaxation cannot automatically be determined by PreconditionSSOR."));
2992
2993 AdditionalData parameters;
2994 parameters.relaxation = 1.0;
2995 parameters.n_iterations = parameters_in.n_iterations;
2996 parameters.preconditioner =
2997 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2998
2999 this->BaseClass::initialize(A, parameters);
3000}
3001
3002
3003
3004//---------------------------------------------------------------------------
3005
3006template <typename MatrixType>
3007inline void
3008PreconditionPSOR<MatrixType>::initialize(
3009 const MatrixType &A,
3010 const std::vector<size_type> &p,
3011 const std::vector<size_type> &ip,
3012 const typename BaseClass::AdditionalData &parameters_in)
3013{
3014 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
3015 Assert(
3016 parameters_in.relaxation != 0.0,
3017 ExcMessage(
3018 "Relaxation cannot automatically be determined by PreconditionPSOR."));
3019
3020 typename BaseClass::AdditionalData parameters;
3021 parameters.relaxation = 1.0;
3022 parameters.n_iterations = parameters_in.n_iterations;
3023 parameters.preconditioner =
3024 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
3025
3026 this->BaseClass::initialize(A, parameters);
3027}
3028
3029
3030template <typename MatrixType>
3031inline void
3032PreconditionPSOR<MatrixType>::initialize(const MatrixType &A,
3033 const AdditionalData &additional_data)
3034{
3035 initialize(A,
3036 additional_data.permutation,
3037 additional_data.inverse_permutation,
3038 additional_data.parameters);
3039}
3040
3041template <typename MatrixType>
3042PreconditionPSOR<MatrixType>::AdditionalData::AdditionalData(
3043 const std::vector<size_type> &permutation,
3044 const std::vector<size_type> &inverse_permutation,
3045 const typename PreconditionPSOR<MatrixType>::BaseClass::AdditionalData
3046 &parameters)
3047 : permutation(permutation)
3048 , inverse_permutation(inverse_permutation)
3049 , parameters(parameters)
3050{}
3051
3052
3053//---------------------------------------------------------------------------
3054
3055
3056template <typename MatrixType, typename VectorType>
3057PreconditionUseMatrix<MatrixType, VectorType>::PreconditionUseMatrix(
3058 const MatrixType &M,
3059 const function_ptr method)
3060 : matrix(M)
3061 , precondition(method)
3062{}
3063
3064
3065
3066template <typename MatrixType, typename VectorType>
3067void
3068PreconditionUseMatrix<MatrixType, VectorType>::vmult(
3069 VectorType &dst,
3070 const VectorType &src) const
3071{
3072 (matrix.*precondition)(dst, src);
3073}
3074
3075//---------------------------------------------------------------------------
3076
3077namespace internal
3078{
3079
3080 template <typename PreconditionerType>
3081 inline EigenvalueAlgorithmAdditionalData<PreconditionerType>::
3082 EigenvalueAlgorithmAdditionalData(
3083 const double smoothing_range,
3084 const unsigned int eig_cg_n_iterations,
3085 const double eig_cg_residual,
3086 const double max_eigenvalue,
3087 const EigenvalueAlgorithm eigenvalue_algorithm,
3088 const double safety_factor)
3089 : smoothing_range(smoothing_range)
3090 , eig_cg_n_iterations(eig_cg_n_iterations)
3091 , eig_cg_residual(eig_cg_residual)
3092 , max_eigenvalue(max_eigenvalue)
3093 , eigenvalue_algorithm(eigenvalue_algorithm)
3094 , safety_factor(safety_factor)
3095 {}
3096
3097
3098
3099 template <typename PreconditionerType>
3100 inline EigenvalueAlgorithmAdditionalData<PreconditionerType> &
3101 EigenvalueAlgorithmAdditionalData<PreconditionerType>::operator=(
3102 const EigenvalueAlgorithmAdditionalData &other_data)
3103 {
3104 smoothing_range = other_data.smoothing_range;
3105 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3106 eig_cg_residual = other_data.eig_cg_residual;
3107 max_eigenvalue = other_data.max_eigenvalue;
3108 preconditioner = other_data.preconditioner;
3109 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3110 safety_factor = other_data.safety_factor;
3111 constraints.copy_from(other_data.constraints);
3112
3113 return *this;
3114 }
3115} // namespace internal
3116
3117template <typename MatrixType, typename PreconditionerType>
3118inline PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData::
3119 AdditionalData(const double relaxation,
3120 const unsigned int n_iterations,
3121 const double smoothing_range,
3122 const unsigned int eig_cg_n_iterations,
3123 const double eig_cg_residual,
3124 const double max_eigenvalue,
3125 const EigenvalueAlgorithm eigenvalue_algorithm,
3126 const double safety_factor)
3127 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3128 smoothing_range,
3129 eig_cg_n_iterations,
3130 eig_cg_residual,
3131 max_eigenvalue,
3132 eigenvalue_algorithm,
3133 safety_factor)
3134 , relaxation(relaxation)
3135 , n_iterations(n_iterations)
3136{}
3137
3138
3139
3140//---------------------------------------------------------------------------
3141
3142namespace internal
3143{
3144 namespace PreconditionChebyshevImplementation
3145 {
3146 // for deal.II vectors, perform updates for Chebyshev preconditioner all
3147 // at once to reduce memory transfer. Here, we select between general
3148 // vectors and deal.II vectors where we expand the loop over the (local)
3149 // size of the vector
3150
3151 // generic part for non-deal.II vectors
3152 template <typename VectorType, typename PreconditionerType>
3153 inline void
3154 vector_updates(const VectorType &rhs,
3155 const PreconditionerType &preconditioner,
3156 const unsigned int iteration_index,
3157 const double factor1,
3158 const double factor2,
3159 VectorType &solution_old,
3160 VectorType &temp_vector1,
3161 VectorType &temp_vector2,
3162 VectorType &solution)
3163 {
3164 if (iteration_index == 0)
3165 {
3166 solution.equ(factor2, rhs);
3167 preconditioner.vmult(solution_old, solution);
3168 }
3169 else if (iteration_index == 1)
3170 {
3171 // compute t = P^{-1} * (b-A*x^{n})
3172 temp_vector1.sadd(-1.0, 1.0, rhs);
3173 preconditioner.vmult(solution_old, temp_vector1);
3174
3175 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
3176 solution_old.sadd(factor2, 1 + factor1, solution);
3177 }
3178 else
3179 {
3180 // compute t = P^{-1} * (b-A*x^{n})
3181 temp_vector1.sadd(-1.0, 1.0, rhs);
3182 preconditioner.vmult(temp_vector2, temp_vector1);
3183
3184 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
3185 solution_old.sadd(-factor1, factor2, temp_vector2);
3186 solution_old.add(1 + factor1, solution);
3187 }
3188
3189 solution.swap(solution_old);
3190 }
3191
3192 // generic part for deal.II vectors
3193 template <
3194 typename Number,
3195 typename PreconditionerType,
3196 std::enable_if_t<
3197 !has_vmult_with_std_functions_for_precondition<
3198 PreconditionerType,
3199 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>>,
3200 int> * = nullptr>
3201 inline void
3202 vector_updates(
3203 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3204 const PreconditionerType &preconditioner,
3205 const unsigned int iteration_index,
3206 const double factor1_,
3207 const double factor2_,
3208 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3209 &solution_old,
3210 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3211 &temp_vector1,
3212 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3213 &temp_vector2,
3214 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
3215 {
3216 const Number factor1 = factor1_;
3217 const Number factor1_plus_1 = 1. + factor1_;
3218 const Number factor2 = factor2_;
3219
3220 if (iteration_index == 0)
3221 {
3222 // compute t = P^{-1} * (b)
3223 preconditioner.vmult(solution_old, rhs);
3224
3225 // compute x^{n+1} = f_2 * t
3226 const auto solution_old_ptr = solution_old.begin();
3227 DEAL_II_OPENMP_SIMD_PRAGMA
3228 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3229 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
3230 }
3231 else if (iteration_index == 1)
3232 {
3233 // compute t = P^{-1} * (b-A*x^{n})
3234 temp_vector1.sadd(-1.0, 1.0, rhs);
3235
3236 preconditioner.vmult(solution_old, temp_vector1);
3237
3238 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
3239 const auto solution_ptr = solution.begin();
3240 const auto solution_old_ptr = solution_old.begin();
3241 DEAL_II_OPENMP_SIMD_PRAGMA
3242 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3243 solution_old_ptr[i] =
3244 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
3245 }
3246 else
3247 {
3248 // compute t = P^{-1} * (b-A*x^{n})
3249 temp_vector1.sadd(-1.0, 1.0, rhs);
3250
3251 preconditioner.vmult(temp_vector2, temp_vector1);
3252
3253 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
3254 const auto solution_ptr = solution.begin();
3255 const auto solution_old_ptr = solution_old.begin();
3256 const auto temp_vector2_ptr = temp_vector2.begin();
3257 DEAL_II_OPENMP_SIMD_PRAGMA
3258 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3259 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3260 factor1 * solution_old_ptr[i] +
3261 temp_vector2_ptr[i] * factor2;
3262 }
3263
3264 solution.swap(solution_old);
3265 }
3266
3267 template <
3268 typename Number,
3269 typename PreconditionerType,
3270 std::enable_if_t<
3271 has_vmult_with_std_functions_for_precondition<
3272 PreconditionerType,
3273 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>>,
3274 int> * = nullptr>
3275 inline void
3276 vector_updates(
3277 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3278 const PreconditionerType &preconditioner,
3279 const unsigned int iteration_index,
3280 const double factor1_,
3281 const double factor2_,
3282 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3283 &solution_old,
3284 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3285 &temp_vector1,
3286 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3287 &temp_vector2,
3288 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
3289 {
3290 const Number factor1 = factor1_;
3291 const Number factor1_plus_1 = 1. + factor1_;
3292 const Number factor2 = factor2_;
3293
3294 const auto rhs_ptr = rhs.begin();
3295 const auto temp_vector1_ptr = temp_vector1.begin();
3296 const auto temp_vector2_ptr = temp_vector2.begin();
3297 const auto solution_ptr = solution.begin();
3298 const auto solution_old_ptr = solution_old.begin();
3299
3300 if (iteration_index == 0)
3301 {
3302 preconditioner.vmult(
3303 temp_vector2,
3304 rhs,
3305 [&](const auto start_range, const auto end_range) {
3306 if (end_range > start_range)
3307 std::memset(temp_vector2_ptr + start_range,
3308 0,
3309 sizeof(Number) * (end_range - start_range));
3310 },
3311 [&](const auto begin, const auto end) {
3312 DEAL_II_OPENMP_SIMD_PRAGMA
3313 for (std::size_t i = begin; i < end; ++i)
3314 temp_vector2_ptr[i] *= factor2;
3315 });
3316 }
3317 else
3318 {
3319 preconditioner.vmult(
3320 temp_vector2,
3321 temp_vector1,
3322 [&](const auto begin, const auto end) {
3323 if (end > begin)
3324 std::memset(temp_vector2_ptr + begin,
3325 0,
3326 sizeof(Number) * (end - begin));
3327
3328 DEAL_II_OPENMP_SIMD_PRAGMA
3329 for (std::size_t i = begin; i < end; ++i)
3330 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3331 },
3332 [&](const auto begin, const auto end) {
3333 if (iteration_index == 1)
3334 {
3335 DEAL_II_OPENMP_SIMD_PRAGMA
3336 for (std::size_t i = begin; i < end; ++i)
3337 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
3338 factor2 * temp_vector2_ptr[i];
3339 }
3340 else
3341 {
3342 DEAL_II_OPENMP_SIMD_PRAGMA
3343 for (std::size_t i = begin; i < end; ++i)
3344 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3345 factor1 * solution_old_ptr[i] +
3346 factor2 * temp_vector2_ptr[i];
3347 }
3348 });
3349 }
3350
3351 solution_old.swap(temp_vector2);
3352 solution_old.swap(solution);
3353 }
3354
3355 // vector updates for device vectors and DiagonalMatrix as preconditioner
3356 template <typename Number>
3357 inline void
3358 vector_updates(
3359 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>
3360 &rhs,
3361 const ::DiagonalMatrix<
3362 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>>
3363 &preconditioner,
3364 const unsigned int iteration_index,
3365 const double factor1_,
3366 const double factor2_,
3367 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>
3368 &solution_old,
3369 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>
3370 &temp_vector1,
3371 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default> &,
3372 LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>
3373 &solution)
3374 {
3375 const Number factor1 = factor1_;
3376 const Number factor1_plus_1 = 1. + factor1_;
3377 const Number factor2 = factor2_;
3378
3379 auto exec = typename ::MemorySpace::Default::kokkos_space::
3380 execution_space{};
3381 Number *sol_data = solution.begin();
3382 const Number *rhs_data = rhs.begin();
3383 const Number *prec_data = preconditioner.get_vector().begin();
3384 Number *sol_old_data = solution_old.begin();
3385 Number *tmp_data = temp_vector1.begin();
3386 using size_type = typename LinearAlgebra::distributed::
3387 Vector<Number, MemorySpace::Default>::size_type;
3388 const size_type size = solution.locally_owned_size();
3389 if (iteration_index == 0)
3390 {
3391 Kokkos::parallel_for(
3392 "::ChebyshevIt0",
3393 Kokkos::RangePolicy<
3394 ::MemorySpace::Default::kokkos_space::execution_space>(
3395 exec, 0, size),
3396 KOKKOS_LAMBDA(size_type i) {
3397 sol_old_data[i] = prec_data[i] * factor2 * rhs_data[i];
3398 });
3399 }
3400 else if (iteration_index == 1)
3401 {
3402 Kokkos::parallel_for(
3403 "::ChebyshevIt1",
3404 Kokkos::RangePolicy<
3405 ::MemorySpace::Default::kokkos_space::execution_space>(
3406 exec, 0, size),
3407 KOKKOS_LAMBDA(size_type i) {
3408 const Number residual = rhs_data[i] - tmp_data[i];
3409 sol_old_data[i] = factor2 * prec_data[i] * residual +
3410 factor1_plus_1 * sol_data[i];
3411 });
3412 }
3413 else
3414 {
3415 Kokkos::parallel_for(
3416 "::ChebyshevItk",
3417 Kokkos::RangePolicy<
3418 ::MemorySpace::Default::kokkos_space::execution_space>(
3419 exec, 0, size),
3420 KOKKOS_LAMBDA(size_type i) {
3421 const Number residual = rhs_data[i] - tmp_data[i];
3422 sol_old_data[i] = factor2 * prec_data[i] * residual -
3423 factor1 * sol_old_data[i] +
3424 factor1_plus_1 * sol_data[i];
3425 });
3426 }
3427 exec.fence();
3428
3429 solution.swap(solution_old);
3430 }
3431
3432 // worker routine for deal.II vectors. Because of vectorization, we need
3433 // to put the loop into an extra structure because the virtual function of
3434 // VectorUpdatesRange prevents the compiler from applying vectorization.
3435 template <typename Number>
3436 struct VectorUpdater
3437 {
3438 VectorUpdater(const Number *rhs,
3439 const Number *matrix_diagonal_inverse,
3440 const unsigned int iteration_index,
3441 const Number factor1,
3442 const Number factor2,
3443 Number *solution_old,
3444 Number *tmp_vector,
3445 Number *solution)
3446 : rhs(rhs)
3447 , matrix_diagonal_inverse(matrix_diagonal_inverse)
3448 , iteration_index(iteration_index)
3449 , factor1(factor1)
3450 , factor2(factor2)
3451 , solution_old(solution_old)
3452 , tmp_vector(tmp_vector)
3453 , solution(solution)
3454 {}
3455
3456 void
3457 apply_to_subrange(const std::size_t begin, const std::size_t end) const
3458 {
3459 // To circumvent a bug in gcc
3460 // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
3461 // copies of the variables factor1 and factor2 and do not check based on
3462 // factor1.
3463 const Number factor1 = this->factor1;
3464 const Number factor1_plus_1 = 1. + this->factor1;
3465 const Number factor2 = this->factor2;
3466 if (iteration_index == 0)
3467 {
3468 DEAL_II_OPENMP_SIMD_PRAGMA
3469 for (std::size_t i = begin; i < end; ++i)
3470 tmp_vector[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
3471 }
3472 else if (iteration_index == 1)
3473 {
3474 // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
3475 DEAL_II_OPENMP_SIMD_PRAGMA
3476 for (std::size_t i = begin; i < end; ++i)
3477 // for efficiency reason, write back to temp_vector that is
3478 // already read (avoid read-for-ownership)
3479 tmp_vector[i] =
3480 factor1_plus_1 * solution[i] +
3481 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3482 }
3483 else
3484 {
3485 // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
3486 // + f_2 * P^{-1} * (b-A*x^{n})
3487 DEAL_II_OPENMP_SIMD_PRAGMA
3488 for (std::size_t i = begin; i < end; ++i)
3489 // for efficiency reason, write back to temp_vector, which is
3490 // already modified during vmult (in best case, the modified
3491 // values are not written back to main memory yet so that
3492 // we do not have to pay additional costs for writing and
3493 // read-for-ownership)
3494 tmp_vector[i] =
3495 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
3496 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3497 }
3498 }
3499
3500 const Number *rhs;
3501 const Number *matrix_diagonal_inverse;
3502 const unsigned int iteration_index;
3503 const Number factor1;
3504 const Number factor2;
3505 mutable Number *solution_old;
3506 mutable Number *tmp_vector;
3507 mutable Number *solution;
3508 };
3509
3510 template <typename Number>
3511 struct VectorUpdatesRange : public ::parallel::ParallelForInteger
3512 {
3513 VectorUpdatesRange(const VectorUpdater<Number> &updater,
3514 const std::size_t size)
3515 : updater(updater)
3516 {
3517 if (size < internal::VectorImplementation::minimum_parallel_grain_size)
3518 VectorUpdatesRange::apply_to_subrange(0, size);
3519 else
3520 apply_parallel(
3521 0,
3522 size,
3523 internal::VectorImplementation::minimum_parallel_grain_size);
3524 }
3525
3526 ~VectorUpdatesRange() override = default;
3527
3528 virtual void
3529 apply_to_subrange(const std::size_t begin,
3530 const std::size_t end) const override
3531 {
3532 updater.apply_to_subrange(begin, end);
3533 }
3534
3535 const VectorUpdater<Number> &updater;
3536 };
3537
3538 // selection for diagonal matrix around deal.II vector
3539 template <typename Number>
3540 inline void
3541 vector_updates(
3542 const ::Vector<Number> &rhs,
3543 const ::DiagonalMatrix<::Vector<Number>> &jacobi,
3544 const unsigned int iteration_index,
3545 const double factor1,
3546 const double factor2,
3547 ::Vector<Number> &solution_old,
3548 ::Vector<Number> &temp_vector1,
3549 ::Vector<Number> &,
3550 ::Vector<Number> &solution)
3551 {
3552 VectorUpdater<Number> upd(rhs.begin(),
3553 jacobi.get_vector().begin(),
3554 iteration_index,
3555 factor1,
3556 factor2,
3557 solution_old.begin(),
3558 temp_vector1.begin(),
3559 solution.begin());
3560 VectorUpdatesRange<Number>(upd, rhs.size());
3561
3562 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3563 solution.swap(temp_vector1);
3564 solution_old.swap(temp_vector1);
3565 }
3566
3567 // selection for diagonal matrix around parallel deal.II vector
3568 template <typename Number>
3569 inline void
3570 vector_updates(
3571 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3572 const ::DiagonalMatrix<
3573 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>> &jacobi,
3574 const unsigned int iteration_index,
3575 const double factor1,
3576 const double factor2,
3577 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3578 &solution_old,
3579 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3580 &temp_vector1,
3581 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &,
3582 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
3583 {
3584 VectorUpdater<Number> upd(rhs.begin(),
3585 jacobi.get_vector().begin(),
3586 iteration_index,
3587 factor1,
3588 factor2,
3589 solution_old.begin(),
3590 temp_vector1.begin(),
3591 solution.begin());
3592 VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
3593
3594 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3595 solution.swap(temp_vector1);
3596 solution_old.swap(temp_vector1);
3597 }
3598
3599 // We need to have a separate declaration for static const members
3600
3601 // general case and the case that the preconditioner can work on
3602 // ranges (covered by vector_updates())
3603 template <
3604 typename MatrixType,
3605 typename VectorType,
3606 typename PreconditionerType,
3607 std::enable_if_t<
3608 !has_vmult_with_std_functions<MatrixType,
3609 VectorType,
3610 PreconditionerType> &&
3611 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3612 VectorType> &&
3613 has_vmult_with_std_functions_for_precondition<MatrixType,
3614 VectorType>),
3615 int> * = nullptr>
3616 inline void
3617 vmult_and_update(const MatrixType &matrix,
3618 const PreconditionerType &preconditioner,
3619 const VectorType &rhs,
3620 const unsigned int iteration_index,
3621 const double factor1,
3622 const double factor2,
3623 VectorType &solution,
3624 VectorType &solution_old,
3625 VectorType &temp_vector1,
3626 VectorType &temp_vector2)
3627 {
3628 if (iteration_index > 0)
3629 matrix.vmult(temp_vector1, solution);
3630 vector_updates(rhs,
3631 preconditioner,
3632 iteration_index,
3633 factor1,
3634 factor2,
3635 solution_old,
3636 temp_vector1,
3637 temp_vector2,
3638 solution);
3639 }
3640
3641 // case that both the operator and the preconditioner can work on
3642 // subranges
3643 template <
3644 typename MatrixType,
3645 typename VectorType,
3646 typename PreconditionerType,
3647 std::enable_if_t<
3648 !has_vmult_with_std_functions<MatrixType,
3649 VectorType,
3650 PreconditionerType> &&
3651 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3652 VectorType> &&
3653 has_vmult_with_std_functions_for_precondition<MatrixType,
3654 VectorType>),
3655 int> * = nullptr>
3656 inline void
3657 vmult_and_update(const MatrixType &matrix,
3658 const PreconditionerType &preconditioner,
3659 const VectorType &rhs,
3660 const unsigned int iteration_index,
3661 const double factor1_,
3662 const double factor2_,
3663 VectorType &solution,
3664 VectorType &solution_old,
3665 VectorType &temp_vector1,
3666 VectorType &temp_vector2)
3667 {
3668 using Number = typename VectorType::value_type;
3669
3670 const Number factor1 = factor1_;
3671 const Number factor1_plus_1 = 1. + factor1_;
3672 const Number factor2 = factor2_;
3673
3674 if (iteration_index == 0)
3675 {
3676 preconditioner.vmult(
3677 temp_vector2,
3678 rhs,
3679 [&](const unsigned int start_range, const unsigned int end_range) {
3680 // zero 'solution' before running the vmult operation
3681 if (end_range > start_range)
3682 std::memset(temp_vector2.begin() + start_range,
3683 0,
3684 sizeof(Number) * (end_range - start_range));
3685 },
3686 [&](const unsigned int start_range, const unsigned int end_range) {
3687 const auto tmp_ptr = temp_vector2.begin();
3688
3689 DEAL_II_OPENMP_SIMD_PRAGMA
3690 for (std::size_t i = start_range; i < end_range; ++i)
3691 tmp_ptr[i] *= factor2;
3692 });
3693 }
3694 else
3695 {
3696 temp_vector1.reinit(rhs, true);
3697 temp_vector2.reinit(rhs, true);
3698
3699 // 1) compute residual (including operator application)
3700 matrix.vmult(
3701 temp_vector1,
3702 solution,
3703 [&](const unsigned int start_range, const unsigned int end_range) {
3704 // zero 'temp_vector1' before running the vmult
3705 // operation
3706 if (end_range > start_range)
3707 std::memset(temp_vector1.begin() + start_range,
3708 0,
3709 sizeof(Number) * (end_range - start_range));
3710 },
3711 [&](const unsigned int start_range, const unsigned int end_range) {
3712 const auto rhs_ptr = rhs.begin();
3713 const auto tmp_ptr = temp_vector1.begin();
3714
3715 DEAL_II_OPENMP_SIMD_PRAGMA
3716 for (std::size_t i = start_range; i < end_range; ++i)
3717 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3718 });
3719
3720 // 2) perform vector updates (including preconditioner application)
3721 preconditioner.vmult(
3722 temp_vector2,
3723 temp_vector1,
3724 [&](const unsigned int start_range, const unsigned int end_range) {
3725 // zero 'temp_vector2' before running the vmult
3726 // operation
3727 if (end_range > start_range)
3728 std::memset(temp_vector2.begin() + start_range,
3729 0,
3730 sizeof(Number) * (end_range - start_range));
3731 },
3732 [&](const unsigned int start_range, const unsigned int end_range) {
3733 const auto solution_ptr = solution.begin();
3734 const auto solution_old_ptr = solution_old.begin();
3735 const auto tmp_ptr = temp_vector2.begin();
3736
3737 if (iteration_index == 1)
3738 {
3739 DEAL_II_OPENMP_SIMD_PRAGMA
3740 for (std::size_t i = start_range; i < end_range; ++i)
3741 tmp_ptr[i] =
3742 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3743 }
3744 else
3745 {
3746 DEAL_II_OPENMP_SIMD_PRAGMA
3747 for (std::size_t i = start_range; i < end_range; ++i)
3748 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3749 factor1 * solution_old_ptr[i] +
3750 factor2 * tmp_ptr[i];
3751 }
3752 });
3753 }
3754 solution.swap(temp_vector2);
3755 solution_old.swap(temp_vector2);
3756 }
3757
3758 // case that the operator can work on subranges and the preconditioner
3759 // is a diagonal
3760 template <typename MatrixType,
3761 typename VectorType,
3762 typename PreconditionerType,
3763 std::enable_if_t<has_vmult_with_std_functions<MatrixType,
3764 VectorType,
3765 PreconditionerType>,
3766 int> * = nullptr>
3767 inline void
3768 vmult_and_update(const MatrixType &matrix,
3769 const PreconditionerType &preconditioner,
3770 const VectorType &rhs,
3771 const unsigned int iteration_index,
3772 const double factor1,
3773 const double factor2,
3774 VectorType &solution,
3775 VectorType &solution_old,
3776 VectorType &temp_vector1,
3777 VectorType &)
3778 {
3779 using Number = typename VectorType::value_type;
3780 VectorUpdater<Number> updater(rhs.begin(),
3781 preconditioner.get_vector().begin(),
3782 iteration_index,
3783 factor1,
3784 factor2,
3785 solution_old.begin(),
3786 temp_vector1.begin(),
3787 solution.begin());
3788 if (iteration_index > 0)
3789 matrix.vmult(
3790 temp_vector1,
3791 solution,
3792 [&](const unsigned int start_range, const unsigned int end_range) {
3793 // zero 'temp_vector1' before running the vmult
3794 // operation
3795 if (end_range > start_range)
3796 std::memset(temp_vector1.begin() + start_range,
3797 0,
3798 sizeof(Number) * (end_range - start_range));
3799 },
3800 [&](const unsigned int start_range, const unsigned int end_range) {
3801 if (end_range > start_range)
3802 updater.apply_to_subrange(start_range, end_range);
3803 });
3804 else
3805 updater.apply_to_subrange(0U, solution.locally_owned_size());
3806
3807 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3808 solution.swap(temp_vector1);
3809 solution_old.swap(temp_vector1);
3810 }
3811
3812 template <typename MatrixType, typename PreconditionerType>
3813 inline void
3814 initialize_preconditioner(
3815 const MatrixType & /*matrix*/,
3816 std::shared_ptr<PreconditionerType> &preconditioner)
3817 {
3818 (void)preconditioner;
3819 AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
3820 }
3821
3822 template <typename MatrixType, typename VectorType>
3823 inline void
3824 initialize_preconditioner(
3825 const MatrixType &matrix,
3826 std::shared_ptr<::DiagonalMatrix<VectorType>> &preconditioner)
3827 {
3828 if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
3829 {
3830 if (preconditioner.get() == nullptr)
3831 preconditioner =
3832 std::make_shared<::DiagonalMatrix<VectorType>>();
3833
3834 Assert(
3835 preconditioner->m() == 0,
3836 ExcMessage(
3837 "Preconditioner appears to be initialized but not sized correctly"));
3838
3839 // This part only works in serial
3840 if (preconditioner->m() != matrix.m())
3841 {
3842 preconditioner->get_vector().reinit(matrix.m());
3843 for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
3844 preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
3845 }
3846 }
3847 }
3848 } // namespace PreconditionChebyshevImplementation
3849} // namespace internal
3850
3851
3852
3853template <typename MatrixType, typename VectorType, typename PreconditionerType>
3854inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3855 AdditionalData::AdditionalData(const unsigned int degree,
3856 const double smoothing_range,
3857 const unsigned int eig_cg_n_iterations,
3858 const double eig_cg_residual,
3859 const double max_eigenvalue,
3860 const EigenvalueAlgorithm eigenvalue_algorithm,
3861 const PolynomialType polynomial_type,
3862 const double safety_factor)
3863 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3864 smoothing_range,
3865 eig_cg_n_iterations,
3866 eig_cg_residual,
3867 max_eigenvalue,
3868 eigenvalue_algorithm,
3869 safety_factor)
3870 , degree(degree)
3871 , polynomial_type(polynomial_type)
3872{}
3873
3874
3875
3876template <typename MatrixType, typename VectorType, typename PreconditionerType>
3877inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3878 PreconditionChebyshev()
3879 : theta(1.)
3880 , delta(1.)
3881 , eigenvalues_are_initialized(false)
3882{
3883 static_assert(
3884 std::is_same_v<size_type, typename VectorType::size_type>,
3885 "PreconditionChebyshev and VectorType must have the same size_type.");
3886}
3887
3888
3889
3890template <typename MatrixType, typename VectorType, typename PreconditionerType>
3891inline void
3892PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::initialize(
3893 const MatrixType &matrix,
3894 const AdditionalData &additional_data)
3895{
3896 matrix_ptr = &matrix;
3897 data = additional_data;
3898 Assert(data.degree > 0,
3899 ExcMessage("The degree of the Chebyshev method must be positive."));
3900 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3901 matrix, data.preconditioner);
3902 eigenvalues_are_initialized = false;
3903}
3904
3905
3906
3907template <typename MatrixType, typename VectorType, typename PreconditionerType>
3908inline void
3909PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::clear()
3910{
3911 eigenvalues_are_initialized = false;
3912 theta = delta = 1.0;
3913 matrix_ptr = nullptr;
3914 {
3915 VectorType empty_vector;
3916 solution_old.reinit(empty_vector);
3917 temp_vector1.reinit(empty_vector);
3918 temp_vector2.reinit(empty_vector);
3919 }
3920 data.preconditioner.reset();
3921}
3922
3923
3924
3925template <typename MatrixType, typename VectorType, typename PreconditionerType>
3926inline typename internal::EigenvalueInformation
3927PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3928 estimate_eigenvalues(const VectorType &src) const
3929{
3930 Assert(eigenvalues_are_initialized == false, ExcInternalError());
3931
3932 solution_old.reinit(src);
3933 temp_vector1.reinit(src, true);
3934
3935 auto info = internal::estimate_eigenvalues<MatrixType>(data,
3936 matrix_ptr,
3937 solution_old,
3938 temp_vector1,
3939 data.degree,
3940 data.safety_factor);
3941
3942 const double alpha = (data.smoothing_range > 1. ?
3943 info.max_eigenvalue_estimate / data.smoothing_range :
3944 std::min(0.9 * info.max_eigenvalue_estimate,
3945 info.min_eigenvalue_estimate));
3946
3947 // in case the user set the degree to invalid unsigned int, we have to
3948 // determine the number of necessary iterations from the Chebyshev error
3949 // estimate, given the target tolerance specified by smoothing_range. This
3950 // estimate is based on the error formula given in section 5.1 of
3951 // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
3952 if (data.degree == numbers::invalid_unsigned_int)
3953 {
3954 // In solver mode, smoothing_range is interpreted as the relative
3955 // target tolerance and must be strictly less than one. Otherwise the
3956 // Chebyshev error formula below evaluates the square root of a
3957 // negative number, producing a NaN that silently gets cast to an
3958 // unsigned int with implementation-defined (typically zero) result.
3959 // This would turn the "solver" into a single damped Jacobi step,
3960 // which is almost certainly not what the user intended.
3961 Assert(data.smoothing_range < 1.,
3962 ExcMessage(
3963 "When using PreconditionChebyshev as a solver (i.e., when "
3964 "AdditionalData::degree == numbers::invalid_unsigned_int), "
3965 "AdditionalData::smoothing_range must be strictly less than "
3966 "one because it is interpreted as the relative target "
3967 "tolerance of the Chebyshev iteration. A value >= 1 is only "
3968 "meaningful when PreconditionChebyshev is used as a smoother "
3969 "with a fixed polynomial degree."));
3970
3971 const double actual_range = info.max_eigenvalue_estimate / alpha;
3972 const double sigma = (1. - std::sqrt(1. / actual_range)) /
3973 (1. + std::sqrt(1. / actual_range));
3974 const double eps = data.smoothing_range;
3975 const_cast<
3976 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(
3977 this)
3978 ->data.degree =
3979 1 + static_cast<unsigned int>(
3980 std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
3981 std::log(1. / sigma));
3982 }
3983
3984 info.degree = data.degree;
3985
3986 const_cast<
3987 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3988 ->delta =
3989 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3990 (info.max_eigenvalue_estimate) :
3991 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3992 const_cast<
3993 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3994 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3995
3996 // We do not need the second temporary vector in case we have a
3997 // DiagonalMatrix as preconditioner and use deal.II's own vectors
3998 using NumberType = typename VectorType::value_type;
3999 if (std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> ==
4000 false ||
4001 (std::is_same_v<VectorType, ::Vector<NumberType>> == false &&
4002 ((std::is_same_v<
4003 VectorType,
4004 LinearAlgebra::distributed::Vector<NumberType, MemorySpace::Host>> ==
4005 false) ||
4006 (std::is_same_v<VectorType,
4007 LinearAlgebra::distributed::
4008 Vector<NumberType, MemorySpace::Default>> == false))))
4009 temp_vector2.reinit(src, true);
4010 else
4011 {
4012 VectorType empty_vector;
4013 temp_vector2.reinit(empty_vector);
4014 }
4015
4016 const_cast<
4017 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
4018 ->eigenvalues_are_initialized = true;
4019
4020 return info;
4021}
4022
4023
4024
4025template <typename MatrixType, typename VectorType, typename PreconditionerType>
4026template <bool do_transpose>
4027inline void
4028PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
4029 apply_internal(const bool zero_out_dst,
4030 const VectorType &rhs,
4031 VectorType &solution) const
4032{
4033 std::scoped_lock lock(mutex);
4034 if (eigenvalues_are_initialized == false)
4035 estimate_eigenvalues(rhs);
4036
4037 if constexpr (do_transpose)
4038 {
4039 matrix_ptr->Tvmult(temp_vector1, solution);
4040 internal::PreconditionChebyshevImplementation::vector_updates(
4041 rhs,
4042 *data.preconditioner,
4043 zero_out_dst ? 0 : 1,
4044 0.,
4045 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4046 (4. / (3. * delta)) :
4047 (1. / theta),
4048 solution_old,
4049 temp_vector1,
4050 temp_vector2,
4051 solution);
4052 }
4053 else
4054 internal::PreconditionChebyshevImplementation::vmult_and_update(
4055 *matrix_ptr,
4056 *data.preconditioner,
4057 rhs,
4058 zero_out_dst ? 0 : 1,
4059 0.,
4060 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4061 (4. / (3. * delta)) :
4062 (1. / theta),
4063 solution,
4064 solution_old,
4065 temp_vector1,
4066 temp_vector2);
4067
4068 // if delta is zero, we do not need to iterate because the updates will be
4069 // zero
4070 if (data.degree < 2 || std::abs(delta) < 1e-40)
4071 return;
4072
4073 double rhok = delta / theta, sigma = theta / delta;
4074 for (unsigned int k = 0; k < data.degree - 1; ++k)
4075 {
4076 double factor1 = 0.0;
4077 double factor2 = 0.0;
4078
4079 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4080 {
4081 factor1 = (2 * k + 1.) / (2 * k + 5.);
4082 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4083 }
4084 else
4085 {
4086 const double rhokp = 1. / (2. * sigma - rhok);
4087 factor1 = rhokp * rhok;
4088 factor2 = 2. * rhokp / delta;
4089 rhok = rhokp;
4090 }
4091
4092 if constexpr (do_transpose)
4093 {
4094 matrix_ptr->Tvmult(temp_vector1, solution);
4095 internal::PreconditionChebyshevImplementation::vector_updates(
4096 rhs,
4097 *data.preconditioner,
4098 k + (zero_out_dst ? 1 : 2),
4099 factor1,
4100 factor2,
4101 solution_old,
4102 temp_vector1,
4103 temp_vector2,
4104 solution);
4105 }
4106 else
4107 {
4108 internal::PreconditionChebyshevImplementation::vmult_and_update(
4109 *matrix_ptr,
4110 *data.preconditioner,
4111 rhs,
4112 k + (zero_out_dst ? 1 : 2),
4113 factor1,
4114 factor2,
4115 solution,
4116 solution_old,
4117 temp_vector1,
4118 temp_vector2);
4119 }
4120 }
4121}
4122
4123
4124
4125template <typename MatrixType, typename VectorType, typename PreconditionerType>
4126inline void
4127PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::vmult(
4128 VectorType &solution,
4129 const VectorType &rhs) const
4130{
4131 apply_internal<false>(/* zero_out_dst */ true, rhs, solution);
4132}
4133
4134
4135
4136template <typename MatrixType, typename VectorType, typename PreconditionerType>
4137inline void
4138PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tvmult(
4139 VectorType &solution,
4140 const VectorType &rhs) const
4141{
4142 apply_internal<true>(/* zero_out_dst */ true, rhs, solution);
4143}
4144
4145
4146
4147template <typename MatrixType, typename VectorType, typename PreconditionerType>
4148inline void
4149PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::step(
4150 VectorType &solution,
4151 const VectorType &rhs) const
4152{
4153 apply_internal<false>(/* zero_out_dst */ false, rhs, solution);
4154}
4155
4156
4157
4158template <typename MatrixType, typename VectorType, typename PreconditionerType>
4159inline void
4160PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tstep(
4161 VectorType &solution,
4162 const VectorType &rhs) const
4163{
4164 apply_internal<true>(/* zero_out_dst */ false, rhs, solution);
4165}
4166
4167
4168
4169template <typename MatrixType, typename VectorType, typename PreconditionerType>
4170inline typename PreconditionChebyshev<MatrixType,
4171 VectorType,
4172 PreconditionerType>::size_type
4173PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::m() const
4174{
4175 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4176 return matrix_ptr->m();
4177}
4178
4179
4180
4181template <typename MatrixType, typename VectorType, typename PreconditionerType>
4182inline typename PreconditionChebyshev<MatrixType,
4183 VectorType,
4184 PreconditionerType>::size_type
4185PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::n() const
4186{
4187 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4188 return matrix_ptr->n();
4189}
4190
4191#endif // DOXYGEN
4192
4193DEAL_II_NAMESPACE_CLOSE
4194
4195#endif
end
*  iterator end()
Definition ALE_Finite_Strain_Plasticity.h:10414
const_iterator
*  *  const_iterator()=default
Definition petsc_matrix_base.h:183
begin
*  *  iterator begin()
Definition ALE_Finite_Strain_Plasticity.h:10413
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:1933
IndexSet::nth_index_in_set
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1982
LinearAlgebra::distributed::Vector::sadd
void sadd(const Number s, const Number a, const Vector< Number, MemorySpace > &V)
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
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::estimate_eigenvalues
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
PreconditionChebyshev::Tstep
void Tstep(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::EigenvalueInformation
internal::EigenvalueInformation EigenvalueInformation
Definition precondition.h:2243
PreconditionChebyshev::size_type
types::global_dof_index size_type
Definition precondition.h:2115
PreconditionChebyshev::apply_internal
void apply_internal(const bool zero_out_dst, const VectorType &src, VectorType &dst) const
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:2277
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::size_type
types::global_dof_index size_type
Definition precondition.h:233
PreconditionIdentity::Tvmult_add
void Tvmult_add(VectorType &, const VectorType &) const
PreconditionJacobi::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1669
PreconditionJacobi::BaseClass
PreconditionRelaxation< MatrixType, PreconditionerType > BaseClass
Definition precondition.h:1663
PreconditionJacobi::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionJacobi::PreconditionerType
internal::PreconditionRelaxation::PreconditionJacobiImpl< MatrixType > PreconditionerType
Definition precondition.h:1661
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::BaseClass
PreconditionRelaxation< MatrixType, PreconditionerType > BaseClass
Definition precondition.h:1841
PreconditionPSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &additional_data)
PreconditionPSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionPSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1839
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::EigenvalueAlgorithm
internal::EigenvalueAlgorithm EigenvalueAlgorithm
Definition precondition.h:589
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< MatrixType, PreconditionerType >::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::EigenvalueInformation
internal::EigenvalueInformation EigenvalueInformation
Definition precondition.h:674
PreconditionRelaxation::m
size_type m() const
PreconditionRelaxation::step
void step(VectorType &x, const VectorType &rhs) const
PreconditionRelaxation< MatrixType, PreconditionerType >::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::BaseClass
PreconditionRelaxation< MatrixType, PreconditionerType > BaseClass
Definition precondition.h:1733
PreconditionSOR::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1739
PreconditionSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1731
PreconditionSSOR< SparseMatrix >::PreconditionerType
internal::PreconditionRelaxation::PreconditionSSORImpl< SparseMatrix > PreconditionerType
Definition precondition.h:1783
PreconditionSSOR< SparseMatrix >::BaseClass
PreconditionRelaxation< SparseMatrix, PreconditionerType > BaseClass
Definition precondition.h:1785
PreconditionSSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionSSOR< SparseMatrix >::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1791
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)
PreconditionUseMatrix::function_ptr
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition precondition.h:512
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::size
virtual size_type size() const override
Vector::swap
virtual void swap(Vector< Number > &v) noexcept
Vector::begin
iterator begin()
DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:213
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
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:138
PlasticityLab::Number
*  *  *  RotationFunction< dim, Number >::RotationFunction Number(dim)
Definition ALE_Finite_Strain_Plasticity.h:8749
internal::VectorImplementation::minimum_parallel_grain_size
unsigned int minimum_parallel_grain_size
Definition parallel.cc:48
internal::is_supported_operation
constexpr bool is_supported_operation
Definition template_constraints.h:159
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:236
std
STL namespace.
std::log
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6647
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6758
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6667
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6725
types::global_dof_index
unsigned int global_dof_index
Definition types.h:92
residual
*  *  residual
Definition nonlinear-heat_transfer_with_AD_NOX.h:1429
observer_pointer.h
solver_cg.h
IsBlockVector::value
static const bool value
Definition block_vector_base.h:85
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)
PreconditionChebyshev::AdditionalData::EigenvalueAlgorithm
internal::EigenvalueAlgorithm EigenvalueAlgorithm
Definition precondition.h:2124
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
Definition symmetric_tensor.h:3134
template_constraints.h
Generated by doxygen 1.16.1