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