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;
1298 Number *dst_ptr = dst.begin();
1299 const Number *src_ptr = src.begin();
1301 preconditioner.vmult(
1304 [&](
const unsigned int start_range,
const unsigned int end_range) {
1306 if (end_range > start_range)
1307 std::memset(dst.begin() + start_range,
1309 sizeof(Number) * (end_range - start_range));
1311 [&](
const unsigned int start_range,
const unsigned int end_range) {
1312 if (relaxation == 1.0)
1315 const auto src_ptr = src.begin();
1316 const auto dst_ptr = dst.begin();
1319 for (std::size_t i = start_range; i < end_range; ++i)
1320 dst_ptr[i] *= relaxation;
1325 tmp.reinit(src,
true);
1331 preconditioner.vmult(
1334 [&](
const unsigned int start_range,
const unsigned int end_range) {
1335 const auto src_ptr = src.begin();
1336 const auto tmp_ptr = tmp.begin();
1338 if (relaxation == 1.0)
1341 for (std::size_t i = start_range; i < end_range; ++i)
1342 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1350 for (std::size_t i = start_range; i < end_range; ++i)
1351 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1354 [&](
const unsigned int,
const unsigned int) {
1364 typename MatrixType,
1365 typename PreconditionerType,
1366 typename VectorType,
1368 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1370 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1373 step_operations(
const MatrixType &A,
1374 const PreconditionerType &preconditioner,
1376 const VectorType &src,
1377 const double relaxation,
1380 const unsigned int i,
1381 const bool transposed)
1384 using Number =
typename VectorType::value_type;
1388 Number *dst_ptr = dst.begin();
1389 const Number *src_ptr = src.begin();
1391 preconditioner.vmult(
1394 [&](
const unsigned int start_range,
const unsigned int end_range) {
1396 if (end_range > start_range)
1397 std::memset(dst.begin() + start_range,
1399 sizeof(Number) * (end_range - start_range));
1401 [&](
const unsigned int start_range,
const unsigned int end_range) {
1402 if (relaxation == 1.0)
1405 const auto src_ptr = src.begin();
1406 const auto dst_ptr = dst.begin();
1409 for (std::size_t i = start_range; i < end_range; ++i)
1410 dst_ptr[i] *= relaxation;
1415 tmp.reinit(src,
true);
1422 [&](
const unsigned int start_range,
const unsigned int end_range) {
1425 if (end_range > start_range)
1426 std::memset(tmp.begin() + start_range,
1428 sizeof(Number) * (end_range - start_range));
1430 [&](
const unsigned int start_range,
const unsigned int end_range) {
1431 const auto src_ptr = src.begin();
1432 const auto tmp_ptr = tmp.begin();
1434 if (relaxation == 1.0)
1437 for (std::size_t i = start_range; i < end_range; ++i)
1438 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1446 for (std::size_t i = start_range; i < end_range; ++i)
1447 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1451 preconditioner.vmult(dst, tmp, [](
const auto,
const auto) {
1459 template <
typename MatrixType,
1460 typename VectorType,
1461 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1462 !has_vmult_with_std_functions<
1466 VectorType> * =
nullptr>
1468 step_operations(
const MatrixType &A,
1469 const ::DiagonalMatrix<VectorType> &preconditioner,
1471 const VectorType &src,
1472 const double relaxation,
1475 const unsigned int i,
1476 const bool transposed)
1478 using Number =
typename VectorType::value_type;
1482 Number *dst_ptr = dst.begin();
1483 const Number *src_ptr = src.begin();
1484 const Number *diag_ptr = preconditioner.get_vector().begin();
1486 if (relaxation == 1.0)
1489 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1490 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1495 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1496 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1501 tmp.reinit(src,
true);
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();
1509 Tvmult(A, tmp, dst);
1513 if (relaxation == 1.0)
1516 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1517 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1522 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1524 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1531 template <
typename MatrixType,
1532 typename VectorType,
1533 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1534 has_vmult_with_std_functions<
1538 VectorType> * =
nullptr>
1540 step_operations(
const MatrixType &A,
1541 const ::DiagonalMatrix<VectorType> &preconditioner,
1543 const VectorType &src,
1544 const double relaxation,
1547 const unsigned int i,
1548 const bool transposed)
1551 using Number =
typename VectorType::value_type;
1555 Number *dst_ptr = dst.begin();
1556 const Number *src_ptr = src.begin();
1557 const Number *diag_ptr = preconditioner.get_vector().begin();
1559 if (relaxation == 1.0)
1562 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1563 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1568 for (
unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1569 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1574 tmp.reinit(src,
true);
1581 [&](
const unsigned int start_range,
const unsigned int end_range) {
1583 if (end_range > start_range)
1584 std::memset(tmp.begin() + start_range,
1586 sizeof(Number) * (end_range - start_range));
1588 [&](
const unsigned int begin,
const unsigned int end) {
1589 const Number *dst_ptr = dst.begin();
1590 const Number *src_ptr = src.begin();
1591 Number *tmp_ptr = tmp.begin();
1592 const Number *diag_ptr = preconditioner.get_vector().begin();
1596 if (relaxation == 1.0)
1599 for (std::size_t i = begin; i < end; ++i)
1601 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1606 for (std::size_t i = begin; i < end; ++i)
1607 tmp_ptr[i] = dst_ptr[i] + relaxation *
1608 (src_ptr[i] - tmp_ptr[i]) *
1650template <
typename MatrixType = SparseMatrix<
double>>
1654 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1657 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1720template <
typename MatrixType = SparseMatrix<
double>>
1724 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1727 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1772template <
typename MatrixType = SparseMatrix<
double>>
1776 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1779 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1828template <
typename MatrixType = SparseMatrix<
double>>
1832 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1835 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1892 const std::vector<size_type> &permutation,
1893 const std::vector<size_type> &inverse_permutation,
2101template <
typename MatrixType = SparseMatrix<
double>,
2102 typename VectorType = Vector<
double>,
2103 typename PreconditionerType = DiagonalMatrix<VectorType>>
2141 const unsigned int degree = 1,
2147 EigenvalueAlgorithm::lanczos,
2196 vmult(VectorType &dst,
const VectorType &src)
const;
2203 Tvmult(VectorType &dst,
const VectorType &src)
const;
2209 step(VectorType &dst,
const VectorType &src)
const;
2215 Tstep(VectorType &dst,
const VectorType &src)
const;
2318 template <
typename VectorType>
2320 set_initial_guess(VectorType &vector)
2322 vector = 1. /
std::sqrt(
static_cast<double>(vector.size()));
2323 if (vector.locally_owned_elements().is_element(0))
2327 template <
typename Number>
2335 for (
unsigned int i = 0; i < vector.
size(); ++i)
2338 const Number mean_value = vector.
mean_value();
2339 vector.
add(-mean_value);
2342 template <
typename Number>
2347 for (
unsigned int block = 0; block < vector.
n_blocks(); ++block)
2348 set_initial_guess(vector.
block(block));
2351 template <
typename Number,
typename MemorySpace>
2368 Kokkos::RangePolicy<
typename MemorySpace::kokkos_space::execution_space,
2369 Kokkos::IndexType<types::global_dof_index>>
2370 policy(0, n_local_elements);
2371 Kokkos::parallel_for(
2372 "::PreconditionChebyshev::set_initial_guess",
2375 values_ptr[i] = (i + first_local_range) % 11;
2377 const Number mean_value = vector.
mean_value();
2378 vector.
add(-mean_value);
2381 struct EigenvalueTracker
2390 std::vector<double>
values;
2397 typename PreconditionerType>
2400 VectorType &eigenvector,
2401 const PreconditionerType &preconditioner,
2402 const unsigned int n_iterations)
2404 typename VectorType::value_type eigenvalue_estimate = 0.;
2405 eigenvector /= eigenvector.l2_norm();
2407 vector1.reinit(eigenvector,
true);
2408 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2409 vector2.reinit(eigenvector,
true);
2410 for (
unsigned int i = 0; i < n_iterations; ++i)
2412 if (!std::is_same_v<PreconditionerType, PreconditionIdentity>)
2414 matrix.vmult(vector2, eigenvector);
2415 preconditioner.vmult(vector1, vector2);
2418 matrix.vmult(vector1, eigenvector);
2420 eigenvalue_estimate = eigenvector * vector1;
2422 vector1 /= vector1.l2_norm();
2423 eigenvector.swap(vector1);
2425 return std::abs(eigenvalue_estimate);
2432 typename PreconditionerType>
2433 EigenvalueInformation
2434 estimate_eigenvalues(
2435 const EigenvalueAlgorithmAdditionalData<PreconditionerType> &data,
2436 const MatrixType *matrix_ptr,
2437 VectorType &solution_old,
2438 VectorType &temp_vector1,
2439 const unsigned int degree)
2443 EigenvalueInformation info{};
2445 if (data.eig_cg_n_iterations > 0)
2447 Assert(data.eig_cg_n_iterations > 2,
2449 "Need to set at least two iterations to find eigenvalues."));
2451 internal::EigenvalueTracker eigenvalue_tracker;
2456 internal::set_initial_guess(temp_vector1);
2457 data.constraints.set_zero(temp_vector1);
2468 solver.connect_eigenvalues_slot(
2469 [&eigenvalue_tracker](
const std::vector<double> &
eigenvalues) {
2473 solver.solve(*matrix_ptr,
2476 *data.preconditioner);
2478 info.cg_iterations = control.last_step();
2480 else if (data.eigenvalue_algorithm ==
2486 ExcMessage(
"Cannot estimate the minimal eigenvalue with the "
2487 "power iteration"));
2489 eigenvalue_tracker.values.push_back(
2492 *data.preconditioner,
2493 data.eig_cg_n_iterations));
2499 if (eigenvalue_tracker.values.empty())
2500 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2503 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2507 info.max_eigenvalue_estimate =
2508 1.2 * eigenvalue_tracker.values.back();
2513 info.max_eigenvalue_estimate = data.max_eigenvalue;
2514 info.min_eigenvalue_estimate =
2515 data.max_eigenvalue / data.smoothing_range;
2529template <
typename MatrixType>
2539template <
typename VectorType>
2548template <
typename VectorType>
2555template <
typename VectorType>
2564template <
typename VectorType>
2596 const double relaxation)
2597 : relaxation(relaxation)
2606 AdditionalData add_data;
2607 relaxation = add_data.relaxation;
2621template <
typename MatrixType>
2624 const MatrixType &matrix,
2634template <
typename VectorType>
2639 std::is_same_v<size_type, typename VectorType::size_type>,
2640 "PreconditionRichardson and VectorType must have the same size_type.");
2647template <
typename VectorType>
2652 std::is_same_v<size_type, typename VectorType::size_type>,
2653 "PreconditionRichardson and VectorType must have the same size_type.");
2658template <
typename VectorType>
2663 std::is_same_v<size_type, typename VectorType::size_type>,
2664 "PreconditionRichardson and VectorType must have the same size_type.");
2671template <
typename VectorType>
2676 std::is_same_v<size_type, typename VectorType::size_type>,
2677 "PreconditionRichardson and VectorType must have the same size_type.");
2698template <
typename MatrixType,
typename PreconditionerType>
2701 const MatrixType &rA,
2702 const AdditionalData ¶meters)
2705 eigenvalues_are_initialized =
false;
2709 this->data = parameters;
2713template <
typename MatrixType,
typename PreconditionerType>
2717 eigenvalues_are_initialized =
false;
2719 data.relaxation = 1.0;
2720 data.preconditioner =
nullptr;
2723template <
typename MatrixType,
typename PreconditionerType>
2732template <
typename MatrixType,
typename PreconditionerType>
2741template <
typename MatrixType,
typename PreconditionerType>
2742template <
typename VectorType>
2746 const VectorType &src)
const
2751 if (eigenvalues_are_initialized ==
false)
2752 estimate_eigenvalues(src);
2756 for (
unsigned int i = 0; i < data.n_iterations; ++i)
2757 internal::PreconditionRelaxation::step_operations(*A,
2758 *data.preconditioner,
2768template <
typename MatrixType,
typename PreconditionerType>
2769template <
typename VectorType>
2773 const VectorType &src)
const
2778 if (eigenvalues_are_initialized ==
false)
2779 estimate_eigenvalues(src);
2783 for (
unsigned int i = 0; i < data.n_iterations; ++i)
2784 internal::PreconditionRelaxation::step_operations(
2785 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i,
true);
2788template <
typename MatrixType,
typename PreconditionerType>
2789template <
typename VectorType>
2793 const VectorType &src)
const
2798 if (eigenvalues_are_initialized ==
false)
2799 estimate_eigenvalues(src);
2803 for (
unsigned int i = 1; i <= data.n_iterations; ++i)
2804 internal::PreconditionRelaxation::step_operations(*A,
2805 *data.preconditioner,
2815template <
typename MatrixType,
typename PreconditionerType>
2816template <
typename VectorType>
2820 const VectorType &src)
const
2825 if (eigenvalues_are_initialized ==
false)
2826 estimate_eigenvalues(src);
2830 for (
unsigned int i = 1; i <= data.n_iterations; ++i)
2831 internal::PreconditionRelaxation::step_operations(
2832 *A, *data.preconditioner, dst, src, data.relaxation, tmp1, tmp2, i,
true);
2835template <
typename MatrixType,
typename PreconditionerType>
2836template <
typename VectorType>
2839 const VectorType &src)
const
2843 EigenvalueInformation info;
2845 if (data.relaxation == 0.0)
2849 solution_old.reinit(src);
2850 temp_vector1.reinit(src,
true);
2852 info = internal::estimate_eigenvalues<MatrixType>(
2853 data, A, solution_old, temp_vector1, data.n_iterations);
2855 const double alpha =
2856 (data.smoothing_range > 1. ?
2857 info.max_eigenvalue_estimate / data.smoothing_range :
2858 std::min(0.9 * info.max_eigenvalue_estimate,
2859 info.min_eigenvalue_estimate));
2862 ->data.relaxation = 2.0 / (alpha + info.max_eigenvalue_estimate);
2866 ->eigenvalues_are_initialized =
true;
2871template <
typename MatrixType,
typename PreconditionerType>
2875 return data.relaxation;
2881template <
typename MatrixType>
2884 const AdditionalData ¶meters_in)
2888 parameters_in.relaxation != 0.0,
2890 "Relaxation cannot automatically be determined by PreconditionJacobi."));
2892 AdditionalData parameters;
2893 parameters.relaxation = 1.0;
2894 parameters.n_iterations = parameters_in.n_iterations;
2895 parameters.preconditioner =
2896 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2898 this->BaseClass::initialize(A, parameters);
2903template <
typename MatrixType>
2906 const AdditionalData ¶meters_in)
2910 parameters_in.relaxation != 0.0,
2912 "Relaxation cannot automatically be determined by PreconditionSOR."));
2914 AdditionalData parameters;
2915 parameters.relaxation = 1.0;
2916 parameters.n_iterations = parameters_in.n_iterations;
2917 parameters.preconditioner =
2918 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2920 this->BaseClass::initialize(A, parameters);
2925template <
typename MatrixType>
2928 const AdditionalData ¶meters_in)
2932 parameters_in.relaxation != 0.0,
2934 "Relaxation cannot automatically be determined by PreconditionSSOR."));
2936 AdditionalData parameters;
2937 parameters.relaxation = 1.0;
2938 parameters.n_iterations = parameters_in.n_iterations;
2939 parameters.preconditioner =
2940 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2942 this->BaseClass::initialize(A, parameters);
2949template <
typename MatrixType>
2952 const MatrixType &A,
2953 const std::vector<size_type> &p,
2954 const std::vector<size_type> &ip,
2955 const typename BaseClass::AdditionalData ¶meters_in)
2959 parameters_in.relaxation != 0.0,
2961 "Relaxation cannot automatically be determined by PreconditionPSOR."));
2963 typename BaseClass::AdditionalData parameters;
2964 parameters.relaxation = 1.0;
2965 parameters.n_iterations = parameters_in.n_iterations;
2966 parameters.preconditioner =
2967 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2969 this->BaseClass::initialize(A, parameters);
2973template <
typename MatrixType>
2976 const AdditionalData &additional_data)
2979 additional_data.permutation,
2980 additional_data.inverse_permutation,
2981 additional_data.parameters);
2984template <
typename MatrixType>
2986 const std::vector<size_type> &permutation,
2987 const std::vector<size_type> &inverse_permutation,
2990 : permutation(permutation)
2991 , inverse_permutation(inverse_permutation)
2992 , parameters(parameters)
2999template <
typename MatrixType,
typename VectorType>
3001 const MatrixType &M,
3002 const function_ptr method)
3004 , precondition(method)
3009template <
typename MatrixType,
typename VectorType>
3013 const VectorType &src)
const
3015 (
matrix.*precondition)(dst, src);
3023 template <
typename PreconditionerType>
3026 const double smoothing_range,
3027 const unsigned int eig_cg_n_iterations,
3028 const double eig_cg_residual,
3029 const double max_eigenvalue,
3031 : smoothing_range(smoothing_range)
3032 , eig_cg_n_iterations(eig_cg_n_iterations)
3033 , eig_cg_residual(eig_cg_residual)
3034 , max_eigenvalue(max_eigenvalue)
3035 , eigenvalue_algorithm(eigenvalue_algorithm)
3040 template <
typename PreconditionerType>
3041 inline EigenvalueAlgorithmAdditionalData<PreconditionerType> &
3042 EigenvalueAlgorithmAdditionalData<PreconditionerType>::operator=(
3043 const EigenvalueAlgorithmAdditionalData &other_data)
3045 smoothing_range = other_data.smoothing_range;
3046 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3047 eig_cg_residual = other_data.eig_cg_residual;
3048 max_eigenvalue = other_data.max_eigenvalue;
3049 preconditioner = other_data.preconditioner;
3050 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3051 constraints.copy_from(other_data.constraints);
3057template <
typename MatrixType,
typename PreconditionerType>
3060 const unsigned int n_iterations,
3061 const double smoothing_range,
3062 const unsigned int eig_cg_n_iterations,
3063 const double eig_cg_residual,
3064 const double max_eigenvalue,
3065 const EigenvalueAlgorithm eigenvalue_algorithm)
3066 :
internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3068 eig_cg_n_iterations,
3071 eigenvalue_algorithm)
3072 , relaxation(relaxation)
3073 , n_iterations(n_iterations)
3082 namespace PreconditionChebyshevImplementation
3090 template <
typename VectorType,
typename PreconditionerType>
3092 vector_updates(
const VectorType &rhs,
3093 const PreconditionerType &preconditioner,
3094 const unsigned int iteration_index,
3095 const double factor1,
3096 const double factor2,
3097 VectorType &solution_old,
3098 VectorType &temp_vector1,
3099 VectorType &temp_vector2,
3100 VectorType &solution)
3102 if (iteration_index == 0)
3104 solution.equ(factor2, rhs);
3105 preconditioner.vmult(solution_old, solution);
3107 else if (iteration_index == 1)
3110 temp_vector1.sadd(-1.0, 1.0, rhs);
3111 preconditioner.vmult(solution_old, temp_vector1);
3114 solution_old.sadd(factor2, 1 + factor1, solution);
3119 temp_vector1.sadd(-1.0, 1.0, rhs);
3120 preconditioner.vmult(temp_vector2, temp_vector1);
3123 solution_old.sadd(-factor1, factor2, temp_vector2);
3124 solution_old.add(1 + factor1, solution);
3127 solution.swap(solution_old);
3133 typename PreconditionerType,
3135 !has_vmult_with_std_functions_for_precondition<
3142 const PreconditionerType &preconditioner,
3143 const unsigned int iteration_index,
3144 const double factor1_,
3145 const double factor2_,
3154 const Number factor1 = factor1_;
3155 const Number factor1_plus_1 = 1. + factor1_;
3156 const Number factor2 = factor2_;
3158 if (iteration_index == 0)
3160 const auto solution_old_ptr = solution_old.
begin();
3163 preconditioner.vmult(solution_old, rhs);
3168 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
3170 else if (iteration_index == 1)
3172 const auto solution_ptr = solution.
begin();
3173 const auto solution_old_ptr = solution_old.
begin();
3176 temp_vector1.
sadd(-1.0, 1.0, rhs);
3178 preconditioner.vmult(solution_old, temp_vector1);
3183 solution_old_ptr[i] =
3184 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
3188 const auto solution_ptr = solution.
begin();
3189 const auto solution_old_ptr = solution_old.
begin();
3190 const auto temp_vector2_ptr = temp_vector2.
begin();
3193 temp_vector1.
sadd(-1.0, 1.0, rhs);
3195 preconditioner.vmult(temp_vector2, temp_vector1);
3200 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3201 factor1 * solution_old_ptr[i] +
3202 temp_vector2_ptr[i] * factor2;
3205 solution.
swap(solution_old);
3210 typename PreconditionerType,
3212 has_vmult_with_std_functions_for_precondition<
3219 const PreconditionerType &preconditioner,
3220 const unsigned int iteration_index,
3221 const double factor1_,
3222 const double factor2_,
3231 const Number factor1 = factor1_;
3232 const Number factor1_plus_1 = 1. + factor1_;
3233 const Number factor2 = factor2_;
3235 const auto rhs_ptr = rhs.
begin();
3236 const auto temp_vector1_ptr = temp_vector1.
begin();
3237 const auto temp_vector2_ptr = temp_vector2.
begin();
3238 const auto solution_ptr = solution.
begin();
3239 const auto solution_old_ptr = solution_old.
begin();
3241 if (iteration_index == 0)
3243 preconditioner.vmult(
3246 [&](
const auto start_range,
const auto end_range) {
3247 if (end_range > start_range)
3248 std::memset(solution.
begin() + start_range,
3250 sizeof(Number) * (end_range - start_range));
3252 [&](
const auto begin,
const auto end) {
3254 for (std::size_t i = begin; i <
end; ++i)
3255 solution_ptr[i] *= factor2;
3260 preconditioner.vmult(
3263 [&](
const auto begin,
const auto end) {
3265 std::memset(temp_vector2.
begin() + begin,
3267 sizeof(Number) * (end - begin));
3270 for (std::size_t i = begin; i <
end; ++i)
3271 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
3273 [&](
const auto begin,
const auto end) {
3274 if (iteration_index == 1)
3277 for (std::size_t i = begin; i <
end; ++i)
3278 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
3279 factor2 * temp_vector2_ptr[i];
3284 for (std::size_t i = begin; i <
end; ++i)
3285 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3286 factor1 * solution_old_ptr[i] +
3287 factor2 * temp_vector2_ptr[i];
3292 if (iteration_index > 0)
3294 solution_old.
swap(temp_vector2);
3295 solution_old.
swap(solution);
3302 template <
typename Number>
3303 struct VectorUpdater
3305 VectorUpdater(
const Number *rhs,
3306 const Number *matrix_diagonal_inverse,
3307 const unsigned int iteration_index,
3308 const Number factor1,
3309 const Number factor2,
3310 Number *solution_old,
3314 , matrix_diagonal_inverse(matrix_diagonal_inverse)
3315 , iteration_index(iteration_index)
3318 , solution_old(solution_old)
3319 , tmp_vector(tmp_vector)
3320 , solution(solution)
3324 apply_to_subrange(
const std::size_t begin,
const std::size_t end)
const
3330 const Number factor1 = this->factor1;
3331 const Number factor1_plus_1 = 1. + this->factor1;
3332 const Number factor2 = this->factor2;
3333 if (iteration_index == 0)
3336 for (std::size_t i = begin; i <
end; ++i)
3337 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
3339 else if (iteration_index == 1)
3343 for (std::size_t i = begin; i <
end; ++i)
3347 factor1_plus_1 * solution[i] +
3348 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3355 for (std::size_t i = begin; i <
end; ++i)
3362 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
3363 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
3368 const Number *matrix_diagonal_inverse;
3369 const unsigned int iteration_index;
3370 const Number factor1;
3371 const Number factor2;
3372 mutable Number *solution_old;
3373 mutable Number *tmp_vector;
3374 mutable Number *solution;
3377 template <
typename Number>
3380 VectorUpdatesRange(
const VectorUpdater<Number> &updater,
3381 const std::size_t size)
3385 VectorUpdatesRange::apply_to_subrange(0, size);
3393 ~VectorUpdatesRange()
override =
default;
3396 apply_to_subrange(
const std::size_t begin,
3397 const std::size_t end)
const override
3399 updater.apply_to_subrange(begin, end);
3402 const VectorUpdater<Number> &updater;
3406 template <
typename Number>
3409 const ::Vector<Number> &rhs,
3411 const unsigned int iteration_index,
3412 const double factor1,
3413 const double factor2,
3419 VectorUpdater<Number> upd(rhs.begin(),
3420 jacobi.get_vector().begin(),
3424 solution_old.
begin(),
3425 temp_vector1.
begin(),
3427 VectorUpdatesRange<Number>(upd, rhs.size());
3430 if (iteration_index == 0)
3437 solution.
swap(temp_vector1);
3438 solution_old.
swap(temp_vector1);
3443 template <
typename Number>
3447 const ::DiagonalMatrix<
3449 const unsigned int iteration_index,
3450 const double factor1,
3451 const double factor2,
3459 VectorUpdater<Number> upd(rhs.
begin(),
3460 jacobi.get_vector().begin(),
3464 solution_old.
begin(),
3465 temp_vector1.
begin(),
3470 if (iteration_index == 0)
3477 solution.
swap(temp_vector1);
3478 solution_old.
swap(temp_vector1);
3489 typename PreconditionerType,
3493 PreconditionerType> &&
3494 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3496 has_vmult_with_std_functions_for_precondition<
MatrixType,
3500 vmult_and_update(
const MatrixType &matrix,
3501 const PreconditionerType &preconditioner,
3502 const VectorType &rhs,
3503 const unsigned int iteration_index,
3504 const double factor1,
3505 const double factor2,
3506 VectorType &solution,
3507 VectorType &solution_old,
3508 VectorType &temp_vector1,
3509 VectorType &temp_vector2)
3511 if (iteration_index > 0)
3512 matrix.vmult(temp_vector1, solution);
3529 typename PreconditionerType,
3533 PreconditionerType> &&
3534 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3536 has_vmult_with_std_functions_for_precondition<
MatrixType,
3540 vmult_and_update(
const MatrixType &matrix,
3541 const PreconditionerType &preconditioner,
3542 const VectorType &rhs,
3543 const unsigned int iteration_index,
3544 const double factor1_,
3545 const double factor2_,
3546 VectorType &solution,
3547 VectorType &solution_old,
3548 VectorType &temp_vector1,
3549 VectorType &temp_vector2)
3551 using Number =
typename VectorType::value_type;
3553 const Number factor1 = factor1_;
3554 const Number factor1_plus_1 = 1. + factor1_;
3555 const Number factor2 = factor2_;
3557 if (iteration_index == 0)
3559 preconditioner.vmult(
3562 [&](
const unsigned int start_range,
const unsigned int end_range) {
3564 if (end_range > start_range)
3565 std::memset(solution.begin() + start_range,
3567 sizeof(Number) * (end_range - start_range));
3569 [&](
const unsigned int start_range,
const unsigned int end_range) {
3570 const auto solution_ptr = solution.begin();
3573 for (std::size_t i = start_range; i < end_range; ++i)
3574 solution_ptr[i] *= factor2;
3579 temp_vector1.reinit(rhs,
true);
3580 temp_vector2.reinit(rhs,
true);
3586 [&](
const unsigned int start_range,
const unsigned int end_range) {
3589 if (end_range > start_range)
3590 std::memset(temp_vector1.begin() + start_range,
3592 sizeof(Number) * (end_range - start_range));
3594 [&](
const unsigned int start_range,
const unsigned int end_range) {
3595 const auto rhs_ptr = rhs.begin();
3596 const auto tmp_ptr = temp_vector1.begin();
3599 for (std::size_t i = start_range; i < end_range; ++i)
3600 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3604 preconditioner.vmult(
3607 [&](
const unsigned int start_range,
const unsigned int end_range) {
3610 if (end_range > start_range)
3611 std::memset(temp_vector2.begin() + start_range,
3613 sizeof(Number) * (end_range - start_range));
3615 [&](
const unsigned int start_range,
const unsigned int end_range) {
3616 const auto solution_ptr = solution.begin();
3617 const auto solution_old_ptr = solution_old.begin();
3618 const auto tmp_ptr = temp_vector2.begin();
3620 if (iteration_index == 1)
3623 for (std::size_t i = start_range; i < end_range; ++i)
3625 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3630 for (std::size_t i = start_range; i < end_range; ++i)
3631 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3632 factor1 * solution_old_ptr[i] +
3633 factor2 * tmp_ptr[i];
3637 solution.swap(temp_vector2);
3638 solution_old.swap(temp_vector2);
3646 typename PreconditionerType,
3647 std::enable_if_t<has_vmult_with_std_functions<
MatrixType,
3649 PreconditionerType>,
3652 vmult_and_update(
const MatrixType &matrix,
3653 const PreconditionerType &preconditioner,
3654 const VectorType &rhs,
3655 const unsigned int iteration_index,
3656 const double factor1,
3657 const double factor2,
3658 VectorType &solution,
3659 VectorType &solution_old,
3660 VectorType &temp_vector1,
3663 using Number =
typename VectorType::value_type;
3664 VectorUpdater<Number> updater(rhs.begin(),
3665 preconditioner.get_vector().begin(),
3669 solution_old.begin(),
3670 temp_vector1.begin(),
3672 if (iteration_index > 0)
3676 [&](
const unsigned int start_range,
const unsigned int end_range) {
3679 if (end_range > start_range)
3680 std::memset(temp_vector1.begin() + start_range,
3682 sizeof(Number) * (end_range - start_range));
3684 [&](
const unsigned int start_range,
const unsigned int end_range) {
3685 if (end_range > start_range)
3686 updater.apply_to_subrange(start_range, end_range);
3689 updater.apply_to_subrange(0U, solution.locally_owned_size());
3692 if (iteration_index == 0)
3699 solution.swap(temp_vector1);
3700 solution_old.swap(temp_vector1);
3704 template <
typename MatrixType,
typename PreconditionerType>
3706 initialize_preconditioner(
3707 const MatrixType &matrix,
3708 std::shared_ptr<PreconditionerType> &preconditioner)
3711 (void)preconditioner;
3715 template <
typename MatrixType,
typename VectorType>
3717 initialize_preconditioner(
3718 const MatrixType &matrix,
3721 if (preconditioner.get() ==
nullptr || preconditioner->m() !=
matrix.m())
3723 if (preconditioner.get() ==
nullptr)
3725 std::make_shared<::DiagonalMatrix<VectorType>>();
3728 preconditioner->m() == 0,
3730 "Preconditioner appears to be initialized but not sized correctly"));
3733 if (preconditioner->m() !=
matrix.m())
3735 preconditioner->get_vector().reinit(
matrix.m());
3736 for (
typename VectorType::size_type i = 0; i <
matrix.m(); ++i)
3737 preconditioner->get_vector()(i) = 1. /
matrix.el(i, i);
3746template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3749 const double smoothing_range,
3750 const unsigned int eig_cg_n_iterations,
3751 const double eig_cg_residual,
3752 const double max_eigenvalue,
3753 const EigenvalueAlgorithm eigenvalue_algorithm,
3754 const PolynomialType polynomial_type)
3755 :
internal::EigenvalueAlgorithmAdditionalData<PreconditionerType>(
3757 eig_cg_n_iterations,
3760 eigenvalue_algorithm)
3762 , polynomial_type(polynomial_type)
3767template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3772 , eigenvalues_are_initialized(false)
3775 std::is_same_v<size_type, typename VectorType::size_type>,
3776 "PreconditionChebyshev and VectorType must have the same size_type.");
3781template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3784 const MatrixType &matrix,
3785 const AdditionalData &additional_data)
3788 data = additional_data;
3790 ExcMessage(
"The degree of the Chebyshev method must be positive."));
3791 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3792 matrix, data.preconditioner);
3793 eigenvalues_are_initialized =
false;
3798template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3802 eigenvalues_are_initialized =
false;
3803 theta = delta = 1.0;
3804 matrix_ptr =
nullptr;
3807 solution_old.reinit(empty_vector);
3808 temp_vector1.reinit(empty_vector);
3809 temp_vector2.reinit(empty_vector);
3811 data.preconditioner.reset();
3816template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3823 solution_old.reinit(src);
3824 temp_vector1.reinit(src,
true);
3826 auto info = internal::estimate_eigenvalues<MatrixType>(
3827 data, matrix_ptr, solution_old, temp_vector1, data.degree);
3829 const double alpha = (data.smoothing_range > 1. ?
3830 info.max_eigenvalue_estimate / data.smoothing_range :
3831 std::min(0.9 * info.max_eigenvalue_estimate,
3832 info.min_eigenvalue_estimate));
3841 const double actual_range = info.max_eigenvalue_estimate / alpha;
3842 const double sigma = (1. -
std::sqrt(1. / actual_range)) /
3844 const double eps = data.smoothing_range;
3849 1 +
static_cast<unsigned int>(
3854 info.degree = data.degree;
3859 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3860 (info.max_eigenvalue_estimate) :
3861 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3864 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3868 using NumberType =
typename VectorType::value_type;
3876 (std::is_same_v<VectorType,
3879 temp_vector2.
reinit(src, true);
3883 temp_vector2.reinit(empty_vector);
3888 ->eigenvalues_are_initialized =
true;
3895template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3898 VectorType &solution,
3899 const VectorType &rhs)
const
3901 std::lock_guard<std::mutex> lock(mutex);
3902 if (eigenvalues_are_initialized ==
false)
3903 estimate_eigenvalues(rhs);
3905 internal::PreconditionChebyshevImplementation::vmult_and_update(
3907 *data.preconditioner,
3911 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3912 (4. / (3. * delta)) :
3921 if (data.degree < 2 ||
std::abs(delta) < 1e-40)
3924 double rhok = delta / theta, sigma = theta / delta;
3925 for (
unsigned int k = 0; k < data.degree - 1; ++k)
3927 double factor1 = 0.0;
3928 double factor2 = 0.0;
3930 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3932 factor1 = (2 * k + 1.) / (2 * k + 5.);
3933 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3937 const double rhokp = 1. / (2. * sigma - rhok);
3938 factor1 = rhokp * rhok;
3939 factor2 = 2. * rhokp / delta;
3943 internal::PreconditionChebyshevImplementation::vmult_and_update(
3945 *data.preconditioner,
3959template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
3962 VectorType &solution,
3963 const VectorType &rhs)
const
3965 std::lock_guard<std::mutex> lock(mutex);
3966 if (eigenvalues_are_initialized ==
false)
3967 estimate_eigenvalues(rhs);
3969 internal::PreconditionChebyshevImplementation::vector_updates(
3971 *data.preconditioner,
3974 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3975 (4. / (3. * delta)) :
3982 if (data.degree < 2 ||
std::abs(delta) < 1e-40)
3985 double rhok = delta / theta, sigma = theta / delta;
3986 for (
unsigned int k = 0; k < data.degree - 1; ++k)
3988 double factor1 = 0.0;
3989 double factor2 = 0.0;
3991 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3993 factor1 = (2 * k + 1.) / (2 * k + 5.);
3994 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3998 const double rhokp = 1. / (2. * sigma - rhok);
3999 factor1 = rhokp * rhok;
4000 factor2 = 2. * rhokp / delta;
4004 matrix_ptr->Tvmult(temp_vector1, solution);
4005 internal::PreconditionChebyshevImplementation::vector_updates(
4007 *data.preconditioner,
4020template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4023 VectorType &solution,
4024 const VectorType &rhs)
const
4026 std::lock_guard<std::mutex> lock(mutex);
4027 if (eigenvalues_are_initialized ==
false)
4028 estimate_eigenvalues(rhs);
4030 internal::PreconditionChebyshevImplementation::vmult_and_update(
4032 *data.preconditioner,
4036 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4037 (4. / (3. * delta)) :
4044 if (data.degree < 2 ||
std::abs(delta) < 1e-40)
4047 double rhok = delta / theta, sigma = theta / delta;
4048 for (
unsigned int k = 0; k < data.degree - 1; ++k)
4050 double factor1 = 0.0;
4051 double factor2 = 0.0;
4053 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4055 factor1 = (2 * k + 1.) / (2 * k + 5.);
4056 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4060 const double rhokp = 1. / (2. * sigma - rhok);
4061 factor1 = rhokp * rhok;
4062 factor2 = 2. * rhokp / delta;
4066 internal::PreconditionChebyshevImplementation::vmult_and_update(
4068 *data.preconditioner,
4082template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4085 VectorType &solution,
4086 const VectorType &rhs)
const
4088 std::lock_guard<std::mutex> lock(mutex);
4089 if (eigenvalues_are_initialized ==
false)
4090 estimate_eigenvalues(rhs);
4092 matrix_ptr->Tvmult(temp_vector1, solution);
4093 internal::PreconditionChebyshevImplementation::vector_updates(
4095 *data.preconditioner,
4098 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
4099 (4. / (3. * delta)) :
4106 if (data.degree < 2 ||
std::abs(delta) < 1e-40)
4109 double rhok = delta / theta, sigma = theta / delta;
4110 for (
unsigned int k = 0; k < data.degree - 1; ++k)
4112 double factor1 = 0.0;
4113 double factor2 = 0.0;
4115 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
4117 factor1 = (2 * k + 1.) / (2 * k + 5.);
4118 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
4122 const double rhokp = 1. / (2. * sigma - rhok);
4123 factor1 = rhokp * rhok;
4124 factor2 = 2. * rhokp / delta;
4128 matrix_ptr->Tvmult(temp_vector1, solution);
4129 internal::PreconditionChebyshevImplementation::vector_updates(
4131 *data.preconditioner,
4144template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4147 PreconditionerType>::size_type
4151 return matrix_ptr->m();
4156template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
4159 PreconditionerType>::size_type
4163 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()
@ 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)
static const 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)