Skip to content

Commit

Permalink
Reorganization of how memory allocated by umfpack is handled.
Browse files Browse the repository at this point in the history 
Step towards letting the user store and resuse a computed (opaque) LU
factorization.
  • Loading branch information
@mlubin
mlubin committed May 13, 2012
1 parent e4ca3a9 commit 8e047d3
Showing 1 changed file with 27 additions and 20 deletions.
47 changes: 27 additions & 20 deletions extras/linalg_suitesparse.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
type MatrixIllConditionedException <: Exception end

# Wrapper for memory allocated by umfpack. Carry along the value and index types.
type UmfpackPtr{Tv<:Union(Float64,Complex128),Ti<:Union(Int64,Int32)}
val::Vector{Ptr{Void}}
end

function _jl_cholmod_transpose{Tv<:Union(Float64,Complex128)}(S::SparseMatrixCSC{Tv})
cm = _jl_cholmod_start()
S = _jl_convert_to_0_based_indexing!(S)
Expand Down Expand Up @@ -49,9 +54,7 @@ function _jl_sparse_lusolve(S, b)
try
symbolic = _jl_umfpack_symbolic(S)
numeric = _jl_umfpack_numeric(S, symbolic)
_jl_umfpack_free_symbolic(S, symbolic)
x = _jl_umfpack_solve(S, b, numeric)
_jl_umfpack_free_numeric(S, numeric)
catch e
S = _jl_convert_to_1_based_indexing!(S)
if is(e,MatrixIllConditionedException)
Expand Down Expand Up @@ -337,29 +340,31 @@ for (f_sym_r, f_sym_c, inttype) in

function _jl_umfpack_symbolic{Tv<:Float64,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti})
# Pointer to store the symbolic factorization returned by UMFPACK
Symbolic = Array(Ptr{Void}, 1)
Symbolic = UmfpackPtr{Tv,Ti}(Array(Ptr{Void},1))
status = ccall(dlsym(_jl_libumfpack, $f_sym_r),
Ti,
(Ti, Ti,
Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
S.m, S.n,
S.colptr, S.rowval, S.nzval, Symbolic, C_NULL, C_NULL)
if status != _jl_UMFPACK_OK; error("Error in symoblic factorization"); end
S.colptr, S.rowval, S.nzval, Symbolic.val, C_NULL, C_NULL)
if status != _jl_UMFPACK_OK; error("Error in symbolic factorization"); end
finalizer(Symbolic,_jl_umfpack_free_symbolic)
return Symbolic
end

function _jl_umfpack_symbolic{Tv<:Complex128,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti})
# Pointer to store the symbolic factorization returned by UMFPACK
Symbolic = Array(Ptr{Void}, 1)
Symbolic = UmfpackPtr{Tv,Ti}(Array(Ptr{Void},1))
status = ccall(dlsym(_jl_libumfpack, $f_sym_c),
Ti,
(Ti, Ti,
Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Void},
Ptr{Float64}, Ptr{Float64}),
S.m, S.n,
S.colptr, S.rowval, real(S.nzval), imag(S.nzval), Symbolic,
S.colptr, S.rowval, real(S.nzval), imag(S.nzval), Symbolic.val,
C_NULL, C_NULL)
if status != _jl_UMFPACK_OK; error("Error in symoblic factorization"); end
if status != _jl_UMFPACK_OK; error("Error in symbolic factorization"); end
finalizer(Symbolic,_jl_umfpack_free_symbolic) # Check: do we need to free if there was an error?
return Symbolic
end

Expand All @@ -373,29 +378,31 @@ for (f_num_r, f_num_c, inttype) in

function _jl_umfpack_numeric{Tv<:Float64,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, Symbolic)
# Pointer to store the numeric factorization returned by UMFPACK
Numeric = Array(Ptr{Void}, 1)
Numeric = UmfpackPtr{Tv,Ti}(Array(Ptr{Void},1))
status = ccall(dlsym(_jl_libumfpack, $f_num_r),
Ti,
(Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Void}, Ptr{Void},
Ptr{Float64}, Ptr{Float64}),
S.colptr, S.rowval, S.nzval, Symbolic[1], Numeric,
S.colptr, S.rowval, S.nzval, Symbolic.val[1], Numeric.val,
C_NULL, C_NULL)
if status > 0; throw(MatrixIllConditionedException); end
if status != _jl_UMFPACK_OK; error("Error in numeric factorization"); end
finalizer(Numeric,_jl_umfpack_free_numeric)
return Numeric
end

