diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl
index b5071b178de100..fe45102d90d778 100644
--- a/stdlib/LinearAlgebra/src/hessenberg.jl
+++ b/stdlib/LinearAlgebra/src/hessenberg.jl
@@ -132,11 +132,11 @@ for T = (:Number, :UniformScaling, :Diagonal)
 end
 
 function *(H::UpperHessenberg, U::UpperOrUnitUpperTriangular)
-    HH = _mulmattri!(_initarray(*, eltype(H), eltype(U), H), H, U)
+    HH = mul!(_initarray(*, eltype(H), eltype(U), H), H, U)
     UpperHessenberg(HH)
 end
 function *(U::UpperOrUnitUpperTriangular, H::UpperHessenberg)
-    HH = _multrimat!(_initarray(*, eltype(U), eltype(H), H), U, H)
+    HH = mul!(_initarray(*, eltype(U), eltype(H), H), U, H)
     UpperHessenberg(HH)
 end
 
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index 295a46f1522a55..873dcb8ec0d156 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -464,9 +464,9 @@ end
 
 # Define `mul!` for (Unit){Upper,Lower}Triangular matrices times a number.
 # be permissive here and require compatibility later in _triscale!
-@inline mul!(A::UpperOrLowerTriangular, B::UpperOrLowerTriangular, C::Number, alpha::Number, beta::Number) =
+@inline mul!(A::AbstractTriangular, B::AbstractTriangular, C::Number, alpha::Number, beta::Number) =
     _triscale!(A, B, C, MulAddMul(alpha, beta))
-@inline mul!(A::UpperOrLowerTriangular, B::Number, C::UpperOrLowerTriangular, alpha::Number, beta::Number) =
+@inline mul!(A::AbstractTriangular, B::Number, C::AbstractTriangular, alpha::Number, beta::Number) =
     _triscale!(A, B, C, MulAddMul(alpha, beta))
 
 function _triscale!(A::UpperTriangular, B::UpperTriangular, c::Number, _add)
