From 7e76526e49dc0d525a9c583bc147c8ae951804d7 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Fri, 5 Feb 2021 10:41:47 +0100
Subject: [PATCH] specialize copyto! and multiplication by numbers for Q from
 qr

This fixes two performance bugs reported in https://github.com/JuliaLang/julia/issues/38972
and https://github.com/JuliaLang/julia/issues/38972 (multiplication of `Q` from `qr` by a
`Diagonal` or `UniformScaling`). In particular, it improves the performance of generating
random orthogonal matrices as described in https://discourse.julialang.org/t/random-orthogonal-matrices/9779/7.
---
 stdlib/LinearAlgebra/src/qr.jl  | 20 ++++++++++++++++++++
 stdlib/LinearAlgebra/test/qr.jl | 32 ++++++++++++++++++++++++++++++++
 2 files changed, 52 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl
index a76577bb63a0d..123787cac9e30 100644
--- a/stdlib/LinearAlgebra/src/qr.jl
+++ b/stdlib/LinearAlgebra/src/qr.jl
@@ -533,6 +533,12 @@ function getindex(Q::AbstractQ, i::Integer, j::Integer)
     return dot(x, lmul!(Q, y))
 end
 
+# specialization avoiding the fallback using slow `getindex`
+function copyto!(dest::AbstractMatrix, src::LinearAlgebra.AbstractQ)
+    copyto!(dest, I)
+    lmul!(src, dest)
+end
+
 ## Multiplication by Q
 ### QB
 lmul!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} =
@@ -590,6 +596,13 @@ function (*)(A::AbstractQ, B::StridedMatrix)
     lmul!(Anew, Bnew)
 end
 
+function (*)(A::AbstractQ, b::Number)
+    TAb = promote_type(eltype(A), typeof(b))
+    dest = similar(A, TAb)
+    copyto!(dest, b*I)
+    lmul!(A, dest)
+end
+
 ### QcB
 lmul!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} =
     (A = adjA.parent; LAPACK.gemqrt!('L','T',A.factors,A.T,B))
@@ -683,6 +696,13 @@ function (*)(A::StridedMatrix, Q::AbstractQ)
     return rmul!(copy_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q))
 end
 
+function (*)(a::Number, B::AbstractQ)
+    TaB = promote_type(typeof(a), eltype(B))
+    dest = similar(B, TaB)
+    copyto!(dest, a*I)
+    rmul!(dest, B)
+end
+
 ### AQc
 rmul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasReal} =
     (B = adjB.parent; LAPACK.gemqrt!('R','T',B.factors,B.T,A))
diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl
index b4e6d383bb262..d9c5c4469c363 100644
--- a/stdlib/LinearAlgebra/test/qr.jl
+++ b/stdlib/LinearAlgebra/test/qr.jl
@@ -322,4 +322,36 @@ end
     end
 end
 
+@testset "QR factorization of Q" begin
+    for T in (Float32, Float64, ComplexF32, ComplexF64)
+        Q1, R1 = qr(randn(T,5,5))
+        Q2, R2 = qr(Q1)
+        @test Q1 ≈ Q2
+        @test R2 ≈ I
+    end
+end
+
+@testset "Generation of orthogonal matrices" begin
+    for T in (Float32, Float64)
+        n = 5
+        Q, R = qr(randn(T,n,n))
+        O = Q * Diagonal(sign.(diag(R)))
+        @test O' * O ≈ I
+    end
+end
+
+@testset "Multiplication of Q by special matrices" begin
+    for T in (Float32, Float64, ComplexF32, ComplexF64)
+        n = 5
+        Q, R = qr(randn(T,n,n))
+        Qmat = Matrix(Q)
+        D = Diagonal(randn(n))
+        @test Q * D ≈ Qmat * D
+        @test D * Q ≈ D * Qmat
+        J = 2*I
+        @test Q * J ≈ Qmat * J
+        @test J * Q ≈ J * Qmat
+    end
+end
+
 end # module TestQR