SRC/dtrsna.f

   1 *> \brief \b DTRSNA
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DTRSNA + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsna.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsna.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsna.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
  22 *                          LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
  23 *                          INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          HOWMNY, JOB
  27 *       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       LOGICAL            SELECT( * )
  31 *       INTEGER            IWORK( * )
  32 *       DOUBLE PRECISION   S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
  33 *      $                   VR( LDVR, * ), WORK( LDWORK, * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> DTRSNA estimates reciprocal condition numbers for specified
  43 *> eigenvalues and/or right eigenvectors of a real upper
  44 *> quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
  45 *> orthogonal).
  46 *>
  47 *> T must be in Schur canonical form (as returned by DHSEQR), that is,
  48 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
  49 *> 2-by-2 diagonal block has its diagonal elements equal and its
  50 *> off-diagonal elements of opposite sign.
  51 *> \endverbatim
  52 *
  53 *  Arguments:
  54 *  ==========
  55 *
  56 *> \param[in] JOB
  57 *> \verbatim
  58 *>          JOB is CHARACTER*1
  59 *>          Specifies whether condition numbers are required for
  60 *>          eigenvalues (S) or eigenvectors (SEP):
  61 *>          = 'E': for eigenvalues only (S);
  62 *>          = 'V': for eigenvectors only (SEP);
  63 *>          = 'B': for both eigenvalues and eigenvectors (S and SEP).
  64 *> \endverbatim
  65 *>
  66 *> \param[in] HOWMNY
  67 *> \verbatim
  68 *>          HOWMNY is CHARACTER*1
  69 *>          = 'A': compute condition numbers for all eigenpairs;
  70 *>          = 'S': compute condition numbers for selected eigenpairs
  71 *>                 specified by the array SELECT.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] SELECT
  75 *> \verbatim
  76 *>          SELECT is LOGICAL array, dimension (N)
  77 *>          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
  78 *>          condition numbers are required. To select condition numbers
  79 *>          for the eigenpair corresponding to a real eigenvalue w(j),
  80 *>          SELECT(j) must be set to .TRUE.. To select condition numbers
  81 *>          corresponding to a complex conjugate pair of eigenvalues w(j)
  82 *>          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
  83 *>          set to .TRUE..
  84 *>          If HOWMNY = 'A', SELECT is not referenced.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] N
  88 *> \verbatim
  89 *>          N is INTEGER
  90 *>          The order of the matrix T. N >= 0.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] T
  94 *> \verbatim
  95 *>          T is DOUBLE PRECISION array, dimension (LDT,N)
  96 *>          The upper quasi-triangular matrix T, in Schur canonical form.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] LDT
 100 *> \verbatim
 101 *>          LDT is INTEGER
 102 *>          The leading dimension of the array T. LDT >= max(1,N).
 103 *> \endverbatim
 104 *>
 105 *> \param[in] VL
 106 *> \verbatim
 107 *>          VL is DOUBLE PRECISION array, dimension (LDVL,M)
 108 *>          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
 109 *>          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
 110 *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
 111 *>          must be stored in consecutive columns of VL, as returned by
 112 *>          DHSEIN or DTREVC.
 113 *>          If JOB = 'V', VL is not referenced.
 114 *> \endverbatim
 115 *>
 116 *> \param[in] LDVL
 117 *> \verbatim
 118 *>          LDVL is INTEGER
 119 *>          The leading dimension of the array VL.
 120 *>          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
 121 *> \endverbatim
 122 *>
 123 *> \param[in] VR
 124 *> \verbatim
 125 *>          VR is DOUBLE PRECISION array, dimension (LDVR,M)
 126 *>          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
 127 *>          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
 128 *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
 129 *>          must be stored in consecutive columns of VR, as returned by
 130 *>          DHSEIN or DTREVC.
 131 *>          If JOB = 'V', VR is not referenced.
 132 *> \endverbatim
 133 *>
 134 *> \param[in] LDVR
 135 *> \verbatim
 136 *>          LDVR is INTEGER
 137 *>          The leading dimension of the array VR.
 138 *>          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
 139 *> \endverbatim
 140 *>
 141 *> \param[out] S
 142 *> \verbatim
 143 *>          S is DOUBLE PRECISION array, dimension (MM)
 144 *>          If JOB = 'E' or 'B', the reciprocal condition numbers of the
 145 *>          selected eigenvalues, stored in consecutive elements of the
 146 *>          array. For a complex conjugate pair of eigenvalues two
 147 *>          consecutive elements of S are set to the same value. Thus
 148 *>          S(j), SEP(j), and the j-th columns of VL and VR all
 149 *>          correspond to the same eigenpair (but not in general the
 150 *>          j-th eigenpair, unless all eigenpairs are selected).
 151 *>          If JOB = 'V', S is not referenced.
 152 *> \endverbatim
 153 *>
 154 *> \param[out] SEP
 155 *> \verbatim
 156 *>          SEP is DOUBLE PRECISION array, dimension (MM)
 157 *>          If JOB = 'V' or 'B', the estimated reciprocal condition
 158 *>          numbers of the selected eigenvectors, stored in consecutive
 159 *>          elements of the array. For a complex eigenvector two
 160 *>          consecutive elements of SEP are set to the same value. If
 161 *>          the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
 162 *>          is set to 0; this can only occur when the true value would be
 163 *>          very small anyway.
 164 *>          If JOB = 'E', SEP is not referenced.
 165 *> \endverbatim
 166 *>
 167 *> \param[in] MM
 168 *> \verbatim
 169 *>          MM is INTEGER
 170 *>          The number of elements in the arrays S (if JOB = 'E' or 'B')
 171 *>           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
 172 *> \endverbatim
 173 *>
 174 *> \param[out] M
 175 *> \verbatim
 176 *>          M is INTEGER
 177 *>          The number of elements of the arrays S and/or SEP actually
 178 *>          used to store the estimated condition numbers.
 179 *>          If HOWMNY = 'A', M is set to N.
 180 *> \endverbatim
 181 *>
 182 *> \param[out] WORK
 183 *> \verbatim
 184 *>          WORK is DOUBLE PRECISION array, dimension (LDWORK,N+6)
 185 *>          If JOB = 'E', WORK is not referenced.
 186 *> \endverbatim
 187 *>
 188 *> \param[in] LDWORK
 189 *> \verbatim
 190 *>          LDWORK is INTEGER
 191 *>          The leading dimension of the array WORK.
 192 *>          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
 193 *> \endverbatim
 194 *>
 195 *> \param[out] IWORK
 196 *> \verbatim
 197 *>          IWORK is INTEGER array, dimension (2*(N-1))
 198 *>          If JOB = 'E', IWORK is not referenced.
 199 *> \endverbatim
 200 *>
 201 *> \param[out] INFO
 202 *> \verbatim
 203 *>          INFO is INTEGER
 204 *>          = 0: successful exit
 205 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 206 *> \endverbatim
 207 *
 208 *  Authors:
 209 *  ========
 210 *
 211 *> \author Univ. of Tennessee
 212 *> \author Univ. of California Berkeley
 213 *> \author Univ. of Colorado Denver
 214 *> \author NAG Ltd.
 215 *
 216 *> \date November 2011
 217 *
 218 *> \ingroup doubleOTHERcomputational
 219 *
 220 *> \par Further Details:
 221 *  =====================
 222 *>
 223 *> \verbatim
 224 *>
 225 *>  The reciprocal of the condition number of an eigenvalue lambda is
 226 *>  defined as
 227 *>
 228 *>          S(lambda) = |v**T*u| / (norm(u)*norm(v))
 229 *>
 230 *>  where u and v are the right and left eigenvectors of T corresponding
 231 *>  to lambda; v**T denotes the transpose of v, and norm(u)
 232 *>  denotes the Euclidean norm. These reciprocal condition numbers always
 233 *>  lie between zero (very badly conditioned) and one (very well
 234 *>  conditioned). If n = 1, S(lambda) is defined to be 1.
 235 *>
 236 *>  An approximate error bound for a computed eigenvalue W(i) is given by
 237 *>
 238 *>                      EPS * norm(T) / S(i)
 239 *>
 240 *>  where EPS is the machine precision.
 241 *>
 242 *>  The reciprocal of the condition number of the right eigenvector u
 243 *>  corresponding to lambda is defined as follows. Suppose
 244 *>
 245 *>              T = ( lambda  c  )
 246 *>                  (   0    T22 )
 247 *>
 248 *>  Then the reciprocal condition number is
 249 *>
 250 *>          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
 251 *>
 252 *>  where sigma-min denotes the smallest singular value. We approximate
 253 *>  the smallest singular value by the reciprocal of an estimate of the
 254 *>  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
 255 *>  defined to be abs(T(1,1)).
 256 *>
 257 *>  An approximate error bound for a computed right eigenvector VR(i)
 258 *>  is given by
 259 *>
 260 *>                      EPS * norm(T) / SEP(i)
 261 *> \endverbatim
 262 *>
 263 *  =====================================================================
 264       SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
 265      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
 266      $                   INFO )
 267 *
 268 *  -- LAPACK computational routine (version 3.4.0) --
 269 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 270 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 271 *     November 2011
 272 *
 273 *     .. Scalar Arguments ..
 274       CHARACTER          HOWMNY, JOB
 275       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
 276 *     ..
 277 *     .. Array Arguments ..
 278       LOGICAL            SELECT( * )
 279       INTEGER            IWORK( * )
 280       DOUBLE PRECISION   S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
 281      $                   VR( LDVR, * ), WORK( LDWORK, * )
 282 *     ..
 283 *
 284 *  =====================================================================
 285 *
 286 *     .. Parameters ..
 287       DOUBLE PRECISION   ZERO, ONE, TWO
 288       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
 289 *     ..
 290 *     .. Local Scalars ..
 291       LOGICAL            PAIR, SOMCON, WANTBH, WANTS, WANTSP
 292       INTEGER            I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
 293       DOUBLE PRECISION   BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
 294      $                   MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
 295 *     ..
 296 *     .. Local Arrays ..
 297       INTEGER            ISAVE( 3 )
 298       DOUBLE PRECISION   DUMMY( 1 )
 299 *     ..
 300 *     .. External Functions ..
 301       LOGICAL            LSAME
 302       DOUBLE PRECISION   DDOT, DLAMCH, DLAPY2, DNRM2
 303       EXTERNAL           LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
 304 *     ..
 305 *     .. External Subroutines ..
 306       EXTERNAL           DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
 307 *     ..
 308 *     .. Intrinsic Functions ..
 309       INTRINSIC          ABS, MAX, SQRT
 310 *     ..
 311 *     .. Executable Statements ..
 312 *
 313 *     Decode and test the input parameters
 314 *
 315       WANTBH = LSAME( JOB, 'B' )
 316       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
 317       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
 318 *
 319       SOMCON = LSAME( HOWMNY, 'S' )
 320 *
 321       INFO = 0
 322       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
 323          INFO = -1
 324       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
 325          INFO = -2
 326       ELSE IF( N.LT.0 ) THEN
 327          INFO = -4
 328       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 329          INFO = -6
 330       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
 331          INFO = -8
 332       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
 333          INFO = -10
 334       ELSE
 335 *
 336 *        Set M to the number of eigenpairs for which condition numbers
 337 *        are required, and test MM.
 338 *
 339          IF( SOMCON ) THEN
 340             M = 0
 341             PAIR = .FALSE.
 342             DO 10 K = 1, N
 343                IF( PAIR ) THEN
 344                   PAIR = .FALSE.
 345                ELSE
 346                   IF( K.LT.N ) THEN
 347                      IF( T( K+1, K ).EQ.ZERO ) THEN
 348                         IF( SELECT( K ) )
 349      $                     M = M + 1
 350                      ELSE
 351                         PAIR = .TRUE.
 352                         IF( SELECT( K ) .OR. SELECT( K+1 ) )
 353      $                     M = M + 2
 354                      END IF
 355                   ELSE
 356                      IF( SELECT( N ) )
 357      $                  M = M + 1
 358                   END IF
 359                END IF
 360    10       CONTINUE
 361          ELSE
 362             M = N
 363          END IF
 364 *
 365          IF( MM.LT.M ) THEN
 366             INFO = -13
 367          ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
 368             INFO = -16
 369          END IF
 370       END IF
 371       IF( INFO.NE.0 ) THEN
 372          CALL XERBLA( 'DTRSNA', -INFO )
 373          RETURN
 374       END IF
 375 *
 376 *     Quick return if possible
 377 *
 378       IF( N.EQ.0 )
 379      $   RETURN
 380 *
 381       IF( N.EQ.1 ) THEN
 382          IF( SOMCON ) THEN
 383             IF( .NOT.SELECT( 1 ) )
 384      $         RETURN
 385          END IF
 386          IF( WANTS )
 387      $      S( 1 ) = ONE
 388          IF( WANTSP )
 389      $      SEP( 1 ) = ABS( T( 1, 1 ) )
 390          RETURN
 391       END IF
 392 *
 393 *     Get machine constants
 394 *
 395       EPS = DLAMCH( 'P' )
 396       SMLNUM = DLAMCH( 'S' ) / EPS
 397       BIGNUM = ONE / SMLNUM
 398       CALL DLABAD( SMLNUM, BIGNUM )
 399 *
 400       KS = 0
 401       PAIR = .FALSE.
 402       DO 60 K = 1, N
 403 *
 404 *        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
 405 *
 406          IF( PAIR ) THEN
 407             PAIR = .FALSE.
 408             GO TO 60
 409          ELSE
 410             IF( K.LT.N )
 411      $         PAIR = T( K+1, K ).NE.ZERO
 412          END IF
 413 *
 414 *        Determine whether condition numbers are required for the k-th
 415 *        eigenpair.
 416 *
 417          IF( SOMCON ) THEN
 418             IF( PAIR ) THEN
 419                IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
 420      $            GO TO 60
 421             ELSE
 422                IF( .NOT.SELECT( K ) )
 423      $            GO TO 60
 424             END IF
 425          END IF
 426 *
 427          KS = KS + 1
 428 *
 429          IF( WANTS ) THEN
 430 *
 431 *           Compute the reciprocal condition number of the k-th
 432 *           eigenvalue.
 433 *
 434             IF( .NOT.PAIR ) THEN
 435 *
 436 *              Real eigenvalue.
 437 *
 438                PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
 439                RNRM = DNRM2( N, VR( 1, KS ), 1 )
 440                LNRM = DNRM2( N, VL( 1, KS ), 1 )
 441                S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
 442             ELSE
 443 *
 444 *              Complex eigenvalue.
 445 *
 446                PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
 447                PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
 448      $                 1 )
 449                PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
 450                PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
 451      $                 1 )
 452                RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
 453      $                DNRM2( N, VR( 1, KS+1 ), 1 ) )
 454                LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
 455      $                DNRM2( N, VL( 1, KS+1 ), 1 ) )
 456                COND = DLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
 457                S( KS ) = COND
 458                S( KS+1 ) = COND
 459             END IF
 460          END IF
 461 *
 462          IF( WANTSP ) THEN
 463 *
 464 *           Estimate the reciprocal condition number of the k-th
 465 *           eigenvector.
 466 *
 467 *           Copy the matrix T to the array WORK and swap the diagonal
 468 *           block beginning at T(k,k) to the (1,1) position.
 469 *
 470             CALL DLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
 471             IFST = K
 472             ILST = 1
 473             CALL DTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
 474      $                   WORK( 1, N+1 ), IERR )
 475 *
 476             IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
 477 *
 478 *              Could not swap because blocks not well separated
 479 *
 480                SCALE = ONE
 481                EST = BIGNUM
 482             ELSE
 483 *
 484 *              Reordering successful
 485 *
 486                IF( WORK( 2, 1 ).EQ.ZERO ) THEN
 487 *
 488 *                 Form C = T22 - lambda*I in WORK(2:N,2:N).
 489 *
 490                   DO 20 I = 2, N
 491                      WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
 492    20             CONTINUE
 493                   N2 = 1
 494                   NN = N - 1
 495                ELSE
 496 *
 497 *                 Triangularize the 2 by 2 block by unitary
 498 *                 transformation U = [  cs   i*ss ]
 499 *                                    [ i*ss   cs  ].
 500 *                 such that the (1,1) position of WORK is complex
 501 *                 eigenvalue lambda with positive imaginary part. (2,2)
 502 *                 position of WORK is the complex eigenvalue lambda
 503 *                 with negative imaginary  part.
 504 *
 505                   MU = SQRT( ABS( WORK( 1, 2 ) ) )*
 506      $                 SQRT( ABS( WORK( 2, 1 ) ) )
 507                   DELTA = DLAPY2( MU, WORK( 2, 1 ) )
 508                   CS = MU / DELTA
 509                   SN = -WORK( 2, 1 ) / DELTA
 510 *
 511 *                 Form
 512 *
 513 *                 C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
 514 *                                          [   mu                     ]
 515 *                                          [         ..               ]
 516 *                                          [             ..           ]
 517 *                                          [                  mu      ]
 518 *                 where C**T is transpose of matrix C,
 519 *                 and RWORK is stored starting in the N+1-st column of
 520 *                 WORK.
 521 *
 522                   DO 30 J = 3, N
 523                      WORK( 2, J ) = CS*WORK( 2, J )
 524                      WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
 525    30             CONTINUE
 526                   WORK( 2, 2 ) = ZERO
 527 *
 528                   WORK( 1, N+1 ) = TWO*MU
 529                   DO 40 I = 2, N - 1
 530                      WORK( I, N+1 ) = SN*WORK( 1, I+1 )
 531    40             CONTINUE
 532                   N2 = 2
 533                   NN = 2*( N-1 )
 534                END IF
 535 *
 536 *              Estimate norm(inv(C**T))
 537 *
 538                EST = ZERO
 539                KASE = 0
 540    50          CONTINUE
 541                CALL DLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
 542      $                      EST, KASE, ISAVE )
 543                IF( KASE.NE.0 ) THEN
 544                   IF( KASE.EQ.1 ) THEN
 545                      IF( N2.EQ.1 ) THEN
 546 *
 547 *                       Real eigenvalue: solve C**T*x = scale*c.
 548 *
 549                         CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
 550      $                               LDWORK, DUMMY, DUMM, SCALE,
 551      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
 552      $                               IERR )
 553                      ELSE
 554 *
 555 *                       Complex eigenvalue: solve
 556 *                       C**T*(p+iq) = scale*(c+id) in real arithmetic.
 557 *
 558                         CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
 559      $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,
 560      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
 561      $                               IERR )
 562                      END IF
 563                   ELSE
 564                      IF( N2.EQ.1 ) THEN
 565 *
 566 *                       Real eigenvalue: solve C*x = scale*c.
 567 *
 568                         CALL DLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
 569      $                               LDWORK, DUMMY, DUMM, SCALE,
 570      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
 571      $                               IERR )
 572                      ELSE
 573 *
 574 *                       Complex eigenvalue: solve
 575 *                       C*(p+iq) = scale*(c+id) in real arithmetic.
 576 *
 577                         CALL DLAQTR( .FALSE., .FALSE., N-1,
 578      $                               WORK( 2, 2 ), LDWORK,
 579      $                               WORK( 1, N+1 ), MU, SCALE,
 580      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
 581      $                               IERR )
 582 *
 583                      END IF
 584                   END IF
 585 *
 586                   GO TO 50
 587                END IF
 588             END IF
 589 *
 590             SEP( KS ) = SCALE / MAX( EST, SMLNUM )
 591             IF( PAIR )
 592      $         SEP( KS+1 ) = SEP( KS )
 593          END IF
 594 *
 595          IF( PAIR )
 596      $      KS = KS + 1
 597 *
 598    60 CONTINUE
 599       RETURN
 600 *
 601 *     End of DTRSNA
 602 *
 603       END