SciML
diff --git a/‎.typos.toml‎
Lines changed: 5 additions & 0 deletions b/‎.typos.toml‎
Lines changed: 5 additions & 0 deletions
diff --git a/‎docs/src/solvers/solvers.md‎
Lines changed: 10 additions & 0 deletions b/‎docs/src/solvers/solvers.md‎
Lines changed: 10 additions & 0 deletions
diff --git a/‎ext/LinearSolveSparseArraysExt.jl‎
Lines changed: 127 additions & 1 deletion b/‎ext/LinearSolveSparseArraysExt.jl‎
Lines changed: 127 additions & 1 deletion
diff --git a/‎src/LinearSolve.jl‎
Lines changed: 2 additions & 0 deletions b/‎src/LinearSolve.jl‎
Lines changed: 2 additions & 0 deletions
diff --git a/‎src/SupernodalLU/NOTICE.md‎
Lines changed: 61 additions & 0 deletions b/‎src/SupernodalLU/NOTICE.md‎
Lines changed: 61 additions & 0 deletions
diff --git a/‎src/SupernodalLU/SupernodalLU.jl‎
Lines changed: 44 additions & 0 deletions b/‎src/SupernodalLU/SupernodalLU.jl‎
Lines changed: 44 additions & 0 deletions
@@ -76,3 +76,8 @@ MTK = "MTK"
 ODE = "ODE"
 PDE = "PDE"
 SDE = "SDE"
+
+# SupernodalLU (src/SupernodalLU) identifiers
+pn = "pn"          # SuiteSparse AMD port variable (amd.jl)
+Tpos = "Tpos"      # transpose-position map in the supernodal factorization
+hashi = "hashi"    # SuiteSparse AMD port variable (amd.jl)
@@ -230,6 +230,7 @@ LinearSolveSparseArraysExt.KLU.klu
 LinearSolveSparseArraysExt.KLU.klu!
 PureKLUFactorization
 PureUMFPACKFactorization