@@ -671,11 +671,50 @@ fillstored!(A::UnitUpperTriangular, x) = (fillband!(A.data, x, 1, size(A,2)-1);
 # BlasFloat routines #
 ######################
 
+# which triangle to use of the underlying data
+uplo_char(::UpperOrUnitUpperTriangular) = 'U'
+uplo_char(::LowerOrUnitLowerTriangular) = 'L'
+uplo_char(::UpperOrUnitUpperTriangular{<:Any,<:AdjOrTrans}) = 'L'
+uplo_char(::LowerOrUnitLowerTriangular{<:Any,<:AdjOrTrans}) = 'U'
+uplo_char(::UpperOrUnitUpperTriangular{<:Any,<:Adjoint{<:Any,<:Transpose}}) = 'U'
+uplo_char(::LowerOrUnitLowerTriangular{<:Any,<:Adjoint{<:Any,<:Transpose}}) = 'L'
+uplo_char(::UpperOrUnitUpperTriangular{<:Any,<:Transpose{<:Any,<:Adjoint}}) = 'U'
+uplo_char(::LowerOrUnitLowerTriangular{<:Any,<:Transpose{<:Any,<:Adjoint}}) = 'L'
+
+isunit_char(::UpperTriangular) = 'N'
+isunit_char(::UnitUpperTriangular) = 'U'
+isunit_char(::LowerTriangular) = 'N'
+isunit_char(::UnitLowerTriangular) = 'U'
+
 lmul!(A::Tridiagonal, B::AbstractTriangular) = A*full!(B)
-mul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVector) = _multrimat!(C, A, B)
-mul!(C::AbstractMatrix, A::AbstractTriangular, B::AbstractMatrix) = _multrimat!(C, A, B)
-mul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractTriangular) = _mulmattri!(C, A, B)
-mul!(C::AbstractMatrix, A::AbstractTriangular, B::AbstractTriangular) = _multrimat!(C, A, B)
+mul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVector) = _trimul!(C, A, B)
+mul!(C::AbstractMatrix, A::AbstractTriangular, B::AbstractMatrix) = _trimul!(C, A, B)
+mul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractTriangular) = _trimul!(C, A, B)
+mul!(C::AbstractMatrix, A::AbstractTriangular, B::AbstractTriangular) = _trimul!(C, A, B)
+
+# generic fallback for AbstractTriangular matrices outside of the four subtypes provided here
+_trimul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVector) =
+    lmul!(A, inplace_adj_or_trans(B)(C, _parent(B)))
+_trimul!(C::AbstractMatrix, A::AbstractTriangular, B::AbstractMatrix) =
+    lmul!(A, inplace_adj_or_trans(B)(C, _parent(B)))
+_trimul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractTriangular) =
+    rmul!(inplace_adj_or_trans(A)(C, _parent(A)), B)
+_trimul!(C::AbstractMatrix, A::AbstractTriangular, B::AbstractTriangular) =
+    lmul!(A, copyto!(C, B))
+# redirect for UpperOrLowerTriangular
+_trimul!(C::AbstractVecOrMat, A::UpperOrLowerTriangular, B::AbstractVector) =
+    generic_trimatmul!(C, uplo_char(A), isunit_char(A), adj_or_trans(parent(A)), _parent(parent(A)), B)
+_trimul!(C::AbstractMatrix, A::UpperOrLowerTriangular, B::AbstractMatrix) =
+    generic_trimatmul!(C, uplo_char(A), isunit_char(A), adj_or_trans(parent(A)), _parent(parent(A)), B)
+_trimul!(C::AbstractMatrix, A::AbstractMatrix, B::UpperOrLowerTriangular) =
+    generic_mattrimul!(C, uplo_char(B), isunit_char(B), adj_or_trans(parent(B)), A, _parent(parent(B)))
+_trimul!(C::AbstractMatrix, A::UpperOrLowerTriangular, B::UpperOrLowerTriangular) =
+    generic_trimatmul!(C, uplo_char(A), isunit_char(A), adj_or_trans(parent(A)), _parent(parent(A)), B)
+# disambiguation with AbstractTriangular
+_trimul!(C::AbstractMatrix, A::UpperOrLowerTriangular, B::AbstractTriangular) =
+    generic_trimatmul!(C, uplo_char(A), isunit_char(A), adj_or_trans(parent(A)), _parent(parent(A)), B)
+_trimul!(C::AbstractMatrix, A::AbstractTriangular, B::UpperOrLowerTriangular) =
+    generic_mattrimul!(C, uplo_char(B), isunit_char(B), adj_or_trans(parent(B)), A, _parent(parent(B)))
 
 for TC in (:AbstractVector, :AbstractMatrix)
     @eval @inline function mul!(C::$TC, A::AbstractTriangular, B::AbstractVector, alpha::Number, beta::Number)
@@ -699,12 +738,6 @@ for (TA, TB) in ((:AbstractTriangular, :AbstractMatrix),
     end
 end
 
-
-# generic fallback for AbstractTriangular matrices outside of the four subtypes provided here
-_multrimat!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) =
-    lmul!(A, inplace_adj_or_trans(B)(C, _parent(B)))
-_mulmattri!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractTriangular) = rmul!(copyto!(C, A), B)
-
 # preserve triangular structure in in-place multiplication
 for (cty, aty, bty) in ((:UpperTriangular, :UpperTriangular, :UpperTriangular),
                         (:UpperTriangular, :UpperTriangular, :UnitUpperTriangular),
@@ -714,8 +747,8 @@ for (cty, aty, bty) in ((:UpperTriangular, :UpperTriangular, :UpperTriangular),
                         (:LowerTriangular, :LowerTriangular, :UnitLowerTriangular),
                         (:LowerTriangular, :UnitLowerTriangular, :LowerTriangular),
                         (:UnitLowerTriangular, :UnitLowerTriangular, :UnitLowerTriangular))
-    @eval function _multrimat!(C::$cty, A::$aty, B::$bty)
-        _multrimat!(parent(C), A, B)
+    @eval function _trimul!(C::$cty, A::$aty, B::$bty)
+        _trimul!(parent(C), A, B)
         return C
     end
 end
@@ -726,16 +759,6 @@ for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'),
                             (:UpperTriangular, 'U', 'N'),
                             (:UnitUpperTriangular, 'U', 'U'))
     @eval begin
