From 3b3976f93c77f343cbc317d0b351d6ca4002a686 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 19 Jan 2015 09:14:55 -0500 Subject: [PATCH 01/10] Make support for CholModSparse more complete. --- base/sparse/cholmod.jl | 19 +++++++++++++++++-- base/sparse/sparsematrix.jl | 6 +++--- 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index ca07b0dbc6118..bcab428f36da4 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -286,13 +286,28 @@ type c_CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes} packed::Cint end -type CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes} +type CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes} <: AbstractSparseMatrix{Tv,Ti} c::c_CholmodSparse{Tv,Ti} colptr0::Vector{Ti} rowval0::Vector{Ti} nzval::Vector{Tv} end +import Base: getindex, size + +size(A::CholmodSparse, i::Integer) = i > 0 ? (i == 1 ? A.c.m : (i == 2 ? A.c.n : 1)) : error("") +size(A::CholmodSparse) = (size(A, 1), size(A, 2)) + +getindex(A::CholmodSparse, i::Integer) = getindex(A, ind2sub(size(A),i)...) +function getindex{T}(A::CholmodSparse{T}, i0::Integer, i1::Integer) + if !(1 <= i0 <= A.c.m && 1 <= i1 <= A.c.n); throw(BoundsError()); end + r1 = int(A.colptr0[i1]);println("r1: $r1") + r2 = int(A.colptr0[i1 + 1] - 1);println("r2: $r2") + (r1 > r2) && return zero(T) + r1 = int(searchsortedfirst(A.rowval0, i0 - 1, r1, r2, Base.Order.Forward));println("r1: $r1") + ((r1 > r2) || (A.rowval0[r1 + 1] + 1 != i0)) ? zero(T) : A.nzval[r1 + 1] +end + type c_CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes} m::Int n::Int @@ -433,7 +448,7 @@ function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(colpt::Vector{Ti}, int32(stype), ityp(Ti), xtyp(Tv), dtyp(Tv), CHOLMOD_FALSE, CHOLMOD_TRUE), - colpt0, rowval, nzval) + colpt0, rowval0, nzval) @isok isvalid(cs) diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index 8deada94d0e92..1f68cb9f5b788 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -881,10 +881,10 @@ getindex(A::SparseMatrixCSC, I::(Integer,Integer)) = getindex(A, I[1], I[2]) function getindex{T}(A::SparseMatrixCSC{T}, i0::Integer, i1::Integer) if !(1 <= i0 <= A.m && 1 <= i1 <= A.n); throw(BoundsError()); end - r1 = int(A.colptr[i1]) - r2 = int(A.colptr[i1+1]-1) + r1 = int(A.colptr[i1]);println("r1: $r1") + r2 = int(A.colptr[i1+1]-1);println("r2: $r2") (r1 > r2) && return zero(T) - r1 = searchsortedfirst(A.rowval, i0, r1, r2, Forward) + r1 = searchsortedfirst(A.rowval, i0, r1, r2, Forward);println("r1: $r1") ((r1 > r2) || (A.rowval[r1] != i0)) ? zero(T) : A.nzval[r1] end From 2ecb6d5c58e8f9007a691881fccf2d14ff9dce6b Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Mon, 19 Jan 2015 16:04:23 -0500 Subject: [PATCH 02/10] Fix \ polyalgorithm Split cholfact and ldltfact functionality Add cholmod_symmetry wrapper. Add issym and ishermitian methods. SuiteSparse version bumped to fix bug in cholmod_symmetry Name cleanup: - cholmod and CHOLMOD avoided within module - use same names as in the C library for the wrapper Use same argument order as the C library Argument types in ccall checked Remove scale methods for Sparse Reorganize tests --- base/linalg/diagonal.jl | 21 +- base/linalg/generic.jl | 2 +- base/linalg/symmetric.jl | 2 +- base/sparse.jl | 19 +- base/sparse/cholmod.jl | 1374 +++++++++++++++------------- base/sparse/cholmod_h.jl | 73 +- base/sparse/linalg.jl | 35 + base/sparse/sparsematrix.jl | 8 +- test/linalg/cholmod.jl | 38 - test/runtests.jl | 10 +- test/sparse/cholmod.jl | 162 ++++ test/{ => sparse}/sparse.jl | 6 + test/{linalg => sparse}/umfpack.jl | 0 test/suitesparse.jl | 127 --- 14 files changed, 1018 insertions(+), 859 deletions(-) delete mode 100644 test/linalg/cholmod.jl create mode 100644 test/sparse/cholmod.jl rename test/{ => sparse}/sparse.jl (99%) rename test/{linalg => sparse}/umfpack.jl (100%) delete mode 100644 test/suitesparse.jl 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..c0567982604d8 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -45,7 +45,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 bcab428f36da4..9dba3dc50a453 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -1,75 +1,65 @@ 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, ctranspose, eltype, getindex, hcat, isvalid, show, size, + sort!, transpose, vcat 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! + full, ishermitian, isposdef!, issym, ldltfact, logdet, norm, scale, scale! + +import Base.SparseMatrix: findnz, sparse -importall Base.SparseMatrix -import Base.SparseMatrix: increment, increment!, decrement, decrement! +export + Dense, + Factor, + Sparse, + Triplet, + etree, + sparse! + +using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, increment!, indtype, decrement, decrement! include("cholmod_h.jl") macro isok(A) - :($A==CHOLMOD_TRUE || throw(CholmodException())) + :($A == TRUE || throw(CHOLMODException(""))) end -const chm_ver = Array(Cint, 3) +const cholmod_ver = Array(Cint, 3) if dlsym(dlopen("libcholmod"), :cholmod_version) != C_NULL - ccall((:cholmod_version, :libcholmod), Cint, (Ptr{Cint},), chm_ver) + ccall((:cholmod_version, :libcholmod), Cint, (Ptr{Cint},), cholmod_ver) else - ccall((:jl_cholmod_version, :libsuitesparse_wrapper), Cint, (Ptr{Cint},), chm_ver) + ccall((:jl_cholmod_version, :libsuitesparse_wrapper), Cint, (Ptr{Cint},), cholmod_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 +const cholmod_com_sz = ccall((:jl_cholmod_common_size,:libsuitesparse_wrapper),Int,()) +const cholmod_com = fill(0xff, cholmod_com_sz) +const cholmod_l_com = fill(0xff, cholmod_com_sz) +## cholmod_com and cholmod_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) +function common(::Type{Int32}) + if isnan(reinterpret(Float64,cholmod_com[1:8])[1]) + @isok ccall((:cholmod_start, :libcholmod), Cint, (Ptr{UInt8},), cholmod_com) end - chm_com + cholmod_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) +function common(::Type{Int64}) + if isnan(reinterpret(Float64,cholmod_l_com[1:8])[1]) + @isok ccall((:cholmod_l_start, :libcholmod), Cint, (Ptr{UInt8},), cholmod_l_com) end - chm_l_com + cholmod_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) +type CHOLMODException <: Exception + msg::AbstractString end - -typealias CHMVTypes Union(Complex128, Float64) -typealias CHMVRealTypes Union(Float64) - -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 +macro cholmod_name(nm,typ) string("cholmod_", eval(typ) == Int64 ? "l_" : "", nm) end -### A way of examining some of the fields in chm_com +### A way of examining some of the fields in cholmod_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 +type Common dbound::Float64 maxrank::Int supernodal_switch::Float64 @@ -92,29 +82,29 @@ type ChmCommon end ### These offsets should be reconfigured to be less error-prone in matches -const chm_com_offsets = Array(Int, length(ChmCommon.types)) +const cholmod_com_offsets = Array(Int, length(Common.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] + Void, (Ptr{Int},), cholmod_com_offsets) +const cholmod_final_ll_inds = (1:4) + cholmod_com_offsets[7] +const cholmod_prt_inds = (1:4) + cholmod_com_offsets[13] +const cholmod_ityp_inds = (1:4) + cholmod_com_offsets[18] ### there must be an easier way but at least this works. -function ChmCommon(aa::Array{UInt8,1}) - typs = ChmCommon.types +function Common(aa::Array{UInt8,1}) + typs = Common.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)) + args = map(i->reinterpret(typs[i], aa[cholmod_com_offsets[i] + (1:sz[i])])[1], 1:length(sz)) + eval(Expr(:call, unshift!(args, :Common), Any)) end -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)]) +function set_cholmod_prt_lev(cm::Array{UInt8}, lev::Integer) # can probably be removed + cm[(1:4) + cholmod_com_offsets[13]] = reinterpret(UInt8, [int32(lev)]) 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 +## Ptr{c_Dense}. The Dense type contains a c_Dense object and other ## fields then ensure the memory pointed to is freed when it should be and not before. -type c_CholmodDense{T<:CHMVTypes} +type c_Dense{T<:VTypes} m::Int n::Int nzmax::Int @@ -125,13 +115,13 @@ type c_CholmodDense{T<:CHMVTypes} dtype::Cint end -type CholmodDense{T<:CHMVTypes} - c::c_CholmodDense +type Dense{T<:VTypes} <: AbstractMatrix{T} + c::c_Dense mat::Matrix{T} end -if (1000chm_ver[1]+chm_ver[2]) >= 2001 # CHOLMOD version 2.1.0 or later - type c_CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} +if (1000cholmod_ver[1]+cholmod_ver[2]) >= 2001 # CHOLMOD version 2.1.0 or later + type c_Factor{Tv<:VTypes,Ti<:ITypes} n::Int minor::Int Perm::Ptr{Ti} @@ -163,8 +153,8 @@ if (1000chm_ver[1]+chm_ver[2]) >= 2001 # CHOLMOD version 2.1.0 or later dtype::Cint end - type CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} - c::c_CholmodFactor{Tv,Ti} + type Factor{Tv<:VTypes,Ti<:ITypes} + c::c_Factor{Tv,Ti} Perm::Vector{Ti} ColCount::Vector{Ti} IPerm::Vector{Ti} @@ -180,7 +170,7 @@ if (1000chm_ver[1]+chm_ver[2]) >= 2001 # CHOLMOD version 2.1.0 or later s::Vector{Ti} end - function CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodFactor{Tv,Ti}}) + function Factor{Tv<:VTypes,Ti<:ITypes}(cp::Ptr{c_Factor{Tv,Ti}}) cfp = unsafe_load(cp) Perm = pointer_to_array(cfp.Perm, (cfp.n,), true) ColCount = pointer_to_array(cfp.ColCount, (cfp.n,), true) @@ -195,13 +185,13 @@ if (1000chm_ver[1]+chm_ver[2]) >= 2001 # CHOLMOD version 2.1.0 or later 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, + cf = Factor{Tv,Ti}(cfp, Perm, ColCount, IPerm, p, i, x, nz, next, prev, super, pi, px, s) c_free(cp) cf end else - type c_CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} + type c_Factor{Tv<:VTypes,Ti<:ITypes} n::Int minor::Int Perm::Ptr{Ti} @@ -232,8 +222,8 @@ else dtype::Cint end - type CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} - c::c_CholmodFactor{Tv,Ti} + type Factor{Tv<:VTypes,Ti<:ITypes} + c::c_Factor{Tv,Ti} Perm::Vector{Ti} ColCount::Vector{Ti} p::Vector{Ti} @@ -248,7 +238,7 @@ else s::Vector{Ti} end - function CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodFactor{Tv,Ti}}) + function Factor{Tv<:VTypes,Ti<:ITypes}(cp::Ptr{c_Factor{Tv,Ti}}) cfp = unsafe_load(cp) Perm = pointer_to_array(cfp.Perm, (cfp.n,), true) ColCount = pointer_to_array(cfp.ColCount, (cfp.n,), true) @@ -262,20 +252,20 @@ else 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, + cf = Factor{Tv,Ti}(cfp, Perm, ColCount, p, i, x, nz, next, prev, super, pi, px, s) c_free(cp) cf end end -type c_CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes} +type c_Sparse{Tv<:VTypes,Ti<:ITypes} m::Csize_t n::Csize_t nzmax::Csize_t ppt::Ptr{Ti} ipt::Ptr{Ti} - nzpt::Ptr{Void} + nzpt::Ptr{Ti} xpt::Ptr{Tv} zpt::Ptr{Void} stype::Cint @@ -286,29 +276,32 @@ type c_CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes} packed::Cint end -type CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes} <: AbstractSparseMatrix{Tv,Ti} - c::c_CholmodSparse{Tv,Ti} +# Corresponds to the exact definition of cholmod_sparse_void in the library. Useful when reading files of unknown type. +type c_SparseVoid + m::Csize_t + n::Csize_t + nzmax::Csize_t + ppt::Ptr{Void} + ipt::Ptr{Void} + nzpt::Ptr{Void} + xpt::Ptr{Void} + zpt::Ptr{Void} + stype::Cint + itype::Cint + xtype::Cint + dtype::Cint + sorted::Cint + packed::Cint +end + +type Sparse{Tv<:VTypes,Ti<:ITypes} <: AbstractSparseMatrix{Tv,Ti} + c::c_Sparse{Tv,Ti} colptr0::Vector{Ti} rowval0::Vector{Ti} nzval::Vector{Tv} end -import Base: getindex, size - -size(A::CholmodSparse, i::Integer) = i > 0 ? (i == 1 ? A.c.m : (i == 2 ? A.c.n : 1)) : error("") -size(A::CholmodSparse) = (size(A, 1), size(A, 2)) - -getindex(A::CholmodSparse, i::Integer) = getindex(A, ind2sub(size(A),i)...) -function getindex{T}(A::CholmodSparse{T}, i0::Integer, i1::Integer) - if !(1 <= i0 <= A.c.m && 1 <= i1 <= A.c.n); throw(BoundsError()); end - r1 = int(A.colptr0[i1]);println("r1: $r1") - r2 = int(A.colptr0[i1 + 1] - 1);println("r2: $r2") - (r1 > r2) && return zero(T) - r1 = int(searchsortedfirst(A.rowval0, i0 - 1, r1, r2, Base.Order.Forward));println("r1: $r1") - ((r1 > r2) || (A.rowval0[r1 + 1] + 1 != i0)) ? zero(T) : A.nzval[r1 + 1] -end - -type c_CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes} +type c_Triplet{Tv<:VTypes,Ti<:ITypes} m::Int n::Int nzmax::Int @@ -323,661 +316,747 @@ type c_CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes} dtype::Cint end -type CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes} - c::c_CholmodTriplet{Tv,Ti} +type Triplet{Tv<:VTypes,Ti<:ITypes} + c::c_Triplet{Tv,Ti} i::Vector{Ti} j::Vector{Ti} x::Vector{Tv} 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 -## 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) -end +eltype{T<:VTypes}(::Type{Dense{T}}) = T +eltype{T<:VTypes}(::Type{Factor{T}}) = T +eltype{T<:VTypes}(::Type{Sparse{T}}) = T -function CholmodDense{T<:CHMVTypes}(c::Ptr{c_CholmodDense{T}}) +## The Dense constructor does not copy the contents, which is generally what you +## want as most uses of Dense objects are read-only. +function Dense{T<:VTypes}(A::VecOrMat{T}) # uses the memory from Julia + m, n = size(A,1), size(A,2) + Dense(c_Dense{T}(m, n, m*n, stride(A, 2), convert(Ptr{T}, A), C_NULL, xtyp(T), dtyp(T)), + length(size(A)) == 2 ? A : reshape(A, (m, n))) +end +function Dense{T<:VTypes}(c::Ptr{c_Dense{T}}) cp = unsafe_load(c) 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)) + val = Dense(cp, pointer_to_array(cp.xpt, (cp.m,cp.n), true)) c_free(c) val 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))) -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))) -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}, +function Sparse{Tv<:VTypes,Ti<:ITypes}(colpt::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}, m::Integer, n::Integer, stype::Signed) bb = colpt[1] - if 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 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 + else # zero based 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)]")) - 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)]")) + # end it = ityp(Ti) - cs = CholmodSparse(c_CholmodSparse{Tv,Ti}(m, n, int(nz), convert(Ptr{Ti}, colpt0), + cs = Sparse(c_Sparse{Tv,Ti}(m, n, int(nz), convert(Ptr{Ti}, colpt0), convert(Ptr{Ti}, rowval0), C_NULL, convert(Ptr{Tv}, nzval), C_NULL, int32(stype), ityp(Ti), xtyp(Tv), dtyp(Tv), - CHOLMOD_FALSE, CHOLMOD_TRUE), + TRUE, TRUE), colpt0, rowval0, nzval) @isok isvalid(cs) - cs = sort!(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) +function Sparse{Tv<:VTypes,Ti<:ITypes}(A::SparseMatrixCSC{Tv,Ti}, stype::Signed) + Sparse(A.colptr, A.rowval, A.nzval, size(A,1), size(A,2), stype) end -function CholmodSparse(A::SparseMatrixCSC) - stype = ishermitian(A) ? 1 : 0 - CholmodSparse(stype > 0 ? triu(A) : A, stype) +function Sparse(A::SparseMatrixCSC) + AC = Sparse(A, 0) + i = symmetry(AC, 0)[1] # Check for Hermitianity + if (isreal(A) && (i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG)) || # real case + (i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG) # complex case + AC.c.stype = -1 + end + AC end -function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(cp::Ptr{c_CholmodSparse{Tv,Ti}}) +Sparse{Ti<:ITypes}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) +Sparse{Ti<:ITypes}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) + + +# Constructs a Sparse object from a pointer to a cholmod_sparse_struct created by CHOLMOD. Julia takes ownership of the memory allocated by CHOLMOD and the pointer can therefore be freed without memory leek. +function Sparse{Tv<:VTypes,Ti<:ITypes}(cp::Ptr{c_Sparse{Tv,Ti}}) csp = unsafe_load(cp) + + # Take control of the memory allocated by CHOLMOD colptr0 = pointer_to_array(csp.ppt, (csp.n + 1,), true) - nnz = int(colptr0[end]) - cms = CholmodSparse{Tv,Ti}(csp, colptr0, - pointer_to_array(csp.ipt, (nnz,), true), - pointer_to_array(csp.xpt, (nnz,), true)) + rowval0 = pointer_to_array(csp.ipt, (csp.nzmax,), true) + nzval = pointer_to_array(csp.xpt, (csp.nzmax,), true) + + cms = Sparse{Tv,Ti}(csp, colptr0, rowval0, nzval) + + # Free memory + csp.packed == 1 && c_free(csp.nzpt) # Julia's Sparse is not using unpacked storage c_free(cp) + cms end -CholmodSparse{Tv<:CHMVTypes}(D::CholmodDense{Tv}) = CholmodSparse(D,1) # default Ti is Int +# Useful when reading in files, but not type stable +function Sparse(p::Ptr{c_SparseVoid}) + s = unsafe_load(p) + + # Check integer type + if s.itype == INT + Ti = Cint + elseif s.itype == INTLONG + throw(CHOLMODException("the value of itype was $s.itype. This combination of integer types shouldn't happen. Please submit a bug report.")) + elseif s.itype == LONG + Ti = SuiteSparse_long + else + throw(CHOLMODException("illegal value of itype: $s.itype")) + end + + # Check for double or single precision + if s.dtype == DOUBLE + Tv = Float64 + elseif s.dtype == SINGLE + Tv = Float32 + else + throw(CHOLMODException("illegal value of dtype: $s.dtype")) + end + + # Check for real or complex + if s.xtype == COMPLEX + Tv = Complex{Tv} + elseif s.xtype != REAL + throw(CHOLMODException("illegal value of xtype: $s.xtype")) + end + + # Take control of the memory allocated by CHOLMOD + colptr0 = pointer_to_array(convert(Ptr{Ti}, s.ppt), (s.n + 1,), true) + rowval0 = pointer_to_array(convert(Ptr{Ti}, s.ipt), (s.nzmax,), true) + nzval = pointer_to_array(convert(Ptr{Tv}, s.xpt), (s.nzmax,), true) + + r = Sparse{Tv,Ti}(c_Sparse{Tv,Ti}(s.m, + s.n, + s.nzmax, + pointer(colptr0), + pointer(rowval0), + C_NULL, + pointer(nzval), + C_NULL, + s.stype, + s.itype, + s.xtype, + s.dtype, + s.sorted, + s.packed), + colptr0, rowval0, nzval) + + # Free memory + s.packed == 1 && c_free(s.nzpt) # Julia's Sparse is not using unpacked storage + c_free(p) + + r +end + +Sparse{Tv<:VTypes}(D::Dense{Tv}) = Sparse(D,1) # default Ti is Int -function CholmodTriplet{Tv<:CHMVTypes,Ti<:CHMITypes}(tp::Ptr{c_CholmodTriplet{Tv,Ti}}) + +# Constructs a Triplet object from a pointer to a cholmod_triplet_struct created by CHOLMOD. Julia takes ownership of the memory allocated by CHOLMOD and the pointer can therefore be freed without memory leek. +function Triplet{Tv<:VTypes,Ti<:ITypes}(tp::Ptr{c_Triplet{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) + ct = Triplet{Tv,Ti}(ctp, i, j, x) + + csp.packed == 1 && c_free(csp.nzpt) # Julia's Sparse is not using unpacked storage c_free(tp) ct 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) +# Dense wrappers +### cholmod_core_h ### +function zeros{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) + Dense(ccall((:cholmod_zeros, :libcholmod), Ptr{c_Dense{T}}, + (Csize_t, Csize_t, Cint, Ptr{UInt8}), + m, n, xtyp(T), common(Int32))) +end +zeros(m::Integer, n::Integer) = zeros(m, n, Float64) + +function ones{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) + Dense(ccall((:cholmod_ones, :libcholmod), Ptr{c_Dense{T}}, + (Csize_t, Csize_t, Cint, Ptr{UInt8}), + m, n, xtyp(T), common(Int32))) +end +ones(m::Integer, n::Integer) = ones(m, n, Float64) + +function eye{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) + Dense(ccall((:cholmod_eye, :libcholmod), Ptr{c_Dense{T}}, + (Csize_t, Csize_t, Cint, Ptr{UInt8}), + m, n, xtyp(T), common(Int32))) +end +eye(m::Integer, n::Integer) = eye(m, n, Float64) +eye(n::Integer) = eye(n, n, Float64) + +function copy_dense{Tv<:VTypes}(A::Dense{Tv}) + Dense(ccall((:cholmod_copy_dense, :libcholmod), Ptr{c_Dense{Tv}}, + (Ptr{c_Dense{Tv}}, Ptr{UInt8}), + &A.c, common(Cint))) +end + +### cholmod_matrixops.h ### +function norm_dense{Tv<:VTypes}(D::Dense{Tv}, p::Integer) + ccall((:cholmod_norm_dense, :libcholmod), Cdouble, + (Ptr{c_Dense{Tv}}, Cint, Ptr{UInt8}), + &D.c, p, common(Int32)) +end + +### cholmod_check.h ### +function check_dense{T<:VTypes}(A::Dense{T}) + bool(ccall((:cholmod_check_dense, :libcholmod), Cint, + (Ptr{c_Dense{T}}, Ptr{UInt8}), + &A.c, common(Cint))) 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))) + # This method is probably not necessary as memory allocated by CHOLMOD is usually overtaken by Julia when cosntructing Julia instances from pointers. + function free_sparse{Tv<:VTypes}(ptr::Ptr{c_Sparse{Tv,$Ti}}) + cm = common($Ti) + @isok ccall((@cholmod_name "free_sparse" $Ti + ,:libcholmod), Cint, + (Ptr{Ptr{c_Sparse{Tv,$Ti}}}, Ptr{UInt8}), + &ptr, cm) 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 + function (*){Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}) + Sparse(ccall((@cholmod_name "ssmult" $Ti + , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Cint, Ptr{UInt8}), + &A.c, &B.c, 0, true, true, common($Ti))) + end - cs = CholmodSparse(aa[2]) + function transpose{Tv<:VTypes}(A::Sparse{Tv,$Ti}, values::Integer) + cm = common($Ti) + Sparse(ccall((@cholmod_name "transpose" $Ti + ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + &A.c, values, cm)) + 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 ssmult{Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, stype::Integer, values::Integer, sorted::Integer) + cm = common($Ti) + Sparse(ccall((@cholmod_name "ssmult" $Ti + ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Cint, Ptr{UInt8}), + &A.c, &B.c, stype, values, sorted, cm)) + end - return sort!(cs) + function aat{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, fset::Vector{$Ti}, mode::Integer) + cm = common($Ti) + Sparse(ccall((@cholmod_name "aat" $Ti + , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{$Ti}, Csize_t, Cint, Ptr{UInt8}), + &A.c, fset, length(fset), mode, cm)) 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 copy{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, stype::Integer, mode::Integer) + cm = common($Ti) + Sparse(ccall((@cholmod_name "copy" $Ti + , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), + &A.c, 0, 1, cm)) 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 sparse_to_dense{Tv<:VTypes}(A::Sparse{Tv,$Ti}) + Dense(ccall((@cholmod_name "sparse_to_dense" $Ti + ,:libcholmod), Ptr{c_Dense{Tv}}, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), + &A.c, common($Ti))) 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 dense_to_sparse{Tv<:VTypes}(D::Dense{Tv},i::$Ti) + Sparse(ccall((@cholmod_name "dense_to_sparse" $Ti + ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Dense{Tv}}, Cint, Ptr{UInt8}), + &D.c, true, common($Ti))) 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 factor_to_sparse{Tv<:VTypes}(L::Factor{Tv,$Ti}) + Sparse(ccall((@cholmod_name "factor_to_sparse" $Ti + ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), + &L.c, common($Ti))) 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 copy_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) + Factor(ccall((@cholmod_name "copy_factor" $Ti + ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, + (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), + &L.c, common($Ti))) 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 change_factor{Tv<:VTypes}(to_xtype::Integer, to_ll::Bool, to_super::Bool, to_packed::Bool, to_monotonic::Bool, L::Factor{Tv,$Ti}) + @isok ccall((@cholmod_name "change_factor" $Ti + ,:libcholmod), Cint, + (Cint, Cint, Cint, Cint, Cint, Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), + to_xtype, to_ll, to_super, to_packed, to_monotonic, L, cm) + L 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 triplet_to_sparse{Tv<:VTypes,Ti<:$Ti}(T::Triplet{Tv,Ti}) + Sparse(ccall((@cholmod_name "triplet_to_sparse" $Ti + ,:libcholmod), Ptr{c_Sparse{Tv,Ti}}, + (Ptr{c_Triplet{Tv,Ti}}, Ptr{UInt8}), + &T.c, 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 sparse_to_triplet{Tv<:VTypes,Ti<:$Ti}(A::Sparse{Tv,Ti}) + Triplet(ccall((@cholmod_name "sparse_to_triplet" $Ti + ,:libcholmod), Ptr{c_Triplet{Tv,Ti}}, + (Ptr{c_Sparse{Tv,Ti}}, Ptr{UInt8}), + &A.c, 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 check_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}) + bool(ccall((@cholmod_name "check_sparse" $Ti + ,:libcholmod), Cint, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), + &A.c, common($Ti))) 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 check_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) + bool(ccall((@cholmod_name "check_factor" $Ti + ,:libcholmod), Cint, + (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), + &L.c, common($Ti))) end - function 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 check_triplet{Tv<:VTypes}(T::Triplet{Tv,$Ti}) + bool(ccall((@cholmod_name "check_triplet" $Ti + ,:libcholmod), Cint, + (Ptr{c_Triplet{Tv,$Ti}}, Ptr{UInt8}), + &T.c, common($Ti))) 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 analyze{Tv<:VTypes}(A::Sparse{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) + Factor(ccall((@cholmod_name "analyze" $Ti + ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), + &A.c, C)) end - function 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 factorize{Tv<:VTypes}(A::Sparse{Tv,$Ti}, F::Factor{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) + @isok ccall((@cholmod_name "factorize" $Ti + ,:libcholmod), Cint, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), + &A.c, &F.c, C) + F end - function 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 factorize_p{Tv<:VTypes}(A::Sparse{Tv,$Ti}, β::Tv, fset::Vector{$Ti}, F::Factor{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) + @isok ccall((@cholmod_name "factorize_p" $Ti + ,:libcholmod), Cint, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{Cdouble}, Ptr{$Ti}, Csize_t, + Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), + &A.c, &β, fset, length(fset), + &F.c, common($Ti)) + F end - 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)) + + ### cholmod_core.h ### + function pack_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) + @isok ccall((@cholmod_name "pack_factor" $Ti + ,:libcholmod), Cint, + (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), + &L.c, common($Ti)) L 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 nnz{Tv<:VTypes}(A::Sparse{Tv,$Ti}) + ccall((@cholmod_name "nnz" $Ti + ,:libcholmod), Int, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), + &A.c, common($Ti)) 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 + + function speye{Tv<:VTypes}(m::Integer, n::Integer, ::Type{Tv}) + 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))) 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 spzeros{Tv<:VTypes}(m::Integer, n::Integer, nzmax::Integer, ::Type{Tv}) + 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))) 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 copy_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) + Factor(ccall((@cholmod_name "copy_factor" $Ti + ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, + (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), + &L.c, common($Ti))) 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)) - A + function copy_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}) + Sparse(ccall((@cholmod_name "copy_sparse" $Ti + ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), + &A.c, common($Ti))) 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")) - 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)) - Y + function copy_triplet{Tv<:VTypes}(T::Triplet{Tv,$Ti}) + Triplet(ccall((@cholmod_name "copy_triplet" $Ti + ,:libcholmod), Ptr{c_Triplet{Tv,$Ti}}, + (Ptr{c_Triplet{Tv,$Ti}}, Ptr{UInt8}), + &T.c, common($Ti))) 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))) + function copy{Tv<:VTypes}(A::Sparse{Tv,$Ti}, stype::Integer, mode::Integer) + Sparse(ccall((@cholmod_name "copy" $Ti + ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), + &A.c, stype, mode, common($Ti))) 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))) + + ### cholmod_check.h ### + function print_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}, name::ASCIIString) + @isok ccall((@cholmod_name "print_sparse" $Ti + ,:libcholmod), Cint, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), + &A.c, name, common($Ti)) + nothing 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))) + function print_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}, name::ASCIIString) + @isok ccall((@cholmod_name "print_factor" $Ti + ,:libcholmod), Cint, + (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), + &L.c, name, common($Ti)) + nothing 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))) + + ### cholmod_matrixops.h ### + function norm_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}, norm::Integer) + norm == 0 || norm == 1 || throw(ArgumentError("norm argument must be either 0 or 1")) + ccall((@cholmod_name "norm_sparse" $Ti + , :libcholmod), Cdouble, + (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + &A.c, norm, common($Ti)) 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))) + + function horzcat{Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) + ccall((@cholmod_name "horzcat" $Ti + , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + &A.c, &B.c, values, common($Ti)) 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))) + + function scale{Tv<:VTypes}(S::Dense{Tv}, scale::Integer, A::Sparse{Tv,$Ti}) + @isok ccall((@cholmod_name "scale" $Ti + ,:libcholmod), Cint, + (Ptr{c_Dense{Tv}}, Cint, Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), + &S.c, scale, &A.c, common($Ti)) + A 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 + + 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 + @isok ccall((@cholmod_name "sdmult" $Ti + ,:libcholmod), Cint, + (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, + Ptr{c_Dense{Tv}}, Ptr{c_Dense{Tv}}, Ptr{UInt8}), + &A.c, transpose, &α, &β, + &X.c, &Y.c, common($Ti)) + Y 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)) + + function vertcat{Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) + ccall((@cholmod_name "vertcat" $Ti + , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + &A.c, &B.c, values, common($Ti)) 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 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.c, option, xmatched, pmatched, + nzoffdiag, nzdiag, common($Ti)) + rv, xmatched[1], pmatched[1], nzoffdiag[1], nzdiag[1] 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)) + + # cholmod_cholesky.h + function etree{Tv<:VTypes}(A::Sparse{Tv,$Ti}) + Parent = Array($Ti, size(A,2)) + @isok ccall((@cholmod_name "etree" $Ti + ,:libcholmod), Cint, + (Ptr{c_Sparse{Tv,$Ti}}, Ptr{$Ti}, Ptr{UInt8}), + &A.c, Parent, common($Ti)) + return Parent end - function solve{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti}, - B::CholmodDense{Tv}, typ::Integer) + + function solve{Tv<:VTypes}(sys::Integer, L::Factor{Tv,$Ti}, B::Dense{Tv}) size(L,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(L,1)) rows, but RHS has $(size(B,1)) rows.")) - 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))) + Dense(ccall((@cholmod_name "solve" $Ti + ,:libcholmod), Ptr{c_Dense{Tv}}, + (Cint, Ptr{c_Factor{Tv,$Ti}}, Ptr{c_Dense{Tv}}, Ptr{UInt8}), + sys, &L.c, &B.c, common($Ti))) end - function solve{Tv<:CHMVTypes}(L::CholmodFactor{Tv,$Ti}, - B::CholmodSparse{Tv,$Ti}, - typ::Integer) + + function spsolve{Tv<:VTypes}(sys::Integer, L::Factor{Tv,$Ti}, B::Sparse{Tv,$Ti}) size(L,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(L,1)) rows, but RHS has $(size(B,1)) rows.")) - 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))) - 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 - 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))) - 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))) - 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)) + Sparse(ccall((@cholmod_name "spsolve" $Ti + ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + (Cint, Ptr{c_Factor{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), + sys, &L.c, &B.c, common($Ti))) 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)) +# Autodetects the types +function read_sparse(file::CFILE) + Sparse(ccall((:cholmod_read_sparse, :libcholmod), Ptr{c_SparseVoid}, + (Ptr{Void}, Ptr{UInt8}), + file.ptr, common(Int32))) +end + +######################### +# High level interfaces # +######################### -function Ac_mul_B{Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::CholmodDense{Tv}) - chm_sdmult(A,true,1.,0.,B) +# Convertion +Dense(A::Sparse) = sparse_to_dense(A) + +Sparse(A::Dense) = dense_to_sparse(A) +Sparse(A::Triplet) = triplet_to_sparse(A) +function Sparse(L::Factor) + if bool(L.c.is_ll) + return factor_to_sparse(L) + end + cm = common($Ti) + Lcll = copy_factor(L) + change_factor(L.c.xtype, true, L.c.is_super == 1, true, true, Lcll, L) + return factor_to_sparse(Lcll) end -function Ac_mul_B{Tv<:CHMVTypes}(A::CholmodSparse{Tv},B::VecOrMat{Tv}) - chm_sdmult(A,true,1.,0.,CholmodDense(B)) +function Sparse(filename::ASCIIString) + f = open(filename) + A = read_sparse(CFILE(f)) + close(f) + A end +Triplet(A::Sparse) = sparse_to_triplet(A) -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)) +show(io::IO, L::Factor) = print_factor(L, "") + +isvalid(A::Dense) = check_dense(A) +isvalid(A::Sparse) = check_sparse(A) +isvalid(A::Triplet) = check_triplet(A) +isvalid(A::Factor) = check_factor(A) + +copy(A::Dense) = copy_dense(A) +copy(A::Sparse) = copy_sparse(A) +copy(A::Triplet) = copy_triplet(A) +copy(A::Factor) = copy_factor(A) + +size(B::Dense) = size(B.mat) +size(B::Dense, d::Integer) = size(B.mat, d) +size(A::Sparse) = (int(A.c.m), int(A.c.n)) +size(A::Sparse, d::Integer) = d > 0 ? (d == 1 ? int(A.c.m) : (d == 2 ? int(A.c.n) : 1)) : BoundsError() +size(L::Factor) = (n = int(L.c.n); (n,n)) +size(L::Factor, d::Integer) = d < 1 ? throw(BoundsError("Dimension $d out of range")) : (d <= 2 ? int(L.c.n) : 1) + +getindex(A::Dense, i::Integer) = A.mat[i] +getindex(A::Dense, i::Integer, j::Integer) = A.mat[i, j] +getindex(A::Sparse, i::Integer) = getindex(A, ind2sub(size(A),i)...) +function getindex{T}(A::Sparse{T}, i0::Integer, i1::Integer) + if !(1 <= i0 <= A.c.m && 1 <= i1 <= A.c.n); throw(BoundsError()); end + r1 = int(A.colptr0[i1] + 1) + r2 = int(A.colptr0[i1 + 1]) + (r1 > r2) && return zero(T) + r1 = int(searchsortedfirst(A.rowval0, i0 - 1, r1, r2, Base.Order.Forward)) + ((r1 > r2) || (A.rowval0[r1] + 1 != i0)) ? zero(T) : A.nzval[r1] end -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) +full(A::Sparse) = Dense(A).mat +full(A::Dense) = A.mat +function sparse(A::Sparse) + 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 Ac_ldiv_B{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::SparseMatrixCSC{Tv,Ti}) - sparse!(solve(L,CholmodSparse(B),CHOLMOD_A)) +function sparse!(A::Sparse) + SparseMatrixCSC(A.c.m, A.c.n, increment!(A.colptr0), increment!(A.rowval0), A.nzval) end +sparse(L::Factor) = sparse!(Sparse(L)) +sparse(D::Dense) = sparse!(Sparse(D)) +sparse(T::Triplet) = sparse!(Sparse(T)) +## Multiplication -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)) -end -function cholfact!{T<:CHMVTypes}(L::CholmodFactor{T},A::SparseMatrixCSC{T},beta::Number) - cholfact!(L,CholmodSparse(A),convert(T,beta)) -end -cholfact!{T<:CHMVTypes}(L::CholmodFactor{T},A::SparseMatrixCSC{T}) = cholfact!(L,CholmodSparse(A)) +(*)(A::Sparse, B::Dense) = sdmult(A, false, 1., 0., B, zeros(size(A, 1), size(B, 2))) +(*)(A::Sparse, B::VecOrMat) = (*)(A, Dense(B)) -chm_analyze(A::SparseMatrixCSC) = chm_analyze(CholmodSparse(A)) +function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti}) + cm = common($Ti) -chm_print(A::CholmodSparse, lev::Integer) = chm_print(A, lev, "") -chm_print(A::CholmodFactor, lev::Integer) = chm_print(L, lev, "") + if !is(A,B) + aa1 = transpose(B, 1) + ## result of ssmult will have stype==0, contain numerical values and be sorted + return ssmult(A, aa1, 0, true, true) + end -function chm_scale!{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T},typ::Integer) - chm_scale!(A,CholmodDense(b),typ) - A + ## 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 + aa1 = copy(A, 0, 1) + return aat(aa1, Ti[], 1) + else + return aat(A, Ti[], 1) + end end -chm_scale{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T},typ::Integer) = chm_scale!(copy(A),b,typ) -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 +function Ac_mul_B{Tv<:VRealTypes}(A::Sparse{Tv}, B::Sparse{Tv}) + aa1 = transpose(A, 1) + if is(A,B) + return A_mul_Bc(aa1, aa1) + end + ## result of ssmult will have stype==0, contain numerical values and be sorted + return ssmult(aa1, B, 0, true, true, cm) +end -chm_spzeros(m::Integer,n::Integer,nzmax::Integer) = chm_spzeros(m,n,nzmax,1.) +A_mul_Bt{Tv<:VRealTypes}(A::Sparse{Tv}, B::Sparse{Tv}) = A_mul_Bc(A, B) # in the unlikely event of writing A*B.' instead of A*B' -function scale!{T<:CHMVTypes}(b::Vector{T}, A::CholmodSparse{T}) - chm_scale!(A,CholmodDense(b),CHOLMOD_ROW) - 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) - A -end -scale{T<:CHMVTypes}(A::CholmodSparse{T},b::Vector{T}) = scale!(copy(A), b) +At_mul_B{Tv<:VRealTypes}(A::Sparse{Tv}, B::Sparse{Tv}) = Ac_mul_B(A,B) # in the unlikely event of writing A.'*B instead of A'*B -norm(A::CholmodSparse) = norm(A,1) +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)) -show(io::IO,L::CholmodFactor) = chm_print(L,int32(4),"") -show(io::IO,A::CholmodSparse) = chm_print(A,int32(4),"") -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) +## Factorization methods -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)) +analyze(A::SparseMatrixCSC) = analyze(Sparse(A)) + +function cholfact(A::Sparse) + cm = common(indtype(A)) + ## may need to change final_asis as well as final_ll + cm[cholmod_final_ll_inds] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + F = analyze(A, cm) + factorize(A, F, cm) + F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) + return F end -function solve{Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti}) - solve(L,B,CHOLMOD_A) + +function ldltfact(A::Sparse) + cm = common(indtype(A)) + ## may need to change final_asis as well as final_ll + cm[cholmod_final_ll_inds] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt + F = analyze(A, cm) + factorize(A, F, cm) + return F end -function (\){Tv<:CHMVTypes,Ti<:CHMITypes}(L::CholmodFactor{Tv,Ti},B::CholmodSparse{Tv,Ti}) - solve(L,B,CHOLMOD_A) + +function cholfact{Tv<:VTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, β::Tv) + cm = common(Ti) + ## may need to change final_asis as well as final_ll + cm[cholmod_final_ll_inds] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + F = analyze(A, cm) + factorize_p(A, β, Ti[], F, cm) + F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) + return F end -function solve{T<:CHMVTypes}(L::CholmodFactor{T},B::Vector{T},typ::Integer) - solve(L,CholmodDense!(B),typ).mat[:] + +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[cholmod_final_ll_inds] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt + F = analyze(A, cm) + factorize_p(A, β, Ti[], F, cm) + F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) + return F end -function solve{T<:CHMVTypes}(L::CholmodFactor{T},B::Matrix{T},typ::Integer) - solve(L,CholmodDense!(B),typ).mat + +cholfact!(L::Factor, A::Sparse) = factorize(A, L) +cholfact!{Tv<:VTypes,Ti<:ITypes}(L::Factor{Tv,Ti}, A::Sparse{Tv,Ti}, β::Tv) = factorize_p(A, β, Ti[], L) + +cholfact{T<:VTypes}(A::Sparse{T}, β::Number) = cholfact(A, convert(T, β)) +cholfact(A::SparseMatrixCSC) = cholfact(Sparse(A)) +cholfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = cholfact(Sparse(A)) +cholfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = cholfact(Sparse(A)) +cholfact!{T<:VTypes}(L::Factor{T}, A::Sparse{T}, β::Number) = cholfact!(L, A, convert(T, β)) +cholfact!{T<:VTypes}(L::Factor{T}, A::SparseMatrixCSC{T}, β::Number) = cholfact!(L, Sparse(A), convert(T, β)) +cholfact!{T<:VTypes}(L::Factor{T}, A::SparseMatrixCSC{T}) = cholfact!(L, Sparse(A)) + +ldltfact(A::SparseMatrixCSC) = ldltfact(Sparse(A)) +ldltfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = ldltfact(Sparse(A)) +ldltfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = ldltfact(Sparse(A)) + +## Solvers + +(\)(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B) +(\)(L::Factor, b::Vector) = reshape(solve(CHOLMOD_A, L, Dense(b)).mat, length(b)) +(\)(L::Factor, B::Matrix) = solve(CHOLMOD_A, L, Dense(B)).mat +(\)(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B) +(\)(L::Factor, B::SparseMatrixCSC) = sparse!(spsolve(CHOLMOD_A, L, Sparse(B))) + +Ac_ldiv_B(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B) +Ac_ldiv_B(L::Factor, B::VecOrMat) = solve(CHOLMOD_A, L, Dense(B)).mat +Ac_ldiv_B(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B) +Ac_ldiv_B(L::Factor, B::SparseMatrixCSC) = sparse!(CHOLMOD_A, spsolve(L, Sparse(B))) + +speye(m::Integer, n::Integer) = speye(m, n, Float64) # default element type is Float64 +speye(n::Integer) = speye(n, n, Float64) # default shape is square + +spzeros(m::Integer, n::Integer, nzmax::Integer) = spzeros(m, n, nzmax, 1.) + +function norm(A::Sparse, p = 1) + if p == 1 + return norm_sparse(A, 1) + elseif p == Inf + return norm_sparse(A, 0) + else + throw(ArgumentError("only 1 and inf norms are supported")) + end end -solve{T<:CHMVTypes}(L::CholmodFactor{T},B::CholmodDense{T}) = solve(L,B,CHOLMOD_A) -function findnz{Tv,Ti}(A::CholmodSparse{Tv,Ti}) +## Julia methods (non-CHOLMOD) +function findnz{Tv,Ti}(A::Sparse{Tv,Ti}) jj = similar(A.rowval0) # expand A.colptr0 to a vector of indices for j in 1:A.c.n, k in (A.colptr0[j]+1):A.colptr0[j+1] jj[k] = j @@ -999,9 +1078,9 @@ function findnz{Tv,Ti}(A::CholmodSparse{Tv,Ti}) (increment!(A.rowval0[ind]), jj[ind], A.nzval[ind]) end -findnz(L::CholmodFactor) = findnz(CholmodSparse(L)) +findnz(L::Factor) = findnz(Sparse(L)) -function diag{Tv}(A::CholmodSparse{Tv}) +function diag{Tv}(A::Sparse{Tv}) minmn = minimum(size(A)) res = zeros(Tv,minmn) cp0 = A.colptr0 @@ -1015,7 +1094,7 @@ function diag{Tv}(A::CholmodSparse{Tv}) res end -function diag{Tv}(L::CholmodFactor{Tv}) +function diag{Tv}(L::Factor{Tv}) res = zeros(Tv,L.c.n) xv = L.x if bool(L.c.is_super) @@ -1045,37 +1124,50 @@ function diag{Tv}(L::CholmodFactor{Tv}) res end -function logdet{Tv,Ti}(L::CholmodFactor{Tv,Ti},sub) +function logdet{Tv,Ti}(L::Factor{Tv,Ti},sub) 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) -end -sparse(L::CholmodFactor) = sparse!(CholmodSparse(L)) -sparse(D::CholmodDense) = sparse!(CholmodSparse(D)) -sparse(T::CholmodTriplet) = sparse!(CholmodSparse(T)) +logdet(L::Factor) = logdet(L, 1:L.c.n) +det(L::Factor) = exp(logdet(L)) -isposdef!{Tv<:CHMVTypes,Ti}(A::SparseMatrixCSC{Tv,Ti}) = ishermitian(A) && cholfact(A).c.minor == size(A,1) +isposdef!{Tv<:VTypes,Ti}(A::SparseMatrixCSC{Tv,Ti}) = ishermitian(A) && cholfact(A).c.minor == size(A,1) -end #module +function issym(A::Sparse) + if isreal(A) + if A.c.stype != 0 + return true + else + i = symmetry(A, 0)[1] + return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG + end + else + if A.c.stype != 0 + return false + else + i = symmetry(A, 0)[1] + return i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG + end + end +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 ishermitian(A::Sparse) + if isreal(A) + if A.c.stype != 0 + return true + else + i = symmetry(A, 0)[1] + return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG + end + else + if A.c.stype != 0 + return true + else + i = symmetry(A, 0)[1] + return i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG + end end - return lufact(A) end + +end #module diff --git a/base/sparse/cholmod_h.jl b/base/sparse/cholmod_h.jl index 52c3572f9f966..71d6ee96efa4a 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,27 @@ 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) 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 1f68cb9f5b788..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? @@ -881,10 +881,10 @@ getindex(A::SparseMatrixCSC, I::(Integer,Integer)) = getindex(A, I[1], I[2]) function getindex{T}(A::SparseMatrixCSC{T}, i0::Integer, i1::Integer) if !(1 <= i0 <= A.m && 1 <= i1 <= A.n); throw(BoundsError()); end - r1 = int(A.colptr[i1]);println("r1: $r1") - r2 = int(A.colptr[i1+1]-1);println("r2: $r2") + r1 = int(A.colptr[i1]) + r2 = int(A.colptr[i1+1]-1) (r1 > r2) && return zero(T) - r1 = searchsortedfirst(A.rowval, i0, r1, r2, Forward);println("r1: $r1") + r1 = searchsortedfirst(A.rowval, i0, r1, r2, Forward) ((r1 > r2) || (A.rowval[r1] != i0)) ? zero(T) : A.nzval[r1] end 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..f8ee3abae1c81 --- /dev/null +++ b/test/sparse/cholmod.jl @@ -0,0 +1,162 @@ +using Base.Test + +using Base.SparseMatrix.CHOLMOD + +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 = Sparse(int32([0,1,2,3,6,9,12,15,18,20,25,30,34,36,39,43,47,52,58,62,67,71,77,84,90,93,95, + 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 = ldltfact(A) # LDL' form +@test isvalid(chma) +x = chma\B +@test_approx_eq x.mat ones(size(x)) + +chma = cholfact(A) # LL' form +@test isvalid(chma) +x = chma\B +@test_approx_eq x.mat ones(size(x)) + +#lp_afiro example +afiro = Sparse(int32([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,25,27,29,33,37, + 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 + +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.parse(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 + +# 9915 +@test speye(2)\speye(2) == eye(2) 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] From c35f6441bfd2d0a81bb8fbf427584d94fb1a631a Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sat, 24 Jan 2015 17:35:57 -0500 Subject: [PATCH 03/10] Remove Triplet support from cholmod.jl --- base/sparse/cholmod.jl | 174 +++++++++++++---------------------------- test/Makefile | 2 +- 2 files changed, 57 insertions(+), 119 deletions(-) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 9dba3dc50a453..02610560672d3 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -13,7 +13,6 @@ export Dense, Factor, Sparse, - Triplet, etree, sparse! @@ -32,22 +31,6 @@ else ccall((:jl_cholmod_version, :libsuitesparse_wrapper), Cint, (Ptr{Cint},), cholmod_ver) end const cholmod_com_sz = ccall((:jl_cholmod_common_size,:libsuitesparse_wrapper),Int,()) -const cholmod_com = fill(0xff, cholmod_com_sz) -const cholmod_l_com = fill(0xff, cholmod_com_sz) -## cholmod_com and cholmod_l_com must be initialized at runtime because they contain pointers -## to functions in libc.so, whose addresses can change -function common(::Type{Int32}) - if isnan(reinterpret(Float64,cholmod_com[1:8])[1]) - @isok ccall((:cholmod_start, :libcholmod), Cint, (Ptr{UInt8},), cholmod_com) - end - cholmod_com -end -function common(::Type{Int64}) - if isnan(reinterpret(Float64,cholmod_l_com[1:8])[1]) - @isok ccall((:cholmod_l_start, :libcholmod), Cint, (Ptr{UInt8},), cholmod_l_com) - end - cholmod_l_com -end type CHOLMODException <: Exception msg::AbstractString @@ -56,6 +39,18 @@ 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 +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 +end + ### A way of examining some of the fields in cholmod_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. @@ -85,20 +80,20 @@ end const cholmod_com_offsets = Array(Int, length(Common.types)) ccall((:jl_cholmod_common_offsets, :libsuitesparse_wrapper), Void, (Ptr{Int},), cholmod_com_offsets) -const cholmod_final_ll_inds = (1:4) + cholmod_com_offsets[7] -const cholmod_prt_inds = (1:4) + cholmod_com_offsets[13] -const cholmod_ityp_inds = (1:4) + cholmod_com_offsets[18] +const common_final_ll = (1:4) + cholmod_com_offsets[7] +const common_print = (1:4) + cholmod_com_offsets[13] +const common_itype = (1:4) + cholmod_com_offsets[18] ### there must be an easier way but at least this works. -function Common(aa::Array{UInt8,1}) - typs = Common.types - sz = map(sizeof, typs) - args = map(i->reinterpret(typs[i], aa[cholmod_com_offsets[i] + (1:sz[i])])[1], 1:length(sz)) - eval(Expr(:call, unshift!(args, :Common), Any)) -end - -function set_cholmod_prt_lev(cm::Array{UInt8}, lev::Integer) # can probably be removed - cm[(1:4) + cholmod_com_offsets[13]] = reinterpret(UInt8, [int32(lev)]) +# function Common(aa::Array{UInt8,1}) +# typs = Common.types +# sz = map(sizeof, typs) +# args = map(i->reinterpret(typs[i], aa[cholmod_com_offsets[i] + (1:sz[i])])[1], 1:length(sz)) +# eval(Expr(:call, unshift!(args, :Common), Any)) +# end + +function set_print_level(cm::Array{UInt8}, lev::Integer) + cm[common_print] = reinterpret(UInt8, [int32(lev)]) end ## cholmod_dense pointers passed to or returned from C functions are of Julia type @@ -276,7 +271,7 @@ type c_Sparse{Tv<:VTypes,Ti<:ITypes} packed::Cint end -# Corresponds to the exact definition of cholmod_sparse_void in the library. Useful when reading files of unknown type. +# Corresponds to the exact definition of cholmod_sparse_struct in the library. Useful when reading matrices of unknown type from files as in cholmod_read_sparse type c_SparseVoid m::Csize_t n::Csize_t @@ -301,32 +296,9 @@ type Sparse{Tv<:VTypes,Ti<:ITypes} <: AbstractSparseMatrix{Tv,Ti} nzval::Vector{Tv} end -type c_Triplet{Tv<:VTypes,Ti<:ITypes} - 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 Triplet{Tv<:VTypes,Ti<:ITypes} - c::c_Triplet{Tv,Ti} - i::Vector{Ti} - j::Vector{Ti} - x::Vector{Tv} -end - - -eltype{T<:VTypes}(::Type{Dense{T}}) = T -eltype{T<:VTypes}(::Type{Factor{T}}) = T -eltype{T<:VTypes}(::Type{Sparse{T}}) = T +eltype{T<:VTypes}(A::Dense{T}) = T +eltype{T<:VTypes}(A::Factor{T}) = T +eltype{T<:VTypes}(A::Sparse{T}) = T ## The Dense constructor does not copy the contents, which is generally what you ## want as most uses of Dense objects are read-only. @@ -485,20 +457,6 @@ end Sparse{Tv<:VTypes}(D::Dense{Tv}) = Sparse(D,1) # default Ti is Int - -# Constructs a Triplet object from a pointer to a cholmod_triplet_struct created by CHOLMOD. Julia takes ownership of the memory allocated by CHOLMOD and the pointer can therefore be freed without memory leek. -function Triplet{Tv<:VTypes,Ti<:ITypes}(tp::Ptr{c_Triplet{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 = Triplet{Tv,Ti}(ctp, i, j, x) - - csp.packed == 1 && c_free(csp.nzpt) # Julia's Sparse is not using unpacked storage - c_free(tp) - ct -end - # Dense wrappers ### cholmod_core_h ### function zeros{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) @@ -629,20 +587,6 @@ for Ti in IndexTypes L end - function triplet_to_sparse{Tv<:VTypes,Ti<:$Ti}(T::Triplet{Tv,Ti}) - Sparse(ccall((@cholmod_name "triplet_to_sparse" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,Ti}}, - (Ptr{c_Triplet{Tv,Ti}}, Ptr{UInt8}), - &T.c, common($Ti))) - end - - function sparse_to_triplet{Tv<:VTypes,Ti<:$Ti}(A::Sparse{Tv,Ti}) - Triplet(ccall((@cholmod_name "sparse_to_triplet" $Ti - ,:libcholmod), Ptr{c_Triplet{Tv,Ti}}, - (Ptr{c_Sparse{Tv,Ti}}, Ptr{UInt8}), - &A.c, common($Ti))) - end - function check_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}) bool(ccall((@cholmod_name "check_sparse" $Ti ,:libcholmod), Cint, @@ -657,13 +601,6 @@ for Ti in IndexTypes &L.c, common($Ti))) end - function check_triplet{Tv<:VTypes}(T::Triplet{Tv,$Ti}) - bool(ccall((@cholmod_name "check_triplet" $Ti - ,:libcholmod), Cint, - (Ptr{c_Triplet{Tv,$Ti}}, Ptr{UInt8}), - &T.c, common($Ti))) - end - function analyze{Tv<:VTypes}(A::Sparse{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) Factor(ccall((@cholmod_name "analyze" $Ti ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, @@ -731,12 +668,6 @@ for Ti in IndexTypes (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), &A.c, common($Ti))) end - function copy_triplet{Tv<:VTypes}(T::Triplet{Tv,$Ti}) - Triplet(ccall((@cholmod_name "copy_triplet" $Ti - ,:libcholmod), Ptr{c_Triplet{Tv,$Ti}}, - (Ptr{c_Triplet{Tv,$Ti}}, Ptr{UInt8}), - &T.c, common($Ti))) - end function copy{Tv<:VTypes}(A::Sparse{Tv,$Ti}, stype::Integer, mode::Integer) Sparse(ccall((@cholmod_name "copy" $Ti ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, @@ -746,17 +677,21 @@ for Ti in IndexTypes ### 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.c, name, common($Ti)) + &A.c, name, cm) nothing end function print_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}, name::ASCIIString) + cm = common($Ti) + set_print_level(cm, 3) @isok ccall((@cholmod_name "print_factor" $Ti ,:libcholmod), Cint, (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), - &L.c, name, common($Ti)) + &L.c, name, cm) nothing end @@ -864,7 +799,6 @@ end Dense(A::Sparse) = sparse_to_dense(A) Sparse(A::Dense) = dense_to_sparse(A) -Sparse(A::Triplet) = triplet_to_sparse(A) function Sparse(L::Factor) if bool(L.c.is_ll) return factor_to_sparse(L) @@ -881,18 +815,14 @@ function Sparse(filename::ASCIIString) A end -Triplet(A::Sparse) = sparse_to_triplet(A) - show(io::IO, L::Factor) = print_factor(L, "") isvalid(A::Dense) = check_dense(A) isvalid(A::Sparse) = check_sparse(A) -isvalid(A::Triplet) = check_triplet(A) isvalid(A::Factor) = check_factor(A) copy(A::Dense) = copy_dense(A) copy(A::Sparse) = copy_sparse(A) -copy(A::Triplet) = copy_triplet(A) copy(A::Factor) = copy_factor(A) size(B::Dense) = size(B.mat) @@ -914,18 +844,26 @@ function getindex{T}(A::Sparse{T}, i0::Integer, i1::Integer) ((r1 > r2) || (A.rowval0[r1] + 1 != i0)) ? zero(T) : A.nzval[r1] end -full(A::Sparse) = Dense(A).mat -full(A::Dense) = A.mat -function sparse(A::Sparse) - 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)) +# Convert invalidates the input matrix +function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, A::Sparse{Tv,Ti}) + B = SparseMatrixCSC(A.c.m, A.c.n, increment!(A.colptr0), increment!(A.rowval0), A.nzval) + if A.c.stype != 0 + return isreal(A) ? Symmetric(B, A.c.stype < 0 ? :L : :U) : Hermitian(B, A.c.stype < 0 ? :L : :U) + end + return B end -function sparse!(A::Sparse) - SparseMatrixCSC(A.c.m, A.c.n, increment!(A.colptr0), increment!(A.rowval0), A.nzval) +sparse!{Tv,Ti}(A::Sparse{Tv,Ti}) = convert(SparseMatrixCSC{Tv,Ti}, A) +sparse(A::Sparse) = sparse!(copy(A)) +sparse(L::Factor) = sparse!(Sparse(L)) +sparse(D::Dense) = sparse!(Sparse(D)) +function full(A::Sparse) + B = Dense(A).mat + if A.c.stype != 0 + return isreal(A) ? Symmetric(AA.c.stype < 0 ? :L : :U) : Hermitian(B, A.c.stype < 0 ? :L : :U) + end + return B end -sparse(L::Factor) = sparse!(Sparse(L)) -sparse(D::Dense) = sparse!(Sparse(D)) -sparse(T::Triplet) = sparse!(Sparse(T)) +full(A::Dense) = A.mat ## Multiplication @@ -976,7 +914,7 @@ analyze(A::SparseMatrixCSC) = analyze(Sparse(A)) function cholfact(A::Sparse) cm = common(indtype(A)) ## may need to change final_asis as well as final_ll - cm[cholmod_final_ll_inds] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt F = analyze(A, cm) factorize(A, F, cm) F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) @@ -986,7 +924,7 @@ end function ldltfact(A::Sparse) cm = common(indtype(A)) ## may need to change final_asis as well as final_ll - cm[cholmod_final_ll_inds] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt + cm[common_final_ll] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt F = analyze(A, cm) factorize(A, F, cm) return F @@ -995,7 +933,7 @@ 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[cholmod_final_ll_inds] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt F = analyze(A, cm) factorize_p(A, β, Ti[], F, cm) F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) @@ -1005,7 +943,7 @@ 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[cholmod_final_ll_inds] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt + cm[common_final_ll] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt F = analyze(A, cm) factorize_p(A, β, Ti[], F, cm) F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) 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 From c9607ffcb45293c21e226039cf154e39544a74d7 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sat, 24 Jan 2015 17:53:58 -0500 Subject: [PATCH 04/10] Use VersionNumber when having version dependent code lines --- base/sparse/cholmod.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 02610560672d3..4379099494344 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -24,12 +24,13 @@ macro isok(A) :($A == TRUE || throw(CHOLMODException(""))) end -const cholmod_ver = Array(Cint, 3) +const version_array = Array(Cint, 3) if dlsym(dlopen("libcholmod"), :cholmod_version) != C_NULL - ccall((:cholmod_version, :libcholmod), Cint, (Ptr{Cint},), cholmod_ver) + ccall((:cholmod_version, :libcholmod), Cint, (Ptr{Cint},), version_array) else - ccall((:jl_cholmod_version, :libsuitesparse_wrapper), Cint, (Ptr{Cint},), cholmod_ver) + 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,()) type CHOLMODException <: Exception @@ -115,7 +116,7 @@ type Dense{T<:VTypes} <: AbstractMatrix{T} mat::Matrix{T} end -if (1000cholmod_ver[1]+cholmod_ver[2]) >= 2001 # CHOLMOD version 2.1.0 or later +if version >= v"2.1.0" # CHOLMOD version 2.1.0 or later type c_Factor{Tv<:VTypes,Ti<:ITypes} n::Int minor::Int @@ -1077,14 +1078,14 @@ function issym(A::Sparse) if A.c.stype != 0 return true else - i = symmetry(A, 0)[1] + i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG end else if A.c.stype != 0 return false else - i = symmetry(A, 0)[1] + i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 return i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG end end @@ -1095,14 +1096,14 @@ function ishermitian(A::Sparse) if A.c.stype != 0 return true else - i = symmetry(A, 0)[1] + i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG end else if A.c.stype != 0 return true else - i = symmetry(A, 0)[1] + i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 return i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG end end From 130fb7745ad48f126717c56ec3ac3c2b992a0643 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sat, 24 Jan 2015 20:28:16 -0500 Subject: [PATCH 05/10] Fix 9915 and getindex for symmetric/Hermitian CHOLMOD.Sparse. Add tests. --- base/sparse/cholmod.jl | 54 +++++++++++++++++++----------------------- test/sparse/cholmod.jl | 25 ++++++++++++++++--- 2 files changed, 47 insertions(+), 32 deletions(-) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 4379099494344..52394096fc74b 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -310,9 +310,6 @@ function Dense{T<:VTypes}(A::VecOrMat{T}) # uses the memory from Julia end function Dense{T<:VTypes}(c::Ptr{c_Dense{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 = Dense(cp, pointer_to_array(cp.xpt, (cp.m,cp.n), true)) c_free(c) @@ -326,15 +323,7 @@ function Sparse{Tv<:VTypes,Ti<:ITypes}(colpt::Vector{Ti}, 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) @@ -343,25 +332,28 @@ function Sparse{Tv<:VTypes,Ti<:ITypes}(colpt::Vector{Ti}, 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)]")) - # end + it = ityp(Ti) - cs = Sparse(c_Sparse{Tv,Ti}(m, n, int(nz), convert(Ptr{Ti}, colpt0), - convert(Ptr{Ti}, rowval0), C_NULL, - convert(Ptr{Tv}, nzval), C_NULL, - int32(stype), ityp(Ti), - xtyp(Tv), dtyp(Tv), - TRUE, TRUE), - colpt0, rowval0, nzval) + cs = Sparse(c_Sparse{Tv,Ti}(m, + n, + int(nz), + convert(Ptr{Ti}, colpt0), + convert(Ptr{Ti}, rowval0), + C_NULL, + convert(Ptr{Tv}, nzval), + C_NULL, + int32(stype), + ityp(Ti), + xtyp(Tv), + dtyp(Tv), + TRUE, + TRUE), + colpt0, + rowval0, + nzval) @isok isvalid(cs) - # cs = sort!(cs) - return cs end function Sparse{Tv<:VTypes,Ti<:ITypes}(A::SparseMatrixCSC{Tv,Ti}, stype::Signed) @@ -837,7 +829,10 @@ getindex(A::Dense, i::Integer) = A.mat[i] getindex(A::Dense, i::Integer, j::Integer) = A.mat[i, j] getindex(A::Sparse, i::Integer) = getindex(A, ind2sub(size(A),i)...) function getindex{T}(A::Sparse{T}, i0::Integer, i1::Integer) - if !(1 <= i0 <= A.c.m && 1 <= i1 <= A.c.n); throw(BoundsError()); end + !(1 <= i0 <= A.c.m && 1 <= i1 <= A.c.n) && throw(BoundsError()) + A.c.stype < 0 && i0 < i1 && return conj(A[i1,i0]) + A.c.stype > 0 && i0 > i1 && return conj(A[i1,i0]) + r1 = int(A.colptr0[i1] + 1) r2 = int(A.colptr0[i1 + 1]) (r1 > r2) && return zero(T) @@ -972,7 +967,8 @@ ldltfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}} (\)(L::Factor, b::Vector) = reshape(solve(CHOLMOD_A, L, Dense(b)).mat, length(b)) (\)(L::Factor, B::Matrix) = solve(CHOLMOD_A, L, Dense(B)).mat (\)(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B) -(\)(L::Factor, B::SparseMatrixCSC) = sparse!(spsolve(CHOLMOD_A, L, Sparse(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) = solve(CHOLMOD_A, L, Dense(B)).mat diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index f8ee3abae1c81..737321a13de8e 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -2,8 +2,6 @@ using Base.Test using Base.SparseMatrix.CHOLMOD -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")) @@ -158,5 +156,26 @@ let # Issue 9160 end end -# 9915 +# 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 = 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("tmp.mtx") == [1 0 0;0 1 0;0 0.5 1] + From e97b2989e65d4526c1027d92057e00fc9dd3b7df Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 25 Jan 2015 11:51:45 -0500 Subject: [PATCH 06/10] Add some more tests and define linear indexing for symmetric matrices --- base/linalg/symmetric.jl | 1 + base/sparse/cholmod.jl | 2 +- test/sparse/cholmod.jl | 3 ++- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index c0567982604d8..18071597db39b 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::Symmetric, i::Integer) = ((q, r) = divrem(i - 1, size(A, 1)); A[r + 1, q + 1]) getindex(A::Symmetric, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.data, i, j) : getindex(A.data, j, i) getindex(A::Hermitian, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.data, i, j) : conj(getindex(A.data, j, i)) full(A::Symmetric) = copytri!(copy(A.data), A.uplo) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 52394096fc74b..8e9afa0293144 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -855,7 +855,7 @@ sparse(D::Dense) = sparse!(Sparse(D)) function full(A::Sparse) B = Dense(A).mat if A.c.stype != 0 - return isreal(A) ? Symmetric(AA.c.stype < 0 ? :L : :U) : Hermitian(B, A.c.stype < 0 ? :L : :U) + return isreal(A) ? Symmetric(B, A.c.stype < 0 ? :L : :U) : Hermitian(B, A.c.stype < 0 ? :L : :U) end return B end diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index 737321a13de8e..f4ff95d67b9d2 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -177,5 +177,6 @@ B = Sparse(SparseMatrixCSC{Float64,Int32}(sprandn(48, 48, 0.1))) # A has Int32 i # 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("tmp.mtx") == [1 0 0;0 1 0;0 0.5 1] +@test full(Sparse("tmp.mtx")) == [1 0 0;0 1 0.5;0 0.5 1] +run(`rm tmp.mtx`) From eca46f6114743b334c9b59653e79cf20f849646d Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sat, 7 Feb 2015 16:45:24 -0500 Subject: [PATCH 07/10] Let CHOLMOD handle memory allocation. Increase test coverage and remove some unused methods. --- base/linalg/symmetric.jl | 2 +- base/sparse/cholmod.jl | 1390 +++++++++++++++++++------------------- test/sparse/cholmod.jl | 280 +++++++- 3 files changed, 951 insertions(+), 721 deletions(-) diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index 18071597db39b..5475c8c9433f0 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -13,7 +13,7 @@ typealias HermOrSym{T,S} Union(Hermitian{T,S}, Symmetric{T,S}) typealias RealHermSymComplexHerm{T<:Real,S} Union(Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}) size(A::HermOrSym, args...) = size(A.data, args...) -getindex(A::Symmetric, i::Integer) = ((q, r) = divrem(i - 1, size(A, 1)); A[r + 1, q + 1]) +getindex(A::HermOrSym, i::Integer) = ((q, r) = divrem(i - 1, size(A, 1)); A[r + 1, q + 1]) getindex(A::Symmetric, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.data, i, j) : getindex(A.data, j, i) getindex(A::Hermitian, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.data, i, j) : conj(getindex(A.data, j, i)) full(A::Symmetric) = copytri!(copy(A.data), A.uplo) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 8e9afa0293144..afbc94c31387d 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -1,23 +1,24 @@ module CHOLMOD -import Base: (*), convert, copy, ctranspose, eltype, getindex, hcat, isvalid, show, size, - sort!, transpose, vcat +import Base: (*), convert, copy, eltype, getindex, show, size import Base.LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B, - A_ldiv_B!, cholfact, cholfact!, copy, det, diag, - full, ishermitian, isposdef!, issym, ldltfact, logdet, norm, scale, scale! + A_ldiv_B!, cholfact, cholfact!, det, diag, ishermitian, isposdef, + issym, isvalid, ldltfact, logdet -import Base.SparseMatrix: findnz, sparse +import Base.SparseMatrix: sparse export Dense, Factor, - Sparse, - etree, - sparse! + Sparse using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, increment!, indtype, decrement, decrement! +######### +# Setup # +######### + include("cholmod_h.jl") macro isok(A) @@ -84,46 +85,84 @@ ccall((:jl_cholmod_common_offsets, :libsuitesparse_wrapper), const common_final_ll = (1:4) + cholmod_com_offsets[7] const common_print = (1:4) + cholmod_com_offsets[13] const common_itype = (1:4) + cholmod_com_offsets[18] - -### there must be an easier way but at least this works. -# function Common(aa::Array{UInt8,1}) -# typs = Common.types -# sz = map(sizeof, typs) -# args = map(i->reinterpret(typs[i], aa[cholmod_com_offsets[i] + (1:sz[i])])[1], 1:length(sz)) -# eval(Expr(:call, unshift!(args, :Common), Any)) -# end +const common_dtype = (1:4) + cholmod_com_offsets[19] function set_print_level(cm::Array{UInt8}, lev::Integer) cm[common_print] = reinterpret(UInt8, [int32(lev)]) end -## cholmod_dense pointers passed to or returned from C functions are of Julia type -## Ptr{c_Dense}. The Dense type contains a c_Dense object and other -## fields then ensure the memory pointed to is freed when it should be and not before. -type c_Dense{T<:VTypes} - m::Int - n::Int - nzmax::Int - lda::Int - xpt::Ptr{T} - zpt::Ptr{Void} +#################### +# Type definitions # +#################### + +# The three core data types for CHOLMOD: Dense, Sparse and Factor. CHOLMOD manages the memory, so the Julia versions only wrap a pointer to a struct. Therefore finalizers should be registret each time a pointer is returned from CHOLMOD. + +# Dense +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 + +type Dense{T<:VTypes} <: DenseMatrix{T} + p::Ptr{C_Dense{T}} +end + +# 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 + +# 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 Dense{T<:VTypes} <: AbstractMatrix{T} - c::c_Dense - mat::Matrix{T} +type Sparse{Tv<:VTypes,Ti<:ITypes} <: AbstractSparseMatrix{Tv,Ti} + p::Ptr{C_Sparse{Tv,Ti}} end +# Factor + if version >= v"2.1.0" # CHOLMOD version 2.1.0 or later - type c_Factor{Tv<:VTypes,Ti<:ITypes} - n::Int - minor::Int + 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} @@ -131,11 +170,11 @@ if version >= v"2.1.0" # CHOLMOD version 2.1.0 or later nz::Ptr{Ti} next::Ptr{Ti} prev::Ptr{Ti} - nsuper::Int - ssize::Int - xsize::Int - maxcsize::Int - maxesize::Int + nsuper::Csize_t + ssize::Csize_t + xsize::Csize_t + maxcsize::Csize_t + maxesize::Csize_t super::Ptr{Ti} pi::Ptr{Ti} px::Ptr{Ti} @@ -148,51 +187,13 @@ if version >= v"2.1.0" # CHOLMOD version 2.1.0 or later xtype::Cint dtype::Cint end - - type Factor{Tv<:VTypes,Ti<:ITypes} - c::c_Factor{Tv,Ti} - Perm::Vector{Ti} - ColCount::Vector{Ti} - IPerm::Vector{Ti} - p::Vector{Ti} - i::Vector{Ti} - x::Vector{Tv} - nz::Vector{Ti} - next::Vector{Ti} - prev::Vector{Ti} - super::Vector{Ti} - pi::Vector{Ti} - px::Vector{Ti} - s::Vector{Ti} - end - - function Factor{Tv<:VTypes,Ti<:ITypes}(cp::Ptr{c_Factor{Tv,Ti}}) - cfp = unsafe_load(cp) - Perm = pointer_to_array(cfp.Perm, (cfp.n,), true) - ColCount = pointer_to_array(cfp.ColCount, (cfp.n,), true) - IPerm = pointer_to_array(cfp.IPerm, (cfp.IPerm == C_NULL ? 0 : cfp.n + 1,), true) - p = pointer_to_array(cfp.p, (cfp.p == C_NULL ? 0 : cfp.n + 1,), true) - i = pointer_to_array(cfp.i, (cfp.i == C_NULL ? 0 : cfp.nzmax,), true) - x = pointer_to_array(cfp.x, (cfp.x == C_NULL ? 0 : max(cfp.nzmax,cfp.xsize),), true) - nz = pointer_to_array(cfp.nz, (cfp.nz == C_NULL ? 0 : cfp.n,), true) - next = pointer_to_array(cfp.next, (cfp.next == C_NULL ? 0 : cfp.n + 2,), true) - prev = pointer_to_array(cfp.prev, (cfp.prev == C_NULL ? 0 : cfp.n + 2,), true) - super = pointer_to_array(cfp.super, (cfp.super == C_NULL ? 0 : cfp.nsuper + 1,), true) - pi = pointer_to_array(cfp.pi, (cfp.pi == C_NULL ? 0 : cfp.nsuper + 1,), true) - px = pointer_to_array(cfp.px, (cfp.px == C_NULL ? 0 : cfp.nsuper + 1,), true) - s = pointer_to_array(cfp.s, (cfp.s == C_NULL ? 0 : cfp.ssize + 1,), true) - cf = Factor{Tv,Ti}(cfp, Perm, ColCount, IPerm, p, i, x, nz, next, prev, - super, pi, px, s) - c_free(cp) - cf - end else - type c_Factor{Tv<:VTypes,Ti<:ITypes} - n::Int - minor::Int + 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} @@ -200,11 +201,11 @@ else nz::Ptr{Ti} next::Ptr{Ti} prev::Ptr{Ti} - nsuper::Int - ssize::Int - xsize::Int - maxcsize::Int - maxesize::Int + nsuper::Csize_t + ssize::Csize_t + xsize::Csize_t + maxcsize::Csize_t + maxesize::Csize_t super::Ptr{Ti} pi::Ptr{Ti} px::Ptr{Ti} @@ -217,522 +218,341 @@ else xtype::Cint dtype::Cint end - - type Factor{Tv<:VTypes,Ti<:ITypes} - c::c_Factor{Tv,Ti} - Perm::Vector{Ti} - ColCount::Vector{Ti} - p::Vector{Ti} - i::Vector{Ti} - x::Vector{Tv} - nz::Vector{Ti} - next::Vector{Ti} - prev::Vector{Ti} - super::Vector{Ti} - pi::Vector{Ti} - px::Vector{Ti} - s::Vector{Ti} - end - - function Factor{Tv<:VTypes,Ti<:ITypes}(cp::Ptr{c_Factor{Tv,Ti}}) - cfp = unsafe_load(cp) - Perm = pointer_to_array(cfp.Perm, (cfp.n,), true) - ColCount = pointer_to_array(cfp.ColCount, (cfp.n,), true) - p = pointer_to_array(cfp.p, (cfp.p == C_NULL ? 0 : cfp.n + 1,), true) - i = pointer_to_array(cfp.i, (cfp.i == C_NULL ? 0 : cfp.nzmax,), true) - x = pointer_to_array(cfp.x, (cfp.x == C_NULL ? 0 : max(cfp.nzmax,cfp.xsize),), true) - nz = pointer_to_array(cfp.nz, (cfp.nz == C_NULL ? 0 : cfp.n,), true) - next = pointer_to_array(cfp.next, (cfp.next == C_NULL ? 0 : cfp.n + 2,), true) - prev = pointer_to_array(cfp.prev, (cfp.prev == C_NULL ? 0 : cfp.n + 2,), true) - super = pointer_to_array(cfp.super, (cfp.super == C_NULL ? 0 : cfp.nsuper + 1,), true) - pi = pointer_to_array(cfp.pi, (cfp.pi == C_NULL ? 0 : cfp.nsuper + 1,), true) - px = pointer_to_array(cfp.px, (cfp.px == C_NULL ? 0 : cfp.nsuper + 1,), true) - s = pointer_to_array(cfp.s, (cfp.s == C_NULL ? 0 : cfp.ssize + 1,), true) - cf = Factor{Tv,Ti}(cfp, Perm, ColCount, p, i, x, nz, next, prev, - super, pi, px, s) - c_free(cp) - cf - end -end - -type c_Sparse{Tv<:VTypes,Ti<:ITypes} - m::Csize_t - n::Csize_t - nzmax::Csize_t - ppt::Ptr{Ti} - ipt::Ptr{Ti} - nzpt::Ptr{Ti} - xpt::Ptr{Tv} - zpt::Ptr{Void} - stype::Cint - itype::Cint - xtype::Cint - dtype::Cint - sorted::Cint - packed::Cint end -# Corresponds to the exact definition of cholmod_sparse_struct in the library. Useful when reading matrices of unknown type from files as in cholmod_read_sparse -type c_SparseVoid - m::Csize_t - n::Csize_t - nzmax::Csize_t - ppt::Ptr{Void} - ipt::Ptr{Void} - nzpt::Ptr{Void} - xpt::Ptr{Void} - zpt::Ptr{Void} - stype::Cint - itype::Cint - xtype::Cint - dtype::Cint - sorted::Cint - packed::Cint +type Factor{Tv,Ti} <: Factorization{Tv} + p::Ptr{C_Factor{Tv,Ti}} end -type Sparse{Tv<:VTypes,Ti<:ITypes} <: AbstractSparseMatrix{Tv,Ti} - c::c_Sparse{Tv,Ti} - colptr0::Vector{Ti} - rowval0::Vector{Ti} - nzval::Vector{Tv} -end - -eltype{T<:VTypes}(A::Dense{T}) = T -eltype{T<:VTypes}(A::Factor{T}) = T -eltype{T<:VTypes}(A::Sparse{T}) = T +################# +# Thin wrappers # +################# -## The Dense constructor does not copy the contents, which is generally what you -## want as most uses of Dense objects are read-only. -function Dense{T<:VTypes}(A::VecOrMat{T}) # uses the memory from Julia - m, n = size(A,1), size(A,2) - Dense(c_Dense{T}(m, n, m*n, stride(A, 2), convert(Ptr{T}, A), C_NULL, xtyp(T), dtyp(T)), - length(size(A)) == 2 ? A : reshape(A, (m, n))) -end -function Dense{T<:VTypes}(c::Ptr{c_Dense{T}}) - cp = unsafe_load(c) - ## the true in the call to pointer_to_array means Julia will free the memory - val = Dense(cp, pointer_to_array(cp.xpt, (cp.m,cp.n), true)) - c_free(c) - val -end +# 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 -function Sparse{Tv<:VTypes,Ti<:ITypes}(colpt::Vector{Ti}, - rowval::Vector{Ti}, - nzval::Vector{Tv}, - m::Integer, - n::Integer, - stype::Signed) - bb = colpt[1] - - if bool(bb) # one-based - colpt0 = decrement(colpt) - rowval0 = decrement(rowval) - else # zero based - colpt0 = colpt - rowval0 = rowval - end - nz = colpt0[end] - - it = ityp(Ti) - cs = Sparse(c_Sparse{Tv,Ti}(m, - n, - int(nz), - convert(Ptr{Ti}, colpt0), - convert(Ptr{Ti}, rowval0), - C_NULL, - convert(Ptr{Tv}, nzval), - C_NULL, - int32(stype), - ityp(Ti), - xtyp(Tv), - dtyp(Tv), - TRUE, - TRUE), - colpt0, - rowval0, - nzval) - - @isok isvalid(cs) - - return cs +### 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 -function Sparse{Tv<:VTypes,Ti<:ITypes}(A::SparseMatrixCSC{Tv,Ti}, stype::Signed) - Sparse(A.colptr, A.rowval, A.nzval, size(A,1), size(A,2), stype) +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 -function Sparse(A::SparseMatrixCSC) - AC = Sparse(A, 0) - i = symmetry(AC, 0)[1] # Check for Hermitianity - if (isreal(A) && (i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG)) || # real case - (i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG) # complex case - AC.c.stype = -1 - end - AC -end -Sparse{Ti<:ITypes}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) -Sparse{Ti<:ITypes}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) +free_dense!{T}(p::Ptr{C_Dense{T}}) = ccall((:cholmod_free_dense, :libcholmod), Cint, (Ptr{Ptr{C_Dense{T}}}, Ptr{Void}), &p, common(Cint)) -# Constructs a Sparse object from a pointer to a cholmod_sparse_struct created by CHOLMOD. Julia takes ownership of the memory allocated by CHOLMOD and the pointer can therefore be freed without memory leek. -function Sparse{Tv<:VTypes,Ti<:ITypes}(cp::Ptr{c_Sparse{Tv,Ti}}) - csp = unsafe_load(cp) - - # Take control of the memory allocated by CHOLMOD - colptr0 = pointer_to_array(csp.ppt, (csp.n + 1,), true) - rowval0 = pointer_to_array(csp.ipt, (csp.nzmax,), true) - nzval = pointer_to_array(csp.xpt, (csp.nzmax,), true) - - cms = Sparse{Tv,Ti}(csp, colptr0, rowval0, nzval) - - # Free memory - csp.packed == 1 && c_free(csp.nzpt) # Julia's Sparse is not using unpacked storage - c_free(cp) - - cms -end -# Useful when reading in files, but not type stable -function Sparse(p::Ptr{c_SparseVoid}) - s = unsafe_load(p) - - # Check integer type - if s.itype == INT - Ti = Cint - elseif s.itype == INTLONG - throw(CHOLMODException("the value of itype was $s.itype. This combination of integer types shouldn't happen. Please submit a bug report.")) - elseif s.itype == LONG - Ti = SuiteSparse_long - else - throw(CHOLMODException("illegal value of itype: $s.itype")) - end - - # Check for double or single precision - if s.dtype == DOUBLE - Tv = Float64 - elseif s.dtype == SINGLE - Tv = Float32 - else - throw(CHOLMODException("illegal value of dtype: $s.dtype")) - end - - # Check for real or complex - if s.xtype == COMPLEX - Tv = Complex{Tv} - elseif s.xtype != REAL - throw(CHOLMODException("illegal value of xtype: $s.xtype")) - end - - # Take control of the memory allocated by CHOLMOD - colptr0 = pointer_to_array(convert(Ptr{Ti}, s.ppt), (s.n + 1,), true) - rowval0 = pointer_to_array(convert(Ptr{Ti}, s.ipt), (s.nzmax,), true) - nzval = pointer_to_array(convert(Ptr{Tv}, s.xpt), (s.nzmax,), true) - - r = Sparse{Tv,Ti}(c_Sparse{Tv,Ti}(s.m, - s.n, - s.nzmax, - pointer(colptr0), - pointer(rowval0), - C_NULL, - pointer(nzval), - C_NULL, - s.stype, - s.itype, - s.xtype, - s.dtype, - s.sorted, - s.packed), - colptr0, rowval0, nzval) - - # Free memory - s.packed == 1 && c_free(s.nzpt) # Julia's Sparse is not using unpacked storage - c_free(p) - - r -end - -Sparse{Tv<:VTypes}(D::Dense{Tv}) = Sparse(D,1) # default Ti is Int - -# Dense wrappers -### cholmod_core_h ### function zeros{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) - Dense(ccall((:cholmod_zeros, :libcholmod), Ptr{c_Dense{T}}, + d = Dense(ccall((:cholmod_zeros, :libcholmod), Ptr{C_Dense{T}}, (Csize_t, Csize_t, Cint, Ptr{UInt8}), - m, n, xtyp(T), common(Int32))) + m, n, xtyp(T), common(Cint))) + finalizer(d, free!) + d end zeros(m::Integer, n::Integer) = zeros(m, n, Float64) function ones{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) - Dense(ccall((:cholmod_ones, :libcholmod), Ptr{c_Dense{T}}, + d = Dense(ccall((:cholmod_ones, :libcholmod), Ptr{C_Dense{T}}, (Csize_t, Csize_t, Cint, Ptr{UInt8}), - m, n, xtyp(T), common(Int32))) + m, n, xtyp(T), common(Cint))) + finalizer(d, free!) + d end ones(m::Integer, n::Integer) = ones(m, n, Float64) function eye{T<:VTypes}(m::Integer, n::Integer, ::Type{T}) - Dense(ccall((:cholmod_eye, :libcholmod), Ptr{c_Dense{T}}, + d = Dense(ccall((:cholmod_eye, :libcholmod), Ptr{C_Dense{T}}, (Csize_t, Csize_t, Cint, Ptr{UInt8}), - m, n, xtyp(T), common(Int32))) + m, n, xtyp(T), common(Cint))) + finalizer(d, free!) + d end eye(m::Integer, n::Integer) = eye(m, n, Float64) eye(n::Integer) = eye(n, n, Float64) function copy_dense{Tv<:VTypes}(A::Dense{Tv}) - Dense(ccall((:cholmod_copy_dense, :libcholmod), Ptr{c_Dense{Tv}}, - (Ptr{c_Dense{Tv}}, Ptr{UInt8}), - &A.c, common(Cint))) + d = Dense(ccall((:cholmod_copy_dense, :libcholmod), Ptr{C_Dense{Tv}}, + (Ptr{C_Dense{Tv}}, Ptr{UInt8}), + A.p, common(Cint))) + finalizer(d, free!) + d end ### cholmod_matrixops.h ### function norm_dense{Tv<:VTypes}(D::Dense{Tv}, p::Integer) + s = unsafe_load(D.p) + if p == 2 + if s.ncol > 1 + throw(ArgumentError("2 norm only supported when matrix has one column")) + end + elseif p != 0 && p != 1 + throw(ArgumentError("second argument must be either 0 (Inf norm), 1, or 2")) + end ccall((:cholmod_norm_dense, :libcholmod), Cdouble, - (Ptr{c_Dense{Tv}}, Cint, Ptr{UInt8}), - &D.c, p, common(Int32)) + (Ptr{C_Dense{Tv}}, Cint, Ptr{UInt8}), + D.p, p, common(Cint)) end ### cholmod_check.h ### function check_dense{T<:VTypes}(A::Dense{T}) bool(ccall((:cholmod_check_dense, :libcholmod), Cint, - (Ptr{c_Dense{T}}, Ptr{UInt8}), - &A.c, common(Cint))) + (Ptr{C_Dense{T}}, Ptr{UInt8}), + A.p, common(Cint))) end # Non-Dense wrappers (which all depend on IType) for Ti in IndexTypes @eval begin - # This method is probably not necessary as memory allocated by CHOLMOD is usually overtaken by Julia when cosntructing Julia instances from pointers. - function free_sparse{Tv<:VTypes}(ptr::Ptr{c_Sparse{Tv,$Ti}}) - cm = common($Ti) - @isok ccall((@cholmod_name "free_sparse" $Ti - ,:libcholmod), Cint, - (Ptr{Ptr{c_Sparse{Tv,$Ti}}}, Ptr{UInt8}), - &ptr, cm) + ### cholmod_core.h ### + function allocate_sparse(nrow::Integer, ncol::Integer, nzmax::Integer, sorted::Bool, packed::Bool, stype::Integer, ::Type{Float64}, ::Type{$Ti}) + s = Sparse(ccall((@cholmod_name("allocate_sparse", $Ti), :libcholmod), Ptr{C_Sparse{Float64,$Ti}}, + (Csize_t, Csize_t, Csize_t, Cint, + Cint, Cint, Cint, Ptr{Void}), + nrow, ncol, nzmax, sorted, + packed, stype, REAL, common($Ti))) + finalizer(s, free!) + s end - - function (*){Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}) - Sparse(ccall((@cholmod_name "ssmult" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Cint, Ptr{UInt8}), - &A.c, &B.c, 0, true, true, common($Ti))) + function allocate_sparse(nrow::Integer, ncol::Integer, nzmax::Integer, sorted::Bool, packed::Bool, stype::Integer, ::Type{Complex{Float64}}, ::Type{$Ti}) + s = Sparse(ccall((@cholmod_name("allocate_sparse", $Ti), :libcholmod), Ptr{C_Sparse{Complex{Float64},$Ti}}, + (Csize_t, Csize_t, Csize_t, Cint, + Cint, Cint, Cint, Ptr{Void}), + nrow, ncol, nzmax, sorted, + packed, stype, COMPLEX, common($Ti))) + finalizer(s, free!) + s end - - function transpose{Tv<:VTypes}(A::Sparse{Tv,$Ti}, values::Integer) - cm = common($Ti) - Sparse(ccall((@cholmod_name "transpose" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, values, cm)) + function free_sparse!{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 ssmult{Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, stype::Integer, values::Integer, sorted::Integer) - cm = common($Ti) - Sparse(ccall((@cholmod_name "ssmult" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Cint, Ptr{UInt8}), - &A.c, &B.c, stype, values, sorted, cm)) + function free_sparse!(ptr::Ptr{C_SparseVoid}) + @isok ccall((@cholmod_name("free_sparse", $Ti), :libcholmod), Cint, + (Ptr{Ptr{C_SparseVoid}}, Ptr{UInt8}), + &ptr, common($Ti)) end - function aat{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, fset::Vector{$Ti}, mode::Integer) - cm = common($Ti) - Sparse(ccall((@cholmod_name "aat" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{$Ti}, Csize_t, Cint, Ptr{UInt8}), - &A.c, fset, length(fset), mode, cm)) + function free_factor!{Tv<:VTypes}(ptr::Ptr{C_Factor{Tv,$Ti}}) + @isok ccall((@cholmod_name("free_factor", $Ti), :libcholmod), Cint, + (Ptr{Ptr{C_Factor{Tv,$Ti}}}, Ptr{Void}), + &ptr, common($Ti)) end - function copy{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, stype::Integer, mode::Integer) - cm = common($Ti) - Sparse(ccall((@cholmod_name "copy" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), - &A.c, 0, 1, cm)) + 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 sparse_to_dense{Tv<:VTypes}(A::Sparse{Tv,$Ti}) - Dense(ccall((@cholmod_name "sparse_to_dense" $Ti - ,:libcholmod), Ptr{c_Dense{Tv}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, common($Ti))) + d = Dense(ccall((@cholmod_name("sparse_to_dense", $Ti),:libcholmod), Ptr{C_Dense{Tv}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, common($Ti))) + finalizer(d, free!) + d end - function dense_to_sparse{Tv<:VTypes}(D::Dense{Tv},i::$Ti) - Sparse(ccall((@cholmod_name "dense_to_sparse" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Dense{Tv}}, Cint, Ptr{UInt8}), - &D.c, true, common($Ti))) + function dense_to_sparse{Tv<:VTypes}(D::Dense{Tv}, ::Type{$Ti}) + s = Sparse(ccall((@cholmod_name("dense_to_sparse", $Ti),:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Dense{Tv}}, Cint, Ptr{UInt8}), + D.p, true, common($Ti))) + finalizer(s, free!) + s end - function factor_to_sparse{Tv<:VTypes}(L::Factor{Tv,$Ti}) - Sparse(ccall((@cholmod_name "factor_to_sparse" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti))) + function factor_to_sparse!{Tv<:VTypes}(F::Factor{Tv,$Ti}) + ss = unsafe_load(F.p) + ss.xtype > PATTERN || throw(CHOLMODException("only numeric factors are supported")) + s = Sparse(ccall((@cholmod_name("factor_to_sparse", $Ti),:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + F.p, common($Ti))) + finalizer(s, free!) + s end - function copy_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) - Factor(ccall((@cholmod_name "copy_factor" $Ti - ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti))) + function change_factor{Tv<:VTypes}(::Type{Float64}, to_ll::Bool, to_super::Bool, to_packed::Bool, to_monotonic::Bool, F::Factor{Tv,$Ti}) + @isok ccall((@cholmod_name("change_factor", $Ti),:libcholmod), Cint, + (Cint, Cint, Cint, Cint, Cint, Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + REAL, to_ll, to_super, to_packed, to_monotonic, F.p, common($Ti)) + Factor{Float64,$Ti}(F.p) end - function change_factor{Tv<:VTypes}(to_xtype::Integer, to_ll::Bool, to_super::Bool, to_packed::Bool, to_monotonic::Bool, L::Factor{Tv,$Ti}) - @isok ccall((@cholmod_name "change_factor" $Ti - ,:libcholmod), Cint, - (Cint, Cint, Cint, Cint, Cint, Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - to_xtype, to_ll, to_super, to_packed, to_monotonic, L, cm) - L + 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 check_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}) - bool(ccall((@cholmod_name "check_sparse" $Ti - ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, common($Ti))) - end - - function check_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) - bool(ccall((@cholmod_name "check_factor" $Ti - ,:libcholmod), Cint, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti))) + bool(ccall((@cholmod_name("check_sparse", $Ti),:libcholmod), Cint, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, common($Ti))) end - function analyze{Tv<:VTypes}(A::Sparse{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) - Factor(ccall((@cholmod_name "analyze" $Ti - ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, C)) - end - - function factorize{Tv<:VTypes}(A::Sparse{Tv,$Ti}, F::Factor{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) - @isok ccall((@cholmod_name "factorize" $Ti - ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &A.c, &F.c, C) - F - end - - function factorize_p{Tv<:VTypes}(A::Sparse{Tv,$Ti}, β::Tv, fset::Vector{$Ti}, F::Factor{Tv,$Ti}, C::Vector{UInt8} = common($Ti)) - @isok ccall((@cholmod_name "factorize_p" $Ti - ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{Cdouble}, Ptr{$Ti}, Csize_t, - Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &A.c, &β, fset, length(fset), - &F.c, common($Ti)) - F - end - - ### cholmod_core.h ### - function pack_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) - @isok ccall((@cholmod_name "pack_factor" $Ti - ,:libcholmod), Cint, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti)) - L + 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 nnz{Tv<:VTypes}(A::Sparse{Tv,$Ti}) - ccall((@cholmod_name "nnz" $Ti - ,:libcholmod), Int, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, common($Ti)) + ccall((@cholmod_name("nnz", $Ti),:libcholmod), Int, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, common($Ti)) end function speye{Tv<:VTypes}(m::Integer, n::Integer, ::Type{Tv}) - Sparse(ccall((@cholmod_name "speye" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + s = Sparse(ccall((@cholmod_name("speye", $Ti), :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, (Csize_t, Csize_t, Cint, Ptr{UInt8}), m, n, xtyp(Tv), common($Ti))) + finalizer(s, free!) + s end function spzeros{Tv<:VTypes}(m::Integer, n::Integer, nzmax::Integer, ::Type{Tv}) - Sparse(ccall((@cholmod_name "spzeros" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, + s = Sparse(ccall((@cholmod_name("spzeros", $Ti), :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, (Csize_t, Csize_t, Csize_t, Cint, Ptr{UInt8}), m, n, nzmax, xtyp(Tv), common($Ti))) + finalizer(s, free!) + s end - function copy_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}) - Factor(ccall((@cholmod_name "copy_factor" $Ti - ,:libcholmod), Ptr{c_Factor{Tv,$Ti}}, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}), - &L.c, common($Ti))) + function transpose{Tv<:VTypes}(A::Sparse{Tv,$Ti}, values::Integer) + s = Sparse(ccall((@cholmod_name("transpose", $Ti),:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + A.p, values, common($Ti))) + finalizer(s, free!) + s + end + + function copy_factor{Tv<:VTypes}(F::Factor{Tv,$Ti}) + f = Factor(ccall((@cholmod_name("copy_factor", $Ti),:libcholmod), Ptr{C_Factor{Tv,$Ti}}, + (Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + F.p, common($Ti))) + finalizer(f, free!) + f end function copy_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}) - Sparse(ccall((@cholmod_name "copy_sparse" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &A.c, common($Ti))) + s = Sparse(ccall((@cholmod_name("copy_sparse", $Ti),:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, common($Ti))) + finalizer(s, free!) + s end - function copy{Tv<:VTypes}(A::Sparse{Tv,$Ti}, stype::Integer, mode::Integer) - Sparse(ccall((@cholmod_name "copy" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), - &A.c, stype, mode, common($Ti))) + function copy{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, stype::Integer, mode::Integer) + s = Sparse(ccall((@cholmod_name("copy", $Ti),:libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Cint, Cint, Ptr{UInt8}), + A.p, stype, mode, common($Ti))) + finalizer(s, free!) + s end ### cholmod_check.h ### 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.c, name, cm) + @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 print_factor{Tv<:VTypes}(L::Factor{Tv,$Ti}, name::ASCIIString) + function print_factor{Tv<:VTypes}(F::Factor{Tv,$Ti}, name::ASCIIString) cm = common($Ti) set_print_level(cm, 3) - @isok ccall((@cholmod_name "print_factor" $Ti - ,:libcholmod), Cint, - (Ptr{c_Factor{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), - &L.c, name, cm) + @isok ccall((@cholmod_name("print_factor", $Ti),:libcholmod), Cint, + (Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}, Ptr{UInt8}), + F.p, name, cm) nothing end ### cholmod_matrixops.h ### + function ssmult{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, stype::Integer, values::Bool, sorted::Bool) + lA = unsafe_load(A.p) + lB = unsafe_load(B.p) + 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 norm_sparse{Tv<:VTypes}(A::Sparse{Tv,$Ti}, norm::Integer) - norm == 0 || norm == 1 || throw(ArgumentError("norm argument must be either 0 or 1")) - ccall((@cholmod_name "norm_sparse" $Ti - , :libcholmod), Cdouble, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, norm, common($Ti)) + 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 horzcat{Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) - ccall((@cholmod_name "horzcat" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, &B.c, values, common($Ti)) + function horzcat{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) + s = Sparse(ccall((@cholmod_name("horzcat", $Ti), :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + A.p, B.p, values, common($Ti))) + finalizer(s, free!) + s end - function scale{Tv<:VTypes}(S::Dense{Tv}, scale::Integer, A::Sparse{Tv,$Ti}) - @isok ccall((@cholmod_name "scale" $Ti - ,:libcholmod), Cint, - (Ptr{c_Dense{Tv}}, Cint, Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - &S.c, scale, &A.c, common($Ti)) + 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 sdmult{Tv<:VTypes}(A::Sparse{Tv,$Ti}, transpose::Bool, α::Tv, β::Tv, X::Dense{Tv}, Y::Dense{Tv}) + 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))")) + if nc != size(X, 1) + throw(DimensionMismatch("incompatible dimensions, $nc and $(size(X,1))")) end - @isok ccall((@cholmod_name "sdmult" $Ti - ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, - Ptr{c_Dense{Tv}}, Ptr{c_Dense{Tv}}, Ptr{UInt8}), - &A.c, transpose, &α, &β, - &X.c, &Y.c, common($Ti)) + @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 vertcat{Tv<:VTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) - ccall((@cholmod_name "vertcat" $Ti - , :libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), - &A.c, &B.c, values, common($Ti)) + function vertcat{Tv<:VRealTypes}(A::Sparse{Tv,$Ti}, B::Sparse{Tv,$Ti}, values::Bool) + s = Sparse(ccall((@cholmod_name("vertcat", $Ti), :libcholmod), Ptr{C_Sparse{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}), + A.p, B.p, values, common($Ti))) + finalizer(s, free!) + s end function symmetry{Tv<:VTypes}(A::Sparse{Tv,$Ti}, option::Integer) @@ -740,75 +560,272 @@ for Ti in IndexTypes 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}, + 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.c, option, xmatched, pmatched, + A.p, option, xmatched, pmatched, nzoffdiag, nzdiag, common($Ti)) rv, xmatched[1], pmatched[1], nzoffdiag[1], nzdiag[1] end # cholmod_cholesky.h - function etree{Tv<:VTypes}(A::Sparse{Tv,$Ti}) - Parent = Array($Ti, size(A,2)) - @isok ccall((@cholmod_name "etree" $Ti - ,:libcholmod), Cint, - (Ptr{c_Sparse{Tv,$Ti}}, Ptr{$Ti}, Ptr{UInt8}), - &A.c, Parent, common($Ti)) - return Parent + # For analyze, factorize and factorize_p, the Common argument must be supplied in order to control if the factorization is LLt or LDLt + function analyze{Tv<:VTypes}(A::Sparse{Tv,$Ti}, cmmn::Vector{UInt8}) + f = Factor(ccall((@cholmod_name("analyze", $Ti),:libcholmod), Ptr{C_Factor{Tv,$Ti}}, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{UInt8}), + A.p, cmmn)) + finalizer(f, free!) + f + end + function factorize!{Tv<:VTypes}(A::Sparse{Tv,$Ti}, F::Factor{Tv,$Ti}, cmmn::Vector{UInt8}) + @isok ccall((@cholmod_name("factorize", $Ti),:libcholmod), Cint, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + A.p, F.p, cmmn) + F + end + function factorize_p!{Tv<:VTypes}(A::Sparse{Tv,$Ti}, β::Tv, fset::Vector{$Ti}, F::Factor{Tv,$Ti}, cmmn::Vector{UInt8}) + @isok ccall((@cholmod_name("factorize_p", $Ti),:libcholmod), Cint, + (Ptr{C_Sparse{Tv,$Ti}}, Ptr{Cdouble}, Ptr{$Ti}, Csize_t, + Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}), + A.p, &β, fset, length(fset), + F.p, cmmn) + F end - function solve{Tv<:VTypes}(sys::Integer, L::Factor{Tv,$Ti}, B::Dense{Tv}) - size(L,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(L,1)) rows, but RHS has $(size(B,1)) rows.")) - Dense(ccall((@cholmod_name "solve" $Ti - ,:libcholmod), Ptr{c_Dense{Tv}}, - (Cint, Ptr{c_Factor{Tv,$Ti}}, Ptr{c_Dense{Tv}}, Ptr{UInt8}), - sys, &L.c, &B.c, common($Ti))) + function solve{Tv<:VTypes}(sys::Integer, F::Factor{Tv,$Ti}, B::Dense{Tv}) + 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 spsolve{Tv<:VTypes}(sys::Integer, L::Factor{Tv,$Ti}, B::Sparse{Tv,$Ti}) - size(L,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(L,1)) rows, but RHS has $(size(B,1)) rows.")) - Sparse(ccall((@cholmod_name "spsolve" $Ti - ,:libcholmod), Ptr{c_Sparse{Tv,$Ti}}, - (Cint, Ptr{c_Factor{Tv,$Ti}}, Ptr{c_Sparse{Tv,$Ti}}, Ptr{UInt8}), - sys, &L.c, &B.c, common($Ti))) + function spsolve{Tv<:VTypes}(sys::Integer, F::Factor{Tv,$Ti}, B::Sparse{Tv,$Ti}) + 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 - end -end -# Autodetects the types -function read_sparse(file::CFILE) - Sparse(ccall((:cholmod_read_sparse, :libcholmod), Ptr{c_SparseVoid}, - (Ptr{Void}, Ptr{UInt8}), - file.ptr, common(Int32))) + # 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 ######################### # High level interfaces # ######################### -# Convertion +# Convertion/construction +function Dense{T<:VTypes}(A::VecOrMat{T}) + d = allocate_dense(size(A, 1), size(A, 2), stride(A, 2), T) + s = unsafe_load(d.p) + unsafe_copy!(s.x, pointer(A), length(A)) + d +end Dense(A::Sparse) = sparse_to_dense(A) -Sparse(A::Dense) = dense_to_sparse(A) -function Sparse(L::Factor) - if bool(L.c.is_ll) - return factor_to_sparse(L) +# This constructior assumes zero based colptr and rowval +function Sparse{Tv<:VTypes,Ti<:ITypes}(m::Integer, n::Integer, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}, stype) + + # check if columns are sorted + iss = true + for i = 2:length(colptr) + if !issorted(sub(rowval, colptr[i - 1] + 1:colptr[i])) + iss = false + break + end + end + + o = allocate_sparse(m, n, length(nzval), iss, true, stype, Tv, Ti) + s = unsafe_load(o.p) + + unsafe_copy!(s.p, pointer(colptr), length(colptr)) + unsafe_copy!(s.i, pointer(rowval), length(rowval)) + unsafe_copy!(s.x, pointer(nzval), length(nzval)) + + @isok check_sparse(o) + + return o + +end +function Sparse{Tv<:VTypes,Ti<:ITypes}(m::Integer, n::Integer, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) + o = Sparse(m, n, colptr, rowval, nzval, 0) + + # check if array is symmetric and change stype if it is + 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) + + return o +end +function Sparse{Tv<:VTypes,Ti<:ITypes}(A::SparseMatrixCSC{Tv,Ti}) + o = Sparse(A, 0) + # check if array is symmetric and change stype if it is + if ishermitian(o) + change_stype!(o, -1) + end + o +end + +Sparse{Ti<:ITypes}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1) +Sparse{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 + + # 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 - cm = common($Ti) - Lcll = copy_factor(L) - change_factor(L.c.xtype, true, L.c.is_super == 1, true, true, Lcll, L) - return factor_to_sparse(Lcll) + + # 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 + + return Sparse(convert(Ptr{C_Sparse{Tv,Ti}}, p)) + end + +Sparse(A::Dense) = dense_to_sparse(A, Cint) # default is Cint which is probably sufficient when converted from dense matrix +Sparse(L::Factor) = factor_to_sparse!(copy(L)) function Sparse(filename::ASCIIString) f = open(filename) - A = read_sparse(CFILE(f)) + A = read_sparse(CFILE(f), Int) close(f) A end -show(io::IO, L::Factor) = print_factor(L, "") +## convertion back to base Julia types +function convert{T<:VTypes}(::Type{Matrix}, D::Dense{T}) + s = unsafe_load(D.p) + a = Array(T, s.nrow, s.ncol) + if s.d == s.nrow + unsafe_copy!(pointer(a), s.x, s.d*s.ncol) + else + for j = 1:s.ncol + for i = 1:s.nrow + a[i,j] = unsafe_load(s.x, i + (j - 1)*s.d) + end + end + end + a +end + +function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, A::Sparse{Tv,Ti}) + s = unsafe_load(A.p) + 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 + +free!(A::Dense) = free_dense!(A.p) +free!(A::Sparse) = free_sparse!(A.p) +free!(F::Factor) = free_factor!(F.p) + +eltype{T<:VTypes}(A::Dense{T}) = T +eltype{T<:VTypes}(A::Factor{T}) = T +eltype{T<:VTypes}(A::Sparse{T}) = T + +function show(io::IO, F::Factor) + s = unsafe_load(F.p) + println(io, typeof(F)) + @printf(io, "type: %12s\n", bool(s.is_ll) ? "LLt" : "LDLt") + @printf(io, "method: %10s\n", bool(s.is_super) ? "supernodal" : "simplicial") + @printf(io, "maxnnz: %10d\n", int(s.nzmax)) + @printf(io, "nnz: %13d\n", nnz(Sparse(F))) +end isvalid(A::Dense) = check_dense(A) isvalid(A::Sparse) = check_sparse(A) @@ -818,59 +835,57 @@ copy(A::Dense) = copy_dense(A) copy(A::Sparse) = copy_sparse(A) copy(A::Factor) = copy_factor(A) -size(B::Dense) = size(B.mat) -size(B::Dense, d::Integer) = size(B.mat, d) -size(A::Sparse) = (int(A.c.m), int(A.c.n)) -size(A::Sparse, d::Integer) = d > 0 ? (d == 1 ? int(A.c.m) : (d == 2 ? int(A.c.n) : 1)) : BoundsError() -size(L::Factor) = (n = int(L.c.n); (n,n)) -size(L::Factor, d::Integer) = d < 1 ? throw(BoundsError("Dimension $d out of range")) : (d <= 2 ? int(L.c.n) : 1) +function size(A::Union(Dense,Sparse)) + s = unsafe_load(A.p) + return (int(s.nrow), int(s.ncol)) +end +function size(F::Factor, i::Integer) + 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 getindex(A::Dense, i::Integer) + s = unsafe_load(A.p) + 0 < i <= s.nrow*s.ncol || throw(BoundsError()) + unsafe_load(s.x, i) +end +function getindex(A::Dense, i::Integer, j::Integer) + s = unsafe_load(A.p) + 0 < i <= s.nrow || throw(BoundsError()) + 0 < j <= s.ncol || throw(BoundsError()) + unsafe_load(s.x, i + (j - 1)*s.d) +end -getindex(A::Dense, i::Integer) = A.mat[i] -getindex(A::Dense, i::Integer, j::Integer) = A.mat[i, j] getindex(A::Sparse, i::Integer) = getindex(A, ind2sub(size(A),i)...) function getindex{T}(A::Sparse{T}, i0::Integer, i1::Integer) - !(1 <= i0 <= A.c.m && 1 <= i1 <= A.c.n) && throw(BoundsError()) - A.c.stype < 0 && i0 < i1 && return conj(A[i1,i0]) - A.c.stype > 0 && i0 > i1 && return conj(A[i1,i0]) + s = unsafe_load(A.p) + !(1 <= i0 <= s.nrow && 1 <= i1 <= s.ncol) && throw(BoundsError()) + s.stype < 0 && i0 < i1 && return conj(A[i1,i0]) + s.stype > 0 && i0 > i1 && return conj(A[i1,i0]) - r1 = int(A.colptr0[i1] + 1) - r2 = int(A.colptr0[i1 + 1]) + r1 = int(unsafe_load(s.p, i1) + 1) + r2 = int(unsafe_load(s.p, i1 + 1)) (r1 > r2) && return zero(T) - r1 = int(searchsortedfirst(A.rowval0, i0 - 1, r1, r2, Base.Order.Forward)) - ((r1 > r2) || (A.rowval0[r1] + 1 != i0)) ? zero(T) : A.nzval[r1] + 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 -# Convert invalidates the input matrix -function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, A::Sparse{Tv,Ti}) - B = SparseMatrixCSC(A.c.m, A.c.n, increment!(A.colptr0), increment!(A.rowval0), A.nzval) - if A.c.stype != 0 - return isreal(A) ? Symmetric(B, A.c.stype < 0 ? :L : :U) : Hermitian(B, A.c.stype < 0 ? :L : :U) - end - return B -end -sparse!{Tv,Ti}(A::Sparse{Tv,Ti}) = convert(SparseMatrixCSC{Tv,Ti}, A) -sparse(A::Sparse) = sparse!(copy(A)) -sparse(L::Factor) = sparse!(Sparse(L)) -sparse(D::Dense) = sparse!(Sparse(D)) -function full(A::Sparse) - B = Dense(A).mat - if A.c.stype != 0 - return isreal(A) ? Symmetric(B, A.c.stype < 0 ? :L : :U) : Hermitian(B, A.c.stype < 0 ? :L : :U) - end - return B -end -full(A::Dense) = A.mat - ## Multiplication - -(*)(A::Sparse, B::Dense) = sdmult(A, false, 1., 0., B, zeros(size(A, 1), size(B, 2))) +(*)(A::Sparse, B::Sparse) = ssmult(A, B, 0, true, true) +(*)(A::Sparse, B::Dense) = sdmult!(A, false, 1., 0., B, zeros(size(A, 1), size(B, 2))) (*)(A::Sparse, B::VecOrMat) = (*)(A, Dense(B)) function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti}) - cm = common($Ti) + cm = common(Ti) if !is(A,B) - aa1 = transpose(B, 1) + aa1 = transpose(B, 2) ## result of ssmult will have stype==0, contain numerical values and be sorted return ssmult(A, aa1, 0, true, true) end @@ -878,61 +893,59 @@ function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti}) ## The A*A' case is handled by cholmod_aat. This routine requires ## A->stype == 0 (storage of upper and lower parts). If neccesary ## the matrix A is first converted to stype == 0 - if A.c.stype != 0 + s = unsafe_load(A.p) + if s.stype != 0 aa1 = copy(A, 0, 1) - return aat(aa1, Ti[], 1) + return aat(aa1, Ti[0:s.ncol-1], 1) else - return aat(A, Ti[], 1) + return aat(A, Ti[0:s.ncol-1], 1) end end -function Ac_mul_B{Tv<:VRealTypes}(A::Sparse{Tv}, B::Sparse{Tv}) - aa1 = transpose(A, 1) +function Ac_mul_B(A::Sparse, B::Sparse) + aa1 = transpose(A, 2) if is(A,B) return A_mul_Bc(aa1, aa1) end ## result of ssmult will have stype==0, contain numerical values and be sorted - return ssmult(aa1, B, 0, true, true, cm) + return ssmult(aa1, B, 0, true, true) end -A_mul_Bt{Tv<:VRealTypes}(A::Sparse{Tv}, B::Sparse{Tv}) = A_mul_Bc(A, B) # in the unlikely event of writing A*B.' instead of A*B' - -At_mul_B{Tv<:VRealTypes}(A::Sparse{Tv}, B::Sparse{Tv}) = Ac_mul_B(A,B) # in the unlikely event of writing A.'*B instead of A'*B - -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::Dense) = sdmult!(A, true, 1., 0., B, zeros(size(A, 2), size(B, 2))) Ac_mul_B(A::Sparse, B::VecOrMat) = Ac_mul_B(A, Dense(B)) ## Factorization methods - -analyze(A::SparseMatrixCSC) = analyze(Sparse(A)) - function cholfact(A::Sparse) + sA = unsafe_load(A.p) + sA.stype == 0 && throw(ArgumentError("sparse matrix is not symmetric/Hermitian")) cm = common(indtype(A)) ## may need to change final_asis as well as final_ll cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt F = analyze(A, cm) - factorize(A, F, cm) - F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) + factorize!(A, 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)) +function cholfact{Tv<:VTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, β::Tv) + cm = common(Ti) ## may need to change final_asis as well as final_ll - cm[common_final_ll] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt + cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt F = analyze(A, cm) - factorize(A, F, cm) + factorize_p!(A, β, Ti[], F, cm) + s = unsafe_load(F.p) + s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor)) return F end -function cholfact{Tv<:VTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, β::Tv) - cm = common(Ti) +function ldltfact(A::Sparse) + cm = common(indtype(A)) ## may need to change final_asis as well as final_ll - cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + cm[common_final_ll] = reinterpret(UInt8, [zero(Cint)]) # Hack! makes it a ldlt F = analyze(A, cm) - factorize_p(A, β, Ti[], F, cm) - F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) + factorize!(A, F, cm) return F end @@ -941,165 +954,132 @@ function ldltfact{Tv<:VTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, β::Tv) ## 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) - F.c.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(F.c.minor)) + 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!(L::Factor, A::Sparse) = factorize(A, L) -cholfact!{Tv<:VTypes,Ti<:ITypes}(L::Factor{Tv,Ti}, A::Sparse{Tv,Ti}, β::Tv) = factorize_p(A, β, Ti[], L) - cholfact{T<:VTypes}(A::Sparse{T}, β::Number) = cholfact(A, convert(T, β)) cholfact(A::SparseMatrixCSC) = cholfact(Sparse(A)) +cholfact(A::SparseMatrixCSC, β::Number) = cholfact(Sparse(A), β) cholfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = cholfact(Sparse(A)) +cholfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}, β::Number) = cholfact(Sparse(A), β) cholfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = cholfact(Sparse(A)) -cholfact!{T<:VTypes}(L::Factor{T}, A::Sparse{T}, β::Number) = cholfact!(L, A, convert(T, β)) -cholfact!{T<:VTypes}(L::Factor{T}, A::SparseMatrixCSC{T}, β::Number) = cholfact!(L, Sparse(A), convert(T, β)) -cholfact!{T<:VTypes}(L::Factor{T}, A::SparseMatrixCSC{T}) = cholfact!(L, Sparse(A)) +cholfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}, β::Number) = cholfact(Sparse(A), β) +ldltfact{T<:VTypes}(A::Sparse{T}, β::Number) = ldltfact(A, convert(T, β)) ldltfact(A::SparseMatrixCSC) = ldltfact(Sparse(A)) +ldltfact(A::SparseMatrixCSC, β::Number) = ldltfact(Sparse(A), β) ldltfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) = ldltfact(Sparse(A)) +ldltfact{Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}, β::Number) = ldltfact(Sparse(A), β) ldltfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}) = ldltfact(Sparse(A)) +ldltfact{Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}, β::Number) = ldltfact(Sparse(A), β) + +function update!{Tv<:VTypes,Ti<:ITypes}(F::Factor{Tv,Ti}, A::Sparse{Tv,Ti}) + cm = common(Ti) + s = unsafe_load(F.p) + if bool(s.is_ll) + cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + end + factorize!(A, F, cm) +end +function update!{Tv<:VTypes,Ti<:ITypes}(F::Factor{Tv,Ti}, A::Sparse{Tv,Ti}, β::Tv) + cm = common(Ti) + s = unsafe_load(F.p) + if bool(s.is_ll) + cm[common_final_ll] = reinterpret(UInt8, [one(Cint)]) # Hack! makes it a llt + end + factorize_p!(A, β, Ti[0:size(F, 1) - 1], F, cm) +end +update!{T<:VTypes}(F::Factor{T}, A::SparseMatrixCSC{T}) = update!(F, Sparse(A)) +update!{T<:VTypes}(F::Factor{T}, A::SparseMatrixCSC{T}, β::Number) = update!(F, Sparse(A), convert(T, β)) ## Solvers (\)(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B) -(\)(L::Factor, b::Vector) = reshape(solve(CHOLMOD_A, L, Dense(b)).mat, length(b)) -(\)(L::Factor, B::Matrix) = solve(CHOLMOD_A, L, Dense(B)).mat +(\)(L::Factor, b::Vector) = reshape(convert(Matrix, solve(CHOLMOD_A, L, Dense(b))), length(b)) +(\)(L::Factor, B::Matrix) = convert(Matrix, solve(CHOLMOD_A, L, Dense(B))) (\)(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B) # When right hand side is sparse, we have to ensure that the rhs is not marked as symmetric. -(\)(L::Factor, B::SparseMatrixCSC) = sparse!(spsolve(CHOLMOD_A, L, Sparse(B, 0))) +(\)(L::Factor, B::SparseMatrixCSC) = sparse(spsolve(CHOLMOD_A, L, Sparse(B, 0))) Ac_ldiv_B(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B) -Ac_ldiv_B(L::Factor, B::VecOrMat) = solve(CHOLMOD_A, L, Dense(B)).mat +Ac_ldiv_B(L::Factor, B::VecOrMat) = convert(Matrix, solve(CHOLMOD_A, L, Dense(B))) Ac_ldiv_B(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B) -Ac_ldiv_B(L::Factor, B::SparseMatrixCSC) = sparse!(CHOLMOD_A, spsolve(L, Sparse(B))) - -speye(m::Integer, n::Integer) = speye(m, n, Float64) # default element type is Float64 -speye(n::Integer) = speye(n, n, Float64) # default shape is square - -spzeros(m::Integer, n::Integer, nzmax::Integer) = spzeros(m, n, nzmax, 1.) - -function norm(A::Sparse, p = 1) - if p == 1 - return norm_sparse(A, 1) - elseif p == Inf - return norm_sparse(A, 0) - else - throw(ArgumentError("only 1 and inf norms are supported")) - end -end - -## Julia methods (non-CHOLMOD) -function findnz{Tv,Ti}(A::Sparse{Tv,Ti}) - jj = similar(A.rowval0) # expand A.colptr0 to a vector of indices - for j in 1:A.c.n, k in (A.colptr0[j]+1):A.colptr0[j+1] - jj[k] = j - end - - ind = similar(A.rowval0) - ipos = 1 - count = 0 - for k in 1:length(A.nzval) - if A.nzval[k] != 0 - ind[ipos] = k - ipos += 1 - count += 1 - else - println("warning: sparse matrix contains explicitly stored zeros") - end - end - ind = ind[1:count] # ind is the indices of nonzeros in A.nzval - (increment!(A.rowval0[ind]), jj[ind], A.nzval[ind]) -end - -findnz(L::Factor) = findnz(Sparse(L)) - -function diag{Tv}(A::Sparse{Tv}) - minmn = minimum(size(A)) - res = zeros(Tv,minmn) - cp0 = A.colptr0 - rv0 = A.rowval0 - anz = A.nzval - for j in 1:minmn, k in (cp0[j]+1):cp0[j+1] - if rv0[k] == j-1 - res[j] += anz[k] - end - end - res -end - -function diag{Tv}(L::Factor{Tv}) - res = zeros(Tv,L.c.n) - xv = L.x - if bool(L.c.is_super) - nec = decrement!(diff(L.super)) # number of excess columns per supernode - dstride = increment!(diff(L.pi)) # stride of diagonal elements (nrow + 1) - px = L.px +Ac_ldiv_B(L::Factor, B::SparseMatrixCSC) = Ac_ldiv_B(L, Sparse(B)) + +## Other convenience methods +function diag{Tv}(F::Factor{Tv}) + f = unsafe_load(F.p) + res = Base.zeros(Tv, int(f.n)) + xv = f.x + if bool(f.is_super) + px = f.px pos = 1 - for i in 1:length(nec) - base = px[i] + 1 - res[pos] = xv[base] + for i in 1:f.nsuper + base = unsafe_load(px, i) + 1 + res[pos] = unsafe_load(xv, base) pos += 1 - for j in 1:nec[i] - res[pos] = xv[base + j*dstride[i]] + for j in 1:unsafe_load(f.super, i + 1) - unsafe_load(f.super, i) - 1 + res[pos] = unsafe_load(xv, base + j*(unsafe_load(f.pi, i + 1) - unsafe_load(f.pi, i) + 1)) pos += 1 end end else - c0 = L.p - r0 = L.i - xv = L.x - for j in 1:L.c.n - jj = c0[j]+1 - assert(r0[jj] == j-1) - res[j] = xv[jj] + c0 = f.p + r0 = f.i + xv = f.x + for j in 1:f.n + jj = unsafe_load(c0, j) + 1 + assert(unsafe_load(r0, jj) == j - 1) + res[j] = unsafe_load(xv, jj) end end res end -function logdet{Tv,Ti}(L::Factor{Tv,Ti},sub) +function logdet{Tv<:VTypes,Ti<:ITypes}(F::Factor{Tv,Ti}) + f = unsafe_load(F.p) res = zero(Tv) - for d in diag(L)[sub] res += log(abs(d)) end - bool(L.c.is_ll) ? 2res : res + for d in diag(F) res += log(abs(d)) end + bool(f.is_ll) ? 2res : res end -logdet(L::Factor) = logdet(L, 1:L.c.n) + det(L::Factor) = exp(logdet(L)) -isposdef!{Tv<:VTypes,Ti}(A::SparseMatrixCSC{Tv,Ti}) = ishermitian(A) && cholfact(A).c.minor == size(A,1) +function isposdef{Tv<:VTypes,Ti}(A::SparseMatrixCSC{Tv,Ti}) + if !ishermitian(A) + return false + end + f = cholfact(A) + s = unsafe_load(f.p) + return s.minor == size(A,1) +end function issym(A::Sparse) - if isreal(A) - if A.c.stype != 0 - return true - else - i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 - return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG - end - else - if A.c.stype != 0 - return false - else - i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 - return i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG - end + s = unsafe_load(A.p) + if s.stype != 0 + return isreal(A) end + i = symmetry(A, ifelse(version >= v"3.0.5", 0, 1))[1] # 0 is faster, but had a bug before 3.0.5 + return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG end function ishermitian(A::Sparse) + s = unsafe_load(A.p) if isreal(A) - if A.c.stype != 0 + if s.stype != 0 return true else - i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 + i = symmetry(A, ifelse(version >= v"3.0.5", 0, 1))[1] # 0 is faster, but had a bug before 3.0.5 return i == MM_SYMMETRIC || i == MM_SYMMETRIC_POSDIAG end else - if A.c.stype != 0 + if s.stype != 0 return true else - i = symmetry(A, ifelse(version >= v"4.4.3", 0, 1))[1] # 0 is faster, but had a bug before 4.4.3 + i = symmetry(A, ifelse(version >= v"3.0.5", 0, 1))[1] # 0 is faster, but had a bug before 3.0.5 return i == MM_HERMITIAN || i == MM_HERMITIAN_POSDIAG end end diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index f4ff95d67b9d2..244719a444f2e 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -33,10 +33,11 @@ using Base.SparseMatrix.CHOLMOD ## residual 1.3e-19 (|Ax-b|/(|A||x|+|b|)) after iterative refinement ## rcond 9.5e-06 -A = Sparse(int32([0,1,2,3,6,9,12,15,18,20,25,30,34,36,39,43,47,52,58,62,67,71,77,84,90,93,95, +A = CHOLMOD.Sparse(48, 48, + int32([0,1,2,3,6,9,12,15,18,20,25,30,34,36,39,43,47,52,58,62,67,71,77,84,90,93,95, 98,103,106,110,115,119,123,130,136,142,146,150,155,161,167,174,182,189,197, 207,215,224]), # zero-based column pointers - int32([0,1,2,1,2,3,0,2,4,0,1,5,0,4,6,1,3,7,2,8,1,3,7,8,9,0,4,6,8,10,5,6,7,11,6,12, + int32([0,1,2,1,2,3,0,2,4,0,1,5,0,4,6,1,3,7,2,8,1,3,7,8,9,0,4,6,8,10,5,6,7,11,6,12, 7,11,13,8,10,13,14,9,13,14,15,8,10,12,14,16,7,11,12,13,16,17,0,12,16,18,1, 5,13,15,19,2,4,14,20,3,13,15,19,20,21,2,4,12,16,18,20,22,1,5,17,18,19,23,0, 5,24,1,25,2,3,26,2,3,25,26,27,4,24,28,0,5,24,29,6,11,24,28,30,7,25,27,31,8, @@ -85,27 +86,29 @@ A = Sparse(int32([0,1,2,3,6,9,12,15,18,20,25,30,34,36,39,43,47,52,58,62,67,71,77 5.0e8,-4.79857058902e6,134990.2747,2.47238730198e9,-1.14928529451e6, 2.29724661236e8,-5.57173510779e7,-833333.333333,-1.25e6,2.5e8,2.39928529451e6, 9.61679848804e8,275828.470683,-5.57173510779e7,1.09411960038e7,2.08333333333e6, - 1.0e8,-2.5e6,140838.195984,-1.09779731332e8,5.31278103775e8], 48, 48, 1) -@test_approx_eq norm(A, Inf) 3.570948074697437e9 -@test_approx_eq norm(A) 3.570948074697437e9 + 1.0e8,-2.5e6,140838.195984,-1.09779731332e8,5.31278103775e8], 1) +@test_approx_eq CHOLMOD.norm_sparse(A, 0) 3.570948074697437e9 +@test_approx_eq CHOLMOD.norm_sparse(A, 1) 3.570948074697437e9 +@test_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.mat ones(size(x)) +@test_approx_eq x ones(size(x)) chma = cholfact(A) # LL' form @test isvalid(chma) x = chma\B -@test_approx_eq x.mat ones(size(x)) +@test_approx_eq x ones(size(x)) #lp_afiro example -afiro = Sparse(int32([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,25,27,29,33,37, +afiro = CHOLMOD.Sparse(27, 51, + int32([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,25,27,29,33,37, 41,45,47,49,51,53,55,57,59,63,65,67,69,71,75,79,83,87,89,91,93,95,97, 99,101,102]), - int32([2,3,6,7,8,9,12,13,16,17,18,19,20,21,22,23,24,25,26,0,1,2,23,0,3,0,21, + int32([2,3,6,7,8,9,12,13,16,17,18,19,20,21,22,23,24,25,26,0,1,2,23,0,3,0,21, 1,25,4,5,6,24,4,5,7,24,4,5,8,24,4,5,9,24,6,20,7,20,8,20,9,20,3,4,4,22, 5,26,10,11,12,21,10,13,10,23,10,20,11,25,14,15,16,22,14,15,17,22,14, 15,18,22,14,15,19,22,16,20,17,20,18,20,19,20,13,15,15,24,14,26,15]), @@ -115,13 +118,15 @@ afiro = Sparse(int32([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,25,27 2.386,-1.0,2.408,-1.0,2.429,1.4,1.0,1.0,-1.0,1.0,1.0,-1.0,-0.43,1.0,0.109, 1.0,-1.0,1.0,-1.0,1.0,-1.0,1.0,1.0,-0.43,1.0,1.0,0.109,-0.43,1.0,1.0,0.108, -0.39,1.0,1.0,0.108,-0.37,1.0,1.0,0.107,-1.0,2.191,-1.0,2.219,-1.0,2.249, - -1.0,2.279,1.4,-1.0,1.0,-1.0,1.0,1.0,1.0], 27, 51, 0) -chmaf = cholfact(afiro) + -1.0,2.279,1.4,-1.0,1.0,-1.0,1.0,1.0,1.0], 0) +afiro2 = CHOLMOD.aat(afiro, Int32[0:50], int32(1)) +CHOLMOD.change_stype!(afiro2, -1) +chmaf = cholfact(afiro2) y = afiro'*ones(size(afiro,1)) sol = chmaf\(afiro*y) # least squares solution @test isvalid(sol) pred = afiro'*sol -@test norm(afiro * (y.mat - pred.mat)) < 1e-8 +@test norm(afiro * (convert(Matrix, y) - convert(Matrix, pred))) < 1e-8 let # Issue 9160 @@ -130,7 +135,7 @@ let # Issue 9160 A = sprand(10,10,0.1) A = convert(SparseMatrixCSC{elty,Ti},A) - cmA = CHOLMOD.parse(A) + cmA = CHOLMOD.Sparse(A) B = sprand(10,10,0.1) B = convert(SparseMatrixCSC{elty,Ti},B) @@ -156,6 +161,7 @@ let # Issue 9160 end end + # Issue #9915 @test speye(2)\speye(2) == eye(2) @@ -172,11 +178,255 @@ ACSC = sprandn(10, 10, 0.3) + I @test ishermitian(Sparse(Hermitian(complex(ACSC), :U))) # test Sparse constructor for c_Sparse{Tv,Ti} input( and Sparse*Sparse) -B = Sparse(SparseMatrixCSC{Float64,Int32}(sprandn(48, 48, 0.1))) # A has Int32 indeces +B = CHOLMOD.Sparse(SparseMatrixCSC{Float64,Int32}(sprandn(48, 48, 0.1))) # A has Int32 indeces @test_approx_eq A*B sparse(A)*sparse(B) # test Sparse constructor for c_SparseVoid (and read_sparse) writedlm("tmp.mtx", ["%%MatrixMarket matrix coordinate real symmetric","3 3 4","1 1 1","2 2 1","3 2 0.5","3 3 1"]) -@test full(Sparse("tmp.mtx")) == [1 0 0;0 1 0.5;0 0.5 1] +@test 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 CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, 2.0)) == CHOLMOD.Sparse(F) + + F = ldltfact(A1pd) + @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) + @test CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd)) == CHOLMOD.Sparse(F) + + F = ldltfact(A1pdSparse, 2) + @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) + @test CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, 2.0)) == CHOLMOD.Sparse(F) + + @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 From fce7b10aa4bd2352f29f2082c55e8f740d0f4fec Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 8 Feb 2015 09:23:28 -0500 Subject: [PATCH 08/10] Remove Common struct and move more of the setup to cholmod_h.jl --- base/sparse/cholmod.jl | 48 +++------------------------------------- base/sparse/cholmod_h.jl | 13 +++++++++++ test/sparse/cholmod.jl | 2 +- 3 files changed, 17 insertions(+), 46 deletions(-) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index afbc94c31387d..19ea4bd9d0c8c 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -21,23 +21,6 @@ using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, incre include("cholmod_h.jl") -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,()) - -type CHOLMODException <: Exception - msg::AbstractString -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 @@ -53,35 +36,10 @@ for Ti in IndexTypes end end -### A way of examining some of the fields in cholmod_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 Common - 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 -end - -### These offsets should be reconfigured to be less error-prone in matches -const cholmod_com_offsets = Array(Int, length(Common.types)) +### 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{Int},), cholmod_com_offsets) + 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] diff --git a/base/sparse/cholmod_h.jl b/base/sparse/cholmod_h.jl index 71d6ee96efa4a..136ade6d9bb90 100644 --- a/base/sparse/cholmod_h.jl +++ b/base/sparse/cholmod_h.jl @@ -67,3 +67,16 @@ 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) +ccall((:cholmod_version, :libcholmod), Cint, (Ptr{Cint},), version_array) +const version = VersionNumber(version_array...) +const cholmod_com_sz = ccall((:jl_cholmod_common_size,:libsuitesparse_wrapper),Int,()) diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index 244719a444f2e..f845da04def5a 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -32,7 +32,7 @@ using Base.SparseMatrix.CHOLMOD ## 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 - +println("version: ", Base.SparseMatrix.CHOLMOD.version) 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, From d3c4c8f4feebbae91e66c5c0a694942c36754ee1 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 8 Feb 2015 11:57:13 -0500 Subject: [PATCH 09/10] Fix checks for symmetry in Sparse constructors --- base/sparse/cholmod.jl | 21 +++++++++++++++------ base/sparse/cholmod_h.jl | 6 +++++- test/sparse/cholmod.jl | 8 ++++---- 3 files changed, 24 insertions(+), 11 deletions(-) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 19ea4bd9d0c8c..6ff6361b71d97 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -53,7 +53,10 @@ end # Type definitions # #################### -# The three core data types for CHOLMOD: Dense, Sparse and Factor. CHOLMOD manages the memory, so the Julia versions only wrap a pointer to a struct. Therefore finalizers should be registret each time a pointer is returned from CHOLMOD. +# The three core data types for CHOLMOD: Dense, Sparse and Factor. CHOLMOD +# manages the memory, so the Julia versions only wrap a pointer to a struct. +# Therefore finalizers should be registret each time a pointer is returned from +# CHOLMOD. # Dense immutable C_Dense{T<:VTypes} @@ -89,7 +92,9 @@ immutable C_Sparse{Tv<:VTypes,Ti<:ITypes} packed::Cint end -# Corresponds to the exact definition of cholmod_sparse_struct in the library. Useful when reading matrices of unknown type from files as in cholmod_read_sparse +# 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 @@ -187,7 +192,8 @@ end ################# # 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 +## 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}) @@ -527,7 +533,8 @@ for Ti in IndexTypes end # 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 + # 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}), @@ -698,7 +705,8 @@ function Sparse(p::Ptr{C_SparseVoid}) end -Sparse(A::Dense) = dense_to_sparse(A, Cint) # default is Cint which is probably sufficient when converted from dense matrix +# 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) @@ -761,7 +769,8 @@ 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 +# 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) diff --git a/base/sparse/cholmod_h.jl b/base/sparse/cholmod_h.jl index 136ade6d9bb90..f97dff5effd6e 100644 --- a/base/sparse/cholmod_h.jl +++ b/base/sparse/cholmod_h.jl @@ -77,6 +77,10 @@ macro isok(A) end const version_array = Array(Cint, 3) -ccall((:cholmod_version, :libcholmod), Cint, (Ptr{Cint},), version_array) +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/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index f845da04def5a..05acf09f1c62d 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -32,7 +32,7 @@ using Base.SparseMatrix.CHOLMOD ## 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 -println("version: ", Base.SparseMatrix.CHOLMOD.version) + 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, @@ -382,15 +382,15 @@ for elty in (Float64, Complex{Float64}) F = cholfact(A1pdSparse, 2) @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) - @test CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, 2.0)) == CHOLMOD.Sparse(F) + @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 CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd)) == CHOLMOD.Sparse(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 F = ldltfact(A1pdSparse, 2) @test isa(CHOLMOD.Sparse(F), CHOLMOD.Sparse{elty, CHOLMOD.SuiteSparse_long}) - @test CHOLMOD.Sparse(CHOLMOD.update!(copy(F), A1pd, 2.0)) == CHOLMOD.Sparse(F) + @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) From c380d394a38d6cbaed6fcc3592c5f1d7277f8fb4 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Wed, 11 Feb 2015 14:23:39 -0500 Subject: [PATCH 10/10] Adapt to new concatenation behavior --- base/sparse/cholmod.jl | 6 +++--- test/sparse/cholmod.jl | 16 ++++++++-------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 6ff6361b71d97..0e3e3d648f857 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -863,9 +863,9 @@ function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti}) s = unsafe_load(A.p) if s.stype != 0 aa1 = copy(A, 0, 1) - return aat(aa1, Ti[0:s.ncol-1], 1) + return aat(aa1, Ti[0:s.ncol-1;], 1) else - return aat(A, Ti[0:s.ncol-1], 1) + return aat(A, Ti[0:s.ncol-1;], 1) end end @@ -957,7 +957,7 @@ function update!{Tv<:VTypes,Ti<:ITypes}(F::Factor{Tv,Ti}, A::Sparse{Tv,Ti}, β:: 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) + 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, β)) diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index 05acf09f1c62d..947b4e4529c9e 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -119,7 +119,7 @@ afiro = CHOLMOD.Sparse(27, 51, 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)) +afiro2 = CHOLMOD.aat(afiro, Int32[0:50;], int32(1)) CHOLMOD.change_stype!(afiro2, -1) chmaf = cholfact(afiro2) y = afiro'*ones(size(afiro,1)) @@ -287,8 +287,8 @@ p = ccall((:cholmod_allocate_sparse, :libcholmod), Ptr{CHOLMOD.C_Sparse{Float64, @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))) + 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) @@ -408,11 +408,11 @@ for elty in (Float64, Complex{Float64}) 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_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_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, 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 @@ -423,8 +423,8 @@ for elty in (Float64, Complex{Float64}) 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_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