Added kron for Diagonal matrices

NEWS.md

@@ -716,6 +716,9 @@ Library improvements
   * Added an optimized method of `vecdot` for taking the Frobenius inner product
     of sparse matrices. ([#27470])
 
+  * Added an optimized method of `kron` for taking the tensor product of two
+    `Diagonal` matrices. ([27581])
+
   * The initial element `v0` in `reduce(op, v0, itr)` has been replaced with an `init`
     optional keyword argument, as in `reduce(op, itr; init=v0)`. Similarly for `foldl`,
     `foldr`, `mapreduce`, `mapfoldl` and `mapfoldr`. ([#27711])

stdlib/LinearAlgebra/src/diagonal.jl

@@ -356,6 +356,16 @@ rdiv!(A::AbstractMatrix{T}, transD::Transpose{<:Any,<:Diagonal{T}}) where {T} =
 (\)(A::Union{QR,QRCompactWY,QRPivoted}, B::Diagonal) =
     invoke(\, Tuple{Union{QR,QRCompactWY,QRPivoted}, AbstractVecOrMat}, A, B)
 
+function kron(A::Diagonal{T1}, B::Diagonal{T2}) where {T1<:Number, T2<:Number}
+    valA = A.diag; nA = length(valA)
+    valB = B.diag; nB = length(valB)
+    valC = Vector{typeof(zero(T1)*zero(T2))}(undef,nA*nB)
+    @inbounds for i = 1:nA, j = 1:nB
+        valC[(i-1)*nB+j] = valA[i] * valB[j]
+    end
+    return Diagonal(valC)
+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

@@ -187,6 +187,11 @@ srand(1)
         DD = copy(D)
         r  = VV * Array(D)'
         @test Array(rmul!(VV, adjoint(DD))) ≈ r
+
+        # kron
+        D3 = Diagonal(convert(Vector{elty}, rand(n÷2)))
+        DM3= Matrix(D3)
+        Matrix(kron(D, D3)) ≈ kron(DM, DM3)
     end
     @testset "triu/tril" begin
         @test istriu(D)