SparseArrays.jl

SparseArrays.jl copied to clipboard

The following errors for cholesky factorizations:

julia> using LinearAlgebra, SparseArrays
julia> A = cholesky(sparse(Matrix(I,10,10))); # sparse cholesky factorization with real elements
julia> z = rand(ComplexF64, 10); # complex RHS vector
julia> A\@views(z[:])
ERROR: InexactError: Float64(0.868342583034311 + 0.07242770884885896im)
Stacktrace:
 [1] Type at ./complex.jl:37 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] unsafe_store! at ./pointer.jl:118 [inlined]
 [4] SuiteSparse.CHOLMOD.Dense{Float64}(::SubArray{Complex{Float64},1,Array{Complex{Float64},1},Tuple{Base.Slice{Base.OneTo{Int64}}},true}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/SuiteSparse/src/cholmod.jl:802
 [5] \(::SuiteSparse.CHOLMOD.Factor{Float64}, ::SubArray{Complex{Float64},1,Array{Complex{Float64},1},Tuple{Base.Slice{Base.OneTo{Int64}}},true}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/SuiteSparse/src/cholmod.jl:1682
 [6] top-level scope at none:0

while the same operation works without the view:

julia> isapprox(A\z, z)
true

I've chased the issue down to being a simple missing definition in SuiteSparse, where we have the following (in cholmod.jl):

(\)(L::Factor{T}, B::Vector{Complex{T}}) where {T<:Float64} = complex.(L\real(B), L\imag(B))
(\)(L::Factor{T}, B::Matrix{Complex{T}}) where {T<:Float64} = complex.(L\real(B), L\imag(B))
(\)(L::Factor{T}, b::StridedVector) where {T<:VTypes} = Vector(L\Dense{T}(b))
(\)(L::Factor{T}, B::StridedMatrix) where {T<:VTypes} = Matrix(L\Dense{T}(B))

So, the analogous definitions for complex StridedVectors and StridedMatrixs are missing. I can submit the PR to base myself, but I wanted to open an issue first in case there is a reason they aren't defined.

Also, this does seem to work for other factorizations, e.g. lu:

julia> B = lu(sparse(Matrix(I,10,10)));
julia> isapprox(B\@views(z[:]), z)
true

My julia version info is following, in case it's relevant:

julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: AMD Phenom(tm) II X6 1100T Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, amdfam10)
Environment:
  JULIA_NUM_THREADS = 24
  JULIA_EDITOR = code

Nice to have a PR here.

I believe CHOLMOD can handle strides, I'll look