16 #ifndef dealii_precondition_h
17 #define dealii_precondition_h
41 template <
typename number>
43 template <
typename number>
49 template <
typename,
typename>
112 template <
typename MatrixType>
120 template <
typename VectorType>
122 vmult(VectorType &,
const VectorType &)
const;
128 template <
typename VectorType>
130 Tvmult(VectorType &,
const VectorType &)
const;
135 template <
typename VectorType>
143 template <
typename VectorType>
241 template <
typename MatrixType>
248 template <
typename VectorType>
250 vmult(VectorType &,
const VectorType &)
const;
256 template <
typename VectorType>
258 Tvmult(VectorType &,
const VectorType &)
const;
262 template <
typename VectorType>
270 template <
typename VectorType>
360 template <
typename MatrixType = SparseMatrix<
double>,
361 typename VectorType = Vector<
double>>
369 const VectorType &)
const;
383 vmult(VectorType &dst,
const VectorType &src)
const;
404 template <
typename MatrixType = SparseMatrix<
double>,
405 typename PreconditionerType = IdentityMatrix>
475 template <
typename VectorType>
477 vmult(VectorType &,
const VectorType &)
const;
483 template <
typename VectorType>
485 Tvmult(VectorType &,
const VectorType &)
const;
490 template <
typename VectorType>
492 step(VectorType &x,
const VectorType &rhs)
const;
497 template <
typename VectorType>
499 Tstep(VectorType &x,
const VectorType &rhs)
const;
532 template <
typename MatrixType,
typename VectorType>
533 using vmult_functions_t = decltype(std::declval<const MatrixType>().vmult(
534 std::declval<VectorType &>(),
535 std::declval<const VectorType &>(),
537 const std::function<
void(
const unsigned int,
const unsigned int)> &>(),
539 const std::function<
void(
const unsigned int,
const unsigned int)> &>()));
541 template <
typename MatrixType,
543 typename PreconditionerType>
544 constexpr
bool has_vmult_with_std_functions =
545 is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
546 std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> &&
547 (std::is_same_v<VectorType,
555 template <
typename MatrixType,
typename VectorType>
556 constexpr
bool has_vmult_with_std_functions_for_precondition =
557 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
561 template <
typename T,
typename VectorType>
562 using Tvmult_t = decltype(std::declval<const T>().Tvmult(
563 std::declval<VectorType &>(),
564 std::declval<const VectorType &>()));
566 template <
typename T,
typename VectorType>
567 constexpr
bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
569 template <
typename T,
typename VectorType>
570 using step_t = decltype(std::declval<const T>().step(
571 std::declval<VectorType &>(),
572 std::declval<const VectorType &>()));
574 template <
typename T,
typename VectorType>
575 constexpr
bool has_step = is_supported_operation<step_t, T, VectorType>;
577 template <
typename T,
typename VectorType>
579 decltype(std::declval<const T>().step(std::declval<VectorType &>(),
580 std::declval<const VectorType &>(),
581 std::declval<const double>()));
583 template <
typename T,
typename VectorType>
584 constexpr
bool has_step_omega =
585 is_supported_operation<step_omega_t, T, VectorType>;
587 template <
typename T,
typename VectorType>
588 using Tstep_t = decltype(std::declval<const T>().Tstep(
589 std::declval<VectorType &>(),
590 std::declval<const VectorType &>()));
592 template <
typename T,
typename VectorType>
593 constexpr
bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
595 template <
typename T,
typename VectorType>
596 using Tstep_omega_t =
597 decltype(std::declval<const T>().Tstep(std::declval<VectorType &>(),
598 std::declval<const VectorType &>(),
599 std::declval<const double>()));
601 template <
typename T,
typename VectorType>
602 constexpr
bool has_Tstep_omega =
603 is_supported_operation<Tstep_omega_t, T, VectorType>;
605 template <
typename T,
typename VectorType>
606 using jacobi_step_t = decltype(std::declval<const T>().Jacobi_step(
607 std::declval<VectorType &>(),
608 std::declval<const VectorType &>(),
609 std::declval<const double>()));
611 template <
typename T,
typename VectorType>
612 constexpr
bool has_jacobi_step =
613 is_supported_operation<jacobi_step_t, T, VectorType>;
615 template <
typename T,
typename VectorType>
616 using SOR_step_t = decltype(std::declval<const T>().SOR_step(
617 std::declval<VectorType &>(),
618 std::declval<const VectorType &>(),
619 std::declval<const double>()));
621 template <
typename T,
typename VectorType>
622 constexpr
bool has_SOR_step =
623 is_supported_operation<SOR_step_t, T, VectorType>;
625 template <
typename T,
typename VectorType>
626 using SSOR_step_t = decltype(std::declval<const T>().SSOR_step(
627 std::declval<VectorType &>(),
628 std::declval<const VectorType &>(),
629 std::declval<const double>()));
631 template <
typename T,
typename VectorType>
632 constexpr
bool has_SSOR_step =
633 is_supported_operation<SSOR_step_t, T, VectorType>;
635 template <
typename MatrixType>
636 class PreconditionJacobiImpl
639 PreconditionJacobiImpl(
const MatrixType &
A,
const double relaxation)
641 , relaxation(relaxation)
644 template <
typename VectorType>
646 vmult(VectorType &dst,
const VectorType &src)
const
648 this->A->precondition_Jacobi(dst, src, this->relaxation);
651 template <
typename VectorType>
653 Tvmult(VectorType &dst,
const VectorType &src)
const
656 this->vmult(dst, src);
659 template <
typename VectorType,
660 std::enable_if_t<has_jacobi_step<MatrixType, VectorType>,
661 MatrixType> * =
nullptr>
663 step(VectorType &dst,
const VectorType &src)
const
665 this->A->Jacobi_step(dst, src, this->relaxation);
668 template <
typename VectorType,
669 std::enable_if_t<!has_jacobi_step<MatrixType, VectorType>,
670 MatrixType> * =
nullptr>
672 step(VectorType &,
const VectorType &)
const
676 "Matrix A does not provide a Jacobi_step() function!"));
679 template <
typename VectorType>
681 Tstep(VectorType &dst,
const VectorType &src)
const
684 this->step(dst, src);
689 const double relaxation;
692 template <
typename MatrixType>
693 class PreconditionSORImpl
696 PreconditionSORImpl(
const MatrixType &
A,
const double relaxation)
698 , relaxation(relaxation)
701 template <
typename VectorType>
703 vmult(VectorType &dst,
const VectorType &src)
const
705 this->A->precondition_SOR(dst, src, this->relaxation);
708 template <
typename VectorType>
710 Tvmult(VectorType &dst,
const VectorType &src)
const
712 this->A->precondition_TSOR(dst, src, this->relaxation);
715 template <
typename VectorType,
716 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
717 MatrixType> * =
nullptr>
719 step(VectorType &dst,
const VectorType &src)
const
721 this->A->SOR_step(dst, src, this->relaxation);
724 template <
typename VectorType,
725 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
726 MatrixType> * =
nullptr>
728 step(VectorType &,
const VectorType &)
const
732 "Matrix A does not provide a SOR_step() function!"));
735 template <
typename VectorType,
736 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
737 MatrixType> * =
nullptr>
739 Tstep(VectorType &dst,
const VectorType &src)
const
741 this->A->TSOR_step(dst, src, this->relaxation);
744 template <
typename VectorType,
745 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
746 MatrixType> * =
nullptr>
748 Tstep(VectorType &,
const VectorType &)
const
752 "Matrix A does not provide a TSOR_step() function!"));
757 const double relaxation;
760 template <
typename MatrixType>
761 class PreconditionSSORImpl
766 PreconditionSSORImpl(
const MatrixType &
A,
const double relaxation)
768 , relaxation(relaxation)
780 pos_right_of_diagonal.resize(n,
static_cast<std::size_t
>(-1));
788 typename MatrixType::value_type>::const_iterator it =
790 for (; it < mat->
end(row); ++it)
791 if (it->column() > row)
793 pos_right_of_diagonal[row] = it - mat->
begin();
798 template <
typename VectorType>
800 vmult(VectorType &dst,
const VectorType &src)
const
802 this->A->precondition_SSOR(dst,
805 pos_right_of_diagonal);
808 template <
typename VectorType>
810 Tvmult(VectorType &dst,
const VectorType &src)
const
812 this->A->precondition_SSOR(dst,
815 pos_right_of_diagonal);
818 template <
typename VectorType,
819 std::enable_if_t<has_SSOR_step<MatrixType, VectorType>,
820 MatrixType> * =
nullptr>
822 step(VectorType &dst,
const VectorType &src)
const
824 this->A->SSOR_step(dst, src, this->relaxation);
827 template <
typename VectorType,
828 std::enable_if_t<!has_SSOR_step<MatrixType, VectorType>,
829 MatrixType> * =
nullptr>
831 step(VectorType &,
const VectorType &)
const
835 "Matrix A does not provide a SSOR_step() function!"));
838 template <
typename VectorType>
840 Tstep(VectorType &dst,
const VectorType &src)
const
843 this->step(dst, src);
848 const double relaxation;
854 std::vector<std::size_t> pos_right_of_diagonal;
857 template <
typename MatrixType>
858 class PreconditionPSORImpl
863 PreconditionPSORImpl(
const MatrixType &
A,
864 const double relaxation,
865 const std::vector<size_type> &permutation,
866 const std::vector<size_type> &inverse_permutation)
868 , relaxation(relaxation)
869 , permutation(permutation)
870 , inverse_permutation(inverse_permutation)
873 template <
typename VectorType>
875 vmult(VectorType &dst,
const VectorType &src)
const
878 this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
881 template <
typename VectorType>
883 Tvmult(VectorType &dst,
const VectorType &src)
const
886 this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
891 const double relaxation;
893 const std::vector<size_type> &permutation;
894 const std::vector<size_type> &inverse_permutation;
897 template <
typename MatrixType,
898 typename PreconditionerType,
900 std::enable_if_t<has_step_omega<PreconditionerType, VectorType>,
901 PreconditionerType> * =
nullptr>
903 step(
const MatrixType &,
904 const PreconditionerType &preconditioner,
906 const VectorType &src,
907 const double relaxation,
911 preconditioner.step(dst, src, relaxation);
916 typename PreconditionerType,
918 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
919 has_step<PreconditionerType, VectorType>,
920 PreconditionerType> * =
nullptr>
922 step(
const MatrixType &,
923 const PreconditionerType &preconditioner,
925 const VectorType &src,
926 const double relaxation,
934 preconditioner.step(dst, src);
939 typename PreconditionerType,
941 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
942 !has_step<PreconditionerType, VectorType>,
943 PreconditionerType> * =
nullptr>
945 step(
const MatrixType &
A,
946 const PreconditionerType &preconditioner,
948 const VectorType &src,
949 const double relaxation,
950 VectorType &residual,
953 residual.reinit(dst,
true);
954 tmp.reinit(dst,
true);
956 A.vmult(residual, dst);
957 residual.sadd(-1.0, 1.0, src);
959 preconditioner.vmult(tmp, residual);
960 dst.add(relaxation, tmp);
963 template <
typename MatrixType,
964 typename PreconditionerType,
966 std::enable_if_t<has_Tstep_omega<PreconditionerType, VectorType>,
967 PreconditionerType> * =
nullptr>
969 Tstep(
const MatrixType &,
970 const PreconditionerType &preconditioner,
972 const VectorType &src,
973 const double relaxation,
977 preconditioner.Tstep(dst, src, relaxation);
982 typename PreconditionerType,
984 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
985 has_Tstep<PreconditionerType, VectorType>,
986 PreconditionerType> * =
nullptr>
988 Tstep(
const MatrixType &,
989 const PreconditionerType &preconditioner,
991 const VectorType &src,
992 const double relaxation,
1000 preconditioner.Tstep(dst, src);
1003 template <
typename MatrixType,
1004 typename VectorType,
1005 std::enable_if_t<has_Tvmult<MatrixType, VectorType>, MatrixType>
1008 Tvmult(
const MatrixType &
A, VectorType &dst,
const VectorType &src)
1013 template <
typename MatrixType,
1014 typename VectorType,
1015 std::enable_if_t<!has_Tvmult<MatrixType, VectorType>, MatrixType>
1018 Tvmult(
const MatrixType &, VectorType &,
const VectorType &)
1021 ExcMessage(
"Matrix A does not provide a Tvmult() function!"));
1025 typename MatrixType,
1026 typename PreconditionerType,
1027 typename VectorType,
1028 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1029 !has_Tstep<PreconditionerType, VectorType>,
1030 PreconditionerType> * =
nullptr>
1032 Tstep(
const MatrixType &
A,
1033 const PreconditionerType &preconditioner,
1035 const VectorType &src,
1036 const double relaxation,
1037 VectorType &residual,
1040 residual.reinit(dst,
true);
1041 tmp.reinit(dst,
true);
1043 Tvmult(
A, residual, dst);
1044 residual.sadd(-1.0, 1.0, src);
1046 Tvmult(preconditioner, tmp, residual);
1047 dst.add(relaxation, tmp);
1051 template <
typename MatrixType,
1052 typename PreconditionerType,
1053 typename VectorType,
1054 std::enable_if_t<!has_vmult_with_std_functions_for_precondition<
1059 step_operations(
const MatrixType &
A,
1060 const PreconditionerType &preconditioner,
1062 const VectorType &src,
1063 const double relaxation,
1066 const unsigned int i,
1067 const bool transposed)
1072 Tvmult(preconditioner, dst, src);
1074 preconditioner.vmult(dst, src);
1076 if (relaxation != 1.0)
1082 Tstep(
A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1084 step(
A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1091 typename MatrixType,
1092 typename PreconditionerType,
1093 typename VectorType,
1095 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1097 !has_vmult_with_std_functions_for_precondition<MatrixType,
1101 step_operations(
const MatrixType &
A,
1102 const PreconditionerType &preconditioner,
1104 const VectorType &src,
1105 const double relaxation,
1108 const unsigned int i,
1109 const bool transposed)
1112 using Number =
typename VectorType::value_type;
1116 Number *dst_ptr = dst.begin();
1117 const Number *src_ptr = src.begin();
1119 preconditioner.vmult(
1122 [&](
const unsigned int start_range,
const unsigned int end_range) {
1124 if (end_range > start_range)
1125 std::memset(dst.begin() + start_range,
1127 sizeof(Number) * (end_range - start_range));
1129 [&](
const unsigned int start_range,
const unsigned int end_range) {
1130 if (relaxation == 1.0)
1133 const auto src_ptr = src.begin();
1134 const auto dst_ptr = dst.begin();
1137 for (std::size_t i = start_range; i < end_range; ++i)
1138 dst_ptr[i] *= relaxation;
1143 tmp.reinit(src,
true);
1149 preconditioner.vmult(
1152 [&](
const unsigned int start_range,
const unsigned int end_range) {
1153 const auto src_ptr = src.begin();
1154 const auto tmp_ptr = tmp.begin();
1156 if (relaxation == 1.0)
1159 for (std::size_t i = start_range; i < end_range; ++i)
1160 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1168 for (std::size_t i = start_range; i < end_range; ++i)
1169 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1172 [&](
const unsigned int,
const unsigned int) {
1182 typename MatrixType,
1183 typename PreconditionerType,
1184 typename VectorType,
1186 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1188 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1191 step_operations(
const MatrixType &
A,
1192 const PreconditionerType &preconditioner,
1194 const VectorType &src,
1195 const double relaxation,
1198 const unsigned int i,
1199 const bool transposed)
1202 using Number =
typename VectorType::value_type;
1206 Number *dst_ptr = dst.begin();
1207 const Number *src_ptr = src.begin();
1209 preconditioner.vmult(
1212 [&](
const unsigned int start_range,
const unsigned int end_range) {
1214 if (end_range > start_range)
1215 std::memset(dst.begin() + start_range,
1217 sizeof(Number) * (end_range - start_range));
1219 [&](
const unsigned int start_range,
const unsigned int end_range) {
1220 if (relaxation == 1.0)
1223 const auto src_ptr = src.begin();
1224 const auto dst_ptr = dst.begin();
1227 for (std::size_t i = start_range; i < end_range; ++i)
1228 dst_ptr[i] *= relaxation;
1233 tmp.reinit(src,
true);
1240 [&](
const unsigned int start_range,
const unsigned int end_range) {
1243 if (end_range > start_range)
1244 std::memset(tmp.begin() + start_range,
1246 sizeof(Number) * (end_range - start_range));
1248 [&](
const unsigned int start_range,
const unsigned int end_range) {
1249 const auto src_ptr = src.begin();
1250 const auto tmp_ptr = tmp.begin();
1252 if (relaxation == 1.0)
1255 for (std::size_t i = start_range; i < end_range; ++i)
1256 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1264 for (std::size_t i = start_range; i < end_range; ++i)
1265 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1269 preconditioner.vmult(dst, tmp, [](
const auto,
const auto) {
1277 template <
typename MatrixType,
1278 typename VectorType,
1279 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1280 !has_vmult_with_std_functions<
1284 VectorType> * =
nullptr>
1286 step_operations(
const MatrixType &
A,
1287 const ::DiagonalMatrix<VectorType> &preconditioner,
1289 const VectorType &src,
1290 const double relaxation,
1293 const unsigned int i,
1294 const bool transposed)
1296 using Number =
typename VectorType::value_type;
1300 Number *dst_ptr = dst.begin();
1301 const Number *src_ptr = src.begin();
1302 const Number *diag_ptr = preconditioner.
get_vector().begin();
1304 if (relaxation == 1.0)
1307 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1308 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1313 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1314 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1319 tmp.reinit(src,
true);
1321 Number *dst_ptr = dst.begin();
1322 const Number *src_ptr = src.begin();
1323 const Number *tmp_ptr = tmp.begin();
1324 const Number *diag_ptr = preconditioner.get_vector().begin();
1327 Tvmult(
A, tmp, dst);
1331 if (relaxation == 1.0)
1334 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1335 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1340 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1342 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1349 template <
typename MatrixType,
1350 typename VectorType,
1351 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1352 has_vmult_with_std_functions<
1356 VectorType> * =
nullptr>
1358 step_operations(
const MatrixType &
A,
1359 const ::DiagonalMatrix<VectorType> &preconditioner,
1361 const VectorType &src,
1362 const double relaxation,
1365 const unsigned int i,
1366 const bool transposed)
1369 using Number =
typename VectorType::value_type;
1373 Number *dst_ptr = dst.begin();
1374 const Number *src_ptr = src.begin();
1375 const Number *diag_ptr = preconditioner.get_vector().begin();
1377 if (relaxation == 1.0)
1380 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1381 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1386 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1387 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1392 tmp.reinit(src,
true);
1399 [&](
const unsigned int start_range,
const unsigned int end_range) {
1401 if (end_range > start_range)
1402 std::memset(tmp.begin() + start_range,
1404 sizeof(Number) * (end_range - start_range));
1406 [&](
const unsigned int begin,
const unsigned int end) {
1407 const Number *dst_ptr = dst.begin();
1408 const Number *src_ptr = src.begin();
1409 Number *tmp_ptr = tmp.begin();
1410 const Number *diag_ptr = preconditioner.get_vector().begin();
1414 if (relaxation == 1.0)
1417 for (std::size_t i =
begin; i <
end; ++i)
1419 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1424 for (std::size_t i =
begin; i <
end; ++i)
1425 tmp_ptr[i] = dst_ptr[i] + relaxation *
1426 (src_ptr[i] - tmp_ptr[i]) *
1468 template <
typename MatrixType = SparseMatrix<
double>>
1472 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1475 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1538 template <
typename MatrixType = SparseMatrix<
double>>
1542 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1545 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1590 template <
typename MatrixType = SparseMatrix<
double>>
1594 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1597 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1646 template <
typename MatrixType = SparseMatrix<
double>>
1650 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1653 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1710 const std::vector<size_type> &permutation,
1711 const std::vector<size_type> &inverse_permutation,
1928 template <
typename MatrixType = SparseMatrix<
double>,
1929 typename VectorType = Vector<
double>,
1930 typename PreconditionerType = DiagonalMatrix<VectorType>>
1988 const unsigned int degree = 1,
2099 vmult(VectorType &dst,
const VectorType &src)
const;
2106 Tvmult(VectorType &dst,
const VectorType &src)
const;
2112 step(VectorType &dst,
const VectorType &src)
const;
2118 Tstep(VectorType &dst,
const VectorType &src)
const;
2186 EigenvalueInformation
2255 template <
typename MatrixType>
2265 template <
typename VectorType>
2274 template <
typename VectorType>
2281 template <
typename VectorType>
2290 template <
typename VectorType>
2322 const double relaxation)
2323 : relaxation(relaxation)
2332 AdditionalData add_data;
2333 relaxation = add_data.relaxation;
2347 template <
typename MatrixType>
2350 const MatrixType &
matrix,
2360 template <
typename VectorType>
2365 std::is_same_v<size_type, typename VectorType::size_type>,
2366 "PreconditionRichardson and VectorType must have the same size_type.");
2373 template <
typename VectorType>
2378 std::is_same_v<size_type, typename VectorType::size_type>,
2379 "PreconditionRichardson and VectorType must have the same size_type.");
2384 template <
typename VectorType>
2389 std::is_same_v<size_type, typename VectorType::size_type>,
2390 "PreconditionRichardson and VectorType must have the same size_type.");
2397 template <
typename VectorType>
2402 std::is_same_v<size_type, typename VectorType::size_type>,
2403 "PreconditionRichardson and VectorType must have the same size_type.");
2424 template <
typename MatrixType,
typename PreconditionerType>
2427 const MatrixType &rA,
2428 const AdditionalData ¶meters)
2431 relaxation = parameters.relaxation;
2435 preconditioner = parameters.preconditioner;
2436 n_iterations = parameters.n_iterations;
2440 template <
typename MatrixType,
typename PreconditionerType>
2445 preconditioner =
nullptr;
2448 template <
typename MatrixType,
typename PreconditionerType>
2457 template <
typename MatrixType,
typename PreconditionerType>
2466 template <
typename MatrixType,
typename PreconditionerType>
2467 template <
typename VectorType>
2471 const VectorType &src)
const
2476 VectorType tmp1, tmp2;
2478 for (
unsigned int i = 0; i < n_iterations; ++i)
2479 internal::PreconditionRelaxation::step_operations(
2480 *
A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i,
false);
2483 template <
typename MatrixType,
typename PreconditionerType>
2484 template <
typename VectorType>
2488 const VectorType &src)
const
2493 VectorType tmp1, tmp2;
2495 for (
unsigned int i = 0; i < n_iterations; ++i)
2496 internal::PreconditionRelaxation::step_operations(
2497 *
A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i,
true);
2500 template <
typename MatrixType,
typename PreconditionerType>
2501 template <
typename VectorType>
2505 const VectorType &src)
const
2510 VectorType tmp1, tmp2;
2512 for (
unsigned int i = 1; i <= n_iterations; ++i)
2513 internal::PreconditionRelaxation::step_operations(
2514 *
A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i,
false);
2517 template <
typename MatrixType,
typename PreconditionerType>
2518 template <
typename VectorType>
2522 const VectorType &src)
const
2527 VectorType tmp1, tmp2;
2529 for (
unsigned int i = 1; i <= n_iterations; ++i)
2530 internal::PreconditionRelaxation::step_operations(
2531 *
A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i,
true);
2536 template <
typename MatrixType>
2539 const AdditionalData ¶meters_in)
2543 AdditionalData parameters;
2544 parameters.relaxation = 1.0;
2545 parameters.n_iterations = parameters_in.n_iterations;
2546 parameters.preconditioner =
2547 std::make_shared<PreconditionerType>(
A, parameters_in.relaxation);
2549 this->BaseClass::initialize(
A, parameters);
2554 template <
typename MatrixType>
2557 const AdditionalData ¶meters_in)
2561 AdditionalData parameters;
2562 parameters.relaxation = 1.0;
2563 parameters.n_iterations = parameters_in.n_iterations;
2564 parameters.preconditioner =
2565 std::make_shared<PreconditionerType>(
A, parameters_in.relaxation);
2567 this->BaseClass::initialize(
A, parameters);
2572 template <
typename MatrixType>
2575 const AdditionalData ¶meters_in)
2579 AdditionalData parameters;
2580 parameters.relaxation = 1.0;
2581 parameters.n_iterations = parameters_in.n_iterations;
2582 parameters.preconditioner =
2583 std::make_shared<PreconditionerType>(
A, parameters_in.relaxation);
2585 this->BaseClass::initialize(
A, parameters);
2592 template <
typename MatrixType>
2595 const MatrixType &
A,
2596 const std::vector<size_type> &p,
2597 const std::vector<size_type> &ip,
2598 const typename BaseClass::AdditionalData ¶meters_in)
2602 typename BaseClass::AdditionalData parameters;
2603 parameters.relaxation = 1.0;
2604 parameters.n_iterations = parameters_in.n_iterations;
2605 parameters.preconditioner =
2606 std::make_shared<PreconditionerType>(
A, parameters_in.relaxation, p, ip);
2608 this->BaseClass::initialize(
A, parameters);
2612 template <
typename MatrixType>
2615 const AdditionalData &additional_data)
2618 additional_data.permutation,
2619 additional_data.inverse_permutation,
2620 additional_data.parameters);
2623 template <
typename MatrixType>
2625 const std::vector<size_type> &permutation,
2626 const std::vector<size_type> &inverse_permutation,
2629 : permutation(permutation)
2630 , inverse_permutation(inverse_permutation)
2631 , parameters(parameters)
2638 template <
typename MatrixType,
typename VectorType>
2640 const MatrixType &M,
2641 const function_ptr method)
2643 , precondition(method)
2648 template <
typename MatrixType,
typename VectorType>
2652 const VectorType &src)
const
2654 (
matrix.*precondition)(dst, src);
2659 template <
typename MatrixType,
typename PreconditionerType>
2661 AdditionalData(
const double relaxation,
const unsigned int n_iterations)
2662 : relaxation(relaxation)
2663 , n_iterations(n_iterations)
2672 namespace PreconditionChebyshevImplementation
2680 template <
typename VectorType,
typename PreconditionerType>
2682 vector_updates(
const VectorType &rhs,
2683 const PreconditionerType &preconditioner,
2684 const unsigned int iteration_index,
2685 const double factor1,
2686 const double factor2,
2687 VectorType &solution_old,
2688 VectorType &temp_vector1,
2689 VectorType &temp_vector2,
2690 VectorType &solution)
2692 if (iteration_index == 0)
2694 solution.equ(factor2, rhs);
2695 preconditioner.vmult(solution_old, solution);
2697 else if (iteration_index == 1)
2700 temp_vector1.sadd(-1.0, 1.0, rhs);
2701 preconditioner.vmult(solution_old, temp_vector1);
2704 solution_old.sadd(factor2, 1 + factor1, solution);
2709 temp_vector1.sadd(-1.0, 1.0, rhs);
2710 preconditioner.vmult(temp_vector2, temp_vector1);
2713 solution_old.sadd(-factor1, factor2, temp_vector2);
2714 solution_old.add(1 + factor1, solution);
2717 solution.swap(solution_old);
2723 typename PreconditionerType,
2725 !has_vmult_with_std_functions_for_precondition<
2732 const PreconditionerType &preconditioner,
2733 const unsigned int iteration_index,
2734 const double factor1_,
2735 const double factor2_,
2744 const Number factor1 = factor1_;
2745 const Number factor1_plus_1 = 1. + factor1_;
2746 const Number factor2 = factor2_;
2748 if (iteration_index == 0)
2750 const auto solution_old_ptr = solution_old.
begin();
2753 preconditioner.vmult(solution_old, rhs);
2758 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
2760 else if (iteration_index == 1)
2762 const auto solution_ptr = solution.
begin();
2763 const auto solution_old_ptr = solution_old.
begin();
2766 temp_vector1.
sadd(-1.0, 1.0, rhs);
2768 preconditioner.vmult(solution_old, temp_vector1);
2773 solution_old_ptr[i] =
2774 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
2778 const auto solution_ptr = solution.
begin();
2779 const auto solution_old_ptr = solution_old.
begin();
2780 const auto temp_vector2_ptr = temp_vector2.
begin();
2783 temp_vector1.
sadd(-1.0, 1.0, rhs);
2785 preconditioner.vmult(temp_vector2, temp_vector1);
2790 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
2791 factor1 * solution_old_ptr[i] +
2792 temp_vector2_ptr[i] * factor2;
2795 solution.
swap(solution_old);
2800 typename PreconditionerType,
2802 has_vmult_with_std_functions_for_precondition<
2809 const PreconditionerType &preconditioner,
2810 const unsigned int iteration_index,
2811 const double factor1_,
2812 const double factor2_,
2821 const Number factor1 = factor1_;
2822 const Number factor1_plus_1 = 1. + factor1_;
2823 const Number factor2 = factor2_;
2825 const auto rhs_ptr = rhs.
begin();
2826 const auto temp_vector1_ptr = temp_vector1.
begin();
2827 const auto temp_vector2_ptr = temp_vector2.
begin();
2828 const auto solution_ptr = solution.
begin();
2829 const auto solution_old_ptr = solution_old.
begin();
2831 if (iteration_index == 0)
2833 preconditioner.vmult(
2836 [&](
const auto start_range,
const auto end_range) {
2837 if (end_range > start_range)
2838 std::memset(solution.
begin() + start_range,
2840 sizeof(Number) * (end_range - start_range));
2842 [&](
const auto begin,
const auto end) {
2844 for (std::size_t i =
begin; i <
end; ++i)
2845 solution_ptr[i] *= factor2;
2850 preconditioner.vmult(
2853 [&](
const auto begin,
const auto end) {
2860 for (std::size_t i =
begin; i <
end; ++i)
2861 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
2863 [&](
const auto begin,
const auto end) {
2864 if (iteration_index == 1)
2867 for (std::size_t i =
begin; i <
end; ++i)
2868 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
2869 factor2 * temp_vector2_ptr[i];
2874 for (std::size_t i =
begin; i <
end; ++i)
2875 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
2876 factor1 * solution_old_ptr[i] +
2877 factor2 * temp_vector2_ptr[i];
2882 if (iteration_index > 0)
2884 solution_old.
swap(temp_vector2);
2885 solution_old.
swap(solution);
2892 template <
typename Number>
2893 struct VectorUpdater
2895 VectorUpdater(
const Number *rhs,
2896 const Number *matrix_diagonal_inverse,
2897 const unsigned int iteration_index,
2898 const Number factor1,
2899 const Number factor2,
2900 Number *solution_old,
2904 , matrix_diagonal_inverse(matrix_diagonal_inverse)
2905 , iteration_index(iteration_index)
2908 , solution_old(solution_old)
2909 , tmp_vector(tmp_vector)
2910 , solution(solution)
2914 apply_to_subrange(
const std::size_t
begin,
const std::size_t
end)
const
2920 const Number factor1 = this->factor1;
2921 const Number factor1_plus_1 = 1. + this->factor1;
2922 const Number factor2 = this->factor2;
2923 if (iteration_index == 0)
2926 for (std::size_t i =
begin; i <
end; ++i)
2927 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
2929 else if (iteration_index == 1)
2933 for (std::size_t i =
begin; i <
end; ++i)
2937 factor1_plus_1 * solution[i] +
2938 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2945 for (std::size_t i =
begin; i <
end; ++i)
2952 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
2953 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2958 const Number *matrix_diagonal_inverse;
2959 const unsigned int iteration_index;
2960 const Number factor1;
2961 const Number factor2;
2962 mutable Number *solution_old;
2963 mutable Number *tmp_vector;
2964 mutable Number *solution;
2967 template <
typename Number>
2970 VectorUpdatesRange(
const VectorUpdater<Number> &updater,
2971 const std::size_t size)
2975 VectorUpdatesRange::apply_to_subrange(0, size);
2983 ~VectorUpdatesRange()
override =
default;
2986 apply_to_subrange(
const std::size_t
begin,
2987 const std::size_t
end)
const override
2989 updater.apply_to_subrange(
begin,
end);
2992 const VectorUpdater<Number> &updater;
2996 template <
typename Number>
2999 const ::Vector<Number> &rhs,
3001 const unsigned int iteration_index,
3002 const double factor1,
3003 const double factor2,
3009 VectorUpdater<Number> upd(rhs.begin(),
3010 jacobi.get_vector().begin(),
3014 solution_old.
begin(),
3015 temp_vector1.
begin(),
3017 VectorUpdatesRange<Number>(upd, rhs.size());
3020 if (iteration_index == 0)
3027 solution.
swap(temp_vector1);
3028 solution_old.
swap(temp_vector1);
3033 template <
typename Number>
3037 const ::DiagonalMatrix<
3039 const unsigned int iteration_index,
3040 const double factor1,
3041 const double factor2,
3049 VectorUpdater<Number> upd(rhs.
begin(),
3050 jacobi.get_vector().begin(),
3054 solution_old.
begin(),
3055 temp_vector1.
begin(),
3060 if (iteration_index == 0)
3067 solution.
swap(temp_vector1);
3068 solution_old.
swap(temp_vector1);
3077 typename MatrixType,
3078 typename VectorType,
3079 typename PreconditionerType,
3081 !has_vmult_with_std_functions<MatrixType,
3083 PreconditionerType> &&
3084 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3086 has_vmult_with_std_functions_for_precondition<MatrixType,
3090 vmult_and_update(
const MatrixType &
matrix,
3091 const PreconditionerType &preconditioner,
3092 const VectorType &rhs,
3093 const unsigned int iteration_index,
3094 const double factor1,
3095 const double factor2,
3096 VectorType &solution,
3097 VectorType &solution_old,
3098 VectorType &temp_vector1,
3099 VectorType &temp_vector2)
3101 if (iteration_index > 0)
3102 matrix.vmult(temp_vector1, solution);
3117 typename MatrixType,
3118 typename VectorType,
3119 typename PreconditionerType,
3121 !has_vmult_with_std_functions<MatrixType,
3123 PreconditionerType> &&
3124 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3126 has_vmult_with_std_functions_for_precondition<MatrixType,
3130 vmult_and_update(
const MatrixType &
matrix,
3131 const PreconditionerType &preconditioner,
3132 const VectorType &rhs,
3133 const unsigned int iteration_index,
3134 const double factor1_,
3135 const double factor2_,
3136 VectorType &solution,
3137 VectorType &solution_old,
3138 VectorType &temp_vector1,
3139 VectorType &temp_vector2)
3141 using Number =
typename VectorType::value_type;
3143 const Number factor1 = factor1_;
3144 const Number factor1_plus_1 = 1. + factor1_;
3145 const Number factor2 = factor2_;
3147 if (iteration_index == 0)
3149 preconditioner.vmult(
3152 [&](
const unsigned int start_range,
const unsigned int end_range) {
3154 if (end_range > start_range)
3155 std::memset(solution.begin() + start_range,
3157 sizeof(Number) * (end_range - start_range));
3159 [&](
const unsigned int start_range,
const unsigned int end_range) {
3160 const auto solution_ptr = solution.begin();
3163 for (std::size_t i = start_range; i < end_range; ++i)
3164 solution_ptr[i] *= factor2;
3169 temp_vector1.reinit(rhs,
true);
3170 temp_vector2.reinit(rhs,
true);
3176 [&](
const unsigned int start_range,
const unsigned int end_range) {
3179 if (end_range > start_range)
3180 std::memset(temp_vector1.begin() + start_range,
3182 sizeof(Number) * (end_range - start_range));
3184 [&](
const unsigned int start_range,
const unsigned int end_range) {
3185 const auto rhs_ptr = rhs.begin();
3186 const auto tmp_ptr = temp_vector1.begin();
3189 for (std::size_t i = start_range; i < end_range; ++i)
3190 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3194 preconditioner.vmult(
3197 [&](
const unsigned int start_range,
const unsigned int end_range) {
3200 if (end_range > start_range)
3201 std::memset(temp_vector2.begin() + start_range,
3203 sizeof(Number) * (end_range - start_range));
3205 [&](
const unsigned int start_range,
const unsigned int end_range) {
3206 const auto solution_ptr = solution.begin();
3207 const auto solution_old_ptr = solution_old.begin();
3208 const auto tmp_ptr = temp_vector2.begin();
3210 if (iteration_index == 1)
3213 for (std::size_t i = start_range; i < end_range; ++i)
3215 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3220 for (std::size_t i = start_range; i < end_range; ++i)
3221 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3222 factor1 * solution_old_ptr[i] +
3223 factor2 * tmp_ptr[i];
3227 solution.swap(temp_vector2);
3228 solution_old.swap(temp_vector2);
3234 template <
typename MatrixType,
3235 typename VectorType,
3236 typename PreconditionerType,
3237 std::enable_if_t<has_vmult_with_std_functions<MatrixType,
3239 PreconditionerType>,
3242 vmult_and_update(
const MatrixType &
matrix,
3243 const PreconditionerType &preconditioner,
3244 const VectorType &rhs,
3245 const unsigned int iteration_index,
3246 const double factor1,
3247 const double factor2,
3248 VectorType &solution,
3249 VectorType &solution_old,
3250 VectorType &temp_vector1,
3253 using Number =
typename VectorType::value_type;
3254 VectorUpdater<Number> updater(rhs.begin(),
3255 preconditioner.get_vector().begin(),
3259 solution_old.begin(),
3260 temp_vector1.begin(),
3262 if (iteration_index > 0)
3266 [&](
const unsigned int start_range,
const unsigned int end_range) {
3269 if (end_range > start_range)
3270 std::memset(temp_vector1.begin() + start_range,
3272 sizeof(Number) * (end_range - start_range));
3274 [&](
const unsigned int start_range,
const unsigned int end_range) {
3275 if (end_range > start_range)
3276 updater.apply_to_subrange(start_range, end_range);
3279 updater.apply_to_subrange(0
U, solution.locally_owned_size());
3282 if (iteration_index == 0)
3289 solution.swap(temp_vector1);
3290 solution_old.swap(temp_vector1);
3294 template <
typename MatrixType,
typename PreconditionerType>
3296 initialize_preconditioner(
3297 const MatrixType &
matrix,
3298 std::shared_ptr<PreconditionerType> &preconditioner)
3301 (void)preconditioner;
3305 template <
typename MatrixType,
typename VectorType>
3307 initialize_preconditioner(
3308 const MatrixType &
matrix,
3311 if (preconditioner.get() ==
nullptr || preconditioner->m() !=
matrix.m())
3313 if (preconditioner.get() ==
nullptr)
3315 std::make_shared<::DiagonalMatrix<VectorType>>();
3318 preconditioner->m() == 0,
3320 "Preconditioner appears to be initialized but not sized correctly"));
3323 if (preconditioner->m() !=
matrix.m())
3325 preconditioner->get_vector().reinit(
matrix.m());
3327 preconditioner->get_vector()(i) = 1. /
matrix.el(i, i);
3332 template <
typename VectorType>
3334 set_initial_guess(VectorType &vector)
3336 vector = 1. /
std::sqrt(
static_cast<double>(vector.size()));
3337 if (vector.locally_owned_elements().is_element(0))
3341 template <
typename Number>
3349 for (
unsigned int i = 0; i < vector.
size(); ++i)
3352 const Number mean_value = vector.
mean_value();
3353 vector.
add(-mean_value);
3356 template <
typename Number>
3361 for (
unsigned int block = 0; block < vector.
n_blocks(); ++block)
3362 set_initial_guess(vector.
block(block));
3365 template <
typename Number,
typename MemorySpace>
3382 Kokkos::RangePolicy<
typename MemorySpace::kokkos_space::execution_space,
3383 Kokkos::IndexType<types::global_dof_index>>
3384 policy(0, n_local_elements);
3386 "::PreconditionChebyshev::set_initial_guess",
3389 values_ptr[i] = (i + first_local_range) % 11;
3391 const Number mean_value = vector.
mean_value();
3392 vector.
add(-mean_value);
3395 struct EigenvalueTracker
3404 std::vector<double>
values;
3409 template <
typename MatrixType,
3410 typename VectorType,
3411 typename PreconditionerType>
3413 power_iteration(
const MatrixType &
matrix,
3414 VectorType &eigenvector,
3415 const PreconditionerType &preconditioner,
3416 const unsigned int n_iterations)
3418 double eigenvalue_estimate = 0.;
3419 eigenvector /= eigenvector.l2_norm();
3420 VectorType vector1, vector2;
3421 vector1.reinit(eigenvector,
true);
3422 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
3423 vector2.reinit(eigenvector,
true);
3424 for (
unsigned int i = 0; i < n_iterations; ++i)
3426 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
3428 matrix.vmult(vector2, eigenvector);
3429 preconditioner.vmult(vector1, vector2);
3432 matrix.vmult(vector1, eigenvector);
3434 eigenvalue_estimate = eigenvector * vector1;
3436 vector1 /= vector1.l2_norm();
3437 eigenvector.swap(vector1);
3439 return eigenvalue_estimate;
3446 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3449 const double smoothing_range,
3450 const unsigned int eig_cg_n_iterations,
3451 const double eig_cg_residual,
3452 const double max_eigenvalue,
3453 const EigenvalueAlgorithm eigenvalue_algorithm,
3454 const PolynomialType polynomial_type)
3456 , smoothing_range(smoothing_range)
3457 , eig_cg_n_iterations(eig_cg_n_iterations)
3458 , eig_cg_residual(eig_cg_residual)
3459 , max_eigenvalue(max_eigenvalue)
3460 , eigenvalue_algorithm(eigenvalue_algorithm)
3461 , polynomial_type(polynomial_type)
3466 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3469 PreconditionerType>::AdditionalData &
3473 degree = other_data.
degree;
3474 smoothing_range = other_data.smoothing_range;
3475 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3476 eig_cg_residual = other_data.eig_cg_residual;
3477 max_eigenvalue = other_data.max_eigenvalue;
3478 preconditioner = other_data.preconditioner;
3479 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3480 polynomial_type = other_data.polynomial_type;
3481 constraints.copy_from(other_data.constraints);
3488 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3493 , eigenvalues_are_initialized(false)
3496 std::is_same_v<size_type, typename VectorType::size_type>,
3497 "PreconditionChebyshev and VectorType must have the same size_type.");
3502 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3505 const MatrixType &
matrix,
3506 const AdditionalData &additional_data)
3509 data = additional_data;
3511 ExcMessage(
"The degree of the Chebyshev method must be positive."));
3512 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3513 matrix, data.preconditioner);
3514 eigenvalues_are_initialized =
false;
3519 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3523 eigenvalues_are_initialized =
false;
3524 theta = delta = 1.0;
3525 matrix_ptr =
nullptr;
3527 VectorType empty_vector;
3528 solution_old.reinit(empty_vector);
3529 temp_vector1.reinit(empty_vector);
3530 temp_vector2.reinit(empty_vector);
3532 data.preconditioner.reset();
3537 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3540 PreconditionerType>::EigenvalueInformation
3548 EigenvalueInformation info{};
3550 solution_old.reinit(src);
3551 temp_vector1.reinit(src,
true);
3553 if (data.eig_cg_n_iterations > 0)
3555 Assert(data.eig_cg_n_iterations > 2,
3557 "Need to set at least two iterations to find eigenvalues."));
3559 internal::PreconditionChebyshevImplementation::EigenvalueTracker
3565 internal::PreconditionChebyshevImplementation::set_initial_guess(
3567 data.constraints.set_zero(temp_vector1);
3569 if (data.eigenvalue_algorithm ==
3570 AdditionalData::EigenvalueAlgorithm::lanczos)
3579 solver.connect_eigenvalues_slot(
3580 [&eigenvalue_tracker](
const std::vector<double> &
eigenvalues) {
3584 solver.solve(*matrix_ptr,
3587 *data.preconditioner);
3589 info.cg_iterations = control.last_step();
3591 else if (data.eigenvalue_algorithm ==
3592 AdditionalData::EigenvalueAlgorithm::power_iteration)
3595 ExcMessage(
"Cannot estimate the minimal eigenvalue with the "
3596 "power iteration"));
3598 eigenvalue_tracker.values.push_back(
3599 internal::PreconditionChebyshevImplementation::power_iteration(
3602 *data.preconditioner,
3603 data.eig_cg_n_iterations));
3609 if (eigenvalue_tracker.values.empty())
3610 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
3613 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
3617 info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
3622 info.max_eigenvalue_estimate = data.max_eigenvalue;
3623 info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
3626 const double alpha = (data.smoothing_range > 1. ?
3627 info.max_eigenvalue_estimate / data.smoothing_range :
3628 std::min(0.9 * info.max_eigenvalue_estimate,
3629 info.min_eigenvalue_estimate));
3638 const double actual_range = info.max_eigenvalue_estimate / alpha;
3639 const double sigma = (1. -
std::sqrt(1. / actual_range)) /
3640 (1. + std::sqrt(1. / actual_range));
3641 const double eps = data.smoothing_range;
3646 1 +
static_cast<unsigned int>(
3647 std::log(1. /
eps + std::sqrt(1. /
eps /
eps - 1.)) /
3648 std::log(1. / sigma));
3651 info.degree = data.degree;
3656 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3657 (info.max_eigenvalue_estimate) :
3658 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3661 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3665 using NumberType =
typename VectorType::value_type;
3673 (std::is_same_v<VectorType,
3674 LinearAlgebra::distributed::
3675 Vector<NumberType, MemorySpace::Default>> ==
false))))
3676 temp_vector2.reinit(src,
true);
3679 VectorType empty_vector;
3680 temp_vector2.reinit(empty_vector);
3685 ->eigenvalues_are_initialized =
true;
3692 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3695 VectorType &solution,
3696 const VectorType &rhs)
const
3698 std::lock_guard<std::mutex> lock(mutex);
3699 if (eigenvalues_are_initialized ==
false)
3700 estimate_eigenvalues(rhs);
3702 internal::PreconditionChebyshevImplementation::vmult_and_update(
3704 *data.preconditioner,
3708 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3709 (4. / (3. * delta)) :
3718 if (data.degree < 2 || std::abs(delta) < 1
e-40)
3721 double rhok = delta / theta, sigma = theta / delta;
3722 for (
unsigned int k = 0; k < data.degree - 1; ++k)
3724 double factor1 = 0.0;
3725 double factor2 = 0.0;
3727 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3729 factor1 = (2 * k + 1.) / (2 * k + 5.);
3730 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3734 const double rhokp = 1. / (2. * sigma - rhok);
3735 factor1 = rhokp * rhok;
3736 factor2 = 2. * rhokp / delta;
3740 internal::PreconditionChebyshevImplementation::vmult_and_update(
3742 *data.preconditioner,
3756 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3759 VectorType &solution,
3760 const VectorType &rhs)
const
3762 std::lock_guard<std::mutex> lock(mutex);
3763 if (eigenvalues_are_initialized ==
false)
3764 estimate_eigenvalues(rhs);
3766 internal::PreconditionChebyshevImplementation::vector_updates(
3768 *data.preconditioner,
3771 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3772 (4. / (3. * delta)) :
3779 if (data.degree < 2 || std::abs(delta) < 1
e-40)
3782 double rhok = delta / theta, sigma = theta / delta;
3783 for (
unsigned int k = 0; k < data.degree - 1; ++k)
3785 double factor1 = 0.0;
3786 double factor2 = 0.0;
3788 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3790 factor1 = (2 * k + 1.) / (2 * k + 5.);
3791 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3795 const double rhokp = 1. / (2. * sigma - rhok);
3796 factor1 = rhokp * rhok;
3797 factor2 = 2. * rhokp / delta;
3801 matrix_ptr->Tvmult(temp_vector1, solution);
3802 internal::PreconditionChebyshevImplementation::vector_updates(
3804 *data.preconditioner,
3817 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3820 VectorType &solution,
3821 const VectorType &rhs)
const
3823 std::lock_guard<std::mutex> lock(mutex);
3824 if (eigenvalues_are_initialized ==
false)
3825 estimate_eigenvalues(rhs);
3827 internal::PreconditionChebyshevImplementation::vmult_and_update(
3829 *data.preconditioner,
3833 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3834 (4. / (3. * delta)) :
3841 if (data.degree < 2 || std::abs(delta) < 1
e-40)
3844 double rhok = delta / theta, sigma = theta / delta;
3845 for (
unsigned int k = 0; k < data.degree - 1; ++k)
3847 double factor1 = 0.0;
3848 double factor2 = 0.0;
3850 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3852 factor1 = (2 * k + 1.) / (2 * k + 5.);
3853 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3857 const double rhokp = 1. / (2. * sigma - rhok);
3858 factor1 = rhokp * rhok;
3859 factor2 = 2. * rhokp / delta;
3863 internal::PreconditionChebyshevImplementation::vmult_and_update(
3865 *data.preconditioner,
3879 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3882 VectorType &solution,
3883 const VectorType &rhs)
const
3885 std::lock_guard<std::mutex> lock(mutex);
3886 if (eigenvalues_are_initialized ==
false)
3887 estimate_eigenvalues(rhs);
3889 matrix_ptr->Tvmult(temp_vector1, solution);
3890 internal::PreconditionChebyshevImplementation::vector_updates(
3892 *data.preconditioner,
3895 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3896 (4. / (3. * delta)) :
3903 if (data.degree < 2 || std::abs(delta) < 1
e-40)
3906 double rhok = delta / theta, sigma = theta / delta;
3907 for (
unsigned int k = 0; k < data.degree - 1; ++k)
3909 double factor1 = 0.0;
3910 double factor2 = 0.0;
3912 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3914 factor1 = (2 * k + 1.) / (2 * k + 5.);
3915 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3919 const double rhokp = 1. / (2. * sigma - rhok);
3920 factor1 = rhokp * rhok;
3921 factor2 = 2. * rhokp / delta;
3925 matrix_ptr->Tvmult(temp_vector1, solution);
3926 internal::PreconditionChebyshevImplementation::vector_updates(
3928 *data.preconditioner,
3941 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3948 return matrix_ptr->m();
3953 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3960 return matrix_ptr->n();
unsigned int n_blocks() const
BlockType & block(const unsigned int i)
VectorType & get_vector()
size_type nth_index_in_set(const size_type local_index) const
void sadd(const Number s, const Number a, const Vector< Number, MemorySpace > &V)
Number mean_value() const
void swap(Vector< Number, MemorySpace > &v)
size_type locally_owned_size() const
::IndexSet locally_owned_elements() const
Number * get_values() const
void Tvmult(VectorType &dst, const VectorType &src) const
void step(VectorType &dst, const VectorType &src) const
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
void Tstep(VectorType &dst, const VectorType &src) const
SmartPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
types::global_dof_index size_type
void vmult(VectorType &dst, const VectorType &src) const
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
bool eigenvalues_are_initialized
void vmult_add(VectorType &, const VectorType &) const
void vmult(VectorType &, const VectorType &) const
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
void Tvmult(VectorType &, const VectorType &) const
types::global_dof_index size_type
void Tvmult_add(VectorType &, const VectorType &) const
typename BaseClass::AdditionalData AdditionalData
internal::PreconditionRelaxation::PreconditionJacobiImpl< MatrixType > PreconditionerType
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
AdditionalData(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData ¶meters=typename BaseClass::AdditionalData())
BaseClass::AdditionalData parameters
const std::vector< size_type > & inverse_permutation
const std::vector< size_type > & permutation
typename BaseClass::size_type size_type
void initialize(const MatrixType &A, const AdditionalData &additional_data)
internal::PreconditionRelaxation::PreconditionPSORImpl< MatrixType > PreconditionerType
void initialize(const MatrixType &A, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData ¶meters=typename BaseClass::AdditionalData())
unsigned int n_iterations
std::shared_ptr< PreconditionerType > preconditioner
AdditionalData(const double relaxation=1., const unsigned int n_iterations=1)
SmartPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
void Tvmult(VectorType &, const VectorType &) const
std::shared_ptr< PreconditionerType > preconditioner
void step(VectorType &x, const VectorType &rhs) const
unsigned int n_iterations
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
void Tstep(VectorType &x, const VectorType &rhs) const
void vmult(VectorType &, const VectorType &) const
types::global_dof_index size_type
AdditionalData(const double relaxation=1.)
types::global_dof_index size_type
void vmult_add(VectorType &, const VectorType &) const
void initialize(const AdditionalData ¶meters)
void initialize(const MatrixType &matrix, const AdditionalData ¶meters)
void vmult(VectorType &, const VectorType &) const
void Tvmult(VectorType &, const VectorType &) const
void Tvmult_add(VectorType &, const VectorType &) const
internal::PreconditionRelaxation::PreconditionSORImpl< MatrixType > PreconditionerType
typename BaseClass::AdditionalData AdditionalData
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
typename BaseClass::AdditionalData AdditionalData
internal::PreconditionRelaxation::PreconditionSSORImpl< MatrixType > PreconditionerType
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
const function_ptr precondition
const MatrixType & matrix
void vmult(VectorType &dst, const VectorType &src) const
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
const_iterator end() const
const_iterator begin() const
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
Number mean_value() const
virtual size_type size() const override
virtual void swap(Vector< Number > &v)
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ matrix
Contents is actually a matrix.
@ eigenvalues
Eigenvalue vector is filled.
types::global_dof_index size_type
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > jacobi(::SymmetricTensor< 2, dim, Number > A)
unsigned int minimum_parallel_grain_size
static const unsigned int invalid_unsigned_int
void parallel_for(Iterator x_begin, Iterator x_end, const Functor &functor, const unsigned int grainsize)
unsigned int global_dof_index
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
PolynomialType polynomial_type
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)
std::shared_ptr< PreconditionerType > preconditioner
AffineConstraints< double > constraints
unsigned int eig_cg_n_iterations
EigenvalueAlgorithm eigenvalue_algorithm
AdditionalData & operator=(const AdditionalData &other_data)