deal.II version GIT relicensing-1941-ga9013bea74 2024-10-08 06:50: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 Subscriptor
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 Subscriptor
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 Subscriptor
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 Subscriptor
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
1296 if (i == 0)
1297 {
1298 Number *dst_ptr = dst.begin();
1299 const Number *src_ptr = src.begin();
1300
1301 preconditioner.vmult(
1302 dst,
1303 src,
1304 [&](const unsigned int start_range, const unsigned int end_range) {
1305 // zero 'dst' before running the vmult operation
1306 if (end_range > start_range)
1307 std::memset(dst.begin() + start_range,
1308 0,
1309 sizeof(Number) * (end_range - start_range));
1310 },
1311 [&](const unsigned int start_range, const unsigned int end_range) {
1312 if (relaxation == 1.0)
1313 return; // nothing to do
1314
1315 const auto src_ptr = src.begin();
1316 const auto dst_ptr = dst.begin();
1317
1318 DEAL_II_OPENMP_SIMD_PRAGMA
1319 for (std::size_t i = start_range; i < end_range; ++i)
1320 dst_ptr[i] *= relaxation;
1321 });
1322 }
1323 else
1324 {
1325 tmp.reinit(src, true);
1326
1327 Assert(transposed == false, ExcNotImplemented());
1328
1329 A.vmult(tmp, dst);
1330
1331 preconditioner.vmult(
1332 dst,
1333 tmp,
1334 [&](const unsigned int start_range, const unsigned int end_range) {
1335 const auto src_ptr = src.begin();
1336 const auto tmp_ptr = tmp.begin();
1337
1338 if (relaxation == 1.0)
1339 {
1340 DEAL_II_OPENMP_SIMD_PRAGMA
1341 for (std::size_t i = start_range; i < end_range; ++i)
1342 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1343 }
1344 else
1345 {
1346 // note: we scale the residual here to be able to add into
1347 // the dst vector, which contains the solution from the last
1348 // iteration
1349 DEAL_II_OPENMP_SIMD_PRAGMA
1350 for (std::size_t i = start_range; i < end_range; ++i)
1351 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1352 }
1353 },
1354 [&](const unsigned int, const unsigned int) {
1355 // nothing to do, since scaling by the relaxation factor
1356 // has been done in the pre operation
1357 });
1358 }
1359 }
1360
1361 // 2) specialized implementation with a preconditioner and a matrix that
1362 // accepts ranges
1363 template <
1364 typename MatrixType,
1365 typename PreconditionerType,
1366 typename VectorType,
1367 std::enable_if_t<
1368 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1369 VectorType> &&
1370 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1371 int> * = nullptr>
1372 void
1373 step_operations(const MatrixType &A,
1374 const PreconditionerType &preconditioner,
1375 VectorType &dst,
1376 const VectorType &src,
1377 const double relaxation,
1378 VectorType &tmp,
1379 VectorType &,
1380 const unsigned int i,
1381 const bool transposed)
1382 {
1383 (void)transposed;
1384 using Number = typename VectorType::value_type;
1385
1386 if (i == 0)
1387 {
1388 Number *dst_ptr = dst.begin();
1389 const Number *src_ptr = src.begin();
1390
1391 preconditioner.vmult(
1392 dst,
1393 src,
1394 [&](const unsigned int start_range, const unsigned int end_range) {
1395 // zero 'dst' before running the vmult operation
1396 if (end_range > start_range)
1397 std::memset(dst.begin() + start_range,
1398 0,
1399 sizeof(Number) * (end_range - start_range));
1400 },
1401 [&](const unsigned int start_range, const unsigned int end_range) {
1402 if (relaxation == 1.0)
1403 return; // nothing to do
1404
1405 const auto src_ptr = src.begin();
1406 const auto dst_ptr = dst.begin();
1407
1408 DEAL_II_OPENMP_SIMD_PRAGMA
1409 for (std::size_t i = start_range; i < end_range; ++i)
1410 dst_ptr[i] *= relaxation;
1411 });
1412 }
1413 else
1414 {
1415 tmp.reinit(src, true);
1416
1417 Assert(transposed == false, ExcNotImplemented());
1418
1419 A.vmult(
1420 tmp,
1421 dst,
1422 [&](const unsigned int start_range, const unsigned int end_range) {
1423 // zero 'tmp' before running the vmult
1424 // operation
1425 if (end_range > start_range)
1426 std::memset(tmp.begin() + start_range,
1427 0,
1428 sizeof(Number) * (end_range - start_range));
1429 },
1430 [&](const unsigned int start_range, const unsigned int end_range) {
1431 const auto src_ptr = src.begin();
1432 const auto tmp_ptr = tmp.begin();
1433
1434 if (relaxation == 1.0)
1435 {
1436 DEAL_II_OPENMP_SIMD_PRAGMA
1437 for (std::size_t i = start_range; i < end_range; ++i)
1438 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1439 }
1440 else
1441 {
1442 // note: we scale the residual here to be able to add into
1443 // the dst vector, which contains the solution from the last
1444 // iteration
1445 DEAL_II_OPENMP_SIMD_PRAGMA
1446 for (std::size_t i = start_range; i < end_range; ++i)
1447 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1448 }
1449 });
1450
1451 preconditioner.vmult(dst, tmp, [](const auto, const auto) {
1452 // note: `dst` vector does not have to be zeroed out
1453 // since we add the result into it
1454 });
1455 }
1456 }
1457
1458 // 3) specialized implementation for inverse-diagonal preconditioner
1459 template <typename MatrixType,
1460 typename VectorType,
1461 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1462 !has_vmult_with_std_functions<
1463 MatrixType,
1464 VectorType,
1465 ::DiagonalMatrix<VectorType>>,
1466 VectorType> * = nullptr>
1467 void
1468 step_operations(const MatrixType &A,
1469 const ::DiagonalMatrix<VectorType> &preconditioner,
1470 VectorType &dst,
1471 const VectorType &src,
1472 const double relaxation,
1473 VectorType &tmp,
1474 VectorType &,
1475 const unsigned int i,
1476 const bool transposed)
1477 {
1478 using Number = typename VectorType::value_type;
1479
1480 if (i == 0)
1481 {
1482 Number *dst_ptr = dst.begin();
1483 const Number *src_ptr = src.begin();
1484 const Number *diag_ptr = preconditioner.get_vector().begin();
1485
1486 if (relaxation == 1.0)
1487 {
1488 DEAL_II_OPENMP_SIMD_PRAGMA
1489 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1490 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1491 }
1492 else
1493 {
1494 DEAL_II_OPENMP_SIMD_PRAGMA
1495 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1496 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1497 }
1498 }
1499 else
1500 {
1501 tmp.reinit(src, true);
1502
1503 Number *dst_ptr = dst.begin();
1504 const Number *src_ptr = src.begin();
1505 const Number *tmp_ptr = tmp.begin();
1506 const Number *diag_ptr = preconditioner.get_vector().begin();
1507
1508 if (transposed)
1509 Tvmult(A, tmp, dst);
1510 else
1511 A.vmult(tmp, dst);
1512
1513 if (relaxation == 1.0)
1514 {
1515 DEAL_II_OPENMP_SIMD_PRAGMA
1516 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1517 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1518 }
1519 else
1520 {
1521 DEAL_II_OPENMP_SIMD_PRAGMA
1522 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1523 dst_ptr[i] +=
1524 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1525 }
1526 }
1527 }
1528
1529 // 4) specialized implementation for inverse-diagonal preconditioner and
1530 // matrix that accepts ranges
1531 template <typename MatrixType,
1532 typename VectorType,
1533 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1534 has_vmult_with_std_functions<
1535 MatrixType,
1536 VectorType,
1537 ::DiagonalMatrix<VectorType>>,
1538 VectorType> * = nullptr>
1539 void
1540 step_operations(const MatrixType &A,
1541 const ::DiagonalMatrix<VectorType> &preconditioner,
1542 VectorType &dst,
1543 const VectorType &src,
1544 const double relaxation,
1545 VectorType &tmp,
1546 VectorType &,
1547 const unsigned int i,
1548 const bool transposed)
1549 {
1550 (void)transposed;
1551 using Number = typename VectorType::value_type;
1552
1553 if (i == 0)
1554 {
1555 Number *dst_ptr = dst.begin();
1556 const Number *src_ptr = src.begin();
1557 const Number *diag_ptr = preconditioner.get_vector().begin();
1558
1559 if (relaxation == 1.0)
1560 {
1561 DEAL_II_OPENMP_SIMD_PRAGMA
1562 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1563 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1564 }
1565 else
1566 {
1567 DEAL_II_OPENMP_SIMD_PRAGMA
1568 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1569 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1570 }
1571 }
1572 else
1573 {
1574 tmp.reinit(src, true);
1575
1576 Assert(transposed == false, ExcNotImplemented());
1577
1578 A.vmult(
1579 tmp,
1580 dst,
1581 [&](const unsigned int start_range, const unsigned int end_range) {
1582 // zero 'tmp' before running the vmult operation
1583 if (end_range > start_range)
1584 std::memset(tmp.begin() + start_range,
1585 0,
1586 sizeof(Number) * (end_range - start_range));
1587 },
1588 [&](const unsigned int begin, const unsigned int end) {
1589 const Number *dst_ptr = dst.begin();
1590 const Number *src_ptr = src.begin();
1591 Number *tmp_ptr = tmp.begin();
1592 const Number *diag_ptr = preconditioner.get_vector().begin();
1593
1594 // for efficiency reason, write back to temp_vector that is
1595 // already read (avoid read-for-ownership)
1596 if (relaxation == 1.0)
1597 {
1598 DEAL_II_OPENMP_SIMD_PRAGMA
1599 for (std::size_t i = begin; i < end; ++i)
1600 tmp_ptr[i] =
1601 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1602 }
1603 else
1604 {
1605 DEAL_II_OPENMP_SIMD_PRAGMA
1606 for (std::size_t i = begin; i < end; ++i)
1607 tmp_ptr[i] = dst_ptr[i] + relaxation *
1608 (src_ptr[i] - tmp_ptr[i]) *
1609 diag_ptr[i];
1610 }
1611 });
1612
1613 tmp.swap(dst);
1614 }
1615 }
1616
1617 } // namespace PreconditionRelaxation
1618} // namespace internal
1619
1620#endif
1621
1622
1623
1650template <typename MatrixType = SparseMatrix<double>>
1651class PreconditionJacobi
1652 : public PreconditionRelaxation<
1653 MatrixType,
1654 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1655{
1656 using PreconditionerType =
1657 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1658 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1659
1660public:
1664 using AdditionalData = typename BaseClass::AdditionalData;
1665
1669 void
1670 initialize(const MatrixType &A,
1671 const AdditionalData &parameters = AdditionalData());
1672};
1673
1674
1720template <typename MatrixType = SparseMatrix<double>>
1721class PreconditionSOR
1722 : public PreconditionRelaxation<
1723 MatrixType,
1724 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1725{
1726 using PreconditionerType =
1727 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1728 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1729
1730public:
1734 using AdditionalData = typename BaseClass::AdditionalData;
1735
1739 void
1740 initialize(const MatrixType &A,
1741 const AdditionalData &parameters = AdditionalData());
1742};
1743
1744
1745
1772template <typename MatrixType = SparseMatrix<double>>
1773class PreconditionSSOR
1774 : public PreconditionRelaxation<
1775 MatrixType,
1776 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1777{
1778 using PreconditionerType =
1779 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1780 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1781
1782public:
1786 using AdditionalData = typename BaseClass::AdditionalData;
1787
1793 void
1794 initialize(const MatrixType &A,
1795 const AdditionalData &parameters = AdditionalData());
1796};
1797
1798
1828template <typename MatrixType = SparseMatrix<double>>
1829class PreconditionPSOR
1830 : public PreconditionRelaxation<
1831 MatrixType,
1832 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1833{
1834 using PreconditionerType =
1835 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1836 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1837
1838public:
1842 using size_type = typename BaseClass::size_type;
1843
1847 class AdditionalData
1848 {
1849 public:
1860 AdditionalData(const std::vector<size_type> &permutation,
1861 const std::vector<size_type> &inverse_permutation,
1862 const typename BaseClass::AdditionalData &parameters =
1863 typename BaseClass::AdditionalData());
1864
1868 const std::vector<size_type> &permutation;
1872 const std::vector<size_type> &inverse_permutation;
1876 typename BaseClass::AdditionalData parameters;
1877 };
1878
1890 void
1891 initialize(const MatrixType &A,
1892 const std::vector<size_type> &permutation,
1893 const std::vector<size_type> &inverse_permutation,
1894 const typename BaseClass::AdditionalData &parameters =
1895 typename BaseClass::AdditionalData());
1896
1907 void
1908 initialize(const MatrixType &A, const AdditionalData &additional_data);
1909};
1910
1911
1912
2100template <typename MatrixType = SparseMatrix<double>,
2101 typename VectorType = Vector<double>,
2102 typename PreconditionerType = DiagonalMatrix<VectorType>>
2103class PreconditionChebyshev : public Subscriptor
2104{
2105public:
2109 using size_type = types::global_dof_index;
2110
2115 struct AdditionalData
2116 : public internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>
2117 {
2118 using EigenvalueAlgorithm = internal::EigenvalueAlgorithm;
2119
2123 enum class PolynomialType
2124 {
2128 first_kind,
2133 fourth_kind
2134 };
2135
2139 AdditionalData(
2140 const unsigned int degree = 1,
2141 const double smoothing_range = 0.,
2142 const unsigned int eig_cg_n_iterations = 8,
2143 const double eig_cg_residual = 1e-2,
2144 const double max_eigenvalue = 1,
2145 const EigenvalueAlgorithm eigenvalue_algorithm =
2146 EigenvalueAlgorithm::lanczos,
2147 const PolynomialType polynomial_type = PolynomialType::first_kind);
2148
2161 unsigned int degree;
2162
2166 PolynomialType polynomial_type;
2167 };
2168
2169
2173 PreconditionChebyshev();
2174
2186 void
2187 initialize(const MatrixType &matrix,
2188 const AdditionalData &additional_data = AdditionalData());
2189
2194 void
2195 vmult(VectorType &dst, const VectorType &src) const;
2196
2201 void
2202 Tvmult(VectorType &dst, const VectorType &src) const;
2203
2207 void
2208 step(VectorType &dst, const VectorType &src) const;
2209
2213 void
2214 Tstep(VectorType &dst, const VectorType &src) const;
2215
2219 void
2220 clear();
2221
2226 size_type
2227 m() const;
2228
2233 size_type
2234 n() const;
2235
2236 using EigenvalueInformation = internal::EigenvalueInformation;
2237
2250 EigenvalueInformation
2251 estimate_eigenvalues(const VectorType &src) const;
2252
2253private:
2257 ObserverPointer<
2258 const MatrixType,
2259 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>>
2260 matrix_ptr;
2261
2265 mutable VectorType solution_old;
2266
2270 mutable VectorType temp_vector1;
2271
2275 mutable VectorType temp_vector2;
2276
2281 AdditionalData data;
2282
2286 double theta;
2287
2292 double delta;
2293
2298 bool eigenvalues_are_initialized;
2299
2304 mutable Threads::Mutex mutex;
2305};
2306
2307
2308
2310/* ---------------------------------- Inline functions ------------------- */
2311
2312#ifndef DOXYGEN
2313
2314
2315namespace internal
2316{
2317 template <typename VectorType>
2318 void
2319 set_initial_guess(VectorType &vector)
2320 {
2321 vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2322 if (vector.locally_owned_elements().is_element(0))
2323 vector(0) = 0.;
2324 }
2325
2326 template <typename Number>
2327 void
2328 set_initial_guess(::Vector<Number> &vector)
2329 {
2330 // Choose a high-frequency mode consisting of numbers between 0 and 1
2331 // that is cheap to compute (cheaper than random numbers) but avoids
2332 // obviously re-occurring numbers in multi-component systems by choosing
2333 // a period of 11
2334 for (unsigned int i = 0; i < vector.size(); ++i)
2335 vector(i) = i % 11;
2336
2337 const Number mean_value = vector.mean_value();
2338 vector.add(-mean_value);
2339 }
2340
2341 template <typename Number>
2342 void
2343 set_initial_guess(
2344 ::LinearAlgebra::distributed::BlockVector<Number> &vector)
2345 {
2346 for (unsigned int block = 0; block < vector.n_blocks(); ++block)
2347 set_initial_guess(vector.block(block));
2348 }
2349
2350 template <typename Number, typename MemorySpace>
2351 void
2352 set_initial_guess(
2353 ::LinearAlgebra::distributed::Vector<Number, MemorySpace> &vector)
2354 {
2355 // Choose a high-frequency mode consisting of numbers between 0 and 1
2356 // that is cheap to compute (cheaper than random numbers) but avoids
2357 // obviously re-occurring numbers in multi-component systems by choosing
2358 // a period of 11.
2359 // Make initial guess robust with respect to number of processors
2360 // by operating on the global index.
2361 types::global_dof_index first_local_range = 0;
2362 if (!vector.locally_owned_elements().is_empty())
2363 first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2364
2365 const auto n_local_elements = vector.locally_owned_size();
2366 Number *values_ptr = vector.get_values();
2367 Kokkos::RangePolicy<typename MemorySpace::kokkos_space::execution_space,
2368 Kokkos::IndexType<types::global_dof_index>>
2369 policy(0, n_local_elements);
2370 Kokkos::parallel_for(
2371 "::PreconditionChebyshev::set_initial_guess",
2372 policy,
2373 KOKKOS_LAMBDA(types::global_dof_index i) {
2374 values_ptr[i] = (i + first_local_range) % 11;
2375 });
2376 const Number mean_value = vector.mean_value();
2377 vector.add(-mean_value);
2378 }
2379
2380 struct EigenvalueTracker
2381 {
2382 public:
2383 void
2384 slot(const std::vector<double> &eigenvalues)
2385 {
2386 values = eigenvalues;
2387 }
2388
2389 std::vector<double> values;
2390 };
2391
2392
2393
2394 template <typename MatrixType,
2395 typename VectorType,
2396 typename PreconditionerType>
2397 double
2398 power_iteration(const MatrixType &matrix,
2399 VectorType &eigenvector,
2400 const PreconditionerType &preconditioner,
2401 const unsigned int n_iterations)
2402 {
2403 typename VectorType::value_type eigenvalue_estimate = 0.;
2404 eigenvector /= eigenvector.l2_norm();
2405 VectorType vector1, vector2;
2406 vector1.reinit(eigenvector, true);
2407 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2408 vector2.reinit(eigenvector, true);
2409 for (unsigned int i = 0; i < n_iterations; ++i)
2410 {
2411 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2412 {
2413 matrix.vmult(vector2, eigenvector);
2414 preconditioner.vmult(vector1, vector2);
2415 }
2416 else
2417 matrix.vmult(vector1, eigenvector);
2418
2419 eigenvalue_estimate = eigenvector * vector1;
2420
2421 vector1 /= vector1.l2_norm();
2422 eigenvector.swap(vector1);
2423 }
2424 return std::abs(eigenvalue_estimate);
2425 }
2426
2427
2428
2429 template <typename MatrixType,
2430 typename VectorType,
2431 typename PreconditionerType>
2432 EigenvalueInformation
2433 estimate_eigenvalues(
2434 const EigenvalueAlgorithmAdditionalData<PreconditionerType> &data,
2435 const MatrixType *matrix_ptr,
2436 VectorType &solution_old,
2437 VectorType &temp_vector1,
2438 const unsigned int degree)
2439 {
2440 Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2441
2442 EigenvalueInformation info{};
2443
2444 if (data.eig_cg_n_iterations > 0)
2445 {
2446 Assert(data.eig_cg_n_iterations > 2,
2447 ExcMessage(
2448 "Need to set at least two iterations to find eigenvalues."));
2449
2450 internal::EigenvalueTracker eigenvalue_tracker;
2451
2452 // set an initial guess that contains some high-frequency parts (to the
2453 // extent possible without knowing the discretization and the numbering)
2454 // to trigger high eigenvalues according to the external function
2455 internal::set_initial_guess(temp_vector1);
2456 data.constraints.set_zero(temp_vector1);
2457
2458 if (data.eigenvalue_algorithm == internal::EigenvalueAlgorithm::lanczos)
2459 {
2460 // set a very strict tolerance to force at least two iterations
2461 IterationNumberControl control(data.eig_cg_n_iterations,
2462 1e-10,
2463 false,
2464 false);
2465
2466 ::SolverCG<VectorType> solver(control);
2467 solver.connect_eigenvalues_slot(
2468 [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2469 eigenvalue_tracker.slot(eigenvalues);
2470 });
2471
2472 solver.solve(*matrix_ptr,
2473 solution_old,
2474 temp_vector1,
2475 *data.preconditioner);
2476
2477 info.cg_iterations = control.last_step();
2478 }
2479 else if (data.eigenvalue_algorithm ==
2480 internal::EigenvalueAlgorithm::power_iteration)
2481 {
2482 (void)degree;
2483
2484 Assert(degree != numbers::invalid_unsigned_int,
2485 ExcMessage("Cannot estimate the minimal eigenvalue with the "
2486 "power iteration"));
2487
2488 eigenvalue_tracker.values.push_back(
2489 internal::power_iteration(*matrix_ptr,
2490 temp_vector1,
2491 *data.preconditioner,
2492 data.eig_cg_n_iterations));
2493 }
2494 else
2495 DEAL_II_NOT_IMPLEMENTED();
2496
2497 // read the eigenvalues from the attached eigenvalue tracker
2498 if (eigenvalue_tracker.values.empty())
2499 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2500 else
2501 {
2502 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2503
2504 // include a safety factor since the CG method will in general not
2505 // be converged
2506 info.max_eigenvalue_estimate =
2507 1.2 * eigenvalue_tracker.values.back();
2508 }
2509 }
2510 else
2511 {
2512 info.max_eigenvalue_estimate = data.max_eigenvalue;
2513 info.min_eigenvalue_estimate =
2514 data.max_eigenvalue / data.smoothing_range;
2515 }
2516
2517 return info;
2518 }
2519} // namespace internal
2520
2521
2522
2523inline PreconditionIdentity::PreconditionIdentity()
2524 : n_rows(0)
2525 , n_columns(0)
2526{}
2527
2528template <typename MatrixType>
2529inline void
2530PreconditionIdentity::initialize(const MatrixType &matrix,
2531 const PreconditionIdentity::AdditionalData &)
2532{
2533 n_rows = matrix.m();
2534 n_columns = matrix.n();
2535}
2536
2537
2538template <typename VectorType>
2539inline void
2540PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
2541{
2542 dst = src;
2543}
2544
2545
2546
2547template <typename VectorType>
2548inline void
2549PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
2550{
2551 dst = src;
2552}
2553
2554template <typename VectorType>
2555inline void
2556PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
2557{
2558 dst += src;
2559}
2560
2561
2562
2563template <typename VectorType>
2564inline void
2565PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
2566{
2567 dst += src;
2568}
2569
2570
2571
2572inline void
2573PreconditionIdentity::clear()
2574{}
2575
2576
2577
2578inline PreconditionIdentity::size_type
2579PreconditionIdentity::m() const
2580{
2581 Assert(n_rows != 0, ExcNotInitialized());
2582 return n_rows;
2583}
2584
2585inline PreconditionIdentity::size_type
2586PreconditionIdentity::n() const
2587{
2588 Assert(n_columns != 0, ExcNotInitialized());
2589 return n_columns;
2590}
2591
2592//---------------------------------------------------------------------------
2593
2594inline PreconditionRichardson::AdditionalData::AdditionalData(
2595 const double relaxation)
2596 : relaxation(relaxation)
2597{}
2598
2599
2600inline PreconditionRichardson::PreconditionRichardson()
2601 : relaxation(0)
2602 , n_rows(0)
2603 , n_columns(0)
2604{
2605 AdditionalData add_data;
2606 relaxation = add_data.relaxation;
2607}
2608
2609
2610
2611inline void
2612PreconditionRichardson::initialize(
2613 const PreconditionRichardson::AdditionalData &parameters)
2614{
2615 relaxation = parameters.relaxation;
2616}
2617
2618
2619
2620template <typename MatrixType>
2621inline void
2622PreconditionRichardson::initialize(
2623 const MatrixType &matrix,
2624 const PreconditionRichardson::AdditionalData &parameters)
2625{
2626 relaxation = parameters.relaxation;
2627 n_rows = matrix.m();
2628 n_columns = matrix.n();
2629}
2630
2631
2632
2633template <typename VectorType>
2634inline void
2635PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2636{
2637 static_assert(
2638 std::is_same_v<size_type, typename VectorType::size_type>,
2639 "PreconditionRichardson and VectorType must have the same size_type.");
2640
2641 dst.equ(relaxation, src);
2642}
2643
2644
2645
2646template <typename VectorType>
2647inline void
2648PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2649{
2650 static_assert(
2651 std::is_same_v<size_type, typename VectorType::size_type>,
2652 "PreconditionRichardson and VectorType must have the same size_type.");
2653
2654 dst.equ(relaxation, src);
2655}
2656
2657template <typename VectorType>
2658inline void
2659PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2660{
2661 static_assert(
2662 std::is_same_v<size_type, typename VectorType::size_type>,
2663 "PreconditionRichardson and VectorType must have the same size_type.");
2664
2665 dst.add(relaxation, src);
2666}
2667
2668
2669
2670template <typename VectorType>
2671inline void
2672PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2673{
2674 static_assert(
2675 std::is_same_v<size_type, typename VectorType::size_type>,
2676 "PreconditionRichardson and VectorType must have the same size_type.");
2677
2678 dst.add(relaxation, src);
2679}
2680
2681inline PreconditionRichardson::size_type
2682PreconditionRichardson::m() const
2683{
2684 Assert(n_rows != 0, ExcNotInitialized());
2685 return n_rows;
2686}
2687
2688inline PreconditionRichardson::size_type
2689PreconditionRichardson::n() const
2690{
2691 Assert(n_columns != 0, ExcNotInitialized());
2692 return n_columns;
2693}
2694
2695//---------------------------------------------------------------------------
2696
2697template <typename MatrixType, typename PreconditionerType>
2698inline void
2699PreconditionRelaxation<MatrixType, PreconditionerType>::initialize(
2700 const MatrixType &rA,
2701 const AdditionalData &parameters)
2702{
2703 A = &rA;
2704 eigenvalues_are_initialized = false;
2705
2706 Assert(parameters.preconditioner, ExcNotInitialized());
2707
2708 this->data = parameters;
2709}
2710
2711
2712template <typename MatrixType, typename PreconditionerType>
2713inline void
2714PreconditionRelaxation<MatrixType, PreconditionerType>::clear()
2715{
2716 eigenvalues_are_initialized = false;
2717 A = nullptr;
2718 data.relaxation = 1.0;
2719 data.preconditioner = nullptr;
2720}
2721
2722template <typename MatrixType, typename PreconditionerType>
2723inline
2724 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2725 PreconditionRelaxation<MatrixType, PreconditionerType>::m() const
2726{
2727 Assert(A != nullptr, ExcNotInitialized());
2728 return A->m();
2729}
2730
2731template <typename MatrixType, typename PreconditionerType>
2732inline
2733 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2734 PreconditionRelaxation<MatrixType, PreconditionerType>::n() const
2735{
2736 Assert(A != nullptr, ExcNotInitialized());
2737 return A->n();
2738}
2739
2740template <typename MatrixType, typename PreconditionerType>
2741template <typename VectorType>
2742inline void
2743PreconditionRelaxation<MatrixType, PreconditionerType>::vmult(
2744 VectorType &dst,
2745 const VectorType &src) const
2746{
2747 Assert(this->A != nullptr, ExcNotInitialized());
2748 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2749
2750 if (eigenvalues_are_initialized == false)
2751 estimate_eigenvalues(src);
2752
2753 VectorType tmp1, tmp2;
2754
2755 for (unsigned int i = 0; i < data.n_iterations; ++i)
2756 internal::PreconditionRelaxation::step_operations(*A,
2757 *data.preconditioner,
2758 dst,
2759 src,
2760 data.relaxation,
2761 tmp1,
2762 tmp2,
2763 i,
2764 false);
2765}
2766
2767template <typename MatrixType, typename PreconditionerType>
2768template <typename VectorType>
2769inline void
2770PreconditionRelaxation<MatrixType, PreconditionerType>::Tvmult(
2771 VectorType &dst,
2772 const VectorType &src) const
2773{
2774 Assert(this->A != nullptr, ExcNotInitialized());
2775 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2776
2777 if (eigenvalues_are_initialized == false)
2778 estimate_eigenvalues(src);
2779
2780 VectorType tmp1, tmp2;
2781
2782 for (unsigned int i = 0; i < data.n_iterations; ++i)
2783 internal::PreconditionRelaxation::step_operations(
2784 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2785}
2786
2787template <typename MatrixType, typename PreconditionerType>
2788template <typename VectorType>
2789inline void
2790PreconditionRelaxation<MatrixType, PreconditionerType>::step(
2791 VectorType &dst,
2792 const VectorType &src) const
2793{
2794 Assert(this->A != nullptr, ExcNotInitialized());
2795 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2796
2797 if (eigenvalues_are_initialized == false)
2798 estimate_eigenvalues(src);
2799
2800 VectorType tmp1, tmp2;
2801
2802 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2803 internal::PreconditionRelaxation::step_operations(*A,
2804 *data.preconditioner,
2805 dst,
2806 src,
2807 data.relaxation,
2808 tmp1,
2809 tmp2,
2810 i,
2811 false);
2812}
2813
2814template <typename MatrixType, typename PreconditionerType>
2815template <typename VectorType>
2816inline void
2817PreconditionRelaxation<MatrixType, PreconditionerType>::Tstep(
2818 VectorType &dst,
2819 const VectorType &src) const
2820{
2821 Assert(this->A != nullptr, ExcNotInitialized());
2822 Assert(data.preconditioner != nullptr, ExcNotInitialized());
2823
2824 if (eigenvalues_are_initialized == false)
2825 estimate_eigenvalues(src);
2826
2827 VectorType tmp1, tmp2;
2828
2829 for (unsigned int i = 1; i <= data.n_iterations; ++i)
2830 internal::PreconditionRelaxation::step_operations(
2831 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i, true);
2832}
2833
2834template <typename MatrixType, typename PreconditionerType>
2835template <typename VectorType>
2836inline internal::EigenvalueInformation
2837PreconditionRelaxation<MatrixType, PreconditionerType>::estimate_eigenvalues(
2838 const VectorType &src) const
2839{
2840 Assert(eigenvalues_are_initialized == false, ExcInternalError());
2841
2842 EigenvalueInformation info;
2843
2844 if (data.relaxation == 0.0)
2845 {
2846 VectorType solution_old, temp_vector1;
2847
2848 solution_old.reinit(src);
2849 temp_vector1.reinit(src, true);
2850
2851 info = internal::estimate_eigenvalues<MatrixType>(
2852 data, A, solution_old, temp_vector1, data.n_iterations);
2853
2854 const double alpha =
2855 (data.smoothing_range > 1. ?
2856 info.max_eigenvalue_estimate / data.smoothing_range :
2857 std::min(0.9 * info.max_eigenvalue_estimate,
2858 info.min_eigenvalue_estimate));
2859
2860 const_cast<PreconditionRelaxation<MatrixType, PreconditionerType> *>(this)
2861 ->data.relaxation = 2.0 / (alpha + info.max_eigenvalue_estimate);
2862 }
2863
2864 const_cast<PreconditionRelaxation<MatrixType, PreconditionerType> *>(this)
2865 ->eigenvalues_are_initialized = true;
2866
2867 return info;
2868}
2869
2870template <typename MatrixType, typename PreconditionerType>
2871double
2872PreconditionRelaxation<MatrixType, PreconditionerType>::get_relaxation() const
2873{
2874 return data.relaxation;
2875}
2876
2877
2878//---------------------------------------------------------------------------
2879
2880template <typename MatrixType>
2881inline void
2882PreconditionJacobi<MatrixType>::initialize(const MatrixType &A,
2883 const AdditionalData &parameters_in)
2884{
2885 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2886 Assert(
2887 parameters_in.relaxation != 0.0,
2888 ExcMessage(
2889 "Relaxation cannot automatically be determined by PreconditionJacobi."));
2890
2891 AdditionalData parameters;
2892 parameters.relaxation = 1.0;
2893 parameters.n_iterations = parameters_in.n_iterations;
2894 parameters.preconditioner =
2895 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2896
2897 this->BaseClass::initialize(A, parameters);
2898}
2899
2900//---------------------------------------------------------------------------
2901
2902template <typename MatrixType>
2903inline void
2904PreconditionSOR<MatrixType>::initialize(const MatrixType &A,
2905 const AdditionalData &parameters_in)
2906{
2907 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2908 Assert(
2909 parameters_in.relaxation != 0.0,
2910 ExcMessage(
2911 "Relaxation cannot automatically be determined by PreconditionSOR."));
2912
2913 AdditionalData parameters;
2914 parameters.relaxation = 1.0;
2915 parameters.n_iterations = parameters_in.n_iterations;
2916 parameters.preconditioner =
2917 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2918
2919 this->BaseClass::initialize(A, parameters);
2920}
2921
2922//---------------------------------------------------------------------------
2923
2924template <typename MatrixType>
2925inline void
2926PreconditionSSOR<MatrixType>::initialize(const MatrixType &A,
2927 const AdditionalData &parameters_in)
2928{
2929 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2930 Assert(
2931 parameters_in.relaxation != 0.0,
2932 ExcMessage(
2933 "Relaxation cannot automatically be determined by PreconditionSSOR."));
2934
2935 AdditionalData parameters;
2936 parameters.relaxation = 1.0;
2937 parameters.n_iterations = parameters_in.n_iterations;
2938 parameters.preconditioner =
2939 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2940
2941 this->BaseClass::initialize(A, parameters);
2942}
2943
2944
2945
2946//---------------------------------------------------------------------------
2947
2948template <typename MatrixType>
2949inline void
2950PreconditionPSOR<MatrixType>::initialize(
2951 const MatrixType &A,
2952 const std::vector<size_type> &p,
2953 const std::vector<size_type> &ip,
2954 const typename BaseClass::AdditionalData &parameters_in)
2955{
2956 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2957 Assert(
2958 parameters_in.relaxation != 0.0,
2959 ExcMessage(
2960 "Relaxation cannot automatically be determined by PreconditionPSOR."));
2961
2962 typename BaseClass::AdditionalData parameters;
2963 parameters.relaxation = 1.0;
2964 parameters.n_iterations = parameters_in.n_iterations;
2965 parameters.preconditioner =
2966 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2967
2968 this->BaseClass::initialize(A, parameters);
2969}
2970
2971
2972template <typename MatrixType>
2973inline void
2974PreconditionPSOR<MatrixType>::initialize(const MatrixType &A,
2975 const AdditionalData &additional_data)
2976{
2977 initialize(A,
2978 additional_data.permutation,
2979 additional_data.inverse_permutation,
2980 additional_data.parameters);
2981}
2982
2983template <typename MatrixType>
2984PreconditionPSOR<MatrixType>::AdditionalData::AdditionalData(
2985 const std::vector<size_type> &permutation,
2986 const std::vector<size_type> &inverse_permutation,
2987 const typename PreconditionPSOR<MatrixType>::BaseClass::AdditionalData
2988 &parameters)
2989 : permutation(permutation)
2990 , inverse_permutation(inverse_permutation)
2991 , parameters(parameters)
2992{}
2993
2994
2995//---------------------------------------------------------------------------
2996
2997
2998template <typename MatrixType, typename VectorType>
2999PreconditionUseMatrix<MatrixType, VectorType>::PreconditionUseMatrix(
3000 const MatrixType &M,
3001 const function_ptr method)
3002 : matrix(M)
3003 , precondition(method)
3004{}
3005
3006
3007
3008template <typename MatrixType, typename VectorType>
3009void
3010PreconditionUseMatrix<MatrixType, VectorType>::vmult(
3011 VectorType &dst,
3012 const VectorType &src) const
3013{
3014 (matrix.*precondition)(dst, src);
3015}
3016
3017//---------------------------------------------------------------------------
3018
3019namespace internal
3020{
3021
3022 template <typename PreconditionerType>
3023 inline EigenvalueAlgorithmAdditionalData<PreconditionerType>::
3024 EigenvalueAlgorithmAdditionalData(
3025 const double smoothing_range,
3026 const unsigned int eig_cg_n_iterations,
3027 const double eig_cg_residual,
3028 const double max_eigenvalue,
3029 const EigenvalueAlgorithm eigenvalue_algorithm)
3030 : smoothing_range(smoothing_range)
3031 , eig_cg_n_iterations(eig_cg_n_iterations)
3032 , eig_cg_residual(eig_cg_residual)
3033 , max_eigenvalue(max_eigenvalue)
3034 , eigenvalue_algorithm(eigenvalue_algorithm)
3035 {}
3036
3037
3038
3039 template <typename PreconditionerType>
3040 inline EigenvalueAlgorithmAdditionalData<PreconditionerType> &
3041 EigenvalueAlgorithmAdditionalData<PreconditionerType>::operator=(
3042 const EigenvalueAlgorithmAdditionalData &other_data)
3043 {
3044 smoothing_range = other_data.smoothing_range;
3045 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3046 eig_cg_residual = other_data.eig_cg_residual;
3047 max_eigenvalue = other_data.max_eigenvalue;
3048 preconditioner = other_data.preconditioner;
3049 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3050 constraints.copy_from(other_data.constraints);
3051
3052 return *this;
3053 }
3054} // namespace internal
3055
3056template <typename MatrixType, typename PreconditionerType>
3057inline PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData::
3058 AdditionalData(const double relaxation,
3059 const unsigned int n_iterations,
3060 const double smoothing_range,
3061 const unsigned int eig_cg_n_iterations,
3062 const double eig_cg_residual,
3063 const double max_eigenvalue,
3064 const EigenvalueAlgorithm eigenvalue_algorithm)
3065 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3066 smoothing_range,
3067 eig_cg_n_iterations,
3068 eig_cg_residual,
3069 max_eigenvalue,
3070 eigenvalue_algorithm)
3071 , relaxation(relaxation)
3072 , n_iterations(n_iterations)
3073{}
3074
3075
3076
3077//---------------------------------------------------------------------------
3078
3079namespace internal
3080{
3081 namespace PreconditionChebyshevImplementation
3082 {
3083 // for deal.II vectors, perform updates for Chebyshev preconditioner all
3084 // at once to reduce memory transfer. Here, we select between general
3085 // vectors and deal.II vectors where we expand the loop over the (local)
3086 // size of the vector
3087
3088 // generic part for non-deal.II vectors
3089 template <typename VectorType, typename PreconditionerType>
3090 inline void
3091 vector_updates(const VectorType &rhs,
3092 const PreconditionerType &preconditioner,
3093 const unsigned int iteration_index,
3094 const double factor1,
3095 const double factor2,
3096 VectorType &solution_old,
3097 VectorType &temp_vector1,
3098 VectorType &temp_vector2,
3099 VectorType &solution)
3100 {
3101 if (iteration_index == 0)
3102 {
3103 solution.equ(factor2, rhs);
3104 preconditioner.vmult(solution_old, solution);
3105 }
3106 else if (iteration_index == 1)
3107 {
3108 // compute t = P^{-1} * (b-A*x^{n})
3109 temp_vector1.sadd(-1.0, 1.0, rhs);
3110 preconditioner.vmult(solution_old, temp_vector1);
3111
3112 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
3113 solution_old.sadd(factor2, 1 + factor1, solution);
3114 }
3115 else
3116 {
3117 // compute t = P^{-1} * (b-A*x^{n})
3118 temp_vector1.sadd(-1.0, 1.0, rhs);
3119 preconditioner.vmult(temp_vector2, temp_vector1);
3120
3121 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
3122 solution_old.sadd(-factor1, factor2, temp_vector2);
3123 solution_old.add(1 + factor1, solution);
3124 }
3125
3126 solution.swap(solution_old);
3127 }
3128
3129 // generic part for deal.II vectors
3130 template <
3131 typename Number,
3132 typename PreconditionerType,
3133 std::enable_if_t<
3134 !has_vmult_with_std_functions_for_precondition<
3135 PreconditionerType,
3136 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>>,
3137 int> * = nullptr>
3138 inline void
3139 vector_updates(
3140 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3141 const PreconditionerType &preconditioner,
3142 const unsigned int iteration_index,
3143 const double factor1_,
3144 const double factor2_,
3145 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3146 &solution_old,
3147 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3148 &temp_vector1,
3149 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3150 &temp_vector2,
3151 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
3152 {
3153 const Number factor1 = factor1_;
3154 const Number factor1_plus_1 = 1. + factor1_;
3155 const Number factor2 = factor2_;
3156
3157 if (iteration_index == 0)
3158 {
3159 const auto solution_old_ptr = solution_old.begin();
3160
3161 // compute t = P^{-1} * (b)
3162 preconditioner.vmult(solution_old, rhs);
3163
3164 // compute x^{n+1} = f_2 * t
3165 DEAL_II_OPENMP_SIMD_PRAGMA
3166 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3167 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
3168 }
3169 else if (iteration_index == 1)
3170 {
3171 const auto solution_ptr = solution.begin();
3172 const auto solution_old_ptr = solution_old.begin();
3173
3174 // compute t = P^{-1} * (b-A*x^{n})
3175 temp_vector1.sadd(-1.0, 1.0, rhs);
3176
3177 preconditioner.vmult(solution_old, temp_vector1);
3178
3179 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
3180 DEAL_II_OPENMP_SIMD_PRAGMA
3181 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3182 solution_old_ptr[i] =
3183 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
3184 }
3185 else
3186 {
3187 const auto solution_ptr = solution.begin();
3188 const auto solution_old_ptr = solution_old.begin();
3189 const auto temp_vector2_ptr = temp_vector2.begin();
3190
3191 // compute t = P^{-1} * (b-A*x^{n})
3192 temp_vector1.sadd(-1.0, 1.0, rhs);
3193
3194 preconditioner.vmult(temp_vector2, temp_vector1);
3195
3196 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
3197 DEAL_II_OPENMP_SIMD_PRAGMA
3198 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
3199 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3200 factor1 * solution_old_ptr[i] +
3201 temp_vector2_ptr[i] * factor2;
3202 }
3203
3204 solution.swap(solution_old);
3205 }
3206
3207 template <
3208 typename Number,
3209 typename PreconditionerType,
3210 std::enable_if_t<
3211 has_vmult_with_std_functions_for_precondition<
3212 PreconditionerType,
3213 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>>,
3214 int> * = nullptr>
3215 inline void
3216 vector_updates(
3217 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3218 const PreconditionerType &preconditioner,
3219 const unsigned int iteration_index,
3220 const double factor1_,
3221 const double factor2_,
3222 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3223 &solution_old,
3224 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3225 &temp_vector1,
3226 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3227 &temp_vector2,
3228 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
3229 {
3230 const Number factor1 = factor1_;
3231 const Number factor1_plus_1 = 1. + factor1_;
3232 const Number factor2 = factor2_;
3233
3234 const auto rhs_ptr = rhs.begin();
3235 const auto temp_vector1_ptr = temp_vector1.begin();
3236 const auto temp_vector2_ptr = temp_vector2.begin();
3237 const auto solution_ptr = solution.begin();
3238 const auto solution_old_ptr = solution_old.begin();
3239
3240 if (iteration_index == 0)
3241 {
3242 preconditioner.vmult(
3243 solution,
3244 rhs,
3245 [&](const auto start_range, const auto end_range) {
3246 if (end_range > start_range)
3247 std::memset(solution.begin() + start_range,
3248 0,
3249 sizeof(Number) * (end_range - start_range));
3250 },
3251 [&](const auto begin, const auto end) {
3252 DEAL_II_OPENMP_SIMD_PRAGMA
3253 for (std::size_t i = begin; i < end; ++i)
3254 solution_ptr[i] *= factor2;
3255 });
3256 }
3257 else
3258 {
3259 preconditioner.vmult(
3260 temp_vector2,
3261 temp_vector1,
3262 [&](const auto begin, const auto end) {
3263 if (end > begin)
3264 std::memset(temp_vector2.begin() + begin,
3265 0,
3266 sizeof(Number) * (end - begin));
3267
3268 DEAL_II_OPENMP_SIMD_PRAGMA
3269 for (std::size_t i = begin; i < end; ++i)
3270 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3271 },
3272 [&](const auto begin, const auto end) {
3273 if (iteration_index == 1)
3274 {
3275 DEAL_II_OPENMP_SIMD_PRAGMA
3276 for (std::size_t i = begin; i < end; ++i)
3277 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
3278 factor2 * temp_vector2_ptr[i];
3279 }
3280 else
3281 {
3282 DEAL_II_OPENMP_SIMD_PRAGMA
3283 for (std::size_t i = begin; i < end; ++i)
3284 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3285 factor1 * solution_old_ptr[i] +
3286 factor2 * temp_vector2_ptr[i];
3287 }
3288 });
3289 }
3290
3291 if (iteration_index > 0)
3292 {
3293 solution_old.swap(temp_vector2);
3294 solution_old.swap(solution);
3295 }
3296 }
3297
3298 // worker routine for deal.II vectors. Because of vectorization, we need
3299 // to put the loop into an extra structure because the virtual function of
3300 // VectorUpdatesRange prevents the compiler from applying vectorization.
3301 template <typename Number>
3302 struct VectorUpdater
3303 {
3304 VectorUpdater(const Number *rhs,
3305 const Number *matrix_diagonal_inverse,
3306 const unsigned int iteration_index,
3307 const Number factor1,
3308 const Number factor2,
3309 Number *solution_old,
3310 Number *tmp_vector,
3311 Number *solution)
3312 : rhs(rhs)
3313 , matrix_diagonal_inverse(matrix_diagonal_inverse)
3314 , iteration_index(iteration_index)
3315 , factor1(factor1)
3316 , factor2(factor2)
3317 , solution_old(solution_old)
3318 , tmp_vector(tmp_vector)
3319 , solution(solution)
3320 {}
3321
3322 void
3323 apply_to_subrange(const std::size_t begin, const std::size_t end) const
3324 {
3325 // To circumvent a bug in gcc
3326 // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
3327 // copies of the variables factor1 and factor2 and do not check based on
3328 // factor1.
3329 const Number factor1 = this->factor1;
3330 const Number factor1_plus_1 = 1. + this->factor1;
3331 const Number factor2 = this->factor2;
3332 if (iteration_index == 0)
3333 {
3334 DEAL_II_OPENMP_SIMD_PRAGMA
3335 for (std::size_t i = begin; i < end; ++i)
3336 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
3337 }
3338 else if (iteration_index == 1)
3339 {
3340 // x^{n+1} = x^{n} + f_1 * x^{n} + 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 that is
3344 // already read (avoid read-for-ownership)
3345 tmp_vector[i] =
3346 factor1_plus_1 * solution[i] +
3347 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3348 }
3349 else
3350 {
3351 // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
3352 // + f_2 * P^{-1} * (b-A*x^{n})
3353 DEAL_II_OPENMP_SIMD_PRAGMA
3354 for (std::size_t i = begin; i < end; ++i)
3355 // for efficiency reason, write back to temp_vector, which is
3356 // already modified during vmult (in best case, the modified
3357 // values are not written back to main memory yet so that
3358 // we do not have to pay additional costs for writing and
3359 // read-for-ownershop)
3360 tmp_vector[i] =
3361 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
3362 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3363 }
3364 }
3365
3366 const Number *rhs;
3367 const Number *matrix_diagonal_inverse;
3368 const unsigned int iteration_index;
3369 const Number factor1;
3370 const Number factor2;
3371 mutable Number *solution_old;
3372 mutable Number *tmp_vector;
3373 mutable Number *solution;
3374 };
3375
3376 template <typename Number>
3377 struct VectorUpdatesRange : public ::parallel::ParallelForInteger
3378 {
3379 VectorUpdatesRange(const VectorUpdater<Number> &updater,
3380 const std::size_t size)
3381 : updater(updater)
3382 {
3383 if (size < internal::VectorImplementation::minimum_parallel_grain_size)
3384 VectorUpdatesRange::apply_to_subrange(0, size);
3385 else
3386 apply_parallel(
3387 0,
3388 size,
3389 internal::VectorImplementation::minimum_parallel_grain_size);
3390 }
3391
3392 ~VectorUpdatesRange() override = default;
3393
3394 virtual void
3395 apply_to_subrange(const std::size_t begin,
3396 const std::size_t end) const override
3397 {
3398 updater.apply_to_subrange(begin, end);
3399 }
3400
3401 const VectorUpdater<Number> &updater;
3402 };
3403
3404 // selection for diagonal matrix around deal.II vector
3405 template <typename Number>
3406 inline void
3407 vector_updates(
3408 const ::Vector<Number> &rhs,
3409 const ::DiagonalMatrix<::Vector<Number>> &jacobi,
3410 const unsigned int iteration_index,
3411 const double factor1,
3412 const double factor2,
3413 ::Vector<Number> &solution_old,
3414 ::Vector<Number> &temp_vector1,
3415 ::Vector<Number> &,
3416 ::Vector<Number> &solution)
3417 {
3418 VectorUpdater<Number> upd(rhs.begin(),
3419 jacobi.get_vector().begin(),
3420 iteration_index,
3421 factor1,
3422 factor2,
3423 solution_old.begin(),
3424 temp_vector1.begin(),
3425 solution.begin());
3426 VectorUpdatesRange<Number>(upd, rhs.size());
3427
3428 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3429 if (iteration_index == 0)
3430 {
3431 // nothing to do here because we can immediately write into the
3432 // solution vector without remembering any of the other vectors
3433 }
3434 else
3435 {
3436 solution.swap(temp_vector1);
3437 solution_old.swap(temp_vector1);
3438 }
3439 }
3440
3441 // selection for diagonal matrix around parallel deal.II vector
3442 template <typename Number>
3443 inline void
3444 vector_updates(
3445 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3446 const ::DiagonalMatrix<
3447 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>> &jacobi,
3448 const unsigned int iteration_index,
3449 const double factor1,
3450 const double factor2,
3451 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3452 &solution_old,
3453 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3454 &temp_vector1,
3455 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &,
3456 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
3457 {
3458 VectorUpdater<Number> upd(rhs.begin(),
3459 jacobi.get_vector().begin(),
3460 iteration_index,
3461 factor1,
3462 factor2,
3463 solution_old.begin(),
3464 temp_vector1.begin(),
3465 solution.begin());
3466 VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
3467
3468 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3469 if (iteration_index == 0)
3470 {
3471 // nothing to do here because we can immediately write into the
3472 // solution vector without remembering any of the other vectors
3473 }
3474 else
3475 {
3476 solution.swap(temp_vector1);
3477 solution_old.swap(temp_vector1);
3478 }
3479 }
3480
3481 // We need to have a separate declaration for static const members
3482
3483 // general case and the case that the preconditioner can work on
3484 // ranges (covered by vector_updates())
3485 template <
3486 typename MatrixType,
3487 typename VectorType,
3488 typename PreconditionerType,
3489 std::enable_if_t<
3490 !has_vmult_with_std_functions<MatrixType,
3491 VectorType,
3492 PreconditionerType> &&
3493 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3494 VectorType> &&
3495 has_vmult_with_std_functions_for_precondition<MatrixType,
3496 VectorType>),
3497 int> * = nullptr>
3498 inline void
3499 vmult_and_update(const MatrixType &matrix,
3500 const PreconditionerType &preconditioner,
3501 const VectorType &rhs,
3502 const unsigned int iteration_index,
3503 const double factor1,
3504 const double factor2,
3505 VectorType &solution,
3506 VectorType &solution_old,
3507 VectorType &temp_vector1,
3508 VectorType &temp_vector2)
3509 {
3510 if (iteration_index > 0)
3511 matrix.vmult(temp_vector1, solution);
3512 vector_updates(rhs,
3513 preconditioner,
3514 iteration_index,
3515 factor1,
3516 factor2,
3517 solution_old,
3518 temp_vector1,
3519 temp_vector2,
3520 solution);
3521 }
3522
3523 // case that both the operator and the preconditioner can work on
3524 // subranges
3525 template <
3526 typename MatrixType,
3527 typename VectorType,
3528 typename PreconditionerType,
3529 std::enable_if_t<
3530 !has_vmult_with_std_functions<MatrixType,
3531 VectorType,
3532 PreconditionerType> &&
3533 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3534 VectorType> &&
3535 has_vmult_with_std_functions_for_precondition<MatrixType,
3536 VectorType>),
3537 int> * = nullptr>
3538 inline void
3539 vmult_and_update(const MatrixType &matrix,
3540 const PreconditionerType &preconditioner,
3541 const VectorType &rhs,
3542 const unsigned int iteration_index,
3543 const double factor1_,
3544 const double factor2_,
3545 VectorType &solution,
3546 VectorType &solution_old,
3547 VectorType &temp_vector1,
3548 VectorType &temp_vector2)
3549 {
3550 using Number = typename VectorType::value_type;
3551
3552 const Number factor1 = factor1_;
3553 const Number factor1_plus_1 = 1. + factor1_;
3554 const Number factor2 = factor2_;
3555
3556 if (iteration_index == 0)
3557 {
3558 preconditioner.vmult(
3559 solution,
3560 rhs,
3561 [&](const unsigned int start_range, const unsigned int end_range) {
3562 // zero 'solution' before running the vmult operation
3563 if (end_range > start_range)
3564 std::memset(solution.begin() + start_range,
3565 0,
3566 sizeof(Number) * (end_range - start_range));
3567 },
3568 [&](const unsigned int start_range, const unsigned int end_range) {
3569 const auto solution_ptr = solution.begin();
3570
3571 DEAL_II_OPENMP_SIMD_PRAGMA
3572 for (std::size_t i = start_range; i < end_range; ++i)
3573 solution_ptr[i] *= factor2;
3574 });
3575 }
3576 else
3577 {
3578 temp_vector1.reinit(rhs, true);
3579 temp_vector2.reinit(rhs, true);
3580
3581 // 1) compute residual (including operator application)
3582 matrix.vmult(
3583 temp_vector1,
3584 solution,
3585 [&](const unsigned int start_range, const unsigned int end_range) {
3586 // zero 'temp_vector1' before running the vmult
3587 // operation
3588 if (end_range > start_range)
3589 std::memset(temp_vector1.begin() + start_range,
3590 0,
3591 sizeof(Number) * (end_range - start_range));
3592 },
3593 [&](const unsigned int start_range, const unsigned int end_range) {
3594 const auto rhs_ptr = rhs.begin();
3595 const auto tmp_ptr = temp_vector1.begin();
3596
3597 DEAL_II_OPENMP_SIMD_PRAGMA
3598 for (std::size_t i = start_range; i < end_range; ++i)
3599 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3600 });
3601
3602 // 2) perform vector updates (including preconditioner application)
3603 preconditioner.vmult(
3604 temp_vector2,
3605 temp_vector1,
3606 [&](const unsigned int start_range, const unsigned int end_range) {
3607 // zero 'temp_vector2' before running the vmult
3608 // operation
3609 if (end_range > start_range)
3610 std::memset(temp_vector2.begin() + start_range,
3611 0,
3612 sizeof(Number) * (end_range - start_range));
3613 },
3614 [&](const unsigned int start_range, const unsigned int end_range) {
3615 const auto solution_ptr = solution.begin();
3616 const auto solution_old_ptr = solution_old.begin();
3617 const auto tmp_ptr = temp_vector2.begin();
3618
3619 if (iteration_index == 1)
3620 {
3621 DEAL_II_OPENMP_SIMD_PRAGMA
3622 for (std::size_t i = start_range; i < end_range; ++i)
3623 tmp_ptr[i] =
3624 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3625 }
3626 else
3627 {
3628 DEAL_II_OPENMP_SIMD_PRAGMA
3629 for (std::size_t i = start_range; i < end_range; ++i)
3630 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3631 factor1 * solution_old_ptr[i] +
3632 factor2 * tmp_ptr[i];
3633 }
3634 });
3635
3636 solution.swap(temp_vector2);
3637 solution_old.swap(temp_vector2);
3638 }
3639 }
3640
3641 // case that the operator can work on subranges and the preconditioner
3642 // is a diagonal
3643 template <typename MatrixType,
3644 typename VectorType,
3645 typename PreconditionerType,
3646 std::enable_if_t<has_vmult_with_std_functions<MatrixType,
3647 VectorType,
3648 PreconditionerType>,
3649 int> * = nullptr>
3650 inline void
3651 vmult_and_update(const MatrixType &matrix,
3652 const PreconditionerType &preconditioner,
3653 const VectorType &rhs,
3654 const unsigned int iteration_index,
3655 const double factor1,
3656 const double factor2,
3657 VectorType &solution,
3658 VectorType &solution_old,
3659 VectorType &temp_vector1,
3660 VectorType &)
3661 {
3662 using Number = typename VectorType::value_type;
3663 VectorUpdater<Number> updater(rhs.begin(),
3664 preconditioner.get_vector().begin(),
3665 iteration_index,
3666 factor1,
3667 factor2,
3668 solution_old.begin(),
3669 temp_vector1.begin(),
3670 solution.begin());
3671 if (iteration_index > 0)
3672 matrix.vmult(
3673 temp_vector1,
3674 solution,
3675 [&](const unsigned int start_range, const unsigned int end_range) {
3676 // zero 'temp_vector1' before running the vmult
3677 // operation
3678 if (end_range > start_range)
3679 std::memset(temp_vector1.begin() + start_range,
3680 0,
3681 sizeof(Number) * (end_range - start_range));
3682 },
3683 [&](const unsigned int start_range, const unsigned int end_range) {
3684 if (end_range > start_range)
3685 updater.apply_to_subrange(start_range, end_range);
3686 });
3687 else
3688 updater.apply_to_subrange(0U, solution.locally_owned_size());
3689
3690 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3691 if (iteration_index == 0)
3692 {
3693 // nothing to do here because we can immediately write into the
3694 // solution vector without remembering any of the other vectors
3695 }
3696 else
3697 {
3698 solution.swap(temp_vector1);
3699 solution_old.swap(temp_vector1);
3700 }
3701 }
3702
3703 template <typename MatrixType, typename PreconditionerType>
3704 inline void
3705 initialize_preconditioner(
3706 const MatrixType &matrix,
3707 std::shared_ptr<PreconditionerType> &preconditioner)
3708 {
3709 (void)matrix;
3710 (void)preconditioner;
3711 AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
3712 }
3713
3714 template <typename MatrixType, typename VectorType>
3715 inline void
3716 initialize_preconditioner(
3717 const MatrixType &matrix,
3718 std::shared_ptr<::DiagonalMatrix<VectorType>> &preconditioner)
3719 {
3720 if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
3721 {
3722 if (preconditioner.get() == nullptr)
3723 preconditioner =
3724 std::make_shared<::DiagonalMatrix<VectorType>>();
3725
3726 Assert(
3727 preconditioner->m() == 0,
3728 ExcMessage(
3729 "Preconditioner appears to be initialized but not sized correctly"));
3730
3731 // This part only works in serial
3732 if (preconditioner->m() != matrix.m())
3733 {
3734 preconditioner->get_vector().reinit(matrix.m());
3735 for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
3736 preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
3737 }
3738 }
3739 }
3740 } // namespace PreconditionChebyshevImplementation
3741} // namespace internal
3742
3743
3744
3745template <typename MatrixType, typename VectorType, typename PreconditionerType>
3746inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3747 AdditionalData::AdditionalData(const unsigned int degree,
3748 const double smoothing_range,
3749 const unsigned int eig_cg_n_iterations,
3750 const double eig_cg_residual,
3751 const double max_eigenvalue,
3752 const EigenvalueAlgorithm eigenvalue_algorithm,
3753 const PolynomialType polynomial_type)
3754 : internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3755 smoothing_range,
3756 eig_cg_n_iterations,
3757 eig_cg_residual,
3758 max_eigenvalue,
3759 eigenvalue_algorithm)
3760 , degree(degree)
3761 , polynomial_type(polynomial_type)
3762{}
3763
3764
3765
3766template <typename MatrixType, typename VectorType, typename PreconditionerType>
3767inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3768 PreconditionChebyshev()
3769 : theta(1.)
3770 , delta(1.)
3771 , eigenvalues_are_initialized(false)
3772{
3773 static_assert(
3774 std::is_same_v<size_type, typename VectorType::size_type>,
3775 "PreconditionChebyshev and VectorType must have the same size_type.");
3776}
3777
3778
3779
3780template <typename MatrixType, typename VectorType, typename PreconditionerType>
3781inline void
3782PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::initialize(
3783 const MatrixType &matrix,
3784 const AdditionalData &additional_data)
3785{
3786 matrix_ptr = &matrix;
3787 data = additional_data;
3788 Assert(data.degree > 0,
3789 ExcMessage("The degree of the Chebyshev method must be positive."));
3790 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3791 matrix, data.preconditioner);
3792 eigenvalues_are_initialized = false;
3793}
3794
3795
3796
3797template <typename MatrixType, typename VectorType, typename PreconditionerType>
3798inline void
3799PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::clear()
3800{
3801 eigenvalues_are_initialized = false;
3802 theta = delta = 1.0;
3803 matrix_ptr = nullptr;
3804 {
3805 VectorType empty_vector;
3806 solution_old.reinit(empty_vector);
3807 temp_vector1.reinit(empty_vector);
3808 temp_vector2.reinit(empty_vector);
3809 }
3810 data.preconditioner.reset();
3811}
3812
3813
3814
3815template <typename MatrixType, typename VectorType, typename PreconditionerType>
3816inline typename internal::EigenvalueInformation
3817PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3818 estimate_eigenvalues(const VectorType &src) const
3819{
3820 Assert(eigenvalues_are_initialized == false, ExcInternalError());
3821
3822 solution_old.reinit(src);
3823 temp_vector1.reinit(src, true);
3824
3825 auto info = internal::estimate_eigenvalues<MatrixType>(
3826 data, matrix_ptr, solution_old, temp_vector1, data.degree);
3827
3828 const double alpha = (data.smoothing_range > 1. ?
3829 info.max_eigenvalue_estimate / data.smoothing_range :
3830 std::min(0.9 * info.max_eigenvalue_estimate,
3831 info.min_eigenvalue_estimate));
3832
3833 // in case the user set the degree to invalid unsigned int, we have to
3834 // determine the number of necessary iterations from the Chebyshev error
3835 // estimate, given the target tolerance specified by smoothing_range. This
3836 // estimate is based on the error formula given in section 5.1 of
3837 // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
3838 if (data.degree == numbers::invalid_unsigned_int)
3839 {
3840 const double actual_range = info.max_eigenvalue_estimate / alpha;
3841 const double sigma = (1. - std::sqrt(1. / actual_range)) /
3842 (1. + std::sqrt(1. / actual_range));
3843 const double eps = data.smoothing_range;
3844 const_cast<
3845 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(
3846 this)
3847 ->data.degree =
3848 1 + static_cast<unsigned int>(
3849 std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
3850 std::log(1. / sigma));
3851 }
3852
3853 info.degree = data.degree;
3854
3855 const_cast<
3856 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3857 ->delta =
3858 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3859 (info.max_eigenvalue_estimate) :
3860 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3861 const_cast<
3862 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3863 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3864
3865 // We do not need the second temporary vector in case we have a
3866 // DiagonalMatrix as preconditioner and use deal.II's own vectors
3867 using NumberType = typename VectorType::value_type;
3868 if (std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> ==
3869 false ||
3870 (std::is_same_v<VectorType, ::Vector<NumberType>> == false &&
3871 ((std::is_same_v<
3872 VectorType,
3873 LinearAlgebra::distributed::Vector<NumberType, MemorySpace::Host>> ==
3874 false) ||
3875 (std::is_same_v<VectorType,
3876 LinearAlgebra::distributed::
3877 Vector<NumberType, MemorySpace::Default>> == false))))
3878 temp_vector2.reinit(src, true);
3879 else
3880 {
3881 VectorType empty_vector;
3882 temp_vector2.reinit(empty_vector);
3883 }
3884
3885 const_cast<
3886 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3887 ->eigenvalues_are_initialized = true;
3888
3889 return info;
3890}
3891
3892
3893
3894template <typename MatrixType, typename VectorType, typename PreconditionerType>
3895inline void
3896PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::vmult(
3897 VectorType &solution,
3898 const VectorType &rhs) const
3899{
3900 std::lock_guard<std::mutex> lock(mutex);
3901 if (eigenvalues_are_initialized == false)
3902 estimate_eigenvalues(rhs);
3903
3904 internal::PreconditionChebyshevImplementation::vmult_and_update(
3905 *matrix_ptr,
3906 *data.preconditioner,
3907 rhs,
3908 0,
3909 0.,
3910 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3911 (4. / (3. * delta)) :
3912 (1. / theta),
3913 solution,
3914 solution_old,
3915 temp_vector1,
3916 temp_vector2);
3917
3918 // if delta is zero, we do not need to iterate because the updates will be
3919 // zero
3920 if (data.degree < 2 || std::abs(delta) < 1e-40)
3921 return;
3922
3923 double rhok = delta / theta, sigma = theta / delta;
3924 for (unsigned int k = 0; k < data.degree - 1; ++k)
3925 {
3926 double factor1 = 0.0;
3927 double factor2 = 0.0;
3928
3929 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3930 {
3931 factor1 = (2 * k + 1.) / (2 * k + 5.);
3932 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3933 }
3934 else
3935 {
3936 const double rhokp = 1. / (2. * sigma - rhok);
3937 factor1 = rhokp * rhok;
3938 factor2 = 2. * rhokp / delta;
3939 rhok = rhokp;
3940 }
3941
3942 internal::PreconditionChebyshevImplementation::vmult_and_update(
3943 *matrix_ptr,
3944 *data.preconditioner,
3945 rhs,
3946 k + 1,
3947 factor1,
3948 factor2,
3949 solution,
3950 solution_old,
3951 temp_vector1,
3952 temp_vector2);
3953 }
3954}
3955
3956
3957
3958template <typename MatrixType, typename VectorType, typename PreconditionerType>
3959inline void
3960PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tvmult(
3961 VectorType &solution,
3962 const VectorType &rhs) const
3963{
3964 std::lock_guard<std::mutex> lock(mutex);
3965 if (eigenvalues_are_initialized == false)
3966 estimate_eigenvalues(rhs);
3967
3968 internal::PreconditionChebyshevImplementation::vector_updates(
3969 rhs,
3970 *data.preconditioner,
3971 0,
3972 0.,
3973 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3974 (4. / (3. * delta)) :
3975 (1. / theta),
3976 solution_old,
3977 temp_vector1,
3978 temp_vector2,
3979 solution);
3980
3981 if (data.degree < 2 || std::abs(delta) < 1e-40)
3982 return;
3983
3984 double rhok = delta / theta, sigma = theta / delta;
3985 for (unsigned int k = 0; k < data.degree - 1; ++k)
3986 {
3987 double factor1 = 0.0;
3988 double factor2 = 0.0;
3989
3990 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3991 {
3992 factor1 = (2 * k + 1.) / (2 * k + 5.);
3993 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3994 }
3995 else
3996 {
3997 const double rhokp = 1. / (2. * sigma - rhok);
3998 factor1 = rhokp * rhok;
3999 factor2 = 2. * rhokp / delta;
4000 rhok = rhokp;
4001 }
4002
4003 matrix_ptr->Tvmult(temp_vector1, solution);
4004 internal::PreconditionChebyshevImplementation::vector_updates(
4005 rhs,
4006 *data.preconditioner,
4007 k + 1,
4008 factor1,
4009 factor2,
4010 solution_old,
4011 temp_vector1,
4012 temp_vector2,
4013 solution);
4014 }
4015}
4016
4017
4018
4019template <typename MatrixType, typename VectorType, typename PreconditionerType>
4020inline void
4021PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::step(
4022 VectorType &solution,
4023 const VectorType &rhs) const
4024{
4025 std::lock_guard<std::mutex> lock(mutex);
4026 if (eigenvalues_are_initialized == false)
4027 estimate_eigenvalues(rhs);
4028
4029 internal::PreconditionChebyshevImplementation::vmult_and_update(
4030 *matrix_ptr,
4031 *data.preconditioner,
4032 rhs,
4033 1,
4034 0.,
4035 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4036 (4. / (3. * delta)) :
4037 (1. / theta),
4038 solution,
4039 solution_old,
4040 temp_vector1,
4041 temp_vector2);
4042
4043 if (data.degree < 2 || std::abs(delta) < 1e-40)
4044 return;
4045
4046 double rhok = delta / theta, sigma = theta / delta;
4047 for (unsigned int k = 0; k < data.degree - 1; ++k)
4048 {
4049 double factor1 = 0.0;
4050 double factor2 = 0.0;
4051
4052 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4053 {
4054 factor1 = (2 * k + 1.) / (2 * k + 5.);
4055 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4056 }
4057 else
4058 {
4059 const double rhokp = 1. / (2. * sigma - rhok);
4060 factor1 = rhokp * rhok;
4061 factor2 = 2. * rhokp / delta;
4062 rhok = rhokp;
4063 }
4064
4065 internal::PreconditionChebyshevImplementation::vmult_and_update(
4066 *matrix_ptr,
4067 *data.preconditioner,
4068 rhs,
4069 k + 2,
4070 factor1,
4071 factor2,
4072 solution,
4073 solution_old,
4074 temp_vector1,
4075 temp_vector2);
4076 }
4077}
4078
4079
4080
4081template <typename MatrixType, typename VectorType, typename PreconditionerType>
4082inline void
4083PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tstep(
4084 VectorType &solution,
4085 const VectorType &rhs) const
4086{
4087 std::lock_guard<std::mutex> lock(mutex);
4088 if (eigenvalues_are_initialized == false)
4089 estimate_eigenvalues(rhs);
4090
4091 matrix_ptr->Tvmult(temp_vector1, solution);
4092 internal::PreconditionChebyshevImplementation::vector_updates(
4093 rhs,
4094 *data.preconditioner,
4095 1,
4096 0.,
4097 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4098 (4. / (3. * delta)) :
4099 (1. / theta),
4100 solution_old,
4101 temp_vector1,
4102 temp_vector2,
4103 solution);
4104
4105 if (data.degree < 2 || std::abs(delta) < 1e-40)
4106 return;
4107
4108 double rhok = delta / theta, sigma = theta / delta;
4109 for (unsigned int k = 0; k < data.degree - 1; ++k)
4110 {
4111 double factor1 = 0.0;
4112 double factor2 = 0.0;
4113
4114 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4115 {
4116 factor1 = (2 * k + 1.) / (2 * k + 5.);
4117 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4118 }
4119 else
4120 {
4121 const double rhokp = 1. / (2. * sigma - rhok);
4122 factor1 = rhokp * rhok;
4123 factor2 = 2. * rhokp / delta;
4124 rhok = rhokp;
4125 }
4126
4127 matrix_ptr->Tvmult(temp_vector1, solution);
4128 internal::PreconditionChebyshevImplementation::vector_updates(
4129 rhs,
4130 *data.preconditioner,
4131 k + 2,
4132 factor1,
4133 factor2,
4134 solution_old,
4135 temp_vector1,
4136 temp_vector2,
4137 solution);
4138 }
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>::m() const
4148{
4149 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4150 return matrix_ptr->m();
4151}
4152
4153
4154
4155template <typename MatrixType, typename VectorType, typename PreconditionerType>
4156inline typename PreconditionChebyshev<MatrixType,
4157 VectorType,
4158 PreconditionerType>::size_type
4159PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::n() const
4160{
4161 Assert(matrix_ptr != nullptr, ExcNotInitialized());
4162 return matrix_ptr->n();
4163}
4164
4165#endif // DOXYGEN
4166
4167DEAL_II_NAMESPACE_CLOSE
4168
4169#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:2260
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:1664
PreconditionJacobi::PreconditionerType
internal::PreconditionRelaxation::PreconditionJacobiImpl< MatrixType > PreconditionerType
Definition precondition.h:1657
PreconditionJacobi::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionPSOR::AdditionalData::AdditionalData
AdditionalData(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
PreconditionPSOR::AdditionalData::parameters
BaseClass::AdditionalData parameters
Definition precondition.h:1876
PreconditionPSOR::AdditionalData::inverse_permutation
const std::vector< size_type > & inverse_permutation
Definition precondition.h:1872
PreconditionPSOR::AdditionalData::permutation
const std::vector< size_type > & permutation
Definition precondition.h:1868
PreconditionPSOR::size_type
typename BaseClass::size_type size_type
Definition precondition.h:1842
PreconditionPSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &additional_data)
PreconditionPSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionPSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1835
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:1727
PreconditionSOR::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1734
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:1786
PreconditionSSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionSSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1779
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:1615
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:1736
identity_matrix.h
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exceptions.h:1811
vector_memory.h
memory_space.h
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