 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, 00002 $ LDVL, VR, LDVR, MM, M, WORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 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 DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ), 00016 $ VR( LDVR, * ), WORK( * ) 00017 * .. 00018 * 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * DTGEVC 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 DGGHRD + DHGEQZ. 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) DOUBLE PRECISION array, dimension (LDS,N) 00078 * The upper quasi-triangular matrix S from a generalized Schur 00079 * factorization, as computed by DHGEQZ. 00080 * 00081 * LDS (input) INTEGER 00082 * The leading dimension of array S. LDS >= max(1,N). 00083 * 00084 * P (input) DOUBLE PRECISION array, dimension (LDP,N) 00085 * The upper triangular matrix P from a generalized Schur 00086 * factorization, as computed by DHGEQZ. 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) DOUBLE PRECISION 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 DHGEQZ). 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) DOUBLE PRECISION 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 DHGEQZ). 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) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE, SAFETY 00215 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, 00216 $ SAFETY = 1.0D+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 DOUBLE PRECISION 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 DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ), 00232 $ SUMP( 2, 2 ) 00233 * .. 00234 * .. External Functions .. 00235 LOGICAL LSAME 00236 DOUBLE PRECISION DLAMCH 00237 EXTERNAL LSAME, DLAMCH 00238 * .. 00239 * .. External Subroutines .. 00240 EXTERNAL DGEMV, DLABAD, DLACPY, DLAG2, DLALN2, 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( 'DTGEVC', -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( 'DTGEVC', -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 = DLAMCH( 'Safe minimum' ) 00365 BIG = ONE / SAFMIN 00366 CALL DLABAD( SAFMIN, BIG ) 00367 ULP = DLAMCH( 'Epsilon' )*DLAMCH( '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 DLAG2( 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 DO 120 JW = 1, NW 00635 DO 110 JA = 1, NA 00636 SUMS( JA, JW ) = ZERO 00637 SUMP( JA, JW ) = ZERO 00638 * 00639 DO 100 JR = JE, J - 1 00640 SUMS( JA, JW ) = SUMS( JA, JW ) + 00641 $ S( JR, J+JA-1 )* 00642 $ WORK( ( JW+1 )*N+JR ) 00643 SUMP( JA, JW ) = SUMP( JA, JW ) + 00644 $ P( JR, J+JA-1 )* 00645 $ WORK( ( JW+1 )*N+JR ) 00646 100 CONTINUE 00647 110 CONTINUE 00648 120 CONTINUE 00649 * 00650 DO 130 JA = 1, NA 00651 IF( ILCPLX ) THEN 00652 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 00653 $ BCOEFR*SUMP( JA, 1 ) - 00654 $ BCOEFI*SUMP( JA, 2 ) 00655 SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) + 00656 $ BCOEFR*SUMP( JA, 2 ) + 00657 $ BCOEFI*SUMP( JA, 1 ) 00658 ELSE 00659 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 00660 $ BCOEFR*SUMP( JA, 1 ) 00661 END IF 00662 130 CONTINUE 00663 * 00664 * T 00665 * Solve ( a A - b B ) y = SUM(,) 00666 * with scaling and perturbation of the denominator 00667 * 00668 CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS, 00669 $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR, 00670 $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP, 00671 $ IINFO ) 00672 IF( SCALE.LT.ONE ) THEN 00673 DO 150 JW = 0, NW - 1 00674 DO 140 JR = JE, J - 1 00675 WORK( ( JW+2 )*N+JR ) = SCALE* 00676 $ WORK( ( JW+2 )*N+JR ) 00677 140 CONTINUE 00678 150 CONTINUE 00679 XMAX = SCALE*XMAX 00680 END IF 00681 XMAX = MAX( XMAX, TEMP ) 00682 160 CONTINUE 00683 * 00684 * Copy eigenvector to VL, back transforming if 00685 * HOWMNY='B'. 00686 * 00687 IEIG = IEIG + 1 00688 IF( ILBACK ) THEN 00689 DO 170 JW = 0, NW - 1 00690 CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL, 00691 $ WORK( ( JW+2 )*N+JE ), 1, ZERO, 00692 $ WORK( ( JW+4 )*N+1 ), 1 ) 00693 170 CONTINUE 00694 CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ), 00695 $ LDVL ) 00696 IBEG = 1 00697 ELSE 00698 CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ), 00699 $ LDVL ) 00700 IBEG = JE 00701 END IF 00702 * 00703 * Scale eigenvector 00704 * 00705 XMAX = ZERO 00706 IF( ILCPLX ) THEN 00707 DO 180 J = IBEG, N 00708 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+ 00709 $ ABS( VL( J, IEIG+1 ) ) ) 00710 180 CONTINUE 00711 ELSE 00712 DO 190 J = IBEG, N 00713 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) ) 00714 190 CONTINUE 00715 END IF 00716 * 00717 IF( XMAX.GT.SAFMIN ) THEN 00718 XSCALE = ONE / XMAX 00719 * 00720 DO 210 JW = 0, NW - 1 00721 DO 200 JR = IBEG, N 00722 VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW ) 00723 200 CONTINUE 00724 210 CONTINUE 00725 END IF 00726 IEIG = IEIG + NW - 1 00727 * 00728 220 CONTINUE 00729 END IF 00730 * 00731 * Right eigenvectors 00732 * 00733 IF( COMPR ) THEN 00734 IEIG = IM + 1 00735 * 00736 * Main loop over eigenvalues 00737 * 00738 ILCPLX = .FALSE. 00739 DO 500 JE = N, 1, -1 00740 * 00741 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 00742 * (b) this would be the second of a complex pair. 00743 * Check for complex eigenvalue, so as to be sure of which 00744 * entry(-ies) of SELECT to look at -- if complex, SELECT(JE) 00745 * or SELECT(JE-1). 00746 * If this is a complex pair, the 2-by-2 diagonal block 00747 * corresponding to the eigenvalue is in rows/columns JE-1:JE 00748 * 00749 IF( ILCPLX ) THEN 00750 ILCPLX = .FALSE. 00751 GO TO 500 00752 END IF 00753 NW = 1 00754 IF( JE.GT.1 ) THEN 00755 IF( S( JE, JE-1 ).NE.ZERO ) THEN 00756 ILCPLX = .TRUE. 00757 NW = 2 00758 END IF 00759 END IF 00760 IF( ILALL ) THEN 00761 ILCOMP = .TRUE. 00762 ELSE IF( ILCPLX ) THEN 00763 ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 ) 00764 ELSE 00765 ILCOMP = SELECT( JE ) 00766 END IF 00767 IF( .NOT.ILCOMP ) 00768 $ GO TO 500 00769 * 00770 * Decide if (a) singular pencil, (b) real eigenvalue, or 00771 * (c) complex eigenvalue. 00772 * 00773 IF( .NOT.ILCPLX ) THEN 00774 IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND. 00775 $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN 00776 * 00777 * Singular matrix pencil -- unit eigenvector 00778 * 00779 IEIG = IEIG - 1 00780 DO 230 JR = 1, N 00781 VR( JR, IEIG ) = ZERO 00782 230 CONTINUE 00783 VR( IEIG, IEIG ) = ONE 00784 GO TO 500 00785 END IF 00786 END IF 00787 * 00788 * Clear vector 00789 * 00790 DO 250 JW = 0, NW - 1 00791 DO 240 JR = 1, N 00792 WORK( ( JW+2 )*N+JR ) = ZERO 00793 240 CONTINUE 00794 250 CONTINUE 00795 * 00796 * Compute coefficients in ( a A - b B ) x = 0 00797 * a is ACOEF 00798 * b is BCOEFR + i*BCOEFI 00799 * 00800 IF( .NOT.ILCPLX ) THEN 00801 * 00802 * Real eigenvalue 00803 * 00804 TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE, 00805 $ ABS( P( JE, JE ) )*BSCALE, SAFMIN ) 00806 SALFAR = ( TEMP*S( JE, JE ) )*ASCALE 00807 SBETA = ( TEMP*P( JE, JE ) )*BSCALE 00808 ACOEF = SBETA*ASCALE 00809 BCOEFR = SALFAR*BSCALE 00810 BCOEFI = ZERO 00811 * 00812 * Scale to avoid underflow 00813 * 00814 SCALE = ONE 00815 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL 00816 LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. 00817 $ SMALL 00818 IF( LSA ) 00819 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 00820 IF( LSB ) 00821 $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* 00822 $ MIN( BNORM, BIG ) ) 00823 IF( LSA .OR. LSB ) THEN 00824 SCALE = MIN( SCALE, ONE / 00825 $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), 00826 $ ABS( BCOEFR ) ) ) ) 00827 IF( LSA ) THEN 00828 ACOEF = ASCALE*( SCALE*SBETA ) 00829 ELSE 00830 ACOEF = SCALE*ACOEF 00831 END IF 00832 IF( LSB ) THEN 00833 BCOEFR = BSCALE*( SCALE*SALFAR ) 00834 ELSE 00835 BCOEFR = SCALE*BCOEFR 00836 END IF 00837 END IF 00838 ACOEFA = ABS( ACOEF ) 00839 BCOEFA = ABS( BCOEFR ) 00840 * 00841 * First component is 1 00842 * 00843 WORK( 2*N+JE ) = ONE 00844 XMAX = ONE 00845 * 00846 * Compute contribution from column JE of A and B to sum 00847 * (See "Further Details", above.) 00848 * 00849 DO 260 JR = 1, JE - 1 00850 WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) - 00851 $ ACOEF*S( JR, JE ) 00852 260 CONTINUE 00853 ELSE 00854 * 00855 * Complex eigenvalue 00856 * 00857 CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP, 00858 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 00859 $ BCOEFI ) 00860 IF( BCOEFI.EQ.ZERO ) THEN 00861 INFO = JE - 1 00862 RETURN 00863 END IF 00864 * 00865 * Scale to avoid over/underflow 00866 * 00867 ACOEFA = ABS( ACOEF ) 00868 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00869 SCALE = ONE 00870 IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) 00871 $ SCALE = ( SAFMIN / ULP ) / ACOEFA 00872 IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) 00873 $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) 00874 IF( SAFMIN*ACOEFA.GT.ASCALE ) 00875 $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) 00876 IF( SAFMIN*BCOEFA.GT.BSCALE ) 00877 $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) 00878 IF( SCALE.NE.ONE ) THEN 00879 ACOEF = SCALE*ACOEF 00880 ACOEFA = ABS( ACOEF ) 00881 BCOEFR = SCALE*BCOEFR 00882 BCOEFI = SCALE*BCOEFI 00883 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00884 END IF 00885 * 00886 * Compute first two components of eigenvector 00887 * and contribution to sums 00888 * 00889 TEMP = ACOEF*S( JE, JE-1 ) 00890 TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE ) 00891 TEMP2I = -BCOEFI*P( JE, JE ) 00892 IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 00893 WORK( 2*N+JE ) = ONE 00894 WORK( 3*N+JE ) = ZERO 00895 WORK( 2*N+JE-1 ) = -TEMP2R / TEMP 00896 WORK( 3*N+JE-1 ) = -TEMP2I / TEMP 00897 ELSE 00898 WORK( 2*N+JE-1 ) = ONE 00899 WORK( 3*N+JE-1 ) = ZERO 00900 TEMP = ACOEF*S( JE-1, JE ) 00901 WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF* 00902 $ S( JE-1, JE-1 ) ) / TEMP 00903 WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP 00904 END IF 00905 * 00906 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 00907 $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) ) 00908 * 00909 * Compute contribution from columns JE and JE-1 00910 * of A and B to the sums. 00911 * 00912 CREALA = ACOEF*WORK( 2*N+JE-1 ) 00913 CIMAGA = ACOEF*WORK( 3*N+JE-1 ) 00914 CREALB = BCOEFR*WORK( 2*N+JE-1 ) - 00915 $ BCOEFI*WORK( 3*N+JE-1 ) 00916 CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) + 00917 $ BCOEFR*WORK( 3*N+JE-1 ) 00918 CRE2A = ACOEF*WORK( 2*N+JE ) 00919 CIM2A = ACOEF*WORK( 3*N+JE ) 00920 CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE ) 00921 CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE ) 00922 DO 270 JR = 1, JE - 2 00923 WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) + 00924 $ CREALB*P( JR, JE-1 ) - 00925 $ CRE2A*S( JR, JE ) + CRE2B*P( JR, JE ) 00926 WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) + 00927 $ CIMAGB*P( JR, JE-1 ) - 00928 $ CIM2A*S( JR, JE ) + CIM2B*P( JR, JE ) 00929 270 CONTINUE 00930 END IF 00931 * 00932 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 00933 * 00934 * Columnwise triangular solve of (a A - b B) x = 0 00935 * 00936 IL2BY2 = .FALSE. 00937 DO 370 J = JE - NW, 1, -1 00938 * 00939 * If a 2-by-2 block, is in position j-1:j, wait until 00940 * next iteration to process it (when it will be j:j+1) 00941 * 00942 IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN 00943 IF( S( J, J-1 ).NE.ZERO ) THEN 00944 IL2BY2 = .TRUE. 00945 GO TO 370 00946 END IF 00947 END IF 00948 BDIAG( 1 ) = P( J, J ) 00949 IF( IL2BY2 ) THEN 00950 NA = 2 00951 BDIAG( 2 ) = P( J+1, J+1 ) 00952 ELSE 00953 NA = 1 00954 END IF 00955 * 00956 * Compute x(j) (and x(j+1), if 2-by-2 block) 00957 * 00958 CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ), 00959 $ LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ), 00960 $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP, 00961 $ IINFO ) 00962 IF( SCALE.LT.ONE ) THEN 00963 * 00964 DO 290 JW = 0, NW - 1 00965 DO 280 JR = 1, JE 00966 WORK( ( JW+2 )*N+JR ) = SCALE* 00967 $ WORK( ( JW+2 )*N+JR ) 00968 280 CONTINUE 00969 290 CONTINUE 00970 END IF 00971 XMAX = MAX( SCALE*XMAX, TEMP ) 00972 * 00973 DO 310 JW = 1, NW 00974 DO 300 JA = 1, NA 00975 WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW ) 00976 300 CONTINUE 00977 310 CONTINUE 00978 * 00979 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling 00980 * 00981 IF( J.GT.1 ) THEN 00982 * 00983 * Check whether scaling is necessary for sum. 00984 * 00985 XSCALE = ONE / MAX( ONE, XMAX ) 00986 TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J ) 00987 IF( IL2BY2 ) 00988 $ TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA* 00989 $ WORK( N+J+1 ) ) 00990 TEMP = MAX( TEMP, ACOEFA, BCOEFA ) 00991 IF( TEMP.GT.BIGNUM*XSCALE ) THEN 00992 * 00993 DO 330 JW = 0, NW - 1 00994 DO 320 JR = 1, JE 00995 WORK( ( JW+2 )*N+JR ) = XSCALE* 00996 $ WORK( ( JW+2 )*N+JR ) 00997 320 CONTINUE 00998 330 CONTINUE 00999 XMAX = XMAX*XSCALE 01000 END IF 01001 * 01002 * Compute the contributions of the off-diagonals of 01003 * column j (and j+1, if 2-by-2 block) of A and B to the 01004 * sums. 01005 * 01006 * 01007 DO 360 JA = 1, NA 01008 IF( ILCPLX ) THEN 01009 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 01010 CIMAGA = ACOEF*WORK( 3*N+J+JA-1 ) 01011 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) - 01012 $ BCOEFI*WORK( 3*N+J+JA-1 ) 01013 CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) + 01014 $ BCOEFR*WORK( 3*N+J+JA-1 ) 01015 DO 340 JR = 1, J - 1 01016 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 01017 $ CREALA*S( JR, J+JA-1 ) + 01018 $ CREALB*P( JR, J+JA-1 ) 01019 WORK( 3*N+JR ) = WORK( 3*N+JR ) - 01020 $ CIMAGA*S( JR, J+JA-1 ) + 01021 $ CIMAGB*P( JR, J+JA-1 ) 01022 340 CONTINUE 01023 ELSE 01024 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 01025 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) 01026 DO 350 JR = 1, J - 1 01027 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 01028 $ CREALA*S( JR, J+JA-1 ) + 01029 $ CREALB*P( JR, J+JA-1 ) 01030 350 CONTINUE 01031 END IF 01032 360 CONTINUE 01033 END IF 01034 * 01035 IL2BY2 = .FALSE. 01036 370 CONTINUE 01037 * 01038 * Copy eigenvector to VR, back transforming if 01039 * HOWMNY='B'. 01040 * 01041 IEIG = IEIG - NW 01042 IF( ILBACK ) THEN 01043 * 01044 DO 410 JW = 0, NW - 1 01045 DO 380 JR = 1, N 01046 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )* 01047 $ VR( JR, 1 ) 01048 380 CONTINUE 01049 * 01050 * A series of compiler directives to defeat 01051 * vectorization for the next loop 01052 * 01053 * 01054 DO 400 JC = 2, JE 01055 DO 390 JR = 1, N 01056 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) + 01057 $ WORK( ( JW+2 )*N+JC )*VR( JR, JC ) 01058 390 CONTINUE 01059 400 CONTINUE 01060 410 CONTINUE 01061 * 01062 DO 430 JW = 0, NW - 1 01063 DO 420 JR = 1, N 01064 VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR ) 01065 420 CONTINUE 01066 430 CONTINUE 01067 * 01068 IEND = N 01069 ELSE 01070 DO 450 JW = 0, NW - 1 01071 DO 440 JR = 1, N 01072 VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR ) 01073 440 CONTINUE 01074 450 CONTINUE 01075 * 01076 IEND = JE 01077 END IF 01078 * 01079 * Scale eigenvector 01080 * 01081 XMAX = ZERO 01082 IF( ILCPLX ) THEN 01083 DO 460 J = 1, IEND 01084 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+ 01085 $ ABS( VR( J, IEIG+1 ) ) ) 01086 460 CONTINUE 01087 ELSE 01088 DO 470 J = 1, IEND 01089 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) ) 01090 470 CONTINUE 01091 END IF 01092 * 01093 IF( XMAX.GT.SAFMIN ) THEN 01094 XSCALE = ONE / XMAX 01095 DO 490 JW = 0, NW - 1 01096 DO 480 JR = 1, IEND 01097 VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW ) 01098 480 CONTINUE 01099 490 CONTINUE 01100 END IF 01101 500 CONTINUE 01102 END IF 01103 * 01104 RETURN 01105 * 01106 * End of DTGEVC 01107 * 01108 END
1.7.3 
