remove copy-allocation on accessing cholesky factors (.L, .U) (#1186)

Author: st--

src/cholesky.jl

@@ -65,7 +65,7 @@ julia> C.U
   ⋅    ⋅    3.0
 
 julia> C.L
-3×3 LowerTriangular{Float64, Matrix{Float64}}:
+3×3 LowerTriangular{Float64, Adjoint{Matrix{Float64}}}:
   2.0   ⋅    ⋅
   6.0  1.0   ⋅
  -8.0  5.0  3.0
@@ -530,7 +530,7 @@ julia> C.U
   ⋅    ⋅    3.0
 
 julia> C.L
-3×3 LowerTriangular{Float64, Matrix{Float64}}:
+3×3 LowerTriangular{Float64, Adjoint{Matrix{Float64}}}:
   2.0   ⋅    ⋅
   6.0  1.0   ⋅
  -8.0  5.0  3.0
@@ -668,14 +668,14 @@ function _choleskyUfactor(Cfactors, Cuplo)
     if Cuplo === 'U'
         return UpperTriangular(Cfactors)
     else
-        return copy(LowerTriangular(Cfactors)')
+        return LowerTriangular(Cfactors)'
     end
 end
 function _choleskyLfactor(Cfactors, Cuplo)
     if Cuplo === 'L'
         return LowerTriangular(Cfactors)
     else
-        return copy(UpperTriangular(Cfactors)')
+        return UpperTriangular(Cfactors)'
     end
 end
 

test/cholesky.jl

@@ -595,18 +595,18 @@ end
     B = cholesky(A)
     B32 = cholesky(Float32.(A))
     @test B isa Cholesky{Float16, Matrix{Float16}}
-    @test B.U isa UpperTriangular{Float16, Matrix{Float16}}
-    @test B.L isa LowerTriangular{Float16, Matrix{Float16}}
-    @test B.UL isa UpperTriangular{Float16, Matrix{Float16}}
+    @test B.U isa UpperTriangular{Float16, <:AbstractMatrix{Float16}}
+    @test B.L isa LowerTriangular{Float16, <:AbstractMatrix{Float16}}
+    @test B.UL isa UpperTriangular{Float16, <:AbstractMatrix{Float16}}
     @test B.U ≈ B32.U
     @test B.L ≈ B32.L
     @test B.UL ≈ B32.UL
     @test Matrix(B) ≈ A
     B = cholesky(A, RowMaximum())
     B32 = cholesky(Float32.(A), RowMaximum())
-    @test B isa CholeskyPivoted{Float16,Matrix{Float16}}
-    @test B.U isa UpperTriangular{Float16, Matrix{Float16}}
-    @test B.L isa LowerTriangular{Float16, Matrix{Float16}}
+    @test B isa CholeskyPivoted{Float16, <:AbstractMatrix{Float16}}
+    @test B.U isa UpperTriangular{Float16, <:AbstractMatrix{Float16}}
+    @test B.L isa LowerTriangular{Float16, <:AbstractMatrix{Float16}}
     @test B.U ≈ B32.U
     @test B.L ≈ B32.L
     @test Matrix(B) ≈ A
@@ -658,4 +658,22 @@ end
     end
 end
 
+@testset "accessing both L and U factors should avoid allocations" begin
+    n = 30
+    A = rand(n, n)
+    Apd = A'A
+    allowed_cost_of_overhead = 32
+    @assert sizeof(Apd) > 4allowed_cost_of_overhead  # ensure that we could positively identify extra copies
+
+    for uplo in (:L, :U)
+        C = Symmetric(Apd, uplo)
+        for val in (NoPivot(), RowMaximum())
+            B = cholesky(C, val)
+            B.L, B.U  # access once to ensure the accessor is compiled already
+            @test (@allocated B.L) <= allowed_cost_of_overhead
+            @test (@allocated B.U) <= allowed_cost_of_overhead
+        end
+    end
+end
+
 end # module TestCholesky