Fix computation of condition number in trsen! (#23521)

* Fix computation of condition number in trsen!

* Add tests for trsen with JOB != 'N'

base/linalg/lapack.jl

@@ -5797,7 +5797,7 @@ for (trexc, trsen, tgsen, elty) in
         #       LOGICAL            SELECT( * )
         #       INTEGER            IWORK( * )
         #       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ), WR( * )
-        function trsen!(compq::Char, job::Char, select::StridedVector{BlasInt},
+        function trsen!(job::Char, compq::Char, select::StridedVector{BlasInt},
                         T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
             chkstride1(T, Q, select)
             n = checksquare(T)
@@ -5812,16 +5812,18 @@ for (trexc, trsen, tgsen, elty) in
             liwork = BlasInt(-1)
             info = Ref{BlasInt}()
             select = convert(Array{BlasInt}, select)
+            s = Ref{$elty}(zero($elty))
+            sep = Ref{$elty}(zero($elty))
             for i = 1:2  # first call returns lwork as work[1] and liwork as iwork[1]
                 ccall((@blasfunc($trsen), liblapack), Void,
                     (Ref{UInt8}, Ref{UInt8}, Ptr{BlasInt}, Ref{BlasInt},
                     Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
-                    Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{Void}, Ptr{Void},
+                    Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ref{$elty},
                     Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
                     Ptr{BlasInt}),
-                    compq, job, select, n,
+                    job, compq, select, n,
                     T, ldt, Q, ldq,
-                    wr, wi, m, C_NULL, C_NULL,
+                    wr, wi, m, s, sep,
                     work, lwork, iwork, liwork,
                     info)
                 chklapackerror(info[])
@@ -5832,7 +5834,7 @@ for (trexc, trsen, tgsen, elty) in
                     resize!(iwork, liwork)
                 end
             end
-            T, Q, iszero(wi) ? wr : complex.(wr, wi)
+            T, Q, iszero(wi) ? wr : complex.(wr, wi), s[], sep[]
         end
         trsen!(select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) =
             trsen!('N', 'V', select, T, Q)
@@ -5906,9 +5908,9 @@ for (trexc, trsen, tgsen, elty) in
     end
 end
 
-for (trexc, trsen, tgsen, elty) in
-    ((:ztrexc_, :ztrsen_, :ztgsen_, :Complex128),
-     (:ctrexc_, :ctrsen_, :ctgsen_, :Complex64))
+for (trexc, trsen, tgsen, elty, relty) in
+    ((:ztrexc_, :ztrsen_, :ztgsen_, :Complex128, :Float64),
+     (:ctrexc_, :ctrsen_, :ctgsen_, :Complex64, :Float32))
     @eval begin
         #      .. Scalar Arguments ..
         #      CHARACTER          COMPQ
@@ -5945,7 +5947,7 @@ for (trexc, trsen, tgsen, elty) in
         #      .. Array Arguments ..
         #      LOGICAL            SELECT( * )
         #      COMPLEX            Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
-        function trsen!(compq::Char, job::Char, select::StridedVector{BlasInt},
+        function trsen!(job::Char, compq::Char, select::StridedVector{BlasInt},
                         T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
             chkstride1(select, T, Q)
             n = checksquare(T)
@@ -5957,16 +5959,18 @@ for (trexc, trsen, tgsen, elty) in
             lwork = BlasInt(-1)
             info = Ref{BlasInt}()
             select = convert(Array{BlasInt}, select)
+            s = Ref{$relty}(zero($relty))
+            sep = Ref{$relty}(zero($relty))
             for i = 1:2  # first call returns lwork as work[1]
                 ccall((@blasfunc($trsen), liblapack), Void,
                     (Ref{UInt8}, Ref{UInt8}, Ptr{BlasInt}, Ref{BlasInt},
                     Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
-                    Ptr{$elty}, Ref{BlasInt}, Ptr{Void}, Ptr{Void},
+                    Ptr{$elty}, Ref{BlasInt}, Ref{$relty}, Ref{$relty},
                     Ptr{$elty}, Ref{BlasInt},
                     Ptr{BlasInt}),
-                    compq, job, select, n,
+                    job, compq, select, n,
                     T, ldt, Q, ldq,
-                    w, m, C_NULL, C_NULL,
+                    w, m, s, sep,
                     work, lwork,
                     info)
                 chklapackerror(info[])
@@ -5975,7 +5979,7 @@ for (trexc, trsen, tgsen, elty) in
                     resize!(work, lwork)
                 end
             end
-            T, Q, w
+            T, Q, w, s[], sep[]
         end
         trsen!(select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) =
             trsen!('N', 'V', select, T, Q)
@@ -6058,7 +6062,7 @@ and `ilst` specify the reordering of the vectors.
 trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::StridedMatrix, Q::StridedMatrix)
 
 """
-    trsen!(compq, job, select, T, Q) -> (T, Q, w)
+    trsen!(compq, job, select, T, Q) -> (T, Q, w, s, sep)
 
 Reorder the Schur factorization of a matrix and optionally finds reciprocal
 condition numbers. If `job = N`, no condition numbers are found. If `job = E`,
@@ -6069,7 +6073,9 @@ found. If `compq = V` the Schur vectors `Q` are updated. If `compq = N`
 the Schur vectors are not modified. `select` determines which
 eigenvalues are in the cluster.
 
-Returns `T`, `Q`, and reordered eigenvalues in `w`.
+Returns `T`, `Q`, reordered eigenvalues in `w`, the condition number of the
+cluster of eigenvalues `s`, and the condition number of the invariant subspace
+`sep`.
 """
 trsen!(compq::Char, job::Char, select::StridedVector{BlasInt}, T::StridedMatrix, Q::StridedMatrix)
 

base/linalg/schur.jl

@@ -142,7 +142,7 @@ ordschur(schur::Schur, select::Union{Vector{Bool},BitVector}) =
 Same as [`ordschur`](@ref) but overwrites the input arguments.
 """
 ordschur!(T::StridedMatrix{Ty}, Z::StridedMatrix{Ty}, select::Union{Vector{Bool},BitVector}) where {Ty<:BlasFloat} =
-    LinAlg.LAPACK.trsen!(convert(Vector{BlasInt}, select), T, Z)
+    LinAlg.LAPACK.trsen!(convert(Vector{BlasInt}, select), T, Z)[1:3]
 
 """
     ordschur(T::StridedMatrix, Z::StridedMatrix, select::Union{Vector{Bool},BitVector}) -> T::StridedMatrix, Z::StridedMatrix, λ::Vector

test/linalg/lapack.jl

@@ -555,22 +555,34 @@ end
     end
 end
 
-@testset "trexc" begin
+@testset "trsen" begin
     @testset for elty in (Float32, Float64, Complex64, Complex128)
-        for c in ('V', 'N')
-            A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4])
-            T,Q,d = schur(A)
-            Base.LinAlg.LAPACK.trsen!('N',c,Array{LinAlg.BlasInt}([0,1,0,0]),T,Q)
-            @test d[1] ≈ T[2,2]
-            @test d[2] ≈ T[1,1]
-            if c == 'V'
-                @test  Q*T*Q' ≈ A
+        for job in ('N', 'E', 'V', 'B')
+            for c in ('V', 'N')
+                A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4])
+                T,Q,d = schur(A)
+                s, sep = Base.LinAlg.LAPACK.trsen!(job,c,Array{LinAlg.BlasInt}([0,1,0,0]),T,Q)[4:5]
+                @test d[1] ≈ T[2,2]
+                @test d[2] ≈ T[1,1]
+                if c == 'V'
+                    @test  Q*T*Q' ≈ A
+                end
+                if job == 'N' || job == 'V'
+                    @test iszero(s)
+                else
+                    @test s ≈ 0.8080423 atol=1e-6
+                end
+                if job == 'N' || job == 'E'
+                    @test iszero(sep)
+                else
+                    @test sep ≈ 2. atol=3e-1
+                end
             end
         end
     end
 end
 
-@testset "trexc and trsen" begin
+@testset "trexc" begin
     @testset for elty in (Float32, Float64, Complex64, Complex128)
         for c in ('V', 'N')
             A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4])