Skip to content

Commit 6b0e706

Browse files
timholytimholy
committed
Support extracting factors from cholfact(::SparseMatrixCSC)
Also supports extracting factors from ldltfact. A wide variety of "unconventional factors," like :PtL, are supported because of the importance of pivoting in sparse factorizations. This also: - changes the behavior of sparse(cholfact(A)) to return the original matrix, just like full(cholfact(A)) for dense matrices *breaking* - Doesn't export transpose(::Sparse, ::Integer) (that seems like an internal interface) - Adds support for ctranspose
1 parent e9fa25b commit 6b0e706

File tree

2 files changed

+289
-10
lines changed

2 files changed

+289
-10
lines changed

‎base/sparse/cholmod.jl

Lines changed: 171 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ import Base: (*), convert, copy, eltype, getindex, show, size,
66
linearindexing, LinearFast, LinearSlow
77

88
import Base.LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B,
9-
cholfact, cholfact!, det, diag, ishermitian, isposdef,
9+
cholfact, det, diag, ishermitian, isposdef,
1010
issym, ldltfact, logdet
1111

1212
import Base.SparseMatrix: sparse, nnz
@@ -132,10 +132,10 @@ end
132132
# Type definitions #
133133
####################
134134

135-
# The three core data types for CHOLMOD: Dense, Sparse and Factor. CHOLMOD
136-
# manages the memory, so the Julia versions only wrap a pointer to a struct.
137-
# Therefore finalizers should be registret each time a pointer is returned from
138-
# CHOLMOD.
135+
# The three core data types for CHOLMOD: Dense, Sparse and Factor.
136+
# CHOLMOD manages the memory, so the Julia versions only wrap a
137+
# pointer to a struct. Therefore finalizers should be registered each
138+
# time a pointer is returned from CHOLMOD.
139139

140140
# Dense
141141
immutable C_Dense{T<:VTypes}
@@ -273,6 +273,28 @@ type Factor{Tv,Ti} <: Factorization{Tv}
273273
p::Ptr{C_Factor{Tv,Ti}}
274274
end
275275

276+
# FactorComponent, for encoding particular factors from a factorization
277+
type FactorComponent{Tv,Ti,S} <: AbstractMatrix{Tv}
278+
F::Factor{Tv,Ti}
279+
280+
function FactorComponent(F::Factor{Tv,Ti})
281+
s = unsafe_load(F.p)
282+
if s.is_ll != 0
283+
S == :L || S == :U || S == :PtL || S == :UP || throw(CHOLMODException(string(S, " not supported for sparse LLt matrices; try :L, :U, :PtL, or :UP")))
284+
else
285+
S == :L || S == :U || S == :PtL || S == :UP ||
286+
S == :D || S == :LD || S == :DU || S == :PtLD || S == :DUP ||
287+
throw(CHOLMODException(string(S, " not supported for sparse LDLt matrices; try :L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, or :DUP")))
288+
end
289+
new(F)
290+
end
291+
end
292+
function FactorComponent{Tv,Ti}(F::Factor{Tv,Ti}, sym::Symbol)
293+
FactorComponent{Tv,Ti,sym}(F)
294+
end
295+
296+
Factor(FC::FactorComponent) = Factor(FC.F)
297+
276298
#################
277299
# Thin wrappers #
278300
#################
@@ -478,7 +500,7 @@ for Ti in IndexTypes
478500
s
479501
end
480502

481-
function transpose{Tv<:VTypes}(A::Sparse{Tv,$Ti}, values::Integer)
503+
function transpose_{Tv<:VTypes}(A::Sparse{Tv,$Ti}, values::Integer)
482504
s = Sparse(ccall((@cholmod_name("transpose", $Ti),:libcholmod), Ptr{C_Sparse{Tv,$Ti}},
483505
(Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}),
484506
A.p, values, common($Ti)))
@@ -637,6 +659,12 @@ for Ti in IndexTypes
637659
finalizer(f, free!)
638660
f
639661
end
662+
function factorize!{Tv<:VTypes}(A::Sparse{Tv,$Ti}, F::Factor{Tv,$Ti}, cmmn::Vector{UInt8})
663+
@isok ccall((@cholmod_name("factorize", $Ti),:libcholmod), Cint,
664+
(Ptr{C_Sparse{Tv,$Ti}}, Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}),
665+
A.p, F.p, cmmn)
666+
F
667+
end
640668
function factorize_p!{Tv<:VTypes}(A::Sparse{Tv,$Ti}, β::Real, F::Factor{Tv,$Ti}, cmmn::Vector{UInt8})
641669
# note that β is passed as a complex number (double beta[2]),
642670
# but the CHOLMOD manual says that only beta[0] (real part) is used
@@ -691,6 +719,13 @@ for Ti in IndexTypes
691719
end
692720
end
693721

