Use modified Gram-Schmidt orthogonalization in GMRES - #600
Use modified Gram-Schmidt orthogonalization in GMRES#600michelebucelli wants to merge 2 commits into
Conversation
There was a problem hiding this comment.
Pull request overview
This pull request updates the GMRES linear solver to use modified Gram–Schmidt (MGS) orthogonalization (extracted into dedicated helper functions) to improve numerical stability and reduce divergence between MPI runs, and also performs a small cleanup in Fourier interpolation evaluation.
Changes:
- Replaced classical Gram–Schmidt with modified Gram–Schmidt in GMRES and refactored the orthogonalization step into shared helper functions.
- Added a “lucky breakdown” handling path based on comparing the subdiagonal Hessenberg entry to machine epsilon-scaled norms.
- Simplified arithmetic in
FourierInterpolationby using compound assignment operators (-=/+=).
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 2 comments.
| File | Description |
|---|---|
Code/Source/linear_solver/gmres.cpp |
Introduces MGS orthogonalization helpers and replaces the previous orthogonalization logic in GMRES variants. |
Code/Source/solver/FourierInterpolation.cpp |
Uses compound assignment operators for trend subtraction and Fourier series accumulation. |
Suppressed comments (2)
Code/Source/linear_solver/gmres.cpp:426
- The new lucky-breakdown handling in orthogonalize_* can set h(i+1,i)=0.0, which makes the subsequent Givens-rotation normalization
tmp = sqrt(h(i,i)^2 + h(i+1,i)^2)potentially become 0.0. That leads to division by zero when computing c(i) and s(i) and can produce NaNs.
orthogonalize_s(lhs, nNo, mynNo, i, u, h);
Code/Source/linear_solver/gmres.cpp:598
- The new lucky-breakdown handling in orthogonalize_* can set h(i+1,i)=0.0, which makes the subsequent Givens-rotation normalization
tmp = sqrt(h(i,i)^2 + h(i+1,i)^2)potentially become 0.0. That leads to division by zero when computing c(i) and s(i) and can produce NaNs.
orthogonalize_v(lhs, dof, nNo, mynNo, i, u, h);
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| for (int j = 0; j <= i; j++) { | ||
| h(j, i) = dot::fsils_dot_v(dof, mynNo, lhs.commu, u.rslice(j), w); | ||
| omp_la::omp_sum_v(dof, nNo, -h(j, i), w, u.rslice(j)); |
|
|
||
| u_slice_1 = u.rslice(i+1); | ||
| omp_la::omp_mul_v(dof, nNo, 1.0/h(i+1,i), u_slice_1); | ||
| orthogonalize_v(lhs, dof, nNo, mynNo, i, u, h); |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #600 +/- ##
==========================================
- Coverage 72.54% 72.54% -0.01%
==========================================
Files 252 252
Lines 39033 39023 -10
Branches 6684 6682 -2
==========================================
- Hits 28318 28308 -10
+ Misses 10480 10479 -1
- Partials 235 236 +1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
I would like to see some tests on real problems to collect data on performance changes, improved convergence, etc. |
|
Moving this to draft, pending the discussion at #599. |
Fixes #599.
Current situation
The GMRES method, as implemented in
gmres.cpp, used unmodified Gram-Schmidt orthogonalization. Multiple experiments indicate that the numerical instability related to that propagates to other parts of the solver. See #599 for additional details.Release Notes
gmresandgmres_v);@todoabout replacingvalue = value + ...withvalue += ...inFourierInterpolation.Documentation
The new functions that wrap the orthogonalization step have Doxygen documentation.
None of these changes is user-facing.
Testing
All automatic tests are still passing.
Code of Conduct & Contributing Guidelines