diff --git a/base/hashing2.jl b/base/hashing2.jl
index 83ef9ea09683f..55d67ca9caef9 100644
--- a/base/hashing2.jl
+++ b/base/hashing2.jl
@@ -178,4 +178,4 @@ function hash(s::Union{String,SubString{String}}, h::UInt)
     # note: use pointer(s) here (see #6058).
     ccall(memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), pointer(s), sizeof(s), h % UInt32) + h
 end
-hash(s::AbstractString, h::UInt) = hash(String(s), h)
\ No newline at end of file
+hash(s::AbstractString, h::UInt) = hash(String(s), h)
diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index c82d6f67d8b9a..4d18d7d5d00fd 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -240,11 +240,9 @@ function char_uplo(uplo::Symbol)
 end
 
 """
-    ldiv!([Y,] A, B) -> Y
+    ldiv!(Y, A, B) -> Y
 
 Compute `A \\ B` in-place and store the result in `Y`, returning the result.
-If only two arguments are passed, then `ldiv!(A, B)` overwrites `B` with
-the result.
 
 The argument `A` should *not* be a matrix.  Rather, instead of matrices it should be a
 factorization object (e.g. produced by [`factorize`](@ref) or [`cholfact`](@ref)).
@@ -256,11 +254,24 @@ control over the factorization of `A`.
 ldiv!(Y, A, B)
 
 """
-    rdiv!([Y,] A, B) -> Y
+    ldiv!(A, B)
 
-Compute `A / B` in-place and store the result in `Y`, returning the result.
-If only two arguments are passed, then `rdiv!(A, B)` overwrites `A` with
-the result.
+Compute `A \\ B` in-place and overwriting `B` to store the result.
+
+The argument `A` should *not* be a matrix.  Rather, instead of matrices it should be a
+factorization object (e.g. produced by [`factorize`](@ref) or [`cholfact`](@ref)).
+The reason for this is that factorization itself is both expensive and typically allocates memory
+(although it can also be done in-place via, e.g., [`lufact!`](@ref)),
+and performance-critical situations requiring `ldiv!` usually also require fine-grained
+control over the factorization of `A`.
+"""
+ldiv!(A, B)
+
+
+"""
+    rdiv!(A, B)
+
+Compute `A / B` in-place and overwriting `A` to store the result.
 
 The argument `B` should *not* be a matrix.  Rather, instead of matrices it should be a
 factorization object (e.g. produced by [`factorize`](@ref) or [`cholfact`](@ref)).
