Lufact SparseMatrixCSC Int32

Specific Domains Numerics
,
1

Hi everyone,

I am very new to Julia, so please let me know if there is info in the documentation that would have avoided the issue. I use the following version:

Julia Version 0.6.1
Commit 0d7248e (2017-10-24 22:15 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, ivybridge)

and get the following error (minimal working example below):

ERROR: LoadError: MethodError: no method matching lufact!(::Hermitian{Float64,SparseMatrixCSC{Float64,Int32}}, ::Type{Val{true}})

which is due to the sparse matrix using Int32 for indexing. If the indices are Int64 everything works just fine. I don’t know if this is a bug or wanted behavior. However, the error message didn’t help me much and I needed quite some time to figure the issue out. I didn’t open an issue on github, because I wasn’t sure if it is appropriate.

Cheers,
Lukas

Code

# Setup (indices, values and r.h.s.)
i1 = [1, 2, 3, 2, 1]
j1 = [1, 2, 3, 1, 2]

i2 = convert(Array{Int32, 1}, i1)
j2 = convert(Array{Int32, 1}, j1)

v = [1.0, 2.0, 3.0, 0.5, 0.5]
b = [1.0, 2.0, 3.0]

# Working Version
A = sparse(i1, j1, v)
println("A has type: ", typeof(A))
x = A \ b
println("x = ", x)

# Non-Working version
A2 = sparse(i2, j2, v)
println("A2 has type: ", typeof(A2))
x2 = A2 \ b
println("x2 = ", x2)

Full Output:

A has type: SparseMatrixCSC{Float64,Int64}
x = [0.571429, 0.857143, 1.0]
A2 has type: SparseMatrixCSC{Float64,Int32}

ERROR: LoadError: MethodError: no method matching lufact!(::Hermitian{Float64,SparseMatrixCSC{Float64,Int32}}, ::Type{Val{true}})
Closest candidates are:
lufact!(!Matched::Union{Base.ReshapedArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T}, DenseArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2}, SubArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T} where N where T, DenseArray}, ::Union{Type{Val{false}}, Type{Val{true}}}) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at linalg/lu.jl:16
lufact!(!Matched::Union{Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T}, DenseArray{T,2}, SubArray{T,2,A,I,L} where L} where I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex, Int64, Range{Int64}},N} where N} where A<:Union{Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{DenseArray, SubArray{T,N,P,I,true} where I<:Tuple{Union{Base.Slice, UnitRange},Vararg{Any,N} where N} where P where N where T} where N where T, DenseArray} where T, ::Union{Type{Val{false}}, Type{Val{true}}}) at linalg/lu.jl:31
lufact!(!Matched::Tridiagonal{T}, ::Union{Type{Val{false}}, Type{Val{true}}}) where T at linalg/lu.jl:321
Stacktrace:
[1] (::Hermitian{Float64,SparseMatrixCSC{Float64,Int32}}, ::Array{Float64,1}) at ./linalg/generic.jl:832
[2] (::SparseMatrixCSC{Float64,Int32}, ::Array{Float64,1}) at ./sparse/linalg.jl:869
[3] include_from_node1(::String) at ./loading.jl:576
[4] include(::String) at ./sysimg.jl:14
[5] process_options(::Base.JLOptions) at ./client.jl:305
[6] _start() at ./client.jl:371
while loading /Users/[private]/sparse_bug_example.jl, in expression starting on line 20

2

x2 = lufact(A2) \ b works so it looks like a missing method, I can take a closer look later. Also, if you put your code inside triple backticks (e.g. ``` code ```) it will format nicer.

4 Likes
3

I’m also having this issue in 0.6.2.

4

What is the preferred solution here? Automatic conversion to Int64 or a better error message?

5

IMO, it should be possible to perform the sparse LU in single precision if the user requests it, without conversion.

1 Like
6

IMO, it should be possible to perform the sparse LU in single precision if the user requests it, without conversion.

It also is. The issue here is A\b tries to Cholesky factorize the matrix because it is symmetric and because of the global common struct used by CHOLMOD we don’t support more than a single integer size for the indices. However, our UMFPACK bindings support Int32 indices so we could fall back on using the LU when the indices are Int32 even though the matrix is symmetric.

7

I wonder if it’s wise to use Cholesky based on numerical symmetry. It may not be what the user wants. Why not let the user specify if their matrix is Symmetric or Hermitian?

8

I agree on this last point

9

Could you elaborate? What is the concern? The user can always ask for a specific factorization, e.g. lufact(A)\b.

1 Like
10

I looked and cholmod_symmetry() checks for actual equality between elements of the lower and upper triangle to determine symmetry. Initially I thought there might be a tolerance, which got me worried, but there isn’t. Still, the O(nnz) cost incurred by this check comes as a surprise to the user and could be quite high for large problems. Perhaps that’s ok for something as generic as sparse backslash. I can’t help but feel cholfact should only accept Symmetric matrices though.