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