trsen.jl

# HEAD FROM lapack.jl
const liblapack = Base.liblapack_name

import Base.LinAlg: BlasFloat, BlasChar, BlasInt, blas_int, chkstride1, chksquare, Schur

# extract to appear in lapack.jl
# 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)

            for i = 1:2
                ccall(($(string(trsen)), liblapack), Void,
                    (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{Int}, 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)

            for i = 1:2
                ccall(($(string(trsen)), liblapack), Void,
                    (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{Int}, 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

ordschur!{Ty<:BlasFloat}(Q::StridedMatrix{Ty}, T::StridedMatrix{Ty}, select::Array{Int}) = Schur(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)

function test()
    n = 10
    areal = randn(n,n)/2
    aimg  = randn(n,n)/2

    for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
        asym = a'+a                  # symmetric indefinite

        ordschura = [1 0; 0 2]
        O=ordschur(schurfact(ordschura), [0 1])
        @test O[:values][1]==2

        if eltya != 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
    end
end

test()

here is the documentation to lapack's dtrsen: http://www.netlib.org/lapack/lapack-3.1.1/html/dtrsen.f.html

According to the documentation lapack's select is a LOGICAL array. I am not sure whether I get the right representation via convert(Array{Bool}, select) :(

now working as expected =)

hi dear Alex,
i am trying to use dtrsen_ as well, but i could not get proper answer, i meas there is not any happen to my output may i ask you please take a look at my code thank you so much.

Eigen::MatrixXd U(8,8) ; //unitary matrix from schur decomposition
Eigen::MatrixXd T(8,8) ;// upper triangular matrix form schur decomposition
U << -0.00553504 , 0.944387 , -0.328788, -0.000122479 , 1.27028e-05, 7.53413e-05, -0.000122381, -3.1325e-06,
0.999984 , 0.0054448 ,-0.00119516 , -4.99167e-07 ,-0.000142322 , 2.38504e-05, -7.77209e-08, -4.64215e-07,
2.02526e-06 ,-0.000794243 , -0.001715 , -0.499969 , 0.0762346 , 0.452152 , -0.734452 , -0.0188134,
-0.000455701, 0.164384 , 0.472173 , -0.00188079 , -0.854124 , 0.143147 ,-0.000459009 , -0.00278575,
-0.00050073 , 0.284744 , 0.817886 ,-0.00325787 , 0.493095 ,-0.0826293 ,0.000271747 , 0.0016084,
-1.91054e-06, 0.00110319 , 0.00306559 , 0.866035 , 0.0440108 , 0.261031 , -0.424006 , -0.0108483,
8.72962e-09 ,-3.63513e-06 ,-1.04414e-05, 4.15912e-08 , 0.139898 , 0.836721, 0.529304, 0.0128374,
3.73269e-11 ,-1.62909e-08 ,-3.8368e-08 ,-1.10561e-05 ,-0.00305771 , 0.00112896 , -0.025222 , 0.999677;

T <<  -146.514  ,   -62510.2    ,  -366637     ,  -9.75076e+07 ,-4.95515e+06 , -2.93897e+07,  4.77396e+07,  1.20219e+06,
       0.200156 ,    -366.366   ,  -2062.17    ,  -549728      ,-27944.4     , -165743     ,  269228     ,  6745.65,
       0        ,    0          ,   248.993    ,   62510.4     , 3152.96     ,  18698.8    , -30371.4    , -862.689,
       0        ,    0          ,  -0.00773466 ,   263.887     , 13.5044     ,  80.1043    , -130.128    , -2.88704,
       0        ,    0          ,   0          ,   0           , 42.4146     ,  257.717    , -425.42     ,  281.112,
       0        ,    0          ,   0          ,   0           ,-9.04028     , -43.1354    ,  61.2901    ,  382.851,
       0        ,    0          ,   0          ,   0           , 0           ,  0          , -6.37213    ,  264.696,
       0        ,    0          ,   0          ,   0           , 0           ,  0          , -2.06069    ,  7.09297;
std::cout << "U_before\n" <<U << std::endl;
std::cout << "T_before\n" <<T << std::endl;

 char job,compq;
int  ldq ,ldt, N=8, M,info , order=N;



job = 'B';
compq = 'V';
int select[N]; 
for(int i =0 ; i < N; i++)
{select[i]=false;}


ldt=N;
ldq =N;
Eigen::VectorXd wr(N);
Eigen::VectorXd wi(N);
int lwork=ldq*8;
Eigen::VectorXi iwork(lwork);
int liwork = std::max(1,M*(N-M));
 Eigen::VectorXd work(lwork);
  double s;                    
double sep; 

LAPACK_dtrsen(&job,&compq,select,&order,U.data(),&ldt,T.data(),&ldq,wr.data(),wi.data(),&M,&s,&sep,work.data(),&lwork,iwork.data(),&liwork,&info);