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>
JuliaLang · May 5, 2023 · e70354d · e70354d
1 parent 1729d40
commit e70354d
Show file tree

Hide file tree

Showing 2 changed files with 34 additions and 5 deletions.
diff --git 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
@@ -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