Harden host-side sparse index width and allocation safety - #2822
Harden host-side sparse index width and allocation safety#2822LwhJesse wants to merge 2 commits into
Conversation
|
After #2816, I think the follow-up direction needs a policy decision before I move #2822 or the resource-cache draft further. I currently have two related follow-ups:
The next steps depend on which CUDA sparse-backend direction the project prefers. I see four practical routes:
My technical preference is route 1, because it gives the cleanest correctness story and the cleanest future GPU-resident Krylov path. If that is too large as the immediate next step, route 2 is a reasonable limited scope for #2822. Could the maintainers confirm which direction they prefer before I mark #2822 ready or move the #2816-based resource-cache draft into the main project? |
|
I think both are going too deep into the software engineering weeds. |
## Proposed Changes This PR adds an end-to-end CUDA linear solve path: the Krylov solvers keep their host control flow, and `CSysVector` operations, the SpMV and the Jacobi preconditioner are dispatched to CUDA kernels when `ENABLE_CUDA=YES`. Transfers are explicit and owned by the object responsible for the data, with no coherency or dirty-flag tracking in `CSysVector` / `CSysMatrix`: - `CSysMatrixVectorProduct` uploads the matrix on construction - the preconditioner uploads its data in `Build()` - `CSysSolve` uploads `b` and `x` and downloads `x` in `HandleTemporariesIn/Out`, which is also where device evaluation is switched on and off - preconditioners without a device implementation (ILU, LU-SGS, Linelet, PaStiX) download the input, apply on the host, and upload the result A 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: - `cuBLAS` for `dot` / `norm`, custom kernels for the block-LDU SpMV, the Jacobi apply, and `CSysVector` expression assignment - `CSysVector` operands 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 in `CSysVectorGPU.cu`, consistent with how `CSysMatrix` is instantiated - device work is issued by a single thread with the OpenMP team synchronized around it, so the GPU path is usable from inside the existing parallel regions; this is internal to the linear algebra layer - the CUDA translation units are only linked into the primal libraries, since they cannot be compiled with the CoDiPack defines; device dispatch is compiled out of the AD builds - adds `LINEAR_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 - CPU vs GPU on inviscid NACA0012, 25 iterations, RTX 4070 Ti SUPER (sm_89), for FGMRES + `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. - Mixed, normal and single precision builds all compile and run; device dispatch confirmed live in each (kernel launch counts track the subspace size). - `OMP_NUM_THREADS=1` and `4` give bit-identical results on the GPU path. - Builds verified: primal, primal + AD + directdiff, `enable-cuda` + `with-omp`, mixed / normal / single precision. - Transfer counts measured per solve: 2 H2D + 1 D2H + 1 matrix upload for Jacobi and `NONE`, plus one download/upload pair per preconditioner application for ILU. Earlier validation of the original design (6 representative cases, `nsys` / `ncu` profiling) predates the rework of the transfer and dispatch model and should be repeated. ## 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>
Harden Host-Side Sparse Index Width and Allocation Safety
Summary
This PR hardens the host-side sparse linear algebra path against integer-width overflow and silent narrowing.
It introduces an explicit 64-bit host sparse-index alias,
and uses it in the host sparse-pattern / CSR metadata path where index width and allocation size matter. It also adds checked multiplication for the main storage-size calculations and checked narrowing at the remaining external boundaries.
This PR does not modify the CUDA backend implementation itself. In particular, it does not change:
Common/src/linear_algebra/CSysMatrixGPU.cuCommon/include/linear_algebra/GPUComms.cuhCUDA matvec correctness remains the responsibility of PR #2816.
Problem
Several sparse linear algebra paths in SU2 currently rely on
unsigned longorintfor sparse metadata, storage-size arithmetic, or boundary conversion.On LP64 systems this often goes unnoticed because
unsigned longis usually 64-bit. On LLP64 systems such as 64-bit Windows,unsigned longis still 32-bit. At sufficiently large local sparse-matrix sizes, that can lead to:nnz * nVar * nEqnFor the main matrix storage, the critical host-side product is:
If that is evaluated in 32-bit arithmetic, the largest representable unsigned value is:
So the first overflowing matrix-scalar count is:
For
double, that means:which is:
32 GiBbinary34.36 GBdecimalThat is the origin of the host-side estimate quoted above.
To connect that number to an actual sparse matrix, take a common square-block case with
nVar = nEqn = 5. Then each nonzero block contributes25matrix scalars, so the first overflowing host-side case is:This is already larger than
2^32 - 1 = 4,294,967,295, so the host-side 32-bit product has crossed the limit. In memory terms, that same case corresponds to:which is again about:
32 GiBbinary34.36 GBdecimalSo the quoted
34.36 GBnumber is not a separate back-of-the-envelope estimate; it is simply the byte size of the first5 x 5double-precision sparse matrix whosennz_blocks * nVar * nEqnproduct no longer fits in 32-bit unsigned arithmetic.Using a rough relation
nnz_blocks ~= N_local * z, withzthe average nonzero-block count per row, this5 x 5threshold corresponds approximately to:17.18 millionlocal points atz = 108.59 millionlocal points atz = 20The same argument gives the corresponding first-overflow cases for larger square blocks:
6 x 6:119,304,648 * 36 = 4,294,967,3287 x 7:87,652,394 * 49 = 4,294,967,306So the threshold is high, but it is still a local sparse-matrix scale that can be reached in large implicit runs.
For reference, the current legacy CUDA kernel in
developcan overflow earlier because it uses a signed 32-bitintmatrix-scalar offset internally. The largest positive signed 32-bit value is:so the first overflowing matrix-scalar count is:
For
double, that means:which is:
16 GiBbinary17.18 GBdecimalThat estimate is specific to the legacy custom CUDA matvec path currently used in
develop.PR #2816 naturally addresses that CUDA-side issue by replacing the old custom matvec kernel, including the legacy signed-
intmatrix-offset path that creates this earlier overflow limit. That is why this PR does not modify the CUDA backend itself.What this PR changes
su2_index_t = std::uint64_tfor the host sparse-index chain that needs it.CSysMatrixsparse metadata,CSysVectorsize/index API,CPastixWrappersparse input handling, and the sparse-pattern-facing geometry accessors required for that propagation.nnz * nVar * nEqnnnz_ilu * nVar * nEqnnPointDomain * nVar * nEqnnumBlk * numVarnumBlkDomain * numVarunsigned long, but range-checks host sparse indices before converting them for the current CUDA upload path.After this change, the targeted host-side sparse-index and allocation path no longer fails by silent wraparound or silent truncation. If a value exceeds the remaining boundary types, the code now fails explicitly through
SU2_MPI::Error(...).Validation
git diff --checkpassed.For numerical validation, I ran the existing six-case correctness harness with a two-way CPU comparison:
develop CPUthis branch CPUCases:
periodic2d_sectorudf_lam_flatplate_sudf_lam_flatplate_mudf_lam_flatplate_ludf_test_11_probes_sudf_test_11_probes_mResult:
this branch CPUmatcheddevelop CPUin all 6 caseshistory.csvmatched exactly in all 6 casesmax_abs_delta = 0.0for every caseFor CUDA, the relevant correctness discussion remains PR #2816, since this branch still uses the legacy pre-#2816 custom CUDA matvec implementation and this PR intentionally does not modify that backend. In particular, the current
develop-side signed-intoffset limit described above is handled naturally by the#2816backend replacement rather than by this PR.Notes on testing
This PR does not add a direct large-memory reproducer for the original host-side overflow. The original failure mode depends on LLP64-style 32-bit
unsigned longarithmetic, and a faithful end-to-end reproduction would require a much larger local sparse structure than is practical for a routine test here.PR Checklist
--warnlevel=3when using meson).pre-commit run --allto format old commits.