`cholesky!` L.data pointer question

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