722+
function get_perm(F::Factor)
723+
s = unsafe_load(F.p)
724+
p = pointer_to_array(s.Perm, s.n, false)
725+
p+1
726+
end
727+
get_perm(FC::FactorComponent) = get_perm(Factor(FC))
728+
694729
#########################
695730
# High level interfaces #
696731
#########################
@@ -879,9 +914,37 @@ function sparse{Ti}(A::Sparse{Complex{Float64},Ti}) # Notice! Cannot be type sta
879914
end
880915
return convert(Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}, A)
881916
end
882-
sparse(L::Factor) = sparse(Sparse(L))
917+
function sparse(F::Factor)
918+
s = unsafe_load(F.p)
919+
if s.is_ll != 0
920+
L = Sparse(F)
921+
A = sparse(L*L')
922+
else
923+
LD = sparse(F[:LD])
924+
L, d = getLd!(LD)
925+
A = scale(L, d)*L'
926+
end
927+
p = get_perm(F)
928+
if p != [1:s.n;]
929+
pinv = Array(Int, length(p))
930+
for k = 1:length(p)
931+
pinv[p[k]] = k
932+
end
933+
A = A[pinv,pinv]
934+
end
935+
A
936+
end
937+
883938
sparse(D::Dense) = sparse(Sparse(D))
884939

940+
function sparse{Tv,Ti}(FC::FactorComponent{Tv,Ti,:L})
941+
F = Factor(FC)
942+
s = unsafe_load(F.p)
943+
s.is_ll != 0 || throw(CHOLMODException("sparse: supported only for :LD on LDLt factorizations"))
944+
sparse(Sparse(F))
945+
end
946+
sparse{Tv,Ti}(FC::FactorComponent{Tv,Ti,:LD}) = sparse(Sparse(Factor(FC)))
947+
885948
# Calculate the offset into the stype field of the cholmod_sparse_struct and
886949
# change the value
887950
let offidx=findfirst(fieldnames(C_Sparse) .== :stype)
@@ -905,14 +968,25 @@ eltype{T<:VTypes}(A::Sparse{T}) = T
905968
nnz(F::Factor) = nnz(Sparse(F))
906969

907970
function show(io::IO, F::Factor)
908-
s = unsafe_load(F.p)
909971
println(io, typeof(F))
972+
showfactor(io, F)
973+
end
974+
975+
function show(io::IO, FC::FactorComponent)
976+
println(io, typeof(FC))
977+
showfactor(io, Factor(FC))
978+
end
979+
980+
function showfactor(io::IO, F::Factor)
981+
s = unsafe_load(F.p)
910982
@printf(io, "type: %12s\n", s.is_ll!=0 ? "LLt" : "LDLt")
911983
@printf(io, "method: %10s\n", s.is_super!=0 ? "supernodal" : "simplicial")
912984
@printf(io, "maxnnz: %10d\n", Int(s.nzmax))
913985
@printf(io, "nnz: %13d\n", nnz(F))
914986
end
915987

988+
989+
916990
isvalid(A::Dense) = check_dense(A)
917991
isvalid(A::Sparse) = check_sparse(A)
918992
isvalid(A::Factor) = check_factor(A)
@@ -936,7 +1010,22 @@ function size(F::Factor, i::Integer)
9361010
return 1
9371011
end
9381012

1013+
9391014
linearindexing(::Dense) = LinearFast()
1015+
1016+
size(FC::FactorComponent, i::Integer) = size(FC.F, i)
1017+
size(FC::FactorComponent) = size(FC.F)
1018+
1019+
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:L}) = FactorComponent{Tv,Ti,:U}(FC.F)
1020+
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:U}) = FactorComponent{Tv,Ti,:L}(FC.F)
1021+
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:PtL}) = FactorComponent{Tv,Ti,:UP}(FC.F)
1022+
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:UP}) = FactorComponent{Tv,Ti,:PtL}(FC.F)
1023+
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:D}) = FC
1024+
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:LD}) = FactorComponent{Tv,Ti,:DU}(FC.F)
1025+
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:DU}) = FactorComponent{Tv,Ti,:LD}(FC.F)
1026+
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:PtLD}) = FactorComponent{Tv,Ti,:DUP}(FC.F)
1027+
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:DUP}) = FactorComponent{Tv,Ti,:PtLD}(FC.F)
1028+
9401029
function getindex(A::Dense, i::Integer)
9411030
s = unsafe_load(A.p)
9421031
0 < i <= s.nrow*s.ncol || throw(BoundsError())
@@ -957,6 +1046,27 @@ function getindex{T}(A::Sparse{T}, i0::Integer, i1::Integer)
9571046
((r1 > r2) || (unsafe_load(s.i, r1) + 1 != i0)) ? zero(T) : unsafe_load(s.x, r1)
9581047
end
9591048

