diff --git a/Project.toml b/Project.toml
index b1af7bf35..7278a431a 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]
@@ -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"
diff --git a/docs/Project.toml b/docs/Project.toml
index b2f6de188..ce9141bec 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 = "4"
 LinearSolveAutotune = "1.1"
 SciMLOperators = "1"
+
+[sources]
+LinearSolve = {path = ".."}
+LinearSolveAutotune = {path = "../lib/LinearSolveAutotune"}
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
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/lib/LinearSolveAutotune/Project.toml b/lib/LinearSolveAutotune/Project.toml
index 3e7c4df57..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"
+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 8e0ff5233..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"
+LinearSolve = "4"
 LinearSolveAutotune = "1"
 SciMLTesting = "1.6"
 Test = "1"
diff --git a/lib/LinearSolvePyAMG/Project.toml b/lib/LinearSolvePyAMG/Project.toml
index 0525312f7..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"
+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 561566f17..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"
+LinearSolve = "4"
 LinearSolvePyAMG = "1"
 SciMLTesting = "1.6"
 Test = "1"
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
 
diff --git a/test/AD/Project.toml b/test/AD/Project.toml
index ce0258f52..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"
+LinearSolve = "4"
 Mooncake = "0.5.15"
 RecursiveFactorization = "0.2.26"
 SafeTestsets = "0.1, 1"
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/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
 
diff --git a/test/GPU/Project.toml b/test/GPU/Project.toml
index 7e5350bfb..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"
+LinearSolve = "4"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1"
 SparseArrays = "1.10"
diff --git a/test/LinearSolveElemental/Project.toml b/test/LinearSolveElemental/Project.toml
index 30b3b3742..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"
+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 7ede803ee..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"
+LinearSolve = "4"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1"
 SparseArrays = "1.10"
diff --git a/test/LinearSolveHYPRE/Project.toml b/test/LinearSolveHYPRE/Project.toml
index d0accd4b7..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"
+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 7200abf7f..7f6d285b9 100644
--- a/test/LinearSolveMUMPS/Project.toml
+++ b/test/LinearSolveMUMPS/Project.toml
@@ -13,7 +13,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+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 307b8da01..2f478ce65 100644
--- a/test/LinearSolvePETSc/Project.toml
+++ b/test/LinearSolvePETSc/Project.toml
@@ -17,7 +17,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+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 5db9db9c6..4891aab50 100644
--- a/test/LinearSolveParU/Project.toml
+++ b/test/LinearSolveParU/Project.toml
@@ -14,7 +14,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+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 447e50a31..267ca098d 100644
--- a/test/LinearSolvePardiso/Project.toml
+++ b/test/LinearSolvePardiso/Project.toml
@@ -13,7 +13,7 @@ LinearSolve = {path = "../.."}
 
 [compat]
 LinearAlgebra = "1.10"
-LinearSolve = "3"
+LinearSolve = "4"
 Pardiso = "1"
 Random = "1.10"
 SafeTestsets = "0.1, 1"
diff --git a/test/LinearSolvePartitionedSolvers/Project.toml b/test/LinearSolvePartitionedSolvers/Project.toml
index a16aa8c6a..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"
+LinearSolve = "4"
 MPI = "0.20"
 PartitionedArrays = "0.5"
 PartitionedSolvers = "0.3"
diff --git a/test/LinearSolveSuperLUDIST/Project.toml b/test/LinearSolveSuperLUDIST/Project.toml
index d17fa3323..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"
+LinearSolve = "4"
 MPI = "0.20"
 SparseArrays = "1.10"
 SuperLUDIST = "1"
diff --git a/test/Trim/Project.toml b/test/Trim/Project.toml
index fa2c79e35..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"
+LinearSolve = "4"
 RecursiveFactorization = "0.2"
 SafeTestsets = "0.1, 1"
 SciMLBase = "2"
diff --git a/test/qa/Project.toml b/test/qa/Project.toml
index 9e34b8fbf..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"
+LinearSolve = "4"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1.6"
 Test = "<0.0.1, 1"
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")