Skip to content

[lapack][cuSOLVER] Implement native LAPACK gaps - #768

Open
zjin-lcf wants to merge 6 commits into
uxlfoundation:developfrom
zjin-lcf:feature/cusolver-native-lapack-gaps
Open

[lapack][cuSOLVER] Implement native LAPACK gaps#768
zjin-lcf wants to merge 6 commits into
uxlfoundation:developfrom
zjin-lcf:feature/cusolver-native-lapack-gaps

Conversation

@zjin-lcf

@zjin-lcf zjin-lcf commented Aug 16, 2026

Copy link
Copy Markdown
Contributor

Summary

  • implement trtrs (buffer and USM) with cublas<t>trsm, and return 0 from trtrs_scratchpad_size since trsm needs no scratchpad. trsm does not detect the exactly singular matrix that LAPACK reports through info, so the diagonal is searched for a zero before the solve is enqueued, which keeps the LAPACK behaviour of leaving b untouched and reporting the one based index of the first zero diagonal element (as a computation_error)
  • implement the grouped (group_count/group_sizes) USM getri_batch with cublas<t>getriBatched, one native call per group, and grow the grouped getri_batch_scratchpad_size query so the out-of-place inverses of the whole batch fit
  • support nrhs > 1 in potrs_batch (buffer strided, USM strided and USM grouped) by solving the columns of B one at a time with cusolverDn<t>potrsBatched, which only handles a single right-hand side
  • report the info of every matrix of a batch through a new lapack_info_check_batch helper that throws lapack::batch_error, fix the USM get_cusolver_devinfo overload that copied a single int regardless of the batch size and did not wait on the asynchronous copy, and give the batched potrs solves a real (zero-initialised) info array instead of nullptr
  • release the temporary device allocations of the strided and grouped potrs_batch, including on failure paths, instead of leaking them
  • leave gerqf/ormrq/unmrq, hetrf and the m < n cases of gebrd/gesvd unimplemented because cuSOLVER and cuBLAS have no native equivalents

Implementation notes

  • cublas<t>getriBatched requires a uniform n and lda within a call, so one call is issued per group. The oneMath 64-bit pivots are converted into a single 32-bit array holding n contiguous pivots per matrix, as expected by cuBLAS, and only the n x n part of each result is copied back from the scratchpad so that the lda padding of a is preserved.
  • For potrs_batch, all right-hand side column pointer arrays are built and uploaded to the device before the first native call. When SYCL_EXT_ONEAPI_ENQUEUE_NATIVE_COMMAND is available the native calls are not synchronised, so rewriting a device pointer array between two calls would race with a call still reading it.
  • cusolverDn<t>potrsBatched only writes its info array when a parameter is invalid, so the array is zero-initialised before the call.
  • A unit diagonal is not referenced by trsm, so trtrs does not search it for zeros and stays fully asynchronous in that case.

Test plan

  • Build the cuSOLVER LAPACK backend with DPC++ and CUDA 13.3.1 (cuSOLVER 12.2.6) on NVIDIA A100
  • Run runtime dispatch LAPACK tests on NVIDIA GPU: 398 tests, 352 passed, 46 skipped, 0 failed (base commit: 320 passed, 78 skipped)
  • Run compile-time dispatch LAPACK tests on NVIDIA GPU: 398 tests, 352 passed, 46 skipped, 0 failed (base commit: 320 passed, 78 skipped)
  • Verify the 32 tests per dispatch mode that stopped being skipped: Trtrs accuracy and dependency, GetriBatchGroup accuracy and dependency, PotrsBatchGroup accuracy (nrhs = {9, 6}), PotrsBatchStride accuracy for both the buffer and the USM API
  • Check the trtrs singularity path, which the test suite does not cover, with a standalone program: a zero third diagonal element gives computation_error with info() == 3 and leaves b unchanged, for both APIs and for real and complex types, while the same matrix with a unit diagonal solves normally
  • Run clang-format 19.1.7 and git diff --check

Portability

The implementation only uses documented cuSOLVER and cuBLAS entry points together with basic SYCL kernels, and keeps the existing SYCL_EXT_ONEAPI_ENQUEUE_NATIVE_COMMAND code paths intact.

Left unimplemented

