SparseArrays.jl

SparseArrays.jl copied to clipboard

The reason seems to be that in the 3-arg ldiv! its 2-arg-version is called, which is missing for sparse Cholesky factorizations.

MWE:

using LinearAlgebra, SparseArrays
A = sprand(100, 100, 0.01) + spdiagm(0 => ones(100))
x = rand(100)
y = similar(x)
ldiv!(y, factorize(A'A), x)

throws

ERROR: MethodError: no method matching ldiv!(::SuiteSparse.CHOLMOD.Factor{Float64}, ::Array{Float64,1})

Asking for the employed method gives

julia> @which ldiv!(y, factorize(A'A), x)
ldiv!(Y::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T, A::Factorization, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) in LinearAlgebra at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/factorization.jl:100

Sorry, I messed something up. I guess this functionality never existed anyway. I wonder how one could avoid the confusion. The docstring of ldiv! emphasizes that the middle argument should be a factorization, but doesn't mention any restriction on the type of factorization. And the exact same call works when A is dense and s.p.d.

I think your first reaction is right. I think the solver code in the CHOLMOD needs adjustment/renaming such that three argument ldiv! becomes supported. Actually, I though it already was. For now, you can work around the issue by explicitly calling lu instead of factorize.

In a sense, you can call the (generic) three argument ldiv! (from LinearAlgebra/src/factorization.jl), it's just that it calls the two argument ldiv! internally, which does not exist for the sparse factorization. So one may need to provide a two-arg ldiv! for sparse factorizations, or add a more specific 3-arg ldiv! with higher precedence. Is the latter what you have in mind?

It was the latter I had in mind.

I started working on this. I realized that both QRsparse and CHOLMOD are missing 3-args ldiv! methods. From what I understand, there is no truly in-place solve method for these two, so pretending there was one by defining, say,

ldiv!(Y, F::Union{SparseQR,CHOLMOD}, X) = copyto!(Y, F\X)

is then probably misleading, because there will be hidden, potentially significant allocations, which the user would not expect. If we choose to not introduce something like the above, maybe we should then at least introduce a more informative error message?

ldiv!(Y, F::Union{SparseQR,CHOLMOD}, X) = error("There is no `ldiv!` in-place solver for the $F factorization. Use \\ instead.")

The error message mentioned in the OP says there is no 2-args ldiv!, but the user called a 3-args ldiv!.

I have the same error with:

hDL = ldlt([mP + (σ * sparse(I, numElementsX, numElementsX)) transpose(mA); mA -ρ¹ * sparse(I, numRowsA, numRowsA)]);
ldiv!(hDL, vV);

Is there any intention to have in place version or is it something which infeasible?
@dkarrasch , I get and agree with your point that a user which uses a ! function expect no allocations. But is it really the case on the Julia eco system? I think it's not.

The real thing is to offer "Buffers" in some functions to indeed prevent any allocation of the function itself.