From 52b76f0949f5a52a1cbc27c0c928a8cc7cc0470d Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 14 Feb 2025 18:27:26 +0530 Subject: [PATCH 01/15] Specialize Diagonal * Adjoint --- src/diagonal.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/diagonal.jl b/src/diagonal.jl index c41cce14..20e5ce10 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -332,6 +332,13 @@ function (*)(D::Diagonal, V::AbstractVector) return D.diag .* V end +function (*)(A::AdjOrTransAbsMat, D::Diagonal) + copy((D' * A')') +end +function (*)(D::Diagonal, A::AdjOrTransAbsMat) + copy((A' * D')') +end + function rmul!(A::AbstractMatrix, D::Diagonal) matmul_size_check(size(A), size(D)) for I in CartesianIndices(A) From c0829bb23f59a2ed7ee98f0660da9302e976f953 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 14 Feb 2025 20:17:25 +0530 Subject: [PATCH 02/15] Get wrapper type from array --- src/diagonal.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 20e5ce10..ad3378b7 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -333,10 +333,12 @@ function (*)(D::Diagonal, V::AbstractVector) end function (*)(A::AdjOrTransAbsMat, D::Diagonal) - copy((D' * A')') + adj = wrapperop(A) + copy(adj(adj(D) * adj(A))) end function (*)(D::Diagonal, A::AdjOrTransAbsMat) - copy((A' * D')') + adj = wrapperop(A) + copy(adj(adj(A) * adj(D))) end function rmul!(A::AbstractMatrix, D::Diagonal) From ea6b1b8af7322b177a22442cc81be7925e893c13 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 14 Feb 2025 18:27:26 +0530 Subject: [PATCH 03/15] Specialize Diagonal * Adjoint --- src/diagonal.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/diagonal.jl b/src/diagonal.jl index c41cce14..20e5ce10 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -332,6 +332,13 @@ function (*)(D::Diagonal, V::AbstractVector) return D.diag .* V end +function (*)(A::AdjOrTransAbsMat, D::Diagonal) + copy((D' * A')') +end +function (*)(D::Diagonal, A::AdjOrTransAbsMat) + copy((A' * D')') +end + function rmul!(A::AbstractMatrix, D::Diagonal) matmul_size_check(size(A), size(D)) for I in CartesianIndices(A) From 7b639c639f89d7fa23fcce366a03ef15b340b340 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 14 Feb 2025 20:17:25 +0530 Subject: [PATCH 04/15] Get wrapper type from array --- src/diagonal.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 20e5ce10..ad3378b7 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -333,10 +333,12 @@ function (*)(D::Diagonal, V::AbstractVector) end function (*)(A::AdjOrTransAbsMat, D::Diagonal) - copy((D' * A')') + adj = wrapperop(A) + copy(adj(adj(D) * adj(A))) end function (*)(D::Diagonal, A::AdjOrTransAbsMat) - copy((A' * D')') + adj = wrapperop(A) + copy(adj(adj(A) * adj(D))) end function rmul!(A::AbstractMatrix, D::Diagonal) From 8879f8c6b01b5657c5ffbf0224a8d1f252f01ce8 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 19 Feb 2025 23:07:46 +0530 Subject: [PATCH 05/15] Add a special path for strided matrices --- src/adjtrans.jl | 4 ++-- src/diagonal.jl | 19 +++++++++++++++++-- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/adjtrans.jl b/src/adjtrans.jl index f7caa9bc..d81aa3ae 100644 --- a/src/adjtrans.jl +++ b/src/adjtrans.jl @@ -319,8 +319,8 @@ const AdjointAbsVec{T} = Adjoint{T,<:AbstractVector} const AdjointAbsMat{T} = Adjoint{T,<:AbstractMatrix} const TransposeAbsVec{T} = Transpose{T,<:AbstractVector} const TransposeAbsMat{T} = Transpose{T,<:AbstractMatrix} -const AdjOrTransAbsVec{T} = AdjOrTrans{T,<:AbstractVector} -const AdjOrTransAbsMat{T} = AdjOrTrans{T,<:AbstractMatrix} +const AdjOrTransAbsVec{T,V<:AbstractVector} = AdjOrTrans{T,V} +const AdjOrTransAbsMat{T,M<:AbstractMatrix} = AdjOrTrans{T,M} # for internal use below wrapperop(_) = identity diff --git a/src/diagonal.jl b/src/diagonal.jl index ad3378b7..6f523a89 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -332,14 +332,29 @@ function (*)(D::Diagonal, V::AbstractVector) return D.diag .* V end -function (*)(A::AdjOrTransAbsMat, D::Diagonal) +function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal) adj = wrapperop(A) copy(adj(adj(D) * adj(A))) end -function (*)(D::Diagonal, A::AdjOrTransAbsMat) +function _diag_adj_mul(A::AdjOrTransAbsMat{<:Any, <:StridedMatrix}, D::Diagonal) + Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) + rmul!(Ac, D) +end +function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) adj = wrapperop(A) copy(adj(adj(A) * adj(D))) end +function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat{<:Any, <:StridedMatrix}) + Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) + lmul!(D, Ac) +end + +function (*)(A::AdjOrTransAbsMat, D::Diagonal) + _diag_adj_mul(A, D) +end +function (*)(D::Diagonal, A::AdjOrTransAbsMat) + _diag_adj_mul(D, A) +end function rmul!(A::AbstractMatrix, D::Diagonal) matmul_size_check(size(A), size(D)) From 7dfba0c824d3b5ec747095d47f0289a1efc580c1 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 20 Feb 2025 09:17:52 +0530 Subject: [PATCH 06/15] Restrict strided methods to `Number` --- src/diagonal.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 6f523a89..a6d41a3c 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -336,7 +336,7 @@ function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal) adj = wrapperop(A) copy(adj(adj(D) * adj(A))) end -function _diag_adj_mul(A::AdjOrTransAbsMat{<:Any, <:StridedMatrix}, D::Diagonal) +function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal) Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) rmul!(Ac, D) end @@ -344,7 +344,7 @@ function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) adj = wrapperop(A) copy(adj(adj(A) * adj(D))) end -function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat{<:Any, <:StridedMatrix}) +function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) lmul!(D, Ac) end From 37f454746afc421d911ed7ebcef196cdc4e84695 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 20 Feb 2025 09:19:21 +0530 Subject: [PATCH 07/15] Restrict `Diagonal` to Number as well in strided methods --- src/diagonal.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index a6d41a3c..0d142996 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -336,7 +336,7 @@ function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal) adj = wrapperop(A) copy(adj(adj(D) * adj(A))) end -function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal) +function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) rmul!(Ac, D) end @@ -344,7 +344,7 @@ function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) adj = wrapperop(A) copy(adj(adj(A) * adj(D))) end -function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) +function _diag_adj_mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) lmul!(D, Ac) end From 640bc84437d6352020ea76af408d08d600c3ae39 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 20 Feb 2025 09:59:29 +0530 Subject: [PATCH 08/15] Use `matprod_dest` to obtain destination --- src/diagonal.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 0d142996..be107420 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -337,7 +337,8 @@ function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal) copy(adj(adj(D) * adj(A))) end function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) - Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) + TS = promote_op(matprod, eltype(A), eltype(D)) + Ac = copyto!(matprod_dest(A, D, TS), A) rmul!(Ac, D) end function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) @@ -345,7 +346,8 @@ function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) copy(adj(adj(A) * adj(D))) end function _diag_adj_mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) - Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) + T = promote_op(matprod, eltype(A), eltype(D)) + Ac = copyto!(matprod_dest(D, A, TS), A) lmul!(D, Ac) end From a51c91d8294c469a365043b398b896b177278de7 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 20 Feb 2025 11:22:22 +0530 Subject: [PATCH 09/15] Call `mul!` instead of `lmul!`/`rmul!` --- src/diagonal.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index be107420..df79de42 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -338,8 +338,8 @@ function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal) end function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) TS = promote_op(matprod, eltype(A), eltype(D)) - Ac = copyto!(matprod_dest(A, D, TS), A) - rmul!(Ac, D) + C = matprod_dest(A, D, TS) + mul!(C, A, D) end function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) adj = wrapperop(A) @@ -347,8 +347,8 @@ function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) end function _diag_adj_mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) T = promote_op(matprod, eltype(A), eltype(D)) - Ac = copyto!(matprod_dest(D, A, TS), A) - lmul!(D, Ac) + C = matprod_dest(A, D, TS) + mul!(C, D, A) end function (*)(A::AdjOrTransAbsMat, D::Diagonal) From 35523844c670a01f7cdc441c5e8c442be6d2d000 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 20 Feb 2025 12:05:05 +0530 Subject: [PATCH 10/15] Restore earlier behavior --- src/diagonal.jl | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index df79de42..5bdd36cc 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -449,10 +449,19 @@ function __muldiag_nonzeroalpha!(out, D::Diagonal, B::UpperOrLowerTriangular, al end @inline function __muldiag_nonzeroalpha_right!(out, A, D::Diagonal, alpha::Number, beta::Number) - @inbounds for j in axes(A, 2) - dja = @stable_muladdmul MulAddMul(alpha,false)(D.diag[j]) - @simd for i in axes(A, 1) - @stable_muladdmul _modify!(MulAddMul(true,beta), A[i,j] * dja, out, (i,j)) + if iszero(beta) + @inbounds for j in axes(A, 2) + dja = D.diag[j] * alpha + @simd for i in axes(A, 1) + out[i,j] = A[i,j] * dja + end + end + else + @inbounds for j in axes(A, 2) + dja = D.diag[j] * alpha + @simd for i in axes(A, 1) + out[i,j] = A[i,j] * dja + out[i,j] * beta + end end end return out From 23d0e50f3cc4b4b129e4f43f3540055e60579d50 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 19 Feb 2025 23:07:46 +0530 Subject: [PATCH 11/15] Add a special path for strided matrices --- src/adjtrans.jl | 4 ++-- src/diagonal.jl | 19 +++++++++++++++++-- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/adjtrans.jl b/src/adjtrans.jl index f7caa9bc..d81aa3ae 100644 --- a/src/adjtrans.jl +++ b/src/adjtrans.jl @@ -319,8 +319,8 @@ const AdjointAbsVec{T} = Adjoint{T,<:AbstractVector} const AdjointAbsMat{T} = Adjoint{T,<:AbstractMatrix} const TransposeAbsVec{T} = Transpose{T,<:AbstractVector} const TransposeAbsMat{T} = Transpose{T,<:AbstractMatrix} -const AdjOrTransAbsVec{T} = AdjOrTrans{T,<:AbstractVector} -const AdjOrTransAbsMat{T} = AdjOrTrans{T,<:AbstractMatrix} +const AdjOrTransAbsVec{T,V<:AbstractVector} = AdjOrTrans{T,V} +const AdjOrTransAbsMat{T,M<:AbstractMatrix} = AdjOrTrans{T,M} # for internal use below wrapperop(_) = identity diff --git a/src/diagonal.jl b/src/diagonal.jl index ad3378b7..0d142996 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -332,14 +332,29 @@ function (*)(D::Diagonal, V::AbstractVector) return D.diag .* V end -function (*)(A::AdjOrTransAbsMat, D::Diagonal) +function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal) adj = wrapperop(A) copy(adj(adj(D) * adj(A))) end -function (*)(D::Diagonal, A::AdjOrTransAbsMat) +function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) + Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) + rmul!(Ac, D) +end +function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) adj = wrapperop(A) copy(adj(adj(A) * adj(D))) end +function _diag_adj_mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) + Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) + lmul!(D, Ac) +end + +function (*)(A::AdjOrTransAbsMat, D::Diagonal) + _diag_adj_mul(A, D) +end +function (*)(D::Diagonal, A::AdjOrTransAbsMat) + _diag_adj_mul(D, A) +end function rmul!(A::AbstractMatrix, D::Diagonal) matmul_size_check(size(A), size(D)) From 35a3df22f4db7ab591d8f8138248b02938f53c00 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 20 Feb 2025 09:59:29 +0530 Subject: [PATCH 12/15] Use `matprod_dest` to obtain destination --- src/diagonal.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 0d142996..bc2508ab 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -337,16 +337,18 @@ function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal) copy(adj(adj(D) * adj(A))) end function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) - Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) - rmul!(Ac, D) + TS = promote_op(matprod, eltype(A), eltype(D)) + C = matprod_dest(A, D, TS) + mul!(C, A, D) end function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) adj = wrapperop(A) copy(adj(adj(A) * adj(D))) end function _diag_adj_mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) - Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D))) - lmul!(D, Ac) + TS = promote_op(matprod, eltype(A), eltype(D)) + C = matprod_dest(D, A, TS) + mul!(C, D, A) end function (*)(A::AdjOrTransAbsMat, D::Diagonal) From c54949c9622b82b4ab0abe6266226f6d5059a247 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 20 Feb 2025 12:49:13 +0100 Subject: [PATCH 13/15] dispatch by `mul` --- src/diagonal.jl | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 793d1f60..9f247162 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -330,26 +330,19 @@ function (*)(D::Diagonal, V::AbstractVector) return D.diag .* V end -function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal) +function mul(A::AdjOrTransAbsMat, D::Diagonal) adj = wrapperop(A) copy(adj(adj(D) * adj(A))) end -function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) - @invoke *(A::AbstractMatrix, D::AbstractMatrix) +function mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) + @invoke mul(A::AbstractMatrix, D::AbstractMatrix) end -function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) +function mul(D::Diagonal, A::AdjOrTransAbsMat) adj = wrapperop(A) copy(adj(adj(A) * adj(D))) end -function _diag_adj_mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) - @invoke *(D::AbstractMatrix, A::AbstractMatrix) -end - -function (*)(A::AdjOrTransAbsMat, D::Diagonal) - _diag_adj_mul(A, D) -end -function (*)(D::Diagonal, A::AdjOrTransAbsMat) - _diag_adj_mul(D, A) +function mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) + @invoke mul(D::AbstractMatrix, A::AbstractMatrix) end function rmul!(A::AbstractMatrix, D::Diagonal) From 33d8a9683e5e4709e90ae090f4cbf085970761c0 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 20 Feb 2025 12:50:09 +0100 Subject: [PATCH 14/15] Revert "dispatch by `mul`" This reverts commit c54949c9622b82b4ab0abe6266226f6d5059a247. --- src/diagonal.jl | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 9f247162..793d1f60 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -330,19 +330,26 @@ function (*)(D::Diagonal, V::AbstractVector) return D.diag .* V end -function mul(A::AdjOrTransAbsMat, D::Diagonal) +function _diag_adj_mul(A::AdjOrTransAbsMat, D::Diagonal) adj = wrapperop(A) copy(adj(adj(D) * adj(A))) end -function mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) - @invoke mul(A::AbstractMatrix, D::AbstractMatrix) +function _diag_adj_mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) + @invoke *(A::AbstractMatrix, D::AbstractMatrix) end -function mul(D::Diagonal, A::AdjOrTransAbsMat) +function _diag_adj_mul(D::Diagonal, A::AdjOrTransAbsMat) adj = wrapperop(A) copy(adj(adj(A) * adj(D))) end -function mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) - @invoke mul(D::AbstractMatrix, A::AbstractMatrix) +function _diag_adj_mul(D::Diagonal{<:Number}, A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}) + @invoke *(D::AbstractMatrix, A::AbstractMatrix) +end + +function (*)(A::AdjOrTransAbsMat, D::Diagonal) + _diag_adj_mul(A, D) +end +function (*)(D::Diagonal, A::AdjOrTransAbsMat) + _diag_adj_mul(D, A) end function rmul!(A::AbstractMatrix, D::Diagonal) From cf947065fd049862f0071c90540809925e839d76 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 20 Feb 2025 19:45:42 +0530 Subject: [PATCH 15/15] Restore `__muldiag_nonzeroalpha_right!` --- src/diagonal.jl | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 793d1f60..5618c313 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -443,19 +443,10 @@ function __muldiag_nonzeroalpha!(out, D::Diagonal, B::UpperOrLowerTriangular, al end @inline function __muldiag_nonzeroalpha_right!(out, A, D::Diagonal, alpha::Number, beta::Number) - if iszero(beta) - @inbounds for j in axes(A, 2) - dja = D.diag[j] * alpha - @simd for i in axes(A, 1) - out[i,j] = A[i,j] * dja - end - end - else - @inbounds for j in axes(A, 2) - dja = D.diag[j] * alpha - @simd for i in axes(A, 1) - out[i,j] = A[i,j] * dja + out[i,j] * beta - end + @inbounds for j in axes(A, 2) + dja = @stable_muladdmul MulAddMul(alpha,false)(D.diag[j]) + @simd for i in axes(A, 1) + @stable_muladdmul _modify!(MulAddMul(true,beta), A[i,j] * dja, out, (i,j)) end end return out