Merge pull request #144 from st--/st/avoid_copy

Avoid copy in `chol_lower`

Author: andreasnoack

src/chol.jl

@@ -1,9 +1,11 @@
-chol_lower(a::Matrix) = cholesky(a).L
-
-# always use the stored cholesky factor, not a copy
+# Accessing a.L directly might involve an extra copy();
+# instead, always use the stored Cholesky factor:
 chol_lower(a::Cholesky) = a.uplo === 'L' ? a.L : a.U'
 chol_upper(a::Cholesky) = a.uplo === 'U' ? a.U : a.L'
 
+# For a dense Matrix, the following allows us to avoid the Adjoint wrapper:
+chol_lower(a::Matrix) = cholesky(Hermitian(a, :L)).L
+
 if HAVE_CHOLMOD
     CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T}
 

test/chol.jl

@@ -0,0 +1,20 @@
+using LinearAlgebra, PDMats
+using PDMats: chol_lower, chol_upper
+
+@testset "chol_lower" begin
+    A = rand(100, 100)
+    C = A'A
+    size_of_one_copy = sizeof(C)
+    @assert size_of_one_copy > 100  # ensure the matrix is large enough that few-byte allocations don't matter
+
+    chol_lower(C)
+    @test (@allocated chol_lower(C)) < 1.05 * size_of_one_copy  # allow 5% overhead
+
+    for uplo in (:L, :U)
+        ch = cholesky(Symmetric(C, uplo))
+        chol_lower(ch)
+        @test (@allocated chol_lower(ch)) < 50  # allow small overhead for wrapper types
+        chol_upper(ch)
+        @test (@allocated chol_upper(ch)) < 50  # allow small overhead for wrapper types
+    end
+end

test/runtests.jl

@@ -1,5 +1,5 @@
 include("testutils.jl")
-tests = ["pdmtypes", "addition", "generics", "kron"]
+tests = ["pdmtypes", "addition", "generics", "kron", "chol"]
 println("Running tests ...")
 
 for t in tests