Add logabsdet method for UmfpackLU (JuliaLang#40716)

Co-authored-by: Daniel Karrasch <[email protected]>

Author: 2 people

stdlib/SuiteSparse/src/umfpack.jl

@@ -6,7 +6,7 @@ export UmfpackLU
 
 import Base: (\), getproperty, show, size
 using LinearAlgebra
-import LinearAlgebra: Factorization, det, lu, lu!, ldiv!
+import LinearAlgebra: Factorization, checksquare, det, logabsdet, lu, lu!, ldiv!
 
 using SparseArrays
 using SparseArrays: getcolptr
@@ -279,6 +279,26 @@ function deserialize(s::AbstractSerializer, t::Type{UmfpackLU{Tv,Ti}}) where {Tv
     return obj
 end
 
+# compute the sign/parity of a permutation
+function _signperm(p)
+    n = length(p)
+    result = 0
+    todo = trues(n)
+    while any(todo)
+        k = findfirst(todo)
+        todo[k] = false
+        result += 1 # increment element count
+        j = p[k]
+        while j != k
+            result += 1 # increment element count
+            todo[j] = false
+            j = p[j]
+        end
+        result += 1 # increment cycle count
+    end
+    return ifelse(isodd(result), -1, 1)
+end
+
 ## Wrappers for UMFPACK functions
 
 # generate the name of the C function according to the value and integer types
@@ -406,6 +426,23 @@ for itype in UmfpackIndexTypes
                         mx, mz, C_NULL, lu.numeric, umf_info)
             complex(mx[], mz[])
         end
+        function logabsdet(F::UmfpackLU{T, $itype}) where {T<:Union{Float64,ComplexF64}} # return log(abs(det)) and sign(det)
+            n = checksquare(F)
+            issuccess(F) || return log(zero(real(T))), zero(T)
+            U = F.U
+            Rs = F.Rs
+            p = F.p
+            q = F.q
+            s = _signperm(p)*_signperm(q)*one(real(T))
+            P = one(T)
+            abs_det = zero(real(T))
+            @inbounds for i in 1:n
+                dg_ii = U[i, i] / Rs[i]
+                P *= sign(dg_ii)
+                abs_det += log(abs(dg_ii))
+            end
+            return abs_det, s * P
+        end
         function umf_lunz(lu::UmfpackLU{Float64,$itype})
             lnz = Ref{$itype}()
             unz = Ref{$itype}()

stdlib/SuiteSparse/test/umfpack.jl

@@ -34,7 +34,11 @@ using SparseArrays: nnz, sparse, sprand, sprandn, SparseMatrixCSC
             @test Rs == lua.Rs
             @test (Diagonal(Rs) * A)[p,q] ≈ L * U
 
-            det(lua) ≈ det(Array(A))
+            @test det(lua) ≈ det(Array(A))
+            logdet_lua, sign_lua = logabsdet(lua)
+            logdet_A, sign_A = logabsdet(Array(A))
+            @test logdet_lua ≈ logdet_A
+            @test sign_lua ≈ sign_A
 
             b = [8., 45., -3., 3., 19.]
             x = lua\b