`cholesky!` L.data pointer question

Internals & Design
1

In the following code, I was surprised that the pointer for L.L.data isn’t the same as the pointer to B, L.UL.data, L.U.data.

using LinearAlgebra, Random, BenchmarkTools

T = Float64
D = 3
rng = Random.Xoshiro(0)

function myfunc!(B, A)
    copy!(B, A)
    L = cholesky!(B)
    return L
end

A = randn(T, D, D)
A = A'*A
B = copy(A)

@btime myfunc!($B, $A)

L = myfunc!(B, A)

pointer_vec = [pointer(L.U.data); pointer(L.L.data); pointer(L.UL.data); pointer(B)]
display(pointer_vec)

My REPL showed:

  162.419 ns (0 allocations: 0 bytes)
4-element Vector{Ptr{Float64}}:
 Ptr{Float64} @0x00007f93a66647a0
 Ptr{Float64} @0x00007f9428ff5b20
 Ptr{Float64} @0x00007f93a66647a0
 Ptr{Float64} @0x00007f93a66647a0

The timing is from @btime, and indicates no allocation. The L.L.data seems to be referencing the transpose of the same matrix as B:

julia> norm(L.L.data - B')
0.0

It is true that a pointer to B transposed can be a different value than B? If so, why does it seem that myfun! avoided (heap) allocation for the new address of L.L.data?

Ultimately, I’m trying to figure out if cholesky! allocates if I pre-allocated all its arguments. AllocCheck.jl seems to pass just fine:

using AllocCheck

@check_allocs testfunc!(a,b) = myfunc!(a,b)
testfunc!(B, A) # no errors when I run this line in the REPL
2

On a Matrix, cholesky! (and cholesky after a copy) checks ishermitian, wraps the input in Hermitian, then passes the wrapped array to a specialized cholesky! method. By default, Hermitian uses the upper triangle. However, you can manually specify which triangle is used in the Hermitian wrapper.

function myfunc!(B, A, UL)
    copy!(B, A)
    L = cholesky!(Hermitian(B, UL)) # specify which half to reference
    return L
end

julia> L = myfunc!(B, A, :U) # stores U
Cholesky{Float64, Matrix{Float64}}
U factor:
3×3 UpperTriangular{Float64, Matrix{Float64}}:
 1.15458  0.836121  -0.935387
  ⋅       2.7127     0.708932
  ⋅        ⋅         0.364093

julia> pointer_vec = [pointer(L.U.data); pointer(L.L.data); pointer(L.UL.data); pointer(B)]
4-element Vector{Ptr{Float64}}:
 Ptr{Float64} @0x00007f486f797920
 Ptr{Float64} @0x00007f486ed4dd00 # .L makes a copy
 Ptr{Float64} @0x00007f486f797920
 Ptr{Float64} @0x00007f486f797920

julia> L = myfunc!(B, A, :L) # stores L
Cholesky{Float64, Matrix{Float64}}
L factor:
3×3 LowerTriangular{Float64, Matrix{Float64}}:
  1.15458    ⋅         ⋅
  0.836121  2.7127     ⋅
 -0.935387  0.708932  0.364093

julia> pointer_vec = [pointer(L.U.data); pointer(L.L.data); pointer(L.UL.data); pointer(B)]
4-element Vector{Ptr{Float64}}:
 Ptr{Float64} @0x00007f486ed4ff20 # .U makes a copy
 Ptr{Float64} @0x00007f486f797920
 Ptr{Float64} @0x00007f486f797920
 Ptr{Float64} @0x00007f486f797920

When the non-stored triangle is requested, .U/.L will make a copy of the proper triangularity. The implementation can be found via @edit L.U and (in v1.11) is

function getproperty(C::Cholesky, d::Symbol)
    Cfactors = getfield(C, :factors)
    Cuplo    = getfield(C, :uplo)
    if d === :U
        return UpperTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors'))
    elseif d === :L
        return LowerTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors'))
    elseif d === :UL
        return (Cuplo === 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors))
    else
        return getfield(C, d)
    end
end

If you want the non-stored triangle without the copy, you can use (e.g., if U is stored), L.U' instead of L.L. I don’t know why the implementation prefers to make a copy rather than do this, but I assume there’s a decent reason.

1 Like