Switch to native julia 0.7 constructions

src/PDMats.jl

@@ -6,7 +6,7 @@ module PDMats
     using LinearAlgebra, SparseArrays
 
     import Base: +, *, \, /, ==, convert, inv, Matrix
-    import LinearAlgebra: logdet, diag, diagm, eigmax, eigmin
+    import LinearAlgebra: logdet, diag, eigmax, eigmin
 
     export
         # Types

src/generics.jl

@@ -40,10 +40,10 @@ unwhiten(a::AbstractPDMat, x::StridedVecOrMat) = unwhiten!(similar(x), a, x)
 
 function quad(a::AbstractPDMat{T}, x::StridedMatrix{S}) where {T<:Real, S<:Real}
     @check_argdims dim(a) == size(x, 1)
-    quad!(Array{promote_type(T, S)}(size(x,2)), a, x)
+    quad!(Array{promote_type(T, S)}(uninitialized, size(x,2)), a, x)
 end
 
 function invquad(a::AbstractPDMat{T}, x::StridedMatrix{S}) where {T<:Real, S<:Real}
     @check_argdims dim(a) == size(x, 1)
-    invquad!(Array{promote_type(T, S)}(size(x,2)), a, x)
+    invquad!(Array{promote_type(T, S)}(uninitialized, size(x,2)), a, x)
 end

src/pdiagmat.jl

@@ -21,7 +21,7 @@ Base.convert(::Type{AbstractArray{T}}, a::PDiagMat) where {T<:Real} = convert(PD
 ### Basics
 
 dim(a::PDiagMat) = a.dim
-Base.Matrix(a::PDiagMat) = diagm(a.diag)
+Base.Matrix(a::PDiagMat) = Matrix(Diagonal(a.diag))
 LinearAlgebra.diag(a::PDiagMat) = copy(a.diag)
 
 
@@ -32,7 +32,7 @@ function pdadd!(r::Matrix, a::Matrix, b::PDiagMat, c)
     if r === a
         _adddiag!(r, b.diag, c)
     else
-        _adddiag!(copy!(r, a), b.diag, c)
+        _adddiag!(copyto!(r, a), b.diag, c)
     end
     return r
 end
@@ -84,28 +84,53 @@ unwhiten!(r::StridedMatrix, a::PDiagMat, x::StridedMatrix) =
 quad(a::PDiagMat, x::StridedVector) = wsumsq(a.diag, x)
 invquad(a::PDiagMat, x::StridedVector) = wsumsq(a.inv_diag, x)
 
-quad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix) = At_mul_B!(r, abs2.(x), a.diag)
-invquad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix) = At_mul_B!(r, abs2.(x), a.inv_diag)
+function quad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix)
+    m, n = size(x)
+    ad = a.diag
+    @check_argdims m == length(ad) && length(r) == n
+    @inbounds for j = 1:n
+        s = zero(promote_type(eltype(ad), eltype(x)))
+        for i in 1:m
+            s += ad[i] * abs2(x[i,j])
+        end
+        r[j] = s
+    end
+    r
+end
+
+function invquad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix)
+    m, n = size(x)
+    ainvd = a.inv_diag
+    @check_argdims m == length(ainvd) && length(r) == n
+    @inbounds for j = 1:n
+        s = zero(promote_type(eltype(ainvd), eltype(x)))
+        for i in 1:m
+            s += ainvd[i] * abs2(x[i,j])
+        end
+        r[j] = s
+    end
+    r
+end
 
 
 ### tri products
 
 function X_A_Xt(a::PDiagMat, x::StridedMatrix)
     z = x .* reshape(sqrt.(a.diag), 1, dim(a))
-    A_mul_Bt(z, z)
+    z * transpose(z)
 end
 
 function Xt_A_X(a::PDiagMat, x::StridedMatrix)
     z = x .* sqrt.(a.diag)
-    At_mul_B(z, z)
+    transpose(z) * z
 end
 
 function X_invA_Xt(a::PDiagMat, x::StridedMatrix)
     z = x .* reshape(sqrt.(a.inv_diag), 1, dim(a))
-    A_mul_Bt(z, z)
+    z * transpose(z)
 end
 
 function Xt_invA_X(a::PDiagMat, x::StridedMatrix)
     z = x .* sqrt.(a.inv_diag)
-    At_mul_B(z, z)
+    transpose(z) * z
 end

src/pdmat.jl

@@ -15,7 +15,7 @@ end
 
 PDMat(mat::Matrix) = PDMat(mat, cholfact(mat))
 PDMat(mat::Symmetric) = PDMat(Matrix(mat))
-PDMat(fac::CholType) = PDMat(Matrix(fac),fac)
+PDMat(fac::CholType) = PDMat(Matrix(fac), fac)
 
 ### Conversion
 Base.convert(::Type{PDMat{T}},         a::PDMat) where {T<:Real} = PDMat(convert(AbstractArray{T}, a.mat))
