[GSoC] Feature: Jacobi preconditioner on GPU and CUDA Unified preconditiong matrix

[ddg93]

Proposed Changes

This draft 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 (current head) 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

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 Add initial end-to-end CUDA FGMRES solver path #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:

• perform further testing on the asynchronous prefetching of the preconditioning matrix;
• extend the CUDA Unified Memory allocation to the CSysVector class
• 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;
• I observe that CSysMatrixVectorProduct object in the FGMRES solver is the Jacobian matrix. If the Jacobian preconditioner is applied, the preconditioning matrix is also calculated from the Jacobian matrix. It would then be convenient to issue the transfer of the Jacobian matrix to the GPU inside the building of the preconditioner and perform the preconditioning matrix calculations on the GPU (sparse matrix inversion based on the current LU decomposition). Then, the preconditioning matrix would already be on the GPU, reducing by one the number of transfers from CPU to GPU. This is not the current case as both the Jacobian matrix and the preconditioning matrix are being transferred to device.
• NEW : benchmark the custom CUDA kernels (e.g. for Linear Combination) against a corresponding execution based on the Abstract Syntax Tree to understand the performance trade;
• perform further testing on the asynchronous prefetching of the preconditioning matrix;

PR Checklist

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• 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.

[bigfooted]

That's great!
Can you run pre-commit to fix the code style/formatting?

[pcarruscag]

I think Areen looked into unified memory a bit.
Do you have some performance measurements? I imagine we save on the data transfer but suffer a little on the kernel performance?
Does unified memory introduce a significant restriction in terms of the required compute capability?

[ddg93]

That's great! Can you run pre-commit to fix the code style/formatting?

Thanks for the heads-up. Updated that commit and the following ones.

[ddg93]

I think Areen looked into unified memory a bit. Do you have some performance measurements? I imagine we save on the data transfer but suffer a little on the kernel performance? Does unified memory introduce a significant restriction in terms of the required compute capability?

For the moment, I've only one performance measurement on my local GPU Quadro P620:

• GPU version based on Unified memory: Avg. s/iter: 0.256129
• CPU version: Avg. s/iter: 0.198681
  This result is clearly not sufficient but shows a decent performance of the GPU version based on Unified memory on a small test-case and on a very modest GPU.
  I intend to measure performance for a GPU version without Unified memory, and also consider larger tests and more advanced hardware (through cloud services).

Preliminary profiling shows indeed that the transfer of the preconditioning matrix is asynchronously launched but then the preconditioning kernel takes longer due to page fault. I intend to better profile what is going on with the CSysVector Unified memory implementation.

Unified memory is available for the compute capability 6.0 or higher. Some of its advanced features seem to depend on the system.

Is there some report by Areen? It could be interesting for me to read it.

[pcarruscag]

Thanks I was wondering about unified vs explicit transfers not CPU vs GPU.
Not sure if Areen registered those tests somewhere.
Based on what nvidia says about unified memory, and the our goal with the linear solver, we should transfer explicitly and avoid complexity of keeping both strategies

[ddg93]

Yes, I understand. I have further developed my PR to cover the whole FGMRES solver so that operations with CSysVectors run on GPU, operators too. In the last commit, I also fell back from Unified to Managed CUDA memory for the CSysVector. In this way, I could profile as you suggested. For the moment I could do this only on my local Quadro P620, I'll be trying to get better GPUs on the cloud. Please find here SU2_ra2822_GPU_MA_vs_UM_nsys_prof.pdf
a pdf report showing the differences between the nsys profiling output for Unified Memory (left column) (second last commit 5fe25b8) and Managed Memory (right column) (last commit 4b43fa0) for one run. UM is surprisingly competitive (should average on multiple runs and run on better GPUs to get more realistic results) and clearly facilitates the development as I did not have to safeguard the CPU manipulation of CSysVectors still on GPU.

[pcarruscag]

Thanks but I want the compile-time version of this that is being worked in https://github.com/su2code/SU2/pull/2825/changes

[ddg93]

Ok. What do you mean by that, exactly? Should I adopt the specific implementation from https://github.com/su2code/SU2/pull/2825/changes or just follow the same approach? Besides the operators kernels, my implementation differs in other parts of the FGMRES solver, especially I introduced custom CUDA kernels rather than CUBLAS.

[pcarruscag]

I will verify that pr does what it should this weekend and merge it since I had already reviewed a few iterations of it.
The multi dot and other custom kernels you have here would be nice additions on top.
Unfortunately both PRs port the jacobi preconditioner, if you have time left in your project and want to take a look at completing the lu-sgs port, that would also be high value.

[pcarruscag]

PR did what it should after I got an LLM to clean it up.
I have a weekly budget left to spend, so I'm getting the AI to try and do the ILU preconditioner.
Let me know if you want to take on the tasks I suggested above, and sorry for the merge conflicts, we can't choose who/what/when contributes to SU2.

[ddg93]

Ok to add the multi dot and other custom kernels as additions on top. I assume this meas updating my PR over the current version of the develop branch. Then, I could propose the CUBLAS or Custom kernels path for the solver based on a pre-processor flag that checks if CUBLAS is available, what do you think?

Ok to work on the lu-sgs porting, I still have some time left with the GSoC so I am happy to continue working. Besides, I am enjoying this so I'll try to stay involved even after the project has ended.
No worries about the conflicts.

[pcarruscag]

yes bring this PR up to date with develop and we can merge the multiDot part, then the lu-sgs can be a second PR