Skip to content

Commit

Permalink
Clarify which args are modified in mul!
Browse files Browse the repository at this point in the history
  • Loading branch information
@simonbyrne @JeffBezanson
simonbyrne authored and JeffBezanson committed Jan 16, 2018
1 parent aea9155 commit 4fd63eb
Show file tree
Hide file tree
Showing 25 changed files with 459 additions and 440 deletions.
202 changes: 105 additions & 97 deletions base/deprecated.jl
View file

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion base/hashing2.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
hash(s::AbstractString, h::UInt) = hash(String(s), h)
75 changes: 34 additions & 41 deletions base/linalg/diagonal.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,107 +151,97 @@ 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
UpperTriangular(B.data)
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}) =
Diagonal(adjoint.(D.parent.diag) .* adjoint.(B.parent.diag))
*(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
Expand Down Expand Up @@ -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"))
Expand Down Expand Up @@ -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} =
Expand All @@ -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} =
Expand Down
33 changes: 10 additions & 23 deletions base/linalg/givens.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
"""
LinAlg.Givens(i1,i2,c,s) -> G
Expand Down Expand Up @@ -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"))
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)))
19 changes: 9 additions & 10 deletions base/linalg/hessenberg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,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
Expand All @@ -82,31 +82,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
27 changes: 19 additions & 8 deletions base/linalg/linalg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,11 +238,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)).
Expand All @@ -254,11 +252,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)).
Expand All @@ -267,7 +278,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)
Expand Down
Loading

0 comments on commit 4fd63eb

Please sign in to comment.