From e70354da556fbf7b08c73110380abc0c34216c1e Mon Sep 17 00:00:00 2001
From: Aravindh Krishnamoorthy <aravindh.krishnamoorthy@fau.de>
Date: Fri, 5 May 2023 08:40:01 +0200
Subject: [PATCH] Fix for #48911: Bugfix for the generic left-divide (ldiv!)
 algorithm with pivoted QR decomposition and wide matrices (#49570)

* Fix for #48911

1. Fix output permutation when A.P is not the identity matrix.
2. Fix the application of the elementary transformation (correct complex conjugation) to the solution as per LAPACK zlarzb https://netlib.org/lapack/explore-html/d5/ddd/zlarzb_8f_source.html
3. Align the final permutation with LAPACK's ZGELSY at https://netlib.org/lapack/explore-html/d8/d6e/zgelsy_8f_source.html
4. Fix a complex conjugate bug in the original solution for #48911, align with LAPACK zlarzb, https://netlib.org/lapack/explore-html/d5/ddd/zlarzb_8f_source.html
5. Improve permutation performance

* Integration fix for #48911: Permutation is taken care by ldiv!(QRPivoted, ...)

* Tests for #48911: LinearAlgebra/test/qr.jl

* Issue #48911: Include original test case mentioned in the issue.

---------

Co-authored-by: Daniel Karrasch <daniel.karrasch@posteo.de>
---
 stdlib/LinearAlgebra/src/qr.jl  |  9 ++++-----
 stdlib/LinearAlgebra/test/qr.jl | 30 ++++++++++++++++++++++++++++++
 2 files changed, 34 insertions(+), 5 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl
index 4e3baa7ca91c4..43d04ac5fa415 100644
--- a/stdlib/LinearAlgebra/src/qr.jl
+++ b/stdlib/LinearAlgebra/src/qr.jl
@@ -614,7 +614,6 @@ ldiv!(A::QRPivoted{T,<:StridedMatrix}, B::AbstractVector{T}) where {T<:BlasFloat
 ldiv!(A::QRPivoted{T,<:StridedMatrix}, B::AbstractMatrix{T}) where {T<:BlasFloat} =
     ldiv!(A, B, min(size(A)...)*eps(real(T)))[1]
 
-
 function _wide_qr_ldiv!(A::QR{T}, B::AbstractMatrix{T}) where T
     m, n = size(A)
     minmn = min(m,n)
@@ -646,14 +645,14 @@ function _wide_qr_ldiv!(A::QR{T}, B::AbstractMatrix{T}) where T
             B[m + 1:mB,1:nB] .= zero(T)
             for j = 1:nB
                 for k = 1:m
-                    vBj = B[k,j]
+                    vBj = B[k,j]'
                     for i = m + 1:n
-                        vBj += B[i,j]*R[k,i]'
+                        vBj += B[i,j]'*R[k,i]'
                     end
                     vBj *= τ[k]
-                    B[k,j] -= vBj
+                    B[k,j] -= vBj'
                     for i = m + 1:n
-                        B[i,j] -= R[k,i]*vBj
+                        B[i,j] -= R[k,i]'*vBj'
                     end
                 end
             end
diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl
index cb15d94d08dcb..6e2e9a7b20603 100644
--- a/stdlib/LinearAlgebra/test/qr.jl
+++ b/stdlib/LinearAlgebra/test/qr.jl
@@ -474,4 +474,34 @@ end
     @test MyIdentity{Float64}()[1,:] == [1.0, 0.0]
 end
 
+@testset "issue #48911" begin
+    # testcase in the original issue
+    # test ldiv!(::QRPivoted, ::AbstractVector)
+    A = Complex{BigFloat}[1+im 1-im]
+    b = Complex{BigFloat}[3+im]
+    x = A\b
+    AF = Complex{Float64}[1+im 1-im]
+    bf = Complex{Float64}[3+im]
+    xf = AF\bf
+    @test x ≈ xf
+
+    # test ldiv!(::QRPivoted, ::AbstractVector)
+    A = Complex{BigFloat}[1+im 2-2im 3+3im; 4-4im 5+5im 6-6im]
+    b = Complex{BigFloat}[1+im; 0]
+    x = A\b
+    AF = Complex{Float64}[1+im 2-2im 3+3im; 4-4im 5+5im 6-6im]
+    bf = Complex{Float64}[1+im; 0]
+    xf = AF\bf
+    @test x ≈ xf
+
+    # test ldiv!(::QRPivoted, ::AbstractMatrix)
+    C = Complex{BigFloat}[1+im 2-2im 3+3im; 4-4im 5+5im 6-6im]
+    D = Complex{BigFloat}[1+im 1-im; 0 0]
+    x = C\D
+    CF = Complex{Float64}[1+im 2-2im 3+3im; 4-4im 5+5im 6-6im]
+    DF = Complex{Float64}[1+im 1-im; 0 0]
+    xf = CF\DF
+    @test x ≈ xf
+end
+
 end # module TestQR