Skip to content

Commit ca8107f

Browse files
st--st--
authored
Merge branch 'master' into st/avoid_copy
2 parents 45a17ae + a0e7191 commit ca8107f

File tree

7 files changed

+146
-122
lines changed

7 files changed

+146
-122
lines changed

‎src/chol.jl‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,7 @@
1-
CholType{T,S<:AbstractMatrix} = Cholesky{T,S}
2-
31
# Accessing a.L directly might involve an extra copy();
42
# instead, always use the stored Cholesky factor:
5-
chol_lower(a::CholType) = a.uplo === 'L' ? a.L : a.U'
6-
chol_upper(a::CholType) = a.uplo === 'U' ? a.U : a.L'
3+
chol_lower(a::Cholesky) = a.uplo === 'L' ? a.L : a.U'
4+
chol_upper(a::Cholesky) = a.uplo === 'U' ? a.U : a.L'
75

86
chol_lower(a::Matrix) = chol_lower(cholesky(a))
97

‎src/deprecates.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,5 @@ using Base: @deprecate
99
@deprecate add_scal(a::Matrix, b::AbstractPDMat, c::Real) pdadd(a, b, c)
1010

1111
@deprecate full(x::AbstractPDMat) Matrix(x)
12+
13+
@deprecate CholType Cholesky

‎src/pdmat.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -4,12 +4,12 @@ Full positive definite matrix together with a Cholesky factorization object.
44
struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T}
55
dim::Int
66
mat::S
7-
chol::CholType{T,S}
7+
chol::Cholesky{T,S}
88

9-
PDMat{T,S}(d::Int,m::AbstractMatrix{T},c::CholType{T,S}) where {T,S} = new{T,S}(d,m,c)
9+
PDMat{T,S}(d::Int,m::AbstractMatrix{T},c::Cholesky{T,S}) where {T,S} = new{T,S}(d,m,c)
1010
end
1111

12-
function PDMat(mat::AbstractMatrix,chol::CholType{T,S}) where {T,S}
12+
function PDMat(mat::AbstractMatrix,chol::Cholesky{T,S}) where {T,S}
1313
d = size(mat, 1)
1414
size(chol, 1) == d ||
1515
throw(DimensionMismatch("Dimensions of mat and chol are inconsistent."))
@@ -18,7 +18,7 @@ end
1818

1919
PDMat(mat::Matrix) = PDMat(mat, cholesky(mat))
2020
PDMat(mat::Symmetric) = PDMat(Matrix(mat))
21-
PDMat(fac::CholType) = PDMat(Matrix(fac), fac)
21+
PDMat(fac::Cholesky) = PDMat(Matrix(fac), fac)
2222

2323
### Conversion
2424
Base.convert(::Type{PDMat{T}}, a::PDMat) where {T<:Real} = PDMat(convert(AbstractArray{T}, a.mat))

‎test/addition.jl‎

Lines changed: 31 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -2,34 +2,35 @@
22

33
using PDMats
44

