15#ifndef dealii_solver_gmres_h
16#define dealii_solver_gmres_h
34#include <boost/signals2/signal.hpp>
50 template <
typename,
typename>
66 namespace SolverGMRESImplementation
75 template <
typename VectorType>
123 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
136 template <
typename Number>
147 const unsigned int max_basis_size,
148 const bool force_reorthogonalization);
172 template <
typename VectorType>
175 const unsigned int n,
177 const unsigned int accumulated_iterations = 0,
178 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal =
179 boost::signals2::signal<
void(
int)>());
284 std::vector<std::pair<double, double>> &rotations,
66 namespace SolverGMRESImplementation {
…}
357template <
typename VectorType = Vector<
double>>
378 const bool right_preconditioning =
false,
379 const bool use_default_residual =
true,
380 const bool force_re_orthogonalization =
false,
381 const bool batched_mode =
false,
383 orthogonalization_strategy =
385 delayed_classical_gram_schmidt);
466 template <
typename MatrixType,
typename PreconditionerType>
470 void solve(const MatrixType &A,
473 const PreconditionerType &preconditioner);
481 boost::signals2::connection
482 connect_condition_number_slot(const
std::function<
void(
double)> &slot,
483 const
bool every_iteration = false);
491 boost::signals2::connection
492 connect_eigenvalues_slot(
493 const
std::function<
void(const
std::vector<
std::complex<
double>> &)> &slot,
494 const
bool every_iteration = false);
503 boost::signals2::connection
504 connect_hessenberg_slot(
506 const
bool every_iteration = true);
514 boost::signals2::connection
515 connect_krylov_space_slot(
517 void(const
internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
525 boost::signals2::connection
526 connect_re_orthogonalization_slot(const
std::function<
void(
int)> &slot);
531 << "The number of temporary vectors you gave (" << arg1
532 << ") is too small. It should be at least 10 for "
533 << "any results, and much more for reasonable ones.");
545 boost::signals2::signal<
void(
double)> condition_number_signal;
551 boost::signals2::signal<
void(
double)> all_condition_numbers_signal;
557 boost::signals2::signal<
void(const
std::vector<
std::complex<
double>> &)>
564 boost::signals2::signal<
void(const
std::vector<
std::complex<
double>> &)>
565 all_eigenvalues_signal;
578 all_hessenberg_signal;
584 boost::signals2::signal<
void(
585 const
internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
592 boost::signals2::signal<
void(
int)> re_orthogonalize_signal;
614 compute_eigs_and_cond(
616 const
unsigned int n,
617 const
boost::signals2::signal<
618 void(const
std::vector<
std::complex<
double>> &)> &eigenvalues_signal,
621 const
boost::signals2::signal<
void(
double)> &cond_signal);
627 internal::SolverGMRESImplementation::ArnoldiProcess<
628 typename VectorType::value_type>
654template <typename VectorType =
Vector<
double>>
669 orthogonalization_strategy =
671 delayed_classical_gram_schmidt)
672 : max_basis_size(max_basis_size)
673 , orthogonalization_strategy(orthogonalization_strategy)
704 template <
typename MatrixType,
typename PreconditionerType>
708 void solve(const MatrixType &A,
711 const PreconditionerType &preconditioner);
723 internal::SolverGMRESImplementation::ArnoldiProcess<
724 typename VectorType::value_type>
734template <
typename VectorType>
737 const unsigned int max_basis_size,
738 const bool right_preconditioning,
739 const bool use_default_residual,
740 const bool force_re_orthogonalization,
741 const bool batched_mode,
743 : max_n_tmp_vectors(0)
744 , max_basis_size(max_basis_size)
745 , right_preconditioning(right_preconditioning)
746 , use_default_residual(use_default_residual)
747 , force_re_orthogonalization(force_re_orthogonalization)
748 , batched_mode(batched_mode)
749 , orthogonalization_strategy(orthogonalization_strategy)
751 Assert(max_basis_size >= 1,
752 ExcMessage(
"SolverGMRES needs at least one vector in the "
758template <
typename VectorType>
762 const AdditionalData &
data)
764 , additional_data(
data)
770template <
typename VectorType>
773 const AdditionalData &
data)
775 , additional_data(
data)
783 namespace SolverGMRESImplementation
785 template <
typename VectorType>
794 template <
typename VectorType>
796 TmpVectors<VectorType>::operator[](
const unsigned int i)
const
806 template <
typename VectorType>
808 TmpVectors<VectorType>::operator()(
const unsigned int i,
809 const VectorType &temp)
812 if (
data[i] ==
nullptr)
815 data[i]->reinit(temp,
true);
822 template <
typename VectorType>
824 TmpVectors<VectorType>::size()
const
826 return (
data.size() > 0 ?
data.size() - 1 : 0);
831 template <
typename VectorType,
typename Enable =
void>
832 struct is_dealii_compatible_vector;
834 template <
typename VectorType>
835 struct is_dealii_compatible_vector<
837 std::enable_if_t<!internal::is_block_vector<VectorType>>>
839 static constexpr bool value =
844 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
849 template <
typename VectorType>
850 struct is_dealii_compatible_vector<
852 std::enable_if_t<internal::is_block_vector<VectorType>>>
854 static constexpr bool value =
856 typename VectorType::BlockType,
859 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
865 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType>
876 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
881 return vector.n_blocks();
887 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType>
890 block(VectorType &vector,
const unsigned int b)
900 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType>
903 block(
const VectorType &vector,
const unsigned int b)
913 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
915 typename VectorType::BlockType &
916 block(VectorType &vector,
const unsigned int b)
918 return vector.block(b);
924 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
926 const typename VectorType::BlockType &
927 block(
const VectorType &vector,
const unsigned int b)
929 return vector.block(b);
934 template <
bool delayed_reorthogonalization,
936 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
939 Tvmult_add(
const unsigned int n,
940 const VectorType &vv,
941 const TmpVectors<VectorType> &orthogonal_vectors,
943 std::vector<const typename VectorType::value_type *> &)
945 for (
unsigned int i = 0; i < n; ++i)
947 h(i) += vv * orthogonal_vectors[i];
948 if (delayed_reorthogonalization)
949 h(n + i) += orthogonal_vectors[i] * orthogonal_vectors[n - 1];
951 if (delayed_reorthogonalization)
958 template <
bool delayed_reorthogonalization,
typename Number>
962 const Number *current_vector,
963 const std::vector<const Number *> &orthogonal_vectors,
968 template <
bool delayed_reorthogonalization,
970 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
974 const unsigned int n,
975 const VectorType &vv,
976 const TmpVectors<VectorType> &orthogonal_vectors,
978 std::vector<const typename VectorType::value_type *> &vector_ptrs)
980 for (
unsigned int b = 0;
b <
n_blocks(vv); ++
b)
982 vector_ptrs.resize(n);
983 for (
unsigned int i = 0; i < n; ++i)
984 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
986 do_Tvmult_add<delayed_reorthogonalization>(n,
988 block(vv, b).
begin(),
989 block(vv, b).
begin(),
999 template <
bool delayed_reorthogonalization,
1001 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1004 subtract_and_norm(
const unsigned int n,
1005 const TmpVectors<VectorType> &orthogonal_vectors,
1008 std::vector<const typename VectorType::value_type *> &)
1013 const_cast<VectorType &
>(orthogonal_vectors[n - 1]);
1014 for (
unsigned int i = 0; i < n - 1; ++i)
1016 if (delayed_reorthogonalization && i + 2 < n)
1017 last_vector.add(-h(n + i), orthogonal_vectors[i]);
1018 vv.add(-h(i), orthogonal_vectors[i]);
1021 if (delayed_reorthogonalization)
1024 last_vector.sadd(1. / h(n + n - 1),
1025 -h(n + n - 2) / h(n + n - 1),
1026 orthogonal_vectors[n - 2]);
1029 const double scaling_factor_vv = h(n + n) > 0.0 ?
1030 1. / (h(n + n - 1) * h(n + n)) :
1031 1. / (h(n + n - 1) * h(n + n - 1));
1032 vv.sadd(scaling_factor_vv,
1033 -h(n - 1) * scaling_factor_vv,
1038 return std::numeric_limits<double>::signaling_NaN();
1042 vv.add_and_dot(-h(n - 1), orthogonal_vectors[n - 1], vv));
1048 template <
bool delayed_reorthogonalization,
typename Number>
1052 const std::vector<const Number *> &orthogonal_vectors,
1054 Number *current_vector);
1058 template <
bool delayed_reorthogonalization,
1060 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1064 const unsigned int n,
1065 const TmpVectors<VectorType> &orthogonal_vectors,
1068 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1070 double norm_vv_temp = 0.0;
1072 for (
unsigned int b = 0;
b <
n_blocks(vv); ++
b)
1074 vector_ptrs.resize(n);
1075 for (
unsigned int i = 0; i < n; ++i)
1076 vector_ptrs[i] = block(orthogonal_vectors[i], b).begin();
1078 norm_vv_temp += do_subtract_and_norm<delayed_reorthogonalization>(
1080 block(vv, b).
end() - block(vv, b).
begin(),
1083 block(vv, b).
begin());
1093 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1097 const unsigned int n,
1099 const TmpVectors<VectorType> &tmp_vectors,
1100 const bool zero_out,
1101 std::vector<const typename VectorType::value_type *> &)
1104 p.equ(h(0), tmp_vectors[0]);
1106 p.add(h(0), tmp_vectors[0]);
1108 for (
unsigned int i = 1; i < n; ++i)
1109 p.add(h(i), tmp_vectors[i]);
1115 template <
typename Number>
1117 do_add(
const unsigned int n_vectors,
1119 const std::vector<const Number *> &tmp_vectors,
1121 const bool zero_out,
1127 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1131 const unsigned int n,
1133 const TmpVectors<VectorType> &tmp_vectors,
1134 const bool zero_out,
1135 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1137 for (
unsigned int b = 0;
b <
n_blocks(p); ++
b)
1139 vector_ptrs.resize(n);
1140 for (
unsigned int i = 0; i < n; ++i)
1141 vector_ptrs[i] = block(tmp_vectors[i], b).begin();
1143 block(p, b).
end() - block(p, b).
begin(),
1147 block(p, b).
begin());
1153 template <
typename Number>
1155 ArnoldiProcess<Number>::initialize(
1157 const unsigned int basis_size,
1158 const bool force_reorthogonalization)
1160 this->orthogonalization_strategy = orthogonalization_strategy;
1161 this->do_reorthogonalization = force_reorthogonalization;
1163 hessenberg_matrix.reinit(basis_size + 1, basis_size);
1164 triangular_matrix.reinit(basis_size + 1, basis_size,
true);
1167 projected_rhs.reinit(basis_size + 1,
true);
1168 givens_rotations.reserve(basis_size);
1170 if (orthogonalization_strategy ==
1173 h.
reinit(2 * basis_size + 3);
1175 h.
reinit(basis_size + 1);
1180 template <
typename Number>
1181 template <
typename VectorType>
1183 ArnoldiProcess<Number>::orthonormalize_nth_vector(
1184 const unsigned int n,
1185 TmpVectors<VectorType> &orthogonal_vectors,
1186 const unsigned int accumulated_iterations,
1187 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal)
1194 double residual_estimate = std::numeric_limits<double>::signaling_NaN();
1197 givens_rotations.clear();
1198 residual_estimate = vv.l2_norm();
1199 if (residual_estimate != 0.)
1200 vv /= residual_estimate;
1201 projected_rhs(0) = residual_estimate;
1203 else if (orthogonalization_strategy ==
1213 const double previous_scaling = n > 0 ? h(n + n - 2) : 1.;
1216 h.reinit(n + n + 1);
1219 Tvmult_add<true>(n, vv, orthogonal_vectors, h, vector_ptrs);
1223 for (
unsigned int i = 0; i < n - 1; ++i)
1224 tmp += h(n + i) * h(n + i);
1225 const double alpha_j = h(n + n - 1) > tmp ?
1228 h(n + n - 1) = alpha_j;
1231 for (
unsigned int i = 0; i < n - 1; ++i)
1232 tmp += h(i) * h(n + i);
1233 h(n - 1) = (h(n - 1) - tmp) / alpha_j;
1238 for (
unsigned int i = 0; i < n - 1; ++i)
1239 hessenberg_matrix(i, n - 2) += h(n + i) * previous_scaling;
1240 hessenberg_matrix(n - 1, n - 2) = alpha_j * previous_scaling;
1242 for (
unsigned int i = 0; i < n; ++i)
1245 for (
unsigned int j = (i == 0 ? 0 : i - 1); j < n - 1; ++j)
1246 sum += hessenberg_matrix(i, j) * h(n + j);
1247 hessenberg_matrix(i, n - 1) = (h(i) -
sum) / alpha_j;
1253 for (
unsigned int i = 0; i < n - 1; ++i)
1255 sum += (2. - 1.) * h(n - 1) * h(n - 1);
1256 hessenberg_matrix(n, n - 1) =
1263 h(n + n) = hessenberg_matrix(n, n - 1);
1264 subtract_and_norm<true>(n, orthogonal_vectors, h, vv, vector_ptrs);
1268 residual_estimate = do_givens_rotation(
1269 true, n - 2, triangular_matrix, givens_rotations, projected_rhs);
1274 double norm_vv = 0.0;
1275 double norm_vv_start = 0;
1276 const bool consider_reorthogonalize =
1277 (do_reorthogonalization ==
false) && (n % 5 == 0);
1278 if (consider_reorthogonalize)
1279 norm_vv_start = vv.l2_norm();
1285 for (
unsigned int c = 0; c < 2; ++c)
1288 if (orthogonalization_strategy ==
1292 double htmp = vv * orthogonal_vectors[0];
1294 for (
unsigned int i = 1; i < n; ++i)
1296 htmp = vv.add_and_dot(-htmp,
1297 orthogonal_vectors[i - 1],
1298 orthogonal_vectors[i]);
1303 vv.add_and_dot(-htmp, orthogonal_vectors[n - 1], vv));
1305 else if (orthogonalization_strategy ==
1309 Tvmult_add<false>(n, vv, orthogonal_vectors, h, vector_ptrs);
1310 norm_vv = subtract_and_norm<false>(
1311 n, orthogonal_vectors, h, vv, vector_ptrs);
1329 if (consider_reorthogonalize)
1332 10. * norm_vv_start *
1334 typename VectorType::value_type>::epsilon()))
1339 do_reorthogonalization =
true;
1340 if (!reorthogonalize_signal.empty())
1341 reorthogonalize_signal(accumulated_iterations);
1345 if (do_reorthogonalization ==
false)
1349 for (
unsigned int i = 0; i < n; ++i)
1350 hessenberg_matrix(i, n - 1) = h(i);
1351 hessenberg_matrix(n, n - 1) = norm_vv;
1358 residual_estimate = do_givens_rotation(
1359 false, n - 1, triangular_matrix, givens_rotations, projected_rhs);
1362 return residual_estimate;
1367 template <
typename Number>
1369 ArnoldiProcess<Number>::do_givens_rotation(
1370 const bool delayed_reorthogonalization,
1373 std::vector<std::pair<double, double>> &rotations,
1381 if (delayed_reorthogonalization)
1386 matrix(0, col) = hessenberg_matrix(0, col);
1388 double H_next = hessenberg_matrix(0, col + 1);
1389 for (
int i = 0; i < col; ++i)
1391 const double c = rotations[i].first;
1392 const double s = rotations[i].second;
1393 const double Hi =
matrix(i, col);
1394 const double Hi1 = hessenberg_matrix(i + 1, col);
1395 H_next = -s * H_next + c * hessenberg_matrix(i + 1, col + 1);
1396 matrix(i, col) = c * Hi + s * Hi1;
1397 matrix(i + 1, col) = -s * Hi + c * Hi1;
1402 const double H_col1 = hessenberg_matrix(col + 1, col);
1403 const double H_col =
matrix(col, col);
1404 const double r = 1. /
std::sqrt(H_col * H_col + H_col1 * H_col1);
1405 rotations.emplace_back(H_col * r, H_col1 * r);
1407 rotations[col].first * H_col + rotations[col].second * H_col1;
1409 rhs(col + 1) = -rotations[col].second * rhs(col);
1410 rhs(col) *= rotations[col].first;
1413 -rotations[col].second * H_next +
1414 rotations[col].first * hessenberg_matrix(col + 1, col + 1);
1417 const double H_last = hessenberg_matrix(col + 2, col + 1);
1418 const double r = 1. /
std::sqrt(H_next * H_next + H_last * H_last);
1419 return std::abs(H_last * r * rhs(col + 1));
1425 matrix(0, col) = hessenberg_matrix(0, col);
1426 for (
int i = 0; i < col; ++i)
1428 const double c = rotations[i].first;
1429 const double s = rotations[i].second;
1430 const double Hi =
matrix(i, col);
1431 const double Hi1 = hessenberg_matrix(i + 1, col);
1432 matrix(i, col) = c * Hi + s * Hi1;
1433 matrix(i + 1, col) = -s * Hi + c * Hi1;
1436 const double Hi =
matrix(col, col);
1437 const double Hi1 = hessenberg_matrix(col + 1, col);
1438 const double r = 1. /
std::sqrt(Hi * Hi + Hi1 * Hi1);
1439 rotations.emplace_back(Hi * r, Hi1 * r);
1441 rotations[col].first * Hi + rotations[col].second * Hi1;
1443 rhs(col + 1) = -rotations[col].second * rhs(col);
1444 rhs(col) *= rotations[col].first;
1452 template <
typename Number>
1454 ArnoldiProcess<Number>::solve_projected_system(
1455 const bool orthogonalization_finished)
1461 unsigned int n = givens_rotations.
size();
1471 if (orthogonalization_strategy ==
1476 if (!orthogonalization_finished)
1478 tmp_triangular_matrix = triangular_matrix;
1479 tmp_rhs = projected_rhs;
1480 std::vector<std::pair<double, double>> tmp_givens_rotations(
1482 do_givens_rotation(
false,
1483 givens_rotations.size(),
1484 tmp_triangular_matrix,
1485 tmp_givens_rotations,
1487 matrix = &tmp_triangular_matrix;
1491 do_givens_rotation(
false,
1492 givens_rotations.size(),
1499 projected_solution.reinit(n);
1500 for (
int i = n - 1; i >= 0; --i)
1502 double s = (*rhs)(i);
1503 for (
unsigned int j = i + 1; j < n; ++j)
1504 s -= projected_solution(j) * (*matrix)(i, j);
1505 projected_solution(i) = s / (*matrix)(i, i);
1509 return projected_solution;
1514 template <
typename Number>
1516 ArnoldiProcess<Number>::get_hessenberg_matrix()
const
1518 return hessenberg_matrix;
1525 complex_less_pred(
const std::complex<double> &x,
1526 const std::complex<double> &y)
1528 return x.real() < y.real() ||
1529 (x.real() == y.real() && x.imag() < y.imag());
1536template <
typename VectorType>
1540 const unsigned int n,
1541 const boost::signals2::signal<
void(
const std::vector<std::complex<double>> &)>
1542 &eigenvalues_signal,
1545 const boost::signals2::signal<
void(
double)> &cond_signal)
1548 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1549 !cond_signal.empty()) &&
1553 for (
unsigned int i = 0; i < n; ++i)
1554 for (
unsigned int j = 0; j < n; ++j)
1555 mat(i, j) = H_orig(i, j);
1556 hessenberg_signal(H_orig);
1558 if (!eigenvalues_signal.empty())
1564 mat_eig.compute_eigenvalues();
1566 for (
unsigned int i = 0; i < mat_eig.n(); ++i)
1571 internal::SolverGMRESImplementation::complex_less_pred);
1576 if (!cond_signal.empty() && (mat.n() > 1))
1579 double condition_number =
1580 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1581 cond_signal(condition_number);
1588template <
typename VectorType>
1590template <
typename MatrixType,
typename PreconditionerType>
1596 const VectorType &b,
1597 const PreconditionerType &preconditioner)
1599 std::unique_ptr<LogStream::Prefix> prefix;
1600 if (!additional_data.batched_mode)
1601 prefix = std::make_unique<LogStream::Prefix>(
"GMRES");
1605 const unsigned int basis_size =
1608 std::max(additional_data.max_n_tmp_vectors, 3u) - 2);
1612 basis_size + 2, this->memory);
1616 unsigned int accumulated_iterations = 0;
1618 const bool do_eigenvalues =
1619 !additional_data.batched_mode &&
1620 (!condition_number_signal.empty() ||
1621 !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1622 !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1623 !all_hessenberg_signal.empty());
1626 double res = std::numeric_limits<double>::lowest();
1630 const bool left_precondition = !additional_data.right_preconditioning;
1635 const bool use_default_residual = additional_data.use_default_residual;
1638 VectorType &p = basis_vectors(basis_size + 1, x);
1644 if (!use_default_residual)
1654 additional_data.force_re_orthogonalization);
1671 if (left_precondition)
1673 if (accumulated_iterations == 0 && x.all_zero())
1674 preconditioner.vmult(v, b);
1679 preconditioner.vmult(v, p);
1684 if (accumulated_iterations == 0 && x.all_zero())
1693 const double norm_v = arnoldi_process.orthonormalize_nth_vector(
1694 0, basis_vectors, accumulated_iterations, re_orthogonalize_signal);
1699 if (use_default_residual)
1702 if (additional_data.batched_mode)
1703 iteration_state = solver_control.
check(accumulated_iterations, res);
1706 this->iteration_status(accumulated_iterations, res, x);
1713 deallog <<
"default_res=" << norm_v << std::endl;
1715 if (left_precondition)
1718 r->sadd(-1., 1., b);
1721 preconditioner.vmult(*r, v);
1724 if (additional_data.batched_mode)
1725 iteration_state = solver_control.
check(accumulated_iterations, res);
1728 this->iteration_status(accumulated_iterations, res, x);
1736 unsigned int inner_iteration = 0;
1737 for (; (inner_iteration < basis_size &&
1741 ++accumulated_iterations;
1743 VectorType &vv = basis_vectors(inner_iteration + 1, x);
1745 if (left_precondition)
1747 A.vmult(p, basis_vectors[inner_iteration]);
1748 preconditioner.vmult(vv, p);
1752 preconditioner.vmult(p, basis_vectors[inner_iteration]);
1757 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1,
1759 accumulated_iterations,
1760 re_orthogonalize_signal);
1762 if (use_default_residual)
1764 if (additional_data.batched_mode)
1766 solver_control.
check(accumulated_iterations, res);
1769 this->iteration_status(accumulated_iterations, res, x);
1773 if (!additional_data.batched_mode)
1774 deallog <<
"default_res=" << res << std::endl;
1778 arnoldi_process.solve_projected_system(
false);
1780 if (left_precondition)
1781 for (
unsigned int i = 0; i < inner_iteration + 1; ++i)
1782 x_->add(projected_solution(i), basis_vectors[i]);
1786 for (
unsigned int i = 0; i < inner_iteration + 1; ++i)
1787 p.add(projected_solution(i), basis_vectors[i]);
1788 preconditioner.vmult(*r, p);
1792 r->sadd(-1., 1., b);
1795 if (left_precondition)
1799 this->iteration_status(accumulated_iterations, res, x);
1803 preconditioner.vmult(*x_, *r);
1804 res = x_->l2_norm();
1806 if (additional_data.batched_mode)
1808 solver_control.
check(accumulated_iterations, res);
1811 this->iteration_status(accumulated_iterations, res, x);
1819 arnoldi_process.solve_projected_system(
true);
1822 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
1824 all_eigenvalues_signal,
1825 all_hessenberg_signal,
1826 condition_number_signal);
1828 if (left_precondition)
1829 ::internal::SolverGMRESImplementation::add(
1835 arnoldi_process.vector_ptrs);
1838 ::internal::SolverGMRESImplementation::add(
1844 arnoldi_process.vector_ptrs);
1845 preconditioner.vmult(v, p);
1853 compute_eigs_and_cond(arnoldi_process.get_hessenberg_matrix(),
1857 condition_number_signal);
1859 if (!additional_data.batched_mode && !krylov_space_signal.empty())
1860 krylov_space_signal(basis_vectors);
1875template <
typename VectorType>
1877boost::signals2::connection
1879 const std::function<
void(
double)> &slot,
1880 const bool every_iteration)
1882 if (every_iteration)
1884 return all_condition_numbers_signal.connect(slot);
1888 return condition_number_signal.connect(slot);
1894template <
typename VectorType>
1897 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
1898 const bool every_iteration)
1900 if (every_iteration)
1902 return all_eigenvalues_signal.connect(slot);
1906 return eigenvalues_signal.connect(slot);
1912template <
typename VectorType>
1916 const bool every_iteration)
1918 if (every_iteration)
1920 return all_hessenberg_signal.connect(slot);
1924 return hessenberg_signal.connect(slot);
1930template <
typename VectorType>
1933 const std::function<
void(
1936 return krylov_space_signal.connect(slot);
1941template <
typename VectorType>
1943boost::signals2::connection
1945 const std::function<
void(
int)> &slot)
1947 return re_orthogonalize_signal.connect(slot);
1952template <
typename VectorType>
1966template <
typename VectorType>
1970 const AdditionalData &
data)
1972 , additional_data(
data)
1977template <
typename VectorType>
1980 const AdditionalData &
data)
1982 , additional_data(
data)
1987template <
typename VectorType>
1989template <
typename MatrixType,
typename PreconditionerType>
1995 const VectorType &b,
1996 const PreconditionerType &preconditioner)
2006 basis_size + 1, this->memory);
2008 basis_size, this->memory);
2012 unsigned int accumulated_iterations = 0;
2020 double res = std::numeric_limits<double>::lowest();
2030 if (accumulated_iterations == 0 && x.all_zero())
2034 A.vmult(v(0, x), x);
2035 v[0].sadd(-1., 1., b);
2038 res = arnoldi_process.orthonormalize_nth_vector(0, v);
2039 iteration_state = this->iteration_status(accumulated_iterations, res, x);
2043 unsigned int inner_iteration = 0;
2044 for (; (inner_iteration < basis_size &&
2048 preconditioner.vmult(z(inner_iteration, x), v[inner_iteration]);
2049 A.vmult(v(inner_iteration + 1, x), z[inner_iteration]);
2052 arnoldi_process.orthonormalize_nth_vector(inner_iteration + 1, v);
2059 this->iteration_status(++accumulated_iterations, res, x);
2065 arnoldi_process.solve_projected_system(
true);
2066 ::internal::SolverGMRESImplementation::add(
2072 arnoldi_process.vector_ptrs);
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())
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
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)
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
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
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)