Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove Thread-local storage from CHOLMOD and SPQR #206

Merged
merged 4 commits into from
Aug 5, 2022
Merged

Remove Thread-local storage from CHOLMOD and SPQR #206

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
12ea3a4
move from global nthreads size vector to task_local_storage
rayegun Aug 4, 2022
9a9de5c
use Ian's workaround for missing nthreads
rayegun Aug 4, 2022
ea193df
use haskey rather than key iterator
rayegun Aug 4, 2022
ed4fd23
type return, use get!
rayegun Aug 4, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 56 additions & 59 deletions src/solvers/cholmod.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,7 @@ macro cholmod_param(kwarg, code)
value = kwarg.args[2]

common_param = # Read `common.param`
Expr(:., :(COMMONS[Threads.threadid()][]), QuoteNode(param))
Expr(:., :(task_local_storage(:cholmod_common)[]), QuoteNode(param))

return quote
default_value = $common_param
Expand All @@ -151,7 +151,18 @@ macro cholmod_param(kwarg, code)
end
end

const COMMONS = Vector{Ref{cholmod_common}}()
function newcommon(; print = 0)
common = finalizer(cholmod_l_finish, Ref(cholmod_common()))
result = cholmod_l_start(common)
@assert result == TRUE "failed to run `cholmod_l_start`!"
common[].print = 0 # no printing from CHOLMOD by default
common[].error_handler = @cfunction(error_handler, Cvoid, (Cint, Cstring, Cint, Cstring))
return common
end

function getcommon()
return get!(task_local_storage(), :cholmod_common, newcommon())::Ref{cholmod_common}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this not be

get!(newcommon, task_local_storage(), :cholmod_common)::Ref{cholmod_common}

? As it is now a new object is created for every call.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The docs say it's collection, key, default right?

Copy link
Contributor

@IanButterworth IanButterworth Aug 6, 2022 •
edited by fredrikekre
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's another method where you can give a function/type as a first arg. Note that Fredrik's suggestion gives the type but doesn't evaluate it

end

const BUILD_VERSION = VersionNumber(CHOLMOD_MAIN_VERSION, CHOLMOD_SUB_VERSION, CHOLMOD_SUBSUB_VERSION)