5-
for T in [Float64,Float32]
6-
7-
printstyled("Testing addition with eltype = $T\n", color=:blue)
8-
M = convert(Array{T,2},[4. -2. -1.; -2. 5. -1.; -1. -1. 6.])
9-
V = convert(Array{T,1},[1.5, 2.5, 2.0])
10-
local X = convert(T,2.0)
11-
12-
pm1 = PDMat(M)
13-
pm2 = PDiagMat(V)
14-
pm3 = ScalMat(3,X)
15-
pm4 = X*I
16-
pm5 = PDSparseMat(sparse(M))
17-
18-
pmats = Any[pm1, pm2, pm3] #, pm5]
19-
20-
for p1 in pmats, p2 in pmats
21-
pr = p1 + p2
22-
@test size(pr) == size(p1)
23-
@test Matrix(pr) ≈ Matrix(p1) + Matrix(p2)
24-
25-
pr = pdadd(p1, p2, convert(T,1.5))
26-
@test size(pr) == size(p1)
27-
@test Matrix(pr) ≈ Matrix(p1) + Matrix(p2) * convert(T,1.5)
28-
end
29-
30-
for p1 in pmats
31-
pr = p1 + pm4
32-
@test size(pr) == size(p1)
33-
@test Matrix(pr) ≈ Matrix(p1) + pm4
34-
end
5+
@testset "addition" begin
6+
for T in (Float64, Float32)
7+
printstyled("Testing addition with eltype = $T\n"; color=:blue)
8+
M = convert(Array{T,2}, [4.0 -2.0 -1.0; -2.0 5.0 -1.0; -1.0 -1.0 6.0])
9+
V = convert(Array{T,1}, [1.5, 2.5, 2.0])
10+
local X = convert(T, 2.0)
11+
12+
pm1 = PDMat(M)
13+
pm2 = PDiagMat(V)
14+
pm3 = ScalMat(3, X)
15+
pm4 = X * I
16+
pm5 = PDSparseMat(sparse(M))
17+
18+
pmats = Any[pm1, pm2, pm3] #, pm5]
19+
20+
for p1 in pmats, p2 in pmats
21+
pr = p1 + p2
22+
@test size(pr) == size(p1)
23+
@test Matrix(pr) ≈ Matrix(p1) + Matrix(p2)
24+
25+
pr = pdadd(p1, p2, convert(T, 1.5))
26+
@test size(pr) == size(p1)
27+
@test Matrix(pr) ≈ Matrix(p1) + Matrix(p2) * convert(T, 1.5)
28+
end
29+
30+
for p1 in pmats
31+
pr = p1 + pm4
32+
@test size(pr) == size(p1)
33+
@test Matrix(pr) ≈ Matrix(p1) + pm4
34+
end
35+
end
3536
end

‎test/generics.jl‎

Lines changed: 16 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -3,22 +3,23 @@
33
using LinearAlgebra, PDMats
44
using Test
55

6-
# test scalar multiplication
7-
printstyled("Testing scalar multiplication\n", color=:blue)
8-
pm1 = PDMat(Matrix(1.0I, 3, 3))
9-
pm2 = PDiagMat(ones(3))
10-
pm3 = ScalMat(3,1)
6+
@testset "scalar multiplication" begin
7+
printstyled("Testing scalar multiplication\n"; color=:blue)
8+
pm1 = PDMat(Matrix(1.0I, 3, 3))
9+
pm2 = PDiagMat(ones(3))
10+
pm3 = ScalMat(3, 1)
1111

12-
pm1a = PDMat(Matrix(3.0I, 3, 3))
13-
pm2a = PDiagMat(3.0 .* ones(3))
14-
pm3a = ScalMat(3, 3)
12+
pm1a = PDMat(Matrix(3.0I, 3, 3))
13+
pm2a = PDiagMat(3.0 .* ones(3))
14+
pm3a = ScalMat(3, 3)
1515

16-
pmats = Any[pm1, pm2, pm3]
17-
pmatsa= Any[pm1a,pm2a,pm3a]
16+
pmats = Any[pm1, pm2, pm3]
17+
pmatsa = Any[pm1a, pm2a, pm3a]
1818

19-
for i in 1:length(pmats)
20-
@test Matrix(3.0 * pmats[i]) == Matrix(pmatsa[i])
21-
@test Matrix(pmats[i] * 3.0) == Matrix(pmatsa[i])
22-
@test Matrix(3 * pmats[i]) == Matrix(pmatsa[i])
23-
@test Matrix(pmats[i] * 3) == Matrix(pmatsa[i])
19+
for i in 1:length(pmats)
20+
@test Matrix(3.0 * pmats[i]) == Matrix(pmatsa[i])
21+
@test Matrix(pmats[i] * 3.0) == Matrix(pmatsa[i])
22+
@test Matrix(3 * pmats[i]) == Matrix(pmatsa[i])
23+
@test Matrix(pmats[i] * 3) == Matrix(pmatsa[i])
24+
end
2425
end

‎test/kron.jl‎

