Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_GMRES.H
Go to the documentation of this file.
1#ifndef AMREX_GMRES_H_
2#define AMREX_GMRES_H_
3#include <AMReX_Config.H>
4
5#include <AMReX_BLProfiler.H>
6#include <AMReX_Print.H>
7#include <AMReX_TableData.H>
8#include <AMReX_Vector.H>
9#include <cmath>
10#include <limits>
11#include <memory>
12
13namespace amrex {
14
78template <typename V, typename M>
79class GMRES
80{
81public:
82
83 using RT = typename M::RT; // double or float
84
86
90 void define (M& linop);
91
101 void solve (V& a_sol, V const& a_rhs, RT a_tol_rel, RT a_tol_abs, int a_its=-1);
102
104 void setVerbose (int v) { m_verbose = v; }
105
107 void setRestartLength (int rl);
108
110 void setMaxIters (int niters) { m_maxiter = niters; }
111
113 [[nodiscard]] int getNumIters () const { return m_its; }
114
116 [[nodiscard]] int getStatus () const { return m_status; }
117
119 [[nodiscard]] RT getResidualNorm () const { return m_res; }
120
121private:
122 void clear ();
124 void cycle (V& a_xx, int& a_status, int& a_itcount, RT& a_rnorm0);
125 void build_solution (V& a_xx, int it);
126 void compute_residual (V& a_rr, V const& a_xx, V const& a_bb);
127
128 bool converged (RT r0, RT r) const;
129
131 void update_hessenberg (int it, bool happyend, RT& res);
132
133 int m_verbose = 0;
134 int m_maxiter = 2000;
135 int m_its = 0;
136 int m_status = -1;
137 int m_restrtlen = 30;
138 RT m_res = std::numeric_limits<RT>::max();
148 std::unique_ptr<V> m_v_tmp_rhs;
149 std::unique_ptr<V> m_v_tmp_lhs;
151 M* m_linop = nullptr;
152};
153
154template <typename V, typename M>
156{
157 allocate_scratch();
158}
159
160template <typename V, typename M>
162{
163 int rs = m_restrtlen;
164
165 m_hh_1d.resize(std::size_t(rs + 2) * (rs + 1));
166 m_hh = Table2D<RT>(m_hh_1d.data(), {0,0}, {rs+1,rs}); // (0:rs+1,0:rs)
167
168 m_hes_1d.resize(std::size_t(rs + 2) * (rs + 1));
169 m_hes = Table2D<RT>(m_hes_1d.data(), {0,0}, {rs+1,rs}); // (0:rs+1,0:rs)
170
171 m_grs.resize(rs + 2);
172 m_cc.resize(rs + 1);
173 m_ss.resize(rs + 1);
174}
175
176template <typename V, typename M>
178{
179 if (m_restrtlen != rl) {
180 m_restrtlen = rl;
181 allocate_scratch();
182 m_vv.clear();
183 }
184}
185
186template <typename V, typename M>
187void GMRES<V,M>::define (M& linop)
188{
189 clear();
190 m_linop = &linop;
191}
192
193template <typename V, typename M>
195{
196 m_its = 0;
197 m_status = -1;
198 m_res = std::numeric_limits<RT>::max();
199 m_v_tmp_rhs.reset();
200 m_v_tmp_lhs.reset();
201 m_vv.clear();
202 m_linop = nullptr;
203}
204
205template <typename V, typename M>
206bool GMRES<V,M>::converged (RT r0, RT r) const
207{
208 return (r < r0*m_rtol) || (r < m_atol);
209}
210
211template <typename V, typename M>
212void GMRES<V,M>::solve (V& a_sol, V const& a_rhs, RT a_tol_rel, RT a_tol_abs, int a_its)
213{
214 BL_PROFILE("GMRES::solve()");
215
216 AMREX_ALWAYS_ASSERT(m_linop != nullptr);
217
218 auto t0 = amrex::second();
219
220 if (m_v_tmp_rhs == nullptr) {
221 m_v_tmp_rhs = std::make_unique<V>(m_linop->makeVecRHS());
222 }
223 if (m_v_tmp_lhs == nullptr) {
224 m_v_tmp_lhs = std::make_unique<V>(m_linop->makeVecLHS());
225 }
226 if (m_vv.empty()) {
227 m_vv.reserve(m_restrtlen+1);
228 for (int i = 0; i < 2; ++i) { // to save space, start with just 2
229 m_vv.emplace_back(m_linop->makeVecRHS());
230 }
231 }
232
233 m_rtol = a_tol_rel;
234 m_atol = a_tol_abs;
235
236 if (a_its < 0) { a_its = m_maxiter; }
237
238 auto rnorm0 = RT(0);
239
240 m_linop->assign(m_vv[0], a_rhs);
241 m_linop->setToZero(a_sol);
242
243 m_its = 0;
244 m_status = -1;
245 cycle(a_sol, m_status, m_its, rnorm0);
246
247 while (m_status == -1 && m_its < a_its) {
248 compute_residual(m_vv[0], a_sol, a_rhs);
249 cycle(a_sol, m_status, m_its, rnorm0);
250 }
251
252 if (m_status == -1 && m_its >= a_its) { m_status = 1; }
253
254 m_v_tmp_rhs.reset();
255 m_v_tmp_lhs.reset();
256 m_vv.clear();
257
258 auto t1 = amrex::second();
259 if (m_verbose > 0) {
260 amrex::Print() << "GMRES: Solve Time = " << t1-t0 << '\n';
261 }
262}
263
264template <typename V, typename M>
265void GMRES<V,M>::cycle (V& a_xx, int& a_status, int& a_itcount, RT& a_rnorm0)
266{
267 BL_PROFILE("GMRES::cycle()");
268
269 m_res = m_linop->norm2(m_vv[0]);
270 m_grs[0] = m_res;
271
272 if (m_res == RT(0.0)) {
273 a_status = 0;
274 return;
275 }
276
277 m_linop->scale(m_vv[0], RT(1.0)/m_res);
278
279 if (a_itcount == 0) { a_rnorm0 = m_res; }
280
281 a_status = converged(a_rnorm0,m_res) ? 0 : -1;
282
283 int it = 0;
284 while (it < m_restrtlen && a_itcount < m_maxiter)
285 {
286 if (m_verbose > 1) {
287 amrex::Print() << "GMRES: iter = " << a_itcount
288 << ", residual = " << m_res << ", " << m_res/a_rnorm0
289 << " (rel.)\n";
290 }
291
292 if (a_status == 0) { break; }
293
294 while (m_vv.size() < it+2) {
295 m_vv.emplace_back(m_linop->makeVecRHS());
296 }
297
298 auto const& vv_it = m_vv[it ];
299 auto & vv_it1 = m_vv[it+1];
300
301 m_linop->precond(*m_v_tmp_lhs, vv_it);
302 m_linop->apply(vv_it1, *m_v_tmp_lhs);
303
304 gram_schmidt_orthogonalization(it);
305
306 auto tt = m_linop->norm2(vv_it1);
307
308 auto const sml = RT((sizeof(RT) == 8) ? 1.e-99 : 1.e-30);
309 bool happyend = (tt < sml);
310 if (!happyend) {
311 m_linop->scale(vv_it1, RT(1.0)/tt);
312 }
313
314 m_hh (it+1,it) = tt;
315 m_hes(it+1,it) = tt;
316
317 update_hessenberg(it, happyend, m_res);
318
319 ++it;
320 ++a_itcount;
321 a_status = converged(a_rnorm0, m_res) ? 0 : -1;
322 if (happyend) { break; }
323 }
324
325 if ((m_verbose > 1) && (a_status != 0 || a_itcount >= m_maxiter)) {
326 amrex::Print() << "GMRES: iter = " << a_itcount
327 << ", residual = " << m_res << ", " << m_res/a_rnorm0
328 << " (rel.)\n";
329 }
330
331 build_solution(a_xx, it-1);
332}
333
334template <typename V, typename M>
336{
337 // Two unmodified Gram-Schmidt Orthogonalization
338
339 BL_PROFILE("GMRES::GramSchmidt");
340
341 auto& vv_1 = m_vv[it+1];
342
343 Vector<RT> lhh(it+1);
344
345 for (int j = 0; j <= it; ++j) {
346 m_hh (j,it) = RT(0.0);
347 m_hes(j,it) = RT(0.0);
348 }
349
350 for (int ncnt = 0; ncnt < 2 ; ++ncnt)
351 {
352 for (int j = 0; j <= it; ++j) {
353 lhh[j] = m_linop->dotProduct(vv_1, m_vv[j]);
354 }
355
356 for (int j = 0; j <= it; ++j) {
357 m_linop->increment(vv_1, m_vv[j], -lhh[j]);
358 m_hh (j,it) += lhh[j];
359 m_hes(j,it) -= lhh[j];
360 }
361 }
362}
363
364template <typename V, typename M>
365void GMRES<V,M>::update_hessenberg (int it, bool happyend, RT& res)
366{
367 BL_PROFILE("GMRES::update_hessenberg()");
368
369 for (int j = 1; j <= it; ++j) {
370 auto tt = m_hh(j-1,it);
371 m_hh(j-1,it) = m_cc[j-1] * tt + m_ss[j-1] * m_hh(j,it);
372 m_hh(j ,it) = m_cc[j-1] * m_hh(j,it) - m_ss[j-1] * tt;
373 }
374
375 if (!happyend)
376 {
377 auto tt = std::sqrt(m_hh(it,it)*m_hh(it,it) + m_hh(it+1,it)*m_hh(it+1,it));
378 m_cc[it] = m_hh(it ,it) / tt;
379 m_ss[it] = m_hh(it+1,it) / tt;
380 m_grs[it+1] = - (m_ss[it] * m_grs[it]);
381 m_grs[it ] = m_cc[it] * m_grs[it];
382 m_hh(it,it) = m_cc[it] * m_hh(it,it) + m_ss[it] * m_hh(it+1,it);
383 res = std::abs(m_grs[it+1]);
384 }
385 else
386 {
387 res = RT(0.0);
388 }
389}
390
391template <typename V, typename M>
392void GMRES<V,M>::build_solution (V& a_xx, int const it)
393{
394 BL_PROFILE("GMRES:build_solution()");
395
396 if (it < 0) { return; }
397
398 if (m_hh(it,it) != RT(0.0)) {
399 m_grs[it] /= m_hh(it,it);
400 } else {
401 m_grs[it] = RT(0.0);
402 }
403
404 for (int ii = 1; ii <= it; ++ii) {
405 int k = it - ii;
406 auto tt = m_grs[k];
407 for (int j = k+1; j <= it; ++j) {
408 tt -= m_hh(k,j) * m_grs[j];
409 }
410 m_grs[k] = tt / m_hh(k,k);
411 }
412
413 m_linop->setToZero(*m_v_tmp_rhs);
414 for (int ii = 0; ii < it+1; ++ii) {
415 m_linop->increment(*m_v_tmp_rhs, m_vv[ii], m_grs[ii]);
416 }
417
418 m_linop->precond(*m_v_tmp_lhs, *m_v_tmp_rhs);
419 m_linop->increment(a_xx, *m_v_tmp_lhs, RT(1.0));
420}
421
422template <typename V, typename M>
423void GMRES<V,M>::compute_residual (V& a_rr, V const& a_xx, V const& a_bb)
424{
425 BL_PROFILE("GMRES::compute_residual()");
426 m_linop->assign(*m_v_tmp_lhs, a_xx);
427 m_linop->apply(*m_v_tmp_rhs, *m_v_tmp_lhs);
428 m_linop->linComb(a_rr, RT(1.0), a_bb, RT(-1.0), *m_v_tmp_rhs);
429}
430
431}
432#endif
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#define AMREX_ALWAYS_ASSERT(EX)
Definition AMReX_BLassert.H:50
GMRES.
Definition AMReX_GMRES.H:80
int getNumIters() const
Gets the number of iterations.
Definition AMReX_GMRES.H:113
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:212
int m_status
Definition AMReX_GMRES.H:136
void build_solution(V &a_xx, int it)
Definition AMReX_GMRES.H:392
void allocate_scratch()
Definition AMReX_GMRES.H:161
int m_verbose
Definition AMReX_GMRES.H:133
void define(M &linop)
Definition AMReX_GMRES.H:187
RT m_rtol
Definition AMReX_GMRES.H:139
bool converged(RT r0, RT r) const
Definition AMReX_GMRES.H:206
Vector< RT > m_hh_1d
Definition AMReX_GMRES.H:141
RT m_atol
Definition AMReX_GMRES.H:140
int m_maxiter
Definition AMReX_GMRES.H:134
void cycle(V &a_xx, int &a_status, int &a_itcount, RT &a_rnorm0)
Definition AMReX_GMRES.H:265
Table2D< RT > m_hes
Definition AMReX_GMRES.H:144
Table2D< RT > m_hh
Definition AMReX_GMRES.H:143
void setVerbose(int v)
Sets verbosity.
Definition AMReX_GMRES.H:104
void compute_residual(V &a_rr, V const &a_xx, V const &a_bb)
Definition AMReX_GMRES.H:423
Vector< RT > m_cc
Definition AMReX_GMRES.H:146
M * m_linop
Definition AMReX_GMRES.H:151
int m_its
Definition AMReX_GMRES.H:135
int m_restrtlen
Definition AMReX_GMRES.H:137
void setMaxIters(int niters)
Sets the max number of iterations.
Definition AMReX_GMRES.H:110
std::unique_ptr< V > m_v_tmp_lhs
Definition AMReX_GMRES.H:149
Vector< RT > m_hes_1d
Definition AMReX_GMRES.H:142
void setRestartLength(int rl)
Sets restart length. The default is 30.
Definition AMReX_GMRES.H:177
void update_hessenberg(int it, bool happyend, RT &res)
Definition AMReX_GMRES.H:365
void gram_schmidt_orthogonalization(int it)
Definition AMReX_GMRES.H:335
RT getResidualNorm() const
Gets the 2-norm of the residual.
Definition AMReX_GMRES.H:119
GMRES()
Definition AMReX_GMRES.H:155
Vector< V > m_vv
Definition AMReX_GMRES.H:150
Vector< RT > m_grs
Definition AMReX_GMRES.H:145
Vector< RT > m_ss
Definition AMReX_GMRES.H:147
void clear()
Definition AMReX_GMRES.H:194
int getStatus() const
Gets the solver status.
Definition AMReX_GMRES.H:116
std::unique_ptr< V > m_v_tmp_rhs
Definition AMReX_GMRES.H:148
RT m_res
Definition AMReX_GMRES.H:138
typename M::RT RT
Definition AMReX_GMRES.H:83
This class provides the user with a few print options.
Definition AMReX_Print.H:35
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:27
Definition AMReX_Amr.cpp:49
double second() noexcept
Definition AMReX_Utility.cpp:922
Definition AMReX_TableData.H:93