diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl index 1d160beb4b116..2a662a63e144a 100644 --- a/base/linalg/diagonal.jl +++ b/base/linalg/diagonal.jl @@ -99,27 +99,28 @@ expm(D::Diagonal) = Diagonal(exp(D.diag)) sqrtm(D::Diagonal) = Diagonal(sqrt(D.diag)) #Linear solver -function \{TD<:Number,TA<:Number}(D::Diagonal{TD}, A::AbstractArray{TA,1}) - m, n = size(A,2)==1 ? (size(A,1),1) : size(A) - m==length(D.diag) || throw(DimensionMismatch()) - (m == 0 || n == 0) && return A - C = Array(typeof(one(TD)/one(TA)),size(A)) +function A_ldiv_B!(D::Diagonal, B::StridedVecOrMat) + m, n = size(B, 1), size(B, 2) + m == length(D.diag) || throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $m rows")) + (m == 0 || n == 0) && return B for j = 1:n for i = 1:m di = D.diag[i] - di==0 && throw(SingularException(i)) - C[i,j] = A[i,j] / di + di == 0 && throw(SingularException(i)) + B[i,j] /= di end end - return C + return B end +\(D::Diagonal, B::StridedMatrix) = scale(1 ./ D.diag, B) +\(D::Diagonal, b::StridedVector) = reshape(scale(1 ./ D.diag, reshape(b, length(b), 1)), length(b)) \(Da::Diagonal, Db::Diagonal) = Diagonal(Db.diag ./ Da.diag) function inv{T}(D::Diagonal{T}) Di = similar(D.diag) for i = 1:length(D.diag) - D.diag[i]==zero(T) && throw(SingularException(i)) - Di[i]=inv(D.diag[i]) + D.diag[i] == zero(T) && throw(SingularException(i)) + Di[i] = inv(D.diag[i]) end Diagonal(Di) end diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl index 795eb3a44dce5..cd9a04e729fd7 100644 --- a/base/linalg/generic.jl +++ b/base/linalg/generic.jl @@ -240,7 +240,7 @@ end function \{TA,TB}(A::AbstractMatrix{TA}, B::AbstractVecOrMat{TB}) TC = typeof(one(TA)/one(TB)) size(A,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(A,1)) rows, but RHS has $(size(B,1)) rows.")) - A_ldiv_B!(factorize(TA == TC ? copy(A) : convert(AbstractMatrix{TC}, A)), TB == TC ? copy(B) : convert(AbstractArray{TC}, B)) + \(factorize(TA == TC ? A : convert(AbstractMatrix{TC}, A)), TB == TC ? copy(B) : convert(AbstractArray{TC}, B)) end \(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b /(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')' diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index 3672f7eb769b1..5475c8c9433f0 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -13,6 +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::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) @@ -45,7 +46,7 @@ A_mul_B!{T<:BlasComplex,S<:AbstractMatrix}(y::StridedMatrix{T}, A::Hermitian{T,S *(A::StridedMatrix, B::HermOrSym) = A*full(B) factorize(A::HermOrSym) = bkfact(A.data, symbol(A.uplo), issym(A)) -\(A::HermOrSym, B::StridedVecOrMat) = \(bkfact(A.data, symbol(A.uplo), issym(A)), B) +\{T,S<:StridedMatrix}(A::HermOrSym{T,S}, B::StridedVecOrMat) = \(bkfact(A.data, symbol(A.uplo), issym(A)), B) inv{T<:BlasFloat,S<:StridedMatrix}(A::Hermitian{T,S}) = Hermitian{T,S}(inv(bkfact(A.data, symbol(A.uplo))), A.uplo) inv{T<:BlasFloat,S<:StridedMatrix}(A::Symmetric{T,S}) = Symmetric{T,S}(inv(bkfact(A.data, symbol(A.uplo), true)), A.uplo) diff --git a/base/sparse.jl b/base/sparse.jl index 3de1216124625..c7eaee447de11 100644 --- a/base/sparse.jl +++ b/base/sparse.jl @@ -1,21 +1,18 @@ -include("sparse/abstractsparse.jl") - module SparseMatrix +using Base: NonTupleType +using Base.Sort: Forward +using Base.LinAlg: AbstractTriangular + importall Base importall Base.LinAlg -import Base.NonTupleType, Base.float, Base.Order, Base.Sort.Forward -import Base.transpose!, Base.ctranspose! -import Base.LinAlg.AbstractTriangular - export SparseMatrixCSC, - blkdiag, dense, diag, diagm, droptol!, dropzeros!, etree, full, - getindex, ishermitian, issparse, issym, istril, istriu, nnz, - setindex!, sparse, sparsevec, spdiagm, speye, spones, - sprand, sprandbool, sprandn, spzeros, symperm, trace, tril, tril!, - triu, triu!, nonzeros, rowvals, nzrange + blkdiag, dense, droptol!, dropzeros!, etree, issparse, nnz, nonzeros, nzrange, + rowvals, sparse, sparsevec, spdiagm, speye, spones, sprand, sprandbool, sprandn, + spzeros, symperm +include("sparse/abstractsparse.jl") include("sparse/sparsematrix.jl") include("sparse/csparse.jl") diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index ca07b0dbc6118..0e3e3d648f857 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -1,143 +1,131 @@ module CHOLMOD -export # types - CholmodDense, - CholmodFactor, - CholmodSparse, - CholmodTriplet, - CholmodDense!, - etree - -import Base: (*), convert, copy, ctranspose, eltype, findnz, 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, isposdef!, logdet, norm, scale, scale! + A_ldiv_B!, cholfact, cholfact!, det, diag, ishermitian, isposdef, + issym, isvalid, ldltfact, logdet + +import Base.SparseMatrix: sparse + +export + Dense, + Factor, + Sparse + +using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, increment!, indtype, decrement, decrement! -importall Base.SparseMatrix -import Base.SparseMatrix: increment, increment!, decrement, decrement! +######### +# Setup # +######### include("cholmod_h.jl") -macro isok(A) - :($A==CHOLMOD_TRUE || throw(CholmodException())) -end +## macro to generate the name of the C function according to the integer type +macro cholmod_name(nm,typ) string("cholmod_", eval(typ) == Int64 ? "l_" : "", nm) end -const chm_ver = Array(Cint, 3) -if dlsym(dlopen("libcholmod"), :cholmod_version) != C_NULL - ccall((:cholmod_version, :libcholmod), Cint, (Ptr{Cint},), chm_ver) -else - ccall((:jl_cholmod_version, :libsuitesparse_wrapper), Cint, (Ptr{Cint},), chm_ver) -end -const chm_com_sz = ccall((:jl_cholmod_common_size,:libsuitesparse_wrapper),Int,()) -const chm_com = fill(0xff, chm_com_sz) -const chm_l_com = fill(0xff, chm_com_sz) -## chm_com and chm_l_com must be initialized at runtime because they contain pointers -## to functions in libc.so, whose addresses can change -function cmn(::Type{Int32}) - if isnan(reinterpret(Float64,chm_com[1:8])[1]) - @isok ccall((:cholmod_start, :libcholmod), Cint, (Ptr{UInt8},), chm_com) - end - chm_com -end -function cmn(::Type{Int64}) - if isnan(reinterpret(Float64,chm_l_com[1:8])[1]) - @isok ccall((:cholmod_l_start, :libcholmod), Cint, (Ptr{UInt8},), chm_l_com) +for Ti in IndexTypes + @eval begin + function common(::Type{$Ti}) + a = fill(0xff, cholmod_com_sz) + @isok ccall((@cholmod_name "start" $Ti + , :libcholmod), Cint, (Ptr{UInt8},), a) + set_print_level(a, 0) # no printing from CHOLMOD by default + a + end end - chm_l_com end -# check the size of SuiteSparse_long -if int(ccall((:jl_cholmod_sizeof_long,:libsuitesparse_wrapper),Csize_t,())) == 4 - const CholmodIndexTypes = (:Int32, ) - typealias CHMITypes Union(Int32) -else - const CholmodIndexTypes = (:Int32, :Int64) - typealias CHMITypes Union(Int32, Int64) +### These offsets are defined in SuiteSparse_wrapper.c +const cholmod_com_offsets = Array(Csize_t, 19) +ccall((:jl_cholmod_common_offsets, :libsuitesparse_wrapper), + Void, (Ptr{Csize_t},), cholmod_com_offsets) +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] +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 +#################### +# Type definitions # +#################### -typealias CHMVTypes Union(Complex128, Float64) -typealias CHMVRealTypes Union(Float64) +# 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. -type CholmodException <: Exception end - -## macro to generate the name of the C function according to the integer type -macro chm_nm(nm,typ) string("cholmod_", eval(typ) == Int64 ? "l_" : "", nm) end - -### A way of examining some of the fields in chm_com -### Probably better to make this a Dict{ASCIIString,Tuple} and -### save the offsets and the lengths and the types. Then the names can be checked. -type ChmCommon - dbound::Float64 - maxrank::Int - supernodal_switch::Float64 - supernodal::Int32 - final_asis::Int32 - final_super::Int32 - final_ll::Int32 - final_pack::Int32 - final_monotonic::Int32 - final_resymbol::Int32 - prefer_zomplex::Int32 # should always be false - prefer_upper::Int32 - print::Int32 # print level. Default: 3 - precise::Int32 # print 16 digits, otherwise 5 - nmethods::Int32 # number of ordering methods - selected::Int32 - postorder::Int32 - itype::Int32 - dtype::Int32 +# Dense +immutable 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 -### These offsets should be reconfigured to be less error-prone in matches -const chm_com_offsets = Array(Int, length(ChmCommon.types)) -ccall((:jl_cholmod_common_offsets, :libsuitesparse_wrapper), - Void, (Ptr{Int},), chm_com_offsets) -const chm_final_ll_inds = (1:4) + chm_com_offsets[7] -const chm_prt_inds = (1:4) + chm_com_offsets[13] -const chm_ityp_inds = (1:4) + chm_com_offsets[18] - -### there must be an easier way but at least this works. -function ChmCommon(aa::Array{UInt8,1}) - typs = ChmCommon.types - sz = map(sizeof, typs) - args = map(i->reinterpret(typs[i], aa[chm_com_offsets[i] + (1:sz[i])])[1], 1:length(sz)) - eval(Expr(:call, unshift!(args, :ChmCommon), Any)) +type Dense{T<:VTypes} <: DenseMatrix{T} + p::Ptr{C_Dense{T}} end -function set_chm_prt_lev(cm::Array{UInt8}, lev::Integer) # can probably be removed - cm[(1:4) + chm_com_offsets[13]] = reinterpret(UInt8, [int32(lev)]) +# Sparse +immutable 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 -## cholmod_dense pointers passed to or returned from C functions are of Julia type -## Ptr{c_CholmodDense}. The CholmodDense type contains a c_CholmodDense object and other -## fields then ensure the memory pointed to is freed when it should be and not before. -type c_CholmodDense{T<:CHMVTypes} - m::Int - n::Int - nzmax::Int - lda::Int - xpt::Ptr{T} - zpt::Ptr{Void} +# 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 +immutable 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 CholmodDense{T<:CHMVTypes} - c::c_CholmodDense - mat::Matrix{T} +type Sparse{Tv<:VTypes,Ti<:ITypes} <: AbstractSparseMatrix{Tv,Ti} + p::Ptr{C_Sparse{Tv,Ti}} end -if (1000chm_ver[1]+chm_ver[2]) >= 2001 # CHOLMOD version 2.1.0 or later - type c_CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} - n::Int - minor::Int +# Factor + +if version >= v"2.1.0" # CHOLMOD version 2.1.0 or later + immutable 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} @@ -145,11 +133,11 @@ if (1000chm_ver[1]+chm_ver[2]) >= 2001 # 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} @@ -162,51 +150,13 @@ if (1000chm_ver[1]+chm_ver[2]) >= 2001 # CHOLMOD version 2.1.0 or later xtype::Cint dtype::Cint end - - type CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} - c::c_CholmodFactor{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 CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodFactor{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 = CholmodFactor{Tv,Ti}(cfp, Perm, ColCount, IPerm, p, i, x, nz, next, prev, - super, pi, px, s) - c_free(cp) - cf - end else - type c_CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} - n::Int - minor::Int + immutable 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} @@ -214,11 +164,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} @@ -231,836 +181,875 @@ else xtype::Cint dtype::Cint end - - type CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} - c::c_CholmodFactor{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 CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodFactor{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 = CholmodFactor{Tv,Ti}(cfp, Perm, ColCount, p, i, x, nz, next, prev, - super, pi, px, s) - c_free(cp) - cf - end end -type c_CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes} - m::Csize_t - n::Csize_t - nzmax::Csize_t - ppt::Ptr{Ti} - ipt::Ptr{Ti} - nzpt::Ptr{Void} - xpt::Ptr{Tv} - zpt::Ptr{Void} - stype::Cint - itype::Cint - xtype::Cint - dtype::Cint - sorted::Cint - packed::Cint +type Factor{Tv,Ti} <: Factorization{Tv} + p::Ptr{C_Factor{Tv,Ti}} end -type CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes} - c::c_CholmodSparse{Tv,Ti} - colptr0::Vector{Ti} - rowval0::Vector{Ti} - nzval::Vector{Tv} +################# +# Thin wrappers # +################# + +# 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 + +### 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 - -type c_CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes} - m::Int - n::Int - nzmax::Int - nnz::Int - i::Ptr{Ti} - j::Ptr{Ti} - x::Ptr{Tv} - z::Ptr{Void} - stype::Cint - itype::Cint - xtype::Cint - dtype::Cint -end - -type CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes} - c::c_CholmodTriplet{Tv,Ti} - i::Vector{Ti} - j::Vector{Ti} - x::Vector{Tv} +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 -eltype{T<:CHMVTypes}(::Type{CholmodDense{T}}) = T -eltype{T<:CHMVTypes}(::Type{CholmodFactor{T}}) = T -eltype{T<:CHMVTypes}(::Type{CholmodSparse{T}}) = T -eltype{T<:CHMVTypes}(::Type{CholmodTriplet{T}}) = T - -## The CholmodDense! constructor does not copy the contents, which is generally what you -## want as most uses of CholmodDense objects are read-only. -function CholmodDense!{T<:CHMVTypes}(aa::VecOrMat{T}) # uses the memory from Julia - m = size(aa,1); n = size(aa,2) - CholmodDense(c_CholmodDense{T}(m, n, m*n, stride(aa,2), convert(Ptr{T}, aa), - C_NULL, xtyp(T), dtyp(T)), - length(size(aa)) == 2 ? aa : reshape(aa, (m,n))) -end +free_dense!{T}(p::Ptr{C_Dense{T}}) = ccall((:cholmod_free_dense, :libcholmod), Cint, (Ptr{Ptr{C_Dense{T}}}, Ptr{Void}), &p, common(Cint)) -## The CholmodDense constructor copies the contents -function CholmodDense{T<:CHMVTypes}(aa::VecOrMat{T}) - m = size(aa,1); n = size(aa,2) - acp = length(size(aa)) == 2 ? copy(aa) : reshape(copy(aa), (m,n)) - CholmodDense(c_CholmodDense{T}(m, n, m*n, stride(aa,2), convert(Ptr{T}, acp), - C_NULL, xtyp(T), dtyp(T)), acp) +function zeros{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) + d = Dense(ccall((:cholmod_zeros, :libcholmod), Ptr{C_Dense{T}}, + (Csize_t, Csize_t, Cint, Ptr{UInt8}), + m, n, xtyp(T), common(Cint))) + finalizer(d, free!) + d end - -function CholmodDense{T<:CHMVTypes}(c::Ptr{c_CholmodDense{T}}) - cp = unsafe_load(c) - if cp.lda != cp.m || cp.nzmax != cp.m * cp.n - throw(DimensionMismatch("overallocated cholmod_dense returned object of size $(cp.m) by $(cp.n) with leading dim $(cp.lda) and nzmax $(cp.nzmax)")) - end - ## the true in the call to pointer_to_array means Julia will free the memory - val = CholmodDense(cp, pointer_to_array(cp.xpt, (cp.m,cp.n), true)) - c_free(c) - val +zeros(m::Integer, n::Integer) = zeros(m, n, Float64) + +function ones{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) + d = Dense(ccall((:cholmod_ones, :libcholmod), Ptr{C_Dense{T}}, + (Csize_t, Csize_t, Cint, Ptr{UInt8}), + m, n, xtyp(T), common(Cint))) + finalizer(d, free!) + d end -CholmodDense!{T<:CHMVTypes}(c::Ptr{c_CholmodDense{T}}) = CholmodDense(c) # no distinction - -function isvalid{T<:CHMVTypes}(cd::CholmodDense{T}) - bool(ccall((:cholmod_check_dense, :libcholmod), Cint, - (Ptr{c_CholmodDense{T}}, Ptr{UInt8}), &cd.c, cmn(Int32))) +ones(m::Integer, n::Integer) = ones(m, n, Float64) + +function eye{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) + d = Dense(ccall((:cholmod_eye, :libcholmod), Ptr{C_Dense{T}}, + (Csize_t, Csize_t, Cint, Ptr{UInt8}), + m, n, xtyp(T), common(Cint))) + finalizer(d, free!) + d end - -function chm_eye{T<:CHMVTypes}(m::Integer, n::Integer, t::T) - CholmodDense(ccall((:cholmod_eye, :libcholmod), Ptr{c_CholmodDense{T}}, - (Int, Int, Cint, Ptr{UInt8}), - m, n,xtyp(T),cmn(Int32))) +eye(m::Integer, n::Integer) = eye(m, n, Float64) +eye(n::Integer) = eye(n, n, Float64) + +function copy_dense{Tv<:VTypes}(A::Dense{Tv}) + 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 -chm_eye(m::Integer, n::Integer) = chm_eye(m, n, 1.) -chm_eye(n::Integer) = chm_eye(n, n, 1.) -function chm_ones{T<:CHMVTypes}(m::Integer, n::Integer, t::T) - CholmodDense(ccall((:cholmod_ones, :libcholmod), Ptr{c_CholmodDense{T}}, - (Int, Int, Cint, Ptr{UInt8}), - m, n, xtyp(T), cmn(Int32))) -end -chm_ones(m::Integer, n::Integer) = chm_ones(m, n, 1.) - -function chm_zeros{T<:CHMVTypes}(m::Integer, n::Integer, t::T) - CholmodDense(ccall((:cholmod_zeros, :libcholmod), Ptr{c_CholmodDense{T}}, - (Int, Int, Cint, Ptr{UInt8}), - m, n, xtyp(T), cmn(Int32))) -end -chm_zeros(m::Integer, n::Integer) = chm_zeros(m, n, 1.) - -function chm_print{T<:CHMVTypes}(cd::CholmodDense{T}, lev::Integer, nm::ASCIIString) - @isok isvalid(cd) - - cm = cmn(Int32) - orig = cm[chm_prt_inds] - cm[chm_prt_inds] = reinterpret(UInt8, [int32(lev)]) - @isok ccall((:cholmod_print_dense, :libcholmod), Cint, - (Ptr{c_CholmodDense{T}}, Ptr{UInt8}, Ptr{UInt8}), - &cd.c, nm, cm) - cm[chm_prt_inds] = orig -end -chm_print(cd::CholmodDense, lev::Integer) = chm_print(cd, lev, "") -chm_print(cd::CholmodDense) = chm_print(cd, int32(4), "") -show(io::IO,cd::CholmodDense) = chm_print(cd, int32(4), "") - -function copy{Tv<:CHMVTypes}(B::CholmodDense{Tv}) - CholmodDense(ccall((:cholmod_copy_dense,:libcholmod), Ptr{c_CholmodDense{Tv}}, - (Ptr{c_CholmodDense{Tv}},Ptr{UInt8}), &B.c, cmn(Int32))) -end - -function norm{Tv<:CHMVTypes}(D::CholmodDense{Tv},p::Real=1) - ccall((:cholmod_norm_dense, :libcholmod), Float64, - (Ptr{c_CholmodDense{Tv}}, Cint, Ptr{UInt8}), - &D.c, p == 1 ? 1 :(p == Inf ? 1 : throw(ArgumentError("p must be 1 or Inf"))),cmn(Int32)) -end - -function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(colpt::Vector{Ti}, - rowval::Vector{Ti}, - nzval::Vector{Tv}, - m::Integer, - n::Integer, - stype::Signed) - bb = colpt[1] - if bb != 0 && bb != 1 throw(DimensionMismatch("colpt[1] is $bb, must be 0 or 1")) end - if any(diff(colpt) .< 0) throw(ArgumentError("elements of colpt must be non-decreasing")) end - if length(colpt) != n + 1 throw(DimensionMismatch("length(colptr) = $(length(colpt)), should be $(n+1)")) end - if bool(bb) # one-based - colpt0 = decrement(colpt) - rowval0 = decrement(rowval) - else - colpt0 = colpt - rowval0 = rowval - end - nz = colpt0[end] - if length(rowval0) != nz || length(nzval) != nz - throw(DimensionMismatch("length(rowval) = $(length(rowval)) and length(nzval) = $(length(nzval)) should be $nz")) - end - if any(rowval0 .< 0) || any(rowval0 .>= m) - throw(ArgumentError("all elements of rowval0 must be in the range [0,$(m-1)]")) +### 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 - it = ityp(Ti) - cs = CholmodSparse(c_CholmodSparse{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), - CHOLMOD_FALSE, CHOLMOD_TRUE), - colpt0, rowval, nzval) - - @isok isvalid(cs) - - cs = sort!(cs) - - return cs -end -function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(A::SparseMatrixCSC{Tv,Ti}, stype::Signed) - CholmodSparse(A.colptr, A.rowval, A.nzval, size(A,1), size(A,2), stype) -end -function CholmodSparse(A::SparseMatrixCSC) - stype = ishermitian(A) ? 1 : 0 - CholmodSparse(stype > 0 ? triu(A) : A, stype) -end -function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodSparse{Tv,Ti}}) - csp = unsafe_load(cp) - colptr0 = pointer_to_array(csp.ppt, (csp.n + 1,), true) - nnz = int(colptr0[end]) - cms = CholmodSparse{Tv,Ti}(csp, colptr0, - pointer_to_array(csp.ipt, (nnz,), true), - pointer_to_array(csp.xpt, (nnz,), true)) - c_free(cp) - cms -end -CholmodSparse{Tv<:CHMVTypes}(D::CholmodDense{Tv}) = CholmodSparse(D,1) # default Ti is Int - -function CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes}(tp::Ptr{c_CholmodTriplet{Tv,Ti}}) - ctp = unsafe_load(tp) - i = pointer_to_array(ctp.i, (ctp.nnz,), true) - j = pointer_to_array(ctp.j, (ctp.nnz,), true) - x = pointer_to_array(ctp.x, (ctp.x == C_NULL ? 0 : ctp.nnz), true) - ct = CholmodTriplet{Tv,Ti}(ctp, i, j, x) - c_free(tp) - ct + ccall((:cholmod_norm_dense, :libcholmod), Cdouble, + (Ptr{C_Dense{Tv}}, Cint, Ptr{UInt8}), + D.p, p, common(Cint)) end -function chm_rdsp(fnm::AbstractString) - fd = ccall(:fopen, Ptr{Void}, (Ptr{UInt8},Ptr{UInt8}), fnm, "r") - res = ccall((:cholmod_read_sparse,:libcholmod), Ptr{c_CholmodSparse{Float64,Cint}}, - (Ptr{Void},Ptr{UInt8}),fd,cmn(Cint)) - ccall(:fclose, Cint, (Ptr{Void},), fd) # should do this in try/finally/end - CholmodSparse(res) +### 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.p, common(Cint))) end -for Ti in CholmodIndexTypes +# Non-Dense wrappers (which all depend on IType) +for Ti in IndexTypes @eval begin - function (*){Tv<:CHMVRealTypes}(A::CholmodSparse{Tv,$Ti}, - B::CholmodSparse{Tv,$Ti}) - CholmodSparse(ccall((@chm_nm "ssmult" $Ti - , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{c_CholmodSparse{Tv,$Ti}}, - Cint,Cint,Cint,Ptr{UInt8}), &A.c,&B.c,0,true,true,cmn($Ti))) + ### 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 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 A_mul_Bc{Tv<:CHMVRealTypes}(A::CholmodSparse{Tv,$Ti}, - B::CholmodSparse{Tv,$Ti}) - cm = cmn($Ti) - aa = Array(Ptr{c_CholmodSparse{Tv,$Ti}}, 2) - if !is(A,B) - aa[1] = ccall((@chm_nm "transpose" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &B.c, 1, cm) - ## result of ssmult will have stype==0, contain numerical values and be sorted - aa[2] = ccall((@chm_nm "ssmult" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, - Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Cint,Cint,Ptr{UInt8}), - &A.c, aa[1], 0, true, true, cm) - @isok ccall((@chm_nm "free_sparse" $Ti - ,:libcholmod), Cint, - (Ptr{Ptr{c_CholmodSparse{Tv,$Ti}}}, Ptr{UInt8}), aa, cm) - return CholmodSparse(aa[2]) - end - - ## 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 - aa[1] = ccall((@chm_nm "copy" $Ti - , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), - &A.c, 0, 1, cm) - - aa[2] = ccall((@chm_nm "aat" $Ti - , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Void}, Int, Cint, Ptr{UInt8}), - aa[1], C_NULL, 0, 1, cm) - - @isok ccall((@chm_nm "free_sparse" $Ti - , :libcholmod), Cint, - (Ptr{Ptr{c_CholmodSparse{Tv,$Ti}}}, Ptr{UInt8}), aa, cm) - - else - aa[2] = ccall((@chm_nm "aat" $Ti - , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Void}, Int, Cint, Ptr{UInt8}), - &A.c, C_NULL, 0, 1, cm) - end - cs = CholmodSparse(aa[2]) + 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 - ## According to the docs conversion to symmetric matrix can be done - ## by changing the stype of the result (stype == 1: upper part). - cs.c.stype = 1 + 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 - return sort!(cs) + function aat{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, fset::Vector{$Ti}, mode::Integer) + 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 Ac_mul_B{Tv<:CHMVRealTypes}(A::CholmodSparse{Tv,$Ti}, - B::CholmodSparse{Tv,$Ti}) - cm = cmn($Ti) - aa = Array(Ptr{c_CholmodSparse{Tv,$Ti}}, 2) - aa[1] = ccall((@chm_nm "transpose" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, 1, cm) - if is(A,B) - Ac = CholmodSparse(aa[1]) - return A_mul_Bc(Ac,Ac) - end - ## result of ssmult will have stype==0, contain numerical values and be sorted - aa[2] = ccall((@chm_nm "ssmult" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, - Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Cint,Cint,Ptr{UInt8}), - aa[1],&B.c,0,true,true,cm) - @isok ccall((@chm_nm "free_sparse" $Ti - , :libcholmod), Cint, - (Ptr{Ptr{c_CholmodSparse{Tv,$Ti}}}, Ptr{UInt8}), aa, cm) - CholmodSparse(aa[2]) + + function sparse_to_dense{Tv<:VTypes}(A::Sparse{Tv,$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 A_mul_Bt{Tv<:CHMVRealTypes}(A::CholmodSparse{Tv,$Ti}, - B::CholmodSparse{Tv,$Ti}) - A_mul_Bc(A,B) # in the unlikely event of writing A*B.' instead of A*B' + 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 At_mul_B{Tv<:CHMVRealTypes}(A::CholmodSparse{Tv,$Ti}, - B::CholmodSparse{Tv,$Ti}) - Ac_mul_B(A,B) # in the unlikely event of writing A.'*B instead of A'*B + + 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 CholmodDense{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - CholmodDense(ccall((@chm_nm "sparse_to_dense" $Ti - ,:libcholmod), Ptr{c_CholmodDense{Tv}}, - (Ptr{c_CholmodSparse{Tv,Ti}},Ptr{UInt8}), - &A.c,cmn($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 CholmodSparse{Tv<:CHMVTypes}(D::CholmodDense{Tv},i::$Ti) - CholmodSparse(ccall((@chm_nm "dense_to_sparse" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodDense{Tv}},Cint,Ptr{UInt8}), - &D.c,true,cmn($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}), + COMPLEX, to_ll, to_super, to_packed, to_monotonic, F.p, common($Ti)) + Factor{Complex{Float64},$Ti}(F.p) end - function CholmodSparse{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti}) - if bool(L.c.is_ll) - return CholmodSparse(ccall((@chm_nm "factor_to_sparse" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{UInt8}), - &L.c,cmn($Ti))) - end - cm = cmn($Ti) - Lcll = ccall((@chm_nm "copy_factor" $Ti - ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}}, - (Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{UInt8}), - &L.c,cm) - @isok ccall((@chm_nm "change_factor" $Ti - ,:libcholmod), Cint, - (Cint,Cint,Cint,Cint,Cint,Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{UInt8}), - L.c.xtype,true,L.c.is_super,true,true,Lcll,cm) - val = CholmodSparse(ccall((@chm_nm "factor_to_sparse" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{UInt8}), - Lcll,cmn($Ti))) - @isok ccall((@chm_nm "free_factor" $Ti - ,:libcholmod), Cint, - (Ptr{Ptr{c_CholmodFactor{Tv,$Ti}}},Ptr{UInt8}), - [Lcll],cm) - val + + 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.p, common($Ti))) end - function CholmodSparse{Tv<:CHMVTypes,Ti<:$Ti}(T::CholmodTriplet{Tv,Ti}) - CholmodSparse(ccall((@chm_nm "triplet_to_sparse" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,Ti}}, - (Ptr{c_CholmodTriplet{Tv,Ti}},Ptr{UInt8}), - &T.c,cmn($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}), + F.p, common($Ti))) end - function CholmodTriplet{Tv<:CHMVTypes,Ti<:$Ti}(A::CholmodSparse{Tv,Ti}) - CholmodTriplet(ccall((@chm_nm "sparse_to_triplet" $Ti - ,:libcholmod), Ptr{c_CholmodTriplet{Tv,Ti}}, - (Ptr{c_CholmodSparse{Tv,Ti}},Ptr{UInt8}), - &A.c,cmn($Ti))) + + function nnz{Tv<:VTypes}(A::Sparse{Tv,$Ti}) + ccall((@cholmod_name("nnz", $Ti),:libcholmod), Int, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, common($Ti)) end - function isvalid{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti}) - bool(ccall((@chm_nm "check_factor" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, cmn($Ti))) + + function speye{Tv<:VTypes}(m::Integer, n::Integer, ::Type{Tv}) + 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 isvalid{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - bool(ccall((@chm_nm "check_sparse" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, cmn($Ti))) + + function spzeros{Tv<:VTypes}(m::Integer, n::Integer, nzmax::Integer, ::Type{Tv}) + 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 isvalid{Tv<:CHMVTypes}(T::CholmodTriplet{Tv,$Ti}) - bool(ccall((@chm_nm "check_triplet" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodTriplet{Tv,$Ti}}, Ptr{UInt8}), - &T.c, cmn($Ti))) + + 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 cholfact{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}, ll::Bool) - cm = cmn($Ti) - ## may need to change final_asis as well as final_ll - if ll cm[chm_final_ll_inds] = reinterpret(UInt8, [one(Cint)]) end - Lpt = ccall((@chm_nm "analyze" $Ti - ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{UInt8}), &A.c, cm) - @isok ccall((@chm_nm "factorize" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodSparse{Tv,$Ti}}, - Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{UInt8}), - &A.c, Lpt, cm) - CholmodFactor(Lpt) + + 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 cholfact{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},beta::Tv,ll::Bool) - cm = cmn($Ti) - ## may need to change final_asis as well as final_ll - if ll cm[chm_final_ll_inds] = reinterpret(UInt8, [one(Cint)]) end - Lpt = ccall((@chm_nm "analyze" $Ti - ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{UInt8}), &A.c, cm) - @isok ccall((@chm_nm "factorize_p" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Tv}, Ptr{Cint}, Csize_t, - Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{UInt8}), - &A.c, &beta, C_NULL, zero(Csize_t), Lpt, cmn($Ti)) - CholmodFactor(Lpt) + function copy_sparse{Tv<:VTypes}(A::Sparse{Tv,$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 chm_analyze{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - CholmodFactor(ccall((@chm_nm "analyze" $Ti - ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{UInt8}), &A.c, cmn($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 - function cholfact!{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti},A::CholmodSparse{Tv,$Ti}) - @isok ccall((@chm_nm "factorize" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodSparse{Tv,$Ti}}, - Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{UInt8}), - &A.c, &L.c, cmn($Ti)) - L + + ### cholmod_check.h ### + function print_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}, name::ASCIIString) + cm = common($Ti) + set_print_level(cm, 3) + @isok ccall((@cholmod_name("print_sparse", $Ti),:libcholmod), Cint, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), + A.p, name, cm) + nothing end - function cholfact!{Tv<:Float64}(L::CholmodFactor{Tv,$Ti},A::CholmodSparse{Tv,$Ti}, - beta::Tv) - @isok ccall((@chm_nm "factorize_p" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{Tv}, Ptr{Cint}, Csize_t, - Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{UInt8}), - &A.c, &beta, C_NULL, zero(Csize_t), &L.c, cmn($Ti)) - L + 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}), + F.p, name, cm) + nothing end - function chm_pack!{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti}) - @isok ccall((@chm_nm "pack_factor" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{UInt8}), - &L.c,cmn($Ti)) - L + + ### 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) + if lA.ncol != lB.nrow + throw(DimensionMismatch("inner matrix dimensions do not fit")) + end + 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 chm_print{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti},lev,nm) - @isok isvalid(L) - - cm = cmn($Ti) - orig = cm[chm_prt_inds] - cm[chm_prt_inds] = reinterpret(UInt8, [int32(lev)]) - @isok ccall((@chm_nm "print_factor" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodFactor{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), - &L.c, nm, cm) - cm[chm_prt_inds] = orig + + function norm_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}, norm::Integer) + if norm != 0 && norm != 1 + throw(ArgumentError("norm argument must be either 0 or 1")) + end + ccall((@cholmod_name("norm_sparse", $Ti), :libcholmod), Cdouble, + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + A.p, norm, common($Ti)) end - function chm_print{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},lev,nm) - @isok isvalid(A) - - cm = cmn($Ti) - orig = cm[chm_prt_inds] - cm[chm_prt_inds] = reinterpret(UInt8, [int32(lev)]) - @isok ccall((@chm_nm "print_sparse" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), - &A.c, nm, cm) - cm[chm_prt_inds] = orig + + 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 chm_scale!{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}, - S::CholmodDense{Tv}, - typ::Integer) - @isok ccall((@chm_nm "scale" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodDense{Tv}},Cint,Ptr{c_CholmodSparse{Tv,$Ti}}, - Ptr{UInt8}), &S.c, typ, &A.c, cmn($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 + if sA.nrow != sA.ncol + throw(DimensionMismatch("matrix must be square")) + elseif 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 + 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.p, scale, A.p, common($Ti)) A end - function chm_sdmult{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}, - trans::Bool, - alpha::Real, - beta::Real, - X::CholmodDense{Tv}) - m,n = size(A) - nc = trans ? m : n - nr = trans ? n : m - if nc != size(X,1) - throw(DimensionMismatch("Incompatible dimensions, $nc and $(size(X,1)), in chm_sdmult")) + + function sdmult!{Tv<:VTypes}(A::Sparse{Tv,$Ti}, transpose::Bool, α::Tv, β::Tv, X::Dense{Tv}, Y::Dense{Tv}) + 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))")) end - aa = float64([alpha, 0.]) - bb = float64([beta, 0.]) - Y = CholmodDense(zeros(Tv,nr,size(X,2))) - @isok ccall((@chm_nm "sdmult" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Ptr{Cdouble},Ptr{Cdouble}, - Ptr{c_CholmodDense{Tv}}, Ptr{c_CholmodDense{Tv}}, Ptr{UInt8}), - &A.c,trans,aa,bb,&X.c,&Y.c,cmn($Ti)) + @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.p, transpose, &α, &β, + X.p, Y.p, common($Ti)) Y end - function chm_speye{Tv<:CHMVTypes,Ti<:$Ti}(m::Ti, n::Ti, x::Tv) - CholmodSparse(ccall((@chm_nm "speye" $Ti - , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Int, Int, Cint, Ptr{UInt8}), - m, n, xtyp(Tv), cmn($Ti))) - end - function chm_spzeros{Tv<:Union(Float64,Complex128)}(m::$Ti, n::$Ti, nzmax::$Ti, x::Tv) - CholmodSparse(ccall((@chm_nm "spzeros" $Ti - , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Int, Int, Int, Cint, Ptr{UInt8}), - m, n, nzmax, xtyp(Tv), cmn($Ti))) - end -## add chm_xtype and chm_pack - function copy{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti}) - CholmodFactor(ccall((@chm_nm "copy_factor" $Ti - ,:libcholmod), Ptr{c_CholmodFactor{Tv,$Ti}}, - (Ptr{c_CholmodFactor{Tv,$Ti}},Ptr{UInt8}), &L.c, cmn($Ti))) - end - function copy{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - CholmodSparse(ccall((@chm_nm "copy_sparse" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{UInt8}), &A.c, cmn($Ti))) - end - function copy{Tv<:CHMVTypes}(T::CholmodTriplet{Tv,$Ti}) - CholmodTriplet(ccall((@chm_nm "copy_triplet" $Ti - ,:libcholmod), Ptr{c_CholmodTriplet{Tv,$Ti}}, - (Ptr{c_CholmodTriplet{Tv,$Ti}},Ptr{UInt8}), &T.c, cmn($Ti))) - end - function ctranspose{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - CholmodSparse(ccall((@chm_nm "transpose" $Ti - ,:libcholmod),Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, 2, cmn($Ti))) - end - function etree{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - tr = Array($Ti,size(A,2)) - @isok ccall((@chm_nm "etree" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{$Ti},Ptr{UInt8}), - &A.c,tr,cmn($Ti)) - tr - end - function hcat{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},B::CholmodSparse{Tv,$Ti}) - ccall((@chm_nm "horzcat" $Ti - , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Ptr{UInt8}), - &A.c,&B.c,true,cmn($Ti)) - end - function nnz{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - ccall((@chm_nm "nnz" $Ti - ,:libcholmod), Int, (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{UInt8}),&A.c,cmn($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 norm{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},p::Real) - ccall((@chm_nm "norm_sparse" $Ti - , :libcholmod), Float64, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c,p == 1 ? 1 :(p == Inf ? 1 : throw(ArgumentError("p must be 1 or Inf"))),cmn($Ti)) + + function symmetry{Tv<:VTypes}(A::Sparse{Tv,$Ti}, option::Integer) + xmatched = Array($Ti, 1) + pmatched = Array($Ti, 1) + nzoffdiag = Array($Ti, 1) + nzdiag = Array($Ti, 1) + rv = ccall((@cholmod_name("symmetry", $Ti), :libcholmod), Cint, + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{$Ti}, Ptr{$Ti}, + Ptr{$Ti}, Ptr{$Ti}, Ptr{UInt8}), + A.p, option, xmatched, pmatched, + nzoffdiag, nzdiag, common($Ti)) + rv, xmatched[1], pmatched[1], nzoffdiag[1], nzdiag[1] end - function solve{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti}, - B::CholmodDense{Tv}, typ::Integer) - 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.")) - CholmodDense(ccall((@chm_nm "solve" $Ti - ,:libcholmod), Ptr{c_CholmodDense{Tv}}, - (Cint, Ptr{c_CholmodFactor{Tv,$Ti}}, - Ptr{c_CholmodDense{Tv}}, Ptr{UInt8}), - typ, &L.c, &B.c, cmn($Ti))) + + # cholmod_cholesky.h + # 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 solve{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti}, - B::CholmodSparse{Tv,$Ti}, - typ::Integer) - 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.")) - CholmodSparse(ccall((@chm_nm "spsolve" $Ti - ,:libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Cint, Ptr{c_CholmodFactor{Tv,$Ti}}, - Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{UInt8}), - typ, &L.c, &B.c, cmn($Ti))) + 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 sort!{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - @isok ccall((@chm_nm "sort" $Ti - ,:libcholmod), Cint, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, cmn($Ti)) - A + 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{Cdouble}, Ptr{$Ti}, Csize_t, + Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + A.p, &β, fset, length(fset), + F.p, cmmn) + F end - function copysym{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - CholmodSparse(ccall((@chm_nm "copy" $Ti - ,:libcholmod),Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Cint,Ptr{UInt8}), - &A.c,0,1,cmn($Ti))) + + function solve{Tv<:VTypes}(sys::Integer, F::Factor{Tv,$Ti}, B::Dense{Tv}) + if 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.")) + end + 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 transpose{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti}) - CholmodSparse(ccall((@chm_nm "transpose" $Ti - ,:libcholmod),Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, 1, cmn($Ti))) + + function spsolve{Tv<:VTypes}(sys::Integer, F::Factor{Tv,$Ti}, B::Sparse{Tv,$Ti}) + if 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.")) + end + 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 - function vcat{Tv<:CHMVTypes}(A::CholmodSparse{Tv,$Ti},B::CholmodSparse{Tv,$Ti}) - ccall((@chm_nm "vertcat" $Ti - , :libcholmod), Ptr{c_CholmodSparse{Tv,$Ti}}, - (Ptr{c_CholmodSparse{Tv,$Ti}},Ptr{c_CholmodSparse{Tv,$Ti}},Cint,Ptr{UInt8}), - &A.c,&B.c,true,cmn($Ti)) + + # Autodetects the types + function read_sparse(file::CFILE, ::Type{$Ti}) + s = Sparse(ccall((@cholmod_name("read_sparse", $Ti), :libcholmod), Ptr{C_SparseVoid}, + (Ptr{Void}, Ptr{UInt8}), + file.ptr, common($Ti))) + finalizer(s, free!) + s end end end -(*){Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::CholmodDense{Tv}) = chm_sdmult(A,false,1.,0.,B) -(*){Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::VecOrMat{Tv}) = chm_sdmult(A,false,1.,0.,CholmodDense(B)) +######################### +# High level interfaces # +######################### -function Ac_mul_B{Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::CholmodDense{Tv}) - chm_sdmult(A,true,1.,0.,B) -end -function Ac_mul_B{Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::VecOrMat{Tv}) - chm_sdmult(A,true,1.,0.,CholmodDense(B)) +# 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) +# 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) -A_ldiv_B!(L::CholmodFactor, B) = L\B # Revisit this to see if allocation can be avoided. It should be possible at least for the right hand side. -(\){T<:CHMVTypes}(L::CholmodFactor{T}, B::CholmodDense{T}) = solve(L, B, CHOLMOD_A) -(\){T<:CHMVTypes}(L::CholmodFactor{T}, b::Vector{T}) = reshape(solve(L, CholmodDense!(b), CHOLMOD_A).mat, length(b)) -(\){T<:CHMVTypes}(L::CholmodFactor{T}, B::Matrix{T}) = solve(L, CholmodDense!(B),CHOLMOD_A).mat -function (\){Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti}) - solve(L,B,CHOLMOD_A) -end -function (\){Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::SparseMatrixCSC{Tv,Ti}) - sparse!(solve(L,CholmodSparse(B, 0),CHOLMOD_A)) -end + # 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 -Ac_ldiv_B{T<:CHMVTypes}(L::CholmodFactor{T},B::CholmodDense{T}) = solve(L,B,CHOLMOD_A) -Ac_ldiv_B{T<:CHMVTypes}(L::CholmodFactor{T},B::VecOrMat{T}) = solve(L,CholmodDense!(B),CHOLMOD_A).mat -function Ac_ldiv_B{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti}) - solve(L,B,CHOLMOD_A) end -function Ac_ldiv_B{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::SparseMatrixCSC{Tv,Ti}) - sparse!(solve(L,CholmodSparse(B),CHOLMOD_A)) +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 + if ishermitian(o) + 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) -cholfact{T<:CHMVTypes}(A::CholmodSparse{T},beta::T) = cholfact(A,beta,false) -cholfact(A::CholmodSparse) = cholfact(A,false) -cholfact(A::SparseMatrixCSC,ll::Bool) = cholfact(CholmodSparse(A),ll) -cholfact(A::SparseMatrixCSC) = cholfact(CholmodSparse(A),false) -function cholfact!{T<:CHMVTypes}(L::CholmodFactor{T},A::CholmodSparse{T},beta::Number) - cholfact!(L,A,convert(T,beta)) + return o end -function cholfact!{T<:CHMVTypes}(L::CholmodFactor{T},A::SparseMatrixCSC{T},beta::Number) - cholfact!(L,CholmodSparse(A),convert(T,beta)) +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 + if ishermitian(o) + change_stype!(o, -1) + end + o end -cholfact!{T<:CHMVTypes}(L::CholmodFactor{T},A::SparseMatrixCSC{T}) = cholfact!(L,CholmodSparse(A)) -chm_analyze(A::SparseMatrixCSC) = chm_analyze(CholmodSparse(A)) +Sparse{Ti<:ITypes}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) +Sparse{Tv<:VTypes,Ti<:ITypes}(A::Hermitian{Tv,SparseMatrixCSC{Tv,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 -chm_print(A::CholmodSparse, lev::Integer) = chm_print(A, lev, "") -chm_print(A::CholmodFactor, lev::Integer) = chm_print(L, lev, "") + # 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 -function chm_scale!{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T},typ::Integer) - chm_scale!(A,CholmodDense(b),typ) - A -end -chm_scale{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T},typ::Integer) = chm_scale!(copy(A),b,typ) + # 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 -chm_speye(m::Integer, n::Integer) = chm_speye(m, n, 1., 1) # default element type is Float32 -chm_speye(n::Integer) = chm_speye(n, n, 1.) # default shape is square + return Sparse(convert(Ptr{C_Sparse{Tv,Ti}}, p)) -chm_spzeros(m::Integer,n::Integer,nzmax::Integer) = chm_spzeros(m,n,nzmax,1.) +end -function scale!{T<:CHMVTypes}(b::Vector{T}, A::CholmodSparse{T}) - chm_scale!(A,CholmodDense(b),CHOLMOD_ROW) +# default is Cint which is probably sufficient when converted from dense matrix +Sparse(A::Dense) = dense_to_sparse(A, Cint) +Sparse(L::Factor) = factor_to_sparse!(copy(L)) +function Sparse(filename::ASCIIString) + f = open(filename) + A = read_sparse(CFILE(f), Int) + close(f) A end -scale{T<:CHMVTypes}(b::Vector{T}, A::CholmodSparse{T}) = scale!(b,copy(A)) -function scale!{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T}) - chm_scale!(A,CholmodDense(b),CHOLMOD_COL) + +## 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) + if s.stype != 0 + throw(ArgumentError("matrix has stype != 0. Convert to matrix with stype == 0 before converting to SparseMatrixCSC")) + end + return 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))) +end +function convert{Ti<:ITypes}(::Type{Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}}, A::Sparse{Float64,Ti}) + s = unsafe_load(A.p) + if !issym(A) + throw(ArgumentError("matrix is not symmetric")) + end + return Symmetric(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))), s.stype > 0 ? :U : :L) +end +function convert{Tv<:VTypes,Ti<:ITypes}(::Type{Hermitian{Tv,SparseMatrixCSC{Tv,Ti}}}, A::Sparse{Tv,Ti}) + s = unsafe_load(A.p) + if !ishermitian(A) + throw(ArgumentError("matrix is not Hermitian")) + end + return Hermitian(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))), s.stype > 0 ? :U : :L) +end +function sparse{Ti}(A::Sparse{Float64,Ti}) # Notice! Cannot be type stable because of stype + s = unsafe_load(A.p) + if s.stype == 0 + return convert(SparseMatrixCSC{Float64,Ti}, A) + end + return convert(Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}, A) +end +function sparse{Ti}(A::Sparse{Complex{Float64},Ti}) # Notice! Cannot be type stable because of stype + s = unsafe_load(A.p) + if s.stype == 0 + return convert(SparseMatrixCSC{Complex{Float64},Ti}, A) + end + return convert(Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}, A) +end +sparse(L::Factor) = sparse(Sparse(L)) +sparse(D::Dense) = sparse(Sparse(D)) + +# 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) A end -scale{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T}) = scale!(copy(A), b) -norm(A::CholmodSparse) = norm(A,1) +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 -show(io::IO,L::CholmodFactor) = chm_print(L,int32(4),"") -show(io::IO,A::CholmodSparse) = chm_print(A,int32(4),"") +isvalid(A::Dense) = check_dense(A) +isvalid(A::Sparse) = check_sparse(A) +isvalid(A::Factor) = check_factor(A) -size(B::CholmodDense) = size(B.mat) -size(B::CholmodDense,d) = size(B.mat,d) -size(A::CholmodSparse) = (int(A.c.m), int(A.c.n)) -function size(A::CholmodSparse, d::Integer) - d == 1 ? A.c.m : (d == 2 ? A.c.n : 1) -end -size(L::CholmodFactor) = (n = int(L.c.n); (n,n)) -size(L::CholmodFactor,d::Integer) = d < 1 ? throw(BoundsError("Dimension $d out of range")) : (d <= 2 ? int(L.c.n) : 1) +copy(A::Dense) = copy_dense(A) +copy(A::Sparse) = copy_sparse(A) +copy(A::Factor) = copy_factor(A) -function solve{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti}, - B::SparseMatrixCSC{Tv,Ti},typ::Integer) - sparse!(solve(L,CholmodSparse(B),typ)) -end -function solve{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::SparseMatrixCSC{Tv,Ti}) - sparse!(solve(L,CholmodSparse(B),CHOLMOD_A)) +function size(A::Union(Dense,Sparse)) + s = unsafe_load(A.p) + return (int(s.nrow), int(s.ncol)) end -function solve{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti}) - solve(L,B,CHOLMOD_A) +function size(F::Factor, i::Integer) + if i < 1 + throw(ArgumentError("dimension must be positive")) + end + s = unsafe_load(F.p) + if i <= 2 + return int(s.n) + end + return 1 end -function (\){Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti}) - solve(L,B,CHOLMOD_A) + +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 solve{T<:CHMVTypes}(L::CholmodFactor{T},B::Vector{T},typ::Integer) - solve(L,CholmodDense!(B),typ).mat[:] +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 -function solve{T<:CHMVTypes}(L::CholmodFactor{T},B::Matrix{T},typ::Integer) - solve(L,CholmodDense!(B),typ).mat + +getindex(A::Sparse, i::Integer) = getindex(A, ind2sub(size(A),i)...) +function getindex{T}(A::Sparse{T}, i0::Integer, i1::Integer) + 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(unsafe_load(s.p, i1) + 1) + r2 = int(unsafe_load(s.p, i1 + 1)) + (r1 > r2) && return zero(T) + 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 -solve{T<:CHMVTypes}(L::CholmodFactor{T},B::CholmodDense{T}) = solve(L,B,CHOLMOD_A) -function findnz{Tv,Ti}(A::CholmodSparse{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 +## 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) + + if !is(A,B) + 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 - 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 + ## 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 + s = unsafe_load(A.p) + if s.stype != 0 + aa1 = copy(A, 0, 1) + return aat(aa1, Ti[0:s.ncol-1;], 1) + else + return aat(A, Ti[0:s.ncol-1;], 1) 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::CholmodFactor) = findnz(CholmodSparse(L)) - -function diag{Tv}(A::CholmodSparse{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 +function Ac_mul_B(A::Sparse, B::Sparse) + aa1 = transpose(A, 2) + if is(A,B) + return A_mul_Bc(aa1, aa1) end - res + ## result of ssmult will have stype==0, contain numerical values and be sorted + return ssmult(aa1, B, 0, true, true) end -function diag{Tv}(L::CholmodFactor{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_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 +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) + 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) + ## 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_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 ldltfact(A::Sparse) + cm = common(indtype(A)) + ## may need to change final_asis as well as final_ll + cm[common_final_ll] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt + F = analyze(A, cm) + factorize!(A, F, cm) + return F +end + +function ldltfact{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 + F = analyze(A, 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 + +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{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(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))) + +Ac_ldiv_B(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B) +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) = 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::CholmodFactor{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 -end -logdet(L::CholmodFactor) = logdet(L, 1:L.c.n) -det(L::CholmodFactor) = exp(logdet(L)) - -full(A::CholmodSparse) = CholmodDense(A).mat -full(A::CholmodDense) = A.mat -function sparse(A::CholmodSparse) - if bool(A.c.stype) return sparse!(copysym(A)) end - SparseMatrixCSC(A.c.m, A.c.n, increment(A.colptr0), increment(A.rowval0), copy(A.nzval)) -end -function sparse!(A::CholmodSparse) - SparseMatrixCSC(A.c.m, A.c.n, increment!(A.colptr0), increment!(A.rowval0), A.nzval) + for d in diag(F) res += log(abs(d)) end + bool(f.is_ll) ? 2res : res end -sparse(L::CholmodFactor) = sparse!(CholmodSparse(L)) -sparse(D::CholmodDense) = sparse!(CholmodSparse(D)) -sparse(T::CholmodTriplet) = sparse!(CholmodSparse(T)) -isposdef!{Tv<:CHMVTypes,Ti}(A::SparseMatrixCSC{Tv,Ti}) = ishermitian(A) && cholfact(A).c.minor == size(A,1) +det(L::Factor) = exp(logdet(L)) -end #module +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 -# placing factorize here for now. Maybe add a new file -function factorize(A::SparseMatrixCSC) - m, n = size(A) - if m == n - Ac = cholfact(A) - Ac.c.minor == m && ishermitian(A) && return Ac +function issym(A::Sparse) + s = unsafe_load(A.p) + if s.stype != 0 + return isreal(A) end - return lufact(A) + 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 s.stype != 0 + return true + else + 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 s.stype != 0 + return true + else + 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 +end + +end #module diff --git a/base/sparse/cholmod_h.jl b/base/sparse/cholmod_h.jl index 52c3572f9f966..f97dff5effd6e 100644 --- a/base/sparse/cholmod_h.jl +++ b/base/sparse/cholmod_h.jl @@ -1,36 +1,37 @@ ## CHOLMOD -const CHOLMOD_TRUE = int32(1) -const CHOLMOD_FALSE = int32(0) +const TRUE = int32(1) +const FALSE = int32(0) ## itype defines the types of integer used: -const CHOLMOD_INT = int32(0) # all integer arrays are int -const CHOLMOD_LONG = int32(2) # all integer arrays are UF_long -ityp(::Type{Int32}) = CHOLMOD_INT -ityp(::Type{Int64}) = CHOLMOD_LONG +const INT = int32(0) # all integer arrays are int +const INTLONG = int32(1) # most are int, some are SuiteSparse_long +const LONG = int32(2) # all integer arrays are SuiteSparse_long +ityp(::Type{Int32}) = INT +ityp(::Type{Int64}) = LONG ## dtype defines what the numerical type is (double or float): -const CHOLMOD_DOUBLE = int32(0) # all numerical values are double -const CHOLMOD_SINGLE = int32(1) # all numerical values are float -dtyp(::Type{Float32}) = CHOLMOD_SINGLE -dtyp(::Type{Float64}) = CHOLMOD_DOUBLE -dtyp(::Type{Complex64}) = CHOLMOD_SINGLE -dtyp(::Type{Complex128}) = CHOLMOD_DOUBLE +const DOUBLE = int32(0) # all numerical values are double +const SINGLE = int32(1) # all numerical values are float +dtyp(::Type{Float32}) = SINGLE +dtyp(::Type{Float64}) = DOUBLE +dtyp(::Type{Complex64}) = SINGLE +dtyp(::Type{Complex128}) = DOUBLE ## xtype defines the kind of numerical values used: -const CHOLMOD_PATTERN = int32(0) # pattern only, no numerical values -const CHOLMOD_REAL = int32(1) # a real matrix -const CHOLMOD_COMPLEX = int32(2) # a complex matrix (ANSI C99 compatible) -const CHOLMOD_ZOMPLEX = int32(3) # a complex matrix (MATLAB compatible) -xtyp(::Type{Float32}) = CHOLMOD_REAL -xtyp(::Type{Float64}) = CHOLMOD_REAL -xtyp(::Type{Complex64}) = CHOLMOD_COMPLEX -xtyp(::Type{Complex128}) = CHOLMOD_COMPLEX +const PATTERN = int32(0) # pattern only, no numerical values +const REAL = int32(1) # a real matrix +const COMPLEX = int32(2) # a complex matrix (ANSI C99 compatible) +const ZOMPLEX = int32(3) # a complex matrix (MATLAB compatible) +xtyp(::Type{Float32}) = REAL +xtyp(::Type{Float64}) = REAL +xtyp(::Type{Complex64}) = COMPLEX +xtyp(::Type{Complex128}) = COMPLEX ## Scaling modes, selected by the scale input parameter: -const CHOLMOD_SCALAR = int32(0) # A = s*A -const CHOLMOD_ROW = int32(1) # A = diag(s)*A -const CHOLMOD_COL = int32(2) # A = A*diag(s) -const CHOLMOD_SYM = int32(3) # A = diag(s)*A*diag(s) +const SCALAR = int32(0) # A = s*A +const ROW = int32(1) # A = diag(s)*A +const COL = int32(2) # A = A*diag(s) +const SYM = int32(3) # A = diag(s)*A*diag(s) ## Types of systems to solve const CHOLMOD_A = int32(0) # solve Ax=b @@ -42,3 +43,44 @@ const CHOLMOD_Lt = int32(5) # solve L'x=b const CHOLMOD_D = int32(6) # solve Dx=b const CHOLMOD_P = int32(7) # permute x=Px const CHOLMOD_Pt = int32(8) # permute x=P'x + +# Symmetry types +const EMPTY =-1 +const MM_RECTANGULAR = 1 +const MM_UNSYMMETRIC = 2 +const MM_SYMMETRIC = 3 +const MM_HERMITIAN = 4 +const MM_SKEW_SYMMETRIC = 5 +const MM_SYMMETRIC_POSDIAG = 6 +const MM_HERMITIAN_POSDIAG = 7 + +# check the size of SuiteSparse_long +if int(ccall((:jl_cholmod_sizeof_long, :libsuitesparse_wrapper),Csize_t,())) == 4 + const SuiteSparse_long = Int32 + const IndexTypes = (:Int32, ) + typealias ITypes Union(Int32) +else + const SuiteSparse_long = Int64 + const IndexTypes = (:Int32, :Int64) + typealias ITypes Union(Int32, Int64) +end + +typealias VTypes Union(Complex128, Float64) +typealias VRealTypes Union(Float64) + +type CHOLMODException <: Exception + msg::AbstractString +end + +macro isok(A) + :($A == TRUE || throw(CHOLMODException(""))) +end + +const version_array = Array(Cint, 3) +if dlsym(dlopen("libcholmod"), :cholmod_version) != C_NULL + ccall((:cholmod_version, :libcholmod), Cint, (Ptr{Cint},), version_array) +else + ccall((:jl_cholmod_version, :libsuitesparse_wrapper), Cint, (Ptr{Cint},), version_array) +end +const version = VersionNumber(version_array...) +const cholmod_com_sz = ccall((:jl_cholmod_common_size,:libsuitesparse_wrapper),Int,()) diff --git a/base/sparse/linalg.jl b/base/sparse/linalg.jl index 1e78e54b847a6..3b5249c313699 100644 --- a/base/sparse/linalg.jl +++ b/base/sparse/linalg.jl @@ -673,6 +673,41 @@ scale{Tv,Ti,T}(A::SparseMatrixCSC{Tv,Ti}, b::Vector{T}) = scale{T,Tv,Ti}(b::Vector{T}, A::SparseMatrixCSC{Tv,Ti}) = scale!(similar(A, promote_type(Tv,T)), b, A) +function factorize(A::SparseMatrixCSC) + m, n = size(A) + if m == n + AC = CHOLMOD.Sparse(A) + if ishermitian(AC) + try + return cholfact(A) + catch e + isa(e, PosDefException) || e + return ldltfact(A) + end + end + return lufact(A) + else + throw(ArgumentError("sparse least squares problems by QR are not handled yet")) + end +end + +function factorize{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) + try + return cholfact(A) + catch e + isa(e, PosDefException) || e + return ldltfact(A) + end +end +function factorize{Ti}(A::Hermitian{Complex{Float64}, SparseMatrixCSC{Complex{Float64},Ti}}) + try + return cholfact(A) + catch e + isa(e, PosDefException) || e + return ldltfact(A) + end +end + chol(A::SparseMatrixCSC) = error("Use cholfact() instead of chol() for sparse matrices.") lu(A::SparseMatrixCSC) = error("Use lufact() instead of lu() for sparse matrices.") eig(A::SparseMatrixCSC) = error("Use eigs() instead of eig() for sparse matrices.") diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index 8deada94d0e92..c21ce36525163 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -22,7 +22,7 @@ nonzeros(S::SparseMatrixCSC) = S.nzval rowvals(S::SparseMatrixCSC) = S.rowval nzrange(S::SparseMatrixCSC, col::Integer) = S.colptr[col]:(S.colptr[col+1]-1) -function Base.showarray(io::IO, S::SparseMatrixCSC; +function showarray(io::IO, S::SparseMatrixCSC; header::Bool=true, limit::Bool=Base._limit_output, rows = Base.tty_size()[1], repr=false) # TODO: repr? diff --git a/test/Makefile b/test/Makefile index 2fe67f6e0ca40..d75a284b0c4d3 100644 --- a/test/Makefile +++ b/test/Makefile @@ -1,7 +1,7 @@ JULIAHOME = $(abspath ..) include ../Make.inc -TESTS = all linalg $(filter-out TestHelpers runtests testdefs,$(subst .jl,,$(wildcard *.jl))) +TESTS = all linalg sparse $(filter-out TestHelpers runtests testdefs,$(subst .jl,,$(wildcard *.jl))) default: all diff --git a/test/linalg/cholmod.jl b/test/linalg/cholmod.jl deleted file mode 100644 index bbe53af41b4bc..0000000000000 --- a/test/linalg/cholmod.jl +++ /dev/null @@ -1,38 +0,0 @@ -using Base.Test - -let # Issue 9160 - const CHOLMOD = Base.SparseMatrix.CHOLMOD - - for Ti in CHOLMOD.CHMITypes.types - for elty in CHOLMOD.CHMVRealTypes.types - - A = sprand(10,10,0.1) - A = convert(SparseMatrixCSC{elty,Ti},A) - cmA = CHOLMOD.CholmodSparse(A) - - B = sprand(10,10,0.1) - B = convert(SparseMatrixCSC{elty,Ti},B) - cmB = CHOLMOD.CholmodSparse(B) - - # Ac_mul_B - @test_approx_eq sparse(cmA'*cmB) A'*B - - # A_mul_Bc - @test_approx_eq sparse(cmA*cmB') A*B' - - # A_mul_Ac - @test_approx_eq sparse(cmA*cmA') A*A' - - # Ac_mul_A - @test_approx_eq sparse(cmA'*cmA) A'*A - - # A_mul_Ac for symmetric A - A = 0.5*(A + A') - cmA = CHOLMOD.CholmodSparse(A) - @test_approx_eq sparse(cmA*cmA') A*A' - end - end -end - -#9915 -@test speye(2)\speye(2) == eye(2) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 8ddcfd4d277f4..f499145432a14 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,7 @@ testnames = [ "subarray", "reduce", "reducedim", "random", "intfuncs", "simdloop", "blas", "fft", "dsp", "sparse", "bitarray", "copy", "math", "fastmath", "functional", "bigint", "sorting", "statistics", "spawn", - "backtrace", "priorityqueue", "arpack", "file", "suitesparse", "version", + "backtrace", "priorityqueue", "arpack", "file", "version", "resolve", "pollfd", "mpfr", "broadcast", "complex", "socket", "floatapprox", "readdlm", "reflection", "regex", "float16", "combinatorics", "sysinfo", "rounding", "ranges", "mod2pi", "euler", "show", @@ -25,10 +25,16 @@ push!(testnames, "parallel") tests = (ARGS==["all"] || isempty(ARGS)) ? testnames : ARGS +if "sparse" in tests + # specifically selected case + filter!(x -> x != "sparse", tests) + prepend!(tests, ["sparse/sparse", "sparse/cholmod", "sparse/umfpack"]) +end + if "linalg" in tests # specifically selected case filter!(x -> x != "linalg", tests) - prepend!(tests, ["linalg1", "linalg2", "linalg3", "linalg4", "linalg/lapack", "linalg/triangular", "linalg/tridiag", "linalg/pinv", "linalg/cholmod", "linalg/umfpack", "linalg/givens"]) + prepend!(tests, ["linalg1", "linalg2", "linalg3", "linalg4", "linalg/lapack", "linalg/triangular", "linalg/tridiag", "linalg/pinv", "linalg/givens"]) end net_required_for = ["socket", "parallel"] diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl new file mode 100644 index 0000000000000..947b4e4529c9e --- /dev/null +++ b/test/sparse/cholmod.jl @@ -0,0 +1,432 @@ +using Base.Test + +using Base.SparseMatrix.CHOLMOD + +# based on deps/SuiteSparse-4.0.2/CHOLMOD/Demo/ + +# chm_rdsp(joinpath(JULIA_HOME, "../../deps/SuiteSparse-4.0.2/CHOLMOD/Demo/Matrix/bcsstk01.tri")) +# because the file may not exist in binary distributions and when a system suitesparse library +# is used + +## Result from C program +## ---------------------------------- cholmod_demo: +## norm (A,inf) = 3.57095e+09 +## norm (A,1) = 3.57095e+09 +## CHOLMOD sparse: A: 48-by-48, nz 224, upper. OK +## CHOLMOD dense: B: 48-by-1, OK +## bnorm 1.97917 +## Analyze: flop 6009 lnz 489 +## Factorizing A +## CHOLMOD factor: L: 48-by-48 simplicial, LDL'. nzmax 489. nz 489 OK +## Ordering: AMD fl/lnz 12.3 lnz/anz 2.2 +## ints in L: 782, doubles in L: 489 +## factor flops 6009 nnz(L) 489 (w/no amalgamation) +## nnz(A*A'): 224 +## flops / nnz(L): 12.3 +## nnz(L) / nnz(A): 2.2 +## analyze cputime: 0.0000 +## factor cputime: 0.0000 mflop: 0.0 +## solve cputime: 0.0000 mflop: 0.0 +## overall cputime: 0.0000 mflop: 0.0 +## peak memory usage: 0 (MB) +## residual 2.5e-19 (|Ax-b|/(|A||x|+|b|)) +## residual 1.3e-19 (|Ax-b|/(|A||x|+|b|)) after iterative refinement +## rcond 9.5e-06 + +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, + 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, + 9,26,32,8,9,25,27,31,32,33,10,24,28,30,32,34,6,11,29,30,31,35,12,17,30,36, + 13,31,35,37,14,15,32,34,38,14,15,33,37,38,39,16,32,34,36,38,40,12,17,31,35, + 36,37,41,12,16,17,18,23,36,40,42,13,14,15,19,37,39,43,13,14,15,20,21,38,43, + 44,13,14,15,20,21,37,39,43,44,45,12,16,17,22,36,40,42,46,12,16,17,18,23,41, + 42,46,47]), + [2.83226851852e6,1.63544753086e6,1.72436728395e6,-2.0e6,-2.08333333333e6, + 1.00333333333e9,1.0e6, -2.77777777778e6,1.0675e9,2.08333333333e6, + 5.55555555555e6,1.53533333333e9,-3333.33333333,-1.0e6,2.83226851852e6, + -6666.66666667,2.0e6,1.63544753086e6,-1.68e6,1.72436728395e6,-2.0e6,4.0e8,2.0e6, + -2.08333333333e6,1.00333333333e9,1.0e6,2.0e8,-1.0e6,-2.77777777778e6,1.0675e9, + -2.0e6,2.08333333333e6,5.55555555555e6,1.53533333333e9,-2.8e6,2.8360994695e6, + -30864.1975309,-5.55555555555e6,1.76741074446e6,-15432.0987654,2.77777777778e6, + 517922.131816,3.89003806848e6,-3.33333333333e6,4.29857058902e6,-2.6349902747e6, + 1.97572063531e9,-2.77777777778e6,3.33333333333e8,-2.14928529451e6, + 2.77777777778e6,1.52734651547e9,5.55555555555e6,6.66666666667e8,2.35916180402e6, + -5.55555555555e6,-1.09779731332e8,1.56411143711e9,-2.8e6,-3333.33333333,1.0e6, + 2.83226851852e6,-30864.1975309,-5.55555555555e6,-6666.66666667,-2.0e6, + 1.63544753086e6,-15432.0987654,2.77777777778e6,-1.68e6,1.72436728395e6, + -3.33333333333e6,2.0e6,4.0e8,-2.0e6,-2.08333333333e6,1.00333333333e9, + -2.77777777778e6,3.33333333333e8,-1.0e6,2.0e8,1.0e6,2.77777777778e6,1.0675e9, + 5.55555555555e6,6.66666666667e8,-2.0e6,2.08333333333e6,-5.55555555555e6, + 1.53533333333e9,-28935.1851852,-2.08333333333e6,60879.6296296,-1.59791666667e6, + 3.37291666667e6,-28935.1851852,2.08333333333e6,2.41171296296e6,-2.08333333333e6, + 1.0e8,-2.5e6,-416666.666667,1.5e9,-833333.333333,1.25e6,5.01833333333e8, + 2.08333333333e6,1.0e8,416666.666667,5.025e8,-28935.1851852,-2.08333333333e6, + -4166.66666667,-1.25e6,3.98587962963e6,-1.59791666667e6,-8333.33333333,2.5e6, + 3.41149691358e6,-28935.1851852,2.08333333333e6,-2.355e6,2.43100308642e6, + -2.08333333333e6,1.0e8,-2.5e6,5.0e8,2.5e6,-416666.666667,1.50416666667e9, + -833333.333333,1.25e6,2.5e8,-1.25e6,-3.47222222222e6,1.33516666667e9, + 2.08333333333e6,1.0e8,-2.5e6,416666.666667,6.94444444444e6,2.16916666667e9, + -28935.1851852,-2.08333333333e6,-3.925e6,3.98587962963e6,-1.59791666667e6, + -38580.2469136,-6.94444444444e6,3.41149691358e6,-28935.1851852,2.08333333333e6, + -19290.1234568,3.47222222222e6,2.43100308642e6,-2.08333333333e6,1.0e8, + -4.16666666667e6,2.5e6,-416666.666667,1.50416666667e9,-833333.333333, + -3.47222222222e6,4.16666666667e8,-1.25e6,3.47222222222e6,1.33516666667e9, + 2.08333333333e6,1.0e8,6.94444444445e6,8.33333333333e8,416666.666667, + -6.94444444445e6,2.16916666667e9,-3830.95098171,1.14928529451e6,-275828.470683, + -28935.1851852,-2.08333333333e6,-4166.66666667,1.25e6,64710.5806113, + -131963.213599,-517922.131816,-2.29857058902e6,-1.59791666667e6,-8333.33333333, + -2.5e6,3.50487988027e6,-517922.131816,-2.16567078453e6,551656.941366, + -28935.1851852,2.08333333333e6,-2.355e6,517922.131816,4.57738374749e6, + 2.29857058902e6,-551656.941367,4.8619365099e8,-2.08333333333e6,1.0e8,2.5e6, + 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], 1) +@test_approx_eq CHOLMOD.norm_sparse(A, 0) 3.570948074697437e9 +@test_approx_eq CHOLMOD.norm_sparse(A, 1) 3.570948074697437e9 +@test_throws ArgumentError CHOLMOD.norm_sparse(A, 2) +@test isvalid(A) + +B = A * ones(size(A,2)) +chma = ldltfact(A) # LDL' form +@test isvalid(chma) +x = chma\B +@test_approx_eq x ones(size(x)) + +chma = cholfact(A) # LL' form +@test isvalid(chma) +x = chma\B +@test_approx_eq x ones(size(x)) + +#lp_afiro example +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, + 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]), + [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, + -1.0,-1.06,1.0,0.301,1.0,-1.0,1.0,-1.0,1.0,1.0,-1.0,-1.06,1.0,0.301,-1.0, + -1.06,1.0,0.313,-1.0,-0.96,1.0,0.313,-1.0,-0.86,1.0,0.326,-1.0,2.364,-1.0, + 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], 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 * (convert(Matrix, y) - convert(Matrix, pred))) < 1e-8 + +let # Issue 9160 + + for Ti in CHOLMOD.ITypes.types + for elty in CHOLMOD.VRealTypes.types + + A = sprand(10,10,0.1) + A = convert(SparseMatrixCSC{elty,Ti},A) + cmA = CHOLMOD.Sparse(A) + + B = sprand(10,10,0.1) + B = convert(SparseMatrixCSC{elty,Ti},B) + cmB = CHOLMOD.Sparse(B) + + # Ac_mul_B + @test_approx_eq sparse(cmA'*cmB) A'*B + + # A_mul_Bc + @test_approx_eq sparse(cmA*cmB') A*B' + + # A_mul_Ac + @test_approx_eq sparse(cmA*cmA') A*A' + + # Ac_mul_A + @test_approx_eq sparse(cmA'*cmA) A'*A + + # A_mul_Ac for symmetric A + A = 0.5*(A + A') + cmA = CHOLMOD.Sparse(A) + @test_approx_eq sparse(cmA*cmA') A*A' + end + end +end + + +# Issue #9915 +@test speye(2)\speye(2) == eye(2) + +# test eltype +@test eltype(Dense(ones(3))) == Float64 +@test eltype(A) == Float64 +@test eltype(chma) == Float64 + +# test Sparse constructor Symmetric and Hermitian input (and issym and ishermitian) +ACSC = sprandn(10, 10, 0.3) + I +@test issym(Sparse(Symmetric(ACSC, :L))) +@test issym(Sparse(Symmetric(ACSC, :U))) +@test ishermitian(Sparse(Hermitian(complex(ACSC), :L))) +@test ishermitian(Sparse(Hermitian(complex(ACSC), :U))) + +# test Sparse constructor for c_Sparse{Tv,Ti} input( and Sparse*Sparse) +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 sparse(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 sparse(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)) +puint = convert(Ptr{Uint32}, p) +unsafe_store!(puint, CHOLMOD.SINGLE, 3*div(sizeof(Csize_t), 4) + 5*div(sizeof(Ptr{Void}), 4) + 4) +@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)) +puint = convert(Ptr{Uint32}, p) +unsafe_store!(puint, 5, 3*div(sizeof(Csize_t), 4) + 5*div(sizeof(Ptr{Void}), 4) + 4) +@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)) +puint = convert(Ptr{Uint32}, p) +unsafe_store!(puint, 3, 3*div(sizeof(Csize_t), 4) + 5*div(sizeof(Ptr{Void}), 4) + 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)) +puint = convert(Ptr{Uint32}, p) +unsafe_store!(puint, CHOLMOD.INTLONG, 3*div(sizeof(Csize_t), 4) + 5*div(sizeof(Ptr{Void}), 4) + 2) +@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)) +puint = convert(Ptr{Uint32}, p) +unsafe_store!(puint, 5, 3*div(sizeof(Csize_t), 4) + 5*div(sizeof(Ptr{Void}), 4) + 2) +@test_throws CHOLMOD.CHOLMODException CHOLMOD.Sparse(p) + +# test that Sparse(Ptr) works for SuiteSparse_long (on 64 bit systems) +if CHOLMOD.SuiteSparse_long == Int64 + 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}) +end + +# 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_throws BoundsError ADense[6, 1] + @test_throws BoundsError ADense[1, 6] + @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_throws BoundsError A1Sparse[6, 1] + @test_throws BoundsError A1Sparse[1, 6] + @test sparse(A1Sparse) == A1 + @test CHOLMOD.sparse(CHOLMOD.Sparse(Hermitian(A1, :L))) == Hermitian(A1, :L) + @test CHOLMOD.sparse(CHOLMOD.Sparse(Hermitian(A1, :U))) == Hermitian(A1, :U) + @test_throws ArgumentError convert(SparseMatrixCSC{elty,Int}, A1pdSparse) + if elty <: Real + @test_throws ArgumentError convert(Symmetric{Float64,SparseMatrixCSC{Float64,Int}}, A1Sparse) + else + @test_throws ArgumentError convert(Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Int}}, A1Sparse) + end + @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_throws DimensionMismatch CHOLMOD.Sparse(A1[:,1:4])*A2Sparse + @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' + + @test_throws DimensionMismatch A1Sparse*CHOLMOD.eye(4, 5, elty) + end + + # Factor + @test_throws ArgumentError cholfact(A1) + @test_throws Base.LinAlg.PosDefException cholfact(A1 + A1' - I) + 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_throws DimensionMismatch F\CHOLMOD.Dense(ones(elty, 4)) + @test_throws DimensionMismatch F\CHOLMOD.Sparse(sparse(ones(elty, 4))) + @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 logdet(F) logdet(full(A1pd)) + @test det(F) == exp(logdet(F)) + let # to test supernodal, we must use a larger matrix + Ftmp = sprandn(100,100,0.1) + Ftmp = Ftmp'Ftmp + I + @test_approx_eq logdet(cholfact(Ftmp)) logdet(full(Ftmp)) + end + @test_approx_eq logdet(ldltfact(A1pd)) logdet(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_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, 2) == 5 + @test size(F, 3) == 1 + @test_throws ArgumentError size(F, 0) + + F = cholfact(A1pdSparse, 2) + @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) + @test_approx_eq CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, 2.0)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality + + F = ldltfact(A1pd) + @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) + @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 + + F = ldltfact(A1pdSparse, 2) + @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) + @test_approx_eq CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, 2.0)) CHOLMOD.Sparse(F) # surprisingly, this can cause small ulp size changes so we cannot test exact equality + + @test isa(CHOLMOD.factor_to_sparse!(F), CHOLMOD.Sparse) + @test_throws CHOLMOD.CHOLMODException CHOLMOD.factor_to_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) + @test_throws DimensionMismatch CHOLMOD.scale!(CHOLMOD.Dense(svec), CHOLMOD.SYM, CHOLMOD.Sparse(A1[:,1:4])) + 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 diff --git a/test/sparse.jl b/test/sparse/sparse.jl similarity index 99% rename from test/sparse.jl rename to test/sparse/sparse.jl index be329e978265f..7a24fdf1420ea 100644 --- a/test/sparse.jl +++ b/test/sparse/sparse.jl @@ -667,3 +667,9 @@ let D = Diagonal(ones(10,10)), @test countnz(sparse(Diagonal(Int[]))) == 0 @test countnz(sparsevec(Diagonal(Int[]))) == 0 end + +# explicit zeros +a = SparseMatrixCSC(2, 2, [1, 3, 5], [1, 2, 1, 2], [1.0, 0.0, 0.0, 1.0]) +@test_approx_eq lufact(a)\[2.0, 3.0] [2.0, 3.0] +@test_approx_eq cholfact(a)\[2.0, 3.0] [2.0, 3.0] + diff --git a/test/linalg/umfpack.jl b/test/sparse/umfpack.jl similarity index 100% rename from test/linalg/umfpack.jl rename to test/sparse/umfpack.jl diff --git a/test/suitesparse.jl b/test/suitesparse.jl deleted file mode 100644 index a999b3c526d88..0000000000000 --- a/test/suitesparse.jl +++ /dev/null @@ -1,127 +0,0 @@ -using Base.SparseMatrix.CHOLMOD - -# based on deps/SuiteSparse-4.0.2/CHOLMOD/Demo/ - -# chm_rdsp(joinpath(JULIA_HOME, "../../deps/SuiteSparse-4.0.2/CHOLMOD/Demo/Matrix/bcsstk01.tri")) -# because the file may not exist in binary distributions and when a system suitesparse library -# is used - -## Result from C program -## ---------------------------------- cholmod_demo: -## norm (A,inf) = 3.57095e+09 -## norm (A,1) = 3.57095e+09 -## CHOLMOD sparse: A: 48-by-48, nz 224, upper. OK -## CHOLMOD dense: B: 48-by-1, OK -## bnorm 1.97917 -## Analyze: flop 6009 lnz 489 -## Factorizing A -## CHOLMOD factor: L: 48-by-48 simplicial, LDL'. nzmax 489. nz 489 OK -## Ordering: AMD fl/lnz 12.3 lnz/anz 2.2 -## ints in L: 782, doubles in L: 489 -## factor flops 6009 nnz(L) 489 (w/no amalgamation) -## nnz(A*A'): 224 -## flops / nnz(L): 12.3 -## nnz(L) / nnz(A): 2.2 -## analyze cputime: 0.0000 -## factor cputime: 0.0000 mflop: 0.0 -## solve cputime: 0.0000 mflop: 0.0 -## overall cputime: 0.0000 mflop: 0.0 -## peak memory usage: 0 (MB) -## residual 2.5e-19 (|Ax-b|/(|A||x|+|b|)) -## residual 1.3e-19 (|Ax-b|/(|A||x|+|b|)) after iterative refinement -## rcond 9.5e-06 - -A = CholmodSparse(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, - 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, - 9,26,32,8,9,25,27,31,32,33,10,24,28,30,32,34,6,11,29,30,31,35,12,17,30,36, - 13,31,35,37,14,15,32,34,38,14,15,33,37,38,39,16,32,34,36,38,40,12,17,31,35, - 36,37,41,12,16,17,18,23,36,40,42,13,14,15,19,37,39,43,13,14,15,20,21,38,43, - 44,13,14,15,20,21,37,39,43,44,45,12,16,17,22,36,40,42,46,12,16,17,18,23,41, - 42,46,47]), - [2.83226851852e6,1.63544753086e6,1.72436728395e6,-2.0e6,-2.08333333333e6, - 1.00333333333e9,1.0e6, -2.77777777778e6,1.0675e9,2.08333333333e6, - 5.55555555555e6,1.53533333333e9,-3333.33333333,-1.0e6,2.83226851852e6, - -6666.66666667,2.0e6,1.63544753086e6,-1.68e6,1.72436728395e6,-2.0e6,4.0e8,2.0e6, - -2.08333333333e6,1.00333333333e9,1.0e6,2.0e8,-1.0e6,-2.77777777778e6,1.0675e9, - -2.0e6,2.08333333333e6,5.55555555555e6,1.53533333333e9,-2.8e6,2.8360994695e6, - -30864.1975309,-5.55555555555e6,1.76741074446e6,-15432.0987654,2.77777777778e6, - 517922.131816,3.89003806848e6,-3.33333333333e6,4.29857058902e6,-2.6349902747e6, - 1.97572063531e9,-2.77777777778e6,3.33333333333e8,-2.14928529451e6, - 2.77777777778e6,1.52734651547e9,5.55555555555e6,6.66666666667e8,2.35916180402e6, - -5.55555555555e6,-1.09779731332e8,1.56411143711e9,-2.8e6,-3333.33333333,1.0e6, - 2.83226851852e6,-30864.1975309,-5.55555555555e6,-6666.66666667,-2.0e6, - 1.63544753086e6,-15432.0987654,2.77777777778e6,-1.68e6,1.72436728395e6, - -3.33333333333e6,2.0e6,4.0e8,-2.0e6,-2.08333333333e6,1.00333333333e9, - -2.77777777778e6,3.33333333333e8,-1.0e6,2.0e8,1.0e6,2.77777777778e6,1.0675e9, - 5.55555555555e6,6.66666666667e8,-2.0e6,2.08333333333e6,-5.55555555555e6, - 1.53533333333e9,-28935.1851852,-2.08333333333e6,60879.6296296,-1.59791666667e6, - 3.37291666667e6,-28935.1851852,2.08333333333e6,2.41171296296e6,-2.08333333333e6, - 1.0e8,-2.5e6,-416666.666667,1.5e9,-833333.333333,1.25e6,5.01833333333e8, - 2.08333333333e6,1.0e8,416666.666667,5.025e8,-28935.1851852,-2.08333333333e6, - -4166.66666667,-1.25e6,3.98587962963e6,-1.59791666667e6,-8333.33333333,2.5e6, - 3.41149691358e6,-28935.1851852,2.08333333333e6,-2.355e6,2.43100308642e6, - -2.08333333333e6,1.0e8,-2.5e6,5.0e8,2.5e6,-416666.666667,1.50416666667e9, - -833333.333333,1.25e6,2.5e8,-1.25e6,-3.47222222222e6,1.33516666667e9, - 2.08333333333e6,1.0e8,-2.5e6,416666.666667,6.94444444444e6,2.16916666667e9, - -28935.1851852,-2.08333333333e6,-3.925e6,3.98587962963e6,-1.59791666667e6, - -38580.2469136,-6.94444444444e6,3.41149691358e6,-28935.1851852,2.08333333333e6, - -19290.1234568,3.47222222222e6,2.43100308642e6,-2.08333333333e6,1.0e8, - -4.16666666667e6,2.5e6,-416666.666667,1.50416666667e9,-833333.333333, - -3.47222222222e6,4.16666666667e8,-1.25e6,3.47222222222e6,1.33516666667e9, - 2.08333333333e6,1.0e8,6.94444444445e6,8.33333333333e8,416666.666667, - -6.94444444445e6,2.16916666667e9,-3830.95098171,1.14928529451e6,-275828.470683, - -28935.1851852,-2.08333333333e6,-4166.66666667,1.25e6,64710.5806113, - -131963.213599,-517922.131816,-2.29857058902e6,-1.59791666667e6,-8333.33333333, - -2.5e6,3.50487988027e6,-517922.131816,-2.16567078453e6,551656.941366, - -28935.1851852,2.08333333333e6,-2.355e6,517922.131816,4.57738374749e6, - 2.29857058902e6,-551656.941367,4.8619365099e8,-2.08333333333e6,1.0e8,2.5e6, - 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 -@test isvalid(A) - -B = A * ones(size(A,2)) -chma = cholfact(A) # LDL' form -@test isvalid(chma) -x = chma\B -@test_approx_eq x.mat ones(size(x)) - -chma = cholfact(A,true) # LL' form -@test isvalid(chma) -x = chma\B -@test_approx_eq x.mat ones(size(x)) - -#lp_afiro example -afiro = CholmodSparse(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, - 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]), - [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, - -1.0,-1.06,1.0,0.301,1.0,-1.0,1.0,-1.0,1.0,1.0,-1.0,-1.06,1.0,0.301,-1.0, - -1.06,1.0,0.313,-1.0,-0.96,1.0,0.313,-1.0,-0.86,1.0,0.326,-1.0,2.364,-1.0, - 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) -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 - -# explicit zeros -a = SparseMatrixCSC(2,2,[1,3,5],[1,2,1,2],[1.0,0.0,0.0,1.0]) -@test_approx_eq lufact(a)\[2.0,3.0] [2.0,3.0] -@test_approx_eq cholfact(a)\[2.0,3.0] [2.0,3.0]