Fix GPU Cholesky cache initialization for non-square matrices

When using DefaultLinearSolver with non-square GPU matrices (e.g., for
least squares problems), the init_cacheval function for CholeskyFactorization
would fail because it tried to compute cholesky(A) on a non-square matrix.

The fix checks assumptions.issq before attempting Cholesky factorization
and returns nothing for non-square matrices, allowing the DefaultLinearSolver
to properly use QRFactorization instead.

Fixes SciML/NonlinearSolve.jl#746

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

Author: ChrisRackauckas

src/factorization.jl

@@ -317,21 +317,24 @@ end
 const PREALLOCATED_QR_ColumnNorm = ArrayInterface.qr_instance(rand(1, 1), ColumnNorm())
 
 function init_cacheval(alg::QRFactorization{ColumnNorm}, A::Matrix{Float64}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     return PREALLOCATED_QR_ColumnNorm
 end
 
 function init_cacheval(
         alg::QRFactorization, A::Union{<:Adjoint, <:Transpose}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     A isa GPUArraysCore.AnyGPUArray && return qr(A)
     return qr(A, alg.pivot)
 end
 
 const PREALLOCATED_QR_NoPivot = ArrayInterface.qr_instance(rand(1, 1))
 
 function init_cacheval(alg::QRFactorization{NoPivot}, A::Matrix{Float64}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     return PREALLOCATED_QR_NoPivot
 end
 
@@ -388,13 +391,18 @@ function init_cacheval(alg::CholeskyFactorization, A::SMatrix{S1, S2}, b, u, Pl,
 end
 
 function init_cacheval(alg::CholeskyFactorization, A::GPUArraysCore.AnyGPUArray, b, u, Pl,
-        Pr, maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
-    cholesky(A; check = false)
+        Pr, maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
+    # Cholesky requires square matrices - return nothing for non-square to avoid errors
+    # during DefaultLinearSolver cache initialization
+    # See https://github.com/SciML/NonlinearSolve.jl/issues/746
+    assumptions.issq ? cholesky(A; check = false) : nothing
 end
 
 function init_cacheval(
         alg::CholeskyFactorization, A::AbstractArray{<:BLASELTYPES}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     if LinearSolve.is_cusparse_csc(A)
         nothing
     elseif LinearSolve.is_cusparse_csr(A) && !LinearSolve.cudss_loaded(A)
@@ -1012,7 +1020,8 @@ const PREALLOCATED_NORMALCHOLESKY_SYMMETRIC = ArrayInterface.cholesky_instance(
     Symmetric(rand(1, 1)), NoPivot())
 
 function init_cacheval(alg::NormalCholeskyFactorization, A::Matrix{Float64}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     return PREALLOCATED_NORMALCHOLESKY_SYMMETRIC
 end
 
@@ -1164,7 +1173,8 @@ function init_cacheval(alg::SparspakFactorization,
 end
 
 function init_cacheval(::SparspakFactorization, ::StaticArray, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     nothing
 end
 
@@ -1190,9 +1200,8 @@ struct CliqueTreesFactorization{A, S} <: AbstractSparseFactorization
             alg::A = nothing,
             snd::S = nothing,
             reuse_symbolic = true,
-            throwerror = true,
-        ) where {A, S}
-
+            throwerror = true
+    ) where {A, S}
         ext = Base.get_extension(@__MODULE__, :LinearSolveCliqueTreesExt)
 
         if throwerror && isnothing(ext)
@@ -1203,30 +1212,36 @@ struct CliqueTreesFactorization{A, S} <: AbstractSparseFactorization
     end
 end
 
-function init_cacheval(::CliqueTreesFactorization, ::Union{AbstractMatrix, Nothing, AbstractSciMLOperator}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+function init_cacheval(::CliqueTreesFactorization,
+        ::Union{AbstractMatrix, Nothing, AbstractSciMLOperator}, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     nothing
 end
 
 function init_cacheval(::CliqueTreesFactorization, ::StaticArray, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     nothing
 end
 
 # Fallback init_cacheval for extension-based algorithms when extensions aren't loaded
 # These return nothing since the actual implementations are in the extensions
 function init_cacheval(::BLISLUFactorization, A, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     nothing
 end
 
 function init_cacheval(::CudaOffloadLUFactorization, A, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     nothing
 end
 
 function init_cacheval(::MetalLUFactorization, A, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions)
+        maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool},
+        assumptions::OperatorAssumptions)
     nothing
 end
 

test/gpu/cuda.jl

@@ -75,7 +75,7 @@ x2 = zero(b);
 prob1 = LinearProblem(A1, b1; u0 = x1)
 prob2 = LinearProblem(A2, b2; u0 = x2)
 
-cache_kwargs = (;abstol = 1e-8, reltol = 1e-8, maxiter = 30)
+cache_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiter = 30)
 
 function test_interface(alg, prob1, prob2)
     A1 = prob1.A
@@ -103,7 +103,8 @@ function test_interface(alg, prob1, prob2)
     return
 end
 
-@testset "$alg" for alg in (CudaOffloadLUFactorization(), CudaOffloadQRFactorization(), NormalCholeskyFactorization())
+@testset "$alg" for alg in (CudaOffloadLUFactorization(), CudaOffloadQRFactorization(),
+    NormalCholeskyFactorization())
     test_interface(alg, prob1, prob2)
 end
 
@@ -171,3 +172,45 @@ if Base.find_package("CUSOLVERRF") !== nothing
         include("cusolverrf.jl")
     end
 end
+
+# Test for non-square GPU matrices (least squares problems)
+# See https://github.com/SciML/NonlinearSolve.jl/issues/746
+@testset "Non-square GPU matrices" begin
+    # Overdetermined system: more rows than columns (4x2)
+    A_rect = cu(Float32[1.0 2.0; 3.0 4.0; 5.0 6.0; 7.0 8.0])
+    b_rect = cu(Float32[1.0, 2.0, 3.0, 4.0])
+
+    prob_rect = LinearProblem(A_rect, b_rect)
+
+    # Test that default solver works (should use QRFactorization)
+    @testset "Default solver for non-square" begin
+        sol = solve(prob_rect)
+        @test sol.alg.alg == LinearSolve.DefaultAlgorithmChoice.QRFactorization
+        # Verify least squares solution
+        @test norm(A_rect * sol.u - b_rect) < norm(b_rect)  # residual should be smaller than b
+    end
+
+    # Test explicit QRFactorization
+    @testset "QRFactorization for non-square" begin
+        sol = solve(prob_rect, QRFactorization())
+        @test norm(A_rect * sol.u - b_rect) < norm(b_rect)
+    end
+
+    # Test NormalCholeskyFactorization (should work via A'*A)
+    @testset "NormalCholeskyFactorization for non-square" begin
+        sol = solve(prob_rect, NormalCholeskyFactorization())
+        @test norm(A_rect * sol.u - b_rect) < norm(b_rect)
+    end
+
+    # Underdetermined system: more columns than rows (2x4)
+    A_under = cu(Float32[1.0 2.0 3.0 4.0; 5.0 6.0 7.0 8.0])
+    b_under = cu(Float32[1.0, 2.0])
+
+    prob_under = LinearProblem(A_under, b_under)
+
+    @testset "Default solver for underdetermined" begin
+        sol = solve(prob_under)
+        # Should still work and give a solution
+        @test norm(A_under * sol.u - b_under) < 1e-4
+    end
+end