STGSEN  - reorder the generalized real Schur decomposition
       of a real matrix pair (A, B) (in terms of  an  orthonormal
       equivalence  trans-  formation Q' * (A, B) * Z), so that a
       selected cluster of eigenvalues  appears  in  the  leading
       diagonal blocks of the upper quasi-triangular matrix A and
       the upper triangular B


SYNOPSIS

       SUBROUTINE STGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A,  LDA,
                          B,  LDB,  ALPHAR, ALPHAI, 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

           REAL           PL, PR

           LOGICAL        SELECT( * )

           INTEGER        IWORK( * )

           REAL           A( LDA, * ), ALPHAI( * ), ALPHAR( *  ),
                          B(  LDB,  *  ), BETA( * ), DIF( * ), Q(
                          LDQ, * ), WORK( * ), Z( LDZ, * )


PURPOSE

       STGSEN reorders the generalized real  Schur  decomposition
       of  a  real matrix pair (A, B) (in terms of an orthonormal
       equivalence trans- formation Q' * (A, B) * Z), so  that  a
       selected  cluster  of  eigenvalues  appears in the leading
       diagonal blocks of the upper quasi-triangular matrix A and
       the  upper  triangular  B.  The leading columns of Q and Z
       form orthonormal bases of the corresponding left and right
       eigen-  spaces  (deflating  subspaces).  (A, B) must be in
       generalized real Schur  canonical  form  (as  returned  by
       SGGES),  i.e.  A is block upper triangular with 1-by-1 and
       2-by-2 diagonal blocks. B is upper triangular.

       STGSEN also computes the generalized eigenvalues

                   w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)

       of the reordered matrix pair (A, B).

       Optionally, STGSEN computes the  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,

       eigenspaces   w.r.t.    the   selected   cluster   in  the
       (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  a  real  eigenvalue   w(j),
               SELECT(j)  must  be set to w(j) and w(j+1), corre­
               sponding  to  a  2-by-2  diagonal  block,   either
               SELECT(j)  or  SELECT(j+1)  or both must be set to
               either  both  included  in  the  cluster  or  both
               excluded.

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

       A       (input/output) REAL array, dimension(LDA,N)
               On  entry,  the  upper  quasi-triangular matrix A,
               with (A, B) in generalized  real  Schur  canonical
               form.   On exit, A is overwritten by the reordered
               matrix A.

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

       B       (input/output) REAL array, dimension(LDB,N)
               On  entry, the upper triangular matrix B, with (A,
               B) in generalized real Schur canonical  form.   On
               exit,  B is overwritten by the reordered matrix B.

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

       ALPHAR  (output) REAL array, dimension (N)
               ALPHAI   (output)  REAL  array, dimension (N) BETA
               (output)  REAL  array,  dimension  (N)  On   exit,
               (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
               be  the  generalized  eigenvalues.   ALPHAR(j)   +
               ALPHAI(j)*i  and BETA(j),j=1,...,N  are the diago­
               nals of the complex Schur form  (S,T)  that  would
               result  if  the 2-by-2 diagonal blocks of the real
               generalized  Schur  form  of  (A,B)  were  further
               reduced  to  triangular form using complex unitary
               transformations.  If ALPHAI(j) is zero,  then  the
               j-th  eigenvalue is real; if positive, then the j-
               th and (j+1)-st eigenvalues are a  complex  conju­
               gate pair, with ALPHAI(j+1) negative.

       Q       (input/output) REAL 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   orthogonal   transformation   matrix  which
               reorder (A, B); The leading M columns  of  Q  form
               orthonormal  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;
               and if WANTQ = .TRUE., LDQ >= N.

       Z       (input/output) REAL 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  orthogonal   transformation   matrix   which
               reorder  (A,  B);  The leading M columns of Z form
               orthonormal 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 eigen- spaces (deflating subspaces). 0 <=  M
               <= N.

               PL,  PR  (output) REAL If IJOB = 1, 4 or 5, PL, PR
               are lower bounds on the reciprocal of the norm  of
               "projections" onto left and right eigenspaces with
               respect to the selected cluster.  0 < PL, PR <= 1.

               or 3, PL and PR are not referenced.

       DIF     (output) REAL 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.  If M = 0
               or N, DIF(1:2) = F-norm([A, B]).  If IJOB =  0  or
               1, DIF is not referenced.

       WORK    (workspace/output) REAL 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 >=  4*N+16.
               If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16,  2*M*(N-
               M)).   If  IJOB  =  3  or  5, LWORK >= MAX(4*N+16,
               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 array, 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+6.  If IJOB = 3  or
               5, LIWORK >= MAX(2*M*(N-M), N+6).

               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-

               reordered.  If requested, 0 is returned in DIF(*),
               PL and PR.


FURTHER DETAILS

       STGSEN  first collects the selected eigenvalues by comput­
       ing orthogonal U and W that move them to the top left cor­
       ner  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
                                     n1  n2    n1  n2

       where N = n1+n2 and U' means the 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 real 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.

       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
       right  eigenspaces  associated  with  (A11,  B11)  may  be
       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))

       erences for more information.

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

       Based on contributions by
          Bo Kagstrom and Peter Poromaa, Department of  Computing
       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