According to the Julia manual, the :L factor obtained from sparse Cholesky
is not useful because the permutation is not applied. The manual says that
the :PtL factor should be used instead. However, there does not seem to be
a documented method to bring the :PtL factor back into Julia:
Here is the trace from Julia 0.4.6:
julia> A = speye(7);
julia> F = cholfact(A);
julia> L = sparse(F[:L]);
julia> PtL = sparse(F[:PtL])
ERROR: MethodError: `size` has no method matching
size(::Base.SparseMatrix.CHOLMOD.Factor{Float64})
Closest candidates are:
size(::Base.SparseMatrix.CHOLMOD.Factor{Tv}, ::Integer)
size(::Any, ::Integer, ::Integer, ::Integer...)
in convert at sparse/sparsematrix.jl:214
in sparse at sparse/sparsematrix.jl:320
I wrote a function to accomplish this task, but it uses
undocumented/unexported functions and is also probably wasteful in its
array creation. (See below.) Does anyone know of a better way?
-- Steve Vavasis
# Compute the permuted cholesky factor of a sparse, symmetric
# positive definite matrix A.
function cholfactPtL(A)
n = size(A,1)
F = cholfact(A)
L0 = sparse(F[:L])
is,js,es = findnz(L0)
p = Base.SparseMatrix.CHOLMOD.get_perm(F)
sparse(p[is], p[js], es, n, n)
end
