Skip to content

Commit 4865eb0

Browse files
ViralBShahViralBShah
committed
Sparse matrices can use Int32 or Int64 indices. The default is to use Int32,
since memory usage doubles with Int64 indices. We need a good way to create sparse matrices with Int64 indices. The index type is essentially decided in constructors, and then percolates through to other functions. UMFPACK can now solve matrices with both, Int32 and Int64 indices.
1 parent d65bfce commit 4865eb0

File tree

2 files changed

+137
-160
lines changed

2 files changed

+137
-160
lines changed

‎j/linalg_suitesparse.j

+56-48
Original file line numberDiff line numberDiff line change
@@ -1,21 +1,21 @@
11
_jl_libSuiteSparse = dlopen("libSuiteSparse")
22

33
## Type of solve
4-
const _jl_UMFPACK_A = int64(0) # Ax=b
5-
const _jl_UMFPACK_At = int64(1) # A'x=b
6-
const _jl_UMFPACK_Aat = int64(2) # A.'x=b
7-
const _jl_UMFPACK_Pt_L = int64(3) # P'Lx=b
8-
const _jl_UMFPACK_L = int64(4) # Lx=b
9-
const _jl_UMFPACK_Lt_P = int64(5) # L'Px=b
10-
const _jl_UMFPACK_Lat_P = int64(6) # L.'Px=b
11-
const _jl_UMFPACK_Lt = int64(7) # L'x=b
12-
const _jl_UMFPACK_Lat = int64(8) # L.'x=b
13-
const _jl_UMFPACK_U_Qt = int64(9) # UQ'x=b
14-
const _jl_UMFPACK_U = int64(10) # Ux=b
15-
const _jl_UMFPACK_Q_Ut = int64(11) # QU'x=b
16-
const _jl_UMFPACK_Q_Uat = int64(12) # QU.'x=b
17-
const _jl_UMFPACK_Ut = int64(13) # U'x=b
18-
const _jl_UMFPACK_Uat = int64(14) # U.'x=b
4+
const _jl_UMFPACK_A = 0 # Ax=b
5+
const _jl_UMFPACK_At = 1 # A'x=b
6+
const _jl_UMFPACK_Aat = 2 # A.'x=b
7+
const _jl_UMFPACK_Pt_L = 3 # P'Lx=b
8+
const _jl_UMFPACK_L = 4 # Lx=b
9+
const _jl_UMFPACK_Lt_P = 5 # L'Px=b
10+
const _jl_UMFPACK_Lat_P = 6 # L.'Px=b
11+
const _jl_UMFPACK_Lt = 7 # L'x=b
12+
const _jl_UMFPACK_Lat = 8 # L.'x=b
13+
const _jl_UMFPACK_U_Qt = 9 # UQ'x=b
14+
const _jl_UMFPACK_U = 10 # Ux=b
15+
const _jl_UMFPACK_Q_Ut = 11 # QU'x=b
16+
const _jl_UMFPACK_Q_Uat = 12 # QU.'x=b
17+
const _jl_UMFPACK_Ut = 13 # U'x=b
18+
const _jl_UMFPACK_Uat = 14 # U.'x=b
1919

2020
## Sizes of Control and Info arrays for returning information from solver
2121
const _jl_UMFPACK_INFO = 90
@@ -53,50 +53,58 @@ function _jl_convert_to_1_based_indexing(S::SparseMatrixCSC)
5353
end
5454

5555
## Wrappers around UMFPACK routines
56-
function _jl_umfpack_symbolic{Tv<:Float64,Ti<:Int64}(S::SparseMatrixCSC{Tv,Ti})
57-
# Pointer to store the symbolic factorization returned by UMFPACK
58-
Symbolic = Array(Ptr{Void}, 1)
59-
status = ccall(dlsym(_jl_libSuiteSparse, :umfpack_dl_symbolic),
60-
Ti,
61-
(Ti, Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Tv}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
62-
S.m, S.n, S.colptr, S.rowval, S.nzval, Symbolic, C_NULL, C_NULL)
63-
return Symbolic
64-
end
6556

66-
function _jl_umfpack_numeric{Tv<:Float64,Ti<:Int64}(S::SparseMatrixCSC{Tv,Ti}, Symbolic)
67-
# Pointer to store the numeric factorization returned by UMFPACK
68-
Numeric = Array(Ptr{Void}, 1)
69-
status = ccall(dlsym(_jl_libSuiteSparse, :umfpack_dl_numeric),
70-
Ti,
71-
(Ptr{Ti}, Ptr{Ti}, Ptr{Tv}, Ptr{Void}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
72-
S.colptr, S.rowval, S.nzval, Symbolic, Numeric, C_NULL, C_NULL)
73-
return Numeric
74-
end
57+
macro _jl_umfpack_macro(eltype, inttype, f_sym, f_num, f_sol, f_symfree, f_numfree)
58+
quote
59+
function _jl_umfpack_symbolic{Tv<:$eltype,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti})
60+
# Pointer to store the symbolic factorization returned by UMFPACK
61+
Symbolic = Array(Ptr{Void}, 1)
62+
status = ccall(dlsym(_jl_libSuiteSparse, $f_sym),
63+
Ti,
64+
(Ti, Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Tv}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
65+
S.m, S.n, S.colptr, S.rowval, S.nzval, Symbolic, C_NULL, C_NULL)
66+
return Symbolic
67+
end
7568

