Deprecate scale! in favor of mul!, mul1!, and mul2!

Author: andreasnoack

Diff for: stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl

@@ -9,7 +9,7 @@ module IterativeEigensolvers
 
 using LinearAlgebra: BlasFloat, BlasInt, SVD, checksquare, mul!,
               UniformScaling, issymmetric, ishermitian,
-              factorize, I, scale!, qr
+              factorize, I, 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 for: stdlib/LinearAlgebra/docs/src/index.md

@@ -354,7 +354,8 @@ LinearAlgebra.tril!
 LinearAlgebra.diagind
 LinearAlgebra.diag
 LinearAlgebra.diagm
-LinearAlgebra.scale!
+LinearAlgebra.mul1!
+LinearAlgebra.mul2!
 LinearAlgebra.rank
 LinearAlgebra.norm
 LinearAlgebra.vecnorm

Diff for: 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 for: 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 for: 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!(X, A, Diagonal(b))
+@deprecate scale!(C::AbstractMatrix, a::AbstractVector, B::AbstractMatrix) mul!(X, Diagonal(a), B)

Diff for: 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 for: 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, b)
 
 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]
 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!(X::AbstractArray, s::Number) = generic_mul1!(X, s)
 
-julia> a = [1 2; 3 4];
+"""
+    mul2!(A, b)
 
-julia> b = [1; 2];
+Scale an array `A` by a scalar `b` overwriting `A` in-place.
 
-julia> scale!(b,a)
+# Examples
+```jldoctest
+julia> a = [1 2; 3 4]
 2×2 Array{Int64,2}:
  1  2
+ 3  4
+
+julia> mul2!(2, a)
+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!(s::Number, X::AbstractArray) = generic_mul2!(s, X)
 
 """
     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 for: 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)