SparseArrays: Speed up right-division by diagonal matrices

[eschnett]

No description provided.

[dkarrasch]

The 2-arg *div!s should work, at least they exist.

[dkarrasch]

S**t, I thought I wrote a long comment, but must have forgotten to push the button. To make the long story short, the left-division is currently working nicely, there is no need for an additional method. The right-division is not too bad in that it doesn't fall back to generic dense methods or whatever, but can be improved still

using SparseArrays, LinearAlgebra, BenchmarkTools, Test

A = sprand(10000,10000,0.1);
D = Diagonal(rand(10000));

function mydiv(A::AbstractMatrix, D::Diagonal)
    rdiv!(typeof(oneunit(eltype(A))/oneunit(eltype(D))).(A), D)
end

@test A/D ≈ mydiv(A, D)
@btime mydiv($A, $D);
@btime $A / $D;
@btime $D \ $A; # for comparison with mydiv

My function mydiv is meant to be the relaxation of

julia/stdlib/LinearAlgebra/src/diagonal.jl

Lines 486 to 487 in ac82513

	(/)(A::Union{StridedMatrix, AbstractTriangular}, D::Diagonal) =
	rdiv!((typeof(oneunit(eltype(D))/oneunit(eltype(A)))).(A), D)

to AbstractMatrixes. I'd suggest to relax that method, remove the newly added methods, and keep the tests, maybe adding further sparsity tests.

[eschnett]

@dkarrasch Do you suggest to provide a generic method for AbstractMatrix that handles division by diagonal matrices? It seems dangerous to apply a point-wise operation to them since we know nothing about their internal storage. We can't guarantee that the result of a division by a diagonal matrix can be stored in their internal representation. For example, this function would fail if called for UniformScaling.

Should we instead relax the existing division operator to include sparse matrices only?

[dkarrasch]

Do you suggest to provide a generic method for AbstractMatrix that handles division by diagonal matrices?

Yes. Look at this one:

julia/stdlib/LinearAlgebra/src/diagonal.jl

Lines 597 to 598 in ac82513

	(\)(D::Diagonal, A::AbstractMatrix) =
	ldiv!(D, (typeof(oneunit(eltype(D))/oneunit(eltype(A)))).(A))

If there is an issue with the underlying ldiv! for that specific array type, we'll get a MethodError, and I think that's just fine. IIRC, UniformScaling is not an AbstractMatrix. And this one is called currently for the left-division of a sparse matrix by a Diagonal.

[eschnett]

@dkarrasch Done.

[dkarrasch]

Ok, my sincere apologies for the mess. Now we've got a method ambiguity. Apparently, there is a subtle reason for the different signatures of ldiv! and rdiv!, for why one can handle AbstractMatrix and the other is more specialized. So I hope we should be fine with adding

function (/)(A::AbstractSparseMatrixCSC, D::Diagonal)
    T = typeof(oneunit(eltype(A))/oneunit(eltype(D)))
    rdiv!(LinearAlgebra.copy_oftype(A, T), D)
end

to SparseArrays (and leaving the one in diagonal.jl unchanged compared to master).

Initially, my main point was that D \ A works just fine, so we don't need to change anything to allow for that. A / D also works, but has a little overhead due to the multiple adjoints, so we can improve that.

[dkarrasch]

Sorry again for the back and forth, this time just a trivial comment on the tests, and I hope with my last suggestion we'll be all set.

[eschnett]

@dkarrasch Done. Thanks for the suggestions; this is now a much better implementation.

[dkarrasch]

@eschnett Ok, the tests look good now, but you still have the method ambiguity due to my own bold suggestion (sorry again). You'll need to undo the change in LinearAlgebra/diagonal.jl, and add #35533 (comment) to SparseArrays/linalg.jl. Then I'd expect this to pass tests.

[eschnett]

Sorry, did not see your final comment earlier.

I didn't spot the method ambiguity when I ran the test suite locally.