Skip to content

Commit

Permalink
cholmod: permutation-related fixes
Browse files Browse the repository at this point in the history 
The old tests were run with a permutation that happened to be its own
inverse, and hence were not sufficiently strict.
  • Loading branch information
@timholy
timholy committed Mar 8, 2015
1 parent 079bcb9 commit 137b19e
Show file tree
Hide file tree
Showing 2 changed files with 99 additions and 73 deletions.
51 changes: 33 additions & 18 deletions base/sparse/cholmod.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,11 +202,11 @@ type FactorComponent{Tv,Ti,S} <: AbstractMatrix{Tv}
function FactorComponent(F::Factor{Tv,Ti})
s = unsafe_load(F.p)
if bool(s.is_ll)
S == :L || S == :U || S == :PL || S == :UPt || throw(CHOLMODException(S, " not supported for sparse LLt matrices; try :L, :U, :PL, or :UPt"))
S == :L || S == :U || S == :PtL || S == :UP || throw(CHOLMODException(string(S, " not supported for sparse LLt matrices; try :L, :U, :PtL, or :UP")))
else
S == :L || S == :U || S == :PL || S == :UPt ||
S == :D || S == :LD || S == :DU || S == :PLD || S == :DUPt ||
throw(CHOLMODException(S, " not supported for sparse LDLt matrices; try :L, :U, :PL, :UPt, :D, :LD, :DU, :PLD, or :DUPt"))
S == :L || S == :U || S == :PtL || S == :UP ||
S == :D || S == :LD || S == :DU || S == :PtLD || S == :DUP ||
throw(CHOLMODException(string(S, " not supported for sparse LDLt matrices; try :L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, or :DUP")))
end
new(F)
end
Expand Down Expand Up @@ -850,7 +850,11 @@ function sparse(F::Factor)
end
p = get_perm(F)
if p != [1:s.n;]
A = A[p,p]
pinv = Array(Int, length(p))
for k = 1:length(p)
pinv[p[k]] = k
end
A = A[pinv,pinv]
end
A
end
Expand Down Expand Up @@ -886,14 +890,25 @@ eltype{T<:VTypes}(A::Factor{T}) = T
eltype{T<:VTypes}(A::Sparse{T}) = T

function show(io::IO, F::Factor)
s = unsafe_load(F.p)
println(io, typeof(F))
showfactor(io, F)
end

function show(io::IO, FC::FactorComponent)
println(io, typeof(FC))
showfactor(io, Factor(FC))
end

function showfactor(io::IO, F::Factor)
s = unsafe_load(F.p)
@printf(io, "type: %12s\n", bool(s.is_ll) ? "LLt" : "LDLt")
@printf(io, "method: %10s\n", bool(s.is_super) ? "supernodal" : "simplicial")
@printf(io, "maxnnz: %10d\n", int(s.nzmax))
@printf(io, "nnz: %13d\n", nnz(Sparse(F)))
end



isvalid(A::Dense) = check_dense(A)
isvalid(A::Sparse) = check_sparse(A)
isvalid(A::Factor) = check_factor(A)
Expand Down Expand Up @@ -922,13 +937,13 @@ size(FC::FactorComponent) = size(FC.F)

ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:L}) = FactorComponent{Tv,Ti,:U}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:U}) = FactorComponent{Tv,Ti,:L}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:PL}) = FactorComponent{Tv,Ti,:UPt}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:UPt}) = FactorComponent{Tv,Ti,:PL}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:PtL}) = FactorComponent{Tv,Ti,:UP}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:UP}) = FactorComponent{Tv,Ti,:PtL}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:D}) = FC
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:LD}) = FactorComponent{Tv,Ti,:DU}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:DU}) = FactorComponent{Tv,Ti,:LD}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:PLD}) = FactorComponent{Tv,Ti,:DUPt}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:DUPt}) = FactorComponent{Tv,Ti,:PLD}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:PtLD}) = FactorComponent{Tv,Ti,:DUP}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:DUP}) = FactorComponent{Tv,Ti,:PtLD}(FC.F)

