Restore dense pinv fallback for inverse Jacobian init

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Author: ChrisRackauckas

lib/NonlinearSolveBase/ext/NonlinearSolveBaseSparseArraysExt.jl

@@ -87,7 +87,7 @@ dense buffer before factorizing.
 """
 function Utils.linsolve_workspace(A::AbstractSparseMatrix)
     dense_A = Matrix(A)
-    workspace, _ = Utils.linsolve_workspace(dense_A)
+    workspace, _ = invoke(Utils.linsolve_workspace, Tuple{AbstractMatrix}, dense_A)
     return workspace, dense_A
 end
 

lib/NonlinearSolveBase/src/utils.jl

@@ -226,13 +226,13 @@ function condition_number(J::AbstractVector)
 end
 condition_number(::Any) = -1
 
-# Explicit inverse-Jacobian initialization for quasi-Newton update rules that store J⁻¹:
-# solve A * X = I with the default linear solver. Scalars, `Diagonal`s, and static
-# matrices take native fast paths (mirroring `construct_linear_solver`'s routing) instead
-# of the linear-solve cache.
+# Explicit inverse-Jacobian initialization for quasi-Newton update rules that store J⁻¹.
+# Dense matrices keep the previous `pinv` fallback semantics for rank-deficient
+# reinitializations; sparse and structured matrices solve A * X = I through LinearSolve.
 linsolve_workspace(A) = nothing, A
 linsolve_workspace(A::Diagonal) = nothing, A
 linsolve_workspace(A::SMatrix) = nothing, A
+linsolve_workspace(A::StridedMatrix) = nothing, A
 function linsolve_workspace(A::AbstractMatrix)
     LinearAlgebra.checksquare(A)
     # the linear solve promotes e.g. integer eltypes; the buffers must hold the promoted
@@ -268,6 +268,25 @@ function linsolve_identity!!(workspace, A::Diagonal)
     return Diagonal(D)
 end
 linsolve_identity!!(workspace, A::AbstractVector) = linsolve_identity!!(workspace, Diagonal(A))
+function linsolve_identity!!(workspace, A::StridedMatrix)
+    LinearAlgebra.checksquare(A)
+    if LinearAlgebra.istriu(A)
+        issingular = any(iszero, @view(A[diagind(A)]))
+        A_ = LinearAlgebra.UpperTriangular(A)
+        !issingular && return LinearAlgebra.triu!(parent(inv(A_)))
+    elseif LinearAlgebra.istril(A)
+        A_ = LinearAlgebra.LowerTriangular(A)
+        issingular = any(iszero, @view(A_[diagind(A_)]))
+        !issingular && return LinearAlgebra.tril!(parent(inv(A_)))
+    else
+        F = LinearAlgebra.lu(A; check = false)
+        if LinearAlgebra.issuccess(F)
+            Ai = LinearAlgebra.inv!(F)
+            return convert(typeof(parent(Ai)), Ai)
+        end
+    end
+    return pinv(A)
+end
 # Static matrices solve `A X = I` through LinearSolve's static fast path directly (NOT
 # `construct_linear_solver`, whose `SMatrix` route is a plain `A \ b` that yields Inf/NaN
 # on a singular Jacobian). The static default's SVD rescue returns the finite min-norm

lib/NonlinearSolveBase/test/linsolve_workspace.jl

@@ -1,5 +1,5 @@
 using NonlinearSolveBase: NonlinearSolveBase, Utils
-using LinearSolve  # linsolve_identity!! routes matrices through the LinearSolve default
+using LinearSolve  # used by sparse and structured linsolve_identity!! workspaces
 using LinearAlgebra
 using Test
 
@@ -20,20 +20,15 @@ using Test
         @test Utils.linsolve_identity!!(workspace, A) ≈ inv(A_orig)
         # nothing workspace (preinverted-init path in QuasiNewton) still works
         @test Utils.linsolve_identity!!(nothing, A) ≈ inv(A_orig)
-        # re-inverting the returned inverse (which aliases the workspace's solution
-        # buffer) recovers A
+        # re-inverting the returned inverse recovers A
         Ai_copy = copy(Utils.linsolve_identity!!(workspace, A))
         @test Utils.linsolve_identity!!(workspace, Ai_copy) ≈ A_orig
         @test Utils.linsolve_identity!!(workspace, Utils.linsolve_identity!!(workspace, A)) ≈
             A_orig
     end
 end
 
-@testset "singular input takes the pivoted-QR rescue" begin
-    # The result is the LinearSolve default algorithm's least-squares generalized
-    # inverse from its singular-LU → pivoted-QR rescue, NOT the SVD `pinv` (an
-    # intentional semantics change: same fitness for quasi-Newton initialization,
-    # different finite matrix).
+@testset "dense singular input matches pinv" begin
     n = 20
     A_singular = let B = rand(n, n)
         B[:, 1] .= 0
@@ -45,32 +40,29 @@ end
     for A in (A_singular, A_singular_triu)
         A_orig = copy(A)
         workspace, _ = Utils.linsolve_workspace(A)
+        @test workspace === nothing
         X = Utils.linsolve_identity!!(workspace, A)
         @test size(X) == size(A)
         @test all(isfinite, X)
-        # X is a least-squares generalized inverse: A * X projects onto range(A)
+        @test X ≈ pinv(A_orig)
         @test A * X * A ≈ A atol = 1.0e-8 * norm(A)
         @test A == A_orig
-        # a subsequent nonsingular solve with the same workspace is unaffected
         B = rand(n, n) + n * I
         @test Utils.linsolve_identity!!(workspace, B) ≈ inv(B)
     end
 end
 
-@testset "workspace reuse does not allocate a matrix copy" begin
+@testset "strided matrices use the native inverse path" begin
     n = 51
     A = rand(n, n) + n * I
-    workspace, _ = Utils.linsolve_workspace(A)
-    Utils.linsolve_identity!!(workspace, A)  # compile
-    # Steady state is the LinearSolve refactorization floor (`lu!` ipiv + wrapper); the
-    # old code copied A (`lu`) and allocated a LAPACK getri work array, both O(n^2).
-    allocs = @allocated Utils.linsolve_identity!!(workspace, A)
-    @test allocs < sizeof(A)
+    workspace, A_ret = Utils.linsolve_workspace(A)
+    @test workspace === nothing && A_ret === A
+    @test Utils.linsolve_identity!!(workspace, A) ≈ inv(A)
 
     A_triu = triu(A)
-    Utils.linsolve_identity!!(workspace, A_triu)  # compile
-    allocs_tri = @allocated Utils.linsolve_identity!!(workspace, A_triu)
-    @test allocs_tri < sizeof(A)
+    workspace, A_triu_ret = Utils.linsolve_workspace(A_triu)
+    @test workspace === nothing && A_triu_ret === A_triu
+    @test Utils.linsolve_identity!!(workspace, A_triu) ≈ inv(A_triu)
 end
 
 @testset "sparse matrices route through the lincache (no sparse pinv)" begin