-        # Vector multiplication
-        lmul!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} =
-            BLAS.trmv!($uploc, 'N', $isunitc, A.data, b)
-
-        # Matrix multiplication
-        lmul!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} =
-            BLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), A.data, B)
-        rmul!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} =
-            BLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), B.data, A)
-
         # Left division
         ldiv!(A::$t{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
             LAPACK.trtrs!($uploc, 'N', $isunitc, A.data, B)
@@ -772,29 +795,6 @@ for (t, uploc, isunitc) in ((:LowerTriangular, 'U', 'N'),
                             (:UpperTriangular, 'L', 'N'),
                             (:UnitUpperTriangular, 'L', 'U'))
     @eval begin
-        # Vector multiplication
-        lmul!(A::$t{<:Any,<:Transpose{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasFloat} =
-            BLAS.trmv!($uploc, 'T', $isunitc, parent(parent(A)), b)
-        lmul!(A::$t{<:Any,<:Adjoint{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasReal} =
-            BLAS.trmv!($uploc, 'T', $isunitc, parent(parent(A)), b)
-        lmul!(A::$t{<:Any,<:Adjoint{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasComplex} =
-            BLAS.trmv!($uploc, 'C', $isunitc, parent(parent(A)), b)
-
-        # Matrix multiplication
-        lmul!(A::$t{<:Any,<:Transpose{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasFloat} =
-            BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), parent(parent(A)), B)
-        lmul!(A::$t{<:Any,<:Adjoint{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasComplex} =
-            BLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), parent(parent(A)), B)
-        lmul!(A::$t{<:Any,<:Adjoint{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasReal} =
-            BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), parent(parent(A)), B)
-
-        rmul!(A::StridedMatrix{T}, B::$t{<:Any,<:Transpose{T,<:StridedMatrix}}) where {T<:BlasFloat} =
-            BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), parent(parent(B)), A)
-        rmul!(A::StridedMatrix{T}, B::$t{<:Any,<:Adjoint{T,<:StridedMatrix}}) where {T<:BlasComplex} =
-            BLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), parent(parent(B)), A)
-        rmul!(A::StridedMatrix{T}, B::$t{<:Any,<:Adjoint{T,<:StridedMatrix}}) where {T<:BlasReal} =
-            BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), parent(parent(B)), A)
-
         # Left division
         ldiv!(A::$t{<:Any,<:Transpose{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
             LAPACK.trtrs!($uploc, 'T', $isunitc, parent(parent(A)), B)
@@ -813,21 +813,17 @@ for (t, uploc, isunitc) in ((:LowerTriangular, 'U', 'N'),
     end
 end
 
+# Vector multiplication
+generic_trimatmul!(c::StridedVector{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, b::AbstractVector{T}) where {T<:BlasFloat} =
+    BLAS.trmv!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, c === b ? c : copyto!(c, b))
+# Matrix multiplication
+generic_trimatmul!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, B::AbstractMatrix{T}) where {T<:BlasFloat} =
+    BLAS.trmm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), A, C === B ? C : copyto!(C, B))
+generic_mattrimul!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function, A::AbstractMatrix{T}, B::StridedMatrix{T}) where {T<:BlasFloat} =
+    BLAS.trmm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), B, C === A ? C : copyto!(C, A))
+
 # redirect back to BLAS
 for t in (:UpperTriangular, :UnitUpperTriangular, :LowerTriangular, :UnitLowerTriangular)
