Add initial end-to-end CUDA FGMRES solver path - #2825
Merged
Conversation
LwhJesse
marked this pull request as ready for review
June 3, 2026 05:20
pcarruscag
reviewed
Jun 13, 2026
LwhJesse
force-pushed
the
gpu/initial-cuda-fgmres
branch
from
June 17, 2026 14:24
8982fcb to
d821ca0
Compare
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
LwhJesse
force-pushed
the
gpu/initial-cuda-fgmres
branch
from
June 18, 2026 05:59
592b302 to
5b01f21
Compare
Move the FGMRES iteration into one shared implementation and select host or CUDA vector-operation backends from the existing solver entry point. Keep the CUDA path GPU-resident for SpMV, Jacobi, dot/norm, and vector updates, with only scalar reductions and final solution synchronization crossing back to the host.
LwhJesse
force-pushed
the
gpu/initial-cuda-fgmres
branch
from
June 18, 2026 11:37
5b01f21 to
a875f56
Compare
pcarruscag
reviewed
Jun 21, 2026
ThreadSanitizer hybrid_regression jobs reported a data race in CSysVector::MarkHostDataModified() on host/device validity flags. Avoid updating those CUDA-only validity flags in non-CUDA builds.
ddg93
added
GSoC
9 tasks
# Conflicts: # Common/include/linear_algebra/CSysMatrix.hpp # Common/src/linear_algebra/CSysMatrix.cpp # Common/src/linear_algebra/CSysMatrixGPU.cu # meson.build
- meson.build: fix CUDA arch for Ada Lovelace GPUs, use find_library for cudart (Ubuntu's distro-packaged CUDA toolkit installs libs outside the path meson's cuda dependency() module searches), skip the AMX-tile intrinsics header nvcc's frontend can't parse, and point nvcc's host compiler at the actual C++ compiler so MPI link flags aren't dropped when nvcc becomes the final linker for CUDA-containing targets. - Common/src/meson.build: stop forwarding CODI_REVERSE_TYPE/ CODI_FORWARD_TYPE into CUDA compilation; the GPU kernels only ever operate on plain floating-point types and nvcc's device-code frontend cannot parse CoDiPack's tape/event-system machinery. - CSysMatrixGPU.cu: fix BlockLDU_SpMV_kernel's row/col index pointer types to match the LDU struct's actual su2uint (uint32_t) storage. - CSysVector/CSysMatrix/CMatrixVectorProduct: the GPU dispatch paths (GPUDot, HtDTransfer, AssignDeviceExpression, GPUMatrixVectorProduct, Jacobi preconditioner GPU hooks) were only gated at runtime, but are only ever instantiated for non-AD scalar types. Add a su2_gpu_capable_v compile-time trait and gate every call site with if constexpr so AD-typed instantiations never reference those symbols and correctly fall back to host computation.
pcarruscag
approved these changes
Aug 2, 2026
pcarruscag
added a commit
that referenced
this pull request
Aug 15, 2026
…nditioner on GPU and CUDA Unified preconditiong matrix) (#2843) ## Proposed Changes This PR introduces improves the GPU implementation of the multiDot product between Krylov vectors used in the FGMRES solver. It moves from a looped instantiation of pairwise dot products to a batched one-time kernel launch. A custom CUDA kernel and the batched gemm using cublas<t>gemm are benchmarked, the latter being eventually selected. ## Original proposal for the GSoC This PR introduces CUDA Unified memory and Managed memory allocation, and memory management for the CSysVector class and the preconditioning matrix inside the Jacobi preconditioner. This allows for a benchmark between the two memory strategies within the FGMRES solver. Also, this draft PR extends the section of GPU execution inside the FGMRES solver Custom CUDA kernels are implemented for the preconditioning matrix, the multi dot product and the linear combination (inside the Modified Gram-Schmidt orthogonalization), and the vector norm calculation. Unary vector-scalar operations based on templates are also offloaded to GPU through a generic kernel. Moreover, abstract Syntax Tree are deployed to offload vector-vector binary operations to the GPU through a generic kernel based on a runtime evaluation of the tree. The solver logic is not modified, and the GPU path is hidden inside the specific methods. ### Unified memory approach 5fe25b8 A custom data() method recovers the CSysVector Unified pointer inside the CUDA logic, cleaning the logic from the double host/device pointers. All memory explicit memory copies are also removed from the CUDA logic, but explicit device synchronizations are introduced around MPI calls and at the end of the CUDA section. ### Managed memory approach 4b43fa0 Operators accessing the CSysVector on Host after GPU operations demand synchronization, which is introduced explicitly in each. ## Performance evaluation (on-going) For the considered test-case (rae2822), the CPU execution expresses an **Avg. s/iter: 0.198681**. The GPU execution with Unified memory expresses an **Avg. s/iter: 0.238953**. The GPU execution with Managed memory expresses an **Avg. s/iter: 0.372046**. The CPU is a Intel(R) Xeon(R) E-2276M CPU @ 2.80GHz with 12 cores. The GPU is a Quadro P620. Tests on more advanced hardware are ongoing. Profiling results comparing the Unified Memory (left column) against the Managed Memory (right column) are available in the attached pdf: [SU2_ra2822_GPU_MA_vs_UM_nsys_prof.pdf](https://github.com/user-attachments/files/30520807/SU2_ra2822_GPU_MA_vs_UM_nsys_prof.pdf) ## Asynchronous pre-fetching: The Jacobi preconditioner calculations are performed on GPU through a new custom CUDA kernel under the preconditioner abstraction. The preconditioning matrix is selected to test CUDA Unified Memory asynchronous prefetching to the GPU. For simplicity, the double CPU/GPU pointer is still maintained in the current logic, although the device pointer reduces to an alias for the Unified Memory pointer when this kind of allocation is adopted. **This strategy introduces a simple context to test the CUDA Unified Memory usage and study the possibility of overlapping memory transfers and calculations without the need to introduce CUDA streams.** Concretely, this PR: - introduces new CUDA Unified Memory allocation methods and asynchronous prefetching; - introduces the GPU logic for the Jacobi preconditioner; - introduces the GPU logic for the multi dot product and the linear combination (Modified Gram-Schmidt orthogonalization); - introduces the GPU logic for the vector norm operation; - introduces the GPU logic for generic scalar-vector unary operations through templates and generic vector/scalar-vector binary operations through Abstract Syntax Tree evaluated at runtime; - finally falls back to CUDA Managed memory as highlighted in the following discussion. This introduces the need for explicity synchronization in all the custom setter/getter methods of the CSysVector. This work is part of my ongoing contribution during the Google Summer of Code 2026 program. ## Validation Validated locally with: - serial CUDA build compilation - serial CPU build compilation - CPU/GPU numerical comparison on 1 representative case (rae2822) tested with LINEAR_SOLVER_PREC=JACOBI with both CUDA Unified and Managed memory approaches. Nsys profiling was performed to confirm the asynchronous prefetching of the CUDA Unified Memory preconditioning matrix on my local GPU. Partial prefetching is observed, although page faults were reported during the preconditioner CUDA kernel, indicating that the calculations were slowed down by the prefetching matrix still being transferred to the GPU. This overlap is expected to largely improve on more modern hardware; tests are ongoing in the cloud. ## Related Work The Jacobi preconditioner kernels come from the PR #2825. ## Observed Issues Some issues were observed during this first period of GSoC: - PR #2825 compiles with CUDA 13.3 but CUSPARSE calls raise an unknown operation at the first matrix-vector product. - the build.meson file has a hard-coded CUDA arch. - the HAVE_MPI flag is not being passed to the nvcc compiler in the develop branch. This raises a linking error if MPI operations are included within the *.cu files. That is not the case in the master branch. ## Next steps I propose to continue working on the following steps: - [X] extend the CUDA Unified Memory allocation to the CSysVector class - [X] evaluate if it might be of interest to extend the CUDA Unified Memory allocation to the CSysMatrix class; Update: discussion with mantainer indicates preference for Managed Memory approach; - [X] benchmark the multiDot approached: looped vs custom kernel batched vs cublas<t>gemm batched: batched approaches are convenient due to just one kernel launch overhead, cublas implementation shows moderate speed-up (0.01s) against the custom CUDA kernel on local P620 but ensures less code complexity, thus it is selected. ## PR Checklist - [X] I am submitting my contribution to the develop branch. - [X] My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson). - [X] My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/). - [X] I used the pre-commit hook to prevent dirty commits and used `pre-commit run --all` to format old commits. - [ ] I have added a test case that demonstrates my contribution, if necessary. - [ ] I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary. --------- Co-authored-by: Pedro Gomes <pcarruscag@gmail.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
4 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Proposed Changes
This PR adds an end-to-end CUDA linear solve path: the Krylov solvers keep their host control flow, and
CSysVectoroperations, the SpMV and the Jacobi preconditioner are dispatched to CUDA kernels whenENABLE_CUDA=YES.Transfers are explicit and owned by the object responsible for the data, with no coherency or dirty-flag tracking in
CSysVector/CSysMatrix:CSysMatrixVectorProductuploads the matrix on constructionBuild()CSysSolveuploadsbandxand downloadsxinHandleTemporariesIn/Out, which is also where device evaluation is switched on and offA solve is therefore fully device resident for Identity and Jacobi preconditioners: 2 uploads and 1 download per linear system, independent of the Krylov subspace size.
Implementation notes:
cuBLASfordot/norm, custom kernels for the block-LDU SpMV, the Jacobi apply, andCSysVectorexpression assignmentCSysVectoroperands are captured in expressions by value, so an arbitrary expression tree is trivially copyable into the assignment kernel; the required expression shapes are explicitly instantiated inCSysVectorGPU.cu, consistent with howCSysMatrixis instantiatedLINEAR_SOLVER_PREC= NONE(identity)Related Work
Follows the review direction in #2822 (show a working end-to-end GPU linear solve before splitting out infrastructure) and the implementation preferences in #2816.
Validation
JACOBI, FGMRES +NONE, FGMRES +ILU(host preconditioner path) and BCGSTAB. Results agree to the printed precision in double, and to ~6 significant figures in mixed precision. BCGSTAB agrees exactly once the linear system is converged (LINEAR_SOLVER_ERROR=1e-10); at loose tolerances the two paths diverge through BCGSTAB's own sensitivity, not a difference in the algebra.OMP_NUM_THREADS=1and4give bit-identical results on the GPU path.enable-cuda+with-omp, mixed / normal / single precision.NONE, plus one download/upload pair per preconditioner application for ILU.Earlier validation of the original design (6 representative cases,
nsys/ncuprofiling) predates the rework of the transfer and dispatch model and should be repeated.PR Checklist
pre-commit run --allto format old commits.