diff --git a/NEWS.md b/NEWS.md index 1540203be012a..61609b4539632 100644 --- a/NEWS.md +++ b/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]) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index e531dbcfdb303..8f5a22098825b 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/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)) diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index b3467621fcc10..637b52bd28453 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/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)