SupernodalLU: route panel triangular solves through TriangularSolve (#1113)

The supernodal factorization applies two BLAS-3 trsms per supernode
against its just-factored diagonal block: the L21 panel on the right by
U11, and the U12 panel on the left by unit-L11. They now go through
overridable hooks whose default is the stdlib triangular solve, with
LinearSolveRecursiveFactorizationExt overriding them with
TriangularSolve.

TriangularSolve is a hard dependency of RecursiveFactorization and is
what RecursiveFactorization uses for its own trsms, so this adds nothing
to anyone's dependency graph: it is declared as a weakdep purely to
trigger the extension, and since RecursiveFactorization pulls it in,
loading RecursiveFactorization alone still activates everything. That is
also exactly the condition under which RFLUFactorization becomes the
dense default, so the sparse solver's panel work and the dense default
run on the same kernels.

Measured refactorization, BLAS pinned to 1 thread:

  poisson2d_512  0.4781 -> 0.3921 s  (-18%)
  poisson3d_32   0.3466 -> 0.3027 s  (-13%)

Residuals unchanged (1.5e-14 / 8.8e-16).

GROUP=Core passes.

Co-authored-by: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

Author: 3 people

ext/LinearSolveRecursiveFactorizationExt.jl

@@ -224,4 +224,31 @@ function LinearSolve._custom_adjoint_factorization_solve(
     return solution[1:n]
 end
 
+# ---- SupernodalLU panel triangular solves ---------------------------------
+# The vendored supernodal sparse LU (src/SupernodalLU) applies two BLAS-3
+# trsms per supernode against its just-factored diagonal block: the L21 panel
+# on the right by U11, and the U12 panel on the left by unit-L11.  Route them
+# through TriangularSolve, which RecursiveFactorization already depends on
+# and uses for its own trsms — so when RFLU is the dense default, the sparse
+# solver's panel work runs on the same kernels.  Measured: recovers the
+# 2D-mesh refactorization gap left by the stdlib trsms.
+const SNLU = LinearSolve.SupernodalLU
+const SNLUTypes = Union{Float32, Float64}
+
+function SNLU._panel_rdiv!(W::Matrix{Tv}, np::Int, len::Int) where {Tv <: SNLUTypes}
+    len > np || return nothing
+    TriangularSolve.rdiv!(
+        view(W, (np + 1):len, 1:np), UpperTriangular(view(W, 1:np, 1:np)), Val(false)
+    )
+    return nothing
+end
+
+function SNLU._panel_ldiv!(W::Matrix{Tv}, np::Int, Z::Matrix{Tv}) where {Tv <: SNLUTypes}
+    isempty(Z) && return nothing
+    TriangularSolve.ldiv!(
+        UnitLowerTriangular(view(W, 1:np, 1:np)), Z, Val(false)
+    )
+    return nothing
+end
+
 end

src/SupernodalLU/numeric.jl

@@ -274,6 +274,23 @@ function _cache_lu!(
     return true
 end
 
+# Panel triangular solves (L21 := L21·U11⁻¹ and U12 := L11⁻¹·U12).  These are
+# BLAS-3 trsms against the just-factored diagonal block, not linear solves in
+# the LinearSolve sense, so they are plain overridable hooks: the defaults
+# below use the stdlib triangular solves, and
+# LinearSolveRecursiveFactorizationExt routes them through TriangularSolve —
+# the same library RecursiveFactorization (and hence the RFLU dense default)
+# already uses internally for its own trsms.
+function _panel_rdiv!(W::Matrix{Tv}, np::Int, len::Int) where {Tv}
+    rdiv!(view(W, (np + 1):len, 1:np), UpperTriangular(view(W, 1:np, 1:np)))
+    return nothing
+end
+
+function _panel_ldiv!(W::Matrix{Tv}, np::Int, Z::Matrix{Tv}) where {Tv}
+    ldiv!(UnitLowerTriangular(view(W, 1:np, 1:np)), Z)
+    return nothing
+end
+
 _lu_from_cacheval(cv::LinearAlgebra.LU) = cv
 _lu_from_cacheval(cv::Tuple) = _lu_from_cacheval(first(cv))
 
@@ -497,8 +514,8 @@ function _process_supernode!(
         end
 
         if nu > 0
-            rdiv!(view(Ws, (np + 1):len, 1:np), UpperTriangular(view(Ws, 1:np, 1:np)))
-            ldiv!(UnitLowerTriangular(view(Ws, 1:np, 1:np)), Zs)
+            _panel_rdiv!(Ws, np, len)
+            _panel_ldiv!(Ws, np, Zs)
         end
 
         for a in 1:np