Batched (matrix) right-hand sides: LinearProblem(A, B) ≡ A \ B (v4.0.0) (#1072)

* Support batched (matrix) right-hand sides in factorization algorithms

`solve(LinearProblem(A, B))` with `B::AbstractMatrix` now computes the
equivalent of `A \ B`: `u0` is initialized as a `size(A, 2) × size(B, 2)`
matrix and the factorization-based algorithms solve all columns against a
single factorization of `A`. This is a breaking change (v4.0.0): previously a
matrix `b` initialized a flattened vector `u` and errored downstream.

- `__init_u0_from_Ab` gains matrix-`b` methods (including SMatrix
  disambiguation).
- Iterative/Krylov algorithms (`KrylovJL`, `SimpleGMRES`, and the other
  `AbstractKrylovSubspaceMethod`s, plus the default algorithm when it selects
  a Krylov method for an operator) throw an informative `ArgumentError` at
  `init` time for matrix `b`; `SimpleLUFactorization` (vector-only workspace)
  does the same.
- The default polyalgorithm's unused Krylov cacheval slots initialize to
  `nothing` for matrix `b` so structured-matrix defaults (Diagonal,
  SymTridiagonal, ...) work.
- `_check_residual_safety` sizes its residual buffer shape-generically.
- Default-algorithm size heuristics use `size(b, 1)` instead of `length(b)`
  so batched right-hand sides don't inflate the apparent problem size (same
  for the sparse `use_klulike_sparse_structure` heuristic).
- MKL/OpenBLAS/AppleAccelerate LU copy solution columns correctly for matrix
  `b` in the rectangular (`m > n`) path.