function getindex(A::Dense, i::Integer)
s = unsafe_load(A.p)
Expand Down Expand Up @@ -1181,13 +1196,13 @@ function (\){T,Ti}(L::FactorComponent{T,Ti,:U}, B::Union(Dense,Sparse))
solve(CHOLMOD_Lt, Factor(L), B)
end
# Solve PLx = b and L'P'x=b where A = P*L*L'*P'
function (\){T,Ti}(L::FactorComponent{T,Ti,:PL}, B::Union(Dense,Sparse))
function (\){T,Ti}(L::FactorComponent{T,Ti,:PtL}, B::Union(Dense,Sparse))
F = Factor(L)
solve(CHOLMOD_L, F, solve(CHOLMOD_Pt, F, B)) # Confusingly, CHOLMOD_Pt solves Px = b
solve(CHOLMOD_L, F, solve(CHOLMOD_P, F, B)) # Confusingly, CHOLMOD_P solves P'x = b
end
function (\){T,Ti}(L::FactorComponent{T,Ti,:UPt}, B::Union(Dense,Sparse))
function (\){T,Ti}(L::FactorComponent{T,Ti,:UP}, B::Union(Dense,Sparse))
F = Factor(L)
solve(CHOLMOD_P, F, solve(CHOLMOD_Lt, F, B))
solve(CHOLMOD_Pt, F, solve(CHOLMOD_Lt, F, B))
end
# Solve various equations for A = L*D*L' and A = P*L*D*L'*P'
function (\){T,Ti}(L::FactorComponent{T,Ti,:D}, B::Union(Dense,Sparse))
Expand All @@ -1199,13 +1214,13 @@ end
function (\){T,Ti}(L::FactorComponent{T,Ti,:DU}, B::Union(Dense,Sparse))
solve(CHOLMOD_DLt, Factor(L), B)
end
function (\){T,Ti}(L::FactorComponent{T,Ti,:PLD}, B::Union(Dense,Sparse))
function (\){T,Ti}(L::FactorComponent{T,Ti,:PtLD}, B::Union(Dense,Sparse))
F = Factor(L)
solve(CHOLMOD_LD, F, solve(CHOLMOD_Pt, F, B))
solve(CHOLMOD_LD, F, solve(CHOLMOD_P, F, B))
end
function (\){T,Ti}(L::FactorComponent{T,Ti,:DUPt}, B::Union(Dense,Sparse))
function (\){T,Ti}(L::FactorComponent{T,Ti,:DUP}, B::Union(Dense,Sparse))
F = Factor(L)
solve(CHOLMOD_P, F, solve(CHOLMOD_DLt, F, B))
solve(CHOLMOD_Pt, F, solve(CHOLMOD_DLt, F, B))
end

function (\)(L::FactorComponent, b::Vector)
Expand Down
121 changes: 66 additions & 55 deletions test/sparsedir/cholmod.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -476,95 +476,106 @@ Fs = cholfact(As, [1:3;])
@test_approx_eq sparse(Fs[:L]) Lf
@test_approx_eq sparse(Fs) As
b = rand(3)
@test_approx_eq Fs\b As\b
@test_approx_eq Fs\b Af\b
@test_approx_eq Fs[:UP]\(Fs[:PtL]\b) Af\b
@test_approx_eq Fs[:L]\b Lf\b
@test_approx_eq Fs[:U]\b Lf'\b
@test_approx_eq Fs[:L]'\b Lf'\b
@test_approx_eq Fs[:U]'\b Lf\b
@test_approx_eq Fs[:PL]\b Lf\b
@test_approx_eq Fs[:UPt]\b Lf'\b
@test_approx_eq Fs[:PL]'\b Lf'\b
@test_approx_eq Fs[:UPt]'\b Lf\b
@test_throws MethodError Fs[:D]
@test_throws MethodError Fs[:LD]
@test_throws MethodError Fs[:DU]
@test_throws MethodError Fs[:PLD]
@test_throws MethodError Fs[:DUPt]
@test_approx_eq Fs[:PtL]\b Lf\b
@test_approx_eq Fs[:UP]\b Lf'\b
@test_approx_eq Fs[:PtL]'\b Lf'\b
@test_approx_eq Fs[:UP]'\b Lf\b
@test_throws CHOLMOD.CHOLMODException Fs[:D]
@test_throws CHOLMOD.CHOLMODException Fs[:LD]
@test_throws CHOLMOD.CHOLMODException Fs[:DU]
@test_throws CHOLMOD.CHOLMODException Fs[:PLD]
@test_throws CHOLMOD.CHOLMODException Fs[:DUPt]

