diff --git a/Project.toml b/Project.toml index 5dc5781..9754928 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" [compat] +BandedMatrices = "0.15, 1" +Random = "<0.0.1, 1" +StaticArrays = "1" +Test = "<0.0.1, 1" julia = "1" [extensions] @@ -15,12 +19,13 @@ PDMatsSparseArraysExt = "SparseArrays" [extras] BandedMatrices = "aae01518-5342-5314-be14-df237901396f" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["BandedMatrices", "SparseArrays", "StaticArrays", "Test"] +test = ["BandedMatrices", "Random", "SparseArrays", "StaticArrays", "Test"] [weakdeps] SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl index e3ae261..e542f0e 100644 --- a/src/pdiagmat.jl +++ b/src/pdiagmat.jl @@ -28,7 +28,7 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDiagMat) where {T<:Real} = convert(PD ### Basics Base.size(a::PDiagMat) = (a.dim, a.dim) -Base.Matrix(a::PDiagMat) = Matrix(Diagonal(a.diag)) +Base.Matrix{T}(a::PDiagMat) where {T} = Matrix{T}(Diagonal(a.diag)) LinearAlgebra.diag(a::PDiagMat) = copy(a.diag) LinearAlgebra.cholesky(a::PDiagMat) = Cholesky(Diagonal(map(sqrt, a.diag)), 'U', 0) @@ -37,11 +37,7 @@ Base.broadcastable(a::PDiagMat) = Base.broadcastable(Diagonal(a.diag)) ### Inheriting from AbstractMatrix -function Base.getindex(a::PDiagMat, i::Integer) - ncol, nrow = fldmod1(i, a.dim) - ncol == nrow ? a.diag[nrow] : zero(eltype(a)) -end -Base.getindex(a::PDiagMat{T},i::Integer,j::Integer) where {T} = i == j ? a.diag[i] : zero(T) +Base.@propagate_inbounds Base.getindex(a::PDiagMat{T}, i::Int, j::Int) where {T} = i == j ? a.diag[i] : zero(T) ### Arithmetics diff --git a/src/pdmat.jl b/src/pdmat.jl index 2e760fe..61e7681 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -46,7 +46,7 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat ### Basics Base.size(a::PDMat) = (a.dim, a.dim) -Base.Matrix(a::PDMat) = Matrix(a.mat) +Base.Matrix{T}(a::PDMat) where {T} = Matrix{T}(a.mat) LinearAlgebra.diag(a::PDMat) = diag(a.mat) LinearAlgebra.cholesky(a::PDMatCholesky) = a.fact @@ -55,8 +55,11 @@ Base.broadcastable(a::PDMat) = Base.broadcastable(a.mat) ### Inheriting from AbstractMatrix -Base.getindex(a::PDMat, i::Int) = getindex(a.mat, i) -Base.getindex(a::PDMat, I::Vararg{Int, N}) where {N} = getindex(a.mat, I...) +Base.IndexStyle(::Type{PDMat{T,S}}) where {T,S} = Base.IndexStyle(S) +# Linear Indexing +Base.@propagate_inbounds Base.getindex(a::PDMat, i::Int) = getindex(a.mat, i) +# Cartesian Indexing +Base.@propagate_inbounds Base.getindex(a::PDMat, I::Vararg{Int, 2}) = getindex(a.mat, I...) ### Arithmetics @@ -84,8 +87,8 @@ end Base.inv(a::PDMat) = PDMat(inv(a.fact)) LinearAlgebra.det(a::PDMat) = det(a.fact) LinearAlgebra.logdet(a::PDMat) = logdet(a.fact) -LinearAlgebra.eigmax(a::PDMat) = eigmax(a.mat) -LinearAlgebra.eigmin(a::PDMat) = eigmin(a.mat) +LinearAlgebra.eigmax(a::PDMat) = eigmax(Symmetric(a.mat)) +LinearAlgebra.eigmin(a::PDMat) = eigmin(Symmetric(a.mat)) Base.sqrt(A::PDMat) = PDMat(sqrt(Hermitian(A.mat))) function Base.kron(A::PDMatCholesky, B::PDMatCholesky) diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index 8753c25..9331671 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -170,21 +170,30 @@ using Test @testset "AbstractPDMat constructors (#136)" begin x = rand(10, 10) - A = x' * x + I + A = Array(Symmetric(x' * x + I)) M = @inferred AbstractPDMat(A) @test M isa PDMat @test cholesky(M) isa Cholesky @test Matrix(M) ≈ A + Mat32 = @inferred Matrix{Float32}(M) + @test eltype(Mat32) == Float32 + @test Mat32 ≈ Float32.(A) M = @inferred AbstractPDMat(cholesky(A)) @test M isa PDMat @test cholesky(M) isa Cholesky @test Matrix(M) ≈ A + Mat32 = @inferred Matrix{Float32}(M) + @test Mat32 isa Matrix{Float32} + @test Mat32 ≈ Float32.(A) M = @inferred AbstractPDMat(Diagonal(A)) @test M isa PDiagMat @test Matrix(M) ≈ Diagonal(A) + Mat32 = @inferred Matrix{Float32}(M) + @test Mat32 isa Matrix{Float32} + @test Mat32 ≈ Float32.(Diagonal(A)) M = @inferred AbstractPDMat(Symmetric(Diagonal(A))) @test M isa PDiagMat @@ -198,6 +207,9 @@ using Test @test M isa PDMat @test cholesky(M) isa CHOLMOD.Factor @test Matrix(M) ≈ A + Mat32 = @inferred Matrix{Float32}(M) + @test Mat32 isa Matrix{Float32} + @test Mat32 ≈ Float32.(A) if VERSION < v"1.6" # inference fails e.g. on Julia 1.0 @@ -213,7 +225,7 @@ using Test @testset "properties and fields" begin for dim in (1, 5, 10) x = rand(dim, dim) - M = PDMat(x' * x + I) + M = PDMat(Array(Symmetric(x' * x + I))) @test fieldnames(typeof(M)) == (:mat, :fact) @test propertynames(M) == (fieldnames(typeof(M))..., :dim) @test getproperty(M, :dim) === dim @@ -238,7 +250,7 @@ using Test if HAVE_CHOLMOD x = sprand(dim, dim, 0.2) - M = PDMat(x' * x + I) + M = PDMat(sparse(Symmetric(x' * x + I))) @test fieldnames(typeof(M)) == (:mat, :fact) @test propertynames(M) == (fieldnames(typeof(M))..., :dim) @test getproperty(M, :dim) === dim @@ -251,14 +263,14 @@ using Test @testset "Incorrect dimensions" begin x = rand(10, 10) - A = x * x' + I + A = Array(Symmetric(x * x' + I)) C = cholesky(A) @test_throws DimensionMismatch PDMat(A[:, 1:(end - 1)], C) @test_throws DimensionMismatch PDMat(A[1:(end - 1), 1:(end - 1)], C) if HAVE_CHOLMOD x = sprand(10, 10, 0.2) - A = x * x' + I + A = sparse(Symmetric(x * x' + I)) C = cholesky(A) @test_throws DimensionMismatch PDMat(A[:, 1:(end - 1)], C) @test_throws DimensionMismatch PDMat(A[1:(end - 1), 1:(end - 1)], C) @@ -269,7 +281,7 @@ using Test # This falls back to the generic method in Julia based on broadcasting dim = 4 x = rand(dim, dim) - A = PDMat(x' * x + I) + A = PDMat(Array(Symmetric(x' * x + I))) @test Base.broadcastable(A) == A.mat B = PDiagMat(rand(dim)) diff --git a/test/specialarrays.jl b/test/specialarrays.jl index 427f9b4..d2e1c2f 100644 --- a/test/specialarrays.jl +++ b/test/specialarrays.jl @@ -4,7 +4,7 @@ using StaticArrays @testset "Special matrix types" begin @testset "StaticArrays" begin # Full matrix - S = (x -> x * x' + I)(@SMatrix(randn(4, 7))) + S = (x -> SMatrix{4,4}(Symmetric(x * x' + I)))(@SMatrix(randn(4, 7))) PDS = PDMat(S) @test PDS isa PDMat{Float64, <:SMatrix{4, 4, Float64}} @test isbits(PDS) diff --git a/test/testutils.jl b/test/testutils.jl index 1163e39..e1abbb2 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -4,7 +4,9 @@ # the implementation of a subtype of AbstractPDMat # -using PDMats, LinearAlgebra, SparseArrays, Test +using PDMats, LinearAlgebra, SparseArrays, Test, Random + +Random.seed!(10) if isdefined(Base, :get_extension) const HAVE_CHOLMOD = isdefined(SparseArrays, :CHOLMOD)