Deprecate scale! in favor of mul!, mul1!, and mul2! (JuliaLang#25701)

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

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 ]

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!
 ```

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,

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)

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)

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

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

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)