# cholfact, with permutation
p = [2,1,3]
p = [2,3,1]
pinv = [3,1,2]
Fs = cholfact(As, p)
@test Fs[:p] == p
Asp = Af[p,p]
Lfp = cholfact(Asp)[:L]
Afp = Af[p,p]
Lfp = cholfact(Afp)[:L]
@test_approx_eq sparse(Fs[:L]) Lfp
@test_approx_eq sparse(Fs) As
b = rand(3)
@test_approx_eq Fs\b As\b
@test_approx_eq Fs[:L]\b[p] Lfp\b[p]
@test_approx_eq Fs[:U]\b[p] Lfp'\b[p]
@test_approx_eq Fs[:L]'\b[p] Lfp'\b[p]
@test_approx_eq Fs[:U]'\b[p] Lfp\b[p]
@test_approx_eq Fs[:PL]\b Lfp\b[p]
@test_approx_eq Fs[:UPt]\b (Lfp'\b)[p]
@test_approx_eq Fs[:PL]'\b (Lfp'\b)[p]
@test_approx_eq Fs[:UPt]'\b Lfp\b[p]
@test_throws MethodError Fs[:D]
@test_throws MethodError Fs[:LD]
@test_throws MethodError Fs[:DU]
@test_throws MethodError Fs[:PLD]
@test_throws MethodError Fs[:DUPt]
@test_approx_eq Fs\b Af\b
@test_approx_eq Fs[:UP]\(Fs[:PtL]\b) Af\b
@test_approx_eq Fs[:L]\b Lfp\b
@test_approx_eq Fs[:U]'\b Lfp\b
@test_approx_eq Fs[:U]\b Lfp'\b
@test_approx_eq Fs[:L]'\b Lfp'\b
@test_approx_eq Fs[:PtL]\b Lfp\b[p]
@test_approx_eq Fs[:UP]\b (Lfp'\b)[pinv]
@test_approx_eq Fs[:PtL]'\b (Lfp'\b)[pinv]
@test_approx_eq Fs[:UP]'\b Lfp\b[p]
@test_throws CHOLMOD.CHOLMODException Fs[:PL]
@test_throws CHOLMOD.CHOLMODException Fs[:UPt]
@test_throws CHOLMOD.CHOLMODException Fs[:D]
@test_throws CHOLMOD.CHOLMODException Fs[:LD]
@test_throws CHOLMOD.CHOLMODException Fs[:DU]
@test_throws CHOLMOD.CHOLMODException Fs[:PLD]
@test_throws CHOLMOD.CHOLMODException Fs[:DUPt]

# ldltfact, no permutation
Fs = ldltfact(As, [1:3;])
@test Fs[:p] == [1:3;]
@test_approx_eq sparse(Fs[:LD]) LDf
@test_approx_eq sparse(Fs) As
b = rand(3)
@test_approx_eq Fs\b As\b
@test_approx_eq Fs\b Af\b
@test_approx_eq Fs[:UP]\(Fs[:PtLD]\b) Af\b
@test_approx_eq Fs[:DUP]\(Fs[:PtL]\b) Af\b
@test_approx_eq Fs[:L]\b L_f\b
@test_approx_eq Fs[:U]\b L_f'\b
@test_approx_eq Fs[:L]'\b L_f'\b
@test_approx_eq Fs[:U]'\b L_f\b
@test_approx_eq Fs[:PL]\b L_f\b
@test_approx_eq Fs[:UPt]\b L_f'\b
@test_approx_eq Fs[:PL]'\b L_f'\b
@test_approx_eq Fs[:UPt]'\b L_f\b
@test_approx_eq Fs[:PtL]\b L_f\b
@test_approx_eq Fs[:UP]\b L_f'\b
@test_approx_eq Fs[:PtL]'\b L_f'\b
@test_approx_eq Fs[:UP]'\b L_f\b
@test_approx_eq Fs[:D]\b D_f\b
@test_approx_eq Fs[:D]'\b D_f\b
@test_approx_eq Fs[:LD]\b D_f\(L_f\b)
@test_approx_eq Fs[:DU]'\b D_f\(L_f\b)
@test_approx_eq Fs[:LD]'\b L_f'\(D_f\b)
@test_approx_eq Fs[:DU]\b L_f'\(D_f\b)
@test_approx_eq Fs[:PLD]\b D_f\(L_f\b)
@test_approx_eq Fs[:DUPt]'\b D_f\(L_f\b)
@test_approx_eq Fs[:PLD]'\b L_f'\(D_f\b)
@test_approx_eq Fs[:DUPt]\b L_f'\(D_f\b)
@test_approx_eq Fs[:PtLD]\b D_f\(L_f\b)
@test_approx_eq Fs[:DUP]'\b D_f\(L_f\b)
@test_approx_eq Fs[:PtLD]'\b L_f'\(D_f\b)
@test_approx_eq Fs[:DUP]\b L_f'\(D_f\b)

