Skip to content

Commit 57fbd18

Browse files
simonbyrneJeffBezanson
simonbyrne
authored and
JeffBezanson
committed
Clarify which args are modified in mul!
1 parent 19429e2 commit 57fbd18

25 files changed

+453
-434
lines changed

‎base/deprecated.jl

Lines changed: 105 additions & 97 deletions
Large diffs are not rendered by default.

‎base/hashing2.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -178,4 +178,4 @@ function hash(s::Union{String,SubString{String}}, h::UInt)
178178
# note: use pointer(s) here (see #6058).
179179
ccall(memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), pointer(s), sizeof(s), h % UInt32) + h
180180
end
181-
hash(s::AbstractString, h::UInt) = hash(String(s), h)
181+
hash(s::AbstractString, h::UInt) = hash(String(s), h)

‎base/linalg/diagonal.jl

Lines changed: 28 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -151,101 +151,91 @@ end
151151
(*)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .* Db.diag)
152152
(*)(D::Diagonal, V::AbstractVector) = D.diag .* V
153153

154-
(*)(A::AbstractTriangular, D::Diagonal) = mul!(copy(A), D)
155-
(*)(D::Diagonal, B::AbstractTriangular) = mul!(D, copy(B))
154+
(*)(A::AbstractTriangular, D::Diagonal) = mul1!(copy(A), D)
155+
(*)(D::Diagonal, B::AbstractTriangular) = mul2!(D, copy(B))
156156

157157
(*)(A::AbstractMatrix, D::Diagonal) =
158158
scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), A, D.diag)
159159
(*)(D::Diagonal, A::AbstractMatrix) =
160160
scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), D.diag, A)
161161

162-
mul!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) = typeof(A)(mul!(A.data, D))
163-
function mul!(A::UnitLowerTriangular, D::Diagonal)
164-
mul!(A.data, D)
162+
163+
mul1!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) = typeof(A)(mul1!(A.data, D))
164+
function mul1!(A::UnitLowerTriangular, D::Diagonal)
165+
mul1!(A.data, D)
165166
for i = 1:size(A, 1)
166167
A.data[i,i] = D.diag[i]
167168
end
168169
LowerTriangular(A.data)
169170
end
170-
function mul!(A::UnitUpperTriangular, D::Diagonal)
171-
mul!(A.data, D)
171+
function mul1!(A::UnitUpperTriangular, D::Diagonal)
172+
mul1!(A.data, D)
172173
for i = 1:size(A, 1)
173174
A.data[i,i] = D.diag[i]
174175
end
175176
UpperTriangular(A.data)
176177
end
177-
function mul!(D::Diagonal, B::UnitLowerTriangular)
178-
mul!(D, B.data)
178+
179+
function mul2!(D::Diagonal, B::UnitLowerTriangular)
180+
mul2!(D, B.data)
179181
for i = 1:size(B, 1)
180182
B.data[i,i] = D.diag[i]
181183
end
182184
LowerTriangular(B.data)
183185
end
184-
function mul!(D::Diagonal, B::UnitUpperTriangular)
185-
mul!(D, B.data)
186+
function mul2!(D::Diagonal, B::UnitUpperTriangular)
187+
mul2!(D, B.data)
186188
for i = 1:size(B, 1)
187189
B.data[i,i] = D.diag[i]
188190
end
189191
UpperTriangular(B.data)
190192
end
191193

192194
*(D::Adjoint{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(adjoint.(D.parent.diag) .* B.diag)
193-
*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) = mul!(copy(A), D)
195+
*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) = mul1!(copy(A), D)
194196
function *(adjA::Adjoint{<:Any,<:AbstractMatrix}, D::Diagonal)
195197
A = adjA.parent
196198
Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
197199
adjoint!(Ac, A)
198-
mul!(Ac, D)
200+
mul1!(Ac, D)
199201
end
200202

201203
*(D::Transpose{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(transpose.(D.parent.diag) .* B.diag)
202-
*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) = mul!(copy(A), D)
204+
*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) = mul1!(copy(A), D)
203205
function *(transA::Transpose{<:Any,<:AbstractMatrix}, D::Diagonal)
204206
A = transA.parent
205207
At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
206208
transpose!(At, A)
207-
mul!(At, D)
209+
mul1!(At, D)
208210
end
209211

