Add GMRES(m) linear solver — scalar, AVX2, NEON, OMP backends (ROADMA…

CMakeLists.txt

@@ -326,6 +326,9 @@ if(BUILD_TESTS)
     add_executable(test_bicgstab tests/math/test_bicgstab.c)
     add_executable(test_bicgstab_avx2 tests/math/test_bicgstab_avx2.c)
     add_executable(test_bicgstab_neon tests/math/test_bicgstab_neon.c)
+    add_executable(test_gmres tests/math/test_gmres.c)
+    add_executable(test_gmres_avx2 tests/math/test_gmres_avx2.c)
+    add_executable(test_gmres_neon tests/math/test_gmres_neon.c)
     add_executable(test_omp_consistency tests/math/test_omp_consistency.c)
     add_executable(test_optimal_omega tests/math/test_optimal_omega.c)
     add_executable(test_solver_breakdown tests/math/test_solver_breakdown.c)
@@ -439,6 +442,9 @@ if(BUILD_TESTS)
     target_link_libraries(test_bicgstab PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_bicgstab_avx2 PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_bicgstab_neon PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
+    target_link_libraries(test_gmres PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
+    target_link_libraries(test_gmres_avx2 PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
+    target_link_libraries(test_gmres_neon PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_omp_consistency PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_optimal_omega PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_solver_breakdown PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
@@ -586,7 +592,11 @@ if(BUILD_TESTS)
     add_test(NAME BiCGSTABTest COMMAND test_bicgstab)
     add_test(NAME BiCGSTAB_AVX2_Consistency COMMAND test_bicgstab_avx2)
     add_test(NAME BiCGSTAB_NEON_Consistency COMMAND test_bicgstab_neon)
+    add_test(NAME GMRESTest COMMAND test_gmres)
+    add_test(NAME GMRES_AVX2_Consistency COMMAND test_gmres_avx2)
+    add_test(NAME GMRES_NEON_Consistency COMMAND test_gmres_neon)
     set_tests_properties(BiCGSTAB_NEON_Consistency PROPERTIES LABELS "cross-arch")
+    set_tests_properties(GMRES_NEON_Consistency PROPERTIES LABELS "cross-arch")
     add_test(NAME OMPConsistencyTest COMMAND test_omp_consistency)
     add_test(NAME OptimalOmegaTest COMMAND test_optimal_omega)
     add_test(NAME SolverBreakdownTest COMMAND test_solver_breakdown)

ROADMAP.md

@@ -52,6 +52,7 @@ The single source of truth for backend gaps. Each algorithm targets scalar (CPU)
 |                     | Red-Black SOR  | done | done     | done     | done     | done |
 |                     | CG / PCG       | done | done     | done     | done     | done |
 |                     | BiCGSTAB       | done | done     | done     | —        | done |
+|                     | GMRES(m)       | done | done     | done     | done     | —    |
 | **Boundary Conds**  | All types      | done | done     | done     | done     | done |
 
 ### Known Limitations
@@ -103,15 +104,22 @@ No-slip, Inlet, Outlet, Symmetry, Moving wall, Time-varying). See CHANGELOG.
 
 ### 1.2 Linear Solvers (P0)
 
-Implemented: Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB (all with SIMD backends; CG is the
-default Poisson solver for projection methods). GPU standalone Jacobi, CG, Red-Black SOR,
-plain SOR, and BiCGSTAB are done and validated vs CPU; `solve_projection_method_gpu` uses
+Implemented: Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB, GMRES(m) (all with SIMD backends;
+CG is the default Poisson solver for projection methods). GPU standalone Jacobi, CG, Red-Black
+SOR, plain SOR, and BiCGSTAB are done and validated vs CPU; `solve_projection_method_gpu` uses
 on-device CG.
 
+GMRES(m) with restart is implemented across scalar, AVX2, NEON, and OMP backends (right-
+preconditioned Jacobi seam; validated vs CG on the SPD Poisson problem, with restart-no-stall
+and cross-backend consistency tests). Since the current linear systems are all the symmetric
+pressure-Poisson operator, GMRES is a forward-looking addition (it becomes load-bearing once
+non-symmetric operators arrive, e.g. implicit advection-diffusion in §1.5).
+
 **Still needed:**
 
-- [ ] GMRES (Generalized Minimal Residual) for non-symmetric systems — scalar, AVX2, NEON,
-      OMP, GPU
+- [x] GMRES (Generalized Minimal Residual) for non-symmetric systems — scalar, AVX2, NEON, OMP
+      (done). GPU variant deferred.
+- [ ] GMRES GPU backend (deferred from the initial GMRES landing)
 - [ ] SSOR (Symmetric SOR) preconditioner
 - [ ] ILU preconditioner
 - [ ] Geometric multigrid

docs/reference/api-reference.md

@@ -408,15 +408,20 @@ void poisson_solver_destroy(poisson_solver_t* solver);
 
 ```c
 typedef enum {
-    POISSON_METHOD_JACOBI = 0,
-    POISSON_METHOD_SOR = 1,
-    POISSON_METHOD_REDBLACK_SOR = 2,
-    POISSON_METHOD_CG = 3,
-    POISSON_METHOD_PCG = 4,
-    POISSON_METHOD_BICGSTAB = 5,
-} poisson_method_t;
+    POISSON_METHOD_JACOBI,        // Jacobi iteration (fully parallelizable)
+    POISSON_METHOD_GAUSS_SEIDEL,  // Gauss-Seidel iteration
+    POISSON_METHOD_SOR,           // Successive Over-Relaxation
+    POISSON_METHOD_REDBLACK_SOR,  // Red-Black SOR (parallelizable)
+    POISSON_METHOD_CG,            // Conjugate Gradient (SPD systems)
+    POISSON_METHOD_BICGSTAB,      // BiCGSTAB (non-symmetric systems)
+    POISSON_METHOD_GMRES,         // Restarted GMRES(m) (scalar/SIMD/OMP)
+    POISSON_METHOD_MULTIGRID      // Multigrid (future)
+} poisson_solver_method_t;
 ```
 
+> Preconditioning is a separate `preconditioner` field on `poisson_solver_params_t`
+> (see below), not a distinct method — CG becomes PCG when a preconditioner is set.
+
 ### Poisson Backends
 
 ```c
@@ -436,10 +441,11 @@ typedef struct {
     double tolerance;           // Relative tolerance (default: 1e-6)
     double absolute_tolerance;  // Absolute tolerance (default: 1e-10)
     int max_iterations;         // Max iterations (default: 5000)
-    double omega;               // SOR relaxation (default: 1.5)
+    double omega;               // SOR relaxation (default: 0 = auto-optimal)
     int check_interval;         // Convergence check interval (default: 1)
-    int verbose;                // Print convergence info (default: 0)
-    poisson_precond_t preconditioner;  // Preconditioner (default: NONE)
+    bool verbose;               // Print convergence info (default: false)
+    poisson_precond_type_t preconditioner;  // Preconditioner (default: POISSON_PRECOND_NONE)
+    int restart;                // GMRES(m) restart length (default: 0 = auto/30)
 } poisson_solver_params_t;
 
 poisson_solver_params_t poisson_solver_params_default(void);

docs/reference/solvers.md

@@ -226,6 +226,35 @@ poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_BICGSTAB,
                                                  POISSON_BACKEND_SCALAR);
 ```
 
+#### 7. GMRES(m)
+
+**Algorithm:** Restarted Generalized Minimal Residual for non-symmetric systems.
+Builds an orthonormal Krylov basis via Arnoldi + modified Gram-Schmidt and
+minimizes the residual over that basis using incremental Givens rotations,
+restarting every `m` inner iterations to bound memory.
+
+**Characteristics:**
+- Handles non-symmetric matrices; minimizes the residual norm at every step
+- Restart length `m` set via `params.restart` (0 = auto, default 30); bounds the
+  Krylov basis to `m+1` grid-sized vectors
+- Right-preconditioned Jacobi seam (`POISSON_PRECOND_JACOBI`); like PCG, Jacobi
+  preconditioning gives no benefit on a uniform grid (constant diagonal)
+- On the symmetric Poisson operator it converges to the same solution as CG (a
+  useful independent cross-check); its real value is future non-symmetric systems
+
+**Backends:** scalar, SIMD (AVX2/NEON), and OpenMP. A CUDA GPU variant is not yet
+implemented.
+
+**Usage:**
+```c
+poisson_solver_params_t params = poisson_solver_params_default();
+params.restart = 30;  // GMRES(30); 0 selects the default
+
+poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_GMRES,
+                                                 POISSON_BACKEND_SCALAR);
+poisson_solver_init(solver, nx, ny, nz, dx, dy, dz, &params);
+```
+
 ### Linear Solver Performance Comparison
 
 **Problem:** 65×65 grid, tolerance = 1e-6

lib/CMakeLists.txt

@@ -128,6 +128,7 @@ set(CFD_SCALAR_SOURCES
     src/solvers/linear/cpu/linear_solver_redblack.c
     src/solvers/linear/cpu/linear_solver_cg.c
     src/solvers/linear/cpu/linear_solver_bicgstab.c
+    src/solvers/linear/cpu/linear_solver_gmres.c
 )
 
 # SIMD sources (AVX2/NEON optimized)
@@ -152,11 +153,13 @@ set(CFD_SIMD_SOURCES
     src/solvers/linear/avx2/linear_solver_cg_avx2.c
     src/solvers/linear/avx2/linear_solver_bicgstab_avx2.c
     src/solvers/linear/avx2/linear_solver_sor_avx2.c
+    src/solvers/linear/avx2/linear_solver_gmres_avx2.c
     src/solvers/linear/neon/linear_solver_jacobi_neon.c
     src/solvers/linear/neon/linear_solver_redblack_neon.c
     src/solvers/linear/neon/linear_solver_cg_neon.c
     src/solvers/linear/neon/linear_solver_bicgstab_neon.c
     src/solvers/linear/neon/linear_solver_sor_neon.c
+    src/solvers/linear/neon/linear_solver_gmres_neon.c
 )
 
 # OpenMP sources
@@ -174,6 +177,7 @@ set(CFD_OMP_SOURCES
     # Linear solvers - OMP
     src/solvers/linear/omp/linear_solver_redblack_omp.c
     src/solvers/linear/omp/linear_solver_cg_omp.c
+    src/solvers/linear/omp/linear_solver_gmres_omp.c
 )
 
 # Note: API sources (solver_registry, simulation_api, output_registry) are in

lib/include/cfd/solvers/poisson_solver.h

@@ -56,7 +56,8 @@ typedef enum {
     POISSON_METHOD_SOR,           /**< Successive Over-Relaxation */
     POISSON_METHOD_REDBLACK_SOR,  /**< Red-Black SOR (parallelizable) */
     POISSON_METHOD_CG,            /**< Conjugate Gradient (for SPD systems) */
-    POISSON_METHOD_BICGSTAB,      /**< BiCGSTAB (future) */
+    POISSON_METHOD_BICGSTAB,      /**< BiCGSTAB (for non-symmetric systems) */
+    POISSON_METHOD_GMRES,         /**< Restarted GMRES(m) (for non-symmetric systems) */
     POISSON_METHOD_MULTIGRID      /**< Multigrid (future) */
 } poisson_solver_method_t;
 
@@ -105,6 +106,7 @@ typedef struct {
     int check_interval;        /**< Check convergence every N iterations (default: 1) */
     bool verbose;              /**< Print iteration progress (default: false) */
     poisson_precond_type_t preconditioner; /**< Preconditioner type (default: NONE) */
+    int restart;               /**< GMRES restart length m (default: 0 = auto/30); ignored by other methods */
 } poisson_solver_params_t;
 
 /**
@@ -391,6 +393,9 @@ CFD_LIBRARY_EXPORT bool poisson_solver_backend_available(poisson_solver_backend_
 #define POISSON_SOLVER_TYPE_BICGSTAB_SCALAR   "bicgstab_scalar"
 #define POISSON_SOLVER_TYPE_BICGSTAB_SIMD     "bicgstab_simd"
 #define POISSON_SOLVER_TYPE_BICGSTAB_GPU      "bicgstab_gpu"
+#define POISSON_SOLVER_TYPE_GMRES_SCALAR      "gmres_scalar"
+#define POISSON_SOLVER_TYPE_GMRES_OMP         "gmres_omp"
+#define POISSON_SOLVER_TYPE_GMRES_SIMD        "gmres_simd"
 
 /* ============================================================================
  * CONVENIENCE API

lib/src/solvers/linear/avx2/linear_solver_gmres_avx2.c

@@ -0,0 +1,76 @@
+/**
+ * @file linear_solver_gmres_avx2.c
+ * @brief Restarted GMRES(m) linear solver - AVX2 + OpenMP backend
+ *
+ * Provides AVX2-specific macro definitions and includes the parameterized SIMD
+ * template. All algorithm logic lives in the template header. Only the O(n)
+ * vector primitives are vectorized; dense Givens/Hessenberg work stays scalar.
+ */
+
+#include "../linear_solver_internal.h"
+#include "cfd/boundary/boundary_conditions.h"
+#include "cfd/core/memory.h"
+#include <math.h>
+#include <string.h>
+#include <limits.h>
+
+#if defined(CFD_HAS_AVX2) && defined(CFD_ENABLE_OPENMP)
+#define GMRES_HAS_AVX2 1
+#include <immintrin.h>
+#include <omp.h>
+#endif
+
+#if defined(GMRES_HAS_AVX2)
+
+//=============================================================================
+// AVX2-SPECIFIC MACRO DEFINITIONS
+//=============================================================================
+
+#define SIMD_SUFFIX avx2
+#define SIMD_VEC __m256d
+#define SIMD_WIDTH 4
+
+#define SIMD_LOAD(ptr) _mm256_loadu_pd(ptr)
+#define SIMD_STORE(ptr, vec) _mm256_storeu_pd(ptr, vec)
+#define SIMD_SET1(val) _mm256_set1_pd(val)
+#define SIMD_SETZERO() _mm256_setzero_pd()
+
+#define SIMD_ADD(a, b) _mm256_add_pd(a, b)
+#define SIMD_SUB(a, b) _mm256_sub_pd(a, b)
+#define SIMD_MUL(a, b) _mm256_mul_pd(a, b)
+#define SIMD_FMA(a, b, c) _mm256_fmadd_pd(a, b, c)
+
+/* AVX2 horizontal sum (4 doubles -> 1 double) - MSVC-compatible */
+static inline double gmres_simd_hsum_avx2(__m256d vec) {
+    __m128d low = _mm256_castpd256_pd128(vec);
+    __m128d high = _mm256_extractf128_pd(vec, 1);
+    __m128d sum128 = _mm_add_pd(low, high);
+    sum128 = _mm_hadd_pd(sum128, sum128);
+    return _mm_cvtsd_f64(sum128);
+}
+
+#define SIMD_HSUM(vec) gmres_simd_hsum_avx2(vec)
+
+#include "../simd_template/linear_solver_gmres_simd_template.h"
+
+#undef SIMD_SUFFIX
+#undef SIMD_VEC
+#undef SIMD_WIDTH
+#undef SIMD_LOAD
+#undef SIMD_STORE
+#undef SIMD_SET1
+#undef SIMD_SETZERO
+#undef SIMD_ADD
+#undef SIMD_SUB
+#undef SIMD_MUL
+#undef SIMD_FMA
+#undef SIMD_HSUM
+
+#else  /* !GMRES_HAS_AVX2 */
+
+/* Stub for platforms without AVX2 */
+poisson_solver_t* create_gmres_avx2_solver(void) {
+    return NULL;
+}
+
+#endif  /* GMRES_HAS_AVX2 */