Skip to content

Commit

Permalink
Optimisation of getproperty for UmfpackLU object (avoid copying irrel…
Browse files Browse the repository at this point in the history 
…evant information) (JuliaLang#40663)
  • Loading branch information
@ElOceanografo
clement-f authored and ElOceanografo committed May 4, 2021
1 parent 9b8b1ad commit 6e1f3e8
Show file tree
Hide file tree
Showing 2 changed files with 231 additions and 73 deletions.
299 changes: 226 additions & 73 deletions stdlib/SuiteSparse/src/umfpack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -451,61 +451,237 @@ for itype in UmfpackIndexTypes
lnz, unz, n_row, n_col, nz_diag, lu.numeric)
(lnz[], unz[], n_row[], n_col[], nz_diag[])
end
function umf_extract(lu::UmfpackLU{Float64,$itype})
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Lp = Vector{$itype}(undef, n_row + 1)
Lj = Vector{$itype}(undef, lnz) # L is returned in CSR (compressed sparse row) format
Lx = Vector{Float64}(undef, lnz)
Up = Vector{$itype}(undef, n_col + 1)
Ui = Vector{$itype}(undef, unz)
Ux = Vector{Float64}(undef, unz)
P = Vector{$itype}(undef, n_row)
Q = Vector{$itype}(undef, n_col)
Rs = Vector{Float64}(undef, n_row)
@isok ccall(($get_num_r,:libumfpack), $itype,
(Ptr{$itype},Ptr{$itype},Ptr{Float64},
Ptr{$itype},Ptr{$itype},Ptr{Float64},
Ptr{$itype},Ptr{$itype},Ptr{Cvoid},
Ref{$itype},Ptr{Float64},Ptr{Cvoid}),
Lp,Lj,Lx,
Up,Ui,Ux,
P, Q, C_NULL,
0, Rs, lu.numeric)
(copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row, increment!(Lp), increment!(Lj), Lx))),
SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up), increment!(Ui), Ux),
increment!(P), increment!(Q), Rs)
function getproperty(lu::UmfpackLU{Float64, $itype}, d::Symbol)
if d === :L
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Lp = Vector{$itype}(undef, n_row + 1)
# L is returned in CSR (compressed sparse row) format
Lj = Vector{$itype}(undef, lnz)
Lx = Vector{Float64}(undef, lnz)
@isok ccall(($get_num_r, :libumfpack), $itype,
(Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
Lp, Lj, Lx,
C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, lu.numeric)
return copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row,
increment!(Lp), increment!(Lj), Lx)))
elseif d === :U
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Up = Vector{$itype}(undef, n_col + 1)
Ui = Vector{$itype}(undef, unz)
Ux = Vector{Float64}(undef, unz)
@isok ccall(($get_num_r, :libumfpack), $itype,
(Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
C_NULL, C_NULL, C_NULL,
Up, Ui, Ux,
C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, lu.numeric)
return SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up),
increment!(Ui), Ux)
elseif d === :p
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
P = Vector{$itype}(undef, n_row)
@isok ccall(($get_num_r, :libumfpack), $itype,
(Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{$itype}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL,
P, C_NULL, C_NULL,
C_NULL, C_NULL, lu.numeric)
return increment!(P)
elseif d === :q
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Q = Vector{$itype}(undef, n_col)
@isok ccall(($get_num_r, :libumfpack), $itype,
(Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{$itype}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL,
C_NULL, Q, C_NULL,
C_NULL, C_NULL, lu.numeric)
return increment!(Q)
elseif d === :Rs
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Rs = Vector{Float64}(undef, n_row)
@isok ccall(($get_num_r, :libumfpack), $itype,
(Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL,
C_NULL, Rs, lu.numeric)
return Rs
elseif d === :(:)
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Lp = Vector{$itype}(undef, n_row + 1)
# L is returned in CSR (compressed sparse row) format
Lj = Vector{$itype}(undef, lnz)
Lx = Vector{Float64}(undef, lnz)
Up = Vector{$itype}(undef, n_col + 1)
Ui = Vector{$itype}(undef, unz)
Ux = Vector{Float64}(undef, unz)
P = Vector{$itype}(undef, n_row)
Q = Vector{$itype}(undef, n_col)
Rs = Vector{Float64}(undef, n_row)
@isok ccall(($get_num_r, :libumfpack), $itype,
(Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
Ptr{$itype}, Ptr{$itype}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
Lp, Lj, Lx,
Up, Ui, Ux,
P, Q, C_NULL,
C_NULL, Rs, lu.numeric)
return (copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row,
increment!(Lp), increment!(Lj),
Lx))),
SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up),
increment!(Ui), Ux),
increment!(P), increment!(Q), Rs)
else
return getfield(lu, d)
end
end
function umf_extract(lu::UmfpackLU{ComplexF64,$itype})
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Lp = Vector{$itype}(undef, n_row + 1)
Lj = Vector{$itype}(undef, lnz) # L is returned in CSR (compressed sparse row) format
Lx = Vector{Float64}(undef, lnz)
Lz = Vector{Float64}(undef, lnz)
Up = Vector{$itype}(undef, n_col + 1)
Ui = Vector{$itype}(undef, unz)
Ux = Vector{Float64}(undef, unz)
Uz = Vector{Float64}(undef, unz)
P = Vector{$itype}(undef, n_row)
Q = Vector{$itype}(undef, n_col)
Rs = Vector{Float64}(undef, n_row)
@isok ccall(($get_num_z,:libumfpack), $itype,
(Ptr{$itype},Ptr{$itype},Ptr{Float64},Ptr{Float64},
Ptr{$itype},Ptr{$itype},Ptr{Float64},Ptr{Float64},
Ptr{$itype},Ptr{$itype},Ptr{Cvoid}, Ptr{Cvoid},
Ref{$itype},Ptr{Float64},Ptr{Cvoid}),
Lp,Lj,Lx,Lz,
Up,Ui,Ux,Uz,
P, Q, C_NULL, C_NULL,
0, Rs, lu.numeric)
(copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row, increment!(Lp), increment!(Lj), complex.(Lx, Lz)))),
SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up), increment!(Ui), complex.(Ux, Uz)),
increment!(P), increment!(Q), Rs)
function getproperty(lu::UmfpackLU{ComplexF64, $itype}, d::Symbol)
if d === :L
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Lp = Vector{$itype}(undef, n_row + 1)
# L is returned in CSR (compressed sparse row) format
Lj = Vector{$itype}(undef, lnz)
Lx = Vector{Float64}(undef, lnz)
Lz = Vector{Float64}(undef, lnz)
@isok ccall(($get_num_z, :libumfpack), $itype,
(Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
Lp, Lj, Lx, Lz,
C_NULL, C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, lu.numeric)
return copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row,
increment!(Lp), increment!(Lj),
complex.(Lx, Lz))))
elseif d === :U
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Up = Vector{$itype}(undef, n_col + 1)
Ui = Vector{$itype}(undef, unz)
Ux = Vector{Float64}(undef, unz)
Uz = Vector{Float64}(undef, unz)
@isok ccall(($get_num_z, :libumfpack), $itype,
(Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
C_NULL, C_NULL, C_NULL, C_NULL,
Up, Ui, Ux, Uz,
C_NULL, C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, lu.numeric)
return SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up),
increment!(Ui), complex.(Ux, Uz))
elseif d === :p
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
P = Vector{$itype}(undef, n_row)
@isok ccall(($get_num_z, :libumfpack), $itype,
(Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{$itype}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
C_NULL, C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL, C_NULL,
P, C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, lu.numeric)
return increment!(P)
elseif d === :q
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Q = Vector{$itype}(undef, n_col)
@isok ccall(($get_num_z, :libumfpack), $itype,
(Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{$itype}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
C_NULL, C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL, C_NULL,
C_NULL, Q, C_NULL, C_NULL,
C_NULL, C_NULL, lu.numeric)
return increment!(Q)
elseif d === :Rs
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Rs = Vector{Float64}(undef, n_row)
@isok ccall(($get_num_z, :libumfpack), $itype,
(Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
C_NULL, C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL, C_NULL,
C_NULL, C_NULL, C_NULL, C_NULL,
C_NULL, Rs, lu.numeric)
return Rs
elseif d === :(:)
umfpack_numeric!(lu) # ensure the numeric decomposition exists
(lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
Lp = Vector{$itype}(undef, n_row + 1)
# L is returned in CSR (compressed sparse row) format
Lj = Vector{$itype}(undef, lnz)
Lx = Vector{Float64}(undef, lnz)
Lz = Vector{Float64}(undef, lnz)
Up = Vector{$itype}(undef, n_col + 1)
Ui = Vector{$itype}(undef, unz)
Ux = Vector{Float64}(undef, unz)
Uz = Vector{Float64}(undef, unz)
P = Vector{$itype}(undef, n_row)
Q = Vector{$itype}(undef, n_col)
Rs = Vector{Float64}(undef, n_row)
@isok ccall(($get_num_z, :libumfpack), $itype,
(Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64},
Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64},
Ptr{$itype}, Ptr{$itype}, Ptr{Cvoid}, Ptr{Cvoid},
Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
Lp, Lj, Lx, Lz,
Up, Ui, Ux, Uz,
P, Q, C_NULL, C_NULL,
C_NULL, Rs, lu.numeric)
return (copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row,
increment!(Lp), increment!(Lj),
complex.(Lx, Lz)))),
SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up),
increment!(Ui), complex.(Ux, Uz)),
increment!(P), increment!(Q), Rs)
else
return getfield(lu, d)
end
end
end
end

