From 6362c920d3323cfeb2ce0700726d6feb0b59c567 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sun, 8 May 2022 17:33:00 +0200 Subject: [PATCH] Add definition of `det` --- Project.toml | 2 +- src/pdiagmat.jl | 1 + src/pdmat.jl | 1 + src/pdsparsemat.jl | 1 + src/scalmat.jl | 1 + test/testutils.jl | 20 +++++++++++++++++--- 6 files changed, 22 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 919d98c..c09c784 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "PDMats" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.8" +version = "0.11.9" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl index 3753d9f..bba7280 100644 --- a/src/pdiagmat.jl +++ b/src/pdiagmat.jl @@ -62,6 +62,7 @@ Base.kron(A::PDiagMat, B::PDiagMat) = PDiagMat( vcat([A.diag[i] * B.diag for i i ### Algebra Base.inv(a::PDiagMat) = PDiagMat(map(inv, a.diag)) +LinearAlgebra.det(a::PDiagMat) = prod(a.diag) function LinearAlgebra.logdet(a::PDiagMat) diag = a.diag return isempty(diag) ? zero(log(zero(eltype(diag)))) : sum(log, diag) diff --git a/src/pdmat.jl b/src/pdmat.jl index 9213072..c15c25c 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -53,6 +53,7 @@ end ### Algebra Base.inv(a::PDMat) = PDMat(inv(a.chol)) +LinearAlgebra.det(a::PDMat) = det(a.chol) LinearAlgebra.logdet(a::PDMat) = logdet(a.chol) LinearAlgebra.eigmax(a::PDMat) = eigmax(a.mat) LinearAlgebra.eigmin(a::PDMat) = eigmin(a.mat) diff --git a/src/pdsparsemat.jl b/src/pdsparsemat.jl index 6a949b7..fa0c6ad 100644 --- a/src/pdsparsemat.jl +++ b/src/pdsparsemat.jl @@ -51,6 +51,7 @@ end ### Algebra Base.inv(a::PDSparseMat{T}) where {T<:Real} = PDMat(inv(a.mat)) +LinearAlgebra.det(a::PDSparseMat) = det(a.chol) LinearAlgebra.logdet(a::PDSparseMat) = logdet(a.chol) ### whiten and unwhiten diff --git a/src/scalmat.jl b/src/scalmat.jl index 3ed50f6..53e771e 100644 --- a/src/scalmat.jl +++ b/src/scalmat.jl @@ -61,6 +61,7 @@ Base.kron(A::ScalMat, B::ScalMat) = ScalMat( dim(A) * dim(B), A.value * B.value ### Algebra Base.inv(a::ScalMat) = ScalMat(a.dim, inv(a.value)) +LinearAlgebra.det(a::ScalMat) = a.value^a.dim LinearAlgebra.logdet(a::ScalMat) = a.dim * log(a.value) LinearAlgebra.eigmax(a::ScalMat) = a.value LinearAlgebra.eigmin(a::ScalMat) = a.value diff --git a/test/testutils.jl b/test/testutils.jl index b2e8a00..8b373c3 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -17,6 +17,7 @@ function test_pdmat(C::AbstractPDMat, Cmat::Matrix; t_cholesky::Bool=true, # whether to test cholesky method t_scale::Bool=true, # whether to test scaling t_add::Bool=true, # whether to test pdadd + t_det::Bool=true, # whether to test det method t_logdet::Bool=true, # whether to test logdet method t_eig::Bool=true, # whether to test eigmax and eigmin t_mul::Bool=true, # whether to test multiplication @@ -36,6 +37,7 @@ function test_pdmat(C::AbstractPDMat, Cmat::Matrix; isa(C, PDMatType) && t_cholesky && pdtest_cholesky(C, Cmat, cmat_eq, verbose) t_scale && pdtest_scale(C, Cmat, verbose) t_add && pdtest_add(C, Cmat, verbose) + t_det && pdtest_det(C, Cmat, verbose) t_logdet && pdtest_logdet(C, Cmat, verbose) t_eig && pdtest_eig(C, Cmat, verbose) @@ -154,12 +156,24 @@ function pdtest_add(C::AbstractPDMat, Cmat::Matrix, verbose::Int) @test Mr ≈ R end +function pdtest_det(C::AbstractPDMat, Cmat::Matrix, verbose::Int) + _pdt(verbose, "det") + @test det(C) ≈ det(Cmat) + + # generic fallback in LinearAlgebra performs LU decomposition + if C isa Union{PDMat,PDiagMat,ScalMat} + @test iszero(@allocated det(C)) + end +end function pdtest_logdet(C::AbstractPDMat, Cmat::Matrix, verbose::Int) _pdt(verbose, "logdet") - # default tolerance in isapprox is different on 0.4. rtol argument can be deleted - # ≈ form used when 0.4 is no longer supported - @test isapprox(logdet(C), logdet(Cmat), rtol=sqrt(max(eps(real(eltype(C))), eps(real(eltype(Cmat)))))) + @test logdet(C) ≈ logdet(Cmat) + + # generic fallback in LinearAlgebra performs LU decomposition + if C isa Union{PDMat,PDiagMat,ScalMat} + @test iszero(@allocated logdet(C)) + end end