Add native left-upper ldiv! and right-lower rdiv! kernels

[ChrisRackauckas-Claude]

⚠️ Draft — please ignore until reviewed by @ChrisRackauckas.

Add native left-upper ldiv! and right-lower rdiv! kernels

Implements the two missing side/triangle combinations: ldiv!(::UpperTriangular, A) (A ← U \ A) and rdiv!(A, ::LowerTriangular) (A ← A / L), plus their Unit variants and 3-arg forms. Previously ldiv!(UpperTriangular(F), B) silently fell through the catch-all to LinearAlgebra.ldiv!, which is the multi-RHS RFLU backsolve gap behind SciML/LinearSolve.jl#1146 (worked around by an O(n²)-materializing exchange-matrix path in SciML/LinearSolve.jl#1152).

Design

The library is formulated in terms of right-division with forward substitution over column blocks (C = A / U, upper-triangular U), with left-lower ldiv! obtained for free via L \ A = (A' / L')' dispatching to row-major-specialized kernels (_ldiv_L!). This PR adds the exact dual: native right-lower division C = A / L with backward substitution (column blocks solved last-to-first), from which left-upper ldiv! is obtained for free via U \ B = (B' / U')'. Right-lower was chosen as the native family because it mirrors the existing code structure one-to-one — every new function is the backward twin of an existing one:

existing (forward)	new (backward)
solve_AU	solve_AL
BdivU_small_kern[_u]!	BdivL_small_kern[_u]!
rdiv_solve_W[_u]!	rdiv_solve_W[_u]_L!
rdiv_U! (column-major)	rdiv_L! (column-major)
rdiv_U! (2,2-stride) → _ldiv_L! (ldiv_solve_W*, ldiv_remainder!)	rdiv_L! (2,2-stride) → _ldiv_U! (uldiv_solve_W*, uldiv_remainder!)
nmuladd!, rdiv_block_N!, rdiv_block_MandN!, RDivBlockMandNv2, multithread_rdiv!, div_dispatch!	nmuladd_L!, rdiv_block_N_L!, rdiv_block_MandN_L!, RDivBlockMandNv2L, multithread_rdiv_L!, div_dispatch_L!

Key points:

