LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
strsna.f
Go to the documentation of this file.
00001 *> \brief \b STRSNA
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STRSNA + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strsna.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strsna.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strsna.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00022 *                          LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
00023 *                          INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          HOWMNY, JOB
00027 *       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       LOGICAL            SELECT( * )
00031 *       INTEGER            IWORK( * )
00032 *       REAL               S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
00033 *      $                   VR( LDVR, * ), WORK( LDWORK, * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> STRSNA estimates reciprocal condition numbers for specified
00043 *> eigenvalues and/or right eigenvectors of a real upper
00044 *> quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
00045 *> orthogonal).
00046 *>
00047 *> T must be in Schur canonical form (as returned by SHSEQR), that is,
00048 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
00049 *> 2-by-2 diagonal block has its diagonal elements equal and its
00050 *> off-diagonal elements of opposite sign.
00051 *> \endverbatim
00052 *
00053 *  Arguments:
00054 *  ==========
00055 *
00056 *> \param[in] JOB
00057 *> \verbatim
00058 *>          JOB is CHARACTER*1
00059 *>          Specifies whether condition numbers are required for
00060 *>          eigenvalues (S) or eigenvectors (SEP):
00061 *>          = 'E': for eigenvalues only (S);
00062 *>          = 'V': for eigenvectors only (SEP);
00063 *>          = 'B': for both eigenvalues and eigenvectors (S and SEP).
00064 *> \endverbatim
00065 *>
00066 *> \param[in] HOWMNY
00067 *> \verbatim
00068 *>          HOWMNY is CHARACTER*1
00069 *>          = 'A': compute condition numbers for all eigenpairs;
00070 *>          = 'S': compute condition numbers for selected eigenpairs
00071 *>                 specified by the array SELECT.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] SELECT
00075 *> \verbatim
00076 *>          SELECT is LOGICAL array, dimension (N)
00077 *>          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
00078 *>          condition numbers are required. To select condition numbers
00079 *>          for the eigenpair corresponding to a real eigenvalue w(j),
00080 *>          SELECT(j) must be set to .TRUE.. To select condition numbers
00081 *>          corresponding to a complex conjugate pair of eigenvalues w(j)
00082 *>          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
00083 *>          set to .TRUE..
00084 *>          If HOWMNY = 'A', SELECT is not referenced.
00085 *> \endverbatim
00086 *>
00087 *> \param[in] N
00088 *> \verbatim
00089 *>          N is INTEGER
00090 *>          The order of the matrix T. N >= 0.
00091 *> \endverbatim
00092 *>
00093 *> \param[in] T
00094 *> \verbatim
00095 *>          T is REAL array, dimension (LDT,N)
00096 *>          The upper quasi-triangular matrix T, in Schur canonical form.
00097 *> \endverbatim
00098 *>
00099 *> \param[in] LDT
00100 *> \verbatim
00101 *>          LDT is INTEGER
00102 *>          The leading dimension of the array T. LDT >= max(1,N).
00103 *> \endverbatim
00104 *>
00105 *> \param[in] VL
00106 *> \verbatim
00107 *>          VL is REAL array, dimension (LDVL,M)
00108 *>          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
00109 *>          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
00110 *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00111 *>          must be stored in consecutive columns of VL, as returned by
00112 *>          SHSEIN or STREVC.
00113 *>          If JOB = 'V', VL is not referenced.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] LDVL
00117 *> \verbatim
00118 *>          LDVL is INTEGER
00119 *>          The leading dimension of the array VL.
00120 *>          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
00121 *> \endverbatim
00122 *>
00123 *> \param[in] VR
00124 *> \verbatim
00125 *>          VR is REAL array, dimension (LDVR,M)
00126 *>          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
00127 *>          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
00128 *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00129 *>          must be stored in consecutive columns of VR, as returned by
00130 *>          SHSEIN or STREVC.
00131 *>          If JOB = 'V', VR is not referenced.
00132 *> \endverbatim
00133 *>
00134 *> \param[in] LDVR
00135 *> \verbatim
00136 *>          LDVR is INTEGER
00137 *>          The leading dimension of the array VR.
00138 *>          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
00139 *> \endverbatim
00140 *>
00141 *> \param[out] S
00142 *> \verbatim
00143 *>          S is REAL array, dimension (MM)
00144 *>          If JOB = 'E' or 'B', the reciprocal condition numbers of the
00145 *>          selected eigenvalues, stored in consecutive elements of the
00146 *>          array. For a complex conjugate pair of eigenvalues two
00147 *>          consecutive elements of S are set to the same value. Thus
00148 *>          S(j), SEP(j), and the j-th columns of VL and VR all
00149 *>          correspond to the same eigenpair (but not in general the
00150 *>          j-th eigenpair, unless all eigenpairs are selected).
00151 *>          If JOB = 'V', S is not referenced.
00152 *> \endverbatim
00153 *>
00154 *> \param[out] SEP
00155 *> \verbatim
00156 *>          SEP is REAL array, dimension (MM)
00157 *>          If JOB = 'V' or 'B', the estimated reciprocal condition
00158 *>          numbers of the selected eigenvectors, stored in consecutive
00159 *>          elements of the array. For a complex eigenvector two
00160 *>          consecutive elements of SEP are set to the same value. If
00161 *>          the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
00162 *>          is set to 0; this can only occur when the true value would be
00163 *>          very small anyway.
00164 *>          If JOB = 'E', SEP is not referenced.
00165 *> \endverbatim
00166 *>
00167 *> \param[in] MM
00168 *> \verbatim
00169 *>          MM is INTEGER
00170 *>          The number of elements in the arrays S (if JOB = 'E' or 'B')
00171 *>           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
00172 *> \endverbatim
00173 *>
00174 *> \param[out] M
00175 *> \verbatim
00176 *>          M is INTEGER
00177 *>          The number of elements of the arrays S and/or SEP actually
00178 *>          used to store the estimated condition numbers.
00179 *>          If HOWMNY = 'A', M is set to N.
00180 *> \endverbatim
00181 *>
00182 *> \param[out] WORK
00183 *> \verbatim
00184 *>          WORK is REAL array, dimension (LDWORK,N+6)
00185 *>          If JOB = 'E', WORK is not referenced.
00186 *> \endverbatim
00187 *>
00188 *> \param[in] LDWORK
00189 *> \verbatim
00190 *>          LDWORK is INTEGER
00191 *>          The leading dimension of the array WORK.
00192 *>          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
00193 *> \endverbatim
00194 *>
00195 *> \param[out] IWORK
00196 *> \verbatim
00197 *>          IWORK is INTEGER array, dimension (2*(N-1))
00198 *>          If JOB = 'E', IWORK is not referenced.
00199 *> \endverbatim
00200 *>
00201 *> \param[out] INFO
00202 *> \verbatim
00203 *>          INFO is INTEGER
00204 *>          = 0: successful exit
00205 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00206 *> \endverbatim
00207 *
00208 *  Authors:
00209 *  ========
00210 *
00211 *> \author Univ. of Tennessee 
00212 *> \author Univ. of California Berkeley 
00213 *> \author Univ. of Colorado Denver 
00214 *> \author NAG Ltd. 
00215 *
00216 *> \date November 2011
00217 *
00218 *> \ingroup realOTHERcomputational
00219 *
00220 *> \par Further Details:
00221 *  =====================
00222 *>
00223 *> \verbatim
00224 *>
00225 *>  The reciprocal of the condition number of an eigenvalue lambda is
00226 *>  defined as
00227 *>
00228 *>          S(lambda) = |v**T*u| / (norm(u)*norm(v))
00229 *>
00230 *>  where u and v are the right and left eigenvectors of T corresponding
00231 *>  to lambda; v**T denotes the transpose of v, and norm(u)
00232 *>  denotes the Euclidean norm. These reciprocal condition numbers always
00233 *>  lie between zero (very badly conditioned) and one (very well
00234 *>  conditioned). If n = 1, S(lambda) is defined to be 1.
00235 *>
00236 *>  An approximate error bound for a computed eigenvalue W(i) is given by
00237 *>
00238 *>                      EPS * norm(T) / S(i)
00239 *>
00240 *>  where EPS is the machine precision.
00241 *>
00242 *>  The reciprocal of the condition number of the right eigenvector u
00243 *>  corresponding to lambda is defined as follows. Suppose
00244 *>
00245 *>              T = ( lambda  c  )
00246 *>                  (   0    T22 )
00247 *>
00248 *>  Then the reciprocal condition number is
00249 *>
00250 *>          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
00251 *>
00252 *>  where sigma-min denotes the smallest singular value. We approximate
00253 *>  the smallest singular value by the reciprocal of an estimate of the
00254 *>  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
00255 *>  defined to be abs(T(1,1)).
00256 *>
00257 *>  An approximate error bound for a computed right eigenvector VR(i)
00258 *>  is given by
00259 *>
00260 *>                      EPS * norm(T) / SEP(i)
00261 *> \endverbatim
00262 *>
00263 *  =====================================================================
00264       SUBROUTINE STRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00265      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
00266      $                   INFO )
00267 *
00268 *  -- LAPACK computational routine (version 3.4.0) --
00269 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00270 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00271 *     November 2011
00272 *
00273 *     .. Scalar Arguments ..
00274       CHARACTER          HOWMNY, JOB
00275       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
00276 *     ..
00277 *     .. Array Arguments ..
00278       LOGICAL            SELECT( * )
00279       INTEGER            IWORK( * )
00280       REAL               S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
00281      $                   VR( LDVR, * ), WORK( LDWORK, * )
00282 *     ..
00283 *
00284 *  =====================================================================
00285 *
00286 *     .. Parameters ..
00287       REAL               ZERO, ONE, TWO
00288       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
00289 *     ..
00290 *     .. Local Scalars ..
00291       LOGICAL            PAIR, SOMCON, WANTBH, WANTS, WANTSP
00292       INTEGER            I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
00293       REAL               BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
00294      $                   MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
00295 *     ..
00296 *     .. Local Arrays ..
00297       INTEGER            ISAVE( 3 )
00298       REAL               DUMMY( 1 )
00299 *     ..
00300 *     .. External Functions ..
00301       LOGICAL            LSAME
00302       REAL               SDOT, SLAMCH, SLAPY2, SNRM2
00303       EXTERNAL           LSAME, SDOT, SLAMCH, SLAPY2, SNRM2
00304 *     ..
00305 *     .. External Subroutines ..
00306       EXTERNAL           SLABAD, SLACN2, SLACPY, SLAQTR, STREXC, XERBLA
00307 *     ..
00308 *     .. Intrinsic Functions ..
00309       INTRINSIC          ABS, MAX, SQRT
00310 *     ..
00311 *     .. Executable Statements ..
00312 *
00313 *     Decode and test the input parameters
00314 *
00315       WANTBH = LSAME( JOB, 'B' )
00316       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00317       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
00318 *
00319       SOMCON = LSAME( HOWMNY, 'S' )
00320 *
00321       INFO = 0
00322       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
00323          INFO = -1
00324       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
00325          INFO = -2
00326       ELSE IF( N.LT.0 ) THEN
00327          INFO = -4
00328       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00329          INFO = -6
00330       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
00331          INFO = -8
00332       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
00333          INFO = -10
00334       ELSE
00335 *
00336 *        Set M to the number of eigenpairs for which condition numbers
00337 *        are required, and test MM.
00338 *
00339          IF( SOMCON ) THEN
00340             M = 0
00341             PAIR = .FALSE.
00342             DO 10 K = 1, N
00343                IF( PAIR ) THEN
00344                   PAIR = .FALSE.
00345                ELSE
00346                   IF( K.LT.N ) THEN
00347                      IF( T( K+1, K ).EQ.ZERO ) THEN
00348                         IF( SELECT( K ) )
00349      $                     M = M + 1
00350                      ELSE
00351                         PAIR = .TRUE.
00352                         IF( SELECT( K ) .OR. SELECT( K+1 ) )
00353      $                     M = M + 2
00354                      END IF
00355                   ELSE
00356                      IF( SELECT( N ) )
00357      $                  M = M + 1
00358                   END IF
00359                END IF
00360    10       CONTINUE
00361          ELSE
00362             M = N
00363          END IF
00364 *
00365          IF( MM.LT.M ) THEN
00366             INFO = -13
00367          ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
00368             INFO = -16
00369          END IF
00370       END IF
00371       IF( INFO.NE.0 ) THEN
00372          CALL XERBLA( 'STRSNA', -INFO )
00373          RETURN
00374       END IF
00375 *
00376 *     Quick return if possible
00377 *
00378       IF( N.EQ.0 )
00379      $   RETURN
00380 *
00381       IF( N.EQ.1 ) THEN
00382          IF( SOMCON ) THEN
00383             IF( .NOT.SELECT( 1 ) )
00384      $         RETURN
00385          END IF
00386          IF( WANTS )
00387      $      S( 1 ) = ONE
00388          IF( WANTSP )
00389      $      SEP( 1 ) = ABS( T( 1, 1 ) )
00390          RETURN
00391       END IF
00392 *
00393 *     Get machine constants
00394 *
00395       EPS = SLAMCH( 'P' )
00396       SMLNUM = SLAMCH( 'S' ) / EPS
00397       BIGNUM = ONE / SMLNUM
00398       CALL SLABAD( SMLNUM, BIGNUM )
00399 *
00400       KS = 0
00401       PAIR = .FALSE.
00402       DO 60 K = 1, N
00403 *
00404 *        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
00405 *
00406          IF( PAIR ) THEN
00407             PAIR = .FALSE.
00408             GO TO 60
00409          ELSE
00410             IF( K.LT.N )
00411      $         PAIR = T( K+1, K ).NE.ZERO
00412          END IF
00413 *
00414 *        Determine whether condition numbers are required for the k-th
00415 *        eigenpair.
00416 *
00417          IF( SOMCON ) THEN
00418             IF( PAIR ) THEN
00419                IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
00420      $            GO TO 60
00421             ELSE
00422                IF( .NOT.SELECT( K ) )
00423      $            GO TO 60
00424             END IF
00425          END IF
00426 *
00427          KS = KS + 1
00428 *
00429          IF( WANTS ) THEN
00430 *
00431 *           Compute the reciprocal condition number of the k-th
00432 *           eigenvalue.
00433 *
00434             IF( .NOT.PAIR ) THEN
00435 *
00436 *              Real eigenvalue.
00437 *
00438                PROD = SDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
00439                RNRM = SNRM2( N, VR( 1, KS ), 1 )
00440                LNRM = SNRM2( N, VL( 1, KS ), 1 )
00441                S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
00442             ELSE
00443 *
00444 *              Complex eigenvalue.
00445 *
00446                PROD1 = SDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
00447                PROD1 = PROD1 + SDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
00448      $                 1 )
00449                PROD2 = SDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
00450                PROD2 = PROD2 - SDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
00451      $                 1 )
00452                RNRM = SLAPY2( SNRM2( N, VR( 1, KS ), 1 ),
00453      $                SNRM2( N, VR( 1, KS+1 ), 1 ) )
00454                LNRM = SLAPY2( SNRM2( N, VL( 1, KS ), 1 ),
00455      $                SNRM2( N, VL( 1, KS+1 ), 1 ) )
00456                COND = SLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
00457                S( KS ) = COND
00458                S( KS+1 ) = COND
00459             END IF
00460          END IF
00461 *
00462          IF( WANTSP ) THEN
00463 *
00464 *           Estimate the reciprocal condition number of the k-th
00465 *           eigenvector.
00466 *
00467 *           Copy the matrix T to the array WORK and swap the diagonal
00468 *           block beginning at T(k,k) to the (1,1) position.
00469 *
00470             CALL SLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
00471             IFST = K
00472             ILST = 1
00473             CALL STREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
00474      $                   WORK( 1, N+1 ), IERR )
00475 *
00476             IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
00477 *
00478 *              Could not swap because blocks not well separated
00479 *
00480                SCALE = ONE
00481                EST = BIGNUM
00482             ELSE
00483 *
00484 *              Reordering successful
00485 *
00486                IF( WORK( 2, 1 ).EQ.ZERO ) THEN
00487 *
00488 *                 Form C = T22 - lambda*I in WORK(2:N,2:N).
00489 *
00490                   DO 20 I = 2, N
00491                      WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
00492    20             CONTINUE
00493                   N2 = 1
00494                   NN = N - 1
00495                ELSE
00496 *
00497 *                 Triangularize the 2 by 2 block by unitary
00498 *                 transformation U = [  cs   i*ss ]
00499 *                                    [ i*ss   cs  ].
00500 *                 such that the (1,1) position of WORK is complex
00501 *                 eigenvalue lambda with positive imaginary part. (2,2)
00502 *                 position of WORK is the complex eigenvalue lambda
00503 *                 with negative imaginary  part.
00504 *
00505                   MU = SQRT( ABS( WORK( 1, 2 ) ) )*
00506      $                 SQRT( ABS( WORK( 2, 1 ) ) )
00507                   DELTA = SLAPY2( MU, WORK( 2, 1 ) )
00508                   CS = MU / DELTA
00509                   SN = -WORK( 2, 1 ) / DELTA
00510 *
00511 *                 Form
00512 *
00513 *                 C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
00514 *                                          [   mu                     ]
00515 *                                          [         ..               ]
00516 *                                          [             ..           ]
00517 *                                          [                  mu      ]
00518 *                 where C**T is transpose of matrix C,
00519 *                 and RWORK is stored starting in the N+1-st column of
00520 *                 WORK.
00521 *
00522                   DO 30 J = 3, N
00523                      WORK( 2, J ) = CS*WORK( 2, J )
00524                      WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
00525    30             CONTINUE
00526                   WORK( 2, 2 ) = ZERO
00527 *
00528                   WORK( 1, N+1 ) = TWO*MU
00529                   DO 40 I = 2, N - 1
00530                      WORK( I, N+1 ) = SN*WORK( 1, I+1 )
00531    40             CONTINUE
00532                   N2 = 2
00533                   NN = 2*( N-1 )
00534                END IF
00535 *
00536 *              Estimate norm(inv(C**T))
00537 *
00538                EST = ZERO
00539                KASE = 0
00540    50          CONTINUE
00541                CALL SLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
00542      $                      EST, KASE, ISAVE )
00543                IF( KASE.NE.0 ) THEN
00544                   IF( KASE.EQ.1 ) THEN
00545                      IF( N2.EQ.1 ) THEN
00546 *
00547 *                       Real eigenvalue: solve C**T*x = scale*c.
00548 *
00549                         CALL SLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
00550      $                               LDWORK, DUMMY, DUMM, SCALE,
00551      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
00552      $                               IERR )
00553                      ELSE
00554 *
00555 *                       Complex eigenvalue: solve
00556 *                       C**T*(p+iq) = scale*(c+id) in real arithmetic.
00557 *
00558                         CALL SLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
00559      $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,
00560      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
00561      $                               IERR )
00562                      END IF
00563                   ELSE
00564                      IF( N2.EQ.1 ) THEN
00565 *
00566 *                       Real eigenvalue: solve C*x = scale*c.
00567 *
00568                         CALL SLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
00569      $                               LDWORK, DUMMY, DUMM, SCALE,
00570      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
00571      $                               IERR )
00572                      ELSE
00573 *
00574 *                       Complex eigenvalue: solve
00575 *                       C*(p+iq) = scale*(c+id) in real arithmetic.
00576 *
00577                         CALL SLAQTR( .FALSE., .FALSE., N-1,
00578      $                               WORK( 2, 2 ), LDWORK,
00579      $                               WORK( 1, N+1 ), MU, SCALE,
00580      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
00581      $                               IERR )
00582 *
00583                      END IF
00584                   END IF
00585 *
00586                   GO TO 50
00587                END IF
00588             END IF
00589 *
00590             SEP( KS ) = SCALE / MAX( EST, SMLNUM )
00591             IF( PAIR )
00592      $         SEP( KS+1 ) = SEP( KS )
00593          END IF
00594 *
00595          IF( PAIR )
00596      $      KS = KS + 1
00597 *
00598    60 CONTINUE
00599       RETURN
00600 *
00601 *     End of STRSNA
00602 *
00603       END
 All Files Functions