210212
*(D::Diagonal, B::Adjoint{<:Any,<:Diagonal}) = Diagonal(D.diag .* adjoint.(B.parent.diag))
211-
*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) = mul!(D, collect(B))
212-
*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; mul!(Array(D), adjoint(Q)))
213+
*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) = mul2!(D, collect(B))
214+
*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; mul1!(Array(D), adjoint(Q)))
213215
function *(D::Diagonal, adjA::Adjoint{<:Any,<:AbstractMatrix})
214216
A = adjA.parent
215217
Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
216218
adjoint!(Ac, A)
217-
mul!(D, Ac)
219+
mul2!(D, Ac)
218220
end
219221

220222
*(D::Diagonal, B::Transpose{<:Any,<:Diagonal}) = Diagonal(D.diag .* transpose.(B.parent.diag))
221-
*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) = mul!(D, copy(B))
223+
*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) = mul2!(D, copy(B))
222224
function *(D::Diagonal, transA::Transpose{<:Any,<:AbstractMatrix})
223225
A = transA.parent
224226
At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
225227
transpose!(At, A)
226-
mul!(D, At)
228+
mul2!(D, At)
227229
end
228230

229231
*(D::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) =
230232
Diagonal(adjoint.(D.parent.diag) .* adjoint.(B.parent.diag))
231233
*(D::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) =
232234
Diagonal(transpose.(D.parent.diag) .* transpose.(B.parent.diag))
233235

