From 410471f0f5cd31312c955af4d1cdae21b7a03241 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreasnoackjensen@gmail.com>
Date: Fri, 27 May 2016 11:53:56 -0400
Subject: [PATCH] Avoid some unintended calls to generic_matmatmul! for special
 matrices

---
 base/linalg/diagonal.jl | 92 ++++++++++++++++++++++++++++++++---------
 base/linalg/lq.jl       | 34 +++++++++++++--
 test/linalg/diagonal.jl | 13 ++++++
 3 files changed, 116 insertions(+), 23 deletions(-)

diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl
index 3d106b99ef186..560ae57257a5d 100644
--- a/base/linalg/diagonal.jl
+++ b/base/linalg/diagonal.jl
@@ -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)
+    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})
diff --git a/base/linalg/lq.jl b/base/linalg/lq.jl
index 8038033770396..b2a28867f9abb 100644
--- a/base/linalg/lq.jl
+++ b/base/linalg/lq.jl
@@ -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)
@@ -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})
@@ -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!))
+    @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)
diff --git a/test/linalg/diagonal.jl b/test/linalg/diagonal.jl
index 2e6a0883785b8..86394dc7ac1a4 100644
--- a/test/linalg/diagonal.jl
+++ b/test/linalg/diagonal.jl
@@ -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'