From 2f6fd1b3c7ef71f3a91b1cce9d9e92724a5c890b Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 8 Dec 2023 10:50:01 +0530 Subject: [PATCH] Typed Matrix constructor and more stable tests --- Project.toml | 7 ++++++- src/pdiagmat.jl | 7 ++----- src/pdmat.jl | 4 ++-- src/pdsparsemat.jl | 4 ++-- test/pdmtypes.jl | 24 ++++++++++++++++++------ test/specialarrays.jl | 2 +- test/testutils.jl | 4 +++- 7 files changed, 34 insertions(+), 18 deletions(-) diff --git a/Project.toml b/Project.toml index d9b4481..f708168 100644 --- a/Project.toml +++ b/Project.toml @@ -8,12 +8,17 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" [compat] +BandedMatrices = "1" +Random = "1" +StaticArrays = "1" +Test = "1" julia = "1" [extras] BandedMatrices = "aae01518-5342-5314-be14-df237901396f" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["BandedMatrices", "StaticArrays", "Test"] +test = ["BandedMatrices", "StaticArrays", "Random", "Test"] diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl index cb43a70..fad2b7e 100644 --- a/src/pdiagmat.jl +++ b/src/pdiagmat.jl @@ -29,6 +29,7 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDiagMat) where {T<:Real} = convert(PD 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 +38,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 abaec69..fc3dca4 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -43,6 +43,7 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat 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::PDMat) = a.chol @@ -51,8 +52,7 @@ 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.@propagate_inbounds Base.getindex(a::PDMat, I::Vararg{Int, 2}) = getindex(a.mat, I...) ### Arithmetics diff --git a/src/pdsparsemat.jl b/src/pdsparsemat.jl index 032e0ed..b1b602e 100644 --- a/src/pdsparsemat.jl +++ b/src/pdsparsemat.jl @@ -46,13 +46,13 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDSparseMat) where {T<:Real} = convert Base.size(a::PDSparseMat) = (a.dim, a.dim) Base.Matrix(a::PDSparseMat) = Matrix(a.mat) +Base.Matrix{T}(a::PDSparseMat) where {T} = Matrix{T}(a.mat) LinearAlgebra.diag(a::PDSparseMat) = diag(a.mat) LinearAlgebra.cholesky(a::PDSparseMat) = a.chol ### Inheriting from AbstractMatrix -Base.getindex(a::PDSparseMat,i::Integer) = getindex(a.mat, i) -Base.getindex(a::PDSparseMat,I::Vararg{Int, N}) where {N} = getindex(a.mat, I...) +Base.getindex(a::PDSparseMat, I::Vararg{Int, 2}) = getindex(a.mat, I...) ### Arithmetics diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index c749cf8..99ce6ce 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -166,19 +166,28 @@ using Test @testset "AbstractPDMat constructors (#136)" begin x = rand(10, 10) - A = x' * x + I + A = Array(Hermitian(x' * x + I)) M = @inferred AbstractPDMat(A) @test M isa PDMat @test Matrix(M) ≈ A + Mat32 = Matrix{Float32}(M) + @test eltype(Mat32) == Float32 + @test Mat32 ≈ Float32.(A) M = @inferred AbstractPDMat(cholesky(A)) @test M isa PDMat @test Matrix(M) ≈ A + Mat32 = Matrix{Float32}(M) + @test eltype(Mat32) == Float32 + @test Mat32 ≈ Float32.(A) M = @inferred AbstractPDMat(Diagonal(A)) @test M isa PDiagMat @test Matrix(M) ≈ Diagonal(A) + Mat32 = Matrix{Float32}(M) + @test eltype(Mat32) == Float32 + @test Mat32 ≈ Float32.(Diagonal(A)) M = @inferred AbstractPDMat(Symmetric(Diagonal(A))) @test M isa PDiagMat @@ -191,6 +200,9 @@ using Test M = @inferred AbstractPDMat(sparse(A)) @test M isa PDSparseMat @test Matrix(M) ≈ A + Mat32 = Matrix{Float32}(M) + @test eltype(Mat32) == Float32 + @test Mat32 ≈ Float32.(A) if VERSION < v"1.6" # inference fails e.g. on Julia 1.0 @@ -205,7 +217,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(Hermitian(x' * x + I))) @test fieldnames(typeof(M)) == (:mat, :chol) @test propertynames(M) == (fieldnames(typeof(M))..., :dim) @test getproperty(M, :dim) === dim @@ -230,7 +242,7 @@ using Test if HAVE_CHOLMOD x = sprand(dim, dim, 0.2) - M = PDSparseMat(x' * x + I) + M = PDSparseMat(sparse(Hermitian(x' * x + I))) @test fieldnames(typeof(M)) == (:mat, :chol) @test propertynames(M) == (fieldnames(typeof(M))..., :dim) @test getproperty(M, :dim) === dim @@ -243,14 +255,14 @@ using Test @testset "Incorrect dimensions" begin x = rand(10, 10) - A = x * x' + I + A = Array(Hermitian(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(Hermitian(x * x' + I)) C = cholesky(A) @test_throws DimensionMismatch PDSparseMat(A[:, 1:(end - 1)], C) @test_throws DimensionMismatch PDSparseMat(A[1:(end - 1), 1:(end - 1)], C) @@ -261,7 +273,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(Hermitian(x' * x + I))) @test Base.broadcastable(A) == A.mat B = PDiagMat(rand(dim)) diff --git a/test/specialarrays.jl b/test/specialarrays.jl index be6cbe8..8e67d77 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}(Hermitian(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 32bd5a9..e3a6e1e 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -4,7 +4,9 @@ # the implementation of a subtype of AbstractPDMat # -using PDMats, SuiteSparse, Test +using PDMats, SuiteSparse, Test, Random + +Random.seed!(10) const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD) const PDMatType = HAVE_CHOLMOD ? Union{PDMat, PDSparseMat, PDiagMat} : Union{PDMat, PDiagMat}