|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE SGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, 00002 $ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, 00003 $ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO ) 00004 * 00005 * -- LAPACK driver routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER BALANC, JOBVL, JOBVR, SENSE 00012 INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N 00013 REAL ABNRM 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IWORK( * ) 00017 REAL A( LDA, * ), RCONDE( * ), RCONDV( * ), 00018 $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ), 00019 $ WI( * ), WORK( * ), WR( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * SGEEVX computes for an N-by-N real nonsymmetric matrix A, the 00026 * eigenvalues and, optionally, the left and/or right eigenvectors. 00027 * 00028 * Optionally also, it computes a balancing transformation to improve 00029 * the conditioning of the eigenvalues and eigenvectors (ILO, IHI, 00030 * SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues 00031 * (RCONDE), and reciprocal condition numbers for the right 00032 * eigenvectors (RCONDV). 00033 * 00034 * The right eigenvector v(j) of A satisfies 00035 * A * v(j) = lambda(j) * v(j) 00036 * where lambda(j) is its eigenvalue. 00037 * The left eigenvector u(j) of A satisfies 00038 * u(j)**H * A = lambda(j) * u(j)**H 00039 * where u(j)**H denotes the conjugate transpose of u(j). 00040 * 00041 * The computed eigenvectors are normalized to have Euclidean norm 00042 * equal to 1 and largest component real. 00043 * 00044 * Balancing a matrix means permuting the rows and columns to make it 00045 * more nearly upper triangular, and applying a diagonal similarity 00046 * transformation D * A * D**(-1), where D is a diagonal matrix, to 00047 * make its rows and columns closer in norm and the condition numbers 00048 * of its eigenvalues and eigenvectors smaller. The computed 00049 * reciprocal condition numbers correspond to the balanced matrix. 00050 * Permuting rows and columns will not change the condition numbers 00051 * (in exact arithmetic) but diagonal scaling will. For further 00052 * explanation of balancing, see section 4.10.2 of the LAPACK 00053 * Users' Guide. 00054 * 00055 * Arguments 00056 * ========= 00057 * 00058 * BALANC (input) CHARACTER*1 00059 * Indicates how the input matrix should be diagonally scaled 00060 * and/or permuted to improve the conditioning of its 00061 * eigenvalues. 00062 * = 'N': Do not diagonally scale or permute; 00063 * = 'P': Perform permutations to make the matrix more nearly 00064 * upper triangular. Do not diagonally scale; 00065 * = 'S': Diagonally scale the matrix, i.e. replace A by 00066 * D*A*D**(-1), where D is a diagonal matrix chosen 00067 * to make the rows and columns of A more equal in 00068 * norm. Do not permute; 00069 * = 'B': Both diagonally scale and permute A. 00070 * 00071 * Computed reciprocal condition numbers will be for the matrix 00072 * after balancing and/or permuting. Permuting does not change 00073 * condition numbers (in exact arithmetic), but balancing does. 00074 * 00075 * JOBVL (input) CHARACTER*1 00076 * = 'N': left eigenvectors of A are not computed; 00077 * = 'V': left eigenvectors of A are computed. 00078 * If SENSE = 'E' or 'B', JOBVL must = 'V'. 00079 * 00080 * JOBVR (input) CHARACTER*1 00081 * = 'N': right eigenvectors of A are not computed; 00082 * = 'V': right eigenvectors of A are computed. 00083 * If SENSE = 'E' or 'B', JOBVR must = 'V'. 00084 * 00085 * SENSE (input) CHARACTER*1 00086 * Determines which reciprocal condition numbers are computed. 00087 * = 'N': None are computed; 00088 * = 'E': Computed for eigenvalues only; 00089 * = 'V': Computed for right eigenvectors only; 00090 * = 'B': Computed for eigenvalues and right eigenvectors. 00091 * 00092 * If SENSE = 'E' or 'B', both left and right eigenvectors 00093 * must also be computed (JOBVL = 'V' and JOBVR = 'V'). 00094 * 00095 * N (input) INTEGER 00096 * The order of the matrix A. N >= 0. 00097 * 00098 * A (input/output) REAL array, dimension (LDA,N) 00099 * On entry, the N-by-N matrix A. 00100 * On exit, A has been overwritten. If JOBVL = 'V' or 00101 * JOBVR = 'V', A contains the real Schur form of the balanced 00102 * version of the input matrix A. 00103 * 00104 * LDA (input) INTEGER 00105 * The leading dimension of the array A. LDA >= max(1,N). 00106 * 00107 * WR (output) REAL array, dimension (N) 00108 * WI (output) REAL array, dimension (N) 00109 * WR and WI contain the real and imaginary parts, 00110 * respectively, of the computed eigenvalues. Complex 00111 * conjugate pairs of eigenvalues will appear consecutively 00112 * with the eigenvalue having the positive imaginary part 00113 * first. 00114 * 00115 * VL (output) REAL array, dimension (LDVL,N) 00116 * If JOBVL = 'V', the left eigenvectors u(j) are stored one 00117 * after another in the columns of VL, in the same order 00118 * as their eigenvalues. 00119 * If JOBVL = 'N', VL is not referenced. 00120 * If the j-th eigenvalue is real, then u(j) = VL(:,j), 00121 * the j-th column of VL. 00122 * If the j-th and (j+1)-st eigenvalues form a complex 00123 * conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and 00124 * u(j+1) = VL(:,j) - i*VL(:,j+1). 00125 * 00126 * LDVL (input) INTEGER 00127 * The leading dimension of the array VL. LDVL >= 1; if 00128 * JOBVL = 'V', LDVL >= N. 00129 * 00130 * VR (output) REAL array, dimension (LDVR,N) 00131 * If JOBVR = 'V', the right eigenvectors v(j) are stored one 00132 * after another in the columns of VR, in the same order 00133 * as their eigenvalues. 00134 * If JOBVR = 'N', VR is not referenced. 00135 * If the j-th eigenvalue is real, then v(j) = VR(:,j), 00136 * the j-th column of VR. 00137 * If the j-th and (j+1)-st eigenvalues form a complex 00138 * conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and 00139 * v(j+1) = VR(:,j) - i*VR(:,j+1). 00140 * 00141 * LDVR (input) INTEGER 00142 * The leading dimension of the array VR. LDVR >= 1, and if 00143 * JOBVR = 'V', LDVR >= N. 00144 * 00145 * ILO (output) INTEGER 00146 * IHI (output) INTEGER 00147 * ILO and IHI are integer values determined when A was 00148 * balanced. The balanced A(i,j) = 0 if I > J and 00149 * J = 1,...,ILO-1 or I = IHI+1,...,N. 00150 * 00151 * SCALE (output) REAL array, dimension (N) 00152 * Details of the permutations and scaling factors applied 00153 * when balancing A. If P(j) is the index of the row and column 00154 * interchanged with row and column j, and D(j) is the scaling 00155 * factor applied to row and column j, then 00156 * SCALE(J) = P(J), for J = 1,...,ILO-1 00157 * = D(J), for J = ILO,...,IHI 00158 * = P(J) for J = IHI+1,...,N. 00159 * The order in which the interchanges are made is N to IHI+1, 00160 * then 1 to ILO-1. 00161 * 00162 * ABNRM (output) REAL 00163 * The one-norm of the balanced matrix (the maximum 00164 * of the sum of absolute values of elements of any column). 00165 * 00166 * RCONDE (output) REAL array, dimension (N) 00167 * RCONDE(j) is the reciprocal condition number of the j-th 00168 * eigenvalue. 00169 * 00170 * RCONDV (output) REAL array, dimension (N) 00171 * RCONDV(j) is the reciprocal condition number of the j-th 00172 * right eigenvector. 00173 * 00174 * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) 00175 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00176 * 00177 * LWORK (input) INTEGER 00178 * The dimension of the array WORK. If SENSE = 'N' or 'E', 00179 * LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V', 00180 * LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6). 00181 * For good performance, LWORK must generally be larger. 00182 * 00183 * If LWORK = -1, then a workspace query is assumed; the routine 00184 * only calculates the optimal size of the WORK array, returns 00185 * this value as the first entry of the WORK array, and no error 00186 * message related to LWORK is issued by XERBLA. 00187 * 00188 * IWORK (workspace) INTEGER array, dimension (2*N-2) 00189 * If SENSE = 'N' or 'E', not referenced. 00190 * 00191 * INFO (output) INTEGER 00192 * = 0: successful exit 00193 * < 0: if INFO = -i, the i-th argument had an illegal value. 00194 * > 0: if INFO = i, the QR algorithm failed to compute all the 00195 * eigenvalues, and no eigenvectors or condition numbers 00196 * have been computed; elements 1:ILO-1 and i+1:N of WR 00197 * and WI contain eigenvalues which have converged. 00198 * 00199 * ===================================================================== 00200 * 00201 * .. Parameters .. 00202 REAL ZERO, ONE 00203 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00204 * .. 00205 * .. Local Scalars .. 00206 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE, 00207 $ WNTSNN, WNTSNV 00208 CHARACTER JOB, SIDE 00209 INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK, 00210 $ MINWRK, NOUT 00211 REAL ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM, 00212 $ SN 00213 * .. 00214 * .. Local Arrays .. 00215 LOGICAL SELECT( 1 ) 00216 REAL DUM( 1 ) 00217 * .. 00218 * .. External Subroutines .. 00219 EXTERNAL SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD, SLACPY, 00220 $ SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC, 00221 $ STRSNA, XERBLA 00222 * .. 00223 * .. External Functions .. 00224 LOGICAL LSAME 00225 INTEGER ILAENV, ISAMAX 00226 REAL SLAMCH, SLANGE, SLAPY2, SNRM2 00227 EXTERNAL LSAME, ILAENV, ISAMAX, SLAMCH, SLANGE, SLAPY2, 00228 $ SNRM2 00229 * .. 00230 * .. Intrinsic Functions .. 00231 INTRINSIC MAX, SQRT 00232 * .. 00233 * .. Executable Statements .. 00234 * 00235 * Test the input arguments 00236 * 00237 INFO = 0 00238 LQUERY = ( LWORK.EQ.-1 ) 00239 WANTVL = LSAME( JOBVL, 'V' ) 00240 WANTVR = LSAME( JOBVR, 'V' ) 00241 WNTSNN = LSAME( SENSE, 'N' ) 00242 WNTSNE = LSAME( SENSE, 'E' ) 00243 WNTSNV = LSAME( SENSE, 'V' ) 00244 WNTSNB = LSAME( SENSE, 'B' ) 00245 IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' ) .OR. 00246 $ LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) THEN 00247 INFO = -1 00248 ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN 00249 INFO = -2 00250 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN 00251 INFO = -3 00252 ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR. 00253 $ ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND. 00254 $ WANTVR ) ) ) THEN 00255 INFO = -4 00256 ELSE IF( N.LT.0 ) THEN 00257 INFO = -5 00258 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00259 INFO = -7 00260 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN 00261 INFO = -11 00262 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN 00263 INFO = -13 00264 END IF 00265 * 00266 * Compute workspace 00267 * (Note: Comments in the code beginning "Workspace:" describe the 00268 * minimal amount of workspace needed at that point in the code, 00269 * as well as the preferred amount for good performance. 00270 * NB refers to the optimal block size for the immediately 00271 * following subroutine, as returned by ILAENV. 00272 * HSWORK refers to the workspace preferred by SHSEQR, as 00273 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 00274 * the worst case.) 00275 * 00276 IF( INFO.EQ.0 ) THEN 00277 IF( N.EQ.0 ) THEN 00278 MINWRK = 1 00279 MAXWRK = 1 00280 ELSE 00281 MAXWRK = N + N*ILAENV( 1, 'SGEHRD', ' ', N, 1, N, 0 ) 00282 * 00283 IF( WANTVL ) THEN 00284 CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL, 00285 $ WORK, -1, INFO ) 00286 ELSE IF( WANTVR ) THEN 00287 CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR, 00288 $ WORK, -1, INFO ) 00289 ELSE 00290 IF( WNTSNN ) THEN 00291 CALL SHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, 00292 $ LDVR, WORK, -1, INFO ) 00293 ELSE 00294 CALL SHSEQR( 'S', 'N', N, 1, N, A, LDA, WR, WI, VR, 00295 $ LDVR, WORK, -1, INFO ) 00296 END IF 00297 END IF 00298 HSWORK = WORK( 1 ) 00299 * 00300 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN 00301 MINWRK = 2*N 00302 IF( .NOT.WNTSNN ) 00303 $ MINWRK = MAX( MINWRK, N*N+6*N ) 00304 MAXWRK = MAX( MAXWRK, HSWORK ) 00305 IF( .NOT.WNTSNN ) 00306 $ MAXWRK = MAX( MAXWRK, N*N + 6*N ) 00307 ELSE 00308 MINWRK = 3*N 00309 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) ) 00310 $ MINWRK = MAX( MINWRK, N*N + 6*N ) 00311 MAXWRK = MAX( MAXWRK, HSWORK ) 00312 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'SORGHR', 00313 $ ' ', N, 1, N, -1 ) ) 00314 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) ) 00315 $ MAXWRK = MAX( MAXWRK, N*N + 6*N ) 00316 MAXWRK = MAX( MAXWRK, 3*N ) 00317 END IF 00318 MAXWRK = MAX( MAXWRK, MINWRK ) 00319 END IF 00320 WORK( 1 ) = MAXWRK 00321 * 00322 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00323 INFO = -21 00324 END IF 00325 END IF 00326 * 00327 IF( INFO.NE.0 ) THEN 00328 CALL XERBLA( 'SGEEVX', -INFO ) 00329 RETURN 00330 ELSE IF( LQUERY ) THEN 00331 RETURN 00332 END IF 00333 * 00334 * Quick return if possible 00335 * 00336 IF( N.EQ.0 ) 00337 $ RETURN 00338 * 00339 * Get machine constants 00340 * 00341 EPS = SLAMCH( 'P' ) 00342 SMLNUM = SLAMCH( 'S' ) 00343 BIGNUM = ONE / SMLNUM 00344 CALL SLABAD( SMLNUM, BIGNUM ) 00345 SMLNUM = SQRT( SMLNUM ) / EPS 00346 BIGNUM = ONE / SMLNUM 00347 * 00348 * Scale A if max element outside range [SMLNUM,BIGNUM] 00349 * 00350 ICOND = 0 00351 ANRM = SLANGE( 'M', N, N, A, LDA, DUM ) 00352 SCALEA = .FALSE. 00353 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00354 SCALEA = .TRUE. 00355 CSCALE = SMLNUM 00356 ELSE IF( ANRM.GT.BIGNUM ) THEN 00357 SCALEA = .TRUE. 00358 CSCALE = BIGNUM 00359 END IF 00360 IF( SCALEA ) 00361 $ CALL SLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 00362 * 00363 * Balance the matrix and compute ABNRM 00364 * 00365 CALL SGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR ) 00366 ABNRM = SLANGE( '1', N, N, A, LDA, DUM ) 00367 IF( SCALEA ) THEN 00368 DUM( 1 ) = ABNRM 00369 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR ) 00370 ABNRM = DUM( 1 ) 00371 END IF 00372 * 00373 * Reduce to upper Hessenberg form 00374 * (Workspace: need 2*N, prefer N+N*NB) 00375 * 00376 ITAU = 1 00377 IWRK = ITAU + N 00378 CALL SGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 00379 $ LWORK-IWRK+1, IERR ) 00380 * 00381 IF( WANTVL ) THEN 00382 * 00383 * Want left eigenvectors 00384 * Copy Householder vectors to VL 00385 * 00386 SIDE = 'L' 00387 CALL SLACPY( 'L', N, N, A, LDA, VL, LDVL ) 00388 * 00389 * Generate orthogonal matrix in VL 00390 * (Workspace: need 2*N-1, prefer N+(N-1)*NB) 00391 * 00392 CALL SORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ), 00393 $ LWORK-IWRK+1, IERR ) 00394 * 00395 * Perform QR iteration, accumulating Schur vectors in VL 00396 * (Workspace: need 1, prefer HSWORK (see comments) ) 00397 * 00398 IWRK = ITAU 00399 CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL, 00400 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00401 * 00402 IF( WANTVR ) THEN 00403 * 00404 * Want left and right eigenvectors 00405 * Copy Schur vectors to VR 00406 * 00407 SIDE = 'B' 00408 CALL SLACPY( 'F', N, N, VL, LDVL, VR, LDVR ) 00409 END IF 00410 * 00411 ELSE IF( WANTVR ) THEN 00412 * 00413 * Want right eigenvectors 00414 * Copy Householder vectors to VR 00415 * 00416 SIDE = 'R' 00417 CALL SLACPY( 'L', N, N, A, LDA, VR, LDVR ) 00418 * 00419 * Generate orthogonal matrix in VR 00420 * (Workspace: need 2*N-1, prefer N+(N-1)*NB) 00421 * 00422 CALL SORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ), 00423 $ LWORK-IWRK+1, IERR ) 00424 * 00425 * Perform QR iteration, accumulating Schur vectors in VR 00426 * (Workspace: need 1, prefer HSWORK (see comments) ) 00427 * 00428 IWRK = ITAU 00429 CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 00430 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00431 * 00432 ELSE 00433 * 00434 * Compute eigenvalues only 00435 * If condition numbers desired, compute Schur form 00436 * 00437 IF( WNTSNN ) THEN 00438 JOB = 'E' 00439 ELSE 00440 JOB = 'S' 00441 END IF 00442 * 00443 * (Workspace: need 1, prefer HSWORK (see comments) ) 00444 * 00445 IWRK = ITAU 00446 CALL SHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 00447 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00448 END IF 00449 * 00450 * If INFO > 0 from SHSEQR, then quit 00451 * 00452 IF( INFO.GT.0 ) 00453 $ GO TO 50 00454 * 00455 IF( WANTVL .OR. WANTVR ) THEN 00456 * 00457 * Compute left and/or right eigenvectors 00458 * (Workspace: need 3*N) 00459 * 00460 CALL STREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 00461 $ N, NOUT, WORK( IWRK ), IERR ) 00462 END IF 00463 * 00464 * Compute condition numbers if desired 00465 * (Workspace: need N*N+6*N unless SENSE = 'E') 00466 * 00467 IF( .NOT.WNTSNN ) THEN 00468 CALL STRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 00469 $ RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, IWORK, 00470 $ ICOND ) 00471 END IF 00472 * 00473 IF( WANTVL ) THEN 00474 * 00475 * Undo balancing of left eigenvectors 00476 * 00477 CALL SGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL, 00478 $ IERR ) 00479 * 00480 * Normalize left eigenvectors and make largest component real 00481 * 00482 DO 20 I = 1, N 00483 IF( WI( I ).EQ.ZERO ) THEN 00484 SCL = ONE / SNRM2( N, VL( 1, I ), 1 ) 00485 CALL SSCAL( N, SCL, VL( 1, I ), 1 ) 00486 ELSE IF( WI( I ).GT.ZERO ) THEN 00487 SCL = ONE / SLAPY2( SNRM2( N, VL( 1, I ), 1 ), 00488 $ SNRM2( N, VL( 1, I+1 ), 1 ) ) 00489 CALL SSCAL( N, SCL, VL( 1, I ), 1 ) 00490 CALL SSCAL( N, SCL, VL( 1, I+1 ), 1 ) 00491 DO 10 K = 1, N 00492 WORK( K ) = VL( K, I )**2 + VL( K, I+1 )**2 00493 10 CONTINUE 00494 K = ISAMAX( N, WORK, 1 ) 00495 CALL SLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R ) 00496 CALL SROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN ) 00497 VL( K, I+1 ) = ZERO 00498 END IF 00499 20 CONTINUE 00500 END IF 00501 * 00502 IF( WANTVR ) THEN 00503 * 00504 * Undo balancing of right eigenvectors 00505 * 00506 CALL SGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR, 00507 $ IERR ) 00508 * 00509 * Normalize right eigenvectors and make largest component real 00510 * 00511 DO 40 I = 1, N 00512 IF( WI( I ).EQ.ZERO ) THEN 00513 SCL = ONE / SNRM2( N, VR( 1, I ), 1 ) 00514 CALL SSCAL( N, SCL, VR( 1, I ), 1 ) 00515 ELSE IF( WI( I ).GT.ZERO ) THEN 00516 SCL = ONE / SLAPY2( SNRM2( N, VR( 1, I ), 1 ), 00517 $ SNRM2( N, VR( 1, I+1 ), 1 ) ) 00518 CALL SSCAL( N, SCL, VR( 1, I ), 1 ) 00519 CALL SSCAL( N, SCL, VR( 1, I+1 ), 1 ) 00520 DO 30 K = 1, N 00521 WORK( K ) = VR( K, I )**2 + VR( K, I+1 )**2 00522 30 CONTINUE 00523 K = ISAMAX( N, WORK, 1 ) 00524 CALL SLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R ) 00525 CALL SROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN ) 00526 VR( K, I+1 ) = ZERO 00527 END IF 00528 40 CONTINUE 00529 END IF 00530 * 00531 * Undo scaling if necessary 00532 * 00533 50 CONTINUE 00534 IF( SCALEA ) THEN 00535 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ), 00536 $ MAX( N-INFO, 1 ), IERR ) 00537 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ), 00538 $ MAX( N-INFO, 1 ), IERR ) 00539 IF( INFO.EQ.0 ) THEN 00540 IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 ) 00541 $ CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N, 00542 $ IERR ) 00543 ELSE 00544 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N, 00545 $ IERR ) 00546 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N, 00547 $ IERR ) 00548 END IF 00549 END IF 00550 * 00551 WORK( 1 ) = MAXWRK 00552 RETURN 00553 * 00554 * End of SGEEVX 00555 * 00556 END
1.7.3 
