RFC: add ordschur function to base/linalg

base/exports.jl

@@ -666,6 +666,8 @@ export
     lyap,
     norm,
     null,
+    ordschur!,
+    ordschur,
     peakflops,
     pinv,
     qr,

base/linalg.jl

@@ -85,6 +85,8 @@ export
     lyap,
     norm,
     null,
+    ordschur!,
+    ordschur,
     peakflops,
     pinv,
     qr,

base/linalg/factorization.jl

@@ -659,6 +659,11 @@ function schur(A::AbstractMatrix)
     SchurF[:T], SchurF[:Z], SchurF[:values]
 end
 
+ordschur!{Ty<:BlasFloat}(Q::StridedMatrix{Ty}, T::StridedMatrix{Ty}, select::Array{Int}) = Schur(LinAlg.LAPACK.trsen!(select, T , Q)...)
+ordschur{Ty<:BlasFloat}(Q::StridedMatrix{Ty}, T::StridedMatrix{Ty}, select::Array{Int}) = ordschur!(copy(Q), copy(T), select)
+ordschur!{Ty<:BlasFloat}(schur::Schur{Ty}, select::Array{Int}) = (res=ordschur!(schur.Z, schur.T, select); schur[:values][:]=res[:values]; res)
+ordschur{Ty<:BlasFloat}(schur::Schur{Ty}, select::Array{Int}) = ordschur(schur.Z, schur.T, select)
+
 immutable GeneralizedSchur{Ty<:BlasFloat} <: Factorization{Ty}
     S::Matrix{Ty}
     T::Matrix{Ty}

base/linalg/lapack.jl

@@ -3508,6 +3508,104 @@ for (gees, gges, elty, relty) in
         end
     end
 end
