Merge pull request #683 from SciML/myb/linsol

Better LU factorization and linear solve

Author: ChrisRackauckas

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...) = lu(x, Val{false}(); kw...)
+LinearAlgebra.lu(x::Array{Num}; check=true, kw...) = sym_lu(x; check=check)
 
 """
 $(TYPEDEF)

src/solve.jl

@@ -9,35 +9,46 @@ 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)
-    L = fill!(Array{Any}(undef, size(A)),0) # TODO: make sparse?
-    for i=1:min(m, n)
-        L[i,i] = 1
-    end
-    U = copy!(Array{Any}(undef, size(A)),A)
-    p = BlasInt[1:m;]
-    for k = 1:m-1
-        _, i = findmin(map(ii->_iszero(U[ii, k]) ? Inf : nterms(U[ii,k]), k:n))
-        i += k - 1
+    F = map(x->x isa Num ? x : Num(x), A)
+    minmn = min(m, n)
+    p = Vector{BlasInt}(undef, minmn)
+    info = 0
+    for k = 1:minmn
+        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
+
         # swap
-        U[k, k:end], U[i, k:end] = U[i, k:end], U[k, k:end]
-        L[k, 1:k-1], L[i, 1:k-1] = L[i, 1:k-1], L[k, 1:k-1]
-        p[k] = i
+        for j in 1:n
+            F[k, j], F[kp, j] = F[kp, j], F[k, j]
+        end
 
-        for j = k+1:m
-            L[j,k] = U[j, k] / U[k, k]
-            U[j,k:m] .= U[j,k:m] .- (L[j,k],) .* U[k,k:m]
+        for i in k+1:m
+            F[i, k] = F[i, k] / F[k, k]
         end
-    end
-    for j=1:m
-        for i=j+1:n
-            U[i,j] = 0
+        for j = k+1:n
+            for i in k+1:m
+                F[i, j] = F[i, j] - F[i, k] * F[k, j]
+            end
         end
     end
-
-    (L, U, LinearAlgebra.ipiv2perm(p, m))
+    check && LinearAlgebra.checknonsingular(info, Val{true}())
+    LU(F, p, convert(BlasInt, info))
 end
 
 # Given a vector of equations and a
@@ -77,46 +88,56 @@ function solve_for(eqs, vars)
 end
 
 function _solve(A::AbstractMatrix, b::AbstractArray)
-    A = SymbolicUtils.simplify.(to_symbolic.(A), polynorm=true)
-    b = SymbolicUtils.simplify.(to_symbolic.(b), polynorm=true)
-    SymbolicUtils.simplify.(ldiv(sym_lu(A), b))
+    A = SymbolicUtils.simplify.(Num.(A), polynorm=true)
+    b = SymbolicUtils.simplify.(Num.(b), polynorm=true)
+    value.(SymbolicUtils.simplify.(sym_lu(A) \ b))
 end
 _solve(a, b) = value(SymbolicUtils.simplify(b/a, polynorm=true))
 
 # ldiv below
 
 _iszero(x::Number) = iszero(x)
 _isone(x::Number) = isone(x)
-_iszero(::Term) = false
-_isone(::Term) = false
-
-function simplifying_dot(x,y)
-    isempty(x) && return 0
-    muls = map(x,y) do xi,yi
-        _isone(xi) ? yi : _isone(yi) ? xi : _iszero(xi) ? 0 : _iszero(yi) ? 0 : xi * yi
+_iszero(::Symbolic) = false
+_isone(::Symbolic) = false
+_iszero(x::Num) = _iszero(value(x))
+_isone(x::Num) = _isone(value(x))
+
+LinearAlgebra.ldiv!(A::UpperTriangular{<:Union{Symbolic,Num}}, b::AbstractVector{<:Union{Symbolic,Num}}, x::AbstractVector{<:Union{Symbolic,Num}} = b) = symsub!(A, b, x)
+function symsub!(A::UpperTriangular, b::AbstractVector, x::AbstractVector = b)
+    LinearAlgebra.require_one_based_indexing(A, b, x)
+    n = size(A, 2)
+    if !(n == length(b) == length(x))
+        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
     end
-
-    reduce(muls) do acc, x
-        _iszero(acc) ? x : _iszero(x) ? acc : acc + x
+    @inbounds for j in n:-1:1
+        _iszero(A.data[j,j]) && throw(SingularException(j))
+        xj = x[j] = b[j] / A.data[j,j]
+        for i in j-1:-1:1
+            sub = _isone(xj) ? A.data[i,j] : A.data[i,j] * xj
+            if !_iszero(sub)
+                b[i] -= sub
+            end
+        end
     end
+    x
 end
 
-function ldiv((L,U,p), b)
-    m, n = size(L)
-    x = Vector{Any}(undef, length(b))
-    b = b[p]
-
-    for i=n:-1:1
-        sub = simplifying_dot(x[i+1:end], U[i,i+1:end])
-        den = U[i,i]
-        x[i] = _iszero(sub) ? b[i] : b[i] - sub
-        x[i] = _isone(den) ? x[i] : _isone(-den) ? -x[i] : x[i] / den
+LinearAlgebra.ldiv!(A::UnitLowerTriangular{<:Union{Symbolic,Num}}, b::AbstractVector{<:Union{Symbolic,Num}}, x::AbstractVector{<:Union{Symbolic,Num}} = b) = symsub!(A, b, x)
+function symsub!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector = b)
+    LinearAlgebra.require_one_based_indexing(A, b, x)
+    n = size(A, 2)
+    if !(n == length(b) == length(x))
+        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
     end
-
-    # unit lower triangular solve first:
-    for i=1:n
-        sub = simplifying_dot(b[1:i-1], L[i, 1:i-1]) # this should be `b` not x
-        x[i] = _iszero(sub) ? x[i] : x[i] - sub
+    @inbounds for j in 1:n
+        xj = x[j] = b[j]
+        for i in j+1:n
+            sub = _isone(xj) ? A.data[i,j] : A.data[i,j] * xj
+            if !_iszero(sub)
+                b[i] -= sub
+            end
+        end
     end
     x
 end

test/operation_overloads.jl

@@ -3,7 +3,7 @@ using LinearAlgebra
 using SparseArrays: sparse
 using Test
 
-@variables a,b,c,d
+@variables a,b,c,d,e,f,g,h,i
 
 # test hashing
 aa = a; # old a
@@ -18,12 +18,22 @@ aa = a; # old a
 @test hash(a+b ~ c+d) == hash(a+b ~ c+d)
 
 # test some matrix operations don't throw errors
-X = [a b;c d]
-det(X)
-lu(X)
+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) == 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