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)
  • Loading branch information
Clément Fauchereau committed Apr 30, 2021
1 parent 3389940 commit 0fbbb20
Show file tree
Hide file tree
Showing 2 changed files with 204 additions and 61 deletions.
260 changes: 199 additions & 61 deletions stdlib/SuiteSparse/src/umfpack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -428,57 +428,208 @@ 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})
function umf_extract(lu::UmfpackLU{Float64,$itype}, d::Symbol)
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)

if d === :L
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
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
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
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
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 === :(:)
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)
end
end
function umf_extract(lu::UmfpackLU{ComplexF64,$itype})
function umf_extract(lu::UmfpackLU{ComplexF64,$itype}, d::Symbol)
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)

if d === :L
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
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
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
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
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 === :(:)
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)
end
end
end
end
Expand Down Expand Up @@ -558,21 +709,8 @@ 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
# Guard the call to umf_extract behind a branch to avoid infinite recursion
return umf_extract(lu, d)
else
getfield(lu, d)
end
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 0fbbb20

Please sign in to comment.