Skip to content

Commit

Permalink
import logabsdetogabsdet(F::UmfpackLU) from future
Browse files Browse the repository at this point in the history
This method is required for the sparse logabsdet rrules, but was
not included in Julia prior to v1.7.
  • Loading branch information
ElOceanografo committed Aug 15, 2023
1 parent 72e6a83 commit 7df39de
Show file tree
Hide file tree
Showing 2 changed files with 35 additions and 0 deletions.
1 change: 1 addition & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SparseInverseSubset = "dc90abb0-5640-4711-901d-7e5b23a2fada"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[compat]
Adapt = "3.4.0"
Expand Down
34 changes: 34 additions & 0 deletions src/rulesets/SparseArrays/sparsematrix.jl
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,40 @@ function rrule(::typeof(findnz), v::AbstractSparseVector)
return (I, V), findnz_pullback
end

if VERSION < v"1.7"
#=
This method for `logabsdet(F::UmfpackLU)` is required to calculate the (log)determinants
of sparse matrices, but was not defined prior to Julia v1.7. In order fo the rrules
for the determinants of sparse matrices below to work, they need to be able to
compute the primals as well, so this import from the future is included. For more
recent versions of Julia, this definition lives in:
julia/stdlib/SuiteSparse/src/umfpack.jl
=#
using SuiteSparse.UMFPACK: _signperm, UmfpackLU

for itype in (:Int32, :Int64)
@eval begin
function LinearAlgebra.logabsdet(F::UmfpackLU{T, $itype}) where {T<:Union{Float64,ComplexF64}}
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
end
end
end


function rrule(::typeof(logabsdet), x::SparseMatrixCSC)
F = cholesky(x)
Expand Down

0 comments on commit 7df39de

Please sign in to comment.