# backward compatibility
umfpack_extract(lu::UmfpackLU) = getproperty(lu, :(:))

function nnz(lu::UmfpackLU)
lnz, unz, = umf_lunz(lu)
return Int(lnz + unz)
Expand Down Expand Up @@ -578,29 +754,6 @@ function _AqldivB_kernel!(X::StridedMatrix{Tb}, lu::UmfpackLU{Float64},
end
end


@inline function getproperty(lu::UmfpackLU, d::Symbol)
if d === :L || d === :U || d === :p || d === :q || d === :Rs || d === :(:)
# Guard the call to umf_extract behaind a branch to avoid infinite recursion
L, U, p, q, Rs = umf_extract(lu)
if d === :L
return L
elseif d === :U
return U
elseif d === :p
return p
elseif d === :q
return q
elseif d === :Rs
return Rs
elseif d === :(:)
return (L, U, p, q, Rs)
end
else
getfield(lu, d)
end
end

for Tv in (:Float64, :ComplexF64), Ti in UmfpackIndexTypes
f = Symbol(umf_nm("free_symbolic", Tv, Ti))
@eval begin
Expand Down
5 changes: 5 additions & 0 deletions stdlib/SuiteSparse/test/umfpack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,11 @@ using SparseArrays: nnz, sparse, sprand, sprandn, SparseMatrixCSC
@test nnz(lua) == 18
@test_throws ErrorException lua.Z
L,U,p,q,Rs = lua.:(:)
@test L == lua.L
@test U == lua.U
@test p == lua.p
@test q == lua.q
@test Rs == lua.Rs
@test (Diagonal(Rs) * A)[p,q] ≈ L * U

det(lua) ≈ det(Array(A))
Expand Down

0 comments on commit 6e1f3e8

Please sign in to comment.