ZTGSEN  - reorder the generalized Schur decomposition of a
       complex matrix pair (A, B) (in terms of an unitary equiva­
       lence  trans-  formation  Q'  *  (A,  B)  *  Z), so that a
       selected cluster of eigenvalues  appears  in  the  leading
       diagonal blocks of the pair (A,B)


SYNOPSIS

       SUBROUTINE ZTGSEN( IJOB,  WANTQ, WANTZ, SELECT, N, A, LDA,
                          B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, M,
                          PL,   PR,   DIF,  WORK,  LWORK,  IWORK,
                          LIWORK, INFO )

           LOGICAL        WANTQ, WANTZ

           INTEGER        IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK,
                          LWORK, M, N

           DOUBLE         PRECISION PL, PR

           LOGICAL        SELECT( * )

           INTEGER        IWORK( * )

           DOUBLE         PRECISION DIF( * )

           COMPLEX*16     A(  LDA,  * ), ALPHA( * ), B( LDB, * ),
                          BETA( * ), Q( LDQ, * ), WORK( *  ),  Z(
                          LDZ, * )


PURPOSE

       ZTGSEN  reorders  the generalized Schur decomposition of a
       complex matrix pair (A, B) (in terms of an unitary equiva­
       lence  trans-  formation  Q'  *  (A,  B)  *  Z), so that a
       selected cluster of eigenvalues  appears  in  the  leading
       diagonal  blocks of the pair (A,B). The leading columns of
       Q and Z form unitary bases of the corresponding  left  and
       right eigenspaces (deflating subspaces). (A, B) must be in
       generalized Schur canonical form, that is,  A  and  B  are
       both upper triangular.

       ZTGSEN also computes the generalized eigenvalues

                w(j)= ALPHA(j) / BETA(j)

       of the reordered matrix pair (A, B).

       Optionally,  the  routine computes estimates of reciprocal
       condition numbers for eigenvalues and  eigenspaces.  These
       are   Difu[(A11,B11),   (A22,B22)]   and   Difl[(A11,B11),
       (A22,B22)], i.e.  the  separation(s)  between  the  matrix
       pairs  (A11,  B11)  and  (A22,B22)  that correspond to the
       selected cluster and the eigenvalues outside the  cluster,
       resp.,  and  norms  of  "projections"  onto left and right

       (1,1)-block.


ARGUMENTS

       IJOB    (input) integer
               Specifies  whether  condition numbers are required
               for the cluster of eigenvalues (PL and PR) or  the
               deflating subspaces (Difu and Difl):
               =0: Only reorder w.r.t. SELECT. No extras.
               =1: Reciprocal of norms of "projections" onto left
               and right eigenspaces w.r.t. the selected  cluster
               (PL  and  PR).  =2: Upper bounds on Difu and Difl.
               F-norm-based estimate
               (DIF(1:2)).
               =3: Estimate of Difu and Difl. 1-norm-based  esti­
               mate
               (DIF(1:2)).   About 5 times as expensive as IJOB =
               2.  =4: Compute PL, PR and DIF (i.e. 0,  1  and  2
               above):  Economic version to get it all.  =5: Com­
               pute PL, PR and DIF (i.e. 0, 1 and 3 above)

       WANTQ   (input) LOGICAL

       WANTZ   (input) LOGICAL

       SELECT  (input) LOGICAL array, dimension (N)
               SELECT specifies the eigenvalues in  the  selected
               cluster.  To  select an eigenvalue w(j), SELECT(j)
               must be set to

       N       (input) INTEGER
               The order of the matrices A and B. N >= 0.

       A       (input/output) COMPLEX*16 array, dimension(LDA,N)
               On entry, the upper triangular matrix A, in gener­
               alized  Schur canonical form.  On exit, A is over­
               written by the reordered matrix A.

       LDA     (input) INTEGER
               The leading dimension  of  the  array  A.  LDA  >=
               max(1,N).

       B       (input/output) COMPLEX*16 array, dimension(LDB,N)
               On entry, the upper triangular matrix B, in gener­
               alized Schur canonical form.  On exit, B is  over­
               written by the reordered matrix B.

       LDB     (input) INTEGER
               The  leading  dimension  of  the  array  B. LDB >=
               max(1,N).

               BETA    (output) COMPLEX*16 array,  dimension  (N)
               The  diagonal  elements  of A and B, respectively,
               when the pair (A,B) has been reduced  to  general­
               ized  Schur  form.  ALPHA(i)/BETA(i) i=1,...,N are
               the generalized eigenvalues.

       Q       (input/output) COMPLEX*16 array, dimension (LDQ,N)
               On  entry,  if  WANTQ  =  .TRUE.,  Q  is an N-by-N
               matrix.  On exit, Q has been postmultiplied by the
               left  unitary  transformation matrix which reorder
               (A, B); The leading M columns of Q form  orthonor­
               mal   bases   for   the  specified  pair  of  left
               eigenspaces (deflating  subspaces).   If  WANTQ  =
               .FALSE., Q is not referenced.

       LDQ     (input) INTEGER
               The  leading  dimension  of the array Q. LDQ >= 1.
               If WANTQ = .TRUE., LDQ >= N.

       Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
               On  entry,  if  WANTZ  =  .TRUE.,  Z  is an N-by-N
               matrix.  On exit, Z has been postmultiplied by the
               left  unitary  transformation matrix which reorder
               (A, B); The leading M columns of Z form  orthonor­
               mal   bases   for   the  specified  pair  of  left
               eigenspaces (deflating  subspaces).   If  WANTZ  =
               .FALSE., Z is not referenced.

       LDZ     (input) INTEGER
               The  leading  dimension  of the array Z. LDZ >= 1.
               If WANTZ = .TRUE., LDZ >= N.

       M       (output) INTEGER
               The dimension of the specified pair  of  left  and
               right eigenspaces, (deflating subspaces) 0 <= M <=
               N.

               PL, PR  (output) DOUBLE PRECISION If IJOB =  1,  4
               or  5,  PL,  PR are lower bounds on the reciprocal
               of the norm of "projections" onto left  and  right
               eigenspace  with  respect to the selected cluster.
               0 < PL, PR <= 1.  If M = 0 or M = N, PL = PR  = 1.
               If IJOB = 0, 2 or 3 PL, PR are not referenced.

       DIF     (output) DOUBLE PRECISION array, dimension (2).
               If IJOB >= 2, DIF(1:2) store the estimates of Difu
               and Difl.
               If IJOB = 2 or 4, DIF(1:2) are F-norm-based  upper
               bounds on
               Difu  and  Difl.  If  IJOB  = 3 or 5, DIF(1:2) are
               1-norm-based estimates of Difu and Difl,  computed
               using  reversed communication with ZLACON.  If M =
               0 or N, DIF(1:2) = F-norm([A, B]).  If IJOB = 0 or


       WORK    (workspace/output) COMPLEX*16 array, dimension
               (LWORK)
               IF IJOB = 0, WORK is not  referenced.   Otherwise,
               on  exit, if INFO = 0, WORK(1) returns the optimal
               LWORK.

       LWORK   (input) INTEGER
               The dimension of the array WORK. LWORK  >=   1  If
               IJOB  = 1, 2 or 4, LWORK >=  2*M*(N-M) If IJOB = 3
               or 5, LWORK >=  4*M*(N-M)

               If LWORK = -1, then a workspace query is  assumed;
               the  routine  only  calculates the optimal size of
               the WORK array, returns this value  as  the  first
               entry  of  the  WORK  array,  and no error message
               related to LWORK is issued by XERBLA.

       IWORK   (workspace/output) INTEGER, dimension (LIWORK)
               IF IJOB = 0, IWORK is not referenced.   Otherwise,
               on exit, if INFO = 0, IWORK(1) returns the optimal
               LIWORK.

       LIWORK  (input) INTEGER
               The dimension of the array IWORK. LIWORK >= 1.  If
               IJOB  =  1, 2 or 4, LIWORK >=  N+2; If IJOB = 3 or
               5, LIWORK >= MAX(N+2, 2*M*(N-M));

               If LIWORK = -1, then a workspace query is assumed;
               the  routine  only  calculates the optimal size of
               the IWORK array, returns this value as  the  first
               entry  of  the  IWORK  array, and no error message
               related to LIWORK is issued by XERBLA.

       INFO    (output) INTEGER
               =0: Successful exit.
               <0: If INFO = -i, the i-th argument had an illegal
               value.
               =1: Reordering of (A, B) failed because the trans­
               formed matrix pair (A, B) would be  too  far  from
               generalized  Schur  form; the problem is very ill-
               conditioned.   (A,  B)  may  have  been  partially
               reordered.  If requested, 0 is returned in DIF(*),
               PL and PR.


FURTHER DETAILS

       ZTGSEN first collects the selected eigenvalues by  comput­
       ing  unitary U and W that move them to the top left corner
       of (A, B). In other words, the  selected  eigenvalues  are
       the eigenvalues of (A11, B11) in

                     U'*(A, B)*W = (A11 A12) (B11 B12) n1
                                   ( 0  A22),( 0  B22) n2


       where N = n1+n2 and U' means the conjugate transpose of U.
       The first n1 columns of U and W span the specified pair of
       left  and  right  eigenspaces (deflating subspaces) of (A,
       B).

       If (A, B) has been  obtained  from  the  generalized  real
       Schur decomposition of a matrix pair (C, D) = Q*(A, B)*Z',
       then the reordered generalized Schur form  of  (C,  D)  is
       given by

                (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',

       and  the  first  n1 columns of Q*U and Z*W span the corre­
       sponding deflating subspaces of (C, D) (Q and Z store  Q*U
       and Z*W, resp.).

       Note  that if the selected eigenvalue is sufficiently ill-
       conditioned, then its value may differ significantly  from
       its value before reordering.

       The  reciprocal  condition  numbers  of the left and right
       eigenspaces spanned by the first n1 columns of U and W (or
       Q*U and Z*W) may be returned in DIF(1:2), corresponding to
       Difu and Difl, resp.

       The Difu and Difl are defined as:

            Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
       and
            Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11,
       B11)],

       where  sigma-min(Zu) is the smallest singular value of the
       (2*n1*n2)-by-(2*n1*n2) matrix

            Zu = [ kron(In2, A11)  -kron(A22', In1) ]
                 [ kron(In2, B11)  -kron(B22', In1) ].

       Here, Inx is the identity matrix of size nx  and  A22'  is
       the  transpose of A22. kron(X, Y) is the Kronecker product
       between the matrices X and Y.

       When DIF(2) is small, small changes in (A,  B)  can  cause
       large  changes  in  the deflating subspace. An approximate
       (asymptotic) bound on the maximum  angular  error  in  the
       computed deflating subspaces is

            EPS * norm((A, B)) / DIF(2),

       where EPS is the machine precision.

       The  reciprocal  norm  of  the  projectors on the left and

       returned  in  PL  and  PR.   They are computed as follows.
       First we compute L and R so that P*(A, B)*Q is block diag­
       onal, where

            P = ( I -L ) n1           Q = ( I R ) n1
                ( 0  I ) n2    and        ( 0 I ) n2
                  n1 n2                    n1 n2

       and  (L,  R)  is the solution to the generalized Sylvester
       equation

            A11*R - L*A22 = -A12
            B11*R - L*B22 = -B12

       Then  PL  =  (F-norm(L)**2+1)**(-1/2)   and   PR   =   (F-
       norm(R)**2+1)**(-1/2).   An approximate (asymptotic) bound
       on the average absolute error of the selected  eigenvalues
       is

            EPS * norm((A, B)) / PL.

       There are also global error bounds which valid for pertur­
       bations up to a certain restriction:  A lower bound (x) on
       the  smallest F-norm(E,F) for which an eigenvalue of (A11,
       B11) may move and coalesce with  an  eigenvalue  of  (A22,
       B22) under perturbation (E,F), (i.e. (A + E, B + F), is

        x                                                       =
       min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).

       An  approximate  bound on x can be computed from DIF(1:2),
       PL and PR.

       If y = ( F-norm(E,F) / x) <= 1,  the  angles  between  the
       perturbed  (L',  R') and unperturbed (L, R) left and right
       deflating subspaces associated with the  selected  cluster
       in the (1,1)-blocks can be bounded as

        max-angle(L,  L')  <= arctan( y * PL / (1 - y * (1 - PL *
       PL)**(1/2))
        max-angle(R, R') <= arctan( y * PR / (1 - y * (1 -  PR  *
       PR)**(1/2))

       See LAPACK User's Guide section 4.11 or the following ref­
       erences for more information.

       Note that if the default method for computing  the  Frobe­
       nius-norm-  based estimate DIF is not wanted (see ZLATDF),
       then the parameter IDIFJB (see below)  should  be  changed
       from  3 to 4 (routine ZLATDF (IJOB = 2 will be used)). See
       ZTGSYL for more details.

       Based on contributions by

       Science,
          Umea University, S-901 87 Umea, Sweden.

       References
       ==========

       [1]  B. Kagstrom; A Direct Method for Reordering Eigenval­
       ues in the
           Generalized Real Schur Form of a Regular  Matrix  Pair
       (A, B), in
           M.S.  Moonen  et  al  (eds),  Linear Algebra for Large
       Scale and
           Real-Time Applications, Kluwer Academic Publ. 1993, pp
       195-218.

       [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with
       Specified
           Eigenvalues of a Regular Matrix Pair (A, B) and Condi­
       tion
           Estimation: Theory, Algorithms and Software, Report
           UMINF  -  94.04, Department of Computing Science, Umea
       University,
           S-901 87 Umea, Sweden, 1994. Also  as  LAPACK  Working
       Note 87.
           To appear in Numerical Algorithms, 1996.

       [3]  B.  Kagstrom  and P. Poromaa, LAPACK-Style Algorithms
       and Software
           for Solving the  Generalized  Sylvester  Equation  and
       Estimating the
           Separation  between Regular Matrix Pairs, Report UMINF
       - 93.23,
           Department  of  Computing  Science,  Umea  University,
       S-901 87 Umea,
           Sweden,  December  1993,  Revised  April 1994, Also as
       LAPACK working
           Note 75. To appear in ACM Trans.  on  Math.  Software,
       Vol 22, No 1,
           1996.


Man(1) output converted with man2html