diff --git a/src/chol.jl b/src/chol.jl index c95ea4c..4da0669 100644 --- a/src/chol.jl +++ b/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} diff --git a/test/chol.jl b/test/chol.jl new file mode 100644 index 0000000..71d4a5e --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 54e50b1..3ac27b1 100644 --- a/test/runtests.jl +++ b/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