|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE STGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, 00002 $ LDVL, VR, LDVR, MM, M, WORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER HOWMNY, SIDE 00011 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N 00012 * .. 00013 * .. Array Arguments .. 00014 LOGICAL SELECT( * ) 00015 REAL P( LDP, * ), S( LDS, * ), VL( LDVL, * ), 00016 $ VR( LDVR, * ), WORK( * ) 00017 * .. 00018 * 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * STGEVC computes some or all of the right and/or left eigenvectors of 00024 * a pair of real matrices (S,P), where S is a quasi-triangular matrix 00025 * and P is upper triangular. Matrix pairs of this type are produced by 00026 * the generalized Schur factorization of a matrix pair (A,B): 00027 * 00028 * A = Q*S*Z**T, B = Q*P*Z**T 00029 * 00030 * as computed by SGGHRD + SHGEQZ. 00031 * 00032 * The right eigenvector x and the left eigenvector y of (S,P) 00033 * corresponding to an eigenvalue w are defined by: 00034 * 00035 * S*x = w*P*x, (y**H)*S = w*(y**H)*P, 00036 * 00037 * where y**H denotes the conjugate tranpose of y. 00038 * The eigenvalues are not input to this routine, but are computed 00039 * directly from the diagonal blocks of S and P. 00040 * 00041 * This routine returns the matrices X and/or Y of right and left 00042 * eigenvectors of (S,P), or the products Z*X and/or Q*Y, 00043 * where Z and Q are input matrices. 00044 * If Q and Z are the orthogonal factors from the generalized Schur 00045 * factorization of a matrix pair (A,B), then Z*X and Q*Y 00046 * are the matrices of right and left eigenvectors of (A,B). 00047 * 00048 * Arguments 00049 * ========= 00050 * 00051 * SIDE (input) CHARACTER*1 00052 * = 'R': compute right eigenvectors only; 00053 * = 'L': compute left eigenvectors only; 00054 * = 'B': compute both right and left eigenvectors. 00055 * 00056 * HOWMNY (input) CHARACTER*1 00057 * = 'A': compute all right and/or left eigenvectors; 00058 * = 'B': compute all right and/or left eigenvectors, 00059 * backtransformed by the matrices in VR and/or VL; 00060 * = 'S': compute selected right and/or left eigenvectors, 00061 * specified by the logical array SELECT. 00062 * 00063 * SELECT (input) LOGICAL array, dimension (N) 00064 * If HOWMNY='S', SELECT specifies the eigenvectors to be 00065 * computed. If w(j) is a real eigenvalue, the corresponding 00066 * real eigenvector is computed if SELECT(j) is .TRUE.. 00067 * If w(j) and w(j+1) are the real and imaginary parts of a 00068 * complex eigenvalue, the corresponding complex eigenvector 00069 * is computed if either SELECT(j) or SELECT(j+1) is .TRUE., 00070 * and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is 00071 * set to .FALSE.. 00072 * Not referenced if HOWMNY = 'A' or 'B'. 00073 * 00074 * N (input) INTEGER 00075 * The order of the matrices S and P. N >= 0. 00076 * 00077 * S (input) REAL array, dimension (LDS,N) 00078 * The upper quasi-triangular matrix S from a generalized Schur 00079 * factorization, as computed by SHGEQZ. 00080 * 00081 * LDS (input) INTEGER 00082 * The leading dimension of array S. LDS >= max(1,N). 00083 * 00084 * P (input) REAL array, dimension (LDP,N) 00085 * The upper triangular matrix P from a generalized Schur 00086 * factorization, as computed by SHGEQZ. 00087 * 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks 00088 * of S must be in positive diagonal form. 00089 * 00090 * LDP (input) INTEGER 00091 * The leading dimension of array P. LDP >= max(1,N). 00092 * 00093 * VL (input/output) REAL array, dimension (LDVL,MM) 00094 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 00095 * contain an N-by-N matrix Q (usually the orthogonal matrix Q 00096 * of left Schur vectors returned by SHGEQZ). 00097 * On exit, if SIDE = 'L' or 'B', VL contains: 00098 * if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P); 00099 * if HOWMNY = 'B', the matrix Q*Y; 00100 * if HOWMNY = 'S', the left eigenvectors of (S,P) specified by 00101 * SELECT, stored consecutively in the columns of 00102 * VL, in the same order as their eigenvalues. 00103 * 00104 * A complex eigenvector corresponding to a complex eigenvalue 00105 * is stored in two consecutive columns, the first holding the 00106 * real part, and the second the imaginary part. 00107 * 00108 * Not referenced if SIDE = 'R'. 00109 * 00110 * LDVL (input) INTEGER 00111 * The leading dimension of array VL. LDVL >= 1, and if 00112 * SIDE = 'L' or 'B', LDVL >= N. 00113 * 00114 * VR (input/output) REAL array, dimension (LDVR,MM) 00115 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 00116 * contain an N-by-N matrix Z (usually the orthogonal matrix Z 00117 * of right Schur vectors returned by SHGEQZ). 00118 * 00119 * On exit, if SIDE = 'R' or 'B', VR contains: 00120 * if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P); 00121 * if HOWMNY = 'B' or 'b', the matrix Z*X; 00122 * if HOWMNY = 'S' or 's', the right eigenvectors of (S,P) 00123 * specified by SELECT, stored consecutively in the 00124 * columns of VR, in the same order as their 00125 * eigenvalues. 00126 * 00127 * A complex eigenvector corresponding to a complex eigenvalue 00128 * is stored in two consecutive columns, the first holding the 00129 * real part and the second the imaginary part. 00130 * 00131 * Not referenced if SIDE = 'L'. 00132 * 00133 * LDVR (input) INTEGER 00134 * The leading dimension of the array VR. LDVR >= 1, and if 00135 * SIDE = 'R' or 'B', LDVR >= N. 00136 * 00137 * MM (input) INTEGER 00138 * The number of columns in the arrays VL and/or VR. MM >= M. 00139 * 00140 * M (output) INTEGER 00141 * The number of columns in the arrays VL and/or VR actually 00142 * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 00143 * is set to N. Each selected real eigenvector occupies one 00144 * column and each selected complex eigenvector occupies two 00145 * columns. 00146 * 00147 * WORK (workspace) REAL array, dimension (6*N) 00148 * 00149 * INFO (output) INTEGER 00150 * = 0: successful exit. 00151 * < 0: if INFO = -i, the i-th argument had an illegal value. 00152 * > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex 00153 * eigenvalue. 00154 * 00155 * Further Details 00156 * =============== 00157 * 00158 * Allocation of workspace: 00159 * ---------- -- --------- 00160 * 00161 * WORK( j ) = 1-norm of j-th column of A, above the diagonal 00162 * WORK( N+j ) = 1-norm of j-th column of B, above the diagonal 00163 * WORK( 2*N+1:3*N ) = real part of eigenvector 00164 * WORK( 3*N+1:4*N ) = imaginary part of eigenvector 00165 * WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector 00166 * WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector 00167 * 00168 * Rowwise vs. columnwise solution methods: 00169 * ------- -- ---------- -------- ------- 00170 * 00171 * Finding a generalized eigenvector consists basically of solving the 00172 * singular triangular system 00173 * 00174 * (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left) 00175 * 00176 * Consider finding the i-th right eigenvector (assume all eigenvalues 00177 * are real). The equation to be solved is: 00178 * n i 00179 * 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1 00180 * k=j k=j 00181 * 00182 * where C = (A - w B) (The components v(i+1:n) are 0.) 00183 * 00184 * The "rowwise" method is: 00185 * 00186 * (1) v(i) := 1 00187 * for j = i-1,. . .,1: 00188 * i 00189 * (2) compute s = - sum C(j,k) v(k) and 00190 * k=j+1 00191 * 00192 * (3) v(j) := s / C(j,j) 00193 * 00194 * Step 2 is sometimes called the "dot product" step, since it is an 00195 * inner product between the j-th row and the portion of the eigenvector 00196 * that has been computed so far. 00197 * 00198 * The "columnwise" method consists basically in doing the sums 00199 * for all the rows in parallel. As each v(j) is computed, the 00200 * contribution of v(j) times the j-th column of C is added to the 00201 * partial sums. Since FORTRAN arrays are stored columnwise, this has 00202 * the advantage that at each step, the elements of C that are accessed 00203 * are adjacent to one another, whereas with the rowwise method, the 00204 * elements accessed at a step are spaced LDS (and LDP) words apart. 00205 * 00206 * When finding left eigenvectors, the matrix in question is the 00207 * transpose of the one in storage, so the rowwise method then 00208 * actually accesses columns of A and B at each step, and so is the 00209 * preferred method. 00210 * 00211 * ===================================================================== 00212 * 00213 * .. Parameters .. 00214 REAL ZERO, ONE, SAFETY 00215 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, 00216 $ SAFETY = 1.0E+2 ) 00217 * .. 00218 * .. Local Scalars .. 00219 LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK, 00220 $ ILBBAD, ILCOMP, ILCPLX, LSA, LSB 00221 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE, 00222 $ J, JA, JC, JE, JR, JW, NA, NW 00223 REAL ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI, 00224 $ BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A, 00225 $ CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA, 00226 $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE, 00227 $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX, 00228 $ XSCALE 00229 * .. 00230 * .. Local Arrays .. 00231 REAL BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ), 00232 $ SUMP( 2, 2 ) 00233 * .. 00234 * .. External Functions .. 00235 LOGICAL LSAME 00236 REAL SLAMCH 00237 EXTERNAL LSAME, SLAMCH 00238 * .. 00239 * .. External Subroutines .. 00240 EXTERNAL SGEMV, SLABAD, SLACPY, SLAG2, SLALN2, XERBLA 00241 * .. 00242 * .. Intrinsic Functions .. 00243 INTRINSIC ABS, MAX, MIN 00244 * .. 00245 * .. Executable Statements .. 00246 * 00247 * Decode and Test the input parameters 00248 * 00249 IF( LSAME( HOWMNY, 'A' ) ) THEN 00250 IHWMNY = 1 00251 ILALL = .TRUE. 00252 ILBACK = .FALSE. 00253 ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN 00254 IHWMNY = 2 00255 ILALL = .FALSE. 00256 ILBACK = .FALSE. 00257 ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN 00258 IHWMNY = 3 00259 ILALL = .TRUE. 00260 ILBACK = .TRUE. 00261 ELSE 00262 IHWMNY = -1 00263 ILALL = .TRUE. 00264 END IF 00265 * 00266 IF( LSAME( SIDE, 'R' ) ) THEN 00267 ISIDE = 1 00268 COMPL = .FALSE. 00269 COMPR = .TRUE. 00270 ELSE IF( LSAME( SIDE, 'L' ) ) THEN 00271 ISIDE = 2 00272 COMPL = .TRUE. 00273 COMPR = .FALSE. 00274 ELSE IF( LSAME( SIDE, 'B' ) ) THEN 00275 ISIDE = 3 00276 COMPL = .TRUE. 00277 COMPR = .TRUE. 00278 ELSE 00279 ISIDE = -1 00280 END IF 00281 * 00282 INFO = 0 00283 IF( ISIDE.LT.0 ) THEN 00284 INFO = -1 00285 ELSE IF( IHWMNY.LT.0 ) THEN 00286 INFO = -2 00287 ELSE IF( N.LT.0 ) THEN 00288 INFO = -4 00289 ELSE IF( LDS.LT.MAX( 1, N ) ) THEN 00290 INFO = -6 00291 ELSE IF( LDP.LT.MAX( 1, N ) ) THEN 00292 INFO = -8 00293 END IF 00294 IF( INFO.NE.0 ) THEN 00295 CALL XERBLA( 'STGEVC', -INFO ) 00296 RETURN 00297 END IF 00298 * 00299 * Count the number of eigenvectors to be computed 00300 * 00301 IF( .NOT.ILALL ) THEN 00302 IM = 0 00303 ILCPLX = .FALSE. 00304 DO 10 J = 1, N 00305 IF( ILCPLX ) THEN 00306 ILCPLX = .FALSE. 00307 GO TO 10 00308 END IF 00309 IF( J.LT.N ) THEN 00310 IF( S( J+1, J ).NE.ZERO ) 00311 $ ILCPLX = .TRUE. 00312 END IF 00313 IF( ILCPLX ) THEN 00314 IF( SELECT( J ) .OR. SELECT( J+1 ) ) 00315 $ IM = IM + 2 00316 ELSE 00317 IF( SELECT( J ) ) 00318 $ IM = IM + 1 00319 END IF 00320 10 CONTINUE 00321 ELSE 00322 IM = N 00323 END IF 00324 * 00325 * Check 2-by-2 diagonal blocks of A, B 00326 * 00327 ILABAD = .FALSE. 00328 ILBBAD = .FALSE. 00329 DO 20 J = 1, N - 1 00330 IF( S( J+1, J ).NE.ZERO ) THEN 00331 IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR. 00332 $ P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE. 00333 IF( J.LT.N-1 ) THEN 00334 IF( S( J+2, J+1 ).NE.ZERO ) 00335 $ ILABAD = .TRUE. 00336 END IF 00337 END IF 00338 20 CONTINUE 00339 * 00340 IF( ILABAD ) THEN 00341 INFO = -5 00342 ELSE IF( ILBBAD ) THEN 00343 INFO = -7 00344 ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN 00345 INFO = -10 00346 ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN 00347 INFO = -12 00348 ELSE IF( MM.LT.IM ) THEN 00349 INFO = -13 00350 END IF 00351 IF( INFO.NE.0 ) THEN 00352 CALL XERBLA( 'STGEVC', -INFO ) 00353 RETURN 00354 END IF 00355 * 00356 * Quick return if possible 00357 * 00358 M = IM 00359 IF( N.EQ.0 ) 00360 $ RETURN 00361 * 00362 * Machine Constants 00363 * 00364 SAFMIN = SLAMCH( 'Safe minimum' ) 00365 BIG = ONE / SAFMIN 00366 CALL SLABAD( SAFMIN, BIG ) 00367 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00368 SMALL = SAFMIN*N / ULP 00369 BIG = ONE / SMALL 00370 BIGNUM = ONE / ( SAFMIN*N ) 00371 * 00372 * Compute the 1-norm of each column of the strictly upper triangular 00373 * part (i.e., excluding all elements belonging to the diagonal 00374 * blocks) of A and B to check for possible overflow in the 00375 * triangular solver. 00376 * 00377 ANORM = ABS( S( 1, 1 ) ) 00378 IF( N.GT.1 ) 00379 $ ANORM = ANORM + ABS( S( 2, 1 ) ) 00380 BNORM = ABS( P( 1, 1 ) ) 00381 WORK( 1 ) = ZERO 00382 WORK( N+1 ) = ZERO 00383 * 00384 DO 50 J = 2, N 00385 TEMP = ZERO 00386 TEMP2 = ZERO 00387 IF( S( J, J-1 ).EQ.ZERO ) THEN 00388 IEND = J - 1 00389 ELSE 00390 IEND = J - 2 00391 END IF 00392 DO 30 I = 1, IEND 00393 TEMP = TEMP + ABS( S( I, J ) ) 00394 TEMP2 = TEMP2 + ABS( P( I, J ) ) 00395 30 CONTINUE 00396 WORK( J ) = TEMP 00397 WORK( N+J ) = TEMP2 00398 DO 40 I = IEND + 1, MIN( J+1, N ) 00399 TEMP = TEMP + ABS( S( I, J ) ) 00400 TEMP2 = TEMP2 + ABS( P( I, J ) ) 00401 40 CONTINUE 00402 ANORM = MAX( ANORM, TEMP ) 00403 BNORM = MAX( BNORM, TEMP2 ) 00404 50 CONTINUE 00405 * 00406 ASCALE = ONE / MAX( ANORM, SAFMIN ) 00407 BSCALE = ONE / MAX( BNORM, SAFMIN ) 00408 * 00409 * Left eigenvectors 00410 * 00411 IF( COMPL ) THEN 00412 IEIG = 0 00413 * 00414 * Main loop over eigenvalues 00415 * 00416 ILCPLX = .FALSE. 00417 DO 220 JE = 1, N 00418 * 00419 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 00420 * (b) this would be the second of a complex pair. 00421 * Check for complex eigenvalue, so as to be sure of which 00422 * entry(-ies) of SELECT to look at. 00423 * 00424 IF( ILCPLX ) THEN 00425 ILCPLX = .FALSE. 00426 GO TO 220 00427 END IF 00428 NW = 1 00429 IF( JE.LT.N ) THEN 00430 IF( S( JE+1, JE ).NE.ZERO ) THEN 00431 ILCPLX = .TRUE. 00432 NW = 2 00433 END IF 00434 END IF 00435 IF( ILALL ) THEN 00436 ILCOMP = .TRUE. 00437 ELSE IF( ILCPLX ) THEN 00438 ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 ) 00439 ELSE 00440 ILCOMP = SELECT( JE ) 00441 END IF 00442 IF( .NOT.ILCOMP ) 00443 $ GO TO 220 00444 * 00445 * Decide if (a) singular pencil, (b) real eigenvalue, or 00446 * (c) complex eigenvalue. 00447 * 00448 IF( .NOT.ILCPLX ) THEN 00449 IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND. 00450 $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN 00451 * 00452 * Singular matrix pencil -- return unit eigenvector 00453 * 00454 IEIG = IEIG + 1 00455 DO 60 JR = 1, N 00456 VL( JR, IEIG ) = ZERO 00457 60 CONTINUE 00458 VL( IEIG, IEIG ) = ONE 00459 GO TO 220 00460 END IF 00461 END IF 00462 * 00463 * Clear vector 00464 * 00465 DO 70 JR = 1, NW*N 00466 WORK( 2*N+JR ) = ZERO 00467 70 CONTINUE 00468 * T 00469 * Compute coefficients in ( a A - b B ) y = 0 00470 * a is ACOEF 00471 * b is BCOEFR + i*BCOEFI 00472 * 00473 IF( .NOT.ILCPLX ) THEN 00474 * 00475 * Real eigenvalue 00476 * 00477 TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE, 00478 $ ABS( P( JE, JE ) )*BSCALE, SAFMIN ) 00479 SALFAR = ( TEMP*S( JE, JE ) )*ASCALE 00480 SBETA = ( TEMP*P( JE, JE ) )*BSCALE 00481 ACOEF = SBETA*ASCALE 00482 BCOEFR = SALFAR*BSCALE 00483 BCOEFI = ZERO 00484 * 00485 * Scale to avoid underflow 00486 * 00487 SCALE = ONE 00488 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL 00489 LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. 00490 $ SMALL 00491 IF( LSA ) 00492 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 00493 IF( LSB ) 00494 $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* 00495 $ MIN( BNORM, BIG ) ) 00496 IF( LSA .OR. LSB ) THEN 00497 SCALE = MIN( SCALE, ONE / 00498 $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), 00499 $ ABS( BCOEFR ) ) ) ) 00500 IF( LSA ) THEN 00501 ACOEF = ASCALE*( SCALE*SBETA ) 00502 ELSE 00503 ACOEF = SCALE*ACOEF 00504 END IF 00505 IF( LSB ) THEN 00506 BCOEFR = BSCALE*( SCALE*SALFAR ) 00507 ELSE 00508 BCOEFR = SCALE*BCOEFR 00509 END IF 00510 END IF 00511 ACOEFA = ABS( ACOEF ) 00512 BCOEFA = ABS( BCOEFR ) 00513 * 00514 * First component is 1 00515 * 00516 WORK( 2*N+JE ) = ONE 00517 XMAX = ONE 00518 ELSE 00519 * 00520 * Complex eigenvalue 00521 * 00522 CALL SLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP, 00523 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 00524 $ BCOEFI ) 00525 BCOEFI = -BCOEFI 00526 IF( BCOEFI.EQ.ZERO ) THEN 00527 INFO = JE 00528 RETURN 00529 END IF 00530 * 00531 * Scale to avoid over/underflow 00532 * 00533 ACOEFA = ABS( ACOEF ) 00534 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00535 SCALE = ONE 00536 IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) 00537 $ SCALE = ( SAFMIN / ULP ) / ACOEFA 00538 IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) 00539 $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) 00540 IF( SAFMIN*ACOEFA.GT.ASCALE ) 00541 $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) 00542 IF( SAFMIN*BCOEFA.GT.BSCALE ) 00543 $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) 00544 IF( SCALE.NE.ONE ) THEN 00545 ACOEF = SCALE*ACOEF 00546 ACOEFA = ABS( ACOEF ) 00547 BCOEFR = SCALE*BCOEFR 00548 BCOEFI = SCALE*BCOEFI 00549 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00550 END IF 00551 * 00552 * Compute first two components of eigenvector 00553 * 00554 TEMP = ACOEF*S( JE+1, JE ) 00555 TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE ) 00556 TEMP2I = -BCOEFI*P( JE, JE ) 00557 IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 00558 WORK( 2*N+JE ) = ONE 00559 WORK( 3*N+JE ) = ZERO 00560 WORK( 2*N+JE+1 ) = -TEMP2R / TEMP 00561 WORK( 3*N+JE+1 ) = -TEMP2I / TEMP 00562 ELSE 00563 WORK( 2*N+JE+1 ) = ONE 00564 WORK( 3*N+JE+1 ) = ZERO 00565 TEMP = ACOEF*S( JE, JE+1 ) 00566 WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF* 00567 $ S( JE+1, JE+1 ) ) / TEMP 00568 WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP 00569 END IF 00570 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 00571 $ ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) ) 00572 END IF 00573 * 00574 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 00575 * 00576 * T 00577 * Triangular solve of (a A - b B) y = 0 00578 * 00579 * T 00580 * (rowwise in (a A - b B) , or columnwise in (a A - b B) ) 00581 * 00582 IL2BY2 = .FALSE. 00583 * 00584 DO 160 J = JE + NW, N 00585 IF( IL2BY2 ) THEN 00586 IL2BY2 = .FALSE. 00587 GO TO 160 00588 END IF 00589 * 00590 NA = 1 00591 BDIAG( 1 ) = P( J, J ) 00592 IF( J.LT.N ) THEN 00593 IF( S( J+1, J ).NE.ZERO ) THEN 00594 IL2BY2 = .TRUE. 00595 BDIAG( 2 ) = P( J+1, J+1 ) 00596 NA = 2 00597 END IF 00598 END IF 00599 * 00600 * Check whether scaling is necessary for dot products 00601 * 00602 XSCALE = ONE / MAX( ONE, XMAX ) 00603 TEMP = MAX( WORK( J ), WORK( N+J ), 00604 $ ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) ) 00605 IF( IL2BY2 ) 00606 $ TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ), 00607 $ ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) ) 00608 IF( TEMP.GT.BIGNUM*XSCALE ) THEN 00609 DO 90 JW = 0, NW - 1 00610 DO 80 JR = JE, J - 1 00611 WORK( ( JW+2 )*N+JR ) = XSCALE* 00612 $ WORK( ( JW+2 )*N+JR ) 00613 80 CONTINUE 00614 90 CONTINUE 00615 XMAX = XMAX*XSCALE 00616 END IF 00617 * 00618 * Compute dot products 00619 * 00620 * j-1 00621 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k) 00622 * k=je 00623 * 00624 * To reduce the op count, this is done as 00625 * 00626 * _ j-1 _ j-1 00627 * a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) ) 00628 * k=je k=je 00629 * 00630 * which may cause underflow problems if A or B are close 00631 * to underflow. (E.g., less than SMALL.) 00632 * 00633 * 00634 * A series of compiler directives to defeat vectorization 00635 * for the next loop 00636 * 00637 *$PL$ CMCHAR=' ' 00638 CDIR$ NEXTSCALAR 00639 C$DIR SCALAR 00640 CDIR$ NEXT SCALAR 00641 CVD$L NOVECTOR 00642 CDEC$ NOVECTOR 00643 CVD$ NOVECTOR 00644 *VDIR NOVECTOR 00645 *VOCL LOOP,SCALAR 00646 CIBM PREFER SCALAR 00647 *$PL$ CMCHAR='*' 00648 * 00649 DO 120 JW = 1, NW 00650 * 00651 *$PL$ CMCHAR=' ' 00652 CDIR$ NEXTSCALAR 00653 C$DIR SCALAR 00654 CDIR$ NEXT SCALAR 00655 CVD$L NOVECTOR 00656 CDEC$ NOVECTOR 00657 CVD$ NOVECTOR 00658 *VDIR NOVECTOR 00659 *VOCL LOOP,SCALAR 00660 CIBM PREFER SCALAR 00661 *$PL$ CMCHAR='*' 00662 * 00663 DO 110 JA = 1, NA 00664 SUMS( JA, JW ) = ZERO 00665 SUMP( JA, JW ) = ZERO 00666 * 00667 DO 100 JR = JE, J - 1 00668 SUMS( JA, JW ) = SUMS( JA, JW ) + 00669 $ S( JR, J+JA-1 )* 00670 $ WORK( ( JW+1 )*N+JR ) 00671 SUMP( JA, JW ) = SUMP( JA, JW ) + 00672 $ P( JR, J+JA-1 )* 00673 $ WORK( ( JW+1 )*N+JR ) 00674 100 CONTINUE 00675 110 CONTINUE 00676 120 CONTINUE 00677 * 00678 *$PL$ CMCHAR=' ' 00679 CDIR$ NEXTSCALAR 00680 C$DIR SCALAR 00681 CDIR$ NEXT SCALAR 00682 CVD$L NOVECTOR 00683 CDEC$ NOVECTOR 00684 CVD$ NOVECTOR 00685 *VDIR NOVECTOR 00686 *VOCL LOOP,SCALAR 00687 CIBM PREFER SCALAR 00688 *$PL$ CMCHAR='*' 00689 * 00690 DO 130 JA = 1, NA 00691 IF( ILCPLX ) THEN 00692 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 00693 $ BCOEFR*SUMP( JA, 1 ) - 00694 $ BCOEFI*SUMP( JA, 2 ) 00695 SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) + 00696 $ BCOEFR*SUMP( JA, 2 ) + 00697 $ BCOEFI*SUMP( JA, 1 ) 00698 ELSE 00699 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 00700 $ BCOEFR*SUMP( JA, 1 ) 00701 END IF 00702 130 CONTINUE 00703 * 00704 * T 00705 * Solve ( a A - b B ) y = SUM(,) 00706 * with scaling and perturbation of the denominator 00707 * 00708 CALL SLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS, 00709 $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR, 00710 $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP, 00711 $ IINFO ) 00712 IF( SCALE.LT.ONE ) THEN 00713 DO 150 JW = 0, NW - 1 00714 DO 140 JR = JE, J - 1 00715 WORK( ( JW+2 )*N+JR ) = SCALE* 00716 $ WORK( ( JW+2 )*N+JR ) 00717 140 CONTINUE 00718 150 CONTINUE 00719 XMAX = SCALE*XMAX 00720 END IF 00721 XMAX = MAX( XMAX, TEMP ) 00722 160 CONTINUE 00723 * 00724 * Copy eigenvector to VL, back transforming if 00725 * HOWMNY='B'. 00726 * 00727 IEIG = IEIG + 1 00728 IF( ILBACK ) THEN 00729 DO 170 JW = 0, NW - 1 00730 CALL SGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL, 00731 $ WORK( ( JW+2 )*N+JE ), 1, ZERO, 00732 $ WORK( ( JW+4 )*N+1 ), 1 ) 00733 170 CONTINUE 00734 CALL SLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ), 00735 $ LDVL ) 00736 IBEG = 1 00737 ELSE 00738 CALL SLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ), 00739 $ LDVL ) 00740 IBEG = JE 00741 END IF 00742 * 00743 * Scale eigenvector 00744 * 00745 XMAX = ZERO 00746 IF( ILCPLX ) THEN 00747 DO 180 J = IBEG, N 00748 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+ 00749 $ ABS( VL( J, IEIG+1 ) ) ) 00750 180 CONTINUE 00751 ELSE 00752 DO 190 J = IBEG, N 00753 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) ) 00754 190 CONTINUE 00755 END IF 00756 * 00757 IF( XMAX.GT.SAFMIN ) THEN 00758 XSCALE = ONE / XMAX 00759 * 00760 DO 210 JW = 0, NW - 1 00761 DO 200 JR = IBEG, N 00762 VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW ) 00763 200 CONTINUE 00764 210 CONTINUE 00765 END IF 00766 IEIG = IEIG + NW - 1 00767 * 00768 220 CONTINUE 00769 END IF 00770 * 00771 * Right eigenvectors 00772 * 00773 IF( COMPR ) THEN 00774 IEIG = IM + 1 00775 * 00776 * Main loop over eigenvalues 00777 * 00778 ILCPLX = .FALSE. 00779 DO 500 JE = N, 1, -1 00780 * 00781 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 00782 * (b) this would be the second of a complex pair. 00783 * Check for complex eigenvalue, so as to be sure of which 00784 * entry(-ies) of SELECT to look at -- if complex, SELECT(JE) 00785 * or SELECT(JE-1). 00786 * If this is a complex pair, the 2-by-2 diagonal block 00787 * corresponding to the eigenvalue is in rows/columns JE-1:JE 00788 * 00789 IF( ILCPLX ) THEN 00790 ILCPLX = .FALSE. 00791 GO TO 500 00792 END IF 00793 NW = 1 00794 IF( JE.GT.1 ) THEN 00795 IF( S( JE, JE-1 ).NE.ZERO ) THEN 00796 ILCPLX = .TRUE. 00797 NW = 2 00798 END IF 00799 END IF 00800 IF( ILALL ) THEN 00801 ILCOMP = .TRUE. 00802 ELSE IF( ILCPLX ) THEN 00803 ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 ) 00804 ELSE 00805 ILCOMP = SELECT( JE ) 00806 END IF 00807 IF( .NOT.ILCOMP ) 00808 $ GO TO 500 00809 * 00810 * Decide if (a) singular pencil, (b) real eigenvalue, or 00811 * (c) complex eigenvalue. 00812 * 00813 IF( .NOT.ILCPLX ) THEN 00814 IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND. 00815 $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN 00816 * 00817 * Singular matrix pencil -- unit eigenvector 00818 * 00819 IEIG = IEIG - 1 00820 DO 230 JR = 1, N 00821 VR( JR, IEIG ) = ZERO 00822 230 CONTINUE 00823 VR( IEIG, IEIG ) = ONE 00824 GO TO 500 00825 END IF 00826 END IF 00827 * 00828 * Clear vector 00829 * 00830 DO 250 JW = 0, NW - 1 00831 DO 240 JR = 1, N 00832 WORK( ( JW+2 )*N+JR ) = ZERO 00833 240 CONTINUE 00834 250 CONTINUE 00835 * 00836 * Compute coefficients in ( a A - b B ) x = 0 00837 * a is ACOEF 00838 * b is BCOEFR + i*BCOEFI 00839 * 00840 IF( .NOT.ILCPLX ) THEN 00841 * 00842 * Real eigenvalue 00843 * 00844 TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE, 00845 $ ABS( P( JE, JE ) )*BSCALE, SAFMIN ) 00846 SALFAR = ( TEMP*S( JE, JE ) )*ASCALE 00847 SBETA = ( TEMP*P( JE, JE ) )*BSCALE 00848 ACOEF = SBETA*ASCALE 00849 BCOEFR = SALFAR*BSCALE 00850 BCOEFI = ZERO 00851 * 00852 * Scale to avoid underflow 00853 * 00854 SCALE = ONE 00855 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL 00856 LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. 00857 $ SMALL 00858 IF( LSA ) 00859 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 00860 IF( LSB ) 00861 $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* 00862 $ MIN( BNORM, BIG ) ) 00863 IF( LSA .OR. LSB ) THEN 00864 SCALE = MIN( SCALE, ONE / 00865 $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), 00866 $ ABS( BCOEFR ) ) ) ) 00867 IF( LSA ) THEN 00868 ACOEF = ASCALE*( SCALE*SBETA ) 00869 ELSE 00870 ACOEF = SCALE*ACOEF 00871 END IF 00872 IF( LSB ) THEN 00873 BCOEFR = BSCALE*( SCALE*SALFAR ) 00874 ELSE 00875 BCOEFR = SCALE*BCOEFR 00876 END IF 00877 END IF 00878 ACOEFA = ABS( ACOEF ) 00879 BCOEFA = ABS( BCOEFR ) 00880 * 00881 * First component is 1 00882 * 00883 WORK( 2*N+JE ) = ONE 00884 XMAX = ONE 00885 * 00886 * Compute contribution from column JE of A and B to sum 00887 * (See "Further Details", above.) 00888 * 00889 DO 260 JR = 1, JE - 1 00890 WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) - 00891 $ ACOEF*S( JR, JE ) 00892 260 CONTINUE 00893 ELSE 00894 * 00895 * Complex eigenvalue 00896 * 00897 CALL SLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP, 00898 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 00899 $ BCOEFI ) 00900 IF( BCOEFI.EQ.ZERO ) THEN 00901 INFO = JE - 1 00902 RETURN 00903 END IF 00904 * 00905 * Scale to avoid over/underflow 00906 * 00907 ACOEFA = ABS( ACOEF ) 00908 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00909 SCALE = ONE 00910 IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) 00911 $ SCALE = ( SAFMIN / ULP ) / ACOEFA 00912 IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) 00913 $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) 00914 IF( SAFMIN*ACOEFA.GT.ASCALE ) 00915 $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) 00916 IF( SAFMIN*BCOEFA.GT.BSCALE ) 00917 $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) 00918 IF( SCALE.NE.ONE ) THEN 00919 ACOEF = SCALE*ACOEF 00920 ACOEFA = ABS( ACOEF ) 00921 BCOEFR = SCALE*BCOEFR 00922 BCOEFI = SCALE*BCOEFI 00923 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00924 END IF 00925 * 00926 * Compute first two components of eigenvector 00927 * and contribution to sums 00928 * 00929 TEMP = ACOEF*S( JE, JE-1 ) 00930 TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE ) 00931 TEMP2I = -BCOEFI*P( JE, JE ) 00932 IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 00933 WORK( 2*N+JE ) = ONE 00934 WORK( 3*N+JE ) = ZERO 00935 WORK( 2*N+JE-1 ) = -TEMP2R / TEMP 00936 WORK( 3*N+JE-1 ) = -TEMP2I / TEMP 00937 ELSE 00938 WORK( 2*N+JE-1 ) = ONE 00939 WORK( 3*N+JE-1 ) = ZERO 00940 TEMP = ACOEF*S( JE-1, JE ) 00941 WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF* 00942 $ S( JE-1, JE-1 ) ) / TEMP 00943 WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP 00944 END IF 00945 * 00946 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 00947 $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) ) 00948 * 00949 * Compute contribution from columns JE and JE-1 00950 * of A and B to the sums. 00951 * 00952 CREALA = ACOEF*WORK( 2*N+JE-1 ) 00953 CIMAGA = ACOEF*WORK( 3*N+JE-1 ) 00954 CREALB = BCOEFR*WORK( 2*N+JE-1 ) - 00955 $ BCOEFI*WORK( 3*N+JE-1 ) 00956 CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) + 00957 $ BCOEFR*WORK( 3*N+JE-1 ) 00958 CRE2A = ACOEF*WORK( 2*N+JE ) 00959 CIM2A = ACOEF*WORK( 3*N+JE ) 00960 CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE ) 00961 CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE ) 00962 DO 270 JR = 1, JE - 2 00963 WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) + 00964 $ CREALB*P( JR, JE-1 ) - 00965 $ CRE2A*S( JR, JE ) + CRE2B*P( JR, JE ) 00966 WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) + 00967 $ CIMAGB*P( JR, JE-1 ) - 00968 $ CIM2A*S( JR, JE ) + CIM2B*P( JR, JE ) 00969 270 CONTINUE 00970 END IF 00971 * 00972 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 00973 * 00974 * Columnwise triangular solve of (a A - b B) x = 0 00975 * 00976 IL2BY2 = .FALSE. 00977 DO 370 J = JE - NW, 1, -1 00978 * 00979 * If a 2-by-2 block, is in position j-1:j, wait until 00980 * next iteration to process it (when it will be j:j+1) 00981 * 00982 IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN 00983 IF( S( J, J-1 ).NE.ZERO ) THEN 00984 IL2BY2 = .TRUE. 00985 GO TO 370 00986 END IF 00987 END IF 00988 BDIAG( 1 ) = P( J, J ) 00989 IF( IL2BY2 ) THEN 00990 NA = 2 00991 BDIAG( 2 ) = P( J+1, J+1 ) 00992 ELSE 00993 NA = 1 00994 END IF 00995 * 00996 * Compute x(j) (and x(j+1), if 2-by-2 block) 00997 * 00998 CALL SLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ), 00999 $ LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ), 01000 $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP, 01001 $ IINFO ) 01002 IF( SCALE.LT.ONE ) THEN 01003 * 01004 DO 290 JW = 0, NW - 1 01005 DO 280 JR = 1, JE 01006 WORK( ( JW+2 )*N+JR ) = SCALE* 01007 $ WORK( ( JW+2 )*N+JR ) 01008 280 CONTINUE 01009 290 CONTINUE 01010 END IF 01011 XMAX = MAX( SCALE*XMAX, TEMP ) 01012 * 01013 DO 310 JW = 1, NW 01014 DO 300 JA = 1, NA 01015 WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW ) 01016 300 CONTINUE 01017 310 CONTINUE 01018 * 01019 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling 01020 * 01021 IF( J.GT.1 ) THEN 01022 * 01023 * Check whether scaling is necessary for sum. 01024 * 01025 XSCALE = ONE / MAX( ONE, XMAX ) 01026 TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J ) 01027 IF( IL2BY2 ) 01028 $ TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA* 01029 $ WORK( N+J+1 ) ) 01030 TEMP = MAX( TEMP, ACOEFA, BCOEFA ) 01031 IF( TEMP.GT.BIGNUM*XSCALE ) THEN 01032 * 01033 DO 330 JW = 0, NW - 1 01034 DO 320 JR = 1, JE 01035 WORK( ( JW+2 )*N+JR ) = XSCALE* 01036 $ WORK( ( JW+2 )*N+JR ) 01037 320 CONTINUE 01038 330 CONTINUE 01039 XMAX = XMAX*XSCALE 01040 END IF 01041 * 01042 * Compute the contributions of the off-diagonals of 01043 * column j (and j+1, if 2-by-2 block) of A and B to the 01044 * sums. 01045 * 01046 * 01047 DO 360 JA = 1, NA 01048 IF( ILCPLX ) THEN 01049 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 01050 CIMAGA = ACOEF*WORK( 3*N+J+JA-1 ) 01051 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) - 01052 $ BCOEFI*WORK( 3*N+J+JA-1 ) 01053 CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) + 01054 $ BCOEFR*WORK( 3*N+J+JA-1 ) 01055 DO 340 JR = 1, J - 1 01056 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 01057 $ CREALA*S( JR, J+JA-1 ) + 01058 $ CREALB*P( JR, J+JA-1 ) 01059 WORK( 3*N+JR ) = WORK( 3*N+JR ) - 01060 $ CIMAGA*S( JR, J+JA-1 ) + 01061 $ CIMAGB*P( JR, J+JA-1 ) 01062 340 CONTINUE 01063 ELSE 01064 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 01065 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) 01066 DO 350 JR = 1, J - 1 01067 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 01068 $ CREALA*S( JR, J+JA-1 ) + 01069 $ CREALB*P( JR, J+JA-1 ) 01070 350 CONTINUE 01071 END IF 01072 360 CONTINUE 01073 END IF 01074 * 01075 IL2BY2 = .FALSE. 01076 370 CONTINUE 01077 * 01078 * Copy eigenvector to VR, back transforming if 01079 * HOWMNY='B'. 01080 * 01081 IEIG = IEIG - NW 01082 IF( ILBACK ) THEN 01083 * 01084 DO 410 JW = 0, NW - 1 01085 DO 380 JR = 1, N 01086 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )* 01087 $ VR( JR, 1 ) 01088 380 CONTINUE 01089 * 01090 * A series of compiler directives to defeat 01091 * vectorization for the next loop 01092 * 01093 * 01094 DO 400 JC = 2, JE 01095 DO 390 JR = 1, N 01096 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) + 01097 $ WORK( ( JW+2 )*N+JC )*VR( JR, JC ) 01098 390 CONTINUE 01099 400 CONTINUE 01100 410 CONTINUE 01101 * 01102 DO 430 JW = 0, NW - 1 01103 DO 420 JR = 1, N 01104 VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR ) 01105 420 CONTINUE 01106 430 CONTINUE 01107 * 01108 IEND = N 01109 ELSE 01110 DO 450 JW = 0, NW - 1 01111 DO 440 JR = 1, N 01112 VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR ) 01113 440 CONTINUE 01114 450 CONTINUE 01115 * 01116 IEND = JE 01117 END IF 01118 * 01119 * Scale eigenvector 01120 * 01121 XMAX = ZERO 01122 IF( ILCPLX ) THEN 01123 DO 460 J = 1, IEND 01124 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+ 01125 $ ABS( VR( J, IEIG+1 ) ) ) 01126 460 CONTINUE 01127 ELSE 01128 DO 470 J = 1, IEND 01129 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) ) 01130 470 CONTINUE 01131 END IF 01132 * 01133 IF( XMAX.GT.SAFMIN ) THEN 01134 XSCALE = ONE / XMAX 01135 DO 490 JW = 0, NW - 1 01136 DO 480 JR = 1, IEND 01137 VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW ) 01138 480 CONTINUE 01139 490 CONTINUE 01140 END IF 01141 500 CONTINUE 01142 END IF 01143 * 01144 RETURN 01145 * 01146 * End of STGEVC 01147 * 01148 END
1.7.3 
