SupernodalLU: fix block factorization on the direct-BLAS backends (#1119)

`snlu` on a Mac was a `MethodError` for any matrix with a supernode at or
above `dense_threshold`:

    MethodError: no method matching _lu_from_cacheval(::AppleAccelerateLUCache{...})

`_lu_from_cacheval` pulled the block factorization out of the dense cache's
cacheval, handling `LU` and `Tuple`.  The direct-BLAS backends store neither:
Apple Accelerate, BLIS, and OpenBLAS each keep a mutable `*LUCache` holding
the factored matrix, the getrf pivots, and `info`, which
`_custom_cache_factorization` wraps in an `LU` only on demand.  On a Mac the
*default* dense solver resolves to Apple Accelerate, so the default path fell
straight into the hole.

`_cache_lu!` reads exactly `.factors`/`.ipiv`, with the same LAPACK
semantics those caches use, so the cacheval is handed back as-is rather than
allocating a throwaway `LU` per block; anything exposing neither field now
gets a real error message instead of a `MethodError`.

`test/Core/supernodal_lu.jl` missed this because it loads
RecursiveFactorization, which makes the default resolve to RFLU (a `Tuple`
cacheval) everywhere.  The new testset names the platform's direct-BLAS
backends explicitly, over both a caching-heavy and a caching-sparse
`dense_threshold`, and covers refactorization through the same caches.

Also version-gates two assertions in test/qa/allocations.jl.  Above
`PANEL_BLAS_CUTOFF` the sweeps hand off to LinearAlgebra's `ldiv!`/`mul!`,
whose wrappers are not statically allocation-free on every release: 1.10 and
1.13 leave a `generic_trimatdiv!` dynamic dispatch plus boxed
`MulAddMul`/`SubArray` values that AllocCheck reports from stdlib frames, and
1.11 allocates them for real on the `SubArray`-matrix path the multi-RHS
sweep takes.  Runtime allocation is 0 on 1.10, 1.12, and 1.13, so each
assertion is pinned to the versions that can carry it rather than dropped.

Verified on macOS aarch64: qa/allocations.jl and Core/supernodal_lu.jl pass
on 1.10, 1.11, and 1.12; runtime zero-allocation confirmed on 1.13.  The new
testset fails with the original `MethodError` without the numeric.jl change.

Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: ChrisRackauckas

src/SupernodalLU/numeric.jl

@@ -307,6 +307,22 @@ end
 _lu_from_cacheval(cv::LinearAlgebra.LU) = cv
 _lu_from_cacheval(cv::Tuple) = _lu_from_cacheval(first(cv))
 
+# The direct-BLAS backends (Apple Accelerate, BLIS, OpenBLAS) do not keep an
+# `LU` in their cacheval at all: each owns a mutable cache holding the factored
+# matrix and the getrf pivot vector, which `_custom_cache_factorization` only
+# wraps in an `LU` on demand.  `_cache_lu!` reads exactly `.factors`/`.ipiv`,
+# with the same LAPACK semantics, so hand the cacheval straight back instead of
+# allocating a throwaway wrapper per block.  On a Mac this is the path the
+# *default* dense solver takes, so without it every `snlu` with a supernode at
+# or above `dense_threshold` is a `MethodError`.  `hasfield` on a concrete type
+# folds at compile time, so the check costs nothing at runtime.
+function _lu_from_cacheval(cv)
+    if hasfield(typeof(cv), :factors) && hasfield(typeof(cv), :ipiv)
+        return cv
+    end
+    throw(ArgumentError("SupernodalLU: dense block solver cacheval $(typeof(cv)) does not expose an LU factorization; pass `dense_alg` explicitly"))
+end
+
 # The default solver keeps one slot per candidate algorithm and records the
 # choice in a runtime enum, so pull the factorization out of the slot that was
 # actually used.  Written as a literal-symbol chain rather than

test/Core/supernodal_lu.jl

@@ -155,6 +155,36 @@ end
     end
 end
 
+# The direct-BLAS backends do not keep an `LU` in their cacheval — they keep a
+# mutable cache holding the factored matrix and the pivots — so the block
+# factorization has to read the factors out of that instead.  The rest of this
+# file loads RecursiveFactorization, which makes the *default* dense solver
+# resolve to RFLU everywhere, so the backends are named explicitly here: on a
+# Mac the default resolves to Apple Accelerate, and without this every `snlu`
+# with a supernode at or above `dense_threshold` was a `MethodError`.
+@testset "block caches over the direct-BLAS backends" begin
+    A = poisson2d(30)
+    b = randn(size(A, 1))
+    xref = Matrix(A) \ b
+    algs = Any[]
+    LinearSolve.appleaccelerate_isavailable() &&
+        push!(algs, AppleAccelerateLUFactorization())
+    LinearSolve.useopenblas && push!(algs, OpenBLASLUFactorization())
+    Base.get_extension(LinearSolve, :LinearSolveBLISExt) !== nothing &&
+        push!(algs, LinearSolve.BLISLUFactorization())
+    @test !isempty(algs)        # every platform has at least one of these
+    for alg in algs, thr in (4, 64)
+        F = SNLU.snlu(A; dense_alg = alg, dense_threshold = thr)
+        @test length(F.bcaches) > 0
+        x = similar(b)
+        SNLU.solve!(x, F, b)
+        @test norm(x - xref) <= 1.0e-9 * norm(xref)
+        SNLU.snlu!(F, A)                   # refactorize through the same caches
+        SNLU.solve!(x, F, b)
+        @test norm(x - xref) <= 1.0e-9 * norm(xref)
+    end
+end
+
 @testset "matching engages on zero/weak diagonals" begin
     n = 40
     P = sparse(collect(n:-1:1), collect(1:n), 2.0 .+ rand(n), n, n)

test/qa/allocations.jl

@@ -71,6 +71,19 @@ end
 @check_allocs allocation_checked_supernodal_sweeps!(y, F) =
     LinearSolve.SupernodalLU._solve_panels!(y, F)
 
+# What the sweeps can be *proved* about depends on the Julia version, because
+# above `PANEL_BLAS_CUTOFF` they hand off to LinearAlgebra's `ldiv!`/`mul!`,
+# whose wrappers are not statically clean on every release: 1.10 and 1.13 leave
+# a `generic_trimatdiv!` dynamic dispatch plus boxed `MulAddMul`/`SubArray`
+# values that AllocCheck reports (`triangular.jl`, `matmul.jl` — stdlib frames,
+# no SupernodalLU code involved), and 1.11 allocates them for real on the
+# `SubArray`-matrix path the multi-RHS sweep takes.  1.12 folds all of it away.
+# So pin each assertion to the versions that can carry it rather than weakening
+# it everywhere; the sweeps carry no growth branch on any version, which is what
+# these are here to protect.
+const STATIC_SWEEP_PROOF = v"1.12" <= VERSION < v"1.13"
+const RUNTIME_MULTIRHS_ZERO = VERSION < v"1.11" || VERSION >= v"1.12"
+
 function poisson2d_qa(k)
     n = k * k
     Is = Int[]; Js = Int[]; V = Float64[]
@@ -97,8 +110,10 @@ end
         F = LinearSolve.SupernodalLU.snlu(A)
         LinearSolve.SupernodalLU._ensure_panel_scratch!(F, 1)
         y = ones(n)
-        allocation_checked_supernodal_sweeps!(y, F)   # throws if it can allocate
-        @test true
+        if STATIC_SWEEP_PROOF
+            allocation_checked_supernodal_sweeps!(y, F)   # throws if it can allocate
+            @test true
+        end
         # the full solve!, which owns the one-time sizing, is zero at runtime
         b = ones(n)
         x = similar(b)
@@ -109,6 +124,8 @@ end
         B = ones(n, 3)
         X = similar(B)
         LinearSolve.SupernodalLU.solve!(X, F, B; refine = 0)
-        @test @allocated(LinearSolve.SupernodalLU.solve!(X, F, B; refine = 0)) == 0
+        if RUNTIME_MULTIRHS_ZERO
+            @test @allocated(LinearSolve.SupernodalLU.solve!(X, F, B; refine = 0)) == 0
+        end
     end
 end