3 #include <AMReX_Config.H>
77 template <
typename V,
typename M>
82 using RT =
typename M::RT;
100 void solve (V& a_sol, V
const& a_rhs,
RT a_tol_rel,
RT a_tol_abs,
int a_its=-1);
123 void cycle (V& a_xx,
int& a_status,
int& a_itcount,
RT& a_rnorm0);
153 template <
typename V,
typename M>
159 template <
typename V,
typename M>
162 int rs = m_restrtlen;
164 m_hh_1d.resize(std::size_t(rs + 2) * (rs + 1));
165 m_hh =
Table2D<RT>(m_hh_1d.data(), {0,0}, {rs+1,rs});
167 m_hes_1d.resize(std::size_t(rs + 2) * (rs + 1));
168 m_hes =
Table2D<RT>(m_hes_1d.data(), {0,0}, {rs+1,rs});
170 m_grs.resize(rs + 2);
175 template <
typename V,
typename M>
178 if (m_restrtlen != rl) {
185 template <
typename V,
typename M>
192 template <
typename V,
typename M>
204 template <
typename V,
typename M>
207 return (
r < r0*m_rtol) || (
r < m_atol);
210 template <
typename V,
typename M>
219 if (m_v_tmp_rhs ==
nullptr) {
220 m_v_tmp_rhs = std::make_unique<V>(m_linop->makeVecRHS());
222 if (m_v_tmp_lhs ==
nullptr) {
223 m_v_tmp_lhs = std::make_unique<V>(m_linop->makeVecLHS());
226 m_vv.reserve(m_restrtlen+1);
227 for (
int i = 0; i < 2; ++i) {
228 m_vv.emplace_back(m_linop->makeVecRHS());
235 if (a_its < 0) { a_its = m_maxiter; }
239 m_linop->assign(m_vv[0], a_rhs);
240 m_linop->setToZero(a_sol);
244 cycle(a_sol, m_status, m_its, rnorm0);
246 while (m_status == -1 && m_its < a_its) {
247 compute_residual(m_vv[0], a_sol, a_rhs);
248 cycle(a_sol, m_status, m_its, rnorm0);
251 if (m_status == -1 && m_its >= a_its) { m_status = 1; }
259 amrex::Print() <<
"GMRES: Solve Time = " << t1-t0 <<
'\n';
263 template <
typename V,
typename M>
268 m_res = m_linop->norm2(m_vv[0]);
271 if (m_res ==
RT(0.0)) {
276 m_linop->scale(m_vv[0],
RT(1.0)/m_res);
278 if (a_itcount == 0) { a_rnorm0 = m_res; }
280 a_status = converged(a_rnorm0,m_res) ? 0 : -1;
283 while (it < m_restrtlen && a_itcount < m_maxiter)
287 <<
", residual = " << m_res <<
", " << m_res/a_rnorm0
291 if (a_status == 0) {
break; }
293 while (m_vv.size() < it+2) {
294 m_vv.emplace_back(m_linop->makeVecRHS());
297 auto const& vv_it = m_vv[it ];
298 auto & vv_it1 = m_vv[it+1];
300 m_linop->precond(*m_v_tmp_lhs, vv_it);
301 m_linop->apply(vv_it1, *m_v_tmp_lhs);
303 gram_schmidt_orthogonalization(it);
305 auto tt = m_linop->norm2(vv_it1);
307 auto const small =
RT((
sizeof(
RT) == 8) ? 1.e-99 : 1.e-30);
308 bool happyend = (tt < small);
310 m_linop->scale(vv_it1,
RT(1.0)/tt);
316 update_hessenberg(it, happyend, m_res);
320 a_status = converged(a_rnorm0, m_res) ? 0 : -1;
321 if (happyend) {
break; }
324 if ((m_verbose > 1) && (a_status != 0 || a_itcount >= m_maxiter)) {
326 <<
", residual = " << m_res <<
", " << m_res/a_rnorm0
330 build_solution(a_xx, it-1);
333 template <
typename V,
typename M>
340 auto& vv_1 = m_vv[it+1];
344 for (
int j = 0; j <= it; ++j) {
345 m_hh (j,it) =
RT(0.0);
346 m_hes(j,it) =
RT(0.0);
349 for (
int ncnt = 0; ncnt < 2 ; ++ncnt)
351 for (
int j = 0; j <= it; ++j) {
352 lhh[j] = m_linop->dotProduct(vv_1, m_vv[j]);
355 for (
int j = 0; j <= it; ++j) {
356 m_linop->increment(vv_1, m_vv[j], -lhh[j]);
357 m_hh (j,it) += lhh[j];
358 m_hes(j,it) -= lhh[j];
363 template <
typename V,
typename M>
368 for (
int j = 1; j <= it; ++j) {
369 auto tt = m_hh(j-1,it);
370 m_hh(j-1,it) = m_cc[j-1] * tt + m_ss[j-1] * m_hh(j,it);
371 m_hh(j ,it) = m_cc[j-1] * m_hh(j,it) - m_ss[j-1] * tt;
376 auto tt =
std::sqrt(m_hh(it,it)*m_hh(it,it) + m_hh(it+1,it)*m_hh(it+1,it));
377 m_cc[it] = m_hh(it ,it) / tt;
378 m_ss[it] = m_hh(it+1,it) / tt;
379 m_grs[it+1] = - (m_ss[it] * m_grs[it]);
380 m_grs[it ] = m_cc[it] * m_grs[it];
381 m_hh(it,it) = m_cc[it] * m_hh(it,it) + m_ss[it] * m_hh(it+1,it);
390 template <
typename V,
typename M>
395 if (it < 0) {
return; }
397 if (m_hh(it,it) !=
RT(0.0)) {
398 m_grs[it] /= m_hh(it,it);
403 for (
int ii = 1; ii <= it; ++ii) {
406 for (
int j = k+1; j <= it; ++j) {
407 tt -= m_hh(k,j) * m_grs[j];
409 m_grs[k] = tt / m_hh(k,k);
412 m_linop->setToZero(*m_v_tmp_rhs);
413 for (
int ii = 0; ii < it+1; ++ii) {
414 m_linop->increment(*m_v_tmp_rhs, m_vv[ii], m_grs[ii]);
417 m_linop->precond(*m_v_tmp_lhs, *m_v_tmp_rhs);
418 m_linop->increment(a_xx, *m_v_tmp_lhs,
RT(1.0));
421 template <
typename V,
typename M>
425 m_linop->assign(*m_v_tmp_lhs, a_xx);
426 m_linop->apply(*m_v_tmp_rhs, *m_v_tmp_lhs);
427 m_linop->linComb(a_rr,
RT(1.0), a_bb,
RT(-1.0), *m_v_tmp_rhs);
#define BL_PROFILE(a)
Definition: AMReX_BLProfiler.H:551
#define AMREX_ALWAYS_ASSERT(EX)
Definition: AMReX_BLassert.H:50
GMRES.
Definition: AMReX_GMRES.H:79
int getNumIters() const
Gets the number of iterations.
Definition: AMReX_GMRES.H:112
void solve(V &a_sol, V const &a_rhs, RT a_tol_rel, RT a_tol_abs, int a_its=-1)
Solve the linear system.
Definition: AMReX_GMRES.H:211
int m_status
Definition: AMReX_GMRES.H:135
void build_solution(V &a_xx, int it)
Definition: AMReX_GMRES.H:391
void allocate_scratch()
Definition: AMReX_GMRES.H:160
int m_verbose
Definition: AMReX_GMRES.H:132
void define(M &linop)
Definition: AMReX_GMRES.H:186
RT m_rtol
Definition: AMReX_GMRES.H:138
bool converged(RT r0, RT r) const
Definition: AMReX_GMRES.H:205
Vector< RT > m_hh_1d
Definition: AMReX_GMRES.H:140
RT m_atol
Definition: AMReX_GMRES.H:139
int m_maxiter
Definition: AMReX_GMRES.H:133
void cycle(V &a_xx, int &a_status, int &a_itcount, RT &a_rnorm0)
Definition: AMReX_GMRES.H:264
Table2D< RT > m_hes
Definition: AMReX_GMRES.H:143
Table2D< RT > m_hh
Definition: AMReX_GMRES.H:142
void setVerbose(int v)
Sets verbosity.
Definition: AMReX_GMRES.H:103
void compute_residual(V &a_rr, V const &a_xx, V const &a_bb)
Definition: AMReX_GMRES.H:422
Vector< RT > m_cc
Definition: AMReX_GMRES.H:145
M * m_linop
Definition: AMReX_GMRES.H:150
int m_its
Definition: AMReX_GMRES.H:134
int m_restrtlen
Definition: AMReX_GMRES.H:136
void setMaxIters(int niters)
Sets the max number of iterations.
Definition: AMReX_GMRES.H:109
std::unique_ptr< V > m_v_tmp_lhs
Definition: AMReX_GMRES.H:148
Vector< RT > m_hes_1d
Definition: AMReX_GMRES.H:141
void setRestartLength(int rl)
Sets restart length. The default is 30.
Definition: AMReX_GMRES.H:176
void update_hessenberg(int it, bool happyend, RT &res)
Definition: AMReX_GMRES.H:364
void gram_schmidt_orthogonalization(int it)
Definition: AMReX_GMRES.H:334
RT getResidualNorm() const
Gets the 2-norm of the residual.
Definition: AMReX_GMRES.H:118
GMRES()
Definition: AMReX_GMRES.H:154
Vector< V > m_vv
Definition: AMReX_GMRES.H:149
Vector< RT > m_grs
Definition: AMReX_GMRES.H:144
Vector< RT > m_ss
Definition: AMReX_GMRES.H:146
void clear()
Definition: AMReX_GMRES.H:193
int getStatus() const
Gets the solver status.
Definition: AMReX_GMRES.H:115
std::unique_ptr< V > m_v_tmp_rhs
Definition: AMReX_GMRES.H:147
RT m_res
Definition: AMReX_GMRES.H:137
typename M::RT RT
Definition: AMReX_GMRES.H:82
This class provides the user with a few print options.
Definition: AMReX_Print.H:35
@ max
Definition: AMReX_ParallelReduce.H:17
static constexpr int M
Definition: AMReX_OpenBC.H:13
Definition: AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Return the absolute value of a complex number.
Definition: AMReX_GpuComplex.H:356
double second() noexcept
Definition: AMReX_Utility.cpp:922
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Return the square root of a complex number.
Definition: AMReX_GpuComplex.H:373