SRC/dtgsna.f

   1 *> \brief \b DTGSNA
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DTGSNA + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsna.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsna.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsna.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
  22 *                          LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
  23 *                          IWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          HOWMNY, JOB
  27 *       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       LOGICAL            SELECT( * )
  31 *       INTEGER            IWORK( * )
  32 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
  33 *      $                   VL( LDVL, * ), VR( LDVR, * ), WORK( * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> DTGSNA estimates reciprocal condition numbers for specified
  43 *> eigenvalues and/or eigenvectors of a matrix pair (A, B) in
  44 *> generalized real Schur canonical form (or of any matrix pair
  45 *> (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
  46 *> Z**T denotes the transpose of Z.
  47 *>
  48 *> (A, B) must be in generalized real Schur form (as returned by DGGES),
  49 *> i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
  50 *> blocks. B is upper triangular.
  51 *>
  52 *> \endverbatim
  53 *
  54 *  Arguments:
  55 *  ==========
  56 *
  57 *> \param[in] JOB
  58 *> \verbatim
  59 *>          JOB is CHARACTER*1
  60 *>          Specifies whether condition numbers are required for
  61 *>          eigenvalues (S) or eigenvectors (DIF):
  62 *>          = 'E': for eigenvalues only (S);
  63 *>          = 'V': for eigenvectors only (DIF);
  64 *>          = 'B': for both eigenvalues and eigenvectors (S and DIF).
  65 *> \endverbatim
  66 *>
  67 *> \param[in] HOWMNY
  68 *> \verbatim
  69 *>          HOWMNY is CHARACTER*1
  70 *>          = 'A': compute condition numbers for all eigenpairs;
  71 *>          = 'S': compute condition numbers for selected eigenpairs
  72 *>                 specified by the array SELECT.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] SELECT
  76 *> \verbatim
  77 *>          SELECT is LOGICAL array, dimension (N)
  78 *>          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
  79 *>          condition numbers are required. To select condition numbers
  80 *>          for the eigenpair corresponding to a real eigenvalue w(j),
  81 *>          SELECT(j) must be set to .TRUE.. To select condition numbers
  82 *>          corresponding to a complex conjugate pair of eigenvalues w(j)
  83 *>          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
  84 *>          set to .TRUE..
  85 *>          If HOWMNY = 'A', SELECT is not referenced.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] N
  89 *> \verbatim
  90 *>          N is INTEGER
  91 *>          The order of the square matrix pair (A, B). N >= 0.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] A
  95 *> \verbatim
  96 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  97 *>          The upper quasi-triangular matrix A in the pair (A,B).
  98 *> \endverbatim
  99 *>
 100 *> \param[in] LDA
 101 *> \verbatim
 102 *>          LDA is INTEGER
 103 *>          The leading dimension of the array A. LDA >= max(1,N).
 104 *> \endverbatim
 105 *>
 106 *> \param[in] B
 107 *> \verbatim
 108 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
 109 *>          The upper triangular matrix B in the pair (A,B).
 110 *> \endverbatim
 111 *>
 112 *> \param[in] LDB
 113 *> \verbatim
 114 *>          LDB is INTEGER
 115 *>          The leading dimension of the array B. LDB >= max(1,N).
 116 *> \endverbatim
 117 *>
 118 *> \param[in] VL
 119 *> \verbatim
 120 *>          VL is DOUBLE PRECISION array, dimension (LDVL,M)
 121 *>          If JOB = 'E' or 'B', VL must contain left eigenvectors of
 122 *>          (A, B), corresponding to the eigenpairs specified by HOWMNY
 123 *>          and SELECT. The eigenvectors must be stored in consecutive
 124 *>          columns of VL, as returned by DTGEVC.
 125 *>          If JOB = 'V', VL is not referenced.
 126 *> \endverbatim
 127 *>
 128 *> \param[in] LDVL
 129 *> \verbatim
 130 *>          LDVL is INTEGER
 131 *>          The leading dimension of the array VL. LDVL >= 1.
 132 *>          If JOB = 'E' or 'B', LDVL >= N.
 133 *> \endverbatim
 134 *>
 135 *> \param[in] VR
 136 *> \verbatim
 137 *>          VR is DOUBLE PRECISION array, dimension (LDVR,M)
 138 *>          If JOB = 'E' or 'B', VR must contain right eigenvectors of
 139 *>          (A, B), corresponding to the eigenpairs specified by HOWMNY
 140 *>          and SELECT. The eigenvectors must be stored in consecutive
 141 *>          columns ov VR, as returned by DTGEVC.
 142 *>          If JOB = 'V', VR is not referenced.
 143 *> \endverbatim
 144 *>
 145 *> \param[in] LDVR
 146 *> \verbatim
 147 *>          LDVR is INTEGER
 148 *>          The leading dimension of the array VR. LDVR >= 1.
 149 *>          If JOB = 'E' or 'B', LDVR >= N.
 150 *> \endverbatim
 151 *>
 152 *> \param[out] S
 153 *> \verbatim
 154 *>          S is DOUBLE PRECISION array, dimension (MM)
 155 *>          If JOB = 'E' or 'B', the reciprocal condition numbers of the
 156 *>          selected eigenvalues, stored in consecutive elements of the
 157 *>          array. For a complex conjugate pair of eigenvalues two
 158 *>          consecutive elements of S are set to the same value. Thus
 159 *>          S(j), DIF(j), and the j-th columns of VL and VR all
 160 *>          correspond to the same eigenpair (but not in general the
 161 *>          j-th eigenpair, unless all eigenpairs are selected).
 162 *>          If JOB = 'V', S is not referenced.
 163 *> \endverbatim
 164 *>
 165 *> \param[out] DIF
 166 *> \verbatim
 167 *>          DIF is DOUBLE PRECISION array, dimension (MM)
 168 *>          If JOB = 'V' or 'B', the estimated reciprocal condition
 169 *>          numbers of the selected eigenvectors, stored in consecutive
 170 *>          elements of the array. For a complex eigenvector two
 171 *>          consecutive elements of DIF are set to the same value. If
 172 *>          the eigenvalues cannot be reordered to compute DIF(j), DIF(j)
 173 *>          is set to 0; this can only occur when the true value would be
 174 *>          very small anyway.
 175 *>          If JOB = 'E', DIF is not referenced.
 176 *> \endverbatim
 177 *>
 178 *> \param[in] MM
 179 *> \verbatim
 180 *>          MM is INTEGER
 181 *>          The number of elements in the arrays S and DIF. MM >= M.
 182 *> \endverbatim
 183 *>
 184 *> \param[out] M
 185 *> \verbatim
 186 *>          M is INTEGER
 187 *>          The number of elements of the arrays S and DIF used to store
 188 *>          the specified condition numbers; for each selected real
 189 *>          eigenvalue one element is used, and for each selected complex
 190 *>          conjugate pair of eigenvalues, two elements are used.
 191 *>          If HOWMNY = 'A', M is set to N.
 192 *> \endverbatim
 193 *>
 194 *> \param[out] WORK
 195 *> \verbatim
 196 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 197 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 198 *> \endverbatim
 199 *>
 200 *> \param[in] LWORK
 201 *> \verbatim
 202 *>          LWORK is INTEGER
 203 *>          The dimension of the array WORK. LWORK >= max(1,N).
 204 *>          If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16.
 205 *>
 206 *>          If LWORK = -1, then a workspace query is assumed; the routine
 207 *>          only calculates the optimal size of the WORK array, returns
 208 *>          this value as the first entry of the WORK array, and no error
 209 *>          message related to LWORK is issued by XERBLA.
 210 *> \endverbatim
 211 *>
 212 *> \param[out] IWORK
 213 *> \verbatim
 214 *>          IWORK is INTEGER array, dimension (N + 6)
 215 *>          If JOB = 'E', IWORK is not referenced.
 216 *> \endverbatim
 217 *>
 218 *> \param[out] INFO
 219 *> \verbatim
 220 *>          INFO is INTEGER
 221 *>          =0: Successful exit
 222 *>          <0: If INFO = -i, the i-th argument had an illegal value
 223 *> \endverbatim
 224 *
 225 *  Authors:
 226 *  ========
 227 *
 228 *> \author Univ. of Tennessee
 229 *> \author Univ. of California Berkeley
 230 *> \author Univ. of Colorado Denver
 231 *> \author NAG Ltd.
 232 *
 233 *> \date November 2011
 234 *
 235 *> \ingroup doubleOTHERcomputational
 236 *
 237 *> \par Further Details:
 238 *  =====================
 239 *>
 240 *> \verbatim
 241 *>
 242 *>  The reciprocal of the condition number of a generalized eigenvalue
 243 *>  w = (a, b) is defined as
 244 *>
 245 *>       S(w) = (|u**TAv|**2 + |u**TBv|**2)**(1/2) / (norm(u)*norm(v))
 246 *>
 247 *>  where u and v are the left and right eigenvectors of (A, B)
 248 *>  corresponding to w; |z| denotes the absolute value of the complex
 249 *>  number, and norm(u) denotes the 2-norm of the vector u.
 250 *>  The pair (a, b) corresponds to an eigenvalue w = a/b (= u**TAv/u**TBv)
 251 *>  of the matrix pair (A, B). If both a and b equal zero, then (A B) is
 252 *>  singular and S(I) = -1 is returned.
 253 *>
 254 *>  An approximate error bound on the chordal distance between the i-th
 255 *>  computed generalized eigenvalue w and the corresponding exact
 256 *>  eigenvalue lambda is
 257 *>
 258 *>       chord(w, lambda) <= EPS * norm(A, B) / S(I)
 259 *>
 260 *>  where EPS is the machine precision.
 261 *>
 262 *>  The reciprocal of the condition number DIF(i) of right eigenvector u
 263 *>  and left eigenvector v corresponding to the generalized eigenvalue w
 264 *>  is defined as follows:
 265 *>
 266 *>  a) If the i-th eigenvalue w = (a,b) is real
 267 *>
 268 *>     Suppose U and V are orthogonal transformations such that
 269 *>
 270 *>              U**T*(A, B)*V  = (S, T) = ( a   *  ) ( b  *  )  1
 271 *>                                        ( 0  S22 ),( 0 T22 )  n-1
 272 *>                                          1  n-1     1 n-1
 273 *>
 274 *>     Then the reciprocal condition number DIF(i) is
 275 *>
 276 *>                Difl((a, b), (S22, T22)) = sigma-min( Zl ),
 277 *>
 278 *>     where sigma-min(Zl) denotes the smallest singular value of the
 279 *>     2(n-1)-by-2(n-1) matrix
 280 *>
 281 *>         Zl = [ kron(a, In-1)  -kron(1, S22) ]
 282 *>              [ kron(b, In-1)  -kron(1, T22) ] .
 283 *>
 284 *>     Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
 285 *>     Kronecker product between the matrices X and Y.
 286 *>
 287 *>     Note that if the default method for computing DIF(i) is wanted
 288 *>     (see DLATDF), then the parameter DIFDRI (see below) should be
 289 *>     changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)).
 290 *>     See DTGSYL for more details.
 291 *>
 292 *>  b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
 293 *>
 294 *>     Suppose U and V are orthogonal transformations such that
 295 *>
 296 *>              U**T*(A, B)*V = (S, T) = ( S11  *   ) ( T11  *  )  2
 297 *>                                       ( 0    S22 ),( 0    T22) n-2
 298 *>                                         2    n-2     2    n-2
 299 *>
 300 *>     and (S11, T11) corresponds to the complex conjugate eigenvalue
 301 *>     pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
 302 *>     that
 303 *>
 304 *>       U1**T*S11*V1 = ( s11 s12 ) and U1**T*T11*V1 = ( t11 t12 )
 305 *>                      (  0  s22 )                    (  0  t22 )
 306 *>
 307 *>     where the generalized eigenvalues w = s11/t11 and
 308 *>     conjg(w) = s22/t22.
 309 *>
 310 *>     Then the reciprocal condition number DIF(i) is bounded by
 311 *>
 312 *>         min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
 313 *>
 314 *>     where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
 315 *>     Z1 is the complex 2-by-2 matrix
 316 *>
 317 *>              Z1 =  [ s11  -s22 ]
 318 *>                    [ t11  -t22 ],
 319 *>
 320 *>     This is done by computing (using real arithmetic) the
 321 *>     roots of the characteristical polynomial det(Z1**T * Z1 - lambda I),
 322 *>     where Z1**T denotes the transpose of Z1 and det(X) denotes
 323 *>     the determinant of X.
 324 *>
 325 *>     and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
 326 *>     upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
 327 *>
 328 *>              Z2 = [ kron(S11**T, In-2)  -kron(I2, S22) ]
 329 *>                   [ kron(T11**T, In-2)  -kron(I2, T22) ]
 330 *>
 331 *>     Note that if the default method for computing DIF is wanted (see
 332 *>     DLATDF), then the parameter DIFDRI (see below) should be changed
 333 *>     from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL
 334 *>     for more details.
 335 *>
 336 *>  For each eigenvalue/vector specified by SELECT, DIF stores a
 337 *>  Frobenius norm-based estimate of Difl.
 338 *>
 339 *>  An approximate error bound for the i-th computed eigenvector VL(i) or
 340 *>  VR(i) is given by
 341 *>
 342 *>             EPS * norm(A, B) / DIF(i).
 343 *>
 344 *>  See ref. [2-3] for more details and further references.
 345 *> \endverbatim
 346 *
 347 *> \par Contributors:
 348 *  ==================
 349 *>
 350 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 351 *>     Umea University, S-901 87 Umea, Sweden.
 352 *
 353 *> \par References:
 354 *  ================
 355 *>
 356 *> \verbatim
 357 *>
 358 *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
 359 *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
 360 *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
 361 *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
 362 *>
 363 *>  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
 364 *>      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
 365 *>      Estimation: Theory, Algorithms and Software,
 366 *>      Report UMINF - 94.04, Department of Computing Science, Umea
 367 *>      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
 368 *>      Note 87. To appear in Numerical Algorithms, 1996.
 369 *>
 370 *>  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
 371 *>      for Solving the Generalized Sylvester Equation and Estimating the
 372 *>      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
 373 *>      Department of Computing Science, Umea University, S-901 87 Umea,
 374 *>      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
 375 *>      Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
 376 *>      No 1, 1996.
 377 *> \endverbatim
 378 *>
 379 *  =====================================================================
 380       SUBROUTINE DTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
 381      $                   LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
 382      $                   IWORK, INFO )
 383 *
 384 *  -- LAPACK computational routine (version 3.4.0) --
 385 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 386 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 387 *     November 2011
 388 *
 389 *     .. Scalar Arguments ..
 390       CHARACTER          HOWMNY, JOB
 391       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
 392 *     ..
 393 *     .. Array Arguments ..
 394       LOGICAL            SELECT( * )
 395       INTEGER            IWORK( * )
 396       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
 397      $                   VL( LDVL, * ), VR( LDVR, * ), WORK( * )
 398 *     ..
 399 *
 400 *  =====================================================================
 401 *
 402 *     .. Parameters ..
 403       INTEGER            DIFDRI
 404       PARAMETER          ( DIFDRI = 3 )
 405       DOUBLE PRECISION   ZERO, ONE, TWO, FOUR
 406       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
 407      $                   FOUR = 4.0D+0 )
 408 *     ..
 409 *     .. Local Scalars ..
 410       LOGICAL            LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
 411       INTEGER            I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
 412       DOUBLE PRECISION   ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
 413      $                   EPS, LNRM, RNRM, ROOT1, ROOT2, SCALE, SMLNUM,
 414      $                   TMPII, TMPIR, TMPRI, TMPRR, UHAV, UHAVI, UHBV,
 415      $                   UHBVI
 416 *     ..
 417 *     .. Local Arrays ..
 418       DOUBLE PRECISION   DUMMY( 1 ), DUMMY1( 1 )
 419 *     ..
 420 *     .. External Functions ..
 421       LOGICAL            LSAME
 422       DOUBLE PRECISION   DDOT, DLAMCH, DLAPY2, DNRM2
 423       EXTERNAL           LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
 424 *     ..
 425 *     .. External Subroutines ..
 426       EXTERNAL           DGEMV, DLACPY, DLAG2, DTGEXC, DTGSYL, XERBLA
 427 *     ..
 428 *     .. Intrinsic Functions ..
 429       INTRINSIC          MAX, MIN, SQRT
 430 *     ..
 431 *     .. Executable Statements ..
 432 *
 433 *     Decode and test the input parameters
 434 *
 435       WANTBH = LSAME( JOB, 'B' )
 436       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
 437       WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
 438 *
 439       SOMCON = LSAME( HOWMNY, 'S' )
 440 *
 441       INFO = 0
 442       LQUERY = ( LWORK.EQ.-1 )
 443 *
 444       IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
 445          INFO = -1
 446       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
 447          INFO = -2
 448       ELSE IF( N.LT.0 ) THEN
 449          INFO = -4
 450       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 451          INFO = -6
 452       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 453          INFO = -8
 454       ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
 455          INFO = -10
 456       ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
 457          INFO = -12
 458       ELSE
 459 *
 460 *        Set M to the number of eigenpairs for which condition numbers
 461 *        are required, and test MM.
 462 *
 463          IF( SOMCON ) THEN
 464             M = 0
 465             PAIR = .FALSE.
 466             DO 10 K = 1, N
 467                IF( PAIR ) THEN
 468                   PAIR = .FALSE.
 469                ELSE
 470                   IF( K.LT.N ) THEN
 471                      IF( A( K+1, K ).EQ.ZERO ) THEN
 472                         IF( SELECT( K ) )
 473      $                     M = M + 1
 474                      ELSE
 475                         PAIR = .TRUE.
 476                         IF( SELECT( K ) .OR. SELECT( K+1 ) )
 477      $                     M = M + 2
 478                      END IF
 479                   ELSE
 480                      IF( SELECT( N ) )
 481      $                  M = M + 1
 482                   END IF
 483                END IF
 484    10       CONTINUE
 485          ELSE
 486             M = N
 487          END IF
 488 *
 489          IF( N.EQ.0 ) THEN
 490             LWMIN = 1
 491          ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
 492             LWMIN = 2*N*( N + 2 ) + 16
 493          ELSE
 494             LWMIN = N
 495          END IF
 496          WORK( 1 ) = LWMIN
 497 *
 498          IF( MM.LT.M ) THEN
 499             INFO = -15
 500          ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 501             INFO = -18
 502          END IF
 503       END IF
 504 *
 505       IF( INFO.NE.0 ) THEN
 506          CALL XERBLA( 'DTGSNA', -INFO )
 507          RETURN
 508       ELSE IF( LQUERY ) THEN
 509          RETURN
 510       END IF
 511 *
 512 *     Quick return if possible
 513 *
 514       IF( N.EQ.0 )
 515      $   RETURN
 516 *
 517 *     Get machine constants
 518 *
 519       EPS = DLAMCH( 'P' )
 520       SMLNUM = DLAMCH( 'S' ) / EPS
 521       KS = 0
 522       PAIR = .FALSE.
 523 *
 524       DO 20 K = 1, N
 525 *
 526 *        Determine whether A(k,k) begins a 1-by-1 or 2-by-2 block.
 527 *
 528          IF( PAIR ) THEN
 529             PAIR = .FALSE.
 530             GO TO 20
 531          ELSE
 532             IF( K.LT.N )
 533      $         PAIR = A( K+1, K ).NE.ZERO
 534          END IF
 535 *
 536 *        Determine whether condition numbers are required for the k-th
 537 *        eigenpair.
 538 *
 539          IF( SOMCON ) THEN
 540             IF( PAIR ) THEN
 541                IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
 542      $            GO TO 20
 543             ELSE
 544                IF( .NOT.SELECT( K ) )
 545      $            GO TO 20
 546             END IF
 547          END IF
 548 *
 549          KS = KS + 1
 550 *
 551          IF( WANTS ) THEN
 552 *
 553 *           Compute the reciprocal condition number of the k-th
 554 *           eigenvalue.
 555 *
 556             IF( PAIR ) THEN
 557 *
 558 *              Complex eigenvalue pair.
 559 *
 560                RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
 561      $                DNRM2( N, VR( 1, KS+1 ), 1 ) )
 562                LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
 563      $                DNRM2( N, VL( 1, KS+1 ), 1 ) )
 564                CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
 565      $                     WORK, 1 )
 566                TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
 567                TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
 568                CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS+1 ), 1,
 569      $                     ZERO, WORK, 1 )
 570                TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
 571                TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
 572                UHAV = TMPRR + TMPII
 573                UHAVI = TMPIR - TMPRI
 574                CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
 575      $                     WORK, 1 )
 576                TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
 577                TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
 578                CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS+1 ), 1,
 579      $                     ZERO, WORK, 1 )
 580                TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
 581                TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
 582                UHBV = TMPRR + TMPII
 583                UHBVI = TMPIR - TMPRI
 584                UHAV = DLAPY2( UHAV, UHAVI )
 585                UHBV = DLAPY2( UHBV, UHBVI )
 586                COND = DLAPY2( UHAV, UHBV )
 587                S( KS ) = COND / ( RNRM*LNRM )
 588                S( KS+1 ) = S( KS )
 589 *
 590             ELSE
 591 *
 592 *              Real eigenvalue.
 593 *
 594                RNRM = DNRM2( N, VR( 1, KS ), 1 )
 595                LNRM = DNRM2( N, VL( 1, KS ), 1 )
 596                CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
 597      $                     WORK, 1 )
 598                UHAV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
 599                CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
 600      $                     WORK, 1 )
 601                UHBV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
 602                COND = DLAPY2( UHAV, UHBV )
 603                IF( COND.EQ.ZERO ) THEN
 604                   S( KS ) = -ONE
 605                ELSE
 606                   S( KS ) = COND / ( RNRM*LNRM )
 607                END IF
 608             END IF
 609          END IF
 610 *
 611          IF( WANTDF ) THEN
 612             IF( N.EQ.1 ) THEN
 613                DIF( KS ) = DLAPY2( A( 1, 1 ), B( 1, 1 ) )
 614                GO TO 20
 615             END IF
 616 *
 617 *           Estimate the reciprocal condition number of the k-th
 618 *           eigenvectors.
 619             IF( PAIR ) THEN
 620 *
 621 *              Copy the  2-by 2 pencil beginning at (A(k,k), B(k, k)).
 622 *              Compute the eigenvalue(s) at position K.
 623 *
 624                WORK( 1 ) = A( K, K )
 625                WORK( 2 ) = A( K+1, K )
 626                WORK( 3 ) = A( K, K+1 )
 627                WORK( 4 ) = A( K+1, K+1 )
 628                WORK( 5 ) = B( K, K )
 629                WORK( 6 ) = B( K+1, K )
 630                WORK( 7 ) = B( K, K+1 )
 631                WORK( 8 ) = B( K+1, K+1 )
 632                CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA,
 633      $                     DUMMY1( 1 ), ALPHAR, DUMMY( 1 ), ALPHAI )
 634                ALPRQT = ONE
 635                C1 = TWO*( ALPHAR*ALPHAR+ALPHAI*ALPHAI+BETA*BETA )
 636                C2 = FOUR*BETA*BETA*ALPHAI*ALPHAI
 637                ROOT1 = C1 + SQRT( C1*C1-4.0D0*C2 )
 638                ROOT2 = C2 / ROOT1
 639                ROOT1 = ROOT1 / TWO
 640                COND = MIN( SQRT( ROOT1 ), SQRT( ROOT2 ) )
 641             END IF
 642 *
 643 *           Copy the matrix (A, B) to the array WORK and swap the
 644 *           diagonal block beginning at A(k,k) to the (1,1) position.
 645 *
 646             CALL DLACPY( 'Full', N, N, A, LDA, WORK, N )
 647             CALL DLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
 648             IFST = K
 649             ILST = 1
 650 *
 651             CALL DTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), N,
 652      $                   DUMMY, 1, DUMMY1, 1, IFST, ILST,
 653      $                   WORK( N*N*2+1 ), LWORK-2*N*N, IERR )
 654 *
 655             IF( IERR.GT.0 ) THEN
 656 *
 657 *              Ill-conditioned problem - swap rejected.
 658 *
 659                DIF( KS ) = ZERO
 660             ELSE
 661 *
 662 *              Reordering successful, solve generalized Sylvester
 663 *              equation for R and L,
 664 *                         A22 * R - L * A11 = A12
 665 *                         B22 * R - L * B11 = B12,
 666 *              and compute estimate of Difl((A11,B11), (A22, B22)).
 667 *
 668                N1 = 1
 669                IF( WORK( 2 ).NE.ZERO )
 670      $            N1 = 2
 671                N2 = N - N1
 672                IF( N2.EQ.0 ) THEN
 673                   DIF( KS ) = COND
 674                ELSE
 675                   I = N*N + 1
 676                   IZ = 2*N*N + 1
 677                   CALL DTGSYL( 'N', DIFDRI, N2, N1, WORK( N*N1+N1+1 ),
 678      $                         N, WORK, N, WORK( N1+1 ), N,
 679      $                         WORK( N*N1+N1+I ), N, WORK( I ), N,
 680      $                         WORK( N1+I ), N, SCALE, DIF( KS ),
 681      $                         WORK( IZ+1 ), LWORK-2*N*N, IWORK, IERR )
 682 *
 683                   IF( PAIR )
 684      $               DIF( KS ) = MIN( MAX( ONE, ALPRQT )*DIF( KS ),
 685      $                           COND )
 686                END IF
 687             END IF
 688             IF( PAIR )
 689      $         DIF( KS+1 ) = DIF( KS )
 690          END IF
 691          IF( PAIR )
 692      $      KS = KS + 1
 693 *
 694    20 CONTINUE
 695       WORK( 1 ) = LWMIN
 696       RETURN
 697 *
 698 *     End of DTGSNA
 699 *
 700       END