LAPACK 3.3.1
Linear Algebra PACKage

dtgsen.f

Go to the documentation of this file.
00001       SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
00002      $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL,
00003      $                   PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO )
00004 *
00005 *  -- LAPACK routine (version 3.3.1) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *  -- April 2011                                                      --
00009 *
00010 *     Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
00011 *
00012 *     .. Scalar Arguments ..
00013       LOGICAL            WANTQ, WANTZ
00014       INTEGER            IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
00015      $                   M, N
00016       DOUBLE PRECISION   PL, PR
00017 *     ..
00018 *     .. Array Arguments ..
00019       LOGICAL            SELECT( * )
00020       INTEGER            IWORK( * )
00021       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00022      $                   B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ),
00023      $                   WORK( * ), Z( LDZ, * )
00024 *     ..
00025 *
00026 *  Purpose
00027 *  =======
00028 *
00029 *  DTGSEN reorders the generalized real Schur decomposition of a real
00030 *  matrix pair (A, B) (in terms of an orthonormal equivalence trans-
00031 *  formation Q**T * (A, B) * Z), so that a selected cluster of eigenvalues
00032 *  appears in the leading diagonal blocks of the upper quasi-triangular
00033 *  matrix A and the upper triangular B. The leading columns of Q and
00034 *  Z form orthonormal bases of the corresponding left and right eigen-
00035 *  spaces (deflating subspaces). (A, B) must be in generalized real
00036 *  Schur canonical form (as returned by DGGES), i.e. A is block upper
00037 *  triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper
00038 *  triangular.
00039 *
00040 *  DTGSEN also computes the generalized eigenvalues
00041 *
00042 *              w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)
00043 *
00044 *  of the reordered matrix pair (A, B).
00045 *
00046 *  Optionally, DTGSEN computes the estimates of reciprocal condition
00047 *  numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
00048 *  (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
00049 *  between the matrix pairs (A11, B11) and (A22,B22) that correspond to
00050 *  the selected cluster and the eigenvalues outside the cluster, resp.,
00051 *  and norms of "projections" onto left and right eigenspaces w.r.t.
00052 *  the selected cluster in the (1,1)-block.
00053 *
00054 *  Arguments
00055 *  =========
00056 *
00057 *  IJOB    (input) INTEGER
00058 *          Specifies whether condition numbers are required for the
00059 *          cluster of eigenvalues (PL and PR) or the deflating subspaces
00060 *          (Difu and Difl):
00061 *           =0: Only reorder w.r.t. SELECT. No extras.
00062 *           =1: Reciprocal of norms of "projections" onto left and right
00063 *               eigenspaces w.r.t. the selected cluster (PL and PR).
00064 *           =2: Upper bounds on Difu and Difl. F-norm-based estimate
00065 *               (DIF(1:2)).
00066 *           =3: Estimate of Difu and Difl. 1-norm-based estimate
00067 *               (DIF(1:2)).
00068 *               About 5 times as expensive as IJOB = 2.
00069 *           =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic
00070 *               version to get it all.
00071 *           =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above)
00072 *
00073 *  WANTQ   (input) LOGICAL
00074 *          .TRUE. : update the left transformation matrix Q;
00075 *          .FALSE.: do not update Q.
00076 *
00077 *  WANTZ   (input) LOGICAL
00078 *          .TRUE. : update the right transformation matrix Z;
00079 *          .FALSE.: do not update Z.
00080 *
00081 *  SELECT  (input) LOGICAL array, dimension (N)
00082 *          SELECT specifies the eigenvalues in the selected cluster.
00083 *          To select a real eigenvalue w(j), SELECT(j) must be set to
00084 *          .TRUE.. To select a complex conjugate pair of eigenvalues
00085 *          w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
00086 *          either SELECT(j) or SELECT(j+1) or both must be set to
00087 *          .TRUE.; a complex conjugate pair of eigenvalues must be
00088 *          either both included in the cluster or both excluded.
00089 *
00090 *  N       (input) INTEGER
00091 *          The order of the matrices A and B. N >= 0.
00092 *
00093 *  A       (input/output) DOUBLE PRECISION array, dimension(LDA,N)
00094 *          On entry, the upper quasi-triangular matrix A, with (A, B) in
00095 *          generalized real Schur canonical form.
00096 *          On exit, A is overwritten by the reordered matrix A.
00097 *
00098 *  LDA     (input) INTEGER
00099 *          The leading dimension of the array A. LDA >= max(1,N).
00100 *
00101 *  B       (input/output) DOUBLE PRECISION array, dimension(LDB,N)
00102 *          On entry, the upper triangular matrix B, with (A, B) in
00103 *          generalized real Schur canonical form.
00104 *          On exit, B is overwritten by the reordered matrix B.
00105 *
00106 *  LDB     (input) INTEGER
00107 *          The leading dimension of the array B. LDB >= max(1,N).
00108 *
00109 *  ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
00110 *  ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
00111 *  BETA    (output) DOUBLE PRECISION array, dimension (N)
00112 *          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
00113 *          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i
00114 *          and BETA(j),j=1,...,N  are the diagonals of the complex Schur
00115 *          form (S,T) that would result if the 2-by-2 diagonal blocks of
00116 *          the real generalized Schur form of (A,B) were further reduced
00117 *          to triangular form using complex unitary transformations.
00118 *          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
00119 *          positive, then the j-th and (j+1)-st eigenvalues are a
00120 *          complex conjugate pair, with ALPHAI(j+1) negative.
00121 *
00122 *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
00123 *          On entry, if WANTQ = .TRUE., Q is an N-by-N matrix.
00124 *          On exit, Q has been postmultiplied by the left orthogonal
00125 *          transformation matrix which reorder (A, B); The leading M
00126 *          columns of Q form orthonormal bases for the specified pair of
00127 *          left eigenspaces (deflating subspaces).
00128 *          If WANTQ = .FALSE., Q is not referenced.
00129 *
00130 *  LDQ     (input) INTEGER
00131 *          The leading dimension of the array Q.  LDQ >= 1;
00132 *          and if WANTQ = .TRUE., LDQ >= N.
00133 *
00134 *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
00135 *          On entry, if WANTZ = .TRUE., Z is an N-by-N matrix.
00136 *          On exit, Z has been postmultiplied by the left orthogonal
00137 *          transformation matrix which reorder (A, B); The leading M
00138 *          columns of Z form orthonormal bases for the specified pair of
00139 *          left eigenspaces (deflating subspaces).
00140 *          If WANTZ = .FALSE., Z is not referenced.
00141 *
00142 *  LDZ     (input) INTEGER
00143 *          The leading dimension of the array Z. LDZ >= 1;
00144 *          If WANTZ = .TRUE., LDZ >= N.
00145 *
00146 *  M       (output) INTEGER
00147 *          The dimension of the specified pair of left and right eigen-
00148 *          spaces (deflating subspaces). 0 <= M <= N.
00149 *
00150 *  PL      (output) DOUBLE PRECISION
00151 *  PR      (output) DOUBLE PRECISION
00152 *          If IJOB = 1, 4 or 5, PL, PR are lower bounds on the
00153 *          reciprocal of the norm of "projections" onto left and right
00154 *          eigenspaces with respect to the selected cluster.
00155 *          0 < PL, PR <= 1.
00156 *          If M = 0 or M = N, PL = PR  = 1.
00157 *          If IJOB = 0, 2 or 3, PL and PR are not referenced.
00158 *
00159 *  DIF     (output) DOUBLE PRECISION array, dimension (2).
00160 *          If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
00161 *          If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
00162 *          Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
00163 *          estimates of Difu and Difl.
00164 *          If M = 0 or N, DIF(1:2) = F-norm([A, B]).
00165 *          If IJOB = 0 or 1, DIF is not referenced.
00166 *
00167 *  WORK    (workspace/output) DOUBLE PRECISION array,
00168 *          dimension (MAX(1,LWORK)) 
00169 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00170 *
00171 *  LWORK   (input) INTEGER
00172 *          The dimension of the array WORK. LWORK >=  4*N+16.
00173 *          If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)).
00174 *          If IJOB = 3 or 5, LWORK >= MAX(4*N+16, 4*M*(N-M)).
00175 *
00176 *          If LWORK = -1, then a workspace query is assumed; the routine
00177 *          only calculates the optimal size of the WORK array, returns
00178 *          this value as the first entry of the WORK array, and no error
00179 *          message related to LWORK is issued by XERBLA.
00180 *
00181 *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
00182 *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
00183 *
00184 *  LIWORK  (input) INTEGER
00185 *          The dimension of the array IWORK. LIWORK >= 1.
00186 *          If IJOB = 1, 2 or 4, LIWORK >=  N+6.
00187 *          If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M), N+6).
00188 *
00189 *          If LIWORK = -1, then a workspace query is assumed; the
00190 *          routine only calculates the optimal size of the IWORK array,
00191 *          returns this value as the first entry of the IWORK array, and
00192 *          no error message related to LIWORK is issued by XERBLA.
00193 *
00194 *  INFO    (output) INTEGER
00195 *            =0: Successful exit.
00196 *            <0: If INFO = -i, the i-th argument had an illegal value.
00197 *            =1: Reordering of (A, B) failed because the transformed
00198 *                matrix pair (A, B) would be too far from generalized
00199 *                Schur form; the problem is very ill-conditioned.
00200 *                (A, B) may have been partially reordered.
00201 *                If requested, 0 is returned in DIF(*), PL and PR.
00202 *
00203 *  Further Details
00204 *  ===============
00205 *
00206 *  DTGSEN first collects the selected eigenvalues by computing
00207 *  orthogonal U and W that move them to the top left corner of (A, B).
00208 *  In other words, the selected eigenvalues are the eigenvalues of
00209 *  (A11, B11) in:
00210 *
00211 *              U**T*(A, B)*W = (A11 A12) (B11 B12) n1
00212 *                              ( 0  A22),( 0  B22) n2
00213 *                                n1  n2    n1  n2
00214 *
00215 *  where N = n1+n2 and U**T means the transpose of U. The first n1 columns
00216 *  of U and W span the specified pair of left and right eigenspaces
00217 *  (deflating subspaces) of (A, B).
00218 *
00219 *  If (A, B) has been obtained from the generalized real Schur
00220 *  decomposition of a matrix pair (C, D) = Q*(A, B)*Z**T, then the
00221 *  reordered generalized real Schur form of (C, D) is given by
00222 *
00223 *           (C, D) = (Q*U)*(U**T*(A, B)*W)*(Z*W)**T,
00224 *
00225 *  and the first n1 columns of Q*U and Z*W span the corresponding
00226 *  deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
00227 *
00228 *  Note that if the selected eigenvalue is sufficiently ill-conditioned,
00229 *  then its value may differ significantly from its value before
00230 *  reordering.
00231 *
00232 *  The reciprocal condition numbers of the left and right eigenspaces
00233 *  spanned by the first n1 columns of U and W (or Q*U and Z*W) may
00234 *  be returned in DIF(1:2), corresponding to Difu and Difl, resp.
00235 *
00236 *  The Difu and Difl are defined as:
00237 *
00238 *       Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
00239 *  and
00240 *       Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
00241 *
00242 *  where sigma-min(Zu) is the smallest singular value of the
00243 *  (2*n1*n2)-by-(2*n1*n2) matrix
00244 *
00245 *       Zu = [ kron(In2, A11)  -kron(A22**T, In1) ]
00246 *            [ kron(In2, B11)  -kron(B22**T, In1) ].
00247 *
00248 *  Here, Inx is the identity matrix of size nx and A22**T is the
00249 *  transpose of A22. kron(X, Y) is the Kronecker product between
00250 *  the matrices X and Y.
00251 *
00252 *  When DIF(2) is small, small changes in (A, B) can cause large changes
00253 *  in the deflating subspace. An approximate (asymptotic) bound on the
00254 *  maximum angular error in the computed deflating subspaces is
00255 *
00256 *       EPS * norm((A, B)) / DIF(2),
00257 *
00258 *  where EPS is the machine precision.
00259 *
00260 *  The reciprocal norm of the projectors on the left and right
00261 *  eigenspaces associated with (A11, B11) may be returned in PL and PR.
00262 *  They are computed as follows. First we compute L and R so that
00263 *  P*(A, B)*Q is block diagonal, where
00264 *
00265 *       P = ( I -L ) n1           Q = ( I R ) n1
00266 *           ( 0  I ) n2    and        ( 0 I ) n2
00267 *             n1 n2                    n1 n2
00268 *
00269 *  and (L, R) is the solution to the generalized Sylvester equation
00270 *
00271 *       A11*R - L*A22 = -A12
00272 *       B11*R - L*B22 = -B12
00273 *
00274 *  Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
00275 *  An approximate (asymptotic) bound on the average absolute error of
00276 *  the selected eigenvalues is
00277 *
00278 *       EPS * norm((A, B)) / PL.
00279 *
00280 *  There are also global error bounds which valid for perturbations up
00281 *  to a certain restriction:  A lower bound (x) on the smallest
00282 *  F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
00283 *  coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
00284 *  (i.e. (A + E, B + F), is
00285 *
00286 *   x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
00287 *
00288 *  An approximate bound on x can be computed from DIF(1:2), PL and PR.
00289 *
00290 *  If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
00291 *  (L', R') and unperturbed (L, R) left and right deflating subspaces
00292 *  associated with the selected cluster in the (1,1)-blocks can be
00293 *  bounded as
00294 *
00295 *   max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
00296 *   max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
00297 *
00298 *  See LAPACK User's Guide section 4.11 or the following references
00299 *  for more information.
00300 *
00301 *  Note that if the default method for computing the Frobenius-norm-
00302 *  based estimate DIF is not wanted (see DLATDF), then the parameter
00303 *  IDIFJB (see below) should be changed from 3 to 4 (routine DLATDF
00304 *  (IJOB = 2 will be used)). See DTGSYL for more details.
00305 *
00306 *  Based on contributions by
00307 *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00308 *     Umea University, S-901 87 Umea, Sweden.
00309 *
00310 *  References
00311 *  ==========
00312 *
00313 *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
00314 *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
00315 *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
00316 *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
00317 *
00318 *  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
00319 *      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
00320 *      Estimation: Theory, Algorithms and Software,
00321 *      Report UMINF - 94.04, Department of Computing Science, Umea
00322 *      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
00323 *      Note 87. To appear in Numerical Algorithms, 1996.
00324 *
00325 *  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
00326 *      for Solving the Generalized Sylvester Equation and Estimating the
00327 *      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
00328 *      Department of Computing Science, Umea University, S-901 87 Umea,
00329 *      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
00330 *      Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
00331 *      1996.
00332 *
00333 *  =====================================================================
00334 *
00335 *     .. Parameters ..
00336       INTEGER            IDIFJB
00337       PARAMETER          ( IDIFJB = 3 )
00338       DOUBLE PRECISION   ZERO, ONE
00339       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00340 *     ..
00341 *     .. Local Scalars ..
00342       LOGICAL            LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2,
00343      $                   WANTP
00344       INTEGER            I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN,
00345      $                   MN2, N1, N2
00346       DOUBLE PRECISION   DSCALE, DSUM, EPS, RDSCAL, SMLNUM
00347 *     ..
00348 *     .. Local Arrays ..
00349       INTEGER            ISAVE( 3 )
00350 *     ..
00351 *     .. External Subroutines ..
00352       EXTERNAL           DLACN2, DLACPY, DLAG2, DLASSQ, DTGEXC, DTGSYL,
00353      $                   XERBLA
00354 *     ..
00355 *     .. External Functions ..
00356       DOUBLE PRECISION   DLAMCH
00357       EXTERNAL           DLAMCH
00358 *     ..
00359 *     .. Intrinsic Functions ..
00360       INTRINSIC          MAX, SIGN, SQRT
00361 *     ..
00362 *     .. Executable Statements ..
00363 *
00364 *     Decode and test the input parameters
00365 *
00366       INFO = 0
00367       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
00368 *
00369       IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN
00370          INFO = -1
00371       ELSE IF( N.LT.0 ) THEN
00372          INFO = -5
00373       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00374          INFO = -7
00375       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00376          INFO = -9
00377       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
00378          INFO = -14
00379       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
00380          INFO = -16
00381       END IF
00382 *
00383       IF( INFO.NE.0 ) THEN
00384          CALL XERBLA( 'DTGSEN', -INFO )
00385          RETURN
00386       END IF
00387 *
00388 *     Get machine constants
00389 *
00390       EPS = DLAMCH( 'P' )
00391       SMLNUM = DLAMCH( 'S' ) / EPS
00392       IERR = 0
00393 *
00394       WANTP = IJOB.EQ.1 .OR. IJOB.GE.4
00395       WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4
00396       WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5
00397       WANTD = WANTD1 .OR. WANTD2
00398 *
00399 *     Set M to the dimension of the specified pair of deflating
00400 *     subspaces.
00401 *
00402       M = 0
00403       PAIR = .FALSE.
00404       DO 10 K = 1, N
00405          IF( PAIR ) THEN
00406             PAIR = .FALSE.
00407          ELSE
00408             IF( K.LT.N ) THEN
00409                IF( A( K+1, K ).EQ.ZERO ) THEN
00410                   IF( SELECT( K ) )
00411      $               M = M + 1
00412                ELSE
00413                   PAIR = .TRUE.
00414                   IF( SELECT( K ) .OR. SELECT( K+1 ) )
00415      $               M = M + 2
00416                END IF
00417             ELSE
00418                IF( SELECT( N ) )
00419      $            M = M + 1
00420             END IF
00421          END IF
00422    10 CONTINUE
00423 *
00424       IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
00425          LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) )
00426          LIWMIN = MAX( 1, N+6 )
00427       ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
00428          LWMIN = MAX( 1, 4*N+16, 4*M*( N-M ) )
00429          LIWMIN = MAX( 1, 2*M*( N-M ), N+6 )
00430       ELSE
00431          LWMIN = MAX( 1, 4*N+16 )
00432          LIWMIN = 1
00433       END IF
00434 *
00435       WORK( 1 ) = LWMIN
00436       IWORK( 1 ) = LIWMIN
00437 *
00438       IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00439          INFO = -22
00440       ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00441          INFO = -24
00442       END IF
00443 *
00444       IF( INFO.NE.0 ) THEN
00445          CALL XERBLA( 'DTGSEN', -INFO )
00446          RETURN
00447       ELSE IF( LQUERY ) THEN
00448          RETURN
00449       END IF
00450 *
00451 *     Quick return if possible.
00452 *
00453       IF( M.EQ.N .OR. M.EQ.0 ) THEN
00454          IF( WANTP ) THEN
00455             PL = ONE
00456             PR = ONE
00457          END IF
00458          IF( WANTD ) THEN
00459             DSCALE = ZERO
00460             DSUM = ONE
00461             DO 20 I = 1, N
00462                CALL DLASSQ( N, A( 1, I ), 1, DSCALE, DSUM )
00463                CALL DLASSQ( N, B( 1, I ), 1, DSCALE, DSUM )
00464    20       CONTINUE
00465             DIF( 1 ) = DSCALE*SQRT( DSUM )
00466             DIF( 2 ) = DIF( 1 )
00467          END IF
00468          GO TO 60
00469       END IF
00470 *
00471 *     Collect the selected blocks at the top-left corner of (A, B).
00472 *
00473       KS = 0
00474       PAIR = .FALSE.
00475       DO 30 K = 1, N
00476          IF( PAIR ) THEN
00477             PAIR = .FALSE.
00478          ELSE
00479 *
00480             SWAP = SELECT( K )
00481             IF( K.LT.N ) THEN
00482                IF( A( K+1, K ).NE.ZERO ) THEN
00483                   PAIR = .TRUE.
00484                   SWAP = SWAP .OR. SELECT( K+1 )
00485                END IF
00486             END IF
00487 *
00488             IF( SWAP ) THEN
00489                KS = KS + 1
00490 *
00491 *              Swap the K-th block to position KS.
00492 *              Perform the reordering of diagonal blocks in (A, B)
00493 *              by orthogonal transformation matrices and update
00494 *              Q and Z accordingly (if requested):
00495 *
00496                KK = K
00497                IF( K.NE.KS )
00498      $            CALL DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00499      $                         Z, LDZ, KK, KS, WORK, LWORK, IERR )
00500 *
00501                IF( IERR.GT.0 ) THEN
00502 *
00503 *                 Swap is rejected: exit.
00504 *
00505                   INFO = 1
00506                   IF( WANTP ) THEN
00507                      PL = ZERO
00508                      PR = ZERO
00509                   END IF
00510                   IF( WANTD ) THEN
00511                      DIF( 1 ) = ZERO
00512                      DIF( 2 ) = ZERO
00513                   END IF
00514                   GO TO 60
00515                END IF
00516 *
00517                IF( PAIR )
00518      $            KS = KS + 1
00519             END IF
00520          END IF
00521    30 CONTINUE
00522       IF( WANTP ) THEN
00523 *
00524 *        Solve generalized Sylvester equation for R and L
00525 *        and compute PL and PR.
00526 *
00527          N1 = M
00528          N2 = N - M
00529          I = N1 + 1
00530          IJB = 0
00531          CALL DLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
00532          CALL DLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
00533      $                N1 )
00534          CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
00535      $                N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
00536      $                DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
00537      $                LWORK-2*N1*N2, IWORK, IERR )
00538 *
00539 *        Estimate the reciprocal of norms of "projections" onto left
00540 *        and right eigenspaces.
00541 *
00542          RDSCAL = ZERO
00543          DSUM = ONE
00544          CALL DLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM )
00545          PL = RDSCAL*SQRT( DSUM )
00546          IF( PL.EQ.ZERO ) THEN
00547             PL = ONE
00548          ELSE
00549             PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) )
00550          END IF
00551          RDSCAL = ZERO
00552          DSUM = ONE
00553          CALL DLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM )
00554          PR = RDSCAL*SQRT( DSUM )
00555          IF( PR.EQ.ZERO ) THEN
00556             PR = ONE
00557          ELSE
00558             PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) )
00559          END IF
00560       END IF
00561 *
00562       IF( WANTD ) THEN
00563 *
00564 *        Compute estimates of Difu and Difl.
00565 *
00566          IF( WANTD1 ) THEN
00567             N1 = M
00568             N2 = N - M
00569             I = N1 + 1
00570             IJB = IDIFJB
00571 *
00572 *           Frobenius norm-based Difu-estimate.
00573 *
00574             CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
00575      $                   N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
00576      $                   N1, DSCALE, DIF( 1 ), WORK( 2*N1*N2+1 ),
00577      $                   LWORK-2*N1*N2, IWORK, IERR )
00578 *
00579 *           Frobenius norm-based Difl-estimate.
00580 *
00581             CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
00582      $                   N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
00583      $                   N2, DSCALE, DIF( 2 ), WORK( 2*N1*N2+1 ),
00584      $                   LWORK-2*N1*N2, IWORK, IERR )
00585          ELSE
00586 *
00587 *
00588 *           Compute 1-norm-based estimates of Difu and Difl using
00589 *           reversed communication with DLACN2. In each step a
00590 *           generalized Sylvester equation or a transposed variant
00591 *           is solved.
00592 *
00593             KASE = 0
00594             N1 = M
00595             N2 = N - M
00596             I = N1 + 1
00597             IJB = 0
00598             MN2 = 2*N1*N2
00599 *
00600 *           1-norm-based estimate of Difu.
00601 *
00602    40       CONTINUE
00603             CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 1 ),
00604      $                   KASE, ISAVE )
00605             IF( KASE.NE.0 ) THEN
00606                IF( KASE.EQ.1 ) THEN
00607 *
00608 *                 Solve generalized Sylvester equation.
00609 *
00610                   CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
00611      $                         WORK, N1, B, LDB, B( I, I ), LDB,
00612      $                         WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
00613      $                         WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
00614      $                         IERR )
00615                ELSE
00616 *
00617 *                 Solve the transposed variant.
00618 *
00619                   CALL DTGSYL( 'T', IJB, N1, N2, A, LDA, A( I, I ), LDA,
00620      $                         WORK, N1, B, LDB, B( I, I ), LDB,
00621      $                         WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
00622      $                         WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
00623      $                         IERR )
00624                END IF
00625                GO TO 40
00626             END IF
00627             DIF( 1 ) = DSCALE / DIF( 1 )
00628 *
00629 *           1-norm-based estimate of Difl.
00630 *
00631    50       CONTINUE
00632             CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 2 ),
00633      $                   KASE, ISAVE )
00634             IF( KASE.NE.0 ) THEN
00635                IF( KASE.EQ.1 ) THEN
00636 *
00637 *                 Solve generalized Sylvester equation.
00638 *
00639                   CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
00640      $                         WORK, N2, B( I, I ), LDB, B, LDB,
00641      $                         WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
00642      $                         WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
00643      $                         IERR )
00644                ELSE
00645 *
00646 *                 Solve the transposed variant.
00647 *
00648                   CALL DTGSYL( 'T', IJB, N2, N1, A( I, I ), LDA, A, LDA,
00649      $                         WORK, N2, B( I, I ), LDB, B, LDB,
00650      $                         WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
00651      $                         WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
00652      $                         IERR )
00653                END IF
00654                GO TO 50
00655             END IF
00656             DIF( 2 ) = DSCALE / DIF( 2 )
00657 *
00658          END IF
00659       END IF
00660 *
00661    60 CONTINUE
00662 *
00663 *     Compute generalized eigenvalues of reordered pair (A, B) and
00664 *     normalize the generalized Schur form.
00665 *
00666       PAIR = .FALSE.
00667       DO 80 K = 1, N
00668          IF( PAIR ) THEN
00669             PAIR = .FALSE.
00670          ELSE
00671 *
00672             IF( K.LT.N ) THEN
00673                IF( A( K+1, K ).NE.ZERO ) THEN
00674                   PAIR = .TRUE.
00675                END IF
00676             END IF
00677 *
00678             IF( PAIR ) THEN
00679 *
00680 *             Compute the eigenvalue(s) at position K.
00681 *
00682                WORK( 1 ) = A( K, K )
00683                WORK( 2 ) = A( K+1, K )
00684                WORK( 3 ) = A( K, K+1 )
00685                WORK( 4 ) = A( K+1, K+1 )
00686                WORK( 5 ) = B( K, K )
00687                WORK( 6 ) = B( K+1, K )
00688                WORK( 7 ) = B( K, K+1 )
00689                WORK( 8 ) = B( K+1, K+1 )
00690                CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA( K ),
00691      $                     BETA( K+1 ), ALPHAR( K ), ALPHAR( K+1 ),
00692      $                     ALPHAI( K ) )
00693                ALPHAI( K+1 ) = -ALPHAI( K )
00694 *
00695             ELSE
00696 *
00697                IF( SIGN( ONE, B( K, K ) ).LT.ZERO ) THEN
00698 *
00699 *                 If B(K,K) is negative, make it positive
00700 *
00701                   DO 70 I = 1, N
00702                      A( K, I ) = -A( K, I )
00703                      B( K, I ) = -B( K, I )
00704                      IF( WANTQ ) Q( I, K ) = -Q( I, K )
00705    70             CONTINUE
00706                END IF
00707 *
00708                ALPHAR( K ) = A( K, K )
00709                ALPHAI( K ) = ZERO
00710                BETA( K ) = B( K, K )
00711 *
00712             END IF
00713          END IF
00714    80 CONTINUE
00715 *
00716       WORK( 1 ) = LWMIN
00717       IWORK( 1 ) = LIWMIN
00718 *
00719       RETURN
00720 *
00721 *     End of DTGSEN
00722 *
00723       END
 All Files Functions