Route RFLU matrix right-hand sides through TriangularSolve (1.5-1.7x multi-RHS solve)#1117
Merged
Merged
Conversation
RecursiveFactorization already uses TriangularSolve inside its own lu!, but the ldiv! that consumes the factorization only does so for the pivotless NotIPIV case (RecursiveFactorization/src/lu.jl); a pivoted LU - which is what RFLUFactorization produces - falls back to LinearAlgebra, i.e. BLAS trsv/trsm. Fixing that inside RecursiveFactorization would be type piracy (it owns neither ldiv! nor LU), so it belongs here: LinearSolve owns solve! for its own algorithm type. Matrix right-hand sides now go through TriangularSolve's blocked kernels. Measured solve-only, factorization already computed, 1 BLAS thread: n nrhs=4 nrhs=8 64 6.31 -> 3.76 us 8.76 -> 5.40 us (1.68x / 1.62x) 128 21.11 -> 12.68 us 26.20 -> 17.22 us (1.66x / 1.52x) 256 86.43 -> 55.11 us 97.09 -> 67.87 us (1.57x / 1.43x) 500 277.11 -> 175.42 us 301.13 -> 225.77 us (1.58x / 1.33x) Vector right-hand sides keep the stdlib path: TriangularSolve has no vector kernel, and reshaping to n x 1 measured 1.09x at n=128 but 0.88x by n=256. A single-column matrix is the same case in disguise, so it is guarded explicitly - with the guard, nrhs=1 and vector results are bit-identical to the stdlib path (agreement exactly 0.0), while multi-RHS keeps 1.5-1.7x. TriangularSolve is added as a weakdep purely to trigger the extension; RecursiveFactorization depends on it, so loading RecursiveFactorization alone still activates everything and no dependency graph changes. New multi-RHS RFLU testset covers the direct and cache-reuse paths. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
ChrisRackauckas
marked this pull request as ready for review
July 26, 2026 10:00
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Independent of the SupernodalLU PRs. Second TriangularSolve PR — the first is #1113 (SupernodalLU panel trsms); this one is the dense RFLU solve path.
Why here and not in RecursiveFactorization
RecursiveFactorization already uses TriangularSolve inside its own
lu!, but theldiv!that consumes the factorization only does so for the pivotlessNotIPIVcase (RecursiveFactorization/src/lu.jl:54). A pivotedLU— which is exactly whatRFLUFactorizationproduces — falls back to LinearAlgebra, i.e. BLAStrsv/trsm.Fixing it inside RecursiveFactorization would be type piracy: it owns neither
ldiv!norLU, so such a method would change behavior for all Julia code holding a standardLU, whether or not it came from RF. LinearSolve ownssolve!for its own algorithm type, so this is the legitimate home.Measured (solve-only, factorization already computed, 1 BLAS thread)
What deliberately does not change
Vector right-hand sides keep the stdlib path — TriangularSolve has no vector kernel, and reshaping to n×1 measured 1.09× at n=128 but 0.88× by n=256. A single-column matrix is the same case in disguise and takes the matrix method, so it is guarded explicitly. With the guard,
nrhs=1and vector results are bit-identical to the stdlib path (agreement exactly0.0), while multi-RHS keeps 1.5–1.7×.TriangularSolve is added as a weakdep purely to trigger the extension; RecursiveFactorization already depends on it, so
using RecursiveFactorizationalone still activates everything and nobody's dependency graph changes.GROUP=Corepasses, with a new multi-RHS RFLU testset covering the direct and cache-reuse paths.🤖 Generated with Claude Code
https://claude.ai/code/session_017rr42T1vFRRT8D7DMAmJzf