Forward Diagonal rmul!/lmul! for adj/trans to parent (#1229)

Currently. `rmul!(A'::Adjoint, D::Diagonal)` loops over
`CartesianIndices(A')`, which is not cache-friendly. In this PR, we use
the fact that `A' = A' * D` may be rewritten as `A = D' * A`, and that
`D' isa Diagonal` by design. The operation therefore becomes
`lmul!(D'::Diagonal, A)`, which is cache-friendly.
On master,
```julia
julia> using LinearAlgebra, Chairmarks

julia> D = Diagonal(rand(4000));

julia> A = rand(size(D)...);

julia> @b (A', D) rmul!(_[1], _[2])
103.276 ms
```
whereas, with this PR,
```julia
julia> @b (A', D) rmul!(_[1], _[2])
9.994 ms
```

Author: jishnub

src/diagonal.jl

@@ -360,6 +360,13 @@ function rmul!(A::AbstractMatrix, D::Diagonal)
     end
     return A
 end
+# A' = A' * D => A = D' * A
+# This uses the fact that D' is a Diagonal
+function rmul!(A::AdjOrTransAbsMat, D::Diagonal)
+    f = wrapperop(A)
+    lmul!(f(D), f(A))
+    A
+end
 # T .= T * D
 function rmul!(T::Tridiagonal, D::Diagonal)
     matmul_size_check(size(T), size(D))
@@ -395,6 +402,13 @@ function lmul!(D::Diagonal, T::Tridiagonal)
     end
     return T
 end
+# A' = D * A' => A = A * D'
+# This uses the fact that D' is a Diagonal
+function lmul!(D::Diagonal, A::AdjOrTransAbsMat)
+    f = wrapperop(A)
+    rmul!(f(A), f(D))
+    A
+end
 
 @inline function __muldiag_nonzeroalpha!(out, D::Diagonal, B, alpha::Number, beta::Number)
     @inbounds for j in axes(B, 2)

test/diagonal.jl

@@ -1257,6 +1257,22 @@ end
     end
 end
 
+@testset "rmul!/lmul! for adj/trans" begin
+    for T in (Float64, ComplexF64)
+        A = rand(T,5,4); B = similar(A)
+        for f in (adjoint, transpose)
+            D = Diagonal(rand(T, size(A,1)))
+            B .= A
+            rmul!(f(B), D)
+            @test f(B) == f(A) * D
+            D = Diagonal(rand(T, size(A,2)))
+            B .= A
+            lmul!(D, f(B))
+            @test f(B) == D * f(A)
+        end
+    end
+end
+
 struct SMatrix1{T} <: AbstractArray{T,2}
     elt::T
 end