- SparseArrays extension: CHOLMOD pre-1.12 `_ldiv!` fallback widened to
  matrices, SPQR gains a matrix `_ldiv!` via `\` (no in-place matrix ldiv! on
  any Julia version), and SparseColumnPivotedQR solves matrix `b`
  column-by-column against the one factorization.

Part of #552

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* Add batched right-hand side tests

Covers dense LU/GenericLU/QR (NoPivot + ColumnNorm)/SVD/Cholesky/
BunchKaufman/NormalCholesky vs A \ B for Float64 and ComplexF64, the default
algorithm (dense, sparse, structured, non-square with k != n), sparse
UMFPACK/KLU/CHOLMOD and the non-square sparse QR default, cache reuse via
`cache.b = B2` and `cache.A = A2`, failure retcodes on singular A (including
the default's QR fallback), residual-safety with matrix B, informative errors
for Krylov methods and SimpleLUFactorization, static arrays, BigFloat, and
single-column-matrix vs vector consistency.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* Document batched right-hand sides

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh

* Widen in-repo LinearSolve compat for v4 and source docs env locally

All sublibrary, test, qa, GPU, Trim, and docs environments pinned
LinearSolve = "3[.x.y]", so every CI job died at resolution against the
4.0.0 bump ("empty intersection between LinearSolve@4.0.0 and project
compatibility"). Widen each to include 4. The docs environment also gains
[sources] entries for LinearSolve and LinearSolveAutotune (the registered
Autotune restricts LinearSolve to 3, so the docs build must use the local
copies, matching the pattern the test environments already use).

Verified locally: docs, lib/LinearSolveAutotune, and test/qa environments
all instantiate against 4.0.0.

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

* Runic-format forwarddiff_overloads.jl

The repo-wide Runic Format Check fails on main because this file landed
unformatted in aabf007 (#1069); format it here alongside the dispatch fix
for the same commit's fallout.

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

* Allow AlgebraicMultigrid 0.5 so the test target resolves against v4

Every registered AlgebraicMultigrid >= 0.6 caps LinearSolve at 3, so the
root test target (which all CI test groups resolve) was unsatisfiable
against the 4.0.0 bump. AMG 0.5.x predates its LinearSolve dependency and
exposes the identical API surface the extension uses
(SmoothedAggregationAMG, RugeStubenAMG, aspreconditioner, ruge_stuben);
the AlgebraicMultigridJL testset passes on 0.5.1 (verified locally,
including the tight-tolerance case). The resolver will pick AMG 2.x again
automatically once AlgebraicMultigrid widens its own compat; the compat
comment marks the entry for removal then.

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

* Retrigger CI (self-hosted runner disconnect on Downgrade Core)

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

* Drop LinearSolve v3 from in-repo lower bounds

All in-repo environments and both registered sublibraries now require
LinearSolve 4: the environments always develop the local copy, and
v4-era releases of LinearSolveAutotune/LinearSolvePyAMG should require
the batched-RHS LinearSolve rather than straddle the major boundary.
Minor-bump both sublibraries (Autotune 1.12.0, PyAMG 1.2.0) for the
dependency floor raise. All environment classes re-verified to resolve.

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

---------

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

Author: 3 people

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearSolve"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "3.87.0"
+version = "4.0.0"
 authors = ["SciML"]
 
 [deps]
@@ -110,7 +110,10 @@ LinearSolveSpecializingFactorizationsExt = "SpecializingFactorizations"
 [compat]
 AMD = "0.5"
 AMDGPU = "1.2, 2"
-AlgebraicMultigrid = "1, 2"
+# 0.5 allowed so the test target resolves against LinearSolve 4 until
+# AlgebraicMultigrid widens its own LinearSolve compat (its 1.x/2.x cap at 3);
+# the ext API is identical and the AMG testset passes on 0.5.1.
+AlgebraicMultigrid = "0.5, 1, 2"
 ArrayInterface = "7.19"
 BandedMatrices = "1.8"
 BlockDiagonals = "0.2"

docs/Project.toml

@@ -6,6 +6,10 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 
 [compat]
 Documenter = "1"
-LinearSolve = "3"
+LinearSolve = "4"
 LinearSolveAutotune = "1.1"
 SciMLOperators = "1"
+
+[sources]
+LinearSolve = {path = ".."}
+LinearSolveAutotune = {path = "../lib/LinearSolveAutotune"}

docs/src/release_notes.md

@@ -1,5 +1,14 @@
 # Release Notes
 
+## v4.0
+
+  - Batched (matrix) right-hand sides are now supported: `solve(LinearProblem(A, B))` with
+    `B::AbstractMatrix` computes the equivalent of `A \ B`, factorizing `A` once and
+    returning `sol.u` as a `size(A, 2) × size(B, 2)` matrix. This is a breaking change:
+    previously a matrix `b` initialized a vector-shaped `u` and generally errored downstream.
+    Batched right-hand sides are supported by the factorization-based algorithms; iterative
+    (Krylov) methods throw an informative `ArgumentError` for matrix `b`.
+
 ## Upcoming Changes
 
   - `CudaOffloadFactorization` has been split into two algorithms:

docs/src/tutorials/linear.md

@@ -36,6 +36,26 @@ pass in an algorithm struct and all wrapped linear solvers are immediately
 available as tweaks to the general algorithm. For more information on the
 available solvers, see [the solvers page](@ref linearsystemsolvers)
 
+## Batched Right-Hand Sides
+
+A matrix right-hand side `B` solves all of its columns against the same `A` at
+once, just like `A \ B`: `solve(LinearProblem(A, B))` factorizes `A` a single
+time and returns `sol.u` as the `size(A, 2) × size(B, 2)` matrix satisfying
+`A * sol.u ≈ B`.
+
+```@example linsys1
+B = rand(4, 2)
+prob = LS.LinearProblem(A, B)
+sol = LS.solve(prob)
+sol.u ≈ A \ B
+```
+
+This works with the caching interface as well (`cache.b = B2; solve!(cache)`
+re-solves all columns against the cached factorization). Batched right-hand
+sides are only supported by factorization-based algorithms; iterative (Krylov)
+methods such as `LS.KrylovJL_GMRES()` accept only vector `b` and throw an
+informative error for matrix `B` — solve column-by-column instead in that case.
+
 ## Sparse and Structured Matrices
 
 There is no difference in the interface for LinearSolve.jl on sparse

ext/LinearSolveSparseArraysExt.jl

@@ -778,12 +778,24 @@ function SciMLBase.solve!(
         cache.isfresh = false
     end
     F = LinearSolve.@get_cacheval(cache, :SparseColumnPivotedQRFactorization)
-    y = ldiv!(cache.u, F, cache.b)
+    y = LinearSolve._ldiv!(cache.u, F, cache.b)
     return SciMLBase.build_linear_solution(
         alg, y, nothing, cache; retcode = ReturnCode.Success
     )
 end
 
+# SparseColumnPivotedQR's ldiv! only accepts vector right-hand sides; batched
+# (matrix) right-hand sides solve column-by-column against the one factorization.
+function LinearSolve._ldiv!(
+        x::AbstractMatrix,
+        F::SCPQR.SparseColumnPivotedQRFactorization, b::AbstractMatrix
+    )
+    for j in axes(b, 2)
+        ldiv!(view(x, :, j), F, view(b, :, j))
+    end
+    return x
+end
+
 # Build a column-pivoted sparse QR factorization for the default sparse-LU
 # singular fallback (`_do_sparse_qr_fallback` in src/default.jl).
 function LinearSolve.sparse_colpivqr_factorize(A)
@@ -874,8 +886,8 @@ end
             x .= A \ b
         end
         function LinearSolve._ldiv!(
-                x::AbstractVector,
-                A::SparseArrays.CHOLMOD.Factor, b::AbstractVector
+                x::AbstractVecOrMat,
+                A::SparseArrays.CHOLMOD.Factor, b::AbstractVecOrMat
             )
             x .= A \ b
         end
@@ -897,6 +909,15 @@ end
         end
     end
 
+    # SPQR has no in-place matrix (batched) ldiv! on any current Julia version,
+    # so route batched right-hand sides through the allocating `\`.
+    function LinearSolve._ldiv!(
+            x::AbstractMatrix,
+            A::SparseArrays.SPQR.QRSparse, b::AbstractMatrix
+        )
+        x .= A \ b
+    end
+
     function LinearSolve._ldiv!(
             ::SVector,
             A::Union{SparseArrays.CHOLMOD.Factor, SparseArrays.SPQR.QRSparse},
@@ -940,11 +961,11 @@ end
 # medium and very sparse / "less structure") favors the pure-Julia KLU-style
 # solvers: PureKLU for LU and SparseColumnPivotedQR for QR. `false` ("more
 # structure") favors the SuiteSparse solvers: UMFPACK for LU and SPQR for QR.
-# The `length(b) <= 1_000` branch is the "small enough" fast path; the second is
+# The `size(b, 1) <= 1_000` branch is the "small enough" fast path; the second is
 # the "medium and sufficiently sparse" rule.
 function LinearSolve.use_klulike_sparse_structure(A::AbstractSparseMatrixCSC, b)
-    return length(b) <= 1_000 ||
-        (length(b) <= 10_000 && length(nonzeros(A)) / length(A) < 2.0e-4)
+    return size(b, 1) <= 1_000 ||
+        (size(b, 1) <= 10_000 && length(nonzeros(A)) / length(A) < 2.0e-4)
 end
 
 @static if Base.USE_GPL_LIBS

lib/LinearSolveAutotune/Project.toml

@@ -1,7 +1,7 @@
 name = "LinearSolveAutotune"
 uuid = "67398393-80e8-4254-b7e4-1b9a36a3c5b6"
 authors = ["SciML"]
-version = "1.11.1"
+version = "1.12.0"
 
 [deps]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
@@ -44,7 +44,7 @@ FastLapackInterface = "2.0.4"
 GitHub = "5"
 LAPACK_jll = "3.12"
 LinearAlgebra = "1"
-LinearSolve = "3.39.2"
+LinearSolve = "4"
 MKL_jll = "2025.2.0"
 Metal = "1.5"
 OpenBLAS_jll = "0.3"

lib/LinearSolveAutotune/test/qa/Project.toml

@@ -13,7 +13,7 @@ LinearSolveAutotune = {path = "../.."}
 [compat]
 Aqua = "0.8"
 JET = "0.9, 0.10, 0.11"
-LinearSolve = "3"
+LinearSolve = "4"
 LinearSolveAutotune = "1"
 SciMLTesting = "1.6"
 Test = "1"

lib/LinearSolvePyAMG/Project.toml

@@ -1,7 +1,7 @@
 name = "LinearSolvePyAMG"
 uuid = "7a56c47d-7ab1-4e99-b0e3-2952e463d64a"
 authors = ["SciML"]
-version = "1.1.1"
+version = "1.2.0"
 
 [deps]
 CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
@@ -17,7 +17,7 @@ LinearSolve = {path = "../.."}
 [compat]
 CondaPkg = "0.2"
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "4"
 Pkg = "1"
 PythonCall = "0.9"
 SafeTestsets = "0.1"

lib/LinearSolvePyAMG/test/qa/Project.toml

@@ -13,7 +13,7 @@ LinearSolvePyAMG = {path = "../.."}
 [compat]
 Aqua = "0.8"
 JET = "0.9, 0.10, 0.11"
-LinearSolve = "3"
+LinearSolve = "4"
 LinearSolvePyAMG = "1"
 SciMLTesting = "1.6"
 Test = "1"

src/appleaccelerate.jl

@@ -382,7 +382,11 @@ function SciMLBase.solve!(
     if m > n
         Bc = copy(cache.b)
         aa_getrs!('N', A.factors, A.ipiv, Bc; info)
-        copyto!(cache.u, 1, Bc, 1, n)
+        if cache.b isa AbstractMatrix
+            copyto!(cache.u, @view(Bc[1:n, :]))
+        else
+            copyto!(cache.u, 1, Bc, 1, n)
+        end
     else
         copyto!(cache.u, cache.b)
         aa_getrs!('N', A.factors, A.ipiv, cache.u; info)
@@ -468,7 +472,11 @@ function SciMLBase.solve!(
     if m > n
         aa_getrs!('N', A_lu.factors, A_lu.ipiv, b_32; info)
         # Convert back to original precision
-        cache.u[1:n] .= Torig.(b_32[1:n])
+        if cache.b isa AbstractMatrix
+            cache.u .= Torig.(@view(b_32[1:n, :]))
+        else
+            cache.u[1:n] .= Torig.(@view(b_32[1:n]))
+        end
     else
         copyto!(u_32, b_32)
         aa_getrs!('N', A_lu.factors, A_lu.ipiv, u_32; info)