diff --git a/src/bidiag.jl b/src/bidiag.jl
index bb5b8830..89e5711c 100644
--- a/src/bidiag.jl
+++ b/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
diff --git a/src/diagonal.jl b/src/diagonal.jl
index 8ded64ad..c4f88276 100644
--- a/src/diagonal.jl
+++ b/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)
diff --git a/src/hessenberg.jl b/src/hessenberg.jl
index ed654c33..056b97e3 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,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)
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 b059f317..11b34a99 100644
--- a/src/symmetric.jl
+++ b/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)
diff --git a/src/symmetriceigen.jl b/src/symmetriceigen.jl
index 68a1b29f..1d92e4f8 100644
--- a/src/symmetriceigen.jl
+++ b/src/symmetriceigen.jl
@@ -6,10 +6,16 @@ eigencopy_oftype(A::Hermitian, S) = Hermitian(copytrito!(similar(parent(A), S, s
 eigencopy_oftype(A::Symmetric, S) = Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo))
 eigencopy_oftype(A::Symmetric{<:Complex}, S) = copyto!(similar(parent(A), S), A)
 
-default_eigen_alg(A) = DivideAndConquer()
+"""
+    default_eigen_alg(A)
+
+Return the default algorithm used to solve the eigensystem `A v = λ v` for a symmetric matrix `A`.
+Defaults to `LinearAlegbra.DivideAndConquer()`, which corresponds to the LAPACK function `LAPACK.syevd!`.
+"""
+default_eigen_alg(@nospecialize(A)) = DivideAndConquer()
 
 # Eigensolvers for symmetric and Hermitian matrices
-function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
+function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
     if alg === DivideAndConquer()
         Eigen(sorteig!(LAPACK.syevd!('V', A.uplo, A.data)..., sortby)...)
     elseif alg === QRIteration()
@@ -22,7 +28,7 @@ function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algo
 end
 
 """
-    eigen(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A)) -> Eigen
+    eigen(A::Union{Hermitian, Symmetric}; alg::LinearAlgebra.Algorithm = LinearAlgebra.default_eigen_alg(A)) -> Eigen
 
 Compute the eigenvalue decomposition of `A`, returning an [`Eigen`](@ref) factorization object `F`
 which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
@@ -45,19 +51,19 @@ The default `alg` used may change in the future.
 
 The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).
 """
-function eigen(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
-    _eigen(A, alg; sortby)
+function eigen(A::RealHermSymComplexHerm; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
+    _eigen(A; alg, sortby)
 end
 
 # we dispatch on the eltype in an internal method to avoid ambiguities
-function _eigen(A::RealHermSymComplexHerm, alg::Algorithm; sortby)
+function _eigen(A::RealHermSymComplexHerm; alg::Algorithm, sortby)
     S = eigtype(eltype(A))
-    eigen!(eigencopy_oftype(A, S), alg; sortby)
+    eigen!(eigencopy_oftype(A, S); alg, sortby)
 end
 
-function _eigen(A::RealHermSymComplexHerm{Float16}, alg::Algorithm; sortby::Union{Function,Nothing}=nothing)
+function _eigen(A::RealHermSymComplexHerm{Float16}; alg::Algorithm, sortby::Union{Function,Nothing}=nothing)
     S = eigtype(eltype(A))
-    E = eigen!(eigencopy_oftype(A, S), alg, sortby=sortby)
+    E = eigen!(eigencopy_oftype(A, S); alg, sortby)
     values = convert(AbstractVector{Float16}, E.values)
     vectors = convert(AbstractMatrix{isreal(E.vectors) ? Float16 : Complex{Float16}}, E.vectors)
     return Eigen(values, vectors)
@@ -114,7 +120,7 @@ function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real)
 end
 
 
-function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
+function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
     vals::Vector{real(eltype(A))} = if alg === DivideAndConquer()
         LAPACK.syevd!('N', A.uplo, A.data)
     elseif alg === QRIteration()
@@ -129,7 +135,7 @@ function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, alg::Al
 end
 
 """
-    eigvals(A::Union{Hermitian, Symmetric}, alg::Algorithm = default_eigen_alg(A))) -> values
+    eigvals(A::Union{Hermitian, Symmetric}; alg::Algorithm = default_eigen_alg(A))) -> values
 
 Return the eigenvalues of `A`.
 
@@ -143,9 +149,9 @@ a comparison of the accuracy and performance of different methods.
 
 The default `alg` used may change in the future.
 """
-function eigvals(A::RealHermSymComplexHerm, alg::Algorithm = default_eigen_alg(A); sortby::Union{Function,Nothing}=nothing)
+function eigvals(A::RealHermSymComplexHerm; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
     S = eigtype(eltype(A))
-    eigvals!(eigencopy_oftype(A, S), alg; sortby)
+    eigvals!(eigencopy_oftype(A, S); alg, sortby)
 end
 
 
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))
diff --git a/test/matmul.jl b/test/matmul.jl
index afa730f1..86c75ae5 100644
--- a/test/matmul.jl
+++ b/test/matmul.jl
@@ -325,6 +325,10 @@ end
             @test 0 == @allocations mul!(C, At, Bt)
         end
         # syrk/herk
+        mul!(C, transpose(A), A)
+        mul!(C, adjoint(A), A)
+        mul!(C, A, transpose(A))
+        mul!(C, A, adjoint(A))
         @test 0 == @allocations mul!(C, transpose(A), A)
         @test 0 == @allocations mul!(C, adjoint(A), A)
         @test 0 == @allocations mul!(C, A, transpose(A))
@@ -334,6 +338,7 @@ end
         Ac = complex(A)
         for t in (identity, adjoint, transpose)
             Bt = t(B)
+            mul!(Cc, Ac, Bt)
             @test 0 == @allocations mul!(Cc, Ac, Bt)
         end
     end
@@ -356,6 +361,9 @@ end
         A = rand(-10:10, n, n)
         B = ones(Float64, n, n)
         C = zeros(Float64, n, n)
+        mul!(C, A, B)
+        mul!(C, A, transpose(B))
+        mul!(C, adjoint(A), B)
         @test 0 == @allocations mul!(C, A, B)
         @test 0 == @allocations mul!(C, A, transpose(B))
         @test 0 == @allocations mul!(C, adjoint(A), B)
diff --git a/test/symmetric.jl b/test/symmetric.jl
index edd3af48..68f0ff47 100644
--- a/test/symmetric.jl
+++ b/test/symmetric.jl
@@ -720,21 +720,6 @@ end
     end
 end
 
-@testset "eigendecomposition Algorithms" begin
-    using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations
-    for T in (Float64, ComplexF64, Float32, ComplexF32)
-        n = 4
-        A = T <: Real ? Symmetric(randn(T, n, n)) : Hermitian(randn(T, n, n))
-        d, v = eigen(A)
-        for alg in (DivideAndConquer(), QRIteration(), RobustRepresentations())
-            @test (@inferred eigvals(A, alg)) ≈ d
-            d2, v2 = @inferred eigen(A, alg)
-            @test d2 ≈ d
-            @test A * v2 ≈ v2 * Diagonal(d2)
-        end
-    end
-end
-
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
 isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl"))
 using .Main.ImmutableArrays
diff --git a/test/symmetriceigen.jl b/test/symmetriceigen.jl
index 71087ae4..7a2b88a5 100644
--- a/test/symmetriceigen.jl
+++ b/test/symmetriceigen.jl
@@ -3,6 +3,7 @@
 module TestSymmetricEigen
 
 using Test, LinearAlgebra
+using LinearAlgebra: DivideAndConquer, QRIteration, RobustRepresentations
 
 @testset "chol-eigen-eigvals" begin
     ## Cholesky decomposition based
@@ -173,7 +174,7 @@ end
     @test D.vectors ≈ D32.vectors
 
     # ensure that different algorithms dispatch correctly
-    λ, V = eigen(C, LinearAlgebra.QRIteration())
+    λ, V = eigen(C; alg=LinearAlgebra.QRIteration())
     @test λ isa Vector{Float16}
     @test C * V ≈ V * Diagonal(λ)
 end
@@ -184,4 +185,18 @@ end
     @test S * v ≈ v * Diagonal(λ)
 end
 
+@testset "eigendecomposition Algorithms" begin
+    for T in (Float64, ComplexF64, Float32, ComplexF32)
+        n = 4
+        A = T <: Real ? Symmetric(randn(T, n, n)) : Hermitian(randn(T, n, n))
+        d, v = eigen(A)
+        for alg in (DivideAndConquer(), QRIteration(), RobustRepresentations())
+            @test (@inferred eigvals(A; alg)) ≈ d
+            d2, v2 = @inferred eigen(A; alg)
+            @test d2 ≈ d
+            @test A * v2 ≈ v2 * Diagonal(d2)
+        end
+    end
+end
+
 end # module TestSymmetricEigen
diff --git a/test/triangular.jl b/test/triangular.jl
index 8eeeb666..eb6814b9 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 = 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