+# Reorder Schur forms
+for (trsen, elty) in
+    ((:dtrsen_,:Float64),
+     (:strsen_,:Float32))
+    @eval begin
+        function trsen!(select::Array{Int}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
+# *     .. Scalar Arguments ..
+#       CHARACTER          COMPQ, JOB
+#       INTEGER            INFO, LDQ, LDT, LIWORK, LWORK, M, N
+#       DOUBLE PRECISION   S, SEP
+# *     ..
+# *     .. Array Arguments ..
+#       LOGICAL            SELECT( * )
+#       INTEGER            IWORK( * )
+#       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ), WR( * )
+            chkstride1(T, Q)
+            n = chksquare(T)
+            ld = max(1, n)
+            wr = similar(T, $elty, n)
+            wi = similar(T, $elty, n)
+            m = sum(select)
+            work = Array($elty, 1)
+            lwork = blas_int(-1)
+            iwork = Array(BlasInt, 1)
+            liwork = blas_int(-1)
+            info = Array(BlasInt, 1)
+            select = convert(Array{BlasInt}, select)
+
+            for i = 1:2
+                ccall(($(string(trsen)), liblapack), Void,
+                    (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt},
+                    Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
+                    Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{Void}, Ptr{Void},
+                    Ptr{$elty}, Ptr  {BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
+                    Ptr{BlasInt}),
+                    &'N', &'V', select, &n,
+                    T, &ld, Q, &ld,
+                    wr, wi, &m, C_NULL, C_NULL,
+                    work, &lwork, iwork, &liwork,
+                    info)
+                @lapackerror
+                if i == 1 # only estimated optimal lwork, liwork
+                    lwork  = blas_int(real(work[1]))
+                    liwork = blas_int(real(iwork[1]))
+                    work   = Array($elty, lwork)
+                    iwork  = Array(BlasInt, liwork)
+                end
+            end
+            T, Q, all(wi .== 0) ? wr : complex(wr, wi)
+        end
+    end
+end
+
+for (trsen, elty) in
+    ((:ztrsen_,:Complex128),
+     (:ctrsen_,:Complex64))
+    @eval begin
+        function trsen!(select::Array{Int}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
+# *     .. Scalar Arguments ..
+#       CHARACTER          COMPQ, JOB
+#       INTEGER            INFO, LDQ, LDT, LWORK, M, N
+#       DOUBLE PRECISION   S, SEP
+# *     ..
+# *     .. Array Arguments ..
+#       LOGICAL            SELECT( * )
+#       COMPLEX            Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
+            chkstride1(T, Q)
+            n = chksquare(T)
+            ld = max(1, n)
+            w = similar(T, $elty, n)
+            m = sum(select)
+            work = Array($elty, 1)
+            lwork = blas_int(-1)
+            info = Array(BlasInt, 1)
+            select = convert(Array{BlasInt}, select)
+
+            for i = 1:2
+                ccall(($(string(trsen)), liblapack), Void,
+                    (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt},
+                    Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
+                    Ptr{$elty}, Ptr{BlasInt}, Ptr{Void}, Ptr{Void},
+                    Ptr{$elty}, Ptr  {BlasInt},
+                    Ptr{BlasInt}),
+                    &'N', &'V', select, &n,
+                    T, &ld, Q, &ld,
+                    w, &m, C_NULL, C_NULL,
+                    work, &lwork,
+                    info)
+                @lapackerror
+                if i == 1 # only estimated optimal lwork, liwork
+                    lwork  = blas_int(real(work[1]))
+                    work   = Array($elty, lwork)
+                end
+            end
+            T, Q, w
+        end
+    end
+end
 
 ### Rectangular full packed format
 

doc/stdlib/linalg.rst

@@ -288,12 +288,28 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
 .. function:: schurfact!(A)
 
-   Computer the Schur factorization of ``A``, overwriting ``A`` in the process. See :func:`schurfact`
+   Computes the Schur factorization of ``A``, overwriting ``A`` in the process. See :func:`schurfact`
 
 .. function:: schur(A) -> Schur[:T], Schur[:Z], Schur[:values]
 
    See :func:`schurfact`
 
+.. function:: ordschur(Q, T, select) -> Schur
+
+   Reorders the Schur factorization of a real matrix ``A=Q*T*Q'`` according to the logical array ``select`` returning a Schur object ``F``. The selected eigenvalues appear in the leading diagonal of ``F[:Schur]`` and the the corresponding leading columns of ``F[:vectors]`` form an orthonormal basis of the corresponding right invariant subspace. A complex conjugate pair of eigenvalues must be either both included or excluded via ``select``. 
+
+.. function:: ordschur!(Q, T, select) -> Schur
+
+   Reorders the Schur factorization of a real matrix ``A=Q*T*Q'``, overwriting ``Q`` and ``T`` in the process. See :func:`ordschur`
+
+.. function:: ordschur(S, select) -> Schur
+
+   Reorders the Schur factorization ``S`` of type ``Schur``.
+
+.. function:: ordschur!(S, select) -> Schur
+
+   Reorders the Schur factorization ``S`` of type ``Schur``, overwriting ``S`` in the process. See :func:`ordschur`
+
 .. function:: schurfact(A, B) -> GeneralizedSchur
 
    Computes the Generalized Schur (or QZ) factorization of the matrices ``A`` and ``B``. The (quasi) triangular Schur factors can be obtained from the ``Schur`` object ``F`` with ``F[:S]`` and ``F[:T]``, the left unitary/orthogonal Schur vectors can be obtained with ``F[:left]`` or ``F[:Q]`` and the right unitary/orthogonal Schur vectors can be obtained with ``F[:right]`` or ``F[:Z]`` such that ``A=F[:left]*F[:S]*F[:right]'`` and ``B=F[:left]*F[:T]*F[:right]'``. The generalized eigenvalues of ``A`` and ``B`` can be obtained with ``F[:alpha]./F[:beta]``.

test/linalg1.jl

@@ -180,6 +180,18 @@ debug && println("Schur")
         @test istriu(f[:Schur]) || iseltype(a,Real)
     end
 
+debug && println("Reorder Schur")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        # use asym for real schur to enforce tridiag structure
+        # avoiding partly selection of conj. eigenvalues
+        ordschura = eltya <: Complex ? a : asym
+        S = schurfact(ordschura)
+        select = rand(range(0,2), n)
+        O = ordschur(S, select)
+        bool(sum(select)) && @test_approx_eq S[:values][find(select)] O[:values][1:sum(select)]
+        @test_approx_eq O[:vectors]*O[:Schur]*O[:vectors]' ordschura
+    end
+
 debug && println("Generalized Schur")
     if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
         f = schurfact(a[1:5,1:5], a[6:10,6:10])