Merge pull request #134 from sethaxen/cholnocopy

Fix #133

Author: andreasnoack

Project.toml

@@ -1,6 +1,6 @@
 name = "PDMats"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.11.0"
+version = "0.11.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/chol.jl

@@ -1,6 +1,10 @@
 CholType{T,S<:AbstractMatrix} = Cholesky{T,S}
 chol_lower(a::Matrix) = cholesky(a).L
 
+# always use the stored cholesky factor, not a copy
+chol_lower(a::CholType) = a.uplo === 'L' ? a.L : a.U'
+chol_upper(a::CholType) = a.uplo === 'U' ? a.U : a.L'
+
 if HAVE_CHOLMOD
     CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T}
 

src/pdmat.jl

@@ -61,22 +61,20 @@ Base.kron(A::PDMat, B::PDMat) = PDMat(kron(A.mat, B.mat), Cholesky(kron(A.chol.U
 ### whiten and unwhiten
 
 function whiten!(r::StridedVecOrMat, a::PDMat, x::StridedVecOrMat)
-    cf = a.chol.U
     v = _rcopy!(r, x)
-    ldiv!(transpose(cf), v)
+    ldiv!(chol_lower(a.chol), v)
 end
 
 function unwhiten!(r::StridedVecOrMat, a::PDMat, x::StridedVecOrMat)
-    cf = a.chol.U
     v = _rcopy!(r, x)
-    lmul!(transpose(cf), v)
+    lmul!(chol_lower(a.chol), v)
 end
 
 
 ### quadratic forms
 
-quad(a::PDMat, x::AbstractVector) = sum(abs2, a.chol.U * x)
-invquad(a::PDMat, x::AbstractVector) = sum(abs2, a.chol.L \ x)
+quad(a::PDMat, x::AbstractVector) = sum(abs2, chol_upper(a.chol) * x)
+invquad(a::PDMat, x::AbstractVector) = sum(abs2, chol_lower(a.chol) \ x)
 
 """
     quad!(r::AbstractArray, a::AbstractPDMat, x::StridedMatrix)
@@ -96,25 +94,21 @@ invquad!(r::AbstractArray, a::PDMat, x::StridedMatrix) = colwise_dot!(r, x, a.ma
 ### tri products
 
 function X_A_Xt(a::PDMat, x::StridedMatrix)
-    cf = a.chol.U
-    z = rmul!(copy(x), transpose(cf))
+    z = rmul!(copy(x), chol_lower(a.chol))
     return z * transpose(z)
 end
 
 function Xt_A_X(a::PDMat, x::StridedMatrix)
-    cf = a.chol.U
-    z = lmul!(cf, copy(x))
+    z = lmul!(chol_upper(a.chol), copy(x))
     return transpose(z) * z
 end
 
 function X_invA_Xt(a::PDMat, x::StridedMatrix)
-    cf = a.chol.U
-    z = rdiv!(copy(x), cf)
+    z = rdiv!(copy(x), chol_upper(a.chol))
     return z * transpose(z)
 end
 
 function Xt_invA_X(a::PDMat, x::StridedMatrix)
-    cf = a.chol.U
-    z = ldiv!(transpose(cf), copy(x))
+    z = ldiv!(chol_lower(a.chol), copy(x))
     return transpose(z) * z
 end

test/pdmtypes.jl

@@ -19,6 +19,8 @@ for T in [Float64, Float32]
     X = convert(T,2.0)
 
     test_pdmat(PDMat(M), M,                        cmat_eq=true, verbose=1) #tests of PDMat
+    cholL = Cholesky(Matrix(transpose(cholesky(M).factors)), 'L', 0)
+    test_pdmat(PDMat(cholL), M,                    cmat_eq=true, verbose=1) #tests of PDMat
     test_pdmat(PDiagMat(V), Matrix(Diagonal(V)),   cmat_eq=true, verbose=1) #tests of PDiagMat
     test_pdmat(ScalMat(3,x), x*Matrix{T}(I, 3, 3), cmat_eq=true, verbose=1) #tests of ScalMat
     test_pdmat(PDSparseMat(sparse(M)), M,          cmat_eq=true, verbose=1, t_eig=false)