@@ -269,7 +280,7 @@ The reason for this is that factorization itself is both expensive and typically
 and performance-critical situations requiring `rdiv!` usually also require fine-grained
 control over the factorization of `B`.
 """
-rdiv!(Y, A, B)
+rdiv!(A, B)
 
 copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A)
 copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A)
diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl
index 6e63a0d353016..3859eb5bdfa86 100644
--- a/stdlib/LinearAlgebra/src/deprecated.jl
+++ b/stdlib/LinearAlgebra/src/deprecated.jl
@@ -727,11 +727,11 @@ export A_ldiv_B!,
 @deprecate A_mul_B!(C::StridedVecOrMat, S::SymTridiagonal, B::StridedVecOrMat)  mul!(C, S, B)
 
 # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/diagonal.jl, to deprecate
-@deprecate A_mul_B!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) mul!(A, D)
-@deprecate A_mul_B!(A::UnitLowerTriangular, D::Diagonal)    mul!(A, D)
-@deprecate A_mul_B!(A::UnitUpperTriangular, D::Diagonal)    mul!(A, D)
-@deprecate A_mul_B!(D::Diagonal, B::UnitLowerTriangular)    mul!(D, B)
-@deprecate A_mul_B!(D::Diagonal, B::UnitUpperTriangular)    mul!(D, B)
+@deprecate A_mul_B!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) mul1!(A, D)
+@deprecate A_mul_B!(A::UnitLowerTriangular, D::Diagonal)    mul1!(A, D)
+@deprecate A_mul_B!(A::UnitUpperTriangular, D::Diagonal)    mul1!(A, D)
+@deprecate A_mul_B!(D::Diagonal, B::UnitLowerTriangular)    mul2!(D, B)
+@deprecate A_mul_B!(D::Diagonal, B::UnitUpperTriangular)    mul2!(D, B)
 @deprecate Ac_mul_B(D::Diagonal, B::Diagonal)           (*)(adjoint(D), B)
 @deprecate Ac_mul_B(A::AbstractTriangular, D::Diagonal) (*)(adjoint(A), D)
 @deprecate Ac_mul_B(A::AbstractMatrix, D::Diagonal)     (*)(adjoint(A), D)
@@ -747,16 +747,19 @@ export A_ldiv_B!,
 @deprecate A_mul_Bt(D::Diagonal, A::AbstractMatrix) (*)(D, transpose(A))
 @deprecate Ac_mul_Bc(D::Diagonal, B::Diagonal)      (*)(adjoint(D), adjoint(B))
 @deprecate At_mul_Bt(D::Diagonal, B::Diagonal)      (*)(transpose(D), transpose(B))
-@deprecate A_mul_B!(A::Diagonal,B::Diagonal)        mul!(A, B)
-@deprecate At_mul_B!(A::Diagonal,B::Diagonal)       mul!(transpose(A), B)
-@deprecate Ac_mul_B!(A::Diagonal,B::Diagonal)       mul!(adjoint(A), B)
-@deprecate A_mul_B!(A::QRPackedQ, D::Diagonal)      mul!(A, D)
-@deprecate A_mul_B!(A::Diagonal,B::AbstractMatrix)      mul!(A, B)
-@deprecate At_mul_B!(A::Diagonal,B::AbstractMatrix)     mul!(transpose(A), B)
-@deprecate Ac_mul_B!(A::Diagonal,B::AbstractMatrix)     mul!(adjoint(A), B)
-@deprecate A_mul_B!(A::AbstractMatrix,B::Diagonal)      mul!(A, B)
-@deprecate A_mul_Bt!(A::AbstractMatrix,B::Diagonal)     mul!(A, transpose(B))
-@deprecate A_mul_Bc!(A::AbstractMatrix,B::Diagonal)     mul!(A, adjoint(B))
+function A_mul_B!(A::Diagonal,B::Diagonal)
+    depwarn("`A_mul_B!(A::Diagonal,B::Diagonal)` should be replaced with `mul1!(A, B)` or `mul2!(A, B)`.", :A_mul_B!)
+    throw(MethodError(A_mul_B!, (A, B)))
+end
+@deprecate At_mul_B!(A::Diagonal,B::Diagonal)       mul2!(transpose(A), B)
+@deprecate Ac_mul_B!(A::Diagonal,B::Diagonal)       mul2!(adjoint(A), B)
+@deprecate A_mul_B!(A::QRPackedQ, D::Diagonal)      mul1!(A, D)
+@deprecate A_mul_B!(A::Diagonal,B::AbstractMatrix)      mul2!(A, B)
+@deprecate At_mul_B!(A::Diagonal,B::AbstractMatrix)     mul2!(transpose(A), B)
+@deprecate Ac_mul_B!(A::Diagonal,B::AbstractMatrix)     mul2!(adjoint(A), B)
+@deprecate A_mul_B!(A::AbstractMatrix,B::Diagonal)      mul1!(A, B)
+@deprecate A_mul_Bt!(A::AbstractMatrix,B::Diagonal)     mul1!(A, transpose(B))
+@deprecate A_mul_Bc!(A::AbstractMatrix,B::Diagonal)     mul1!(A, adjoint(B))
 @deprecate A_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector)       mul!(out, A, in)
 @deprecate Ac_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector)      mul!(out, adjoint(A), in)
 @deprecate At_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector)      mul!(out, transpose(A), in)
@@ -780,7 +783,7 @@ export A_ldiv_B!,
 @deprecate A_ldiv_B!(D::Diagonal, B::StridedVecOrMat)   ldiv!(D, B)
 
 # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/special.jl, to deprecate
-@deprecate A_mul_Bc!(A::AbstractTriangular, B::Union{QRCompactWYQ,QRPackedQ})   mul!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::AbstractTriangular, B::Union{QRCompactWYQ,QRPackedQ})   mul1!(A, adjoint(B))
 @deprecate A_mul_Bc(A::AbstractTriangular, B::Union{QRCompactWYQ,QRPackedQ})    (*)(A, adjoint(B))
 
 # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/bunchkaufman.jl, to deprecate
@@ -805,10 +808,10 @@ export A_ldiv_B!,
 @deprecate At_ldiv_B(F::Factorization, B)   (\)(transpose(F), B)
 
 # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/hessenberg.jl, to deprecate
-@deprecate A_mul_B!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat}  mul!(Q, X)
-@deprecate A_mul_B!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat}    mul!(X, Q)
-@deprecate Ac_mul_B!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} mul!(adjoint(Q), X)
-@deprecate A_mul_Bc!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat}   mul!(X, adjoint(Q))
+@deprecate A_mul_B!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat}  mul2!(Q, X)
+@deprecate A_mul_B!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat}    mul1!(X, Q)
+@deprecate Ac_mul_B!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} mul2!(adjoint(Q), X)
+@deprecate A_mul_Bc!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat}   mul1!(X, adjoint(Q))
 @deprecate Ac_mul_B(Q::HessenbergQ{T}, X::StridedVecOrMat{S}) where {T,S}   (*)(adjoint(Q), X)
 @deprecate A_mul_Bc(X::StridedVecOrMat{S}, Q::HessenbergQ{T}) where {T,S}   (*)(X, adjoint(Q))
 
@@ -858,43 +861,43 @@ export A_ldiv_B!,
 @deprecate Ac_ldiv_B!(A::LU{T,Tridiagonal{T,V}}, B::AbstractVecOrMat) where {T,V}   ldiv!(adjoint(A), B)
 
 # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/lq.jl, to deprecate
-@deprecate A_mul_B!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat}   mul!(A, B)
-@deprecate A_mul_B!(A::LQ{T}, B::QR{T}) where {T<:BlasFloat}    mul!(A, B)
-@deprecate A_mul_B!(A::QR{T}, B::LQ{T}) where {T<:BlasFloat}    mul!(A, B)
-@deprecate A_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat}    mul!(A, B)
-@deprecate Ac_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasReal}    mul!(adjoint(A), B)
-@deprecate Ac_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasComplex} mul!(adjoint(A), B)
+@deprecate A_mul_B!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat}   mul2!(A, B)
+@deprecate A_mul_B!(A::LQ{T}, B::QR{T}) where {T<:BlasFloat}    A*B
+@deprecate A_mul_B!(A::QR{T}, B::LQ{T}) where {T<:BlasFloat}    A*B
+@deprecate A_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat}    mul2!(A, B)
+@deprecate Ac_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasReal}    mul2!(adjoint(A), B)
+@deprecate Ac_mul_B!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasComplex} mul2!(adjoint(A), B)
 @deprecate Ac_mul_B(A::LQPackedQ, B::StridedVecOrMat)   (*)(adjoint(A), B)
 @deprecate A_mul_Bc(A::LQPackedQ, B::StridedVecOrMat)   (*)(A, adjoint(B))
 @deprecate Ac_mul_Bc(A::LQPackedQ, B::StridedVecOrMat)  (*)(adjoint(A), adjoint(B))
-@deprecate A_mul_B!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat}      mul!(A, B)
-@deprecate A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasReal}      mul!(A, adjoint(B))
-@deprecate A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasComplex}   mul!(A, adjoint(B))
+@deprecate A_mul_B!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat}      mul1!(A, B)
+@deprecate A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasReal}      mul1!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasComplex}   mul1!(A, adjoint(B))
 @deprecate A_mul_Bc(A::StridedVecOrMat, Q::LQPackedQ)   (*)(A, adjoint(Q))
 @deprecate Ac_mul_Bc(A::StridedMatrix, Q::LQPackedQ)    (*)(adjoint(A), adjoint(Q))
 @deprecate Ac_mul_B(A::StridedMatrix, Q::LQPackedQ)     (*)(adjoint(A), Q)
 @deprecate A_ldiv_B!(A::LQ{T}, B::StridedVecOrMat{T}) where {T} ldiv!(A, B)
 
 # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/qr.jl, to deprecate
-@deprecate A_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} mul!(A, B)
-@deprecate A_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix}    mul!(A, B)
-@deprecate A_mul_B!(A::QRPackedQ, B::AbstractVecOrMat)  mul!(A, B)
-@deprecate Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix}      mul!(adjoint(A), B)
-@deprecate Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix}   mul!(adjoint(A), B)
-@deprecate Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix}     mul!(adjoint(A), B)
-@deprecate Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix}  mul!(adjoint(A), B)
-@deprecate Ac_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) mul!(adjoint(A), B)
+@deprecate A_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} mul2!(A, B)
+@deprecate A_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix}    mul2!(A, B)
+@deprecate A_mul_B!(A::QRPackedQ, B::AbstractVecOrMat)  mul2!(A, B)
+@deprecate Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix}      mul2!(adjoint(A), B)
+@deprecate Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix}   mul2!(adjoint(A), B)
+@deprecate Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix}     mul2!(adjoint(A), B)
+@deprecate Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix}  mul2!(adjoint(A), B)
+@deprecate Ac_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) mul2!(adjoint(A), B)
 @deprecate Ac_mul_B(Q::AbstractQ, B::StridedVecOrMat)   (*)(adjoint(Q), B)
 @deprecate A_mul_Bc(Q::AbstractQ, B::StridedVecOrMat)   (*)(Q, adjoint(B))
 @deprecate Ac_mul_Bc(Q::AbstractQ, B::StridedVecOrMat)  (*)(adjoint(Q), adjoint(B))
-@deprecate A_mul_B!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix}  mul!(A, B)
-@deprecate A_mul_B!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix}     mul!(A, B)
-@deprecate A_mul_B!(A::StridedMatrix,Q::QRPackedQ)  mul!(A, Q)
-@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasReal}     mul!(A, adjoint(B))
-@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasComplex}  mul!(A, adjoint(B))
-@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasReal}    mul!(A, adjoint(B))
-@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasComplex} mul!(A, adjoint(B))
-@deprecate A_mul_Bc!(A::StridedMatrix,Q::QRPackedQ)     mul!(A, adjoint(Q))
+@deprecate A_mul_B!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix}  mul1!(A, B)
+@deprecate A_mul_B!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix}     mul1!(A, B)
+@deprecate A_mul_B!(A::StridedMatrix,Q::QRPackedQ)  mul1!(A, Q)
+@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasReal}     mul1!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasComplex}  mul1!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasReal}    mul1!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasComplex} mul1!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::StridedMatrix,Q::QRPackedQ)     mul1!(A, adjoint(Q))
 @deprecate A_mul_Bc(A::StridedMatrix, B::AbstractQ)     (*)(A, adjoint(B))
 @deprecate A_mul_Bc(rowvec::RowVector, B::AbstractQ)    (*)(rowvec, adjoint(B))
 @deprecate Ac_mul_B(A::StridedVecOrMat, Q::AbstractQ)   (*)(adjoint(A), Q)
@@ -987,20 +990,20 @@ export A_ldiv_B!,
 @deprecate At_mul_Bt(A::AbstractTriangular, B::AbstractTriangular)  (*)(transpose(A), transpose(B))
 @deprecate At_mul_Bt(A::AbstractTriangular, B::AbstractMatrix)      (*)(transpose(A), transpose(B))
 @deprecate At_mul_Bt(A::AbstractMatrix, B::AbstractTriangular)      (*)(transpose(A), transpose(B))
-@deprecate A_mul_Bc!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular})     mul!(A, adjoint(B))
-@deprecate A_mul_Bc!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular})     mul!(A, adjoint(B))
-@deprecate A_mul_Bt!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular})     mul!(A, transpose(B))
-@deprecate A_mul_Bt!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular})     mul!(A, transpose(B))
+@deprecate A_mul_Bc!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular})     mul1!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular})     mul1!(A, adjoint(B))
+@deprecate A_mul_Bt!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular})     mul1!(A, transpose(B))
+@deprecate A_mul_Bt!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular})     mul1!(A, transpose(B))
 @deprecate A_rdiv_Bc!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular})    rdiv!(A, adjoint(B))
 @deprecate A_rdiv_Bc!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular})    rdiv!(A, adjoint(B))
 @deprecate A_rdiv_Bt!(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular})    rdiv!(A, transpose(B))
 @deprecate A_rdiv_Bt!(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular})    rdiv!(A, transpose(B))
 @deprecate A_rdiv_B!(A::UpperTriangular, B::Union{UpperTriangular,UnitUpperTriangular})     rdiv!(A, B)
 @deprecate A_rdiv_B!(A::LowerTriangular, B::Union{LowerTriangular,UnitLowerTriangular})     rdiv!(A, B)
-@deprecate Ac_mul_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular)     mul!(adjoint(A), B)
-@deprecate Ac_mul_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular)     mul!(adjoint(A), B)
-@deprecate At_mul_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular)     mul!(transpose(A), B)
-@deprecate At_mul_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular)     mul!(transpose(A), B)
+@deprecate Ac_mul_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular)     mul2!(adjoint(A), B)
+@deprecate Ac_mul_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular)     mul2!(adjoint(A), B)
+@deprecate At_mul_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular)     mul2!(transpose(A), B)
+@deprecate At_mul_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular)     mul2!(transpose(A), B)
 @deprecate Ac_ldiv_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular)    ldiv!(adjoint(A), B)
 @deprecate Ac_ldiv_B!(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular)    ldiv!(adjoint(A), B)
 @deprecate At_ldiv_B!(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular)    ldiv!(transpose(A), B)
@@ -1025,30 +1028,30 @@ export A_ldiv_B!,
 @deprecate At_ldiv_B!(A::UpperTriangular, b::AbstractVector, x::AbstractVector = b)     ldiv!(transpose(A), b, x)
 @deprecate At_ldiv_B!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector = b) ldiv!(transpose(A), b, x)
 @deprecate At_ldiv_B!(A::LowerTriangular, b::AbstractVector, x::AbstractVector = b)     ldiv!(transpose(A), b, x)
-@deprecate A_mul_Bt!(A::StridedMatrix, B::UnitLowerTriangular)  mul!(A, transpose(B))
-@deprecate A_mul_Bt!(A::StridedMatrix, B::LowerTriangular)      mul!(A, transpose(B))
-@deprecate A_mul_Bt!(A::StridedMatrix, B::UnitUpperTriangular)  mul!(A, transpose(B))
-@deprecate A_mul_Bt!(A::StridedMatrix, B::UpperTriangular)      mul!(A, transpose(B))
-@deprecate A_mul_Bc!(A::StridedMatrix, B::UnitLowerTriangular)  mul!(A, adjoint(B))
-@deprecate A_mul_Bc!(A::StridedMatrix, B::LowerTriangular)      mul!(A, adjoint(B))
-@deprecate A_mul_Bc!(A::StridedMatrix, B::UnitUpperTriangular)  mul!(A, adjoint(B))
-@deprecate A_mul_Bc!(A::StridedMatrix, B::UpperTriangular)      mul!(A, adjoint(B))
-@deprecate A_mul_B!(A::StridedMatrix, B::UnitLowerTriangular)   mul!(A, B)
-@deprecate A_mul_B!(A::StridedMatrix, B::LowerTriangular)       mul!(A, B)
-@deprecate A_mul_B!(A::StridedMatrix, B::UnitUpperTriangular)   mul!(A, B)
-@deprecate A_mul_B!(A::StridedMatrix, B::UpperTriangular)       mul!(A, B)
-@deprecate At_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat)    mul!(transpose(A), B)
-@deprecate At_mul_B!(A::LowerTriangular, B::StridedVecOrMat)        mul!(transpose(A), B)
-@deprecate At_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat)    mul!(transpose(A), B)
-@deprecate At_mul_B!(A::UpperTriangular, B::StridedVecOrMat)        mul!(transpose(A), B)
-@deprecate Ac_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat)    mul!(adjoint(A), B)
-@deprecate Ac_mul_B!(A::LowerTriangular, B::StridedVecOrMat)        mul!(adjoint(A), B)
-@deprecate Ac_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat)    mul!(adjoint(A), B)
-@deprecate Ac_mul_B!(A::UpperTriangular, B::StridedVecOrMat)    mul!(adjoint(A), B)
-@deprecate A_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) mul!(A, B)
-@deprecate A_mul_B!(A::LowerTriangular, B::StridedVecOrMat)     mul!(A, B)
-@deprecate A_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) mul!(A, B)
-@deprecate A_mul_B!(A::UpperTriangular, B::StridedVecOrMat)     mul!(A, B)
+@deprecate A_mul_Bt!(A::StridedMatrix, B::UnitLowerTriangular)  mul1!(A, transpose(B))
+@deprecate A_mul_Bt!(A::StridedMatrix, B::LowerTriangular)      mul1!(A, transpose(B))
+@deprecate A_mul_Bt!(A::StridedMatrix, B::UnitUpperTriangular)  mul1!(A, transpose(B))
+@deprecate A_mul_Bt!(A::StridedMatrix, B::UpperTriangular)      mul1!(A, transpose(B))
+@deprecate A_mul_Bc!(A::StridedMatrix, B::UnitLowerTriangular)  mul1!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::StridedMatrix, B::LowerTriangular)      mul1!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::StridedMatrix, B::UnitUpperTriangular)  mul1!(A, adjoint(B))
+@deprecate A_mul_Bc!(A::StridedMatrix, B::UpperTriangular)      mul1!(A, adjoint(B))
+@deprecate A_mul_B!(A::StridedMatrix, B::UnitLowerTriangular)   mul1!(A, B)
+@deprecate A_mul_B!(A::StridedMatrix, B::LowerTriangular)       mul1!(A, B)
+@deprecate A_mul_B!(A::StridedMatrix, B::UnitUpperTriangular)   mul1!(A, B)
+@deprecate A_mul_B!(A::StridedMatrix, B::UpperTriangular)       mul1!(A, B)
+@deprecate At_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat)    mul2!(transpose(A), B)
+@deprecate At_mul_B!(A::LowerTriangular, B::StridedVecOrMat)        mul2!(transpose(A), B)
+@deprecate At_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat)    mul2!(transpose(A), B)
+@deprecate At_mul_B!(A::UpperTriangular, B::StridedVecOrMat)        mul2!(transpose(A), B)
+@deprecate Ac_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat)    mul2!(adjoint(A), B)
+@deprecate Ac_mul_B!(A::LowerTriangular, B::StridedVecOrMat)        mul2!(adjoint(A), B)
+@deprecate Ac_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat)    mul2!(adjoint(A), B)
+@deprecate Ac_mul_B!(A::UpperTriangular, B::StridedVecOrMat)    mul2!(adjoint(A), B)
+@deprecate A_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat) mul2!(A, B)
+@deprecate A_mul_B!(A::LowerTriangular, B::StridedVecOrMat)     mul2!(A, B)
+@deprecate A_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat) mul2!(A, B)
+@deprecate A_mul_B!(A::UpperTriangular, B::StridedVecOrMat)     mul2!(A, B)
 @deprecate A_mul_B!(C::AbstractVector  , A::AbstractTriangular, B::AbstractVector)      mul!(C, A, B)
 @deprecate A_mul_B!(C::AbstractMatrix  , A::AbstractTriangular, B::AbstractVecOrMat)    mul!(C, A, B)
 @deprecate A_mul_B!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat)    mul!(C, A, B)
@@ -1058,7 +1061,7 @@ export A_ldiv_B!,
 @deprecate At_mul_B!(C::AbstractVector  , A::AbstractTriangular, B::AbstractVector)     mul!(C, transpose(A), B)
 @deprecate At_mul_B!(C::AbstractMatrix  , A::AbstractTriangular, B::AbstractVecOrMat)   mul!(C, transpose(A), B)
 @deprecate At_mul_B!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat)   mul!(C, transpose(A), B)
-@deprecate A_mul_B!(A::Tridiagonal, B::AbstractTriangular)  mul!(A, B)
+@deprecate A_mul_B!(A::Tridiagonal, B::AbstractTriangular)  mul2!(A, B)
 @deprecate A_mul_B!(C::AbstractMatrix, A::AbstractTriangular, B::Tridiagonal)   mul!(C, A, B)
 @deprecate A_mul_B!(C::AbstractMatrix, A::Tridiagonal, B::AbstractTriangular)   mul!(C, A, B)
 @deprecate A_mul_Bt!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat)   mul!(C, A, transpose(B))
@@ -1110,22 +1113,22 @@ for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'),
                             (:UnitUpperTriangular, 'U', 'U'))
     @eval begin
         # Vector multiplication
-        @deprecate A_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat}     mul!(A, b)
-        @deprecate At_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat}    mul!(transpose(A), b)
-        @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasReal}     mul!(adjoint(A), b)
-        @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasComplex}  mul!(adjoint(A), b)
+        @deprecate A_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat}     mul2!(A, b)
+        @deprecate At_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat}    mul2!(transpose(A), b)
+        @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasReal}     mul2!(adjoint(A), b)
+        @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasComplex}  mul2!(adjoint(A), b)
 
         # Matrix multiplication
-        @deprecate A_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat}     mul!(A, B)
-        @deprecate A_mul_B!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat}     mul!(A, B)
+        @deprecate A_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat}     mul2!(A, B)
+        @deprecate A_mul_B!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat}     mul2!(A, B)
 
-        @deprecate At_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat}       mul!(transpose(A), B)
-        @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasComplex}     mul!(adjoint(A), B)
-        @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasReal}        mul!(adjoint(A), B)
+        @deprecate At_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat}       mul2!(transpose(A), B)
+        @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasComplex}     mul2!(adjoint(A), B)
+        @deprecate Ac_mul_B!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasReal}        mul2!(adjoint(A), B)
 
-        @deprecate A_mul_Bt!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat}    mul!(A, transpose(B))
-        @deprecate A_mul_Bc!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasComplex}  mul!(A, adjoint(B))
-        @deprecate A_mul_Bc!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasReal}     mul!(A, adjoint(B))
+        @deprecate A_mul_Bt!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat}    mul1!(A, transpose(B))
+        @deprecate A_mul_Bc!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasComplex}  mul1!(A, adjoint(B))
+        @deprecate A_mul_Bc!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasReal}     mul1!(A, adjoint(B))
 
         # Left division
         @deprecate A_ldiv_B!(A::$t{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat}  ldiv!(A, B)
@@ -1183,12 +1186,12 @@ end
 @deprecate A_mul_Bt(mat::AbstractMatrix, rowvec::RowVector)     (*)(mat, transpose(rowvec))
 
 # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/linalg/givens.jl, to deprecate
-@deprecate A_mul_Bc!(A::AbstractMatrix, R::Rotation)    mul!(A, adjoint(R))
-@deprecate A_mul_B!(R::Rotation, A::AbstractMatrix)     mul!(R, A)
-@deprecate A_mul_B!(G::Givens, R::Rotation)             mul!(G, R)
-@deprecate A_mul_Bc!(A::AbstractMatrix, G::Givens)      mul!(A, adjoint(G))
-@deprecate A_mul_B!(G::Givens, A::AbstractVecOrMat)     mul!(G, A)
-@deprecate A_mul_B!(G1::Givens, G2::Givens)             mul!(G1, G2)
+@deprecate A_mul_Bc!(A::AbstractMatrix, R::Rotation)    mul1!(A, adjoint(R))
+@deprecate A_mul_B!(R::Rotation, A::AbstractMatrix)     mul2!(R, A)
+@deprecate A_mul_B!(G::Givens, R::Rotation)             mul2!(G, R)
+@deprecate A_mul_Bc!(A::AbstractMatrix, G::Givens)      mul1!(A, adjoint(G))
+@deprecate A_mul_B!(G::Givens, A::AbstractVecOrMat)     mul2!(G, A)
+@deprecate A_mul_B!(G1::Givens, G2::Givens)             G1 * G2
 @deprecate A_mul_Bc(A::AbstractVecOrMat{T}, R::AbstractRotation{S}) where {T,S}     (*)(A, adjoint(R))
 
 
@@ -1254,7 +1257,10 @@ norm(tv::RowVector) = norm(rvtranspose(tv))
 *(A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:RowVector}) = A * rvtranspose(B.parent)
 *(A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:RowVector}) = A * rvadjoint(B.parent)
 
-
+@deprecate *(A::LQ,B::QR) A*Matrix(B)
+@deprecate *(A::QR,B::LQ) A*Matrix(B)
+@deprecate *(A::Adjoint{<:Any,<:LQ}, B::LQ) A*Matrix(B)
+@deprecate *(A::LQ, B::Adjoint{<:Any,<:LQ}) A*Matrix(B)
 
 # PR #25184. Use getproperty instead of getindex for Factorizations
 function getindex(F::Factorization, s::Symbol)
diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl
index 8e3f06e1000d1..e90820b949eed 100644
--- a/stdlib/LinearAlgebra/src/diagonal.jl
+++ b/stdlib/LinearAlgebra/src/diagonal.jl
@@ -151,38 +151,40 @@ end
 (*)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .* Db.diag)
 (*)(D::Diagonal, V::AbstractVector) = D.diag .* V
 
-(*)(A::AbstractTriangular, D::Diagonal) = mul!(copy(A), D)
-(*)(D::Diagonal, B::AbstractTriangular) = mul!(D, copy(B))
+(*)(A::AbstractTriangular, D::Diagonal) = mul1!(copy(A), D)
+(*)(D::Diagonal, B::AbstractTriangular) = mul2!(D, copy(B))
 
 (*)(A::AbstractMatrix, D::Diagonal) =
     scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), A, D.diag)
 (*)(D::Diagonal, A::AbstractMatrix) =
     scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), D.diag, A)
 
-mul!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) = typeof(A)(mul!(A.data, D))
-function mul!(A::UnitLowerTriangular, D::Diagonal)
-    mul!(A.data, D)
+
+mul1!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) = typeof(A)(mul1!(A.data, D))
+function mul1!(A::UnitLowerTriangular, D::Diagonal)
+    mul1!(A.data, D)
     for i = 1:size(A, 1)
         A.data[i,i] = D.diag[i]
     end
     LowerTriangular(A.data)
 end
-function mul!(A::UnitUpperTriangular, D::Diagonal)
-    mul!(A.data, D)
+function mul1!(A::UnitUpperTriangular, D::Diagonal)
+    mul1!(A.data, D)
     for i = 1:size(A, 1)
         A.data[i,i] = D.diag[i]
     end
     UpperTriangular(A.data)
 end
-function mul!(D::Diagonal, B::UnitLowerTriangular)
-    mul!(D, B.data)
+
+function mul2!(D::Diagonal, B::UnitLowerTriangular)
+    mul2!(D, B.data)
     for i = 1:size(B, 1)
         B.data[i,i] = D.diag[i]
     end
     LowerTriangular(B.data)
 end
-function mul!(D::Diagonal, B::UnitUpperTriangular)
-    mul!(D, B.data)
+function mul2!(D::Diagonal, B::UnitUpperTriangular)
+    mul2!(D, B.data)
     for i = 1:size(B, 1)
         B.data[i,i] = D.diag[i]
     end
@@ -190,40 +192,40 @@ function mul!(D::Diagonal, B::UnitUpperTriangular)
 end
 
 *(D::Adjoint{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(adjoint.(D.parent.diag) .* B.diag)
-*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) = mul!(copy(A), D)
+*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) = mul1!(copy(A), D)
 function *(adjA::Adjoint{<:Any,<:AbstractMatrix}, D::Diagonal)
     A = adjA.parent
     Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
     adjoint!(Ac, A)
-    mul!(Ac, D)
+    mul1!(Ac, D)
 end
 
 *(D::Transpose{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(transpose.(D.parent.diag) .* B.diag)
-*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) = mul!(copy(A), D)
+*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) = mul1!(copy(A), D)
 function *(transA::Transpose{<:Any,<:AbstractMatrix}, D::Diagonal)
     A = transA.parent
     At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
     transpose!(At, A)
-    mul!(At, D)
+    mul1!(At, D)
 end
 
 *(D::Diagonal, B::Adjoint{<:Any,<:Diagonal}) = Diagonal(D.diag .* adjoint.(B.parent.diag))
-*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) = mul!(D, collect(B))
-*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; mul!(Array(D), adjoint(Q)))
+*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) = mul2!(D, collect(B))
+*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; mul1!(Array(D), adjoint(Q)))
 function *(D::Diagonal, adjA::Adjoint{<:Any,<:AbstractMatrix})
     A = adjA.parent
     Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
     adjoint!(Ac, A)
-    mul!(D, Ac)
+    mul2!(D, Ac)
 end
 
 *(D::Diagonal, B::Transpose{<:Any,<:Diagonal}) = Diagonal(D.diag .* transpose.(B.parent.diag))
-*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) = mul!(D, copy(B))
+*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) = mul2!(D, copy(B))
 function *(D::Diagonal, transA::Transpose{<:Any,<:AbstractMatrix})
     A = transA.parent
     At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
     transpose!(At, A)
-    mul!(D, At)
+    mul2!(D, At)
 end
 
 *(D::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) =
@@ -231,27 +233,15 @@ end
 *(D::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) =
     Diagonal(transpose.(D.parent.diag) .* transpose.(B.parent.diag))
 
-mul!(A::Diagonal, B::Diagonal) = throw(MethodError(mul!, (A, B)))
-mul!(A::QRPackedQ, D::Diagonal) = throw(MethodError(mul!, (A, D)))
-mul!(A::QRPackedQ, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::QRPackedQ, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Diagonal) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Diagonal, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Diagonal, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Transpose{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Transpose{<:Any,<:Diagonal}, B::Diagonal) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:Diagonal}, B::Diagonal) = throw(MethodError(mul!, (A, B)))
-mul!(A::Diagonal, B::AbstractMatrix)  = scale!(A.diag, B)
-mul!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = adjA.parent; scale!(conj(A.diag),B))
-mul!(transA::Transpose{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = transA.parent; scale!(A.diag,B))
-mul!(A::AbstractMatrix, B::Diagonal)  = scale!(A,B.diag)
-mul!(A::AbstractMatrix, adjB::Adjoint{<:Any,<:Diagonal}) = (B = adjB.parent; scale!(A,conj(B.diag)))
-mul!(A::AbstractMatrix, transB::Transpose{<:Any,<:Diagonal}) = (B = transB.parent; scale!(A,B.diag))
+mul1!(A::Diagonal,B::Diagonal) = Diagonal(A.diag .*= B.diag)
+mul2!(A::Diagonal,B::Diagonal) = Diagonal(B.diag .= A.diag .* B.diag)
+
+mul2!(A::Diagonal, B::AbstractMatrix)  = scale!(A.diag, B)
+mul2!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = adjA.parent; scale!(conj(A.diag),B))
+mul2!(transA::Transpose{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = transA.parent; scale!(A.diag,B))
+mul1!(A::AbstractMatrix, B::Diagonal)  = scale!(A,B.diag)
+mul1!(A::AbstractMatrix, adjB::Adjoint{<:Any,<:Diagonal}) = (B = adjB.parent; scale!(A,conj(B.diag)))
+mul1!(A::AbstractMatrix, transB::Transpose{<:Any,<:Diagonal}) = (B = transB.parent; scale!(A,B.diag))
 
 # Get ambiguous method if try to unify AbstractVector/AbstractMatrix here using AbstractVecOrMat
 mul!(out::AbstractVector, A::Diagonal, in::AbstractVector) = out .= A.diag .* in
@@ -287,6 +277,7 @@ mul!(C::AbstractMatrix, A::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Rea
 
 
 (/)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag ./ Db.diag)
+
 function ldiv!(D::Diagonal{T}, v::AbstractVector{T}) where {T}
     if length(v) != length(D.diag)
         throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $(length(v)) rows"))
@@ -316,6 +307,7 @@ function ldiv!(D::Diagonal{T}, V::AbstractMatrix{T}) where {T}
     V
 end
 
+
 ldiv!(adjD::Adjoint{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} =
     (D = adjD.parent; ldiv!(conj(D), B))
 ldiv!(transD::Transpose{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} =
@@ -339,6 +331,7 @@ function rdiv!(A::AbstractMatrix{T}, D::Diagonal{T}) where {T}
     A
 end
 
+
 rdiv!(A::AbstractMatrix{T}, adjD::Adjoint{<:Any,<:Diagonal{T}}) where {T} =
     (D = adjD.parent; rdiv!(A, conj(D)))
 rdiv!(A::AbstractMatrix{T}, transD::Transpose{<:Any,<:Diagonal{T}}) where {T} =
diff --git a/stdlib/LinearAlgebra/src/givens.jl b/stdlib/LinearAlgebra/src/givens.jl
index 4c84aac293974..9b508d5e1677c 100644
--- a/stdlib/LinearAlgebra/src/givens.jl
+++ b/stdlib/LinearAlgebra/src/givens.jl
@@ -7,14 +7,14 @@ transpose(R::AbstractRotation) = error("transpose not implemented for $(typeof(R
 
 function *(R::AbstractRotation{T}, A::AbstractVecOrMat{S}) where {T,S}
     TS = typeof(zero(T)*zero(S) + zero(T)*zero(S))
-    mul!(convert(AbstractRotation{TS}, R), TS == S ? copy(A) : convert(AbstractArray{TS}, A))
+    mul2!(convert(AbstractRotation{TS}, R), TS == S ? copy(A) : convert(AbstractArray{TS}, A))
 end
 *(A::AbstractVector, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)
 *(A::AbstractMatrix, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)
 function _absvecormat_mul_adjrot(A::AbstractVecOrMat{T}, adjR::Adjoint{<:Any,<:AbstractRotation{S}}) where {T,S}
     R = adjR.parent
     TS = typeof(zero(T)*zero(S) + zero(T)*zero(S))
-    mul!(TS == T ? copy(A) : convert(AbstractArray{TS}, A), adjoint(convert(AbstractRotation{TS}, R)))
+    mul1!(TS == T ? copy(A) : convert(AbstractArray{TS}, A), adjoint(convert(AbstractRotation{TS}, R)))
 end
 """
     LinearAlgebra.Givens(i1,i2,c,s) -> G