+SupernodalLUFactorization
 UMFPACKFactorization
 SparseColumnPivotedQRFactorization
 ```
@@ -241,6 +242,15 @@ SparseColumnPivotedQRFactorization
     patterns, replacing the SuiteSparse-backed `KLUFactorization` in the default
     polyalgorithm.
 
+!!! note
+    
+    `SupernodalLUFactorization` is a pure-Julia implementation of the supernodal
+    left–right-looking sparse LU method of Schenk & Gärtner (vendored
+    self-contained in `src/SupernodalLU`, no binary dependency). It is the
+    strongest choice for "more structured" (PDE-mesh-like) sparse systems,
+    where it outperforms both `UMFPACKFactorization` and `KLUFactorization`.
+    It is not yet wired into the default polyalgorithm; request it explicitly.
+
 !!! note
 
     `SparseColumnPivotedQRFactorization` is a pure-Julia, rank-revealing
 
@@ -7,7 +7,7 @@ using LinearSolve: LinearSolve, BLASELTYPES, pattern_changed, ArrayInterface,
     NormalCholeskyFactorization,
     OperatorAssumptions, LinearVerbosity,
     QRFactorization, RFLUFactorization, UMFPACKFactorization,
-    SparseColumnPivotedQRFactorization, solve
+    SparseColumnPivotedQRFactorization, SupernodalLUFactorization, solve
 using SciMLOperators: AbstractSciMLOperator, has_concretization
 using ArrayInterface: ArrayInterface
 using LinearAlgebra: LinearAlgebra, I, Hermitian, Symmetric, cholesky, ldiv!, lu, lu!
@@ -31,6 +31,10 @@ include("../src/KLU/klu.jl")
 # PureKLU (pure-Julia, no SuiteSparse) is a hard dependency and the default
 # sparse LU; the SuiteSparse `KLUFactorization` above is unchanged.
 import PureKLU
+# SupernodalLU (pure-Julia supernodal left-right-looking LU, Schenk-Gärtner
+# method) is vendored in src/SupernodalLU: the BLAS-3 sparse LU for
+# structured systems.
+const SNLU = LinearSolve.SupernodalLU
 # SparseColumnPivotedQR (pure-Julia, rank-revealing column-pivoted sparse QR) is a
 # hard dependency: the default sparse QR and the singular-LU fallback.
 import SparseColumnPivotedQR
@@ -680,6 +684,128 @@ function SciMLBase.solve!(
     end
 end
 
+# --- SupernodalLU: pure-Julia supernodal left-right-looking LU (Schenk-Gärtner) ---
+# The BLAS-3 sparse LU for structured (PDE-mesh-like) systems (vendored in src/SupernodalLU).
+
+function LinearSolve.init_cacheval(
+        alg::SupernodalLUFactorization, A::AbstractArray, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    return nothing
+end
+
+function LinearSolve.init_cacheval(
+        alg::SupernodalLUFactorization, A::LinearSolve.GPUArraysCore.AnyGPUArray, b, u,
+        Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    return nothing
+end
+
+const PREALLOCATED_SUPERNODAL = SNLU.snlu(
+    SparseMatrixCSC{Float64, Int64}(0, 0, [Int64(1)], Int64[], Float64[])
+)
+
+function LinearSolve.init_cacheval(
+        alg::SupernodalLUFactorization, A::AbstractSparseArray{Float64, Int64}, b, u,
+        Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    return PREALLOCATED_SUPERNODAL
+end
+
+# SupernodalLU is pure Julia and factors any `Number` element type; the empty
+# cacheval carries the correct types so `pplu!`/`pplu` dispatch is concrete.
+function LinearSolve.init_cacheval(
+        alg::SupernodalLUFactorization, A::AbstractSparseArray{T, Ti}, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    ) where {T <: Number, Ti <: Integer}
+    return SNLU.snlu(
+        SparseMatrixCSC{T, Ti}(0, 0, [one(Ti)], Ti[], T[])
+    )
+end
+
+function LinearSolve.init_cacheval(
+        alg::SupernodalLUFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    if has_concretization(A)
+        return LinearSolve.init_cacheval(
+            alg, convert(AbstractMatrix, A), b, u, Pl, Pr,
+            maxiters, abstol, reltol, verbose, assumptions
+        )
+    else
+        nothing
+    end
+end
+
+function LinearSolve.pattern_changed(
+        F::SNLU.SupernodalLUFactor, A::SparseArrays.AbstractSparseMatrixCSC
+    )
+    Aold = F.A
+    return getcolptr(Aold) != getcolptr(A) || rowvals(Aold) != rowvals(A)
+end
+
+function SciMLBase.solve!(
+        cache::LinearSolve.LinearCache, alg::SupernodalLUFactorization; kwargs...
+    )
+    A = cache.A
+    A = LinearSolve.reduce_operand!(cache.sparse_reduction, A)
+    A = convert(AbstractMatrix, A)
+    if cache.isfresh
+        cacheval = LinearSolve.@get_cacheval(cache, :SupernodalLUFactorization)
+        As = SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))
+        if alg.reuse_symbolic && size(cacheval) == size(As) &&
+                nnz(cacheval.A) == nnz(As) &&
+                !(alg.check_pattern && pattern_changed(cacheval, As))
+            # numeric-only refactorization: reuses the analysis, matching, and
+            # all numeric storage (allocation-free)
+            fact = SNLU.snlu!(cacheval, As)
+        else
+            # `check = false`: static pivoting never aborts — numerically
+            # singular systems surface through the finiteness check below.
+            fact = SNLU.snlu(
+                As; ordering = alg.ordering, matching = alg.matching,
+                eps_pivot = alg.eps_pivot, threaded = alg.threaded,
+                dense_alg = alg.dense_alg, check = false
+            )
+        end
+        cache.cacheval = fact
+        cache.isfresh = false
+    end
+    F = LinearSolve.@get_cacheval(cache, :SupernodalLUFactorization)
+    y = SNLU.solve!(cache.u, F, cache.b)
+    # Static pivoting never aborts: a numerically singular system factors with
+    # perturbed pivots and produces a finite but meaningless solution. When
+    # pivots were perturbed (rare), verify the residual (one sparse mat-vec)
+    # so singularity surfaces as `Infeasible` instead of a silent `Success`.
+    ok = all(isfinite, y)
+    if ok && SNLU.nperturbed(F) > 0
+        r = F.A * y
+        r .-= cache.b
+        bn = LinearAlgebra.norm(cache.b)
+        ok = LinearAlgebra.norm(r) <= 1.0e-6 * max(bn, floatmin(real(eltype(r))))
+    end
+    return if ok
+        SciMLBase.build_linear_solution(
+            alg, y, nothing, nothing; retcode = ReturnCode.Success
+        )
+    else
+        @SciMLMessage(
+            "Solver produced a non-finite or inaccurate solution; matrix is likely singular",
+            cache.verbose, :solver_failure
+        )
+        SciMLBase.build_linear_solution(
+            alg, cache.u, nothing, nothing; retcode = ReturnCode.Infeasible
+        )
+    end
+end
+
 # --- SparseColumnPivotedQR: pure-Julia rank-revealing column-pivoted sparse QR ---
 # The default sparse QR (non-square sparse systems) and the singular-LU fallback.
 
 
@@ -409,6 +409,7 @@ function defaultalg_symbol end
 
 include("verbosity.jl")
 include("blas_logging.jl")
+include("SupernodalLU/SupernodalLU.jl")
 include("generic_lufact.jl")
 include("eigenvalue.jl")
 include("common.jl")
@@ -524,6 +525,7 @@ export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
     RFLUFactorization, ButterflyFactorization,
     NormalCholeskyFactorization, NormalBunchKaufmanFactorization,
     UMFPACKFactorization, KLUFactorization, PureKLUFactorization,
+    SupernodalLUFactorization,
     PureUMFPACKFactorization, SparseColumnPivotedQRFactorization, FastLUFactorization,
     FastQRFactorization,
     SparspakFactorization, DiagonalFactorization, CholeskyFactorization,
 
@@ -0,0 +1,61 @@
+SupernodalLU (LinearSolve src/SupernodalLU) — code lineage and licensing notice
+==================================================
+
+Summary: MIT for this folder, BSD-3-Clause for the vendored `amd.jl`.
+No GPL/LGPL-licensed code and no proprietary code is included.  The full
+license texts are in LICENSE.  Per-component provenance:
+
+1. The PARDISO method itself (supernodal left-looking LU on the symmetrized
+   pattern, restricted-block pivoting with static perturbation, weighted
+   matching preprocessing) is implemented from the published literature:
+
+     - O. Schenk, K. Gärtner: "Solving unsymmetric sparse systems of linear
+       equations with PARDISO", FGCS 20(3), 2004.
+     - O. Schenk, K. Gärtner: "On fast factorization pivoting methods for
+       sparse symmetric indefinite systems", ETNA 23, 2006.
+
+   PARDISO (pardiso-project.org) and Intel MKL PARDISO are closed-source;
+   no source code of either exists publicly and none was used.
+
+2. `amd.jl` is a pure-Julia port of SuiteSparse AMD, vendored from
+   PureKLU.jl.  SuiteSparse AMD is BSD-3-Clause, (c) 1996-2025 Timothy A.
+   Davis, Patrick R. Amestoy, and Iain S. Duff; the file remains under
+   BSD-3-Clause (text in LICENSE).
+
+3. `src/matching.jl` implements maximum-weight (max-product) bipartite
+   matching with dual-variable scalings from the published algorithm
+   descriptions (I. S. Duff, J. Koster, SIMAX 22(4), 2001; M. Olschowka,
+   A. Neumaier, Linear Algebra Appl. 240, 1996), as shortest augmenting
+   paths with potentials (Jonker–Volgenant style).  HSL MC64 is proprietary
+   and its source was not used.
+
+4. `src/symbolic.jl` implements standard published symbolic-analysis
+   algorithms — the elimination tree via union-find with path halving
+   (J. W. H. Liu, "The role of elimination trees in sparse factorization",
+   SIMAX 11(1), 1990), L-structure prediction by row subtrees (J. R.
+   Gilbert, E. Ng, B. W. Peyton, SIMAX 15(1), 1994), and relaxed supernode
+   amalgamation in the style of Ashcraft–Grimes with the relaxation
+   constants published for CHOLMOD (Y. Chen, T. A. Davis, W. W. Hager,
+   S. Rajamanickam, ACM TOMS 35(3), 2008).  The implementations were
+   written from the algorithm descriptions and deliberately use different
+   formulations from the CSparse/SuiteSparse (LGPL/GPL) implementations of
+   the same algorithms: union-find etree instead of ancestor stamping,
+   cursor-stack DFS over counting-sorted child arrays instead of
+   destructive sibling lists, row-subtree column construction instead of
+   child-set merging, and single bucket-pass pattern permutation instead
+   of per-column sorting.
+
+5. The banded-detection acceptance thresholds in `symbolic.jl` follow the
+   heuristic introduced in PureKLU.jl by the same copyright holder.
+
+6. Everything else (`numeric.jl`, `solve.jl`, `interface.jl`,
+   `ordering.jl` apart from the AMD calls, and the
+   RecursiveFactorization/TriangularSolve kernel overrides in
+   LinearSolveRecursiveFactorizationExt) is original work of the same
+   authors, MIT-licensed.
+
+For comparison, the related packages carry different licenses because they
+ARE ports: PureKLU.jl is LGPL-2.1+ (KLU/BTF translation) and PureUMFPACK.jl
+is GPL-2.0+ (UMFPACK/CSparse translation).  This folder is not a
+translation of any existing codebase, which is why MIT (+ BSD-3 for the
+vendored `amd.jl`) is the appropriate license, matching LinearSolve's.
@@ -0,0 +1,44 @@
+# SPDX-FileCopyrightText: 2026 Chris Rackauckas <accounts@chrisrackauckas.com> and contributors
+# SPDX-License-Identifier: MIT
+#
+# SupernodalLU — a pure-Julia implementation of the supernodal
+# left-right-looking sparse LU method of O. Schenk and K. Gärtner
+# (FGCS 20(3), 2004; ETNA 23, 2006): supernodal BLAS-3 LU on the symmetrized
+# nonzero pattern of A + Aᵀ, pivoting restricted to the supernode diagonal
+# block with static pivot perturbation (ε‖A‖) compensated by iterative
+# refinement, and maximum-weight matching + scaling preprocessing.
+# Implemented from the papers; see NOTICE.md for the full per-component
+# lineage (no code from the proprietary PARDISO library, which popularized
+# the method, nor from Intel MKL PARDISO or HSL; `amd.jl` is a BSD-3-Clause
+# SuiteSparse AMD port).
+#
+# Internal to LinearSolve: the public surface is `SupernodalLUFactorization`.
+# Entry points here are `snlu` (analyze + factor), `snlu!` (refactorize, same
+# pattern), `solve!`, and the `snlu_symbolic` analysis.
+# Each supernode's dense diagonal block is factored through LinearSolve's own
+# dense `init`/`solve!` machinery (a per-supernode dense `LinearCache` over a
+# view of the panel), so the sparse solver automatically uses whatever
+# LinearSolve's dense default resolves to — RFLUFactorization when
+# RecursiveFactorization is loaded, MKL/LAPACK/generic otherwise — with no
+# backend-specific wiring here.  Static pivoting is preserved by an
+# optimistic fast path: accept the dense factorization iff every pivot
+# clears the ε‖A‖ threshold, else restore the block and rerun the built-in
+# static-perturbation kernel.
+
+module SupernodalLU
+
+using SparseArrays: SparseMatrixCSC, sparse, getcolptr, rowvals, nonzeros, nnz
+using LinearAlgebra: LinearAlgebra, SingularException, UpperTriangular,
+    UnitLowerTriangular, ldiv!, rdiv!, mul!, norm, BLAS
+using ..LinearSolve: LinearSolve
+using SciMLBase: SciMLBase
+
+include("amd.jl")          # vendored BSD-3 SuiteSparse AMD port (module AMD)
+include("ordering.jl")
+include("symbolic.jl")
+include("matching.jl")
+include("numeric.jl")
+include("solve.jl")
+include("interface.jl")
+
+end # module