From 3ec69dbf8bd9403d421a2c36f90d26dc6d908643 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sat, 7 Feb 2015 16:45:24 -0500 Subject: [PATCH] Let CHOLMOD handle memory allocation. Increase test coverage and remove some unused methods. --- base/linalg/symmetric.jl | 2 +- base/sparse/cholmod.jl | 1333 +++++++++++++++++++------------------- test/sparse/cholmod.jl | 246 ++++++- 3 files changed, 899 insertions(+), 682 deletions(-) diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index aae84206ea6fa..770989812c8ec 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -13,7 +13,7 @@ typealias HermOrSym{T,S} Union(Hermitian{T,S}, Symmetric{T,S}) typealias RealHermSymComplexHerm{T<:Real,S} Union(Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}) size(A::HermOrSym, args...) = size(A.data, args...) -getindex(A::Symmetric, i::Integer) = ((q, r) = divrem(i - 1, size(A, 1)); A[r + 1, q + 1]) +getindex(A::HermOrSym, i::Integer) = ((q, r) = divrem(i - 1, size(A, 1)); A[r + 1, q + 1]) getindex(A::Symmetric, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.data, i, j) : getindex(A.data, j, i) getindex(A::Hermitian, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.data, i, j) : conj(getindex(A.data, j, i)) full(A::Symmetric) = copytri!(copy(A.data), A.uplo) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 8e9afa0293144..b657e1da3f66d 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -1,23 +1,24 @@ module CHOLMOD -import Base: (*), convert, copy, ctranspose, eltype, getindex, hcat, isvalid, show, size, - sort!, transpose, vcat +import Base: (*), convert, copy, eltype, getindex, show, size import Base.LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B, - A_ldiv_B!, cholfact, cholfact!, copy, det, diag, - full, ishermitian, isposdef!, issym, ldltfact, logdet, norm, scale, scale! + A_ldiv_B!, cholfact, cholfact!, det, diag, full, ishermitian, isposdef, + issym, isvalid, ldltfact, logdet -import Base.SparseMatrix: findnz, sparse +import Base.SparseMatrix: sparse export Dense, Factor, - Sparse, - etree, - sparse! + Sparse using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, increment!, indtype, decrement, decrement! +######### +# Setup # +######### + include("cholmod_h.jl") macro isok(A) @@ -84,46 +85,84 @@ ccall((:jl_cholmod_common_offsets, :libsuitesparse_wrapper), const common_final_ll = (1:4) + cholmod_com_offsets[7] const common_print = (1:4) + cholmod_com_offsets[13] const common_itype = (1:4) + cholmod_com_offsets[18] - -### there must be an easier way but at least this works. -# function Common(aa::Array{UInt8,1}) -# typs = Common.types -# sz = map(sizeof, typs) -# args = map(i->reinterpret(typs[i], aa[cholmod_com_offsets[i] + (1:sz[i])])[1], 1:length(sz)) -# eval(Expr(:call, unshift!(args, :Common), Any)) -# end +const common_dtype = (1:4) + cholmod_com_offsets[19] function set_print_level(cm::Array{UInt8}, lev::Integer) cm[common_print] = reinterpret(UInt8, [int32(lev)]) end -## cholmod_dense pointers passed to or returned from C functions are of Julia type -## Ptr{c_Dense}. The Dense type contains a c_Dense object and other -## fields then ensure the memory pointed to is freed when it should be and not before. -type c_Dense{T<:VTypes} - m::Int - n::Int - nzmax::Int - lda::Int - xpt::Ptr{T} - zpt::Ptr{Void} +#################### +# Type definitions # +#################### + +# The three core data types for CHOLMOD: Dense, Sparse and Factor. CHOLMOD manages the memory, so the Julia versions only wrap a pointer to a struct. Therefore finalizers should be registret each time a pointer is returned from CHOLMOD. + +# Dense +type C_Dense{T<:VTypes} + nrow::Csize_t + ncol::Csize_t + nzmax::Csize_t + d::Csize_t + x::Ptr{T} + z::Ptr{Void} xtype::Cint dtype::Cint end -type Dense{T<:VTypes} <: AbstractMatrix{T} - c::c_Dense - mat::Matrix{T} +type Dense{T<:VTypes} <: DenseMatrix{T} + p::Ptr{C_Dense{T}} +end + +# Sparse +type C_Sparse{Tv<:VTypes,Ti<:ITypes} + nrow::Csize_t + ncol::Csize_t + nzmax::Csize_t + p::Ptr{Ti} + i::Ptr{Ti} + nz::Ptr{Ti} + x::Ptr{Tv} + z::Ptr{Void} + stype::Cint + itype::Cint + xtype::Cint + dtype::Cint + sorted::Cint + packed::Cint end +# Corresponds to the exact definition of cholmod_sparse_struct in the library. Useful when reading matrices of unknown type from files as in cholmod_read_sparse +type C_SparseVoid + nrow::Csize_t + ncol::Csize_t + nzmax::Csize_t + p::Ptr{Void} + i::Ptr{Void} + nz::Ptr{Void} + x::Ptr{Void} + z::Ptr{Void} + stype::Cint + itype::Cint + xtype::Cint + dtype::Cint + sorted::Cint + packed::Cint +end + +type Sparse{Tv<:VTypes,Ti<:ITypes} <: AbstractSparseMatrix{Tv,Ti} + p::Ptr{C_Sparse{Tv,Ti}} +end + +# Factor + if version >= v"2.1.0" # CHOLMOD version 2.1.0 or later - type c_Factor{Tv<:VTypes,Ti<:ITypes} - n::Int - minor::Int + type C_Factor{Tv<:VTypes,Ti<:ITypes} + n::Csize_t + minor::Csize_t Perm::Ptr{Ti} ColCount::Ptr{Ti} IPerm::Ptr{Ti} # this pointer was added in verison 2.1.0 - nzmax::Int + nzmax::Csize_t p::Ptr{Ti} i::Ptr{Ti} x::Ptr{Tv} @@ -131,11 +170,11 @@ if version >= v"2.1.0" # CHOLMOD version 2.1.0 or later nz::Ptr{Ti} next::Ptr{Ti} prev::Ptr{Ti} - nsuper::Int - ssize::Int - xsize::Int - maxcsize::Int - maxesize::Int + nsuper::Csize_t + ssize::Csize_t + xsize::Csize_t + maxcsize::Csize_t + maxesize::Csize_t super::Ptr{Ti} pi::Ptr{Ti} px::Ptr{Ti} @@ -148,51 +187,13 @@ if version >= v"2.1.0" # CHOLMOD version 2.1.0 or later xtype::Cint dtype::Cint end - - type Factor{Tv<:VTypes,Ti<:ITypes} - c::c_Factor{Tv,Ti} - Perm::Vector{Ti} - ColCount::Vector{Ti} - IPerm::Vector{Ti} - p::Vector{Ti} - i::Vector{Ti} - x::Vector{Tv} - nz::Vector{Ti} - next::Vector{Ti} - prev::Vector{Ti} - super::Vector{Ti} - pi::Vector{Ti} - px::Vector{Ti} - s::Vector{Ti} - end - - function Factor{Tv<:VTypes,Ti<:ITypes}(cp::Ptr{c_Factor{Tv,Ti}}) - cfp = unsafe_load(cp) - Perm = pointer_to_array(cfp.Perm, (cfp.n,), true) - ColCount = pointer_to_array(cfp.ColCount, (cfp.n,), true) - IPerm = pointer_to_array(cfp.IPerm, (cfp.IPerm == C_NULL ? 0 : cfp.n + 1,), true) - p = pointer_to_array(cfp.p, (cfp.p == C_NULL ? 0 : cfp.n + 1,), true) - i = pointer_to_array(cfp.i, (cfp.i == C_NULL ? 0 : cfp.nzmax,), true) - x = pointer_to_array(cfp.x, (cfp.x == C_NULL ? 0 : max(cfp.nzmax,cfp.xsize),), true) - nz = pointer_to_array(cfp.nz, (cfp.nz == C_NULL ? 0 : cfp.n,), true) - next = pointer_to_array(cfp.next, (cfp.next == C_NULL ? 0 : cfp.n + 2,), true) - prev = pointer_to_array(cfp.prev, (cfp.prev == C_NULL ? 0 : cfp.n + 2,), true) - super = pointer_to_array(cfp.super, (cfp.super == C_NULL ? 0 : cfp.nsuper + 1,), true) - pi = pointer_to_array(cfp.pi, (cfp.pi == C_NULL ? 0 : cfp.nsuper + 1,), true) - px = pointer_to_array(cfp.px, (cfp.px == C_NULL ? 0 : cfp.nsuper + 1,), true) - s = pointer_to_array(cfp.s, (cfp.s == C_NULL ? 0 : cfp.ssize + 1,), true) - cf = Factor{Tv,Ti}(cfp, Perm, ColCount, IPerm, p, i, x, nz, next, prev, - super, pi, px, s) - c_free(cp) - cf - end else - type c_Factor{Tv<:VTypes,Ti<:ITypes} - n::Int - minor::Int + type C_Factor{Tv<:VTypes,Ti<:ITypes} + n::Csize_t + minor::Csize_t Perm::Ptr{Ti} ColCount::Ptr{Ti} - nzmax::Int + nzmax::Csize_t p::Ptr{Ti} i::Ptr{Ti} x::Ptr{Tv} @@ -200,11 +201,11 @@ else nz::Ptr{Ti} next::Ptr{Ti} prev::Ptr{Ti} - nsuper::Int - ssize::Int - xsize::Int - maxcsize::Int - maxesize::Int + nsuper::Csize_t + ssize::Csize_t + xsize::Csize_t + maxcsize::Csize_t + maxesize::Csize_t super::Ptr{Ti} pi::Ptr{Ti} px::Ptr{Ti} @@ -217,455 +218,273 @@ else xtype::Cint dtype::Cint end - - type Factor{Tv<:VTypes,Ti<:ITypes} - c::c_Factor{Tv,Ti} - Perm::Vector{Ti} - ColCount::Vector{Ti} - p::Vector{Ti} - i::Vector{Ti} - x::Vector{Tv} - nz::Vector{Ti} - next::Vector{Ti} - prev::Vector{Ti} - super::Vector{Ti} - pi::Vector{Ti} - px::Vector{Ti} - s::Vector{Ti} - end - - function Factor{Tv<:VTypes,Ti<:ITypes}(cp::Ptr{c_Factor{Tv,Ti}}) - cfp = unsafe_load(cp) - Perm = pointer_to_array(cfp.Perm, (cfp.n,), true) - ColCount = pointer_to_array(cfp.ColCount, (cfp.n,), true) - p = pointer_to_array(cfp.p, (cfp.p == C_NULL ? 0 : cfp.n + 1,), true) - i = pointer_to_array(cfp.i, (cfp.i == C_NULL ? 0 : cfp.nzmax,), true) - x = pointer_to_array(cfp.x, (cfp.x == C_NULL ? 0 : max(cfp.nzmax,cfp.xsize),), true) - nz = pointer_to_array(cfp.nz, (cfp.nz == C_NULL ? 0 : cfp.n,), true) - next = pointer_to_array(cfp.next, (cfp.next == C_NULL ? 0 : cfp.n + 2,), true) - prev = pointer_to_array(cfp.prev, (cfp.prev == C_NULL ? 0 : cfp.n + 2,), true) - super = pointer_to_array(cfp.super, (cfp.super == C_NULL ? 0 : cfp.nsuper + 1,), true) - pi = pointer_to_array(cfp.pi, (cfp.pi == C_NULL ? 0 : cfp.nsuper + 1,), true) - px = pointer_to_array(cfp.px, (cfp.px == C_NULL ? 0 : cfp.nsuper + 1,), true) - s = pointer_to_array(cfp.s, (cfp.s == C_NULL ? 0 : cfp.ssize + 1,), true) - cf = Factor{Tv,Ti}(cfp, Perm, ColCount, p, i, x, nz, next, prev, - super, pi, px, s) - c_free(cp) - cf - end -end - -type c_Sparse{Tv<:VTypes,Ti<:ITypes} - m::Csize_t - n::Csize_t - nzmax::Csize_t - ppt::Ptr{Ti} - ipt::Ptr{Ti} - nzpt::Ptr{Ti} - xpt::Ptr{Tv} - zpt::Ptr{Void} - stype::Cint - itype::Cint - xtype::Cint - dtype::Cint - sorted::Cint - packed::Cint -end - -# Corresponds to the exact definition of cholmod_sparse_struct in the library. Useful when reading matrices of unknown type from files as in cholmod_read_sparse -type c_SparseVoid - m::Csize_t - n::Csize_t - nzmax::Csize_t - ppt::Ptr{Void} - ipt::Ptr{Void} - nzpt::Ptr{Void} - xpt::Ptr{Void} - zpt::Ptr{Void} - stype::Cint - itype::Cint - xtype::Cint - dtype::Cint - sorted::Cint - packed::Cint end -type Sparse{Tv<:VTypes,Ti<:ITypes} <: AbstractSparseMatrix{Tv,Ti} - c::c_Sparse{Tv,Ti} - colptr0::Vector{Ti} - rowval0::Vector{Ti} - nzval::Vector{Tv} +type Factor{Tv,Ti} <: Factorization{Tv} + p::Ptr{C_Factor{Tv,Ti}} end -eltype{T<:VTypes}(A::Dense{T}) = T -eltype{T<:VTypes}(A::Factor{T}) = T -eltype{T<:VTypes}(A::Sparse{T}) = T +################# +# Thin wrappers # +################# -## The Dense constructor does not copy the contents, which is generally what you -## want as most uses of Dense objects are read-only. -function Dense{T<:VTypes}(A::VecOrMat{T}) # uses the memory from Julia - m, n = size(A,1), size(A,2) - Dense(c_Dense{T}(m, n, m*n, stride(A, 2), convert(Ptr{T}, A), C_NULL, xtyp(T), dtyp(T)), - length(size(A)) == 2 ? A : reshape(A, (m, n))) -end -function Dense{T<:VTypes}(c::Ptr{c_Dense{T}}) - cp = unsafe_load(c) - ## the true in the call to pointer_to_array means Julia will free the memory - val = Dense(cp, pointer_to_array(cp.xpt, (cp.m,cp.n), true)) - c_free(c) - val -end - -function Sparse{Tv<:VTypes,Ti<:ITypes}(colpt::Vector{Ti}, - rowval::Vector{Ti}, - nzval::Vector{Tv}, - m::Integer, - n::Integer, - stype::Signed) - bb = colpt[1] - - if bool(bb) # one-based - colpt0 = decrement(colpt) - rowval0 = decrement(rowval) - else # zero based - colpt0 = colpt - rowval0 = rowval - end - nz = colpt0[end] - - it = ityp(Ti) - cs = Sparse(c_Sparse{Tv,Ti}(m, - n, - int(nz), - convert(Ptr{Ti}, colpt0), - convert(Ptr{Ti}, rowval0), - C_NULL, - convert(Ptr{Tv}, nzval), - C_NULL, - int32(stype), - ityp(Ti), - xtyp(Tv), - dtyp(Tv), - TRUE, - TRUE), - colpt0, - rowval0, - nzval) - - @isok isvalid(cs) - - return cs -end -function Sparse{Tv<:VTypes,Ti<:ITypes}(A::SparseMatrixCSC{Tv,Ti}, stype::Signed) - Sparse(A.colptr, A.rowval, A.nzval, size(A,1), size(A,2), stype) -end -function Sparse(A::SparseMatrixCSC) - AC = Sparse(A, 0) - i = symmetry(AC, 0)[1] # Check for Hermitianity - if (isreal(A) && (i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG)) || # real case - (i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG) # complex case - AC.c.stype = -1 - end - AC -end -Sparse{Ti<:ITypes}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) -Sparse{Ti<:ITypes}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) - - -# Constructs a Sparse object from a pointer to a cholmod_sparse_struct created by CHOLMOD. Julia takes ownership of the memory allocated by CHOLMOD and the pointer can therefore be freed without memory leek. -function Sparse{Tv<:VTypes,Ti<:ITypes}(cp::Ptr{c_Sparse{Tv,Ti}}) - csp = unsafe_load(cp) - - # Take control of the memory allocated by CHOLMOD - colptr0 = pointer_to_array(csp.ppt, (csp.n + 1,), true) - rowval0 = pointer_to_array(csp.ipt, (csp.nzmax,), true) - nzval = pointer_to_array(csp.xpt, (csp.nzmax,), true) - - cms = Sparse{Tv,Ti}(csp, colptr0, rowval0, nzval) - - # Free memory - csp.packed == 1 && c_free(csp.nzpt) # Julia's Sparse is not using unpacked storage - c_free(cp) +# Dense wrappers +## Note! Integer type defaults to Cint, but this is actually not necessary, but making this a choice would require another type parameter in the Dense type - cms +### cholmod_core_h ### +function allocate_dense(nrow::Integer, ncol::Integer, d::Integer, ::Type{Float64}) + d = Dense(ccall((:cholmod_allocate_dense, :libcholmod), Ptr{C_Dense{Float64}}, + (Csize_t, Csize_t, Csize_t, Cint, Ptr{Void}), + nrow, ncol, d, REAL, common(Cint))) + finalizer(d, free!) + d end -# Useful when reading in files, but not type stable -function Sparse(p::Ptr{c_SparseVoid}) - s = unsafe_load(p) - - # Check integer type - if s.itype == INT - Ti = Cint - elseif s.itype == INTLONG - throw(CHOLMODException("the value of itype was $s.itype. This combination of integer types shouldn't happen. Please submit a bug report.")) - elseif s.itype == LONG - Ti = SuiteSparse_long - else - throw(CHOLMODException("illegal value of itype: $s.itype")) - end - - # Check for double or single precision - if s.dtype == DOUBLE - Tv = Float64 - elseif s.dtype == SINGLE - Tv = Float32 - else - throw(CHOLMODException("illegal value of dtype: $s.dtype")) - end - - # Check for real or complex - if s.xtype == COMPLEX - Tv = Complex{Tv} - elseif s.xtype != REAL - throw(CHOLMODException("illegal value of xtype: $s.xtype")) - end - - # Take control of the memory allocated by CHOLMOD - colptr0 = pointer_to_array(convert(Ptr{Ti}, s.ppt), (s.n + 1,), true) - rowval0 = pointer_to_array(convert(Ptr{Ti}, s.ipt), (s.nzmax,), true) - nzval = pointer_to_array(convert(Ptr{Tv}, s.xpt), (s.nzmax,), true) - - r = Sparse{Tv,Ti}(c_Sparse{Tv,Ti}(s.m, - s.n, - s.nzmax, - pointer(colptr0), - pointer(rowval0), - C_NULL, - pointer(nzval), - C_NULL, - s.stype, - s.itype, - s.xtype, - s.dtype, - s.sorted, - s.packed), - colptr0, rowval0, nzval) - - # Free memory - s.packed == 1 && c_free(s.nzpt) # Julia's Sparse is not using unpacked storage - c_free(p) - - r +function allocate_dense(nrow::Integer, ncol::Integer, d::Integer, ::Type{Complex{Float64}}) + d = Dense(ccall((:cholmod_allocate_dense, :libcholmod), Ptr{C_Dense{Complex{Float64}}}, + (Csize_t, Csize_t, Csize_t, Cint, Ptr{Void}), + nrow, ncol, d, COMPLEX, common(Cint))) + finalizer(d, free!) + d end -Sparse{Tv<:VTypes}(D::Dense{Tv}) = Sparse(D,1) # default Ti is Int +free_dense{T}(p::Ptr{C_Dense{T}}) = ccall((:cholmod_free_dense, :libcholmod), Cint, (Ptr{Ptr{C_Dense{T}}}, Ptr{Void}), &p, common(Cint)) -# Dense wrappers -### cholmod_core_h ### function zeros{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) - Dense(ccall((:cholmod_zeros, :libcholmod), Ptr{c_Dense{T}}, + d = Dense(ccall((:cholmod_zeros, :libcholmod), Ptr{C_Dense{T}}, (Csize_t, Csize_t, Cint, Ptr{UInt8}), - m, n, xtyp(T), common(Int32))) + m, n, xtyp(T), common(Cint))) + finalizer(d, free!) + d end zeros(m::Integer, n::Integer) = zeros(m, n, Float64) function ones{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) - Dense(ccall((:cholmod_ones, :libcholmod), Ptr{c_Dense{T}}, + d = Dense(ccall((:cholmod_ones, :libcholmod), Ptr{C_Dense{T}}, (Csize_t, Csize_t, Cint, Ptr{UInt8}), - m, n, xtyp(T), common(Int32))) + m, n, xtyp(T), common(Cint))) + finalizer(d, free!) + d end ones(m::Integer, n::Integer) = ones(m, n, Float64) function eye{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) - Dense(ccall((:cholmod_eye, :libcholmod), Ptr{c_Dense{T}}, + d = Dense(ccall((:cholmod_eye, :libcholmod), Ptr{C_Dense{T}}, (Csize_t, Csize_t, Cint, Ptr{UInt8}), - m, n, xtyp(T), common(Int32))) + m, n, xtyp(T), common(Cint))) + finalizer(d, free!) + d end eye(m::Integer, n::Integer) = eye(m, n, Float64) eye(n::Integer) = eye(n, n, Float64) function copy_dense{Tv<:VTypes}(A::Dense{Tv}) - Dense(ccall((:cholmod_copy_dense, :libcholmod), Ptr{c_Dense{Tv}}, - (Ptr{c_Dense{Tv}}, Ptr{UInt8}), - &A.c, common(Cint))) + d = Dense(ccall((:cholmod_copy_dense, :libcholmod), Ptr{C_Dense{Tv}}, + (Ptr{C_Dense{Tv}}, Ptr{UInt8}), + A.p, common(Cint))) + finalizer(d, free!) + d end ### cholmod_matrixops.h ### function norm_dense{Tv<:VTypes}(D::Dense{Tv}, p::Integer) + s = unsafe_load(D.p) + if p == 2 + if s.ncol > 1 + throw(ArgumentError("2 norm only supported when matrix has one column")) + end + elseif p != 0 && p != 1 + throw(ArgumentError("second argument must be either 0 (Inf norm), 1, or 2")) + end ccall((:cholmod_norm_dense, :libcholmod), Cdouble, - (Ptr{c_Dense{Tv}}, Cint, Ptr{UInt8}), - &D.c, p, common(Int32)) + (Ptr{C_Dense{Tv}}, Cint, Ptr{UInt8}), + D.p, p, common(Cint)) end ### cholmod_check.h ### function check_dense{T<:VTypes}(A::Dense{T}) bool(ccall((:cholmod_check_dense, :libcholmod), Cint, - (Ptr{c_Dense{T}}, Ptr{UInt8}), - &A.c, common(Cint))) + (Ptr{C_Dense{T}}, Ptr{UInt8}), + A.p, common(Cint))) end # Non-Dense wrappers (which all depend on IType) for Ti in IndexTypes @eval begin - # This method is probably not necessary as memory allocated by CHOLMOD is usually overtaken by Julia when cosntructing Julia instances from pointers. - function free_sparse{Tv<:VTypes}(ptr::Ptr{c_Sparse{Tv,$Ti}}) - cm = common($Ti) - @isok ccall((@cholmod_name "free_sparse" $Ti - ,:libcholmod), Cint, - (Ptr{Ptr{c_Sparse{Tv,$Ti}}}, Ptr{UInt8}), - &ptr, cm) + ### cholmod_core.h ### + function allocate_sparse(nrow::Integer, ncol::Integer, nzmax::Integer, sorted::Bool, packed::Bool, stype::Integer, ::Type{Float64}, ::Type{$Ti}) + s = Sparse(ccall((@cholmod_name "allocate_sparse" $Ti + , :libcholmod), Ptr{C_Sparse{Float64,$Ti}}, + (Csize_t, Csize_t, Csize_t, Cint, + Cint, Cint, Cint, Ptr{Void}), + nrow, ncol, nzmax, sorted, + packed, stype, REAL, common($Ti))) + finalizer(s, free!) + s end - - function (*){Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}) - Sparse(ccall((@cholmod_name "ssmult" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Cint, Ptr{UInt8}), - &A.c, &B.c, 0, true, true, common($Ti))) + function allocate_sparse(nrow::Integer, ncol::Integer, nzmax::Integer, sorted::Bool, packed::Bool, stype::Integer, ::Type{Complex{Float64}}, ::Type{$Ti}) + s = Sparse(ccall((@cholmod_name "allocate_sparse" $Ti + , :libcholmod), Ptr{C_Sparse{Complex{Float64},$Ti}}, + (Csize_t, Csize_t, Csize_t, Cint, + Cint, Cint, Cint, Ptr{Void}), + nrow, ncol, nzmax, sorted, + packed, stype, COMPLEX, common($Ti))) + finalizer(s, free!) + s + end + function free_sparse{Tv<:VTypes}(ptr::Ptr{C_Sparse{Tv,$Ti}}) + @isok ccall((@cholmod_name "free_sparse" $Ti + , :libcholmod), Cint, + (Ptr{Ptr{C_Sparse{Tv,$Ti}}}, Ptr{UInt8}), + &ptr, common($Ti)) end - function transpose{Tv<:VTypes}(A::Sparse{Tv,$Ti}, values::Integer) - cm = common($Ti) - Sparse(ccall((@cholmod_name "transpose" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, values, cm)) + function free_sparse(ptr::Ptr{C_SparseVoid}) + @isok ccall((@cholmod_name "free_sparse" $Ti + , :libcholmod), Cint, + (Ptr{Ptr{C_SparseVoid}}, Ptr{UInt8}), + &ptr, common($Ti)) end - function ssmult{Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, stype::Integer, values::Integer, sorted::Integer) - cm = common($Ti) - Sparse(ccall((@cholmod_name "ssmult" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Cint, Ptr{UInt8}), - &A.c, &B.c, stype, values, sorted, cm)) + function free_factor{Tv<:VTypes}(ptr::Ptr{C_Factor{Tv,$Ti}}) + @isok ccall((@cholmod_name "free_factor" $Ti + , :libcholmod), Cint, + (Ptr{Ptr{C_Factor{Tv,$Ti}}}, Ptr{Void}), + &ptr, common($Ti)) end function aat{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, fset::Vector{$Ti}, mode::Integer) - cm = common($Ti) - Sparse(ccall((@cholmod_name "aat" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{$Ti}, Csize_t, Cint, Ptr{UInt8}), - &A.c, fset, length(fset), mode, cm)) + s = Sparse(ccall((@cholmod_name "aat" $Ti + , :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{$Ti}, Csize_t, Cint, Ptr{UInt8}), + A.p, fset, length(fset), mode, common($Ti))) + finalizer(s, free!) + s end function copy{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, stype::Integer, mode::Integer) - cm = common($Ti) - Sparse(ccall((@cholmod_name "copy" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), - &A.c, 0, 1, cm)) + s = Sparse(ccall((@cholmod_name "copy" $Ti + , :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), + A.p, 0, 1, common($Ti))) + finalizer(s, free!) + s end function sparse_to_dense{Tv<:VTypes}(A::Sparse{Tv,$Ti}) - Dense(ccall((@cholmod_name "sparse_to_dense" $Ti - ,:libcholmod), Ptr{c_Dense{Tv}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, common($Ti))) + d = Dense(ccall((@cholmod_name "sparse_to_dense" $Ti + ,:libcholmod), Ptr{C_Dense{Tv}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, common($Ti))) + finalizer(d, free!) + d end - function dense_to_sparse{Tv<:VTypes}(D::Dense{Tv},i::$Ti) - Sparse(ccall((@cholmod_name "dense_to_sparse" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Dense{Tv}}, Cint, Ptr{UInt8}), - &D.c, true, common($Ti))) + function dense_to_sparse{Tv<:VTypes}(D::Dense{Tv}, ::Type{$Ti}) + s = Sparse(ccall((@cholmod_name "dense_to_sparse" $Ti + ,:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Dense{Tv}}, Cint, Ptr{UInt8}), + D.p, true, common($Ti))) + finalizer(s, free!) + s end - function factor_to_sparse{Tv<:VTypes}(L::Factor{Tv,$Ti}) - Sparse(ccall((@cholmod_name "factor_to_sparse" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti))) + function factor_to_sparse!{Tv<:VTypes}(F::Factor{Tv,$Ti}) + ss = unsafe_load(F.p) + ss.xtype > PATTERN || throw(CHOLMODException("only numeric factors are supported")) + s = Sparse(ccall((@cholmod_name "factor_to_sparse" $Ti + ,:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + F.p, common($Ti))) + finalizer(s, free!) + s end - function copy_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) - Factor(ccall((@cholmod_name "copy_factor" $Ti - ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti))) + function change_factor{Tv<:VTypes}(::Type{Float64}, to_ll::Bool, to_super::Bool, to_packed::Bool, to_monotonic::Bool, F::Factor{Tv,$Ti}) + @isok ccall((@cholmod_name "change_factor" $Ti + ,:libcholmod), Cint, + (Cint, Cint, Cint, Cint, Cint, Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + REAL, to_ll, to_super, to_packed, to_monotonic, F.p, common($Ti)) + Factor{Float64,$Ti}(F.p) end - function change_factor{Tv<:VTypes}(to_xtype::Integer, to_ll::Bool, to_super::Bool, to_packed::Bool, to_monotonic::Bool, L::Factor{Tv,$Ti}) + function change_factor{Tv<:VTypes}(::Type{Complex{Float64}}, to_ll::Bool, to_super::Bool, to_packed::Bool, to_monotonic::Bool, F::Factor{Tv,$Ti}) @isok ccall((@cholmod_name "change_factor" $Ti ,:libcholmod), Cint, - (Cint, Cint, Cint, Cint, Cint, Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - to_xtype, to_ll, to_super, to_packed, to_monotonic, L, cm) - L + (Cint, Cint, Cint, Cint, Cint, Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + COMPLEX, to_ll, to_super, to_packed, to_monotonic, F.p, common($Ti)) + Factor{Complex{Float64},$Ti}(F.p) end function check_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}) bool(ccall((@cholmod_name "check_sparse" $Ti ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, common($Ti))) + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, common($Ti))) end - function check_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) + function check_factor{Tv<:VTypes}(F::Factor{Tv,$Ti}) bool(ccall((@cholmod_name "check_factor" $Ti ,:libcholmod), Cint, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti))) - end - - function analyze{Tv<:VTypes}(A::Sparse{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) - Factor(ccall((@cholmod_name "analyze" $Ti - ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, C)) - end - - function factorize{Tv<:VTypes}(A::Sparse{Tv,$Ti}, F::Factor{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) - @isok ccall((@cholmod_name "factorize" $Ti - ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &A.c, &F.c, C) - F - end - - function factorize_p{Tv<:VTypes}(A::Sparse{Tv,$Ti}, β::Tv, fset::Vector{$Ti}, F::Factor{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) - @isok ccall((@cholmod_name "factorize_p" $Ti - ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{Cdouble}, Ptr{$Ti}, Csize_t, - Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &A.c, &β, fset, length(fset), - &F.c, common($Ti)) - F - end - - ### cholmod_core.h ### - function pack_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) - @isok ccall((@cholmod_name "pack_factor" $Ti - ,:libcholmod), Cint, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti)) - L + (Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + F.p, common($Ti))) end function nnz{Tv<:VTypes}(A::Sparse{Tv,$Ti}) ccall((@cholmod_name "nnz" $Ti ,:libcholmod), Int, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, common($Ti)) + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, common($Ti)) end function speye{Tv<:VTypes}(m::Integer, n::Integer, ::Type{Tv}) - Sparse(ccall((@cholmod_name "speye" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + s = Sparse(ccall((@cholmod_name "speye" $Ti + , :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, (Csize_t, Csize_t, Cint, Ptr{UInt8}), m, n, xtyp(Tv), common($Ti))) + finalizer(s, free!) + s end function spzeros{Tv<:VTypes}(m::Integer, n::Integer, nzmax::Integer, ::Type{Tv}) - Sparse(ccall((@cholmod_name "spzeros" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + s = Sparse(ccall((@cholmod_name "spzeros" $Ti + , :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, (Csize_t, Csize_t, Csize_t, Cint, Ptr{UInt8}), m, n, nzmax, xtyp(Tv), common($Ti))) + finalizer(s, free!) + s + end + + function transpose{Tv<:VTypes}(A::Sparse{Tv,$Ti}, values::Integer) + s = Sparse(ccall((@cholmod_name "transpose" $Ti + ,:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + A.p, values, common($Ti))) + finalizer(s, free!) + s end - function copy_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) - Factor(ccall((@cholmod_name "copy_factor" $Ti - ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti))) + function copy_factor{Tv<:VTypes}(F::Factor{Tv,$Ti}) + f = Factor(ccall((@cholmod_name "copy_factor" $Ti + ,:libcholmod), Ptr{C_Factor{Tv,$Ti}}, + (Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + F.p, common($Ti))) + finalizer(f, free!) + f end function copy_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}) - Sparse(ccall((@cholmod_name "copy_sparse" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, common($Ti))) + s = Sparse(ccall((@cholmod_name "copy_sparse" $Ti + ,:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, common($Ti))) + finalizer(s, free!) + s end - function copy{Tv<:VTypes}(A::Sparse{Tv,$Ti}, stype::Integer, mode::Integer) - Sparse(ccall((@cholmod_name "copy" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), - &A.c, stype, mode, common($Ti))) + function copy{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, stype::Integer, mode::Integer) + s = Sparse(ccall((@cholmod_name "copy" $Ti + ,:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), + A.p, stype, mode, common($Ti))) + finalizer(s, free!) + s end ### cholmod_check.h ### @@ -674,41 +493,70 @@ for Ti in IndexTypes set_print_level(cm, 3) @isok ccall((@cholmod_name "print_sparse" $Ti ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), - &A.c, name, cm) + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), + A.p, name, cm) nothing end - function print_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}, name::ASCIIString) + function print_factor{Tv<:VTypes}(F::Factor{Tv,$Ti}, name::ASCIIString) cm = common($Ti) set_print_level(cm, 3) @isok ccall((@cholmod_name "print_factor" $Ti ,:libcholmod), Cint, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), - &L.c, name, cm) + (Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), + F.p, name, cm) nothing end ### cholmod_matrixops.h ### + function ssmult{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, stype::Integer, values::Bool, sorted::Bool) + lA = unsafe_load(A.p) + lB = unsafe_load(B.p) + lA.ncol == lB.nrow || throw(DimensionMismatch("inner matrix dimensions do not fit")) + s = Sparse(ccall((@cholmod_name "ssmult" $Ti + ,:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{C_Sparse{Tv,$Ti}}, Cint, Cint, Cint, Ptr{UInt8}), + A.p, B.p, stype, values, sorted, common($Ti))) + finalizer(s, free!) + s + end + function norm_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}, norm::Integer) norm == 0 || norm == 1 || throw(ArgumentError("norm argument must be either 0 or 1")) ccall((@cholmod_name "norm_sparse" $Ti , :libcholmod), Cdouble, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, norm, common($Ti)) + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + A.p, norm, common($Ti)) end - function horzcat{Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) - ccall((@cholmod_name "horzcat" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, &B.c, values, common($Ti)) + function horzcat{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) + s = Sparse(ccall((@cholmod_name "horzcat" $Ti + , :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + A.p, B.p, values, common($Ti))) + finalizer(s, free!) + s end - function scale{Tv<:VTypes}(S::Dense{Tv}, scale::Integer, A::Sparse{Tv,$Ti}) + function scale{Tv<:VRealTypes}(S::Dense{Tv}, scale::Integer, A::Sparse{Tv,$Ti}) + sS = unsafe_load(S.p) + sA = unsafe_load(A.p) + sS.ncol == 1 || sS.nrow == 1 || throw(DimensionMismatch("first argument must be a vector")) + if scale == SCALAR + sS.nrow == 1 || throw(DimensionMismatch("scaling argument must have length one")) + elseif scale == ROW + sS.nrow*sS.ncol == sA.nrow || throw(DimensionMismatch("scaling vector has length $(sS.nrow*sS.ncol), but matrix has $(sA.nrow) rows.")) + elseif scale == COL + sS.nrow*sS.ncol == sA.ncol || throw(DimensionMismatch("scaling vector has length $(sS.nrow*sS.ncol), but matrix has $(sA.ncol) columns")) + elseif scale == SYM + sA.nrow == sA.ncol || throw(DimensionMismatch("matrix must be square")) + sS.nrow*sS.ncol == sA.nrow || throw(DimensionMismatch("scaling vector has length $(sS.nrow*sS.ncol), but matrix has $(sA.ncol) columns and rows")) + end + + sA = unsafe_load(A.p) @isok ccall((@cholmod_name "scale" $Ti ,:libcholmod), Cint, - (Ptr{c_Dense{Tv}}, Cint, Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &S.c, scale, &A.c, common($Ti)) + (Ptr{C_Dense{Tv}}, Cint, Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + S.p, scale, A.p, common($Ti)) A end @@ -716,23 +564,25 @@ for Ti in IndexTypes m, n = size(A) nc = transpose ? m : n nr = transpose ? n : m - if nc != size(X,1) - throw(DimensionMismatch("Incompatible dimensions, $nc and $(size(X,1))")) + if nc != size(X, 1) + throw(DimensionMismatch("incompatible dimensions, $nc and $(size(X,1))")) end @isok ccall((@cholmod_name "sdmult" $Ti ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, - Ptr{c_Dense{Tv}}, Ptr{c_Dense{Tv}}, Ptr{UInt8}), - &A.c, transpose, &α, &β, - &X.c, &Y.c, common($Ti)) + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, + Ptr{C_Dense{Tv}}, Ptr{C_Dense{Tv}}, Ptr{UInt8}), + A.p, transpose, &α, &β, + X.p, Y.p, common($Ti)) Y end - function vertcat{Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) - ccall((@cholmod_name "vertcat" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, &B.c, values, common($Ti)) + function vertcat{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) + s = Sparse(ccall((@cholmod_name "vertcat" $Ti + , :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + A.p, B.p, values, common($Ti))) + finalizer(s, free!) + s end function symmetry{Tv<:VTypes}(A::Sparse{Tv,$Ti}, option::Integer) @@ -742,65 +592,192 @@ for Ti in IndexTypes nzdiag = Array($Ti, 1) rv = ccall((@cholmod_name "symmetry" $Ti , :libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{$Ti}, Ptr{$Ti}, + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{$Ti}, Ptr{$Ti}, Ptr{$Ti}, Ptr{$Ti}, Ptr{UInt8}), - &A.c, option, xmatched, pmatched, + A.p, option, xmatched, pmatched, nzoffdiag, nzdiag, common($Ti)) rv, xmatched[1], pmatched[1], nzoffdiag[1], nzdiag[1] end # cholmod_cholesky.h - function etree{Tv<:VTypes}(A::Sparse{Tv,$Ti}) - Parent = Array($Ti, size(A,2)) - @isok ccall((@cholmod_name "etree" $Ti + # For analyze, factorize and factorize_p, the Common argument must be supplied in order to control if the factorization is LLt or LDLt + function analyze{Tv<:VTypes}(A::Sparse{Tv,$Ti}, cmmn::Vector{UInt8}) + f = Factor(ccall((@cholmod_name "analyze" $Ti + ,:libcholmod), Ptr{C_Factor{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, cmmn)) + finalizer(f, free!) + f + end + function factorize{Tv<:VTypes}(A::Sparse{Tv,$Ti}, F::Factor{Tv,$Ti}, cmmn::Vector{UInt8}) + @isok ccall((@cholmod_name "factorize" $Ti + ,:libcholmod), Cint, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + A.p, F.p, cmmn) + F + end + function factorize_p{Tv<:VTypes}(A::Sparse{Tv,$Ti}, β::Tv, fset::Vector{$Ti}, F::Factor{Tv,$Ti}, cmmn::Vector{UInt8}) + @isok ccall((@cholmod_name "factorize_p" $Ti ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{$Ti}, Ptr{UInt8}), - &A.c, Parent, common($Ti)) - return Parent + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{Cdouble}, Ptr{$Ti}, Csize_t, + Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + A.p, &β, fset, length(fset), + F.p, cmmn) + F end - function solve{Tv<:VTypes}(sys::Integer, L::Factor{Tv,$Ti}, B::Dense{Tv}) - size(L,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(L,1)) rows, but RHS has $(size(B,1)) rows.")) - Dense(ccall((@cholmod_name "solve" $Ti - ,:libcholmod), Ptr{c_Dense{Tv}}, - (Cint, Ptr{c_Factor{Tv,$Ti}}, Ptr{c_Dense{Tv}}, Ptr{UInt8}), - sys, &L.c, &B.c, common($Ti))) + function solve{Tv<:VTypes}(sys::Integer, F::Factor{Tv,$Ti}, B::Dense{Tv}) + 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.")) + d = Dense(ccall((@cholmod_name "solve" $Ti + ,:libcholmod), Ptr{C_Dense{Tv}}, + (Cint, Ptr{C_Factor{Tv,$Ti}}, Ptr{C_Dense{Tv}}, Ptr{UInt8}), + sys, F.p, B.p, common($Ti))) + finalizer(d, free!) + d end - function spsolve{Tv<:VTypes}(sys::Integer, L::Factor{Tv,$Ti}, B::Sparse{Tv,$Ti}) - size(L,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(L,1)) rows, but RHS has $(size(B,1)) rows.")) - Sparse(ccall((@cholmod_name "spsolve" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Cint, Ptr{c_Factor{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - sys, &L.c, &B.c, common($Ti))) + function spsolve{Tv<:VTypes}(sys::Integer, F::Factor{Tv,$Ti}, B::Sparse{Tv,$Ti}) + size(F,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(L,1)) rows, but RHS has $(size(B,1)) rows.")) + s = Sparse(ccall((@cholmod_name "spsolve" $Ti + ,:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Cint, Ptr{C_Factor{Tv,$Ti}}, Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + sys, F.p, B.p, common($Ti))) + finalizer(s, free!) + s end end end # Autodetects the types function read_sparse(file::CFILE) - Sparse(ccall((:cholmod_read_sparse, :libcholmod), Ptr{c_SparseVoid}, + s = Sparse(ccall((:cholmod_read_sparse, :libcholmod), Ptr{C_SparseVoid}, (Ptr{Void}, Ptr{UInt8}), file.ptr, common(Int32))) + finalizer(s, free!) + s end ######################### # High level interfaces # ######################### -# Convertion +# Convertion/construction +function Dense{T<:VTypes}(A::VecOrMat{T}) + d = allocate_dense(size(A, 1), size(A, 2), stride(A, 2), T) + s = unsafe_load(d.p) + unsafe_copy!(s.x, pointer(A), length(A)) + d +end Dense(A::Sparse) = sparse_to_dense(A) -Sparse(A::Dense) = dense_to_sparse(A) -function Sparse(L::Factor) - if bool(L.c.is_ll) - return factor_to_sparse(L) +# This constructior assumes zero based colptr and rowval +function Sparse{Tv<:VTypes,Ti<:ITypes}(m::Integer, n::Integer, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}, stype) + + # check if columns are sorted + iss = true + for i = 2:length(colptr) + if !issorted(sub(rowval, colptr[i - 1] + 1:colptr[i])) + iss = false + break + end + end + + o = allocate_sparse(m, n, length(nzval), iss, true, stype, Tv, Ti) + s = unsafe_load(o.p) + + unsafe_copy!(s.p, pointer(colptr), length(colptr)) + unsafe_copy!(s.i, pointer(rowval), length(rowval)) + unsafe_copy!(s.x, pointer(nzval), length(nzval)) + + @isok check_sparse(o) + + return o + +end +function Sparse{Tv<:VTypes,Ti<:ITypes}(m::Integer, n::Integer, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) + o = Sparse(m, n, colptr, rowval, nzval, 0) + + # check if array is symmetric and change stype if it is + i = symmetry(o, 0)[1] # Check for Hermitianity + if (isreal(nzval) && (i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG)) || # real case + (i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG) # complex case + change_stype!(o, -1) + end + o +end + +function Sparse{Tv<:VTypes,Ti<:ITypes}(A::SparseMatrixCSC{Tv,Ti}, stype::Integer) + o = allocate_sparse(A.m, A.n, length(A.nzval), true, true, stype, Tv, Ti) + s = unsafe_load(o.p) + for i = 1:length(A.colptr) + unsafe_store!(s.p, A.colptr[i] - 1, i) + end + for i = 1:length(A.rowval) + unsafe_store!(s.i, A.rowval[i] - 1, i) + end + unsafe_copy!(s.x, pointer(A.nzval), length(A.nzval)) + + @isok check_sparse(o) + + return o +end +function Sparse{Tv<:VTypes,Ti<:ITypes}(A::SparseMatrixCSC{Tv,Ti}) + o = Sparse(A, 0) + # check if array is symmetric and change stype if it is + i = symmetry(o, 0)[1] # Check for Hermitianity + if (isreal(A) && (i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG)) || # real case + (i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG) # complex case + change_stype!(o, -1) + end + o +end + +Sparse{Ti<:ITypes}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) +Sparse{Ti<:ITypes}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) + +# Useful when reading in files, but not type stable +function Sparse(p::Ptr{C_SparseVoid}) + s = unsafe_load(p) + + # Check integer type + if s.itype == INT + Ti = Cint + elseif s.itype == INTLONG + free_sparse(p) + throw(CHOLMODException("the value of itype was $s.itype. This combination of integer types shouldn't happen. Please submit a bug report.")) + elseif s.itype == LONG + Ti = SuiteSparse_long + else + free_sparse(p) + throw(CHOLMODException("illegal value of itype: $s.itype")) + end + + # Check for double or single precision + if s.dtype == DOUBLE + Tv = Float64 + elseif s.dtype == SINGLE + # Tv = Float32 # this should be supported at some point + free_sparse(p) + throw(CHOLMODException("single precision not supported yet")) + else + free_sparse(p) + throw(CHOLMODException("illegal value of dtype: $s.dtype")) + end + + # Check for real or complex + if s.xtype == COMPLEX + Tv = Complex{Tv} + elseif s.xtype != REAL + free_sparse(p) + throw(CHOLMODException("illegal value of xtype: $s.xtype")) end - cm = common($Ti) - Lcll = copy_factor(L) - change_factor(L.c.xtype, true, L.c.is_super == 1, true, true, Lcll, L) - return factor_to_sparse(Lcll) + + return Sparse(convert(Ptr{C_Sparse{Tv,Ti}}, p)) + end + +Sparse(A::Dense) = dense_to_sparse(A, Cint) # default is Cint which is probably sufficient when converted from dense matrix +Sparse(L::Factor) = factor_to_sparse!(copy(L)) function Sparse(filename::ASCIIString) f = open(filename) A = read_sparse(CFILE(f)) @@ -808,7 +785,66 @@ function Sparse(filename::ASCIIString) A end -show(io::IO, L::Factor) = print_factor(L, "") +## convertion back to base Julia types +function convert{T<:VTypes}(::Type{Matrix}, D::Dense{T}) + s = unsafe_load(D.p) + a = Array(T, s.nrow, s.ncol) + if s.d == s.nrow + unsafe_copy!(pointer(a), s.x, s.d*s.ncol) + else + for j = 1:s.ncol + for i = 1:s.nrow + a[i,j] = unsafe_load(s.x, i + (j - 1)*s.d) + end + end + end + a +end + +function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, A::Sparse{Tv,Ti}) + s = unsafe_load(A.p) + B = SparseMatrixCSC(s.nrow, s.ncol, increment(pointer_to_array(s.p, (s.ncol + 1,), false)), increment(pointer_to_array(s.i, (s.nzmax,), false)), copy(pointer_to_array(s.x, (s.nzmax,), false))) + if s.stype != 0 + return isreal(A) ? Symmetric(B, s.stype < 0 ? :L : :U) : Hermitian(B, s.stype < 0 ? :L : :U) + end + return B +end +sparse{Tv,Ti}(A::Sparse{Tv,Ti}) = convert(SparseMatrixCSC{Tv,Ti}, A) +sparse(L::Factor) = sparse(Sparse(L)) +sparse(D::Dense) = sparse(Sparse(D)) +function full(A::Sparse) + B = convert(Matrix, Dense(A)) + s = unsafe_load(A.p) + if s.stype != 0 + return isreal(A) ? Symmetric(B, s.stype < 0 ? :L : :U) : Hermitian(B, s.stype < 0 ? :L : :U) + end + return B +end +full(A::Dense) = convert(Matrix, A) + +# Calculate the offset into the stype field of the cholmod_sparse_struct and change the value +function change_stype!(A::Sparse, i::Integer) + offset = fieldoffsets(C_Sparse)[names(C_Sparse) .== :stype][1] + unsafe_store!(convert(Ptr{Cint}, A.p), i, div(offset, 4) + 1) + nothing +end + +free!(A::Dense) = free_dense(A.p) +free!(A::Sparse) = free_sparse(A.p) +free!(F::Factor) = free_factor(F.p) + +eltype{T<:VTypes}(A::Dense{T}) = T +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)) + @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) @@ -818,59 +854,55 @@ copy(A::Dense) = copy_dense(A) copy(A::Sparse) = copy_sparse(A) copy(A::Factor) = copy_factor(A) -size(B::Dense) = size(B.mat) -size(B::Dense, d::Integer) = size(B.mat, d) -size(A::Sparse) = (int(A.c.m), int(A.c.n)) -size(A::Sparse, d::Integer) = d > 0 ? (d == 1 ? int(A.c.m) : (d == 2 ? int(A.c.n) : 1)) : BoundsError() -size(L::Factor) = (n = int(L.c.n); (n,n)) -size(L::Factor, d::Integer) = d < 1 ? throw(BoundsError("Dimension $d out of range")) : (d <= 2 ? int(L.c.n) : 1) +function size(A::Union(Dense,Sparse)) + s = unsafe_load(A.p) + return (int(s.nrow), int(s.ncol)) +end +function size(F::Factor, i::Integer) + i < 1 && error("arraysize: dimension out of range") + s = unsafe_load(F.p) + if i <= 2 + return int(s.n) + end + return 1 +end + +function getindex(A::Dense, i::Integer) + s = unsafe_load(A.p) + 0 < i <= s.nrow*s.ncol || throw(BoundsError()) + unsafe_load(s.x, i) +end +function getindex(A::Dense, i::Integer, j::Integer) + s = unsafe_load(A.p) + 0 < i <= s.nrow || throw(BoundsError()) + 0 < j <= s.ncol || throw(BoundsError()) + unsafe_load(s.x, i + (j - 1)*s.d) +end -getindex(A::Dense, i::Integer) = A.mat[i] -getindex(A::Dense, i::Integer, j::Integer) = A.mat[i, j] getindex(A::Sparse, i::Integer) = getindex(A, ind2sub(size(A),i)...) function getindex{T}(A::Sparse{T}, i0::Integer, i1::Integer) - !(1 <= i0 <= A.c.m && 1 <= i1 <= A.c.n) && throw(BoundsError()) - A.c.stype < 0 && i0 < i1 && return conj(A[i1,i0]) - A.c.stype > 0 && i0 > i1 && return conj(A[i1,i0]) + s = unsafe_load(A.p) + !(1 <= i0 <= s.nrow && 1 <= i1 <= s.ncol) && throw(BoundsError()) + s.stype < 0 && i0 < i1 && return conj(A[i1,i0]) + s.stype > 0 && i0 > i1 && return conj(A[i1,i0]) - r1 = int(A.colptr0[i1] + 1) - r2 = int(A.colptr0[i1 + 1]) + r1 = int(unsafe_load(s.p, i1) + 1) + r2 = int(unsafe_load(s.p, i1 + 1)) (r1 > r2) && return zero(T) - r1 = int(searchsortedfirst(A.rowval0, i0 - 1, r1, r2, Base.Order.Forward)) - ((r1 > r2) || (A.rowval0[r1] + 1 != i0)) ? zero(T) : A.nzval[r1] -end - -# Convert invalidates the input matrix -function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, A::Sparse{Tv,Ti}) - B = SparseMatrixCSC(A.c.m, A.c.n, increment!(A.colptr0), increment!(A.rowval0), A.nzval) - if A.c.stype != 0 - return isreal(A) ? Symmetric(B, A.c.stype < 0 ? :L : :U) : Hermitian(B, A.c.stype < 0 ? :L : :U) - end - return B -end -sparse!{Tv,Ti}(A::Sparse{Tv,Ti}) = convert(SparseMatrixCSC{Tv,Ti}, A) -sparse(A::Sparse) = sparse!(copy(A)) -sparse(L::Factor) = sparse!(Sparse(L)) -sparse(D::Dense) = sparse!(Sparse(D)) -function full(A::Sparse) - B = Dense(A).mat - if A.c.stype != 0 - return isreal(A) ? Symmetric(B, A.c.stype < 0 ? :L : :U) : Hermitian(B, A.c.stype < 0 ? :L : :U) - end - return B + r1 = int(searchsortedfirst(pointer_to_array(s.i, (s.nzmax,), false), i0 - 1, r1, r2, Base.Order.Forward)) + ((r1 > r2) || (unsafe_load(s.i, r1) + 1 != i0)) ? zero(T) : unsafe_load(s.x, r1) end -full(A::Dense) = A.mat ## Multiplication - +(*)(A::Sparse, B::Sparse) = ssmult(A, B, 0, true, true) (*)(A::Sparse, B::Dense) = sdmult(A, false, 1., 0., B, zeros(size(A, 1), size(B, 2))) (*)(A::Sparse, B::VecOrMat) = (*)(A, Dense(B)) function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti}) - cm = common($Ti) + cm = common(Ti) if !is(A,B) - aa1 = transpose(B, 1) + aa1 = transpose(B, 2) ## result of ssmult will have stype==0, contain numerical values and be sorted return ssmult(A, aa1, 0, true, true) end @@ -878,61 +910,63 @@ function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti}) ## The A*A' case is handled by cholmod_aat. This routine requires ## A->stype == 0 (storage of upper and lower parts). If neccesary ## the matrix A is first converted to stype == 0 - if A.c.stype != 0 + s = unsafe_load(A.p) + if s.stype != 0 aa1 = copy(A, 0, 1) - return aat(aa1, Ti[], 1) + return aat(aa1, Ti[0:s.ncol-1], 1) else - return aat(A, Ti[], 1) + return aat(A, Ti[0:s.ncol-1], 1) end end -function Ac_mul_B{Tv<:VRealTypes}(A::Sparse{Tv}, B::Sparse{Tv}) - aa1 = transpose(A, 1) +function Ac_mul_B(A::Sparse, B::Sparse) + aa1 = transpose(A, 2) if is(A,B) return A_mul_Bc(aa1, aa1) end ## result of ssmult will have stype==0, contain numerical values and be sorted - return ssmult(aa1, B, 0, true, true, cm) + return ssmult(aa1, B, 0, true, true) end -A_mul_Bt{Tv<:VRealTypes}(A::Sparse{Tv}, B::Sparse{Tv}) = A_mul_Bc(A, B) # in the unlikely event of writing A*B.' instead of A*B' +A_mul_Bt(A::Sparse, B::Sparse) = A_mul_Bc(A, B) # in the unlikely event of writing A*B.' instead of A*B' -At_mul_B{Tv<:VRealTypes}(A::Sparse{Tv}, B::Sparse{Tv}) = Ac_mul_B(A,B) # in the unlikely event of writing A.'*B instead of A'*B +At_mul_B(A::Sparse, B::Sparse) = Ac_mul_B(A,B) # in the unlikely event of writing A.'*B instead of A'*B Ac_mul_B(A::Sparse, B::Dense) = sdmult(A, true, 1., 0., B, zeros(size(A, 2), size(B, 2))) Ac_mul_B(A::Sparse, B::VecOrMat) = Ac_mul_B(A, Dense(B)) ## Factorization methods - -analyze(A::SparseMatrixCSC) = analyze(Sparse(A)) - function cholfact(A::Sparse) + sA = unsafe_load(A.p) + sA.stype == 0 && throw(ArgumentError("sparse matrix is not symmetric/Hermitian")) cm = common(indtype(A)) ## may need to change final_asis as well as final_ll cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt F = analyze(A, cm) factorize(A, F, cm) - F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) + s = unsafe_load(F.p) + s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor)) return F end -function ldltfact(A::Sparse) - cm = common(indtype(A)) +function cholfact{Tv<:VTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, β::Tv) + cm = common(Ti) ## may need to change final_asis as well as final_ll - cm[common_final_ll] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt + cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt F = analyze(A, cm) - factorize(A, F, cm) + factorize_p(A, β, Ti[], F, cm) + s = unsafe_load(F.p) + s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor)) return F end -function cholfact{Tv<:VTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, β::Tv) - cm = common(Ti) +function ldltfact(A::Sparse) + cm = common(indtype(A)) ## may need to change final_asis as well as final_ll - cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + cm[common_final_ll] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt F = analyze(A, cm) - factorize_p(A, β, Ti[], F, cm) - F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) + factorize(A, F, cm) return F end @@ -942,164 +976,131 @@ function ldltfact{Tv<:VTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, β::Tv) cm[common_final_ll] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt F = analyze(A, cm) factorize_p(A, β, Ti[], F, cm) - F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) + s = unsafe_load(F.p) + s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor)) return F end -cholfact!(L::Factor, A::Sparse) = factorize(A, L) -cholfact!{Tv<:VTypes,Ti<:ITypes}(L::Factor{Tv,Ti}, A::Sparse{Tv,Ti}, β::Tv) = factorize_p(A, β, Ti[], L) - cholfact{T<:VTypes}(A::Sparse{T}, β::Number) = cholfact(A, convert(T, β)) cholfact(A::SparseMatrixCSC) = cholfact(Sparse(A)) +cholfact(A::SparseMatrixCSC, β::Number) = cholfact(Sparse(A), β) cholfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = cholfact(Sparse(A)) +cholfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}, β::Number) = cholfact(Sparse(A), β) cholfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = cholfact(Sparse(A)) -cholfact!{T<:VTypes}(L::Factor{T}, A::Sparse{T}, β::Number) = cholfact!(L, A, convert(T, β)) -cholfact!{T<:VTypes}(L::Factor{T}, A::SparseMatrixCSC{T}, β::Number) = cholfact!(L, Sparse(A), convert(T, β)) -cholfact!{T<:VTypes}(L::Factor{T}, A::SparseMatrixCSC{T}) = cholfact!(L, Sparse(A)) +cholfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}, β::Number) = cholfact(Sparse(A), β) +ldltfact{T<:VTypes}(A::Sparse{T}, β::Number) = ldltfact(A, convert(T, β)) ldltfact(A::SparseMatrixCSC) = ldltfact(Sparse(A)) +ldltfact(A::SparseMatrixCSC, β::Number) = ldltfact(Sparse(A), β) ldltfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = ldltfact(Sparse(A)) +ldltfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}, β::Number) = ldltfact(Sparse(A), β) ldltfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = ldltfact(Sparse(A)) +ldltfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}, β::Number) = ldltfact(Sparse(A), β) + +function update!{Tv<:VTypes,Ti<:ITypes}(F::Factor{Tv,Ti}, A::Sparse{Tv,Ti}) + cm = common(Ti) + s = unsafe_load(F.p) + if bool(s.is_ll) + cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + end + factorize(A, F, cm) +end +function update!{Tv<:VTypes,Ti<:ITypes}(F::Factor{Tv,Ti}, A::Sparse{Tv,Ti}, β::Tv) + cm = common(Ti) + s = unsafe_load(F.p) + if bool(s.is_ll) + cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + end + factorize_p(A, β, Ti[0:size(F, 1) - 1], F, cm) +end +update!{T<:VTypes}(F::Factor{T}, A::SparseMatrixCSC{T}) = update!(F, Sparse(A)) +update!{T<:VTypes}(F::Factor{T}, A::SparseMatrixCSC{T}, β::Number) = update!(F, Sparse(A), convert(T, β)) ## Solvers (\)(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B) -(\)(L::Factor, b::Vector) = reshape(solve(CHOLMOD_A, L, Dense(b)).mat, length(b)) -(\)(L::Factor, B::Matrix) = solve(CHOLMOD_A, L, Dense(B)).mat +(\)(L::Factor, b::Vector) = reshape(convert(Matrix, solve(CHOLMOD_A, L, Dense(b))), length(b)) +(\)(L::Factor, B::Matrix) = convert(Matrix, solve(CHOLMOD_A, L, Dense(B))) (\)(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B) # When right hand side is sparse, we have to ensure that the rhs is not marked as symmetric. -(\)(L::Factor, B::SparseMatrixCSC) = sparse!(spsolve(CHOLMOD_A, L, Sparse(B, 0))) +(\)(L::Factor, B::SparseMatrixCSC) = sparse(spsolve(CHOLMOD_A, L, Sparse(B, 0))) Ac_ldiv_B(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B) -Ac_ldiv_B(L::Factor, B::VecOrMat) = solve(CHOLMOD_A, L, Dense(B)).mat +Ac_ldiv_B(L::Factor, B::VecOrMat) = convert(Matrix, solve(CHOLMOD_A, L, Dense(B))) Ac_ldiv_B(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B) -Ac_ldiv_B(L::Factor, B::SparseMatrixCSC) = sparse!(CHOLMOD_A, spsolve(L, Sparse(B))) - -speye(m::Integer, n::Integer) = speye(m, n, Float64) # default element type is Float64 -speye(n::Integer) = speye(n, n, Float64) # default shape is square - -spzeros(m::Integer, n::Integer, nzmax::Integer) = spzeros(m, n, nzmax, 1.) - -function norm(A::Sparse, p = 1) - if p == 1 - return norm_sparse(A, 1) - elseif p == Inf - return norm_sparse(A, 0) - else - throw(ArgumentError("only 1 and inf norms are supported")) - end -end - -## Julia methods (non-CHOLMOD) -function findnz{Tv,Ti}(A::Sparse{Tv,Ti}) - jj = similar(A.rowval0) # expand A.colptr0 to a vector of indices - for j in 1:A.c.n, k in (A.colptr0[j]+1):A.colptr0[j+1] - jj[k] = j - end - - ind = similar(A.rowval0) - ipos = 1 - count = 0 - for k in 1:length(A.nzval) - if A.nzval[k] != 0 - ind[ipos] = k - ipos += 1 - count += 1 - else - println("warning: sparse matrix contains explicitly stored zeros") - end - end - ind = ind[1:count] # ind is the indices of nonzeros in A.nzval - (increment!(A.rowval0[ind]), jj[ind], A.nzval[ind]) -end - -findnz(L::Factor) = findnz(Sparse(L)) - -function diag{Tv}(A::Sparse{Tv}) - minmn = minimum(size(A)) - res = zeros(Tv,minmn) - cp0 = A.colptr0 - rv0 = A.rowval0 - anz = A.nzval - for j in 1:minmn, k in (cp0[j]+1):cp0[j+1] - if rv0[k] == j-1 - res[j] += anz[k] - end - end - res -end - -function diag{Tv}(L::Factor{Tv}) - res = zeros(Tv,L.c.n) - xv = L.x - if bool(L.c.is_super) - nec = decrement!(diff(L.super)) # number of excess columns per supernode - dstride = increment!(diff(L.pi)) # stride of diagonal elements (nrow + 1) - px = L.px +Ac_ldiv_B(L::Factor, B::SparseMatrixCSC) = Ac_ldiv_B(L, Sparse(B)) + +## Other convenience methods +function diag{Tv}(F::Factor{Tv}) + f = unsafe_load(F.p) + res = Base.zeros(Tv, int(f.n)) + xv = f.x + if bool(f.is_super) + px = f.px pos = 1 - for i in 1:length(nec) - base = px[i] + 1 - res[pos] = xv[base] + for i in 1:f.nsuper + base = unsafe_load(px, i) + 1 + res[pos] = unsafe_load(xv, base) pos += 1 - for j in 1:nec[i] - res[pos] = xv[base + j*dstride[i]] + for j in 1:unsafe_load(f.super, i + 1) - unsafe_load(f.super, i) - 1 + res[pos] = unsafe_load(xv, base + j*(unsafe_load(f.pi, i + 1) - unsafe_load(f.pi, i) + 1)) pos += 1 end end else - c0 = L.p - r0 = L.i - xv = L.x - for j in 1:L.c.n - jj = c0[j]+1 - assert(r0[jj] == j-1) - res[j] = xv[jj] + c0 = f.p + r0 = f.i + xv = f.x + for j in 1:f.n + jj = unsafe_load(c0, j) + 1 + assert(unsafe_load(r0, jj) == j - 1) + res[j] = unsafe_load(xv, jj) end end res end -function logdet{Tv,Ti}(L::Factor{Tv,Ti},sub) +function logdet{Tv<:VTypes,Ti<:ITypes}(F::Factor{Tv,Ti}) + f = unsafe_load(F.p) res = zero(Tv) - for d in diag(L)[sub] res += log(abs(d)) end - bool(L.c.is_ll) ? 2res : res + for d in diag(F) res += log(abs(d)) end + bool(f.is_ll) ? 2res : res end -logdet(L::Factor) = logdet(L, 1:L.c.n) + det(L::Factor) = exp(logdet(L)) -isposdef!{Tv<:VTypes,Ti}(A::SparseMatrixCSC{Tv,Ti}) = ishermitian(A) && cholfact(A).c.minor == size(A,1) +function isposdef{Tv<:VTypes,Ti}(A::SparseMatrixCSC{Tv,Ti}) + if !ishermitian(A) + return false + end + f = cholfact(A) + s = unsafe_load(f.p) + return s.minor == size(A,1) +end function issym(A::Sparse) - if isreal(A) - if A.c.stype != 0 - return true - else - i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 - return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG - end - else - if A.c.stype != 0 - return false - else - i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 - return i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG - end + s = unsafe_load(A.p) + if s.stype != 0 + return isreal(A) end + i = symmetry(A, ifelse(version >= v"3.0.5", 0, 1))[1] # 0 is faster, but had a bug before 3.0.5 + return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG end function ishermitian(A::Sparse) + s = unsafe_load(A.p) if isreal(A) - if A.c.stype != 0 + if s.stype != 0 return true else - i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 + i = symmetry(A, ifelse(version >= v"3.0.5", 0, 1))[1] # 0 is faster, but had a bug before 3.0.5 return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG end else - if A.c.stype != 0 + if s.stype != 0 return true else - i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 + i = symmetry(A, ifelse(version >= v"3.0.5", 0, 1))[1] # 0 is faster, but had a bug before 3.0.5 return i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG end end diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index f4ff95d67b9d2..0641090dc821b 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -33,10 +33,11 @@ using Base.SparseMatrix.CHOLMOD ## residual 1.3e-19 (|Ax-b|/(|A||x|+|b|)) after iterative refinement ## rcond 9.5e-06 -A = Sparse(int32([0,1,2,3,6,9,12,15,18,20,25,30,34,36,39,43,47,52,58,62,67,71,77,84,90,93,95, +A = CHOLMOD.Sparse(48, 48, + int32([0,1,2,3,6,9,12,15,18,20,25,30,34,36,39,43,47,52,58,62,67,71,77,84,90,93,95, 98,103,106,110,115,119,123,130,136,142,146,150,155,161,167,174,182,189,197, 207,215,224]), # zero-based column pointers - int32([0,1,2,1,2,3,0,2,4,0,1,5,0,4,6,1,3,7,2,8,1,3,7,8,9,0,4,6,8,10,5,6,7,11,6,12, + int32([0,1,2,1,2,3,0,2,4,0,1,5,0,4,6,1,3,7,2,8,1,3,7,8,9,0,4,6,8,10,5,6,7,11,6,12, 7,11,13,8,10,13,14,9,13,14,15,8,10,12,14,16,7,11,12,13,16,17,0,12,16,18,1, 5,13,15,19,2,4,14,20,3,13,15,19,20,21,2,4,12,16,18,20,22,1,5,17,18,19,23,0, 5,24,1,25,2,3,26,2,3,25,26,27,4,24,28,0,5,24,29,6,11,24,28,30,7,25,27,31,8, @@ -85,27 +86,28 @@ A = Sparse(int32([0,1,2,3,6,9,12,15,18,20,25,30,34,36,39,43,47,52,58,62,67,71,77 5.0e8,-4.79857058902e6,134990.2747,2.47238730198e9,-1.14928529451e6, 2.29724661236e8,-5.57173510779e7,-833333.333333,-1.25e6,2.5e8,2.39928529451e6, 9.61679848804e8,275828.470683,-5.57173510779e7,1.09411960038e7,2.08333333333e6, - 1.0e8,-2.5e6,140838.195984,-1.09779731332e8,5.31278103775e8], 48, 48, 1) -@test_approx_eq norm(A, Inf) 3.570948074697437e9 -@test_approx_eq norm(A) 3.570948074697437e9 + 1.0e8,-2.5e6,140838.195984,-1.09779731332e8,5.31278103775e8], 1) +@test_approx_eq CHOLMOD.norm_sparse(A, 0) 3.570948074697437e9 +@test_approx_eq CHOLMOD.norm_sparse(A, 1) 3.570948074697437e9 @test isvalid(A) B = A * ones(size(A,2)) chma = ldltfact(A) # LDL' form @test isvalid(chma) x = chma\B -@test_approx_eq x.mat ones(size(x)) +@test_approx_eq x ones(size(x)) chma = cholfact(A) # LL' form @test isvalid(chma) x = chma\B -@test_approx_eq x.mat ones(size(x)) +@test_approx_eq x ones(size(x)) #lp_afiro example -afiro = Sparse(int32([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,25,27,29,33,37, +afiro = CHOLMOD.Sparse(27, 51, + int32([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,25,27,29,33,37, 41,45,47,49,51,53,55,57,59,63,65,67,69,71,75,79,83,87,89,91,93,95,97, 99,101,102]), - int32([2,3,6,7,8,9,12,13,16,17,18,19,20,21,22,23,24,25,26,0,1,2,23,0,3,0,21, + int32([2,3,6,7,8,9,12,13,16,17,18,19,20,21,22,23,24,25,26,0,1,2,23,0,3,0,21, 1,25,4,5,6,24,4,5,7,24,4,5,8,24,4,5,9,24,6,20,7,20,8,20,9,20,3,4,4,22, 5,26,10,11,12,21,10,13,10,23,10,20,11,25,14,15,16,22,14,15,17,22,14, 15,18,22,14,15,19,22,16,20,17,20,18,20,19,20,13,15,15,24,14,26,15]), @@ -115,13 +117,15 @@ afiro = Sparse(int32([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,25,27 2.386,-1.0,2.408,-1.0,2.429,1.4,1.0,1.0,-1.0,1.0,1.0,-1.0,-0.43,1.0,0.109, 1.0,-1.0,1.0,-1.0,1.0,-1.0,1.0,1.0,-0.43,1.0,1.0,0.109,-0.43,1.0,1.0,0.108, -0.39,1.0,1.0,0.108,-0.37,1.0,1.0,0.107,-1.0,2.191,-1.0,2.219,-1.0,2.249, - -1.0,2.279,1.4,-1.0,1.0,-1.0,1.0,1.0,1.0], 27, 51, 0) -chmaf = cholfact(afiro) + -1.0,2.279,1.4,-1.0,1.0,-1.0,1.0,1.0,1.0], 0) +afiro2 = CHOLMOD.aat(afiro, Int32[0:50], int32(1)) +CHOLMOD.change_stype!(afiro2, -1) +chmaf = cholfact(afiro2) y = afiro'*ones(size(afiro,1)) sol = chmaf\(afiro*y) # least squares solution @test isvalid(sol) pred = afiro'*sol -@test norm(afiro * (y.mat - pred.mat)) < 1e-8 +@test norm(afiro * (convert(Matrix, y) - convert(Matrix, pred))) < 1e-8 let # Issue 9160 @@ -130,7 +134,7 @@ let # Issue 9160 A = sprand(10,10,0.1) A = convert(SparseMatrixCSC{elty,Ti},A) - cmA = CHOLMOD.parse(A) + cmA = CHOLMOD.Sparse(A) B = sprand(10,10,0.1) B = convert(SparseMatrixCSC{elty,Ti},B) @@ -156,6 +160,7 @@ let # Issue 9160 end end + # Issue #9915 @test speye(2)\speye(2) == eye(2) @@ -172,11 +177,222 @@ ACSC = sprandn(10, 10, 0.3) + I @test ishermitian(Sparse(Hermitian(complex(ACSC), :U))) # test Sparse constructor for c_Sparse{Tv,Ti} input( and Sparse*Sparse) -B = Sparse(SparseMatrixCSC{Float64,Int32}(sprandn(48, 48, 0.1))) # A has Int32 indeces +B = CHOLMOD.Sparse(SparseMatrixCSC{Float64,Int32}(sprandn(48, 48, 0.1))) # A has Int32 indeces @test_approx_eq A*B sparse(A)*sparse(B) # test Sparse constructor for c_SparseVoid (and read_sparse) writedlm("tmp.mtx", ["%%MatrixMarket matrix coordinate real symmetric","3 3 4","1 1 1","2 2 1","3 2 0.5","3 3 1"]) -@test full(Sparse("tmp.mtx")) == [1 0 0;0 1 0.5;0 0.5 1] +@test full(CHOLMOD.Sparse("tmp.mtx")) == [1 0 0;0 1 0.5;0 0.5 1] +run(`rm tmp.mtx`) +writedlm("tmp.mtx", ["%%MatrixMarket matrix coordinate complex Hermitian","3 3 4","1 1 1.0 0.0","2 2 1.0 0.0","3 2 0.5 0.5","3 3 1.0 0.0"]) +@test full(CHOLMOD.Sparse("tmp.mtx")) == [1 0 0;0 1 0.5-0.5im;0 0.5+0.5im 1] run(`rm tmp.mtx`) +# test that Sparse(Ptr) constructor throws the right places +## The struct pointer must be constructed by the library constructor and then modified afterwards to checks that the method throws +### illegal dtype (for now but should be supprted at some point) +p = ccall((:cholmod_allocate_sparse, :libcholmod), Ptr{CHOLMOD.C_SparseVoid}, + (Csize_t, Csize_t, Csize_t, Cint, Cint, Cint, Cint, Ptr{Void}), + 1, 1, 1, true, true, 0, CHOLMOD.REAL, CHOLMOD.common(Cint)) +pa = pointer_to_array(convert(Ptr{Uint32}, p), 8*2 + 6, false) +pa[8*2 + 4] = CHOLMOD.SINGLE +@test_throws CHOLMOD.CHOLMODException CHOLMOD.Sparse(p) + +### illegal dtype +p = ccall((:cholmod_allocate_sparse, :libcholmod), Ptr{CHOLMOD.C_SparseVoid}, + (Csize_t, Csize_t, Csize_t, Cint, Cint, Cint, Cint, Ptr{Void}), + 1, 1, 1, true, true, 0, CHOLMOD.REAL, CHOLMOD.common(Cint)) +pa = pointer_to_array(convert(Ptr{Uint32}, p), 8*2 + 6, false) +pa[8*2 + 4] = 5 +@test_throws CHOLMOD.CHOLMODException CHOLMOD.Sparse(p) + +### illegal xtype +p = ccall((:cholmod_allocate_sparse, :libcholmod), Ptr{CHOLMOD.C_SparseVoid}, + (Csize_t, Csize_t, Csize_t, Cint, Cint, Cint, Cint, Ptr{Void}), + 1, 1, 1, true, true, 0, CHOLMOD.REAL, CHOLMOD.common(Cint)) +pa = pointer_to_array(convert(Ptr{Uint32}, p), 8*2 + 6, false) +pa[8*2 + 3] = 3 +@test_throws CHOLMOD.CHOLMODException CHOLMOD.Sparse(p) + +### illegal itype +p = ccall((:cholmod_allocate_sparse, :libcholmod), Ptr{CHOLMOD.C_SparseVoid}, + (Csize_t, Csize_t, Csize_t, Cint, Cint, Cint, Cint, Ptr{Void}), + 1, 1, 1, true, true, 0, CHOLMOD.REAL, CHOLMOD.common(Cint)) +pa = pointer_to_array(convert(Ptr{Uint32}, p), 8*2 + 6, false) +pa[8*2 + 2] = CHOLMOD.INTLONG +@test_throws CHOLMOD.CHOLMODException CHOLMOD.Sparse(p) + +### illegal itype +p = ccall((:cholmod_allocate_sparse, :libcholmod), Ptr{CHOLMOD.C_SparseVoid}, + (Csize_t, Csize_t, Csize_t, Cint, Cint, Cint, Cint, Ptr{Void}), + 1, 1, 1, true, true, 0, CHOLMOD.REAL, CHOLMOD.common(Cint)) +pa = pointer_to_array(convert(Ptr{Uint32}, p), 8*2 + 6, false) +pa[8*2 + 2] = 5 +@test_throws CHOLMOD.CHOLMODException CHOLMOD.Sparse(p) + +# test that Sparse(Ptr) works for SuiteSparse_long +p = ccall((:cholmod_l_allocate_sparse, :libcholmod), Ptr{CHOLMOD.C_SparseVoid}, + (Csize_t, Csize_t, Csize_t, Cint, Cint, Cint, Cint, Ptr{Void}), + 1, 1, 1, true, true, 0, CHOLMOD.REAL, CHOLMOD.common(CHOLMOD.SuiteSparse_long)) +@test isa(CHOLMOD.Sparse(p), CHOLMOD.Sparse{Float64,CHOLMOD.SuiteSparse_long}) + +# Test Dense wrappers (only Float64 supported a present) + +## High level interfact +for elty in (Float64, Complex{Float64}) + if elty == Float64 + A = randn(5, 5) + b = randn(5) + else + A = complex(randn(5, 5), randn(5, 5)) + b = complex(randn(5), randn(5)) + end + ADense = CHOLMOD.Dense(A) + bDense = CHOLMOD.Dense(b) + + @test copy(ADense) == ADense + @test_approx_eq CHOLMOD.norm_dense(ADense, 1) norm(A, 1) + @test_approx_eq CHOLMOD.norm_dense(ADense, 0) norm(A, Inf) + @test_throws ArgumentError CHOLMOD.norm_dense(ADense, 2) + @test_throws ArgumentError CHOLMOD.norm_dense(ADense, 3) + + @test_approx_eq CHOLMOD.norm_dense(bDense, 2) norm(b) + @test CHOLMOD.check_dense(bDense) + + AA = CHOLMOD.eye(3) + unsafe_store!(convert(Ptr{Csize_t}, AA.p), 2, 1) # change size, but not stride, of Dense + @test convert(Matrix, AA) == eye(2, 3) +end + +## Low level interface +@test isa(CHOLMOD.zeros(3, 3, Float64), CHOLMOD.Dense{Float64}) +@test isa(CHOLMOD.zeros(3, 3), CHOLMOD.Dense{Float64}) +@test isa(CHOLMOD.zeros(3, 3, Float64), CHOLMOD.Dense{Float64}) +@test isa(CHOLMOD.ones(3, 3), CHOLMOD.Dense{Float64}) +@test isa(CHOLMOD.eye(3, 4, Float64), CHOLMOD.Dense{Float64}) +@test isa(CHOLMOD.eye(3, 4), CHOLMOD.Dense{Float64}) +@test isa(CHOLMOD.eye(3), CHOLMOD.Dense{Float64}) +@test isa(CHOLMOD.copy_dense(CHOLMOD.eye(3)), CHOLMOD.Dense{Float64}) + +# Test Sparse and Factor +## test free_sparse +p = ccall((:cholmod_allocate_sparse, :libcholmod), Ptr{CHOLMOD.C_Sparse{Float64,Cint}}, + (Csize_t, Csize_t, Csize_t, Cint, Cint, Cint, Cint, Ptr{Void}), + 1, 1, 1, true, true, 0, CHOLMOD.REAL, CHOLMOD.common(Cint)) +@test CHOLMOD.free_sparse(p) + +for elty in (Float64, Complex{Float64}) + A1 = sparse([1:5, 1], [1:5, 2], elty == Float64 ? randn(6) : complex(randn(6), randn(6))) + A2 = sparse([1:5, 1], [1:5, 2], elty == Float64 ? randn(6) : complex(randn(6), randn(6))) + A1pd = A1'A1 + A1Sparse = CHOLMOD.Sparse(A1) + A2Sparse = CHOLMOD.Sparse(A2) + A1pdSparse = CHOLMOD.Sparse( + A1pd.m, + A1pd.n, + Base.SparseMatrix.decrement(A1pd.colptr), + Base.SparseMatrix.decrement(A1pd.rowval), + A1pd.nzval) + + ## High level interface + @test isa(CHOLMOD.Sparse(3, 3, [0,1,3,4], [0,2,1,2], ones(4)), CHOLMOD.Sparse) # Sparse doesn't require columns to be sorted + @test full(A1Sparse) == A1 + @test copy(A1Sparse) == A1Sparse + @test size(A1Sparse, 3) == 1 + if elty <: Real # multiplcation only defined for real matrices in CHOLMOD + @test_approx_eq A1Sparse*A2Sparse A1*A2 + @test_approx_eq A1Sparse'A2Sparse A1'A2 + @test_approx_eq A1Sparse*A2Sparse' A1*A2' + + @test_approx_eq A1Sparse*A1Sparse A1*A1 + @test_approx_eq A1Sparse'A1Sparse A1'A1 + @test_approx_eq A1Sparse*A1Sparse' A1*A1' + + @test_approx_eq A1pdSparse*A1pdSparse A1pd*A1pd + @test_approx_eq A1pdSparse'A1pdSparse A1pd'A1pd + @test_approx_eq A1pdSparse*A1pdSparse' A1pd*A1pd' + end + + # Factor + F = cholfact(A1pd) + tmp = IOBuffer() + show(tmp, F) + @test tmp.size > 0 + @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) + @test_approx_eq F\CHOLMOD.Sparse(sparse(ones(elty, 5))) A1pd\ones(5) + @test_approx_eq F'\ones(elty, 5) full(A1pd)'\ones(5) + @test_approx_eq F'\sparse(ones(elty, 5)) full(A1pd)'\ones(5) + @test_approx_eq det(F) det(full(A1pd)) + @test isposdef(A1pd) + @test !isposdef(A1) + if elty <: Real + @test CHOLMOD.issym(Sparse(A1pd, 0)) + @test CHOLMOD.Sparse(cholfact(Symmetric(A1pd, :L))) == CHOLMOD.Sparse(cholfact(A1pd)) + @test CHOLMOD.Sparse(cholfact(Symmetric(A1pd, :L), 2)) == CHOLMOD.Sparse(cholfact(A1pd, 2)) + @test CHOLMOD.Sparse(ldltfact(Symmetric(A1pd, :L))) == CHOLMOD.Sparse(ldltfact(A1pd)) + @test CHOLMOD.Sparse(ldltfact(Symmetric(A1pd, :L), 2)) == CHOLMOD.Sparse(ldltfact(A1pd, 2)) + else + @test !CHOLMOD.issym(Sparse(A1pd, 0)) + @test CHOLMOD.ishermitian(Sparse(A1pd, 0)) + @test CHOLMOD.Sparse(cholfact(Hermitian(A1pd, :L))) == CHOLMOD.Sparse(cholfact(A1pd)) + @test CHOLMOD.Sparse(cholfact(Hermitian(A1pd, :L), 2)) == CHOLMOD.Sparse(cholfact(A1pd, 2)) + @test CHOLMOD.Sparse(ldltfact(Hermitian(A1pd, :L))) == CHOLMOD.Sparse(ldltfact(A1pd)) + @test CHOLMOD.Sparse(ldltfact(Hermitian(A1pd, :L), 2)) == CHOLMOD.Sparse(ldltfact(A1pd, 2)) + end + + + ### update! + F = cholfact(A1pd) + CHOLMOD.change_factor(elty, false, false, true, true, F) + @test unsafe_load(F.p).is_ll == 0 + CHOLMOD.change_factor(elty, true, false, true, true, F) + @test unsafe_load(F.p).is_ll == 1 + @test_approx_eq CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality + @test size(F, 3) == 1 + + F = cholfact(A1pdSparse, 2) + @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) + @test CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, 2.0)) == CHOLMOD.Sparse(F) + + F = ldltfact(A1pd) + @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) + @test CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd)) == CHOLMOD.Sparse(F) + + F = ldltfact(A1pdSparse, 2) + @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) + @test CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, 2.0)) == CHOLMOD.Sparse(F) + + ## Low level interface + @test CHOLMOD.nnz(A1Sparse) == nnz(A1) + @test CHOLMOD.speye(5, 5, elty) == eye(elty, 5, 5) + @test CHOLMOD.spzeros(5, 5, 5, elty) == zeros(elty, 5, 5) + if elty <: Real + @test CHOLMOD.copy(A1Sparse, 0, 1) == A1Sparse + @test CHOLMOD.horzcat(A1Sparse, A2Sparse, true) == [A1 A2] + @test CHOLMOD.vertcat(A1Sparse, A2Sparse, true) == [A1; A2] + svec = ones(elty, 1) + @test CHOLMOD.scale(CHOLMOD.Dense(svec), CHOLMOD.SCALAR, A1Sparse) == A1Sparse + svec = ones(elty, 5) + @test_throws DimensionMismatch CHOLMOD.scale(CHOLMOD.Dense(svec), CHOLMOD.SCALAR, A1Sparse) + @test CHOLMOD.scale(CHOLMOD.Dense(svec), CHOLMOD.ROW, A1Sparse) == A1Sparse + @test_throws DimensionMismatch CHOLMOD.scale(CHOLMOD.Dense([svec, 1]), CHOLMOD.ROW, A1Sparse) + @test CHOLMOD.scale(CHOLMOD.Dense(svec), CHOLMOD.COL, A1Sparse) == A1Sparse + @test_throws DimensionMismatch CHOLMOD.scale(CHOLMOD.Dense([svec, 1]), CHOLMOD.COL, A1Sparse) + @test CHOLMOD.scale(CHOLMOD.Dense(svec), CHOLMOD.SYM, A1Sparse) == A1Sparse + @test_throws DimensionMismatch CHOLMOD.scale(CHOLMOD.Dense([svec, 1]), CHOLMOD.SYM, A1Sparse) + else + @test_throws MethodError CHOLMOD.copy(A1Sparse, 0, 1) == A1Sparse + @test_throws MethodError CHOLMOD.horzcat(A1Sparse, A2Sparse, true) == [A1 A2] + @test_throws MethodError CHOLMOD.vertcat(A1Sparse, A2Sparse, true) == [A1; A2] + end + + + if elty <: Real + @test_approx_eq CHOLMOD.ssmult(A1Sparse, A2Sparse, 0, true, true) A1*A2 + @test_approx_eq CHOLMOD.aat(A1Sparse, [0:size(A1,2)-1], 1) A1*A1' + @test_approx_eq CHOLMOD.aat(A1Sparse, [0:1], 1) A1[:,1:2]*A1[:,1:2]' + @test CHOLMOD.copy(A1Sparse, 0, 1) == A1Sparse + end + + @test CHOLMOD.Sparse(CHOLMOD.Dense(A1Sparse)) == A1Sparse +end