Merge branch 'master' into dw/unify_pdmat_pdsparsemat

Project.toml

@@ -8,19 +8,24 @@ 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]
 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"

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
 

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)

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))

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)

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)