From 17a127e470a123a8b9d52613819fdd2cdbe51a41 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Fri, 3 Jul 2026 11:42:52 -0400
Subject: [PATCH 1/8] Support batched (matrix) right-hand sides in
 factorization algorithms
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`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 https://github.com/SciML/LinearSolve.jl/issues/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
---
 Project.toml                      |  2 +-
 ext/LinearSolveSparseArraysExt.jl | 33 +++++++++++++++----
 src/appleaccelerate.jl            | 12 +++++--
 src/common.jl                     | 54 +++++++++++++++++++++++++++++++
 src/default.jl                    | 10 +++---
 src/factorization.jl              |  6 ++--
 src/iterative_wrappers.jl         | 13 ++++++++
 src/mkl.jl                        | 12 +++++--
 src/openblas.jl                   | 12 +++++--
 src/simplelu.jl                   |  6 ++++
 10 files changed, 141 insertions(+), 19 deletions(-)

diff --git a/Project.toml b/Project.toml
index b1af7bf35..54c830227 100644
--- a/Project.toml
+++ b/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]
diff --git a/ext/LinearSolveSparseArraysExt.jl b/ext/LinearSolveSparseArraysExt.jl
index 923739705..379f9e52f 100644
--- a/ext/LinearSolveSparseArraysExt.jl
+++ b/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
diff --git a/src/appleaccelerate.jl b/src/appleaccelerate.jl
index ae86c5483..c533ede85 100644
--- a/src/appleaccelerate.jl
+++ b/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)
diff --git a/src/common.jl b/src/common.jl
index 4bf3af799..fe5084147 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -405,6 +405,9 @@ same element type as `b` and sized to match the number of columns in `A`.
 
 ## Returns
 A zero-initialized vector of size `(size(A, 2),)` with element type matching `b`.
+For a matrix (batched) right-hand side `b` of size `(size(A, 1), k)`, returns a
+zero-initialized matrix of size `(size(A, 2), k)` so that each column of `u0`
+corresponds to a column of `b`.
 
 ## Specializations
 - For static matrices (`SMatrix`): Returns a static vector (`SVector`)
@@ -415,7 +418,56 @@ function __init_u0_from_Ab(A, b)
     fill!(u0, false)
     return u0
 end
+function __init_u0_from_Ab(A, b::AbstractMatrix)
+    u0 = similar(b, size(A, 2), size(b, 2))
+    fill!(u0, false)
+    return u0
+end
 __init_u0_from_Ab(::SMatrix{S1, S2}, b) where {S1, S2} = zeros(SVector{S2, eltype(b)})
+function __init_u0_from_Ab(::SMatrix{S1, S2}, b::AbstractMatrix) where {S1, S2}
+    u0 = similar(b, S2, size(b, 2))
+    fill!(u0, false)
+    return u0
+end
+function __init_u0_from_Ab(
+        ::SMatrix{S1, S2}, ::SMatrix{S1b, S2b, Tb}
+    ) where {S1, S2, S1b, S2b, Tb}
+    return zeros(SMatrix{S2, S2b, Tb})
+end
+
+"""
+    _check_batched_rhs_support(alg, b)
+
+Throw an informative `ArgumentError` at `init` time when a matrix (batched)
+right-hand side `b` is used with an algorithm that only supports vector `b`
+(Krylov subspace / iterative methods). Factorization-based algorithms support
+matrix `b` and pass through the generic no-op fallback.
+"""
+_check_batched_rhs_support(alg, b) = nothing
+function _check_batched_rhs_support(alg::AbstractKrylovSubspaceMethod, b::AbstractMatrix)
+    throw(
+        ArgumentError(
+            "Batched (matrix) right-hand sides are only supported by factorization " *
+                "algorithms; $(nameof(typeof(alg))) supports only vector `b`. Solve " *
+                "column-by-column or use a factorization algorithm (e.g. `LUFactorization()`)."
+        )
+    )
+end
+function _check_batched_rhs_support(alg::DefaultLinearSolver, b::AbstractMatrix)
+    if alg.alg === DefaultAlgorithmChoice.KrylovJL_GMRES ||
+            alg.alg === DefaultAlgorithmChoice.KrylovJL_CRAIGMR ||
+            alg.alg === DefaultAlgorithmChoice.KrylovJL_LSMR
+        throw(
+            ArgumentError(
+                "Batched (matrix) right-hand sides are only supported by factorization " *
+                    "algorithms; the default algorithm selected the Krylov method " *
+                    "$(alg.alg) for this operator, which supports only vector `b`. " *
+                    "Solve column-by-column or use a factorization algorithm."
+            )
+        )
+    end
+    return nothing
+end
 
 function SciMLBase.init(prob::LinearProblem, alg::SciMLLinearSolveAlgorithm, args...; kwargs...)
     return __init(prob, alg, args...; kwargs...)
@@ -506,6 +558,8 @@ function __init(
         copy(b)
     end
 
+    _check_batched_rhs_support(alg, b)
+
     u0_ = u0 !== nothing ? u0 : __init_u0_from_Ab(A, b)
 
     # Guard against type mismatch for user-specified reltol/abstol
diff --git a/src/default.jl b/src/default.jl
index 92025b2a3..69b27f7f3 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -354,11 +354,13 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool})
                 )
 
                 # Small matrix override - always use GenericLUFactorization for tiny problems
-                if length(b) <= 10
+                # `size(b, 1)` (not `length(b)`) so batched (matrix) right-hand
+                # sides don't inflate the apparent problem size.
+                if size(b, 1) <= 10
                     DefaultAlgorithmChoice.GenericLUFactorization
                 else
                     # Check if autotune preferences exist for larger matrices
-                    matrix_size = length(b)
+                    matrix_size = size(b, 1)
                     eltype_A = A === nothing ? Nothing : eltype(A)
                     tuned_alg = get_tuned_algorithm(eltype_A, eltype(b), matrix_size)
 
@@ -369,8 +371,8 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool})
                             eltype(b) <: Union{Float32, Float64, ComplexF32, ComplexF64}
                         DefaultAlgorithmChoice.AppleAccelerateLUFactorization
                     elseif (
-                            length(b) <= 100 || (isopenblas() && length(b) <= 500) ||
-                                (usemkl && length(b) <= 200)
+                            size(b, 1) <= 100 || (isopenblas() && size(b, 1) <= 500) ||
+                                (usemkl && size(b, 1) <= 200)
                         ) &&
                             (
                             A === nothing ? eltype(b) <: Union{Float32, Float64} :
diff --git a/src/factorization.jl b/src/factorization.jl
index 6cc596255..488ac0960 100644
--- a/src/factorization.jl
+++ b/src/factorization.jl
@@ -122,8 +122,10 @@ function _check_residual_safety(cache::LinearCache, alg, A_original, y)
     b = cache.b
     if cache.alg isa DefaultLinearSolver
         buf = cache.cacheval.residual_buf
-        if length(buf) != length(b)
-            resize!(buf, length(b))
+        if size(buf) != size(b)
+            # `resize!` only applies to vectors; matrix (batched) b just allocates.
+            buf = buf isa Vector && b isa AbstractVector ? resize!(buf, length(b)) :
+                similar(b)
         end
     else
         buf = similar(b)
diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl
index 547b7010e..dd203ca22 100644
--- a/src/iterative_wrappers.jl
+++ b/src/iterative_wrappers.jl
@@ -276,6 +276,19 @@ function init_cacheval(
     return nothing
 end
 
+# Krylov workspaces only support vector right-hand sides. Batched (matrix) `b`
+# with a Krylov algorithm errors informatively at `init` time
+# (`_check_batched_rhs_support`); this method only exists so the default
+# polyalgorithm can still initialize its (unused) Krylov cacheval slots when a
+# factorization algorithm is chosen for a batched problem.
+function init_cacheval(
+        alg::LinearSolve.KrylovJL, A, b::AbstractMatrix, u, Pl, Pr,
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        ::LinearSolve.OperatorAssumptions; zeroinit = true
+    )
+    return nothing
+end
+
 function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...)
     if cache.precsisfresh && !isnothing(alg.precs)
         Pl, Pr = alg.precs(cache.A, cache.p)
diff --git a/src/mkl.jl b/src/mkl.jl
index e794342b7..604e10bb3 100644
--- a/src/mkl.jl
+++ b/src/mkl.jl
@@ -386,7 +386,11 @@ function SciMLBase.solve!(
     if m > n
         Bc = copy(cache.b)
         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)
         getrs!('N', A.factors, A.ipiv, cache.u; info)
@@ -469,7 +473,11 @@ function SciMLBase.solve!(
     if m > n
         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)
         getrs!('N', A_lu.factors, A_lu.ipiv, u_32; info)
diff --git a/src/openblas.jl b/src/openblas.jl
index 2ca7ec0fd..83e0b1fa1 100644
--- a/src/openblas.jl
+++ b/src/openblas.jl
@@ -392,7 +392,11 @@ function SciMLBase.solve!(
     if m > n
         Bc = copy(cache.b)
         openblas_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)
         openblas_getrs!('N', A.factors, A.ipiv, cache.u; info)
@@ -476,7 +480,11 @@ function SciMLBase.solve!(
     if m > n
         openblas_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)
         openblas_getrs!('N', A_lu.factors, A_lu.ipiv, u_32; info)
diff --git a/src/simplelu.jl b/src/simplelu.jl
index fdfbfcccb..4ebef9357 100644
--- a/src/simplelu.jl
+++ b/src/simplelu.jl
@@ -224,6 +224,12 @@ function init_cacheval(
         alg::SimpleLUFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol,
         reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
     )
+    b isa AbstractMatrix && throw(
+        ArgumentError(
+            "SimpleLUFactorization supports only vector right-hand sides. Use " *
+                "`LUFactorization()` or `GenericLUFactorization()` for batched (matrix) `b`."
+        )
+    )
     return LUSolver(convert(AbstractMatrix, A))
 end
 

From b4299e679e112bb881b1d298a3d06bb9650e2d83 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Fri, 3 Jul 2026 11:43:02 -0400
Subject: [PATCH 2/8] 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
---
 test/Core/batch.jl | 180 +++++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl   |   1 +
 2 files changed, 181 insertions(+)
 create mode 100644 test/Core/batch.jl

diff --git a/test/Core/batch.jl b/test/Core/batch.jl
new file mode 100644
index 000000000..ab25ad7b9
--- /dev/null
+++ b/test/Core/batch.jl
@@ -0,0 +1,180 @@
+using LinearSolve, LinearAlgebra, SparseArrays, StaticArrays, Test
+using Random
+
+Random.seed!(1234)
+
+const n = 12
+const k = 4
+
+@testset "Dense factorizations vs A \\ B" begin
+    for T in (Float64, ComplexF64)
+        A = rand(T, n, n) + n * I
+        B = rand(T, n, k)
+        Xref = A \ B
+        algs = Any[
+            LUFactorization(), GenericLUFactorization(),
+            QRFactorization(), QRFactorization(ColumnNorm()),
+            SVDFactorization(),
+        ]
+        LinearSolve.usemkl && push!(algs, MKLLUFactorization())
+        LinearSolve.useopenblas && push!(algs, OpenBLASLUFactorization())
+        LinearSolve.appleaccelerate_isavailable() &&
+            push!(algs, AppleAccelerateLUFactorization())
+        @testset "$T $(nameof(typeof(alg)))" for alg in algs
+            sol = solve(LinearProblem(A, B), alg)
+            @test SciMLBase.successful_retcode(sol)
+            @test size(sol.u) == (n, k)
+            @test sol.u ≈ Xref
+        end
+
+        # SPD / symmetric algorithms
+        Aspd = Matrix(Hermitian(A' * A + n * I))
+        Awrap = T <: Complex ? Hermitian(Aspd) : Symmetric(Aspd)
+        sol = solve(LinearProblem(Awrap, B), CholeskyFactorization())
+        @test sol.u ≈ Awrap \ B
+
+        if T <: Real
+            sol = solve(LinearProblem(Symmetric(Aspd), B), BunchKaufmanFactorization())
+            @test sol.u ≈ Symmetric(Aspd) \ B
+
+            sol = solve(LinearProblem(A, B), NormalCholeskyFactorization())
+            @test sol.u ≈ Xref rtol = 1.0e-6
+        end
+
+        # Default algorithm
+        sol = solve(LinearProblem(A, B))
+        @test size(sol.u) == (n, k)
+        @test sol.u ≈ Xref
+    end
+end
+
+@testset "Residual safety check with matrix B" begin
+    A = rand(n, n) + n * I
+    B = rand(n, k)
+    sol = solve(LinearProblem(A, B), LUFactorization(residualsafety = true))
+    @test SciMLBase.successful_retcode(sol)
+    @test sol.u ≈ A \ B
+end
+
+@testset "Generic element types (GenericLUFactorization)" begin
+    A = big.(rand(n, n)) + n * I
+    B = big.(rand(n, k))
+    sol = solve(LinearProblem(A, B), GenericLUFactorization())
+    @test sol.u ≈ A \ B
+end
+
+@testset "Non-square (least squares / minimum norm), k != n" begin
+    m2, n2, k2 = 20, 10, 3
+    A = rand(m2, n2)
+    B = rand(m2, k2)
+    for alg in (QRFactorization(), SVDFactorization(), nothing)
+        sol = solve(LinearProblem(A, B), alg)
+        @test size(sol.u) == (n2, k2)
+        @test sol.u ≈ A \ B rtol = 1.0e-8
+    end
+end
+
+@testset "Structured matrices via default algorithm" begin
+    B = rand(n, k)
+
+    A = Diagonal(rand(n) .+ 1)
+    sol = solve(LinearProblem(A, B))
+    @test sol.u ≈ Matrix(A) \ B
+
+    A = Tridiagonal(rand(n - 1), rand(n) .+ 4, rand(n - 1))
+    sol = solve(LinearProblem(A, B))
+    @test sol.u ≈ Matrix(A) \ B
+
+    A = SymTridiagonal(fill(4.0, n), fill(1.0, n - 1))
+    sol = solve(LinearProblem(A, B))
+    @test sol.u ≈ Matrix(A) \ B
+end
+
+@testset "Sparse A" begin
+    Random.seed!(42)
+    As = sprand(30, 30, 0.3) + 30I
+    B = rand(30, k)
+    Xref = Matrix(As) \ B
+    # Default (PureKLU slot) and explicit sparse LU algorithms
+    for alg in (nothing, UMFPACKFactorization(), KLUFactorization())
+        sol = solve(LinearProblem(As, B), alg)
+        @test size(sol.u) == (30, k)
+        @test sol.u ≈ Xref
+    end
+
+    # Sparse Cholesky (CHOLMOD)
+    Asym = sparse(Symmetric(As + As'))
+    sol = solve(LinearProblem(Symmetric(Asym), B), CHOLMODFactorization())
+    @test sol.u ≈ Matrix(Asym) \ B
+
+    # Non-square sparse default (column-pivoted sparse QR), k != n
+    Ans = sprand(40, 20, 0.5)
+    Bns = rand(40, k)
+    sol = solve(LinearProblem(Ans, Bns))
+    @test size(sol.u) == (20, k)
+    @test sol.u ≈ Matrix(Ans) \ Bns rtol = 1.0e-8
+end
+
+@testset "Cache reuse: new b and new A" begin
+    A1 = rand(n, n) + n * I
+    A2 = rand(n, n) + n * I
+    B1 = rand(n, k)
+    B2 = rand(n, k)
+    for alg in (LUFactorization(), QRFactorization(), nothing)
+        cache = init(LinearProblem(A1, B1), alg)
+        @test solve!(cache).u ≈ A1 \ B1
+        cache.b = B2
+        @test solve!(cache).u ≈ A1 \ B2
+        # `cache.A = ...` aliases (in-place factorizations mutate it), so hand the
+        # cache its own copy and compare against the pristine A2.
+        cache.A = copy(A2)
+        @test solve!(cache).u ≈ A2 \ B2
+    end
+end
+
+@testset "Singular A returns failure retcode without throwing" begin
+    A = zeros(n, n)
+    B = rand(n, k)
+    sol = solve(LinearProblem(A, B), LUFactorization())
+    @test !SciMLBase.successful_retcode(sol)
+
+    # Default polyalgorithm on a rank-deficient matrix must not throw
+    # (QR fallback path with matrix B)
+    A = rand(n, n)
+    A[:, 1] .= A[:, 2]
+    sol = solve(LinearProblem(A, B))
+    @test sol.u isa Matrix
+end
+
+@testset "Iterative algorithms error informatively" begin
+    A = rand(n, n) + n * I
+    B = rand(n, k)
+    @test_throws ArgumentError solve(LinearProblem(A, B), KrylovJL_GMRES())
+    @test_throws ArgumentError solve(LinearProblem(A, B), SimpleGMRES())
+    @test_throws ArgumentError solve(LinearProblem(A, B), KrylovJL_CG())
+    @test_throws ArgumentError solve(LinearProblem(A, B), SimpleLUFactorization())
+    err = try
+        solve(LinearProblem(A, B), KrylovJL_GMRES())
+    catch e
+        e
+    end
+    @test occursin("Batched", err.msg)
+    @test occursin("KrylovJL", err.msg)
+end
+
+@testset "Static arrays" begin
+    A = (@SMatrix rand(4, 4)) + 4 * I
+    B = @SMatrix rand(4, 2)
+    sol = solve(LinearProblem(A, B))
+    @test sol.u ≈ A \ B
+end
+
+@testset "Single-column matrix B matches vector b" begin
+    A = rand(n, n) + n * I
+    b = rand(n)
+    B = reshape(copy(b), n, 1)
+    solvec = solve(LinearProblem(A, b), LUFactorization())
+    solmat = solve(LinearProblem(A, B), LUFactorization())
+    @test size(solmat.u) == (n, 1)
+    @test vec(solmat.u) ≈ solvec.u
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 118343b96..cd7817b9d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -62,6 +62,7 @@ else
         default = "All",
         core = function ()
             @time @safetestset "Basic Tests" include("Core/basictests.jl")
+            @time @safetestset "Batched RHS" include("Core/batch.jl")
             @time @safetestset "Return codes" include("Core/retcodes.jl")
             @time @safetestset "Re-solve" include("Core/resolve.jl")
             @time @safetestset "Zero Initialization Tests" include("Core/zeroinittests.jl")

From 843656c2d240ba6c95a7396c87fd1aef449ac146 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Fri, 3 Jul 2026 11:43:02 -0400
Subject: [PATCH 3/8] 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
---
 docs/src/release_notes.md    |  9 +++++++++
 docs/src/tutorials/linear.md | 20 ++++++++++++++++++++
 2 files changed, 29 insertions(+)

diff --git a/docs/src/release_notes.md b/docs/src/release_notes.md
index 6aa85f6ef..47987c702 100644
--- a/docs/src/release_notes.md
+++ b/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:
diff --git a/docs/src/tutorials/linear.md b/docs/src/tutorials/linear.md
index 5cffc444b..b305a0042 100644
--- a/docs/src/tutorials/linear.md
+++ b/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

From 9c30317a8f2e5d87816de1d0b460a1ac11d584a4 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Fri, 3 Jul 2026 14:12:05 -0400
Subject: [PATCH 4/8] 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>
---
 docs/Project.toml                               | 6 +++++-
 lib/LinearSolveAutotune/Project.toml            | 2 +-
 lib/LinearSolveAutotune/test/qa/Project.toml    | 2 +-
 lib/LinearSolvePyAMG/Project.toml               | 2 +-
 lib/LinearSolvePyAMG/test/qa/Project.toml       | 2 +-
 test/AD/Project.toml                            | 2 +-
 test/GPU/Project.toml                           | 2 +-
 test/LinearSolveElemental/Project.toml          | 2 +-
 test/LinearSolveGinkgo/Project.toml             | 2 +-
 test/LinearSolveHYPRE/Project.toml              | 2 +-
 test/LinearSolveMUMPS/Project.toml              | 2 +-
 test/LinearSolvePETSc/Project.toml              | 2 +-
 test/LinearSolveParU/Project.toml               | 2 +-
 test/LinearSolvePardiso/Project.toml            | 2 +-
 test/LinearSolvePartitionedSolvers/Project.toml | 2 +-
 test/LinearSolveSuperLUDIST/Project.toml        | 2 +-
 test/Trim/Project.toml                          | 2 +-
 test/qa/Project.toml                            | 2 +-
 18 files changed, 22 insertions(+), 18 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index b2f6de188..c1756e034 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,6 +6,10 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 
 [compat]
 Documenter = "1"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 LinearSolveAutotune = "1.1"
 SciMLOperators = "1"
+
+[sources]
+LinearSolve = {path = ".."}
+LinearSolveAutotune = {path = "../lib/LinearSolveAutotune"}
diff --git a/lib/LinearSolveAutotune/Project.toml b/lib/LinearSolveAutotune/Project.toml
index 3e7c4df57..da909d28e 100644
--- a/lib/LinearSolveAutotune/Project.toml
+++ b/lib/LinearSolveAutotune/Project.toml
@@ -44,7 +44,7 @@ FastLapackInterface = "2.0.4"
 GitHub = "5"
 LAPACK_jll = "3.12"
 LinearAlgebra = "1"
-LinearSolve = "3.39.2"
+LinearSolve = "3.39.2, 4"
 MKL_jll = "2025.2.0"
 Metal = "1.5"
 OpenBLAS_jll = "0.3"
diff --git a/lib/LinearSolveAutotune/test/qa/Project.toml b/lib/LinearSolveAutotune/test/qa/Project.toml
index 8e0ff5233..568312340 100644
--- a/lib/LinearSolveAutotune/test/qa/Project.toml
+++ b/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 = "3, 4"
 LinearSolveAutotune = "1"
 SciMLTesting = "1.6"
 Test = "1"
diff --git a/lib/LinearSolvePyAMG/Project.toml b/lib/LinearSolvePyAMG/Project.toml
index 0525312f7..5eb36042c 100644
--- a/lib/LinearSolvePyAMG/Project.toml
+++ b/lib/LinearSolvePyAMG/Project.toml
@@ -17,7 +17,7 @@ LinearSolve = {path = "../.."}
 [compat]
 CondaPkg = "0.2"
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 Pkg = "1"
 PythonCall = "0.9"
 SafeTestsets = "0.1"
diff --git a/lib/LinearSolvePyAMG/test/qa/Project.toml b/lib/LinearSolvePyAMG/test/qa/Project.toml
index 561566f17..bd46f93a3 100644
--- a/lib/LinearSolvePyAMG/test/qa/Project.toml
+++ b/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 = "3, 4"
 LinearSolvePyAMG = "1"
 SciMLTesting = "1.6"
 Test = "1"
diff --git a/test/AD/Project.toml b/test/AD/Project.toml
index ce0258f52..09c27eeb6 100644
--- a/test/AD/Project.toml
+++ b/test/AD/Project.toml
@@ -27,7 +27,7 @@ ForwardDiff = "0.10.38, 1"
 InteractiveUtils = "1.10"
 JET = "0.9, 0.11"
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 Mooncake = "0.5.15"
 RecursiveFactorization = "0.2.26"
 SafeTestsets = "0.1, 1"
diff --git a/test/GPU/Project.toml b/test/GPU/Project.toml
index 7e5350bfb..c6787b093 100644
--- a/test/GPU/Project.toml
+++ b/test/GPU/Project.toml
@@ -19,7 +19,7 @@ LinearSolve = {path = "/home/crackauc/sandbox/tmp_20260531_035946_10202/dg/fr_Li
 BlockDiagonals = "0.2"
 CUDSS = "0.7"
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1"
 SparseArrays = "1.10"
diff --git a/test/LinearSolveElemental/Project.toml b/test/LinearSolveElemental/Project.toml
index 30b3b3742..aa69233a1 100644
--- a/test/LinearSolveElemental/Project.toml
+++ b/test/LinearSolveElemental/Project.toml
@@ -14,7 +14,7 @@ LinearSolve = {path = "../.."}
 [compat]
 Elemental = "0.6.1"
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 Random = "1.10"
 SafeTestsets = "0.1, 1"
 SciMLBase = "2.148, 3"
diff --git a/test/LinearSolveGinkgo/Project.toml b/test/LinearSolveGinkgo/Project.toml
index 7ede803ee..bcf5caff8 100644
--- a/test/LinearSolveGinkgo/Project.toml
+++ b/test/LinearSolveGinkgo/Project.toml
@@ -13,7 +13,7 @@ LinearSolve = {path = "../.."}
 [compat]
 Ginkgo = "1"
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1"
 SparseArrays = "1.10"
diff --git a/test/LinearSolveHYPRE/Project.toml b/test/LinearSolveHYPRE/Project.toml
index d0accd4b7..3159c9ad2 100644
--- a/test/LinearSolveHYPRE/Project.toml
+++ b/test/LinearSolveHYPRE/Project.toml
@@ -15,7 +15,7 @@ LinearSolve = {path = "../.."}
 [compat]
 HYPRE = "1.7"
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 MPI = "0.20"
 Random = "1.10"
 SafeTestsets = "0.1, 1"
diff --git a/test/LinearSolveMUMPS/Project.toml b/test/LinearSolveMUMPS/Project.toml
index 7200abf7f..41bf13147 100644
--- a/test/LinearSolveMUMPS/Project.toml
+++ b/test/LinearSolveMUMPS/Project.toml
@@ -13,7 +13,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 MPI = "0.20"
 MUMPS = "1.4"
 SafeTestsets = "0.1, 1"
diff --git a/test/LinearSolvePETSc/Project.toml b/test/LinearSolvePETSc/Project.toml
index 307b8da01..12b23c3be 100644
--- a/test/LinearSolvePETSc/Project.toml
+++ b/test/LinearSolvePETSc/Project.toml
@@ -17,7 +17,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 MPI = "0.20"
 PETSc = "0.4.10"
 PartitionedArrays = "0.5"
diff --git a/test/LinearSolveParU/Project.toml b/test/LinearSolveParU/Project.toml
index 5db9db9c6..3622b2bbb 100644
--- a/test/LinearSolveParU/Project.toml
+++ b/test/LinearSolveParU/Project.toml
@@ -14,7 +14,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 ParU_jll = "1"
 Random = "1.10"
 SafeTestsets = "0.1, 1"
diff --git a/test/LinearSolvePardiso/Project.toml b/test/LinearSolvePardiso/Project.toml
index 447e50a31..d4d060000 100644
--- a/test/LinearSolvePardiso/Project.toml
+++ b/test/LinearSolvePardiso/Project.toml
@@ -13,7 +13,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 Pardiso = "1"
 Random = "1.10"
 SafeTestsets = "0.1, 1"
diff --git a/test/LinearSolvePartitionedSolvers/Project.toml b/test/LinearSolvePartitionedSolvers/Project.toml
index a16aa8c6a..da455830a 100644
--- a/test/LinearSolvePartitionedSolvers/Project.toml
+++ b/test/LinearSolvePartitionedSolvers/Project.toml
@@ -13,7 +13,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 LinearSolve = {path = "../.."}
 
 [compat]
-LinearSolve = "3"
+LinearSolve = "3, 4"
 MPI = "0.20"
 PartitionedArrays = "0.5"
 PartitionedSolvers = "0.3"
diff --git a/test/LinearSolveSuperLUDIST/Project.toml b/test/LinearSolveSuperLUDIST/Project.toml
index d17fa3323..c16708e74 100644
--- a/test/LinearSolveSuperLUDIST/Project.toml
+++ b/test/LinearSolveSuperLUDIST/Project.toml
@@ -9,7 +9,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 LinearSolve = {path = "../.."}
 
 [compat]
-LinearSolve = "3"
+LinearSolve = "3, 4"
 MPI = "0.20"
 SparseArrays = "1.10"
 SuperLUDIST = "1"
diff --git a/test/Trim/Project.toml b/test/Trim/Project.toml
index fa2c79e35..f90d14d6b 100644
--- a/test/Trim/Project.toml
+++ b/test/Trim/Project.toml
@@ -22,7 +22,7 @@ JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 test = ["JET"]
 
 [compat]
-LinearSolve = "3"
+LinearSolve = "3, 4"
 RecursiveFactorization = "0.2"
 SafeTestsets = "0.1, 1"
 SciMLBase = "2"
diff --git a/test/qa/Project.toml b/test/qa/Project.toml
index 9e34b8fbf..efa1307e1 100644
--- a/test/qa/Project.toml
+++ b/test/qa/Project.toml
@@ -12,7 +12,7 @@ LinearSolve = {path = "../.."}
 [compat]
 Aqua = "0.8"
 JET = "0.9, 0.11"
-LinearSolve = "3"
+LinearSolve = "3, 4"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1.6"
 Test = "<0.0.1, 1"

From 0c5e242d159db3ebf926ae6dae5639f495f61f48 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Fri, 3 Jul 2026 14:12:09 -0400
Subject: [PATCH 5/8] 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>
---
 test/Core/forwarddiff_overloads.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/Core/forwarddiff_overloads.jl b/test/Core/forwarddiff_overloads.jl
index 5578b1132..2f68e04f3 100644
--- a/test/Core/forwarddiff_overloads.jl
+++ b/test/Core/forwarddiff_overloads.jl
@@ -293,8 +293,8 @@ prob = LinearProblem(sparse(plain_A), b)
     for nchunk in (1, 2, 3)
         bd = [
             ForwardDiff.Dual{Nothing, Float64, nchunk}(
-                Float64(i), ForwardDiff.Partials(ntuple(k -> sin(i + k), nchunk))
-            ) for i in 1:5
+                    Float64(i), ForwardDiff.Partials(ntuple(k -> sin(i + k), nchunk))
+                ) for i in 1:5
         ]
         cache = LinearSolve.__init(LinearProblem(Asp, bd), PureKLUFactorization())
         @test eltype(cache.A) == Float64                 # A not promoted
@@ -303,9 +303,9 @@ prob = LinearProblem(sparse(plain_A), b)
         @test isapprox(ForwardDiff.value.(u), ForwardDiff.value.(uref); rtol = 1.0e-10)
         @test all(
             isapprox(
-                ForwardDiff.partials(u[i], j), ForwardDiff.partials(uref[i], j);
-                rtol = 1.0e-8, atol = 1.0e-12
-            ) for i in 1:5, j in 1:nchunk
+                    ForwardDiff.partials(u[i], j), ForwardDiff.partials(uref[i], j);
+                    rtol = 1.0e-8, atol = 1.0e-12
+                ) for i in 1:5, j in 1:nchunk
         )
     end
 

From d64bad89dee98d4b2e62256e82ba33997b459c4c Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Fri, 3 Jul 2026 14:41:40 -0400
Subject: [PATCH 6/8] 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>
---
 Project.toml | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 54c830227..7278a431a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -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"

From 46228a6f4cb61db95ba8642112f545be64b6ccb6 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Fri, 3 Jul 2026 15:01:40 -0400
Subject: [PATCH 7/8] 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>

From 7f0d6c480f157799d503e6f736e583298dd44f08 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Sat, 4 Jul 2026 16:35:36 -0400
Subject: [PATCH 8/8] 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>
---
 docs/Project.toml                               | 2 +-
 lib/LinearSolveAutotune/Project.toml            | 4 ++--
 lib/LinearSolveAutotune/test/qa/Project.toml    | 2 +-
 lib/LinearSolvePyAMG/Project.toml               | 4 ++--
 lib/LinearSolvePyAMG/test/qa/Project.toml       | 2 +-
 test/AD/Project.toml                            | 2 +-
 test/GPU/Project.toml                           | 2 +-
 test/LinearSolveElemental/Project.toml          | 2 +-
 test/LinearSolveGinkgo/Project.toml             | 2 +-
 test/LinearSolveHYPRE/Project.toml              | 2 +-
 test/LinearSolveMUMPS/Project.toml              | 2 +-
 test/LinearSolvePETSc/Project.toml              | 2 +-
 test/LinearSolveParU/Project.toml               | 2 +-
 test/LinearSolvePardiso/Project.toml            | 2 +-
 test/LinearSolvePartitionedSolvers/Project.toml | 2 +-
 test/LinearSolveSuperLUDIST/Project.toml        | 2 +-
 test/Trim/Project.toml                          | 2 +-
 test/qa/Project.toml                            | 2 +-
 18 files changed, 20 insertions(+), 20 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index c1756e034..ce9141bec 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,7 +6,7 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 
 [compat]
 Documenter = "1"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 LinearSolveAutotune = "1.1"
 SciMLOperators = "1"
 
diff --git a/lib/LinearSolveAutotune/Project.toml b/lib/LinearSolveAutotune/Project.toml
index da909d28e..6054ce65d 100644
--- a/lib/LinearSolveAutotune/Project.toml
+++ b/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, 4"
+LinearSolve = "4"
 MKL_jll = "2025.2.0"
 Metal = "1.5"
 OpenBLAS_jll = "0.3"
diff --git a/lib/LinearSolveAutotune/test/qa/Project.toml b/lib/LinearSolveAutotune/test/qa/Project.toml
index 568312340..9083a6594 100644
--- a/lib/LinearSolveAutotune/test/qa/Project.toml
+++ b/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, 4"
+LinearSolve = "4"
 LinearSolveAutotune = "1"
 SciMLTesting = "1.6"
 Test = "1"
diff --git a/lib/LinearSolvePyAMG/Project.toml b/lib/LinearSolvePyAMG/Project.toml
index 5eb36042c..e8cb5a04b 100644
--- a/lib/LinearSolvePyAMG/Project.toml
+++ b/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, 4"
+LinearSolve = "4"
 Pkg = "1"
 PythonCall = "0.9"
 SafeTestsets = "0.1"
diff --git a/lib/LinearSolvePyAMG/test/qa/Project.toml b/lib/LinearSolvePyAMG/test/qa/Project.toml
index bd46f93a3..0933d0eb5 100644
--- a/lib/LinearSolvePyAMG/test/qa/Project.toml
+++ b/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, 4"
+LinearSolve = "4"
 LinearSolvePyAMG = "1"
 SciMLTesting = "1.6"
 Test = "1"
diff --git a/test/AD/Project.toml b/test/AD/Project.toml
index 09c27eeb6..80ccdf35a 100644
--- a/test/AD/Project.toml
+++ b/test/AD/Project.toml
@@ -27,7 +27,7 @@ ForwardDiff = "0.10.38, 1"
 InteractiveUtils = "1.10"
 JET = "0.9, 0.11"
 LinearAlgebra = "1.10"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 Mooncake = "0.5.15"
 RecursiveFactorization = "0.2.26"
 SafeTestsets = "0.1, 1"
diff --git a/test/GPU/Project.toml b/test/GPU/Project.toml
index c6787b093..27f05c058 100644
--- a/test/GPU/Project.toml
+++ b/test/GPU/Project.toml
@@ -19,7 +19,7 @@ LinearSolve = {path = "/home/crackauc/sandbox/tmp_20260531_035946_10202/dg/fr_Li
 BlockDiagonals = "0.2"
 CUDSS = "0.7"
 LinearAlgebra = "1.10"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1"
 SparseArrays = "1.10"
diff --git a/test/LinearSolveElemental/Project.toml b/test/LinearSolveElemental/Project.toml
index aa69233a1..47d76b43f 100644
--- a/test/LinearSolveElemental/Project.toml
+++ b/test/LinearSolveElemental/Project.toml
@@ -14,7 +14,7 @@ LinearSolve = {path = "../.."}
 [compat]
 Elemental = "0.6.1"
 LinearAlgebra = "1.10"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 Random = "1.10"
 SafeTestsets = "0.1, 1"
 SciMLBase = "2.148, 3"
diff --git a/test/LinearSolveGinkgo/Project.toml b/test/LinearSolveGinkgo/Project.toml
index bcf5caff8..2c3afe028 100644
--- a/test/LinearSolveGinkgo/Project.toml
+++ b/test/LinearSolveGinkgo/Project.toml
@@ -13,7 +13,7 @@ LinearSolve = {path = "../.."}
 [compat]
 Ginkgo = "1"
 LinearAlgebra = "1.10"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1"
 SparseArrays = "1.10"
diff --git a/test/LinearSolveHYPRE/Project.toml b/test/LinearSolveHYPRE/Project.toml
index 3159c9ad2..c9a527f68 100644
--- a/test/LinearSolveHYPRE/Project.toml
+++ b/test/LinearSolveHYPRE/Project.toml
@@ -15,7 +15,7 @@ LinearSolve = {path = "../.."}
 [compat]
 HYPRE = "1.7"
 LinearAlgebra = "1.10"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 MPI = "0.20"
 Random = "1.10"
 SafeTestsets = "0.1, 1"
diff --git a/test/LinearSolveMUMPS/Project.toml b/test/LinearSolveMUMPS/Project.toml
index 41bf13147..7f6d285b9 100644
--- a/test/LinearSolveMUMPS/Project.toml
+++ b/test/LinearSolveMUMPS/Project.toml
@@ -13,7 +13,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 MPI = "0.20"
 MUMPS = "1.4"
 SafeTestsets = "0.1, 1"
diff --git a/test/LinearSolvePETSc/Project.toml b/test/LinearSolvePETSc/Project.toml
index 12b23c3be..2f478ce65 100644
--- a/test/LinearSolvePETSc/Project.toml
+++ b/test/LinearSolvePETSc/Project.toml
@@ -17,7 +17,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 MPI = "0.20"
 PETSc = "0.4.10"
 PartitionedArrays = "0.5"
diff --git a/test/LinearSolveParU/Project.toml b/test/LinearSolveParU/Project.toml
index 3622b2bbb..4891aab50 100644
--- a/test/LinearSolveParU/Project.toml
+++ b/test/LinearSolveParU/Project.toml
@@ -14,7 +14,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 ParU_jll = "1"
 Random = "1.10"
 SafeTestsets = "0.1, 1"
diff --git a/test/LinearSolvePardiso/Project.toml b/test/LinearSolvePardiso/Project.toml
index d4d060000..267ca098d 100644
--- a/test/LinearSolvePardiso/Project.toml
+++ b/test/LinearSolvePardiso/Project.toml
@@ -13,7 +13,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 Pardiso = "1"
 Random = "1.10"
 SafeTestsets = "0.1, 1"
diff --git a/test/LinearSolvePartitionedSolvers/Project.toml b/test/LinearSolvePartitionedSolvers/Project.toml
index da455830a..5dd4533dc 100644
--- a/test/LinearSolvePartitionedSolvers/Project.toml
+++ b/test/LinearSolvePartitionedSolvers/Project.toml
@@ -13,7 +13,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 LinearSolve = {path = "../.."}
 
 [compat]
-LinearSolve = "3, 4"
+LinearSolve = "4"
 MPI = "0.20"
 PartitionedArrays = "0.5"
 PartitionedSolvers = "0.3"
diff --git a/test/LinearSolveSuperLUDIST/Project.toml b/test/LinearSolveSuperLUDIST/Project.toml
index c16708e74..bab0d4e45 100644
--- a/test/LinearSolveSuperLUDIST/Project.toml
+++ b/test/LinearSolveSuperLUDIST/Project.toml
@@ -9,7 +9,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 LinearSolve = {path = "../.."}
 
 [compat]
-LinearSolve = "3, 4"
+LinearSolve = "4"
 MPI = "0.20"
 SparseArrays = "1.10"
 SuperLUDIST = "1"
diff --git a/test/Trim/Project.toml b/test/Trim/Project.toml
index f90d14d6b..9345e31f1 100644
--- a/test/Trim/Project.toml
+++ b/test/Trim/Project.toml
@@ -22,7 +22,7 @@ JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 test = ["JET"]
 
 [compat]
-LinearSolve = "3, 4"
+LinearSolve = "4"
 RecursiveFactorization = "0.2"
 SafeTestsets = "0.1, 1"
 SciMLBase = "2"
diff --git a/test/qa/Project.toml b/test/qa/Project.toml
index efa1307e1..b240834c5 100644
--- a/test/qa/Project.toml
+++ b/test/qa/Project.toml
@@ -12,7 +12,7 @@ LinearSolve = {path = "../.."}
 [compat]
 Aqua = "0.8"
 JET = "0.9, 0.11"
-LinearSolve = "3, 4"
+LinearSolve = "4"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1.6"
 Test = "<0.0.1, 1"