From 223d61a7e41cf156824c084c5a44a5829540fc3d Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 4 Jul 2026 15:53:04 +0300 Subject: [PATCH 1/5] =?UTF-8?q?Add=20GMRES(m)=20linear=20solver=20?= =?UTF-8?q?=E2=80=94=20scalar,=20AVX2,=20NEON,=20OMP=20backends=20(ROADMAP?= =?UTF-8?q?=20=C2=A71.2)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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). --- CMakeLists.txt | 10 + ROADMAP.md | 18 +- docs/reference/api-reference.md | 2 + docs/reference/solvers.md | 29 + lib/CMakeLists.txt | 4 + lib/include/cfd/solvers/poisson_solver.h | 7 +- .../linear/avx2/linear_solver_gmres_avx2.c | 76 ++ .../solvers/linear/cpu/linear_solver_gmres.c | 576 +++++++++++++++ lib/src/solvers/linear/linear_solver.c | 15 + .../solvers/linear/linear_solver_internal.h | 27 + .../linear/neon/linear_solver_gmres_neon.c | 68 ++ .../linear/omp/linear_solver_gmres_omp.c | 519 ++++++++++++++ .../linear/simd/linear_solver_simd_dispatch.c | 23 + .../linear_solver_gmres_simd_template.h | 662 ++++++++++++++++++ tests/math/test_gmres.c | 632 +++++++++++++++++ tests/math/test_gmres_avx2.c | 186 +++++ tests/math/test_gmres_neon.c | 183 +++++ tests/math/test_omp_consistency.c | 99 +++ 18 files changed, 3130 insertions(+), 6 deletions(-) create mode 100644 lib/src/solvers/linear/avx2/linear_solver_gmres_avx2.c create mode 100644 lib/src/solvers/linear/cpu/linear_solver_gmres.c create mode 100644 lib/src/solvers/linear/neon/linear_solver_gmres_neon.c create mode 100644 lib/src/solvers/linear/omp/linear_solver_gmres_omp.c create mode 100644 lib/src/solvers/linear/simd_template/linear_solver_gmres_simd_template.h create mode 100644 tests/math/test_gmres.c create mode 100644 tests/math/test_gmres_avx2.c create mode 100644 tests/math/test_gmres_neon.c diff --git a/CMakeLists.txt b/CMakeLists.txt index c70c56b6..4899c727 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -325,6 +325,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) @@ -437,6 +440,9 @@ if(BUILD_TESTS) target_link_libraries(test_bicgstab PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_bicgstab_avx2 PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_bicgstab_neon PRIVATE CFD::Library unity $<$>:m>) + target_link_libraries(test_gmres PRIVATE CFD::Library unity $<$>:m>) + target_link_libraries(test_gmres_avx2 PRIVATE CFD::Library unity $<$>:m>) + target_link_libraries(test_gmres_neon PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_omp_consistency PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_optimal_omega PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_solver_breakdown PRIVATE CFD::Library unity $<$>:m>) @@ -583,7 +589,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) diff --git a/ROADMAP.md b/ROADMAP.md index 0e2abeef..4166a74f 100644 --- a/ROADMAP.md +++ b/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,14 +104,21 @@ 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, and -BiCGSTAB are done and validated vs CPU; `solve_projection_method_gpu` uses on-device CG. +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, 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 diff --git a/docs/reference/api-reference.md b/docs/reference/api-reference.md index 6f10c225..582791a9 100644 --- a/docs/reference/api-reference.md +++ b/docs/reference/api-reference.md @@ -414,6 +414,7 @@ typedef enum { POISSON_METHOD_CG = 3, POISSON_METHOD_PCG = 4, POISSON_METHOD_BICGSTAB = 5, + POISSON_METHOD_GMRES = 6, // Restarted GMRES(m) (scalar/SIMD/OMP) } poisson_method_t; ``` @@ -440,6 +441,7 @@ typedef struct { int check_interval; // Convergence check interval (default: 1) int verbose; // Print convergence info (default: 0) poisson_precond_t preconditioner; // Preconditioner (default: NONE) + int restart; // GMRES(m) restart length (default: 0 = auto/30) } poisson_solver_params_t; poisson_solver_params_t poisson_solver_params_default(void); diff --git a/docs/reference/solvers.md b/docs/reference/solvers.md index f523e0de..0c62caf6 100644 --- a/docs/reference/solvers.md +++ b/docs/reference/solvers.md @@ -224,6 +224,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, ¶ms); +``` + ### Linear Solver Performance Comparison **Problem:** 65×65 grid, tolerance = 1e-6 diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 153f05bd..2e7a3409 100644 --- a/lib/CMakeLists.txt +++ b/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 diff --git a/lib/include/cfd/solvers/poisson_solver.h b/lib/include/cfd/solvers/poisson_solver.h index 93e6f4a3..b7887969 100644 --- a/lib/include/cfd/solvers/poisson_solver.h +++ b/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 diff --git a/lib/src/solvers/linear/avx2/linear_solver_gmres_avx2.c b/lib/src/solvers/linear/avx2/linear_solver_gmres_avx2.c new file mode 100644 index 00000000..d9daa13d --- /dev/null +++ b/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 +#include +#include + +#if defined(CFD_HAS_AVX2) && defined(CFD_ENABLE_OPENMP) +#define GMRES_HAS_AVX2 1 +#include +#include +#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 */ diff --git a/lib/src/solvers/linear/cpu/linear_solver_gmres.c b/lib/src/solvers/linear/cpu/linear_solver_gmres.c new file mode 100644 index 00000000..650cafe2 --- /dev/null +++ b/lib/src/solvers/linear/cpu/linear_solver_gmres.c @@ -0,0 +1,576 @@ +/** + * @file linear_solver_gmres.c + * @brief Restarted GMRES(m) solver - scalar CPU implementation + * + * Generalized Minimal Residual method for solving Ax = b. Unlike CG, GMRES does + * not require A to be symmetric positive definite, so it is the Krylov method of + * choice for non-symmetric systems (e.g. implicit advection-diffusion). For the + * pressure-Poisson operator used here (A = -nabla^2, SPD) GMRES converges to the + * same solution as CG and serves as an independent cross-check. + * + * GMRES characteristics: + * - Minimizes the residual norm over the Krylov subspace at every inner step. + * - Restarted GMRES(m) caps the basis at m+1 vectors to bound memory; the outer + * loop restarts from the current iterate until convergence or max_iterations. + * - Each inner iteration: 1 matrix-vector product, (j+1) dot products + axpys for + * modified Gram-Schmidt, plus O(m) scalar work for Givens rotations. + * + * Algorithm (restarted GMRES(m), right-preconditioned when enabled): + * outer restart: + * r = b - A x; beta = ||r||; v_0 = r / beta; g = (beta, 0, ..., 0) + * inner j = 0..m-1 (Arnoldi + modified Gram-Schmidt): + * w = A M^{-1} v_j (w = A v_j when unpreconditioned) + * for i <= j: H[i,j] = (w, v_i); w -= H[i,j] v_i + * H[j+1,j] = ||w||; v_{j+1} = w / H[j+1,j] + * apply stored Givens rotations to column j; build a new one to zero + * H[j+1,j]; rotate g. Residual estimate = |g[j+1]|. + * solve upper-triangular H y = g; x += M^{-1}(V y) (x += V y unpreconditioned) + * + * Right (not left) preconditioning is used so the cheap Givens estimate |g[k]| + * equals the true unpreconditioned residual ||b - A x||, keeping the convergence + * check, reported statistics, and the CG cross-check consistent. See the plan and + * the PCG notes: on a uniform grid the Laplacian diagonal is constant, so Jacobi + * preconditioning is a scalar multiply and does not reduce iteration count. + */ + +#include "../linear_solver_internal.h" + +#include "cfd/core/indexing.h" +#include "cfd/core/logging.h" +#include "cfd/core/memory.h" + +#include +#include + +/* ============================================================================ + * GMRES CONTEXT + * ============================================================================ */ + +typedef struct { + double dx2; /* dx^2 */ + double dy2; /* dy^2 */ + double inv_dz2; /* 1/dz^2 (0 for 2D) */ + double diag_inv; /* Jacobi preconditioner: 1/(2/dx2 + 2/dy2 + 2*inv_dz2) */ + + size_t stride_z; /* nx*ny for 3D, 0 for 2D */ + size_t k_start; /* first interior k index */ + size_t k_end; /* one-past-last interior k index */ + + int m; /* restart length (Krylov subspace dimension) */ + size_t n; /* total field size nx*ny*nz */ + + /* Large O(n) working vectors (allocated during init) */ + double* Vblock; /* (m+1) Krylov basis vectors, V[j] = Vblock + j*n */ + double* w; /* A*v_j scratch */ + double* Mv; /* M^{-1}*v scratch (NULL if no precond) */ + + /* Small dense O(m^2)/O(m) arrays (scalar in every backend) */ + double* H; /* upper-Hessenberg (m+1) x m, column-major H[i + j*(m+1)] */ + double* cs; /* Givens cosines, length m */ + double* sn; /* Givens sines, length m */ + double* g; /* rotated least-squares RHS, length m+1 */ + double* y; /* back-substitution solution, length m */ + + int use_precond; /* Flag: is preconditioner enabled? */ + int initialized; +} gmres_context_t; + +/* ============================================================================ + * O(n) VECTOR PRIMITIVES (interior points only) + * + * dot_product / axpy / apply_laplacian / compute_residual / copy_vector / + * apply_jacobi_precond mirror linear_solver_cg.c exactly. scale_vector, + * vector_norm, and linear_combination are the GMRES-specific additions. + * ============================================================================ */ + +/** Dot product over interior points */ +static double dot_product(const double* a, const double* b, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + double sum = 0.0; + for (size_t k = k_start; k < k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = (k * stride_z) + (IDX_2D(i, j, nx)); + sum += a[idx] * b[idx]; + } + } + } + return sum; +} + +/** y = y + alpha * x (interior points only) */ +static void axpy(double alpha, const double* x, double* y, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + for (size_t k = k_start; k < k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = (k * stride_z) + (IDX_2D(i, j, nx)); + y[idx] += alpha * x[idx]; + } + } + } +} + +/** x = alpha * x (interior points only) */ +static void scale_vector(double alpha, double* x, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + for (size_t k = k_start; k < k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = (k * stride_z) + (IDX_2D(i, j, nx)); + x[idx] *= alpha; + } + } + } +} + +/** Euclidean norm over interior points: sqrt(dot(x,x)) */ +static double vector_norm(const double* x, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + return sqrt(dot_product(x, x, nx, ny, k_start, k_end, stride_z)); +} + +/** + * Apply negative Laplacian operator: Ap = -nabla^2(p) (matches CG's A = -nabla^2). + */ +static void apply_laplacian(const double* p, double* Ap, + size_t nx, size_t ny, + double dx2, double dy2, double inv_dz2, + size_t k_start, size_t k_end, size_t stride_z) { + double dx2_inv = 1.0 / dx2; + double dy2_inv = 1.0 / dy2; + + for (size_t k = k_start; k < k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = (k * stride_z) + (IDX_2D(i, j, nx)); + + double laplacian = + ((p[idx + 1] - (2.0 * p[idx]) + p[idx - 1]) * dx2_inv) + + ((p[idx + (nx)] - (2.0 * p[idx]) + p[idx - (nx)]) * dy2_inv) + + ((p[idx + (stride_z)] + p[idx - (stride_z)] - (2.0 * p[idx])) * inv_dz2); + Ap[idx] = -laplacian; + } + } + } +} + +/** + * Compute residual r = b - A x = -rhs + nabla^2 x (matches CG convention). + */ +static void compute_residual(const double* x, const double* rhs, double* r, + size_t nx, size_t ny, + double dx2, double dy2, double inv_dz2, + size_t k_start, size_t k_end, size_t stride_z) { + double dx2_inv = 1.0 / dx2; + double dy2_inv = 1.0 / dy2; + + for (size_t k = k_start; k < k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = (k * stride_z) + (IDX_2D(i, j, nx)); + + double laplacian = + ((x[idx + 1] - (2.0 * x[idx]) + x[idx - 1]) * dx2_inv) + + ((x[idx + (nx)] - (2.0 * x[idx]) + x[idx - (nx)]) * dy2_inv) + + ((x[idx + (stride_z)] + x[idx - (stride_z)] - (2.0 * x[idx])) * inv_dz2); + + r[idx] = -rhs[idx] + laplacian; + } + } + } +} + +/** Copy dst = src (interior points only) */ +static void copy_vector(const double* src, double* dst, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + for (size_t k = k_start; k < k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = (k * stride_z) + (IDX_2D(i, j, nx)); + dst[idx] = src[idx]; + } + } + } +} + +/** Apply Jacobi preconditioner: z = M^{-1} r = diag_inv * r (interior points) */ +static void apply_jacobi_precond(const double* r, double* z, + size_t nx, size_t ny, + double diag_inv, + size_t k_start, size_t k_end, size_t stride_z) { + for (size_t k = k_start; k < k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = (k * stride_z) + (IDX_2D(i, j, nx)); + z[idx] = diag_inv * r[idx]; + } + } + } +} + +/** + * Accumulate x += sum_{j=0..k-1} y[j] * V[j] over interior points. + * Simplest correct form is k axpy calls; V holds the Krylov basis. + */ +static void linear_combination(double* x, const double* Vblock, size_t n, + const double* y, int k, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + for (int jj = 0; jj < k; jj++) { + axpy(y[jj], Vblock + (size_t)jj * n, x, nx, ny, k_start, k_end, stride_z); + } +} + +/* ============================================================================ + * DENSE SCALAR HELPERS (small O(m^2)/O(m) work — scalar in every backend) + * ============================================================================ */ + +/** + * Apply previously stored Givens rotations to column j of H, then build and apply + * a new rotation that zeros H[j+1,j], and rotate the RHS vector g accordingly. + * H is column-major, (m+1) x m: element (row i, col j) is H[i + j*(m+1)]. + * After the call, |g[j+1]| is the GMRES residual estimate for the current basis. + */ +static void apply_givens_and_update_hessenberg(double* H, double* cs, double* sn, + double* g, int j, int m) { + const int lead = m + 1; /* column stride */ + + /* Apply rotations 0..j-1 to the new column j */ + for (int i = 0; i < j; i++) { + double temp = cs[i] * H[i + j * lead] + sn[i] * H[(i + 1) + j * lead]; + H[(i + 1) + j * lead] = -sn[i] * H[i + j * lead] + cs[i] * H[(i + 1) + j * lead]; + H[i + j * lead] = temp; + } + + /* Compute new rotation to zero H[j+1,j] */ + double h_jj = H[j + j * lead]; + double h_j1j = H[(j + 1) + j * lead]; + double denom = sqrt(h_jj * h_jj + h_j1j * h_j1j); + if (denom < GMRES_BREAKDOWN_THRESHOLD) { + /* Degenerate column (both entries ~0): identity rotation. */ + cs[j] = 1.0; + sn[j] = 0.0; + } else { + cs[j] = h_jj / denom; + sn[j] = h_j1j / denom; + } + + /* Apply the new rotation to H (H[j,j] becomes denom, H[j+1,j] becomes 0) */ + H[j + j * lead] = cs[j] * h_jj + sn[j] * h_j1j; + H[(j + 1) + j * lead] = 0.0; + + /* Rotate g */ + double g_j = g[j]; + g[j] = cs[j] * g_j; + g[j + 1] = -sn[j] * g_j; +} + +/** + * Solve the k x k upper-triangular top block of H for y: H(0:k,0:k) y = g(0:k). + * H column-major, (m+1) x m. + */ +static void back_substitution(const double* H, const double* g, double* y, + int k, int m) { + const int lead = m + 1; + for (int i = k - 1; i >= 0; i--) { + double sum = g[i]; + for (int l = i + 1; l < k; l++) { + sum -= H[i + l * lead] * y[l]; + } + double diag = H[i + i * lead]; + y[i] = (fabs(diag) > GMRES_BREAKDOWN_THRESHOLD) ? (sum / diag) : 0.0; + } +} + +/* ============================================================================ + * GMRES SCALAR IMPLEMENTATION + * ============================================================================ */ + +static void gmres_scalar_destroy(poisson_solver_t* solver) { + if (solver && solver->context) { + gmres_context_t* ctx = (gmres_context_t*)solver->context; + cfd_free(ctx->Vblock); + cfd_free(ctx->w); + cfd_free(ctx->Mv); + cfd_free(ctx->H); + cfd_free(ctx->cs); + cfd_free(ctx->sn); + cfd_free(ctx->g); + cfd_free(ctx->y); + cfd_free(ctx); + solver->context = NULL; + } +} + +static cfd_status_t gmres_scalar_init( + poisson_solver_t* solver, + size_t nx, size_t ny, size_t nz, + double dx, double dy, double dz, + const poisson_solver_params_t* params) +{ + gmres_context_t* ctx = (gmres_context_t*)cfd_calloc(1, sizeof(gmres_context_t)); + if (!ctx) { + return CFD_ERROR_NOMEM; + } + + ctx->dx2 = dx * dx; + ctx->dy2 = dy * dy; + ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz); + poisson_solver_compute_3d_bounds(nz, nx, ny, &ctx->stride_z, &ctx->k_start, &ctx->k_end); + + ctx->diag_inv = 1.0 / (2.0 / ctx->dx2 + 2.0 / ctx->dy2 + 2.0 * ctx->inv_dz2); + ctx->use_precond = (params && params->preconditioner == POISSON_PRECOND_JACOBI); + + /* Resolve restart length m (0 = auto), clamp to >= 1 */ + int m = (params && params->restart > 0) ? params->restart : GMRES_DEFAULT_RESTART; + if (m < 1) { + m = 1; + } + ctx->m = m; + + size_t n = nx * ny * nz; + ctx->n = n; + + /* Large O(n) arrays */ + ctx->Vblock = (double*)cfd_calloc((size_t)(m + 1) * n, sizeof(double)); + ctx->w = (double*)cfd_calloc(n, sizeof(double)); + ctx->Mv = ctx->use_precond ? (double*)cfd_calloc(n, sizeof(double)) : NULL; + + /* Small dense arrays */ + ctx->H = (double*)cfd_calloc((size_t)(m + 1) * m, sizeof(double)); + ctx->cs = (double*)cfd_calloc(m, sizeof(double)); + ctx->sn = (double*)cfd_calloc(m, sizeof(double)); + ctx->g = (double*)cfd_calloc((size_t)m + 1, sizeof(double)); + ctx->y = (double*)cfd_calloc(m, sizeof(double)); + + int precond_ok = (!ctx->use_precond) || (ctx->Mv != NULL); + if (!ctx->Vblock || !ctx->w || !ctx->H || !ctx->cs || !ctx->sn || + !ctx->g || !ctx->y || !precond_ok) { + solver->context = ctx; + gmres_scalar_destroy(solver); + return CFD_ERROR_NOMEM; + } + + ctx->initialized = 1; + solver->context = ctx; + return CFD_SUCCESS; +} + +static cfd_status_t gmres_scalar_solve( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + poisson_solver_stats_t* stats) +{ + (void)x_temp; /* GMRES uses its own basis storage */ + + gmres_context_t* ctx = (gmres_context_t*)solver->context; + size_t nx = solver->nx; + size_t ny = solver->ny; + double dx2 = ctx->dx2; + double dy2 = ctx->dy2; + double inv_dz2 = ctx->inv_dz2; + size_t stride_z = ctx->stride_z; + size_t k_start = ctx->k_start; + size_t k_end = ctx->k_end; + + int m = ctx->m; + size_t n = ctx->n; + const int lead = m + 1; + + double* Vblock = ctx->Vblock; + double* w = ctx->w; + double* Mv = ctx->Mv; + double* H = ctx->H; + double* cs = ctx->cs; + double* sn = ctx->sn; + double* g = ctx->g; + double* y = ctx->y; + int use_precond = ctx->use_precond; + double diag_inv = ctx->diag_inv; + + poisson_solver_params_t* params = &solver->params; + double start_time = poisson_solver_get_time_ms(); + + /* Apply initial boundary conditions */ + poisson_solver_apply_bc(solver, x); + + /* Initial residual r_0 = b - A x_0 into V[0] */ + compute_residual(x, rhs, Vblock, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + double beta = vector_norm(Vblock, nx, ny, k_start, k_end, stride_z); + double initial_res = beta; + + if (stats) { + stats->initial_residual = initial_res; + } + + double tolerance = params->tolerance * initial_res; + if (tolerance < params->absolute_tolerance) { + tolerance = params->absolute_tolerance; + } + + if (initial_res < params->absolute_tolerance) { + poisson_solver_apply_bc(solver, x); + if (stats) { + stats->status = POISSON_CONVERGED; + stats->iterations = 0; + stats->final_residual = initial_res; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return CFD_SUCCESS; + } + + int total_inner = 0; + int converged = 0; + double final_res = initial_res; + + while (total_inner < params->max_iterations && !converged) { + /* V[0] and beta hold the current TRUE residual (fresh on first pass, + * recomputed at the end of each restart below). */ + if (beta < tolerance || beta < params->absolute_tolerance) { + converged = 1; + final_res = beta; + break; + } + + /* v_0 = r / beta */ + scale_vector(1.0 / beta, Vblock, nx, ny, k_start, k_end, stride_z); + + /* Reset least-squares RHS: g = (beta, 0, ..., 0) */ + memset(g, 0, (size_t)(m + 1) * sizeof(double)); + g[0] = beta; + + int k = 0; + for (int j = 0; j < m; j++) { + double* v_j = Vblock + (size_t)j * n; + + /* Arnoldi: w = A M^{-1} v_j (w = A v_j unpreconditioned) */ + if (use_precond) { + apply_jacobi_precond(v_j, Mv, nx, ny, diag_inv, k_start, k_end, stride_z); + apply_laplacian(Mv, w, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + } else { + apply_laplacian(v_j, w, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + } + + /* Modified Gram-Schmidt against v_0..v_j */ + for (int i = 0; i <= j; i++) { + double* v_i = Vblock + (size_t)i * n; + double hij = dot_product(w, v_i, nx, ny, k_start, k_end, stride_z); + H[i + j * lead] = hij; + axpy(-hij, v_i, w, nx, ny, k_start, k_end, stride_z); + } + + double hnorm = vector_norm(w, nx, ny, k_start, k_end, stride_z); + H[(j + 1) + j * lead] = hnorm; + + /* Happy breakdown: solution lies in the current Krylov subspace. */ + int happy = (hnorm < GMRES_BREAKDOWN_THRESHOLD); + if (!happy) { + scale_vector(1.0 / hnorm, w, nx, ny, k_start, k_end, stride_z); + copy_vector(w, Vblock + (size_t)(j + 1) * n, nx, ny, k_start, k_end, stride_z); + } + + apply_givens_and_update_hessenberg(H, cs, sn, g, j, m); + total_inner++; + k = j + 1; + + /* Cheap Givens residual estimate: only decides when to END the inner + * cycle early (to update x and recheck the TRUE residual). It does NOT + * decide convergence, so drift in the estimate cannot cause a false + * positive. */ + double resid_est = fabs(g[j + 1]); + if (params->verbose) { + CFD_LOG_DEBUG("poisson", "GMRES inner %d (total %d): residual est = %.6e", + j, total_inner, resid_est); + } + if (resid_est < tolerance || happy || total_inner >= params->max_iterations) { + break; + } + } + + /* Solve the least-squares problem and update x */ + back_substitution(H, g, y, k, m); + if (use_precond) { + memset(w, 0, n * sizeof(double)); + linear_combination(w, Vblock, n, y, k, nx, ny, k_start, k_end, stride_z); + apply_jacobi_precond(w, Mv, nx, ny, diag_inv, k_start, k_end, stride_z); + axpy(1.0, Mv, x, nx, ny, k_start, k_end, stride_z); + } else { + linear_combination(x, Vblock, n, y, k, nx, ny, k_start, k_end, stride_z); + } + + /* Recompute the TRUE residual and decide convergence on it (prevents the + * cheap Givens estimate from drifting below the real residual). This is + * the residual reported to the caller, measured with the same fixed + * boundary values used throughout the solve. */ + compute_residual(x, rhs, Vblock, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + beta = vector_norm(Vblock, nx, ny, k_start, k_end, stride_z); + final_res = beta; + if (beta < tolerance || beta < params->absolute_tolerance) { + converged = 1; + } + } + + /* Apply final boundary conditions for output (mirrors CG). The reported + * residual is the solve-consistent one measured above, NOT recomputed here. */ + poisson_solver_apply_bc(solver, x); + + if (stats) { + stats->iterations = total_inner; + stats->final_residual = final_res; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + stats->status = converged ? POISSON_CONVERGED : POISSON_MAX_ITER; + } + + return converged ? CFD_SUCCESS : CFD_ERROR_MAX_ITER; +} + +/** + * Single iteration is not well-defined for GMRES (it maintains an internal + * Krylov basis). Mirror CG: return the current residual only. + */ +static cfd_status_t gmres_scalar_iterate( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + double* residual) +{ + (void)x_temp; + if (residual) { + *residual = poisson_solver_compute_residual(solver, x, rhs); + } + return CFD_SUCCESS; +} + +/* ============================================================================ + * FACTORY FUNCTION + * ============================================================================ */ + +poisson_solver_t* create_gmres_scalar_solver(void) { + poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t)); + if (!solver) { + return NULL; + } + + solver->name = POISSON_SOLVER_TYPE_GMRES_SCALAR; + solver->description = "Restarted GMRES(m) (scalar CPU)"; + solver->method = POISSON_METHOD_GMRES; + solver->backend = POISSON_BACKEND_SCALAR; + solver->params = poisson_solver_params_default(); + + solver->init = gmres_scalar_init; + solver->destroy = gmres_scalar_destroy; + solver->solve = gmres_scalar_solve; /* GMRES uses custom solve loop */ + solver->iterate = gmres_scalar_iterate; + solver->apply_bc = NULL; /* Use default Neumann */ + + return solver; +} diff --git a/lib/src/solvers/linear/linear_solver.c b/lib/src/solvers/linear/linear_solver.c index d5744102..2f30e52b 100644 --- a/lib/src/solvers/linear/linear_solver.c +++ b/lib/src/solvers/linear/linear_solver.c @@ -43,6 +43,7 @@ poisson_solver_params_t poisson_solver_params_default(void) { params.check_interval = 1; params.verbose = false; params.preconditioner = POISSON_PRECOND_NONE; + params.restart = 0; /* 0 = auto (GMRES_DEFAULT_RESTART); ignored by non-GMRES methods */ return params; } @@ -233,6 +234,20 @@ poisson_solver_t* poisson_solver_create( return NULL; /* Requested backend not available for BiCGSTAB */ } + case POISSON_METHOD_GMRES: + switch (backend) { + case POISSON_BACKEND_SIMD: + return create_gmres_simd_solver(); +#ifdef CFD_ENABLE_OPENMP + case POISSON_BACKEND_OMP: + return create_gmres_omp_solver(); +#endif + case POISSON_BACKEND_SCALAR: + return create_gmres_scalar_solver(); + default: + return NULL; /* Requested backend not available for GMRES */ + } + case POISSON_METHOD_MULTIGRID: /* Not yet implemented */ return NULL; diff --git a/lib/src/solvers/linear/linear_solver_internal.h b/lib/src/solvers/linear/linear_solver_internal.h index eb733668..7664405f 100644 --- a/lib/src/solvers/linear/linear_solver_internal.h +++ b/lib/src/solvers/linear/linear_solver_internal.h @@ -61,6 +61,14 @@ poisson_solver_t* create_redblack_gpu_solver(void); poisson_solver_t* create_bicgstab_scalar_solver(void); poisson_solver_t* create_bicgstab_simd_solver(void); +/* GMRES solvers (restarted GMRES(m) for non-symmetric systems) */ +poisson_solver_t* create_gmres_scalar_solver(void); +poisson_solver_t* create_gmres_simd_solver(void); + +#ifdef CFD_ENABLE_OPENMP +poisson_solver_t* create_gmres_omp_solver(void); +#endif + /* ============================================================================ * CG ALGORITHM CONSTANTS * ============================================================================ */ @@ -105,6 +113,25 @@ poisson_solver_t* create_bicgstab_simd_solver(void); */ #define BICGSTAB_BREAKDOWN_THRESHOLD 1e-30 +/* ============================================================================ + * GMRES ALGORITHM CONSTANTS + * ============================================================================ */ + +/** + * Default restart length m for GMRES(m) when params.restart <= 0. + * The Krylov basis holds m+1 grid-sized vectors, so this bounds memory. + */ +#define GMRES_DEFAULT_RESTART 30 + +/** + * Threshold for detecting GMRES happy (lucky) breakdown. + * When the Arnoldi subdiagonal H[j+1,j] (a vector norm) falls below this, the + * exact solution already lies in the current Krylov subspace: this signals + * CONVERGENCE, not failure. Larger than CG_BREAKDOWN_THRESHOLD because + * H[j+1,j] is an sqrt-scaled norm rather than a squared dot product. + */ +#define GMRES_BREAKDOWN_THRESHOLD 1e-12 + /** * Convert size_t to int for OpenMP loop bounds. * OpenMP requires int loop variables, but grid dimensions are size_t. diff --git a/lib/src/solvers/linear/neon/linear_solver_gmres_neon.c b/lib/src/solvers/linear/neon/linear_solver_gmres_neon.c new file mode 100644 index 00000000..f907e005 --- /dev/null +++ b/lib/src/solvers/linear/neon/linear_solver_gmres_neon.c @@ -0,0 +1,68 @@ +/** + * @file linear_solver_gmres_neon.c + * @brief Restarted GMRES(m) linear solver - ARM NEON + OpenMP backend + * + * Provides NEON-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 +#include +#include + +#if (defined(__aarch64__) || defined(_M_ARM64) || defined(__ARM_NEON) || defined(__ARM_NEON__)) && defined(CFD_ENABLE_OPENMP) +#define GMRES_HAS_NEON 1 +#include +#include +#endif + +#if defined(GMRES_HAS_NEON) + +//============================================================================= +// NEON-SPECIFIC MACRO DEFINITIONS +//============================================================================= + +#define SIMD_SUFFIX neon +#define SIMD_VEC float64x2_t +#define SIMD_WIDTH 2 + +#define SIMD_LOAD(ptr) vld1q_f64(ptr) +#define SIMD_STORE(ptr, vec) vst1q_f64(ptr, vec) +#define SIMD_SET1(val) vdupq_n_f64(val) +#define SIMD_SETZERO() vdupq_n_f64(0.0) + +#define SIMD_ADD(a, b) vaddq_f64(a, b) +#define SIMD_SUB(a, b) vsubq_f64(a, b) +#define SIMD_MUL(a, b) vmulq_f64(a, b) +#define SIMD_FMA(a, b, c) vfmaq_f64(c, a, b) + +/* NEON horizontal sum (2 doubles -> 1 double) */ +#define SIMD_HSUM(vec) (vgetq_lane_f64(vec, 0) + vgetq_lane_f64(vec, 1)) + +#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_NEON */ + +/* Stub for platforms without NEON */ +poisson_solver_t* create_gmres_neon_solver(void) { + return NULL; +} + +#endif /* GMRES_HAS_NEON */ diff --git a/lib/src/solvers/linear/omp/linear_solver_gmres_omp.c b/lib/src/solvers/linear/omp/linear_solver_gmres_omp.c new file mode 100644 index 00000000..3937c7e0 --- /dev/null +++ b/lib/src/solvers/linear/omp/linear_solver_gmres_omp.c @@ -0,0 +1,519 @@ +/** + * @file linear_solver_gmres_omp.c + * @brief Restarted GMRES(m) solver - OpenMP parallelized implementation + * + * Same algorithm as scalar GMRES (linear_solver_gmres.c) but with the O(n) + * vector primitives OpenMP-parallelized. The dense Givens/Hessenberg and + * back-substitution work is O(m^2) and stays scalar, identical to the scalar + * and SIMD backends, so results are consistent across backends. + * + * Operator/sign convention matches CG and scalar GMRES: A = -nabla^2, b = -rhs. + */ + +#include "../linear_solver_internal.h" + +#include "cfd/core/indexing.h" +#include "cfd/core/logging.h" +#include "cfd/core/memory.h" + +#include +#include + +#ifdef CFD_ENABLE_OPENMP + +#include + +/* ============================================================================ + * GMRES OMP CONTEXT + * ============================================================================ */ + +typedef struct { + double dx2; + double dy2; + double inv_dz2; + double diag_inv; + + size_t stride_z; + size_t k_start; + size_t k_end; + + int m; + size_t n; + + double* Vblock; + double* w; + double* Mv; + + double* H; + double* cs; + double* sn; + double* g; + double* y; + + int use_precond; + int initialized; +} gmres_omp_context_t; + +/* ============================================================================ + * OMP-PARALLELIZED O(n) PRIMITIVES (interior points only) + * ============================================================================ */ + +static inline int size_to_int(size_t val) { + if (val > (size_t)INT_MAX) { + return INT_MAX; + } + return (int)val; +} + +static double dot_product_omp(const double* a, const double* b, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + double sum = 0.0; + int ny_int = size_to_int(ny); + int nx_int = size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) reduction(+:sum) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + sum += a[idx] * b[idx]; + } + } + } + return sum; +} + +static double vector_norm_omp(const double* x, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + return sqrt(dot_product_omp(x, x, nx, ny, k_start, k_end, stride_z)); +} + +static void axpy_omp(double alpha, const double* x, double* y, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = size_to_int(ny); + int nx_int = size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + y[idx] += alpha * x[idx]; + } + } + } +} + +static void scale_vector_omp(double alpha, double* x, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = size_to_int(ny); + int nx_int = size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + x[idx] *= alpha; + } + } + } +} + +static void apply_laplacian_omp(const double* p, double* Ap, + size_t nx, size_t ny, + double dx2, double dy2, double inv_dz2, + size_t k_start, size_t k_end, size_t stride_z) { + double dx2_inv = 1.0 / dx2; + double dy2_inv = 1.0 / dy2; + int ny_int = size_to_int(ny); + int nx_int = size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + double laplacian = + (p[idx + 1] - 2.0 * p[idx] + p[idx - 1]) * dx2_inv + + (p[idx + nx] - 2.0 * p[idx] + p[idx - nx]) * dy2_inv + + (p[idx + stride_z] + p[idx - stride_z] - 2.0 * p[idx]) * inv_dz2; + Ap[idx] = -laplacian; + } + } + } +} + +static void compute_residual_omp(const double* x, const double* rhs, double* r, + size_t nx, size_t ny, + double dx2, double dy2, double inv_dz2, + size_t k_start, size_t k_end, size_t stride_z) { + double dx2_inv = 1.0 / dx2; + double dy2_inv = 1.0 / dy2; + int ny_int = size_to_int(ny); + int nx_int = size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + double laplacian = + (x[idx + 1] - 2.0 * x[idx] + x[idx - 1]) * dx2_inv + + (x[idx + nx] - 2.0 * x[idx] + x[idx - nx]) * dy2_inv + + (x[idx + stride_z] + x[idx - stride_z] - 2.0 * x[idx]) * inv_dz2; + r[idx] = -rhs[idx] + laplacian; + } + } + } +} + +static void copy_vector_omp(const double* src, double* dst, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = size_to_int(ny); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + size_t row_start = k * stride_z + (size_t)j * nx; + memcpy(&dst[row_start + 1], &src[row_start + 1], (nx - 2) * sizeof(double)); + } + } +} + +static void apply_jacobi_precond_omp(const double* r, double* z, + size_t nx, size_t ny, + double diag_inv, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = size_to_int(ny); + int nx_int = size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + z[idx] = diag_inv * r[idx]; + } + } + } +} + +/* ============================================================================ + * DENSE SCALAR HELPERS (identical to scalar backend; not parallelized) + * ============================================================================ */ + +static void gmres_givens_omp(double* H, double* cs, double* sn, + double* g, int j, int m) { + const int lead = m + 1; + for (int i = 0; i < j; i++) { + double temp = cs[i] * H[i + j * lead] + sn[i] * H[(i + 1) + j * lead]; + H[(i + 1) + j * lead] = -sn[i] * H[i + j * lead] + cs[i] * H[(i + 1) + j * lead]; + H[i + j * lead] = temp; + } + + double h_jj = H[j + j * lead]; + double h_j1j = H[(j + 1) + j * lead]; + double denom = sqrt(h_jj * h_jj + h_j1j * h_j1j); + if (denom < GMRES_BREAKDOWN_THRESHOLD) { + cs[j] = 1.0; + sn[j] = 0.0; + } else { + cs[j] = h_jj / denom; + sn[j] = h_j1j / denom; + } + + H[j + j * lead] = cs[j] * h_jj + sn[j] * h_j1j; + H[(j + 1) + j * lead] = 0.0; + + double g_j = g[j]; + g[j] = cs[j] * g_j; + g[j + 1] = -sn[j] * g_j; +} + +static void gmres_backsub_omp(const double* H, const double* g, double* y, + int k, int m) { + const int lead = m + 1; + for (int i = k - 1; i >= 0; i--) { + double sum = g[i]; + for (int l = i + 1; l < k; l++) { + sum -= H[i + l * lead] * y[l]; + } + double diag = H[i + i * lead]; + y[i] = (fabs(diag) > GMRES_BREAKDOWN_THRESHOLD) ? (sum / diag) : 0.0; + } +} + +/* ============================================================================ + * GMRES OMP IMPLEMENTATION + * ============================================================================ */ + +static void gmres_omp_destroy(poisson_solver_t* solver) { + if (solver && solver->context) { + gmres_omp_context_t* ctx = (gmres_omp_context_t*)solver->context; + cfd_free(ctx->Vblock); + cfd_free(ctx->w); + cfd_free(ctx->Mv); + cfd_free(ctx->H); + cfd_free(ctx->cs); + cfd_free(ctx->sn); + cfd_free(ctx->g); + cfd_free(ctx->y); + cfd_free(ctx); + solver->context = NULL; + } +} + +static cfd_status_t gmres_omp_init( + poisson_solver_t* solver, + size_t nx, size_t ny, size_t nz, + double dx, double dy, double dz, + const poisson_solver_params_t* params) +{ + gmres_omp_context_t* ctx = (gmres_omp_context_t*)cfd_calloc(1, sizeof(gmres_omp_context_t)); + if (!ctx) { + return CFD_ERROR_NOMEM; + } + + ctx->dx2 = dx * dx; + ctx->dy2 = dy * dy; + ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz); + poisson_solver_compute_3d_bounds(nz, nx, ny, &ctx->stride_z, &ctx->k_start, &ctx->k_end); + ctx->diag_inv = 1.0 / (2.0 / ctx->dx2 + 2.0 / ctx->dy2 + 2.0 * ctx->inv_dz2); + ctx->use_precond = (params && params->preconditioner == POISSON_PRECOND_JACOBI); + + int m = (params && params->restart > 0) ? params->restart : GMRES_DEFAULT_RESTART; + if (m < 1) m = 1; + ctx->m = m; + + size_t n = nx * ny * nz; + ctx->n = n; + + ctx->Vblock = (double*)cfd_calloc((size_t)(m + 1) * n, sizeof(double)); + ctx->w = (double*)cfd_calloc(n, sizeof(double)); + ctx->Mv = ctx->use_precond ? (double*)cfd_calloc(n, sizeof(double)) : NULL; + + ctx->H = (double*)cfd_calloc((size_t)(m + 1) * m, sizeof(double)); + ctx->cs = (double*)cfd_calloc(m, sizeof(double)); + ctx->sn = (double*)cfd_calloc(m, sizeof(double)); + ctx->g = (double*)cfd_calloc((size_t)m + 1, sizeof(double)); + ctx->y = (double*)cfd_calloc(m, sizeof(double)); + + int precond_ok = (!ctx->use_precond) || (ctx->Mv != NULL); + if (!ctx->Vblock || !ctx->w || !ctx->H || !ctx->cs || !ctx->sn || + !ctx->g || !ctx->y || !precond_ok) { + solver->context = ctx; + gmres_omp_destroy(solver); + return CFD_ERROR_NOMEM; + } + + ctx->initialized = 1; + solver->context = ctx; + return CFD_SUCCESS; +} + +static cfd_status_t gmres_omp_solve( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + poisson_solver_stats_t* stats) +{ + (void)x_temp; + + gmres_omp_context_t* ctx = (gmres_omp_context_t*)solver->context; + size_t nx = solver->nx; + size_t ny = solver->ny; + double dx2 = ctx->dx2; + double dy2 = ctx->dy2; + double inv_dz2 = ctx->inv_dz2; + size_t stride_z = ctx->stride_z; + size_t k_start = ctx->k_start; + size_t k_end = ctx->k_end; + + int m = ctx->m; + size_t n = ctx->n; + const int lead = m + 1; + + double* Vblock = ctx->Vblock; + double* w = ctx->w; + double* Mv = ctx->Mv; + double* H = ctx->H; + double* cs = ctx->cs; + double* sn = ctx->sn; + double* g = ctx->g; + double* y = ctx->y; + int use_precond = ctx->use_precond; + double diag_inv = ctx->diag_inv; + + poisson_solver_params_t* params = &solver->params; + double start_time = poisson_solver_get_time_ms(); + + poisson_solver_apply_bc(solver, x); + + compute_residual_omp(x, rhs, Vblock, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + double beta = vector_norm_omp(Vblock, nx, ny, k_start, k_end, stride_z); + double initial_res = beta; + + if (stats) { + stats->initial_residual = initial_res; + } + + double tolerance = params->tolerance * initial_res; + if (tolerance < params->absolute_tolerance) { + tolerance = params->absolute_tolerance; + } + + if (initial_res < params->absolute_tolerance) { + poisson_solver_apply_bc(solver, x); + if (stats) { + stats->status = POISSON_CONVERGED; + stats->iterations = 0; + stats->final_residual = initial_res; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return CFD_SUCCESS; + } + + int total_inner = 0; + int converged = 0; + double final_res = initial_res; + + while (total_inner < params->max_iterations && !converged) { + if (beta < tolerance || beta < params->absolute_tolerance) { + converged = 1; + final_res = beta; + break; + } + + scale_vector_omp(1.0 / beta, Vblock, nx, ny, k_start, k_end, stride_z); + memset(g, 0, (size_t)(m + 1) * sizeof(double)); + g[0] = beta; + + int k = 0; + for (int j = 0; j < m; j++) { + double* v_j = Vblock + (size_t)j * n; + + if (use_precond) { + apply_jacobi_precond_omp(v_j, Mv, nx, ny, diag_inv, k_start, k_end, stride_z); + apply_laplacian_omp(Mv, w, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + } else { + apply_laplacian_omp(v_j, w, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + } + + for (int i = 0; i <= j; i++) { + double* v_i = Vblock + (size_t)i * n; + double hij = dot_product_omp(w, v_i, nx, ny, k_start, k_end, stride_z); + H[i + j * lead] = hij; + axpy_omp(-hij, v_i, w, nx, ny, k_start, k_end, stride_z); + } + + double hnorm = vector_norm_omp(w, nx, ny, k_start, k_end, stride_z); + H[(j + 1) + j * lead] = hnorm; + + int happy = (hnorm < GMRES_BREAKDOWN_THRESHOLD); + if (!happy) { + scale_vector_omp(1.0 / hnorm, w, nx, ny, k_start, k_end, stride_z); + copy_vector_omp(w, Vblock + (size_t)(j + 1) * n, nx, ny, k_start, k_end, stride_z); + } + + gmres_givens_omp(H, cs, sn, g, j, m); + total_inner++; + k = j + 1; + + double resid_est = fabs(g[j + 1]); + if (params->verbose) { + CFD_LOG_DEBUG("poisson", "GMRES-OMP inner %d (total %d): residual est = %.6e", + j, total_inner, resid_est); + } + if (resid_est < tolerance || happy || total_inner >= params->max_iterations) { + break; + } + } + + gmres_backsub_omp(H, g, y, k, m); + if (use_precond) { + memset(w, 0, n * sizeof(double)); + for (int jj = 0; jj < k; jj++) { + axpy_omp(y[jj], Vblock + (size_t)jj * n, w, nx, ny, k_start, k_end, stride_z); + } + apply_jacobi_precond_omp(w, Mv, nx, ny, diag_inv, k_start, k_end, stride_z); + axpy_omp(1.0, Mv, x, nx, ny, k_start, k_end, stride_z); + } else { + for (int jj = 0; jj < k; jj++) { + axpy_omp(y[jj], Vblock + (size_t)jj * n, x, nx, ny, k_start, k_end, stride_z); + } + } + + compute_residual_omp(x, rhs, Vblock, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + beta = vector_norm_omp(Vblock, nx, ny, k_start, k_end, stride_z); + final_res = beta; + if (beta < tolerance || beta < params->absolute_tolerance) { + converged = 1; + } + } + + poisson_solver_apply_bc(solver, x); + + if (stats) { + stats->iterations = total_inner; + stats->final_residual = final_res; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + stats->status = converged ? POISSON_CONVERGED : POISSON_MAX_ITER; + } + + return converged ? CFD_SUCCESS : CFD_ERROR_MAX_ITER; +} + +static cfd_status_t gmres_omp_iterate( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + double* residual) +{ + (void)x_temp; + if (residual) { + *residual = poisson_solver_compute_residual(solver, x, rhs); + } + return CFD_SUCCESS; +} + +/* ============================================================================ + * FACTORY FUNCTION + * ============================================================================ */ + +poisson_solver_t* create_gmres_omp_solver(void) { + poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t)); + if (!solver) { + return NULL; + } + + solver->name = POISSON_SOLVER_TYPE_GMRES_OMP; + solver->description = "Restarted GMRES(m) (OpenMP)"; + solver->method = POISSON_METHOD_GMRES; + solver->backend = POISSON_BACKEND_OMP; + solver->params = poisson_solver_params_default(); + + solver->init = gmres_omp_init; + solver->destroy = gmres_omp_destroy; + solver->solve = gmres_omp_solve; + solver->iterate = gmres_omp_iterate; + solver->apply_bc = NULL; + + return solver; +} + +#endif /* CFD_ENABLE_OPENMP */ diff --git a/lib/src/solvers/linear/simd/linear_solver_simd_dispatch.c b/lib/src/solvers/linear/simd/linear_solver_simd_dispatch.c index 25668320..782c9671 100644 --- a/lib/src/solvers/linear/simd/linear_solver_simd_dispatch.c +++ b/lib/src/solvers/linear/simd/linear_solver_simd_dispatch.c @@ -48,6 +48,7 @@ extern poisson_solver_t* create_redblack_avx2_solver(void); extern poisson_solver_t* create_cg_avx2_solver(void); extern poisson_solver_t* create_bicgstab_avx2_solver(void); extern poisson_solver_t* create_sor_avx2_solver(void); +extern poisson_solver_t* create_gmres_avx2_solver(void); /* NEON implementations (ARM64) */ extern poisson_solver_t* create_jacobi_neon_solver(void); @@ -55,6 +56,7 @@ extern poisson_solver_t* create_redblack_neon_solver(void); extern poisson_solver_t* create_cg_neon_solver(void); extern poisson_solver_t* create_bicgstab_neon_solver(void); extern poisson_solver_t* create_sor_neon_solver(void); +extern poisson_solver_t* create_gmres_neon_solver(void); /* ============================================================================ * SIMD BACKEND AVAILABILITY (Runtime + Compile-time detection) @@ -217,3 +219,24 @@ poisson_solver_t* create_sor_simd_solver(void) { log_no_simd_available("SOR"); return NULL; } + +/* ============================================================================ + * GMRES SIMD DISPATCHER + * ============================================================================ */ + +poisson_solver_t* create_gmres_simd_solver(void) { + cfd_simd_arch_t arch = cfd_detect_simd_arch(); + + /* Check compile-time AND runtime availability before calling factory */ + if (HAS_AVX2_IMPL && arch == CFD_SIMD_AVX2) { + return create_gmres_avx2_solver(); + } + + if (HAS_NEON_IMPL && arch == CFD_SIMD_NEON) { + return create_gmres_neon_solver(); + } + + /* No SIMD backend available - report error and return NULL (no fallback) */ + log_no_simd_available("GMRES"); + return NULL; +} diff --git a/lib/src/solvers/linear/simd_template/linear_solver_gmres_simd_template.h b/lib/src/solvers/linear/simd_template/linear_solver_gmres_simd_template.h new file mode 100644 index 00000000..4266119d --- /dev/null +++ b/lib/src/solvers/linear/simd_template/linear_solver_gmres_simd_template.h @@ -0,0 +1,662 @@ +/** + * @file linear_solver_gmres_simd_template.h + * @brief Restarted GMRES(m) SIMD implementation template (AVX2/NEON parameterized) + * + * Parameterized header implementing GMRES(m) with architecture-agnostic SIMD + * macros. Included twice (AVX2, NEON) with different macro definitions. Only the + * O(n) vector primitives are vectorized; the small dense Givens/Hessenberg and + * back-substitution work is plain scalar (identical to the scalar backend), so + * SIMD and scalar results stay consistent to rounding. + * + * The operator/sign convention matches linear_solver_gmres.c and CG exactly: + * A = -nabla^2, b = -rhs. Supports 2D (nz==1) and 3D via branch-free stencil + * (stride_z=0, inv_dz2=0.0 for 2D). + * + * REQUIRED MACROS (see linear_solver_bicgstab_simd_template.h for the full list): + * SIMD_SUFFIX, SIMD_VEC, SIMD_WIDTH, SIMD_LOAD, SIMD_STORE, SIMD_SET1, + * SIMD_SETZERO, SIMD_ADD, SIMD_SUB, SIMD_MUL, SIMD_FMA, SIMD_HSUM. + */ + +#include "cfd/core/indexing.h" + +#include +#include + +#ifndef SIMD_SUFFIX +#error "SIMD_SUFFIX must be defined before including linear_solver_gmres_simd_template.h" +#endif +#ifndef SIMD_VEC +#error "SIMD_VEC must be defined before including linear_solver_gmres_simd_template.h" +#endif +#ifndef SIMD_WIDTH +#error "SIMD_WIDTH must be defined before including linear_solver_gmres_simd_template.h" +#endif + +//============================================================================= +// TOKEN PASTING MACROS +//============================================================================= + +#define CONCAT_IMPL(a, b) a##_##b +#define CONCAT(a, b) CONCAT_IMPL(a, b) +#define SIMD_FUNC(name) CONCAT(name, SIMD_SUFFIX) + +/* Factory function name: create_gmres__solver */ +#define FACTORY_NAME_IMPL(suffix) create_gmres_##suffix##_solver +#define FACTORY_NAME(suffix) FACTORY_NAME_IMPL(suffix) + +//============================================================================= +// CONTEXT STRUCTURE +//============================================================================= + +#define gmres_simd_context_t SIMD_FUNC(gmres_context_t) + +typedef struct { + double dx2; + double dy2; + double inv_dz2; + double diag_inv; + + size_t stride_z; + size_t k_start; + size_t k_end; + + int m; /* restart length */ + size_t n; /* total field size */ + + /* Precomputed SIMD constant vectors */ + SIMD_VEC dx2_inv_vec; + SIMD_VEC dy2_inv_vec; + SIMD_VEC dz2_inv_vec; + SIMD_VEC two_vec; + + /* Large O(n) working arrays */ + double* Vblock; /* (m+1) Krylov basis vectors */ + double* w; + double* Mv; /* NULL if no precond */ + + /* Small dense O(m^2)/O(m) arrays (scalar) */ + double* H; + double* cs; + double* sn; + double* g; + double* y; + + int use_precond; + int initialized; +} gmres_simd_context_t; + +//============================================================================= +// SIMD O(n) PRIMITIVES (interior points only) +//============================================================================= + +/** Dot product over interior points (SIMD + OpenMP reduction) */ +static inline double SIMD_FUNC(gmres_dot)(const double* a, const double* b, + size_t nx, size_t ny, + size_t k_start, size_t k_end, + size_t stride_z) { + double sum = 0.0; + int ny_int = bicgstab_size_to_int(ny); + if (ny_int == 0) return 0.0; + + for (size_t k = k_start; k < k_end; k++) { + int jj; + #pragma omp parallel for reduction(+:sum) schedule(static) + for (jj = 1; jj < ny_int - 1; jj++) { + size_t j = (size_t)jj; + double row_sum = 0.0; + size_t i = 1; + SIMD_VEC acc = SIMD_SETZERO(); + + for (; i + SIMD_WIDTH - 1 < nx - 1; i += SIMD_WIDTH) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + SIMD_VEC va = SIMD_LOAD(&a[idx]); + SIMD_VEC vb = SIMD_LOAD(&b[idx]); + acc = SIMD_FMA(va, vb, acc); + } + row_sum += SIMD_HSUM(acc); + + for (; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + row_sum += a[idx] * b[idx]; + } + sum += row_sum; + } + } + return sum; +} + +/** Euclidean norm over interior points */ +static inline double SIMD_FUNC(gmres_norm)(const double* x, + size_t nx, size_t ny, + size_t k_start, size_t k_end, + size_t stride_z) { + return sqrt(SIMD_FUNC(gmres_dot)(x, x, nx, ny, k_start, k_end, stride_z)); +} + +/** y = y + alpha*x (SIMD FMA) */ +static inline void SIMD_FUNC(gmres_axpy)(double alpha, const double* x, double* y, + size_t nx, size_t ny, + size_t k_start, size_t k_end, + size_t stride_z) { + SIMD_VEC alpha_vec = SIMD_SET1(alpha); + int ny_int = bicgstab_size_to_int(ny); + if (ny_int == 0) return; + + for (size_t k = k_start; k < k_end; k++) { + int jj; + #pragma omp parallel for schedule(static) + for (jj = 1; jj < ny_int - 1; jj++) { + size_t j = (size_t)jj; + size_t i = 1; + for (; i + SIMD_WIDTH - 1 < nx - 1; i += SIMD_WIDTH) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + SIMD_VEC vx = SIMD_LOAD(&x[idx]); + SIMD_VEC vy = SIMD_LOAD(&y[idx]); + vy = SIMD_FMA(alpha_vec, vx, vy); + SIMD_STORE(&y[idx], vy); + } + for (; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + y[idx] += alpha * x[idx]; + } + } + } +} + +/** x = alpha*x (SIMD) */ +static inline void SIMD_FUNC(gmres_scale)(double alpha, double* x, + size_t nx, size_t ny, + size_t k_start, size_t k_end, + size_t stride_z) { + SIMD_VEC alpha_vec = SIMD_SET1(alpha); + int ny_int = bicgstab_size_to_int(ny); + if (ny_int == 0) return; + + for (size_t k = k_start; k < k_end; k++) { + int jj; + #pragma omp parallel for schedule(static) + for (jj = 1; jj < ny_int - 1; jj++) { + size_t j = (size_t)jj; + size_t i = 1; + for (; i + SIMD_WIDTH - 1 < nx - 1; i += SIMD_WIDTH) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + SIMD_VEC vx = SIMD_LOAD(&x[idx]); + SIMD_STORE(&x[idx], SIMD_MUL(alpha_vec, vx)); + } + for (; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + x[idx] *= alpha; + } + } + } +} + +/** dst = src (interior points) */ +static inline void SIMD_FUNC(gmres_copy)(const double* src, double* dst, + size_t nx, size_t ny, + size_t k_start, size_t k_end, + size_t stride_z) { + int ny_int = bicgstab_size_to_int(ny); + if (ny_int == 0) return; + + for (size_t k = k_start; k < k_end; k++) { + int jj; + #pragma omp parallel for schedule(static) + for (jj = 1; jj < ny_int - 1; jj++) { + size_t j = (size_t)jj; + size_t i = 1; + for (; i + SIMD_WIDTH - 1 < nx - 1; i += SIMD_WIDTH) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + SIMD_STORE(&dst[idx], SIMD_LOAD(&src[idx])); + } + for (; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + dst[idx] = src[idx]; + } + } + } +} + +/** Apply negative Laplacian: Ap = -nabla^2(p) = A p (SIMD stencil) */ +static inline void SIMD_FUNC(gmres_apply_A)(const double* p, double* Ap, + size_t nx, size_t ny, + const gmres_simd_context_t* ctx) { + SIMD_VEC dx2_inv = ctx->dx2_inv_vec; + SIMD_VEC dy2_inv = ctx->dy2_inv_vec; + SIMD_VEC dz2_inv = ctx->dz2_inv_vec; + SIMD_VEC two_vec = ctx->two_vec; + SIMD_VEC zero = SIMD_SETZERO(); + size_t stride_z = ctx->stride_z; + int ny_int = bicgstab_size_to_int(ny); + if (ny_int == 0) return; + + for (size_t k = ctx->k_start; k < ctx->k_end; k++) { + int jj; + #pragma omp parallel for schedule(static) + for (jj = 1; jj < ny_int - 1; jj++) { + size_t j = (size_t)jj; + size_t i = 1; + for (; i + SIMD_WIDTH - 1 < nx - 1; i += SIMD_WIDTH) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + SIMD_VEC p_c = SIMD_LOAD(&p[idx]); + SIMD_VEC p_w = SIMD_LOAD(&p[idx - 1]); + SIMD_VEC p_e = SIMD_LOAD(&p[idx + 1]); + SIMD_VEC p_s = SIMD_LOAD(&p[idx - nx]); + SIMD_VEC p_n = SIMD_LOAD(&p[idx + nx]); + SIMD_VEC p_zm = SIMD_LOAD(&p[idx - stride_z]); + SIMD_VEC p_zp = SIMD_LOAD(&p[idx + stride_z]); + + SIMD_VEC two_center = SIMD_MUL(two_vec, p_c); + SIMD_VEC d2x = SIMD_MUL(SIMD_SUB(SIMD_ADD(p_e, p_w), two_center), dx2_inv); + SIMD_VEC d2y = SIMD_MUL(SIMD_SUB(SIMD_ADD(p_n, p_s), two_center), dy2_inv); + SIMD_VEC d2z = SIMD_MUL(SIMD_SUB(SIMD_ADD(p_zp, p_zm), two_center), dz2_inv); + SIMD_VEC laplacian = SIMD_ADD(SIMD_ADD(d2x, d2y), d2z); + /* Ap = -laplacian */ + SIMD_STORE(&Ap[idx], SIMD_SUB(zero, laplacian)); + } + for (; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + double lap = (p[idx + 1] - 2.0 * p[idx] + p[idx - 1]) / ctx->dx2 + + (p[idx + nx] - 2.0 * p[idx] + p[idx - nx]) / ctx->dy2 + + (p[idx + stride_z] + p[idx - stride_z] - 2.0 * p[idx]) * ctx->inv_dz2; + Ap[idx] = -lap; + } + } + } +} + +/** Residual r = b - A x = -rhs + nabla^2 x (SIMD, fused) */ +static inline void SIMD_FUNC(gmres_residual)(const double* x, const double* rhs, + double* r, size_t nx, size_t ny, + const gmres_simd_context_t* ctx) { + SIMD_VEC dx2_inv = ctx->dx2_inv_vec; + SIMD_VEC dy2_inv = ctx->dy2_inv_vec; + SIMD_VEC dz2_inv = ctx->dz2_inv_vec; + SIMD_VEC two_vec = ctx->two_vec; + size_t stride_z = ctx->stride_z; + int ny_int = bicgstab_size_to_int(ny); + if (ny_int == 0) return; + + for (size_t k = ctx->k_start; k < ctx->k_end; k++) { + int jj; + #pragma omp parallel for schedule(static) + for (jj = 1; jj < ny_int - 1; jj++) { + size_t j = (size_t)jj; + size_t i = 1; + for (; i + SIMD_WIDTH - 1 < nx - 1; i += SIMD_WIDTH) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + SIMD_VEC c = SIMD_LOAD(&x[idx]); + SIMD_VEC l = SIMD_LOAD(&x[idx - 1]); + SIMD_VEC rr = SIMD_LOAD(&x[idx + 1]); + SIMD_VEC d = SIMD_LOAD(&x[idx - nx]); + SIMD_VEC u = SIMD_LOAD(&x[idx + nx]); + SIMD_VEC bk = SIMD_LOAD(&x[idx - stride_z]); + SIMD_VEC fr = SIMD_LOAD(&x[idx + stride_z]); + + SIMD_VEC two_center = SIMD_MUL(two_vec, c); + SIMD_VEC d2x = SIMD_MUL(SIMD_SUB(SIMD_ADD(l, rr), two_center), dx2_inv); + SIMD_VEC d2y = SIMD_MUL(SIMD_SUB(SIMD_ADD(d, u), two_center), dy2_inv); + SIMD_VEC d2z = SIMD_MUL(SIMD_SUB(SIMD_ADD(bk, fr), two_center), dz2_inv); + SIMD_VEC lap = SIMD_ADD(SIMD_ADD(d2x, d2y), d2z); + /* r = -rhs + lap */ + SIMD_VEC vrhs = SIMD_LOAD(&rhs[idx]); + SIMD_STORE(&r[idx], SIMD_SUB(lap, vrhs)); + } + for (; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + double lap = (x[idx + 1] - 2.0 * x[idx] + x[idx - 1]) / ctx->dx2 + + (x[idx + nx] - 2.0 * x[idx] + x[idx - nx]) / ctx->dy2 + + (x[idx + stride_z] + x[idx - stride_z] - 2.0 * x[idx]) * ctx->inv_dz2; + r[idx] = -rhs[idx] + lap; + } + } + } +} + +/** z = diag_inv * r (Jacobi preconditioner, SIMD) */ +static inline void SIMD_FUNC(gmres_precond)(const double* r, double* z, + double diag_inv, + size_t nx, size_t ny, + size_t k_start, size_t k_end, + size_t stride_z) { + SIMD_VEC diag_vec = SIMD_SET1(diag_inv); + int ny_int = bicgstab_size_to_int(ny); + if (ny_int == 0) return; + + for (size_t k = k_start; k < k_end; k++) { + int jj; + #pragma omp parallel for schedule(static) + for (jj = 1; jj < ny_int - 1; jj++) { + size_t j = (size_t)jj; + size_t i = 1; + for (; i + SIMD_WIDTH - 1 < nx - 1; i += SIMD_WIDTH) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + SIMD_STORE(&z[idx], SIMD_MUL(diag_vec, SIMD_LOAD(&r[idx]))); + } + for (; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + z[idx] = diag_inv * r[idx]; + } + } + } +} + +//============================================================================= +// DENSE SCALAR HELPERS (identical to scalar backend; not vectorized) +//============================================================================= + +static void SIMD_FUNC(gmres_givens)(double* H, double* cs, double* sn, + double* g, int j, int m) { + const int lead = m + 1; + for (int i = 0; i < j; i++) { + double temp = cs[i] * H[i + j * lead] + sn[i] * H[(i + 1) + j * lead]; + H[(i + 1) + j * lead] = -sn[i] * H[i + j * lead] + cs[i] * H[(i + 1) + j * lead]; + H[i + j * lead] = temp; + } + + double h_jj = H[j + j * lead]; + double h_j1j = H[(j + 1) + j * lead]; + double denom = sqrt(h_jj * h_jj + h_j1j * h_j1j); + if (denom < GMRES_BREAKDOWN_THRESHOLD) { + cs[j] = 1.0; + sn[j] = 0.0; + } else { + cs[j] = h_jj / denom; + sn[j] = h_j1j / denom; + } + + H[j + j * lead] = cs[j] * h_jj + sn[j] * h_j1j; + H[(j + 1) + j * lead] = 0.0; + + double g_j = g[j]; + g[j] = cs[j] * g_j; + g[j + 1] = -sn[j] * g_j; +} + +static void SIMD_FUNC(gmres_backsub)(const double* H, const double* g, double* y, + int k, int m) { + const int lead = m + 1; + for (int i = k - 1; i >= 0; i--) { + double sum = g[i]; + for (int l = i + 1; l < k; l++) { + sum -= H[i + l * lead] * y[l]; + } + double diag = H[i + i * lead]; + y[i] = (fabs(diag) > GMRES_BREAKDOWN_THRESHOLD) ? (sum / diag) : 0.0; + } +} + +//============================================================================= +// INIT / DESTROY +//============================================================================= + +static void SIMD_FUNC(gmres_destroy)(poisson_solver_t* solver) { + if (!solver || !solver->context) return; + gmres_simd_context_t* ctx = (gmres_simd_context_t*)(solver->context); + cfd_aligned_free(ctx->Vblock); + cfd_aligned_free(ctx->w); + cfd_aligned_free(ctx->Mv); + cfd_free(ctx->H); + cfd_free(ctx->cs); + cfd_free(ctx->sn); + cfd_free(ctx->g); + cfd_free(ctx->y); + cfd_aligned_free(ctx); + solver->context = NULL; +} + +static cfd_status_t SIMD_FUNC(gmres_init)( + poisson_solver_t* solver, + size_t nx, size_t ny, size_t nz, + double dx, double dy, double dz, + const poisson_solver_params_t* params) { + + gmres_simd_context_t* ctx = (gmres_simd_context_t*)cfd_aligned_calloc( + 1, sizeof(gmres_simd_context_t)); + if (!ctx) { + return CFD_ERROR_NOMEM; + } + + ctx->dx2 = dx * dx; + ctx->dy2 = dy * dy; + ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz); + poisson_solver_compute_3d_bounds(nz, nx, ny, &ctx->stride_z, &ctx->k_start, &ctx->k_end); + ctx->diag_inv = 1.0 / (2.0 / ctx->dx2 + 2.0 / ctx->dy2 + 2.0 * ctx->inv_dz2); + ctx->use_precond = (params && params->preconditioner == POISSON_PRECOND_JACOBI); + + int m = (params && params->restart > 0) ? params->restart : GMRES_DEFAULT_RESTART; + if (m < 1) m = 1; + ctx->m = m; + + size_t n = nx * ny * nz; + ctx->n = n; + + ctx->Vblock = (double*)cfd_aligned_calloc((size_t)(m + 1) * n, sizeof(double)); + ctx->w = (double*)cfd_aligned_calloc(n, sizeof(double)); + ctx->Mv = ctx->use_precond ? (double*)cfd_aligned_calloc(n, sizeof(double)) : NULL; + + ctx->H = (double*)cfd_calloc((size_t)(m + 1) * m, sizeof(double)); + ctx->cs = (double*)cfd_calloc(m, sizeof(double)); + ctx->sn = (double*)cfd_calloc(m, sizeof(double)); + ctx->g = (double*)cfd_calloc((size_t)m + 1, sizeof(double)); + ctx->y = (double*)cfd_calloc(m, sizeof(double)); + + int precond_ok = (!ctx->use_precond) || (ctx->Mv != NULL); + if (!ctx->Vblock || !ctx->w || !ctx->H || !ctx->cs || !ctx->sn || + !ctx->g || !ctx->y || !precond_ok) { + solver->context = ctx; + SIMD_FUNC(gmres_destroy)(solver); + return CFD_ERROR_NOMEM; + } + + ctx->dx2_inv_vec = SIMD_SET1(1.0 / ctx->dx2); + ctx->dy2_inv_vec = SIMD_SET1(1.0 / ctx->dy2); + ctx->dz2_inv_vec = SIMD_SET1(ctx->inv_dz2); + ctx->two_vec = SIMD_SET1(2.0); + + ctx->initialized = 1; + solver->context = ctx; + return CFD_SUCCESS; +} + +//============================================================================= +// SOLVE +//============================================================================= + +static cfd_status_t SIMD_FUNC(gmres_solve)( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + poisson_solver_stats_t* stats) { + (void)x_temp; + + gmres_simd_context_t* ctx = (gmres_simd_context_t*)(solver->context); + if (!ctx || !ctx->initialized) { + return CFD_ERROR; + } + + size_t nx = solver->nx; + size_t ny = solver->ny; + size_t stride_z = ctx->stride_z; + size_t k_start = ctx->k_start; + size_t k_end = ctx->k_end; + + int m = ctx->m; + size_t n = ctx->n; + const int lead = m + 1; + + double* Vblock = ctx->Vblock; + double* w = ctx->w; + double* Mv = ctx->Mv; + double* H = ctx->H; + double* cs = ctx->cs; + double* sn = ctx->sn; + double* g = ctx->g; + double* y = ctx->y; + int use_precond = ctx->use_precond; + double diag_inv = ctx->diag_inv; + + poisson_solver_params_t* params = &solver->params; + double start_time = poisson_solver_get_time_ms(); + + poisson_solver_apply_bc(solver, x); + + SIMD_FUNC(gmres_residual)(x, rhs, Vblock, nx, ny, ctx); + double beta = SIMD_FUNC(gmres_norm)(Vblock, nx, ny, k_start, k_end, stride_z); + double initial_res = beta; + + if (stats) { + stats->initial_residual = initial_res; + } + + double tolerance = params->tolerance * initial_res; + if (tolerance < params->absolute_tolerance) { + tolerance = params->absolute_tolerance; + } + + if (initial_res < params->absolute_tolerance) { + poisson_solver_apply_bc(solver, x); + if (stats) { + stats->status = POISSON_CONVERGED; + stats->iterations = 0; + stats->final_residual = initial_res; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return CFD_SUCCESS; + } + + int total_inner = 0; + int converged = 0; + double final_res = initial_res; + + while (total_inner < params->max_iterations && !converged) { + if (beta < tolerance || beta < params->absolute_tolerance) { + converged = 1; + final_res = beta; + break; + } + + SIMD_FUNC(gmres_scale)(1.0 / beta, Vblock, nx, ny, k_start, k_end, stride_z); + memset(g, 0, (size_t)(m + 1) * sizeof(double)); + g[0] = beta; + + int k = 0; + for (int j = 0; j < m; j++) { + double* v_j = Vblock + (size_t)j * n; + + if (use_precond) { + SIMD_FUNC(gmres_precond)(v_j, Mv, diag_inv, nx, ny, k_start, k_end, stride_z); + SIMD_FUNC(gmres_apply_A)(Mv, w, nx, ny, ctx); + } else { + SIMD_FUNC(gmres_apply_A)(v_j, w, nx, ny, ctx); + } + + for (int i = 0; i <= j; i++) { + double* v_i = Vblock + (size_t)i * n; + double hij = SIMD_FUNC(gmres_dot)(w, v_i, nx, ny, k_start, k_end, stride_z); + H[i + j * lead] = hij; + SIMD_FUNC(gmres_axpy)(-hij, v_i, w, nx, ny, k_start, k_end, stride_z); + } + + double hnorm = SIMD_FUNC(gmres_norm)(w, nx, ny, k_start, k_end, stride_z); + H[(j + 1) + j * lead] = hnorm; + + int happy = (hnorm < GMRES_BREAKDOWN_THRESHOLD); + if (!happy) { + SIMD_FUNC(gmres_scale)(1.0 / hnorm, w, nx, ny, k_start, k_end, stride_z); + SIMD_FUNC(gmres_copy)(w, Vblock + (size_t)(j + 1) * n, nx, ny, k_start, k_end, stride_z); + } + + SIMD_FUNC(gmres_givens)(H, cs, sn, g, j, m); + total_inner++; + k = j + 1; + + double resid_est = fabs(g[j + 1]); + if (resid_est < tolerance || happy || total_inner >= params->max_iterations) { + break; + } + } + + SIMD_FUNC(gmres_backsub)(H, g, y, k, m); + if (use_precond) { + memset(w, 0, n * sizeof(double)); + for (int jj = 0; jj < k; jj++) { + SIMD_FUNC(gmres_axpy)(y[jj], Vblock + (size_t)jj * n, w, + nx, ny, k_start, k_end, stride_z); + } + SIMD_FUNC(gmres_precond)(w, Mv, diag_inv, nx, ny, k_start, k_end, stride_z); + SIMD_FUNC(gmres_axpy)(1.0, Mv, x, nx, ny, k_start, k_end, stride_z); + } else { + for (int jj = 0; jj < k; jj++) { + SIMD_FUNC(gmres_axpy)(y[jj], Vblock + (size_t)jj * n, x, + nx, ny, k_start, k_end, stride_z); + } + } + + SIMD_FUNC(gmres_residual)(x, rhs, Vblock, nx, ny, ctx); + beta = SIMD_FUNC(gmres_norm)(Vblock, nx, ny, k_start, k_end, stride_z); + final_res = beta; + if (beta < tolerance || beta < params->absolute_tolerance) { + converged = 1; + } + } + + poisson_solver_apply_bc(solver, x); + + if (stats) { + stats->iterations = total_inner; + stats->final_residual = final_res; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + stats->status = converged ? POISSON_CONVERGED : POISSON_MAX_ITER; + } + + return converged ? CFD_SUCCESS : CFD_ERROR_MAX_ITER; +} + +static cfd_status_t SIMD_FUNC(gmres_iterate)( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + double* residual) { + (void)x_temp; + if (residual) { + *residual = poisson_solver_compute_residual(solver, x, rhs); + } + return CFD_SUCCESS; +} + +//============================================================================= +// FACTORY +//============================================================================= + +poisson_solver_t* FACTORY_NAME(SIMD_SUFFIX)(void) { + poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t)); + if (!solver) return NULL; + + solver->name = POISSON_SOLVER_TYPE_GMRES_SIMD; + solver->description = "Restarted GMRES(m) (SIMD + OpenMP)"; + solver->method = POISSON_METHOD_GMRES; + solver->backend = POISSON_BACKEND_SIMD; + solver->params = poisson_solver_params_default(); + + solver->init = SIMD_FUNC(gmres_init); + solver->destroy = SIMD_FUNC(gmres_destroy); + solver->solve = SIMD_FUNC(gmres_solve); + solver->iterate = SIMD_FUNC(gmres_iterate); + solver->apply_bc = NULL; + + return solver; +} + +//============================================================================= +// CLEANUP MACROS +//============================================================================= + +#undef gmres_simd_context_t +#undef CONCAT_IMPL +#undef CONCAT +#undef SIMD_FUNC +#undef FACTORY_NAME_IMPL +#undef FACTORY_NAME diff --git a/tests/math/test_gmres.c b/tests/math/test_gmres.c new file mode 100644 index 00000000..cb6bb753 --- /dev/null +++ b/tests/math/test_gmres.c @@ -0,0 +1,632 @@ +/** + * @file test_gmres.c + * @brief Restarted GMRES(m) solver unit tests + * + * Tests the GMRES solver on the Poisson equation. GMRES is designed for + * non-symmetric systems but must also solve the symmetric Poisson operator + * used here, converging to the same solution as CG (a strong cross-check). + * + * Tests cover: + * - Creation / initialization / restart-parameter plumbing + * - Convergence on zero RHS (happy/early-return path) + * - Convergence on sinusoidal RHS (Neumann compatible) + * - Comparison with CG (max-norm and L2-norm) on the same SPD problem + * - Manufactured Dirichlet solution accuracy (O(h^2)) + * - Hessenberg estimate vs independently recomputed true residual + * - Restart-no-stall: small m and default m reach the same solution + * - Jacobi preconditioner path converges (no iteration-count assertion) + * - Error handling (unsupported backends, NULL destroy) + */ + +#include "unity.h" +#include "cfd/solvers/poisson_solver.h" +#include "cfd/core/memory.h" +#include "cfd/core/indexing.h" +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +/* ============================================================================ + * TEST PARAMETERS + * ============================================================================ */ + +#define DOMAIN_XMIN 0.0 +#define DOMAIN_XMAX 1.0 +#define DOMAIN_YMIN 0.0 +#define DOMAIN_YMAX 1.0 + +#define NX_SMALL 17 +#define NY_SMALL 17 +#define NX_MEDIUM 33 +#define NY_MEDIUM 33 + +#define TOLERANCE 1e-6 +#define ABS_TOLERANCE 1e-10 +#define MAX_ITERATIONS 2000 + +/* ============================================================================ + * TEST FIXTURES + * ============================================================================ */ + +void setUp(void) {} +void tearDown(void) {} + +/* ============================================================================ + * HELPER FUNCTIONS (mirror test_bicgstab.c) + * ============================================================================ */ + +static double* create_field(size_t nx, size_t ny) { + double* field = (double*)cfd_malloc(nx * ny * sizeof(double)); + if (field) { + for (size_t i = 0; i < nx * ny; i++) { + field[i] = 0.0; + } + } + return field; +} + +/** Sinusoidal RHS f = cos(2πx)cos(2πy) with interior mean removed (Neumann compatible) */ +static void init_sinusoidal_rhs(double* rhs, size_t nx, size_t ny, + double dx, double dy) { + for (size_t j = 0; j < ny; j++) { + double y = DOMAIN_YMIN + j * dy; + for (size_t i = 0; i < nx; i++) { + double x = DOMAIN_XMIN + i * dx; + rhs[IDX_2D(i, j, nx)] = cos(2.0 * M_PI * x) * cos(2.0 * M_PI * y); + } + } + + double interior_sum = 0.0; + size_t interior_count = 0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + interior_sum += rhs[IDX_2D(i, j, nx)]; + interior_count++; + } + } + double interior_mean = interior_sum / (double)interior_count; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + rhs[IDX_2D(i, j, nx)] -= interior_mean; + } + } +} + +/** RMS difference between two fields over interior points */ +static double compute_l2_error(const double* a, const double* b, + size_t nx, size_t ny) { + double sum = 0.0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + double diff = a[IDX_2D(i, j, nx)] - b[IDX_2D(i, j, nx)]; + sum += diff * diff; + } + } + return sqrt(sum / ((nx - 2) * (ny - 2))); +} + +/** RHS for manufactured p = sin(πx)sin(πy): ∇²p = -2π² sin(πx)sin(πy) */ +static void init_dirichlet_rhs(double* rhs, size_t nx, size_t ny, + double dx, double dy) { + (void)dx; (void)dy; + double coeff = -2.0 * M_PI * M_PI; + for (size_t j = 0; j < ny; j++) { + double y = DOMAIN_YMIN + j * (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + for (size_t i = 0; i < nx; i++) { + double x = DOMAIN_XMIN + i * (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + rhs[IDX_2D(i, j, nx)] = coeff * sin(M_PI * x) * sin(M_PI * y); + } + } +} + +static void compute_exact_solution(double* exact, size_t nx, size_t ny) { + for (size_t j = 0; j < ny; j++) { + double y = DOMAIN_YMIN + j * (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + for (size_t i = 0; i < nx; i++) { + double x = DOMAIN_XMIN + i * (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + exact[IDX_2D(i, j, nx)] = sin(M_PI * x) * sin(M_PI * y); + } + } +} + +static void remove_interior_mean(double* field, size_t nx, size_t ny) { + double sum = 0.0; + size_t count = 0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + sum += field[IDX_2D(i, j, nx)]; + count++; + } + } + double mean = sum / count; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + field[IDX_2D(i, j, nx)] -= mean; + } + } +} + +/** + * Independently recompute the true residual max-norm ||rhs - ∇²x||_inf. + * + * Measured over the DEEP interior (indices 2..n-3), skipping the first interior + * ring. The solver applies output boundary conditions to x after the solve, which + * rewrites the boundary cells; that changes the stencil residual only at the ring + * i==1 / j==1. The deep interior residual is exactly the solve-consistent residual + * the solver reports, independently recomputed here without the Hessenberg estimate. + */ +static double true_residual_inf(const double* x, const double* rhs, + size_t nx, size_t ny, double dx, double dy) { + double dx2 = dx * dx; + double dy2 = dy * dy; + double max_r = 0.0; + for (size_t j = 2; j < ny - 2; j++) { + for (size_t i = 2; i < nx - 2; i++) { + size_t idx = IDX_2D(i, j, nx); + double lap = (x[idx + 1] - 2.0 * x[idx] + x[idx - 1]) / dx2 + + (x[idx + nx] - 2.0 * x[idx] + x[idx - nx]) / dy2; + double r = fabs(lap - rhs[idx]); + if (r > max_r) max_r = r; + } + } + return max_r; +} + +/* ============================================================================ + * BASIC TESTS + * ============================================================================ */ + +void test_gmres_create(void) { + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + + TEST_ASSERT_NOT_NULL(solver); + TEST_ASSERT_EQUAL(POISSON_METHOD_GMRES, solver->method); + TEST_ASSERT_EQUAL(POISSON_BACKEND_SCALAR, solver->backend); + TEST_ASSERT_EQUAL_STRING("gmres_scalar", solver->name); + + poisson_solver_destroy(solver); +} + +void test_gmres_init(void) { + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + + size_t nx = NX_SMALL, ny = NY_SMALL; + double dx = (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + double dy = (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + + cfd_status_t status = poisson_solver_init(solver, nx, ny, 1, dx, dy, 0.0, NULL); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(nx, solver->nx); + TEST_ASSERT_EQUAL(ny, solver->ny); + + poisson_solver_destroy(solver); +} + +/* ============================================================================ + * CONVERGENCE TESTS + * ============================================================================ */ + +/** + * Zero RHS with non-zero initial guess. The initial residual is tiny (the guess + * is a constant, whose discrete Laplacian is zero on the interior), so this + * exercises the already-converged early-return / happy path. Must NOT error. + */ +void test_gmres_zero_rhs(void) { + size_t nx = NX_SMALL, ny = NY_SMALL; + double dx = (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + double dy = (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + + double* p = create_field(nx, ny); + double* rhs = create_field(nx, ny); + TEST_ASSERT_NOT_NULL(p); + TEST_ASSERT_NOT_NULL(rhs); + + for (size_t i = 0; i < nx * ny; i++) { + p[i] = 1.0; + } + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = MAX_ITERATIONS; + + cfd_status_t status = poisson_solver_init(solver, nx, ny, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + status = poisson_solver_solve(solver, p, NULL, rhs, &stats); + + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + TEST_ASSERT_TRUE(stats.iterations < MAX_ITERATIONS); + + poisson_solver_destroy(solver); + cfd_free(p); + cfd_free(rhs); +} + +void test_gmres_sinusoidal_rhs(void) { + size_t nx = NX_MEDIUM, ny = NY_MEDIUM; + double dx = (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + double dy = (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + + double* p = create_field(nx, ny); + double* rhs = create_field(nx, ny); + TEST_ASSERT_NOT_NULL(p); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, nx, ny, dx, dy); + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = MAX_ITERATIONS; + + cfd_status_t status = poisson_solver_init(solver, nx, ny, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + status = poisson_solver_solve(solver, p, NULL, rhs, &stats); + + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + + double relative_tol = stats.initial_residual * params.tolerance; + TEST_ASSERT_TRUE(stats.final_residual < relative_tol || + stats.final_residual < params.absolute_tolerance); + + poisson_solver_destroy(solver); + cfd_free(p); + cfd_free(rhs); +} + +/** + * GMRES stats.final_residual must match an independently recomputed true + * residual. Right preconditioning is chosen precisely so the reported residual + * is the true unpreconditioned one. Here (no precond) they must agree closely. + */ +void test_gmres_residual_matches_true(void) { + size_t nx = NX_MEDIUM, ny = NY_MEDIUM; + double dx = (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + double dy = (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + + double* p = create_field(nx, ny); + double* rhs = create_field(nx, ny); + TEST_ASSERT_NOT_NULL(p); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, nx, ny, dx, dy); + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = MAX_ITERATIONS; + poisson_solver_init(solver, nx, ny, 1, dx, dy, 0.0, ¶ms); + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + cfd_status_t status = poisson_solver_solve(solver, p, NULL, rhs, &stats); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + + /* Independently recomputed max-norm residual must be small... */ + double true_res = true_residual_inf(p, rhs, nx, ny, dx, dy); + double relative_tol = stats.initial_residual * params.tolerance; + TEST_ASSERT_TRUE(true_res < 10.0 * relative_tol || + true_res < 1e-6); + + poisson_solver_destroy(solver); + cfd_free(p); + cfd_free(rhs); +} + +/** + * GMRES vs CG on the SPD Poisson problem: same solution up to a Neumann + * constant offset. Max-norm difference cross-check. + */ +void test_gmres_vs_cg(void) { + size_t nx = NX_SMALL, ny = NY_SMALL; + double dx = (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + double dy = (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + + double* p_gmres = create_field(nx, ny); + double* p_cg = create_field(nx, ny); + double* rhs = create_field(nx, ny); + TEST_ASSERT_NOT_NULL(p_gmres); + TEST_ASSERT_NOT_NULL(p_cg); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, nx, ny, dx, dy); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = MAX_ITERATIONS; + + poisson_solver_t* gmres = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(gmres); + poisson_solver_init(gmres, nx, ny, 1, dx, dy, 0.0, ¶ms); + poisson_solver_stats_t stats_gmres = poisson_solver_stats_default(); + cfd_status_t status = poisson_solver_solve(gmres, p_gmres, NULL, rhs, &stats_gmres); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_t* cg = poisson_solver_create( + POISSON_METHOD_CG, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(cg); + poisson_solver_init(cg, nx, ny, 1, dx, dy, 0.0, ¶ms); + poisson_solver_stats_t stats_cg = poisson_solver_stats_default(); + status = poisson_solver_solve(cg, p_cg, NULL, rhs, &stats_cg); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_gmres.status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_cg.status); + + remove_interior_mean(p_gmres, nx, ny); + remove_interior_mean(p_cg, nx, ny); + + double max_diff = 0.0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + double diff = fabs(p_gmres[IDX_2D(i, j, nx)] - p_cg[IDX_2D(i, j, nx)]); + if (diff > max_diff) max_diff = diff; + } + } + + TEST_ASSERT_TRUE(max_diff < 1e-4); + + poisson_solver_destroy(gmres); + poisson_solver_destroy(cg); + cfd_free(p_gmres); + cfd_free(p_cg); + cfd_free(rhs); +} + +/** GMRES vs CG L2 (RMS) cross-check on the medium grid */ +void test_gmres_vs_cg_l2(void) { + size_t nx = NX_MEDIUM, ny = NY_MEDIUM; + double dx = (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + double dy = (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + + double* p_gmres = create_field(nx, ny); + double* p_cg = create_field(nx, ny); + double* rhs = create_field(nx, ny); + TEST_ASSERT_NOT_NULL(p_gmres); + TEST_ASSERT_NOT_NULL(p_cg); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, nx, ny, dx, dy); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = MAX_ITERATIONS; + + poisson_solver_t* gmres = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + poisson_solver_init(gmres, nx, ny, 1, dx, dy, 0.0, ¶ms); + poisson_solver_stats_t s1 = poisson_solver_stats_default(); + poisson_solver_solve(gmres, p_gmres, NULL, rhs, &s1); + + poisson_solver_t* cg = poisson_solver_create( + POISSON_METHOD_CG, POISSON_BACKEND_SCALAR); + poisson_solver_init(cg, nx, ny, 1, dx, dy, 0.0, ¶ms); + poisson_solver_stats_t s2 = poisson_solver_stats_default(); + poisson_solver_solve(cg, p_cg, NULL, rhs, &s2); + + remove_interior_mean(p_gmres, nx, ny); + remove_interior_mean(p_cg, nx, ny); + + double l2_error = compute_l2_error(p_gmres, p_cg, nx, ny); + TEST_ASSERT_TRUE(l2_error < 1e-5); + + poisson_solver_destroy(gmres); + poisson_solver_destroy(cg); + cfd_free(p_gmres); + cfd_free(p_cg); + cfd_free(rhs); +} + +/** Manufactured Dirichlet solution p = sin(πx)sin(πy): O(h²) accuracy */ +void test_gmres_dirichlet(void) { + size_t nx = NX_MEDIUM, ny = NY_MEDIUM; + double dx = (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + double dy = (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + + double* p = create_field(nx, ny); + double* rhs = create_field(nx, ny); + double* exact = create_field(nx, ny); + TEST_ASSERT_NOT_NULL(p); + TEST_ASSERT_NOT_NULL(rhs); + TEST_ASSERT_NOT_NULL(exact); + + init_dirichlet_rhs(rhs, nx, ny, dx, dy); + compute_exact_solution(exact, nx, ny); + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = MAX_ITERATIONS; + + cfd_status_t status = poisson_solver_init(solver, nx, ny, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + status = poisson_solver_solve(solver, p, NULL, rhs, &stats); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + + remove_interior_mean(p, nx, ny); + remove_interior_mean(exact, nx, ny); + + double l2_error = compute_l2_error(p, exact, nx, ny); + TEST_ASSERT_TRUE(l2_error < 0.01); + + poisson_solver_destroy(solver); + cfd_free(p); + cfd_free(rhs); + cfd_free(exact); +} + +/** + * Restart must not stall: GMRES(5) (many restarts) and GMRES(30, default) must + * reach the same mean-removed solution, and the small-m run must finish in a + * finite number of iterations below the cap. + */ +void test_gmres_restart_no_stall(void) { + size_t nx = NX_MEDIUM, ny = NY_MEDIUM; + double dx = (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + double dy = (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + + double* p_small = create_field(nx, ny); + double* p_big = create_field(nx, ny); + double* rhs = create_field(nx, ny); + TEST_ASSERT_NOT_NULL(p_small); + TEST_ASSERT_NOT_NULL(p_big); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, nx, ny, dx, dy); + + /* Small restart m = 5 */ + poisson_solver_params_t params_small = poisson_solver_params_default(); + params_small.tolerance = TOLERANCE; + params_small.max_iterations = MAX_ITERATIONS; + params_small.restart = 5; + + poisson_solver_t* g_small = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + poisson_solver_init(g_small, nx, ny, 1, dx, dy, 0.0, ¶ms_small); + poisson_solver_stats_t stats_small = poisson_solver_stats_default(); + cfd_status_t status = poisson_solver_solve(g_small, p_small, NULL, rhs, &stats_small); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_small.status); + TEST_ASSERT_TRUE(stats_small.iterations < MAX_ITERATIONS); + + /* Default restart (m = 30) */ + poisson_solver_params_t params_big = poisson_solver_params_default(); + params_big.tolerance = TOLERANCE; + params_big.max_iterations = MAX_ITERATIONS; + params_big.restart = 0; /* auto -> 30 */ + + poisson_solver_t* g_big = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + poisson_solver_init(g_big, nx, ny, 1, dx, dy, 0.0, ¶ms_big); + poisson_solver_stats_t stats_big = poisson_solver_stats_default(); + status = poisson_solver_solve(g_big, p_big, NULL, rhs, &stats_big); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_big.status); + + remove_interior_mean(p_small, nx, ny); + remove_interior_mean(p_big, nx, ny); + + double l2_error = compute_l2_error(p_small, p_big, nx, ny); + TEST_ASSERT_TRUE(l2_error < 1e-6); + + poisson_solver_destroy(g_small); + poisson_solver_destroy(g_big); + cfd_free(p_small); + cfd_free(p_big); + cfd_free(rhs); +} + +/** + * Jacobi preconditioner path must converge. On a uniform grid the Laplacian + * diagonal is constant, so Jacobi preconditioning is a scalar multiply and does + * NOT reduce iteration count — we assert convergence only, not fewer iterations. + */ +void test_gmres_jacobi_precond(void) { + size_t nx = NX_MEDIUM, ny = NY_MEDIUM; + double dx = (DOMAIN_XMAX - DOMAIN_XMIN) / (nx - 1); + double dy = (DOMAIN_YMAX - DOMAIN_YMIN) / (ny - 1); + + double* p = create_field(nx, ny); + double* rhs = create_field(nx, ny); + TEST_ASSERT_NOT_NULL(p); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, nx, ny, dx, dy); + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = MAX_ITERATIONS; + params.preconditioner = POISSON_PRECOND_JACOBI; + + cfd_status_t status = poisson_solver_init(solver, nx, ny, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + status = poisson_solver_solve(solver, p, NULL, rhs, &stats); + + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + + double relative_tol = stats.initial_residual * params.tolerance; + TEST_ASSERT_TRUE(stats.final_residual < relative_tol || + stats.final_residual < params.absolute_tolerance); + + poisson_solver_destroy(solver); + cfd_free(p); + cfd_free(rhs); +} + +/* ============================================================================ + * ERROR HANDLING TESTS + * ============================================================================ */ + +void test_gmres_unsupported_backend(void) { + /* GPU backend: create() returns the factory solver whenever the library was + * built with CUDA (device check happens in init). GMRES has no GPU backend + * yet, so it must return NULL regardless. */ + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_GPU); + TEST_ASSERT_NULL(solver); +} + +void test_gmres_destroy_null(void) { + poisson_solver_destroy(NULL); +} + +/* ============================================================================ + * MAIN + * ============================================================================ */ + +int main(void) { + UNITY_BEGIN(); + + RUN_TEST(test_gmres_create); + RUN_TEST(test_gmres_init); + + RUN_TEST(test_gmres_zero_rhs); + RUN_TEST(test_gmres_sinusoidal_rhs); + RUN_TEST(test_gmres_residual_matches_true); + RUN_TEST(test_gmres_vs_cg); + RUN_TEST(test_gmres_vs_cg_l2); + RUN_TEST(test_gmres_dirichlet); + RUN_TEST(test_gmres_restart_no_stall); + RUN_TEST(test_gmres_jacobi_precond); + + RUN_TEST(test_gmres_unsupported_backend); + RUN_TEST(test_gmres_destroy_null); + + return UNITY_END(); +} diff --git a/tests/math/test_gmres_avx2.c b/tests/math/test_gmres_avx2.c new file mode 100644 index 00000000..27a01c33 --- /dev/null +++ b/tests/math/test_gmres_avx2.c @@ -0,0 +1,186 @@ +/** + * @file test_gmres_avx2.c + * @brief Consistency test: GMRES AVX2 vs Scalar + * + * Verifies the AVX2-optimized GMRES(m) solver produces numerically consistent + * results with the scalar reference: + * - L2 difference < 1e-7 between solutions + * - Iteration counts within ±2 (SIMD rounding accumulation) + * + * The dense Givens/Hessenberg/back-substitution work is byte-identical scalar in + * both backends, so the only divergence source is the vectorized O(n) primitives. + */ + +#include "unity.h" +#include "cfd/solvers/poisson_solver.h" +#include "cfd/core/memory.h" +#include "cfd/core/cpu_features.h" +#include "cfd/core/indexing.h" +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#define NX 33 +#define NY 33 +#define XMIN 0.0 +#define XMAX 1.0 +#define YMIN 0.0 +#define YMAX 1.0 +#define TOLERANCE 1e-6 + +void setUp(void) {} +void tearDown(void) {} + +static void init_sinusoidal_rhs(double* rhs, size_t nx, size_t ny, + double dx, double dy) { + for (size_t j = 0; j < ny; j++) { + double y = YMIN + j * dy; + for (size_t i = 0; i < nx; i++) { + double x = XMIN + i * dx; + rhs[IDX_2D(i, j, nx)] = cos(2.0 * M_PI * x) * cos(2.0 * M_PI * y); + } + } + + double interior_sum = 0.0; + size_t interior_count = 0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + interior_sum += rhs[IDX_2D(i, j, nx)]; + interior_count++; + } + } + if (interior_count > 0) { + double interior_mean = interior_sum / (double)interior_count; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + rhs[IDX_2D(i, j, nx)] -= interior_mean; + } + } + } + for (size_t i = 0; i < nx; i++) { + rhs[i] = 0.0; + rhs[(ny - 1) * nx + i] = 0.0; + } + for (size_t j = 0; j < ny; j++) { + rhs[j * nx] = 0.0; + rhs[j * nx + (nx - 1)] = 0.0; + } +} + +void test_gmres_avx2_scalar_consistency(void) { + double dx = (XMAX - XMIN) / (NX - 1); + double dy = (YMAX - YMIN) / (NY - 1); + size_t n = NX * NY; + + double* x_scalar = (double*)cfd_calloc(n, sizeof(double)); + double* x_avx2 = (double*)cfd_calloc(n, sizeof(double)); + double* x_temp = (double*)cfd_calloc(n, sizeof(double)); + double* rhs = (double*)cfd_calloc(n, sizeof(double)); + TEST_ASSERT_NOT_NULL(x_scalar); + TEST_ASSERT_NOT_NULL(x_avx2); + TEST_ASSERT_NOT_NULL(x_temp); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, NX, NY, dx, dy); + + poisson_solver_t* solver_scalar = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver_scalar); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = 5000; + + cfd_status_t status = poisson_solver_init(solver_scalar, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_scalar = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_scalar, x_scalar, x_temp, rhs, &stats_scalar); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_scalar.status); + + /* GMRES SIMD backends require OpenMP - skip if not available */ +#ifndef CFD_ENABLE_OPENMP + cfd_free(x_scalar); + cfd_free(x_avx2); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("GMRES SIMD requires OpenMP (not enabled in this build)"); + return; +#endif + + bool simd_available = poisson_solver_backend_available(POISSON_BACKEND_SIMD); + + cfd_simd_arch_t arch = cfd_detect_simd_arch(); + if (simd_available && arch != CFD_SIMD_AVX2) { + cfd_free(x_scalar); + cfd_free(x_avx2); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("AVX2 test skipped: platform uses different SIMD architecture"); + return; + } + + poisson_solver_t* solver_avx2 = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SIMD); + + if (!solver_avx2) { + cfd_free(x_scalar); + cfd_free(x_avx2); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + if (simd_available) { + TEST_FAIL_MESSAGE("SIMD backend available but GMRES SIMD solver creation failed"); + } else { + TEST_IGNORE_MESSAGE("SIMD backend not available on this platform"); + } + return; + } + + status = poisson_solver_init(solver_avx2, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_avx2 = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_avx2, x_avx2, x_temp, rhs, &stats_avx2); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_avx2.status); + + double l2_diff = 0.0; + size_t count = 0; + for (size_t j = 1; j < NY - 1; j++) { + for (size_t i = 1; i < NX - 1; i++) { + size_t idx = IDX_2D(i, j, NX); + double diff = x_scalar[idx] - x_avx2[idx]; + l2_diff += diff * diff; + count++; + } + } + l2_diff = sqrt(l2_diff / count); + + /* SIMD FMA + reduction ordering causes minor rounding vs scalar; 1e-7 stays + * an order of magnitude below the 1e-6 solver tolerance while leaving + * cross-platform margin (mirrors the BiCGSTAB consistency test). */ + TEST_ASSERT_DOUBLE_WITHIN(1.0e-7, 0.0, l2_diff); + + int iter_diff = abs((int)stats_scalar.iterations - (int)stats_avx2.iterations); + TEST_ASSERT_LESS_OR_EQUAL(2, iter_diff); + + poisson_solver_destroy(solver_scalar); + poisson_solver_destroy(solver_avx2); + cfd_free(x_scalar); + cfd_free(x_avx2); + cfd_free(x_temp); + cfd_free(rhs); +} + +int main(void) { + UNITY_BEGIN(); + RUN_TEST(test_gmres_avx2_scalar_consistency); + return UNITY_END(); +} diff --git a/tests/math/test_gmres_neon.c b/tests/math/test_gmres_neon.c new file mode 100644 index 00000000..8e6e18e6 --- /dev/null +++ b/tests/math/test_gmres_neon.c @@ -0,0 +1,183 @@ +/** + * @file test_gmres_neon.c + * @brief Consistency test: GMRES ARM NEON vs Scalar + * + * Verifies the NEON-optimized GMRES(m) solver produces numerically consistent + * results with the scalar reference: + * - L2 difference < 1e-7 between solutions + * - Iteration counts within ±2 (SIMD rounding accumulation) + * + * Skips gracefully on non-ARM platforms. The dense Givens/Hessenberg/back-sub + * work is byte-identical scalar in both backends. + */ + +#include "unity.h" +#include "cfd/solvers/poisson_solver.h" +#include "cfd/core/memory.h" +#include "cfd/core/cpu_features.h" +#include "cfd/core/indexing.h" +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#define NX 33 +#define NY 33 +#define XMIN 0.0 +#define XMAX 1.0 +#define YMIN 0.0 +#define YMAX 1.0 +#define TOLERANCE 1e-6 + +void setUp(void) {} +void tearDown(void) {} + +static void init_sinusoidal_rhs(double* rhs, size_t nx, size_t ny, + double dx, double dy) { + for (size_t j = 0; j < ny; j++) { + double y = YMIN + j * dy; + for (size_t i = 0; i < nx; i++) { + double x = XMIN + i * dx; + rhs[IDX_2D(i, j, nx)] = cos(2.0 * M_PI * x) * cos(2.0 * M_PI * y); + } + } + + double interior_sum = 0.0; + size_t interior_count = 0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + interior_sum += rhs[IDX_2D(i, j, nx)]; + interior_count++; + } + } + if (interior_count > 0) { + double interior_mean = interior_sum / (double)interior_count; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + rhs[IDX_2D(i, j, nx)] -= interior_mean; + } + } + } + for (size_t i = 0; i < nx; i++) { + rhs[i] = 0.0; + rhs[(ny - 1) * nx + i] = 0.0; + } + for (size_t j = 0; j < ny; j++) { + rhs[j * nx] = 0.0; + rhs[j * nx + (nx - 1)] = 0.0; + } +} + +void test_gmres_neon_scalar_consistency(void) { + double dx = (XMAX - XMIN) / (NX - 1); + double dy = (YMAX - YMIN) / (NY - 1); + size_t n = NX * NY; + + double* x_scalar = (double*)cfd_calloc(n, sizeof(double)); + double* x_neon = (double*)cfd_calloc(n, sizeof(double)); + double* x_temp = (double*)cfd_calloc(n, sizeof(double)); + double* rhs = (double*)cfd_calloc(n, sizeof(double)); + TEST_ASSERT_NOT_NULL(x_scalar); + TEST_ASSERT_NOT_NULL(x_neon); + TEST_ASSERT_NOT_NULL(x_temp); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, NX, NY, dx, dy); + + poisson_solver_t* solver_scalar = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver_scalar); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = 5000; + + cfd_status_t status = poisson_solver_init(solver_scalar, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_scalar = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_scalar, x_scalar, x_temp, rhs, &stats_scalar); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_scalar.status); + + /* GMRES SIMD backends require OpenMP - skip if not available */ +#ifndef CFD_ENABLE_OPENMP + cfd_free(x_scalar); + cfd_free(x_neon); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("GMRES SIMD requires OpenMP (not enabled in this build)"); + return; +#endif + + bool simd_available = poisson_solver_backend_available(POISSON_BACKEND_SIMD); + + cfd_simd_arch_t arch = cfd_detect_simd_arch(); + if (simd_available && arch != CFD_SIMD_NEON) { + cfd_free(x_scalar); + cfd_free(x_neon); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("NEON test skipped: platform uses different SIMD architecture"); + return; + } + + poisson_solver_t* solver_neon = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SIMD); + + if (!solver_neon) { + cfd_free(x_scalar); + cfd_free(x_neon); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + if (simd_available) { + TEST_FAIL_MESSAGE("SIMD backend available but GMRES SIMD solver creation failed"); + } else { + TEST_IGNORE_MESSAGE("SIMD backend not available on this platform"); + } + return; + } + + status = poisson_solver_init(solver_neon, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_neon = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_neon, x_neon, x_temp, rhs, &stats_neon); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_neon.status); + + double l2_diff = 0.0; + size_t count = 0; + for (size_t j = 1; j < NY - 1; j++) { + for (size_t i = 1; i < NX - 1; i++) { + size_t idx = IDX_2D(i, j, NX); + double diff = x_scalar[idx] - x_neon[idx]; + l2_diff += diff * diff; + count++; + } + } + l2_diff = sqrt(l2_diff / count); + + TEST_ASSERT_DOUBLE_WITHIN(1.0e-7, 0.0, l2_diff); + + int iter_diff = abs((int)stats_scalar.iterations - (int)stats_neon.iterations); + TEST_ASSERT_LESS_OR_EQUAL(2, iter_diff); + + poisson_solver_destroy(solver_scalar); + poisson_solver_destroy(solver_neon); + cfd_free(x_scalar); + cfd_free(x_neon); + cfd_free(x_temp); + cfd_free(rhs); +} + +int main(void) { + UNITY_BEGIN(); + RUN_TEST(test_gmres_neon_scalar_consistency); + return UNITY_END(); +} diff --git a/tests/math/test_omp_consistency.c b/tests/math/test_omp_consistency.c index 4cca9427..76c816e7 100644 --- a/tests/math/test_omp_consistency.c +++ b/tests/math/test_omp_consistency.c @@ -299,9 +299,108 @@ void test_redblack_omp_vs_scalar(void) { cfd_free(rhs); } +/** + * Test: GMRES OMP vs Scalar Consistency + * + * Solves the same Poisson problem with both scalar and OMP GMRES solvers, then + * verifies the solutions are numerically identical. The dense Givens/Hessenberg + * work is byte-identical scalar in both backends, so only the parallel-reduction + * dot products introduce rounding differences. + */ +void test_gmres_omp_vs_scalar(void) { + double dx = (XMAX - XMIN) / (NX - 1); + double dy = (YMAX - YMIN) / (NY - 1); + size_t n = NX * NY; + + double* x_scalar = (double*)cfd_calloc(n, sizeof(double)); + double* x_omp = (double*)cfd_calloc(n, sizeof(double)); + double* x_temp = (double*)cfd_calloc(n, sizeof(double)); + double* rhs = (double*)cfd_calloc(n, sizeof(double)); + TEST_ASSERT_NOT_NULL(x_scalar); + TEST_ASSERT_NOT_NULL(x_omp); + TEST_ASSERT_NOT_NULL(x_temp); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, NX, NY, dx, dy); + + poisson_solver_t* solver_scalar = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver_scalar); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = 1000; + + cfd_status_t status = poisson_solver_init(solver_scalar, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_scalar = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_scalar, x_scalar, x_temp, rhs, &stats_scalar); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_scalar.status); + + if (!poisson_solver_backend_available(POISSON_BACKEND_OMP)) { + cfd_free(x_scalar); + cfd_free(x_omp); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("OMP backend not available on this platform"); + return; + } + + poisson_solver_t* solver_omp = poisson_solver_create( + POISSON_METHOD_GMRES, POISSON_BACKEND_OMP); + + if (!solver_omp) { + cfd_free(x_scalar); + cfd_free(x_omp); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("OMP GMRES solver creation failed; backend may be unavailable"); + return; + } + + status = poisson_solver_init(solver_omp, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_omp = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_omp, x_omp, x_temp, rhs, &stats_omp); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_omp.status); + + double l2_diff = 0.0; + size_t count = 0; + for (size_t j = 1; j < NY - 1; j++) { + for (size_t i = 1; i < NX - 1; i++) { + size_t idx = IDX_2D(i, j, NX); + double diff = x_scalar[idx] - x_omp[idx]; + l2_diff += diff * diff; + count++; + } + } + l2_diff = sqrt(l2_diff / count); + + /* GMRES OMP and scalar share byte-identical dense math; only the + * parallel-reduction dot products differ, so agreement is excellent. */ + TEST_ASSERT_DOUBLE_WITHIN(1.0e-9, 0.0, l2_diff); + + int iter_diff = abs((int)stats_scalar.iterations - (int)stats_omp.iterations); + TEST_ASSERT_LESS_OR_EQUAL(2, iter_diff); + + poisson_solver_destroy(solver_scalar); + poisson_solver_destroy(solver_omp); + cfd_free(x_scalar); + cfd_free(x_omp); + cfd_free(x_temp); + cfd_free(rhs); +} + int main(void) { UNITY_BEGIN(); RUN_TEST(test_cg_omp_vs_scalar); RUN_TEST(test_redblack_omp_vs_scalar); + RUN_TEST(test_gmres_omp_vs_scalar); return UNITY_END(); } From 5b0f8d89ed2977fd2fd64ab15a0e789b4a8f6d24 Mon Sep 17 00:00:00 2001 From: shaia Date: Mon, 6 Jul 2026 20:37:56 +0300 Subject: [PATCH 2/5] Clarify GMRES residual test docstring: max-norm smallness check, not norm-for-norm stats match MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- tests/math/test_gmres.c | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/tests/math/test_gmres.c b/tests/math/test_gmres.c index cb6bb753..d4ef222e 100644 --- a/tests/math/test_gmres.c +++ b/tests/math/test_gmres.c @@ -294,9 +294,12 @@ void test_gmres_sinusoidal_rhs(void) { } /** - * GMRES stats.final_residual must match an independently recomputed true - * residual. Right preconditioning is chosen precisely so the reported residual - * is the true unpreconditioned one. Here (no precond) they must agree closely. + * After GMRES converges, an independently recomputed true residual must be + * small. This cross-checks that the solver's reported convergence corresponds + * to a genuinely small unpreconditioned residual (no precond here). Note the + * independent check uses the max-norm, whereas stats.final_residual is the + * 2-norm the solver measures internally, so this validates smallness rather + * than an exact norm-for-norm match. */ void test_gmres_residual_matches_true(void) { size_t nx = NX_MEDIUM, ny = NY_MEDIUM; From ef4c1061b17862eb16c32e3378f422d274aedc74 Mon Sep 17 00:00:00 2001 From: shaia Date: Mon, 6 Jul 2026 20:39:16 +0300 Subject: [PATCH 3/5] Reword GMRES OMP consistency docstring: "consistent" not "identical" 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. --- tests/math/test_omp_consistency.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/math/test_omp_consistency.c b/tests/math/test_omp_consistency.c index 76c816e7..4fbacb3d 100644 --- a/tests/math/test_omp_consistency.c +++ b/tests/math/test_omp_consistency.c @@ -303,9 +303,10 @@ void test_redblack_omp_vs_scalar(void) { * Test: GMRES OMP vs Scalar Consistency * * Solves the same Poisson problem with both scalar and OMP GMRES solvers, then - * verifies the solutions are numerically identical. The dense Givens/Hessenberg - * work is byte-identical scalar in both backends, so only the parallel-reduction - * dot products introduce rounding differences. + * verifies the solutions are numerically consistent (RMS L2 difference within + * 1e-9). The dense Givens/Hessenberg work is byte-identical scalar in both + * backends, so only the parallel-reduction dot products introduce rounding + * differences — the solutions are consistent, not bit-for-bit identical. */ void test_gmres_omp_vs_scalar(void) { double dx = (XMAX - XMIN) / (NX - 1); From b33919dfd70b97abdefcedf1582bb87613c2eb7d Mon Sep 17 00:00:00 2001 From: shaia Date: Tue, 7 Jul 2026 07:34:31 +0300 Subject: [PATCH 4/5] docs: mirror poisson_solver_method_t enum in API reference 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. --- docs/reference/api-reference.md | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/docs/reference/api-reference.md b/docs/reference/api-reference.md index 582791a9..2790c6c6 100644 --- a/docs/reference/api-reference.md +++ b/docs/reference/api-reference.md @@ -408,15 +408,19 @@ 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_GMRES = 6, // Restarted GMRES(m) (scalar/SIMD/OMP) -} 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 From b21d0feb2b5ffcac7ec47ced11fa22b880726e30 Mon Sep 17 00:00:00 2001 From: shaia Date: Tue, 7 Jul 2026 07:34:46 +0300 Subject: [PATCH 5/5] docs: sync poisson_solver_params_t snippet with public header 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. --- docs/reference/api-reference.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/reference/api-reference.md b/docs/reference/api-reference.md index 2790c6c6..088eb19d 100644 --- a/docs/reference/api-reference.md +++ b/docs/reference/api-reference.md @@ -441,10 +441,10 @@ 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;