Skip to content

Avoid some unintended calls to generic_matmatmul! for special matrices #16615

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Jun 3, 2016
Merged

Avoid some unintended calls to generic_matmatmul! for special matrices #16615

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
92 changes: 73 additions & 19 deletions base/linalg/diagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,31 +104,85 @@ end
/{T<:Number}(D::Diagonal, x::T) = Diagonal(D.diag / x)
*(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .* Db.diag)
*(D::Diagonal, V::AbstractVector) = D.diag .* V
# To avoid ambiguity in the definitions below
for uplo in (:LowerTriangular, :UpperTriangular)
@eval begin
(*)(A::$uplo, D::Diagonal) = $uplo(A.data * D)

function (*)(A::$(Symbol(:Unit, uplo)), D::Diagonal)
B = A.data * D
for i = 1:size(A, 1)
B[i,i] = D.diag[i]
end
return $uplo(B)
end
end
end
(*)(A::AbstractTriangular, D::Diagonal) = error("this method should never be reached")
(*)(D::Diagonal, A::AbstractTriangular) = error("this method should never be reached")

(*)(A::AbstractTriangular, D::Diagonal) = A_mul_B!(copy(A), D)
(*)(D::Diagonal, B::AbstractTriangular) = A_mul_B!(D, copy(B))

(*)(A::AbstractMatrix, D::Diagonal) =
scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), A, D.diag)
(*)(D::Diagonal, A::AbstractMatrix) =
scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), D.diag, A)

A_mul_B!(A::Diagonal,B::AbstractMatrix) = scale!(A.diag,B)
At_mul_B!(A::Diagonal,B::AbstractMatrix)= scale!(A.diag,B)
Ac_mul_B!(A::Diagonal,B::AbstractMatrix)= scale!(conj(A.diag),B)
A_mul_B!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) =
typeof(A)(A_mul_B!(A.data, D))
function A_mul_B!(A::UnitLowerTriangular, D::Diagonal)
A_mul_B!(A.data, D)
for i = 1:size(A, 1)
A.data[i,i] = D.diag[i]
end
LowerTriangular(A.data)
end
function A_mul_B!(A::UnitUpperTriangular, D::Diagonal)
A_mul_B!(A.data, D)
for i = 1:size(A, 1)
A.data[i,i] = D.diag[i]
end
UpperTriangular(A.data)
end
function A_mul_B!(D::Diagonal, B::UnitLowerTriangular)
A_mul_B!(D, B.data)
for i = 1:size(B, 1)
B.data[i,i] = D.diag[i]
end
LowerTriangular(B.data)
end
function A_mul_B!(D::Diagonal, B::UnitUpperTriangular)
A_mul_B!(D, B.data)
for i = 1:size(B, 1)
B.data[i,i] = D.diag[i]
end
UpperTriangular(B.data)
end

Ac_mul_B(A::AbstractTriangular, D::Diagonal) = A_mul_B!(ctranspose(A), D)
function Ac_mul_B(A::AbstractMatrix, D::Diagonal)
Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
ctranspose!(Ac, A)
A_mul_B!(Ac, D)
end

At_mul_B(A::AbstractTriangular, D::Diagonal) = A_mul_B!(transpose(A), D)
function At_mul_B(A::AbstractMatrix, D::Diagonal)
Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
transpose!(Ac, A)
A_mul_B!(Ac, D)
end

A_mul_Bc(D::Diagonal, B::AbstractTriangular) = A_mul_B!(D, ctranspose(B))
A_mul_Bc(D::Diagonal, Q::Union{Base.LinAlg.QRCompactWYQ,Base.LinAlg.QRPackedQ}) = A_mul_Bc!(Array(D), Q)
function A_mul_Bc(D::Diagonal, A::AbstractMatrix)
Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
ctranspose!(Ac, A)
A_mul_B!(D, Ac)
end

A_mul_Bt(D::Diagonal, B::AbstractTriangular) = A_mul_B!(D, transpose(B))
function A_mul_Bt(D::Diagonal, A::AbstractMatrix)
Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
ctranspose!(Ac, A)
Copy link
Contributor

@tkelman tkelman May 27, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

transpose!

A_mul_B!(D, Ac)
end

A_mul_B!(A::Diagonal,B::Diagonal) = throw(MethodError(A_mul_B!, Tuple{Diagonal,Diagonal}))
At_mul_B!(A::Diagonal,B::Diagonal) = throw(MethodError(At_mul_B!, Tuple{Diagonal,Diagonal}))
Ac_mul_B!(A::Diagonal,B::Diagonal) = throw(MethodError(Ac_mul_B!, Tuple{Diagonal,Diagonal}))
A_mul_B!(A::Base.LinAlg.QRPackedQ, D::Diagonal) = throw(MethodError(A_mul_B!, Tuple{Diagonal,Diagonal}))
A_mul_B!(A::Diagonal,B::AbstractMatrix) = scale!(A.diag,B)
At_mul_B!(A::Diagonal,B::AbstractMatrix) = scale!(A.diag,B)
Ac_mul_B!(A::Diagonal,B::AbstractMatrix) = scale!(conj(A.diag),B)
A_mul_B!(A::AbstractMatrix,B::Diagonal) = scale!(A,B.diag)
A_mul_Bt!(A::AbstractMatrix,B::Diagonal) = scale!(A,B.diag)
A_mul_Bc!(A::AbstractMatrix,B::Diagonal) = scale!(A,conj(B.diag))