Lines changed: 17 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -5,31 +5,33 @@ using LinearAlgebra: LinearAlgebra
55
_randPDMat(T, n) = (X = randn(T, n, n); PDMat(X * X' + LinearAlgebra.I))
66
_randPDiagMat(T, n) = PDiagMat(rand(T, n))
77
_randScalMat(T, n) = ScalMat(n, rand(T))
8+
89
function _pd_compare(A::AbstractPDMat, B::AbstractPDMat)
910
@test dim(A) == dim(B)
1011
@test Matrix(A) ≈ Matrix(B)
1112
@test cholesky(A).L ≈ cholesky(B).L
1213
@test cholesky(A).U ≈ cholesky(B).U
13-
nothing
1414
end
15+
1516
function _pd_kron_compare(A::AbstractPDMat, B::AbstractPDMat)
16-
PDAkB1 = kron(A, B)
17-
PDAkB2 = PDMat( kron( Matrix(A), Matrix(B) ) )
18-
_pd_compare(PDAkB1, PDAkB2)
19-
nothing
17+
PDAkB_kron = kron(A, B)
18+
PDAkB_dense = PDMat( kron( Matrix(A), Matrix(B) ) )
19+
_pd_compare(PDAkB_kron, PDAkB_dense)
2020
end
2121

2222
n = 4
2323
m = 7
2424

25-
for T in [Float32, Float64]
26-
_pd_kron_compare( _randPDMat(T, n), _randPDMat(T, m) )
27-
_pd_kron_compare( _randPDiagMat(T, n), _randPDiagMat(T, m) )
28-
_pd_kron_compare( _randScalMat(T, n), _randScalMat(T, m) )
29-
_pd_kron_compare( _randPDMat(T, n), _randPDiagMat(T, m) )
30-
_pd_kron_compare( _randPDiagMat(T, m), _randPDMat(T, n) )
31-
_pd_kron_compare( _randPDMat(T, n), _randScalMat(T, m) )
32-
_pd_kron_compare( _randScalMat(T, m), _randPDMat(T, n) )
33-
_pd_kron_compare( _randPDiagMat(T, n), _randScalMat(T, m) )
34-
_pd_kron_compare( _randScalMat(T, m), _randPDiagMat(T, n) )
25+
@testset "Kronecker product" begin
26+
for T in [Float32, Float64]
27+
_pd_kron_compare( _randPDMat(T, n), _randPDMat(T, m) )
28+
_pd_kron_compare( _randPDiagMat(T, n), _randPDiagMat(T, m) )
29+
_pd_kron_compare( _randScalMat(T, n), _randScalMat(T, m) )
30+
_pd_kron_compare( _randPDMat(T, n), _randPDiagMat(T, m) )
31+
_pd_kron_compare( _randPDiagMat(T, m), _randPDMat(T, n) )
32+
_pd_kron_compare( _randPDMat(T, n), _randScalMat(T, m) )
33+
_pd_kron_compare( _randScalMat(T, m), _randPDMat(T, n) )
34+
_pd_kron_compare( _randPDiagMat(T, n), _randScalMat(T, m) )
35+
_pd_kron_compare( _randScalMat(T, m), _randPDiagMat(T, n) )
36+
end
3537
end

‎test/pdmtypes.jl‎

Lines changed: 74 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -1,58 +1,78 @@
1-
# test pd matrix types
21
using LinearAlgebra, PDMats, SparseArrays, SuiteSparse
32
using Test
43

5-
for T in [Float64, Float32]
6-
#test that all external constructors are accessible
7-
m = Matrix{T}(I, 2, 2)
8-
@test PDMat(m, cholesky(m)).mat == PDMat(Symmetric(m)).mat == PDMat(m).mat == PDMat(cholesky(m)).mat
9-
d = ones(T,2)
10-
@test PDiagMat(d,d).inv_diag == PDiagMat(d).inv_diag
11-
x = one(T)
12-
@test ScalMat(2,x,x).inv_value == ScalMat(2,x).inv_value
13-
s = SparseMatrixCSC{T}(I, 2, 2)
14-
@test PDSparseMat(s, cholesky(s)).mat == PDSparseMat(s).mat == PDSparseMat(cholesky(s)).mat
15-
16-
#test the functionality
17-
M = convert(Array{T,2}, [4. -2. -1.; -2. 5. -1.; -1. -1. 6.])
18-
V = convert(Array{T,1}, [1.5, 2.5, 2.0])
19-
X = convert(T,2.0)
20-
21-
test_pdmat(PDMat(M), M, cmat_eq=true, verbose=1) #tests of PDMat
22-
cholL = Cholesky(Matrix(transpose(cholesky(M).factors)), 'L', 0)
23-
test_pdmat(PDMat(cholL), M, cmat_eq=true, verbose=1) #tests of PDMat
24-
test_pdmat(PDiagMat(V), Matrix(Diagonal(V)), cmat_eq=true, verbose=1) #tests of PDiagMat
25-
test_pdmat(ScalMat(3,x), x*Matrix{T}(I, 3, 3), cmat_eq=true, verbose=1) #tests of ScalMat
26-
test_pdmat(PDSparseMat(sparse(M)), M, cmat_eq=true, verbose=1, t_eig=false)
27-
end
4+
@testset "pd matrix types" begin
5+
for T in [Float64, Float32]
6+
@testset "test that all external constructors are accessible" begin
7+
m = Matrix{T}(I, 2, 2)
8+
@test PDMat(m, cholesky(m)).mat == PDMat(Symmetric(m)).mat == PDMat(m).mat == PDMat(cholesky(m)).mat
9+
d = ones(T,2)
10+
@test PDiagMat(d,d).inv_diag == PDiagMat(d).inv_diag
11+
x = one(T)
12+
@test ScalMat(2,x,x).inv_value == ScalMat(2,x).inv_value
13+
s = SparseMatrixCSC{T}(I, 2, 2)
14+
@test PDSparseMat(s, cholesky(s)).mat == PDSparseMat(s).mat == PDSparseMat(cholesky(s)).mat
15+
end
16+
17+
@testset "test the functionality" begin
18+
M = convert(Array{T,2}, [4. -2. -1.; -2. 5. -1.; -1. -1. 6.])
19+
V = convert(Array{T,1}, [1.5, 2.5, 2.0])
20+
X = convert(T,2.0)
21+
22+
@testset "PDMat from Matrix" begin
23+
test_pdmat(PDMat(M), M, cmat_eq=true, verbose=1)
24+
end
25+
@testset "PDMat from Cholesky" begin
26+
cholL = Cholesky(Matrix(transpose(cholesky(M).factors)), 'L', 0)
27+
test_pdmat(PDMat(cholL), M, cmat_eq=true, verbose=1)
28+
end
29+
@testset "PDiagMat" begin
30+
test_pdmat(PDiagMat(V), Matrix(Diagonal(V)), cmat_eq=true, verbose=1)
31+
end
32+
@testset "ScalMat" begin
33+
test_pdmat(ScalMat(3,X), X*Matrix{T}(I, 3, 3), cmat_eq=true, verbose=1)
34+
end
35+
@testset "PDSparseMat" begin
36+
test_pdmat(PDSparseMat(sparse(M)), M, cmat_eq=true, verbose=1, t_eig=false)
37+
end
38+
end
39+
end
40+
41+
@testset "zero-dimensional matrices" begin
42+
Z = zeros(0, 0)
43+
test_pdmat(PDMat(Z), Z; t_eig=false)
44+
test_pdmat(PDiagMat(diag(Z)), Z; t_eig=false)
45+
end
2846

29-
m = Matrix{Float32}(I, 2, 2)
30-
@test convert(PDMat{Float64}, PDMat(m)).mat == PDMat(convert(Array{Float64}, m)).mat
31-
@test convert(AbstractArray{Float64}, PDMat(m)).mat == PDMat(convert(Array{Float64}, m)).mat
32-
m = ones(Float32,2)
33-
@test convert(PDiagMat{Float64}, PDiagMat(m)).diag == PDiagMat(convert(Array{Float64}, m)).diag
34-
@test convert(AbstractArray{Float64}, PDiagMat(m)).diag == PDiagMat(convert(Array{Float64}, m)).diag
35-
x = one(Float32); d = 4
36-
@test convert(ScalMat{Float64}, ScalMat(d, x)).value == ScalMat(d, convert(Float64, x)).value
37-
@test convert(AbstractArray{Float64}, ScalMat(d, x)).value == ScalMat(d, convert(Float64, x)).value
38-
s = SparseMatrixCSC{Float32}(I, 2, 2)
39-
@test convert(PDSparseMat{Float64}, PDSparseMat(s)).mat == PDSparseMat(convert(SparseMatrixCSC{Float64}, s)).mat
40-
41-
Z = zeros(0, 0)
42-
test_pdmat(PDMat(Z), Z; t_eig=false)
43-
test_pdmat(PDiagMat(diag(Z)), Z; t_eig=false)
44-
45-
# no-op conversion with correct eltype (#101)
46-
X = PDMat((Y->Y'Y)(randn(Float32, 4, 4)))
47-
@test convert(AbstractArray{Float32}, X) === X
48-
@test convert(AbstractArray{Float64}, X) !== X
49-
50-
# type stability of whiten! and unwhiten!
51-
a = PDMat([1 0.5; 0.5 1])
52-
@inferred whiten!(ones(2), a, ones(2))
53-
@inferred unwhiten!(ones(2), a, ones(2))
54-
@inferred whiten(a, ones(2))
55-
@inferred unwhiten(a, ones(2))
56-
57-
# convert Matrix type to the same Cholesky type (#117)
58-
@test PDMat([1 0; 0 1]) == [1.0 0.0; 0.0 1.0]
47+
@testset "float type conversions" begin
48+
m = Matrix{Float32}(I, 2, 2)
49+
@test convert(PDMat{Float64}, PDMat(m)).mat == PDMat(convert(Array{Float64}, m)).mat
50+
@test convert(AbstractArray{Float64}, PDMat(m)).mat == PDMat(convert(Array{Float64}, m)).mat
51+
m = ones(Float32,2)
52+
@test convert(PDiagMat{Float64}, PDiagMat(m)).diag == PDiagMat(convert(Array{Float64}, m)).diag
53+
@test convert(AbstractArray{Float64}, PDiagMat(m)).diag == PDiagMat(convert(Array{Float64}, m)).diag
54+
x = one(Float32); d = 4
55+
@test convert(ScalMat{Float64}, ScalMat(d, x)).value == ScalMat(d, convert(Float64, x)).value
56+
@test convert(AbstractArray{Float64}, ScalMat(d, x)).value == ScalMat(d, convert(Float64, x)).value
57+
s = SparseMatrixCSC{Float32}(I, 2, 2)
58+
@test convert(PDSparseMat{Float64}, PDSparseMat(s)).mat == PDSparseMat(convert(SparseMatrixCSC{Float64}, s)).mat
59+
end
60+
61+
@testset "no-op conversion with correct eltype (#101)" begin
62+
X = PDMat((Y->Y'Y)(randn(Float32, 4, 4)))
63+
@test convert(AbstractArray{Float32}, X) === X
64+
@test convert(AbstractArray{Float64}, X) !== X
65+
end
66+
67+
@testset "type stability of whiten! and unwhiten!" begin
68+
a = PDMat([1 0.5; 0.5 1])
69+
@inferred whiten!(ones(2), a, ones(2))
70+
@inferred unwhiten!(ones(2), a, ones(2))
71+
@inferred whiten(a, ones(2))
72+
@inferred unwhiten(a, ones(2))
73+
end
74+
75+
@testset "convert Matrix type to the same Cholesky type (#117)" begin
76+
@test PDMat([1 0; 0 1]) == [1.0 0.0; 0.0 1.0]
77+
end
78+
end

0 commit comments

Comments
 (0)