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.