@@ -325,10 +325,7 @@ function getindex(G::Givens, i::Integer, j::Integer)
     end
 end
 
-
-mul!(G1::Givens, G2::Givens) = error("Operation not supported. Consider *")
-
-function mul!(G::Givens, A::AbstractVecOrMat)
+function mul2!(G::Givens, A::AbstractVecOrMat)
     m, n = size(A, 1), size(A, 2)
     if G.i2 > m
         throw(DimensionMismatch("column indices for rotation are outside the matrix"))
@@ -340,7 +337,7 @@ function mul!(G::Givens, A::AbstractVecOrMat)
     end
     return A
 end
-function mul!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens})
+function mul1!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens})
     G = adjG.parent
     m, n = size(A, 1), size(A, 2)
     if G.i2 > n
@@ -353,20 +350,21 @@ function mul!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens})
     end
     return A
 end
-function mul!(G::Givens, R::Rotation)
+
+function mul2!(G::Givens, R::Rotation)
     push!(R.rotations, G)
     return R
 end
-function mul!(R::Rotation, A::AbstractMatrix)
+function mul2!(R::Rotation, A::AbstractMatrix)
     @inbounds for i = 1:length(R.rotations)
-        mul!(R.rotations[i], A)
+        mul2!(R.rotations[i], A)
     end
     return A
 end
-function mul!(A::AbstractMatrix, adjR::Adjoint{<:Any,<:Rotation})
+function mul1!(A::AbstractMatrix, adjR::Adjoint{<:Any,<:Rotation})
     R = adjR.parent
     @inbounds for i = 1:length(R.rotations)
-        mul!(A, adjoint(R.rotations[i]))
+        mul1!(A, adjoint(R.rotations[i]))
     end
     return A
 end
@@ -388,14 +386,3 @@ end
 # dismabiguation methods: *(Diag/AbsTri, Adj of AbstractRotation)
 *(A::Diagonal, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B)
 *(A::AbstractTriangular, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B)
-# moar disambiguation
-mul!(A::QRPackedQ, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))
-mul!(A::QRPackedQ, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Diagonal, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Diagonal, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Transpose{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
-mul!(A::Transpose{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))
diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl
index 0d2c0795efd08..18a098f1fb8f7 100644
--- a/stdlib/LinearAlgebra/src/hessenberg.jl
+++ b/stdlib/LinearAlgebra/src/hessenberg.jl
@@ -73,7 +73,7 @@ function getindex(A::HessenbergQ, i::Integer, j::Integer)
     x[i] = 1
     y = zeros(eltype(A), size(A, 2))
     y[j] = 1
-    return dot(x, mul!(A, y))
+    return dot(x, mul2!(A, y))
 end
 
 ## reconstruct the original matrix
@@ -84,31 +84,30 @@ AbstractArray(F::Hessenberg) = AbstractMatrix(F)
 Matrix(F::Hessenberg) = Array(AbstractArray(F))
 Array(F::Hessenberg) = Matrix(F)
 
-mul!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
+mul2!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
     LAPACK.ormhr!('L', 'N', 1, size(Q.factors, 1), Q.factors, Q.τ, X)
-mul!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} =
+mul1!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} =
     LAPACK.ormhr!('R', 'N', 1, size(Q.factors, 1), Q.factors, Q.τ, X)
