diff --git a/NEWS.md b/NEWS.md
index 3db283a9a3b2c..2e90e85584bd4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -991,6 +991,8 @@ Deprecated or removed
 
   * `Base.@gc_preserve` has been deprecated in favor of `GC.@preserve` ([#25616]).
 
+  * `scale!` has been deprecated in favor of `mul!`, `mul1!`, and `mul2!` ([#25701]).
+
   * `DateTime()`, `Date()`, and `Time()` have been deprecated, instead use `DateTime(1)`, `Date(1)`
     and `Time(0)` respectively ([#23724]).
 
diff --git a/stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl b/stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl
index a9138bf66fb3e..83efd53258309 100644
--- a/stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl
+++ b/stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl
@@ -7,9 +7,9 @@ Arnoldi and Lanczos iteration for computing eigenvalues
 """
 module IterativeEigensolvers
 
-using LinearAlgebra: BlasFloat, BlasInt, SVD, checksquare, mul!,
-              UniformScaling, issymmetric, ishermitian,
-              factorize, I, scale!, qr
+using LinearAlgebra: BlasFloat, BlasInt, Diagonal, I, SVD, UniformScaling,
+                     checksquare, factorize,ishermitian, issymmetric, mul!,
+                     mul1!, qr
 import LinearAlgebra
 
 export eigs, svds
@@ -317,10 +317,10 @@ function _svds(X; nsv::Int = 6, ritzvec::Bool = true, tol::Float64 = 0.0, maxite
         # left_sv  = sqrt(2) * ex[2][ 1:size(X,1),     ind ] .* sign.(ex[1][ind]')
         if size(X, 1) >= size(X, 2)
             V = ex[2]
-            U = qr(scale!(X*V, inv.(svals)))[1]
+            U = qr(mul1!(X*V, Diagonal(inv.(svals))))[1]
         else
             U = ex[2]
-            V = qr(scale!(X'U, inv.(svals)))[1]
+            V = qr(mul1!(X'U, Diagonal(inv.(svals))))[1]
         end
 
         # right_sv = sqrt(2) * ex[2][ size(X,1)+1:end, ind ]
diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md
index 2bdd8028c9579..4ab2278fd1a06 100644
--- a/stdlib/LinearAlgebra/docs/src/index.md
+++ b/stdlib/LinearAlgebra/docs/src/index.md
@@ -354,7 +354,6 @@ LinearAlgebra.tril!
 LinearAlgebra.diagind
 LinearAlgebra.diag
 LinearAlgebra.diagm
-LinearAlgebra.scale!
 LinearAlgebra.rank
 LinearAlgebra.norm
 LinearAlgebra.vecnorm
@@ -430,6 +429,8 @@ below (e.g. `mul!`) according to the usual Julia convention.
 
 ```@docs
 LinearAlgebra.mul!
+LinearAlgebra.mul1!
+LinearAlgebra.mul2!
 LinearAlgebra.ldiv!
 LinearAlgebra.rdiv!
 ```
diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index 4d18d7d5d00fd..94697f6a00b4f 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -98,6 +98,7 @@ export
     istril,
     istriu,
     kron,
+    ldiv!,
     ldltfact!,
     ldltfact,
     linreg,
@@ -107,6 +108,9 @@ export
     lufact,
     lufact!,
     lyap,
+    mul!,
+    mul1!,
+    mul2!,
     norm,
     normalize,
     normalize!,
@@ -122,7 +126,7 @@ export
     lqfact!,
     lqfact,
     rank,
-    scale!,
+    rdiv!,
     schur,
     schurfact!,
     schurfact,
diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl
index 8748dd308a3e5..d838e347fadf2 100644
--- a/stdlib/LinearAlgebra/src/dense.jl
+++ b/stdlib/LinearAlgebra/src/dense.jl
@@ -14,21 +14,21 @@ const NRM2_CUTOFF = 32
 # This constant should ideally be determined by the actual CPU cache size
 const ISONE_CUTOFF = 2^21 # 2M
 
-function scale!(X::Array{T}, s::T) where T<:BlasFloat
+function mul1!(X::Array{T}, s::T) where T<:BlasFloat
     s == 0 && return fill!(X, zero(T))
     s == 1 && return X
     if length(X) < SCAL_CUTOFF
-        generic_scale!(X, s)
+        generic_mul1!(X, s)
     else
         BLAS.scal!(length(X), s, X, 1)
     end
     X
 end
 
-scale!(s::T, X::Array{T}) where {T<:BlasFloat} = scale!(X, s)
+mul2!(s::T, X::Array{T}) where {T<:BlasFloat} = mul1!(X, s)
 
-scale!(X::Array{T}, s::Number) where {T<:BlasFloat} = scale!(X, convert(T, s))
-function scale!(X::Array{T}, s::Real) where T<:BlasComplex
+mul1!(X::Array{T}, s::Number) where {T<:BlasFloat} = mul1!(X, convert(T, s))
+function mul1!(X::Array{T}, s::Real) where T<:BlasComplex
     R = typeof(real(zero(T)))
     GC.@preserve X BLAS.scal!(2*length(X), convert(R,s), convert(Ptr{R},pointer(X)), 1)
     X
@@ -1402,7 +1402,7 @@ function sylvester(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedMatrix{T})
 
     D = -(adjoint(QA) * (C*QB))
     Y, scale = LAPACK.trsyl!('N','N', RA, RB, D)
-    scale!(QA*(Y * adjoint(QB)), inv(scale))
+    mul1!(QA*(Y * adjoint(QB)), inv(scale))
 end
 sylvester(A::StridedMatrix{T}, B::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:Integer} = sylvester(float(A), float(B), float(C))
 
@@ -1445,7 +1445,7 @@ function lyap(A::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:BlasFloat}
 
     D = -(adjoint(Q) * (C*Q))
     Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, D)
-    scale!(Q*(Y * adjoint(Q)), inv(scale))
+    mul1!(Q*(Y * adjoint(Q)), inv(scale))
 end
 lyap(A::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:Integer} = lyap(float(A), float(C))
 lyap(a::T, c::T) where {T<:Number} = -c/(2a)
diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl
index 3859eb5bdfa86..34092c9ef3a9f 100644
--- a/stdlib/LinearAlgebra/src/deprecated.jl
+++ b/stdlib/LinearAlgebra/src/deprecated.jl
@@ -1268,3 +1268,11 @@ function getindex(F::Factorization, s::Symbol)
     return getproperty(F, s)
 end
 @deprecate getq(F::Factorization) F.Q
+
+# Deprecate scaling
+@deprecate scale!(A::AbstractArray, b::Number)                             mul1!(A, b)
+@deprecate scale!(a::Number, B::AbstractArray)                             mul2!(a, B)
+@deprecate scale!(A::AbstractMatrix, b::AbstractVector)                    mul1!(A, Diagonal(b))
+@deprecate scale!(a::AbstractVector, B::AbstractMatrix)                    mul2!(Diagonal(a), B)
+@deprecate scale!(C::AbstractMatrix, A::AbstractMatrix, b::AbstractVector) mul!(C, A, Diagonal(b))
+@deprecate scale!(C::AbstractMatrix, a::AbstractVector, B::AbstractMatrix) mul!(C, Diagonal(a), B)
diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl
index e90820b949eed..beb4921780eb6 100644
--- a/stdlib/LinearAlgebra/src/diagonal.jl
+++ b/stdlib/LinearAlgebra/src/diagonal.jl
@@ -155,10 +155,19 @@ end
 (*)(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)
+    mul!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), A, D)
 (*)(D::Diagonal, A::AbstractMatrix) =
-    scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), D.diag, A)
+    mul!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), D, A)
 
+function mul1!(A::AbstractMatrix, D::Diagonal)
+    A .= A .* transpose(D.diag)
+    return A
+end
+
+function mul2!(D::Diagonal, B::AbstractMatrix)
+    B .= D.diag .* B
+    return B
+end
 
 mul1!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) = typeof(A)(mul1!(A.data, D))
 function mul1!(A::UnitLowerTriangular, D::Diagonal)
@@ -233,15 +242,26 @@ end
 *(D::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) =
     Diagonal(transpose.(D.parent.diag) .* transpose.(B.parent.diag))
 
-mul1!(A::Diagonal,B::Diagonal) = Diagonal(A.diag .*= B.diag)
-mul2!(A::Diagonal,B::Diagonal) = Diagonal(B.diag .= A.diag .* 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))
+function mul2!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix)
+    A = adjA.parent
+    return mul2!(conj(A.diag), B)
+end
+function mul2!(transA::Transpose{<:Any,<:Diagonal}, B::AbstractMatrix)
+    A = transA.parent
+    return mul2!(A.diag, B)
+end
+
+function mul1!(A::AbstractMatrix, adjB::Adjoint{<:Any,<:Diagonal})
+    B = adjB.parent
+    return mul1!(A, conj(B.diag))
+end
+function mul1!(A::AbstractMatrix, transB::Transpose{<:Any,<:Diagonal})
+    B = transB.parent
+    return mul1!(A, B.diag)
+end
 
 # Get ambiguous method if try to unify AbstractVector/AbstractMatrix here using AbstractVecOrMat
 mul!(out::AbstractVector, A::Diagonal, in::AbstractVector) = out .= A.diag .* in
diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index f1f7f7b63dd5d..792fbbc39e044 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -4,21 +4,21 @@
 
 # For better performance when input and output are the same array
 # See https://github.com/JuliaLang/julia/issues/8415#issuecomment-56608729
-function generic_scale!(X::AbstractArray, s::Number)
+function generic_mul1!(X::AbstractArray, s::Number)
     @simd for I in eachindex(X)
         @inbounds X[I] *= s
     end
     X
 end
 
-function generic_scale!(s::Number, X::AbstractArray)
+function generic_mul2!(s::Number, X::AbstractArray)
     @simd for I in eachindex(X)
         @inbounds X[I] = s*X[I]
     end
     X
 end
 
-function generic_scale!(C::AbstractArray, X::AbstractArray, s::Number)
+function generic_mul!(C::AbstractArray, X::AbstractArray, s::Number)
     if _length(C) != _length(X)
         throw(DimensionMismatch("first array has length $(_length(C)) which does not match the length of the second, $(_length(X))."))
     end
@@ -28,7 +28,7 @@ function generic_scale!(C::AbstractArray, X::AbstractArray, s::Number)
     C
 end
 
-function generic_scale!(C::AbstractArray, s::Number, X::AbstractArray)
+function generic_mul!(C::AbstractArray, s::Number, X::AbstractArray)
     if _length(C) != _length(X)
         throw(DimensionMismatch("first array has length $(_length(C)) which does not
 match the length of the second, $(_length(X))."))
@@ -39,50 +39,48 @@ match the length of the second, $(_length(X))."))
     C
 end
 
-scale!(C::AbstractArray, s::Number, X::AbstractArray) = generic_scale!(C, X, s)
-scale!(C::AbstractArray, X::AbstractArray, s::Number) = generic_scale!(C, s, X)
+mul!(C::AbstractArray, s::Number, X::AbstractArray) = generic_mul!(C, X, s)
+mul!(C::AbstractArray, X::AbstractArray, s::Number) = generic_mul!(C, s, X)
 
 """
-    scale!(A, b)
-    scale!(b, A)
+    mul1!(A::AbstractArray, b::Number)
 
 Scale an array `A` by a scalar `b` overwriting `A` in-place.
 
-If `A` is a matrix and `b` is a vector, then `scale!(A,b)` scales each column `i` of `A` by
-`b[i]` (similar to `A*Diagonal(b)`), while `scale!(b,A)` scales each row `i` of `A` by `b[i]`
-(similar to `Diagonal(b)*A`), again operating in-place on `A`. An `InexactError` exception is
-thrown if the scaling produces a number not representable by the element type of `A`,
-e.g. for integer types.
-
 # Examples
 ```jldoctest
-julia> a = [1 2; 3 4]
+julia> A = [1 2; 3 4]
 2×2 Array{Int64,2}:
  1  2
  3  4
 
-julia> b = [1; 2]
-2-element Array{Int64,1}:
- 1
- 2
-
-julia> scale!(a,b)
+julia> mul1!(A, 2)
 2×2 Array{Int64,2}:
- 1  4
- 3  8
+ 2  4
+ 6  8
+```
+"""
+mul1!(A::AbstractArray, b::Number) = generic_mul1!(A, b)
 
-julia> a = [1 2; 3 4];
+"""
+    mul2!(a::Number, B::AbstractArray)
 
-julia> b = [1; 2];
+Scale an array `B` by a scalar `a` overwriting `B` in-place.
 
-julia> scale!(b,a)
+# Examples
+```jldoctest
+julia> B = [1 2; 3 4]
 2×2 Array{Int64,2}:
  1  2
+ 3  4
+
+julia> mul2!(2, B)
+2×2 Array{Int64,2}:
+ 2  4
  6  8
 ```
 """
-scale!(X::AbstractArray, s::Number) = generic_scale!(X, s)
-scale!(s::Number, X::AbstractArray) = generic_scale!(s, X)
+mul2!(a::Number, B::AbstractArray) = generic_mul2!(a, B)
 
 """
     cross(x, y)
@@ -1185,22 +1183,6 @@ function linreg(x::AbstractVector, y::AbstractVector)
     return (a, b)
 end
 
-# multiply by diagonal matrix as vector
-#diagmm!(C::AbstractMatrix, A::AbstractMatrix, b::AbstractVector)
-
-#diagmm!(C::AbstractMatrix, b::AbstractVector, A::AbstractMatrix)
-
-scale!(A::AbstractMatrix, b::AbstractVector) = scale!(A,A,b)
-scale!(b::AbstractVector, A::AbstractMatrix) = scale!(A,b,A)
-
-#diagmm(A::AbstractMatrix, b::AbstractVector)
-#diagmm(b::AbstractVector, A::AbstractMatrix)
-
-#^(A::AbstractMatrix, p::Number)
-
-#findmax(a::AbstractArray)
-#findmin(a::AbstractArray)
-
 """
     peakflops(n::Integer=2000; parallel::Bool=false)
 
@@ -1457,12 +1439,12 @@ end
 
     if nrm ≥ δ # Safe to multiply with inverse
         invnrm = inv(nrm)
-        scale!(v, invnrm)
+        mul1!(v, invnrm)
 
     else # scale elements to avoid overflow
         εδ = eps(one(nrm))/δ
-        scale!(v, εδ)
-        scale!(v, inv(nrm*εδ))
+        mul1!(v, εδ)
+        mul1!(v, inv(nrm*εδ))
     end
 
     v
diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl
index 72f056263bb47..16ca79e14b09a 100644
--- a/stdlib/LinearAlgebra/src/matmul.jl
+++ b/stdlib/LinearAlgebra/src/matmul.jl
@@ -4,38 +4,6 @@
 
 matprod(x, y) = x*y + x*y
 
-# multiply by diagonal matrix as vector
-function scale!(C::AbstractMatrix, A::AbstractMatrix, b::AbstractVector)
-    m, n = size(A)
-    if size(A) != size(C)
-        throw(DimensionMismatch("size of A, $(size(A)), does not match size of C, $(size(C))"))
-    end
-    if n != length(b)
-        throw(DimensionMismatch("second dimension of A, $n, does not match length of b, $(length(b))"))
-    end
-    @inbounds for j = 1:n
-        bj = b[j]
-        for i = 1:m
-            C[i,j] = A[i,j]*bj
-        end
-    end
-    C
-end
-
-function scale!(C::AbstractMatrix, b::AbstractVector, A::AbstractMatrix)
-    m, n = size(A)
-    if size(A) != size(C)
-        throw(DimensionMismatch("size of A, $(size(A)), does not match size of C, $(size(C))"))
-    end
-    if m != length(b)
-        throw(DimensionMismatch("first dimension of A, $m, does not match length of b, $(length(b))"))
-    end
-    @inbounds for j = 1:n, i = 1:m
-        C[i,j] = A[i,j]*b[i]
-    end
-    C
-end
-
 # Dot products
 
 vecdot(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasReal} = BLAS.dot(x, y)
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index bb82bd18fcd34..550b0a01bf2c8 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -400,32 +400,85 @@ function copyto!(A::T, B::T) where T<:Union{LowerTriangular,UnitLowerTriangular}
     return A
 end
 
-function scale!(A::UpperTriangular, B::Union{UpperTriangular,UnitUpperTriangular}, c::Number)
+function mul!(A::UpperTriangular, B::UpperTriangular, c::Number)
     n = checksquare(B)
     for j = 1:n
-        if isa(B, UnitUpperTriangular)
-            @inbounds A[j,j] = c
+        for i = 1:j
+            @inbounds A[i,j] = B[i,j] * c
         end
-        for i = 1:(isa(B, UnitUpperTriangular) ? j-1 : j)
+    end
+    return A
+end
+function mul!(A::UpperTriangular, c::Number, B::UpperTriangular)
+    n = checksquare(B)
+    for j = 1:n
+        for i = 1:j
             @inbounds A[i,j] = c * B[i,j]
         end
     end
     return A
 end
-function scale!(A::LowerTriangular, B::Union{LowerTriangular,UnitLowerTriangular}, c::Number)
+function mul!(A::UpperTriangular, B::UnitUpperTriangular, c::Number)
+    n = checksquare(B)
+    for j = 1:n
+        @inbounds A[j,j] = c
+        for i = 1:(j - 1)
+            @inbounds A[i,j] = B[i,j] * c
+        end
+    end
+    return A
+end
+function mul!(A::UpperTriangular, c::Number, B::UnitUpperTriangular)
+    n = checksquare(B)
+    for j = 1:n
+        @inbounds A[j,j] = c
+        for i = 1:(j - 1)
+            @inbounds A[i,j] = c * B[i,j]
+        end
+    end
+    return A
+end
+function mul!(A::LowerTriangular, B::LowerTriangular, c::Number)
+    n = checksquare(B)
+    for j = 1:n
+        for i = j:n
+            @inbounds A[i,j] = B[i,j] * c
+        end
+    end
+    return A
+end
+function mul!(A::LowerTriangular, c::Number, B::LowerTriangular)
+    n = checksquare(B)
+    for j = 1:n
+        for i = j:n
+            @inbounds A[i,j] = c * B[i,j]
+        end
+    end
+    return A
+end
+function mul!(A::LowerTriangular, B::UnitLowerTriangular, c::Number)
     n = checksquare(B)
     for j = 1:n
-        if isa(B, UnitLowerTriangular)
             @inbounds A[j,j] = c
+        for i = (j + 1):n
+            @inbounds A[i,j] = B[i,j] * c
         end
-        for i = (isa(B, UnitLowerTriangular) ? j+1 : j):n
+    end
+    return A
+end
+function mul!(A::LowerTriangular, c::Number, B::UnitLowerTriangular)
+    n = checksquare(B)
+    for j = 1:n
+            @inbounds A[j,j] = c
+        for i = (j + 1):n
             @inbounds A[i,j] = c * B[i,j]
         end
     end
     return A
 end
-scale!(A::Union{UpperTriangular,LowerTriangular}, c::Number) = scale!(A,A,c)
-scale!(c::Number, A::Union{UpperTriangular,LowerTriangular}) = scale!(A,c)
+
+mul1!(A::Union{UpperTriangular,LowerTriangular}, c::Number) = mul!(A, A, c)
+mul2!(c::Number, A::Union{UpperTriangular,LowerTriangular}) = mul!(A, c, A)
 
 fillstored!(A::LowerTriangular, x)     = (fillband!(A.data, x, 1-size(A,1), 0); A)
 fillstored!(A::UnitLowerTriangular, x) = (fillband!(A.data, x, 1-size(A,1), -1); A)
@@ -1902,7 +1955,7 @@ function powm!(A0::UpperTriangular{<:BlasFloat}, p::Real)
     end
 
     normA0 = norm(A0, 1)
-    scale!(A0, 1/normA0)
+    mul1!(A0, 1/normA0)
 
     theta = [1.53e-5, 2.25e-3, 1.92e-2, 6.08e-2, 1.25e-1, 2.03e-1, 2.84e-1]
     n = checksquare(A0)
@@ -1927,7 +1980,7 @@ function powm!(A0::UpperTriangular{<:BlasFloat}, p::Real)
             @inbounds S[i, i] = S[i, i] + 1
         end
         copyto!(Stmp, S)
-        scale!(S, A, c)
+        mul!(S, A, c)
         ldiv!(Stmp, S.data)
 
         c = (p - j) / (j4 - 2)
@@ -1935,14 +1988,14 @@ function powm!(A0::UpperTriangular{<:BlasFloat}, p::Real)
             @inbounds S[i, i] = S[i, i] + 1
         end
         copyto!(Stmp, S)
-        scale!(S, A, c)
+        mul!(S, A, c)
         ldiv!(Stmp, S.data)
     end
     for i = 1:n
         S[i, i] = S[i, i] + 1
     end
     copyto!(Stmp, S)
-    scale!(S, A, -p)
+    mul!(S, A, -p)
     ldiv!(Stmp, S.data)
     for i = 1:n
         @inbounds S[i, i] = S[i, i] + 1
@@ -1954,7 +2007,7 @@ function powm!(A0::UpperTriangular{<:BlasFloat}, p::Real)
         copyto!(S, Stmp)
         blockpower!(A0, S, p/(2^(s-m)))
     end
-    scale!(S, normA0^p)
+    mul1!(S, normA0^p)
     return S
 end
 powm(A::LowerTriangular, p::Real) = copy(transpose(powm(copy(transpose(A)), p::Real)))
@@ -2119,7 +2172,7 @@ function log(A0::UpperTriangular{T}) where T<:BlasFloat
     end
 
     # Scale back
-    scale!(2^s, Y)
+    mul2!(2^s, Y)
 
     # Compute accurate diagonal and superdiagonal of log(T)
     for k = 1:n-1
diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl
index fcaa3b4a04e3b..f2dbdfa3b74ba 100644
--- a/stdlib/LinearAlgebra/src/uniformscaling.jl
+++ b/stdlib/LinearAlgebra/src/uniformscaling.jl
@@ -195,15 +195,16 @@ end
 *(J::UniformScaling, x::Number) = UniformScaling(J.λ*x)
 
 /(J1::UniformScaling, J2::UniformScaling) = J2.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ/J2.λ)
-/(J::UniformScaling, A::AbstractMatrix) = scale!(J.λ, inv(A))
+/(J::UniformScaling, A::AbstractMatrix) = mul2!(J.λ, inv(A))
 /(A::AbstractMatrix, J::UniformScaling) = J.λ == 0 ? throw(SingularException(1)) : A/J.λ
 
 /(J::UniformScaling, x::Number) = UniformScaling(J.λ/x)
 
 \(J1::UniformScaling, J2::UniformScaling) = J1.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ\J2.λ)
-\(A::Union{Bidiagonal{T},AbstractTriangular{T}}, J::UniformScaling) where {T<:Number} = scale!(inv(A), J.λ)
+\(A::Union{Bidiagonal{T},AbstractTriangular{T}}, J::UniformScaling) where {T<:Number} =
+    mul1!(inv(A), J.λ)
 \(J::UniformScaling, A::AbstractVecOrMat) = J.λ == 0 ? throw(SingularException(1)) : J.λ\A
-\(A::AbstractMatrix, J::UniformScaling) = scale!(inv(A), J.λ)
+\(A::AbstractMatrix, J::UniformScaling) = mul1!(inv(A), J.λ)
 
 \(x::Number, J::UniformScaling) = UniformScaling(x\J.λ)
 
diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl
index d17c25fbd4e11..8bc2d9403cdad 100644
--- a/stdlib/LinearAlgebra/test/diagonal.jl
+++ b/stdlib/LinearAlgebra/test/diagonal.jl
@@ -105,7 +105,7 @@ srand(1)
                 @test ldiv!(transpose(D), copy(U)) ≈ DM\U atol=atol_three
                 @test ldiv!(adjoint(conj(D)), copy(U)) ≈ DM\U atol=atol_three
                 Uc = copy(U')
-                target = scale!(Uc, inv.(D.diag))
+                target = mul1!(Uc, Diagonal(inv.(D.diag)))
                 @test rdiv!(Uc, D) ≈ target atol=atol_three
                 @test_throws DimensionMismatch rdiv!(Matrix{elty}(I, n-1, n-1), D)
                 @test_throws SingularException rdiv!(Uc, Diagonal(fill!(similar(D.diag), 0)))
diff --git a/stdlib/LinearAlgebra/test/generic.jl b/stdlib/LinearAlgebra/test/generic.jl
index 048246392366e..d961181547290 100644
--- a/stdlib/LinearAlgebra/test/generic.jl
+++ b/stdlib/LinearAlgebra/test/generic.jl
@@ -6,7 +6,7 @@ using Test, LinearAlgebra, Random
 import Base: -, *, /, \
 
 # A custom Quaternion type with minimal defined interface and methods.
-# Used to test scale and scale! methods to show non-commutativity.
+# Used to test mul and mul! methods to show non-commutativity.
 struct Quaternion{T<:Real} <: Number
     s::T
     v1::T
@@ -182,37 +182,37 @@ end
     aa = reshape([1.:6;], (2,3))
     for a in (aa, view(aa, 1:2, 1:2))
         am, an = size(a)
-        @testset "2-argument version of scale!" begin
-            @test scale!(copy(a), 5.) == a*5
-            @test scale!(5., copy(a)) == a*5
+        @testset "Scaling with mul1! and mul2" begin
+            @test mul1!(copy(a), 5.) == a*5
+            @test mul2!(5., copy(a)) == a*5
             b = randn(LinearAlgebra.SCAL_CUTOFF) # make sure we try BLAS path
             subB = view(b, :, :)
-            @test scale!(copy(b), 5.) == b*5
-            @test scale!(copy(subB), 5.) == subB*5
-            @test scale!([1.; 2.], copy(a)) == a.*[1; 2]
-            @test scale!([1; 2], copy(a)) == a.*[1; 2]
-            @test scale!(copy(a), Vector(1.:an)) == a.*Vector(1:an)'
-            @test scale!(copy(a), Vector(1:an)) == a.*Vector(1:an)'
-            @test_throws DimensionMismatch scale!(Vector{Float64}(uninitialized,am+1), a)
-            @test_throws DimensionMismatch scale!(a, Vector{Float64}(uninitialized,an+1))
+            @test mul1!(copy(b), 5.) == b*5
+            @test mul1!(copy(subB), 5.) == subB*5
+            @test mul2!(Diagonal([1.; 2.]), copy(a)) == a.*[1; 2]
+            @test mul2!(Diagonal([1; 2]), copy(a)) == a.*[1; 2]
+            @test mul1!(copy(a), Diagonal(1.:an)) == a.*Vector(1:an)'
+            @test mul1!(copy(a), Diagonal(1:an)) == a.*Vector(1:an)'
+            @test_throws DimensionMismatch mul2!(Diagonal(Vector{Float64}(uninitialized,am+1)), a)
+            @test_throws DimensionMismatch mul1!(a, Diagonal(Vector{Float64}(uninitialized,an+1)))
         end
 
-        @testset "3-argument version of scale!" begin
-            @test scale!(similar(a), 5., a) == a*5
-            @test scale!(similar(a), a, 5.) == a*5
-            @test scale!(similar(a), [1.; 2.], a) == a.*[1; 2]
-            @test scale!(similar(a), [1; 2], a) == a.*[1; 2]
-            @test_throws DimensionMismatch scale!(similar(a), Vector{Float64}(uninitialized, am+1), a)
-            @test_throws DimensionMismatch scale!(Matrix{Float64}(uninitialized, 3, 2), a, Vector{Float64}(uninitialized, an+1))
-            @test_throws DimensionMismatch scale!(similar(a), a, Vector{Float64}(uninitialized, an+1))
-            @test scale!(similar(a), a, Vector(1.:an)) == a.*Vector(1:an)'
-            @test scale!(similar(a), a, Vector(1:an)) == a.*Vector(1:an)'
+        @testset "Scaling with 3-argument mul!" begin
+            @test mul!(similar(a), 5., a) == a*5
+            @test mul!(similar(a), a, 5.) == a*5
+            @test mul!(similar(a), Diagonal([1.; 2.]), a) == a.*[1; 2]
+            @test mul!(similar(a), Diagonal([1; 2]), a)   == a.*[1; 2]
+            @test_throws DimensionMismatch mul!(similar(a), Diagonal(Vector{Float64}(uninitialized, am+1)), a)
+            @test_throws DimensionMismatch mul!(Matrix{Float64}(uninitialized, 3, 2), a, Diagonal(Vector{Float64}(uninitialized, an+1)))
+            @test_throws DimensionMismatch mul!(similar(a), a, Diagonal(Vector{Float64}(uninitialized, an+1)))
+            @test mul!(similar(a), a, Diagonal(1.:an)) == a.*Vector(1:an)'
+            @test mul!(similar(a), a, Diagonal(1:an))  == a.*Vector(1:an)'
         end
     end
 end
 
 @testset "scale real matrix by complex type" begin
-    @test_throws InexactError scale!([1.0], 2.0im)
+    @test_throws InexactError mul1!([1.0], 2.0im)
     @test isequal([1.0] * 2.0im,             Complex{Float64}[2.0im])
     @test isequal(2.0im * [1.0],             Complex{Float64}[2.0im])
     @test isequal(Float32[1.0] * 2.0f0im,    Complex{Float32}[2.0im])
@@ -223,11 +223,10 @@ end
     @test isequal(BigFloat[1.0] * 2.0im,     Complex{BigFloat}[2.0im])
     @test isequal(BigFloat[1.0] * 2.0f0im,   Complex{BigFloat}[2.0im])
 end
-@testset "scale and scale! for non-commutative multiplication" begin
+@testset "* and mul! for non-commutative scaling" begin
     q = Quaternion(0.44567, 0.755871, 0.882548, 0.423612)
     qmat = [Quaternion(0.015007, 0.355067, 0.418645, 0.318373)]
-    @test scale!(q, copy(qmat)) != scale!(copy(qmat), q)
-    ## Test * because it doesn't dispatch to scale!
+    @test mul2!(q, copy(qmat)) != mul1!(copy(qmat), q)
     @test q*qmat ≉ qmat*q
     @test conj(q*qmat) ≈ conj(qmat)*conj(q)
     @test q * (q \ qmat) ≈ qmat ≈ (qmat / q) * q
diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl
index 03d65760d2210..ce3dc20c2bcae 100644
--- a/stdlib/LinearAlgebra/test/matmul.jl
+++ b/stdlib/LinearAlgebra/test/matmul.jl
@@ -207,10 +207,10 @@ end
     end
 end
 
-@testset "scale!" begin
+@testset "mul! (scaling)" begin
     A5x5, b5, C5x6 = Array{Float64}.(uninitialized,((5,5), 5, (5,6)))
     for A in (A5x5, view(A5x5, :, :)), b in (b5,  view(b5, :)), C in (C5x6, view(C5x6, :, :))
-        @test_throws DimensionMismatch scale!(A, b, C)
+        @test_throws DimensionMismatch mul!(A, Diagonal(b), C)
     end
 end
 
diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl
index 9a7974e0f9318..0ef12cd0ad61b 100644
--- a/stdlib/LinearAlgebra/test/triangular.jl
+++ b/stdlib/LinearAlgebra/test/triangular.jl
@@ -195,25 +195,25 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo
             ci = cr * im
             if elty1 <: Real
                 A1tmp = copy(A1)
-                scale!(A1tmp,cr)
+                mul1!(A1tmp, cr)
                 @test A1tmp == cr*A1
                 A1tmp = copy(A1)
-                scale!(cr,A1tmp)
+                mul2!(cr, A1tmp)
                 @test A1tmp == cr*A1
                 A1tmp = copy(A1)
                 A2tmp = unitt(A1)
-                scale!(A1tmp,A2tmp,cr)
+                mul!(A1tmp, A2tmp, cr)
                 @test A1tmp == cr * A2tmp
             else
                 A1tmp = copy(A1)
-                scale!(A1tmp,ci)
+                mul1!(A1tmp, ci)
                 @test A1tmp == ci*A1
                 A1tmp = copy(A1)
-                scale!(ci,A1tmp)
+                mul2!(ci, A1tmp)
                 @test A1tmp == ci*A1
                 A1tmp = copy(A1)
                 A2tmp = unitt(A1)
-                scale!(A1tmp,A2tmp,ci)
+                mul!(A1tmp, A2tmp, ci)
                 @test A1tmp == ci * A2tmp
             end
         end
diff --git a/stdlib/SparseArrays/src/SparseArrays.jl b/stdlib/SparseArrays/src/SparseArrays.jl
index 3d34c67d91cc5..714a4ee6c8af8 100644
--- a/stdlib/SparseArrays/src/SparseArrays.jl
+++ b/stdlib/SparseArrays/src/SparseArrays.jl
@@ -14,7 +14,7 @@ using LinearAlgebra
 import Base: +, -, *, \, /, &, |, xor, ==
 import LinearAlgebra: mul!, ldiv!, rdiv!, chol, adjoint!, diag, diff, dot, eig,
     issymmetric, istril, istriu, lu, trace, transpose!, tril!, triu!,
-    vecnorm, cond, diagm, factorize, ishermitian, norm, scale!, tril, triu
+    vecnorm, cond, diagm, factorize, ishermitian, norm, mul1!, mul2!, tril, triu
 
 import Base: @get!, acos, acosd, acot, acotd, acsch, asech, asin, asind, asinh,
     atan, atand, atanh, broadcast!, conj!, cos, cosc, cosd, cosh, cospi, cot,
diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl
index becba99b8c889..fc7e052913774 100644
--- a/stdlib/SparseArrays/src/linalg.jl
+++ b/stdlib/SparseArrays/src/linalg.jl
@@ -36,7 +36,7 @@ function mul!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C:
     nzv = A.nzval
     rv = A.rowval
     if β != 1
-        β != 0 ? scale!(C, β) : fill!(C, zero(eltype(C)))
+        β != 0 ? mul1!(C, β) : fill!(C, zero(eltype(C)))
     end
     for k = 1:size(C, 2)
         for col = 1:A.n
@@ -61,7 +61,7 @@ function mul!(α::Number, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, B::StridedVecO
     nzv = A.nzval
     rv = A.rowval
     if β != 1
-        β != 0 ? scale!(C, β) : fill!(C, zero(eltype(C)))
+        β != 0 ? mul1!(C, β) : fill!(C, zero(eltype(C)))
     end
     for k = 1:size(C, 2)
         for col = 1:A.n
@@ -87,7 +87,7 @@ function mul!(α::Number, transA::Transpose{<:Any,<:SparseMatrixCSC}, B::Strided
     nzv = A.nzval
     rv = A.rowval
     if β != 1
-        β != 0 ? scale!(C, β) : fill!(C, zero(eltype(C)))
+        β != 0 ? mul1!(C, β) : fill!(C, zero(eltype(C)))
     end
     for k = 1:size(C, 2)
         for col = 1:A.n
@@ -128,11 +128,11 @@ end
 
 function (*)(D::Diagonal, A::SparseMatrixCSC)
     T = Base.promote_op(*, eltype(D), eltype(A))
-    scale!(LinearAlgebra.copy_oftype(A, T), D.diag, A)
+    mul!(LinearAlgebra.copy_oftype(A, T), D, A)
 end
 function (*)(A::SparseMatrixCSC, D::Diagonal)
     T = Base.promote_op(*, eltype(D), eltype(A))
-    scale!(LinearAlgebra.copy_oftype(A, T), A, D.diag)
+    mul!(LinearAlgebra.copy_oftype(A, T), A, D)
 end
 
 # Sparse matrix multiplication as described in [Gustavson, 1978]:
@@ -628,7 +628,7 @@ function normestinv(A::SparseMatrixCSC{T}, t::Integer = min(2,maximum(size(A))))
             end
         end
     end
-    scale!(X, 1 ./ n)
+    mul1!(X, inv(n))
 
     iter = 0
     local est
@@ -868,8 +868,9 @@ function copyinds!(C::SparseMatrixCSC, A::SparseMatrixCSC)
 end
 
 # multiply by diagonal matrix as vector
-function scale!(C::SparseMatrixCSC, A::SparseMatrixCSC, b::Vector)
+function mul!(C::SparseMatrixCSC, A::SparseMatrixCSC, D::Diagonal{<:Vector})
     m, n = size(A)
+    b    = D.diag
     (n==length(b) && size(A)==size(C)) || throw(DimensionMismatch())
     copyinds!(C, A)
     Cnzval = C.nzval
@@ -881,8 +882,9 @@ function scale!(C::SparseMatrixCSC, A::SparseMatrixCSC, b::Vector)
     C
 end
 
-function scale!(C::SparseMatrixCSC, b::Vector, A::SparseMatrixCSC)
+function mul!(C::SparseMatrixCSC, D::Diagonal{<:Vector}, A::SparseMatrixCSC)
     m, n = size(A)
+    b    = D.diag
     (m==length(b) && size(A)==size(C)) || throw(DimensionMismatch())
     copyinds!(C, A)
     Cnzval = C.nzval
@@ -895,24 +897,30 @@ function scale!(C::SparseMatrixCSC, b::Vector, A::SparseMatrixCSC)
     C
 end
 
-function scale!(C::SparseMatrixCSC, A::SparseMatrixCSC, b::Number)
+function mul!(C::SparseMatrixCSC, A::SparseMatrixCSC, b::Number)
     size(A)==size(C) || throw(DimensionMismatch())
     copyinds!(C, A)
     resize!(C.nzval, length(A.nzval))
-    scale!(C.nzval, A.nzval, b)
+    mul!(C.nzval, A.nzval, b)
     C
 end
 
-function scale!(C::SparseMatrixCSC, b::Number, A::SparseMatrixCSC)
+function mul!(C::SparseMatrixCSC, b::Number, A::SparseMatrixCSC)
     size(A)==size(C) || throw(DimensionMismatch())
     copyinds!(C, A)
     resize!(C.nzval, length(A.nzval))
-    scale!(C.nzval, b, A.nzval)
+    mul!(C.nzval, b, A.nzval)
     C
 end
 
-scale!(A::SparseMatrixCSC, b::Number) = (scale!(A.nzval, b); A)
-scale!(b::Number, A::SparseMatrixCSC) = (scale!(b, A.nzval); A)
+function mul1!(A::SparseMatrixCSC, b::Number)
+    mul1!(A.nzval, b)
+    return A
+end
+function mul2!(b::Number, A::SparseMatrixCSC)
+    mul2!(b, A.nzval)
+    return A
+end
 
 function \(A::SparseMatrixCSC, B::AbstractVecOrMat)
     m, n = size(A)
@@ -1024,5 +1032,5 @@ function Base.cov(X::SparseMatrixCSC, vardim::Int=1; corrected::Bool=true)
     end
 
     # scale with the sample size n or the corrected sample size n - 1
-    return scale!(out, inv(n - corrected))
+    return mul1!(out, inv(n - corrected))
 end
diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl
index 84cfb1844ccfc..bd4a7dfef2fea 100644
--- a/stdlib/SparseArrays/src/sparsevector.jl
+++ b/stdlib/SparseArrays/src/sparsevector.jl
@@ -1470,12 +1470,24 @@ function LinearAlgebra.axpy!(a::Number, x::SparseVectorUnion, y::AbstractVector)
 end
 
 
-# scale
+# scaling
 
-scale!(x::SparseVectorUnion, a::Real)    = (scale!(nonzeros(x), a); x)
-scale!(x::SparseVectorUnion, a::Complex) = (scale!(nonzeros(x), a); x)
-scale!(a::Real, x::SparseVectorUnion)    = (scale!(nonzeros(x), a); x)
-scale!(a::Complex, x::SparseVectorUnion) = (scale!(nonzeros(x), a); x)
+function mul1!(x::SparseVectorUnion, a::Real)
+    mul1!(nonzeros(x), a)
+    return x
+end
+function mul1!(x::SparseVectorUnion, a::Complex)
+    mul1!(nonzeros(x), a)
+    return x
+end
+function mul2!(a::Real, x::SparseVectorUnion)
+    mul1!(nonzeros(x), a)
+    return x
+end
+function mul2!(a::Complex, x::SparseVectorUnion)
+    mul1!(nonzeros(x), a)
+    return x
+end
 
 (*)(x::SparseVectorUnion, a::Number) = SparseVector(length(x), copy(nonzeroinds(x)), nonzeros(x) * a)
 (*)(a::Number, x::SparseVectorUnion) = SparseVector(length(x), copy(nonzeroinds(x)), a * nonzeros(x))
@@ -1576,7 +1588,7 @@ function mul!(α::Number, A::StridedMatrix, x::AbstractSparseVector, β::Number,
     length(x) == n && length(y) == m || throw(DimensionMismatch())
     m == 0 && return y
     if β != one(β)
-        β == zero(β) ? fill!(y, zero(eltype(y))) : scale!(y, β)
+        β == zero(β) ? fill!(y, zero(eltype(y))) : mul1!(y, β)
     end
     α == zero(α) && return y
 
@@ -1615,7 +1627,7 @@ function mul!(α::Number, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractS
     length(x) == m && length(y) == n || throw(DimensionMismatch())
     n == 0 && return y
     if β != one(β)
-        β == zero(β) ? fill!(y, zero(eltype(y))) : scale!(y, β)
+        β == zero(β) ? fill!(y, zero(eltype(y))) : mul1!(y, β)
     end
     α == zero(α) && return y
 
@@ -1673,7 +1685,7 @@ function mul!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Numbe
     length(x) == n && length(y) == m || throw(DimensionMismatch())
     m == 0 && return y
     if β != one(β)
-        β == zero(β) ? fill!(y, zero(eltype(y))) : scale!(y, β)
+        β == zero(β) ? fill!(y, zero(eltype(y))) : mul1!(y, β)
     end
     α == zero(α) && return y
 
@@ -1717,7 +1729,7 @@ function _At_or_Ac_mul_B!(tfun::Function,
     length(x) == m && length(y) == n || throw(DimensionMismatch())
     n == 0 && return y
     if β != one(β)
-        β == zero(β) ? fill!(y, zero(eltype(y))) : scale!(y, β)
+        β == zero(β) ? fill!(y, zero(eltype(y))) : mul1!(y, β)
     end
     α == zero(α) && return y
 
diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl
index 3bf048920dba5..10eaf1efe9c51 100644
--- a/stdlib/SparseArrays/test/sparse.jl
+++ b/stdlib/SparseArrays/test/sparse.jl
@@ -322,31 +322,31 @@ sA = sprandn(3, 7, 0.5)
 sC = similar(sA)
 dA = Array(sA)
 
-@testset "scale and scale!" begin
+@testset "scaling with * and mul!, mul1!, and mul2!" begin
     b = randn(7)
     @test dA * Diagonal(b) == sA * Diagonal(b)
-    @test dA * Diagonal(b) == scale!(sC, sA, b)
-    @test dA * Diagonal(b) == scale!(copy(sA), b)
+    @test dA * Diagonal(b) == mul!(sC, sA, Diagonal(b))
+    @test dA * Diagonal(b) == mul1!(copy(sA), Diagonal(b))
     b = randn(3)
     @test Diagonal(b) * dA == Diagonal(b) * sA
-    @test Diagonal(b) * dA == scale!(sC, b, sA)
-    @test Diagonal(b) * dA == scale!(b, copy(sA))
+    @test Diagonal(b) * dA == mul!(sC, Diagonal(b), sA)
+    @test Diagonal(b) * dA == mul2!(Diagonal(b), copy(sA))
 
     @test dA * 0.5            == sA * 0.5
-    @test dA * 0.5            == scale!(sC, sA, 0.5)
-    @test dA * 0.5            == scale!(copy(sA), 0.5)
+    @test dA * 0.5            == mul!(sC, sA, 0.5)
+    @test dA * 0.5            == mul1!(copy(sA), 0.5)
     @test 0.5 * dA            == 0.5 * sA
-    @test 0.5 * dA            == scale!(sC, sA, 0.5)
-    @test 0.5 * dA            == scale!(0.5, copy(sA))
-    @test scale!(sC, 0.5, sA) == scale!(sC, sA, 0.5)
+    @test 0.5 * dA            == mul!(sC, sA, 0.5)
+    @test 0.5 * dA            == mul2!(0.5, copy(sA))
+    @test mul!(sC, 0.5, sA)   == mul!(sC, sA, 0.5)
 
-    @testset "inverse scale!" begin
+    @testset "inverse scaling with mul!" begin
         bi = inv.(b)
         dAt = copy(transpose(dA))
         sAt = copy(transpose(sA))
-        @test scale!(copy(dAt), bi) ≈ rdiv!(copy(sAt), Diagonal(b))
-        @test scale!(copy(dAt), bi) ≈ rdiv!(copy(sAt), transpose(Diagonal(b)))
-        @test scale!(copy(dAt), conj(bi)) ≈ rdiv!(copy(sAt), adjoint(Diagonal(b)))
+        @test mul1!(copy(dAt), Diagonal(bi)) ≈ rdiv!(copy(sAt), Diagonal(b))
+        @test mul1!(copy(dAt), Diagonal(bi)) ≈ rdiv!(copy(sAt), transpose(Diagonal(b)))
+        @test mul1!(copy(dAt), Diagonal(conj(bi))) ≈ rdiv!(copy(sAt), adjoint(Diagonal(b)))
         @test_throws DimensionMismatch rdiv!(copy(sAt), Diagonal(fill(1., length(b)+1)))
         @test_throws LinearAlgebra.SingularException rdiv!(copy(sAt), Diagonal(zeros(length(b))))
     end
diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl
index 46cf1105cb31f..1b1bbf051c35d 100644
--- a/stdlib/SparseArrays/test/sparsevector.jl
+++ b/stdlib/SparseArrays/test/sparsevector.jl
@@ -778,16 +778,16 @@ end
             @test exact_equal(x / α, SparseVector(x.n, x.nzind, x.nzval / α))
 
             xc = copy(x)
-            @test scale!(xc, α) === xc
+            @test mul1!(xc, α) === xc
             @test exact_equal(xc, sx)
             xc = copy(x)
-            @test scale!(α, xc) === xc
+            @test mul2!(α, xc) === xc
             @test exact_equal(xc, sx)
             xc = copy(x)
-            @test scale!(xc, complex(α, 0.0)) === xc
+            @test mul1!(xc, complex(α, 0.0)) === xc
             @test exact_equal(xc, sx)
             xc = copy(x)
-            @test scale!(complex(α, 0.0), xc) === xc
+            @test mul2!(complex(α, 0.0), xc) === xc
             @test exact_equal(xc, sx)
         end
 
@@ -1254,7 +1254,7 @@ end
         Aj, Ajview = A[:, j], view(A, :, j)
         @test norm(Aj)          == norm(Ajview)
         @test dot(Aj, copy(Aj)) == dot(Ajview, Aj) # don't alias since it takes a different code path
-        @test scale!(Aj, 0.1)   == scale!(Ajview, 0.1)
+        @test mul1!(Aj, 0.1)    == mul1!(Ajview, 0.1)
         @test Aj*0.1            == Ajview*0.1
         @test 0.1*Aj            == 0.1*Ajview
         @test Aj/0.1            == Ajview/0.1