1049+
function getindex(F::Factor, sym::Symbol)
1050+
sym == :p && return get_perm(F)
1051+
FactorComponent(F, sym)
1052+
end
1053+
1054+
function getLd!(S::SparseMatrixCSC)
1055+
d = Array(eltype(S), size(S, 1))
1056+
fill!(d, 0)
1057+
col = 1
1058+
for k = 1:length(S.nzval)
1059+
while k >= S.colptr[col+1]
1060+
col += 1
1061+
end
1062+
if S.rowval[k] == col
1063+
d[col] = S.nzval[k]
1064+
S.nzval[k] = 1
1065+
end
1066+
end
1067+
S, d
1068+
end
1069+
9601070
## Multiplication
9611071
(*)(A::Sparse, B::Sparse) = ssmult(A, B, 0, true, true)
9621072
(*)(A::Sparse, B::Dense) = sdmult!(A, false, 1., 0., B, zeros(size(A, 1), size(B, 2)))
@@ -966,7 +1076,7 @@ function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti})
9661076
cm = common(Ti)
9671077

9681078
if !is(A,B)
969-
aa1 = transpose(B, 2)
1079+
aa1 = transpose_(B, 2)
9701080
## result of ssmult will have stype==0, contain numerical values and be sorted
9711081
return ssmult(A, aa1, 0, true, true)
9721082
end
@@ -984,7 +1094,7 @@ function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti})
9841094
end
9851095

9861096
function Ac_mul_B(A::Sparse, B::Sparse)
987-
aa1 = transpose(A, 2)
1097+
aa1 = transpose_(A, 2)
9881098
if is(A,B)
9891099
return A_mul_Bc(aa1, aa1)
9901100
end
@@ -1071,6 +1181,57 @@ update!{T<:VTypes}(F::Factor{T}, A::SparseMatrixCSC{T}; kws...) = update!(F, Spa
10711181

10721182
## Solvers
10731183

1184+
for (T, f) in ((:Dense, :solve), (:Sparse, :spsolve))
1185+
@eval begin
1186+
# Solve Lx = b and L'x=b where A = L*L'
1187+
function (\){T,Ti}(L::FactorComponent{T,Ti,:L}, B::$T)
1188+
($f)(CHOLMOD_L, Factor(L), B)
1189+
end
1190+
function (\){T,Ti}(L::FactorComponent{T,Ti,:U}, B::$T)
1191+
($f)(CHOLMOD_Lt, Factor(L), B)
1192+
end
1193+
# Solve PLx = b and L'P'x=b where A = P*L*L'*P'
1194+
function (\){T,Ti}(L::FactorComponent{T,Ti,:PtL}, B::$T)
1195+
F = Factor(L)
1196+
($f)(CHOLMOD_L, F, ($f)(CHOLMOD_P, F, B)) # Confusingly, CHOLMOD_P solves P'x = b
1197+
end
1198+
function (\){T,Ti}(L::FactorComponent{T,Ti,:UP}, B::$T)
1199+
F = Factor(L)
1200+
($f)(CHOLMOD_Pt, F, ($f)(CHOLMOD_Lt, F, B))
1201+
end
1202+
# Solve various equations for A = L*D*L' and A = P*L*D*L'*P'
1203+
function (\){T,Ti}(L::FactorComponent{T,Ti,:D}, B::$T)
1204+
($f)(CHOLMOD_D, Factor(L), B)
1205+
end
1206+
function (\){T,Ti}(L::FactorComponent{T,Ti,:LD}, B::$T)
1207+
($f)(CHOLMOD_LD, Factor(L), B)
1208+
end
1209+
function (\){T,Ti}(L::FactorComponent{T,Ti,:DU}, B::$T)
1210+
($f)(CHOLMOD_DLt, Factor(L), B)
1211+
end
1212+
function (\){T,Ti}(L::FactorComponent{T,Ti,:PtLD}, B::$T)
1213+
F = Factor(L)
1214+
($f)(CHOLMOD_LD, F, ($f)(CHOLMOD_P, F, B))
1215+
end
1216+
function (\){T,Ti}(L::FactorComponent{T,Ti,:DUP}, B::$T)
1217+
F = Factor(L)
1218+
($f)(CHOLMOD_Pt, F, ($f)(CHOLMOD_DLt, F, B))
1219+
end
1220+
end
1221+
end
1222+
1223+
function (\)(L::FactorComponent, b::Vector)
1224+
reshape(convert(Matrix, L\Dense(b)), length(b))
1225+
end
1226+
function (\)(L::FactorComponent, B::Matrix)
1227+
convert(Matrix, L\Dense(B))
1228+
end
1229+
function (\)(L::FactorComponent, B::SparseMatrixCSC)
1230+
sparse(L\Sparse(B,0))
1231+
end
1232+
1233+
Ac_ldiv_B(L::FactorComponent, B) = ctranspose(L)\B
1234+
10741235
(\)(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B)
10751236
(\)(L::Factor, b::Vector) = reshape(convert(Matrix, solve(CHOLMOD_A, L, Dense(b))), length(b))
10761237
(\)(L::Factor, B::Matrix) = convert(Matrix, solve(CHOLMOD_A, L, Dense(B)))

0 commit comments

Comments
 (0)