Implement methods for kron with a Diagonal argument (#30232)

* This implements methods for kron(::Diagonal, ::AbstractMatrix) and
  kron(::AbstractMatrix, ::Diagonal). These methods are typically
  twice as fast as the fallback methods when the other argument
  is a Matrix.

* The test for the existing method for kron(::Diagonal, ::Diagonal)
  was missing the @test macro, probably inadvertently. I added the
  macro and checked that the test indeed passes.

stdlib/LinearAlgebra/src/diagonal.jl

@@ -395,6 +395,44 @@ function kron(A::Diagonal{T1}, B::Diagonal{T2}) where {T1<:Number, T2<:Number}
     return Diagonal(valC)
 end
 
+function kron(A::Diagonal{T}, B::AbstractMatrix{S}) where {T<:Number, S<:Number}
+    @assert ! Base.has_offset_axes(B)
+    (mA, nA) = size(A); (mB, nB) = size(B)
+    R = zeros(Base.promote_op(*, T, S), mA * mB, nA * nB)
+    m = 1
+    for j = 1:nA
+        A_jj = A[j,j]
+        for k = 1:nB
+            for l = 1:mB
+                R[m] = A_jj * B[l,k]
+                m += 1
+            end
+            m += (nA - 1) * mB
+        end
+        m += mB
+    end
+    return R
+end
+
+function kron(A::AbstractMatrix{T}, B::Diagonal{S}) where {T<:Number, S<:Number}
+    @assert ! has_offset_axes(A)
+    (mA, nA) = size(A); (mB, nB) = size(B)
+    R = zeros(promote_op(*, T, S), mA * mB, nA * nB)
+    m = 1
+    for j = 1:nA
+        for l = 1:mB
+            Bll = B[l,l]
+            for k = 1:mA
+                R[m] = A[k,j] * Bll
+                m += nB
+            end
+            m += 1
+        end
+        m -= nB
+    end
+    return R
+end
+
 conj(D::Diagonal) = Diagonal(conj(D.diag))
 transpose(D::Diagonal{<:Number}) = D
 transpose(D::Diagonal) = Diagonal(transpose.(D.diag))

stdlib/LinearAlgebra/test/diagonal.jl

@@ -199,7 +199,10 @@ Random.seed!(1)
         # kron
         D3 = Diagonal(convert(Vector{elty}, rand(n÷2)))
         DM3= Matrix(D3)
-        Matrix(kron(D, D3)) ≈ kron(DM, DM3)
+        @test Matrix(kron(D, D3)) ≈ kron(DM, DM3)
+        M4 = rand(elty, n÷2, n÷2)
+        @test kron(D3, M4) ≈ kron(DM3, M4)
+        @test kron(M4, D3) ≈ kron(M4, DM3)
     end
     @testset "iszero, isone, triu, tril" begin
         Dzero = Diagonal(zeros(elty, 10))