234-
mul!(A::Diagonal, B::Diagonal) = throw(MethodError(mul!, (A, B)))
235-
mul!(A::QRPackedQ, D::Diagonal) = throw(MethodError(mul!, (A, D)))
236-
mul!(A::QRPackedQ, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
237-
mul!(A::QRPackedQ, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
238-
mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Diagonal) = throw(MethodError(mul!, (A, B)))
239-
mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
240-
mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
241-
mul!(A::Diagonal, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
242-
mul!(A::Diagonal, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
243-
mul!(A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
244-
mul!(A::Adjoint{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
245-
mul!(A::Transpose{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
246-
mul!(A::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) = throw(MethodError(mul!, (A, B)))
247-
mul!(A::Transpose{<:Any,<:Diagonal}, B::Diagonal) = throw(MethodError(mul!, (A, B)))
248-
mul!(A::Adjoint{<:Any,<:Diagonal}, B::Diagonal) = throw(MethodError(mul!, (A, B)))
236+
mul1!(A::Diagonal,B::Diagonal) = Diagonal(A.diag .*= B.diag)
237+
mul2!(A::Diagonal,B::Diagonal) = Diagonal(B.diag .= A.diag .* B.diag)
238+
249239
mul!(A::Diagonal, B::AbstractMatrix) = scale!(A.diag, B)
250240
mul!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = adjA.parent; scale!(conj(A.diag),B))
251241
mul!(transA::Transpose{<:Any,<:Diagonal}, B::AbstractMatrix) = (A = transA.parent; scale!(A.diag,B))
@@ -287,6 +277,7 @@ mul!(C::AbstractMatrix, A::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Rea
287277

288278

289279
(/)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag ./ Db.diag)
280+
290281
function ldiv!(D::Diagonal{T}, v::AbstractVector{T}) where {T}
291282
if length(v) != length(D.diag)
292283
throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $(length(v)) rows"))
@@ -316,6 +307,7 @@ function ldiv!(D::Diagonal{T}, V::AbstractMatrix{T}) where {T}
316307
V
317308
end
318309

310+
319311
ldiv!(adjD::Adjoint{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} =
320312
(D = adjD.parent; ldiv!(conj(D), B))
321313
ldiv!(transD::Transpose{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} =
@@ -339,6 +331,7 @@ function rdiv!(A::AbstractMatrix{T}, D::Diagonal{T}) where {T}
339331
A
340332
end
341333

334+
342335
rdiv!(A::AbstractMatrix{T}, adjD::Adjoint{<:Any,<:Diagonal{T}}) where {T} =
343336
(D = adjD.parent; rdiv!(A, conj(D)))
344337
rdiv!(A::AbstractMatrix{T}, transD::Transpose{<:Any,<:Diagonal{T}}) where {T} =

‎base/linalg/givens.jl

Lines changed: 10 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -7,14 +7,14 @@ transpose(R::AbstractRotation) = error("transpose not implemented for $(typeof(R
77

88
function *(R::AbstractRotation{T}, A::AbstractVecOrMat{S}) where {T,S}
99
TS = typeof(zero(T)*zero(S) + zero(T)*zero(S))
10-
mul!(convert(AbstractRotation{TS}, R), TS == S ? copy(A) : convert(AbstractArray{TS}, A))
10+
mul2!(convert(AbstractRotation{TS}, R), TS == S ? copy(A) : convert(AbstractArray{TS}, A))
1111
end
1212
*(A::AbstractVector, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)
1313
*(A::AbstractMatrix, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)
1414
function _absvecormat_mul_adjrot(A::AbstractVecOrMat{T}, adjR::Adjoint{<:Any,<:AbstractRotation{S}}) where {T,S}
1515
R = adjR.parent
1616
TS = typeof(zero(T)*zero(S) + zero(T)*zero(S))
17-
mul!(TS == T ? copy(A) : convert(AbstractArray{TS}, A), adjoint(convert(AbstractRotation{TS}, R)))
17+
mul1!(TS == T ? copy(A) : convert(AbstractArray{TS}, A), adjoint(convert(AbstractRotation{TS}, R)))
1818
end
1919
"""
2020
LinAlg.Givens(i1,i2,c,s) -> G
@@ -325,10 +325,7 @@ function getindex(G::Givens, i::Integer, j::Integer)
325325
end
326326
end
327327

328-
329-
mul!(G1::Givens, G2::Givens) = error("Operation not supported. Consider *")
330-
331-
function mul!(G::Givens, A::AbstractVecOrMat)
328+
function mul2!(G::Givens, A::AbstractVecOrMat)
332329
m, n = size(A, 1), size(A, 2)
333330
if G.i2 > m
334331
throw(DimensionMismatch("column indices for rotation are outside the matrix"))
@@ -340,7 +337,7 @@ function mul!(G::Givens, A::AbstractVecOrMat)
340337
end
341338
return A
342339
end
343-
function mul!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens})
340+
function mul1!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens})
344341
G = adjG.parent
345342
m, n = size(A, 1), size(A, 2)
346343
if G.i2 > n
@@ -353,20 +350,21 @@ function mul!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens})
353350
end
354351
return A
355352
end
356-
function mul!(G::Givens, R::Rotation)
353+
354+
function mul2!(G::Givens, R::Rotation)
357355
push!(R.rotations, G)
358356
return R
359357
end
360-
function mul!(R::Rotation, A::AbstractMatrix)
358+
function mul2!(R::Rotation, A::AbstractMatrix)
361359
@inbounds for i = 1:length(R.rotations)
362-
mul!(R.rotations[i], A)
360+
mul2!(R.rotations[i], A)
363361
end
364362
return A
365363
end
366-
function mul!(A::AbstractMatrix, adjR::Adjoint{<:Any,<:Rotation})
364+
function mul1!(A::AbstractMatrix, adjR::Adjoint{<:Any,<:Rotation})
367365
R = adjR.parent
368366
@inbounds for i = 1:length(R.rotations)
369-
mul!(A, adjoint(R.rotations[i]))
367+
mul1!(A, adjoint(R.rotations[i]))
370368
end
371369
return A
372370
end
@@ -388,14 +386,3 @@ end
388386
# dismabiguation methods: *(Diag/AbsTri, Adj of AbstractRotation)
389387
*(A::Diagonal, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B)
390388
*(A::AbstractTriangular, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B)
391-
# moar disambiguation
392-
mul!(A::QRPackedQ, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))
393-
mul!(A::QRPackedQ, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
394-
mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))
395-
mul!(A::Adjoint{<:Any,<:QRPackedQ}, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
396-
mul!(A::Diagonal, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))
397-
mul!(A::Diagonal, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
398-
mul!(A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))
399-
mul!(A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
400-
mul!(A::Transpose{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Rotation}) = throw(MethodError(mul!, (A, B)))
401-
mul!(A::Transpose{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Givens}) = throw(MethodError(mul!, (A, B)))

‎base/linalg/hessenberg.jl

Lines changed: 9 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,7 @@ function getindex(A::HessenbergQ, i::Integer, j::Integer)
7171
x[i] = 1
7272
y = zeros(eltype(A), size(A, 2))
7373
y[j] = 1
74-
return dot(x, mul!(A, y))
74+
return dot(x, mul2!(A, y))
7575
end
7676

7777
## reconstruct the original matrix
@@ -82,31 +82,30 @@ AbstractArray(F::Hessenberg) = AbstractMatrix(F)
8282
Matrix(F::Hessenberg) = Array(AbstractArray(F))
8383
Array(F::Hessenberg) = Matrix(F)
8484

85-
mul!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
85+
mul2!(Q::HessenbergQ{T}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
8686
LAPACK.ormhr!('L', 'N', 1, size(Q.factors, 1), Q.factors, Q.τ, X)
87-
mul!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} =
87+
mul1!(X::StridedMatrix{T}, Q::HessenbergQ{T}) where {T<:BlasFloat} =
8888
LAPACK.ormhr!('R', 'N', 1, size(Q.factors, 1), Q.factors, Q.τ, X)
89-
mul!(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
89+
mul2!(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
9090
(Q = adjQ.parent; LAPACK.ormhr!('L', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X))
91-
mul!(X::StridedMatrix{T}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T<:BlasFloat} =
91+
mul1!(X::StridedMatrix{T}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T<:BlasFloat} =
9292
(Q = adjQ.parent; LAPACK.ormhr!('R', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X))
9393

94-
9594
function (*)(Q::HessenbergQ{T}, X::StridedVecOrMat{S}) where {T,S}
9695
TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
97-
return mul!(Q, copy_oftype(X, TT))
96+
return mul2!(Q, copy_oftype(X, TT))
9897
end
9998
function (*)(X::StridedVecOrMat{S}, Q::HessenbergQ{T}) where {T,S}
10099
TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
101-
return mul!(copy_oftype(X, TT), Q)
100+
return mul1!(copy_oftype(X, TT), Q)
102101
end
103102
function *(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{S}) where {T,S}
104103
Q = adjQ.parent
105104
TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
106-
return mul!(adjoint(Q), copy_oftype(X, TT))
105+
return mul2!(adjoint(Q), copy_oftype(X, TT))
107106
end
108107
function *(X::StridedVecOrMat{S}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T,S}
109108
Q = adjQ.parent
110109
TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
111-
return mul!(copy_oftype(X, TT), adjoint(Q))
110+
return mul1!(copy_oftype(X, TT), adjoint(Q))
112111
end

‎base/linalg/linalg.jl

Lines changed: 19 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -238,11 +238,9 @@ function char_uplo(uplo::Symbol)
238238
end
239239

240240
"""
241-
ldiv!([Y,] A, B) -> Y
241+
ldiv!(Y, A, B) -> Y
242242
243243
Compute `A \\ B` in-place and store the result in `Y`, returning the result.
244-
If only two arguments are passed, then `ldiv!(A, B)` overwrites `B` with
245-
the result.
246244
247245
The argument `A` should *not* be a matrix. Rather, instead of matrices it should be a
248246
factorization object (e.g. produced by [`factorize`](@ref) or [`cholfact`](@ref)).
@@ -254,11 +252,24 @@ control over the factorization of `A`.
254252
ldiv!(Y, A, B)
255253

256254
"""
257-
rdiv!([Y,] A, B) -> Y
255+
ldiv!(A, B)
258256
259-
Compute `A / B` in-place and store the result in `Y`, returning the result.
260-
If only two arguments are passed, then `rdiv!(A, B)` overwrites `A` with
261-
the result.
257+
Compute `A \\ B` in-place and overwriting `B` to store the result.
258+
259+
The argument `A` should *not* be a matrix. Rather, instead of matrices it should be a
260+
factorization object (e.g. produced by [`factorize`](@ref) or [`cholfact`](@ref)).
261+
The reason for this is that factorization itself is both expensive and typically allocates memory
262+
(although it can also be done in-place via, e.g., [`lufact!`](@ref)),
263+
and performance-critical situations requiring `ldiv!` usually also require fine-grained
264+
control over the factorization of `A`.
265+
"""
266+
ldiv!(A, B)
267+
268+
269+
"""
270+
rdiv!(A, B)
271+
272+
Compute `A / B` in-place and overwriting `A` to store the result.
262273
263274
The argument `B` should *not* be a matrix. Rather, instead of matrices it should be a
264275
factorization object (e.g. produced by [`factorize`](@ref) or [`cholfact`](@ref)).
@@ -267,7 +278,7 @@ The reason for this is that factorization itself is both expensive and typically
267278
and performance-critical situations requiring `rdiv!` usually also require fine-grained
268279
control over the factorization of `B`.
269280
"""
270-
rdiv!(Y, A, B)
281+
rdiv!(A, B)
271282

272283
copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A)
273284
copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A)

0 commit comments

Comments
 (0)