Expand Down Expand Up @@ -218,20 +229,6 @@ function __init__()
"""
end

### Initiate CHOLMOD
### common controls the type of factorization and keeps pointers
### to temporary memory. We need to manage a copy for each thread.
nt = Threads.nthreads()
resize!(COMMONS, nt)
errorhandler = @cfunction(error_handler, Cvoid, (Cint, Cstring, Cint, Cstring))
for i in 1:nt
COMMONS[i] = Ref(cholmod_common())
result = cholmod_l_start(COMMONS[i])
@assert result == TRUE "failed to run `cholmod_l_start`!"
COMMONS[i][].print = 0 # no printing from CHOLMOD by default
COMMONS[i][].error_handler = errorhandler
end

# Register gc tracked allocator if CHOLMOD is new enough
if current_version >= v"3.0.0"
cnfg = cglobal((:SuiteSparse_config, :libsuitesparseconfig), Ptr{Cvoid})
Expand Down Expand Up @@ -390,35 +387,35 @@ Factor(FC::FactorComponent) = FC.F

# Dense wrappers
function allocate_dense(m::Integer, n::Integer, d::Integer, ::Type{Tv}) where {Tv<:VTypes}
Dense{Tv}(cholmod_l_allocate_dense(m, n, d, xtyp(Tv), COMMONS[Threads.threadid()]))
Dense{Tv}(cholmod_l_allocate_dense(m, n, d, xtyp(Tv), getcommon()))
end

function free!(p::Ptr{cholmod_dense})
cholmod_l_free_dense(Ref(p), COMMONS[Threads.threadid()]) == TRUE
cholmod_l_free_dense(Ref(p), getcommon()) == TRUE
end

function zeros(m::Integer, n::Integer, ::Type{Tv}) where Tv<:VTypes
Dense{Tv}(cholmod_l_zeros(m, n, xtyp(Tv), COMMONS[Threads.threadid()]))
Dense{Tv}(cholmod_l_zeros(m, n, xtyp(Tv), getcommon()))
end
zeros(m::Integer, n::Integer) = zeros(m, n, Float64)

function ones(m::Integer, n::Integer, ::Type{Tv}) where Tv<:VTypes
Dense{Tv}(cholmod_l_ones(m, n, xtyp(Tv), COMMONS[Threads.threadid()]))
Dense{Tv}(cholmod_l_ones(m, n, xtyp(Tv), getcommon()))
end
ones(m::Integer, n::Integer) = ones(m, n, Float64)

function eye(m::Integer, n::Integer, ::Type{Tv}) where Tv<:VTypes
Dense{Tv}(cholmod_l_eye(m, n, xtyp(Tv), COMMONS[Threads.threadid()]))
Dense{Tv}(cholmod_l_eye(m, n, xtyp(Tv), getcommon()))
end
eye(m::Integer, n::Integer) = eye(m, n, Float64)
eye(n::Integer) = eye(n, n, Float64)

function copy(A::Dense{Tv}) where Tv<:VTypes
Dense{Tv}(cholmod_l_copy_dense(A, COMMONS[Threads.threadid()]))
Dense{Tv}(cholmod_l_copy_dense(A, getcommon()))
end

function sort!(S::Sparse{Tv}) where Tv<:VTypes
cholmod_l_sort(S, COMMONS[Threads.threadid()])
cholmod_l_sort(S, getcommon())
return S
end

Expand All @@ -431,93 +428,93 @@ function norm_dense(D::Dense{Tv}, p::Integer) where Tv<:VTypes
elseif p != 0 && p != 1
throw(ArgumentError("second argument must be either 0 (Inf norm), 1, or 2"))
end
cholmod_l_norm_dense(D, p, COMMONS[Threads.threadid()])
cholmod_l_norm_dense(D, p, getcommon())
end

function check_dense(A::Dense{Tv}) where Tv<:VTypes
cholmod_l_check_dense(pointer(A), COMMONS[Threads.threadid()]) != 0
cholmod_l_check_dense(pointer(A), getcommon()) != 0
end

# Non-Dense wrappers
function allocate_sparse(nrow::Integer, ncol::Integer, nzmax::Integer,
sorted::Bool, packed::Bool, stype::Integer, ::Type{Tv}) where {Tv<:VTypes}
Sparse{Tv}(cholmod_l_allocate_sparse(nrow, ncol, nzmax, sorted, packed, stype,
xtyp(Tv), COMMONS[Threads.threadid()]))
xtyp(Tv), getcommon()))
end

function free!(ptr::Ptr{cholmod_sparse})
cholmod_l_free_sparse(Ref(ptr), COMMONS[Threads.threadid()]) == TRUE
cholmod_l_free_sparse(Ref(ptr), getcommon()) == TRUE
end

function free!(ptr::Ptr{cholmod_factor})
# Warning! Important that finalizer doesn't modify the global Common struct.
cholmod_l_free_factor(Ref(ptr), COMMONS[Threads.threadid()]) == TRUE
cholmod_l_free_factor(Ref(ptr), getcommon()) == TRUE
end

function aat(A::Sparse{Tv}, fset::Vector{SuiteSparse_long}, mode::Integer) where Tv<:VRealTypes
Sparse{Tv}(cholmod_l_aat(A, fset, length(fset), mode, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_aat(A, fset, length(fset), mode, getcommon()))
end

function sparse_to_dense(A::Sparse{Tv}) where Tv<:VTypes
Dense{Tv}(cholmod_l_sparse_to_dense(A, COMMONS[Threads.threadid()]))
Dense{Tv}(cholmod_l_sparse_to_dense(A, getcommon()))
end
function dense_to_sparse(D::Dense{Tv}, ::Type{SuiteSparse_long}) where Tv<:VTypes
Sparse{Tv}(cholmod_l_dense_to_sparse(D, true, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_dense_to_sparse(D, true, getcommon()))
end

function factor_to_sparse!(F::Factor{Tv}) where Tv<:VTypes
ss = unsafe_load(pointer(F))
ss.xtype == CHOLMOD_PATTERN && throw(CHOLMODException("only numeric factors are supported"))
Sparse{Tv}(cholmod_l_factor_to_sparse(F, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_factor_to_sparse(F, getcommon()))
end

function change_factor!(F::Factor{Tv}, to_ll::Bool, to_super::Bool, to_packed::Bool,
to_monotonic::Bool) where Tv<:VTypes
cholmod_l_change_factor(xtyp(Tv), to_ll, to_super, to_packed, to_monotonic, F, COMMONS[Threads.threadid()]) == TRUE
cholmod_l_change_factor(xtyp(Tv), to_ll, to_super, to_packed, to_monotonic, F, getcommon()) == TRUE
end

function check_sparse(A::Sparse{Tv}) where Tv<:VTypes
cholmod_l_check_sparse(A, COMMONS[Threads.threadid()]) != 0
cholmod_l_check_sparse(A, getcommon()) != 0
end

function check_factor(F::Factor{Tv}) where Tv<:VTypes
cholmod_l_check_factor(F, COMMONS[Threads.threadid()]) != 0
cholmod_l_check_factor(F, getcommon()) != 0
end

nnz(A::Sparse{<:VTypes}) = cholmod_l_nnz(A, COMMONS[Threads.threadid()])
nnz(A::Sparse{<:VTypes}) = cholmod_l_nnz(A, getcommon())

function speye(m::Integer, n::Integer, ::Type{Tv}) where Tv<:VTypes
Sparse{Tv}(cholmod_l_speye(m, n, xtyp(Tv), COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_speye(m, n, xtyp(Tv), getcommon()))
end

function spzeros(m::Integer, n::Integer, nzmax::Integer, ::Type{Tv}) where Tv<:VTypes
Sparse{Tv}(cholmod_l_spzeros(m, n, nzmax, xtyp(Tv), COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_spzeros(m, n, nzmax, xtyp(Tv), getcommon()))
end

function transpose_(A::Sparse{Tv}, values::Integer) where Tv<:VTypes
Sparse{Tv}(cholmod_l_transpose(A, values, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_transpose(A, values, getcommon()))
end

function copy(F::Factor{Tv}) where Tv<:VTypes
Factor{Tv}(cholmod_l_copy_factor(F, COMMONS[Threads.threadid()]))
Factor{Tv}(cholmod_l_copy_factor(F, getcommon()))
end
function copy(A::Sparse{Tv}) where Tv<:VTypes
Sparse{Tv}(cholmod_l_copy_sparse(A, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_copy_sparse(A, getcommon()))
end
function copy(A::Sparse{Tv}, stype::Integer, mode::Integer) where Tv<:VRealTypes
Sparse{Tv}(cholmod_l_copy(A, stype, mode, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_copy(A, stype, mode, getcommon()))
end

function print_sparse(A::Sparse{Tv}, name::String) where Tv<:VTypes
isascii(name) || error("non-ASCII name: $name")
@cholmod_param print = 3 begin
cholmod_l_print_sparse(A, name, COMMONS[Threads.threadid()])
cholmod_l_print_sparse(A, name, getcommon())
end
nothing
end
function print_factor(F::Factor{Tv}, name::String) where Tv<:VTypes
@cholmod_param print = 3 begin
cholmod_l_print_factor(F, name, COMMONS[Threads.threadid()])
cholmod_l_print_factor(F, name, getcommon())
end
nothing
end
Expand All @@ -529,18 +526,18 @@ function ssmult(A::Sparse{Tv}, B::Sparse{Tv}, stype::Integer,
if lA.ncol != lB.nrow
throw(DimensionMismatch("inner matrix dimensions do not fit"))
end
Sparse{Tv}(cholmod_l_ssmult(A, B, stype, values, sorted, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_ssmult(A, B, stype, values, sorted, getcommon()))
end

function norm_sparse(A::Sparse{Tv}, norm::Integer) where Tv<:VTypes
if norm != 0 && norm != 1
throw(ArgumentError("norm argument must be either 0 or 1"))
end
cholmod_l_norm_sparse(A, norm, COMMONS[Threads.threadid()])
cholmod_l_norm_sparse(A, norm, getcommon())
end

function horzcat(A::Sparse{Tv}, B::Sparse{Tv}, values::Bool) where Tv<:VRealTypes
Sparse{Tv}(cholmod_l_horzcat(A, B, values, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_horzcat(A, B, values, getcommon()))
end

function scale!(S::Dense{Tv}, scale::Integer, A::Sparse{Tv}) where Tv<:VRealTypes
Expand All @@ -567,7 +564,7 @@ function scale!(S::Dense{Tv}, scale::Integer, A::Sparse{Tv}) where Tv<:VRealType
end

sA = unsafe_load(pointer(A))
cholmod_l_scale(S, scale, A, COMMONS[Threads.threadid()])
cholmod_l_scale(S, scale, A, getcommon())
A
end

Expand All @@ -579,12 +576,12 @@ function sdmult!(A::Sparse{Tv}, transpose::Bool,
if nc != size(X, 1)
throw(DimensionMismatch("incompatible dimensions, $nc and $(size(X,1))"))
end
cholmod_l_sdmult(A, transpose, Ref(α), Ref(β), X, Y, COMMONS[Threads.threadid()])
cholmod_l_sdmult(A, transpose, Ref(α), Ref(β), X, Y, getcommon())
Y
end

function vertcat(A::Sparse{Tv}, B::Sparse{Tv}, values::Bool) where Tv<:VRealTypes
Sparse{Tv}(cholmod_l_vertcat(A, B, values, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_vertcat(A, B, values, getcommon()))
end

function symmetry(A::Sparse{Tv}, option::Integer) where Tv<:VTypes
Expand All @@ -593,27 +590,27 @@ function symmetry(A::Sparse{Tv}, option::Integer) where Tv<:VTypes
nzoffdiag = Ref{SuiteSparse_long}()
nzdiag = Ref{SuiteSparse_long}()
rv = cholmod_l_symmetry(A, option, xmatched, pmatched,
nzoffdiag, nzdiag, COMMONS[Threads.threadid()])
nzoffdiag, nzdiag, getcommon())
rv, xmatched[], pmatched[], nzoffdiag[], nzdiag[]
end

# For analyze, analyze_p, and factorize_p!, the Common argument must be
# supplied in order to control if the factorization is LLt or LDLt
function analyze(A::Sparse{Tv}) where Tv<:VTypes
Factor{Tv}(cholmod_l_analyze(A, COMMONS[Threads.threadid()]))
Factor{Tv}(cholmod_l_analyze(A, getcommon()))
end
function analyze_p(A::Sparse{Tv}, perm::Vector{SuiteSparse_long}) where Tv<:VTypes
length(perm) != size(A,1) && throw(BoundsError())
Factor{Tv}(cholmod_l_analyze_p(A, perm, C_NULL, 0, COMMONS[Threads.threadid()]))
Factor{Tv}(cholmod_l_analyze_p(A, perm, C_NULL, 0, getcommon()))
end
function factorize!(A::Sparse{Tv}, F::Factor{Tv}) where Tv<:VTypes
cholmod_l_factorize(A, F, COMMONS[Threads.threadid()])
cholmod_l_factorize(A, F, getcommon())
F
end
function factorize_p!(A::Sparse{Tv}, β::Real, F::Factor{Tv}) where Tv<:VTypes
# note that β is passed as a complex number (double beta[2]),
# but the CHOLMOD manual says that only beta[0] (real part) is used
cholmod_l_factorize_p(A, Ref{Cdouble}(β), C_NULL, 0, F, COMMONS[Threads.threadid()])
cholmod_l_factorize_p(A, Ref{Cdouble}(β), C_NULL, 0, F, getcommon())
F
end

Expand All @@ -630,20 +627,20 @@ function solve(sys::Integer, F::Factor{Tv}, B::Dense{Tv}) where Tv<:VTypes
throw(LinearAlgebra.ZeroPivotException(s.minor))
end
end
Dense{Tv}(cholmod_l_solve(sys, F, B, COMMONS[Threads.threadid()]))
Dense{Tv}(cholmod_l_solve(sys, F, B, getcommon()))
end

function spsolve(sys::Integer, F::Factor{Tv}, B::Sparse{Tv}) where Tv<:VTypes
if size(F,1) != size(B,1)
throw(DimensionMismatch("LHS and RHS should have the same number of rows. " *
"LHS has $(size(F,1)) rows, but RHS has $(size(B,1)) rows."))
end
Sparse{Tv}(cholmod_l_spsolve(sys, F, B, COMMONS[Threads.threadid()]))
Sparse{Tv}(cholmod_l_spsolve(sys, F, B, getcommon()))
end

# Autodetects the types
function read_sparse(file::Libc.FILE, ::Type{SuiteSparse_long})
Sparse(cholmod_l_read_sparse(file.ptr, COMMONS[Threads.threadid()]))
Sparse(cholmod_l_read_sparse(file.ptr, getcommon()))
end

function read_sparse(file::IO, T)
Expand Down Expand Up @@ -1418,7 +1415,7 @@ function lowrankupdowndate!(F::Factor{Tv}, C::Sparse{Tv}, update::Cint) where Tv
if lF.n != lC.nrow
throw(DimensionMismatch("matrix dimensions do not fit"))
end
cholmod_l_updown(update, C, F, COMMONS[Threads.threadid()])
cholmod_l_updown(update, C, F, getcommon())
F
end

Expand Down
6 changes: 3 additions & 3 deletions src/solvers/spqr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ function _qr!(ordering::Integer, tol::Real, econ::Integer, getCTX::Integer,
H, # m-by-nh Householder vectors
HPinv, # size m row permutation
HTau, # 1-by-nh Householder coefficients
CHOLMOD.COMMONS[Threads.threadid()]) # /* workspace and parameters */
CHOLMOD.getcommon()) # /* workspace and parameters */

if rnk < 0
error("Sparse QR factorization failed")
Expand All @@ -83,7 +83,7 @@ function _qr!(ordering::Integer, tol::Real, econ::Integer, getCTX::Integer,
# Free memory allocated by SPQR. This call will make sure that the
# correct deallocator function is called and that the memory count in
# the common struct is updated
cholmod_l_free(n, sizeof(CHOLMOD.SuiteSparse_long), e, CHOLMOD.COMMONS[Threads.threadid()])
cholmod_l_free(n, sizeof(CHOLMOD.SuiteSparse_long), e, CHOLMOD.getcommon())
end
hpinv = HPinv[]
if hpinv == C_NULL
Expand All @@ -96,7 +96,7 @@ function _qr!(ordering::Integer, tol::Real, econ::Integer, getCTX::Integer,
# Free memory allocated by SPQR. This call will make sure that the
# correct deallocator function is called and that the memory count in
# the common struct is updated
cholmod_l_free(m, sizeof(CHOLMOD.SuiteSparse_long), hpinv, CHOLMOD.COMMONS[Threads.threadid()])
cholmod_l_free(m, sizeof(CHOLMOD.SuiteSparse_long), hpinv, CHOLMOD.getcommon())
end

return rnk, _E, _HPinv
Expand Down
Loading