16 #ifndef dealii_solver_gmres_h
17 #define dealii_solver_gmres_h
48 template <
typename,
typename>
64 namespace SolverGMRESImplementation
74 template <
typename VectorType>
122 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
191 template <
class VectorType = Vector<
double>>
296 template <
typename MatrixType,
typename PreconditionerType>
300 const VectorType &
b,
301 const PreconditionerType &preconditioner);
309 boost::signals2::connection
311 const bool every_iteration =
false);
319 boost::signals2::connection
321 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
322 const bool every_iteration =
false);
331 boost::signals2::connection
334 const bool every_iteration =
true);
342 boost::signals2::connection
353 boost::signals2::connection
359 <<
"The number of temporary vectors you gave (" << arg1
360 <<
") is too small. It should be at least 10 for "
361 <<
"any results, and much more for reasonable ones.");
385 boost::signals2::signal<void(
const std::vector<std::complex<double>> &)>
392 boost::signals2::signal<void(
const std::vector<std::complex<double>> &)>
412 boost::signals2::signal<void(
455 const unsigned int dim,
456 const boost::signals2::signal<
460 const boost::signals2::signal<
void(
double)> &cond_signal);
510 template <
class VectorType = Vector<
double>>
549 template <
typename MatrixType,
typename PreconditionerType>
553 const VectorType &
b,
554 const PreconditionerType &preconditioner);
580 namespace SolverGMRESImplementation
582 template <
class VectorType>
591 template <
class VectorType>
603 template <
class VectorType>
606 const VectorType & temp)
609 if (data[i] ==
nullptr)
612 data[i]->reinit(temp,
true);
619 template <
class VectorType>
623 return (data.size() > 0 ? data.size() - 1 : 0);
630 complex_less_pred(
const std::complex<double> &x,
631 const std::complex<double> &y)
633 return x.real() < y.real() ||
634 (x.real() == y.real() && x.imag() < y.imag());
640 solve_triangular(
const unsigned int dim,
645 for (
int i = dim - 1; i >= 0; --i)
648 for (
unsigned int j = i + 1; j < dim; ++j)
649 s -= solution(j) * H(i, j);
650 solution(i) = s / H(i, i);
659 template <
class VectorType>
661 const unsigned int max_n_tmp_vectors,
662 const bool right_preconditioning,
663 const bool use_default_residual,
664 const bool force_re_orthogonalization,
665 const bool batched_mode,
666 const OrthogonalizationStrategy orthogonalization_strategy)
667 : max_n_tmp_vectors(max_n_tmp_vectors)
668 , right_preconditioning(right_preconditioning)
669 , use_default_residual(use_default_residual)
670 , force_re_orthogonalization(force_re_orthogonalization)
671 , batched_mode(batched_mode)
672 , orthogonalization_strategy(orthogonalization_strategy)
674 Assert(3 <= max_n_tmp_vectors,
675 ExcMessage(
"SolverGMRES needs at least three "
676 "temporary vectors."));
681 template <
class VectorType>
684 const AdditionalData & data)
686 , additional_data(data)
692 template <
class VectorType>
694 const AdditionalData &data)
696 , additional_data(data)
702 template <
class VectorType>
710 for (
int i = 0; i < col; ++i)
712 const double s = si(i);
713 const double c = ci(i);
714 const double dummy = h(i);
715 h(i) = c * dummy + s * h(i + 1);
716 h(i + 1) = -s * dummy + c * h(i + 1);
719 const double r = 1. /
std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
720 si(col) = h(col + 1) * r;
721 ci(col) = h(col) * r;
722 h(col) = ci(col) * h(col) + si(col) * h(col + 1);
723 b(col + 1) = -si(col) *
b(col);
731 namespace SolverGMRESImplementation
733 template <
class VectorType>
735 Tvmult_add(
const unsigned int dim,
736 const VectorType & vv,
738 & orthogonal_vectors,
741 for (
unsigned int i = 0; i < dim; ++i)
742 h[i] += vv * orthogonal_vectors[i];
747 template <
class Number>
750 const unsigned int dim,
754 & orthogonal_vectors,
762 static constexpr
unsigned int n_lanes =
766 for (
unsigned int d = 0;
d < dim; ++
d)
772 ++c, j += n_lanes * 4)
773 for (
unsigned int i = 0; i < dim; ++i)
776 for (
unsigned int k = 0; k < 4; ++k)
777 vvec[k].load(vv.
begin() + j + k * n_lanes);
779 for (
unsigned int k = 0; k < 4; ++k)
782 temp.
load(orthogonal_vectors[i].
begin() + j + k * n_lanes);
783 hs[i] += temp * vvec[k];
789 for (
unsigned int i = 0; i < dim; ++i)
793 temp.
load(orthogonal_vectors[i].
begin() + j);
794 hs[i] += temp * vvec;
797 for (
unsigned int i = 0; i < dim; ++i)
798 for (
unsigned int v = 0; v < n_lanes; ++v)
805 for (
unsigned int i = 0; i < dim; ++i)
806 h(i) += orthogonal_vectors[i].local_element(j) * vv.
local_element(j);
813 template <
class VectorType>
816 const unsigned int dim,
818 & orthogonal_vectors,
824 for (
unsigned int i = 0; i < dim; ++i)
825 vv.add(-h(i), orthogonal_vectors[i]);
827 return std::sqrt(vv.add_and_dot(-h(dim), orthogonal_vectors[dim], vv));
832 template <
class Number>
835 const unsigned int dim,
838 & orthogonal_vectors,
844 double norm_vv_temp = 0.0;
853 for (
unsigned int k = 0; k < 4; ++k)
854 temp[k].load(vv.
begin() + j + k * n_lanes);
856 for (
unsigned int i = 0; i < dim; ++i)
858 const double factor = h(i);
859 for (
unsigned int k = 0; k < 4; ++k)
862 vec.
load(orthogonal_vectors[i].
begin() + j + k * n_lanes);
863 temp[k] -= factor * vec;
867 for (
unsigned int k = 0; k < 4; ++k)
868 temp[k].store(vv.
begin() + j + k * n_lanes);
870 norm_vv_temp_vectorized += (temp[0] * temp[0] + temp[1] * temp[1]) +
871 (temp[2] * temp[2] + temp[3] * temp[3]);
880 for (
unsigned int i = 0; i < dim; ++i)
883 vec.
load(orthogonal_vectors[i].
begin() + j);
889 norm_vv_temp_vectorized += temp * temp;
892 for (
unsigned int v = 0; v < n_lanes; ++v)
893 norm_vv_temp += norm_vv_temp_vectorized[v];
898 for (
unsigned int i = 0; i < dim; ++i)
899 temp -= h(i) * orthogonal_vectors[i].local_element(j);
902 norm_vv_temp += temp * temp;
909 template <
class VectorType>
911 sadd_and_norm(VectorType & v,
912 const double factor_a,
914 const double factor_b)
916 v.sadd(factor_a, factor_b,
b);
921 template <
class Number>
925 const double factor_a,
927 const double factor_b)
946 template <
class VectorType>
949 const unsigned int dim,
956 p.equ(h(0), tmp_vectors[0]);
958 p.add(h(0), tmp_vectors[0]);
960 for (
unsigned int i = 1; i < dim; ++i)
961 p.add(h(i), tmp_vectors[i]);
966 template <
class Number>
969 const unsigned int dim,
979 for (
unsigned int i = 0; i < dim; ++i)
980 temp += tmp_vectors[i].local_element(j) * h(i);
998 template <
class VectorType>
1000 iterated_gram_schmidt(
1002 OrthogonalizationStrategy orthogonalization_strategy,
1004 & orthogonal_vectors,
1005 const unsigned int dim,
1006 const unsigned int accumulated_iterations,
1009 bool & reorthogonalize,
1010 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal =
1011 boost::signals2::signal<
void(
int)>())
1014 const unsigned int inner_iteration = dim - 1;
1017 double norm_vv_start = 0;
1018 const bool consider_reorthogonalize =
1019 (reorthogonalize ==
false) && (inner_iteration % 5 == 4);
1020 if (consider_reorthogonalize)
1021 norm_vv_start = vv.l2_norm();
1023 for (
unsigned int i = 0; i < dim; ++i)
1026 for (
unsigned int c = 0; c < 2;
1030 double norm_vv = 0.0;
1032 if (orthogonalization_strategy ==
1034 OrthogonalizationStrategy::modified_gram_schmidt)
1036 double htmp = vv * orthogonal_vectors[0];
1038 for (
unsigned int i = 1; i < dim; ++i)
1040 htmp = vv.add_and_dot(-htmp,
1041 orthogonal_vectors[i - 1],
1042 orthogonal_vectors[i]);
1047 vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
1049 else if (orthogonalization_strategy ==
1051 OrthogonalizationStrategy::classical_gram_schmidt)
1053 Tvmult_add(dim, vv, orthogonal_vectors, h);
1054 norm_vv = substract_and_norm(dim, orthogonal_vectors, h, vv);
1072 if (consider_reorthogonalize)
1075 10. * norm_vv_start *
1076 std::sqrt(std::numeric_limits<
1077 typename VectorType::value_type>::
epsilon()))
1082 reorthogonalize =
true;
1083 if (!reorthogonalize_signal.empty())
1084 reorthogonalize_signal(accumulated_iterations);
1088 if (reorthogonalize ==
false)
1101 template <
class VectorType>
1105 const unsigned int dim,
1106 const boost::signals2::signal<
void(
const std::vector<std::complex<double>> &)>
1107 &eigenvalues_signal,
1109 & hessenberg_signal,
1110 const boost::signals2::signal<
void(
double)> &cond_signal)
1113 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1114 !cond_signal.empty()) &&
1118 for (
unsigned int i = 0; i < dim; ++i)
1119 for (
unsigned int j = 0; j < dim; ++j)
1120 mat(i, j) = H_orig(i, j);
1121 hessenberg_signal(H_orig);
1123 if (!eigenvalues_signal.empty())
1129 mat_eig.compute_eigenvalues();
1130 std::vector<std::complex<double>>
eigenvalues(dim);
1131 for (
unsigned int i = 0; i < mat_eig.n(); ++i)
1136 internal::SolverGMRESImplementation::complex_less_pred);
1141 if (!cond_signal.empty() && (mat.n() > 1))
1144 double condition_number =
1145 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
1146 cond_signal(condition_number);
1153 template <
class VectorType>
1154 template <
typename MatrixType,
typename PreconditionerType>
1158 const VectorType &
b,
1159 const PreconditionerType &preconditioner)
1165 std::unique_ptr<LogStream::Prefix> prefix;
1167 prefix = std::make_unique<LogStream::Prefix>(
"GMRES");
1171 const unsigned int n_tmp_vectors =
1176 n_tmp_vectors, this->memory);
1181 unsigned int accumulated_iterations = 0;
1183 const bool do_eigenvalues =
1185 (!condition_number_signal.empty() ||
1186 !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() ||
1187 !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
1188 !all_hessenberg_signal.empty());
1193 H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
1196 H.reinit(n_tmp_vectors, n_tmp_vectors - 1,
true);
1199 gamma.reinit(n_tmp_vectors);
1200 ci.
reinit(n_tmp_vectors - 1);
1201 si.
reinit(n_tmp_vectors - 1);
1202 h.
reinit(n_tmp_vectors - 1);
1204 unsigned int dim = 0;
1207 double last_res = std::numeric_limits<double>::lowest();
1219 VectorType &v = tmp_vectors(0, x);
1220 VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
1226 std::unique_ptr<::Vector<double>> gamma_;
1227 if (!use_default_residual)
1234 gamma_ = std::make_unique<::Vector<double>>(
gamma.size());
1246 h.
reinit(n_tmp_vectors - 1);
1250 if (left_precondition)
1254 preconditioner.vmult(v, p);
1260 rho = ::internal::SolverGMRESImplementation::sadd_and_norm(v,
1269 if (use_default_residual)
1273 iteration_state = solver_control.
check(accumulated_iterations, rho);
1276 this->iteration_status(accumulated_iterations, rho, x);
1283 deallog <<
"default_res=" << rho << std::endl;
1285 if (left_precondition)
1288 r->sadd(-1., 1.,
b);
1291 preconditioner.vmult(*r, v);
1293 double res = r->l2_norm();
1296 iteration_state = solver_control.
check(accumulated_iterations, rho);
1299 this->iteration_status(accumulated_iterations, res, x);
1312 for (
unsigned int inner_iteration = 0;
1313 ((inner_iteration < n_tmp_vectors - 2) &&
1317 ++accumulated_iterations;
1319 VectorType &vv = tmp_vectors(inner_iteration + 1, x);
1321 if (left_precondition)
1323 A.vmult(p, tmp_vectors[inner_iteration]);
1324 preconditioner.vmult(vv, p);
1328 preconditioner.vmult(p, tmp_vectors[inner_iteration]);
1332 dim = inner_iteration + 1;
1335 internal::SolverGMRESImplementation::iterated_gram_schmidt(
1339 accumulated_iterations,
1343 re_orthogonalize_signal);
1344 h(inner_iteration + 1) = s;
1354 for (
unsigned int i = 0; i < dim + 1; ++i)
1355 H_orig(i, inner_iteration) = h(i);
1361 for (
unsigned int i = 0; i < dim; ++i)
1362 H(i, inner_iteration) = h(i);
1367 if (use_default_residual)
1372 solver_control.
check(accumulated_iterations, rho);
1375 this->iteration_status(accumulated_iterations, rho, x);
1380 deallog <<
"default_res=" << rho << std::endl;
1384 internal::SolverGMRESImplementation::solve_triangular(dim,
1389 if (left_precondition)
1390 for (
unsigned int i = 0; i < dim; ++i)
1391 x_->add(h(i), tmp_vectors[i]);
1395 for (
unsigned int i = 0; i < dim; ++i)
1396 p.add(h(i), tmp_vectors[i]);
1397 preconditioner.vmult(*r, p);
1401 r->sadd(-1., 1.,
b);
1403 if (left_precondition)
1405 const double res = r->l2_norm();
1409 this->iteration_status(accumulated_iterations, res, x);
1413 preconditioner.vmult(*x_, *r);
1414 const double preconditioned_res = x_->l2_norm();
1415 last_res = preconditioned_res;
1419 solver_control.
check(accumulated_iterations, rho);
1422 this->iteration_status(accumulated_iterations,
1431 internal::SolverGMRESImplementation::solve_triangular(dim, H,
gamma, h);
1434 compute_eigs_and_cond(H_orig,
1436 all_eigenvalues_signal,
1437 all_hessenberg_signal,
1438 condition_number_signal);
1440 if (left_precondition)
1441 ::internal::SolverGMRESImplementation::add(
1442 x, dim, h, tmp_vectors,
false);
1445 ::internal::SolverGMRESImplementation::add(
1446 p, dim, h, tmp_vectors,
true);
1447 preconditioner.vmult(v, p);
1456 compute_eigs_and_cond(H_orig,
1460 condition_number_signal);
1462 if (!additional_data.
batched_mode && !krylov_space_signal.empty())
1463 krylov_space_signal(tmp_vectors);
1472 template <
class VectorType>
1473 boost::signals2::connection
1475 const std::function<
void(
double)> &slot,
1476 const bool every_iteration)
1478 if (every_iteration)
1480 return all_condition_numbers_signal.connect(slot);
1484 return condition_number_signal.connect(slot);
1490 template <
class VectorType>
1491 boost::signals2::connection
1493 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
1494 const bool every_iteration)
1496 if (every_iteration)
1498 return all_eigenvalues_signal.connect(slot);
1502 return eigenvalues_signal.connect(slot);
1508 template <
class VectorType>
1509 boost::signals2::connection
1512 const bool every_iteration)
1514 if (every_iteration)
1516 return all_hessenberg_signal.connect(slot);
1520 return hessenberg_signal.connect(slot);
1526 template <
class VectorType>
1527 boost::signals2::connection
1529 const std::function<
void(
1532 return krylov_space_signal.connect(slot);
1537 template <
class VectorType>
1538 boost::signals2::connection
1540 const std::function<
void(
int)> &slot)
1542 return re_orthogonalize_signal.connect(slot);
1547 template <
class VectorType>
1560 template <
class VectorType>
1563 const AdditionalData & data)
1565 , additional_data(data)
1570 template <
class VectorType>
1572 const AdditionalData &data)
1574 , additional_data(data)
1579 template <
class VectorType>
1580 template <
typename MatrixType,
typename PreconditionerType>
1584 const VectorType &
b,
1585 const PreconditionerType &preconditioner)
1591 const unsigned int basis_size = additional_data.max_basis_size;
1595 basis_size, this->memory);
1597 basis_size, this->memory);
1601 unsigned int accumulated_iterations = 0;
1604 H.reinit(basis_size + 1, basis_size);
1611 double res = std::numeric_limits<double>::lowest();
1618 aux->sadd(-1., 1.,
b);
1620 double beta = aux->l2_norm();
1622 iteration_state = this->iteration_status(accumulated_iterations, res, x);
1626 H.reinit(basis_size + 1, basis_size);
1629 for (
unsigned int j = 0; j < basis_size; ++j)
1632 v(j, x).equ(1. / a, *aux);
1637 preconditioner.vmult(z(j, x), v[j]);
1638 A.vmult(*aux, z[j]);
1641 H(0, j) = *aux * v[0];
1642 for (
unsigned int i = 1; i <= j; ++i)
1643 H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1644 H(j + 1, j) = a =
std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1650 H1.reinit(j + 1, j);
1651 projected_rhs.
reinit(j + 1);
1653 projected_rhs(0) = beta;
1661 res = house.least_squares(y, projected_rhs);
1663 this->iteration_status(++accumulated_iterations, res, x);
1670 for (
unsigned int j = 0; j < y.
size(); ++j)
Number local_element(const size_type local_index) const
size_type locally_owned_size() const
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())
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
static void compute_eigs_and_cond(const FullMatrix< double > &H_orig, const unsigned int dim, 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)
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
SolverControl & solver_control
SolverGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
void givens_rotation(Vector< double > &h, Vector< double > &b, Vector< double > &ci, Vector< double > &si, int col) const
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())
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< std::complex< double >> &)> &slot, const bool every_iteration=false)
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::signal< void(int)> re_orthogonalize_signal
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
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)
VectorMemory< VectorType > & mem
std::vector< typename VectorMemory< VectorType >::Pointer > data
VectorType & operator()(const unsigned int i, const VectorType &temp)
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
unsigned int size() const
VectorType & operator[](const unsigned int i) const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
static ::ExceptionBase & ExcTooFewTmpVectors(int arg1)
#define AssertIndexRange(index, range)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Expression fabs(const Expression &x)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
T sum(const T &t, const MPI_Comm &mpi_communicator)
long double gamma(const unsigned int n)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
AdditionalData(const unsigned int max_basis_size=30)
unsigned int max_basis_size
bool right_preconditioning
OrthogonalizationStrategy
bool force_re_orthogonalization
OrthogonalizationStrategy orthogonalization_strategy
AdditionalData(const unsigned int max_n_tmp_vectors=30, const bool right_preconditioning=false, const bool use_default_residual=true, const bool force_re_orthogonalization=false, const bool batched_mode=false, const OrthogonalizationStrategy orthogonalization_strategy=OrthogonalizationStrategy::modified_gram_schmidt)
bool use_default_residual
unsigned int max_n_tmp_vectors