diff --git a/Project.toml b/Project.toml index f24757b..9302756 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "PDMats" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.37" +version = "0.11.38" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/pdsparsemat.jl b/src/pdsparsemat.jl index de0638f..c315f27 100644 --- a/src/pdsparsemat.jl +++ b/src/pdsparsemat.jl @@ -1,23 +1,39 @@ +const AbstractSparseLike{T<:Real, S<:AbstractSparseMatrix{T}} = Union{S, Symmetric{T,S}, Hermitian{T,S}} +# We can't define `convert(AbstractSparseLike{T}, mat)` here because of piracy +convert_abstractsparselike(::Type{T}, mat::AbstractSparseMatrix) where T<:Real = convert(AbstractSparseMatrix{T}, mat) +convert_abstractsparselike(::Type{T}, mat::SparseMatrixCSC) where T<:Real = convert(SparseMatrixCSC{T}, mat) # the AbstractSparseMatrix method doesn't work +convert_abstractsparselike(::Type{T}, mat::Symmetric{S,<:AbstractSparseMatrix{S}} where S) where T<:Real = Symmetric(convert_abstractsparselike(T, mat.data), Symbol(mat.uplo)) +convert_abstractsparselike(::Type{T}, mat::Hermitian{S,<:AbstractSparseMatrix{S}} where S) where T<:Real = Hermitian(convert_abstractsparselike(T, mat.data), Symbol(mat.uplo)) + """ Sparse positive definite matrix together with a Cholesky factorization object. """ -struct PDSparseMat{T <: Real, S <: AbstractSparseMatrix} <: AbstractPDMat{T} +struct PDSparseMat{T <: Real, S <: AbstractSparseLike{T}, C <: CholTypeSparse} <: AbstractPDMat{T} mat::S - chol::CholTypeSparse - - PDSparseMat{T, S}(m::AbstractSparseMatrix{T}, c::CholTypeSparse) where {T, S} = - new{T, S}(m, c) #add {T} to CholTypeSparse argument once #14076 is implemented + chol::C + # Note: cholesky(::SparseMatrixCSC{Float32}) returns CholTypeSparse{Float64}, so we can't enforce C <: CholTypeSparse{T} + + function PDSparseMat{T, S, C}(m::AbstractSparseLike{T}, c::CholTypeSparse) where {T, S, C} + d = LinearAlgebra.checksquare(m) + size(c, 1) == d || + throw(DimensionMismatch("Dimensions of mat and chol are inconsistent.")) + return new{T, S, C}(m, c) + end end -@deprecate PDSparseMat{T, S}(d::Int, m::AbstractSparseMatrix{T}, c::CholTypeSparse) where {T, S} PDSparseMat{T, S}(m, c) +@deprecate PDSparseMat{T, S}(d::Int, m::AbstractSparseMatrix{T}, c::CholTypeSparse) where {T, S} PDSparseMat{T, S, typeof(c)}(m, c) -function PDSparseMat(mat::AbstractSparseMatrix, chol::CholTypeSparse) - d = LinearAlgebra.checksquare(mat) - size(chol, 1) == d || - throw(DimensionMismatch("Dimensions of mat and chol are inconsistent.")) - return PDSparseMat{eltype(mat), typeof(mat)}(mat, chol) +PDSparseMat{T, S}(mat::AbstractSparseLike, chol::CholTypeSparse) where {T <: Real, S <: AbstractSparseLike{T}} = PDSparseMat{T, S, typeof(chol)}(mat, chol) +PDSparseMat{T}(mat::AbstractSparseLike{T}, chol::CholTypeSparse) where {T <: Real} = PDSparseMat{T, typeof(mat)}(mat, chol) +function PDSparseMat{T}(mat::AbstractSparseLike, chol::CholTypeSparse) where {T <: Real} + mat = convert_abstractsparselike(T, mat)::AbstractSparseLike{T} # prevent StackOverflowError if the eltype conversion lies + return PDSparseMat{T}(mat, chol) end +PDSparseMat(mat::AbstractSparseLike{T}, chol::CholTypeSparse) where {T <: Real} = PDSparseMat{T}(mat, chol) + +PDSparseMat{T}(mat::AbstractSparseLike{T,S}) where {T<:Real,S<:SparseMatrixCSC{T}} = PDSparseMat{T}(mat, cholesky(mat)) +PDSparseMat{T}(mat::AbstractSparseLike) where T<:Real = PDSparseMat{T}(convert_abstractsparselike(T, mat)) -PDSparseMat(mat::SparseMatrixCSC) = PDSparseMat(mat, cholesky(mat)) +PDSparseMat(mat::AbstractSparseLike{T,S}) where {T<:Real,S<:SparseMatrixCSC{T}} = PDSparseMat{T}(mat) PDSparseMat(fac::CholTypeSparse) = PDSparseMat(sparse(fac), fac) function Base.getproperty(a::PDSparseMat, s::Symbol) diff --git a/src/utils.jl b/src/utils.jl index f589815..7ba161d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -7,7 +7,7 @@ macro check_argdims(cond) end end -function _addscal!(r::Matrix, a::Matrix, b::Union{Matrix, SparseMatrixCSC}, c::Real) +function _addscal!(r::Matrix, a::Matrix, b::Union{M, Symmetric{T, <:M}, Hermitian{T, <:M}}, c::Real) where {T<:Real, M <: Union{Matrix, SparseMatrixCSC}} if c == one(c) for i in eachindex(a) @inbounds r[i] = a[i] + b[i] diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index 6420467..f01973a 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -16,6 +16,8 @@ using Test @testset "test the functionality" begin M = convert(Array{T, 2}, [4.0 -2.0 -1.0; -2.0 5.0 -1.0; -1.0 -1.0 6.0]) + MS = Symmetric(convert(Matrix{T}, [4.0 NaN NaN; -2.0 5.0 NaN; -1.0 -1.0 6.0]), :L) + MH = Hermitian(convert(Matrix{T}, [4.0 NaN NaN; -2.0 5.0 NaN; -1.0 -1.0 6.0]), :L) V = convert(Array{T, 1}, [1.5, 2.5, 2.0]) X = convert(T, 2.0) f64M = Float64.(M) @@ -25,6 +27,8 @@ using Test test_pdmat(PDMat(M), M, cmat_eq = true, verbose = 1) test_pdmat(PDMat{Float64}(M), f64M, cmat_eq = true, verbose = 1) test_pdmat(PDMat{Float64, Matrix{Float64}}(M), f64M, cmat_eq = true, verbose = 1) + test_pdmat(PDMat(MS), M, cmat_eq = true, verbose = 1) + test_pdmat(PDMat(MH), M, cmat_eq = true, verbose = 1) @test_throws TypeError PDMat{Float32, Matrix{Float64}}(M) end @testset "PDMat from PDMat" begin @@ -57,6 +61,16 @@ using Test end @testset "PDSparseMat" begin test_pdmat(PDSparseMat(sparse(M)), M, cmat_eq = true, verbose = 1, t_eig = false) + test_pdmat(PDSparseMat{Float64}(sparse(M)), Float64.(M), cmat_eq = true, verbose = 1, t_eig = false) + test_pdmat(PDSparseMat{Float64}(sparse(M), cholesky(sparse(M))), Float64.(M), cmat_eq = true, verbose = 1, t_eig = false) + sparseMS = Symmetric(sparse(MS.data), :L) + test_pdmat(PDSparseMat(sparseMS), M, cmat_eq = true, verbose = 1, t_eig = false) + test_pdmat(PDSparseMat{Float64}(sparseMS), Float64.(M), cmat_eq = true, verbose = 1, t_eig = false) + test_pdmat(PDSparseMat{Float64}(sparseMS, cholesky(sparseMS)), Float64.(M), cmat_eq = true, verbose = 1, t_eig = false) + sparseMH = Hermitian(sparse(MH.data), :L) + test_pdmat(PDSparseMat(sparseMH), M, cmat_eq = true, verbose = 1, t_eig = false) + test_pdmat(PDSparseMat{Float64}(sparseMH), Float64.(M), cmat_eq = true, verbose = 1, t_eig = false) + test_pdmat(PDSparseMat{Float64}(sparseMH, cholesky(sparseMH)), Float64.(M), cmat_eq = true, verbose = 1, t_eig = false) end end @@ -222,7 +236,9 @@ using Test @test Mat32 isa Matrix{Float32} @test Mat32 ≈ Float32.(A) - M = @inferred AbstractPDMat(cholesky(sparse(A))) + # sparse(::SparseArrays.CHOLMOD.Factor) is inferred to return a Union{Symmetric{T, <:SparseMatrixCSC{T}}, SparseMatrixCSC{T}} + # Now that we support Symmetric directly, the following is no longer inferrable + M = AbstractPDMat(cholesky(sparse(A))) @test M isa PDSparseMat @test Matrix(M) ≈ A end