# ldltfact, with permutation
p = [2,1,3]
Fs = ldltfact(As, p)
@test Fs[:p] == p
@test_approx_eq sparse(Fs) As
b = rand(3)
Asp = As[p,p]
LDp = sparse(ldltfact(Asp)[:LD])
LDp = sparse(ldltfact(Asp, [1,2,3])[:LD])
# LDp = sparse(Fs[:LD])
Lp, dp = Base.SparseMatrix.CHOLMOD.getLd!(copy(LDp))
Dp = spdiagm(dp)
@test_approx_eq Fs\b As\b
@test_approx_eq Fs[:L]\b[p] Lp\b[p]
@test_approx_eq Fs[:U]\b[p] Lp'\b[p]
@test_approx_eq Fs[:L]'\b[p] Lp'\b[p]
@test_approx_eq Fs[:U]'\b[p] Lp\b[p]
@test_approx_eq Fs[:PL]\b Lp\b[p]
@test_approx_eq Fs[:UPt]\b (Lp'\b)[p]
@test_approx_eq Fs[:PL]'\b (Lp'\b)[p]
@test_approx_eq Fs[:UPt]'\b Lp\b[p]
@test_approx_eq Fs[:LD]\b[p] Dp\(Lp\b[p])
@test_approx_eq Fs[:DU]'\b[p] Dp\(Lp\b[p])
@test_approx_eq Fs[:LD]'\b[p] Lp'\(Dp\b[p])
@test_approx_eq Fs[:DU]\b[p] Lp'\(Dp\b[p])
@test_approx_eq Fs[:PLD]\b Dp\(Lp\b[p])
@test_approx_eq Fs[:DUPt]'\b Dp\(Lp\b[p])
@test_approx_eq Fs[:PLD]'\b (Lp'\(Dp\b))[p]
@test_approx_eq Fs[:DUPt]\b (Lp'\(Dp\b))[p]
@test_approx_eq Fs\b Af\b
@test_approx_eq Fs[:UP]\(Fs[:PtLD]\b) Af\b
@test_approx_eq Fs[:DUP]\(Fs[:PtL]\b) Af\b
@test_approx_eq Fs[:L]\b Lp\b
@test_approx_eq Fs[:U]\b Lp'\b
@test_approx_eq Fs[:L]'\b Lp'\b
@test_approx_eq Fs[:U]'\b Lp\b
@test_approx_eq Fs[:PtL]\b Lp\b[p]
@test_approx_eq Fs[:UP]\b (Lp'\b)[pinv]
@test_approx_eq Fs[:PtL]'\b (Lp'\b)[pinv]
@test_approx_eq Fs[:UP]'\b Lp\b[p]
@test_approx_eq Fs[:LD]\b Dp\(Lp\b)
@test_approx_eq Fs[:DU]'\b Dp\(Lp\b)
@test_approx_eq Fs[:LD]'\b Lp'\(Dp\b)
@test_approx_eq Fs[:DU]\b Lp'\(Dp\b)
@test_approx_eq Fs[:PtLD]\b Dp\(Lp\b[p])
@test_approx_eq Fs[:DUP]'\b Dp\(Lp\b[p])
@test_approx_eq Fs[:PtLD]'\b (Lp'\(Dp\b))[pinv]
@test_approx_eq Fs[:DUP]\b (Lp'\(Dp\b))[pinv]
@test_throws CHOLMOD.CHOLMODException Fs[:DUPt]
@test_throws CHOLMOD.CHOLMODException Fs[:PLD]

0 comments on commit 137b19e

Please sign in to comment.