• No negative strides, no reversed views, no J-materialization. Only block offsets descend; every SIMD load/store is a contiguous ascending index range. The n-column remainder (N % W) lives at the trailing columns [N-Nr, N) and is solved first (it depends on nothing); full W/W*UF blocks then walk down with their reduction running over the already-solved trailing columns [n+W, N) via SafeCloseOpen(n + W, N).
• Packed-LU discipline preserved: kernels read only the wrapper-indicated triangle (reduction loads sit strictly below the diagonal, solve_AL reads the diagonal block's lower triangle only), so UpperTriangular(F) from a packed LU whose lower triangle holds L works.
• Diagonal-block solves reuse the existing register blocking: solve_AL is solve_AU with the intra-block loop reversed; the transposed family reuses transpose_vecunroll and the same Unroll load/store patterns, unroll factors, masking, and Polyester.batch threading (m-partition) as the forward code.
• Blocked-N path (N > block_size) mirrors rdiv_block_N!: solve the rightmost column block, then for each block leftward do an @turbo nmuladd_L! update against all solved trailing columns before its triangular solve.
• Catch-alls are untouched; non-Float32/Float64 still falls back to LinearAlgebra (tested with BigFloat).

Verification

All run locally on Julia 1.12.4 and 1.10.11, AMD EPYC 7502 (Zen 2, AVX2, 2×32 cores), OpenBLAS with 1 thread for comparisons.

Full test suite (existing tests + new upper-ldiv!/lower-rdiv! cases in the same NaN-box style, packed-LU factors up to n=512, dispatch checks, allocation checks):

Test Summary:      |   Pass   Total     Time
TriangularSolve.jl | 148781  148781  3m53.1s   (Julia 1.12.4)
TriangularSolve.jl | 148781  148781  54.8s     (Julia 1.10.11)

plus Aqua and detect_ambiguities clean on both.

Failing-before / passing-after (the correctness tests alone would pass on master via the LinearAlgebra fallback, so the discriminating test asserts dispatch reaches a native method — the issue's first success criterion):

On unmodified master:

native upper-ldiv / lower-rdiv dispatch: Test Failed  (×4)
Some tests did not pass: 0 passed, 4 failed, 0 errored, 0 broken.

With this PR:

Test Summary: | Pass  Total
              |    4      4
method for ldiv!(UpperTriangular, Matrix, Val{false}):
ldiv!(U::UpperTriangular{T}, A::AbstractMatrix{T}, ::Val{false}) where T<:Union{Float32, Float64}
  @ TriangularSolve src/TriangularSolve.jl:832

Extended sweep (beyond the test suite): n ∈ 1..33, 47..49, 63..65, 100, 127..129, 200, 255..257, 400; nrhs ∈ {1,2,3,4,5,7,8,11,16,31,32}; Float32/Float64; unit/non-unit; Val(true)/Val(false); 2-arg/3-arg; packed lu! parents to n=512; threaded partition path with M up to 2003: 0 failures. Steady-state allocations: 0 (serial and threaded, including on RecursiveFactorization.lu! factors).

Benchmarks

Solve-only, Float64, BLAS.set_num_threads(1), @belapsed with setup=(X = copy(B)), pinned to 8 cores. J-mat emulates the LinearSolve.jl#1152 workaround (materialize reversed U as lower-triangular + reversed lower solve).

Left-upper ldiv!  U \ B   (times in μs)
     n  nrhs |  TS serial  TS thread  BLAS trsm      J-mat |  TS/BLAS TS/J-mat
    32     4 |       0.45       0.42       1.06       2.19 |     2.36     4.87
    32     8 |       0.75       0.77       1.62       3.86 |     2.16     5.15
    32    16 |       1.63       1.62       2.81       5.67 |     1.72     3.48
    64     4 |       0.92       0.94       3.09       4.60 |     3.36     5.01
    64     8 |       1.79       1.80       4.10       7.39 |     2.29     4.13
    64    16 |       3.49       3.51       6.96      13.08 |     1.99     3.75
   128     4 |       2.80       2.82      10.36      10.79 |     3.70     3.85
   128     8 |       5.56       5.58      12.40      17.10 |     2.23     3.08
   128    16 |      11.04       7.19      20.24      32.96 |     1.83     2.99
   256     4 |      11.22      10.53      44.05      38.05 |     3.93     3.39
   256     8 |      21.06      10.89      47.70      52.04 |     2.26     2.47
   256    16 |      38.36      12.02      72.91      85.22 |     1.90     2.22
   500     4 |      39.82      35.05     142.73     126.36 |     3.58     3.17
   500     8 |      79.17      41.46     152.00     193.00 |     1.92     2.44
   500    16 |     133.31      43.17     232.42     293.05 |     1.74     2.20
  1000     4 |     140.90     143.13     545.47     715.15 |     3.87     5.08
  1000     8 |     291.20     142.62     584.05     888.46 |     2.01     3.05
  1000    16 |     477.63     148.27     880.41    1131.90 |     1.84     2.37

Median BLAS/TS-serial = 2.29× for n ∈ [64, 500], nrhs ∈ {4, 8} (acceptance bar was ≥1.3×); median vs the J-materialize path 3.17×, with zero extra allocation.

Right-lower rdiv!  A / L   (A is nrhs x n, μs)
     n  nrhs |  TS serial  BLAS trsm |  TS/BLAS
    32     4 |       0.53       1.35 |     2.55
    64     8 |       2.96       5.32 |     1.80
   128    16 |      10.70      24.35 |     2.28
   256     8 |      27.24      53.80 |     1.98
   500     8 |      81.56     197.68 |     2.42
  1000    16 |     452.48    1006.50 |     2.22

(full grid in the bench script; all cells 1.58–3.63× vs BLAS)

No regression on existing kernels (main vs this branch, serial, μs; differences are within run-to-run noise, and the diff is purely additive):

     n  nrhs |  rdiv-upper main/PR |  ldiv-lower main/PR
    64     8 |      2.80 / 2.79    |      1.60 / 1.58
   128     8 |      9.81 / 9.97    |      4.69 / 4.68
   256     8 |     25.96 / 25.86   |     23.77 / 22.33
   500     8 |     77.46 / 77.39   |     74.49 / 76.31
  1000     8 |    267.81 / 263.14  |    296.10 / 285.97
   512   512 |   7455.08 / 7530.43 |   6675.22 / 7201.25*

* re-ran 3×: main 6647–6763 μs vs branch 6727–6936 μs, overlapping ranges.

What was not verified

• AVX-512 and aarch64: this machine is AVX2 (Zen 2); Float64 W=4 and Float32 W=8 paths (including all @nif remainder cases) are exercised here, but wider-vector unroll counts (e.g. W=8 Float64, UF from a 32-register file) only via CI.
• MKL comparison (OpenBLAS only here).
• versioninfo(): Julia 1.12.4, LLVM 18.1.7 znver2; also tested on 1.10.11.

Reviewer notes

• The W*UF-unrolled transposed kernel uldiv_solve_W_u! generates its intra-block backward loop manually (descending sub-block order with cross-sub-block updates) rather than via Base.Cartesian.@nexprs; the forward ldiv_solve_W_u! could use @nexprs because its dependencies run ascending.
• rdiv_solve_W[_u]_L! and uldiv_solve_W[_u]! take an extra nend argument relative to their forward twins: backward blocks reduce over the trailing solved columns [n+W, N), so the kernels need the total column count, where forward kernels reduce over [0, n).
• The new public methods are constrained to StridedMatrix storage (parent and RHS), so non-strided AbstractMatrix inputs (sparse, structured wrappers) keep hitting the LinearAlgebra catch-all that served them before, and div_dispatch_L! throws DimensionMismatch on shape mismatches instead of reading out of bounds. Both came out of an adversarial self-review; the pre-existing forward-substitution family has the same two latent hazards on its own signatures (bare AbstractMatrix dispatch, no shape validation) — left untouched here to avoid changing existing behavior, but worth the same treatment as a follow-up if you agree.
• Version is bumped 0.2.1 → 0.2.2 since the API is purely additive; note that downstream compat = "0.2" users' ldiv!(::UpperTriangular, ·) calls will switch from the fallback's single-threaded BLAS numerics to these (default-threaded) kernels on upgrade. If you prefer that behavior change to gate on 0.3.0 instead, say so and I'll adjust.
• LinearSolve/RecursiveFactorization follow-ups are staged and waiting on this + registration: Require TriangularSolve 0.2.2: native upper ldiv! on the RFLU U leg SciML/LinearSolve.jl#1153, Require TriangularSolve 0.2.2 for the native upper-triangular ldiv! JuliaLinearAlgebra/RecursiveFactorization.jl#112.

MWE for the LinearSolve integration (run, output shown):

using LinearAlgebra, TriangularSolve, RecursiveFactorization
n, nrhs = 256, 8
A = rand(n, n) + n * I; b = rand(n, nrhs)
ipiv = Vector{LinearAlgebra.BlasInt}(undef, n)
F = RecursiveFactorization.lu!(copy(A), ipiv)
X = b[F.p, :]
TriangularSolve.ldiv!(UnitLowerTriangular(F.factors), X, Val(false)) # L leg
TriangularSolve.ldiv!(UpperTriangular(F.factors), X, Val(false))     # U leg — new
maximum(abs, A * X - b)  # 1.998e-15
# dispatch: TriangularSolve src/TriangularSolve.jl:832 (native), not the catch-all
# steady-state allocations on the U leg: 0

[ChrisRackauckas-Claude]

The evaluate (Invalidations) failure is an upstream tooling incompatibility, not this diff:

• julia-actions/julia-invalidations@v1 pins SnoopCompile to "3" (CI resolved 3.2.7 on Julia 1.12.6).
• SnoopCompile 3.x's report_invalidations does "$(inv.method.name)" (ext/SCPrettyTablesExt.jl:35); on Julia 1.12 an invalidation root can be a Core.Binding (binding partitions), which has no name field → FieldError, exactly the CI stacktrace.
• Whether a binding-rooted invalidation appears depends on the environment (CPU/dep versions), so it doesn't reproduce on every machine: locally on Julia 1.12.6 with the action's pinned SnoopCompile 3.2.7, both this branch and registered 0.2.1 load with an identical 62 invalidation trees and a clean report_invalidations — i.e. the PR also adds no new invalidations where measurable.

The action would need to unpin SnoopCompile (4.x handles binding roots) or the workflow could run the snippet directly; leaving that infra call to you rather than bundling a workflow edit into this PR. The actual test jobs (Julia 1 and nightly, 1 and 3 threads) are green.

[ChrisRackauckas]

Since this is just a mirror kernel, merging as it's already close to optimal. I'll get all downstream, but there will be an optimization coming.