Add singularity checks for lu, method for issuccess. (#766)

This makes the implementation of `lu` conform to the API of
LinearAlgebra: check by default, optionally skip check.

Tests and methods that use `lu` are modified accordingly.

Author: tpapp

src/StaticArrays.jl

@@ -17,7 +17,7 @@ import Base: sqrt, exp, log
 using LinearAlgebra
 import LinearAlgebra: transpose, adjoint, dot, eigvals, eigen, lyap, tr,
                       kron, diag, norm, dot, diagm, lu, svd, svdvals,
-                      factorize, ishermitian, issymmetric, isposdef, normalize,
+                      factorize, ishermitian, issymmetric, isposdef, issuccess, normalize,
                       normalize!, Eigen, det, logdet, cross, diff, qr, \
 using LinearAlgebra: checksquare
 

src/det.jl

@@ -47,8 +47,8 @@ end
     if prod(S) ≤ 14*14
         quote
             @_inline_meta
-            LUp = lu(A)
-            det(LUp.U)*_parity(LUp.p)
+            _, U, p = _lu(A, Val(true), false)
+            det(UpperTriangular(U))*_parity(p)
         end
     else
         :(@_inline_meta; det(Matrix(A)))
@@ -62,7 +62,7 @@ end
     if prod(S) ≤ 14*14
         quote
             @_inline_meta
-            LUp = lu(A)
+            LUp = lu(A; check = false)
             d, s = logabsdet(LUp.U)
             d + log(s*_parity(LUp.p))
         end

src/lu.jl

@@ -30,25 +30,54 @@ function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, F::LU)
 end
 
 # LU decomposition
-function lu(A::StaticMatrix, pivot::Union{Val{false},Val{true}}=Val(true))
-    L, U, p = _lu(A, pivot)
+function lu(A::StaticMatrix, pivot::Union{Val{false},Val{true}}=Val(true); check = true)
+    L, U, p = _lu(A, pivot, check)
     LU(L, U, p)
 end
 
 # For the square version, return explicit lower and upper triangular matrices.
 # We would do this for the rectangular case too, but Base doesn't support that.
-function lu(A::StaticMatrix{N,N}, pivot::Union{Val{false},Val{true}}=Val(true)) where {N}
-    L, U, p = _lu(A, pivot)
+function lu(A::StaticMatrix{N,N}, pivot::Union{Val{false},Val{true}}=Val(true);
+            check = true) where {N}
+    L, U, p = _lu(A, pivot, check)
     LU(LowerTriangular(L), UpperTriangular(U), p)
 end
 
-@generated function _lu(A::StaticMatrix{M,N,T}, pivot) where {M,N,T}
+# location of the first zero on the diagonal, 0 when not found
+function _first_zero_on_diagonal(A::StaticMatrix{M,N,T}) where {M,N,T}
+    if @generated
+        quote
+            $(map(i -> :(A[$i, $i] == zero(T) && return $i), 1:min(M, N))...)
+            0
+        end
+    else
+        for i in 1:min(M, N)
+            A[i, i] == 0 && return i
+        end
+        0
+    end
+end
+
+function _first_zero_on_diagonal(A::LinearAlgebra.AbstractTriangular{<:Any,<:StaticMatrix})
+    _first_zero_on_diagonal(A.data)
+end
+
+issuccess(F::LU) = _first_zero_on_diagonal(F.U) == 0
+
+@generated function _lu(A::StaticMatrix{M,N,T}, pivot, check) where {M,N,T}
     if M*N ≤ 14*14
-        :(__lu(A, pivot))
+        quote
+            L, U, P = __lu(A, pivot)
+            if check
+                i = _first_zero_on_diagonal(U)
+                i == 0 || throw(SingularException(i))
+            end
+            L, U, P
+        end
     else
         quote
             # call through to Base to avoid excessive time spent on type inference for large matrices
-            f = lu(Matrix(A), pivot; check = false)
+            f = lu(Matrix(A), pivot; check = check)
             # Trick to get the output eltype - can't rely on the result of f.L as
             # it's not type inferrable.
             T2 = arithmetic_closure(T)

test/lu.jl

@@ -11,7 +11,7 @@ end
 
 @testset "LU decomposition ($m×$n, pivot=$pivot)" for pivot in (true, false), m in [0:4..., 15], n in [0:4..., 15]
     a = SMatrix{m,n,Int}(1:(m*n))
-    l, u, p = @inferred(lu(a, Val{pivot}()))
+    l, u, p = @inferred(lu(a, Val{pivot}(); check = false))
 
     # expected types
     @test p isa SVector{m,Int}
@@ -57,5 +57,11 @@ end
     # test if / and \ work with lu:
     @test a\b_col ≈ a_lu\b_col
     @test b_line/a ≈ b_line/a_lu
+end
 
+@testset "LU singularity check" for m in [2, 3, 20], n in [2, 3, 20]
+    # NOTE: large dimensions test fallback to LinearAlgebra.lu
+    A = ones(SMatrix{m,n})
+    @test_throws SingularException lu(A)
+    @test !issuccess(lu(A; check = false))
 end