15#ifndef dealii_precondition_h
16#define dealii_precondition_h
33#include <Kokkos_Core.hpp>
41template <
typename number>
43template <
typename number>
49 template <
typename,
typename>
51 template <
typename,
typename>
128 template <
typename PreconditionerType>
249 template <
typename MatrixType>
257 template <
typename VectorType>
259 vmult(VectorType &,
const VectorType &)
const;
265 template <
typename VectorType>
267 Tvmult(VectorType &,
const VectorType &)
const;
272 template <
typename VectorType>
280 template <
typename VectorType>
378 template <
typename MatrixType>
385 template <
typename VectorType>
387 vmult(VectorType &,
const VectorType &)
const;
393 template <
typename VectorType>
395 Tvmult(VectorType &,
const VectorType &)
const;
399 template <
typename VectorType>
407 template <
typename VectorType>
497template <
typename MatrixType = SparseMatrix<
double>,
498 typename VectorType = Vector<
double>>
506 const VectorType &) const;
520 vmult(VectorType &dst,
const VectorType &src)
const;
565template <
typename MatrixType = SparseMatrix<
double>,
566 typename PreconditionerType = IdentityMatrix>
594 EigenvalueAlgorithm::lanczos);
640 template <
typename VectorType>
642 vmult(VectorType &,
const VectorType &)
const;
648 template <
typename VectorType>
650 Tvmult(VectorType &,
const VectorType &)
const;
655 template <
typename VectorType>
657 step(VectorType &x,
const VectorType &rhs)
const;
662 template <
typename VectorType>
664 Tstep(VectorType &x,
const VectorType &rhs)
const;
671 template <
typename VectorType>
716 template <
typename MatrixType,
typename VectorType>
717 using vmult_functions_t =
decltype(std::declval<const MatrixType>().vmult(
718 std::declval<VectorType &>(),
719 std::declval<const VectorType &>(),
721 const std::function<
void(
const unsigned int,
const unsigned int)> &>(),
723 const std::function<
void(
const unsigned int,
const unsigned int)> &>()));
725 template <
typename MatrixType,
727 typename PreconditionerType>
728 constexpr bool has_vmult_with_std_functions =
729 is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
730 std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> &&
731 (std::is_same_v<VectorType,
739 template <
typename MatrixType,
typename VectorType>
740 constexpr bool has_vmult_with_std_functions_for_precondition =
741 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
745 template <
typename T,
typename VectorType>
746 using Tvmult_t =
decltype(std::declval<const T>().Tvmult(
747 std::declval<VectorType &>(),
748 std::declval<const VectorType &>()));
750 template <
typename T,
typename VectorType>
751 constexpr bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
753 template <
typename T,
typename VectorType>
754 using step_t =
decltype(std::declval<const T>().step(
755 std::declval<VectorType &>(),
756 std::declval<const VectorType &>()));
758 template <
typename T,
typename VectorType>
759 constexpr bool has_step = is_supported_operation<step_t, T, VectorType>;
761 template <
typename T,
typename VectorType>
763 decltype(std::declval<const T>().step(std::declval<VectorType &>(),
764 std::declval<const VectorType &>(),
765 std::declval<const double>()));
767 template <
typename T,
typename VectorType>
768 constexpr bool has_step_omega =
769 is_supported_operation<step_omega_t, T, VectorType>;
771 template <
typename T,
typename VectorType>
772 using Tstep_t =
decltype(std::declval<const T>().Tstep(
773 std::declval<VectorType &>(),
774 std::declval<const VectorType &>()));
776 template <
typename T,
typename VectorType>
777 constexpr bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
779 template <
typename T,
typename VectorType>
780 using Tstep_omega_t =
781 decltype(std::declval<const T>().Tstep(std::declval<VectorType &>(),
782 std::declval<const VectorType &>(),
783 std::declval<const double>()));
785 template <
typename T,
typename VectorType>
786 constexpr bool has_Tstep_omega =
787 is_supported_operation<Tstep_omega_t, T, VectorType>;
789 template <
typename T,
typename VectorType>
790 using jacobi_step_t =
decltype(std::declval<const T>().Jacobi_step(
791 std::declval<VectorType &>(),
792 std::declval<const VectorType &>(),
793 std::declval<const double>()));
795 template <
typename T,
typename VectorType>
796 constexpr bool has_jacobi_step =
797 is_supported_operation<jacobi_step_t, T, VectorType>;
799 template <
typename T,
typename VectorType>
800 using SOR_step_t =
decltype(std::declval<const T>().SOR_step(
801 std::declval<VectorType &>(),
802 std::declval<const VectorType &>(),
803 std::declval<const double>()));
805 template <
typename T,
typename VectorType>
806 constexpr bool has_SOR_step =
807 is_supported_operation<SOR_step_t, T, VectorType>;
809 template <
typename T,
typename VectorType>
810 using SSOR_step_t =
decltype(std::declval<const T>().SSOR_step(
811 std::declval<VectorType &>(),
812 std::declval<const VectorType &>(),
813 std::declval<const double>()));
815 template <
typename T,
typename VectorType>
816 constexpr bool has_SSOR_step =
817 is_supported_operation<SSOR_step_t, T, VectorType>;
819 template <
typename MatrixType>
820 class PreconditionJacobiImpl
823 PreconditionJacobiImpl(
const MatrixType &A,
const double relaxation)
825 , relaxation(relaxation)
828 template <
typename VectorType>
830 vmult(VectorType &dst,
const VectorType &src)
const
832 this->A->precondition_Jacobi(dst, src, this->relaxation);
835 template <
typename VectorType>
837 Tvmult(VectorType &dst,
const VectorType &src)
const
840 this->vmult(dst, src);
843 template <
typename VectorType,
844 std::enable_if_t<has_jacobi_step<MatrixType, VectorType>,
845 MatrixType> * =
nullptr>
847 step(VectorType &dst,
const VectorType &src)
const
849 this->A->Jacobi_step(dst, src, this->relaxation);
852 template <
typename VectorType,
853 std::enable_if_t<!has_jacobi_step<MatrixType, VectorType>,
854 MatrixType> * =
nullptr>
856 step(VectorType &,
const VectorType &)
const
860 "Matrix A does not provide a Jacobi_step() function!"));
863 template <
typename VectorType>
865 Tstep(VectorType &dst,
const VectorType &src)
const
868 this->step(dst, src);
873 const double relaxation;
876 template <
typename MatrixType>
877 class PreconditionSORImpl
880 PreconditionSORImpl(
const MatrixType &A,
const double relaxation)
882 , relaxation(relaxation)
885 template <
typename VectorType>
887 vmult(VectorType &dst,
const VectorType &src)
const
889 this->A->precondition_SOR(dst, src, this->relaxation);
892 template <
typename VectorType>
894 Tvmult(VectorType &dst,
const VectorType &src)
const
896 this->A->precondition_TSOR(dst, src, this->relaxation);
899 template <
typename VectorType,
900 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
901 MatrixType> * =
nullptr>
903 step(VectorType &dst,
const VectorType &src)
const
905 this->A->SOR_step(dst, src, this->relaxation);
908 template <
typename VectorType,
909 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
910 MatrixType> * =
nullptr>
912 step(VectorType &,
const VectorType &)
const
916 "Matrix A does not provide a SOR_step() function!"));
919 template <
typename VectorType,
920 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
921 MatrixType> * =
nullptr>
923 Tstep(VectorType &dst,
const VectorType &src)
const
925 this->A->TSOR_step(dst, src, this->relaxation);
928 template <
typename VectorType,
929 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
930 MatrixType> * =
nullptr>
932 Tstep(VectorType &,
const VectorType &)
const
936 "Matrix A does not provide a TSOR_step() function!"));
941 const double relaxation;
944 template <
typename MatrixType>
945 class PreconditionSSORImpl
948 using size_type =
typename MatrixType::size_type;
950 PreconditionSSORImpl(
const MatrixType &A,
const double relaxation)
952 , relaxation(relaxation)
963 const size_type n = this->A->n();
964 pos_right_of_diagonal.resize(n,
static_cast<std::size_t
>(-1));
965 for (size_type row = 0; row < n; ++row)
972 typename MatrixType::value_type>::const_iterator it =
974 for (; it < mat->
end(row); ++it)
975 if (it->column() > row)
977 pos_right_of_diagonal[row] = it - mat->
begin();
982 template <
typename VectorType>
984 vmult(VectorType &dst,
const VectorType &src)
const
986 this->A->precondition_SSOR(dst,
989 pos_right_of_diagonal);
992 template <
typename VectorType>
994 Tvmult(VectorType &dst,
const VectorType &src)
const
996 this->A->precondition_SSOR(dst,
999 pos_right_of_diagonal);
1002 template <
typename VectorType,
1003 std::enable_if_t<has_SSOR_step<MatrixType, VectorType>,
1004 MatrixType> * =
nullptr>
1006 step(VectorType &dst,
const VectorType &src)
const
1008 this->A->SSOR_step(dst, src, this->relaxation);
1011 template <
typename VectorType,
1012 std::enable_if_t<!has_SSOR_step<MatrixType, VectorType>,
1013 MatrixType> * =
nullptr>
1015 step(VectorType &,
const VectorType &)
const
1019 "Matrix A does not provide a SSOR_step() function!"));
1022 template <
typename VectorType>
1024 Tstep(VectorType &dst,
const VectorType &src)
const
1027 this->step(dst, src);
1032 const double relaxation;
1038 std::vector<std::size_t> pos_right_of_diagonal;
1041 template <
typename MatrixType>
1042 class PreconditionPSORImpl
1045 using size_type =
typename MatrixType::size_type;
1047 PreconditionPSORImpl(
const MatrixType &A,
1048 const double relaxation,
1049 const std::vector<size_type> &permutation,
1050 const std::vector<size_type> &inverse_permutation)
1052 , relaxation(relaxation)
1053 , permutation(permutation)
1054 , inverse_permutation(inverse_permutation)
1057 template <
typename VectorType>
1059 vmult(VectorType &dst,
const VectorType &src)
const
1062 this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
1065 template <
typename VectorType>
1067 Tvmult(VectorType &dst,
const VectorType &src)
const
1070 this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
1075 const double relaxation;
1077 const std::vector<size_type> &permutation;
1078 const std::vector<size_type> &inverse_permutation;
1081 template <
typename MatrixType,
1082 typename PreconditionerType,
1083 typename VectorType,
1084 std::enable_if_t<has_step_omega<PreconditionerType, VectorType>,
1085 PreconditionerType> * =
nullptr>
1087 step(
const MatrixType &,
1088 const PreconditionerType &preconditioner,
1090 const VectorType &src,
1091 const double relaxation,
1095 preconditioner.step(dst, src, relaxation);
1099 typename MatrixType,
1100 typename PreconditionerType,
1101 typename VectorType,
1102 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1103 has_step<PreconditionerType, VectorType>,
1104 PreconditionerType> * =
nullptr>
1106 step(
const MatrixType &,
1107 const PreconditionerType &preconditioner,
1109 const VectorType &src,
1110 const double relaxation,
1118 preconditioner.step(dst, src);
1122 typename MatrixType,
1123 typename PreconditionerType,
1124 typename VectorType,
1125 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1126 !has_step<PreconditionerType, VectorType>,
1127 PreconditionerType> * =
nullptr>
1129 step(
const MatrixType &A,
1130 const PreconditionerType &preconditioner,
1132 const VectorType &src,
1133 const double relaxation,
1134 VectorType &residual,
1137 residual.reinit(dst,
true);
1138 tmp.reinit(dst,
true);
1140 A.vmult(residual, dst);
1141 residual.sadd(-1.0, 1.0, src);
1143 preconditioner.vmult(tmp, residual);
1144 dst.add(relaxation, tmp);
1147 template <
typename MatrixType,
1148 typename PreconditionerType,
1149 typename VectorType,
1150 std::enable_if_t<has_Tstep_omega<PreconditionerType, VectorType>,
1151 PreconditionerType> * =
nullptr>
1153 Tstep(
const MatrixType &,
1154 const PreconditionerType &preconditioner,
1156 const VectorType &src,
1157 const double relaxation,
1161 preconditioner.Tstep(dst, src, relaxation);
1165 typename MatrixType,
1166 typename PreconditionerType,
1167 typename VectorType,
1168 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1169 has_Tstep<PreconditionerType, VectorType>,
1170 PreconditionerType> * =
nullptr>
1172 Tstep(
const MatrixType &,
1173 const PreconditionerType &preconditioner,
1175 const VectorType &src,
1176 const double relaxation,
1184 preconditioner.Tstep(dst, src);
1187 template <
typename MatrixType,
1188 typename VectorType,
1189 std::enable_if_t<has_Tvmult<MatrixType, VectorType>, MatrixType>
1192 Tvmult(
const MatrixType &A, VectorType &dst,
const VectorType &src)
1197 template <
typename MatrixType,
1198 typename VectorType,
1199 std::enable_if_t<!has_Tvmult<MatrixType, VectorType>, MatrixType>
1202 Tvmult(
const MatrixType &, VectorType &,
const VectorType &)
1205 ExcMessage(
"Matrix A does not provide a Tvmult() function!"));
1209 typename MatrixType,
1210 typename PreconditionerType,
1211 typename VectorType,
1212 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1213 !has_Tstep<PreconditionerType, VectorType>,
1214 PreconditionerType> * =
nullptr>
1216 Tstep(
const MatrixType &A,
1217 const PreconditionerType &preconditioner,
1219 const VectorType &src,
1220 const double relaxation,
1221 VectorType &residual,
1224 residual.reinit(dst,
true);
1225 tmp.reinit(dst,
true);
1227 Tvmult(A, residual, dst);
1228 residual.sadd(-1.0, 1.0, src);
1230 Tvmult(preconditioner, tmp, residual);
1231 dst.add(relaxation, tmp);
1235 template <
typename MatrixType,
1236 typename PreconditionerType,
1237 typename VectorType,
1238 std::enable_if_t<!has_vmult_with_std_functions_for_precondition<
1243 step_operations(
const MatrixType &A,
1244 const PreconditionerType &preconditioner,
1246 const VectorType &src,
1247 const double relaxation,
1250 const unsigned int i,
1251 const bool transposed)
1256 Tvmult(preconditioner, dst, src);
1258 preconditioner.vmult(dst, src);
1260 if (relaxation != 1.0)
1266 Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1268 step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1275 typename MatrixType,
1276 typename PreconditionerType,
1277 typename VectorType,
1279 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1281 !has_vmult_with_std_functions_for_precondition<MatrixType,
1285 step_operations(
const MatrixType &A,
1286 const PreconditionerType &preconditioner,
1288 const VectorType &src,
1289 const double relaxation,
1292 const unsigned int i,
1293 const bool transposed)
1296 using Number =
typename VectorType::value_type;
1297 Number *dst_ptr = dst.begin();
1298 const Number *src_ptr = src.begin();
1302 preconditioner.vmult(
1305 [&](
const unsigned int start_range,
const unsigned int end_range) {
1307 if (end_range > start_range)
1308 std::memset(dst_ptr + start_range,
1310 sizeof(Number) * (end_range - start_range));
1312 [&](
const unsigned int start_range,
const unsigned int end_range) {
1313 if (relaxation == 1.0)
1317 for (std::size_t i = start_range; i < end_range; ++i)
1318 dst_ptr[i] *= relaxation;
1323 tmp.reinit(src,
true);
1329 preconditioner.vmult(
1332 [&](
const unsigned int start_range,
const unsigned int end_range) {
1333 const auto tmp_ptr = tmp.begin();
1335 if (relaxation == 1.0)
1338 for (std::size_t i = start_range; i < end_range; ++i)
1339 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1347 for (std::size_t i = start_range; i < end_range; ++i)
1348 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1351 [&](
const unsigned int,
const unsigned int) {
1361 typename MatrixType,
1362 typename PreconditionerType,
1363 typename VectorType,
1365 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1367 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1370 step_operations(
const MatrixType &A,
1371 const PreconditionerType &preconditioner,
1373 const VectorType &src,
1374 const double relaxation,
1377 const unsigned int i,
1378 const bool transposed)
1381 using Number =
typename VectorType::value_type;
1383 Number *dst_ptr = dst.begin();
1384 const Number *src_ptr = src.begin();
1388 preconditioner.vmult(
1391 [&](
const unsigned int start_range,
const unsigned int end_range) {
1393 if (end_range > start_range)
1394 std::memset(dst_ptr + start_range,
1396 sizeof(Number) * (end_range - start_range));
1398 [&](
const unsigned int start_range,
const unsigned int end_range) {
1399 if (relaxation == 1.0)
1403 for (std::size_t i = start_range; i < end_range; ++i)
1404 dst_ptr[i] *= relaxation;
1409 tmp.reinit(src,
true);
1410 const auto tmp_ptr = tmp.begin();
1417 [&](
const unsigned int start_range,
const unsigned int end_range) {
1420 if (end_range > start_range)
1421 std::memset(tmp_ptr + start_range,
1423 sizeof(Number) * (end_range - start_range));
1425 [&](
const unsigned int start_range,
const unsigned int end_range) {
1426 if (relaxation == 1.0)
1429 for (std::size_t i = start_range; i < end_range; ++i)
1430 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1438 for (std::size_t i = start_range; i < end_range; ++i)
1439 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1443 preconditioner.vmult(dst, tmp, [](
const auto,
const auto) {
1452 typename MatrixType,
1453 typename VectorType,
1460 !has_vmult_with_std_functions<MatrixType,
1463 VectorType> * =
nullptr>
1465 step_operations(
const MatrixType &A,
1466 const ::DiagonalMatrix<VectorType> &preconditioner,
1468 const VectorType &src,
1469 const double relaxation,
1472 const unsigned int i,
1473 const bool transposed)
1475 using Number =
typename VectorType::value_type;
1479 Number *dst_ptr = dst.begin();
1480 const Number *src_ptr = src.begin();
1481 const Number *diag_ptr = preconditioner.get_vector().begin();
1483 if (relaxation == 1.0)
1486 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1487 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1492 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1493 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1498 tmp.reinit(src,
true);
1501 Tvmult(A, tmp, dst);
1505 Number *dst_ptr = dst.begin();
1506 const Number *src_ptr = src.begin();
1507 const Number *tmp_ptr = tmp.begin();
1508 const Number *diag_ptr = preconditioner.get_vector().begin();
1510 if (relaxation == 1.0)
1513 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1514 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1519 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1521 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1528 template <
typename MatrixType,
1529 typename VectorType,
1530 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1531 has_vmult_with_std_functions<
1535 VectorType> * =
nullptr>
1537 step_operations(
const MatrixType &A,
1538 const ::DiagonalMatrix<VectorType> &preconditioner,
1540 const VectorType &src,
1541 const double relaxation,
1544 const unsigned int i,
1545 const bool transposed)
1548 using Number =
typename VectorType::value_type;
1552 Number *dst_ptr = dst.begin();
1553 const Number *src_ptr = src.begin();
1554 const Number *diag_ptr = preconditioner.get_vector().begin();
1556 if (relaxation == 1.0)
1559 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1560 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1565 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1566 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1571 tmp.reinit(src,
true);
1578 [&](
const unsigned int start_range,
const unsigned int end_range) {
1580 if (end_range > start_range)
1581 std::memset(tmp.begin() + start_range,
1583 sizeof(Number) * (end_range - start_range));
1585 [&](
const unsigned int begin,
const unsigned int end) {
1586 const Number *dst_ptr = dst.begin();
1587 const Number *src_ptr = src.begin();
1588 Number *tmp_ptr = tmp.begin();
1589 const Number *diag_ptr = preconditioner.get_vector().begin();
1593 if (relaxation == 1.0)
1596 for (std::size_t i = begin; i < end; ++i)
1598 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1603 for (std::size_t i = begin; i < end; ++i)
1604 tmp_ptr[i] = dst_ptr[i] + relaxation *
1605 (src_ptr[i] - tmp_ptr[i]) *
1647template <
typename MatrixType = SparseMatrix<
double>>
1651 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1654 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1717template <
typename MatrixType = SparseMatrix<
double>>
1721 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1724 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1769template <
typename MatrixType = SparseMatrix<
double>>
1773 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1776 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1825template <
typename MatrixType = SparseMatrix<
double>>
1829 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1832 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1889 const std::vector<size_type> &permutation,
1890 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,
typename MemorySpace>
2345 for (
unsigned int block = 0; block < vector.
n_blocks(); ++block)
2346 set_initial_guess(vector.
block(block));
2349 template <
typename Number,
typename MemorySpace>
2366 Kokkos::RangePolicy<
typename MemorySpace::kokkos_space::execution_space,
2367 Kokkos::IndexType<types::global_dof_index>>
2368 policy(0, n_local_elements);
2369 Kokkos::parallel_for(
2370 "::PreconditionChebyshev::set_initial_guess",
2373 values_ptr[i] = (i + first_local_range) % 11;
2375 const Number mean_value = vector.
mean_value();
2376 vector.
add(-mean_value);
2379 struct EigenvalueTracker
2388 std::vector<double>
values;
2395 typename PreconditionerType>
2398 VectorType &eigenvector,
2399 const PreconditionerType &preconditioner,
2400 const unsigned int n_iterations)
2402 typename VectorType::value_type eigenvalue_estimate = 0.;
2403 eigenvector /= eigenvector.l2_norm();
2405 vector1.reinit(eigenvector,
true);
2406 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2407 vector2.reinit(eigenvector,
true);
2408 for (
unsigned int i = 0; i < n_iterations; ++i)
2410 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2412 matrix.vmult(vector2, eigenvector);
2413 preconditioner.vmult(vector1, vector2);
2416 matrix.vmult(vector1, eigenvector);
2418 eigenvalue_estimate = eigenvector * vector1;
2420 vector1 /= vector1.l2_norm();
2421 eigenvector.swap(vector1);
2423 return std::abs(eigenvalue_estimate);
2430 typename PreconditionerType>
2431 EigenvalueInformation
2432 estimate_eigenvalues(
2433 const EigenvalueAlgorithmAdditionalData<PreconditionerType> &
data,
2434 const MatrixType *matrix_ptr,
2435 VectorType &solution_old,
2436 VectorType &temp_vector1,
2437 const unsigned int degree)
2441 EigenvalueInformation info{};
2443 if (
data.eig_cg_n_iterations > 0)
2447 "Need to set at least two iterations to find eigenvalues."));
2449 internal::EigenvalueTracker eigenvalue_tracker;
2454 internal::set_initial_guess(temp_vector1);
2455 data.constraints.set_zero(temp_vector1);
2466 solver.connect_eigenvalues_slot(
2467 [&eigenvalue_tracker](
const std::vector<double> &
eigenvalues) {
2471 solver.solve(*matrix_ptr,
2474 *
data.preconditioner);
2476 info.cg_iterations = control.last_step();
2478 else if (
data.eigenvalue_algorithm ==
2484 ExcMessage(
"Cannot estimate the minimal eigenvalue with the "
2485 "power iteration"));
2487 eigenvalue_tracker.values.push_back(
2490 *
data.preconditioner,
2491 data.eig_cg_n_iterations));
2497 if (eigenvalue_tracker.values.empty())
2498 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2501 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2505 info.max_eigenvalue_estimate =
2506 1.2 * eigenvalue_tracker.values.back();
2511 info.max_eigenvalue_estimate =
data.max_eigenvalue;
2512 info.min_eigenvalue_estimate =
2513 data.max_eigenvalue /
data.smoothing_range;
2527template <
typename MatrixType>
2537template <
typename VectorType>
2546template <
typename VectorType>
2553template <
typename VectorType>
2562template <
typename VectorType>
2594 const double relaxation)
2595 : relaxation(relaxation)
2604 AdditionalData add_data;
2605 relaxation = add_data.relaxation;
2619template <
typename MatrixType>
2622 const MatrixType &matrix,
2632template <
typename VectorType>
2637 std::is_same_v<size_type, typename VectorType::size_type>,
2638 "PreconditionRichardson and VectorType must have the same size_type.");
2645template <
typename VectorType>
2650 std::is_same_v<size_type, typename VectorType::size_type>,
2651 "PreconditionRichardson and VectorType must have the same size_type.");
2656template <
typename VectorType>
2661 std::is_same_v<size_type, typename VectorType::size_type>,
2662 "PreconditionRichardson and VectorType must have the same size_type.");
2669template <
typename VectorType>
2674 std::is_same_v<size_type, typename VectorType::size_type>,
2675 "PreconditionRichardson and VectorType must have the same size_type.");
2696template <
typename MatrixType,
typename PreconditionerType>
2699 const MatrixType &rA,
2700 const AdditionalData ¶meters)
2703 eigenvalues_are_initialized =
false;
2707 this->data = parameters;
2711template <
typename MatrixType,
typename PreconditionerType>
2715 eigenvalues_are_initialized =
false;
2717 data.relaxation = 1.0;
2718 data.preconditioner =
nullptr;
2721template <
typename MatrixType,
typename PreconditionerType>
2730template <
typename MatrixType,
typename PreconditionerType>
2739template <
typename MatrixType,
typename PreconditionerType>
2740template <
typename VectorType>
2744 const VectorType &src)
const
2749 if (eigenvalues_are_initialized ==
false)
2750 estimate_eigenvalues(src);
2754 for (
unsigned int i = 0; i <
data.n_iterations; ++i)
2755 internal::PreconditionRelaxation::step_operations(*A,
2756 *
data.preconditioner,
2766template <
typename MatrixType,
typename PreconditionerType>
2767template <
typename VectorType>
2771 const VectorType &src)
const
2776 if (eigenvalues_are_initialized ==
false)
2777 estimate_eigenvalues(src);
2781 for (
unsigned int i = 0; i <
data.n_iterations; ++i)
2782 internal::PreconditionRelaxation::step_operations(
2783 *A, *
data.preconditioner, dst, src,
data.relaxation, tmp1, tmp2, i,
true);
2786template <
typename MatrixType,
typename PreconditionerType>
2787template <
typename VectorType>
2791 const VectorType &src)
const
2796 if (eigenvalues_are_initialized ==
false)
2797 estimate_eigenvalues(src);
2801 for (
unsigned int i = 1; i <=
data.n_iterations; ++i)
2802 internal::PreconditionRelaxation::step_operations(*A,
2803 *
data.preconditioner,
2813template <
typename MatrixType,
typename PreconditionerType>
2814template <
typename VectorType>
2818 const VectorType &src)
const
2823 if (eigenvalues_are_initialized ==
false)
2824 estimate_eigenvalues(src);
2828 for (
unsigned int i = 1; i <=
data.n_iterations; ++i)
2829 internal::PreconditionRelaxation::step_operations(
2830 *A, *
data.preconditioner, dst, src,
data.relaxation, tmp1, tmp2, i,
true);
2833template <
typename MatrixType,
typename PreconditionerType>
2834template <
typename VectorType>
2837 const VectorType &src)
const
2841 EigenvalueInformation info;
2843 if (
data.relaxation == 0.0)
2847 solution_old.reinit(src);
2848 temp_vector1.reinit(src,
true);
2850 info = internal::estimate_eigenvalues<MatrixType>(
2851 data, A, solution_old, temp_vector1,
data.n_iterations);
2853 const double alpha =
2854 (
data.smoothing_range > 1. ?
2855 info.max_eigenvalue_estimate /
data.smoothing_range :
2856 std::min(0.9 * info.max_eigenvalue_estimate,
2857 info.min_eigenvalue_estimate));
2860 ->
data.relaxation = 2.0 / (alpha + info.max_eigenvalue_estimate);
2864 ->eigenvalues_are_initialized =
true;
2869template <
typename MatrixType,
typename PreconditionerType>
2873 return data.relaxation;
2879template <
typename MatrixType>
2882 const AdditionalData ¶meters_in)
2886 parameters_in.relaxation != 0.0,
2888 "Relaxation cannot automatically be determined by PreconditionJacobi."));
2890 AdditionalData parameters;
2891 parameters.relaxation = 1.0;
2892 parameters.n_iterations = parameters_in.n_iterations;
2893 parameters.preconditioner =
2894 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2896 this->BaseClass::initialize(A, parameters);
2901template <
typename MatrixType>
2904 const AdditionalData ¶meters_in)
2908 parameters_in.relaxation != 0.0,
2910 "Relaxation cannot automatically be determined by PreconditionSOR."));
2912 AdditionalData parameters;
2913 parameters.relaxation = 1.0;
2914 parameters.n_iterations = parameters_in.n_iterations;
2915 parameters.preconditioner =
2916 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2918 this->BaseClass::initialize(A, parameters);
2923template <
typename MatrixType>
2926 const AdditionalData ¶meters_in)
2930 parameters_in.relaxation != 0.0,
2932 "Relaxation cannot automatically be determined by PreconditionSSOR."));
2934 AdditionalData parameters;
2935 parameters.relaxation = 1.0;
2936 parameters.n_iterations = parameters_in.n_iterations;
2937 parameters.preconditioner =
2938 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2940 this->BaseClass::initialize(A, parameters);
2947template <
typename MatrixType>
2950 const MatrixType &A,
2951 const std::vector<size_type> &p,
2952 const std::vector<size_type> &ip,
2953 const typename BaseClass::AdditionalData ¶meters_in)
2957 parameters_in.relaxation != 0.0,
2959 "Relaxation cannot automatically be determined by PreconditionPSOR."));
2961 typename BaseClass::AdditionalData parameters;
2962 parameters.relaxation = 1.0;
2963 parameters.n_iterations = parameters_in.n_iterations;
2964 parameters.preconditioner =
2965 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2967 this->BaseClass::initialize(A, parameters);
2971template <
typename MatrixType>
2974 const AdditionalData &additional_data)
2977 additional_data.permutation,
2978 additional_data.inverse_permutation,
2979 additional_data.parameters);
2982template <
typename MatrixType>
2984 const std::vector<size_type> &permutation,
2985 const std::vector<size_type> &inverse_permutation,
2988 : permutation(permutation)
2989 , inverse_permutation(inverse_permutation)
2990 , parameters(parameters)
2997template <
typename MatrixType,
typename VectorType>
2999 const MatrixType &M,
3000 const function_ptr method)
3002 , precondition(method)
3007template <
typename MatrixType,
typename VectorType>
3011 const VectorType &src)
const
3013 (
matrix.*precondition)(dst, src);
3021 template <
typename PreconditionerType>
3024 const double smoothing_range,
3025 const unsigned int eig_cg_n_iterations,
3026 const double eig_cg_residual,
3027 const double max_eigenvalue,
3029 : smoothing_range(smoothing_range)
3030 , eig_cg_n_iterations(eig_cg_n_iterations)
3031 , eig_cg_residual(eig_cg_residual)
3032 , max_eigenvalue(max_eigenvalue)
3033 , eigenvalue_algorithm(eigenvalue_algorithm)
3038 template <
typename PreconditionerType>
3039 inline EigenvalueAlgorithmAdditionalData<PreconditionerType> &
3040 EigenvalueAlgorithmAdditionalData<PreconditionerType>::operator=(
3041 const EigenvalueAlgorithmAdditionalData &other_data)
3043 smoothing_range = other_data.smoothing_range;
3044 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3045 eig_cg_residual = other_data.eig_cg_residual;
3046 max_eigenvalue = other_data.max_eigenvalue;
3047 preconditioner = other_data.preconditioner;
3048 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3049 constraints.copy_from(other_data.constraints);
3055template <
typename MatrixType,
typename PreconditionerType>
3058 const unsigned int n_iterations,
3059 const double smoothing_range,
3060 const unsigned int eig_cg_n_iterations,
3061 const double eig_cg_residual,
3062 const double max_eigenvalue,
3063 const EigenvalueAlgorithm eigenvalue_algorithm)
3064 :
internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3066 eig_cg_n_iterations,
3069 eigenvalue_algorithm)
3070 , relaxation(relaxation)
3071 , n_iterations(n_iterations)
3080 namespace PreconditionChebyshevImplementation
3088 template <
typename VectorType,
typename PreconditionerType>
3090 vector_updates(
const VectorType &rhs,
3091 const PreconditionerType &preconditioner,
3092 const unsigned int iteration_index,
3093 const double factor1,
3094 const double factor2,
3095 VectorType &solution_old,
3096 VectorType &temp_vector1,
3097 VectorType &temp_vector2,
3098 VectorType &solution)
3100 if (iteration_index == 0)
3102 solution.equ(factor2, rhs);
3103 preconditioner.vmult(solution_old, solution);
3105 else if (iteration_index == 1)
3108 temp_vector1.sadd(-1.0, 1.0, rhs);
3109 preconditioner.vmult(solution_old, temp_vector1);
3112 solution_old.sadd(factor2, 1 + factor1, solution);
3117 temp_vector1.sadd(-1.0, 1.0, rhs);
3118 preconditioner.vmult(temp_vector2, temp_vector1);
3121 solution_old.sadd(-factor1, factor2, temp_vector2);
3122 solution_old.add(1 + factor1, solution);
3125 solution.swap(solution_old);
3131 typename PreconditionerType,
3133 !has_vmult_with_std_functions_for_precondition<
3140 const PreconditionerType &preconditioner,
3141 const unsigned int iteration_index,
3142 const double factor1_,
3143 const double factor2_,
3152 const Number factor1 = factor1_;
3153 const Number factor1_plus_1 = 1. + factor1_;
3154 const Number factor2 = factor2_;
3156 if (iteration_index == 0)
3159 preconditioner.vmult(solution_old, rhs);
3162 const auto solution_old_ptr = solution_old.
begin();
3165 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
3167 else if (iteration_index == 1)
3170 temp_vector1.
sadd(-1.0, 1.0, rhs);
3172 preconditioner.vmult(solution_old, temp_vector1);
3175 const auto solution_ptr = solution.
begin();
3176 const auto solution_old_ptr = solution_old.
begin();
3179 solution_old_ptr[i] =
3180 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
3185 temp_vector1.
sadd(-1.0, 1.0, rhs);
3187 preconditioner.vmult(temp_vector2, temp_vector1);
3190 const auto solution_ptr = solution.
begin();
3191 const auto solution_old_ptr = solution_old.
begin();
3192 const auto temp_vector2_ptr = temp_vector2.
begin();
3195 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3196 factor1 * solution_old_ptr[i] +
3197 temp_vector2_ptr[i] * factor2;
3200 solution.
swap(solution_old);
3205 typename PreconditionerType,
3207 has_vmult_with_std_functions_for_precondition<
3214 const PreconditionerType &preconditioner,
3215 const unsigned int iteration_index,
3216 const double factor1_,
3217 const double factor2_,
3226 const Number factor1 = factor1_;
3227 const Number factor1_plus_1 = 1. + factor1_;
3228 const Number factor2 = factor2_;
3230 const auto rhs_ptr = rhs.
begin();
3231 const auto temp_vector1_ptr = temp_vector1.
begin();
3232 const auto temp_vector2_ptr = temp_vector2.
begin();
3233 const auto solution_ptr = solution.
begin();
3234 const auto solution_old_ptr = solution_old.
begin();
3236 if (iteration_index == 0)
3238 preconditioner.vmult(
3241 [&](
const auto start_range,
const auto end_range) {
3242 if (end_range > start_range)
3243 std::memset(solution_ptr + start_range,
3245 sizeof(Number) * (end_range - start_range));
3247 [&](
const auto begin,
const auto end) {
3249 for (std::size_t i = begin; i <
end; ++i)
3250 solution_ptr[i] *= factor2;
3255 preconditioner.vmult(
3258 [&](
const auto begin,
const auto end) {
3260 std::memset(temp_vector2_ptr + begin,
3262 sizeof(Number) * (end - begin));
3265 for (std::size_t i = begin; i <
end; ++i)
3266 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3268 [&](
const auto begin,
const auto end) {
3269 if (iteration_index == 1)
3272 for (std::size_t i = begin; i <
end; ++i)
3273 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
3274 factor2 * temp_vector2_ptr[i];
3279 for (std::size_t i = begin; i <
end; ++i)
3280 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3281 factor1 * solution_old_ptr[i] +
3282 factor2 * temp_vector2_ptr[i];
3287 if (iteration_index > 0)
3289 solution_old.
swap(temp_vector2);
3290 solution_old.
swap(solution);
3297 template <
typename Number>
3298 struct VectorUpdater
3300 VectorUpdater(
const Number *rhs,
3301 const Number *matrix_diagonal_inverse,
3302 const unsigned int iteration_index,
3303 const Number factor1,
3304 const Number factor2,
3305 Number *solution_old,
3309 , matrix_diagonal_inverse(matrix_diagonal_inverse)
3310 , iteration_index(iteration_index)
3313 , solution_old(solution_old)
3314 , tmp_vector(tmp_vector)
3315 , solution(solution)
3319 apply_to_subrange(
const std::size_t begin,
const std::size_t end)
const
3325 const Number factor1 = this->factor1;
3326 const Number factor1_plus_1 = 1. + this->factor1;
3327 const Number factor2 = this->factor2;
3328 if (iteration_index == 0)
3331 for (std::size_t i = begin; i <
end; ++i)
3332 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
3334 else if (iteration_index == 1)
3338 for (std::size_t i = begin; i <
end; ++i)
3342 factor1_plus_1 * solution[i] +
3343 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3350 for (std::size_t i = begin; i <
end; ++i)
3357 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
3358 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3363 const Number *matrix_diagonal_inverse;
3364 const unsigned int iteration_index;
3365 const Number factor1;
3366 const Number factor2;
3367 mutable Number *solution_old;
3368 mutable Number *tmp_vector;
3369 mutable Number *solution;
3372 template <
typename Number>
3375 VectorUpdatesRange(
const VectorUpdater<Number> &updater,
3376 const std::size_t
size)
3380 VectorUpdatesRange::apply_to_subrange(0,
size);
3388 ~VectorUpdatesRange()
override =
default;
3391 apply_to_subrange(
const std::size_t begin,
3392 const std::size_t end)
const override
3394 updater.apply_to_subrange(begin, end);
3397 const VectorUpdater<Number> &updater;
3401 template <
typename Number>
3404 const ::Vector<Number> &rhs,
3406 const unsigned int iteration_index,
3407 const double factor1,
3408 const double factor2,
3414 VectorUpdater<Number> upd(rhs.begin(),
3415 jacobi.get_vector().begin(),
3419 solution_old.
begin(),
3420 temp_vector1.
begin(),
3422 VectorUpdatesRange<Number>(upd, rhs.size());
3425 if (iteration_index == 0)
3432 solution.
swap(temp_vector1);
3433 solution_old.
swap(temp_vector1);
3438 template <
typename Number>
3442 const ::DiagonalMatrix<
3444 const unsigned int iteration_index,
3445 const double factor1,
3446 const double factor2,
3454 VectorUpdater<Number> upd(rhs.
begin(),
3455 jacobi.get_vector().begin(),
3459 solution_old.
begin(),
3460 temp_vector1.
begin(),
3465 if (iteration_index == 0)
3472 solution.
swap(temp_vector1);
3473 solution_old.
swap(temp_vector1);
3484 typename PreconditionerType,
3488 PreconditionerType> &&
3489 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3491 has_vmult_with_std_functions_for_precondition<
MatrixType,
3495 vmult_and_update(
const MatrixType &matrix,
3496 const PreconditionerType &preconditioner,
3497 const VectorType &rhs,
3498 const unsigned int iteration_index,
3499 const double factor1,
3500 const double factor2,
3501 VectorType &solution,
3502 VectorType &solution_old,
3503 VectorType &temp_vector1,
3504 VectorType &temp_vector2)
3506 if (iteration_index > 0)
3507 matrix.vmult(temp_vector1, solution);
3524 typename PreconditionerType,
3528 PreconditionerType> &&
3529 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3531 has_vmult_with_std_functions_for_precondition<
MatrixType,
3535 vmult_and_update(
const MatrixType &matrix,
3536 const PreconditionerType &preconditioner,
3537 const VectorType &rhs,
3538 const unsigned int iteration_index,
3539 const double factor1_,
3540 const double factor2_,
3541 VectorType &solution,
3542 VectorType &solution_old,
3543 VectorType &temp_vector1,
3544 VectorType &temp_vector2)
3546 using Number =
typename VectorType::value_type;
3548 const Number factor1 = factor1_;
3549 const Number factor1_plus_1 = 1. + factor1_;
3550 const Number factor2 = factor2_;
3552 if (iteration_index == 0)
3554 preconditioner.vmult(
3557 [&](
const unsigned int start_range,
const unsigned int end_range) {
3559 if (end_range > start_range)
3560 std::memset(solution.begin() + start_range,
3562 sizeof(Number) * (end_range - start_range));
3564 [&](
const unsigned int start_range,
const unsigned int end_range) {
3565 const auto solution_ptr = solution.begin();
3568 for (std::size_t i = start_range; i < end_range; ++i)
3569 solution_ptr[i] *= factor2;
3574 temp_vector1.reinit(rhs,
true);
3575 temp_vector2.reinit(rhs,
true);
3581 [&](
const unsigned int start_range,
const unsigned int end_range) {
3584 if (end_range > start_range)
3585 std::memset(temp_vector1.begin() + start_range,
3587 sizeof(Number) * (end_range - start_range));
3589 [&](
const unsigned int start_range,
const unsigned int end_range) {
3590 const auto rhs_ptr = rhs.begin();
3591 const auto tmp_ptr = temp_vector1.begin();
3594 for (std::size_t i = start_range; i < end_range; ++i)
3595 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3599 preconditioner.vmult(
3602 [&](
const unsigned int start_range,
const unsigned int end_range) {
3605 if (end_range > start_range)
3606 std::memset(temp_vector2.begin() + start_range,
3608 sizeof(Number) * (end_range - start_range));
3610 [&](
const unsigned int start_range,
const unsigned int end_range) {
3611 const auto solution_ptr = solution.begin();
3612 const auto solution_old_ptr = solution_old.begin();
3613 const auto tmp_ptr = temp_vector2.begin();
3615 if (iteration_index == 1)
3618 for (std::size_t i = start_range; i < end_range; ++i)
3620 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3625 for (std::size_t i = start_range; i < end_range; ++i)
3626 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3627 factor1 * solution_old_ptr[i] +
3628 factor2 * tmp_ptr[i];
3632 solution.swap(temp_vector2);
3633 solution_old.swap(temp_vector2);
3641 typename PreconditionerType,
3642 std::enable_if_t<has_vmult_with_std_functions<
MatrixType,
3644 PreconditionerType>,
3647 vmult_and_update(
const MatrixType &matrix,
3648 const PreconditionerType &preconditioner,
3649 const VectorType &rhs,
3650 const unsigned int iteration_index,
3651 const double factor1,
3652 const double factor2,
3653 VectorType &solution,
3654 VectorType &solution_old,
3655 VectorType &temp_vector1,
3658 using Number =
typename VectorType::value_type;
3659 VectorUpdater<Number> updater(rhs.begin(),
3660 preconditioner.get_vector().begin(),
3664 solution_old.begin(),
3665 temp_vector1.begin(),
3667 if (iteration_index > 0)
3671 [&](
const unsigned int start_range,
const unsigned int end_range) {
3674 if (end_range > start_range)
3675 std::memset(temp_vector1.begin() + start_range,
3677 sizeof(Number) * (end_range - start_range));
3679 [&](
const unsigned int start_range,
const unsigned int end_range) {
3680 if (end_range > start_range)
3681 updater.apply_to_subrange(start_range, end_range);
3684 updater.apply_to_subrange(0U, solution.locally_owned_size());
3687 if (iteration_index == 0)
3694 solution.swap(temp_vector1);
3695 solution_old.swap(temp_vector1);
3699 template <
typename MatrixType,
typename PreconditionerType>
3701 initialize_preconditioner(
3702 const MatrixType &matrix,
3703 std::shared_ptr<PreconditionerType> &preconditioner)
3706 (void)preconditioner;
3710 template <
typename MatrixType,
typename VectorType>
3712 initialize_preconditioner(
3713 const MatrixType &matrix,
3716 if (preconditioner.get() ==
nullptr || preconditioner->m() !=
matrix.m())
3718 if (preconditioner.get() ==
nullptr)
3720 std::make_shared<::DiagonalMatrix<VectorType>>();
3723 preconditioner->m() == 0,
3725 "Preconditioner appears to be initialized but not sized correctly"));
3728 if (preconditioner->m() !=
matrix.m())
3730 preconditioner->get_vector().reinit(
matrix.m());
3731 for (
typename VectorType::size_type i = 0; i <
matrix.m(); ++i)
3732 preconditioner->get_vector()(i) = 1. /
matrix.el(i, i);
3741template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3744 const double smoothing_range,
3745 const unsigned int eig_cg_n_iterations,
3746 const double eig_cg_residual,
3747 const double max_eigenvalue,
3748 const EigenvalueAlgorithm eigenvalue_algorithm,
3749 const PolynomialType polynomial_type)
3750 :
internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3752 eig_cg_n_iterations,
3755 eigenvalue_algorithm)
3757 , polynomial_type(polynomial_type)
3762template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3767 , eigenvalues_are_initialized(false)
3770 std::is_same_v<size_type, typename VectorType::size_type>,
3771 "PreconditionChebyshev and VectorType must have the same size_type.");
3776template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3779 const MatrixType &matrix,
3780 const AdditionalData &additional_data)
3783 data = additional_data;
3785 ExcMessage(
"The degree of the Chebyshev method must be positive."));
3786 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3787 matrix,
data.preconditioner);
3788 eigenvalues_are_initialized =
false;
3793template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3797 eigenvalues_are_initialized =
false;
3798 theta = delta = 1.0;
3799 matrix_ptr =
nullptr;
3802 solution_old.reinit(empty_vector);
3803 temp_vector1.reinit(empty_vector);
3804 temp_vector2.reinit(empty_vector);
3806 data.preconditioner.reset();
3811template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3818 solution_old.reinit(src);
3819 temp_vector1.reinit(src,
true);
3821 auto info = internal::estimate_eigenvalues<MatrixType>(
3822 data, matrix_ptr, solution_old, temp_vector1,
data.degree);
3824 const double alpha = (
data.smoothing_range > 1. ?
3825 info.max_eigenvalue_estimate /
data.smoothing_range :
3826 std::min(0.9 * info.max_eigenvalue_estimate,
3827 info.min_eigenvalue_estimate));
3836 const double actual_range = info.max_eigenvalue_estimate / alpha;
3837 const double sigma = (1. -
std::sqrt(1. / actual_range)) /
3839 const double eps =
data.smoothing_range;
3844 1 +
static_cast<unsigned int>(
3849 info.degree =
data.degree;
3854 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3855 (info.max_eigenvalue_estimate) :
3856 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3859 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3863 using NumberType =
typename VectorType::value_type;
3871 (std::is_same_v<VectorType,
3874 temp_vector2.
reinit(src, true);
3878 temp_vector2.reinit(empty_vector);
3883 ->eigenvalues_are_initialized =
true;
3890template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3893 VectorType &solution,
3894 const VectorType &rhs)
const
3896 std::lock_guard<std::mutex> lock(mutex);
3897 if (eigenvalues_are_initialized ==
false)
3898 estimate_eigenvalues(rhs);
3900 internal::PreconditionChebyshevImplementation::vmult_and_update(
3902 *
data.preconditioner,
3906 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3907 (4. / (3. * delta)) :
3919 double rhok = delta / theta, sigma = theta / delta;
3920 for (
unsigned int k = 0; k <
data.degree - 1; ++k)
3922 double factor1 = 0.0;
3923 double factor2 = 0.0;
3925 if (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3927 factor1 = (2 * k + 1.) / (2 * k + 5.);
3928 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3932 const double rhokp = 1. / (2. * sigma - rhok);
3933 factor1 = rhokp * rhok;
3934 factor2 = 2. * rhokp / delta;
3938 internal::PreconditionChebyshevImplementation::vmult_and_update(
3940 *
data.preconditioner,
3954template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3957 VectorType &solution,
3958 const VectorType &rhs)
const
3960 std::lock_guard<std::mutex> lock(mutex);
3961 if (eigenvalues_are_initialized ==
false)
3962 estimate_eigenvalues(rhs);
3964 internal::PreconditionChebyshevImplementation::vector_updates(
3966 *
data.preconditioner,
3969 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3970 (4. / (3. * delta)) :
3980 double rhok = delta / theta, sigma = theta / delta;
3981 for (
unsigned int k = 0; k <
data.degree - 1; ++k)
3983 double factor1 = 0.0;
3984 double factor2 = 0.0;
3986 if (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3988 factor1 = (2 * k + 1.) / (2 * k + 5.);
3989 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3993 const double rhokp = 1. / (2. * sigma - rhok);
3994 factor1 = rhokp * rhok;
3995 factor2 = 2. * rhokp / delta;
3999 matrix_ptr->Tvmult(temp_vector1, solution);
4000 internal::PreconditionChebyshevImplementation::vector_updates(
4002 *
data.preconditioner,
4015template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4018 VectorType &solution,
4019 const VectorType &rhs)
const
4021 std::lock_guard<std::mutex> lock(mutex);
4022 if (eigenvalues_are_initialized ==
false)
4023 estimate_eigenvalues(rhs);
4025 internal::PreconditionChebyshevImplementation::vmult_and_update(
4027 *
data.preconditioner,
4031 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4032 (4. / (3. * delta)) :
4042 double rhok = delta / theta, sigma = theta / delta;
4043 for (
unsigned int k = 0; k <
data.degree - 1; ++k)
4045 double factor1 = 0.0;
4046 double factor2 = 0.0;
4048 if (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4050 factor1 = (2 * k + 1.) / (2 * k + 5.);
4051 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4055 const double rhokp = 1. / (2. * sigma - rhok);
4056 factor1 = rhokp * rhok;
4057 factor2 = 2. * rhokp / delta;
4061 internal::PreconditionChebyshevImplementation::vmult_and_update(
4063 *
data.preconditioner,
4077template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4080 VectorType &solution,
4081 const VectorType &rhs)
const
4083 std::lock_guard<std::mutex> lock(mutex);
4084 if (eigenvalues_are_initialized ==
false)
4085 estimate_eigenvalues(rhs);
4087 matrix_ptr->Tvmult(temp_vector1, solution);
4088 internal::PreconditionChebyshevImplementation::vector_updates(
4090 *
data.preconditioner,
4093 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4094 (4. / (3. * delta)) :
4104 double rhok = delta / theta, sigma = theta / delta;
4105 for (
unsigned int k = 0; k <
data.degree - 1; ++k)
4107 double factor1 = 0.0;
4108 double factor2 = 0.0;
4110 if (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4112 factor1 = (2 * k + 1.) / (2 * k + 5.);
4113 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4117 const double rhokp = 1. / (2. * sigma - rhok);
4118 factor1 = rhokp * rhok;
4119 factor2 = 2. * rhokp / delta;
4123 matrix_ptr->Tvmult(temp_vector1, solution);
4124 internal::PreconditionChebyshevImplementation::vector_updates(
4126 *
data.preconditioner,
4139template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4142 PreconditionerType>::size_type
4146 return matrix_ptr->m();
4151template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4154 PreconditionerType>::size_type
4158 return matrix_ptr->n();
unsigned int n_blocks() const
BlockType & block(const unsigned int i)
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
size_type locally_owned_size() const
void swap(Vector< Number, MemorySpace > &v) noexcept
::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
void vmult(VectorType &dst, const VectorType &src) const
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
bool eigenvalues_are_initialized
ObserverPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
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
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
ObserverPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
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) noexcept
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
@ matrix
Contents is actually a matrix.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
unsigned int minimum_parallel_grain_size
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr unsigned int invalid_unsigned_int
::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)