Add wrapper for LAPACK _trexc methods (reordering of the Schur form)

Author: mfasi

base/linalg/lapack.jl

@@ -3855,10 +3855,37 @@ for (gees, gges, elty, relty) in
     end
 end
 # Reorder Schur forms
-for (trsen, tgsen, elty) in
-    ((:dtrsen_, :dtgsen_, :Float64),
-     (:strsen_, :stgsen_, :Float32))
+for (trexc, trsen, tgsen, elty) in
+    ((:dtrexc_, :dtrsen_, :dtgsen_, :Float64),
+     (:strexc_, :strsen_, :stgsen_, :Float32))
     @eval begin
+        trexc!(ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) = trexc!('V', ifst, ilst, T, Q)
+        function trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
+# *     .. Scalar Arguments ..
+#       CHARACTER          COMPQ
+#       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
+# *     ..
+# *     .. Array Arguments ..
+#       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
+            chkstride1(T, Q)
+            n = chksquare(T)
+            ldt = max(1, stride(T, 2))
+            ldq = max(1, stride(Q, 2))
+            work = Array($elty, n)
+            info = Array(BlasInt, 1)
+
+            ccall(($(blasfunc(trexc)), liblapack), Void,
+                  (Ptr{UInt8},  Ptr{BlasInt},
+                   Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
+                   Ptr{BlasInt}, Ptr{BlasInt},
+                   Ptr{$elty}, Ptr{BlasInt}),
+                  &compq, &n,
+                  T, &ldt, Q, &ldq,
+                  &ifst, &ilst,
+                  work, info)
+            @lapackerror
+            T, Q
+        end
         function trsen!(select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
 # *     .. Scalar Arguments ..
 #       CHARACTER          COMPQ, JOB
@@ -3975,10 +4002,36 @@ for (trsen, tgsen, elty) in
     end
 end
 
-for (trsen, tgsen, elty) in
-    ((:ztrsen_, :ztgsen_, :Complex128),
-     (:ctrsen_, :ctgsen_, :Complex64))
+for (trexc, trsen, tgsen, elty) in
+    ((:ztrexc_, :ztrsen_, :ztgsen_, :Complex128),
+     (:ctrexc_, :ctrsen_, :ctgsen_, :Complex64))
     @eval begin
+        trexc!(ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) = trexc!('V', ifst, ilst, T, Q)
+        function trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
+# *     .. Scalar Arguments ..
+#       CHARACTER          COMPQ
+#       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
+# *     ..
+# *     .. Array Arguments ..
+#       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
+            chkstride1(T, Q)
+            n = chksquare(T)
+            ldt = max(1, stride(T, 2))
+            ldq = max(1, stride(Q, 2))
+            info = Array(BlasInt, 1)
+
+            ccall(($(blasfunc(trexc)), liblapack), Void,
+                  (Ptr{UInt8},  Ptr{BlasInt},
+                   Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
+                   Ptr{BlasInt}, Ptr{BlasInt},
+                   Ptr{BlasInt}),
+                  &compq, &n,
+                  T, &ldt, Q, &ldq,
+                  &ifst, &ilst,
+                  info)
+            @lapackerror
+            T, Q
+        end
         function trsen!(select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
 # *     .. Scalar Arguments ..
 #       CHARACTER          COMPQ, JOB

test/linalg/lapack.jl

@@ -427,6 +427,20 @@ for elty in (Float32, Float64, Complex64, Complex128)
     @test_throws DimensionMismatch LAPACK.laic1!(1,rand(elty,10),real(rand(elty)),rand(elty,11),rand(elty))
 end
 
+#trexc
+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.trexc!(c,1,2,T,Q)
+        @test d[1] ≈ T[2,2]
+        @test d[2] ≈ T[1,1]
+        if (c == 'V')
+            @test  Q * T * Q' ≈ A
+        end
+    end
+end
+
 # Test our own linear algebra functionaly against LAPACK
 for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
     for nn in (5,10,15)