Julia: Unable to solve asymmetric sparse system with sparse rhs

[jecs]

I'm trying to solve a sparse system of equations A\b, where

A::SparseMatrixCSC{Complex{Float64},Int64}
b::SparseMatrixCSC{Complex{Float64},Int64}

I'm getting the following error:

ERROR: LoadError: MethodError: `A_ldiv_B!` has no method matching A_ldiv_B!(::Base.SparseMatrix.UMFPACK.UmfpackLU{Complex{Float64},Int64}, ::SparseMatrixCSC{Complex{Float64},Int64})
Closest candidates are:
  A_ldiv_B!{T}(::Base.LinAlg.LU{T,Tridiagonal{T}}, ::Union{AbstractArray{T,1},AbstractArray{T,2}})
  A_ldiv_B!{T}(::Diagonal{T}, ::AbstractArray{T,2})
  A_ldiv_B!(::Union{Base.LinAlg.AbstractTriangular{T,S<:AbstractArray{T,2}},Bidiagonal{T}}, ::AbstractArray{T,2})
  ...
 in \ at ./linalg/factorization.jl:25
....

Version Info:

Julia Version 0.4.7
Commit ae26b25 (2016-09-18 16:17 UTC)
Platform Info:
  System: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU           E5430  @ 2.66GHz
  WORD_SIZE: 64
  BLAS: libopenblas (NO_LAPACKE DYNAMIC_ARCH NO_AFFINITY Penryn)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.8

[jecs]

Update: Interestingly, I can solve the system when I cast b to Array{Complex{64},1}.

[andreasnoack]

It is way more common to solve for a dense vector. In many cases, the solution will also be dense even though both A and b are sparse. In 0.5 this works for symmetric and Hermitian matrices but not in the general case. This is because the underlying library supports this in the symmetric/Hermitian case (CHOLMOD) but not in the general case (UMFPACK) so I don't think we'll be able to support this any time soon.

Update: A solution could be to automatically promote the sparse vector to a dense.

[nalimilan]

So just print an error message explaining this?

[tkelman]

Didn't sparse triangular solves with sparse rhs get implemented at some point, probably by @Sacha0? Maybe we could get whatever information is needed out of UMFPACK's data structure to use the Julia solve routines?

[Sacha0]

Didn't sparse triangular solves with sparse rhs get implemented at some point, probably by @Sacha0?

Are you thinking of #14447? (That code does not perform proper sparsetriangular-sparse solves.) Best!

[fredrikekre]

Promote SparseVector to dense and throw a more descriptive error for SparseMatrix rhs, since promoting that could potentially be a memory problem?

[ViralBShah]

Works.

julia> x = ones(10) + im*ones(10);

julia> sparse(Diagonal(x)) \ x
10-element Array{Complex{Float64},1}:
 1.0 - 0.0im
 1.0 - 0.0im
 1.0 - 0.0im
 1.0 - 0.0im
 1.0 - 0.0im
 1.0 - 0.0im
 1.0 - 0.0im
 1.0 - 0.0im
 1.0 - 0.0im
 1.0 - 0.0im

[jecs]

@ViralBShah I'm afraid that you did not properly test the issue. This pertains to sparse b and asymmetric and sparse A.

julia> b=sparse((1+im)ones(10,1));

julia> A=Diagonal(b[:]);

julia> A=A+√(eps())*sparse(max.(randn(size(A)),0.));

julia> typeof(b)
SparseMatrixCSC{Complex{Float64},Int64}

julia> typeof(A)
SparseMatrixCSC{Complex{Float64},Int64}

julia> issymmetric(A)
false

julia> A\b
ERROR: MethodError: no method matching ldiv!(::SuiteSparse.UMFPACK.UmfpackLU{Complex{Float64},Int64}, ::SparseMatrixCSC{Complex{Float64},Int64})

[ViralBShah]

Thank you for pointing that out. I did indeed test the wrong thing. The easy thing to do is to add a dispatch case to convert the rhs to dense.

[thudso]

While I appreciate that there is an easy workaround to this issue, it would be great to have this fixed or even to have a warning message implemented: the current error is not particularly transparent.

[ViralBShah]

Closing in favour of #34052