-mul!(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
+mul2!(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
     (Q = adjQ.parent; LAPACK.ormhr!('L', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X))
-mul!(X::StridedMatrix{T}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T<:BlasFloat} =
+mul1!(X::StridedMatrix{T}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T<:BlasFloat} =
     (Q = adjQ.parent; LAPACK.ormhr!('R', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X))
 
-
 function (*)(Q::HessenbergQ{T}, X::StridedVecOrMat{S}) where {T,S}
     TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
-    return mul!(Q, copy_oftype(X, TT))
+    return mul2!(Q, copy_oftype(X, TT))
 end
 function (*)(X::StridedVecOrMat{S}, Q::HessenbergQ{T}) where {T,S}
     TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
-    return mul!(copy_oftype(X, TT), Q)
+    return mul1!(copy_oftype(X, TT), Q)
 end
 function *(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{S}) where {T,S}
     Q = adjQ.parent
     TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
-    return mul!(adjoint(Q), copy_oftype(X, TT))
+    return mul2!(adjoint(Q), copy_oftype(X, TT))
 end
 function *(X::StridedVecOrMat{S}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T,S}
     Q = adjQ.parent
     TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
-    return mul!(copy_oftype(X, TT), adjoint(Q))
+    return mul1!(copy_oftype(X, TT), adjoint(Q))
 end
diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl
index 254edae9d6036..7197c3cfb9d8b 100644
--- a/stdlib/LinearAlgebra/src/lq.jl
+++ b/stdlib/LinearAlgebra/src/lq.jl
@@ -57,7 +57,7 @@ function lq(A::Union{Number,AbstractMatrix}; full::Bool = false, thin::Union{Boo
     end
     F = lqfact(A)
     L, Q = F.L, F.Q
-    return L, !full ? Array(Q) : mul!(Q, Matrix{eltype(Q)}(I, size(Q.factors, 2), size(Q.factors, 2)))
+    return L, !full ? Array(Q) : mul2!(Q, Matrix{eltype(Q)}(I, size(Q.factors, 2), size(Q.factors, 2)))
 end
 
 copy(A::LQ) = LQ(copy(A.factors), copy(A.τ))
@@ -88,7 +88,7 @@ end
 Base.propertynames(F::LQ, private::Bool=false) = append!([:L,:Q], private ? fieldnames(typeof(F)) : Symbol[])
 
 getindex(A::LQPackedQ, i::Integer, j::Integer) =
-    mul!(A, setindex!(zeros(eltype(A), size(A, 2)), 1, j))[i]
+    mul2!(A, setindex!(zeros(eltype(A), size(A, 2)), 1, j))[i]
 
 function show(io::IO, C::LQ)
     println(io, "$(typeof(C)) with factors L and Q:")
@@ -122,47 +122,34 @@ end
 
 
 ## Multiplication by LQ
-mul!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
-    A.L * LAPACK.ormlq!('L', 'N', A.factors, A.τ, B)
-mul!(A::LQ{T}, B::QR{T}) where {T<:BlasFloat} =
-    A.L * LAPACK.ormlq!('L', 'N', A.factors, A.τ, Matrix(B))
-mul!(A::QR{T}, B::LQ{T}) where {T<:BlasFloat} =
-    mul!(zeros(eltype(A), size(A)), Matrix(A), Matrix(B))
+mul2!(A::LQ, B::StridedVecOrMat) =
+    mul2!(LowerTriangular(A.L), mul2!(A.Q, B))
 function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB}
     TAB = promote_type(TA, TB)
-    mul!(Factorization{TAB}(A), copy_oftype(B, TAB))
+    mul2!(Factorization{TAB}(A), copy_oftype(B, TAB))
 end
-function *(A::LQ{TA},B::QR{TB}) where {TA,TB}
-    TAB = promote_type(TA, TB)
-    mul!(Factorization{TAB}(A), Factorization{TAB}(B))
-end
-function *(A::QR{TA},B::LQ{TB}) where {TA,TB}
-    TAB = promote_type(TA, TB)
-    mul!(Factorization{TAB}(A), Factorization{TAB}(B))
-end
-*(A::Adjoint{<:Any,<:LQ}, B::LQ) = copy(A) * B
-*(A::LQ, B::Adjoint{<:Any,<:LQ}) = A * copy(B)
 
 ## Multiplication by Q
 ### QB
-mul!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B)
+mul2!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B)
 function (*)(A::LQPackedQ, B::StridedVecOrMat)
     TAB = promote_type(eltype(A), eltype(B))
-    mul!(AbstractMatrix{TAB}(A), copy_oftype(B, TAB))
+    mul2!(AbstractMatrix{TAB}(A), copy_oftype(B, TAB))
 end
 
 ### QcB
-mul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} =
+mul2!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} =
     (A = adjA.parent; LAPACK.ormlq!('L','T',A.factors,A.τ,B))
-mul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} =
+mul2!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} =
     (A = adjA.parent; LAPACK.ormlq!('L','C',A.factors,A.τ,B))
+
 function *(adjA::Adjoint{<:Any,<:LQPackedQ}, B::StridedVecOrMat)
     A = adjA.parent
     TAB = promote_type(eltype(A), eltype(B))
     if size(B,1) == size(A.factors,2)
-        mul!(adjoint(AbstractMatrix{TAB}(A)), copy_oftype(B, TAB))
+        mul2!(adjoint(AbstractMatrix{TAB}(A)), copy_oftype(B, TAB))
     elseif size(B,1) == size(A.factors,1)
-        mul!(adjoint(AbstractMatrix{TAB}(A)), [B; zeros(TAB, size(A.factors, 2) - size(A.factors, 1), size(B, 2))])
+        mul2!(adjoint(AbstractMatrix{TAB}(A)), [B; zeros(TAB, size(A.factors, 2) - size(A.factors, 1), size(B, 2))])
     else
         throw(DimensionMismatch("first dimension of B, $(size(B,1)), must equal one of the dimensions of A, $(size(A))"))
     end
@@ -174,14 +161,14 @@ function *(A::LQPackedQ, adjB::Adjoint{<:Any,<:StridedVecOrMat})
     TAB = promote_type(eltype(A), eltype(B))
     BB = similar(B, TAB, (size(B, 2), size(B, 1)))
     adjoint!(BB, B)
-    return mul!(A, BB)
+    return mul2!(A, BB)
 end
 function *(adjA::Adjoint{<:Any,<:LQPackedQ}, adjB::Adjoint{<:Any,<:StridedVecOrMat})
     A, B = adjA.parent, adjB.parent
     TAB = promote_type(eltype(A), eltype(B))
     BB = similar(B, TAB, (size(B, 2), size(B, 1)))
     adjoint!(BB, B)
-    return mul!(adjoint(A), BB)
+    return mul2!(adjoint(A), BB)
 end
 
 # in-place right-application of LQPackedQs
@@ -189,11 +176,11 @@ end
 # match the number of columns (nQ) of the LQPackedQ (Q) (necessary for in-place
 # operation, and the underlying LAPACK routine (ormlq) treats the implicit Q
 # as its (nQ-by-nQ) square form)
-mul!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} =
+mul1!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} =
     LAPACK.ormlq!('R', 'N', B.factors, B.τ, A)
-mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasReal} =
+mul1!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasReal} =
     (B = adjB.parent; LAPACK.ormlq!('R', 'T', B.factors, B.τ, A))
-mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasComplex} =
+mul1!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasComplex} =
     (B = adjB.parent; LAPACK.ormlq!('R', 'C', B.factors, B.τ, A))
 
 # out-of-place right application of LQPackedQs
@@ -211,13 +198,13 @@ mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasCom
 function *(A::StridedVecOrMat, adjQ::Adjoint{<:Any,<:LQPackedQ})
     Q = adjQ.parent
     TR = promote_type(eltype(A), eltype(Q))
-    return mul!(copy_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q)))
+    return mul1!(copy_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q)))
 end
 function *(adjA::Adjoint{<:Any,<:StridedMatrix}, adjQ::Adjoint{<:Any,<:LQPackedQ})
     A, Q = adjA.parent, adjQ.parent
     TR = promote_type(eltype(A), eltype(Q))
     C = adjoint!(similar(A, TR, reverse(size(A))), A)
-    return mul!(C, adjoint(AbstractMatrix{TR}(Q)))
+    return mul1!(C, adjoint(AbstractMatrix{TR}(Q)))
 end
 #
 # (2) the inner dimension in the multiplication is the LQPackedQ's first dimension.
@@ -242,7 +229,7 @@ function *(A::StridedVecOrMat, Q::LQPackedQ)
     else
         _rightappdimmismatch("columns")
     end
-    return mul!(C, AbstractMatrix{TR}(Q))
+    return mul1!(C, AbstractMatrix{TR}(Q))
 end
 function *(adjA::Adjoint{<:Any,<:StridedMatrix}, Q::LQPackedQ)
     A = adjA.parent
@@ -255,7 +242,7 @@ function *(adjA::Adjoint{<:Any,<:StridedMatrix}, Q::LQPackedQ)
     else
         _rightappdimmismatch("rows")
     end
