18#ifndef ITERATIVESOLVER_HPP
19#define ITERATIVESOLVER_HPP
24template<
typename T,
typename data_t>
concept HasEvaluate =
requires(
T* t,
const Col<data_t>& x) { { t->evaluate(x) } -> std::convertible_to<Col<data_t>>; };
27 constexpr sp_d auto ZERO = data_t(0);
28 constexpr sp_d auto ONE = data_t(1);
32 auto generate_rotation = [](
const data_t dx,
const data_t dy, data_t& cs, data_t& sn) ->
void {
37 else if(std::fabs(dy) > std::fabs(dx)) {
38 const data_t fraction = dx / dy;
39 sn =
ONE / std::sqrt(
ONE + fraction * fraction);
43 const data_t fraction = dy / dx;
44 cs =
ONE / std::sqrt(
ONE + fraction * fraction);
49 auto apply_rotation = [](data_t& dx, data_t& dy,
const data_t cs,
const data_t sn) ->
void {
50 const data_t factor = cs * dx + sn * dy;
51 dy = cs * dy - sn * dx;
55 if(x.empty()) x = conditioner->apply(b);
56 else x.zeros(arma::size(b));
58 const auto mp = setting.
restart + 1;
60 Mat<data_t> hessenberg(mp, setting.
restart, fill::zeros);
63 data_t beta, residual;
64 Col<data_t> s(mp, fill::none), cs(mp, fill::none), sn(mp, fill::none), r;
66 auto norm_b =
arma::norm(conditioner->apply(b));
69 auto stop_criterion = [&] {
70 residual = (beta =
arma::norm(r = conditioner->apply(b - system->
evaluate(x)))) / norm_b;
71 suanpan_debug(
"GMRES iterative solver residual: {:.5E}.\n", residual);
80 Mat<data_t> v(b.n_rows, mp, fill::none);
82 auto update = [&](
const int k) -> Col<data_t> {
83 Col<data_t> y = s.head(k + 1llu);
85 for(
auto i = k; i >= 0; --i) {
86 y(i) /= hessenberg(i, i);
87 y.head(i) -= hessenberg.col(i).head(i) * y(i);
90 return v.head_cols(k + 1llu) * y;
98 for(
auto i = 0, j = 1; i < setting.
restart && counter <= setting.
max_iteration; ++i, ++j, ++counter) {
99 auto w = conditioner->apply(system->
evaluate(v.col(i)));
100 for(
auto k = 0; k <= i; ++k) w -= (hessenberg(k, i) = arma::dot(w, v.col(k))) * v.col(k);
101 v.col(j) = w / (hessenberg(j, i) =
arma::norm(w));
103 for(
auto k = 0; k < i; ++k) apply_rotation(hessenberg(k, i), hessenberg(k + 1llu, i), cs(k), sn(k));
105 generate_rotation(hessenberg(i, i), hessenberg(j, i), cs(i), sn(i));
106 apply_rotation(hessenberg(i, i), hessenberg(j, i), cs(i), sn(i));
107 apply_rotation(s(i), s(j), cs(i), sn(i));
109 residual = std::fabs(s(j)) / norm_b;
110 suanpan_debug(
"GMRES iterative solver residual: {:.5E}.\n", residual);
119 x += update(setting.
restart - 1);
128 constexpr sp_d auto ZERO = data_t(0);
129 constexpr sp_d auto ONE = data_t(1);
136 if(x.empty()) x = conditioner->apply(b);
137 else x.zeros(arma::size(b));
139 Col<data_t> r = b - system->
evaluate(x);
140 const auto initial_r = r;
143 suanpan_debug(
"BiCGSTAB iterative solver residual: {:.5E}.\n", residual);
154 const auto rho = arma::dot(initial_r, r);
162 else p = r + rho / pre_rho * alpha / omega * (p - omega * v);
164 const auto phat = conditioner->apply(p);
166 alpha = rho / arma::dot(initial_r, v);
167 const Col<data_t> s = r - alpha * v;
177 const auto shat = conditioner->apply(s);
178 const Col<data_t> t = system->
evaluate(shat);
179 omega = arma::dot(t, s) / arma::dot(t, t);
180 x += alpha * phat + omega * shat;
int BiCGSTAB(const System *system, Col< data_t > &x, const Col< data_t > &b, SolverSetting< data_t > &setting)
Definition: IterativeSolver.hpp:127
int GMRES(const System *system, Col< data_t > &x, const Col< data_t > &b, SolverSetting< data_t > &setting)
Definition: IterativeSolver.hpp:26
Definition: IterativeSolver.hpp:24
Definition: suanPan.h:330
std::enable_if_t<!std::numeric_limits< T >::is_integer, bool > approx_equal(T x, T y, int ulp=2)
Definition: utility.h:60
double norm(const vec &)
Definition: tensor.cpp:370
Definition: SolverSetting.hpp:40
data_t tolerance
Definition: SolverSetting.hpp:43
int restart
Definition: SolverSetting.hpp:41
int max_iteration
Definition: SolverSetting.hpp:42
Preconditioner< data_t > * preconditioner
Definition: SolverSetting.hpp:48
Definition: TestSolver.h:6
vec evaluate(const vec &x) const
Definition: TestSolver.h:12
#define suanpan_debug(...)
Definition: suanPan.h:307
constexpr auto SUANPAN_SUCCESS
Definition: suanPan.h:172
constexpr auto SUANPAN_FAIL
Definition: suanPan.h:173