[GSoC] Feature: MultiDot and Linear Combination CUDA kernels (Originally: 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.

[ddg93]

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

For the sake of completeness, I report a benchmark between Unified Memory (left column, commit 5fe25b8 ) and Managed Memory (right column, commit 4b43fa0 ) here: : SU2_rae2822_UM_vs_MA_NVIDIARTX2000AdaGene.pdf. I was able to run this on an RTX 2000 for the rae 2822 grid with 24576 elements.

[pcarruscag]

Good direction — fusing the n*m dots into one kernel removes develop's launch-bound loop of separate cuBLAS dot calls.

Suggest dropping LinearCombinationGPU and the CSysSolve.cpp hunk. You measured no speedup, which is expected since the expression templates already fuse the combination into a single kernel. It also writes v.nElmDomain where the expression path writes nElm, so the halos of x/r go stale and the FGMRES restart is wrong on 2+ ranks, and it skips the n == 0 zeroing that LinearCombinationImpl does at the top. Both fixable, but with no perf win, removing is simpler.

Blockers in multiDot:

1. The if (VecExpr::UseDeviceExpressions()) guard is gone and the host path moved under #else. A CUDA build running ENABLE_CUDA= NO now reads device pointers that were never uploaded, with no host fallback compiled at all. The same removal makes the su2_gpu_capable_v else-branch fire unconditionally, so CUDA+AD builds abort. Please keep the runtime guard and leave the host path compiled, as in CSysMatrix.cpp:818 and ApplyPreconditionerOnHost.

2. multiDotGPU adds an SU2_OMP_MASTER and a trailing SU2_OMP_BARRIER inside the caller's BEGIN_SU2_DEVICE_REGION, which is already barrier + omp master. A barrier inside master is non-conforming and hangs at OMP_NUM_THREADS > 1 — the macro's own comment says these regions must not be nested. Drop both; only the master thread runs that body, so the statics don't need protecting either.

3. template class CSysVector<su2double>; is unguarded (develop uses #ifdef CODI_REVERSE_TYPE). In a mixed-precision build su2double == passivedouble, so it duplicates the explicit member instantiations below it. Instantiate just multiDotGPU with the neighbouring guards — the new USE_MIXED_PRECISION guard can then go away.

Determinism: atomicAdd across blocks makes V^T W non-reproducible; the host code reduces in a fixed order on purpose. Since this is a tall-skinny GEMM, worth benchmarking cublas<t>gemmBatched before committing to the custom kernel — likely faster and deterministic.

Minor: 1024-thread blocks mean 8 KB shared memory per block for a reduction (128-256 is more typical), and the tree reduction assumes a power-of-two size — worth a static_assert. The device statics are never freed, su2matrix local(n, m) heap-allocates per call, and multiDot/multiDotGPU keep two separate shared statics with a copy between them. The #ifndef SWIG change is unrelated and hides PassiveValue from pySU2.i — please split it out or explain it. Note that .cu files aren't covered by the clang-format hook (types_or: [c++, c]), so the Allman braces / 4-space indent have to be matched by hand.

Testing before this leaves draft: OMP_NUM_THREADS > 1, a CUDA build running ENABLE_CUDA= NO, a mixed-precision build, and 2+ MPI ranks.

[pcarruscag]

Sorry for the previous huge comment from Claude, I had it clean up, the findings look valid.