/(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag ./ Db.diag )
function A_ldiv_B!{T}(D::Diagonal{T}, v::AbstractVector{T})
Expand Down
34 changes: 30 additions & 4 deletions base/linalg/lq.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,16 +106,16 @@ end
## Multiplication by Q
### QB
A_mul_B!{T<:BlasFloat}(A::LQPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormlq!('L','N',A.factors,A.τ,B)
function *{TA,TB}(A::LQPackedQ{TA},B::StridedVecOrMat{TB})
TAB = promote_type(TA, TB)
function (*)(A::LQPackedQ,B::StridedVecOrMat)
TAB = promote_type(eltype(A), eltype(B))
A_mul_B!(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB))
end

### QcB
Ac_mul_B!{T<:BlasReal}(A::LQPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormlq!('L','T',A.factors,A.τ,B)
Ac_mul_B!{T<:BlasComplex}(A::LQPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormlq!('L','C',A.factors,A.τ,B)
function Ac_mul_B{TA,TB}(A::LQPackedQ{TA}, B::StridedVecOrMat{TB})
TAB = promote_type(TA,TB)
function Ac_mul_B(A::LQPackedQ, B::StridedVecOrMat)
TAB = promote_type(eltype(A), eltype(B))
if size(B,1) == size(A.factors,2)
Ac_mul_B!(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB))
elseif size(B,1) == size(A.factors,1)
Expand All @@ -125,6 +125,19 @@ function Ac_mul_B{TA,TB}(A::LQPackedQ{TA}, B::StridedVecOrMat{TB})
end
end

### QBc/QcBc
for (f1, f2) in ((:A_mul_Bc, :A_mul_B!),
(:Ac_mul_Bc, :Ac_mul_B!))
@eval begin
function ($f1)(A::LQPackedQ, B::StridedVecOrMat)
TAB = promote_type(eltype(A), eltype(B))
BB = similar(B, TAB, (size(B, 2), size(B, 1)))
ctranspose!(BB, B)
return ($f2)(A, BB)
end
end
end

### AQ
A_mul_B!{T<:BlasFloat}(A::StridedMatrix{T}, B::LQPackedQ{T}) = LAPACK.ormlq!('R', 'N', B.factors, B.τ, A)
function *{TA,TB}(A::StridedMatrix{TA},B::LQPackedQ{TB})
Expand All @@ -146,6 +159,19 @@ function A_mul_Bc{TA<:Number,TB<:Number}( A::StridedVecOrMat{TA}, B::LQPackedQ{T
A_mul_Bc!(copy_oftype(A, TAB), convert(AbstractMatrix{TAB},(B)))
end

### AcQ/AcQc
for (f1, f2) in ((:Ac_mul_B, :A_mul_B!),
(:Ac_mul_Bc, :A_mul_Bc!))
Copy link
Contributor

@tkelman tkelman May 27, 2016 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

mutating ones should also be for Ac?

and won't mutate A. which are they supposed to modify?

Copy link
Member Author

@andreasnoack andreasnoack May 27, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, because I explicitly transpose A into new memory which will become the return array.

Copy link
Contributor

@tkelman tkelman May 27, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah right I think I was reading this as a loop of 4 definitions, sorry

@eval begin
function ($f1)(A::StridedMatrix, B::LQPackedQ)
TAB = promote_type(eltype(A), eltype(B))
AA = similar(A, TAB, (size(A, 2), size(A, 1)))
ctranspose!(AA, A)
return ($f2)(AA, B)
end
end
end

function \{TA,Tb}(A::LQ{TA}, b::StridedVector{Tb})
S = promote_type(TA,Tb)
m = checksquare(A)
Expand Down
13 changes: 13 additions & 0 deletions test/linalg/diagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,3 +246,16 @@ end

# Issue 15401
@test eye(5) \ Diagonal(ones(5)) == eye(5)

# Triangular and Diagonal
for T in (LowerTriangular(randn(5,5)), LinAlg.UnitLowerTriangular(randn(5,5)))
D = Diagonal(randn(5))
@test T'D == Array(T)'*Array(D)
@test T.'D == Array(T).'*Array(D)
@test D*T' == Array(D)*Array(T)'
@test D*T.' == Array(D)*Array(T).'
end

# Diagonal and Q
Q = qrfact(randn(5,5))[:Q]
@test D*Q' == Array(D)*Q'