15#ifndef dealii_solver_gmres_h
16#define dealii_solver_gmres_h
48 template <
typename,
typename>
64 namespace SolverGMRESImplementation
73 template <
typename VectorType>
121 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
144 const unsigned int max_basis_size,
145 const bool force_reorthogonalization);
169 template <
typename VectorType>
172 const unsigned int n,
174 const unsigned int accumulated_iterations = 0,
175 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal =
176 boost::signals2::signal<
void(
int)>());
276 std::vector<std::pair<double, double>> &rotations,
349template <
typename VectorType = Vector<
double>>
376 delayed_classical_gram_schmidt);
457 template <
typename MatrixType,
typename PreconditionerType>
462 const PreconditionerType &preconditioner);
470 boost::signals2::connection
472 const bool every_iteration =
false);
480 boost::signals2::connection
482 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
483 const bool every_iteration =
false);
492 boost::signals2::connection
495 const bool every_iteration =
true);
503 boost::signals2::connection
514 boost::signals2::connection
520 <<
"The number of temporary vectors you gave (" << arg1
521 <<
") is too small. It should be at least 10 for "
522 <<
"any results, and much more for reasonable ones.");
546 boost::signals2::signal<void(
const std::vector<std::complex<double>> &)>
553 boost::signals2::signal<void(
const std::vector<std::complex<double>> &)>
573 boost::signals2::signal<void(
605 const unsigned int n,
606 const boost::signals2::signal<
610 const boost::signals2::signal<
void(
double)> &cond_signal);
641template <
typename VectorType = Vector<
double>>
690 template <
typename MatrixType,
typename PreconditionerType>
695 const PreconditionerType &preconditioner);
716template <
typename VectorType>
718 const unsigned int max_basis_size,
719 const bool right_preconditioning,
720 const bool use_default_residual,
721 const bool force_re_orthogonalization,
722 const bool batched_mode,
724 : max_n_tmp_vectors(0)
725 , max_basis_size(max_basis_size)
726 , right_preconditioning(right_preconditioning)
727 , use_default_residual(use_default_residual)
728 , force_re_orthogonalization(force_re_orthogonalization)
729 , batched_mode(batched_mode)
730 , orthogonalization_strategy(orthogonalization_strategy)
732 Assert(max_basis_size >= 1,
733 ExcMessage(
"SolverGMRES needs at least one vector in the "
739template <
typename VectorType>
742 const AdditionalData &data)
744 , additional_data(data)
750template <
typename VectorType>
752 const AdditionalData &data)
754 , additional_data(data)
762 namespace SolverGMRESImplementation
764 template <
typename VectorType>
773 template <
typename VectorType>
775 TmpVectors<VectorType>::operator[](
const unsigned int i)
const
785 template <
typename VectorType>
787 TmpVectors<VectorType>::operator()(
const unsigned int i,
788 const VectorType &temp)
791 if (data[i] ==
nullptr)
794 data[i]->reinit(temp,
true);
801 template <
typename VectorType>
803 TmpVectors<VectorType>::size()
const
805 return (data.size() > 0 ? data.size() - 1 : 0);
810 template <
typename VectorType,
typename Enable =
void>
811 struct is_dealii_compatible_distributed_vector;
813 template <
typename VectorType>
814 struct is_dealii_compatible_distributed_vector<
816 std::enable_if_t<!internal::is_block_vector<VectorType>>>
818 static constexpr bool value = std::is_same_v<
826 template <
typename VectorType>
827 struct is_dealii_compatible_distributed_vector<
829 std::enable_if_t<internal::is_block_vector<VectorType>>>
831 static constexpr bool value = std::is_same_v<
832 typename VectorType::BlockType,
839 template <
typename VectorType,
840 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
850 template <
typename VectorType,
851 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
856 return vector.n_blocks();
861 template <
typename VectorType,
862 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
865 block(VectorType &vector,
const unsigned int b)
874 template <
typename VectorType,
875 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType>
878 block(
const VectorType &vector,
const unsigned int b)
887 template <
typename VectorType,
888 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
890 typename VectorType::BlockType &
891 block(VectorType &vector,
const unsigned int b)
893 return vector.block(b);
898 template <
typename VectorType,
899 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
901 const typename VectorType::BlockType &
902 block(
const VectorType &vector,
const unsigned int b)
904 return vector.block(b);
909 template <
bool delayed_reorthogonalization,
912 !is_dealii_compatible_distributed_vector<VectorType>::value,
913 VectorType> * =
nullptr>
915 Tvmult_add(
const unsigned int n,
916 const VectorType &vv,
917 const TmpVectors<VectorType> &orthogonal_vectors,
920 for (
unsigned int i = 0; i < n; ++i)
922 h(i) += vv * orthogonal_vectors[i];
923 if (delayed_reorthogonalization)
924 h(n + i) += orthogonal_vectors[i] * orthogonal_vectors[n - 1];
926 if (delayed_reorthogonalization)
932 template <
bool delayed_reorthogonalization,
935 is_dealii_compatible_distributed_vector<VectorType>::value,
936 VectorType> * =
nullptr>
938 Tvmult_add(
const unsigned int n,
939 const VectorType &vv,
940 const TmpVectors<VectorType> &orthogonal_vectors,
943 for (
unsigned int b = 0;
b <
n_blocks(vv); ++
b)
950 static constexpr unsigned int n_lanes =
954 for (
unsigned int i = 0; i < n; ++i)
957 correct[delayed_reorthogonalization ? 129 : 1];
958 if (delayed_reorthogonalization)
959 for (
unsigned int i = 0; i < n + 1; ++i)
964 constexpr unsigned int inner_batch_size =
965 delayed_reorthogonalization ? 6 : 12;
967 for (; c < block(vv, b).locally_owned_size() / n_lanes /
969 ++c, j += n_lanes * inner_batch_size)
972 for (
unsigned int k = 0; k < inner_batch_size; ++k)
973 vvec[k].load(block(vv, b).begin() + j + k * n_lanes);
975 for (
unsigned int k = 0; k < inner_batch_size; ++k)
977 block(orthogonal_vectors[n - 1], b).begin() + j +
982 last_vector[0] * vvec[0];
984 last_vector[0] * last_vector[0];
986 for (
unsigned int k = 1; k < inner_batch_size; ++k)
988 local_sum_0 += last_vector[k] * vvec[k];
989 if (delayed_reorthogonalization)
991 local_sum_1 += last_vector[k] * last_vector[k];
992 local_sum_2 += vvec[k] * vvec[k];
995 hs[n - 1] += local_sum_0;
996 if (delayed_reorthogonalization)
998 correct[n - 1] += local_sum_1;
999 correct[n] += local_sum_2;
1003 for (
unsigned int i = 0; i < n - 1; ++i)
1009 temp.
load(block(orthogonal_vectors[i], b).begin() + j);
1012 delayed_reorthogonalization ? temp * last_vector[0] :
1014 for (
unsigned int k = 1; k < inner_batch_size; ++k)
1016 temp.
load(block(orthogonal_vectors[i], b).begin() +
1018 local_sum_0 += temp * vvec[k];
1019 if (delayed_reorthogonalization)
1020 local_sum_1 += temp * last_vector[k];
1022 hs[i] += local_sum_0;
1023 if (delayed_reorthogonalization)
1024 correct[i] += local_sum_1;
1028 c *= inner_batch_size;
1029 for (; c < block(vv, b).locally_owned_size() / n_lanes;
1033 vvec.
load(block(vv, b).begin() + j);
1034 last_vector.
load(block(orthogonal_vectors[n - 1], b).begin() +
1036 hs[n - 1] += last_vector * vvec;
1037 if (delayed_reorthogonalization)
1039 correct[n - 1] += last_vector * last_vector;
1040 correct[n] += vvec * vvec;
1043 for (
unsigned int i = 0; i < n - 1; ++i)
1046 temp.
load(block(orthogonal_vectors[i], b).begin() + j);
1047 hs[i] += temp * vvec;
1048 if (delayed_reorthogonalization)
1049 correct[i] += temp * last_vector;
1053 for (
unsigned int i = 0; i < n; ++i)
1055 h(i) += hs[i].
sum();
1056 if (delayed_reorthogonalization)
1057 h(i + n) += correct[i].
sum();
1059 if (delayed_reorthogonalization)
1060 h(n + n) += correct[n].
sum();
1065 for (; j < block(vv, b).locally_owned_size(); ++j)
1067 const double vvec = block(vv, b).local_element(j);
1068 const double last_vector =
1069 block(orthogonal_vectors[n - 1], b).local_element(j);
1070 h(n - 1) += last_vector * vvec;
1071 if (delayed_reorthogonalization)
1073 h(n + n - 1) += last_vector * last_vector;
1074 h(n + n) += vvec * vvec;
1076 for (
unsigned int i = 0; i < n - 1; ++i)
1079 block(orthogonal_vectors[i], b).local_element(j);
1080 h(i) += temp * vvec;
1081 if (delayed_reorthogonalization)
1082 h(n + i) += temp * last_vector;
1092 template <
bool delayed_reorthogonalization,
1093 typename VectorType,
1095 !is_dealii_compatible_distributed_vector<VectorType>::value,
1096 VectorType> * =
nullptr>
1098 subtract_and_norm(
const unsigned int n,
1099 const TmpVectors<VectorType> &orthogonal_vectors,
1105 VectorType &last_vector =
1106 const_cast<VectorType &
>(orthogonal_vectors[n - 1]);
1107 for (
unsigned int i = 0; i < n - 1; ++i)
1109 if (delayed_reorthogonalization && i + 2 < n)
1110 last_vector.add(-h(n + i), orthogonal_vectors[i]);
1111 vv.add(-h(i), orthogonal_vectors[i]);
1114 if (delayed_reorthogonalization)
1117 last_vector.sadd(1. / h(n + n - 1),
1118 -h(n + n - 2) / h(n + n - 1),
1119 orthogonal_vectors[n - 2]);
1122 const double scaling_factor_vv = h(n + n) > 0.0 ?
1123 1. / (h(n + n - 1) * h(n + n)) :
1124 1. / (h(n + n - 1) * h(n + n - 1));
1125 vv.sadd(scaling_factor_vv,
1126 -h(n - 1) * scaling_factor_vv,
1131 return std::numeric_limits<double>::signaling_NaN();
1135 vv.add_and_dot(-h(n - 1), orthogonal_vectors[n - 1], vv));
1140 template <
bool delayed_reorthogonalization,
1141 typename VectorType,
1143 is_dealii_compatible_distributed_vector<VectorType>::value,
1144 VectorType> * =
nullptr>
1146 subtract_and_norm(
const unsigned int n,
1147 const TmpVectors<VectorType> &orthogonal_vectors,
1153 double norm_vv_temp = 0.0;
1154 VectorType &last_vector =
1155 const_cast<VectorType &
>(orthogonal_vectors[n - 1]);
1156 const double inverse_norm_previous =
1157 delayed_reorthogonalization ? 1. / h(n + n - 1) : 0.;
1158 const double scaling_factor_vv =
1159 delayed_reorthogonalization ?
1160 (h(n + n) > 0.0 ? inverse_norm_previous / h(n + n) :
1161 inverse_norm_previous / h(n + n - 1)) :
1164 for (
unsigned int b = 0;
b <
n_blocks(vv); ++
b)
1168 constexpr unsigned int inner_batch_size =
1169 delayed_reorthogonalization ? 6 : 12;
1174 block(vv, b).locally_owned_size() / n_lanes / inner_batch_size;
1175 ++c, j += n_lanes * inner_batch_size)
1180 const double last_factor = h(n - 1);
1181 for (
unsigned int k = 0; k < inner_batch_size; ++k)
1183 temp[k].
load(block(vv, b).
begin() + j + k * n_lanes);
1184 last_vec[k].
load(block(last_vector, b).
begin() + j +
1186 if (!delayed_reorthogonalization)
1187 temp[k] -= last_factor * last_vec[k];
1190 for (
unsigned int i = 0; i < n - 1; ++i)
1192 const double factor = h(i);
1193 const double correction_factor =
1194 (delayed_reorthogonalization ? h(n + i) : 0.0);
1195 for (
unsigned int k = 0; k < inner_batch_size; ++k)
1198 vec.
load(block(orthogonal_vectors[i], b).
begin() + j +
1200 temp[k] -= factor * vec;
1201 if (delayed_reorthogonalization)
1202 last_vec[k] -= correction_factor * vec;
1206 if (delayed_reorthogonalization)
1207 for (
unsigned int k = 0; k < inner_batch_size; ++k)
1209 last_vec[k] = last_vec[k] * inverse_norm_previous;
1210 last_vec[k].
store(block(last_vector, b).
begin() + j +
1212 temp[k] -= last_factor * last_vec[k];
1213 temp[k] = temp[k] * scaling_factor_vv;
1214 temp[k].
store(block(vv, b).
begin() + j + k * n_lanes);
1217 for (
unsigned int k = 0; k < inner_batch_size; ++k)
1219 temp[k].
store(block(vv, b).
begin() + j + k * n_lanes);
1220 norm_vv_temp_vectorized += temp[k] * temp[k];
1224 c *= inner_batch_size;
1225 for (; c < block(vv, b).locally_owned_size() / n_lanes;
1229 temp.
load(block(vv, b).begin() + j);
1230 last_vec.
load(block(last_vector, b).begin() + j);
1231 if (!delayed_reorthogonalization)
1232 temp -= h(n - 1) * last_vec;
1234 for (
unsigned int i = 0; i < n - 1; ++i)
1237 vec.
load(block(orthogonal_vectors[i], b).
begin() + j);
1239 if (delayed_reorthogonalization)
1240 last_vec -= h(n + i) * vec;
1243 if (delayed_reorthogonalization)
1245 last_vec = last_vec * inverse_norm_previous;
1246 last_vec.
store(block(last_vector, b).
begin() + j);
1247 temp -= h(n - 1) * last_vec;
1248 temp = temp * scaling_factor_vv;
1254 norm_vv_temp_vectorized += temp * temp;
1258 if (!delayed_reorthogonalization)
1259 norm_vv_temp += norm_vv_temp_vectorized.
sum();
1261 for (; j < block(vv, b).locally_owned_size(); ++j)
1263 double temp = block(vv, b).local_element(j);
1264 double last_vec = block(last_vector, b).local_element(j);
1265 if (delayed_reorthogonalization)
1267 for (
unsigned int i = 0; i < n - 1; ++i)
1270 block(orthogonal_vectors[i], b).local_element(j);
1272 last_vec -= h(n + i) * vec;
1274 last_vec *= inverse_norm_previous;
1275 block(last_vector, b).local_element(j) = last_vec;
1276 temp -= h(n - 1) * last_vec;
1277 temp *= scaling_factor_vv;
1281 temp -= h(n - 1) * last_vec;
1282 for (
unsigned int i = 0; i < n - 1; ++i)
1284 h(i) * block(orthogonal_vectors[i], b).local_element(j);
1285 norm_vv_temp += temp * temp;
1287 block(vv, b).local_element(j) = temp;
1297 template <
typename VectorType,
1299 !is_dealii_compatible_distributed_vector<VectorType>::value,
1300 VectorType> * =
nullptr>
1303 const unsigned int n,
1305 const TmpVectors<VectorType> &tmp_vectors,
1306 const bool zero_out)
1309 p.equ(h(0), tmp_vectors[0]);
1311 p.add(h(0), tmp_vectors[0]);
1313 for (
unsigned int i = 1; i < n; ++i)
1314 p.add(h(i), tmp_vectors[i]);
1319 template <
typename VectorType,
1321 is_dealii_compatible_distributed_vector<VectorType>::value,
1322 VectorType> * =
nullptr>
1325 const unsigned int n,
1327 const TmpVectors<VectorType> &tmp_vectors,
1328 const bool zero_out)
1330 for (
unsigned int b = 0;
b <
n_blocks(p); ++
b)
1331 for (
unsigned int j = 0; j < block(p, b).locally_owned_size(); ++j)
1333 double temp = zero_out ? 0 : block(p, b).local_element(j);
1334 for (
unsigned int i = 0; i < n; ++i)
1335 temp += block(tmp_vectors[i], b).local_element(j) * h(i);
1336 block(p, b).local_element(j) = temp;
1343 ArnoldiProcess::initialize(
1345 const unsigned int basis_size,
1346 const bool force_reorthogonalization)
1348 this->orthogonalization_strategy = orthogonalization_strategy;
1349 this->do_reorthogonalization = force_reorthogonalization;
1351 hessenberg_matrix.reinit(basis_size + 1, basis_size);
1352 triangular_matrix.reinit(basis_size + 1, basis_size,
true);
1355 projected_rhs.reinit(basis_size + 1,
true);
1356 givens_rotations.reserve(basis_size);
1358 if (orthogonalization_strategy ==
1361 h.
reinit(2 * basis_size + 3);
1363 h.
reinit(basis_size + 1);
1368 template <
typename VectorType>
1370 ArnoldiProcess::orthonormalize_nth_vector(
1371 const unsigned int n,
1372 TmpVectors<VectorType> &orthogonal_vectors,
1373 const unsigned int accumulated_iterations,
1374 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal)
1379 VectorType &vv = orthogonal_vectors[n];
1381 double residual_estimate = std::numeric_limits<double>::signaling_NaN();
1384 givens_rotations.clear();
1385 residual_estimate = vv.l2_norm();
1386 if (residual_estimate != 0.)
1387 vv /= residual_estimate;
1388 projected_rhs(0) = residual_estimate;
1390 else if (orthogonalization_strategy ==
1400 const double previous_scaling = n > 0 ? h(n + n - 2) : 1.;
1403 h.reinit(n + n + 1);
1406 Tvmult_add<true>(n, vv, orthogonal_vectors, h);
1410 for (
unsigned int i = 0; i < n - 1; ++i)
1411 tmp += h(n + i) * h(n + i);
1412 const double alpha_j = h(n + n - 1) > tmp ?
1415 h(n + n - 1) = alpha_j;
1418 for (
unsigned int i = 0; i < n - 1; ++i)
1419 tmp += h(i) * h(n + i);
1420 h(n - 1) = (h(n - 1) - tmp) / alpha_j;
1425 for (
unsigned int i = 0; i < n - 1; ++i)
1426 hessenberg_matrix(i, n - 2) += h(n + i) * previous_scaling;
1427 hessenberg_matrix(n - 1, n - 2) = alpha_j * previous_scaling;
1429 for (
unsigned int i = 0; i < n; ++i)
1432 for (
unsigned int j = (i == 0 ? 0 : i - 1); j < n - 1; ++j)
1433 sum += hessenberg_matrix(i, j) * h(n + j);
1434 hessenberg_matrix(i, n - 1) = (h(i) -
sum) / alpha_j;
1440 for (
unsigned int i = 0; i < n - 1; ++i)
1442 sum += (2. - 1.) * h(n - 1) * h(n - 1);
1443 hessenberg_matrix(n, n - 1) =
1450 h(n + n) = hessenberg_matrix(n, n - 1);
1451 subtract_and_norm<true>(n, orthogonal_vectors, h, vv);
1455 residual_estimate = do_givens_rotation(
1456 true, n - 2, triangular_matrix, givens_rotations, projected_rhs);
1461 double norm_vv = 0.0;
1462 double norm_vv_start = 0;
1463 const bool consider_reorthogonalize =
1464 (do_reorthogonalization ==
false) && (n % 5 == 0);
1465 if (consider_reorthogonalize)
1466 norm_vv_start = vv.l2_norm();
1472 for (
unsigned int c = 0; c < 2; ++c)
1475 if (orthogonalization_strategy ==
1479 double htmp = vv * orthogonal_vectors[0];
1481 for (
unsigned int i = 1; i < n; ++i)
1483 htmp = vv.add_and_dot(-htmp,
1484 orthogonal_vectors[i - 1],
1485 orthogonal_vectors[i]);
1490 vv.add_and_dot(-htmp, orthogonal_vectors[n - 1], vv));
1492 else if (orthogonalization_strategy ==
1496 Tvmult_add<false>(n, vv, orthogonal_vectors, h);
1498 subtract_and_norm<false>(n, orthogonal_vectors, h, vv);
1516 if (consider_reorthogonalize)
1519 10. * norm_vv_start *
1521 typename VectorType::value_type>::epsilon()))
1526 do_reorthogonalization =
true;
1527 if (!reorthogonalize_signal.empty())
1528 reorthogonalize_signal(accumulated_iterations);
1532 if (do_reorthogonalization ==
false)
1536 for (
unsigned int i = 0; i < n; ++i)
1537 hessenberg_matrix(i, n - 1) = h(i);
1538 hessenberg_matrix(n, n - 1) = norm_vv;
1545 residual_estimate = do_givens_rotation(
1546 false, n - 1, triangular_matrix, givens_rotations, projected_rhs);
1549 return residual_estimate;
1555 ArnoldiProcess::do_givens_rotation(
1556 const bool delayed_reorthogonalization,
1559 std::vector<std::pair<double, double>> &rotations,
1567 if (delayed_reorthogonalization)
1572 matrix(0, col) = hessenberg_matrix(0, col);
1574 double H_next = hessenberg_matrix(0, col + 1);
1575 for (
int i = 0; i < col; ++i)
1577 const double c = rotations[i].first;
1578 const double s = rotations[i].second;
1579 const double Hi =
matrix(i, col);
1580 const double Hi1 = hessenberg_matrix(i + 1, col);
1581 H_next = -s * H_next + c * hessenberg_matrix(i + 1, col + 1);
1582 matrix(i, col) = c * Hi + s * Hi1;
1583 matrix(i + 1, col) = -s * Hi + c * Hi1;
1588 const double H_col1 = hessenberg_matrix(col + 1, col);
1589 const double H_col =
matrix(col, col);
1590 const double r = 1. /
std::sqrt(H_col * H_col + H_col1 * H_col1);
1591 rotations.emplace_back(H_col * r, H_col1 * r);
1593 rotations[col].first * H_col + rotations[col].second * H_col1;
1595 rhs(col + 1) = -rotations[col].second * rhs(col);
1596 rhs(col) *= rotations[col].first;
1599 -rotations[col].second * H_next +
1600 rotations[col].first * hessenberg_matrix(col + 1, col + 1);
1603 const double H_last = hessenberg_matrix(col + 2, col + 1);
1604 const double r = 1. /
std::sqrt(H_next * H_next + H_last * H_last);
1605 return std::abs(H_last * r * rhs(col + 1));
1611 matrix(0, col) = hessenberg_matrix(0, col);
1612 for (
int i = 0; i < col; ++i)
1614 const double c = rotations[i].first;
1615 const double s = rotations[i].second;
1616 const double Hi =
matrix(i, col);
1617 const double Hi1 = hessenberg_matrix(i + 1, col);
1618 matrix(i, col) = c * Hi + s * Hi1;
1619 matrix(i + 1, col) = -s * Hi + c * Hi1;
1622 const double Hi =
matrix(col, col);
1623 const double Hi1 = hessenberg_matrix(col + 1, col);
1624 const double r = 1. /
std::sqrt(Hi * Hi + Hi1 * Hi1);
1625 rotations.emplace_back(Hi * r, Hi1 * r);
1627 rotations[col].first * Hi + rotations[col].second * Hi1;
1629 rhs(col + 1) = -rotations[col].second * rhs(col);
1630 rhs(col) *= rotations[col].first;
1639 ArnoldiProcess::solve_projected_system(
1640 const bool orthogonalization_finished)
1646 unsigned int n = givens_rotations.
size();
1655 if (orthogonalization_strategy ==
1660 if (!orthogonalization_finished)
1662 tmp_triangular_matrix = triangular_matrix;
1663 tmp_rhs = projected_rhs;
1664 std::vector<std::pair<double, double>> tmp_givens_rotations(
1666 do_givens_rotation(
false,
1667 givens_rotations.size(),
1668 tmp_triangular_matrix,
1669 tmp_givens_rotations,
1671 matrix = &tmp_triangular_matrix;
1675 do_givens_rotation(
false,
1676 givens_rotations.size(),
1683 projected_solution.reinit(n);
1684 for (
int i = n - 1; i >= 0; --i)
1686 double s = (*rhs)(i);
1687 for (
unsigned int j = i + 1; j < n; ++j)
1688 s -= projected_solution(j) * (*matrix)(i, j);
1689 projected_solution(i) = s / (*matrix)(i, i);
1693 return projected_solution;
1699 ArnoldiProcess::get_hessenberg_matrix()
const
1701 return hessenberg_matrix;
1708 complex_less_pred(
const std::complex<double> &x,
1709 const std::complex<double> &y)
1711 return x.real() < y.real() ||
1712 (x.real() == y.real() && x.imag() < y.imag());
1719template <
typename VectorType>
1723 const unsigned int n,
1724 const boost::signals2::signal<
void(
const std::vector<std::complex<double>> &)>
1725 &eigenvalues_signal,
1728 const boost::signals2::signal<
void(
double)> &cond_signal)
1731 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1732 !cond_signal.empty()) &&
1736 for (
unsigned int i = 0; i < n; ++i)
1737 for (
unsigned int j = 0; j < n; ++j)
1738 mat(i, j) = H_orig(i, j);
1739 hessenberg_signal(H_orig);
1741 if (!eigenvalues_signal.empty())
1747 mat_eig.compute_eigenvalues();
1749 for (
unsigned int i = 0; i < mat_eig.n(); ++i)
1754 internal::SolverGMRESImplementation::complex_less_pred);
1759 if (!cond_signal.empty() && (mat.n() > 1))
1762 double condition_number =
1763 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1764 cond_signal(condition_number);
1771template <
typename VectorType>
1772template <
typename MatrixType,
typename PreconditionerType>
1776 const VectorType &b,
1777 const PreconditionerType &preconditioner)
1779 std::unique_ptr<LogStream::Prefix> prefix;
1781 prefix = std::make_unique<LogStream::Prefix>(
"GMRES");
1785 const unsigned int basis_size =
1792 basis_size + 2, this->memory);
1796 unsigned int accumulated_iterations = 0;
1798 const bool do_eigenvalues =
1800 (!condition_number_signal.empty() ||
1801 !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1802 !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1803 !all_hessenberg_signal.empty());
1806 double res = std::numeric_limits<double>::lowest();
1818 VectorType &
p = basis_vectors(basis_size + 1, x);
1824 if (!use_default_residual)
1842 VectorType &v = basis_vectors(0, x);
1851 if (left_precondition)
1853 if (accumulated_iterations == 0 && x.all_zero())
1854 preconditioner.vmult(v, b);
1859 preconditioner.vmult(v, p);
1864 if (accumulated_iterations == 0 && x.all_zero())
1873 const double norm_v =
1874 arnoldi_process.orthonormalize_nth_vector(0,
1876 accumulated_iterations);
1881 if (use_default_residual)
1885 iteration_state = solver_control.
check(accumulated_iterations, res);
1888 this->iteration_status(accumulated_iterations, res, x);
1895 deallog <<
"default_res=" << norm_v << std::endl;
1897 if (left_precondition)
1900 r->sadd(-1., 1., b);
1903 preconditioner.vmult(*r, v);
1907 iteration_state = solver_control.
check(accumulated_iterations, res);
1910 this->iteration_status(accumulated_iterations, res, x);
1918 unsigned int inner_iteration = 0;
1919 for (; (inner_iteration < basis_size &&
1923 ++accumulated_iterations;
1925 VectorType &vv = basis_vectors(inner_iteration + 1, x);
1927 if (left_precondition)
1929 A.vmult(p, basis_vectors[inner_iteration]);
1930 preconditioner.vmult(vv, p);
1934 preconditioner.vmult(p, basis_vectors[inner_iteration]);
1939 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1,
1941 accumulated_iterations,
1942 re_orthogonalize_signal);
1944 if (use_default_residual)
1948 solver_control.
check(accumulated_iterations, res);
1951 this->iteration_status(accumulated_iterations, res, x);
1956 deallog <<
"default_res=" << res << std::endl;
1960 arnoldi_process.solve_projected_system(
false);
1962 if (left_precondition)
1963 for (
unsigned int i = 0; i < inner_iteration + 1; ++i)
1964 x_->add(projected_solution(i), basis_vectors[i]);
1968 for (
unsigned int i = 0; i < inner_iteration + 1; ++i)
1969 p.add(projected_solution(i), basis_vectors[i]);
1970 preconditioner.vmult(*r, p);
1974 r->sadd(-1., 1., b);
1977 if (left_precondition)
1981 this->iteration_status(accumulated_iterations, res, x);
1985 preconditioner.vmult(*x_, *r);
1986 res = x_->l2_norm();
1990 solver_control.
check(accumulated_iterations, res);
1993 this->iteration_status(accumulated_iterations, res, x);
2001 arnoldi_process.solve_projected_system(
true);
2004 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
2006 all_eigenvalues_signal,
2007 all_hessenberg_signal,
2008 condition_number_signal);
2010 if (left_precondition)
2011 ::internal::SolverGMRESImplementation::add(
2012 x, inner_iteration, projected_solution, basis_vectors,
false);
2015 ::internal::SolverGMRESImplementation::add(
2016 p, inner_iteration, projected_solution, basis_vectors,
true);
2017 preconditioner.vmult(v, p);
2025 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
2029 condition_number_signal);
2031 if (!additional_data.
batched_mode && !krylov_space_signal.empty())
2032 krylov_space_signal(basis_vectors);
2047template <
typename VectorType>
2048boost::signals2::connection
2050 const std::function<
void(
double)> &slot,
2051 const bool every_iteration)
2053 if (every_iteration)
2055 return all_condition_numbers_signal.connect(slot);
2059 return condition_number_signal.connect(slot);
2065template <
typename VectorType>
2066boost::signals2::connection
2068 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
2069 const bool every_iteration)
2071 if (every_iteration)
2073 return all_eigenvalues_signal.connect(slot);
2077 return eigenvalues_signal.connect(slot);
2083template <
typename VectorType>
2084boost::signals2::connection
2087 const bool every_iteration)
2089 if (every_iteration)
2091 return all_hessenberg_signal.connect(slot);
2095 return hessenberg_signal.connect(slot);
2101template <
typename VectorType>
2102boost::signals2::connection
2104 const std::function<
void(
2107 return krylov_space_signal.connect(slot);
2112template <
typename VectorType>
2113boost::signals2::connection
2115 const std::function<
void(
int)> &slot)
2117 return re_orthogonalize_signal.connect(slot);
2122template <
typename VectorType>
2135template <
typename VectorType>
2138 const AdditionalData &data)
2140 , additional_data(data)
2145template <
typename VectorType>
2147 const AdditionalData &data)
2149 , additional_data(data)
2154template <
typename VectorType>
2155template <
typename MatrixType,
typename PreconditionerType>
2159 const VectorType &b,
2160 const PreconditionerType &preconditioner)
2170 basis_size + 1, this->memory);
2172 basis_size, this->memory);
2176 unsigned int accumulated_iterations = 0;
2184 double res = std::numeric_limits<double>::lowest();
2194 if (accumulated_iterations == 0 && x.all_zero())
2198 A.vmult(v(0, x), x);
2199 v[0].sadd(-1., 1., b);
2202 res = arnoldi_process.orthonormalize_nth_vector(0, v);
2203 iteration_state = this->iteration_status(accumulated_iterations, res, x);
2207 unsigned int inner_iteration = 0;
2208 for (; (inner_iteration < basis_size &&
2212 preconditioner.vmult(z(inner_iteration, x), v[inner_iteration]);
2213 A.vmult(v(inner_iteration + 1, x), z[inner_iteration]);
2216 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1, v);
2223 this->iteration_status(++accumulated_iterations, res, x);
2229 arnoldi_process.solve_projected_system(
true);
2230 ::internal::SolverGMRESImplementation::add(
2231 x, inner_iteration, projected_solution, z,
false);
virtual State check(const unsigned int step, const double check_value)
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
internal::SolverGMRESImplementation::ArnoldiProcess arnoldi_process
SolverFGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
SolverFGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
boost::signals2::signal< void(const FullMatrix< double > &)> hessenberg_signal
AdditionalData additional_data
boost::signals2::signal< void(double)> condition_number_signal
SolverGMRES(const SolverGMRES< VectorType > &)=delete
boost::signals2::signal< void(const std::vector< std::complex< double > > &)> eigenvalues_signal
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
boost::signals2::connection connect_krylov_space_slot(const std::function< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> &slot)
boost::signals2::signal< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> krylov_space_signal
boost::signals2::signal< void(const FullMatrix< double > &)> all_hessenberg_signal
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
SolverControl & solver_control
static void compute_eigs_and_cond(const FullMatrix< double > &H_orig, const unsigned int n, const boost::signals2::signal< void(const std::vector< std::complex< double > > &)> &eigenvalues_signal, const boost::signals2::signal< void(const FullMatrix< double > &)> &hessenberg_signal, const boost::signals2::signal< void(double)> &cond_signal)
internal::SolverGMRESImplementation::ArnoldiProcess arnoldi_process
SolverGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
boost::signals2::connection connect_hessenberg_slot(const std::function< void(const FullMatrix< double > &)> &slot, const bool every_iteration=true)
boost::signals2::signal< void(const std::vector< std::complex< double > > &)> all_eigenvalues_signal
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
virtual double criterion()
boost::signals2::signal< void(double)> all_condition_numbers_signal
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< std::complex< double > > &)> &slot, const bool every_iteration=false)
boost::signals2::signal< void(int)> re_orthogonalize_signal
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
static constexpr std::size_t size()
void store(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
void initialize(const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy, const unsigned int max_basis_size, const bool force_reorthogonalization)
double orthonormalize_nth_vector(const unsigned int n, TmpVectors< VectorType > &orthogonal_vectors, const unsigned int accumulated_iterations=0, const boost::signals2::signal< void(int)> &reorthogonalize_signal=boost::signals2::signal< void(int)>())
Vector< double > projected_rhs
FullMatrix< double > hessenberg_matrix
Vector< double > projected_solution
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
std::vector< std::pair< double, double > > givens_rotations
bool do_reorthogonalization
const Vector< double > & solve_projected_system(const bool orthogonalization_finished)
FullMatrix< double > triangular_matrix
double do_givens_rotation(const bool delayed_reorthogonalization, const int col, FullMatrix< double > &matrix, std::vector< std::pair< double, double > > &rotations, Vector< double > &rhs)
const FullMatrix< double > & get_hessenberg_matrix() const
VectorMemory< VectorType > & mem
std::vector< typename VectorMemory< VectorType >::Pointer > data
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
unsigned int size() const
VectorType & operator[](const unsigned int i) const
VectorType & operator()(const unsigned int i, const VectorType &temp)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcTooFewTmpVectors(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
@ matrix
Contents is actually a matrix.
OrthogonalizationStrategy
@ delayed_classical_gram_schmidt
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
T sum(const T &t, const MPI_Comm mpi_communicator)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > max(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 max_basis_size
AdditionalData(const unsigned int max_basis_size=30, const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy=LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt)
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
AdditionalData(const unsigned int max_basis_size=30, const bool right_preconditioning=false, const bool use_default_residual=true, const bool force_re_orthogonalization=false, const bool batched_mode=false, const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy=LinearAlgebra::OrthogonalizationStrategy::delayed_classical_gram_schmidt)
bool right_preconditioning
bool force_re_orthogonalization
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
bool use_default_residual
unsigned int max_basis_size
unsigned int max_n_tmp_vectors
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)