15#ifndef dealii_precondition_h
16#define dealii_precondition_h
39template <
typename number>
41template <
typename number>
47 template <
typename,
typename>
126 template <
typename PreconditionerType>
247 template <
typename MatrixType>
255 template <
typename VectorType>
257 vmult(VectorType &,
const VectorType &)
const;
263 template <
typename VectorType>
265 Tvmult(VectorType &,
const VectorType &)
const;
270 template <
typename VectorType>
278 template <
typename VectorType>
376 template <
typename MatrixType>
383 template <
typename VectorType>
385 vmult(VectorType &,
const VectorType &)
const;
391 template <
typename VectorType>
393 Tvmult(VectorType &,
const VectorType &)
const;
397 template <
typename VectorType>
405 template <
typename VectorType>
495template <
typename MatrixType = SparseMatrix<
double>,
496 typename VectorType = Vector<
double>>
504 const VectorType &) const;
518 vmult(VectorType &dst,
const VectorType &src)
const;
563template <
typename MatrixType = SparseMatrix<
double>,
564 typename PreconditionerType = IdentityMatrix>
592 EigenvalueAlgorithm::lanczos);
638 template <
typename VectorType>
640 vmult(VectorType &,
const VectorType &)
const;
646 template <
typename VectorType>
648 Tvmult(VectorType &,
const VectorType &)
const;
653 template <
typename VectorType>
655 step(VectorType &x,
const VectorType &rhs)
const;
660 template <
typename VectorType>
662 Tstep(VectorType &x,
const VectorType &rhs)
const;
669 template <
typename VectorType>
714 template <
typename MatrixType,
typename VectorType>
715 using vmult_functions_t =
decltype(std::declval<const MatrixType>().vmult(
716 std::declval<VectorType &>(),
717 std::declval<const VectorType &>(),
719 const std::function<
void(
const unsigned int,
const unsigned int)> &>(),
721 const std::function<
void(
const unsigned int,
const unsigned int)> &>()));
723 template <
typename MatrixType,
725 typename PreconditionerType>
726 constexpr bool has_vmult_with_std_functions =
727 is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
728 std::is_same_v<PreconditionerType, ::DiagonalMatrix<VectorType>> &&
729 (std::is_same_v<VectorType,
737 template <
typename MatrixType,
typename VectorType>
738 constexpr bool has_vmult_with_std_functions_for_precondition =
739 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
743 template <
typename T,
typename VectorType>
744 using Tvmult_t =
decltype(std::declval<const T>().Tvmult(
745 std::declval<VectorType &>(),
746 std::declval<const VectorType &>()));
748 template <
typename T,
typename VectorType>
749 constexpr bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
751 template <
typename T,
typename VectorType>
752 using step_t =
decltype(std::declval<const T>().step(
753 std::declval<VectorType &>(),
754 std::declval<const VectorType &>()));
756 template <
typename T,
typename VectorType>
757 constexpr bool has_step = is_supported_operation<step_t, T, VectorType>;
759 template <
typename T,
typename VectorType>
761 decltype(std::declval<const T>().step(std::declval<VectorType &>(),
762 std::declval<const VectorType &>(),
763 std::declval<const double>()));
765 template <
typename T,
typename VectorType>
766 constexpr bool has_step_omega =
767 is_supported_operation<step_omega_t, T, VectorType>;
769 template <
typename T,
typename VectorType>
770 using Tstep_t =
decltype(std::declval<const T>().Tstep(
771 std::declval<VectorType &>(),
772 std::declval<const VectorType &>()));
774 template <
typename T,
typename VectorType>
775 constexpr bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
777 template <
typename T,
typename VectorType>
778 using Tstep_omega_t =
779 decltype(std::declval<const T>().Tstep(std::declval<VectorType &>(),
780 std::declval<const VectorType &>(),
781 std::declval<const double>()));
783 template <
typename T,
typename VectorType>
784 constexpr bool has_Tstep_omega =
785 is_supported_operation<Tstep_omega_t, T, VectorType>;
787 template <
typename T,
typename VectorType>
788 using jacobi_step_t =
decltype(std::declval<const T>().Jacobi_step(
789 std::declval<VectorType &>(),
790 std::declval<const VectorType &>(),
791 std::declval<const double>()));
793 template <
typename T,
typename VectorType>
794 constexpr bool has_jacobi_step =
795 is_supported_operation<jacobi_step_t, T, VectorType>;
797 template <
typename T,
typename VectorType>
798 using SOR_step_t =
decltype(std::declval<const T>().SOR_step(
799 std::declval<VectorType &>(),
800 std::declval<const VectorType &>(),
801 std::declval<const double>()));
803 template <
typename T,
typename VectorType>
804 constexpr bool has_SOR_step =
805 is_supported_operation<SOR_step_t, T, VectorType>;
807 template <
typename T,
typename VectorType>
808 using SSOR_step_t =
decltype(std::declval<const T>().SSOR_step(
809 std::declval<VectorType &>(),
810 std::declval<const VectorType &>(),
811 std::declval<const double>()));
813 template <
typename T,
typename VectorType>
814 constexpr bool has_SSOR_step =
815 is_supported_operation<SSOR_step_t, T, VectorType>;
817 template <
typename MatrixType>
818 class PreconditionJacobiImpl
821 PreconditionJacobiImpl(
const MatrixType &A,
const double relaxation)
823 , relaxation(relaxation)
826 template <
typename VectorType>
828 vmult(VectorType &dst,
const VectorType &src)
const
830 this->A->precondition_Jacobi(dst, src, this->relaxation);
833 template <
typename VectorType>
835 Tvmult(VectorType &dst,
const VectorType &src)
const
838 this->vmult(dst, src);
841 template <
typename VectorType,
842 std::enable_if_t<has_jacobi_step<MatrixType, VectorType>,
843 MatrixType> * =
nullptr>
845 step(VectorType &dst,
const VectorType &src)
const
847 this->A->Jacobi_step(dst, src, this->relaxation);
850 template <
typename VectorType,
851 std::enable_if_t<!has_jacobi_step<MatrixType, VectorType>,
852 MatrixType> * =
nullptr>
854 step(VectorType &,
const VectorType &)
const
858 "Matrix A does not provide a Jacobi_step() function!"));
861 template <
typename VectorType>
863 Tstep(VectorType &dst,
const VectorType &src)
const
866 this->step(dst, src);
871 const double relaxation;
874 template <
typename MatrixType>
875 class PreconditionSORImpl
878 PreconditionSORImpl(
const MatrixType &A,
const double relaxation)
880 , relaxation(relaxation)
883 template <
typename VectorType>
885 vmult(VectorType &dst,
const VectorType &src)
const
887 this->A->precondition_SOR(dst, src, this->relaxation);
890 template <
typename VectorType>
892 Tvmult(VectorType &dst,
const VectorType &src)
const
894 this->A->precondition_TSOR(dst, src, this->relaxation);
897 template <
typename VectorType,
898 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
899 MatrixType> * =
nullptr>
901 step(VectorType &dst,
const VectorType &src)
const
903 this->A->SOR_step(dst, src, this->relaxation);
906 template <
typename VectorType,
907 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
908 MatrixType> * =
nullptr>
910 step(VectorType &,
const VectorType &)
const
914 "Matrix A does not provide a SOR_step() function!"));
917 template <
typename VectorType,
918 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
919 MatrixType> * =
nullptr>
921 Tstep(VectorType &dst,
const VectorType &src)
const
923 this->A->TSOR_step(dst, src, this->relaxation);
926 template <
typename VectorType,
927 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
928 MatrixType> * =
nullptr>
930 Tstep(VectorType &,
const VectorType &)
const
934 "Matrix A does not provide a TSOR_step() function!"));
939 const double relaxation;
942 template <
typename MatrixType>
943 class PreconditionSSORImpl
946 using size_type =
typename MatrixType::size_type;
948 PreconditionSSORImpl(
const MatrixType &A,
const double relaxation)
950 , relaxation(relaxation)
961 const size_type n = this->A->n();
962 pos_right_of_diagonal.resize(n,
static_cast<std::size_t
>(-1));
963 for (size_type row = 0; row < n; ++row)
970 typename MatrixType::value_type>::const_iterator it =
972 for (; it < mat->
end(row); ++it)
973 if (it->column() > row)
975 pos_right_of_diagonal[row] = it - mat->
begin();
980 template <
typename VectorType>
982 vmult(VectorType &dst,
const VectorType &src)
const
984 this->A->precondition_SSOR(dst,
987 pos_right_of_diagonal);
990 template <
typename VectorType>
992 Tvmult(VectorType &dst,
const VectorType &src)
const
994 this->A->precondition_SSOR(dst,
997 pos_right_of_diagonal);
1000 template <
typename VectorType,
1001 std::enable_if_t<has_SSOR_step<MatrixType, VectorType>,
1002 MatrixType> * =
nullptr>
1004 step(VectorType &dst,
const VectorType &src)
const
1006 this->A->SSOR_step(dst, src, this->relaxation);
1009 template <
typename VectorType,
1010 std::enable_if_t<!has_SSOR_step<MatrixType, VectorType>,
1011 MatrixType> * =
nullptr>
1013 step(VectorType &,
const VectorType &)
const
1017 "Matrix A does not provide a SSOR_step() function!"));
1020 template <
typename VectorType>
1022 Tstep(VectorType &dst,
const VectorType &src)
const
1025 this->step(dst, src);
1030 const double relaxation;
1036 std::vector<std::size_t> pos_right_of_diagonal;
1039 template <
typename MatrixType>
1040 class PreconditionPSORImpl
1043 using size_type =
typename MatrixType::size_type;
1045 PreconditionPSORImpl(
const MatrixType &A,
1046 const double relaxation,
1047 const std::vector<size_type> &permutation,
1048 const std::vector<size_type> &inverse_permutation)
1050 , relaxation(relaxation)
1051 , permutation(permutation)
1052 , inverse_permutation(inverse_permutation)
1055 template <
typename VectorType>
1057 vmult(VectorType &dst,
const VectorType &src)
const
1060 this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
1063 template <
typename VectorType>
1065 Tvmult(VectorType &dst,
const VectorType &src)
const
1068 this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
1073 const double relaxation;
1075 const std::vector<size_type> &permutation;
1076 const std::vector<size_type> &inverse_permutation;
1079 template <
typename MatrixType,
1080 typename PreconditionerType,
1081 typename VectorType,
1082 std::enable_if_t<has_step_omega<PreconditionerType, VectorType>,
1083 PreconditionerType> * =
nullptr>
1085 step(
const MatrixType &,
1086 const PreconditionerType &preconditioner,
1088 const VectorType &src,
1089 const double relaxation,
1093 preconditioner.step(dst, src, relaxation);
1097 typename MatrixType,
1098 typename PreconditionerType,
1099 typename VectorType,
1100 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1101 has_step<PreconditionerType, VectorType>,
1102 PreconditionerType> * =
nullptr>
1104 step(
const MatrixType &,
1105 const PreconditionerType &preconditioner,
1107 const VectorType &src,
1108 const double relaxation,
1116 preconditioner.step(dst, src);
1120 typename MatrixType,
1121 typename PreconditionerType,
1122 typename VectorType,
1123 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
1124 !has_step<PreconditionerType, VectorType>,
1125 PreconditionerType> * =
nullptr>
1127 step(
const MatrixType &A,
1128 const PreconditionerType &preconditioner,
1130 const VectorType &src,
1131 const double relaxation,
1132 VectorType &residual,
1135 residual.reinit(dst,
true);
1136 tmp.reinit(dst,
true);
1138 A.vmult(residual, dst);
1139 residual.sadd(-1.0, 1.0, src);
1141 preconditioner.vmult(tmp, residual);
1142 dst.add(relaxation, tmp);
1145 template <
typename MatrixType,
1146 typename PreconditionerType,
1147 typename VectorType,
1148 std::enable_if_t<has_Tstep_omega<PreconditionerType, VectorType>,
1149 PreconditionerType> * =
nullptr>
1151 Tstep(
const MatrixType &,
1152 const PreconditionerType &preconditioner,
1154 const VectorType &src,
1155 const double relaxation,
1159 preconditioner.Tstep(dst, src, relaxation);
1163 typename MatrixType,
1164 typename PreconditionerType,
1165 typename VectorType,
1166 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1167 has_Tstep<PreconditionerType, VectorType>,
1168 PreconditionerType> * =
nullptr>
1170 Tstep(
const MatrixType &,
1171 const PreconditionerType &preconditioner,
1173 const VectorType &src,
1174 const double relaxation,
1182 preconditioner.Tstep(dst, src);
1185 template <
typename MatrixType,
1186 typename VectorType,
1187 std::enable_if_t<has_Tvmult<MatrixType, VectorType>, MatrixType>
1190 Tvmult(
const MatrixType &A, VectorType &dst,
const VectorType &src)
1195 template <
typename MatrixType,
1196 typename VectorType,
1197 std::enable_if_t<!has_Tvmult<MatrixType, VectorType>, MatrixType>
1200 Tvmult(
const MatrixType &, VectorType &,
const VectorType &)
1203 ExcMessage(
"Matrix A does not provide a Tvmult() function!"));
1207 typename MatrixType,
1208 typename PreconditionerType,
1209 typename VectorType,
1210 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1211 !has_Tstep<PreconditionerType, VectorType>,
1212 PreconditionerType> * =
nullptr>
1214 Tstep(
const MatrixType &A,
1215 const PreconditionerType &preconditioner,
1217 const VectorType &src,
1218 const double relaxation,
1219 VectorType &residual,
1222 residual.reinit(dst,
true);
1223 tmp.reinit(dst,
true);
1225 Tvmult(A, residual, dst);
1226 residual.sadd(-1.0, 1.0, src);
1228 Tvmult(preconditioner, tmp, residual);
1229 dst.add(relaxation, tmp);
1233 template <
typename MatrixType,
1234 typename PreconditionerType,
1235 typename VectorType,
1236 std::enable_if_t<!has_vmult_with_std_functions_for_precondition<
1241 step_operations(
const MatrixType &A,
1242 const PreconditionerType &preconditioner,
1244 const VectorType &src,
1245 const double relaxation,
1248 const unsigned int i,
1249 const bool transposed)
1254 Tvmult(preconditioner, dst, src);
1256 preconditioner.vmult(dst, src);
1258 if (relaxation != 1.0)
1264 Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1266 step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1273 typename MatrixType,
1274 typename PreconditionerType,
1275 typename VectorType,
1277 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1279 !has_vmult_with_std_functions_for_precondition<MatrixType,
1283 step_operations(
const MatrixType &A,
1284 const PreconditionerType &preconditioner,
1286 const VectorType &src,
1287 const double relaxation,
1290 const unsigned int i,
1291 const bool transposed)
1294 using Number =
typename VectorType::value_type;
1295 Number *dst_ptr = dst.begin();
1296 const Number *src_ptr = src.begin();
1300 preconditioner.vmult(
1303 [&](
const unsigned int start_range,
const unsigned int end_range) {
1305 if (end_range > start_range)
1306 std::memset(dst_ptr + start_range,
1308 sizeof(Number) * (end_range - start_range));
1310 [&](
const unsigned int start_range,
const unsigned int end_range) {
1311 if (relaxation == 1.0)
1315 for (std::size_t i = start_range; i < end_range; ++i)
1316 dst_ptr[i] *= relaxation;
1321 tmp.reinit(src,
true);
1327 preconditioner.vmult(
1330 [&](
const unsigned int start_range,
const unsigned int end_range) {
1331 const auto tmp_ptr = tmp.begin();
1333 if (relaxation == 1.0)
1336 for (std::size_t i = start_range; i < end_range; ++i)
1337 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1345 for (std::size_t i = start_range; i < end_range; ++i)
1346 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1349 [&](
const unsigned int,
const unsigned int) {
1359 typename MatrixType,
1360 typename PreconditionerType,
1361 typename VectorType,
1363 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1365 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1368 step_operations(
const MatrixType &A,
1369 const PreconditionerType &preconditioner,
1371 const VectorType &src,
1372 const double relaxation,
1375 const unsigned int i,
1376 const bool transposed)
1379 using Number =
typename VectorType::value_type;
1381 Number *dst_ptr = dst.begin();
1382 const Number *src_ptr = src.begin();
1386 preconditioner.vmult(
1389 [&](
const unsigned int start_range,
const unsigned int end_range) {
1391 if (end_range > start_range)
1392 std::memset(dst_ptr + start_range,
1394 sizeof(Number) * (end_range - start_range));
1396 [&](
const unsigned int start_range,
const unsigned int end_range) {
1397 if (relaxation == 1.0)
1401 for (std::size_t i = start_range; i < end_range; ++i)
1402 dst_ptr[i] *= relaxation;
1407 tmp.reinit(src,
true);
1408 const auto tmp_ptr = tmp.begin();
1415 [&](
const unsigned int start_range,
const unsigned int end_range) {
1418 if (end_range > start_range)
1419 std::memset(tmp_ptr + start_range,
1421 sizeof(Number) * (end_range - start_range));
1423 [&](
const unsigned int start_range,
const unsigned int end_range) {
1424 if (relaxation == 1.0)
1427 for (std::size_t i = start_range; i < end_range; ++i)
1428 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1436 for (std::size_t i = start_range; i < end_range; ++i)
1437 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1441 preconditioner.vmult(dst, tmp, [](
const auto,
const auto) {
1450 typename MatrixType,
1451 typename VectorType,
1458 !has_vmult_with_std_functions<MatrixType,
1461 VectorType> * =
nullptr>
1463 step_operations(
const MatrixType &A,
1464 const ::DiagonalMatrix<VectorType> &preconditioner,
1466 const VectorType &src,
1467 const double relaxation,
1470 const unsigned int i,
1471 const bool transposed)
1473 using Number =
typename VectorType::value_type;
1477 Number *dst_ptr = dst.begin();
1478 const Number *src_ptr = src.begin();
1479 const Number *diag_ptr = preconditioner.get_vector().begin();
1481 if (relaxation == 1.0)
1484 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1485 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1490 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1491 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1496 tmp.reinit(src,
true);
1499 Tvmult(A, tmp, dst);
1503 Number *dst_ptr = dst.begin();
1504 const Number *src_ptr = src.begin();
1505 const Number *tmp_ptr = tmp.begin();
1506 const Number *diag_ptr = preconditioner.get_vector().begin();
1508 if (relaxation == 1.0)
1511 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1512 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1517 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1519 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1526 template <
typename MatrixType,
1527 typename VectorType,
1528 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1529 has_vmult_with_std_functions<
1533 VectorType> * =
nullptr>
1535 step_operations(
const MatrixType &A,
1536 const ::DiagonalMatrix<VectorType> &preconditioner,
1538 const VectorType &src,
1539 const double relaxation,
1542 const unsigned int i,
1543 const bool transposed)
1546 using Number =
typename VectorType::value_type;
1550 Number *dst_ptr = dst.begin();
1551 const Number *src_ptr = src.begin();
1552 const Number *diag_ptr = preconditioner.get_vector().begin();
1554 if (relaxation == 1.0)
1557 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1558 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1563 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1564 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1569 tmp.reinit(src,
true);
1576 [&](
const unsigned int start_range,
const unsigned int end_range) {
1578 if (end_range > start_range)
1579 std::memset(tmp.begin() + start_range,
1581 sizeof(Number) * (end_range - start_range));
1583 [&](
const unsigned int begin,
const unsigned int end) {
1584 const Number *dst_ptr = dst.begin();
1585 const Number *src_ptr = src.begin();
1586 Number *tmp_ptr = tmp.begin();
1587 const Number *diag_ptr = preconditioner.get_vector().begin();
1591 if (relaxation == 1.0)
1594 for (std::size_t i = begin; i < end; ++i)
1596 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1601 for (std::size_t i = begin; i < end; ++i)
1602 tmp_ptr[i] = dst_ptr[i] + relaxation *
1603 (src_ptr[i] - tmp_ptr[i]) *
1645template <
typename MatrixType = SparseMatrix<
double>>
1649 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1652 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1715template <
typename MatrixType = SparseMatrix<
double>>
1719 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1722 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1767template <
typename MatrixType = SparseMatrix<
double>>
1771 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1774 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1823template <
typename MatrixType = SparseMatrix<
double>>
1827 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1830 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1887 const std::vector<size_type> &permutation,
1888 const std::vector<size_type> &inverse_permutation,
2096template <
typename MatrixType = SparseMatrix<
double>,
2097 typename VectorType = Vector<
double>,
2098 typename PreconditionerType = DiagonalMatrix<VectorType>>
2136 const unsigned int degree = 1,
2142 EigenvalueAlgorithm::lanczos,
2191 vmult(VectorType &dst,
const VectorType &src)
const;
2198 Tvmult(VectorType &dst,
const VectorType &src)
const;
2204 step(VectorType &dst,
const VectorType &src)
const;
2210 Tstep(VectorType &dst,
const VectorType &src)
const;
2313 template <
typename VectorType>
2315 set_initial_guess(VectorType &vector)
2317 vector = 1. /
std::sqrt(
static_cast<double>(vector.size()));
2318 if (vector.locally_owned_elements().is_element(0))
2322 template <
typename Number>
2330 for (
unsigned int i = 0; i < vector.
size(); ++i)
2333 const Number mean_value = vector.
mean_value();
2334 vector.
add(-mean_value);
2337 template <
typename Number>
2342 for (
unsigned int block = 0; block < vector.
n_blocks(); ++block)
2343 set_initial_guess(vector.
block(block));
2346 template <
typename Number,
typename MemorySpace>
2363 Kokkos::RangePolicy<
typename MemorySpace::kokkos_space::execution_space,
2364 Kokkos::IndexType<types::global_dof_index>>
2365 policy(0, n_local_elements);
2366 Kokkos::parallel_for(
2367 "::PreconditionChebyshev::set_initial_guess",
2370 values_ptr[i] = (i + first_local_range) % 11;
2372 const Number mean_value = vector.
mean_value();
2373 vector.
add(-mean_value);
2376 struct EigenvalueTracker
2385 std::vector<double>
values;
2392 typename PreconditionerType>
2395 VectorType &eigenvector,
2396 const PreconditionerType &preconditioner,
2397 const unsigned int n_iterations)
2399 typename VectorType::value_type eigenvalue_estimate = 0.;
2400 eigenvector /= eigenvector.l2_norm();
2402 vector1.reinit(eigenvector,
true);
2403 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2404 vector2.reinit(eigenvector,
true);
2405 for (
unsigned int i = 0; i < n_iterations; ++i)
2407 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2409 matrix.vmult(vector2, eigenvector);
2410 preconditioner.vmult(vector1, vector2);
2413 matrix.vmult(vector1, eigenvector);
2415 eigenvalue_estimate = eigenvector * vector1;
2417 vector1 /= vector1.l2_norm();
2418 eigenvector.swap(vector1);
2420 return std::abs(eigenvalue_estimate);
2427 typename PreconditionerType>
2428 EigenvalueInformation
2429 estimate_eigenvalues(
2430 const EigenvalueAlgorithmAdditionalData<PreconditionerType> &
data,
2431 const MatrixType *matrix_ptr,
2432 VectorType &solution_old,
2433 VectorType &temp_vector1,
2434 const unsigned int degree)
2438 EigenvalueInformation info{};
2440 if (
data.eig_cg_n_iterations > 0)
2444 "Need to set at least two iterations to find eigenvalues."));
2446 internal::EigenvalueTracker eigenvalue_tracker;
2451 internal::set_initial_guess(temp_vector1);
2452 data.constraints.set_zero(temp_vector1);
2463 solver.connect_eigenvalues_slot(
2464 [&eigenvalue_tracker](
const std::vector<double> &
eigenvalues) {
2468 solver.solve(*matrix_ptr,
2471 *
data.preconditioner);
2473 info.cg_iterations = control.last_step();
2475 else if (
data.eigenvalue_algorithm ==
2481 ExcMessage(
"Cannot estimate the minimal eigenvalue with the "
2482 "power iteration"));
2484 eigenvalue_tracker.values.push_back(
2487 *
data.preconditioner,
2488 data.eig_cg_n_iterations));
2494 if (eigenvalue_tracker.values.empty())
2495 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2498 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2502 info.max_eigenvalue_estimate =
2503 1.2 * eigenvalue_tracker.values.back();
2508 info.max_eigenvalue_estimate =
data.max_eigenvalue;
2509 info.min_eigenvalue_estimate =
2510 data.max_eigenvalue /
data.smoothing_range;
2524template <
typename MatrixType>
2534template <
typename VectorType>
2543template <
typename VectorType>
2550template <
typename VectorType>
2559template <
typename VectorType>
2591 const double relaxation)
2592 : relaxation(relaxation)
2601 AdditionalData add_data;
2602 relaxation = add_data.relaxation;
2616template <
typename MatrixType>
2619 const MatrixType &matrix,
2629template <
typename VectorType>
2634 std::is_same_v<size_type, typename VectorType::size_type>,
2635 "PreconditionRichardson and VectorType must have the same size_type.");
2642template <
typename VectorType>
2647 std::is_same_v<size_type, typename VectorType::size_type>,
2648 "PreconditionRichardson and VectorType must have the same size_type.");
2653template <
typename VectorType>
2658 std::is_same_v<size_type, typename VectorType::size_type>,
2659 "PreconditionRichardson and VectorType must have the same size_type.");
2666template <
typename VectorType>
2671 std::is_same_v<size_type, typename VectorType::size_type>,
2672 "PreconditionRichardson and VectorType must have the same size_type.");
2693template <
typename MatrixType,
typename PreconditionerType>
2696 const MatrixType &rA,
2697 const AdditionalData ¶meters)
2700 eigenvalues_are_initialized =
false;
2704 this->data = parameters;
2708template <
typename MatrixType,
typename PreconditionerType>
2712 eigenvalues_are_initialized =
false;
2714 data.relaxation = 1.0;
2715 data.preconditioner =
nullptr;
2718template <
typename MatrixType,
typename PreconditionerType>
2727template <
typename MatrixType,
typename PreconditionerType>
2736template <
typename MatrixType,
typename PreconditionerType>
2737template <
typename VectorType>
2741 const VectorType &src)
const
2746 if (eigenvalues_are_initialized ==
false)
2747 estimate_eigenvalues(src);
2751 for (
unsigned int i = 0; i <
data.n_iterations; ++i)
2752 internal::PreconditionRelaxation::step_operations(*A,
2753 *
data.preconditioner,
2763template <
typename MatrixType,
typename PreconditionerType>
2764template <
typename VectorType>
2768 const VectorType &src)
const
2773 if (eigenvalues_are_initialized ==
false)
2774 estimate_eigenvalues(src);
2778 for (
unsigned int i = 0; i <
data.n_iterations; ++i)
2779 internal::PreconditionRelaxation::step_operations(
2780 *A, *
data.preconditioner, dst, src,
data.relaxation, tmp1, tmp2, i,
true);
2783template <
typename MatrixType,
typename PreconditionerType>
2784template <
typename VectorType>
2788 const VectorType &src)
const
2793 if (eigenvalues_are_initialized ==
false)
2794 estimate_eigenvalues(src);
2798 for (
unsigned int i = 1; i <=
data.n_iterations; ++i)
2799 internal::PreconditionRelaxation::step_operations(*A,
2800 *
data.preconditioner,
2810template <
typename MatrixType,
typename PreconditionerType>
2811template <
typename VectorType>
2815 const VectorType &src)
const
2820 if (eigenvalues_are_initialized ==
false)
2821 estimate_eigenvalues(src);
2825 for (
unsigned int i = 1; i <=
data.n_iterations; ++i)
2826 internal::PreconditionRelaxation::step_operations(
2827 *A, *
data.preconditioner, dst, src,
data.relaxation, tmp1, tmp2, i,
true);
2830template <
typename MatrixType,
typename PreconditionerType>
2831template <
typename VectorType>
2834 const VectorType &src)
const
2838 EigenvalueInformation info;
2840 if (
data.relaxation == 0.0)
2844 solution_old.reinit(src);
2845 temp_vector1.reinit(src,
true);
2847 info = internal::estimate_eigenvalues<MatrixType>(
2848 data, A, solution_old, temp_vector1,
data.n_iterations);
2850 const double alpha =
2851 (
data.smoothing_range > 1. ?
2852 info.max_eigenvalue_estimate /
data.smoothing_range :
2853 std::min(0.9 * info.max_eigenvalue_estimate,
2854 info.min_eigenvalue_estimate));
2857 ->
data.relaxation = 2.0 / (alpha + info.max_eigenvalue_estimate);
2861 ->eigenvalues_are_initialized =
true;
2866template <
typename MatrixType,
typename PreconditionerType>
2870 return data.relaxation;
2876template <
typename MatrixType>
2879 const AdditionalData ¶meters_in)
2883 parameters_in.relaxation != 0.0,
2885 "Relaxation cannot automatically be determined by PreconditionJacobi."));
2887 AdditionalData parameters;
2888 parameters.relaxation = 1.0;
2889 parameters.n_iterations = parameters_in.n_iterations;
2890 parameters.preconditioner =
2891 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2893 this->BaseClass::initialize(A, parameters);
2898template <
typename MatrixType>
2901 const AdditionalData ¶meters_in)
2905 parameters_in.relaxation != 0.0,
2907 "Relaxation cannot automatically be determined by PreconditionSOR."));
2909 AdditionalData parameters;
2910 parameters.relaxation = 1.0;
2911 parameters.n_iterations = parameters_in.n_iterations;
2912 parameters.preconditioner =
2913 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2915 this->BaseClass::initialize(A, parameters);
2920template <
typename MatrixType>
2923 const AdditionalData ¶meters_in)
2927 parameters_in.relaxation != 0.0,
2929 "Relaxation cannot automatically be determined by PreconditionSSOR."));
2931 AdditionalData parameters;
2932 parameters.relaxation = 1.0;
2933 parameters.n_iterations = parameters_in.n_iterations;
2934 parameters.preconditioner =
2935 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2937 this->BaseClass::initialize(A, parameters);
2944template <
typename MatrixType>
2947 const MatrixType &A,
2948 const std::vector<size_type> &p,
2949 const std::vector<size_type> &ip,
2950 const typename BaseClass::AdditionalData ¶meters_in)
2954 parameters_in.relaxation != 0.0,
2956 "Relaxation cannot automatically be determined by PreconditionPSOR."));
2958 typename BaseClass::AdditionalData parameters;
2959 parameters.relaxation = 1.0;
2960 parameters.n_iterations = parameters_in.n_iterations;
2961 parameters.preconditioner =
2962 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2964 this->BaseClass::initialize(A, parameters);
2968template <
typename MatrixType>
2971 const AdditionalData &additional_data)
2974 additional_data.permutation,
2975 additional_data.inverse_permutation,
2976 additional_data.parameters);
2979template <
typename MatrixType>
2981 const std::vector<size_type> &permutation,
2982 const std::vector<size_type> &inverse_permutation,
2985 : permutation(permutation)
2986 , inverse_permutation(inverse_permutation)
2987 , parameters(parameters)
2994template <
typename MatrixType,
typename VectorType>
2996 const MatrixType &M,
2997 const function_ptr method)
2999 , precondition(method)
3004template <
typename MatrixType,
typename VectorType>
3008 const VectorType &src)
const
3010 (
matrix.*precondition)(dst, src);
3018 template <
typename PreconditionerType>
3021 const double smoothing_range,
3022 const unsigned int eig_cg_n_iterations,
3023 const double eig_cg_residual,
3024 const double max_eigenvalue,
3026 : smoothing_range(smoothing_range)
3027 , eig_cg_n_iterations(eig_cg_n_iterations)
3028 , eig_cg_residual(eig_cg_residual)
3029 , max_eigenvalue(max_eigenvalue)
3030 , eigenvalue_algorithm(eigenvalue_algorithm)
3035 template <
typename PreconditionerType>
3036 inline EigenvalueAlgorithmAdditionalData<PreconditionerType> &
3037 EigenvalueAlgorithmAdditionalData<PreconditionerType>::operator=(
3038 const EigenvalueAlgorithmAdditionalData &other_data)
3040 smoothing_range = other_data.smoothing_range;
3041 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3042 eig_cg_residual = other_data.eig_cg_residual;
3043 max_eigenvalue = other_data.max_eigenvalue;
3044 preconditioner = other_data.preconditioner;
3045 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3046 constraints.copy_from(other_data.constraints);
3052template <
typename MatrixType,
typename PreconditionerType>
3055 const unsigned int n_iterations,
3056 const double smoothing_range,
3057 const unsigned int eig_cg_n_iterations,
3058 const double eig_cg_residual,
3059 const double max_eigenvalue,
3060 const EigenvalueAlgorithm eigenvalue_algorithm)
3061 :
internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3063 eig_cg_n_iterations,
3066 eigenvalue_algorithm)
3067 , relaxation(relaxation)
3068 , n_iterations(n_iterations)
3077 namespace PreconditionChebyshevImplementation
3085 template <
typename VectorType,
typename PreconditionerType>
3087 vector_updates(
const VectorType &rhs,
3088 const PreconditionerType &preconditioner,
3089 const unsigned int iteration_index,
3090 const double factor1,
3091 const double factor2,
3092 VectorType &solution_old,
3093 VectorType &temp_vector1,
3094 VectorType &temp_vector2,
3095 VectorType &solution)
3097 if (iteration_index == 0)
3099 solution.equ(factor2, rhs);
3100 preconditioner.vmult(solution_old, solution);
3102 else if (iteration_index == 1)
3105 temp_vector1.sadd(-1.0, 1.0, rhs);
3106 preconditioner.vmult(solution_old, temp_vector1);
3109 solution_old.sadd(factor2, 1 + factor1, solution);
3114 temp_vector1.sadd(-1.0, 1.0, rhs);
3115 preconditioner.vmult(temp_vector2, temp_vector1);
3118 solution_old.sadd(-factor1, factor2, temp_vector2);
3119 solution_old.add(1 + factor1, solution);
3122 solution.swap(solution_old);
3128 typename PreconditionerType,
3130 !has_vmult_with_std_functions_for_precondition<
3137 const PreconditionerType &preconditioner,
3138 const unsigned int iteration_index,
3139 const double factor1_,
3140 const double factor2_,
3149 const Number factor1 = factor1_;
3150 const Number factor1_plus_1 = 1. + factor1_;
3151 const Number factor2 = factor2_;
3153 if (iteration_index == 0)
3156 preconditioner.vmult(solution_old, rhs);
3159 const auto solution_old_ptr = solution_old.
begin();
3162 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
3164 else if (iteration_index == 1)
3167 temp_vector1.
sadd(-1.0, 1.0, rhs);
3169 preconditioner.vmult(solution_old, temp_vector1);
3172 const auto solution_ptr = solution.
begin();
3173 const auto solution_old_ptr = solution_old.
begin();
3176 solution_old_ptr[i] =
3177 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
3182 temp_vector1.
sadd(-1.0, 1.0, rhs);
3184 preconditioner.vmult(temp_vector2, temp_vector1);
3187 const auto solution_ptr = solution.
begin();
3188 const auto solution_old_ptr = solution_old.
begin();
3189 const auto temp_vector2_ptr = temp_vector2.
begin();
3192 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3193 factor1 * solution_old_ptr[i] +
3194 temp_vector2_ptr[i] * factor2;
3197 solution.
swap(solution_old);
3202 typename PreconditionerType,
3204 has_vmult_with_std_functions_for_precondition<
3211 const PreconditionerType &preconditioner,
3212 const unsigned int iteration_index,
3213 const double factor1_,
3214 const double factor2_,
3223 const Number factor1 = factor1_;
3224 const Number factor1_plus_1 = 1. + factor1_;
3225 const Number factor2 = factor2_;
3227 const auto rhs_ptr = rhs.
begin();
3228 const auto temp_vector1_ptr = temp_vector1.
begin();
3229 const auto temp_vector2_ptr = temp_vector2.
begin();
3230 const auto solution_ptr = solution.
begin();
3231 const auto solution_old_ptr = solution_old.
begin();
3233 if (iteration_index == 0)
3235 preconditioner.vmult(
3238 [&](
const auto start_range,
const auto end_range) {
3239 if (end_range > start_range)
3240 std::memset(solution_ptr + start_range,
3242 sizeof(Number) * (end_range - start_range));
3244 [&](
const auto begin,
const auto end) {
3246 for (std::size_t i = begin; i <
end; ++i)
3247 solution_ptr[i] *= factor2;
3252 preconditioner.vmult(
3255 [&](
const auto begin,
const auto end) {
3257 std::memset(temp_vector2_ptr + begin,
3259 sizeof(Number) * (end - begin));
3262 for (std::size_t i = begin; i <
end; ++i)
3263 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3265 [&](
const auto begin,
const auto end) {
3266 if (iteration_index == 1)
3269 for (std::size_t i = begin; i <
end; ++i)
3270 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
3271 factor2 * temp_vector2_ptr[i];
3276 for (std::size_t i = begin; i <
end; ++i)
3277 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3278 factor1 * solution_old_ptr[i] +
3279 factor2 * temp_vector2_ptr[i];
3284 if (iteration_index > 0)
3286 solution_old.
swap(temp_vector2);
3287 solution_old.
swap(solution);
3294 template <
typename Number>
3295 struct VectorUpdater
3297 VectorUpdater(
const Number *rhs,
3298 const Number *matrix_diagonal_inverse,
3299 const unsigned int iteration_index,
3300 const Number factor1,
3301 const Number factor2,
3302 Number *solution_old,
3306 , matrix_diagonal_inverse(matrix_diagonal_inverse)
3307 , iteration_index(iteration_index)
3310 , solution_old(solution_old)
3311 , tmp_vector(tmp_vector)
3312 , solution(solution)
3316 apply_to_subrange(
const std::size_t begin,
const std::size_t end)
const
3322 const Number factor1 = this->factor1;
3323 const Number factor1_plus_1 = 1. + this->factor1;
3324 const Number factor2 = this->factor2;
3325 if (iteration_index == 0)
3328 for (std::size_t i = begin; i <
end; ++i)
3329 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
3331 else if (iteration_index == 1)
3335 for (std::size_t i = begin; i <
end; ++i)
3339 factor1_plus_1 * solution[i] +
3340 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3347 for (std::size_t i = begin; i <
end; ++i)
3354 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
3355 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3360 const Number *matrix_diagonal_inverse;
3361 const unsigned int iteration_index;
3362 const Number factor1;
3363 const Number factor2;
3364 mutable Number *solution_old;
3365 mutable Number *tmp_vector;
3366 mutable Number *solution;
3369 template <
typename Number>
3372 VectorUpdatesRange(
const VectorUpdater<Number> &updater,
3373 const std::size_t
size)
3377 VectorUpdatesRange::apply_to_subrange(0,
size);
3385 ~VectorUpdatesRange()
override =
default;
3388 apply_to_subrange(
const std::size_t begin,
3389 const std::size_t end)
const override
3391 updater.apply_to_subrange(begin, end);
3394 const VectorUpdater<Number> &updater;
3398 template <
typename Number>
3401 const ::Vector<Number> &rhs,
3403 const unsigned int iteration_index,
3404 const double factor1,
3405 const double factor2,
3411 VectorUpdater<Number> upd(rhs.begin(),
3412 jacobi.get_vector().begin(),
3416 solution_old.
begin(),
3417 temp_vector1.
begin(),
3419 VectorUpdatesRange<Number>(upd, rhs.size());
3422 if (iteration_index == 0)
3429 solution.
swap(temp_vector1);
3430 solution_old.
swap(temp_vector1);
3435 template <
typename Number>
3439 const ::DiagonalMatrix<
3441 const unsigned int iteration_index,
3442 const double factor1,
3443 const double factor2,
3451 VectorUpdater<Number> upd(rhs.
begin(),
3452 jacobi.get_vector().begin(),
3456 solution_old.
begin(),
3457 temp_vector1.
begin(),
3462 if (iteration_index == 0)
3469 solution.
swap(temp_vector1);
3470 solution_old.
swap(temp_vector1);
3481 typename PreconditionerType,
3485 PreconditionerType> &&
3486 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3488 has_vmult_with_std_functions_for_precondition<
MatrixType,
3492 vmult_and_update(
const MatrixType &matrix,
3493 const PreconditionerType &preconditioner,
3494 const VectorType &rhs,
3495 const unsigned int iteration_index,
3496 const double factor1,
3497 const double factor2,
3498 VectorType &solution,
3499 VectorType &solution_old,
3500 VectorType &temp_vector1,
3501 VectorType &temp_vector2)
3503 if (iteration_index > 0)
3504 matrix.vmult(temp_vector1, solution);
3521 typename PreconditionerType,
3525 PreconditionerType> &&
3526 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3528 has_vmult_with_std_functions_for_precondition<
MatrixType,
3532 vmult_and_update(
const MatrixType &matrix,
3533 const PreconditionerType &preconditioner,
3534 const VectorType &rhs,
3535 const unsigned int iteration_index,
3536 const double factor1_,
3537 const double factor2_,
3538 VectorType &solution,
3539 VectorType &solution_old,
3540 VectorType &temp_vector1,
3541 VectorType &temp_vector2)
3543 using Number =
typename VectorType::value_type;
3545 const Number factor1 = factor1_;
3546 const Number factor1_plus_1 = 1. + factor1_;
3547 const Number factor2 = factor2_;
3549 if (iteration_index == 0)
3551 preconditioner.vmult(
3554 [&](
const unsigned int start_range,
const unsigned int end_range) {
3556 if (end_range > start_range)
3557 std::memset(solution.begin() + start_range,
3559 sizeof(Number) * (end_range - start_range));
3561 [&](
const unsigned int start_range,
const unsigned int end_range) {
3562 const auto solution_ptr = solution.begin();
3565 for (std::size_t i = start_range; i < end_range; ++i)
3566 solution_ptr[i] *= factor2;
3571 temp_vector1.reinit(rhs,
true);
3572 temp_vector2.reinit(rhs,
true);
3578 [&](
const unsigned int start_range,
const unsigned int end_range) {
3581 if (end_range > start_range)
3582 std::memset(temp_vector1.begin() + start_range,
3584 sizeof(Number) * (end_range - start_range));
3586 [&](
const unsigned int start_range,
const unsigned int end_range) {
3587 const auto rhs_ptr = rhs.begin();
3588 const auto tmp_ptr = temp_vector1.begin();
3591 for (std::size_t i = start_range; i < end_range; ++i)
3592 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3596 preconditioner.vmult(
3599 [&](
const unsigned int start_range,
const unsigned int end_range) {
3602 if (end_range > start_range)
3603 std::memset(temp_vector2.begin() + start_range,
3605 sizeof(Number) * (end_range - start_range));
3607 [&](
const unsigned int start_range,
const unsigned int end_range) {
3608 const auto solution_ptr = solution.begin();
3609 const auto solution_old_ptr = solution_old.begin();
3610 const auto tmp_ptr = temp_vector2.begin();
3612 if (iteration_index == 1)
3615 for (std::size_t i = start_range; i < end_range; ++i)
3617 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3622 for (std::size_t i = start_range; i < end_range; ++i)
3623 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3624 factor1 * solution_old_ptr[i] +
3625 factor2 * tmp_ptr[i];
3629 solution.swap(temp_vector2);
3630 solution_old.swap(temp_vector2);
3638 typename PreconditionerType,
3639 std::enable_if_t<has_vmult_with_std_functions<
MatrixType,
3641 PreconditionerType>,
3644 vmult_and_update(
const MatrixType &matrix,
3645 const PreconditionerType &preconditioner,
3646 const VectorType &rhs,
3647 const unsigned int iteration_index,
3648 const double factor1,
3649 const double factor2,
3650 VectorType &solution,
3651 VectorType &solution_old,
3652 VectorType &temp_vector1,
3655 using Number =
typename VectorType::value_type;
3656 VectorUpdater<Number> updater(rhs.begin(),
3657 preconditioner.get_vector().begin(),
3661 solution_old.begin(),
3662 temp_vector1.begin(),
3664 if (iteration_index > 0)
3668 [&](
const unsigned int start_range,
const unsigned int end_range) {
3671 if (end_range > start_range)
3672 std::memset(temp_vector1.begin() + start_range,
3674 sizeof(Number) * (end_range - start_range));
3676 [&](
const unsigned int start_range,
const unsigned int end_range) {
3677 if (end_range > start_range)
3678 updater.apply_to_subrange(start_range, end_range);
3681 updater.apply_to_subrange(0U, solution.locally_owned_size());
3684 if (iteration_index == 0)
3691 solution.swap(temp_vector1);
3692 solution_old.swap(temp_vector1);
3696 template <
typename MatrixType,
typename PreconditionerType>
3698 initialize_preconditioner(
3699 const MatrixType &matrix,
3700 std::shared_ptr<PreconditionerType> &preconditioner)
3703 (void)preconditioner;
3707 template <
typename MatrixType,
typename VectorType>
3709 initialize_preconditioner(
3710 const MatrixType &matrix,
3713 if (preconditioner.get() ==
nullptr || preconditioner->m() !=
matrix.m())
3715 if (preconditioner.get() ==
nullptr)
3717 std::make_shared<::DiagonalMatrix<VectorType>>();
3720 preconditioner->m() == 0,
3722 "Preconditioner appears to be initialized but not sized correctly"));
3725 if (preconditioner->m() !=
matrix.m())
3727 preconditioner->get_vector().reinit(
matrix.m());
3728 for (
typename VectorType::size_type i = 0; i <
matrix.m(); ++i)
3729 preconditioner->get_vector()(i) = 1. /
matrix.el(i, i);
3738template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3741 const double smoothing_range,
3742 const unsigned int eig_cg_n_iterations,
3743 const double eig_cg_residual,
3744 const double max_eigenvalue,
3745 const EigenvalueAlgorithm eigenvalue_algorithm,
3746 const PolynomialType polynomial_type)
3747 :
internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3749 eig_cg_n_iterations,
3752 eigenvalue_algorithm)
3754 , polynomial_type(polynomial_type)
3759template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3764 , eigenvalues_are_initialized(false)
3767 std::is_same_v<size_type, typename VectorType::size_type>,
3768 "PreconditionChebyshev and VectorType must have the same size_type.");
3773template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3776 const MatrixType &matrix,
3777 const AdditionalData &additional_data)
3780 data = additional_data;
3782 ExcMessage(
"The degree of the Chebyshev method must be positive."));
3783 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3784 matrix,
data.preconditioner);
3785 eigenvalues_are_initialized =
false;
3790template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3794 eigenvalues_are_initialized =
false;
3795 theta = delta = 1.0;
3796 matrix_ptr =
nullptr;
3799 solution_old.reinit(empty_vector);
3800 temp_vector1.reinit(empty_vector);
3801 temp_vector2.reinit(empty_vector);
3803 data.preconditioner.reset();
3808template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3815 solution_old.reinit(src);
3816 temp_vector1.reinit(src,
true);
3818 auto info = internal::estimate_eigenvalues<MatrixType>(
3819 data, matrix_ptr, solution_old, temp_vector1,
data.degree);
3821 const double alpha = (
data.smoothing_range > 1. ?
3822 info.max_eigenvalue_estimate /
data.smoothing_range :
3823 std::min(0.9 * info.max_eigenvalue_estimate,
3824 info.min_eigenvalue_estimate));
3833 const double actual_range = info.max_eigenvalue_estimate / alpha;
3834 const double sigma = (1. -
std::sqrt(1. / actual_range)) /
3836 const double eps =
data.smoothing_range;
3841 1 +
static_cast<unsigned int>(
3846 info.degree =
data.degree;
3851 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3852 (info.max_eigenvalue_estimate) :
3853 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3856 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3860 using NumberType =
typename VectorType::value_type;
3868 (std::is_same_v<VectorType,
3871 temp_vector2.
reinit(src, true);
3875 temp_vector2.reinit(empty_vector);
3880 ->eigenvalues_are_initialized =
true;
3887template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3890 VectorType &solution,
3891 const VectorType &rhs)
const
3893 std::lock_guard<std::mutex> lock(mutex);
3894 if (eigenvalues_are_initialized ==
false)
3895 estimate_eigenvalues(rhs);
3897 internal::PreconditionChebyshevImplementation::vmult_and_update(
3899 *
data.preconditioner,
3903 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3904 (4. / (3. * delta)) :
3916 double rhok = delta / theta, sigma = theta / delta;
3917 for (
unsigned int k = 0; k <
data.degree - 1; ++k)
3919 double factor1 = 0.0;
3920 double factor2 = 0.0;
3922 if (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3924 factor1 = (2 * k + 1.) / (2 * k + 5.);
3925 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3929 const double rhokp = 1. / (2. * sigma - rhok);
3930 factor1 = rhokp * rhok;
3931 factor2 = 2. * rhokp / delta;
3935 internal::PreconditionChebyshevImplementation::vmult_and_update(
3937 *
data.preconditioner,
3951template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3954 VectorType &solution,
3955 const VectorType &rhs)
const
3957 std::lock_guard<std::mutex> lock(mutex);
3958 if (eigenvalues_are_initialized ==
false)
3959 estimate_eigenvalues(rhs);
3961 internal::PreconditionChebyshevImplementation::vector_updates(
3963 *
data.preconditioner,
3966 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3967 (4. / (3. * delta)) :
3977 double rhok = delta / theta, sigma = theta / delta;
3978 for (
unsigned int k = 0; k <
data.degree - 1; ++k)
3980 double factor1 = 0.0;
3981 double factor2 = 0.0;
3983 if (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3985 factor1 = (2 * k + 1.) / (2 * k + 5.);
3986 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3990 const double rhokp = 1. / (2. * sigma - rhok);
3991 factor1 = rhokp * rhok;
3992 factor2 = 2. * rhokp / delta;
3996 matrix_ptr->Tvmult(temp_vector1, solution);
3997 internal::PreconditionChebyshevImplementation::vector_updates(
3999 *
data.preconditioner,
4012template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4015 VectorType &solution,
4016 const VectorType &rhs)
const
4018 std::lock_guard<std::mutex> lock(mutex);
4019 if (eigenvalues_are_initialized ==
false)
4020 estimate_eigenvalues(rhs);
4022 internal::PreconditionChebyshevImplementation::vmult_and_update(
4024 *
data.preconditioner,
4028 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4029 (4. / (3. * delta)) :
4039 double rhok = delta / theta, sigma = theta / delta;
4040 for (
unsigned int k = 0; k <
data.degree - 1; ++k)
4042 double factor1 = 0.0;
4043 double factor2 = 0.0;
4045 if (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4047 factor1 = (2 * k + 1.) / (2 * k + 5.);
4048 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4052 const double rhokp = 1. / (2. * sigma - rhok);
4053 factor1 = rhokp * rhok;
4054 factor2 = 2. * rhokp / delta;
4058 internal::PreconditionChebyshevImplementation::vmult_and_update(
4060 *
data.preconditioner,
4074template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4077 VectorType &solution,
4078 const VectorType &rhs)
const
4080 std::lock_guard<std::mutex> lock(mutex);
4081 if (eigenvalues_are_initialized ==
false)
4082 estimate_eigenvalues(rhs);
4084 matrix_ptr->Tvmult(temp_vector1, solution);
4085 internal::PreconditionChebyshevImplementation::vector_updates(
4087 *
data.preconditioner,
4090 (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4091 (4. / (3. * delta)) :
4101 double rhok = delta / theta, sigma = theta / delta;
4102 for (
unsigned int k = 0; k <
data.degree - 1; ++k)
4104 double factor1 = 0.0;
4105 double factor2 = 0.0;
4107 if (
data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4109 factor1 = (2 * k + 1.) / (2 * k + 5.);
4110 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4114 const double rhokp = 1. / (2. * sigma - rhok);
4115 factor1 = rhokp * rhok;
4116 factor2 = 2. * rhokp / delta;
4120 matrix_ptr->Tvmult(temp_vector1, solution);
4121 internal::PreconditionChebyshevImplementation::vector_updates(
4123 *
data.preconditioner,
4136template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4139 PreconditionerType>::size_type
4143 return matrix_ptr->m();
4148template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4151 PreconditionerType>::size_type
4155 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
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()
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)