@@ -51,15 +51,13 @@ LinearAlgebra.eigmin(a::PDMat) = eigmin(a.mat)
 ### whiten and unwhiten
 
 function whiten!(r::StridedVecOrMat, a::PDMat, x::StridedVecOrMat)
-    cf = a.chol[:UL]
-    istriu(cf) ? Ac_ldiv_B!(cf, _rcopy!(r, x)) : A_ldiv_B!(cf, _rcopy!(r, x))
-    return r
+    cf = a.chol.UL
+    ldiv!(istriu(cf) ? transpose(cf) : cf, _rcopy!(r, x))
 end
 
 function unwhiten!(r::StridedVecOrMat, a::PDMat, x::StridedVecOrMat)
-    cf = a.chol[:UL]
-    istriu(cf) ? Ac_mul_B!(cf, _rcopy!(r, x)) : A_mul_B!(cf, _rcopy!(r, x))
-    return r
+    cf = a.chol.UL
+    lmul!(istriu(cf) ? transpose(cf) : cf, _rcopy!(r, x))
 end
 
 
@@ -76,28 +74,25 @@ invquad!(r::AbstractArray, a::PDMat, x::StridedMatrix) = colwise_dot!(r, x, a.ma
 
 function X_A_Xt(a::PDMat, x::StridedMatrix)
     z = copy(x)
-    cf = a.chol[:UL]
-    istriu(cf) ? A_mul_Bc!(z, cf) : A_mul_B!(z, cf)
-    A_mul_Bt(z, z)
+    cf = a.chol.UL
+    rmul!(z, istriu(cf) ? transpose(cf) : cf)
+    z * transpose(z)
 end
 
 function Xt_A_X(a::PDMat, x::StridedMatrix)
-    z = copy(x)
-    cf = a.chol[:UL]
-    istriu(cf) ? A_mul_B!(cf, z) : Ac_mul_B!(cf, z)
-    At_mul_B(z, z)
+    cf = a.chol.UL
+    z = lmul!(istriu(cf) ? cf : transpose(cf), copy(x))
+    transpose(z) * z
 end
 
 function X_invA_Xt(a::PDMat, x::StridedMatrix)
-    z = copy(x)
-    cf = a.chol[:UL]
-    istriu(cf) ? A_rdiv_B!(z, cf) : A_rdiv_Bc!(z, cf)
-    A_mul_Bt(z, z)
+    cf = a.chol.UL
+    z = rdiv!(copy(x), istriu(cf) ? cf : transpose(cf))
+    z * transpose(z)
 end
 
 function Xt_invA_X(a::PDMat, x::StridedMatrix)
-    z = copy(x)
-    cf = a.chol[:UL]
-    istriu(cf) ? Ac_ldiv_B!(cf, z) : A_ldiv_B!(cf, z)
-    At_mul_B(z, z)
+    cf = a.chol.UL
+    z = ldiv!(istriu(cf) ? transpose(cf) : cf, copy(x))
+    transpose(z) * z
 end

src/scalmat.jl

@@ -15,8 +15,8 @@ convert(::Type{AbstractArray{T}}, a::ScalMat) where {T<:Real} = convert(ScalMat{
 ### Basics
 
 dim(a::ScalMat) = a.dim
-full(a::ScalMat) = diagm(fill(a.value, a.dim))
-diag(a::ScalMat) = fill(a.value, a.dim)
+Base.Matrix(a::ScalMat) = Matrix(Diagonal(fill(a.value, a.dim)))
+LinearAlgebra.diag(a::ScalMat) = fill(a.value, a.dim)
 
 
 ### Arithmetics
@@ -26,7 +26,7 @@ function pdadd!(r::Matrix, a::Matrix, b::ScalMat, c)
     if r === a
         _adddiag!(r, b.value * c)
     else
-        _adddiag!(copy!(r, a), b.value * c)
+        _adddiag!(copyto!(r, a), b.value * c)
     end
     return r
 end
@@ -39,30 +39,22 @@ end
 
 ### Algebra
 
-inv(a::ScalMat) = ScalMat(a.dim, a.inv_value, a.value)
-logdet(a::ScalMat) = a.dim * log(a.value)
-eigmax(a::ScalMat) = a.value
-eigmin(a::ScalMat) = a.value
+Base.inv(a::ScalMat) = ScalMat(a.dim, a.inv_value, a.value)
+LinearAlgebra.logdet(a::ScalMat) = a.dim * log(a.value)
+LinearAlgebra.eigmax(a::ScalMat) = a.value
+LinearAlgebra.eigmin(a::ScalMat) = a.value
 
 
 ### whiten and unwhiten
 
 function whiten!(r::StridedVecOrMat, a::ScalMat, x::StridedVecOrMat)
     @check_argdims dim(a) == size(x, 1)
-    c = sqrt(a.inv_value)
-    for i = 1:length(x)
-        @inbounds r[i] = x[i] * c
-    end
-    return r
+    mul!(r, x, sqrt(a.inv_value))
 end
 
 function unwhiten!(r::StridedVecOrMat, a::ScalMat, x::StridedVecOrMat)
     @check_argdims dim(a) == size(x, 1)
-    c = sqrt(a.value)
-    for i = 1:length(x)
-        @inbounds r[i] = x[i] * c
-    end
-    return r
+    mul!(r, x, sqrt(a.value))
 end
 
 
@@ -79,20 +71,20 @@ invquad!(r::AbstractArray, a::ScalMat, x::StridedMatrix) = colwise_sumsq!(r, x,
 
 function X_A_Xt(a::ScalMat, x::StridedMatrix)
     @check_argdims dim(a) == size(x, 2)
-    gemm('N', 'T', a.value, x, x)
+    lmul!(a.value, x * transpose(x))
 end
 
 function Xt_A_X(a::ScalMat, x::StridedMatrix)
     @check_argdims dim(a) == size(x, 1)
-    gemm('T', 'N', a.value, x, x)
+    lmul!(a.value, transpose(x) * x)
 end
 
 function X_invA_Xt(a::ScalMat, x::StridedMatrix)
     @check_argdims dim(a) == size(x, 2)
-    gemm('N', 'T', a.inv_value, x, x)
+    lmul!(a.inv_value, x * transpose(x))
 end
 
 function Xt_invA_X(a::ScalMat, x::StridedMatrix)
     @check_argdims dim(a) == size(x, 1)
-    gemm('T', 'N', a.inv_value, x, x)
+    lmul!(a.inv_value, transpose(x) * x)
 end

src/testutils.jl

@@ -4,8 +4,7 @@
 #       the implementation of a subtype of AbstractPDMat
 #
 
-import Compat.Test: @test
-using Compat: view
+using Test: @test
 
 ## driver function
 function test_pdmat(C::AbstractPDMat, Cmat::Matrix;

src/utils.jl

@@ -7,7 +7,7 @@ macro check_argdims(cond)
     end
 end
 
-_rcopy!(r::StridedVecOrMat, x::StridedVecOrMat) = (r === x || copy!(r, x); r)
+_rcopy!(r::StridedVecOrMat, x::StridedVecOrMat) = (r === x || copyto!(r, x); r)
 
 
 function _addscal!(r::Matrix, a::Matrix, b::Union{Matrix, SparseMatrixCSC}, c::Real)

test/addition.jl

@@ -22,16 +22,16 @@ for T in [Float64,Float32]
   for p1 in pmats, p2 in pmats
       pr = p1 + p2
       @test size(pr) == size(p1)
-      @test full(pr) ≈ full(p1) + full(p2)
+      @test Matrix(pr) ≈ Matrix(p1) + Matrix(p2)
 
       pr = pdadd(p1, p2, convert(T,1.5))
       @test size(pr) == size(p1)
-      @test full(pr) ≈ full(p1) + full(p2) * convert(T,1.5)
+      @test Matrix(pr) ≈ Matrix(p1) + Matrix(p2) * convert(T,1.5)
   end
 
   for p1 in pmats
         pr = p1 + pm4
         @test size(pr) == size(p1)
-        @test full(pr) ≈ full(p1) + pm4
+        @test Matrix(pr) ≈ Matrix(p1) + pm4
   end
 end

test/generics.jl

@@ -17,8 +17,8 @@ pmats = Any[pm1, pm2, pm3]
 pmatsa= Any[pm1a,pm2a,pm3a]
 
 for i in 1:length(pmats)
-    @test full(3.0 * pmats[i]) == full(pmatsa[i])
-    @test full(pmats[i] * 3.0) == full(pmatsa[i])
-    @test full(3 * pmats[i])   == full(pmatsa[i])
-    @test full(pmats[i] * 3)   == full(pmatsa[i])
+    @test Matrix(3.0 * pmats[i]) == Matrix(pmatsa[i])
+    @test Matrix(pmats[i] * 3.0) == Matrix(pmatsa[i])
+    @test Matrix(3 * pmats[i])   == Matrix(pmatsa[i])
+    @test Matrix(pmats[i] * 3)   == Matrix(pmatsa[i])
 end

test/pdmtypes.jl

@@ -23,7 +23,7 @@ for T in [Float64, Float32]
     X = convert(T,2.0)
 
     call_test_pdmat(PDMat(M), M) #tests of PDMat
-    call_test_pdmat(PDiagMat(V), diagm(V)) #tests of PDiagMat
+    call_test_pdmat(PDiagMat(V), Matrix(Diagonal(V))) #tests of PDiagMat
     call_test_pdmat(ScalMat(3,x), x*Matrix{T}(I, 3, 3)) #tests of ScalMat
 #=
     call_test_pdmat(PDSparseMat(sparse(M)), M)