Routine Reason
gerqf, ormrq, unmrq neither cuSOLVER nor cuBLAS provides an RQ factorization; left unimplemented for rocSOLVER as well
hetrf oneMath's hetrf is Hermitian, while cuSOLVER only offers the complex symmetric cusolverDn<C,Z>sytrf
gebrd, gesvd with m < n pre-existing documented cuSOLVER limitations

zjin-lcf and others added 4 commits August 16, 2026 07:14
The USM overload of get_cusolver_devinfo copied a single int regardless of
the number of matrices queried and did not wait on the asynchronous copy, so
getrf_batch silently only checked the first matrix of a batch.

Add lapack_info_check_batch, which collects the info value of every matrix
and throws a lapack::batch_error listing all the failing ones, as done by the
rocSOLVER backend.

Co-authored-by: Cursor <cursoragent@cursor.com>
cuSOLVER has no getri, but the strided batch already inverts through
cublas<t>getriBatched. Reuse it for the group API by issuing one call per
group, as cublas requires a uniform n and lda within a call.

The pivots of a group are converted to a single 32-bit array holding n
contiguous pivots per matrix, and the out of place inverses are computed into
the scratchpad before their n x n part is copied back into a. The scratchpad
query now reports room for one output matrix per matrix of the batch instead
of a single matrix.

Co-authored-by: Cursor <cursoragent@cursor.com>
cusolverDnXpotrsBatched only solves for a single right hand side, but solving
the columns of B one at a time is equivalent, so the strided and group
batches no longer report nrhs > 1 as unimplemented. The pointers of all the
columns are uploaded before the first call because the native calls are not
synchronised and would otherwise race with the device array being rewritten.

The batched solves now also pass a real info array, zero initialised because
cuSOLVER only writes it when a parameter is invalid, and the temporary device
allocations of the strided and group batches are released instead of leaked.

Co-authored-by: Cursor <cursoragent@cursor.com>
cuSOLVER has no triangular solve, but cublas<t>trsm solves the same system
with alpha = 1, so trtrs, its USM overload and its scratchpad query no longer
report unimplemented. The query returns 0 because trsm needs no scratchpad.

trsm divides by the diagonal without detecting the exactly singular matrix
that LAPACK reports through info, so the diagonal is searched for a zero
before the solve is enqueued. That keeps the behaviour of LAPACK, which
leaves b untouched and reports the one based index of the first zero diagonal
element, here as a computation_error. A unit diagonal is not referenced by
trsm and is not searched.

Co-authored-by: Cursor <cursoragent@cursor.com>
@zjin-lcf
zjin-lcf requested a review from a team as a code owner August 16, 2026 16:33
zjin-lcf and others added 2 commits August 16, 2026 13:04
The QR based cusolverDn<t>gesvd only factorises m >= n, but the Jacobi
cusolverDn<t>gesvdj has no such restriction, so wide matrices are dispatched
to it instead of reporting unimplemented. gesvdj always computes both vector
sets and returns V rather than V**H, so U and V are staged in the scratchpad
behind the native workspace and copied out transposed and conjugated, which
also covers the jobs that overwrite a. econ is only cleared when vt must hold
all n rows of V**H, and the two jobs that overwrite a are rejected as LAPACK
does. The scratchpad query adds the staging area to the native size.

The zero dimension case returns without referencing s, u or vt, matching
netlib LAPACK, instead of calling cuSOLVER with a zero leading dimension.

The accuracy test now covers every valid job combination for a wide matrix,
reconstructs a from the min(m,n) computed vectors rather than assuming U is
m x n, and copies the vectors out of a before the checks so that the jobs
overwriting a are verified too.

Co-authored-by: Cursor <cursoragent@cursor.com>
The wide job combinations were added as 7 x 11 matrices, which rocSOLVER
factorises to a residual of roughly 80 eps. That is well inside the tolerance
of the LAPACK gesvd test, which accepts a residual of 50 n eps, but
rel_mat_err_check scales by min(m,n) and accepts 10 of those, so 7 x 11 leaves
no headroom in double precision and the reconstruction check fails on that
backend. The bound is equally tight for small tall matrices, so this is not
specific to the wide path.

Use the 24 x 30 shape of the neighbouring case instead, which covers the same
job combinations with the leading dimensions already padded.

Co-authored-by: Cursor <cursoragent@cursor.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant