From 0f1c9f2f4fd8f85088009737b82d9463b45aa6be Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 7 Oct 2026 21:11:33 +0000 Subject: [PATCH 1/3] Fix: view-free G, zero-pivot info, analysis-time pair semantics, overflow guards, docstrings (#84) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01AvwMxoGmJm3eSmpzxAy72A --- TASKS.md | 2 + src/matrix.jl | 10 ++++- src/numeric/assembly.jl | 8 ++-- src/numeric/factorize.jl | 30 ++++++++++--- src/numeric/ldlt.jl | 85 ++++++++++++++++++++++++++----------- src/numeric/ldlt_c.jl | 3 +- src/numeric/lu.jl | 17 +++++--- src/numeric/storage.jl | 6 ++- src/options.jl | 15 +++++-- src/reference/ldlt.jl | 32 +++++++++++--- src/reference/lu.jl | 18 +++++--- src/solve/refinement.jl | 18 ++++---- src/solve/sweeps.jl | 8 ++-- src/solver.jl | 48 +++++++++++++++------ src/symbolic/layout.jl | 9 +++- src/symbolic/maps.jl | 5 ++- src/symbolic/ordering.jl | 5 +++ src/symbolic/pairs.jl | 5 ++- src/symbolic/pattern.jl | 10 ++--- test/ported/cudss_solver.jl | 32 ++++++++++---- test/test_api.jl | 22 +++++++++- test/test_matrix.jl | 5 +++ test/test_numeric_ldlt.jl | 35 +++++++++++++++ test/test_numeric_lu.jl | 43 +++++++++++++++++-- test/test_symbolic_etree.jl | 7 ++- 25 files changed, 369 insertions(+), 109 deletions(-) diff --git a/TASKS.md b/TASKS.md index 62d161d..4eac47a 100644 --- a/TASKS.md +++ b/TASKS.md @@ -2718,6 +2718,8 @@ for no gain; the matching-based pairs of #67 must do better there. Decide in 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. +#86: compare the device and the reference on the K2 dumps after scaling (inertia, nperturbed, piv); close #86 if +scaling removed the divergence, else 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 diff --git a/src/matrix.jl b/src/matrix.jl index 4983444..bfd720f 100644 --- a/src/matrix.jl +++ b/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 diff --git a/src/numeric/assembly.jl b/src/numeric/assembly.jl index cb7589f..cceb08d 100644 --- a/src/numeric/assembly.jl +++ b/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] diff --git a/src/numeric/factorize.jl b/src/numeric/factorize.jl index 4865d4d..f784977 100644 --- a/src/numeric/factorize.jl +++ b/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 diff --git a/src/numeric/ldlt.jl b/src/numeric/ldlt.jl index 043874d..b8ebc39 100644 --- a/src/numeric/ldlt.jl +++ b/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,6 +1336,7 @@ 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 @@ -1327,6 +1344,10 @@ 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 """ @@ -1334,15 +1355,18 @@ end 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,7 +1451,7 @@ 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 @@ -1435,15 +1459,28 @@ 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 diff --git a/src/numeric/ldlt_c.jl b/src/numeric/ldlt_c.jl index 12f3105..b6d3f5c 100644 --- a/src/numeric/ldlt_c.jl +++ b/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 diff --git a/src/numeric/lu.jl b/src/numeric/lu.jl index c9b3c7a..275e3bc 100644 --- a/src/numeric/lu.jl +++ b/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 diff --git a/src/numeric/storage.jl b/src/numeric/storage.jl index dcc1847..ecf85e7 100644 --- a/src/numeric/storage.jl +++ b/src/numeric/storage.jl @@ -197,7 +197,8 @@ statistics and their totals, status vector, pivot order, pivot kinds, pivot sign `"G"` the `Uᵀ` panels and update stack) of `symbolic`'s [`Layout`](@ref) for element type `T` on the KernelAbstractions `backend`, for a uniform batch of `nbatch` matrices (every length times `nbatch`, all members active), and build its [`NumericPlan`](@ref). -This is the only allocation of the numeric phase. +This is the only allocation of the numeric phase. Raises [`InvalidValueError`](@ref) when the +panels or the regime-C workspace of the batch do not fit the index type `INT` of `symbolic`. """ function allocate_numeric(S::Symbolic{INT}, ::Type{T}, backend::KernelAbstractions.Backend = KernelAbstractions.CPU(); nbatch::Integer = 1) where {INT, T} @@ -207,6 +208,9 @@ function allocate_numeric(S::Symbolic{INT}, ::Type{T}, backend::KernelAbstractio S.n < typemax(Int32) || throw(InvalidValueError("n = $(S.n) does not fit the Int32 status vector")) nb * L.factor_len < typemax(INT) && nb * (nsupernodes(S) + 1) < typemax(Int32) || throw(InvalidValueError("a batch of $nb factors of $(L.factor_len) entries does not fit the index type $INT")) + nb * L.work_len < typemax(INT) || + throw(InvalidValueError("a batch of $nb regime-C workspaces of $(L.work_len) entries does not fit the " * + "index type $INT; use Int64 indices")) factor = KernelAbstractions.zeros(backend, T, nb * L.factor_len) d = KernelAbstractions.zeros(backend, T, nb * L.d_len) stack = KernelAbstractions.zeros(backend, T, nb * L.stack_len) diff --git a/src/options.jl b/src/options.jl index 2b5d928..c20e946 100644 --- a/src/options.jl +++ b/src/options.jl @@ -109,7 +109,7 @@ defaults below. Keyword arguments are applied through [`setparam!`](@ref), so | `ir_n_steps` | `0` | no refinement (cuDSS parity) | | `ir_tol` | `0.0` | no early exit | | `pivot_type` | `PIVOT_AUTO` | Bunch–Kaufman for `S`/`H`, none for `SPD`/`HPD` | -| `pivot_threshold` | `0.01` | pivot acceptance threshold | +| `pivot_threshold` | `0.01` | pivot acceptance threshold; with `pivot_pairs = "default"` also the partner threshold of the 2×2 pairs chosen at analysis (a value set after `"analysis"` changes the in-front test, not the pairs already chosen) | | `pivot_epsilon` | `nothing` | [`default_pivot_epsilon`](@ref)`(T)` | | `pivot_epsilon_alg` | `PIVOT_EPSILON_DEFAULT` | | | `max_lu_nnz` | `-1` | no limit (any negative value) | @@ -129,7 +129,7 @@ defaults below. Keyword arguments are applied through [`setparam!`](@ref), so | `factor_precision` | `nothing` | factors in the input precision | | `amalgamation` | `(max_width = 32, zero_fraction = 0.25, min_width = 8)` | | | `schedule` | `SCHEDULE_AUTO` | | -| `pivot_pairs` | `PIVOT_PAIRS_DEFAULT` | 2×2 pivot pairs in the analysis of `"S"`/`"H"` (structurally zero pivots; `"all"`: every candidate; with matching `"algo5"`/`"algo6"`: from the cycles of the matching) | +| `pivot_pairs` | `PIVOT_PAIRS_DEFAULT` | 2×2 pivot pairs in the analysis of `"S"`/`"H"` (structurally zero pivots; `"all"`: every candidate; with matching `"algo5"`/`"algo6"`: from the cycles of the matching); decided from the values present at analysis: an all-zero `nzval` gives no pairs, an undefined one arbitrary pairs | | `pivot_pair_tolerance` | `1e-6` ([`PIVOT_PAIR_TOLERANCE`](@ref)) | relative diagonal size below which a row is a 2×2 candidate | | `regime_c_width` | `64` | fronts wider than this go to regime C (vendor dense calls) | | `regime_c_rows` | `512` | fronts with more rows than this go to regime C | @@ -541,7 +541,16 @@ Accepted values: or `"none"`; ignored for the other structures, with `user_perm` and with the natural ordering; with `matching_alg = "algo5"`/`"algo6"` the pairs come from the cycles of the - scaled symmetric matching instead ([`matching_pairs`](@ref)); + scaled symmetric matching instead ([`matching_pairs`](@ref)). + Pairs are decided from the values present at `"analysis"` (the first batch + member): an all-zero `nzval` gives no pairs and an undefined one arbitrary pairs, + so run `"analysis"` after the first assembly of the matrix (MadNLP: after the + first KKT assembly); `"all"` reads the values as well (its candidates and + partners come from them); +* `"pivot_threshold"` also shapes the ordering with `pivot_pairs = "default"`: + partners are accepted at analysis down to `pivot_threshold · maxₖ |aᵢₖ|` + ([`zero_pivot_pairs!`](@ref)); a value set after `"analysis"` changes the + in-front pivot test but not the pairs already chosen; * `"pivot_pair_tolerance"` (beyond cuDSS): a finite real `≥ 0`, the relative diagonal size `τ` below which a row is a 2×2 candidate (`|aᵢᵢ| ≤ τ maxⱼ≠ᵢ |aᵢⱼ|`, on the scaled matrix with matching; default [`PIVOT_PAIR_TOLERANCE`](@ref)); diff --git a/src/reference/ldlt.jl b/src/reference/ldlt.jl index 24cb859..1b1f06d 100644 --- a/src/reference/ldlt.jl +++ b/src/reference/ldlt.jl @@ -109,8 +109,14 @@ Fills `numeric.factor` (unit-lower panels), `numeric.d`, `numeric.piv`, nperturbed, n2x2, info`; `nzero` counts pivots that were exactly zero before their perturbation; `npos`/`nneg` count the signs of D after perturbation, a 2×2 block with negative determinant counting one of each, and stay `0` for complex -symmetric matrices, which have no inertia). Returns `info = 0`: with static -perturbation the factorization always completes. +symmetric matrices, which have no inertia). An exactly zero 1×1 pivot is +always tiny; with `ε = 0` its perturbation leaves it zero and the +factorization fails there: the front's `info` statistic is the failed local +column (assembled order, the first in pivot order), the following divisions +are left to IEEE arithmetic (the factor is unusable), and the return value is +the original column of the failed pivot with the smallest factor column, as for +the Cholesky oracle. Otherwise returns `info = 0`: with static perturbation +`ε > 0` the factorization always completes. """ function ref_ldlt!(N::Numeric{T, Vector{T}}, S::Symbolic, nzval::AbstractVector; opts::Options = Options()) where {T} length(nzval) == S.nnz || @@ -172,7 +178,7 @@ function ref_ldlt!(N::Numeric{T, Vector{T}}, S::Symbolic, nzval::AbstractVector; N.d[g] = dd[k] N.d[n + g] = ee[k] end - for q in 1:5 + for q in 1:FRONT_STATS_FIELDS N.stats[(s - 1) * FRONT_STATS_FIELDS + q] = st[q] end if m > 0 @@ -192,7 +198,19 @@ function ref_ldlt!(N::Numeric{T, Vector{T}}, S::Symbolic, nzval::AbstractVector; end end _host_totals!(N) - return 0 + return _host_info(N, S) +end + +# the original column of the failed pivot with the smallest factor column (`super_ptr[s] + info - 1`, +# `info` = the failed assembled local column of front `s`; 0 = none), as the device reduction +function _host_info(N::Numeric, S::Symbolic) + sp = S.partition + best = 0 + for s in 1:nsupernodes(sp) + fi = Int(N.stats[s * FRONT_STATS_FIELDS]) + fi > 0 && (best = best == 0 ? sp.super_ptr[s] + fi - 1 : min(best, sp.super_ptr[s] + fi - 1)) + end + return best == 0 ? 0 : sp.perm[best] end # the totals of the per-front statistics (host storage; `reduce_stats!` on the device) @@ -369,7 +387,7 @@ function _ldlt_front!(F::Matrix{T}, w::Int, prm::_LDLTParams, signof) where {T} kind = zeros(Int8, w) dd = zeros(T, w) ee = zeros(T, w) - st = zeros(Int, 5) + st = zeros(Int, FRONT_STATS_FIELDS) # npos, nneg, nzero, nperturbed, n2x2, info k = 1 while k <= w c, r = _choose_pivot(F, k, w, prm) @@ -377,11 +395,13 @@ function _ldlt_front!(F::Matrix{T}, w::Int, prm::_LDLTParams, signof) where {T} _swap_front!(F, k, c, w, lp) d = prm.herm ? T(real(F[k, k])) : F[k, k] kind[k] = PIVOT_KIND_1X1 - if !(abs(d) >= prm.eps) # tiny (or NaN): perturb + if !(abs(d) >= prm.eps) || iszero(d) # tiny (or NaN, or exactly zero): perturb iszero(d) && (st[STAT_NZERO] += 1) d = _perturbation_sign(d, signof(lp[k]), prm.herm) * prm.eps kind[k] = PIVOT_KIND_PERTURBED st[STAT_NPERTURBED] += 1 + # ε = 0 leaves the pivot zero: the factorization fails at this (assembled) column + iszero(d) && st[STAT_INFO] == 0 && (st[STAT_INFO] = lp[k]) end F[k, k] = d dd[k] = d diff --git a/src/reference/lu.jl b/src/reference/lu.jl index 07a812c..9716e34 100644 --- a/src/reference/lu.jl +++ b/src/reference/lu.jl @@ -50,7 +50,7 @@ function _lu_params(S::Symbolic, ::Type{T}, opts::Options, nz::AbstractVector) w end """ - ref_lu!(numeric::Numeric{T, Vector{T}}, symbolic, nzval; opts = Options()) -> 0 + ref_lu!(numeric::Numeric{T, Vector{T}}, symbolic, nzval; opts = Options()) -> info Reference multifrontal `L D U` of the general matrix whose stored values are `nzval` (structure `"G"`, same CSR pattern and index base as the analysis) @@ -75,8 +75,11 @@ Fills `numeric.factor` (unit-lower `L` panels), `numeric.ufactor` (unit-lower `Uᵀ` panels), `numeric.d`, `numeric.piv` (local row order), `numeric.pivot_kind` (`PIVOT_KIND_1X1` or `PIVOT_KIND_PERTURBED`) and the per-front `numeric.stats` (`nzero`, `nperturbed`; `npos`, `nneg` and `n2x2` stay 0: a general matrix has -no inertia). Returns 0: with static perturbation the factorization always -completes. +no inertia). An exactly zero pivot is always tiny; with `ε = 0` it stays zero +and the factorization fails there: the front's `info` statistic is the failed +local column and the return value the original column of the failed pivot with +the smallest factor column (as [`ref_ldlt!`](@ref)). Otherwise returns 0: with +static perturbation `ε > 0` the factorization always completes. """ function ref_lu!(N::Numeric{T, Vector{T}}, S::Symbolic, nzval::AbstractVector; opts::Options = Options()) where {T} length(nzval) == S.nnz || @@ -133,7 +136,7 @@ function ref_lu!(N::Numeric{T, Vector{T}}, S::Symbolic, nzval::AbstractVector; o N.pivot_kind[g] = kind[k] N.d[g] = dd[k] end - for q in 1:5 + for q in 1:FRONT_STATS_FIELDS N.stats[(s - 1) * FRONT_STATS_FIELDS + q] = st[q] end if m > 0 @@ -158,7 +161,7 @@ function ref_lu!(N::Numeric{T, Vector{T}}, S::Symbolic, nzval::AbstractVector; o end end _host_totals!(N) - return 0 + return _host_info(N, S) end ref_lu!(N::Numeric, S::Symbolic, A::CSR; kwargs...) = ref_lu!(N, S, vec(A.nzval); kwargs...) @@ -173,7 +176,7 @@ function _lu_front!(F::AbstractMatrix{T}, w::Int, prm::_LUParams) where {T} lp = collect(1:w) kind = fill(PIVOT_KIND_1X1, w) dd = zeros(T, w) - st = zeros(Int64, 5) + st = zeros(Int64, FRONT_STATS_FIELDS) # nzero, nperturbed (and the zeros), info for k in 1:w r = prm.pivot ? _lu_choose_row(F, k, w, prm.u) : k if r != k @@ -183,11 +186,12 @@ function _lu_front!(F::AbstractMatrix{T}, w::Int, prm::_LUParams) where {T} lp[k], lp[r] = lp[r], lp[k] end d = F[k, k] - if !(abs(d) >= prm.eps) # tiny (or NaN): perturb + if !(abs(d) >= prm.eps) || iszero(d) # tiny (or NaN, or exactly zero): perturb iszero(d) && (st[STAT_NZERO] += 1) d = _perturbation_sign(d, 0, !(T <: Complex)) * prm.eps kind[k] = PIVOT_KIND_PERTURBED st[STAT_NPERTURBED] += 1 + iszero(d) && st[STAT_INFO] == 0 && (st[STAT_INFO] = k) # ε = 0: failed column end F[k, k] = d dd[k] = d diff --git a/src/solve/refinement.jl b/src/solve/refinement.jl index 5cf92a4..edcb61e 100644 --- a/src/solve/refinement.jl +++ b/src/solve/refinement.jl @@ -366,7 +366,7 @@ function refine!(X::AbstractVecOrMat, B::AbstractVecOrMat, W::RefinementWorkspac max_rhs(W) >= ncols && max_rhs(ws) >= ncols || throw(DimensionMismatch("the refinement workspace holds $(max_rhs(W)) right-hand sides, need $ncols")) bm = batch_map(N; nrhs) - pin, sin, pout, sout = solve_io(S, scaling, transpose_matrix) + pin, scin, pout, sout = solve_io(S, scaling, transpose_matrix) steps = 0 progress[] = 0 for _ in 1:nsteps @@ -378,7 +378,7 @@ function refine!(X::AbstractVecOrMat, B::AbstractVecOrMat, W::RefinementWorkspac _log(LOG_DEBUG, () -> "refinement: step $steps, relative residual $rel") rel <= tol && break end - permute_rhs!(ws.Y, W.R, pin; conjugate, scale = sin) + permute_rhs!(ws.Y, W.R, pin; conjugate, scale = scin) forward_sweep!(ws, S, N; nrhs, deterministic, transpose = transpose_matrix) diagonal_sweep!(ws, S, N; nrhs) backward_sweep!(ws, S, N; nrhs, transpose = transpose_matrix) @@ -444,7 +444,8 @@ RefinementOperator(W::RefinementWorkspace, nzval::AbstractVector, n::Integer, nc conjugate::Bool) = RefinementOperator(W, nzval, Int(n), Int(ncols), bm, conjugate, 0) Base.size(op::RefinementOperator) = (op.n * op.ncols, op.n * op.ncols) -Base.size(op::RefinementOperator, d::Integer) = d <= 2 ? op.n * op.ncols : 1 +Base.size(op::RefinementOperator, d::Integer) = + d < 1 ? throw(BoundsError(size(op), d)) : d <= 2 ? op.n * op.ncols : 1 Base.eltype(::RefinementOperator{T}) where {T} = T @kernel function _spmv_kernel!(y, rowptr, colval, src, nzval_all, x, n, ncols, bm, roff, ::Val{CJ}) where {CJ} @@ -492,7 +493,7 @@ permutations and scalings of a [`SolveScaling`](@ref) `scaling` ([`solve_io`](@r `interrupt` first ([`InterruptedError`](@ref)); asynchronous otherwise. The right preconditioner of FGMRES-IR. """ -struct FactorPreconditioner{T, WS <: SolveWorkspace{T}, SY <: Symbolic, NU <: Numeric} +struct FactorPreconditioner{T, WS <: SolveWorkspace{T}, SY <: Symbolic, NU <: Numeric, SC} ws::WS S::SY N::NU @@ -502,7 +503,7 @@ struct FactorPreconditioner{T, WS <: SolveWorkspace{T}, SY <: Symbolic, NU <: Nu deterministic::Bool interrupt::Union{Nothing, Threads.Atomic{Bool}} transpose::Bool - scaling::Any + scaling::SC # a SolveScaling, or nothing end FactorPreconditioner(ws::SolveWorkspace, S::Symbolic, N::Numeric, nrhs::Integer, ncols::Integer, conjugate::Bool, @@ -510,13 +511,14 @@ FactorPreconditioner(ws::SolveWorkspace, S::Symbolic, N::Numeric, nrhs::Integer, FactorPreconditioner(ws, S, N, Int(nrhs), Int(ncols), conjugate, deterministic, interrupt, transpose, scaling) Base.size(P::FactorPreconditioner) = (P.S.n * P.ncols, P.S.n * P.ncols) -Base.size(P::FactorPreconditioner, d::Integer) = d <= 2 ? P.S.n * P.ncols : 1 +Base.size(P::FactorPreconditioner, d::Integer) = + d < 1 ? throw(BoundsError(size(P), d)) : d <= 2 ? P.S.n * P.ncols : 1 Base.eltype(::FactorPreconditioner{T}) where {T} = T function LinearAlgebra.mul!(y::AbstractVector, P::FactorPreconditioner, x::AbstractVector) _poll_interrupt(P.interrupt) - pin, sin, pout, sout = solve_io(P.S, P.scaling, P.transpose) - permute_rhs!(P.ws.Y, x, pin; conjugate = P.conjugate, scale = sin) + pin, scin, pout, sout = solve_io(P.S, P.scaling, P.transpose) + permute_rhs!(P.ws.Y, x, pin; conjugate = P.conjugate, scale = scin) forward_sweep!(P.ws, P.S, P.N; nrhs = P.nrhs, deterministic = P.deterministic, transpose = P.transpose) diagonal_sweep!(P.ws, P.S, P.N; nrhs = P.nrhs) backward_sweep!(P.ws, P.S, P.N; nrhs = P.nrhs, transpose = P.transpose) diff --git a/src/solve/sweeps.jl b/src/solve/sweeps.jl index c096320..26716bc 100644 --- a/src/solve/sweeps.jl +++ b/src/solve/sweeps.jl @@ -867,7 +867,8 @@ function _bwd_dense!(ws::SolveWorkspace{T}, S::Symbolic, N::Numeric, s::Int, nrh end """ - forward_sweep!(ws, symbolic, numeric; nrhs = max_rhs(ws), deterministic = false, impl = :auto) -> ws + forward_sweep!(ws, symbolic, numeric; nrhs = max_rhs(ws), deterministic = false, impl = :auto, + transpose = false) -> ws Forward sweep `L Z = Y` in place on the first `nrhs` columns of `ws.Y` (the permuted right-hand side, [`permute_rhs!`](@ref)), with the factor of @@ -950,7 +951,7 @@ function diagonal_sweep!(ws::SolveWorkspace, S::Symbolic, N::Numeric; nrhs::Inte end """ - backward_sweep!(ws, symbolic, numeric; nrhs = max_rhs(ws), impl = :auto) -> ws + backward_sweep!(ws, symbolic, numeric; nrhs = max_rhs(ws), impl = :auto, transpose = false) -> ws Backward sweep `Lᴴ X = W` (`Lᵀ` for complex symmetric LDLᵀ) in place on the first `nrhs` columns of `ws.Y`, the forward plan in reverse: per step the regime-C-path fronts (gather kernel, @@ -985,7 +986,8 @@ function backward_sweep!(ws::SolveWorkspace, S::Symbolic, N::Numeric; nrhs::Inte end """ - sweep_solve!(X, ws, symbolic, numeric, B; transposed = false, deterministic = false, impl = :auto) -> X + sweep_solve!(X, ws, symbolic, numeric, B; transposed = false, deterministic = false, impl = :auto, + transpose = false) -> X Solve `A X = B` with the factor `P A Pᵀ = L Lᴴ` (or the LDLᵀ/LDLᴴ factor with its local pivot orders) of [`factorize!`](@ref): diff --git a/src/solver.jl b/src/solver.jl index 46da173..f5fa126 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -109,7 +109,9 @@ Sparse direct solver handle (≅ `CudssSolver`), PLAN §3.1–§3.2. DirectSolver(A::CuSparseMatrixCSC, structure, view; index = 'O') # CUDA extension `structure` is `"G"`, `"S"`, `"H"`, `"SPD"` or `"HPD"` and `view` is `'L'`, -`'U'` or `'F'` (the triangle of the matrix that is read), as in CUDSS.jl; +`'U'` or `'F'` (the triangle of the matrix that is read; `'F'` reads the lower +triangle of a symmetric or Hermitian structure, and for `"G"` the view is ignored +and every stored entry is read, as in cuDSS), as in CUDSS.jl; `index` (`'O'`/`'Z'`) is the base of the row pointers and column indices. The CSR arrays are wrapped without copies and live on a KernelAbstractions backend (the solver's `backend`), except for a `SparseMatrixCSC`, which is converted to @@ -123,7 +125,7 @@ right-hand side and the solution (as for `solve_mode`, see Implemented at this point: structures `"SPD"` (real `T`), `"HPD"`, `"S"` (LDLᵀ; complex symmetric for complex `T`), `"H"` (LDLᴴ) and `"G"` (LU on the symmetric pattern of `A + Aᵀ` with row pivoting inside the fully-summed block -of each front, `P_r P A Pᵀ = L D U`; view `'F'` only), single matrices +of each front, `P_r P A Pᵀ = L D U`; any view, ignored), single matrices and uniform batches (see below), the phases `"reordering"`, `"symbolic_factorization"`, `"analysis"`, `"factorization"`, `"refactorization"`, `"solve"`, `"solve_fwd_perm"`, `"solve_fwd"`, `"solve_diag"`, `"solve_bwd"`, @@ -254,7 +256,8 @@ DirectSolver(A::SparseMatrixCSC, structure, view; index = INDEX_ONE) = DirectSolver(CSR(A; index), structure, view) Base.size(solver::DirectSolver) = (solver.A.nrows, solver.A.nrows) -Base.size(solver::DirectSolver, d::Integer) = d <= 2 ? solver.A.nrows : 1 +Base.size(solver::DirectSolver, d::Integer) = + d < 1 ? throw(BoundsError(size(solver), d)) : d <= 2 ? solver.A.nrows : 1 const _STAGE_NAMES = ("created", "reordered", "analyzed", "factorized") @@ -343,8 +346,9 @@ Execute `phase` on `solver` (≅ `cudss(phase, solver, x, b)`), PLAN §1.2: * `"analysis"`: both; * `"factorization"`: numeric factorization with the current values of the matrix (needs the analysis); sets `"info"` (`0` or, for `"SPD"`/`"HPD"`, the - original column of the first non-positive pivot; LDLᵀ/LDLᴴ perturbs tiny - pivots instead and always reports `0`, see `"npivots"`) and clears + original column of the first non-positive pivot; LDLᵀ/LDLᴴ and LU perturb tiny + pivots instead and report `0`, see `"npivots"`, except for an exactly zero + pivot with `pivot_epsilon = 0`, whose original column they report) and clears `fresh_factorization`; * `"refactorization"`: the same with the analysis and storage reused, `"info"` reset first (needs a previous `"factorization"`); @@ -882,7 +886,9 @@ flags 0/1, or `nothing`) select the members of the next phases; the index and the mask are checked against the batch size when a phase runs. `"pivot_sign"` (PLAN §1.7: a vector of `n` entries in `(-1, 0, 1)`, host or device, or `nothing`) is checked against the size of the matrix; the next -`"factorization"`/`"refactorization"` copies it to the device. Setting +`"factorization"`/`"refactorization"` copies it to the device; it steers the +perturbation of LDLᵀ/LDLᴴ only and is ignored for `"G"` (LU perturbs with the +sign of the pivot) and the Cholesky structures. Setting `"ir_n_steps"` (the number of refinement steps requested) makes [`getparam`](@ref) report that value again until the next solve. `"schur_matrix"` (Schur complement mode, after `"reordering"`/`"analysis"`) @@ -936,7 +942,7 @@ The data parameters computed by the solver: | name | value | available after | | --- | --- | --- | -| `"info"` | `Int`: `0`, or the original column (1-based) of the first non-positive pivot | always | +| `"info"` | `Int`: `0`, or the original column (1-based) of the first failed pivot (Cholesky: non-positive; LDLᵀ/LDLᴴ and LU: exactly zero with `pivot_epsilon = 0`, the other tiny pivots are perturbed) | always | | `"ubatch_mask"` | the `"ubatch_mask"` set (or `nothing`) | always | | `"ir_n_steps"` | `Int`: refinement steps performed by the last `"solve"`/`"solve_refinement"` (≤ the configured `ir_n_steps`, fewer when `ir_tol` stopped early); the configured value before the first solve and after `setparam!(solver, "ir_n_steps", k)` | always | | `"lu_nnz"` | `Int64`: nonzeros of `L` (diagonal included, amalgamation zeros excluded) | analysis | @@ -944,20 +950,23 @@ The data parameters computed by the solver: | `"nsuperpanels"` | `Int`: supernodes after amalgamation | analysis | | `"memory_estimates"` | `Vector{Int64}` (16 entries, see [`memory_estimates`](@ref)) | analysis | | `"perm_reorder_row"`, `"perm_reorder_col"` | `Vector{Int}`: the fill-reducing permutation, 1-based (`perm[k]` = original index of the `k`-th pivot) | reordering | -| `"perm_row"`, `"perm_col"` | `Vector{Int}`: the final permutation of the factor (= the reordering for Cholesky and LDLᵀ/LDLᴴ); LU (`"G"`) after a factorization: `perm_row[k]` is the original row of factor row `k` (the reordering composed with the local row interchanges of batch member 1), `perm_col` the reordering (with matching: composed with the matching, `perm_matching[reordering]`), so `A[perm_row, perm_col] = L D U` (with matching `Dr A[perm_row, perm_col] Dc = L D U`, the scalings in that order) | analysis | +| `"perm_row"`, `"perm_col"` | `Vector{Int}`: the final permutation of the factor (= the reordering for Cholesky and LDLᵀ/LDLᴴ); LU (`"G"`) after a factorization: `perm_row[k]` is the original row of factor row `k` (the reordering composed with the local row interchanges of batch member 1), `perm_col` the reordering (with matching: composed with the matching, `perm_matching[reordering]`), so `A[perm_row, perm_col] = L D U` (with matching `Dr[perm_row] A[perm_row, perm_col] Dc[perm_col] = L D U`, the scalings `Dr`, `Dc` of `"scale_row"`/`"scale_col"` as diagonal matrices in the original numbering). For a CSC input (`CuSparseMatrixCSC`, [`csr_of_transpose`](@ref)) `A` is the stored CSR matrix, the transpose of the matrix given | analysis | | `"perm_matching"` | `Vector{Int}`: the matching permutation ([`Matching`](@ref)): row `i` is matched to column `perm_matching[i]` (`"G"`: `A[:, perm_matching]` has the matched entries on its diagonal; symmetric structures: only its cycles are used, for the 2×2 pivot pairs) | reordering (`matching_alg ≠ "default"`) | | `"scale_row"`, `"scale_col"` | `Vector{real(T)}`: the row and column scaling factors `Dr`, `Dc` of the factored matrix, in the original numbering (ones for `"algo1"`–`"algo4"`; equal for the symmetric structures) | reordering (`matching_alg ≠ "default"`) | | `"diag"` | vector of `T` on the solver's backend: the diagonal of `L` (Cholesky), of `D` (LDLᵀ/LDLᴴ; for a 2×2 block its two diagonal entries) or of `D` in `L D U` (LU, = the diagonal of `U` in `L U`) in factor order; Schur complement mode: `1` on the (unfactored) Schur block | factorization | -| `"npivots"` | `INT`: perturbed pivots (LDLᵀ/LDLᴴ, LU; `0` for Cholesky) | factorization | -| `"inertia"` | `Tuple{INT, INT}`: `(npos, nneg)` of D (after perturbation, so the inertia of `A + E`; read it with `"npivots"`); `(0, 0)` for complex symmetric `"S"` and for `"G"`; Cholesky: `(number of positive pivots, 0)` | factorization | -| `"pivot_stats"` | `NamedTuple` `(npos, nneg, nzero, nperturbed, n2x2)` of `Int64` (PLAN §1.7); in Schur complement mode `"npivots"`, `"inertia"` and `"pivot_stats"` count the factored block `A₁₁` only | factorization | +| `"npivots"` | `INT`: perturbed pivots (LDLᵀ/LDLᴴ, LU; `0` for Cholesky); valid when `"info"` is `0` (after a failed front the counts of the fronts above it are stale) | factorization | +| `"inertia"` | `Tuple{INT, INT}`: `(npos, nneg)` of D (after perturbation, so the inertia of `A + E`; read it with `"npivots"`); `(0, 0)` for complex symmetric `"S"` and for `"G"`; Cholesky: `(number of positive pivots, 0)`; valid when `"info"` is `0`, as `"npivots"` | factorization | +| `"pivot_stats"` | `NamedTuple` `(npos, nneg, nzero, nperturbed, n2x2)` of `Int64` (PLAN §1.7); in Schur complement mode `"npivots"`, `"inertia"` and `"pivot_stats"` count the factored block `A₁₁` only; valid when `"info"` is `0`, as `"npivots"` | factorization | | `"schur_shape"` | `(nrows, ncols, nnz)` of `Int64`: the size `ns` of the Schur complement and the nonzeros of its symbolic pattern ([`schur_pattern`](@ref)); for the symmetric structures one triangle with the diagonal (as the sparse export of one triangle) | analysis (`schur_mode = 1`) | | `"schur_matrix"` | the Schur complement `S = A₂₂ − A₂₁ A₁₁⁻¹ A₁₂` of the rows and columns of `"user_schur_indices"` (in their increasing order), written into the destination registered with [`setparam!`](@ref) and returned; without one, a new dense `ns × ns` matrix on the solver's backend. Symmetric structures: a dense destination gets the full matrix, a sparse one the triangle of its view | factorization (`schur_mode = 1`) | Uniform batch (`nbatch > 1`): `"info"`, `"npivots"`, `"inertia"` and `"pivot_stats"` are vectors with one entry per batch member, `"diag"` is the `n · nbatch` vector of the members' diagonals one after the other; the -analysis outputs are shared by the members. `"perm_row"` of an LU batch is +analysis outputs are shared by the members. A factorization restricted by +`"ubatch_index"`/`"ubatch_mask"` updates the active members only: the others +keep the `"info"`, `"npivots"`, `"inertia"` and `"pivot_stats"` of their last +factorization. `"perm_row"` of an LU batch is batch member 1's row order, by design: cuDSS returns the vector outputs of a uniform batch (`"diag"`) for member 1 only and does not implement `"perm_row"` at all (cuDSS 0.8). `"diag"` is the deliberate exception and @@ -1024,7 +1033,20 @@ function getparam(solver::DirectSolver, name::AbstractString) return copy(sp.perm) # perm_row, perm_col end -_memory_estimates(solver::DirectSolver{T, INT}) where {T, INT} = memory_estimates(solver.host_symbolic, T, INT) +# the estimates of the analysis plus, with matching, its device buffers: the scaled values (every member), the +# entry weights, the two scale vectors and, for "G", the column permutation +function _memory_estimates(solver::DirectSolver{T, INT}) where {T, INT} + est = memory_estimates(solver.host_symbolic, T, INT) + m = solver.host_matching + m === nothing && return est + n = solver.A.nrows + extra = Int64(length(solver.A.nzval)) * sizeof(T) + Int64(nnz(solver.A) + 2 * n) * sizeof(real(T)) + + (m.symmetric ? 0 : Int64(n) * sizeof(INT)) + for k in (1, 2, 5) + est[k] += extra + end + return est +end """ getparam!(buffer, solver::DirectSolver, name::String) -> buffer diff --git a/src/symbolic/layout.jl b/src/symbolic/layout.jl index 325550c..7a8d76a 100644 --- a/src/symbolic/layout.jl +++ b/src/symbolic/layout.jl @@ -182,10 +182,15 @@ ldlt_c_work_len(f::Integer, w::Integer, cb::Bool, nb::Integer = LDLT_C_NB) = (cb ? (Int(f) - Int(w))^2 : 0) + 2 * Int(f) * (nb + 1) + 3 * Int(f) """ -Smallest width class and row class of a regime-B bin that the LDLᵀ/LDLᴴ factorization runs on the blocked -regime-C path ([`ldlt_blocked_path`](@ref)) instead of the fused front kernel. +Smallest width class of a regime-B bin that the LDLᵀ/LDLᴴ factorization runs on the blocked regime-C path +([`ldlt_blocked_path`](@ref)) instead of the fused front kernel (together with [`LDLT_BLOCKED_MIN_FCLASS`](@ref)). """ const LDLT_BLOCKED_MIN_WCLASS = 32 + +""" +Smallest row class of a regime-B bin that the LDLᵀ/LDLᴴ factorization runs on the blocked regime-C path +([`ldlt_blocked_path`](@ref)) instead of the fused front kernel (together with [`LDLT_BLOCKED_MIN_WCLASS`](@ref)). +""" const LDLT_BLOCKED_MIN_FCLASS = 256 """ diff --git a/src/symbolic/maps.jl b/src/symbolic/maps.jl index b4a989c..de0e46f 100644 --- a/src/symbolic/maps.jl +++ b/src/symbolic/maps.jl @@ -312,7 +312,10 @@ element type `T` and device indices `INT` (default: the index type of | 12 | largest regime-A local memory in use (per workgroup, [`subtree_local_bytes`](@ref) of its class) | | 13–16 | 0 (reserved) | -Slots 1–6 follow cuDSS's `CUDSS_DATA_MEMORY_ESTIMATES`. +Slots 1–6 follow cuDSS's `CUDSS_DATA_MEMORY_ESTIMATES`. `getparam(solver, "memory_estimates")` +adds to slots 1, 2 and 5 the device buffers of the matching when there is one (`matching_alg`): the +scaled values (`nnz · nbatch` of `T`), the entry weights (`nnz`) and the two scale vectors (`n` +each) of `real(T)`, and for `"G"` the column permutation (`n` of `INT`). """ function memory_estimates(S::Symbolic{INT0}, ::Type{T}, ::Type{INT} = INT0) where {INT0, T, INT} est = zeros(Int64, 16) diff --git a/src/symbolic/ordering.jl b/src/symbolic/ordering.jl index d5d6763..adac1e7 100644 --- a/src/symbolic/ordering.jl +++ b/src/symbolic/ordering.jl @@ -166,6 +166,11 @@ in the analysis when [`pairs_enabled`](@ref): graph compressed by the pairs is re-ordered, until a round adds no pair (at most [`PIVOT_PAIRS_MAX_ROUNDS`](@ref)). Each algorithm keeps its own pairs. +Both come from the values of the matrix at analysis: an all-zero `nzval` gives +no pairs and an undefined one arbitrary pairs, so the analysis must run after the +first assembly of the values (MadNLP: after the first KKT assembly), for +`"default"` and `"all"` alike. + AMD, MMD and ND then order the graph compressed by the pairs ([`compressed_pattern`](@ref): a pair is one vertex with the union adjacency), and the order is expanded with each partner right before its candidate diff --git a/src/symbolic/pairs.jl b/src/symbolic/pairs.jl index 7c7ae5f..c2ff822 100644 --- a/src/symbolic/pairs.jl +++ b/src/symbolic/pairs.jl @@ -339,7 +339,10 @@ end The 2×2 pivot pair input of [`compute_ordering`](@ref) for the analysis: nothing when [`pairs_enabled`](@ref) is false; [`pivot_candidates`](@ref) for `pivot_pairs = "default"` (pairs for the structurally zero pivots); the fixed [`pivot_pairs`](@ref) for -`"all"`. Reads `nzval` on the host once. +`"all"`. Reads `nzval` on the host once. Pairs are decided from the values present at `"analysis"` (the first batch +member): an all-zero `nzval` gives no pairs and an undefined one arbitrary pairs, +so run `"analysis"` after the first assembly of the matrix (MadNLP: after the +first KKT assembly); `"all"` reads the values as well. """ function analysis_pairs(P::SymmetricPattern, rowptr, colval, nzval, n::Integer, structure, opts::Options; view = VIEW_FULL, index = INDEX_ONE) diff --git a/src/symbolic/pattern.jl b/src/symbolic/pattern.jl index de50b10..59d97c1 100644 --- a/src/symbolic/pattern.jl +++ b/src/symbolic/pattern.jl @@ -32,7 +32,7 @@ others build the pattern from a CSR pattern, following PLAN §2.3 step 1: * `structure` `"S"`, `"H"`, `"SPD"`, `"HPD"` (or a [`Structure`](@ref)): only the triangle selected by `view` is read (`'L'`: `col ≤ row`, `'U'`: `col ≥ row`; `'F'` reads the lower triangle, as cuDSS does) and mirrored; -* `"G"`: the pattern of `A + Aᵀ`; `view` must be `'F'`; +* `"G"`: the pattern of `A + Aᵀ`; the view is ignored (every stored entry is read, as in cuDSS); * duplicates are dropped, `index` (`'O'`/`'Z'`) is the base of `rowptr`/`colval`, non-square matrices and out-of-range column indices raise [`InvalidValueError`](@ref). """ @@ -106,11 +106,9 @@ _is_hermitian(s::Structure) = s == STRUCTURE_HERMITIAN || s == STRUCTURE_HPD # Which stored entries (r, c) are read, and whether the mirrored entry is generated. # For symmetric structures 'F' reads the lower triangle, as in cuDSS (PLAN §1.1). function _entry_filter(structure::Structure, view::MatrixView) - if structure == STRUCTURE_GENERAL - view == VIEW_FULL || - throw(InvalidValueError("structure \"G\" needs view 'F', got '$(convert(Char, view))'")) - return (r, c) -> true - end + # "G" ignores the view and reads every stored entry, as in cuDSS ("If the + # accompanying matrix type is CUDSS_MTYPE_GENERAL, the matrix view is ignored"). + structure == STRUCTURE_GENERAL && return (r, c) -> true view == VIEW_UPPER && return (r, c) -> c >= r return (r, c) -> c <= r end diff --git a/test/ported/cudss_solver.jl b/test/ported/cudss_solver.jl index e8c1e02..caf5b61 100644 --- a/test/ported/cudss_solver.jl +++ b/test/ported/cudss_solver.jl @@ -7,9 +7,10 @@ # be rejected at analysis); the matrix is SPD/HPD, so `"inertia"` is `(n, 0)` # for real "S" and "H" (complex symmetric "S": `(0, 0)`, no inertia) and # `"diag"` is real positive except for complex "S" (and has a positive real -# part for "G", whose complex `D` is not exactly real); "G" reads the full -# matrix: views 'L'/'U' are rejected at analysis (`InvalidValueError`; cuDSS -# ignores the view of a general matrix); +# part for "G", whose complex `D` is not exactly real); "G" ignores the view +# and reads every stored entry, as cuDSS does (#84): it is given the full +# matrix with each view, and analysis, factorization and solve give the same +# solution for 'L', 'U' and 'F'; # * `matching_alg = "algo6"` is set before the analysis, as in CUDSS.jl (T21): # `"perm_matching"` is checked to be a permutation and `"scale_row"`/`"scale_col"` # to be positive; the inertia stays `(n, 0)` under matching (cuDSS 0.8 reports @@ -186,14 +187,27 @@ end @testset "structure = $structure" for structure in (T <: Real ? ("G", "S", "H", "SPD", "HPD") : ("G", "S", "H", "HPD")) @testset "view = $view" for view in ('L', 'U', 'F') - A_gpu = api_matrix(backend, triangle_view(A_cpu, view), INT) - if structure == "G" && view != 'F' - @test_throws InvalidValueError execute!("analysis", DirectSolver(A_gpu, structure, view), nothing, - nothing) - continue - end + A_gpu = api_matrix(backend, structure == "G" ? A_cpu : triangle_view(A_cpu, view), INT) ported_solver(backend, T, INT, A_gpu, structure, view, n) end + if structure == "G" + @testset "\"G\" ignores the view" begin + b_cpu = rand(T, n) + xs = map(('L', 'U', 'F')) do view + solver = DirectSolver(api_matrix(backend, A_cpu, INT), "G", view) + x_gpu = to_device(backend, zeros(T, n)) + b_gpu = to_device(backend, b_cpu) + for phase in ("analysis", "factorization", "solve") + execute!(phase, solver, x_gpu, b_gpu) + end + @test getparam(solver, "info") == 0 + Array(x_gpu) + end + @test xs[1] == xs[3] + @test xs[2] == xs[3] + @test relres(A_cpu, xs[3], b_cpu) <= tol(T) + end + end end if T <: Complex solver = DirectSolver(api_matrix(backend, A_cpu, INT), "SPD", 'F') diff --git a/test/test_api.jl b/test/test_api.jl index b7388d1..99e39f4 100644 --- a/test/test_api.jl +++ b/test/test_api.jl @@ -29,6 +29,7 @@ end @test solver isa AbstractDirectSolver{T, INT} @test solver isa LinearAlgebra.Factorization{T} @test size(solver) == (60, 60) && size(solver, 1) == 60 + @test size(solver, 3) == 1 && thrown(() -> size(solver, 0)) isa BoundsError @test solver.fresh_factorization @test solver.backend == backend @test getparam(solver, "info") == 0 @@ -177,8 +178,11 @@ end setparam!(solver, "solve_alg", "algo1") @test thrown(() -> execute!("solve", solver, x, b)) isa NotSupportedError setparam!(solver, "solve_alg", "default") - s3 = DirectSolver(api_matrix(backend, tril(A)), "G", 'L') # LU (T19) reads the full matrix: view 'F' only - @test thrown(() -> execute!("analysis", s3, x, b)) isa InvalidValueError + s3 = DirectSolver(api_matrix(backend, A), "G", 'L') # "G" ignores the view (cuDSS, #84) + for phase in ("analysis", "factorization", "solve") + execute!(phase, s3, x, b) + end + @test getparam(s3, "info") == 0 && relres(A, to_host(x), to_host(b)) <= tol(T) if T <: Complex s4 = DirectSolver(api_matrix(backend, tril(A)), "SPD", 'L') @test thrown(() -> execute!("analysis", s4, x, b)) isa InvalidValueError @@ -206,6 +210,20 @@ end est = getparam(solver, "memory_estimates") @test est isa Vector{Int64} && length(est) == 16 && est == SDS.memory_estimates(S, T, INT) @test est[1] > 0 && est[2] >= est[1] + # with matching, its device buffers count in slots 1, 2 and 5 (#97): scaled values, weights, two scale vectors, + # and for "G" the column permutation + for structure in ("G", "S") + ms = DirectSolver(api_matrix(backend, structure == "G" ? A : tril(A), INT), structure, + structure == "G" ? 'F' : 'L') + setparam!(ms, "matching_alg", "algo5") + execute!("analysis", ms, nothing, nothing) + nzm = nnz(structure == "G" ? A : tril(A)) + extra = nzm * sizeof(T) + (nzm + 2 * n) * sizeof(real(T)) + (structure == "G" ? n * sizeof(INT) : 0) + mref = SDS.memory_estimates(ms.host_symbolic, T, INT) + mm = getparam(ms, "memory_estimates") + @test mm[[1, 2, 5]] == mref[[1, 2, 5]] .+ extra + @test mm[[3, 4, 6:16...]] == mref[[3, 4, 6:16...]] + end # diag: the diagonal of the Cholesky factor of A[perm, perm] dref = diag(cholesky(Hermitian(Matrix(A[perm, perm]), :L)).L) d = getparam(solver, "diag") diff --git a/test/test_matrix.jl b/test/test_matrix.jl index 82c9f2f..7aad711 100644 --- a/test/test_matrix.jl +++ b/test/test_matrix.jl @@ -43,6 +43,11 @@ end @test SparseMatrixCSC(B) == transpose(A) # plain transpose, also for complex T @test logical_matrix(B) == A @test occursin("transposed", sprint(show, B)) + # spare capacity in rowval/nzval (allowed by SparseMatrixCSC) is refused, not read as entries (#31) + for Ag in (SparseMatrixCSC(size(A)..., copy(A.colptr), [A.rowval; INT(1)], copy(A.nzval)), + SparseMatrixCSC(size(A)..., copy(A.colptr), copy(A.rowval), [A.nzval; zero(T)])) + @test thrown(() -> csr_of_transpose(Ag)) isa InvalidValueError + end end @testset "nbatch ($T)" for T in ELTYPES diff --git a/test/test_numeric_ldlt.jl b/test/test_numeric_ldlt.jl index 40bd1b4..c0fe2c1 100644 --- a/test/test_numeric_ldlt.jl +++ b/test/test_numeric_ldlt.jl @@ -82,6 +82,41 @@ end end end +@testset "pivot_epsilon = 0: a zero pivot fails ($(backend_name(backend)), $T)" for backend in BACKENDS, T in ELTYPES + Random.seed!(666) + n, j = 60, 23 + F = SDS.FRONT_STATS_FIELDS + # an exactly zero pivot is tiny whatever ε; with ε = 0 it stays zero and info is its original column (#65) + for (rname, rkw) in (("default", (;)), ("regime-C root", (factorization_alg = "algo2", subtree_budgets = Int[]))), + stored_zero in (false, true) + A = singular_block_matrix(T, n, j; stored_zero) + opts = Options(; pivot_epsilon = 0.0, rkw...) + S, Nr, info_ref, C = reference_ldlt(A; opts) + @test info_ref == j + Sd = SDS.adapt(backend, S, Int32) + Nd = SDS.allocate_numeric(Sd, T, backend) + @test SDS.factorize!(Nd, Sd, to_device(backend, C.nzval); opts) == j + Nh = SDS.host_numeric(Nd) + @test Nh.pivot_kind == Nr.pivot_kind && Nh.piv == Nr.piv + @test d_error(Nh, Nr) <= panel_tol(T) + k = findfirst(==(SDS.PIVOT_KIND_PERTURBED), Nh.pivot_kind) + @test k !== nothing && iszero(Nh.d[k]) && iszero(Nr.d[k]) + @test to_host(Nd.stats)[F:F:end] == Nr.stats[F:F:end] + @test count(!iszero, Nr.stats[F:F:end]) == 1 + st = SDS.pivot_totals(Nd) + @test st.nzero == 1 && st.nperturbed == 1 + end + # the public API reports it in "info"; ε > 0 again completes + solver = DirectSolver(api_matrix(backend, tril(singular_block_matrix(T, n, j))), sym_structure(T), 'L') + setparam!(solver, "pivot_epsilon", 0.0) + execute!("analysis", solver, nothing, nothing) + execute!("factorization", solver, nothing, nothing) + @test getparam(solver, "info") == j + setparam!(solver, "pivot_epsilon", nothing) + execute!("refactorization", solver, nothing, nothing) + @test getparam(solver, "info") == 0 +end + @testset "perturbation and pivot_sign ($(backend_name(backend)), $T)" for backend in BACKENDS, T in ELTYPES n, j = 60, 23 R = real(T) diff --git a/test/test_numeric_lu.jl b/test/test_numeric_lu.jl index 73ea9e5..f894236 100644 --- a/test/test_numeric_lu.jl +++ b/test/test_numeric_lu.jl @@ -93,6 +93,37 @@ end @test lu_error(A, S, Nh) <= tol(T) end +@testset "pivot_epsilon = 0: a zero pivot fails ($(backend_name(backend)), $T)" for backend in BACKENDS, T in ELTYPES + Random.seed!(666) + n, j = 60, 23 + F = SDS.FRONT_STATS_FIELDS + # an exactly zero pivot is tiny whatever ε; with ε = 0 it stays zero and info is its original column (#65) + for (rname, rkw) in (("default", (;)), ("regime-C root", (factorization_alg = "algo2", subtree_budgets = Int[]))), + stored_zero in (false, true) + A = singular_block_matrix(T, n, j; stored_zero) + opts = Options(; pivot_epsilon = 0.0, rkw...) + C = SDS.CSR(A) + S = SDS.symbolic_analysis(C, "G", 'F'; opts) + Nr = SDS.allocate_numeric(S, T) + @test SDS.ref_lu!(Nr, S, C.nzval; opts) == j + Sd = SDS.adapt(backend, S, Int32) + Nd = SDS.allocate_numeric(Sd, T, backend) + @test SDS.factorize!(Nd, Sd, to_device(backend, C.nzval); opts) == j + Nh = SDS.host_numeric(Nd) + @test Nh.pivot_kind == Nr.pivot_kind && Nh.piv == Nr.piv + @test d_error(Nh, Nr) <= panel_tol(T) + @test to_host(Nd.stats)[F:F:end] == Nr.stats[F:F:end] + @test count(!iszero, Nr.stats[F:F:end]) == 1 + st = SDS.pivot_totals(Nd) + @test st.nzero == 1 && st.nperturbed == 1 + end + solver = DirectSolver(api_matrix(backend, singular_block_matrix(T, n, j)), "G", 'F') + setparam!(solver, "pivot_epsilon", 0.0) + execute!("analysis", solver, nothing, nothing) + execute!("factorization", solver, nothing, nothing) + @test getparam(solver, "info") == j +end + @testset "perturbation ($(backend_name(backend)), $T)" for backend in BACKENDS, T in ELTYPES n, j = 60, 23 R = real(T) @@ -327,9 +358,15 @@ end execute!("solve", cs, x, b) @test relres(M, to_host(x), to_host(b)) <= tol(T) end - # views other than 'F' are refused for "G" - @test thrown(() -> execute!("analysis", DirectSolver(api_matrix(backend, A, INT), "G", 'L'), nothing, nothing)) isa - InvalidValueError + # "G" ignores the view (cuDSS, #84): 'L' and 'U' read the full matrix, as 'F' + for view in ('L', 'U') + vs = DirectSolver(api_matrix(backend, A, INT), "G", view) + setparam!(vs, "pivot_threshold", 1.0) + for phase in ("analysis", "factorization", "solve") + execute!(phase, vs, x, b) + end + @test getparam(vs, "info") == 0 && relres(A, to_host(x), to_host(b)) <= tol(T) + end # the LinearAlgebra layer F = lu(SDS.CSR(SparseMatrixCSC{T, INT}(A))) @test F.structure == SDS.STRUCTURE_GENERAL && getparam(F, "ir_n_steps") == 2 diff --git a/test/test_symbolic_etree.jl b/test/test_symbolic_etree.jl index 7afa9bc..40e35c2 100644 --- a/test/test_symbolic_etree.jl +++ b/test/test_symbolic_etree.jl @@ -181,6 +181,12 @@ end for index in ('O', 'Z') PG = SDS.SymmetricPattern(SDS.CSR(G; index), "G") @test SparseMatrixCSC(PG) == offdiag_pattern(Gp + Gp') + # "G" ignores the view (cuDSS, #84): every stored entry is read + @test SDS.SymmetricPattern(SDS.CSR(G; index), "G"; view = 'L') == PG + @test SDS.SymmetricPattern(SDS.CSR(G; index), "G"; view = 'U') == PG + FL = SDS.full_pattern_map(SDS.CSR(G; index), "G"; view = 'L') + FF = SDS.full_pattern_map(SDS.CSR(G; index), "G"; view = 'F') + @test all(getfield(FL, f) == getfield(FF, f) for f in fieldnames(SDS.FullPatternMap)) end end @@ -197,7 +203,6 @@ end R = sprand(4, 5, 0.5) @test thrown(() -> SDS.SymmetricPattern(SDS.CSR(R), "G")) isa InvalidValueError @test thrown(() -> SDS.full_pattern_map(SDS.CSR(R), "G")) isa InvalidValueError - @test thrown(() -> SDS.SymmetricPattern(SDS.CSR(laplacian2d(3, 3)), "G"; view = 'L')) isa InvalidValueError @test thrown(() -> SDS.SymmetricPattern([1, 2, 3], [1, 4], 2, "S")) isa InvalidValueError @test thrown(() -> SDS.SymmetricPattern([1, 2, 3], [1, 0], 2, "S")) isa InvalidValueError @test thrown(() -> SDS.SymmetricPattern([1, 2], [1, 1], 2, "S")) isa InvalidValueError From be091df2c17c578f3682df08ff90c1c1c3b7cc73 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 7 Oct 2026 21:25:13 +0000 Subject: [PATCH 2/3] Fix: csr_of_transpose test builds the spare capacity after construction Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01AvwMxoGmJm3eSmpzxAy72A --- test/test_matrix.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/test/test_matrix.jl b/test/test_matrix.jl index 7aad711..699a192 100644 --- a/test/test_matrix.jl +++ b/test/test_matrix.jl @@ -44,10 +44,12 @@ end @test logical_matrix(B) == A @test occursin("transposed", sprint(show, B)) # spare capacity in rowval/nzval (allowed by SparseMatrixCSC) is refused, not read as entries (#31) - for Ag in (SparseMatrixCSC(size(A)..., copy(A.colptr), [A.rowval; INT(1)], copy(A.nzval)), - SparseMatrixCSC(size(A)..., copy(A.colptr), copy(A.rowval), [A.nzval; zero(T)])) - @test thrown(() -> csr_of_transpose(Ag)) isa InvalidValueError - end + # (the constructor checks the lengths, so the spare entries are pushed afterwards) + Ar, Av = copy(A), copy(A) + push!(Ar.rowval, INT(1)) + push!(Av.nzval, zero(T)) + @test thrown(() -> csr_of_transpose(Ar)) isa InvalidValueError + @test thrown(() -> csr_of_transpose(Av)) isa InvalidValueError end @testset "nbatch ($T)" for T in ELTYPES From 2ef4a26753281e806edb7a7dedc19c76b20a8878 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 7 Oct 2026 22:03:52 +0000 Subject: [PATCH 3/3] Fix: deterministic solve for the bitwise "G" view check; #86 owner note format Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01AvwMxoGmJm3eSmpzxAy72A --- TASKS.md | 6 ++++-- test/ported/cudss_solver.jl | 3 +++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/TASKS.md b/TASKS.md index 4eac47a..82ece82 100644 --- a/TASKS.md +++ b/TASKS.md @@ -2718,8 +2718,10 @@ for no gain; the matching-based pairs of #67 must do better there. Decide in 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. -#86: compare the device and the reference on the K2 dumps after scaling (inertia, nperturbed, piv); close #86 if -scaling removed the divergence, else as documented. + +**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 diff --git a/test/ported/cudss_solver.jl b/test/ported/cudss_solver.jl index caf5b61..594afb7 100644 --- a/test/ported/cudss_solver.jl +++ b/test/ported/cudss_solver.jl @@ -195,6 +195,9 @@ end b_cpu = rand(T, n) xs = map(('L', 'U', 'F')) do view solver = DirectSolver(api_matrix(backend, A_cpu, INT), "G", view) + # the atomic forward sweep (default) adds in launch order on GPUs; bitwise + # equality needs the deterministic solve (#36) + setparam!(solver, "deterministic_mode", 1) x_gpu = to_device(backend, zeros(T, n)) b_gpu = to_device(backend, b_cpu) for phase in ("analysis", "factorization", "solve")