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,
172 template <
typename VectorType>
175 const unsigned int n,
179 boost::signals2::signal<
void(
int)>());
284 std::vector<std::pair<double, double>> &
rotations,
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>
473 const PreconditionerType &preconditioner);
492 connect_eigenvalues_slot(
504 connect_hessenberg_slot(
515 connect_krylov_space_slot(
517 void(
const internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
526 connect_re_orthogonalization_slot(
const std::function<
void(
int)> &
slot);
565 all_eigenvalues_signal;
578 all_hessenberg_signal;
585 const internal::SolverGMRESImplementation::TmpVectors<VectorType> &)>
614 compute_eigs_and_cond(
616 const unsigned int n,
618 void(
const std::vector<
std::complex<
double>> &)> &eigenvalues_signal,
627 internal::SolverGMRESImplementation::ArnoldiProcess<
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>
711 const PreconditionerType &preconditioner);
723 internal::SolverGMRESImplementation::ArnoldiProcess<
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)
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>
834 template <
typename VectorType>
839 static constexpr bool value =
844 std::is_same_v<VectorType, Vector<typename VectorType::value_type>>;
849 template <
typename 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,
943 std::vector<const typename VectorType::value_type *> &)
945 for (
unsigned int i = 0; i < n; ++i)
948 if (delayed_reorthogonalization)
951 if (delayed_reorthogonalization)
958 template <
bool delayed_reorthogonalization,
typename Number>
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,
978 std::vector<const typename VectorType::value_type *> &vector_ptrs)
982 vector_ptrs.resize(n);
983 for (
unsigned int i = 0; i < n; ++i)
999 template <
bool delayed_reorthogonalization,
1001 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1008 std::vector<const typename VectorType::value_type *> &)
1014 for (
unsigned int i = 0; i < n - 1; ++i)
1016 if (delayed_reorthogonalization && i + 2 < n)
1021 if (delayed_reorthogonalization)
1025 -h(n + n - 2) / h(n + n - 1),
1030 1. / (h(n + n - 1) * h(n + n)) :
1031 1. / (h(n + n - 1) * h(n + n - 1));
1038 return std::numeric_limits<double>::signaling_NaN();
1048 template <
bool delayed_reorthogonalization,
typename Number>
1058 template <
bool delayed_reorthogonalization,
1060 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1064 const unsigned int n,
1068 std::vector<const typename VectorType::value_type *> &vector_ptrs)
1074 vector_ptrs.resize(n);
1075 for (
unsigned int i = 0; i < n; ++i)
1093 std::enable_if_t<!is_dealii_compatible_vector<VectorType>::value,
1097 const unsigned int n,
1101 std::vector<const typename VectorType::value_type *> &)
1108 for (
unsigned int i = 1; i < n; ++i)
1115 template <
typename Number>
1117 do_add(
const unsigned int n_vectors,
1127 std::enable_if_t<is_dealii_compatible_vector<VectorType>::value,
1131 const unsigned int n,
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(
1160 this->orthogonalization_strategy = orthogonalization_strategy;
1170 if (orthogonalization_strategy ==
1180 template <
typename Number>
1181 template <
typename VectorType>
1183 ArnoldiProcess<Number>::orthonormalize_nth_vector(
1184 const unsigned int n,
1197 givens_rotations.clear();
1203 else if (orthogonalization_strategy ==
1216 h.reinit(n + n + 1);
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 ?
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)
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);
1269 true, n - 2, triangular_matrix, givens_rotations, projected_rhs);
1277 (do_reorthogonalization ==
false) && (n % 5 == 0);
1285 for (
unsigned int c = 0; c < 2; ++c)
1288 if (orthogonalization_strategy ==
1294 for (
unsigned int i = 1; i < n; ++i)
1305 else if (orthogonalization_strategy ==
1334 typename VectorType::value_type>::epsilon()))
1339 do_reorthogonalization =
true;
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;
1359 false, n - 1, triangular_matrix, givens_rotations, projected_rhs);
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)
1394 const double Hi1 = hessenberg_matrix(i + 1, col);
1395 H_next = -s *
H_next + c * hessenberg_matrix(i + 1, col + 1);
1402 const double H_col1 = hessenberg_matrix(col + 1, col);
1414 rotations[col].first * hessenberg_matrix(col + 1, col + 1);
1417 const double H_last = hessenberg_matrix(col + 2, col + 1);
1425 matrix(0, col) = hessenberg_matrix(0, col);
1426 for (
int i = 0; i < col; ++i)
1431 const double Hi1 = hessenberg_matrix(i + 1, col);
1436 const double Hi =
matrix(col, col);
1437 const double Hi1 = hessenberg_matrix(col + 1, col);
1452 template <
typename Number>
1454 ArnoldiProcess<Number>::solve_projected_system(
1461 unsigned int n = givens_rotations.size();
1471 if (orthogonalization_strategy ==
1482 do_givens_rotation(
false,
1483 givens_rotations.size(),
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;
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() ||
1553 for (
unsigned int i = 0; i < n; ++i)
1554 for (
unsigned int j = 0;
j < n; ++
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);
1580 mat.singular_value(0) /
mat.singular_value(
mat.n() - 1);
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");
1608 std::max(additional_data.max_n_tmp_vectors, 3u) - 2);
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();
1635 const bool use_default_residual = additional_data.use_default_residual;
1644 if (!use_default_residual)
1654 additional_data.force_re_orthogonalization);
1674 preconditioner.vmult(v, b);
1679 preconditioner.vmult(v, p);
1693 const double norm_v = arnoldi_process.orthonormalize_nth_vector(
1699 if (use_default_residual)
1702 if (additional_data.batched_mode)
1718 r->sadd(-1., 1., b);
1721 preconditioner.vmult(*
r, v);
1724 if (additional_data.batched_mode)
1748 preconditioner.vmult(
vv, p);
1760 re_orthogonalize_signal);
1762 if (use_default_residual)
1764 if (additional_data.batched_mode)
1773 if (!additional_data.batched_mode)
1774 deallog <<
"default_res=" <<
res << std::endl;
1778 arnoldi_process.solve_projected_system(
false);
1788 preconditioner.vmult(*
r, p);
1792 r->sadd(-1., 1., b);
1803 preconditioner.vmult(*
x_, *
r);
1804 res =
x_->l2_norm();
1806 if (additional_data.batched_mode)
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);
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())
1875template <
typename VectorType>
1877boost::signals2::connection
1879 const std::function<
void(
double)> &
slot,
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,
1902 return all_eigenvalues_signal.connect(
slot);
1906 return eigenvalues_signal.connect(
slot);
1912template <
typename VectorType>
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)
2020 double res = std::numeric_limits<double>::lowest();
2034 A.vmult(v(0,
x),
x);
2035 v[0].sadd(-1., 1., b);
2038 res = arnoldi_process.orthonormalize_nth_vector(0, v);
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)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
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
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()
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)