-    @eval _multrimat!(C::StridedVecOrMat{T}, A::$t{T,<:StridedMatrix}, B::AbstractVecOrMat{T}) where {T<:BlasFloat} =
-        lmul!(A, copyto!(C, B))
-    @eval _multrimat!(C::StridedVecOrMat{T}, A::$t{<:Any,<:Adjoint{T,<:StridedMatrix}}, B::AbstractVecOrMat{T}) where {T<:BlasFloat} =
-        lmul!(A, copyto!(C, B))
-    @eval _multrimat!(C::StridedVecOrMat{T}, A::$t{<:Any,<:Transpose{T,<:StridedMatrix}}, B::AbstractVecOrMat{T}) where {T<:BlasFloat} =
-        lmul!(A, copyto!(C, B))
-    @eval _mulmattri!(C::StridedMatrix{T}, A::AbstractMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} =
-        rmul!(copyto!(C, A), B)
-    @eval _mulmattri!(C::StridedMatrix{T}, A::AbstractMatrix{T}, B::$t{<:Any,<:Adjoint{T,<:StridedMatrix}}) where {T<:BlasFloat} =
-        rmul!(copyto!(C, A), B)
-    @eval _mulmattri!(C::StridedMatrix{T}, A::AbstractMatrix{T}, B::$t{<:Any,<:Transpose{T,<:StridedMatrix}}) where {T<:BlasFloat} =
-        rmul!(copyto!(C, A), B)
-
     @eval ldiv!(C::StridedVecOrMat{T}, A::$t{T,<:StridedMatrix}, B::AbstractVecOrMat{T}) where {T<:BlasFloat} =
         ldiv!(A, copyto!(C, B))
     @eval ldiv!(C::StridedVecOrMat{T}, A::$t{<:Any,<:Adjoint{T,<:StridedMatrix}}, B::AbstractVecOrMat{T}) where {T<:BlasFloat} =
@@ -886,6 +882,9 @@ end
 # Generic routines #
 ####################
 
+lmul!(A::UpperOrLowerTriangular, B::AbstractVecOrMat) = @inline _trimul!(B, A, B)
+rmul!(A::AbstractMatrix, B::UpperOrLowerTriangular)   = @inline _trimul!(A, A, B)
+
 for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
                    (LowerTriangular, UnitLowerTriangular))
     @eval begin
@@ -930,17 +929,11 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
             end
             $t(B)
         end
-
-        lmul!(A::$t, B::AbstractVecOrMat)     = @inline _multrimat!(B, A, B)
-        lmul!(A::$unitt, B::AbstractVecOrMat) = @inline _multrimat!(B, A, B)
-
-        rmul!(A::AbstractMatrix, B::$t)     = @inline _mulmattri!(A, A, B)
-        rmul!(A::AbstractMatrix, B::$unitt) = @inline _mulmattri!(A, A, B)
     end
 end
 
 ## Generic triangular multiplication
-function _multrimat!(C::AbstractVecOrMat, A::UpperTriangular, B::AbstractVecOrMat)
+function generic_trimatmul!(C::AbstractVecOrMat, uploc, isunitc, tfun, A::AbstractMatrix, B::AbstractVecOrMat)
     require_one_based_indexing(C, A, B)
     m, n = size(B, 1), size(B, 2)
     N = size(A, 1)
