15#ifndef dealii_precondition_h
16#define dealii_precondition_h
40template <
typename number>
42template <
typename number>
48 template <
typename,
typename>
127 template <
typename PreconditionerType>
248 template <
typename MatrixType>
256 template <
typename VectorType>
258 vmult(VectorType &,
const VectorType &)
const;
264 template <
typename VectorType>
266 Tvmult(VectorType &,
const VectorType &)
const;
271 template <
typename VectorType>
279 template <
typename VectorType>
377 template <
typename MatrixType>
384 template <
typename VectorType>
386 vmult(VectorType &,
const VectorType &)
const;
392 template <
typename VectorType>
394 Tvmult(VectorType &,
const VectorType &)
const;
398 template <
typename VectorType>
406 template <
typename VectorType>
496template <
typename MatrixType = SparseMatrix<
double>,
497 typename VectorType = Vector<
double>>
505 const VectorType &) const;
519 vmult(VectorType &dst,
const VectorType &src)
const;
551template <
typename MatrixType = SparseMatrix<
double>,
552 typename PreconditionerType = IdentityMatrix>
580 EigenvalueAlgorithm::lanczos);
626 template <
typename VectorType>
628 vmult(VectorType &,
const VectorType &)
const;
634 template <
typename VectorType>
636 Tvmult(VectorType &,
const VectorType &)
const;
641 template <
typename VectorType>
643 step(VectorType &x,
const VectorType &rhs)
const;
648 template <
typename VectorType>
650 Tstep(VectorType &x,
const VectorType &rhs)
const;
657 template <
typename VectorType>
702 template <
typename MatrixType,
typename VectorType>
703 using vmult_functions_t =
decltype(std::declval<const MatrixType>().vmult(
704 std::declval<VectorType &>(),
705 std::declval<const VectorType &>(),
707 const std::function<
void(
const unsigned int,
const unsigned int)> &>(),
709 const std::function<
void(
const unsigned int,
const unsigned int)> &>()));
711 template <
typename MatrixType,
713 typename PreconditionerType>
714 constexpr bool has_vmult_with_std_functions =
715 is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
716 std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> &&
717 (std::is_same_v<VectorType,
725 template <
typename MatrixType,
typename VectorType>
726 constexpr bool has_vmult_with_std_functions_for_precondition =
727 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
731 template <
typename T,
typename VectorType>
732 using Tvmult_t =
decltype(std::declval<const T>().Tvmult(
733 std::declval<VectorType &>(),
734 std::declval<const VectorType &>()));
736 template <
typename T,
typename VectorType>
737 constexpr bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
739 template <
typename T,
typename VectorType>
740 using step_t =
decltype(std::declval<const T>().step(
741 std::declval<VectorType &>(),
742 std::declval<const VectorType &>()));
744 template <
typename T,
typename VectorType>
745 constexpr bool has_step = is_supported_operation<step_t, T, VectorType>;
747 template <
typename T,
typename VectorType>
749 decltype(std::declval<const T>().step(std::declval<VectorType &>(),
750 std::declval<const VectorType &>(),
751 std::declval<const double>()));
753 template <
typename T,
typename VectorType>
754 constexpr bool has_step_omega =
755 is_supported_operation<step_omega_t, T, VectorType>;
757 template <
typename T,
typename VectorType>
758 using Tstep_t =
decltype(std::declval<const T>().Tstep(
759 std::declval<VectorType &>(),
760 std::declval<const VectorType &>()));
762 template <
typename T,
typename VectorType>
763 constexpr bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
765 template <
typename T,
typename VectorType>
766 using Tstep_omega_t =
767 decltype(std::declval<const T>().Tstep(std::declval<VectorType &>(),
768 std::declval<const VectorType &>(),
769 std::declval<const double>()));
771 template <
typename T,
typename VectorType>
772 constexpr bool has_Tstep_omega =
773 is_supported_operation<Tstep_omega_t, T, VectorType>;
775 template <
typename T,
typename VectorType>
776 using jacobi_step_t =
decltype(std::declval<const T>().Jacobi_step(
777 std::declval<VectorType &>(),
778 std::declval<const VectorType &>(),
779 std::declval<const double>()));
781 template <
typename T,
typename VectorType>
782 constexpr bool has_jacobi_step =
783 is_supported_operation<jacobi_step_t, T, VectorType>;
785 template <
typename T,
typename VectorType>
786 using SOR_step_t =
decltype(std::declval<const T>().SOR_step(
787 std::declval<VectorType &>(),
788 std::declval<const VectorType &>(),
789 std::declval<const double>()));
791 template <
typename T,
typename VectorType>
792 constexpr bool has_SOR_step =
793 is_supported_operation<SOR_step_t, T, VectorType>;
795 template <
typename T,
typename VectorType>
796 using SSOR_step_t =
decltype(std::declval<const T>().SSOR_step(
797 std::declval<VectorType &>(),
798 std::declval<const VectorType &>(),
799 std::declval<const double>()));
801 template <
typename T,
typename VectorType>
802 constexpr bool has_SSOR_step =
803 is_supported_operation<SSOR_step_t, T, VectorType>;
805 template <
typename MatrixType>
806 class PreconditionJacobiImpl
809 PreconditionJacobiImpl(
const MatrixType &A,
const double relaxation)
811 , relaxation(relaxation)
814 template <
typename VectorType>
816 vmult(VectorType &dst,
const VectorType &src)
const
818 this->A->precondition_Jacobi(dst, src, this->relaxation);
821 template <
typename VectorType>
823 Tvmult(VectorType &dst,
const VectorType &src)
const
826 this->vmult(dst, src);
829 template <
typename VectorType,
830 std::enable_if_t<has_jacobi_step<MatrixType, VectorType>,
831 MatrixType> * =
nullptr>
833 step(VectorType &dst,
const VectorType &src)
const
835 this->A->Jacobi_step(dst, src, this->relaxation);
838 template <
typename VectorType,
839 std::enable_if_t<!has_jacobi_step<MatrixType, VectorType>,
840 MatrixType> * =
nullptr>
842 step(VectorType &,
const VectorType &)
const
846 "Matrix A does not provide a Jacobi_step() function!"));
849 template <
typename VectorType>
851 Tstep(VectorType &dst,
const VectorType &src)
const
854 this->step(dst, src);
859 const double relaxation;
862 template <
typename MatrixType>
863 class PreconditionSORImpl
866 PreconditionSORImpl(
const MatrixType &A,
const double relaxation)
868 , relaxation(relaxation)
871 template <
typename VectorType>
873 vmult(VectorType &dst,
const VectorType &src)
const
875 this->A->precondition_SOR(dst, src, this->relaxation);
878 template <
typename VectorType>
880 Tvmult(VectorType &dst,
const VectorType &src)
const
882 this->A->precondition_TSOR(dst, src, this->relaxation);
885 template <
typename VectorType,
886 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
887 MatrixType> * =
nullptr>
889 step(VectorType &dst,
const VectorType &src)
const
891 this->A->SOR_step(dst, src, this->relaxation);
894 template <
typename VectorType,
895 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
896 MatrixType> * =
nullptr>
898 step(VectorType &,
const VectorType &)
const
902 "Matrix A does not provide a SOR_step() function!"));
905 template <
typename VectorType,
906 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
907 MatrixType> * =
nullptr>
909 Tstep(VectorType &dst,
const VectorType &src)
const
911 this->A->TSOR_step(dst, src, this->relaxation);
914 template <
typename VectorType,
915 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
916 MatrixType> * =
nullptr>
918 Tstep(VectorType &,
const VectorType &)
const
922 "Matrix A does not provide a TSOR_step() function!"));
927 const double relaxation;
930 template <
typename MatrixType>
931 class PreconditionSSORImpl
934 using size_type =
typename MatrixType::size_type;
936 PreconditionSSORImpl(
const MatrixType &A,
const double relaxation)
938 , relaxation(relaxation)
949 const size_type n = this->A->n();
950 pos_right_of_diagonal.resize(n,
static_cast<std::size_t
>(-1));
951 for (size_type row = 0; row < n; ++row)
958 typename MatrixType::value_type>::const_iterator it =
960 for (; it < mat->
end(row); ++it)
961 if (it->column() > row)
963 pos_right_of_diagonal[row] = it - mat->
begin();
968 template <
typename VectorType>
970 vmult(VectorType &dst,
const VectorType &src)
const
972 this->A->precondition_SSOR(dst,
975 pos_right_of_diagonal);
978 template <
typename VectorType>
980 Tvmult(VectorType &dst,
const VectorType &src)
const
982 this->A->precondition_SSOR(dst,
985 pos_right_of_diagonal);
988 template <
typename VectorType,
989 std::enable_if_t<has_SSOR_step<MatrixType, VectorType>,
990 MatrixType> * =
nullptr>
992 step(VectorType &dst,
const VectorType &src)
const
994 this->A->SSOR_step(dst, src, this->relaxation);
997 template <
typename VectorType,
998 std::enable_if_t<!has_SSOR_step<MatrixType, VectorType>,
999 MatrixType> * =
nullptr>
1001 step(VectorType &,
const VectorType &)
const
1005 "Matrix A does not provide a SSOR_step() function!"));
1008 template <
typename VectorType>
1010 Tstep(VectorType &dst,
const VectorType &src)
const
1013 this->step(dst, src);
1018 const double relaxation;
1024 std::vector<std::size_t> pos_right_of_diagonal;
1027 template <
typename MatrixType>
1028 class PreconditionPSORImpl
1031 using size_type =
typename MatrixType::size_type;
1033 PreconditionPSORImpl(
const MatrixType &A,
1034 const double relaxation,
1035 const std::vector<size_type> &permutation,
1036 const std::vector<size_type> &inverse_permutation)
1038 , relaxation(relaxation)
1039 , permutation(permutation)
1040 , inverse_permutation(inverse_permutation)
1043 template <
typename VectorType>
1045 vmult(VectorType &dst,
const VectorType &src)
const
1048 this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
1051 template <
typename VectorType>
1053 Tvmult(VectorType &dst,
const VectorType &src)
const
1056 this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
1061 const double relaxation;
1063 const std::vector<size_type> &permutation;
1064 const std::vector<size_type> &inverse_permutation;
1067 template <
typename MatrixType,
1068 typename PreconditionerType,
1069 typename VectorType,
1070 std::enable_if_t<has_step_omega<PreconditionerType, VectorType>,
1071 PreconditionerType> * =
nullptr>
1073 step(
const MatrixType &,
1074 const PreconditionerType &preconditioner,
1076 const VectorType &src,
1077 const double relaxation,
1081 preconditioner.step(dst, src, relaxation);
1085 typename MatrixType,
1086 typename PreconditionerType,
1087 typename VectorType,
1088 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1089 has_step<PreconditionerType, VectorType>,
1090 PreconditionerType> * =
nullptr>
1092 step(
const MatrixType &,
1093 const PreconditionerType &preconditioner,
1095 const VectorType &src,
1096 const double relaxation,
1104 preconditioner.step(dst, src);
1108 typename MatrixType,
1109 typename PreconditionerType,
1110 typename VectorType,
1111 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1112 !has_step<PreconditionerType, VectorType>,
1113 PreconditionerType> * =
nullptr>
1115 step(
const MatrixType &A,
1116 const PreconditionerType &preconditioner,
1118 const VectorType &src,
1119 const double relaxation,
1120 VectorType &residual,
1123 residual.reinit(dst,
true);
1124 tmp.reinit(dst,
true);
1126 A.vmult(residual, dst);
1127 residual.sadd(-1.0, 1.0, src);
1129 preconditioner.vmult(tmp, residual);
1130 dst.add(relaxation, tmp);
1133 template <
typename MatrixType,
1134 typename PreconditionerType,
1135 typename VectorType,
1136 std::enable_if_t<has_Tstep_omega<PreconditionerType, VectorType>,
1137 PreconditionerType> * =
nullptr>
1139 Tstep(
const MatrixType &,
1140 const PreconditionerType &preconditioner,
1142 const VectorType &src,
1143 const double relaxation,
1147 preconditioner.Tstep(dst, src, relaxation);
1151 typename MatrixType,
1152 typename PreconditionerType,
1153 typename VectorType,
1154 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1155 has_Tstep<PreconditionerType, VectorType>,
1156 PreconditionerType> * =
nullptr>
1158 Tstep(
const MatrixType &,
1159 const PreconditionerType &preconditioner,
1161 const VectorType &src,
1162 const double relaxation,
1170 preconditioner.Tstep(dst, src);
1173 template <
typename MatrixType,
1174 typename VectorType,
1175 std::enable_if_t<has_Tvmult<MatrixType, VectorType>, MatrixType>
1178 Tvmult(
const MatrixType &A, VectorType &dst,
const VectorType &src)
1183 template <
typename MatrixType,
1184 typename VectorType,
1185 std::enable_if_t<!has_Tvmult<MatrixType, VectorType>, MatrixType>
1188 Tvmult(
const MatrixType &, VectorType &,
const VectorType &)
1191 ExcMessage(
"Matrix A does not provide a Tvmult() function!"));
1195 typename MatrixType,
1196 typename PreconditionerType,
1197 typename VectorType,
1198 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1199 !has_Tstep<PreconditionerType, VectorType>,
1200 PreconditionerType> * =
nullptr>
1202 Tstep(
const MatrixType &A,
1203 const PreconditionerType &preconditioner,
1205 const VectorType &src,
1206 const double relaxation,
1207 VectorType &residual,
1210 residual.reinit(dst,
true);
1211 tmp.reinit(dst,
true);
1213 Tvmult(A, residual, dst);
1214 residual.sadd(-1.0, 1.0, src);
1216 Tvmult(preconditioner, tmp, residual);
1217 dst.add(relaxation, tmp);
1221 template <
typename MatrixType,
1222 typename PreconditionerType,
1223 typename VectorType,
1224 std::enable_if_t<!has_vmult_with_std_functions_for_precondition<
1229 step_operations(
const MatrixType &A,
1230 const PreconditionerType &preconditioner,
1232 const VectorType &src,
1233 const double relaxation,
1236 const unsigned int i,
1237 const bool transposed)
1242 Tvmult(preconditioner, dst, src);
1244 preconditioner.vmult(dst, src);
1246 if (relaxation != 1.0)
1252 Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1254 step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1261 typename MatrixType,
1262 typename PreconditionerType,
1263 typename VectorType,
1265 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1267 !has_vmult_with_std_functions_for_precondition<MatrixType,
1271 step_operations(
const MatrixType &A,
1272 const PreconditionerType &preconditioner,
1274 const VectorType &src,
1275 const double relaxation,
1278 const unsigned int i,
1279 const bool transposed)
1282 using Number =
typename VectorType::value_type;
1286 Number *dst_ptr = dst.begin();
1287 const Number *src_ptr = src.begin();
1289 preconditioner.vmult(
1292 [&](
const unsigned int start_range,
const unsigned int end_range) {
1294 if (end_range > start_range)
1295 std::memset(dst.begin() + start_range,
1297 sizeof(Number) * (end_range - start_range));
1299 [&](
const unsigned int start_range,
const unsigned int end_range) {
1300 if (relaxation == 1.0)
1303 const auto src_ptr = src.begin();
1304 const auto dst_ptr = dst.begin();
1307 for (std::size_t i = start_range; i < end_range; ++i)
1308 dst_ptr[i] *= relaxation;
1313 tmp.reinit(src,
true);
1319 preconditioner.vmult(
1322 [&](
const unsigned int start_range,
const unsigned int end_range) {
1323 const auto src_ptr = src.begin();
1324 const auto tmp_ptr = tmp.begin();
1326 if (relaxation == 1.0)
1329 for (std::size_t i = start_range; i < end_range; ++i)
1330 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1338 for (std::size_t i = start_range; i < end_range; ++i)
1339 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1342 [&](
const unsigned int,
const unsigned int) {
1352 typename MatrixType,
1353 typename PreconditionerType,
1354 typename VectorType,
1356 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1358 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1361 step_operations(
const MatrixType &A,
1362 const PreconditionerType &preconditioner,
1364 const VectorType &src,
1365 const double relaxation,
1368 const unsigned int i,
1369 const bool transposed)
1372 using Number =
typename VectorType::value_type;
1376 Number *dst_ptr = dst.begin();
1377 const Number *src_ptr = src.begin();
1379 preconditioner.vmult(
1382 [&](
const unsigned int start_range,
const unsigned int end_range) {
1384 if (end_range > start_range)
1385 std::memset(dst.begin() + start_range,
1387 sizeof(Number) * (end_range - start_range));
1389 [&](
const unsigned int start_range,
const unsigned int end_range) {
1390 if (relaxation == 1.0)
1393 const auto src_ptr = src.begin();
1394 const auto dst_ptr = dst.begin();
1397 for (std::size_t i = start_range; i < end_range; ++i)
1398 dst_ptr[i] *= relaxation;
1403 tmp.reinit(src,
true);
1410 [&](
const unsigned int start_range,
const unsigned int end_range) {
1413 if (end_range > start_range)
1414 std::memset(tmp.begin() + start_range,
1416 sizeof(Number) * (end_range - start_range));
1418 [&](
const unsigned int start_range,
const unsigned int end_range) {
1419 const auto src_ptr = src.begin();
1420 const auto tmp_ptr = tmp.begin();
1422 if (relaxation == 1.0)
1425 for (std::size_t i = start_range; i < end_range; ++i)
1426 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1434 for (std::size_t i = start_range; i < end_range; ++i)
1435 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1439 preconditioner.vmult(dst, tmp, [](
const auto,
const auto) {
1447 template <
typename MatrixType,
1448 typename VectorType,
1449 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1450 !has_vmult_with_std_functions<
1454 VectorType> * =
nullptr>
1456 step_operations(
const MatrixType &A,
1457 const ::DiagonalMatrix<VectorType> &preconditioner,
1459 const VectorType &src,
1460 const double relaxation,
1463 const unsigned int i,
1464 const bool transposed)
1466 using Number =
typename VectorType::value_type;
1470 Number *dst_ptr = dst.begin();
1471 const Number *src_ptr = src.begin();
1472 const Number *diag_ptr = preconditioner.
get_vector().begin();
1474 if (relaxation == 1.0)
1477 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1478 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1483 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1484 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1489 tmp.reinit(src,
true);
1491 Number *dst_ptr = dst.begin();
1492 const Number *src_ptr = src.begin();
1493 const Number *tmp_ptr = tmp.begin();
1494 const Number *diag_ptr = preconditioner.get_vector().begin();
1497 Tvmult(A, tmp, dst);
1501 if (relaxation == 1.0)
1504 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1505 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1510 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1512 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1519 template <
typename MatrixType,
1520 typename VectorType,
1521 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1522 has_vmult_with_std_functions<
1526 VectorType> * =
nullptr>
1528 step_operations(
const MatrixType &A,
1529 const ::DiagonalMatrix<VectorType> &preconditioner,
1531 const VectorType &src,
1532 const double relaxation,
1535 const unsigned int i,
1536 const bool transposed)
1539 using Number =
typename VectorType::value_type;
1543 Number *dst_ptr = dst.begin();
1544 const Number *src_ptr = src.begin();
1545 const Number *diag_ptr = preconditioner.get_vector().begin();
1547 if (relaxation == 1.0)
1550 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1551 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1556 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1557 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1562 tmp.reinit(src,
true);
1569 [&](
const unsigned int start_range,
const unsigned int end_range) {
1571 if (end_range > start_range)
1572 std::memset(tmp.begin() + start_range,
1574 sizeof(Number) * (end_range - start_range));
1576 [&](
const unsigned int begin,
const unsigned int end) {
1577 const Number *dst_ptr = dst.begin();
1578 const Number *src_ptr = src.begin();
1579 Number *tmp_ptr = tmp.begin();
1580 const Number *diag_ptr = preconditioner.get_vector().begin();
1584 if (relaxation == 1.0)
1587 for (std::size_t i = begin; i < end; ++i)
1589 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1594 for (std::size_t i = begin; i < end; ++i)
1595 tmp_ptr[i] = dst_ptr[i] + relaxation *
1596 (src_ptr[i] - tmp_ptr[i]) *
1638template <
typename MatrixType = SparseMatrix<
double>>
1642 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1645 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1708template <
typename MatrixType = SparseMatrix<
double>>
1712 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1715 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1760template <
typename MatrixType = SparseMatrix<
double>>
1764 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1767 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1816template <
typename MatrixType = SparseMatrix<
double>>
1820 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1823 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1880 const std::vector<size_type> &permutation,
1881 const std::vector<size_type> &inverse_permutation,
2098template <
typename MatrixType = SparseMatrix<
double>,
2099 typename VectorType = Vector<
double>,
2100 typename PreconditionerType = DiagonalMatrix<VectorType>>
2138 const unsigned int degree = 1,
2144 EigenvalueAlgorithm::lanczos,
2193 vmult(VectorType &dst,
const VectorType &src)
const;
2200 Tvmult(VectorType &dst,
const VectorType &src)
const;
2206 step(VectorType &dst,
const VectorType &src)
const;
2212 Tstep(VectorType &dst,
const VectorType &src)
const;
2315 template <
typename VectorType>
2317 set_initial_guess(VectorType &vector)
2319 vector = 1. /
std::sqrt(
static_cast<double>(vector.size()));
2320 if (vector.locally_owned_elements().is_element(0))
2324 template <
typename Number>
2332 for (
unsigned int i = 0; i < vector.
size(); ++i)
2335 const Number mean_value = vector.
mean_value();
2336 vector.
add(-mean_value);
2339 template <
typename Number>
2344 for (
unsigned int block = 0; block < vector.
n_blocks(); ++block)
2345 set_initial_guess(vector.
block(block));
2348 template <
typename Number,
typename MemorySpace>
2365 Kokkos::RangePolicy<
typename MemorySpace::kokkos_space::execution_space,
2366 Kokkos::IndexType<types::global_dof_index>>
2367 policy(0, n_local_elements);
2368 Kokkos::parallel_for(
2369 "::PreconditionChebyshev::set_initial_guess",
2372 values_ptr[i] = (i + first_local_range) % 11;
2374 const Number mean_value = vector.
mean_value();
2375 vector.
add(-mean_value);
2378 struct EigenvalueTracker
2387 std::vector<double>
values;
2392 template <
typename MatrixType,
2393 typename VectorType,
2394 typename PreconditionerType>
2397 VectorType &eigenvector,
2398 const PreconditionerType &preconditioner,
2399 const unsigned int n_iterations)
2401 typename VectorType::value_type eigenvalue_estimate = 0.;
2402 eigenvector /= eigenvector.l2_norm();
2403 VectorType vector1, vector2;
2404 vector1.reinit(eigenvector,
true);
2405 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2406 vector2.reinit(eigenvector,
true);
2407 for (
unsigned int i = 0; i < n_iterations; ++i)
2409 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2411 matrix.vmult(vector2, eigenvector);
2412 preconditioner.vmult(vector1, vector2);
2415 matrix.vmult(vector1, eigenvector);
2417 eigenvalue_estimate = eigenvector * vector1;
2419 vector1 /= vector1.l2_norm();
2420 eigenvector.swap(vector1);
2422 return std::abs(eigenvalue_estimate);
2427 template <
typename MatrixType,
2428 typename VectorType,
2429 typename PreconditionerType>
2430 EigenvalueInformation
2431 estimate_eigenvalues(
2432 const EigenvalueAlgorithmAdditionalData<PreconditionerType> &data,
2433 const MatrixType *matrix_ptr,
2434 VectorType &solution_old,
2435 VectorType &temp_vector1,
2436 const unsigned int degree)
2440 EigenvalueInformation info{};
2442 if (data.eig_cg_n_iterations > 0)
2444 Assert(data.eig_cg_n_iterations > 2,
2446 "Need to set at least two iterations to find eigenvalues."));
2448 internal::EigenvalueTracker eigenvalue_tracker;
2453 internal::set_initial_guess(temp_vector1);
2454 data.constraints.set_zero(temp_vector1);
2465 solver.connect_eigenvalues_slot(
2466 [&eigenvalue_tracker](
const std::vector<double> &
eigenvalues) {
2470 solver.solve(*matrix_ptr,
2473 *data.preconditioner);
2475 info.cg_iterations = control.last_step();
2477 else if (data.eigenvalue_algorithm ==
2483 ExcMessage(
"Cannot estimate the minimal eigenvalue with the "
2484 "power iteration"));
2486 eigenvalue_tracker.values.push_back(
2489 *data.preconditioner,
2490 data.eig_cg_n_iterations));
2496 if (eigenvalue_tracker.values.empty())
2497 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2500 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2504 info.max_eigenvalue_estimate =
2505 1.2 * eigenvalue_tracker.values.back();
2510 info.max_eigenvalue_estimate = data.max_eigenvalue;
2511 info.min_eigenvalue_estimate =
2512 data.max_eigenvalue / data.smoothing_range;
2526template <
typename MatrixType>
2536template <
typename VectorType>
2545template <
typename VectorType>
2552template <
typename VectorType>
2561template <
typename VectorType>
2593 const double relaxation)
2594 : relaxation(relaxation)
2603 AdditionalData add_data;
2604 relaxation = add_data.relaxation;
2618template <
typename MatrixType>
2621 const MatrixType &matrix,
2631template <
typename VectorType>
2636 std::is_same_v<size_type, typename VectorType::size_type>,
2637 "PreconditionRichardson and VectorType must have the same size_type.");
2644template <
typename VectorType>
2649 std::is_same_v<size_type, typename VectorType::size_type>,
2650 "PreconditionRichardson and VectorType must have the same size_type.");
2655template <
typename VectorType>
2660 std::is_same_v<size_type, typename VectorType::size_type>,
2661 "PreconditionRichardson and VectorType must have the same size_type.");
2668template <
typename VectorType>
2673 std::is_same_v<size_type, typename VectorType::size_type>,
2674 "PreconditionRichardson and VectorType must have the same size_type.");
2695template <
typename MatrixType,
typename PreconditionerType>
2698 const MatrixType &rA,
2699 const AdditionalData ¶meters)
2702 eigenvalues_are_initialized =
false;
2706 this->data = parameters;
2710template <
typename MatrixType,
typename PreconditionerType>
2714 eigenvalues_are_initialized =
false;
2716 data.relaxation = 1.0;
2717 data.preconditioner =
nullptr;
2720template <
typename MatrixType,
typename PreconditionerType>
2729template <
typename MatrixType,
typename PreconditionerType>
2738template <
typename MatrixType,
typename PreconditionerType>
2739template <
typename VectorType>
2743 const VectorType &src)
const
2748 if (eigenvalues_are_initialized ==
false)
2749 estimate_eigenvalues(src);
2751 VectorType tmp1, tmp2;
2753 for (
unsigned int i = 0; i < data.n_iterations; ++i)
2754 internal::PreconditionRelaxation::step_operations(*A,
2755 *data.preconditioner,
2765template <
typename MatrixType,
typename PreconditionerType>
2766template <
typename VectorType>
2770 const VectorType &src)
const
2775 if (eigenvalues_are_initialized ==
false)
2776 estimate_eigenvalues(src);
2778 VectorType tmp1, tmp2;
2780 for (
unsigned int i = 0; i < data.n_iterations; ++i)
2781 internal::PreconditionRelaxation::step_operations(
2782 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i,
true);
2785template <
typename MatrixType,
typename PreconditionerType>
2786template <
typename VectorType>
2790 const VectorType &src)
const
2795 if (eigenvalues_are_initialized ==
false)
2796 estimate_eigenvalues(src);
2798 VectorType tmp1, tmp2;
2800 for (
unsigned int i = 1; i <= data.n_iterations; ++i)
2801 internal::PreconditionRelaxation::step_operations(*A,
2802 *data.preconditioner,
2812template <
typename MatrixType,
typename PreconditionerType>
2813template <
typename VectorType>
2817 const VectorType &src)
const
2822 if (eigenvalues_are_initialized ==
false)
2823 estimate_eigenvalues(src);
2825 VectorType tmp1, tmp2;
2827 for (
unsigned int i = 1; i <= data.n_iterations; ++i)
2828 internal::PreconditionRelaxation::step_operations(
2829 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i,
true);
2832template <
typename MatrixType,
typename PreconditionerType>
2833template <
typename VectorType>
2836 const VectorType &src)
const
2840 EigenvalueInformation info;
2842 if (data.relaxation == 0.0)
2844 VectorType solution_old, temp_vector1;
2846 solution_old.reinit(src);
2847 temp_vector1.reinit(src,
true);
2849 info = internal::estimate_eigenvalues<MatrixType>(
2850 data, A, solution_old, temp_vector1, data.n_iterations);
2852 const double alpha =
2853 (data.smoothing_range > 1. ?
2854 info.max_eigenvalue_estimate / data.smoothing_range :
2855 std::min(0.9 * info.max_eigenvalue_estimate,
2856 info.min_eigenvalue_estimate));
2859 ->data.relaxation = 2.0 / (alpha + info.max_eigenvalue_estimate);
2863 ->eigenvalues_are_initialized =
true;
2868template <
typename MatrixType,
typename PreconditionerType>
2872 return data.relaxation;
2878template <
typename MatrixType>
2881 const AdditionalData ¶meters_in)
2885 parameters_in.relaxation != 0.0,
2887 "Relaxation cannot automatically be determined by PreconditionJacobi."));
2889 AdditionalData parameters;
2890 parameters.relaxation = 1.0;
2891 parameters.n_iterations = parameters_in.n_iterations;
2892 parameters.preconditioner =
2893 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2895 this->BaseClass::initialize(A, parameters);
2900template <
typename MatrixType>
2903 const AdditionalData ¶meters_in)
2907 parameters_in.relaxation != 0.0,
2909 "Relaxation cannot automatically be determined by PreconditionSOR."));
2911 AdditionalData parameters;
2912 parameters.relaxation = 1.0;
2913 parameters.n_iterations = parameters_in.n_iterations;
2914 parameters.preconditioner =
2915 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2917 this->BaseClass::initialize(A, parameters);
2922template <
typename MatrixType>
2925 const AdditionalData ¶meters_in)
2929 parameters_in.relaxation != 0.0,
2931 "Relaxation cannot automatically be determined by PreconditionSSOR."));
2933 AdditionalData parameters;
2934 parameters.relaxation = 1.0;
2935 parameters.n_iterations = parameters_in.n_iterations;
2936 parameters.preconditioner =
2937 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2939 this->BaseClass::initialize(A, parameters);
2946template <
typename MatrixType>
2949 const MatrixType &A,
2950 const std::vector<size_type> &p,
2951 const std::vector<size_type> &ip,
2952 const typename BaseClass::AdditionalData ¶meters_in)
2956 parameters_in.relaxation != 0.0,
2958 "Relaxation cannot automatically be determined by PreconditionPSOR."));
2960 typename BaseClass::AdditionalData parameters;
2961 parameters.relaxation = 1.0;
2962 parameters.n_iterations = parameters_in.n_iterations;
2963 parameters.preconditioner =
2964 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2966 this->BaseClass::initialize(A, parameters);
2970template <
typename MatrixType>
2973 const AdditionalData &additional_data)
2976 additional_data.permutation,
2977 additional_data.inverse_permutation,
2978 additional_data.parameters);
2981template <
typename MatrixType>
2983 const std::vector<size_type> &permutation,
2984 const std::vector<size_type> &inverse_permutation,
2987 : permutation(permutation)
2988 , inverse_permutation(inverse_permutation)
2989 , parameters(parameters)
2996template <
typename MatrixType,
typename VectorType>
2998 const MatrixType &M,
2999 const function_ptr method)
3001 , precondition(method)
3006template <
typename MatrixType,
typename VectorType>
3010 const VectorType &src)
const
3012 (
matrix.*precondition)(dst, src);
3020 template <
typename PreconditionerType>
3023 const double smoothing_range,
3024 const unsigned int eig_cg_n_iterations,
3025 const double eig_cg_residual,
3026 const double max_eigenvalue,
3028 : smoothing_range(smoothing_range)
3029 , eig_cg_n_iterations(eig_cg_n_iterations)
3030 , eig_cg_residual(eig_cg_residual)
3031 , max_eigenvalue(max_eigenvalue)
3032 , eigenvalue_algorithm(eigenvalue_algorithm)
3037 template <
typename PreconditionerType>
3038 inline EigenvalueAlgorithmAdditionalData<PreconditionerType> &
3039 EigenvalueAlgorithmAdditionalData<PreconditionerType>::operator=(
3040 const EigenvalueAlgorithmAdditionalData &other_data)
3042 smoothing_range = other_data.smoothing_range;
3043 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3044 eig_cg_residual = other_data.eig_cg_residual;
3045 max_eigenvalue = other_data.max_eigenvalue;
3046 preconditioner = other_data.preconditioner;
3047 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3048 constraints.copy_from(other_data.constraints);
3054template <
typename MatrixType,
typename PreconditionerType>
3057 const unsigned int n_iterations,
3058 const double smoothing_range,
3059 const unsigned int eig_cg_n_iterations,
3060 const double eig_cg_residual,
3061 const double max_eigenvalue,
3062 const EigenvalueAlgorithm eigenvalue_algorithm)
3063 :
internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3065 eig_cg_n_iterations,
3068 eigenvalue_algorithm)
3069 , relaxation(relaxation)
3070 , n_iterations(n_iterations)
3079 namespace PreconditionChebyshevImplementation
3087 template <
typename VectorType,
typename PreconditionerType>
3089 vector_updates(
const VectorType &rhs,
3090 const PreconditionerType &preconditioner,
3091 const unsigned int iteration_index,
3092 const double factor1,
3093 const double factor2,
3094 VectorType &solution_old,
3095 VectorType &temp_vector1,
3096 VectorType &temp_vector2,
3097 VectorType &solution)
3099 if (iteration_index == 0)
3101 solution.equ(factor2, rhs);
3102 preconditioner.vmult(solution_old, solution);
3104 else if (iteration_index == 1)
3107 temp_vector1.sadd(-1.0, 1.0, rhs);
3108 preconditioner.vmult(solution_old, temp_vector1);
3111 solution_old.sadd(factor2, 1 + factor1, solution);
3116 temp_vector1.sadd(-1.0, 1.0, rhs);
3117 preconditioner.vmult(temp_vector2, temp_vector1);
3120 solution_old.sadd(-factor1, factor2, temp_vector2);
3121 solution_old.add(1 + factor1, solution);
3124 solution.swap(solution_old);
3130 typename PreconditionerType,
3132 !has_vmult_with_std_functions_for_precondition<
3139 const PreconditionerType &preconditioner,
3140 const unsigned int iteration_index,
3141 const double factor1_,
3142 const double factor2_,
3151 const Number factor1 = factor1_;
3152 const Number factor1_plus_1 = 1. + factor1_;
3153 const Number factor2 = factor2_;
3155 if (iteration_index == 0)
3157 const auto solution_old_ptr = solution_old.
begin();
3160 preconditioner.vmult(solution_old, rhs);
3165 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
3167 else if (iteration_index == 1)
3169 const auto solution_ptr = solution.
begin();
3170 const auto solution_old_ptr = solution_old.
begin();
3173 temp_vector1.
sadd(-1.0, 1.0, rhs);
3175 preconditioner.vmult(solution_old, temp_vector1);
3180 solution_old_ptr[i] =
3181 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
3185 const auto solution_ptr = solution.
begin();
3186 const auto solution_old_ptr = solution_old.
begin();
3187 const auto temp_vector2_ptr = temp_vector2.
begin();
3190 temp_vector1.
sadd(-1.0, 1.0, rhs);
3192 preconditioner.vmult(temp_vector2, temp_vector1);
3197 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3198 factor1 * solution_old_ptr[i] +
3199 temp_vector2_ptr[i] * factor2;
3202 solution.
swap(solution_old);
3207 typename PreconditionerType,
3209 has_vmult_with_std_functions_for_precondition<
3216 const PreconditionerType &preconditioner,
3217 const unsigned int iteration_index,
3218 const double factor1_,
3219 const double factor2_,
3228 const Number factor1 = factor1_;
3229 const Number factor1_plus_1 = 1. + factor1_;
3230 const Number factor2 = factor2_;
3232 const auto rhs_ptr = rhs.
begin();
3233 const auto temp_vector1_ptr = temp_vector1.
begin();
3234 const auto temp_vector2_ptr = temp_vector2.
begin();
3235 const auto solution_ptr = solution.
begin();
3236 const auto solution_old_ptr = solution_old.
begin();
3238 if (iteration_index == 0)
3240 preconditioner.vmult(
3243 [&](
const auto start_range,
const auto end_range) {
3244 if (end_range > start_range)
3245 std::memset(solution.
begin() + start_range,
3247 sizeof(Number) * (end_range - start_range));
3249 [&](
const auto begin,
const auto end) {
3251 for (std::size_t i = begin; i <
end; ++i)
3252 solution_ptr[i] *= factor2;
3257 preconditioner.vmult(
3260 [&](
const auto begin,
const auto end) {
3262 std::memset(temp_vector2.
begin() + begin,
3264 sizeof(Number) * (end - begin));
3267 for (std::size_t i = begin; i <
end; ++i)
3268 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3270 [&](
const auto begin,
const auto end) {
3271 if (iteration_index == 1)
3274 for (std::size_t i = begin; i <
end; ++i)
3275 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
3276 factor2 * temp_vector2_ptr[i];
3281 for (std::size_t i = begin; i <
end; ++i)
3282 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3283 factor1 * solution_old_ptr[i] +
3284 factor2 * temp_vector2_ptr[i];
3289 if (iteration_index > 0)
3291 solution_old.
swap(temp_vector2);
3292 solution_old.
swap(solution);
3299 template <
typename Number>
3300 struct VectorUpdater
3302 VectorUpdater(
const Number *rhs,
3303 const Number *matrix_diagonal_inverse,
3304 const unsigned int iteration_index,
3305 const Number factor1,
3306 const Number factor2,
3307 Number *solution_old,
3311 , matrix_diagonal_inverse(matrix_diagonal_inverse)
3312 , iteration_index(iteration_index)
3315 , solution_old(solution_old)
3316 , tmp_vector(tmp_vector)
3317 , solution(solution)
3321 apply_to_subrange(
const std::size_t begin,
const std::size_t end)
const
3327 const Number factor1 = this->factor1;
3328 const Number factor1_plus_1 = 1. + this->factor1;
3329 const Number factor2 = this->factor2;
3330 if (iteration_index == 0)
3333 for (std::size_t i = begin; i <
end; ++i)
3334 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
3336 else if (iteration_index == 1)
3340 for (std::size_t i = begin; i <
end; ++i)
3344 factor1_plus_1 * solution[i] +
3345 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3352 for (std::size_t i = begin; i <
end; ++i)
3359 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
3360 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3365 const Number *matrix_diagonal_inverse;
3366 const unsigned int iteration_index;
3367 const Number factor1;
3368 const Number factor2;
3369 mutable Number *solution_old;
3370 mutable Number *tmp_vector;
3371 mutable Number *solution;
3374 template <
typename Number>
3377 VectorUpdatesRange(
const VectorUpdater<Number> &updater,
3378 const std::size_t size)
3382 VectorUpdatesRange::apply_to_subrange(0, size);
3390 ~VectorUpdatesRange()
override =
default;
3393 apply_to_subrange(
const std::size_t begin,
3394 const std::size_t end)
const override
3396 updater.apply_to_subrange(begin, end);
3399 const VectorUpdater<Number> &updater;
3403 template <
typename Number>
3406 const ::Vector<Number> &rhs,
3408 const unsigned int iteration_index,
3409 const double factor1,
3410 const double factor2,
3416 VectorUpdater<Number> upd(rhs.begin(),
3417 jacobi.get_vector().begin(),
3421 solution_old.
begin(),
3422 temp_vector1.
begin(),
3424 VectorUpdatesRange<Number>(upd, rhs.size());
3427 if (iteration_index == 0)
3434 solution.
swap(temp_vector1);
3435 solution_old.
swap(temp_vector1);
3440 template <
typename Number>
3444 const ::DiagonalMatrix<
3446 const unsigned int iteration_index,
3447 const double factor1,
3448 const double factor2,
3456 VectorUpdater<Number> upd(rhs.
begin(),
3457 jacobi.get_vector().begin(),
3461 solution_old.
begin(),
3462 temp_vector1.
begin(),
3467 if (iteration_index == 0)
3474 solution.
swap(temp_vector1);
3475 solution_old.
swap(temp_vector1);
3484 typename MatrixType,
3485 typename VectorType,
3486 typename PreconditionerType,
3488 !has_vmult_with_std_functions<MatrixType,
3490 PreconditionerType> &&
3491 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3493 has_vmult_with_std_functions_for_precondition<MatrixType,
3497 vmult_and_update(
const MatrixType &matrix,
3498 const PreconditionerType &preconditioner,
3499 const VectorType &rhs,
3500 const unsigned int iteration_index,
3501 const double factor1,
3502 const double factor2,
3503 VectorType &solution,
3504 VectorType &solution_old,
3505 VectorType &temp_vector1,
3506 VectorType &temp_vector2)
3508 if (iteration_index > 0)
3509 matrix.vmult(temp_vector1, solution);
3524 typename MatrixType,
3525 typename VectorType,
3526 typename PreconditionerType,
3528 !has_vmult_with_std_functions<MatrixType,
3530 PreconditionerType> &&
3531 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3533 has_vmult_with_std_functions_for_precondition<MatrixType,
3537 vmult_and_update(
const MatrixType &matrix,
3538 const PreconditionerType &preconditioner,
3539 const VectorType &rhs,
3540 const unsigned int iteration_index,
3541 const double factor1_,
3542 const double factor2_,
3543 VectorType &solution,
3544 VectorType &solution_old,
3545 VectorType &temp_vector1,
3546 VectorType &temp_vector2)
3548 using Number =
typename VectorType::value_type;
3550 const Number factor1 = factor1_;
3551 const Number factor1_plus_1 = 1. + factor1_;
3552 const Number factor2 = factor2_;
3554 if (iteration_index == 0)
3556 preconditioner.vmult(
3559 [&](
const unsigned int start_range,
const unsigned int end_range) {
3561 if (end_range > start_range)
3562 std::memset(solution.begin() + start_range,
3564 sizeof(Number) * (end_range - start_range));
3566 [&](
const unsigned int start_range,
const unsigned int end_range) {
3567 const auto solution_ptr = solution.begin();
3570 for (std::size_t i = start_range; i < end_range; ++i)
3571 solution_ptr[i] *= factor2;
3576 temp_vector1.reinit(rhs,
true);
3577 temp_vector2.reinit(rhs,
true);
3583 [&](
const unsigned int start_range,
const unsigned int end_range) {
3586 if (end_range > start_range)
3587 std::memset(temp_vector1.begin() + start_range,
3589 sizeof(Number) * (end_range - start_range));
3591 [&](
const unsigned int start_range,
const unsigned int end_range) {
3592 const auto rhs_ptr = rhs.begin();
3593 const auto tmp_ptr = temp_vector1.begin();
3596 for (std::size_t i = start_range; i < end_range; ++i)
3597 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3601 preconditioner.vmult(
3604 [&](
const unsigned int start_range,
const unsigned int end_range) {
3607 if (end_range > start_range)
3608 std::memset(temp_vector2.begin() + start_range,
3610 sizeof(Number) * (end_range - start_range));
3612 [&](
const unsigned int start_range,
const unsigned int end_range) {
3613 const auto solution_ptr = solution.begin();
3614 const auto solution_old_ptr = solution_old.begin();
3615 const auto tmp_ptr = temp_vector2.begin();
3617 if (iteration_index == 1)
3620 for (std::size_t i = start_range; i < end_range; ++i)
3622 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3627 for (std::size_t i = start_range; i < end_range; ++i)
3628 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3629 factor1 * solution_old_ptr[i] +
3630 factor2 * tmp_ptr[i];
3634 solution.swap(temp_vector2);
3635 solution_old.swap(temp_vector2);
3641 template <
typename MatrixType,
3642 typename VectorType,
3643 typename PreconditionerType,
3644 std::enable_if_t<has_vmult_with_std_functions<MatrixType,
3646 PreconditionerType>,
3649 vmult_and_update(
const MatrixType &matrix,
3650 const PreconditionerType &preconditioner,
3651 const VectorType &rhs,
3652 const unsigned int iteration_index,
3653 const double factor1,
3654 const double factor2,
3655 VectorType &solution,
3656 VectorType &solution_old,
3657 VectorType &temp_vector1,
3660 using Number =
typename VectorType::value_type;
3661 VectorUpdater<Number> updater(rhs.begin(),
3662 preconditioner.get_vector().begin(),
3666 solution_old.begin(),
3667 temp_vector1.begin(),
3669 if (iteration_index > 0)
3673 [&](
const unsigned int start_range,
const unsigned int end_range) {
3676 if (end_range > start_range)
3677 std::memset(temp_vector1.begin() + start_range,
3679 sizeof(Number) * (end_range - start_range));
3681 [&](
const unsigned int start_range,
const unsigned int end_range) {
3682 if (end_range > start_range)
3683 updater.apply_to_subrange(start_range, end_range);
3686 updater.apply_to_subrange(0U, solution.locally_owned_size());
3689 if (iteration_index == 0)
3696 solution.swap(temp_vector1);
3697 solution_old.swap(temp_vector1);
3701 template <
typename MatrixType,
typename PreconditionerType>
3703 initialize_preconditioner(
3704 const MatrixType &matrix,
3705 std::shared_ptr<PreconditionerType> &preconditioner)
3708 (void)preconditioner;
3712 template <
typename MatrixType,
typename VectorType>
3714 initialize_preconditioner(
3715 const MatrixType &matrix,
3718 if (preconditioner.get() ==
nullptr || preconditioner->m() !=
matrix.m())
3720 if (preconditioner.get() ==
nullptr)
3722 std::make_shared<::DiagonalMatrix<VectorType>>();
3725 preconditioner->m() == 0,
3727 "Preconditioner appears to be initialized but not sized correctly"));
3730 if (preconditioner->m() !=
matrix.m())
3732 preconditioner->get_vector().reinit(
matrix.m());
3733 for (
typename VectorType::size_type i = 0; i <
matrix.m(); ++i)
3734 preconditioner->get_vector()(i) = 1. /
matrix.el(i, i);
3743template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3746 const double smoothing_range,
3747 const unsigned int eig_cg_n_iterations,
3748 const double eig_cg_residual,
3749 const double max_eigenvalue,
3750 const EigenvalueAlgorithm eigenvalue_algorithm,
3751 const PolynomialType polynomial_type)
3752 :
internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3754 eig_cg_n_iterations,
3757 eigenvalue_algorithm)
3759 , polynomial_type(polynomial_type)
3764template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3769 , eigenvalues_are_initialized(false)
3772 std::is_same_v<size_type, typename VectorType::size_type>,
3773 "PreconditionChebyshev and VectorType must have the same size_type.");
3778template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3781 const MatrixType &matrix,
3782 const AdditionalData &additional_data)
3785 data = additional_data;
3787 ExcMessage(
"The degree of the Chebyshev method must be positive."));
3788 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3789 matrix, data.preconditioner);
3790 eigenvalues_are_initialized =
false;
3795template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3799 eigenvalues_are_initialized =
false;
3800 theta = delta = 1.0;
3801 matrix_ptr =
nullptr;
3803 VectorType empty_vector;
3804 solution_old.reinit(empty_vector);
3805 temp_vector1.reinit(empty_vector);
3806 temp_vector2.reinit(empty_vector);
3808 data.preconditioner.reset();
3813template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3820 solution_old.reinit(src);
3821 temp_vector1.reinit(src,
true);
3823 auto info = internal::estimate_eigenvalues<MatrixType>(
3824 data, matrix_ptr, solution_old, temp_vector1, data.degree);
3826 const double alpha = (data.smoothing_range > 1. ?
3827 info.max_eigenvalue_estimate / data.smoothing_range :
3828 std::min(0.9 * info.max_eigenvalue_estimate,
3829 info.min_eigenvalue_estimate));
3838 const double actual_range = info.max_eigenvalue_estimate / alpha;
3839 const double sigma = (1. -
std::sqrt(1. / actual_range)) /
3841 const double eps = data.smoothing_range;
3846 1 +
static_cast<unsigned int>(
3851 info.degree = data.degree;
3856 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3857 (info.max_eigenvalue_estimate) :
3858 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3861 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3865 using NumberType =
typename VectorType::value_type;
3873 (std::is_same_v<VectorType,
3876 temp_vector2.
reinit(src, true);
3879 VectorType empty_vector;
3880 temp_vector2.reinit(empty_vector);
3885 ->eigenvalues_are_initialized =
true;
3892template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3895 VectorType &solution,
3896 const VectorType &rhs)
const
3898 std::lock_guard<std::mutex> lock(mutex);
3899 if (eigenvalues_are_initialized ==
false)
3900 estimate_eigenvalues(rhs);
3902 internal::PreconditionChebyshevImplementation::vmult_and_update(
3904 *data.preconditioner,
3908 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3909 (4. / (3. * delta)) :
3918 if (data.degree < 2 ||
std::abs(delta) < 1e-40)
3921 double rhok = delta / theta, sigma = theta / delta;
3922 for (
unsigned int k = 0; k < data.degree - 1; ++k)
3924 double factor1 = 0.0;
3925 double factor2 = 0.0;
3927 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3929 factor1 = (2 * k + 1.) / (2 * k + 5.);
3930 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3934 const double rhokp = 1. / (2. * sigma - rhok);
3935 factor1 = rhokp * rhok;
3936 factor2 = 2. * rhokp / delta;
3940 internal::PreconditionChebyshevImplementation::vmult_and_update(
3942 *data.preconditioner,
3956template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3959 VectorType &solution,
3960 const VectorType &rhs)
const
3962 std::lock_guard<std::mutex> lock(mutex);
3963 if (eigenvalues_are_initialized ==
false)
3964 estimate_eigenvalues(rhs);
3966 internal::PreconditionChebyshevImplementation::vector_updates(
3968 *data.preconditioner,
3971 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3972 (4. / (3. * delta)) :
3979 if (data.degree < 2 ||
std::abs(delta) < 1e-40)
3982 double rhok = delta / theta, sigma = theta / delta;
3983 for (
unsigned int k = 0; k < data.degree - 1; ++k)
3985 double factor1 = 0.0;
3986 double factor2 = 0.0;
3988 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3990 factor1 = (2 * k + 1.) / (2 * k + 5.);
3991 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3995 const double rhokp = 1. / (2. * sigma - rhok);
3996 factor1 = rhokp * rhok;
3997 factor2 = 2. * rhokp / delta;
4001 matrix_ptr->Tvmult(temp_vector1, solution);
4002 internal::PreconditionChebyshevImplementation::vector_updates(
4004 *data.preconditioner,
4017template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4020 VectorType &solution,
4021 const VectorType &rhs)
const
4023 std::lock_guard<std::mutex> lock(mutex);
4024 if (eigenvalues_are_initialized ==
false)
4025 estimate_eigenvalues(rhs);
4027 internal::PreconditionChebyshevImplementation::vmult_and_update(
4029 *data.preconditioner,
4033 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4034 (4. / (3. * delta)) :
4041 if (data.degree < 2 ||
std::abs(delta) < 1e-40)
4044 double rhok = delta / theta, sigma = theta / delta;
4045 for (
unsigned int k = 0; k < data.degree - 1; ++k)
4047 double factor1 = 0.0;
4048 double factor2 = 0.0;
4050 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4052 factor1 = (2 * k + 1.) / (2 * k + 5.);
4053 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4057 const double rhokp = 1. / (2. * sigma - rhok);
4058 factor1 = rhokp * rhok;
4059 factor2 = 2. * rhokp / delta;
4063 internal::PreconditionChebyshevImplementation::vmult_and_update(
4065 *data.preconditioner,
4079template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4082 VectorType &solution,
4083 const VectorType &rhs)
const
4085 std::lock_guard<std::mutex> lock(mutex);
4086 if (eigenvalues_are_initialized ==
false)
4087 estimate_eigenvalues(rhs);
4089 matrix_ptr->Tvmult(temp_vector1, solution);
4090 internal::PreconditionChebyshevImplementation::vector_updates(
4092 *data.preconditioner,
4095 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4096 (4. / (3. * delta)) :
4103 if (data.degree < 2 ||
std::abs(delta) < 1e-40)
4106 double rhok = delta / theta, sigma = theta / delta;
4107 for (
unsigned int k = 0; k < data.degree - 1; ++k)
4109 double factor1 = 0.0;
4110 double factor2 = 0.0;
4112 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4114 factor1 = (2 * k + 1.) / (2 * k + 5.);
4115 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4119 const double rhokp = 1. / (2. * sigma - rhok);
4120 factor1 = rhokp * rhok;
4121 factor2 = 2. * rhokp / delta;
4125 matrix_ptr->Tvmult(temp_vector1, solution);
4126 internal::PreconditionChebyshevImplementation::vector_updates(
4128 *data.preconditioner,
4141template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4144 PreconditionerType>::size_type
4148 return matrix_ptr->m();
4153template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4156 PreconditionerType>::size_type
4160 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
Number * get_values() const
void swap(Vector< Number, MemorySpace > &v)
size_type locally_owned_size() const
::IndexSet locally_owned_elements() 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
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
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
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)
double get_relaxation() const
SmartPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
void Tvmult(VectorType &, const VectorType &) const
std::shared_ptr< PreconditionerType > preconditioner
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
void step(VectorType &x, const VectorType &rhs) const
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
bool eigenvalues_are_initialized
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
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 & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
@ matrix
Contents is actually a matrix.
unsigned int minimum_parallel_grain_size
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static const unsigned int invalid_unsigned_int
const Iterator const std_cxx20::type_identity_t< Iterator > & end
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
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
EigenvalueAlgorithmAdditionalData< PreconditionerType > & operator=(const EigenvalueAlgorithmAdditionalData< PreconditionerType > &other_data)
::AffineConstraints< double > constraints
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)
unsigned int eig_cg_n_iterations
EigenvalueAlgorithm eigenvalue_algorithm
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)