T17: report

TASKS.md

@@ -2272,7 +2272,7 @@ repeat it once T21 has landed. The bar is cuDSS with `"matching_alg" =
 
 ---
 
-## T17 — Uniform batch (v1 cut line)   `[ ]`
+## T17 — Uniform batch (v1 cut line)   `[!]`
 
 **Reads**: PLAN §1.6, §3.5, §1.3 (`ubatch_*`), §1.4 (`ubatch_mask`);
 `../CUDSS.jl/docs/src/uniform_batch.md`, `test/test_uniform_batch_cudss.jl`.
@@ -2303,6 +2303,79 @@ repeat it once T21 has landed. The bar is cuDSS with `"matching_alg" =
 **This completes v1 (PLAN M0–M6).** The owner re-evaluates the plan before the
 tasks below are refined.
 
+### Report
+
+- Status: [!] (done on the KA CPU backend; CUDA from CI; deviations below)
+- What was built:
+  - `src/numeric/batch.jl` (new): `BatchMap` (active members, `nact`, `nbatch`, the launch's `first`, `nrhs` per
+    member; Adapt-able, `members = nothing` for single-matrix helper launches), member decoding of the flattened
+    group index (`_bm_node`/`_bm_slot`/`_bm_cmember`/`_bm_ucol`), member views (`_mview` member-major,
+    `_iview` interleaved), `MemberPanels` (the panel pointers of member `k`, passed to the kernels in place of
+    `front_ptr`), `panel_offset`, `active_members` (`ubatch_index` 0-based ∩ `ubatch_mask`), `member_runs!`.
+  - Layout (PLAN §3.5): factor panels interleaved per front (the `nb` panels of supernode `s` contiguous, so a
+    regime-C front of a run of consecutive members is one `f×w×count` strided batch); `info` interleaved per
+    entry (batched `potrf` statuses contiguous); `d`, `piv`, `pivot_kind`, update stack, regime-C workspace,
+    `stats`, `totals`, `aux`, `nzval` and the right-hand sides member after member. `nb = 1` is exactly the old
+    layout and the old code paths (regime C, solve dense path).
+  - `Numeric` gains `members` (device `Int32`) and `nbatch`; `NumericPlan` gains `members_host`, `nact`, `runs`,
+    `vendor_ptrs`. `allocate_numeric(S, T, backend; nbatch)`, `set_members!`, `batch_map`, `member_info!`,
+    `member_numeric` (host single-member copy), `panel(S, N, s, k)`, `pivot_totals(N, k)`.
+  - Every numeric kernel (zero/scatter/extend-add, `pack_add`, fused regime-B Cholesky, regime-A Cholesky, both
+    LDLᵀ kernels, `cholesky_stats`, `reduce_stats`, `abs_max`) and every solve kernel (permute/unpermute,
+    forward/backward sweeps, pivot order, diagonal, residual, residual norms, add-correction, `"diag"`) takes the
+    batch index as one more grid dimension, flattened into the 1-D ndrange (as the right-hand sides already were;
+    2-D ndranges allocate on the KA CPU backend). The `BatchMap` replaces the `first` argument, so no kernel goes
+    over the 32-argument specialization limit; the launches stay allocation-free (CPU budgets unchanged).
+  - Regime C with `nb > 1`: per run of consecutive active members, strided-batched `potrf` (`potrf_batched_info!`,
+    new in the dense interface: vendor or KA `ka_potrf!`), `trsm_strided_batched!`, `gemm_strided_batched!`
+    (`F21 F21ᴴ` into `count` workspace blocks) and a batched `pack_add!` (forces the herk real diagonal);
+    `impl = :generic` falls back to the single-matrix calls per member. Vendor `potrf`/`trsm` run on member pointer
+    arrays built once at allocation (`vendor_batch_pointers`, `vendor_potrf_batched_ptrs!`,
+    `vendor_trsm_batched_ptrs!`, new vendor functions; CUDA extension implements them, the existing allocating
+    `vendor_*_batched!` now delegate to them). Solve regime C: strided-batched `gemm`; `trsm` batched on the KA
+    path, member by member on the vendor path (see deviations).
+  - `DirectSolver`: `nbatch` deduced from `length(nzval) ÷ nnz` or an `nnz × nbatch` matrix; `ubatch_size` must be
+    0 or equal (`InvalidValueError`); `info` is a per-member vector (`getparam` returns `Int` for `nb = 1`, a vector
+    otherwise; `setparam!` accepts either); `npivots`, `inertia`, `pivot_stats` per member; `"diag"` = members'
+    diagonals one after the other; `ubatch_index`/`ubatch_mask` restrict factorization and solve phases (inactive
+    members' factors, `info` and solution columns untouched). Right-hand sides: strided vectors, `n × (nrhs nb)`
+    matrices, `n × nrhs × nb` arrays, `MatrixDescriptor(T, n[, nrhs]; nbatch)` (`update!` now also accepts the
+    `n × (ncols nbatch)` matrix of a batched descriptor, as `CudssMatrix(T, n; nbatch)` + `CuMatrix` in CUDSS.jl).
+    Refinement per member (residual with each member's values). LinearAlgebra layer: auto-detection through
+    `DirectSolver`, `ldiv!`/`\` on 3-D arrays, `check` reports the info vector; `logabsdet` of a batch raises
+    `NotSupportedError`.
+  - Tests: `test/test_ubatch.jl` (new), `test/ported/cudss_uniform_batch.jl` (new: LDLᴴ and Cholesky parts of
+    CUDSS.jl's `test_uniform_batch_cudss.jl`, both APIs, strided/non-strided, views L/U/F; LU skipped, T19);
+    shared helpers `batch_members` (test/matrices.jl), `batch_csr`, `batch_relres` (test/utils.jl), `api_csr`,
+    `api_batch_matrix` (test/backends.jl).
+- Tests: `SDS_TEST_GPU=0 julia --project=. -e 'using Pkg; Pkg.test()'` (CPU): 62074 pass, 0 fail, 1 broken (the
+  T16 `@test_broken`), ≈ 14 min; the counts combine the full run (62050 pass, 24 failures of one wrong test
+  assertion in `test_ubatch`, member 1 of a fresh `batch_members` draw is `A` itself) and the rerun of
+  `test_ubatch` after fixing that assertion (1092 pass) — `test_ported`, `test_solve`, `test_refinement` were
+  rerun on the final code as well (9655 pass, 1 broken). `test_ubatch` alone: 1092 pass. CUDA/AMDGPU: pending
+  CI on the PR (the CUDA pointer-array bindings are untested locally).
+- Measurements: KA CPU backend, 1 thread, `laplacian2d(40, 40)` (`n = 1600`), 64 members: batched
+  refactorization 209 ms / solve 36 ms vs 64 single solvers 212 ms / 39 ms — the same arithmetic, as expected on
+  the CPU; the batch's gain is the launch count on GPUs (one launch per group for all members instead of 64).
+- Deviations from PLAN.md / this task:
+  - Solve phase, regime C: cuBLAS has no strided `trsm`, and pointer arrays of `Y` depend on `nrhs` (would be a
+    device allocation per solve), so the vendor path solves `trsm` member by member (the `gemm` stays strided
+    batched; the KA path is batched). The numeric phase uses the prebuilt pointer arrays as asked.
+  - Row-major (`transposed`) right-hand sides of a batch raise `NotSupportedError`; `logabsdet` of a batch too.
+  - The 2×2 pivot pairs of `"S"`/`"H"` (T14 analysis, value based) are computed from the first member's values.
+  - `test/test_api.jl`: three assertions that encoded "not implemented until T17" (`DirectSolver` on batched
+    values, a 3-D right-hand side, `ubatch_size = 2`) now assert the implemented behaviour (a solved batch,
+    `DimensionMismatch` for 4 columns against 1, `InvalidValueError` for a size mismatch).
+- Open issues / follow-ups:
+  - The CUDA bindings `vendor_batch_pointers`/`vendor_*_batched_ptrs!` and the batched regime C on CUDA are
+    exercised only by CI. T23 (AMDGPU/oneAPI/Metal) needs the same three bindings, or falls back to the
+    allocating batched calls / KA (handled: a missing binding leaves `vendor_ptrs` empty).
+  - The GitHub App token expired before the last push: the final commits are local and published by the
+    workflow's "Publish the branch" step.
+- Suggested plan changes:
+  - PLAN §3.5: record the layout choice (panels interleaved per front, `info` interleaved, the rest member-major)
+    and the solve-phase vendor `trsm` per member.
+
 ---
 
 ## Post-v1 tasks (to be refined after the v1 evaluation)

ext/SparseDirectSolverCUDAExt.jl

@@ -3,7 +3,8 @@ module SparseDirectSolverCUDAExt
 # CUDA support: CSR adapters (T02) and the vendor dense bindings of the dense
 # layer (T03, PLAN §2.6): cuBLAS gemm/syrk/herk/trsm/strided-batched gemm/
 # batched trsm, cuSOLVER potrf (also with a device info)/getrf/sytrf/potrfBatched,
-# cuBLAS getrfBatched.
+# cuBLAS getrfBatched; the batched trsm/potrf also on member pointers built once
+# (uniform batch, T17).
 
 using SparseDirectSolver
 using SparseDirectSolver: CSR, DirectSolver, INDEX_ONE, INDEX_ZERO, InvalidValueError
@@ -225,24 +226,36 @@ function _batch_pointers(A::StridedCuArray{T, 3}) where {T}
     return CuArray([base + (i - 1) * s for i in 1:size(A, 3)])
 end
 
+# built once per regime-C front of a uniform batch (T17): the numeric phase then allocates nothing
+SDS.vendor_batch_pointers(A::StridedCuArray{T, 3}) where {T <: CuBlasT} = _batch_pointers(A)
+
 for (fname, fname_64, elty) in ((:cublasStrsmBatched, :cublasStrsmBatched_64, :Float32),
                                 (:cublasDtrsmBatched, :cublasDtrsmBatched_64, :Float64),
                                 (:cublasCtrsmBatched, :cublasCtrsmBatched_64, :ComplexF32),
                                 (:cublasZtrsmBatched, :cublasZtrsmBatched_64, :ComplexF64))
+    @eval function SDS.vendor_trsm_batched_ptrs!(side::Char, uplo::Char, trans::Char, diag::Char, α, m::Integer,
+                                                 n::Integer, Aptrs::StridedCuArray{CuPtr{$elty}, 1}, lda::Integer,
+                                                 Bptrs::StridedCuArray{CuPtr{$elty}, 1}, ldb::Integer, count::Integer)
+        count == 0 && return Bptrs
+        a = _blas_scalar($elty, α)
+        if cuBLAS.version() >= v"12.0"
+            cuBLAS.$fname_64(cuBLAS.handle(), side, uplo, trans, diag, m, n, a, Aptrs, max(1, lda), Bptrs, max(1, ldb),
+                             count)
+        else
+            cuBLAS.$fname(cuBLAS.handle(), side, uplo, trans, diag, m, n, a, Aptrs, max(1, lda), Bptrs, max(1, ldb),
+                          count)
+        end
+        return Bptrs
+    end
     @eval function SDS.vendor_trsm_batched!(side::Char, uplo::Char, trans::Char, diag::Char, α,
                                             A::StridedCuArray{$elty, 3}, B::StridedCuArray{$elty, 3})
         m, n, nb = size(B)
         nb == 0 && return B
-        lda, ldb = max(1, stride(A, 2)), max(1, stride(B, 2))
         GC.@preserve A B begin
             Aptrs = _batch_pointers(A)
             Bptrs = _batch_pointers(B)
-            a = _blas_scalar($elty, α)
-            if cuBLAS.version() >= v"12.0"
-                cuBLAS.$fname_64(cuBLAS.handle(), side, uplo, trans, diag, m, n, a, Aptrs, lda, Bptrs, ldb, nb)
-            else
-                cuBLAS.$fname(cuBLAS.handle(), side, uplo, trans, diag, m, n, a, Aptrs, lda, Bptrs, ldb, nb)
-            end
+            SDS.vendor_trsm_batched_ptrs!(side, uplo, trans, diag, α, m, n, Aptrs, stride(A, 2), Bptrs, stride(B, 2),
+                                          nb)
             CUDACore.unsafe_free!(Aptrs)
             CUDACore.unsafe_free!(Bptrs)
         end
@@ -252,13 +265,19 @@ end
 
 for (fname, elty) in ((:cusolverDnSpotrfBatched, :Float32), (:cusolverDnDpotrfBatched, :Float64),
                       (:cusolverDnCpotrfBatched, :ComplexF32), (:cusolverDnZpotrfBatched, :ComplexF64))
+    @eval function SDS.vendor_potrf_batched_ptrs!(uplo::Char, n::Integer, Aptrs::StridedCuArray{CuPtr{$elty}, 1},
+                                                  lda::Integer, info::StridedCuArray{Cint, 1}, count::Integer)
+        count == 0 && return info
+        cuSOLVER.$fname(cuSOLVER.dense_handle(), uplo, n, Aptrs, max(1, lda), info, count)
+        return info
+    end
     @eval function SDS.vendor_potrf_batched!(uplo::Char, A::StridedCuArray{$elty, 3}, info::CuVector{Cint})
         n = _check_square_batch(A)
         nb = size(A, 3)
         nb == 0 && return info
         GC.@preserve A begin
             Aptrs = _batch_pointers(A)
-            cuSOLVER.$fname(cuSOLVER.dense_handle(), uplo, n, Aptrs, max(1, stride(A, 2)), info, nb)
+            SDS.vendor_potrf_batched_ptrs!(uplo, n, Aptrs, stride(A, 2), info, nb)
             CUDACore.unsafe_free!(Aptrs)
         end
         return info

src/SparseDirectSolver.jl

@@ -52,6 +52,7 @@ include("symbolic/layout.jl")
 include("symbolic/maps.jl")
 
 # numeric storage, CPU reference multifrontal factorization (PLAN §7, the oracle)
+include("numeric/batch.jl")
 include("numeric/storage.jl")
 include("reference/cholesky.jl")
 include("reference/ldlt.jl")

src/dense/fallback/potrf.jl

@@ -115,17 +115,18 @@ end
 
 @kernel function _ka_chol_check_kernel!(info, A, n, idx)
     li = @index(Local, Linear)
+    b = @index(Group, Linear)
     if li == 1
-        @inbounds if info[idx] == 0
+        @inbounds if info[idx + b - 1] == 0
             st = Int32(0)
             for j in 1:n
-                d = _get(A, j, j, 1)
+                d = _get(A, j, j, b)
                 if !(isfinite(real(d)) && real(d) > 0 && iszero(imag(d)))
                     st = Int32(j)
                     break
                 end
             end
-            info[idx] = st
+            info[idx + b - 1] = st
         end
     end
 end
@@ -145,3 +146,19 @@ function ka_chol_check_info!(info::AbstractVector{Int32}, idx::Integer, L::Abstr
     kernel!(info, L, size(L, 1), Int(idx); ndrange = 1)
     return info
 end
+
+"""
+    ka_chol_check_info!(info, offset, L::AbstractArray{T,3}) -> info
+
+Batched form: the check of `ka_chol_check_info!(info, offset + b, L[:, :, b])`
+for every member `b` of the 3-D strided batch `L`, one launch.
+"""
+function ka_chol_check_info!(info::AbstractVector{Int32}, offset::Integer, L::AbstractArray{<:Any, 3})
+    nb = size(L, 3)
+    0 <= offset && offset + nb <= length(info) ||
+        throw(DimensionMismatch("info has length $(length(info)) < offset $offset + batch count $nb"))
+    nb == 0 && return info
+    kernel! = _ka_chol_check_kernel!(KernelAbstractions.get_backend(L), 1)
+    kernel!(info, L, size(L, 1), Int(offset) + 1; ndrange = nb)
+    return info
+end

src/dense/interface.jl

@@ -25,6 +25,7 @@ const DENSE_OPS = (
     laswp = (generic = nothing, vendor = nothing),
     gemm_strided_batched = (generic = nothing, vendor = :vendor_gemm_strided_batched),
     trsm_strided_batched = (generic = nothing, vendor = :vendor_trsm_batched),
+    potrf_batched = (generic = nothing, vendor = :vendor_potrf_batched),
 )
 
 # DENSE_OPS as a concretely typed table (`:none` = no such path), for `select_impl`
@@ -74,7 +75,8 @@ end
     dense_impls(op::Symbol, backend, T) -> Vector{Symbol}
 
 The implementations of dense op `op` (`:gemm`, `:syrk`, `:herk`, `:trsm`,
-`:potrf`, `:getrf`, `:laswp`, `:gemm_strided_batched`, `:trsm_strided_batched`)
+`:potrf`, `:getrf`, `:laswp`, `:gemm_strided_batched`, `:trsm_strided_batched`,
+`:potrf_batched`)
 available for element type `T` on `backend` according to
 [`capabilities`](@ref), in `:auto` preference order (`:vendor`, `:generic`,
 `:ka`). `:ka` is always present.
@@ -382,6 +384,38 @@ function trsm_strided_batched!(side, uplo, trans, diag, α, A::AbstractArray{<:A
     return B
 end
 
+"""
+    potrf_batched_info!(uplo, A, info, offset = 0; impl = :auto) -> info
+
+Uniform strided batch of [`potrf_info!`](@ref): Cholesky of the `uplo`
+triangle of every member `A[:, :, b]` of the 3-D array `A` in place, its
+status (0, or the first non-positive pivot column) into `info[offset + b]`.
+Implementations: `:vendor` (`vendor_potrf_batched!`, validated on the device
+as `potrf_info!`) and `:ka` ([`ka_potrf!`](@ref), one workgroup per member).
+Asynchronous.
+"""
+function potrf_batched_info!(uplo, A::AbstractArray{<:Any, 3}, info::AbstractVector{Int32}, offset::Integer = 0;
+                             impl::Symbol = :auto)
+    ul = _uplo_char(uplo)
+    size(A, 1) == size(A, 2) || throw(DimensionMismatch("potrf: members are $(size(A, 1))×$(size(A, 2)), not square"))
+    0 <= offset && offset + size(A, 3) <= length(info) ||
+        throw(DimensionMismatch("info has length $(length(info)) < offset $offset + batch count $(size(A, 3))"))
+    size(A, 3) == 0 && return info
+    return _potrf_batched_impl!(select_impl(:potrf_batched, A, impl), A, info, offset, ul)
+end
+
+# `potrf_batched_info!` with an already resolved implementation `p` (no capability lookup)
+function _potrf_batched_impl!(p::Symbol, A::AbstractArray{<:Any, 3}, info::AbstractVector{Int32}, offset::Integer,
+                              ul::Char = 'L')
+    if p === :vendor
+        vendor_potrf_batched!(ul, A, view(info, (offset + 1):(offset + size(A, 3))))
+        ka_chol_check_info!(info, offset, A)
+    else
+        ka_potrf!(ul, A, info; offset)
+    end
+    return info
+end
+
 """
     strided_batch(buf, offset, m, n, stride, count) -> AbstractArray{T,3}
 

src/dense/vendor.jl

@@ -9,7 +9,8 @@
 
 const VENDOR_FUNCTIONS = (:vendor_gemm!, :vendor_syrk!, :vendor_herk!, :vendor_trsm!, :vendor_potrf!,
                           :vendor_potrf_info!, :vendor_getrf!, :vendor_sytrf!, :vendor_gemm_strided_batched!,
-                          :vendor_trsm_batched!, :vendor_potrf_batched!, :vendor_getrf_batched!)
+                          :vendor_trsm_batched!, :vendor_potrf_batched!, :vendor_getrf_batched!,
+                          :vendor_batch_pointers, :vendor_trsm_batched_ptrs!, :vendor_potrf_batched_ptrs!)
 
 for f in VENDOR_FUNCTIONS
     @eval $f(args...) = throw(NotSupportedError(string($(string(f)), " has no vendor binding for argument types ",
@@ -29,6 +30,13 @@ end
     vendor_trsm_batched!(side, uplo, trans, diag, α, A, B) -> B          # 3-D arrays
     vendor_potrf_batched!(uplo, A, info) -> info                         # 3-D A, Int32 info
     vendor_getrf_batched!(A, ipiv, info) -> info                         # 3-D A, Int32 ipiv matrix
+    vendor_batch_pointers(A) -> ptrs                                     # device vector of the member pointers of 3-D A
+    vendor_trsm_batched_ptrs!(side, uplo, trans, diag, α, m, n, Aptrs, lda, Bptrs, ldb, count)
+    vendor_potrf_batched_ptrs!(uplo, n, Aptrs, lda, info, count) -> info # Int32 info, count entries
+
+The `_ptrs` variants take device vectors of member pointers built once by
+`vendor_batch_pointers` (the uniform-batch regime C, T17: no allocation in the
+numeric phase); `vendor_trsm_batched!`/`vendor_potrf_batched!` build them per call.
 
 Vendor dense routines behind the `impl = :vendor` path of the dense interface.
 The generic methods throw `NotSupportedError`; backends add methods for their