Julia: Solving with sparse cholesky returns two-dimensional array

[clason]

In v1.0, using the Cholesky decomposition of a sparse array to solve a linear system for the transpose results in a two-dimensional array instead of a vector:

using LinearAlgebra,SparseArrays
N = 4
A = sparse(I,N,N)
C = cholesky(A)
b = randn(N)
C'\b

returns

4×1 Array{Float64,2}:
...

By contrast, C\b (correctly) returns

4-element Array{Float64,1}:
...

(The corresponding test for Array(I,N,N) throws a ERROR: MethodError: no method matching adjoint(::Cholesky{Float64,Array{Float64,2}}).)

[StefanKarpinski]

@Sacha0, @andreasnoack?

[andreasnoack]

This is a bug and should be fixed.

[ViralBShah]

Yes, we must fix this. The culprit is probably in here:

\(adjL::Adjoint{<:Any,<:Factor}, B::Dense) = (L = adjL.parent; solve(CHOLMOD_A, L, B))
\(adjL::Adjoint{<:Any,<:Factor}, B::VecOrMat) = (L = adjL.parent; Matrix(solve(CHOLMOD_A, L, Dense(B))))

[raghav9-97]

I would like to work on this if nobody is currently working on it.

[ViralBShah]

Go for it.

[raghav9-97]

I think the problem is explicit conversion to Matrix in \(adjL::Adjoint{<:Any,<:Factor}, B::VecOrMat) = (L = adjL.parent; Matrix(solve(CHOLMOD_A, L, Dense(B))))
Is that correct?

[ViralBShah]

I believe the matrix and the vector cases need to be separated out.

[raghav9-97]

\(adjL::Adjoint{<:Any,<:Factor}, B::Vector) = (L = adjL.parent; Vector(solve(CHOLMOD_A, L, Dense(B))))
\(adjL::Adjoint{<:Any,<:Factor}, B::Matrix) = (L = adjL.parent; Matrix(solve(CHOLMOD_A, L, Dense(B))))

Something like this??

[ViralBShah]

Probably. It needs to be tried out and all tests should pass.

[andreasnoack]

Something like this??

Yes that should work. While you are at it, please allow that B can be StridedVector/StridedMatrix. Also, it's slightly nicer to use lower case for vectors (but still upper case for matrices). Finally, please also split the function over multiple lines instead of using ;.