-    return mul!(C, AbstractMatrix{TR}(Q))
+    return mul1!(C, AbstractMatrix{TR}(Q))
 end
 _rightappdimmismatch(rowsorcols) =
     throw(DimensionMismatch(string("the number of $(rowsorcols) of the matrix on the left ",
@@ -291,6 +278,6 @@ end
 
 
 function ldiv!(A::LQ{T}, B::StridedVecOrMat{T}) where T
-    mul!(adjoint(A.Q), ldiv!(LowerTriangular(A.L),B))
+    mul2!(adjoint(A.Q), ldiv!(LowerTriangular(A.L),B))
     return B
 end
diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl
index 3f3bf8143db3b..ae2579715ddf5 100644
--- a/stdlib/LinearAlgebra/src/matmul.jl
+++ b/stdlib/LinearAlgebra/src/matmul.jl
@@ -190,12 +190,18 @@ julia> Y
 mul!(C::AbstractMatrix, A::AbstractVecOrMat, B::AbstractVecOrMat) = generic_matmatmul!(C, 'N', 'N', A, B)
 
 """
-    mul!(A, B)
+    mul1!(A, B)
 
-Calculate the matrix-matrix product ``AB``, overwriting one of `A` or `B` (but not both),
-and return the result (the overwritten argument).
+Calculate the matrix-matrix product ``AB``, overwriting `A`, and return the result.
 """
-mul!(A, B)
+mul1!(A, B)
+
+"""
+    mul2!(A, B)
+
+Calculate the matrix-matrix product ``AB``, overwriting `B`, and return the result.
+"""
+mul2!(A, B)
 
 function *(transA::Transpose{<:Any,<:AbstractMatrix}, B::AbstractMatrix)
     A = transA.parent
diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl
index dae11f197942b..2fdc2473e7c80 100644
--- a/stdlib/LinearAlgebra/src/qr.jl
+++ b/stdlib/LinearAlgebra/src/qr.jl
@@ -329,13 +329,13 @@ function _qr(A::Union{Number,AbstractMatrix}, ::Val{false}; full::Bool = false)
     F = qrfact(A, Val(false))
     Q, R = F.Q, F.R
     sQf1 = size(Q.factors, 1)
-    return (!full ? Array(Q) : mul!(Q, Matrix{eltype(Q)}(I, sQf1, sQf1))), R
+    return (!full ? Array(Q) : mul2!(Q, Matrix{eltype(Q)}(I, sQf1, sQf1))), R
 end
 function _qr(A::Union{Number, AbstractMatrix}, ::Val{true}; full::Bool = false)
     F = qrfact(A, Val(true))
     Q, R, p = F.Q, F.R, F.p
     sQf1 = size(Q.factors, 1)
-    return (!full ? Array(Q) : mul!(Q, Matrix{eltype(Q)}(I, sQf1, sQf1))), R, p
+    return (!full ? Array(Q) : mul2!(Q, Matrix{eltype(Q)}(I, sQf1, sQf1))), R, p
 end
 
 """
@@ -506,7 +506,7 @@ AbstractMatrix{T}(Q::QRPackedQ) where {T} = QRPackedQ{T}(Q)
 QRCompactWYQ{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ(convert(AbstractMatrix{S}, Q.factors), convert(AbstractMatrix{S}, Q.T))
 AbstractMatrix{S}(Q::QRCompactWYQ{S}) where {S} = Q
 AbstractMatrix{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ{S}(Q)
-Matrix(A::AbstractQ{T}) where {T} = mul!(A, Matrix{T}(I, size(A.factors, 1), min(size(A.factors)...)))
+Matrix(A::AbstractQ{T}) where {T} = mul2!(A, Matrix{T}(I, size(A.factors, 1), min(size(A.factors)...)))
 Array(A::AbstractQ) = Matrix(A)
 
 size(A::Union{QR,QRCompactWY,QRPivoted}, dim::Integer) = size(getfield(A, :factors), dim)
@@ -520,16 +520,16 @@ function getindex(A::AbstractQ, i::Integer, j::Integer)
     x[i] = 1
     y = zeros(eltype(A), size(A, 2))
     y[j] = 1
-    return dot(x, mul!(A, y))
+    return dot(x, mul2!(A, y))
 end
 
 ## Multiplication by Q
 ### QB
-mul!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} =
+mul2!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} =
     LAPACK.gemqrt!('L','N',A.factors,A.T,B)
-mul!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} =
+mul2!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} =
     LAPACK.ormqr!('L','N',A.factors,A.τ,B)
-function mul!(A::QRPackedQ, B::AbstractVecOrMat)
+function mul2!(A::QRPackedQ, B::AbstractVecOrMat)
     mA, nA = size(A.factors)
     mB, nB = size(B,1), size(B,2)
     if mA != mB
@@ -564,7 +564,7 @@ function (*)(A::AbstractQ, b::StridedVector)
     else
         throw(DimensionMismatch("vector must have length either $(size(A.factors, 1)) or $(size(A.factors, 2))"))
     end
-    mul!(Anew, bnew)
+    mul2!(Anew, bnew)
 end
 function (*)(A::AbstractQ, B::StridedMatrix)
     TAB = promote_type(eltype(A), eltype(B))
@@ -576,19 +576,19 @@ function (*)(A::AbstractQ, B::StridedMatrix)
     else
         throw(DimensionMismatch("first dimension of matrix must have size either $(size(A.factors, 1)) or $(size(A.factors, 2))"))
     end
-    mul!(Anew, Bnew)
+    mul2!(Anew, Bnew)
 end
 
 ### QcB
-mul!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} =
+mul2!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} =
     (A = adjA.parent; LAPACK.gemqrt!('L','T',A.factors,A.T,B))
-mul!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} =
+mul2!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} =
     (A = adjA.parent; LAPACK.gemqrt!('L','C',A.factors,A.T,B))
-mul!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} =
+mul2!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} =
     (A = adjA.parent; LAPACK.ormqr!('L','T',A.factors,A.τ,B))
-mul!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} =
+mul2!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} =
     (A = adjA.parent; LAPACK.ormqr!('L','C',A.factors,A.τ,B))
-function mul!(adjA::Adjoint{<:Any,<:QRPackedQ}, B::AbstractVecOrMat)
+function mul2!(adjA::Adjoint{<:Any,<:QRPackedQ}, B::AbstractVecOrMat)
     A = adjA.parent
     mA, nA = size(A.factors)
     mB, nB = size(B,1), size(B,2)
@@ -616,7 +616,7 @@ end
 function *(adjQ::Adjoint{<:Any,<:AbstractQ}, B::StridedVecOrMat)
     Q = adjQ.parent
     TQB = promote_type(eltype(Q), eltype(B))
-    return mul!(adjoint(convert(AbstractMatrix{TQB}, Q)), copy_oftype(B, TQB))
+    return mul2!(adjoint(convert(AbstractMatrix{TQB}, Q)), copy_oftype(B, TQB))
 end
 
 ### QBc/QcBc
@@ -625,22 +625,22 @@ function *(Q::AbstractQ, adjB::Adjoint{<:Any,<:StridedVecOrMat})
     TQB = promote_type(eltype(Q), eltype(B))
     Bc = similar(B, TQB, (size(B, 2), size(B, 1)))
     adjoint!(Bc, B)
-    return mul!(convert(AbstractMatrix{TQB}, Q), Bc)
+    return mul2!(convert(AbstractMatrix{TQB}, Q), Bc)
 end
 function *(adjQ::Adjoint{<:Any,<:AbstractQ}, adjB::Adjoint{<:Any,<:StridedVecOrMat})
     Q, B = adjQ.parent, adjB.parent
     TQB = promote_type(eltype(Q), eltype(B))
     Bc = similar(B, TQB, (size(B, 2), size(B, 1)))
     adjoint!(Bc, B)
-    return mul!(adjoint(convert(AbstractMatrix{TQB}, Q)), Bc)
+    return mul2!(adjoint(convert(AbstractMatrix{TQB}, Q)), Bc)
 end
 
 ### AQ
-mul!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} =
+mul1!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} =
     LAPACK.gemqrt!('R','N', B.factors, B.T, A)
-mul!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} =
+mul1!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} =
     LAPACK.ormqr!('R', 'N', B.factors, B.τ, A)
-function mul!(A::StridedMatrix,Q::QRPackedQ)
+function mul1!(A::StridedMatrix,Q::QRPackedQ)
     mQ, nQ = size(Q.factors)
     mA, nA = size(A,1), size(A,2)
     if nA != mQ
@@ -667,19 +667,20 @@ end
 
 function (*)(A::StridedMatrix, Q::AbstractQ)
     TAQ = promote_type(eltype(A), eltype(Q))
-    return mul!(copy_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q))
+
+    return mul1!(copy_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q))
 end
 
 ### AQc
-mul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasReal} =
+mul1!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasReal} =
     (B = adjB.parent; LAPACK.gemqrt!('R','T',B.factors,B.T,A))
-mul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasComplex} =
+mul1!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasComplex} =
     (B = adjB.parent; LAPACK.gemqrt!('R','C',B.factors,B.T,A))
-mul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasReal} =
+mul1!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasReal} =
     (B = adjB.parent; LAPACK.ormqr!('R','T',B.factors,B.τ,A))
-mul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasComplex} =
+mul1!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasComplex} =
     (B = adjB.parent; LAPACK.ormqr!('R','C',B.factors,B.τ,A))
-function mul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ})
+function mul1!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ})
     Q = adjQ.parent
     mQ, nQ = size(Q.factors)
     mA, nA = size(A,1), size(A,2)
@@ -711,9 +712,9 @@ function *(A::StridedMatrix, adjB::Adjoint{<:Any,<:AbstractQ})
     if size(A,2) == size(B.factors, 1)
         AA = similar(A, TAB, size(A))
         copyto!(AA, A)
-        return mul!(AA, adjoint(BB))
+        return mul1!(AA, adjoint(BB))
     elseif size(A,2) == size(B.factors,2)
-        return mul!([A zeros(TAB, size(A, 1), size(B.factors, 1) - size(B.factors, 2))], adjoint(BB))
+        return mul1!([A zeros(TAB, size(A, 1), size(B.factors, 1) - size(B.factors, 2))], adjoint(BB))
     else
         throw(DimensionMismatch("matrix A has dimensions $(size(A)) but matrix B has dimensions $(size(B))"))
     end
@@ -727,20 +728,20 @@ function *(adjA::Adjoint{<:Any,<:StridedVecOrMat}, Q::AbstractQ)
     TAQ = promote_type(eltype(A), eltype(Q))
     Ac = similar(A, TAQ, (size(A, 2), size(A, 1)))
     adjoint!(Ac, A)
-    return mul!(Ac, convert(AbstractMatrix{TAQ}, Q))
+    return mul1!(Ac, convert(AbstractMatrix{TAQ}, Q))
 end
 function *(adjA::Adjoint{<:Any,<:StridedVecOrMat}, adjQ::Adjoint{<:Any,<:AbstractQ})
     A, Q = adjA.parent, adjQ.parent
     TAQ = promote_type(eltype(A), eltype(Q))
     Ac = similar(A, TAQ, (size(A, 2), size(A, 1)))
     adjoint!(Ac, A)
-    return mul!(Ac, adjoint(convert(AbstractMatrix{TAQ}, Q)))
+    return mul1!(Ac, adjoint(convert(AbstractMatrix{TAQ}, Q)))
 end
 
 ldiv!(A::QRCompactWY{T}, b::StridedVector{T}) where {T<:BlasFloat} =
-    (ldiv!(UpperTriangular(A.R), view(mul!(adjoint(A.Q), b), 1:size(A, 2))); b)
+    (ldiv!(UpperTriangular(A.R), view(mul2!(adjoint(A.Q), b), 1:size(A, 2))); b)
 ldiv!(A::QRCompactWY{T}, B::StridedMatrix{T}) where {T<:BlasFloat} =
-    (ldiv!(UpperTriangular(A.R), view(mul!(adjoint(A.Q), B), 1:size(A, 2), 1:size(B, 2))); B)
+    (ldiv!(UpperTriangular(A.R), view(mul2!(adjoint(A.Q), B), 1:size(A, 2), 1:size(B, 2))); B)
 
 # Julia implementation similar to xgelsy
 function ldiv!(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real) where T<:BlasFloat
@@ -772,7 +773,7 @@ function ldiv!(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real) where T<:BlasF
         rnk += 1
     end
     C, τ = LAPACK.tzrzf!(A.factors[1:rnk,:])
-    ldiv!(UpperTriangular(C[1:rnk,1:rnk]),view(mul!(adjoint(A.Q), view(B, 1:mA, 1:nrhs)), 1:rnk, 1:nrhs))
+    ldiv!(UpperTriangular(C[1:rnk,1:rnk]),view(mul2!(adjoint(A.Q), view(B, 1:mA, 1:nrhs)), 1:rnk, 1:nrhs))
     B[rnk+1:end,:] = zero(T)
     LAPACK.ormrz!('L', eltype(B)<:Complex ? 'C' : 'T', C, τ, view(B,1:nA,1:nrhs))
     B[1:nA,:] = view(B, 1:nA, :)[invperm(A.p),:]
@@ -786,7 +787,7 @@ function ldiv!(A::QR{T}, B::StridedMatrix{T}) where T
     m, n = size(A)
     minmn = min(m,n)
     mB, nB = size(B)
-    mul!(adjoint(A.Q), view(B, 1:m, :))
+    mul2!(adjoint(A.Q), view(B, 1:m, :))
     R = A.R
     @inbounds begin
         if n > m # minimum norm solution
@@ -859,9 +860,7 @@ function (\)(A::Union{QR{TA},QRCompactWY{TA},QRPivoted{TA}}, B::AbstractVecOrMat
 
     X = _zeros(S, B, n)
     X[1:size(B, 1), :] = B
-
     ldiv!(AA, X)
-
     return _cut_B(X, 1:n)
 end
 
diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl
index 2f3c32a3887b7..06ec00d34b312 100644
--- a/stdlib/LinearAlgebra/src/special.jl
+++ b/stdlib/LinearAlgebra/src/special.jl
@@ -118,8 +118,8 @@ for op in (:+, :-)
     end
 end
 
-mul!(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
-    (B = adjB.parent; mul!(full!(A), adjoint(B)))
+mul1!(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
+    (B = adjB.parent; mul1!(full!(A), adjoint(B)))
 *(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
     (B = adjB.parent; *(copyto!(similar(parent(A)), A), adjoint(B)))
 
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index 0f69a192d00c0..bb82bd18fcd34 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -457,33 +457,35 @@ fillstored!(A::UnitUpperTriangular, x) = (fillband!(A.data, x, 1, size(A,2)-1);
 # BlasFloat routines #
 ######################
 
-mul!(A::Tridiagonal, B::AbstractTriangular) = A*full!(B)
+mul2!(A::Tridiagonal, B::AbstractTriangular) = A*full!(B) # is this necessary?
+
 mul!(C::AbstractMatrix, A::AbstractTriangular, B::Tridiagonal) = mul!(C, copyto!(similar(parent(A)), A), B)
 mul!(C::AbstractMatrix, A::Tridiagonal, B::AbstractTriangular) = mul!(C, A, copyto!(similar(parent(B)), B))
 mul!(C::AbstractVector, A::AbstractTriangular, transB::Transpose{<:Any,<:AbstractVecOrMat}) =
-    (B = transB.parent; mul!(A, transpose!(C, B)))
+    (B = transB.parent; mul2!(A, transpose!(C, B)))
 mul!(C::AbstractMatrix, A::AbstractTriangular, transB::Transpose{<:Any,<:AbstractVecOrMat}) =
-    (B = transB.parent; mul!(A, transpose!(C, B)))
+    (B = transB.parent; mul2!(A, transpose!(C, B)))
 mul!(C::AbstractMatrix, A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractVecOrMat}) =
-    (B = adjB.parent; mul!(A, adjoint!(C, B)))
+    (B = adjB.parent; mul2!(A, adjoint!(C, B)))
 mul!(C::AbstractVecOrMat, A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractVecOrMat}) =
-    (B = adjB.parent; mul!(A, adjoint!(C, B)))
+    (B = adjB.parent; mul2!(A, adjoint!(C, B)))
+
 # The three methods for each op are neceesary to avoid ambiguities with definitions in matmul.jl
-mul!(C::AbstractVector  , A::AbstractTriangular, B::AbstractVector)   = mul!(A, copyto!(C, B))
-mul!(C::AbstractMatrix  , A::AbstractTriangular, B::AbstractVecOrMat) = mul!(A, copyto!(C, B))
-mul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) = mul!(A, copyto!(C, B))
+mul!(C::AbstractVector  , A::AbstractTriangular, B::AbstractVector)   = mul2!(A, copyto!(C, B))
+mul!(C::AbstractMatrix  , A::AbstractTriangular, B::AbstractVecOrMat) = mul2!(A, copyto!(C, B))
+mul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) = mul2!(A, copyto!(C, B))
 mul!(C::AbstractVector  , adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractVector) =
-    (A = adjA.parent; mul!(adjoint(A), copyto!(C, B)))
+    (A = adjA.parent; mul2!(adjoint(A), copyto!(C, B)))
 mul!(C::AbstractMatrix  , adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) =
-    (A = adjA.parent; mul!(adjoint(A), copyto!(C, B)))
+    (A = adjA.parent; mul2!(adjoint(A), copyto!(C, B)))
 mul!(C::AbstractVecOrMat, adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) =
-    (A = adjA.parent; mul!(adjoint(A), copyto!(C, B)))
+    (A = adjA.parent; mul2!(adjoint(A), copyto!(C, B)))
 mul!(C::AbstractVector  , transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractVector) =
-    (A = transA.parent; mul!(transpose(A), copyto!(C, B)))
+    (A = transA.parent; mul2!(transpose(A), copyto!(C, B)))
 mul!(C::AbstractMatrix  , transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) =
-    (A = transA.parent; mul!(transpose(A), copyto!(C, B)))
+    (A = transA.parent; mul2!(transpose(A), copyto!(C, B)))
 mul!(C::AbstractVecOrMat, transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractVecOrMat) =
-    (A = transA.parent; mul!(transpose(A), copyto!(C, B)))
+    (A = transA.parent; mul2!(transpose(A), copyto!(C, B)))
 mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B))
 mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B))
 mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = mul!(C, A, copy(B))
@@ -491,40 +493,39 @@ mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<
 mul!(C::AbstractVector, A::Adjoint{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B)))
 mul!(C::AbstractVector, A::Transpose{<:Any,<:AbstractTriangular}, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B)))
 
-
 for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'),
                             (:UnitLowerTriangular, 'L', 'U'),
                             (:UpperTriangular, 'U', 'N'),
                             (:UnitUpperTriangular, 'U', 'U'))
     @eval begin
         # Vector multiplication
-        mul!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} =
+        mul2!(A::$t{T,<:StridedMatrix}, b::StridedVector{T}) where {T<:BlasFloat} =
             BLAS.trmv!($uploc, 'N', $isunitc, A.data, b)
-        mul!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasFloat} =
+        mul2!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasFloat} =
             (A = transA.parent; BLAS.trmv!($uploc, 'T', $isunitc, A.data, b))
-        mul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasReal} =
+        mul2!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasReal} =
             (A = adjA.parent; BLAS.trmv!($uploc, 'T', $isunitc, A.data, b))
-        mul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasComplex} =
+        mul2!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, b::StridedVector{T}) where {T<:BlasComplex} =
             (A = adjA.parent; BLAS.trmv!($uploc, 'C', $isunitc, A.data, b))
 
         # Matrix multiplication
-        mul!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} =
+        mul2!(A::$t{T,<:StridedMatrix}, B::StridedMatrix{T}) where {T<:BlasFloat} =
             BLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), A.data, B)
-        mul!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} =
+        mul1!(A::StridedMatrix{T}, B::$t{T,<:StridedMatrix}) where {T<:BlasFloat} =
             BLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), B.data, A)
 
-        mul!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasFloat} =
+        mul2!(transA::Transpose{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasFloat} =
             (A = transA.parent; BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B))
-        mul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasComplex} =
+        mul2!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasComplex} =
             (A = adjA.parent; BLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), A.data, B))
-        mul!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasReal} =
+        mul2!(adjA::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}, B::StridedMatrix{T}) where {T<:BlasReal} =
             (A = adjA.parent; BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B))
 
-        mul!(A::StridedMatrix{T}, transB::Transpose{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasFloat} =
+        mul1!(A::StridedMatrix{T}, transB::Transpose{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasFloat} =
             (B = transB.parent; BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A))
-        mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasComplex} =
+        mul1!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasComplex} =
             (B = adjB.parent; BLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), B.data, A))
-        mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasReal} =
+        mul1!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:$t{T,<:StridedMatrix}}) where {T<:BlasReal} =
             (B = adjB.parent; BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A))
 
         # Left division
@@ -660,7 +661,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
 end
 
 ## Generic triangular multiplication
-function mul!(A::UpperTriangular, B::StridedVecOrMat)
+function mul2!(A::UpperTriangular, B::StridedVecOrMat)
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
         throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
@@ -676,7 +677,8 @@ function mul!(A::UpperTriangular, B::StridedVecOrMat)
     end
     B
 end
-function mul!(A::UnitUpperTriangular, B::StridedVecOrMat)
+
+function mul2!(A::UnitUpperTriangular, B::StridedVecOrMat)
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
         throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
@@ -693,7 +695,7 @@ function mul!(A::UnitUpperTriangular, B::StridedVecOrMat)
     B
 end
 
-function mul!(A::LowerTriangular, B::StridedVecOrMat)
+function mul2!(A::LowerTriangular, B::StridedVecOrMat)
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
         throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
@@ -709,7 +711,7 @@ function mul!(A::LowerTriangular, B::StridedVecOrMat)
     end
     B
 end
-function mul!(A::UnitLowerTriangular, B::StridedVecOrMat)
+function mul2!(A::UnitLowerTriangular, B::StridedVecOrMat)
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
         throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
@@ -726,7 +728,7 @@ function mul!(A::UnitLowerTriangular, B::StridedVecOrMat)
     B
 end
 
-function mul!(adjA::Adjoint{<:Any,<:UpperTriangular}, B::StridedVecOrMat)
+function mul2!(adjA::Adjoint{<:Any,<:UpperTriangular}, B::StridedVecOrMat)
     A = adjA.parent
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
@@ -743,7 +745,8 @@ function mul!(adjA::Adjoint{<:Any,<:UpperTriangular}, B::StridedVecOrMat)
     end
     B
 end
-function mul!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat)
+
+function mul2!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat)
     A = adjA.parent
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
@@ -761,7 +764,7 @@ function mul!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat)
     B
 end
 
-function mul!(adjA::Adjoint{<:Any,<:LowerTriangular}, B::StridedVecOrMat)
+function mul2!(adjA::Adjoint{<:Any,<:LowerTriangular}, B::StridedVecOrMat)
     A = adjA.parent
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
@@ -778,7 +781,7 @@ function mul!(adjA::Adjoint{<:Any,<:LowerTriangular}, B::StridedVecOrMat)
     end
     B
 end
-function mul!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat)
+function mul2!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat)
     A = adjA.parent
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
@@ -796,7 +799,7 @@ function mul!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat)
     B
 end
 
-function mul!(transA::Transpose{<:Any,<:UpperTriangular}, B::StridedVecOrMat)
+function mul2!(transA::Transpose{<:Any,<:UpperTriangular}, B::StridedVecOrMat)
     A = transA.parent
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
@@ -813,7 +816,7 @@ function mul!(transA::Transpose{<:Any,<:UpperTriangular}, B::StridedVecOrMat)
     end
     B
 end
-function mul!(transA::Transpose{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat)
+function mul2!(transA::Transpose{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat)
     A = transA.parent
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
@@ -831,7 +834,7 @@ function mul!(transA::Transpose{<:Any,<:UnitUpperTriangular}, B::StridedVecOrMat
     B
 end
 
-function mul!(transA::Transpose{<:Any,<:LowerTriangular}, B::StridedVecOrMat)
+function mul2!(transA::Transpose{<:Any,<:LowerTriangular}, B::StridedVecOrMat)
     A = transA.parent
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
@@ -848,7 +851,7 @@ function mul!(transA::Transpose{<:Any,<:LowerTriangular}, B::StridedVecOrMat)
     end
     B
 end
-function mul!(transA::Transpose{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat)
+function mul2!(transA::Transpose{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat)
     A = transA.parent
     m, n = size(B, 1), size(B, 2)
     if m != size(A, 1)
@@ -866,7 +869,7 @@ function mul!(transA::Transpose{<:Any,<:UnitLowerTriangular}, B::StridedVecOrMat
     B
 end
 
-function mul!(A::StridedMatrix, B::UpperTriangular)
+function mul1!(A::StridedMatrix, B::UpperTriangular)
     m, n = size(A)
     if size(B, 1) != n
         throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
@@ -882,7 +885,7 @@ function mul!(A::StridedMatrix, B::UpperTriangular)
     end
     A
 end
-function mul!(A::StridedMatrix, B::UnitUpperTriangular)
+function mul1!(A::StridedMatrix, B::UnitUpperTriangular)
     m, n = size(A)
     if size(B, 1) != n
         throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
@@ -899,7 +902,7 @@ function mul!(A::StridedMatrix, B::UnitUpperTriangular)
     A
 end
 
-function mul!(A::StridedMatrix, B::LowerTriangular)
+function mul1!(A::StridedMatrix, B::LowerTriangular)
     m, n = size(A)
     if size(B, 1) != n
         throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
@@ -915,7 +918,7 @@ function mul!(A::StridedMatrix, B::LowerTriangular)
     end
     A
 end
-function mul!(A::StridedMatrix, B::UnitLowerTriangular)
+function mul1!(A::StridedMatrix, B::UnitLowerTriangular)
     m, n = size(A)
     if size(B, 1) != n
         throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
@@ -932,7 +935,7 @@ function mul!(A::StridedMatrix, B::UnitLowerTriangular)
     A
 end
 
-function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UpperTriangular})
+function mul1!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UpperTriangular})
     B = adjB.parent
     m, n = size(A)
     if size(B, 1) != n
@@ -949,7 +952,8 @@ function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UpperTriangular})
     end
     A
 end
-function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitUpperTriangular})
+
+function mul1!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitUpperTriangular})
     B = adjB.parent
     m, n = size(A)
     if size(B, 1) != n
@@ -967,7 +971,7 @@ function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitUpperTriangular})
     A
 end
 
-function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:LowerTriangular})
+function mul1!(A::StridedMatrix, adjB::Adjoint{<:Any,<:LowerTriangular})
     B = adjB.parent
     m, n = size(A)
     if size(B, 1) != n
@@ -984,7 +988,8 @@ function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:LowerTriangular})
     end
     A
 end
-function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitLowerTriangular})
+
+function mul1!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitLowerTriangular})
     B = adjB.parent
     m, n = size(A)
     if size(B, 1) != n
@@ -1002,7 +1007,7 @@ function mul!(A::StridedMatrix, adjB::Adjoint{<:Any,<:UnitLowerTriangular})
     A
 end
 
-function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UpperTriangular})
+function mul1!(A::StridedMatrix, transB::Transpose{<:Any,<:UpperTriangular})
     B = transB.parent
     m, n = size(A)
     if size(B, 1) != n
@@ -1019,7 +1024,7 @@ function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UpperTriangular})
     end
     A
 end
-function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitUpperTriangular})
+function mul1!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitUpperTriangular})
     B = transB.parent
     m, n = size(A)
     if size(B, 1) != n
@@ -1037,7 +1042,7 @@ function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitUpperTriangular})
     A
 end
 
-function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:LowerTriangular})
+function mul1!(A::StridedMatrix, transB::Transpose{<:Any,<:LowerTriangular})
     B = transB.parent
     m, n = size(A)
     if size(B, 1) != n
@@ -1054,7 +1059,8 @@ function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:LowerTriangular})
     end
     A
 end
-function mul!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitLowerTriangular})
+
+function mul1!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitLowerTriangular})
     B = transB.parent
     m, n = size(A)
     if size(B, 1) != n
@@ -1136,7 +1142,7 @@ function naivesub!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector
 end
 # in the following transpose and conjugate transpose naive substitution variants,
 # accumulating in z rather than b[j] significantly improves performance as of Dec 2015
-function ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector = b)
+function ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector)
     A = transA.parent
     n = size(A, 1)
     if !(n == length(b) == length(x))
@@ -1152,7 +1158,9 @@ function ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector, x:
     end
     x
 end
-function ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector = b)
+ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector) = ldiv!(transA, b, b)
+
+function ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector)
     A = transA.parent
     n = size(A, 1)
     if !(n == length(b) == length(x))
@@ -1167,7 +1175,9 @@ function ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector
     end
     x
 end
-function ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector = b)
+ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector) = ldiv!(transA, b, b)
+
+function ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector)
     A = transA.parent
     n = size(A, 1)
     if !(n == length(b) == length(x))
@@ -1183,7 +1193,9 @@ function ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector, x:
     end
     x
 end
-function ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector = b)
+ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector) = ldiv!(transA, b, b)
+
+function ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector)
     A = transA.parent
     n = size(A, 1)
     if !(n == length(b) == length(x))
@@ -1198,7 +1210,9 @@ function ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector
     end
     x
 end
-function ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector = b)
+ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector) = ldiv!(transA, b, b)
+
+function ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector)
     A = adjA.parent
     n = size(A, 1)
     if !(n == length(b) == length(x))
@@ -1214,7 +1228,9 @@ function ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::Abs
     end
     x
 end
-function ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector = b)
+ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector) = ldiv!(adjA, b, b)
+
+function ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector)
     A = adjA.parent
     n = size(A, 1)
     if !(n == length(b) == length(x))
@@ -1229,7 +1245,9 @@ function ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x:
     end
     x
 end
-function ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector = b)
+ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector) = ldiv!(adjA, b, b)
+
+function ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector)
     A = adjA.parent
     n = size(A, 1)
     if !(n == length(b) == length(x))
@@ -1245,7 +1263,9 @@ function ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector, x::Abs
     end
     x
 end
-function ldiv!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector = b)
+ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector) = ldiv!(adjA, b, b)
+
+function ldiv!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector)
     A = adjA.parent
     n = size(A, 1)
     if !(n == length(b) == length(x))
@@ -1260,6 +1280,7 @@ function ldiv!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x:
     end
     x
 end
+ldiv!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, b::AbstractVector) = ldiv!(adjA, b, b)
 
 function rdiv!(A::StridedMatrix, B::UpperTriangular)
     m, n = size(A)
@@ -1467,14 +1488,14 @@ function rdiv!(A::StridedMatrix, transB::Transpose{<:Any,<:UnitLowerTriangular})
     A
 end
 
-mul!(adjA::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) =
-    (A = adjA.parent; UpperTriangular(mul!(adjoint(A), triu!(B.data))))
-mul!(adjA::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) =
-    (A = adjA.parent; LowerTriangular(mul!(adjoint(A), tril!(B.data))))
-mul!(transA::Transpose{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) =
-    (A = transA.parent; UpperTriangular(mul!(transpose(A), triu!(B.data))))
-mul!(transA::Transpose{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) =
-    (A = transA.parent; LowerTriangular(mul!(transpose(A), tril!(B.data))))
+mul2!(adjA::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) =
+    (A = adjA.parent; UpperTriangular(mul2!(adjoint(A), triu!(B.data))))
+mul2!(adjA::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) =
+    (A = adjA.parent; LowerTriangular(mul2!(adjoint(A), tril!(B.data))))
+mul2!(transA::Transpose{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) =
+    (A = transA.parent; UpperTriangular(mul2!(transpose(A), triu!(B.data))))
+mul2!(transA::Transpose{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) =
+    (A = transA.parent; LowerTriangular(mul2!(transpose(A), tril!(B.data))))
 ldiv!(adjA::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}, B::UpperTriangular) =
     (A = adjA.parent; UpperTriangular(ldiv!(adjoint(A), triu!(B.data))))
 ldiv!(adjA::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}, B::LowerTriangular) =
@@ -1489,14 +1510,14 @@ rdiv!(A::UpperTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) =
 rdiv!(A::LowerTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) =
     LowerTriangular(rdiv!(tril!(A.data), B))
 
-mul!(A::UpperTriangular, adjB::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) =
-    (B = adjB.parent; UpperTriangular(mul!(triu!(A.data), adjoint(B))))
-mul!(A::LowerTriangular, adjB::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) =
-    (B = adjB.parent; LowerTriangular(mul!(tril!(A.data), adjoint(B))))
-mul!(A::UpperTriangular, transB::Transpose{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) =
-    (B = transB.parent; UpperTriangular(mul!(triu!(A.data), transpose(B))))
-mul!(A::LowerTriangular, transB::Transpose{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) =
-    (B = transB.parent; LowerTriangular(mul!(tril!(A.data), transpose(B))))
+mul1!(A::UpperTriangular, adjB::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) =
+    (B = adjB.parent; UpperTriangular(mul1!(triu!(A.data), adjoint(B))))
+mul1!(A::LowerTriangular, adjB::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) =
+    (B = adjB.parent; LowerTriangular(mul1!(tril!(A.data), adjoint(B))))
+mul1!(A::UpperTriangular, transB::Transpose{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) =
+    (B = transB.parent; UpperTriangular(mul1!(triu!(A.data), transpose(B))))
+mul1!(A::LowerTriangular, transB::Transpose{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) =
+    (B = transB.parent; LowerTriangular(mul1!(tril!(A.data), transpose(B))))
 rdiv!(A::UpperTriangular, adjB::Adjoint{<:Any,<:Union{LowerTriangular,UnitLowerTriangular}}) =
     (B = adjB.parent; UpperTriangular(rdiv!(triu!(A.data), adjoint(B))))
 rdiv!(A::LowerTriangular, adjB::Adjoint{<:Any,<:Union{UpperTriangular,UnitUpperTriangular}}) =
@@ -1515,47 +1536,47 @@ rdiv!(A::LowerTriangular, transB::Transpose{<:Any,<:Union{UpperTriangular,UnitUp
 
 ## Some Triangular-Triangular cases. We might want to write taylored methods
 ## for these cases, but I'm not sure it is worth it.
-(*)(A::Union{Tridiagonal,SymTridiagonal}, B::AbstractTriangular) = mul!(Matrix(A), B)
+(*)(A::Union{Tridiagonal,SymTridiagonal}, B::AbstractTriangular) = mul1!(Matrix(A), B)
 
-for (f1, f2) in ((:*, :mul!), (:\, :ldiv!))
+for (f, f2!) in ((:*, :mul2!), (:\, :ldiv!))
     @eval begin
-        function ($f1)(A::LowerTriangular, B::LowerTriangular)
-            TAB = typeof(($f1)(zero(eltype(A)), zero(eltype(B))) +
-                         ($f1)(zero(eltype(A)), zero(eltype(B))))
+        function ($f)(A::LowerTriangular, B::LowerTriangular)
+            TAB = typeof(($f)(zero(eltype(A)), zero(eltype(B))) +
+                         ($f)(zero(eltype(A)), zero(eltype(B))))
             BB = similar(B, TAB, size(B))
             copyto!(BB, B)
-            return LowerTriangular($f2(convert(AbstractMatrix{TAB}, A), BB))
+            return LowerTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
         end
 
-        function $(f1)(A::UnitLowerTriangular, B::LowerTriangular)
+        function $(f)(A::UnitLowerTriangular, B::LowerTriangular)
             TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) +
                          (*)(zero(eltype(A)), zero(eltype(B))))
             BB = similar(B, TAB, size(B))
             copyto!(BB, B)
-            return LowerTriangular($f2(convert(AbstractMatrix{TAB}, A), BB))
+            return LowerTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
         end
 
-        function ($f1)(A::UpperTriangular, B::UpperTriangular)
-            TAB = typeof(($f1)(zero(eltype(A)), zero(eltype(B))) +
-                         ($f1)(zero(eltype(A)), zero(eltype(B))))
+        function ($f)(A::UpperTriangular, B::UpperTriangular)
+            TAB = typeof(($f)(zero(eltype(A)), zero(eltype(B))) +
+                         ($f)(zero(eltype(A)), zero(eltype(B))))
             BB = similar(B, TAB, size(B))
             copyto!(BB, B)
-            return UpperTriangular($f2(convert(AbstractMatrix{TAB}, A), BB))
+            return UpperTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
         end
 
-        function ($f1)(A::UnitUpperTriangular, B::UpperTriangular)
+        function ($f)(A::UnitUpperTriangular, B::UpperTriangular)
             TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) +
                          (*)(zero(eltype(A)), zero(eltype(B))))
             BB = similar(B, TAB, size(B))
             copyto!(BB, B)
-            return UpperTriangular($f2(convert(AbstractMatrix{TAB}, A), BB))
+            return UpperTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
         end
     end
 end
 
 for (ipop, op, xformtype, xformop) in (
-        (:mul!, :*, :Adjoint, :adjoint),
-        (:mul!, :*, :Transpose, :transpose),
+        (:mul2!, :*, :Adjoint, :adjoint),
+        (:mul2!, :*, :Transpose, :transpose),
         (:ldiv!, :\, :Adjoint, :adjoint),
         (:ldiv!, :\, :Transpose, :transpose))
     @eval begin
@@ -1627,8 +1648,8 @@ function (/)(A::UpperTriangular, B::UnitUpperTriangular)
 end
 
 for (ipop, op, xformtype, xformop) in (
-        (:mul!, :*, :Adjoint, :adjoint),
-        (:mul!, :*, :Transpose, :transpose),
+        (:mul1!, :*, :Adjoint, :adjoint),
+        (:mul1!, :*, :Transpose, :transpose),
         (:rdiv!, :/, :Adjoint, :adjoint),
         (:rdiv!, :/, :Transpose, :transpose))
     @eval begin
@@ -1671,26 +1692,25 @@ for (ipop, op, xformtype, xformop) in (
 end
 
 ## The general promotion methods
-
 function *(A::AbstractTriangular, B::AbstractTriangular)
     TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
     BB = similar(B, TAB, size(B))
     copyto!(BB, B)
-    mul!(convert(AbstractArray{TAB}, A), BB)
+    mul2!(convert(AbstractArray{TAB}, A), BB)
 end
 function *(adjA::Adjoint{<:Any,<:AbstractTriangular}, B::AbstractTriangular)
     A = adjA.parent
     TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
     BB = similar(B, TAB, size(B))
     copyto!(BB, B)
-    mul!(adjoint(convert(AbstractArray{TAB}, A)), BB)
+    mul2!(adjoint(convert(AbstractArray{TAB}, A)), BB)
 end
 function *(transA::Transpose{<:Any,<:AbstractTriangular}, B::AbstractTriangular)
     A = transA.parent
     TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
     BB = similar(B, TAB, size(B))
     copyto!(BB, B)
-    mul!(transpose(convert(AbstractArray{TAB}, A)), BB)
+    mul2!(transpose(convert(AbstractArray{TAB}, A)), BB)
 end
 
 function *(A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractTriangular})
@@ -1698,14 +1718,14 @@ function *(A::AbstractTriangular, adjB::Adjoint{<:Any,<:AbstractTriangular})
     TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
     AA = similar(A, TAB, size(A))
     copyto!(AA, A)
-    mul!(AA, adjoint(convert(AbstractArray{TAB}, B)))
+    mul1!(AA, adjoint(convert(AbstractArray{TAB}, B)))
 end
 function *(A::AbstractTriangular, transB::Transpose{<:Any,<:AbstractTriangular})
     B = transB.parent
     TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
     AA = similar(A, TAB, size(A))
     copyto!(AA, A)
-    mul!(AA, transpose(convert(AbstractArray{TAB}, B)))
+    mul1!(AA, transpose(convert(AbstractArray{TAB}, B)))
 end
 
 for mat in (:AbstractVector, :AbstractMatrix)
@@ -1715,21 +1735,21 @@ for mat in (:AbstractVector, :AbstractMatrix)
             TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
             BB = similar(B, TAB, size(B))
             copyto!(BB, B)
-            mul!(convert(AbstractArray{TAB}, A), BB)
+            mul2!(convert(AbstractArray{TAB}, A), BB)
         end
         function *(adjA::Adjoint{<:Any,<:AbstractTriangular}, B::$mat)
             A = adjA.parent
             TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
             BB = similar(B, TAB, size(B))
             copyto!(BB, B)
-            mul!(adjoint(convert(AbstractArray{TAB}, A)), BB)
+            mul2!(adjoint(convert(AbstractArray{TAB}, A)), BB)
         end
         function *(transA::Transpose{<:Any,<:AbstractTriangular}, B::$mat)
             A = transA.parent
             TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
             BB = similar(B, TAB, size(B))
             copyto!(BB, B)
-            mul!(transpose(convert(AbstractArray{TAB}, A)), BB)
+            mul2!(transpose(convert(AbstractArray{TAB}, A)), BB)
         end
     end
     ### Left division with triangle to the left hence rhs cannot be transposed. No quotients.
@@ -1831,21 +1851,21 @@ function *(A::AbstractMatrix, B::AbstractTriangular)
     TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
     AA = similar(A, TAB, size(A))
     copyto!(AA, A)
-    mul!(AA, convert(AbstractArray{TAB}, B))
+    mul1!(AA, convert(AbstractArray{TAB}, B))
 end
 function *(A::AbstractMatrix, adjB::Adjoint{<:Any,<:AbstractTriangular})
     B = adjB.parent
     TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
     AA = similar(A, TAB, size(A))
     copyto!(AA, A)
-    mul!(AA, adjoint(convert(AbstractArray{TAB}, B)))
+    mul1!(AA, adjoint(convert(AbstractArray{TAB}, B)))
 end
 function *(A::AbstractMatrix, transB::Transpose{<:Any,<:AbstractTriangular})
     B = transB.parent
     TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
     AA = similar(A, TAB, size(A))
     copyto!(AA, A)
-    mul!(AA, transpose(convert(AbstractArray{TAB}, B)))
+    mul1!(AA, transpose(convert(AbstractArray{TAB}, B)))
 end
 # ambiguity resolution with definitions in linalg/rowvector.jl
 *(v::AdjointAbsVec, A::AbstractTriangular) = adjoint(adjoint(A) * v.parent)
diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl
index 6c467d293f1f0..e0df205b0ca3f 100644
--- a/stdlib/LinearAlgebra/test/diagonal.jl
+++ b/stdlib/LinearAlgebra/test/diagonal.jl
@@ -3,7 +3,7 @@
 module TestDiagonal
 
 using Test, LinearAlgebra, SparseArrays, Random
-using LinearAlgebra: mul!, ldiv!, rdiv!, BlasFloat, BlasComplex, SingularException
+using LinearAlgebra: mul!, mul1!, mul2!, ldiv!, rdiv!, BlasFloat, BlasComplex, SingularException
 
 n=12 #Size of matrix problem to test
 srand(1)
@@ -147,9 +147,9 @@ srand(1)
             if relty <: BlasFloat
                 b = rand(elty,n,n)
                 b = sparse(b)
-                @test mul!(copy(D), copy(b)) ≈ Array(D)*Array(b)
-                @test mul!(transpose(copy(D)), copy(b)) ≈ transpose(Array(D))*Array(b)
-                @test mul!(adjoint(copy(D)), copy(b)) ≈ Array(D)'*Array(b)
+                @test mul2!(copy(D), copy(b)) ≈ Array(D)*Array(b)
+                @test mul2!(transpose(copy(D)), copy(b)) ≈ transpose(Array(D))*Array(b)
+                @test mul2!(adjoint(copy(D)), copy(b)) ≈ Array(D)'*Array(b)
             end
         end
 
@@ -178,13 +178,13 @@ srand(1)
         VV = Array(D)
         DD = copy(D)
         r  = VV * Matrix(D)
-        @test Array(mul!(VV, DD)) ≈ r ≈ Array(D)*Array(D)
+        @test Array(mul1!(VV, DD)) ≈ r ≈ Array(D)*Array(D)
         DD = copy(D)
         r  = VV * transpose(Array(D))
-        @test Array(mul!(VV, transpose(DD))) ≈ r
+        @test Array(mul1!(VV, transpose(DD))) ≈ r
         DD = copy(D)
         r  = VV * Array(D)'
-        @test Array(mul!(VV, adjoint(DD))) ≈ r
+        @test Array(mul1!(VV, adjoint(DD))) ≈ r
     end
     @testset "triu/tril" begin
         @test istriu(D)
@@ -348,9 +348,12 @@ end
 end
 
 let D1 = Diagonal(rand(5)), D2 = Diagonal(rand(5))
-    @test_throws MethodError mul!(D1,D2)
-    @test_throws MethodError mul!(transpose(D1),D2)
-    @test_throws MethodError mul!(adjoint(D1),D2)
+    @test LinAlg.mul1!(copy(D1),D2) == D1*D2
+    @test LinAlg.mul2!(D1,copy(D2)) == D1*D2
+    @test LinAlg.mul1!(copy(D1),transpose(D2)) == D1*transpose(D2)
+    @test LinAlg.mul2!(transpose(D1),copy(D2)) == transpose(D1)*D2
+    @test LinAlg.mul1!(copy(D1),adjoint(D2)) == D1*adjoint(D2)
+    @test LinAlg.mul2!(adjoint(D1),copy(D2)) == adjoint(D1)*D2
 end
 
 @testset "multiplication of QR Q-factor and Diagonal (#16615 spot test)" begin
@@ -358,7 +361,7 @@ end
     Q = qrfact(randn(5, 5)).Q
     @test D * Q' == Array(D) * Q'
     Q = qrfact(randn(5, 5), Val(true)).Q
-    @test_throws MethodError mul!(Q, D)
+    @test_throws ArgumentError mul2!(Q, D)
 end
 
 @testset "block diagonal matrices" begin
diff --git a/stdlib/LinearAlgebra/test/givens.jl b/stdlib/LinearAlgebra/test/givens.jl
index fd722d8aa4868..df15b0a088926 100644
--- a/stdlib/LinearAlgebra/test/givens.jl
+++ b/stdlib/LinearAlgebra/test/givens.jl
@@ -3,6 +3,7 @@
 module TestGivens
 
 using Test, LinearAlgebra, Random
+using LinearAlgebra: mul1!, mul2!
 
 # Test givens rotations
 @testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
@@ -17,11 +18,11 @@ using Test, LinearAlgebra, Random
         for j = 1:8
             for i = j+2:10
                 G, _ = givens(A, j+1, i, j)
-                mul!(G, A)
-                mul!(A, adjoint(G))
-                mul!(G, R)
+                mul2!(G, A)
+                mul1!(A, adjoint(G))
+                mul2!(G, R)
 
-                @test mul!(G,Matrix{elty}(I, 10, 10)) == [G[i,j] for i=1:10,j=1:10]
+                @test mul2!(G,Matrix{elty}(I, 10, 10)) == [G[i,j] for i=1:10,j=1:10]
 
                 @testset "transposes" begin
                     @test copy(G')*G*Matrix(elty(1)I, 10, 10) ≈ Matrix(I, 10, 10)
@@ -34,8 +35,8 @@ using Test, LinearAlgebra, Random
         @test_throws ArgumentError givens(A, 3, 3, 2)
         @test_throws ArgumentError givens(one(elty),zero(elty),2,2)
         G, _ = givens(one(elty),zero(elty),11,12)
-        @test_throws DimensionMismatch mul!(G, A)
-        @test_throws DimensionMismatch mul!(A, adjoint(G))
+        @test_throws DimensionMismatch mul2!(G, A)
+        @test_throws DimensionMismatch mul1!(A, adjoint(G))
         @test abs.(A) ≈ abs.(hessfact(Ac).H)
         @test norm(R*Matrix{elty}(I, 10, 10)) ≈ one(elty)
 
diff --git a/stdlib/LinearAlgebra/test/lq.jl b/stdlib/LinearAlgebra/test/lq.jl
index dc63e77a0441c..8b51b496f09f9 100644
--- a/stdlib/LinearAlgebra/test/lq.jl
+++ b/stdlib/LinearAlgebra/test/lq.jl
@@ -3,7 +3,7 @@
 module TestLQ
 
 using Test, LinearAlgebra, Random
-using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, mul!
+using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, mul1!, mul2!
 
 n = 10
 
@@ -21,7 +21,7 @@ breal = randn(n,2)/2
 bimg  = randn(n,2)/2
 
 # helper functions to unambiguously recover explicit forms of an LQPackedQ
-squareQ(Q::LinearAlgebra.LQPackedQ) = (n = size(Q.factors, 2); mul!(Q, Matrix{eltype(Q)}(I, n, n)))
+squareQ(Q::LinearAlgebra.LQPackedQ) = (n = size(Q.factors, 2); mul2!(Q, Matrix{eltype(Q)}(I, n, n)))
 rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q)
 
 @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64)
@@ -54,8 +54,6 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q)
                     @test size(lqa.Q,3) == 1
                     @test_throws ErrorException lqa.Z
                     @test Array(copy(adjoint(lqa))) ≈ a'
-                    @test lqa * lqa' ≈ a * a'
-                    @test lqa' * lqa ≈ a' * a
                     @test q*squareQ(q)' ≈ Matrix(I, n, n)
                     @test l*q ≈ a
                     @test Array(lqa) ≈ a
@@ -99,9 +97,9 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q)
             l,q = lqa.L, lqa.Q
             @test rectangularQ(q)*rectangularQ(q)' ≈ Matrix(I, n1, n1)
             @test squareQ(q)'*squareQ(q) ≈ Matrix(I, n1, n1)
-            @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1),q)
-            @test mul!(adjoint(q), rectangularQ(q)) ≈ Matrix(I, n1, n1)
-            @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1), adjoint(q))
+            @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1),q)
+            @test mul2!(adjoint(q), rectangularQ(q)) ≈ Matrix(I, n1, n1)
+            @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1), adjoint(q))
             @test_throws BoundsError size(q,-1)
         end
     end
@@ -149,7 +147,7 @@ end
     function getqs(F::LinearAlgebra.LQ)
         implicitQ = F.Q
         sq = size(implicitQ.factors, 2)
-        explicitQ = mul!(implicitQ, Matrix{eltype(implicitQ)}(I, sq, sq))
+        explicitQ = mul2!(implicitQ, Matrix{eltype(implicitQ)}(I, sq, sq))
         return implicitQ, explicitQ
     end
 
@@ -190,7 +188,7 @@ end
 @testset "postmultiplication with / right-application of LQPackedQ (#23779)" begin
     function getqs(F::LinearAlgebra.LQ)
         implicitQ = F.Q
-        explicitQ = mul!(implicitQ, Matrix{eltype(implicitQ)}(I, size(implicitQ)...))
+        explicitQ = mul2!(implicitQ, Matrix{eltype(implicitQ)}(I, size(implicitQ)...))
         return implicitQ, explicitQ
     end
     # for any shape m-by-n of LQ-factored matrix, where Q is an LQPackedQ
diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl
index b3ace0f81089c..99385030b556d 100644
--- a/stdlib/LinearAlgebra/test/qr.jl
+++ b/stdlib/LinearAlgebra/test/qr.jl
@@ -3,7 +3,7 @@
 module TestQR
 
 using Test, LinearAlgebra, Random
-using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, mul!
+using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, mul1!, mul2!
 
 n = 10
 
@@ -21,7 +21,7 @@ breal = randn(n,2)/2
 bimg  = randn(n,2)/2
 
 # helper functions to unambiguously recover explicit forms of an implicit QR Q
-squareQ(Q::LinearAlgebra.AbstractQ) = (sq = size(Q.factors, 1); mul!(Q, Matrix{eltype(Q)}(I, sq, sq)))
+squareQ(Q::LinearAlgebra.AbstractQ) = (sq = size(Q.factors, 1); mul2!(Q, Matrix{eltype(Q)}(I, sq, sq)))
 rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q)
 
 @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
@@ -136,20 +136,20 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q)
                 a = raw_a
                 qrpa = factorize(a[:,1:n1])
                 q, r = qrpa.Q, qrpa.R
-                @test mul!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n)
-                @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1),q)
-                @test mul!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n)
-                @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1), adjoint(q))
+                @test mul1!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n)
+                @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1),q)
+                @test mul1!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n)
+                @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1), adjoint(q))
                 @test_throws BoundsError size(q,-1)
-                @test_throws DimensionMismatch LinearAlgebra.mul!(q,zeros(eltya,n1+1))
-                @test_throws DimensionMismatch LinearAlgebra.mul!(adjoint(q), zeros(eltya,n1+1))
+                @test_throws DimensionMismatch LinearAlgebra.mul2!(q,zeros(eltya,n1+1))
+                @test_throws DimensionMismatch LinearAlgebra.mul2!(adjoint(q), zeros(eltya,n1+1))
 
                 qra = qrfact(a[:,1:n1], Val(false))
                 q, r = qra.Q, qra.R
-                @test mul!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n)
-                @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1),q)
-                @test mul!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n)
-                @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n+1),adjoint(q))
+                @test mul1!(copy(squareQ(q)'), q) ≈ Matrix(I, n, n)
+                @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1),q)
+                @test mul1!(squareQ(q), adjoint(q)) ≈ Matrix(I, n, n)
+                @test_throws DimensionMismatch mul1!(Matrix{eltya}(I, n+1, n+1),adjoint(q))
                 @test_throws BoundsError size(q,-1)
                 @test_throws DimensionMismatch q * Matrix{Int8}(I, n+4, n+4)
             end
diff --git a/stdlib/LinearAlgebra/test/special.jl b/stdlib/LinearAlgebra/test/special.jl
index 8799c298f5e3e..9563ce13507ca 100644
--- a/stdlib/LinearAlgebra/test/special.jl
+++ b/stdlib/LinearAlgebra/test/special.jl
@@ -3,7 +3,7 @@
 module TestSpecial
 
 using Test, LinearAlgebra, SparseArrays, Random
-using LinearAlgebra: mul!
+using LinearAlgebra: mul1!
 
 n= 10 #Size of matrix to test
 srand(1)
@@ -118,10 +118,10 @@ end
         b = rand(n,n)
         qrb = qrfact(b,Val(true))
         @test *(atri, adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q'
-        @test mul!(copy(atri), adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q'
+        @test mul1!(copy(atri), adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q'
         qrb = qrfact(b,Val(false))
         @test *(atri, adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q'
-        @test mul!(copy(atri), adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q'
+        @test mul1!(copy(atri), adjoint(qrb.Q)) ≈ Matrix(atri) * qrb.Q'
     end
 end
 
diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl
index 08469d27fd589..c2e6420af3734 100644
--- a/stdlib/LinearAlgebra/test/triangular.jl
+++ b/stdlib/LinearAlgebra/test/triangular.jl
@@ -6,7 +6,7 @@ debug = false
 using Test, LinearAlgebra, SparseArrays, Random
 using LinearAlgebra: BlasFloat, errorbounds, full!, naivesub!, transpose!,
     UnitUpperTriangular, UnitLowerTriangular,
-    mul!, rdiv!
+    mul!, rdiv!, mul1!, mul2!
 
 debug && println("Triangular matrices")
 
@@ -322,7 +322,7 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo
 
             if !(eltyB in (BigFloat, Complex{BigFloat})) # rand does not support BigFloat and Complex{BigFloat} as of Dec 2015
                 Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1))
-                @test mul!(Tri,copy(A1)) ≈ Tri*Matrix(A1)
+                @test mul2!(Tri,copy(A1)) ≈ Tri*Matrix(A1)
             end
 
             # Triangular-dense Matrix/vector multiplication
@@ -360,12 +360,12 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo
             end
             #error handling
             Ann, Bmm, bm = A1, Matrix{eltyB}(uninitialized, n+1, n+1), Vector{eltyB}(uninitialized, n+1)
-            @test_throws DimensionMismatch mul!(Ann, bm)
-            @test_throws DimensionMismatch mul!(Bmm, Ann)
-            @test_throws DimensionMismatch mul!(transpose(Ann), bm)
-            @test_throws DimensionMismatch mul!(adjoint(Ann), bm)
-            @test_throws DimensionMismatch mul!(Bmm, adjoint(Ann))
-            @test_throws DimensionMismatch mul!(Bmm, transpose(Ann))
+            @test_throws DimensionMismatch mul2!(Ann, bm)
+            @test_throws DimensionMismatch mul1!(Bmm, Ann)
+            @test_throws DimensionMismatch mul2!(transpose(Ann), bm)
+            @test_throws DimensionMismatch mul2!(adjoint(Ann), bm)
+            @test_throws DimensionMismatch mul1!(Bmm, adjoint(Ann))
+            @test_throws DimensionMismatch mul1!(Bmm, transpose(Ann))
 
             # ... and division
             @test A1\B[:,1] ≈ Matrix(A1)\B[:,1]
diff --git a/stdlib/SparseArrays/src/deprecated.jl b/stdlib/SparseArrays/src/deprecated.jl
index 829d6c9401a3d..a6691f3ffd80e 100644
--- a/stdlib/SparseArrays/src/deprecated.jl
+++ b/stdlib/SparseArrays/src/deprecated.jl
@@ -222,6 +222,11 @@ end
 import Base: asyncmap
 @deprecate asyncmap(f, s::AbstractSparseArray...; kwargs...) sparse(asyncmap(f, map(Array, s)...; kwargs...))
 
+# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from linalg.jl, to deprecate
+import Base: A_mul_B!, Ac_mul_B, Ac_mul_B!, At_mul_B, At_mul_B!
+import Base: A_mul_Bc, A_mul_Bt, Ac_mul_Bc, At_mul_Bt
+import Base: At_ldiv_B, Ac_ldiv_B, A_ldiv_B!
+import Base.LinAlg: At_ldiv_B!, Ac_ldiv_B!, A_rdiv_B!, A_rdiv_Bc!
 
 # END 0.7 deprecations
 
diff --git a/stdlib/SuiteSparse/src/deprecated.jl b/stdlib/SuiteSparse/src/deprecated.jl
index 2c36db610d028..14bdd94e7c2ac 100644
--- a/stdlib/SuiteSparse/src/deprecated.jl
+++ b/stdlib/SuiteSparse/src/deprecated.jl
@@ -30,12 +30,12 @@ end
         LinearAlgebra.ldiv!(X, transpose(lu), B)
     LinearAlgebra.Ac_ldiv_B!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
         LinearAlgebra.ldiv!(X, adjoint(lu), B)
-    LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(lu, B)
-    LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(transpose(lu), B)
-    LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(adjoint(lu), B)
-    LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(lu, B)
-    LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(transpose(lu), B)
-    LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(adjoint(lu), B)
+    LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(B, lu, copy(B))
+    LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(B, transpose(lu), copy(B))
+    LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = LinearAlgebra.ldiv!(B, adjoint(lu), copy(B))
+    LinearAlgebra.A_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(B, lu, copy(B))
+    LinearAlgebra.At_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(B, transpose(lu), copy(B))
+    LinearAlgebra.Ac_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = LinearAlgebra.ldiv!(B, adjoint(lu), copy(B))
 end
 
 # A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from src/spqr.jl, to deprecate
diff --git a/stdlib/SuiteSparse/src/spqr.jl b/stdlib/SuiteSparse/src/spqr.jl
index 742342d7efa34..4245c264d8da6 100644
--- a/stdlib/SuiteSparse/src/spqr.jl
+++ b/stdlib/SuiteSparse/src/spqr.jl
@@ -198,7 +198,7 @@ LinearAlgebra.qrfact(A::SparseMatrixCSC; tol = _default_tol(A)) = qrfact(A, Val{
 
 LinearAlgebra.qr(A::SparseMatrixCSC; tol = _default_tol(A)) = qr(A, Val{true}, tol = tol)
 
-function LinearAlgebra.mul!(Q::QRSparseQ, A::StridedVecOrMat)
+function LinearAlgebra.mul2!(Q::QRSparseQ, A::StridedVecOrMat)
     if size(A, 1) != size(Q, 1)
         throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))"))
     end
@@ -213,7 +213,7 @@ function LinearAlgebra.mul!(Q::QRSparseQ, A::StridedVecOrMat)
     return A
 end
 
-function LinearAlgebra.mul!(A::StridedMatrix, Q::QRSparseQ)
+function LinearAlgebra.mul1!(A::StridedMatrix, Q::QRSparseQ)
     if size(A, 2) != size(Q, 1)
         throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))"))
     end
@@ -227,7 +227,7 @@ function LinearAlgebra.mul!(A::StridedMatrix, Q::QRSparseQ)
     return A
 end
 
-function LinearAlgebra.mul!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMat)
+function LinearAlgebra.mul2!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMat)
     Q = adjQ.parent
     if size(A, 1) != size(Q, 1)
         throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))"))
@@ -243,7 +243,7 @@ function LinearAlgebra.mul!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMat
     return A
 end
 
-function LinearAlgebra.mul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRSparseQ})
+function LinearAlgebra.mul1!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRSparseQ})
     Q = adjQ.parent
     if size(A, 2) != size(Q, 1)
         throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))"))
@@ -375,7 +375,7 @@ function _ldiv_basic(F::QRSparse, B::StridedVecOrMat)
     X0 = view(X, 1:size(B, 1), :)
 
     # Apply Q' to B
-    LinearAlgebra.mul!(adjoint(F.Q), X0)
+    LinearAlgebra.mul2!(adjoint(F.Q), X0)
 
     # Zero out to get basic solution
     X[rnk + 1:end, :] = 0
diff --git a/stdlib/SuiteSparse/src/umfpack.jl b/stdlib/SuiteSparse/src/umfpack.jl
index 9cc32f2261e23..74ff6a6cfbb0f 100644
--- a/stdlib/SuiteSparse/src/umfpack.jl
+++ b/stdlib/SuiteSparse/src/umfpack.jl
@@ -388,6 +388,9 @@ function nnz(lu::UmfpackLU)
 end
 
 ### Solve with Factorization
+
+import Base.LinAlg.ldiv!
+
 ldiv!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
     ldiv!(B, lu, copy(B))
 ldiv!(translu::Transpose{T,<:UmfpackLU{T}}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
diff --git a/stdlib/SuiteSparse/test/spqr.jl b/stdlib/SuiteSparse/test/spqr.jl
index e59e1ff241026..c36545e173cd2 100644
--- a/stdlib/SuiteSparse/test/spqr.jl
+++ b/stdlib/SuiteSparse/test/spqr.jl
@@ -2,7 +2,7 @@
 
 using SuiteSparse.SPQR
 using SuiteSparse.CHOLMOD
-using LinearAlgebra: mul!, Adjoint, Transpose
+using LinearAlgebra: mul1!, mul2!, Adjoint, Transpose
 
 @testset "Sparse QR" begin
 m, n = 100, 10
@@ -43,10 +43,10 @@ nn = 100
         @test norm(R0[n + 1:end, :], 1) < 1e-12
 
         offsizeA = Matrix{Float64}(I, m+1, m+1)
-        @test_throws DimensionMismatch mul!(Q, offsizeA)
-        @test_throws DimensionMismatch mul!(adjoint(Q), offsizeA)
-        @test_throws DimensionMismatch mul!(offsizeA, Q)
-        @test_throws DimensionMismatch mul!(offsizeA, adjoint(Q))
+        @test_throws DimensionMismatch mul2!(Q, offsizeA)
+        @test_throws DimensionMismatch mul2!(adjoint(Q), offsizeA)
+        @test_throws DimensionMismatch mul1!(offsizeA, Q)
+        @test_throws DimensionMismatch mul1!(offsizeA, adjoint(Q))
     end
 
     @testset "element type of B: $eltyB" for eltyB in (Int, Float64, Complex{Float64})
diff --git a/stdlib/SuiteSparse/test/umfpack.jl b/stdlib/SuiteSparse/test/umfpack.jl
index 1e48075937895..9b4f965f6c361 100644
--- a/stdlib/SuiteSparse/test/umfpack.jl
+++ b/stdlib/SuiteSparse/test/umfpack.jl
@@ -32,7 +32,7 @@
 
             @test A*x ≈ b
             z = complex.(b)
-            x = SuiteSparse.ldiv!(lua, z)
+            x = Base.LinAlg.ldiv!(lua, z)
             @test x ≈ float([1:5;])
             @test z === x
             y = similar(z)
@@ -47,11 +47,11 @@
 
             @test A'*x ≈ b
             z = complex.(b)
-            x = SuiteSparse.ldiv!(adjoint(lua), z)
+            x = Base.LinAlg.ldiv!(adjoint(lua), z)
             @test x ≈ float([1:5;])
             @test x === z
             y = similar(x)
-            SuiteSparse.ldiv!(y, adjoint(lua), complex.(b))
+            Base.LinAlg.ldiv!(y, adjoint(lua), complex.(b))
             @test y ≈ x
 
             @test A'*x ≈ b
@@ -59,10 +59,10 @@
             @test x ≈ float([1:5;])
 
             @test transpose(A) * x ≈ b
-            x = SuiteSparse.ldiv!(transpose(lua), complex.(b))
+            x = Base.LinAlg.ldiv!(transpose(lua), complex.(b))
             @test x ≈ float([1:5;])
             y = similar(x)
-            SuiteSparse.ldiv!(y, transpose(lua), complex.(b))
+            Base.LinAlg.ldiv!(y, transpose(lua), complex.(b))
             @test y ≈ x
 
             @test transpose(A) * x ≈ b