Add check for singularity

Author: YingboMa

src/ModelingToolkit.jl

@@ -123,7 +123,7 @@ Base.convert(::Type{<:Array{Num}}, x::AbstractArray) = map(Num, x)
 Base.convert(::Type{<:Array{Num}}, x::AbstractArray{Num}) = x
 Base.convert(::Type{Sym}, x::Num) = value(x) isa Sym ? value(x) : error("cannot convert $x to Sym")
 
-LinearAlgebra.lu(x::Array{Num}; kw...) = sym_lu(x)
+LinearAlgebra.lu(x::Array{Num}, ::Union{Val{false}, Val{true}}; check=true, kw...) = sym_lu(x; check=check)
 
 """
 $(TYPEDEF)

src/solve.jl

@@ -9,24 +9,35 @@ function nterms(t)
 end
 # Soft pivoted
 # Note: we call this function with a matrix of Union{SymbolicUtils.Symbolic, Any}
-# It should work as-is with Operation type too.
-function sym_lu(A)
+function sym_lu(A; check=true)
+    SINGULAR = typemax(Int)
     m, n = size(A)
     F = map(x->x isa Num ? x : Num(x), A)
     minmn = min(m, n)
     p = Vector{BlasInt}(undef, minmn)
-    info = zero(BlasInt)
+    info = 0
     for k = 1:minmn
-        val, i = findmin(map(ii->_iszero(F[ii, k]) ? Inf : nterms(F[ii,k]), k:n))
-        if !(val isa Symbolic) && (val isa Number) && val == Inf && iszero(info)
+        kp = k
+        amin = SINGULAR
+        for i in k:m
+            absi = _iszero(F[i, k]) ? SINGULAR : nterms(F[i,k])
+            if absi < amin
+                kp = i
+                amin = absi
+            end
+        end
+
+        p[k] = kp
+
+        if amin == SINGULAR && !(amin isa Symbolic) && (amin isa Number) && iszero(info)
             info = k
         end
-        i += k - 1
+
         # swap
         for j in 1:n
-            F[k, j], F[i, j] = F[i, j], F[k, j]
+            F[k, j], F[kp, j] = F[kp, j], F[k, j]
         end
-        p[k] = i
+
         for i in k+1:m
             F[i, k] = F[i, k] / F[k, k]
         end
@@ -36,7 +47,8 @@ function sym_lu(A)
             end
         end
     end
-    LU(F, p, info)
+    check && LinearAlgebra.checknonsingular(info, Val{true}())
+    LU(F, p, convert(BlasInt, info))
 end
 
 # Given a vector of equations and a

test/operation_overloads.jl

@@ -21,13 +21,19 @@ aa = a; # old a
 X = [0 b c; d e f; g h i]
 @test iszero(simplify(det(X) - ((b * f * g) + (c * d * h) - (b * d * i) - (c * e * g)), polynorm=true))
 F = lu(X)
+@test F.p == [2, 1, 3]
 R = simplify.(F.L * F.U - X[F.p, :], polynorm=true)
 @test iszero(R)
-@test simplify.(F \ X, polynorm=true) == I
+@test simplify.(F \ X) == I
 @test ModelingToolkit._solve(X, X) == I
 inv(X)
 qr(X)
 
+X2 = [0 b c; 0 0 0; 0 h 0]
+@test_throws SingularException lu(X2)
+F2 = lu(X2, check=false)
+@test F2.info == 1
+
 # test operations with sparse arrays and Operations
 # note `isequal` instead of `==` because `==` would give another Operation