DTGSEN  - 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 DTGSEN( 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

           DOUBLE         PRECISION PL, PR

           LOGICAL        SELECT( * )

           INTEGER        IWORK( * )

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


PURPOSE

       DTGSEN  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
       DGGES), i.e. A is block upper triangular with  1-by-1  and
       2-by-2 diagonal blocks. B is upper triangular.

       DTGSEN also computes the generalized eigenvalues

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

       of the reordered matrix pair (A, B).

       Optionally,  DTGSEN  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

       resp.,  and  norms  of  "projections"  onto left and right
       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) DOUBLE PRECISION array, dimen­
               sion(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) DOUBLE PRECISION array, dimen­
               sion(LDB,N)
               On entry, the upper triangular matrix B, with  (A,

               exit, B is overwritten by the reordered matrix  B.

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

       ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
               ALPHAI  (output) DOUBLE PRECISION array, dimension
               (N)   BETA     (output)  DOUBLE  PRECISION  array,
               dimension    (N)    On    exit,    (ALPHAR(j)    +
               ALPHAI(j)*i)/BETA(j),  j=1,...,N, will be the gen­
               eralized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i and
               BETA(j),j=1,...,N   are  the diagonals of the com­
               plex 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 trian­
               gular  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  conjugate  pair,  with
               ALPHAI(j+1) negative.

       Q       (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION 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


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

       WORK     (workspace/output)  DOUBLE  PRECISION  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.

               =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

       DTGSEN 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

                 [ 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
       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


        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  DLATDF),
       then  the  parameter  IDIFJB (see below) should be changed
       from 3 to 4 (routine DLATDF (IJOB = 2 will be used)).  See
       DTGSYL 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

       Vol 22, No 1,
           1996.


Man(1) output converted with man2html