From cf9163fd57eece7d2435f2bddca139af0cb69565 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 20 Feb 2025 12:41:22 +0100 Subject: [PATCH 1/6] Make use of `mul` indirection --- src/bidiag.jl | 8 ++++---- src/diagonal.jl | 6 +++--- src/hessenberg.jl | 8 +++++--- src/special.jl | 4 ++-- src/symmetric.jl | 6 +++--- src/triangular.jl | 2 +- 6 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/bidiag.jl b/src/bidiag.jl index 90385ab7..8698fc54 100644 --- a/src/bidiag.jl +++ b/src/bidiag.jl @@ -1184,25 +1184,25 @@ function _dibimul!(C::Bidiagonal, A::Diagonal, B::Bidiagonal, _add) C end -function *(A::UpperOrUnitUpperTriangular, B::Bidiagonal) +function mul(A::UpperOrUnitUpperTriangular, B::Bidiagonal) TS = promote_op(matprod, eltype(A), eltype(B)) C = mul!(similar(A, TS, size(A)), A, B) return B.uplo == 'U' ? UpperTriangular(C) : C end -function *(A::LowerOrUnitLowerTriangular, B::Bidiagonal) +function mul(A::LowerOrUnitLowerTriangular, B::Bidiagonal) TS = promote_op(matprod, eltype(A), eltype(B)) C = mul!(similar(A, TS, size(A)), A, B) return B.uplo == 'L' ? LowerTriangular(C) : C end -function *(A::Bidiagonal, B::UpperOrUnitUpperTriangular) +function mul(A::Bidiagonal, B::UpperOrUnitUpperTriangular) TS = promote_op(matprod, eltype(A), eltype(B)) C = mul!(similar(B, TS, size(B)), A, B) return A.uplo == 'U' ? UpperTriangular(C) : C end -function *(A::Bidiagonal, B::LowerOrUnitLowerTriangular) +function mul(A::Bidiagonal, B::LowerOrUnitLowerTriangular) TS = promote_op(matprod, eltype(A), eltype(B)) C = mul!(similar(B, TS, size(B)), A, B) return A.uplo == 'L' ? LowerTriangular(C) : C diff --git a/src/diagonal.jl b/src/diagonal.jl index 27ff3b1f..cb373184 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -320,7 +320,7 @@ Base.literal_pow(::typeof(^), D::Diagonal, valp::Val) = Diagonal(Base.literal_pow.(^, D.diag, valp)) # for speed Base.literal_pow(::typeof(^), D::Diagonal, ::Val{-1}) = inv(D) # for disambiguation -function (*)(Da::Diagonal, Db::Diagonal) +function mul(Da::Diagonal, Db::Diagonal) matmul_size_check(size(Da), size(Db)) return Diagonal(Da.diag .* Db.diag) end @@ -1018,8 +1018,8 @@ end /(u::AdjointAbsVec, D::Diagonal) = (D' \ u')' /(u::TransposeAbsVec, D::Diagonal) = transpose(transpose(D) \ transpose(u)) # disambiguation methods: Call unoptimized version for user defined AbstractTriangular. -*(A::AbstractTriangular, D::Diagonal) = @invoke *(A::AbstractMatrix, D::Diagonal) -*(D::Diagonal, A::AbstractTriangular) = @invoke *(D::Diagonal, A::AbstractMatrix) +mul(A::AbstractTriangular, D::Diagonal) = @invoke mul(A::AbstractMatrix, D::Diagonal) +mul(D::Diagonal, A::AbstractTriangular) = @invoke mul(D::Diagonal, A::AbstractMatrix) _opnorm1(A::Diagonal) = maximum(norm(x) for x in A.diag) _opnormInf(A::Diagonal) = maximum(norm(x) for x in A.diag) diff --git a/src/hessenberg.jl b/src/hessenberg.jl index ed654c33..0dc9f17d 100644 --- a/src/hessenberg.jl +++ b/src/hessenberg.jl @@ -137,7 +137,7 @@ for T = (:UniformScaling, :Diagonal, :Bidiagonal, :Tridiagonal, :SymTridiagonal, end end -for T = (:Number, :UniformScaling, :Diagonal) +for T = (:Number, :UniformScaling) @eval begin *(H::UpperHessenberg, x::$T) = UpperHessenberg(H.data * x) *(x::$T, H::UpperHessenberg) = UpperHessenberg(x * H.data) @@ -146,11 +146,13 @@ for T = (:Number, :UniformScaling, :Diagonal) end end -function *(H::UpperHessenberg, U::UpperOrUnitUpperTriangular) +mul(H::UpperHessenberg, D::Diagonal) = UpperHessenberg(H.data * D) +mul(D::Diagonal, H::UpperHessenberg) = UpperHessenberg(D * H.data) +function mul(H::UpperHessenberg, U::UpperOrUnitUpperTriangular) HH = mul!(matprod_dest(H, U, promote_op(matprod, eltype(H), eltype(U))), H, U) UpperHessenberg(HH) end -function *(U::UpperOrUnitUpperTriangular, H::UpperHessenberg) +function mul(U::UpperOrUnitUpperTriangular, H::UpperHessenberg) HH = mul!(matprod_dest(U, H, promote_op(matprod, eltype(U), eltype(H))), U, H) UpperHessenberg(HH) end diff --git a/src/special.jl b/src/special.jl index 636fd8d1..58863858 100644 --- a/src/special.jl +++ b/src/special.jl @@ -120,12 +120,12 @@ _mul!(C::AbstractMatrix, A::AbstractTriangular, B::BandedMatrix, alpha::Number, _mul!(C::AbstractMatrix, A::BandedMatrix, B::AbstractTriangular, alpha::Number, beta::Number) = @stable_muladdmul _mul!(C, A, B, MulAddMul(alpha, beta)) -function *(H::UpperHessenberg, B::Bidiagonal) +function mul(H::UpperHessenberg, B::Bidiagonal) T = promote_op(matprod, eltype(H), eltype(B)) A = mul!(similar(H, T, size(H)), H, B) return B.uplo == 'U' ? UpperHessenberg(A) : A end -function *(B::Bidiagonal, H::UpperHessenberg) +function mul(B::Bidiagonal, H::UpperHessenberg) T = promote_op(matprod, eltype(B), eltype(H)) A = mul!(similar(H, T, size(H)), B, H) return B.uplo == 'U' ? UpperHessenberg(A) : A diff --git a/src/symmetric.jl b/src/symmetric.jl index 8e92109f..6c06dec1 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -702,15 +702,15 @@ for f in (:+, :-) end end -*(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B) +mul(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B) # catch a few potential BLAS-cases -function *(A::HermOrSym{<:BlasFloat,<:StridedMatrix}, B::AdjOrTrans{<:BlasFloat,<:StridedMatrix}) +function mul(A::HermOrSym{<:BlasFloat,<:StridedMatrix}, B::AdjOrTrans{<:BlasFloat,<:StridedMatrix}) T = promote_type(eltype(A), eltype(B)) mul!(similar(B, T, (size(A, 1), size(B, 2))), convert(AbstractMatrix{T}, A), copy_oftype(B, T)) # make sure the AdjOrTrans wrapper is resolved end -function *(A::AdjOrTrans{<:BlasFloat,<:StridedMatrix}, B::HermOrSym{<:BlasFloat,<:StridedMatrix}) +function mul(A::AdjOrTrans{<:BlasFloat,<:StridedMatrix}, B::HermOrSym{<:BlasFloat,<:StridedMatrix}) T = promote_type(eltype(A), eltype(B)) mul!(similar(B, T, (size(A, 1), size(B, 2))), copy_oftype(A, T), # make sure the AdjOrTrans wrapper is resolved diff --git a/src/triangular.jl b/src/triangular.jl index 3faa12da..4a03cff6 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -1886,7 +1886,7 @@ end ## Some Triangular-Triangular cases. We might want to write tailored methods ## for these cases, but I'm not sure it is worth it. -for f in (:*, :\) +for f in (:mul, :\) @eval begin ($f)(A::LowerTriangular, B::LowerTriangular) = LowerTriangular(@invoke $f(A::LowerTriangular, B::AbstractMatrix)) From 6e332dec9fc62439a1772a22bcf6614838005226 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 20 Feb 2025 14:35:16 +0100 Subject: [PATCH 2/6] fix --- src/hessenberg.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/hessenberg.jl b/src/hessenberg.jl index 0dc9f17d..056b97e3 100644 --- a/src/hessenberg.jl +++ b/src/hessenberg.jl @@ -157,10 +157,12 @@ function mul(U::UpperOrUnitUpperTriangular, H::UpperHessenberg) UpperHessenberg(HH) end +/(H::UpperHessenberg, D::Diagonal) = UpperHessenberg(H.data / D) function /(H::UpperHessenberg, U::UpperTriangular) HH = _rdiv!(matprod_dest(H, U, promote_op(/, eltype(H), eltype(U))), H, U) UpperHessenberg(HH) end +\(D::Diagonal, H::UpperHessenberg) = UpperHessenberg(D \ H.data) function /(H::UpperHessenberg, U::UnitUpperTriangular) HH = _rdiv!(matprod_dest(H, U, promote_op(/, eltype(H), eltype(U))), H, U) UpperHessenberg(HH) From ea046289e7acd75529c599cf9e435ec32842392c Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 21 Feb 2025 10:36:39 +0100 Subject: [PATCH 3/6] switch to `mul` in diagonal.jl --- src/diagonal.jl | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 69cd1fd8..f1515926 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -330,28 +330,21 @@ 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}) +function mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) @invoke *(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}) +function 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) matmul_size_check(size(A), size(D)) for I in CartesianIndices(A) From 915b8eb2fbba48936194ecf2350a9e8e03a6a6bb Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 27 Feb 2025 11:15:03 +0100 Subject: [PATCH 4/6] Apply suggestions from code review --- src/diagonal.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 61b69323..cef55429 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -335,14 +335,14 @@ function mul(A::AdjOrTransAbsMat, D::Diagonal) copy(adj(adj(D) * adj(A))) end function mul(A::AdjOrTransAbsMat{<:Number, <:StridedMatrix}, D::Diagonal{<:Number}) - @invoke *(A::AbstractMatrix, D::AbstractMatrix) + @invoke mul(A::AbstractMatrix, D::AbstractMatrix) end function 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 *(D::AbstractMatrix, A::AbstractMatrix) + @invoke mul(D::AbstractMatrix, A::AbstractMatrix) end function rmul!(A::AbstractMatrix, D::Diagonal) From e8bcbad8410ed74efc57497963d28c6c129fc475 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 28 Feb 2025 10:31:37 +0100 Subject: [PATCH 5/6] remove methods, but test abstracttriangular-diag-mul --- src/diagonal.jl | 3 --- test/triangular.jl | 5 +++++ 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index cef55429..6a0652a3 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -1046,9 +1046,6 @@ end *(x::TransposeAbsVec, D::Diagonal, y::AbstractVector) = _mapreduce_prod(*, x, D, y) /(u::AdjointAbsVec, D::Diagonal) = (D' \ u')' /(u::TransposeAbsVec, D::Diagonal) = transpose(transpose(D) \ transpose(u)) -# disambiguation methods: Call unoptimized version for user defined AbstractTriangular. -mul(A::AbstractTriangular, D::Diagonal) = @invoke mul(A::AbstractMatrix, D::Diagonal) -mul(D::Diagonal, A::AbstractTriangular) = @invoke mul(D::Diagonal, A::AbstractMatrix) _opnorm1(A::Diagonal) = maximum(norm(x) for x in A.diag) _opnormInf(A::Diagonal) = maximum(norm(x) for x in A.diag) diff --git a/test/triangular.jl b/test/triangular.jl index 8eeeb666..a0ef5947 100644 --- a/test/triangular.jl +++ b/test/triangular.jl @@ -703,6 +703,7 @@ end @testset "(l/r)mul! and (l/r)div! for generic triangular" begin @testset for T in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular) M = MyTriangular(T(rand(4,4))) + D = randn(4) A = rand(4,4) Ac = similar(A) @testset "lmul!" begin @@ -725,6 +726,10 @@ end rdiv!(Ac, M) @test Ac ≈ A / M end + @testset "diagonal mul" begin + @test D * M ≈ D * M.data + @test M * D ≈ M.data * D + end end end From f7f4deb4392c58cb37611cef9632fa69d7a260d3 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 28 Feb 2025 13:29:00 +0100 Subject: [PATCH 6/6] Update test/triangular.jl --- test/triangular.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/triangular.jl b/test/triangular.jl index a0ef5947..eb6814b9 100644 --- a/test/triangular.jl +++ b/test/triangular.jl @@ -703,7 +703,7 @@ end @testset "(l/r)mul! and (l/r)div! for generic triangular" begin @testset for T in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular) M = MyTriangular(T(rand(4,4))) - D = randn(4) + D = Diagonal(randn(4)) A = rand(4,4) Ac = similar(A) @testset "lmul!" begin