function _jl_umfpack_numeric{Tv<:Complex128,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, Symbolic)
# Pointer to store the numeric factorization returned by UMFPACK
Numeric = Array(Ptr{Void}, 1)
Numeric = UmfpackPtr{Tv,Ti}(Array(Ptr{Void},1))
status = ccall(dlsym(_jl_libumfpack, $f_num_c),
Ti,
(Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Void},
Ptr{Float64}, Ptr{Float64}),
S.colptr, S.rowval, real(S.nzval), imag(S.nzval), Symbolic[1], Numeric,
S.colptr, S.rowval, real(S.nzval), imag(S.nzval), Symbolic.val[1], Numeric.val,
C_NULL, C_NULL)
if status > 0; throw(MatrixIllConditionedException); end
if status != _jl_UMFPACK_OK; error("Error in numeric factorization"); end
finalizer(Numeric,_jl_umfpack_free_numeric)
return Numeric
end

Expand All @@ -408,28 +415,28 @@ for (f_sol_r, f_sol_c, inttype) in
@eval begin

function _jl_umfpack_solve{Tv<:Float64,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti},
b::Vector{Tv}, Numeric)
b::Vector{Tv}, Numeric::UmfpackPtr{Tv,Ti})
x = similar(b)
status = ccall(dlsym(_jl_libumfpack, $f_sol_r),
Ti,
(Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64},
Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
_jl_UMFPACK_A, S.colptr, S.rowval, S.nzval,
x, b, Numeric[1], C_NULL, C_NULL)
x, b, Numeric.val[1], C_NULL, C_NULL)
if status != _jl_UMFPACK_OK; error("Error in solve"); end
return x
end

function _jl_umfpack_solve{Tv<:Complex128,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti},
b::Vector{Tv}, Numeric)
b::Vector{Tv}, Numeric::UmfpackPtr{Tv,Ti})
xr = similar(b, Float64)
xi = similar(b, Float64)
status = ccall(dlsym(_jl_libumfpack, $f_sol_c),
Ti,
(Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
_jl_UMFPACK_A, S.colptr, S.rowval, real(S.nzval), imag(S.nzval),
xr, xi, real(b), imag(b), Numeric[1], C_NULL, C_NULL)
xr, xi, real(b), imag(b), Numeric.val[1], C_NULL, C_NULL)
if status != _jl_UMFPACK_OK; error("Error in solve"); end
return complex(xr,xi)
end
Expand All @@ -444,11 +451,11 @@ for (f_symfree, f_numfree, elty, inttype) in
("umfpack_zl_free_symbolic","umfpack_zl_free_numeric",:Complex128,:Int64))
@eval begin

_jl_umfpack_free_symbolic{Tv<:$elty,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, Symbolic) =
ccall(dlsym(_jl_libumfpack, $f_symfree), Void, (Ptr{Void},), Symbolic)
_jl_umfpack_free_symbolic{Tv<:$elty,Ti<:$inttype}(Symbolic::UmfpackPtr{Tv,Ti}) =
ccall(dlsym(_jl_libumfpack, $f_symfree), Void, (Ptr{Void},), Symbolic.val)

_jl_umfpack_free_numeric{Tv<:$elty,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, Numeric) =
ccall(dlsym(_jl_libumfpack, $f_numfree), Void, (Ptr{Void},), Numeric)
_jl_umfpack_free_numeric{Tv<:$elty,Ti<:$inttype}(Numeric::UmfpackPtr{Tv,Ti}) =
ccall(dlsym(_jl_libumfpack, $f_numfree), Void, (Ptr{Void},), Numeric.val)

end
end

0 comments on commit 8e047d3

Please sign in to comment.