Generic path in ldiv! for Diagonal (#1261)

The performance improvement in `ldiv!` is best seen using
`SparseArrays`, as this seems to lack a specialized `ldiv!` method
currently:
```julia
julia> D = Diagonal(rand(3000));

julia> S = sprand(size(D,1), 0.01);

julia> @Btime ldiv!($D, $S);
  29.620 μs (0 allocations: 0 bytes) # master
  15.388 μs (22 allocations: 154.47 KiB) # This PR
```
This is because the broadcasting method is taken. Unfortunately, in this
example, the in-place broadcasting for a sparse array is allocating, but
this should be allocation-free in general.

Author: jishnub

src/diagonal.jl

@@ -612,10 +612,10 @@ function _rdiv!(B::AbstractVecOrMat, A::AbstractVecOrMat, D::Diagonal)
     if (k = length(dd)) != n
         throw(DimensionMismatch(lazy"left hand side has $n columns but D is $k by $k"))
     end
-    @inbounds for j in 1:n
+    @inbounds for j in axes(A,2)
         ddj = dd[j]
         iszero(ddj) && throw(SingularException(j))
-        for i in 1:m
+        for i in axes(A,1)
             B[i, j] = A[i, j] / ddj
         end
     end
@@ -640,7 +640,17 @@ function ldiv!(B::AbstractVecOrMat, D::Diagonal, A::AbstractVecOrMat)
     (m, n) == (m′, n′) || throw(DimensionMismatch(lazy"expect output to be $m by $n, but got $m′ by $n′"))
     j = findfirst(iszero, D.diag)
     isnothing(j) || throw(SingularException(j))
-    @inbounds for j = 1:n, i = 1:m
+    _ldiv_Diagonal_loop!(B, D, A)
+    B
+end
+function _ldiv_Diagonal_loop!(B::AbstractVecOrMat, D::Diagonal, A::AbstractVecOrMat)
+    dd = D.diag
+    @. B = dd \ A
+    B
+end
+function _ldiv_Diagonal_loop!(B::StridedVecOrMat, D::Diagonal, A::StridedVecOrMat)
+    dd = D.diag
+    @inbounds for j in axes(A,2), i in axes(A,1)
         B[i, j] = dd[i] \ A[i, j]
     end
     B