@@ -951,41 +944,104 @@ function _multrimat!(C::AbstractVecOrMat, A::UpperTriangular, B::AbstractVecOrMa
     if mc != N || nc != n
         throw(DimensionMismatch("output has dimensions ($mc,$nc), should have ($N,$n)"))
     end
-    @inbounds for j in 1:n
-        for i in 1:m
-            Cij = A.data[i,i] * B[i,j]
-            for k in i + 1:m
-                Cij += A.data[i,k] * B[k,j]
+    if uploc == 'U'
+        if isunitc == 'N'
+            if tfun === identity
+                @inbounds for j in 1:n
+                    for i in 1:m
+                        Cij = A[i,i] * B[i,j]
+                        for k in i + 1:m
+                            Cij += A[i,k] * B[k,j]
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            else # tfun in (transpose, adjoint)
+                @inbounds for j in 1:n
+                    for i in m:-1:1
+                        Cij = tfun(A[i,i]) * B[i,j]
+                        for k in 1:i - 1
+                            Cij += tfun(A[k,i]) * B[k,j]
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            end
+        else # unitc == 'U'
+            if tfun === identity
+                @inbounds for j in 1:n
+                    for i in 1:m
+                        Cij = oneunit(eltype(A)) * B[i,j]
+                        for k in i + 1:m
+                            Cij += A[i,k] * B[k,j]
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            else # tfun in (transpose, adjoint)
+                @inbounds for j in 1:n
+                    for i in m:-1:1
+                        Cij = oneunit(eltype(A)) * B[i,j]
+                        for k in 1:i - 1
+                            Cij += tfun(A[k,i]) * B[k,j]
+                        end
+                        C[i,j] = Cij
+                    end
+                end
             end
-            C[i,j] = Cij
         end
-    end
-    return C
-end
-function _multrimat!(C::AbstractVecOrMat, A::UnitUpperTriangular, B::AbstractVecOrMat)
-    require_one_based_indexing(C, A, B)
-    m, n = size(B, 1), size(B, 2)
-    N = size(A, 1)
-    if m != N
-        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
-    end
-
-    mc, nc = size(C, 1), size(C, 2)
-    if mc != N || nc != n
-        throw(DimensionMismatch("output has dimensions ($mc,$nc), should have ($N,$n)"))
-    end
-    @inbounds for j in 1:n
-        for i in 1:m
-            Cij = oneunit(eltype(A)) * B[i,j]
-            for k in i + 1:m
-                Cij += A.data[i,k] * B[k,j]
+    else # uploc == 'L'
+        if isunitc == 'N'
+            if tfun === identity
+                @inbounds for j in 1:n
+                    for i in m:-1:1
+                        Cij = A[i,i] * B[i,j]
+                        for k in 1:i - 1
+                            Cij += A[i,k] * B[k,j]
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            else # tfun in (transpose, adjoint)
+                @inbounds for j in 1:n
+                    for i in 1:m
+                        Cij = tfun(A[i,i]) * B[i,j]
+                        for k in i + 1:m
+                            Cij += tfun(A[k,i]) * B[k,j]
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            end
+        else # isunitc == 'U'
+            if tfun === identity
+                @inbounds for j in 1:n
+                    for i in m:-1:1
+                        Cij = oneunit(eltype(A)) * B[i,j]
+                        for k in 1:i - 1
+                            Cij += A[i,k] * B[k,j]
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            else # tfun in (transpose, adjoint)
+                @inbounds for j in 1:n
+                    for i in 1:m
+                        Cij = oneunit(eltype(A)) * B[i,j]
+                        for k in i + 1:m
+                            Cij += tfun(A[k,i]) * B[k,j]
+                        end
+                        C[i,j] = Cij
+                    end
+                end
             end
-            C[i,j] = Cij
         end
     end
     return C
 end
-function _multrimat!(C::AbstractVecOrMat, A::LowerTriangular, B::AbstractVecOrMat)
+# Conjugate cases transpose(adjoint(A)) and adjoint(transpose(A))
+function generic_trimatmul!(C::AbstractVecOrMat, uploc, isunitc, _, xA::AdjOrTrans, B::AbstractVecOrMat)
+    A = parent(xA)
     require_one_based_indexing(C, A, B)
     m, n = size(B, 1), size(B, 2)
     N = size(A, 1)
@@ -996,41 +1052,55 @@ function _multrimat!(C::AbstractVecOrMat, A::LowerTriangular, B::AbstractVecOrMa
     if mc != N || nc != n
         throw(DimensionMismatch("output has dimensions ($mc,$nc), should have ($N,$n)"))
     end
-    @inbounds for j in 1:n
-        for i in m:-1:1
-            Cij = A.data[i,i] * B[i,j]
-            for k in 1:i - 1
-                Cij += A.data[i,k] * B[k,j]
+    if uploc == 'U'
+        if isunitc == 'N'
+            @inbounds for j in 1:n
+                for i in 1:m
+                    Cij = conj(A[i,i]) * B[i,j]
+                    for k in i + 1:m
+                        Cij += conj(A[i,k]) * B[k,j]
+                    end
+                    C[i,j] = Cij
+                end
+            end
+        else # unitc == 'U'
+            @inbounds for j in 1:n
+                for i in 1:m
+                    Cij = oneunit(eltype(A)) * B[i,j]
+                    for k in i + 1:m
+                        Cij += conj(A[i,k]) * B[k,j]
+                    end
+                    C[i,j] = Cij
+                end
             end
-            C[i,j] = Cij
         end
-    end
-    return C
-end
-function _multrimat!(C::AbstractVecOrMat, A::UnitLowerTriangular, B::AbstractVecOrMat)
-    require_one_based_indexing(C, A, B)
-    m, n = size(B, 1), size(B, 2)
-    N = size(A, 1)
-    if m != N
-        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
-    end
-    mc, nc = size(C, 1), size(C, 2)
-    if mc != N || nc != n
-        throw(DimensionMismatch("output has dimensions ($mc,$nc), should have ($N,$n)"))
-    end
-    @inbounds for j in 1:n
-        for i in m:-1:1
-            Cij = oneunit(eltype(A)) * B[i,j]
-            for k in 1:i - 1
-                Cij += A.data[i,k] * B[k,j]
+    else # uploc == 'L'
+        if isunitc == 'N'
+            @inbounds for j in 1:n
+                for i in m:-1:1
+                    Cij = conj(A[i,i]) * B[i,j]
+                    for k in 1:i - 1
+                        Cij += conj(A[i,k]) * B[k,j]
+                    end
+                    C[i,j] = Cij
+                end
+            end
+        else # isunitc == 'U'
+            @inbounds for j in 1:n
+                for i in m:-1:1
+                    Cij = oneunit(eltype(A)) * B[i,j]
+                    for k in 1:i - 1
+                        Cij += conj(A[i,k]) * B[k,j]
+                    end
+                    C[i,j] = Cij
+                end
             end
-            C[i,j] = Cij
         end
     end
     return C
 end
 
-function _mulmattri!(C::AbstractMatrix, A::AbstractMatrix, B::UpperTriangular)
+function generic_mattrimul!(C::AbstractMatrix, uploc, isunitc, tfun, A::AbstractMatrix, B::AbstractMatrix)
     require_one_based_indexing(C, A, B)
     m, n = size(A, 1), size(A, 2)
     N = size(B, 1)
@@ -1041,40 +1111,104 @@ function _mulmattri!(C::AbstractMatrix, A::AbstractMatrix, B::UpperTriangular)
     if mc != m || nc != N
         throw(DimensionMismatch("output has dimensions ($mc,$nc), should have ($m,$N)"))
     end
-    @inbounds for i in 1:m
-        for j in n:-1:1
-            Cij = A[i,j] * B.data[j,j]
-            for k in 1:j - 1
-                Cij += A[i,k] * B.data[k,j]
+    if uploc == 'U'
+        if isunitc == 'N'
+            if tfun === identity
+                @inbounds for i in 1:m
+                    for j in n:-1:1
+                        Cij = A[i,j] * tfun(B[j,j])
+                        for k in 1:j - 1
+                            Cij += A[i,k] * tfun(B[k,j])
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            else # tfun in (transpose, adjoint)
+                @inbounds for i in 1:m
+                    for j in 1:n
+                        Cij = A[i,j] * tfun(B[j,j])
+                        for k in j + 1:n
+                            Cij += A[i,k] * tfun(B[j,k])
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            end
+        else # isunitc == 'U'
+            if tfun === identity
+                @inbounds for i in 1:m
+                    for j in n:-1:1
+                        Cij = A[i,j] * oneunit(eltype(B))
+                        for k in 1:j - 1
+                            Cij += A[i,k] * tfun(B[k,j])
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            else # tfun in (transpose, adjoint)
+                @inbounds for i in 1:m
+                    for j in 1:n
+                        Cij = A[i,j] * oneunit(eltype(B))
+                        for k in j + 1:n
+                            Cij += A[i,k] * tfun(B[j,k])
+                        end
+                        C[i,j] = Cij
+                    end
+                end
             end
-            C[i,j] = Cij
         end
-    end
-    return C
-end
-function _mulmattri!(C::AbstractMatrix, A::AbstractMatrix, B::UnitUpperTriangular)
-    require_one_based_indexing(C, A, B)
-    m, n = size(A, 1), size(A, 2)
-    N = size(B, 1)
-    if n != N
-        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $N"))
-    end
-    mc, nc = size(C, 1), size(C, 2)
-    if mc != m || nc != N
-        throw(DimensionMismatch("output has dimensions ($mc,$nc), should have ($m,$N)"))
-    end
-    @inbounds for i in 1:m
-        for j in n:-1:1
-            Cij = A[i,j] * oneunit(eltype(B))
-            for k in 1:j - 1
-                Cij += A[i,k] * B.data[k,j]
+    else # uploc == 'L'
+        if isunitc == 'N'
+            if tfun === identity
+                @inbounds for i in 1:m
+                    for j in 1:n
+                        Cij = A[i,j] * tfun(B[j,j])
+                        for k in j + 1:n
+                            Cij += A[i,k] * tfun(B[k,j])
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            else # tfun in (transpose, adjoint)
+                @inbounds for i in 1:m
+                    for j in n:-1:1
+                        Cij = A[i,j] * tfun(B[j,j])
+                        for k in 1:j - 1
+                            Cij += A[i,k] * tfun(B[j,k])
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            end
+        else # unitc == 'U'
+            if tfun === identity
+                @inbounds for i in 1:m
+                    for j in 1:n
+                        Cij = A[i,j] * oneunit(eltype(B))
+                        for k in j + 1:n
+                            Cij += A[i,k] * tfun(B[k,j])
+                        end
+                        C[i,j] = Cij
+                    end
+                end
+            else # tfun in (transpose, adjoint)
+                @inbounds for i in 1:m
+                    for j in n:-1:1
+                        Cij = A[i,j] * oneunit(eltype(B))
+                        for k in 1:j - 1
+                            Cij += A[i,k] * tfun(B[j,k])
+                        end
+                        C[i,j] = Cij
+                    end
+                end
             end
-            C[i,j] = Cij
         end
     end
     return C
 end
-function _mulmattri!(C::AbstractMatrix, A::AbstractMatrix, B::LowerTriangular)
+# Conjugate cases
+function generic_mattrimul!(C::AbstractMatrix, uploc, isunitc, _, A::AbstractMatrix, xB::AdjOrTrans)
+    B = parent(xB)
     require_one_based_indexing(C, A, B)
     m, n = size(A, 1), size(A, 2)
     N = size(B, 1)
@@ -1085,35 +1219,49 @@ function _mulmattri!(C::AbstractMatrix, A::AbstractMatrix, B::LowerTriangular)
     if mc != m || nc != N
         throw(DimensionMismatch("output has dimensions ($mc,$nc), should have ($m,$N)"))
     end
-    @inbounds for i in 1:m
-        for j in 1:n
-            Cij = A[i,j] * B.data[j,j]
-            for k in j + 1:n
-                Cij += A[i,k] * B.data[k,j]
+    if uploc == 'U'
+        if isunitc == 'N'
+            @inbounds for i in 1:m
+                for j in n:-1:1
+                    Cij = A[i,j] * conj(B[j,j])
+                    for k in 1:j - 1
+                        Cij += A[i,k] * conj(B[k,j])
+                    end
+                    C[i,j] = Cij
+                end
+            end
+        else # isunitc == 'U'
+            @inbounds for i in 1:m
+                for j in n:-1:1
+                    Cij = A[i,j] * oneunit(eltype(B))
+                    for k in 1:j - 1
+                        Cij += A[i,k] * conj(B[k,j])
+                    end
+                    C[i,j] = Cij
+                end
             end
-            C[i,j] = Cij
         end
-    end
-    return C
-end
-function _mulmattri!(C::AbstractMatrix, A::AbstractMatrix, B::UnitLowerTriangular)
-    require_one_based_indexing(C, A, B)
-    m, n = size(A, 1), size(A, 2)
-    N = size(B, 1)
-    if n != N
-        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $N"))
-    end
-    mc, nc = size(C, 1), size(C, 2)
-    if mc != m || nc != N
-        throw(DimensionMismatch("output has dimensions ($mc,$nc), should have ($m,$N)"))
-    end
-    @inbounds for i in 1:m
-        for j in 1:n
-            Cij = A[i,j] * oneunit(eltype(B))
-            for k in j + 1:n
-                Cij += A[i,k] * B.data[k,j]
+    else # uploc == 'L'
+        if isunitc == 'N'
+            @inbounds for i in 1:m
+                for j in 1:n
+                    Cij = A[i,j] * conj(B[j,j])
+                    for k in j + 1:n
+                        Cij += A[i,k] * conj(B[k,j])
+                    end
+                    C[i,j] = Cij
+                end
+            end
+        else # unitc == 'U'
+            @inbounds for i in 1:m
+                for j in 1:n
+                    Cij = A[i,j] * oneunit(eltype(B))
+                    for k in j + 1:n
+                        Cij += A[i,k] * conj(B[k,j])
+                    end
+                    C[i,j] = Cij
+                end
             end
-            C[i,j] = Cij
         end
     end
     return C
@@ -1394,11 +1542,7 @@ _inner_type_promotion(op, ::Type{TA}, ::Type{TB}) where {TA,TB} =
 ## The general promotion methods
 function *(A::AbstractTriangular, B::AbstractTriangular)
     TAB = _init_eltype(*, eltype(A), eltype(B))
-    if TAB <: BlasFloat
-        lmul!(convert(AbstractArray{TAB}, A), copy_similar(B, TAB))
-    else
-        mul!(similar(B, TAB, size(B)), A, B)
-    end
+    mul!(similar(B, TAB, size(B)), A, B)
 end
 
 for mat in (:AbstractVector, :AbstractMatrix)
@@ -1406,11 +1550,7 @@ for mat in (:AbstractVector, :AbstractMatrix)
     @eval function *(A::AbstractTriangular, B::$mat)
         require_one_based_indexing(B)
         TAB = _init_eltype(*, eltype(A), eltype(B))
-        if TAB <: BlasFloat
-            lmul!(convert(AbstractArray{TAB}, A), copy_similar(B, TAB))
-        else
-            mul!(similar(B, TAB, size(B)), A, B)
-        end
+        mul!(similar(B, TAB, size(B)), A, B)
     end
     ### Left division with triangle to the left hence rhs cannot be transposed. No quotients.
     @eval function \(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::$mat)
@@ -1458,11 +1598,7 @@ end
 function *(A::AbstractMatrix, B::AbstractTriangular)
     require_one_based_indexing(A)
     TAB = _init_eltype(*, eltype(A), eltype(B))
-    if TAB <: BlasFloat
-        rmul!(copy_similar(A, TAB), convert(AbstractArray{TAB}, B))
-    else
-        mul!(similar(A, TAB, size(A)), A, B)
-    end
+    mul!(similar(A, TAB, size(A)), A, B)
 end
 # ambiguity resolution with definitions in matmul.jl
 *(v::AdjointAbsVec, A::AbstractTriangular) = adjoint(adjoint(A) * v.parent)