![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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