15#ifndef dealii_solver_gmres_h
16#define dealii_solver_gmres_h
34#include <boost/signals2/signal.hpp>
51 template <
typename,
typename>
67 namespace SolverGMRESImplementation
76 template <
typename VectorType>
124 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
138 template <
typename Number>
149 const unsigned int max_basis_size,
150 const bool force_reorthogonalization);
174 template <
typename VectorType>
177 const unsigned int n,
179 const unsigned int accumulated_iterations = 0,
180 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal =
181 boost::signals2::signal<
void(
int)>());
286 std::vector<std::pair<double, double>> &rotations,
413template <
typename VectorType = Vector<
double>>
434 const bool right_preconditioning =
false,
435 const bool use_default_residual =
true,
436 const bool force_re_orthogonalization =
false,
437 const bool batched_mode =
false,
439 orthogonalization_strategy =
441 delayed_classical_gram_schmidt);
522 template <
typename MatrixType,
typename PreconditionerType>
526 void solve(const MatrixType &A,
529 const PreconditionerType &preconditioner);
537 boost::signals2::connection
538 connect_condition_number_slot(const
std::function<
void(
double)> &slot,
539 const
bool every_iteration = false);
547 boost::signals2::connection
548 connect_eigenvalues_slot(
549 const
std::function<
void(const
std::vector<
std::complex<
double>> &)> &slot,
550 const
bool every_iteration = false);
559 boost::signals2::connection
560 connect_hessenberg_slot(
562 const
bool every_iteration = true);
570 boost::signals2::connection
571 connect_krylov_space_slot(
573 void(const
internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
581 boost::signals2::connection
582 connect_re_orthogonalization_slot(const
std::function<
void(
int)> &slot);
587 << "The number of temporary vectors you gave (" << arg1
588 << ") is too small. It should be at least 10 for "
589 << "any results, and much more for reasonable ones.");
601 boost::signals2::signal<
void(
double)> condition_number_signal;
607 boost::signals2::signal<
void(
double)> all_condition_numbers_signal;
613 boost::signals2::signal<
void(const
std::vector<
std::complex<
double>> &)>
620 boost::signals2::signal<
void(const
std::vector<
std::complex<
double>> &)>
621 all_eigenvalues_signal;
634 all_hessenberg_signal;
640 boost::signals2::signal<
void(
641 const
internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
648 boost::signals2::signal<
void(
int)> re_orthogonalize_signal;
670 compute_eigs_and_cond(
672 const
unsigned int n,
673 const
boost::signals2::signal<
674 void(const
std::vector<
std::complex<
double>> &)> &eigenvalues_signal,
677 const
boost::signals2::signal<
void(
double)> &cond_signal);
683 internal::SolverGMRESImplementation::ArnoldiProcess<
684 typename VectorType::value_type>
752template <typename VectorType =
Vector<
double>>
767 orthogonalization_strategy =
769 delayed_classical_gram_schmidt,
770 const
bool use_truncated_mpgmres_strategy = true)
771 : max_basis_size(max_basis_size)
772 , orthogonalization_strategy(orthogonalization_strategy)
773 , use_truncated_mpgmres_strategy(use_truncated_mpgmres_strategy)
815 template <
typename MatrixType,
typename... PreconditionerTypes>
819 void solve(const MatrixType &A,
822 const PreconditionerTypes &...preconditioners);
841 template <
typename MatrixType,
typename... PreconditionerTypes>
847 const PreconditionerTypes &...preconditioners);
860 typename VectorType::value_type>
891template <
typename VectorType = Vector<
double>>
905 const unsigned int max_basis_size = 30,
908 delayed_classical_gram_schmidt)
909 : max_basis_size(max_basis_size)
910 , orthogonalization_strategy(orthogonalization_strategy)
942 template <
typename MatrixType,
typename... PreconditionerTypes>
946 void solve(const MatrixType &A,
949 const PreconditionerTypes &...preconditioners);
959template <
typename VectorType>
962 const unsigned int max_basis_size,
963 const bool right_preconditioning,
964 const bool use_default_residual,
965 const bool force_re_orthogonalization,
966 const bool batched_mode,
968 : max_n_tmp_vectors(0)
969 , max_basis_size(max_basis_size)
970 , right_preconditioning(right_preconditioning)
971 , use_default_residual(use_default_residual)
972 , force_re_orthogonalization(force_re_orthogonalization)
973 , batched_mode(batched_mode)
974 , orthogonalization_strategy(orthogonalization_strategy)
976 Assert(max_basis_size >= 1,
977 ExcMessage(
"SolverGMRES needs at least one vector in the "
983template <
typename VectorType>
987 const AdditionalData &
data)
989 , additional_data(
data)
995template <
typename VectorType>
998 const AdditionalData &
data)
1000 , additional_data(
data)
1001 , solver_control(cn)
1008 namespace SolverGMRESImplementation
1010 template <
typename VectorType>
1019 template <
typename VectorType>
1021 TmpVectors<VectorType>::operator[](
const unsigned int i)
const
1031 template <
typename VectorType>
1033 TmpVectors<VectorType>::operator()(
const unsigned int i,
1034 const VectorType &temp)
1037 if (
data[i] ==
nullptr)
1040 data[i]->reinit(temp,
true);
1047 template <
typename VectorType>
1049 TmpVectors<VectorType>::size()
const
1051 return (
data.size() > 0 ?
data.size() - 1 : 0);
1056 template <
typename VectorType,
typename Enable =
void>
1057 struct is_dealii_compatible_vector;
1059 template <
typename VectorType>
1060 struct is_dealii_compatible_vector<
1062 std::enable_if_t<!internal::is_block_vector<VectorType>>>
1064 static constexpr bool value =
1069 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
1074 template <
typename VectorType>
1075 struct is_dealii_compatible_vector<
1077 std::enable_if_t<internal::is_block_vector<VectorType>>>
1079 static constexpr bool value =
1081 typename VectorType::BlockType,
1084 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
1090 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType>
1101 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
1106 return vector.n_blocks();
1112 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType>
1115 block(VectorType &vector,
const unsigned int b)
1124 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType>
1127 block(
const VectorType &vector,
const unsigned int b)
1136 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
1138 typename VectorType::BlockType &
1139 block(VectorType &vector,
const unsigned int b)
1141 return vector.block(b);
1147 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
1149 const typename VectorType::BlockType &
1150 block(
const VectorType &vector,
const unsigned int b)
1152 return vector.block(b);
1157 template <
bool delayed_reorthogonalization,
1159 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1162 Tvmult_add(
const unsigned int n,
1163 const VectorType &vv,
1164 const TmpVectors<VectorType> &orthogonal_vectors,
1166 std::vector<const typename VectorType::value_type *> &)
1168 for (
unsigned int i = 0; i < n; ++i)
1170 h(i) += vv * orthogonal_vectors[i];
1171 if (delayed_reorthogonalization)
1172 h(n + i) += orthogonal_vectors[i] * orthogonal_vectors[n - 1];
1174 if (delayed_reorthogonalization)
1175 h(n + n) += vv * vv;
1181 template <
bool delayed_reorthogonalization,
typename Number>
1185 const Number *current_vector,
1186 const std::vector<const Number *> &orthogonal_vectors,
1191 template <
bool delayed_reorthogonalization,
1193 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1197 const unsigned int n,
1198 const VectorType &vv,
1199 const TmpVectors<VectorType> &orthogonal_vectors,
1201 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1203 for (
unsigned int b = 0;
b <
n_blocks(vv); ++
b)
1205 vector_ptrs.resize(n);
1206 for (
unsigned int i = 0; i < n; ++i)
1207 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
1209 do_Tvmult_add<delayed_reorthogonalization>(n,
1210 block(vv, b).
end() -
1211 block(vv, b).
begin(),
1212 block(vv, b).
begin(),
1222 template <
bool delayed_reorthogonalization,
1224 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1227 subtract_and_norm(
const unsigned int n,
1228 const TmpVectors<VectorType> &orthogonal_vectors,
1231 std::vector<const typename VectorType::value_type *> &)
1236 const_cast<VectorType &
>(orthogonal_vectors[n - 1]);
1237 for (
unsigned int i = 0; i < n - 1; ++i)
1239 if (delayed_reorthogonalization && i + 2 < n)
1240 last_vector.add(-h(n + i), orthogonal_vectors[i]);
1241 vv.add(-h(i), orthogonal_vectors[i]);
1244 if (delayed_reorthogonalization)
1247 last_vector.sadd(1. / h(n + n - 1),
1248 -h(n + n - 2) / h(n + n - 1),
1249 orthogonal_vectors[n - 2]);
1252 const double scaling_factor_vv = h(n + n) > 0.0 ?
1253 1. / (h(n + n - 1) * h(n + n)) :
1254 1. / (h(n + n - 1) * h(n + n - 1));
1255 vv.sadd(scaling_factor_vv,
1256 -h(n - 1) * scaling_factor_vv,
1261 return std::numeric_limits<double>::signaling_NaN();
1265 vv.add_and_dot(-h(n - 1), orthogonal_vectors[n - 1], vv));
1271 template <
bool delayed_reorthogonalization,
typename Number>
1275 const std::vector<const Number *> &orthogonal_vectors,
1277 Number *current_vector);
1281 template <
bool delayed_reorthogonalization,
1283 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1287 const unsigned int n,
1288 const TmpVectors<VectorType> &orthogonal_vectors,
1291 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1293 double norm_vv_temp = 0.0;
1295 for (
unsigned int b = 0;
b <
n_blocks(vv); ++
b)
1297 vector_ptrs.resize(n);
1298 for (
unsigned int i = 0; i < n; ++i)
1299 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
1301 norm_vv_temp += do_subtract_and_norm<delayed_reorthogonalization>(
1303 block(vv, b).
end() - block(vv, b).
begin(),
1306 block(vv, b).
begin());
1316 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1320 const unsigned int n,
1322 const TmpVectors<VectorType> &tmp_vectors,
1323 const bool zero_out,
1324 std::vector<const typename VectorType::value_type *> &)
1327 p.equ(h(0), tmp_vectors[0]);
1329 p.add(h(0), tmp_vectors[0]);
1331 for (
unsigned int i = 1; i < n; ++i)
1332 p.add(h(i), tmp_vectors[i]);
1338 template <
typename Number>
1340 do_add(
const unsigned int n_vectors,
1342 const std::vector<const Number *> &tmp_vectors,
1344 const bool zero_out,
1350 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1354 const unsigned int n,
1356 const TmpVectors<VectorType> &tmp_vectors,
1357 const bool zero_out,
1358 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1360 for (
unsigned int b = 0;
b <
n_blocks(p); ++
b)
1362 vector_ptrs.resize(n);
1363 for (
unsigned int i = 0; i < n; ++i)
1364 vector_ptrs[i] = block(tmp_vectors[i], b).begin();
1366 block(p, b).
end() - block(p, b).
begin(),
1370 block(p, b).
begin());
1376 template <
typename Number>
1378 ArnoldiProcess<Number>::initialize(
1380 const unsigned int basis_size,
1381 const bool force_reorthogonalization)
1383 this->orthogonalization_strategy = orthogonalization_strategy;
1384 this->do_reorthogonalization = force_reorthogonalization;
1386 hessenberg_matrix.reinit(basis_size + 1, basis_size);
1387 triangular_matrix.reinit(basis_size + 1, basis_size,
true);
1390 projected_rhs.reinit(basis_size + 1,
true);
1391 givens_rotations.reserve(basis_size);
1393 if (orthogonalization_strategy ==
1396 h.
reinit(2 * basis_size + 3);
1398 h.
reinit(basis_size + 1);
1403 template <
typename Number>
1404 template <
typename VectorType>
1406 ArnoldiProcess<Number>::orthonormalize_nth_vector(
1407 const unsigned int n,
1408 TmpVectors<VectorType> &orthogonal_vectors,
1409 const unsigned int accumulated_iterations,
1410 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal)
1417 double residual_estimate = std::numeric_limits<double>::signaling_NaN();
1420 givens_rotations.clear();
1421 residual_estimate = vv.l2_norm();
1422 if (residual_estimate != 0.)
1423 vv /= residual_estimate;
1424 projected_rhs(0) = residual_estimate;
1426 else if (orthogonalization_strategy ==
1436 const double previous_scaling = n > 0 ? h(n + n - 2) : 1.;
1439 h.reinit(n + n + 1);
1442 Tvmult_add<true>(n, vv, orthogonal_vectors, h, vector_ptrs);
1446 for (
unsigned int i = 0; i < n - 1; ++i)
1447 tmp += h(n + i) * h(n + i);
1448 const double alpha_j = h(n + n - 1) > tmp ?
1451 h(n + n - 1) = alpha_j;
1454 for (
unsigned int i = 0; i < n - 1; ++i)
1455 tmp += h(i) * h(n + i);
1456 h(n - 1) = (h(n - 1) - tmp) / alpha_j;
1461 for (
unsigned int i = 0; i < n - 1; ++i)
1462 hessenberg_matrix(i, n - 2) += h(n + i) * previous_scaling;
1463 hessenberg_matrix(n - 1, n - 2) = alpha_j * previous_scaling;
1465 for (
unsigned int i = 0; i < n; ++i)
1468 for (
unsigned int j = (i == 0 ? 0 : i - 1); j < n - 1; ++j)
1469 sum += hessenberg_matrix(i, j) * h(n + j);
1470 hessenberg_matrix(i, n - 1) = (h(i) -
sum) / alpha_j;
1476 for (
unsigned int i = 0; i < n - 1; ++i)
1478 sum += (2. - 1.) * h(n - 1) * h(n - 1);
1479 hessenberg_matrix(n, n - 1) =
1486 h(n + n) = hessenberg_matrix(n, n - 1);
1487 subtract_and_norm<true>(n, orthogonal_vectors, h, vv, vector_ptrs);
1491 residual_estimate = do_givens_rotation(
1492 true, n - 2, triangular_matrix, givens_rotations, projected_rhs);
1497 double norm_vv = 0.0;
1498 double norm_vv_start = 0;
1499 const bool consider_reorthogonalize =
1500 (do_reorthogonalization ==
false) && (n % 5 == 0);
1501 if (consider_reorthogonalize)
1502 norm_vv_start = vv.l2_norm();
1508 for (
unsigned int c = 0; c < 2; ++c)
1511 if (orthogonalization_strategy ==
1515 double htmp = vv * orthogonal_vectors[0];
1517 for (
unsigned int i = 1; i < n; ++i)
1519 htmp = vv.add_and_dot(-htmp,
1520 orthogonal_vectors[i - 1],
1521 orthogonal_vectors[i]);
1526 vv.add_and_dot(-htmp, orthogonal_vectors[n - 1], vv));
1528 else if (orthogonalization_strategy ==
1532 Tvmult_add<false>(n, vv, orthogonal_vectors, h, vector_ptrs);
1533 norm_vv = subtract_and_norm<false>(
1534 n, orthogonal_vectors, h, vv, vector_ptrs);
1552 if (consider_reorthogonalize)
1555 10. * norm_vv_start *
1557 typename VectorType::value_type>::epsilon()))
1562 do_reorthogonalization =
true;
1563 if (!reorthogonalize_signal.empty())
1564 reorthogonalize_signal(accumulated_iterations);
1568 if (do_reorthogonalization ==
false)
1572 for (
unsigned int i = 0; i < n; ++i)
1573 hessenberg_matrix(i, n - 1) = h(i);
1574 hessenberg_matrix(n, n - 1) = norm_vv;
1581 residual_estimate = do_givens_rotation(
1582 false, n - 1, triangular_matrix, givens_rotations, projected_rhs);
1585 return residual_estimate;
1590 template <
typename Number>
1592 ArnoldiProcess<Number>::do_givens_rotation(
1593 const bool delayed_reorthogonalization,
1596 std::vector<std::pair<double, double>> &rotations,
1604 if (delayed_reorthogonalization)
1609 matrix(0, col) = hessenberg_matrix(0, col);
1611 double H_next = hessenberg_matrix(0, col + 1);
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 H_next = -s * H_next + c * hessenberg_matrix(i + 1, col + 1);
1619 matrix(i, col) = c * Hi + s * Hi1;
1620 matrix(i + 1, col) = -s * Hi + c * Hi1;
1625 const double H_col1 = hessenberg_matrix(col + 1, col);
1626 const double H_col =
matrix(col, col);
1627 const double r = 1. /
std::sqrt(H_col * H_col + H_col1 * H_col1);
1628 rotations.emplace_back(H_col * r, H_col1 * r);
1630 rotations[col].first * H_col + rotations[col].second * H_col1;
1632 rhs(col + 1) = -rotations[col].second * rhs(col);
1633 rhs(col) *= rotations[col].first;
1636 -rotations[col].second * H_next +
1637 rotations[col].first * hessenberg_matrix(col + 1, col + 1);
1640 const double H_last = hessenberg_matrix(col + 2, col + 1);
1641 const double r = 1. /
std::sqrt(H_next * H_next + H_last * H_last);
1642 return std::abs(H_last * r * rhs(col + 1));
1648 matrix(0, col) = hessenberg_matrix(0, col);
1649 for (
int i = 0; i < col; ++i)
1651 const double c = rotations[i].first;
1652 const double s = rotations[i].second;
1653 const double Hi =
matrix(i, col);
1654 const double Hi1 = hessenberg_matrix(i + 1, col);
1655 matrix(i, col) = c * Hi + s * Hi1;
1656 matrix(i + 1, col) = -s * Hi + c * Hi1;
1659 const double Hi =
matrix(col, col);
1660 const double Hi1 = hessenberg_matrix(col + 1, col);
1661 const double r = 1. /
std::sqrt(Hi * Hi + Hi1 * Hi1);
1662 rotations.emplace_back(Hi * r, Hi1 * r);
1664 rotations[col].first * Hi + rotations[col].second * Hi1;
1666 rhs(col + 1) = -rotations[col].second * rhs(col);
1667 rhs(col) *= rotations[col].first;
1675 template <
typename Number>
1677 ArnoldiProcess<Number>::solve_projected_system(
1678 const bool orthogonalization_finished)
1684 unsigned int n = givens_rotations.
size();
1694 if (orthogonalization_strategy ==
1699 if (!orthogonalization_finished)
1701 tmp_triangular_matrix = triangular_matrix;
1702 tmp_rhs = projected_rhs;
1703 std::vector<std::pair<double, double>> tmp_givens_rotations(
1705 do_givens_rotation(
false,
1706 givens_rotations.size(),
1707 tmp_triangular_matrix,
1708 tmp_givens_rotations,
1710 matrix = &tmp_triangular_matrix;
1714 do_givens_rotation(
false,
1715 givens_rotations.size(),
1722 projected_solution.reinit(n);
1723 for (
int i = n - 1; i >= 0; --i)
1725 double s = (*rhs)(i);
1726 for (
unsigned int j = i + 1; j < n; ++j)
1727 s -= projected_solution(j) * (*matrix)(i, j);
1728 projected_solution(i) = s / (*matrix)(i, i);
1732 return projected_solution;
1737 template <
typename Number>
1739 ArnoldiProcess<Number>::get_hessenberg_matrix()
const
1741 return hessenberg_matrix;
1748 complex_less_pred(
const std::complex<double> &x,
1749 const std::complex<double> &y)
1751 return x.real() < y.real() ||
1752 (x.real() == y.real() && x.imag() < y.imag());
1759template <
typename VectorType>
1763 const unsigned int n,
1764 const boost::signals2::signal<
void(
const std::vector<std::complex<double>> &)>
1765 &eigenvalues_signal,
1768 const boost::signals2::signal<
void(
double)> &cond_signal)
1771 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1772 !cond_signal.empty()) &&
1776 for (
unsigned int i = 0; i < n; ++i)
1777 for (
unsigned int j = 0; j < n; ++j)
1778 mat(i, j) = H_orig(i, j);
1779 hessenberg_signal(H_orig);
1781 if (!eigenvalues_signal.empty())
1787 mat_eig.compute_eigenvalues();
1789 for (
unsigned int i = 0; i < mat_eig.n(); ++i)
1794 internal::SolverGMRESImplementation::complex_less_pred);
1799 if (!cond_signal.empty() && (mat.n() > 1))
1802 double condition_number =
1803 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1804 cond_signal(condition_number);
1811template <
typename VectorType>
1813template <
typename MatrixType,
typename PreconditionerType>
1819 const VectorType &b,
1820 const PreconditionerType &preconditioner)
1822 std::unique_ptr<LogStream::Prefix> prefix;
1824 prefix = std::make_unique<LogStream::Prefix>(
"GMRES");
1828 const unsigned int basis_size =
1835 basis_size + 2, this->memory);
1839 unsigned int accumulated_iterations = 0;
1841 const bool do_eigenvalues =
1843 (!condition_number_signal.empty() ||
1844 !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1845 !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1846 !all_hessenberg_signal.empty());
1849 double res = std::numeric_limits<double>::lowest();
1861 VectorType &p = basis_vectors(basis_size + 1, x);
1867 if (!use_default_residual)
1894 if (left_precondition)
1896 if (accumulated_iterations == 0 && x.all_zero())
1897 preconditioner.vmult(v, b);
1902 preconditioner.vmult(v, p);
1907 if (accumulated_iterations == 0 && x.all_zero())
1916 const double norm_v = arnoldi_process.orthonormalize_nth_vector(
1917 0, basis_vectors, accumulated_iterations, re_orthogonalize_signal);
1922 if (use_default_residual)
1926 iteration_state = solver_control.
check(accumulated_iterations, res);
1929 this->iteration_status(accumulated_iterations, res, x);
1936 deallog <<
"default_res=" << norm_v << std::endl;
1938 if (left_precondition)
1941 r->sadd(-1., 1., b);
1944 preconditioner.vmult(*r, v);
1948 iteration_state = solver_control.
check(accumulated_iterations, res);
1951 this->iteration_status(accumulated_iterations, res, x);
1959 unsigned int inner_iteration = 0;
1960 for (; (inner_iteration < basis_size &&
1964 ++accumulated_iterations;
1966 VectorType &vv = basis_vectors(inner_iteration + 1, x);
1968 if (left_precondition)
1970 A.vmult(p, basis_vectors[inner_iteration]);
1971 preconditioner.vmult(vv, p);
1975 preconditioner.vmult(p, basis_vectors[inner_iteration]);
1980 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1,
1982 accumulated_iterations,
1983 re_orthogonalize_signal);
1985 if (use_default_residual)
1989 solver_control.
check(accumulated_iterations, res);
1992 this->iteration_status(accumulated_iterations, res, x);
1997 deallog <<
"default_res=" << res << std::endl;
2001 arnoldi_process.solve_projected_system(
false);
2003 if (left_precondition)
2004 for (
unsigned int i = 0; i < inner_iteration + 1; ++i)
2005 x_->add(projected_solution(i), basis_vectors[i]);
2009 for (
unsigned int i = 0; i < inner_iteration + 1; ++i)
2010 p.add(projected_solution(i), basis_vectors[i]);
2011 preconditioner.vmult(*r, p);
2015 r->sadd(-1., 1., b);
2018 if (left_precondition)
2022 this->iteration_status(accumulated_iterations, res, x);
2026 preconditioner.vmult(*x_, *r);
2027 res = x_->l2_norm();
2031 solver_control.
check(accumulated_iterations, res);
2034 this->iteration_status(accumulated_iterations, res, x);
2042 arnoldi_process.solve_projected_system(
true);
2045 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
2047 all_eigenvalues_signal,
2048 all_hessenberg_signal,
2049 condition_number_signal);
2051 if (left_precondition)
2052 ::internal::SolverGMRESImplementation::add(
2058 arnoldi_process.vector_ptrs);
2061 ::internal::SolverGMRESImplementation::add(
2067 arnoldi_process.vector_ptrs);
2068 preconditioner.vmult(v, p);
2076 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
2080 condition_number_signal);
2082 if (!additional_data.
batched_mode && !krylov_space_signal.empty())
2083 krylov_space_signal(basis_vectors);
2098template <
typename VectorType>
2100boost::signals2::connection
2102 const std::function<
void(
double)> &slot,
2103 const bool every_iteration)
2105 if (every_iteration)
2107 return all_condition_numbers_signal.connect(slot);
2111 return condition_number_signal.connect(slot);
2117template <
typename VectorType>
2120 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
2121 const bool every_iteration)
2123 if (every_iteration)
2125 return all_eigenvalues_signal.connect(slot);
2129 return eigenvalues_signal.connect(slot);
2135template <
typename VectorType>
2139 const bool every_iteration)
2141 if (every_iteration)
2143 return all_hessenberg_signal.connect(slot);
2147 return hessenberg_signal.connect(slot);
2153template <
typename VectorType>
2156 const std::function<
void(
2159 return krylov_space_signal.connect(slot);
2164template <
typename VectorType>
2166boost::signals2::connection
2168 const std::function<
void(
int)> &slot)
2170 return re_orthogonalize_signal.connect(slot);
2175template <
typename VectorType>
2191template <
typename VectorType>
2195 const AdditionalData &
data)
2197 , additional_data(
data)
2202template <
typename VectorType>
2205 const AdditionalData &
data)
2207 , additional_data(
data)
2212template <
typename VectorType>
2214template <
typename MatrixType,
typename... PreconditionerTypes>
2219 const MatrixType &A,
2221 const VectorType &b,
2222 const PreconditionerTypes &...preconditioners)
2226 if (additional_data.use_truncated_mpgmres_strategy)
2228 A, x, b, IndexingStrategy::truncated_mpgmres, preconditioners...);
2231 A, x, b, IndexingStrategy::full_mpgmres, preconditioners...);
2236template <
typename VectorType>
2238template <
typename MatrixType,
typename... PreconditionerTypes>
2240 const MatrixType &A,
2242 const VectorType &b,
2243 const IndexingStrategy &indexing_strategy,
2244 const PreconditionerTypes &...preconditioners)
2246 constexpr std::size_t n_preconditioners =
sizeof...(PreconditionerTypes);
2251 const auto apply_nth_preconditioner = [&](
unsigned int n,
2257 [[maybe_unused]]
bool preconditioner_called =
false;
2259 const auto call_matching_preconditioner = [&](
const auto &preconditioner) {
2263 preconditioner_called =
true;
2264 preconditioner.vmult(dst, src);
2269 (call_matching_preconditioner(preconditioners), ...);
2273 std::size_t current_index = 0;
2279 const auto preconditioner_vmult = [&](
auto &dst,
const auto &src) {
2281 if (n_preconditioners == 0)
2285 apply_nth_preconditioner(current_index, dst, src);
2286 current_index = (current_index + 1) % n_preconditioners;
2293 const auto previous_vector_index =
2294 [n_preconditioners, indexing_strategy](
unsigned int i) ->
unsigned int {
2297 if (n_preconditioners == 0)
2302 switch (indexing_strategy)
2304 case IndexingStrategy::fgmres:
2307 case IndexingStrategy::full_mpgmres:
2309 return i / n_preconditioners;
2310 case IndexingStrategy::truncated_mpgmres:
2312 return (1 + i >= n_preconditioners) ? (1 + i - n_preconditioners) : 0;
2325 basis_size + 1, this->memory);
2327 basis_size, this->memory);
2331 unsigned int accumulated_iterations = 0;
2339 double res = std::numeric_limits<double>::lowest();
2349 if (accumulated_iterations == 0 && x.all_zero())
2353 A.vmult(v(0, x), x);
2354 v[0].sadd(-1., 1., b);
2357 res = arnoldi_process.orthonormalize_nth_vector(0, v);
2358 iteration_state = this->iteration_status(accumulated_iterations, res, x);
2362 unsigned int inner_iteration = 0;
2363 for (; (inner_iteration < basis_size &&
2367 preconditioner_vmult(z(inner_iteration, x),
2368 v[previous_vector_index(inner_iteration)]);
2369 A.vmult(v(inner_iteration + 1, x), z[inner_iteration]);
2372 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1, v);
2379 this->iteration_status(++accumulated_iterations, res, x);
2385 arnoldi_process.solve_projected_system(
true);
2386 ::internal::SolverGMRESImplementation::add(
2392 arnoldi_process.vector_ptrs);
2408template <
typename VectorType>
2412 const AdditionalData &
data)
2417 data.max_basis_size,
2418 data.orthogonalization_strategy,
2424template <
typename VectorType>
2427 const AdditionalData &
data)
2431 data.max_basis_size,
2432 data.orthogonalization_strategy,
2438template <
typename VectorType>
2440template <
typename MatrixType,
typename... PreconditionerTypes>
2445 const MatrixType &A,
2447 const VectorType &b,
2448 const PreconditionerTypes &...preconditioners)
virtual State check(const unsigned int step, const double check_value)
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
SolverFGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
SolverFGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverGMRES(const SolverGMRES< VectorType > &)=delete
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::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
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)
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)
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::connection connect_eigenvalues_slot(const std::function< void(const std::vector< std::complex< double > > &)> &slot, const bool every_iteration=false)
SolverMPGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
internal::SolverGMRESImplementation::ArnoldiProcess< typename VectorType::value_type > arnoldi_process
SolverMPGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
void solve_internal(const MatrixType &A, VectorType &x, const VectorType &b, const IndexingStrategy &indexing_strategy, const PreconditionerTypes &...preconditioners)
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
const FullMatrix< double > & get_hessenberg_matrix() const
Vector< double > projected_rhs
bool do_reorthogonalization
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
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)>())
const Vector< double > & solve_projected_system(const bool orthogonalization_finished)
void initialize(const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy, const unsigned int max_basis_size, const bool force_reorthogonalization)
Vector< double > projected_solution
FullMatrix< double > hessenberg_matrix
std::vector< const Number * > vector_ptrs
FullMatrix< double > triangular_matrix
std::vector< std::pair< double, double > > givens_rotations
double do_givens_rotation(const bool delayed_reorthogonalization, const int col, FullMatrix< double > &matrix, std::vector< std::pair< double, double > > &rotations, Vector< double > &rhs)
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_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ASSERT_UNREACHABLE()
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)
std::vector< index_type > data
types::global_dof_index locally_owned_size
@ matrix
Contents is actually a matrix.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
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)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
void do_Tvmult_add(const unsigned int n_vectors, const std::size_t locally_owned_size, const Number *current_vector, const std::vector< const Number * > &orthogonal_vectors, Vector< double > &h)
double do_subtract_and_norm(const unsigned int n_vectors, const std::size_t locally_owned_size, const std::vector< const Number * > &orthogonal_vectors, const Vector< double > &h, Number *current_vector)
void do_add(const unsigned int n_vectors, const std::size_t locally_owned_size, const std::vector< const Number * > &tmp_vectors, const Vector< double > &h, const bool zero_out, Number *output)
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 > &)
AdditionalData(const unsigned int max_basis_size=30, const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy=LinearAlgebra::OrthogonalizationStrategy::delayed_classical_gram_schmidt)
unsigned int max_basis_size
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
AdditionalData(const unsigned int max_basis_size=30, const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy=LinearAlgebra::OrthogonalizationStrategy::delayed_classical_gram_schmidt, const bool use_truncated_mpgmres_strategy=true)
LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy
bool use_truncated_mpgmres_strategy
unsigned int max_basis_size
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)