Fix: view-free G, zero-pivot info, analysis-time pair semantics, overflow guards, docstrings (#84)

TASKS.md

@@ -2719,6 +2719,10 @@ the Report whether the a posteriori pivoting of M13 (T27) must move before
 the MadNLP integration; the cuDSS evidence says matching is enough. Closes
 #71.
 
+**Owner note (issue #86)**: compare the device and the reference on the K2
+dumps after scaling (inertia, `nperturbed`, `piv`); close #86 if scaling
+removed the divergence, else close it as documented.
+
 Host MC64-style matching (job 5 first), `perm_matching`, `scale_row/col`,
 composition with the ordering; test: on a badly scaled unsymmetric matrix the
 LU with matching has `nperturbed == 0` where the LU without has `> 0`;

src/matrix.jl

@@ -96,10 +96,18 @@ end
 Zero-copy reinterpretation of the CSC matrix `A` (`m × n`) as the CSR matrix of
 `transpose(A)` (`n × m`): `A.colptr` becomes `rowptr`, `A.rowval` becomes
 `colval`, the arrays are shared and `transposed = true` records that the caller
-means `A`. This is how MadNLP passes `colPtr`/`rowVal` (PLAN §1.1).
+means `A`. This is how MadNLP passes `colPtr`/`rowVal` (PLAN §1.1). A solver
+on such a matrix reads it with the view flipped (`'L'` ↔ `'U'`; for `"G"` the
+view is ignored). Raises [`InvalidValueError`](@ref) when `A.rowval` or
+`A.nzval` is longer than `nnz(A)` (a `SparseMatrixCSC` may keep spare capacity,
+which the shared CSR arrays cannot describe).
 """
 function csr_of_transpose(A::SparseMatrixCSC)
     m, n = size(A)
+    nz = Int(A.colptr[end]) - 1
+    (length(A.rowval) == nz && length(A.nzval) == nz) ||
+        throw(InvalidValueError("csr_of_transpose: rowval/nzval have $(length(A.rowval))/$(length(A.nzval)) entries, " *
+                                "nnz(A) = $nz; `resize!` the buffers to nnz(A) first"))
     return CSR(A.colptr, A.rowval, A.nzval, n, m; index = INDEX_ONE, transposed = true)
 end
 

src/numeric/assembly.jl

@@ -24,7 +24,7 @@ const ASSEMBLY_WORKGROUP = 256
         end
         c = cb_ptr[s]
         if c > 0
-            m = front_nrows[s] - front_ncols[s]
+            m = Int(front_nrows[s]) - Int(front_ncols[s])  # packed sizes in Int: m(m+1)/2 overflows Int32 first
             for q in (c + li - 1):WG:(c + m * (m + 1) ÷ 2 - 1)
                 stack[q] = z
             end
@@ -133,9 +133,9 @@ end
             c = child_list[kc]
             cb = cb_ptr[c]
             if cb > 0
-                mc = front_nrows[c] - front_ncols[c]
-                fp = front_nrows[s]
-                wp = front_ncols[s]
+                mc = Int(front_nrows[c]) - Int(front_ncols[c])   # packed offsets in Int (`_packed` of a block that
+                fp = Int(front_nrows[s])                         # fits INT can overflow INT for m > 46340)
+                wp = Int(front_ncols[s])
                 mp = fp - wp
                 pp = front_ptr[s]
                 cp = cb_ptr[s]

src/numeric/factorize.jl

@@ -199,7 +199,10 @@ _check_ldlt_eltype(::Type{T}) where {T} = T <: LinearAlgebra.BlasFloat ||
     factorize!(numeric, symbolic, nzval; impl = :auto, opts = Options()) -> info::Int
 
 Structures `"S"`/`"H"`: [`factorize_ldlt!`](@ref) with the pivoting options
-`opts` (`impl` is not used); structure `"G"`: [`factorize_lu!`](@ref) likewise. Structures `"SPD"`/`"HPD"`: multifrontal Cholesky `P A Pᵀ = L Lᴴ` on the device: `nzval` are the stored
+`opts` (`impl` selects its regime-C GEMMs); structure `"G"`:
+[`factorize_lu!`](@ref) (`impl` is not used). Both return `info` (nonzero only
+for a zero pivot with `pivot_epsilon = 0`). Structures `"SPD"`/`"HPD"`:
+multifrontal Cholesky `P A Pᵀ = L Lᴴ` on the device: `nzval` are the stored
 values of A (same pattern, view and index base as the analysis) on the backend
 of `numeric`, and `symbolic` has its maps on that backend
 (`adapt(backend, symbolic, INT)`). For every launch group of the schedule,
@@ -274,8 +277,21 @@ function factorize_cholesky!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector;
     end
     assemble_schur!(N, S, nzval)
     cholesky_stats!(N, S)
+    return _numeric_info!(N, S, true)
+end
+
+# the reduced `info` of the active members into `plan.info_host` (`read`: copied from the device, the phase's only
+# host synchronization; else 0, a factorization that cannot fail) and the first nonzero one as an original column
+function _numeric_info!(N::Numeric, S::Symbolic, read::Bool)
+    plan = N.plan
     nb = N.nbatch
-    copyto!(plan.info_host, 1, N.info, nsupernodes(S) * nb + 1, nb)  # the phase's only host synchronization
+    if read
+        copyto!(plan.info_host, 1, N.info, nsupernodes(S) * nb + 1, nb)
+    else
+        for j in 1:plan.nact[]
+            plan.info_host[plan.members_host[j]] = 0
+        end
+    end
     nb == 1 && return _original_column(S, plan.info_host[1])
     for j in 1:plan.nact[]
         k = Int(plan.members_host[j])
@@ -291,16 +307,16 @@ _original_column(S::Symbolic, k::Integer) = k == 0 ? 0 : S.partition.perm[k]
     member_info!(info, numeric, symbolic) -> info
 
 Write the `info` of the last [`factorize!`](@ref) of every active batch member
-`k` into `info[k]` (the original column of the member's first non-positive
-pivot, 0 = success; inactive members keep their entries). Host data only: the
-statuses were read by `factorize!` (Cholesky) or are 0 (LDLᵀ/LDLᴴ, LU).
+`k` into `info[k]` (the original column of the member's first failed pivot:
+non-positive for Cholesky, zero with `pivot_epsilon = 0` for LDLᵀ/LDLᴴ and LU;
+0 = success; inactive members keep their entries). Host data only: the
+statuses were read by `factorize!`.
 """
 function member_info!(info::AbstractVector{<:Integer}, N::Numeric, S::Symbolic)
     plan = N.plan
-    chol = !_is_ldlt_structure(S.structure) && S.structure != STRUCTURE_GENERAL
     for j in 1:plan.nact[]
         k = Int(plan.members_host[j])
-        info[k] = chol ? _original_column(S, plan.info_host[k]) : 0
+        info[k] = _original_column(S, plan.info_host[k])
     end
     return info
 end

src/numeric/ldlt.jl

@@ -8,11 +8,16 @@
 # * Regime A: `subtree_ldlt_kernel!`, one workgroup per subtree, the serial
 #   stack of packed fronts and contribution blocks in `@localmem` (as the
 #   Cholesky kernel of `src/numeric/subtree.jl`).
-# * Regimes B and C: `front_ldlt_kernel!`, one workgroup per front of a launch
+# * Regime B: `front_ldlt_kernel!`, one workgroup per front of a launch
 #   group, fused assembly (zero, scatter A, owner-pull extend-add), then the
-#   pivoted factorization of the panel (regime B: `F₁₁` staged in `@localmem`,
-#   the rows below it in global memory; regime C: all in global memory) and the
-#   update of the front's packed contribution block on the update stack.
+#   pivoted factorization of the panel and the update of the front's packed
+#   contribution block on the update stack. The kernel's `Val(W)`: `W > 0` (width
+#   classes above `_LT_GLOBAL_MAX_W`) stages `F₁₁` as a `W×W` triangle in
+#   `@localmem`, the rows below it in global memory; `W = 0` (the narrow classes
+#   `≤ _LT_GLOBAL_MAX_W`) keeps the whole panel in global memory.
+# * Regime C, and the wide regime-B bins of `ldlt_blocked_path`: the blocked
+#   path of `src/numeric/ldlt_c.jl` (panel kernel plus GEMMs). In the launch plan,
+#   `group_width[k] = 0` marks a group on that path (not the kernel's `W = 0`).
 #
 # Per pivot step (a loop over the `w` columns of the front with a uniform trip
 # count; a 2×2 pivot leaves the last iterations idle): the workgroup chooses the
@@ -392,11 +397,13 @@ end
             x = _fget(F, c, c)
             dk = H ? T(real(x)) : x
             kind = PIVOT_KIND_1X1
-            if !(abs(dk) >= ε)                        # tiny (or NaN): perturb
+            if !(abs(dk) >= ε) || iszero(dk)          # tiny (or NaN, or exactly zero): perturb
                 iszero(dk) && (ctl[_LT_STAT + STAT_NZERO] += one(IT))
                 dk = _perturbation_sign(dk, Int(psign[perm[piv[c0 + c - 1]]]), H) * ε
                 kind = PIVOT_KIND_PERTURBED
                 ctl[_LT_STAT + STAT_NPERTURBED] += one(IT)
+                # ε = 0 leaves the pivot zero: the front fails at this assembled local column
+                iszero(dk) && ctl[_ST_STATUS] == 0 && (ctl[_ST_STATUS] = (Int(piv[c0 + c - 1]) - c0 + 1) % IT)
             end
             d[g] = dk
             d[n + g] = zero(T)
@@ -573,6 +580,7 @@ end
         for q in 1:5
             ctl[_LT_STAT + q] = zero(IT)
         end
+        ctl[_ST_STATUS] = zero(IT)
     end
     return nothing
 end
@@ -925,8 +933,8 @@ end
         for q in 1:5
             stats[base + q] = Int64(ctl[_LT_STAT + q])
         end
-        stats[base + STAT_INFO] = 0
-        info[s] = Int32(0)
+        stats[base + STAT_INFO] = Int64(ctl[_ST_STATUS])
+        info[s] = Int32(ctl[_ST_STATUS])
     end
     return nothing
 end
@@ -965,14 +973,15 @@ end
                                     front_ncols, cb_ptr, child_ptr, child_list, relind_ptr, relind, maxchild, prm,
                                     Val(W), Val(WG); ndrange = WG * count * bm.nact)
 
-LDLᵀ/LDLᴴ kernel of regimes B and C: workgroup `G` takes front
+LDLᵀ/LDLᴴ kernel of regime B: workgroup `G` takes front
 `s = nodes[bm.first + g - 1]` of batch member `k` (`g`, `k` from the
 [`BatchMap`](@ref) `bm`; the per-member arrays are the member's), zeroes and assembles it (A through the `amap`, then
 the `maxchild` children's packed contribution blocks in `child_list` order),
 factors its `w` fully-summed columns in place with in-block
-Bunch–Kaufman pivoting (`W = 0`, regime C: in the panel; regime B, width
-class `W ≥ w`: `F₁₁` staged as a packed `W×W` lower triangle in `@localmem`,
-the rows below it in the panel; the pivot search is a workgroup reduction), threshold acceptance and perturbation (`prm`, a
+Bunch–Kaufman pivoting (`W = 0`, the narrow width classes `≤ _LT_GLOBAL_MAX_W`:
+in the panel in global memory; width class `W ≥ w` above that: `F₁₁` staged
+as a packed `W×W` lower triangle in `@localmem`, the rows below it in the panel;
+the pivot search is a workgroup reduction), threshold acceptance and perturbation (`prm`, a
 `_LDLTDevice{R, HERM}`; `HERM`: Hermitian or real symmetric, else complex symmetric),
 writes D, the local pivot order `piv`, the pivot kinds and the front's
 statistics, updates its packed contribution block on the update stack
@@ -1296,11 +1305,12 @@ function abs_max!(aux::AbstractVector, nzval::AbstractVector, bm::BatchMap = sin
     return aux
 end
 
-@kernel function _reduce_stats_kernel!(totals, stats, ns, bm, ::Val{WG}, ::Val{NF},
+@kernel function _reduce_stats_kernel!(totals, info_all, stats, super_ptr, ns, bm, ::Val{WG}, ::Val{NF},
                                        ::Val{LOG2WG}) where {WG, NF, LOG2WG}
     li = @index(Local, Linear)
     G = @index(Group, Linear)
     acc = @localmem Int64 (WG * NF,)
+    best = @localmem Int64 (WG,)
     @inbounds begin
         st = _mview(stats, _bm_gmember(bm, G), bm.nbatch)
         for q in 1:NF
@@ -1311,6 +1321,12 @@ end
             end
             acc[(q - 1) * WG + li] = a
         end
+        m = typemax(Int64)
+        for s in li:WG:ns
+            fi = st[s * NF]
+            fi > 0 && (m = min(m, Int64(super_ptr[s]) + fi - 1))
+        end
+        best[li] = m
     end
     @synchronize
     for lev in 1:LOG2WG
@@ -1320,29 +1336,37 @@ end
                 for q in 1:NF
                     acc[(q - 1) * WG + li] += acc[(q - 1) * WG + li + h]
                 end
+                best[li] = min(best[li], best[li + h])
             end
         end
         @synchronize
     end
     if li <= NF
         @inbounds _mview(totals, _bm_gmember(bm, G), bm.nbatch)[li] = acc[(li - 1) * WG + 1]
     end
+    if li == 1
+        @inbounds _iview(info_all, _bm_gmember(bm, G), bm.nbatch)[ns + 1] =
+            best[1] == typemax(Int64) ? Int32(0) : Int32(best[1])
+    end
 end
 
 """
     reduce_stats!(numeric, symbolic) -> numeric
 
 Sum the per-front statistics `numeric.stats` into `numeric.totals`
 (`npos, nneg, nzero, nperturbed, n2x2`, and the number of fronts with a failed
-pivot) of every active batch member: one workgroup per member, `@localmem`
+pivot) of every active batch member, and the smallest failed factor column
+(`super_ptr[s] + info - 1` over the fronts with a nonzero `info` statistic,
+0 = none) into `numeric.info[ns + 1]`: one workgroup per member, `@localmem`
 tree reduction, no atomics. Asynchronous; [`pivot_totals`](@ref) reads the
 result.
 """
 function reduce_stats!(N::Numeric, S::Symbolic)
     WG = STATS_WORKGROUP
     bm = batch_map(N)
-    _reduce_stats_kernel!(KernelAbstractions.get_backend(N.stats), WG)(N.totals, N.stats, nsupernodes(S), bm,
-                                                                       Val(WG), Val(FRONT_STATS_FIELDS),
+    _reduce_stats_kernel!(KernelAbstractions.get_backend(N.stats), WG)(N.totals, N.info, N.stats, S.super_ptr,
+                                                                       nsupernodes(S), bm, Val(WG),
+                                                                       Val(FRONT_STATS_FIELDS),
                                                                        Val(_ilog2(WG)); ndrange = WG * bm.nact)
     return N
 end
@@ -1427,23 +1451,36 @@ function _factorize_ldlt_herm!(N::Numeric{T}, S::Symbolic, nzval, opts, prm, gim
 end
 
 """
-    factorize_ldlt!(numeric, symbolic, nzval; opts = Options()) -> 0
+    factorize_ldlt!(numeric, symbolic, nzval; impl = :auto, opts = Options(), nb = LDLT_C_NB) -> info::Int
 
 Multifrontal `P A Pᵀ = L D Lᴴ` (structure `"H"`, or `"S"` with real `T`) or
 `L D Lᵀ` (complex symmetric `"S"`) on the device, with the in-front pivoting
 and the static perturbation of the reference [`ref_ldlt!`](@ref) (same pivot
 choice, same storage): `opts.pivot_type`, `pivot_threshold`, `pivot_epsilon`,
 `pivot_epsilon_alg` and `pivot_sign` (copied to `numeric.psign`; a length
 other than `n` raises [`InvalidValueError`](@ref)). Regime-A groups are one
-[`subtree_ldlt_kernel!`](@ref) launch per budget class; every regime-B and
-regime-C launch group is one [`front_ldlt_kernel!`](@ref) launch (fused
-assembly and factorization, one workgroup per front, `F₁₁` in local memory for
-regime-B width classes; no vendor calls, see the T15 report); then [`reduce_stats!`](@ref). `opts.user_host_interrupt` is polled
+[`subtree_ldlt_kernel!`](@ref) launch per budget class; a regime-B launch
+group is one [`front_ldlt_kernel!`](@ref) launch (fused assembly and
+factorization, one workgroup per front, `F₁₁` in local memory above the narrow
+width classes); regime-C groups and the wide regime-B bins of
+[`ldlt_blocked_path`](@ref) take the blocked path of `src/numeric/ldlt_c.jl`
+(`panel_ldlt_kernel!` per block of `nb` pivot columns, then GEMMs through the
+dense interface: `impl = :auto` means vendor GEMMs when the analysis chose
+vendor calls for regime C, else `:ka`; other values as in
+[`select_impl`](@ref); `nb < LDLT_C_NB` is for tests). Then
+[`reduce_stats!`](@ref). `opts.user_host_interrupt` is polled
 before every launch group ([`InterruptedError`](@ref)). With `pivot_epsilon_alg = "algo1"`
 [`abs_max!`](@ref) computes the scale first. Fills `numeric.factor`
-(unit-lower panels), `d`, `piv`, `pivot_kind`, `stats` and `totals`. The
-factorization always completes (`info = 0`); the phase allocates nothing on
-the device, never synchronizes with the host, and is deterministic.
+(unit-lower panels), `d`, `piv`, `pivot_kind`, `stats` and `totals`.
+
+Returns `info`: an exactly zero 1×1 pivot is always perturbed, and with an
+effective `ε = 0` it stays zero; the factorization then fails there (the
+following divisions are left to IEEE arithmetic) and `info` is the original
+column of the failed pivot with the smallest factor column, as for
+[`ref_ldlt!`](@ref). With `ε > 0` the factorization always completes
+(`info = 0`). The phase allocates nothing on the device and is deterministic;
+it synchronizes with the host only to read `info` when the factorization can
+fail (`pivot_epsilon = 0`, or `pivot_epsilon_alg = "algo1"`), otherwise never.
 """
 function factorize_ldlt!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector; impl::Symbol = :auto,
                          opts::Options = Options(), nb::Integer = LDLT_C_NB) where {T}
@@ -1472,5 +1509,5 @@ function factorize_ldlt!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector; impl
     end
     assemble_schur!(N, S, nzval)
     reduce_stats!(N, S)
-    return 0
+    return _numeric_info!(N, S, iszero(prm.eps) || prm.scaled)
 end

src/numeric/ldlt_c.jl

@@ -355,7 +355,8 @@ end
                 x = Int64(ctl[_LT_STAT + q])
                 stats[base + q] = b == 1 ? x : stats[base + q] + x
             end
-            stats[base + STAT_INFO] = 0
+            fi = Int64(ctl[_ST_STATUS])                   # the first failed column over the blocks
+            stats[base + STAT_INFO] = b == 1 || stats[base + STAT_INFO] == 0 ? fi : stats[base + STAT_INFO]
         end
     end
     return nothing

src/numeric/lu.jl

@@ -121,11 +121,12 @@ end
         ε = _lu_eps(p, aux)
         dk = _fget(Fl, k, k)
         kind = PIVOT_KIND_1X1
-        if !(abs(dk) >= ε)                                # tiny (or NaN): perturb
+        if !(abs(dk) >= ε) || iszero(dk)                  # tiny (or NaN, or exactly zero): perturb
             iszero(dk) && (ctl[_LT_STAT + STAT_NZERO] += one(IT))
             dk = _perturbation_sign(dk, 0, !(T <: Complex)) * ε
             kind = PIVOT_KIND_PERTURBED
             ctl[_LT_STAT + STAT_NPERTURBED] += one(IT)
+            iszero(dk) && ctl[_ST_STATUS] == 0 && (ctl[_ST_STATUS] = k % IT)   # ε = 0: the front fails here
         end
         _fset!(Fl, k, k, dk)
         d[g] = dk
@@ -534,7 +535,7 @@ function _launch_subtrees_lu!(N::Numeric{T}, S::Symbolic, nzval, first, count, p
 end
 
 """
-    factorize_lu!(numeric, symbolic, nzval; opts = Options()) -> 0
+    factorize_lu!(numeric, symbolic, nzval; opts = Options()) -> info::Int
 
 Multifrontal `P_r P A Pᵀ = L D U` of a general matrix (structure `"G"`) on the
 device, with the in-front row pivoting and the static perturbation of the
@@ -548,9 +549,13 @@ workgroup per front; no vendor calls); then [`reduce_stats!`](@ref).
 ([`InterruptedError`](@ref)). With `pivot_epsilon_alg = "algo1"`
 [`abs_max!`](@ref) computes the scale first. Fills `numeric.factor` (`L`),
 `numeric.ufactor` (`Uᵀ`), `d`, `piv` (local row order), `pivot_kind`, `stats`
-and `totals`. The factorization always completes (`info = 0`); the phase
-allocates nothing on the device, never synchronizes with the host, and is
-deterministic. Uniform batch as [`factorize!`](@ref).
+and `totals`. Returns `info` as [`factorize_ldlt!`](@ref): an exactly zero
+pivot stays zero with an effective `ε = 0` and the factorization fails there
+(`info` = the original column, as [`ref_lu!`](@ref)); with `ε > 0` it always
+completes (`info = 0`). The phase allocates nothing on the device, is
+deterministic, and synchronizes with the host only to read `info` when the
+factorization can fail (`pivot_epsilon = 0`, or `pivot_epsilon_alg = "algo1"`).
+Uniform batch as [`factorize!`](@ref).
 """
 function factorize_lu!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector; opts::Options = Options()) where {T}
     _check_numeric(N, S, nzval)
@@ -572,5 +577,5 @@ function factorize_lu!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector; opts::
     end
     assemble_schur!(N, S, nzval)
     reduce_stats!(N, S)
-    return 0
+    return _numeric_info!(N, S, iszero(prm.eps) || prm.scaled)
 end