Skip to content

Commit

Permalink
Deprecate scale! in favor of mul!, mul1!, and mul2! (#25701)
Browse files Browse the repository at this point in the history
  • Loading branch information
andreasnoack authored and JeffBezanson committed Jan 24, 2018
1 parent 27407e0 commit cf6303a
Show file tree
Hide file tree
Showing 20 changed files with 258 additions and 200 deletions.
2 changes: 2 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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]).

Expand Down
10 changes: 5 additions & 5 deletions stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 ]
Expand Down
3 changes: 2 additions & 1 deletion stdlib/LinearAlgebra/docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -354,7 +354,6 @@ LinearAlgebra.tril!
LinearAlgebra.diagind
LinearAlgebra.diag
LinearAlgebra.diagm
LinearAlgebra.scale!
LinearAlgebra.rank
LinearAlgebra.norm
LinearAlgebra.vecnorm
Expand Down Expand Up @@ -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!
```
Expand Down
6 changes: 5 additions & 1 deletion stdlib/LinearAlgebra/src/LinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ export
istril,
istriu,
kron,
ldiv!,
ldltfact!,
ldltfact,
linreg,
Expand All @@ -107,6 +108,9 @@ export
lufact,
lufact!,
lyap,
mul!,
mul1!,
mul2!,
norm,
normalize,
normalize!,
Expand All @@ -122,7 +126,7 @@ export
lqfact!,
lqfact,
rank,
scale!,
rdiv!,
schur,
schurfact!,
schurfact,
Expand Down
14 changes: 7 additions & 7 deletions stdlib/LinearAlgebra/src/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))

Expand Down Expand Up @@ -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)
8 changes: 8 additions & 0 deletions stdlib/LinearAlgebra/src/deprecated.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
40 changes: 30 additions & 10 deletions stdlib/LinearAlgebra/src/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
76 changes: 29 additions & 47 deletions stdlib/LinearAlgebra/src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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))."))
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
32 changes: 0 additions & 32 deletions stdlib/LinearAlgebra/src/matmul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading

0 comments on commit cf6303a

Please sign in to comment.