76-
function _jl_umfpack_solve{Tv<:Float64,Ti<:Int64}(S::SparseMatrixCSC{Tv,Ti}, b::Vector{Tv}, Numeric)
77-
x = similar(b)
78-
status = ccall(dlsym(_jl_libSuiteSparse, :umfpack_dl_solve),
79-
Ti,
80-
(Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Tv}, Ptr{Tv}, Ptr{Tv}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
81-
_jl_UMFPACK_A, S.colptr, S.rowval, S.nzval, x, b, Numeric, C_NULL, C_NULL)
82-
return x
83-
end
69+
function _jl_umfpack_numeric{Tv<:$eltype,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, Symbolic)
70+
# Pointer to store the numeric factorization returned by UMFPACK
71+
Numeric = Array(Ptr{Void}, 1)
72+
status = ccall(dlsym(_jl_libSuiteSparse, $f_num),
73+
Ti,
74+
(Ptr{Ti}, Ptr{Ti}, Ptr{Tv}, Ptr{Void}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
75+
S.colptr, S.rowval, S.nzval, Symbolic[1], Numeric, C_NULL, C_NULL)
76+
return Numeric
77+
end
8478

85-
_jl_umfpack_free_symbolic{Tv<:Float64,Ti<:Int64}(S::SparseMatrixCSC{Tv,Ti}, Symbolic) =
86-
ccall(dlsym(_jl_libSuiteSparse, :umfpack_dl_free_symbolic), Void, (Ptr{Void},), Symbolic)
79+
function _jl_umfpack_solve{Tv<:$eltype,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, b::Vector{Tv}, Numeric)
80+
x = similar(b)
81+
status = ccall(dlsym(_jl_libSuiteSparse, $f_sol),
82+
Ti,
83+
(Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Tv}, Ptr{Tv}, Ptr{Tv}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
84+
convert(Ti, _jl_UMFPACK_A), S.colptr, S.rowval, S.nzval, x, b, Numeric[1], C_NULL, C_NULL)
85+
return x
86+
end
8787

88-
_jl_umfpack_free_numeric{Tv<:Float64,Ti<:Int64}(S::SparseMatrixCSC{Tv,Ti}, Numeric) =
89-
ccall(dlsym(_jl_libSuiteSparse, :umfpack_dl_free_numeric), Void, (Ptr{Void},), Numeric)
88+
_jl_umfpack_free_symbolic{Tv<:$eltype,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, Symbolic) =
89+
ccall(dlsym(_jl_libSuiteSparse, $f_symfree), Void, (Ptr{Void},), Symbolic)
90+
91+
_jl_umfpack_free_numeric{Tv<:$eltype,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, Numeric) =
92+
ccall(dlsym(_jl_libSuiteSparse, $f_numfree), Void, (Ptr{Void},), Numeric)
93+
94+
end # quote
95+
end # macro
9096

97+
@_jl_umfpack_macro Float64 Int64 :umfpack_dl_symbolic :umfpack_dl_numeric :umfpack_dl_solve :umfpack_dl_free_symbolic :umfpack_dl_free_numeric
98+
@_jl_umfpack_macro Float64 Int32 :umfpack_di_symbolic :umfpack_di_numeric :umfpack_di_solve :umfpack_di_free_symbolic :umfpack_di_free_numeric
9199

92100
## User-callable functions
93-
function (\){Tv<:Float64,Ti<:Int64}(S::SparseMatrixCSC{Tv,Ti}, b::Vector{Tv})
101+
function (\){Tv<:Float64,Ti<:Union(Int64,Int32)}(S::SparseMatrixCSC{Tv,Ti}, b::Vector{Tv})
94102
S = _jl_convert_to_0_based_indexing(S)
95103

96104
symbolic = _jl_umfpack_symbolic(S)
97-
numeric = _jl_umfpack_numeric(S, symbolic[1])
105+
numeric = _jl_umfpack_numeric(S, symbolic)
98106
_jl_umfpack_free_symbolic(S, symbolic)
99-
x = _jl_umfpack_solve(S, b, numeric[1])
107+
x = _jl_umfpack_solve(S, b, numeric)
100108
_jl_umfpack_free_numeric(S, numeric)
101109

102110
S = _jl_convert_to_1_based_indexing(S)

0 commit comments

Comments
 (0)