Make use of mul indirection (#1215)

Author: dkarrasch

src/bidiag.jl

@@ -1183,25 +1183,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

src/diagonal.jl

@@ -322,7 +322,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
@@ -332,26 +332,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)
@@ -1037,9 +1030,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.
-*(A::AbstractTriangular, D::Diagonal) = @invoke *(A::AbstractMatrix, D::Diagonal)
-*(D::Diagonal, A::AbstractTriangular) = @invoke *(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)

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,19 +146,23 @@ 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
 
+/(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)

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

src/symmetric.jl

@@ -702,7 +702,20 @@ 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 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 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
+            convert(AbstractMatrix{T}, B))
+end
 
 function dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
     require_one_based_indexing(x, y)

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))

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 = Diagonal(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