Use pinv if inv fails

Author: avik-pal

src/broyden.jl

@@ -155,7 +155,7 @@ function perform_step!(cache::GeneralBroydenCache{iip, IJ, UR}) where {iip, IJ,
     T = eltype(cache.u)
 
     if IJ === :true_jacobian && cache.stats.nsteps == 0
-        cache.J⁻¹ = inv(jacobian!!(cache.J⁻¹, cache)) # This allocates
+        cache.J⁻¹ = __safe_inv(jacobian!!(cache.J⁻¹, cache)) # This allocates
     end
 
     @bb cache.du = cache.J⁻¹ × vec(cache.fu)
@@ -179,7 +179,7 @@ function perform_step!(cache::GeneralBroydenCache{iip, IJ, UR}) where {iip, IJ,
             return nothing
         end
         if IJ === :true_jacobian
-            cache.J⁻¹ = inv(jacobian!!(cache.J⁻¹, cache))
+            cache.J⁻¹ = __safe_inv(jacobian!!(cache.J⁻¹, cache))
         else
             cache.inv_alpha = __initial_inv_alpha(cache.inv_alpha, cache.alpha_initial,
                 cache.u, cache.fu, cache.internalnorm)

src/utils.jl

@@ -526,3 +526,24 @@ end
 @inline __diag(x::AbstractMatrix) = diag(x)
 @inline __diag(x::AbstractVector) = x
 @inline __diag(x::Number) = x
+
+# Safe Inverse: Try to use `inv` but if lu fails use `pinv`
+function __safe_inv(A::StridedMatrix{T}) where {T}
+    LinearAlgebra.checksquare(A)
+    if istriu(A)
+        A_ = UpperTriangular(A)
+        issingular = any(iszero, @view(A_[diagind(A_)]))
+        !issingular && return triu!(parent(inv(A_)))
+    elseif istril(A)
+        A_ = LowerTriangular(A)
+        issingular = any(iszero, @view(A_[diagind(A_)]))
+        !issingular && return tril!(parent(inv(A_)))
+    else
+        F = lu(A; check = false)
+        if issuccess(F)
+            Ai = LinearAlgebra.inv!(F)
+            return convert(typeof(parent(Ai)), Ai)
+        end
+    end
+    return pinv(A)
+end