Add GMRES(m) linear solver — scalar, AVX2, NEON, OMP backends (ROADMA…#198
Merged
Conversation
Contributor
There was a problem hiding this comment.
Pull request overview
This PR adds a restarted GMRES(m) linear solver to the Poisson solver suite, including scalar, OpenMP, and SIMD (AVX2/NEON) implementations, along with new unit/consistency tests and documentation updates. This expands the library’s Krylov solver coverage to include a method suitable for future non-symmetric operators while cross-validating against existing CG behavior on the current SPD Poisson operator.
Changes:
- Implement restarted GMRES(m) across scalar, OpenMP, and SIMD (AVX2/NEON) backends and wire it into solver creation/dispatch.
- Add new GMRES-focused unit tests plus cross-backend consistency tests (scalar vs SIMD, scalar vs OMP), and hook them into CTest.
- Update public API/docs/roadmap to reflect the new method and the new
params.restartknob.
Reviewed changes
Copilot reviewed 18 out of 18 changed files in this pull request and generated 2 comments.
Show a summary per file
| File | Description |
|---|---|
| tests/math/test_omp_consistency.c | Adds OMP vs scalar GMRES consistency test and registers it in the existing OMP consistency test runner. |
| tests/math/test_gmres.c | New comprehensive GMRES(m) unit test suite (creation/init, convergence, CG cross-checks, restart behavior, preconditioning seam, error handling). |
| tests/math/test_gmres_neon.c | New NEON-vs-scalar numerical consistency test for GMRES SIMD backend. |
| tests/math/test_gmres_avx2.c | New AVX2-vs-scalar numerical consistency test for GMRES SIMD backend. |
| ROADMAP.md | Updates solver backend matrix and roadmap text to mark GMRES(m) implemented and clarify GPU SOR status. |
| lib/src/solvers/linear/simd/linear_solver_simd_dispatch.c | Extends SIMD dispatcher to route GMRES to AVX2 or NEON implementations at runtime. |
| lib/src/solvers/linear/simd_template/linear_solver_gmres_simd_template.h | Introduces the parameterized SIMD GMRES(m) implementation used by AVX2 and NEON backends. |
| lib/src/solvers/linear/omp/linear_solver_gmres_omp.c | Adds OpenMP-parallel GMRES(m) backend (parallel O(n) primitives; dense pieces remain scalar). |
| lib/src/solvers/linear/neon/linear_solver_gmres_neon.c | Adds NEON macro wrapper and includes the GMRES SIMD template (with OpenMP requirement). |
| lib/src/solvers/linear/avx2/linear_solver_gmres_avx2.c | Adds AVX2 macro wrapper and includes the GMRES SIMD template (with OpenMP requirement). |
| lib/src/solvers/linear/cpu/linear_solver_gmres.c | Adds scalar reference GMRES(m) implementation (right-preconditioned Jacobi seam, restart loop). |
| lib/src/solvers/linear/linear_solver.c | Adds params.restart defaulting and GMRES factory wiring in poisson_solver_create(). |
| lib/src/solvers/linear/linear_solver_internal.h | Declares GMRES factories and introduces GMRES constants (default restart, breakdown threshold). |
| lib/include/cfd/solvers/poisson_solver.h | Extends public API with POISSON_METHOD_GMRES, params.restart, and GMRES solver type strings. |
| lib/CMakeLists.txt | Adds GMRES sources to scalar/SIMD/OMP build source lists. |
| docs/reference/solvers.md | Documents GMRES(m) usage, parameters, and supported backends. |
| docs/reference/api-reference.md | Updates reference docs to include GMRES method and the restart parameter. |
| CMakeLists.txt | Adds GMRES tests to the build and registers them with CTest. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
shaia
force-pushed
the
feature/gmres-linear-solver
branch
from
d04469c to
70f7c65
Compare
…P §1.2) Restarted GMRES(m) completes the Krylov solver suite and provides a non-symmetric-capable solver for future implicit integrators (§1.5). On the current SPD Poisson operator it cross-validates against CG to within 1e-4 max-diff. Algorithm choices that prevented silent bugs: - Right preconditioning so the Hessenberg residual estimate equals the true unpreconditioned residual, keeping stats.final_residual and the CG cross-check consistent. - Convergence decided on the recomputed true residual after each restart, not the cheap Givens estimate, which drifts below the real residual on finer grids due to loss of MGS orthogonality. - Dense O(m²) Givens/Hessenberg/back-sub kept as plain scalar in every backend (AVX2, NEON, OMP) for byte-consistent results across backends. New sources: linear_solver_gmres.c (scalar), linear_solver_gmres_simd_template.h (vectorizes only O(n) primitives), linear_solver_gmres_avx2.c, _neon.c, _omp.c. New tests: test_gmres.c (12 cases), test_gmres_avx2.c, test_gmres_neon.c, plus GMRES OMP case in test_omp_consistency.c. All 75 fast-suite tests pass. GPU backend deferred (noted in ROADMAP).
…norm-for-norm stats match The docstring claimed the test validates stats.final_residual against a recomputed true residual, but the independent check uses the max-norm while the reported residual is a 2-norm — so it verifies smallness, not an exact norm match. Reword to match what the test actually asserts.
The test allows an RMS L2 difference up to 1e-9 (parallel-reduction rounding), so the solutions are numerically consistent, not bit-for-bit identical. Fix the docstring to reflect the actual tolerance-bounded assertion.
The Poisson methods snippet used a nonexistent type name (poisson_method_t), fabricated fixed numeric values, listed a POISSON_METHOD_PCG that does not exist, and omitted GAUSS_SEIDEL and MULTIGRID. Mirror the real enum from poisson_solver.h and note that preconditioning is a separate params field, so consumers get an accurate, copy-pasteable reference.
verbose is a bool (not int), the preconditioner field type is poisson_precond_type_t (not poisson_precond_t), and the omega default is 0 (auto-optimal), not 1.5. Mirror poisson_solver.h so the struct is copy-pasteable into real code.
shaia
force-pushed
the
feature/gmres-linear-solver
branch
from
70f7c65 to
b21d0fe
Compare
…r-solver # Conflicts: # ROADMAP.md
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
…P §1.2)
Restarted GMRES(m) completes the Krylov solver suite and provides a non-symmetric-capable solver for future implicit integrators (§1.5). On the current SPD Poisson operator it cross-validates against CG to within 1e-4 max-diff.
Algorithm choices that prevented silent bugs:
New sources: linear_solver_gmres.c (scalar), linear_solver_gmres_simd_template.h (vectorizes only O(n) primitives), linear_solver_gmres_avx2.c, _neon.c, _omp.c. New tests: test_gmres.c (12 cases), test_gmres_avx2.c, test_gmres_neon.c, plus GMRES OMP case in test_omp_consistency.c. All 75 fast-suite tests pass. GPU backend deferred (noted in ROADMAP).