![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, 2: $ LDVL, VR, LDVR, MM, M, WORK, INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * .. Scalar Arguments .. 10: CHARACTER HOWMNY, SIDE 11: INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N 12: * .. 13: * .. Array Arguments .. 14: LOGICAL SELECT( * ) 15: DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ), 16: $ VR( LDVR, * ), WORK( * ) 17: * .. 18: * 19: * 20: * Purpose 21: * ======= 22: * 23: * DTGEVC computes some or all of the right and/or left eigenvectors of 24: * a pair of real matrices (S,P), where S is a quasi-triangular matrix 25: * and P is upper triangular. Matrix pairs of this type are produced by 26: * the generalized Schur factorization of a matrix pair (A,B): 27: * 28: * A = Q*S*Z**T, B = Q*P*Z**T 29: * 30: * as computed by DGGHRD + DHGEQZ. 31: * 32: * The right eigenvector x and the left eigenvector y of (S,P) 33: * corresponding to an eigenvalue w are defined by: 34: * 35: * S*x = w*P*x, (y**H)*S = w*(y**H)*P, 36: * 37: * where y**H denotes the conjugate tranpose of y. 38: * The eigenvalues are not input to this routine, but are computed 39: * directly from the diagonal blocks of S and P. 40: * 41: * This routine returns the matrices X and/or Y of right and left 42: * eigenvectors of (S,P), or the products Z*X and/or Q*Y, 43: * where Z and Q are input matrices. 44: * If Q and Z are the orthogonal factors from the generalized Schur 45: * factorization of a matrix pair (A,B), then Z*X and Q*Y 46: * are the matrices of right and left eigenvectors of (A,B). 47: * 48: * Arguments 49: * ========= 50: * 51: * SIDE (input) CHARACTER*1 52: * = 'R': compute right eigenvectors only; 53: * = 'L': compute left eigenvectors only; 54: * = 'B': compute both right and left eigenvectors. 55: * 56: * HOWMNY (input) CHARACTER*1 57: * = 'A': compute all right and/or left eigenvectors; 58: * = 'B': compute all right and/or left eigenvectors, 59: * backtransformed by the matrices in VR and/or VL; 60: * = 'S': compute selected right and/or left eigenvectors, 61: * specified by the logical array SELECT. 62: * 63: * SELECT (input) LOGICAL array, dimension (N) 64: * If HOWMNY='S', SELECT specifies the eigenvectors to be 65: * computed. If w(j) is a real eigenvalue, the corresponding 66: * real eigenvector is computed if SELECT(j) is .TRUE.. 67: * If w(j) and w(j+1) are the real and imaginary parts of a 68: * complex eigenvalue, the corresponding complex eigenvector 69: * is computed if either SELECT(j) or SELECT(j+1) is .TRUE., 70: * and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is 71: * set to .FALSE.. 72: * Not referenced if HOWMNY = 'A' or 'B'. 73: * 74: * N (input) INTEGER 75: * The order of the matrices S and P. N >= 0. 76: * 77: * S (input) DOUBLE PRECISION array, dimension (LDS,N) 78: * The upper quasi-triangular matrix S from a generalized Schur 79: * factorization, as computed by DHGEQZ. 80: * 81: * LDS (input) INTEGER 82: * The leading dimension of array S. LDS >= max(1,N). 83: * 84: * P (input) DOUBLE PRECISION array, dimension (LDP,N) 85: * The upper triangular matrix P from a generalized Schur 86: * factorization, as computed by DHGEQZ. 87: * 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks 88: * of S must be in positive diagonal form. 89: * 90: * LDP (input) INTEGER 91: * The leading dimension of array P. LDP >= max(1,N). 92: * 93: * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM) 94: * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 95: * contain an N-by-N matrix Q (usually the orthogonal matrix Q 96: * of left Schur vectors returned by DHGEQZ). 97: * On exit, if SIDE = 'L' or 'B', VL contains: 98: * if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P); 99: * if HOWMNY = 'B', the matrix Q*Y; 100: * if HOWMNY = 'S', the left eigenvectors of (S,P) specified by 101: * SELECT, stored consecutively in the columns of 102: * VL, in the same order as their eigenvalues. 103: * 104: * A complex eigenvector corresponding to a complex eigenvalue 105: * is stored in two consecutive columns, the first holding the 106: * real part, and the second the imaginary part. 107: * 108: * Not referenced if SIDE = 'R'. 109: * 110: * LDVL (input) INTEGER 111: * The leading dimension of array VL. LDVL >= 1, and if 112: * SIDE = 'L' or 'B', LDVL >= N. 113: * 114: * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) 115: * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 116: * contain an N-by-N matrix Z (usually the orthogonal matrix Z 117: * of right Schur vectors returned by DHGEQZ). 118: * 119: * On exit, if SIDE = 'R' or 'B', VR contains: 120: * if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P); 121: * if HOWMNY = 'B' or 'b', the matrix Z*X; 122: * if HOWMNY = 'S' or 's', the right eigenvectors of (S,P) 123: * specified by SELECT, stored consecutively in the 124: * columns of VR, in the same order as their 125: * eigenvalues. 126: * 127: * A complex eigenvector corresponding to a complex eigenvalue 128: * is stored in two consecutive columns, the first holding the 129: * real part and the second the imaginary part. 130: * 131: * Not referenced if SIDE = 'L'. 132: * 133: * LDVR (input) INTEGER 134: * The leading dimension of the array VR. LDVR >= 1, and if 135: * SIDE = 'R' or 'B', LDVR >= N. 136: * 137: * MM (input) INTEGER 138: * The number of columns in the arrays VL and/or VR. MM >= M. 139: * 140: * M (output) INTEGER 141: * The number of columns in the arrays VL and/or VR actually 142: * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 143: * is set to N. Each selected real eigenvector occupies one 144: * column and each selected complex eigenvector occupies two 145: * columns. 146: * 147: * WORK (workspace) DOUBLE PRECISION array, dimension (6*N) 148: * 149: * INFO (output) INTEGER 150: * = 0: successful exit. 151: * < 0: if INFO = -i, the i-th argument had an illegal value. 152: * > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex 153: * eigenvalue. 154: * 155: * Further Details 156: * =============== 157: * 158: * Allocation of workspace: 159: * ---------- -- --------- 160: * 161: * WORK( j ) = 1-norm of j-th column of A, above the diagonal 162: * WORK( N+j ) = 1-norm of j-th column of B, above the diagonal 163: * WORK( 2*N+1:3*N ) = real part of eigenvector 164: * WORK( 3*N+1:4*N ) = imaginary part of eigenvector 165: * WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector 166: * WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector 167: * 168: * Rowwise vs. columnwise solution methods: 169: * ------- -- ---------- -------- ------- 170: * 171: * Finding a generalized eigenvector consists basically of solving the 172: * singular triangular system 173: * 174: * (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left) 175: * 176: * Consider finding the i-th right eigenvector (assume all eigenvalues 177: * are real). The equation to be solved is: 178: * n i 179: * 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1 180: * k=j k=j 181: * 182: * where C = (A - w B) (The components v(i+1:n) are 0.) 183: * 184: * The "rowwise" method is: 185: * 186: * (1) v(i) := 1 187: * for j = i-1,. . .,1: 188: * i 189: * (2) compute s = - sum C(j,k) v(k) and 190: * k=j+1 191: * 192: * (3) v(j) := s / C(j,j) 193: * 194: * Step 2 is sometimes called the "dot product" step, since it is an 195: * inner product between the j-th row and the portion of the eigenvector 196: * that has been computed so far. 197: * 198: * The "columnwise" method consists basically in doing the sums 199: * for all the rows in parallel. As each v(j) is computed, the 200: * contribution of v(j) times the j-th column of C is added to the 201: * partial sums. Since FORTRAN arrays are stored columnwise, this has 202: * the advantage that at each step, the elements of C that are accessed 203: * are adjacent to one another, whereas with the rowwise method, the 204: * elements accessed at a step are spaced LDS (and LDP) words apart. 205: * 206: * When finding left eigenvectors, the matrix in question is the 207: * transpose of the one in storage, so the rowwise method then 208: * actually accesses columns of A and B at each step, and so is the 209: * preferred method. 210: * 211: * ===================================================================== 212: * 213: * .. Parameters .. 214: DOUBLE PRECISION ZERO, ONE, SAFETY 215: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, 216: $ SAFETY = 1.0D+2 ) 217: * .. 218: * .. Local Scalars .. 219: LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK, 220: $ ILBBAD, ILCOMP, ILCPLX, LSA, LSB 221: INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE, 222: $ J, JA, JC, JE, JR, JW, NA, NW 223: DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI, 224: $ BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A, 225: $ CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA, 226: $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE, 227: $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX, 228: $ XSCALE 229: * .. 230: * .. Local Arrays .. 231: DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ), 232: $ SUMP( 2, 2 ) 233: * .. 234: * .. External Functions .. 235: LOGICAL LSAME 236: DOUBLE PRECISION DLAMCH 237: EXTERNAL LSAME, DLAMCH 238: * .. 239: * .. External Subroutines .. 240: EXTERNAL DGEMV, DLABAD, DLACPY, DLAG2, DLALN2, XERBLA 241: * .. 242: * .. Intrinsic Functions .. 243: INTRINSIC ABS, MAX, MIN 244: * .. 245: * .. Executable Statements .. 246: * 247: * Decode and Test the input parameters 248: * 249: IF( LSAME( HOWMNY, 'A' ) ) THEN 250: IHWMNY = 1 251: ILALL = .TRUE. 252: ILBACK = .FALSE. 253: ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN 254: IHWMNY = 2 255: ILALL = .FALSE. 256: ILBACK = .FALSE. 257: ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN 258: IHWMNY = 3 259: ILALL = .TRUE. 260: ILBACK = .TRUE. 261: ELSE 262: IHWMNY = -1 263: ILALL = .TRUE. 264: END IF 265: * 266: IF( LSAME( SIDE, 'R' ) ) THEN 267: ISIDE = 1 268: COMPL = .FALSE. 269: COMPR = .TRUE. 270: ELSE IF( LSAME( SIDE, 'L' ) ) THEN 271: ISIDE = 2 272: COMPL = .TRUE. 273: COMPR = .FALSE. 274: ELSE IF( LSAME( SIDE, 'B' ) ) THEN 275: ISIDE = 3 276: COMPL = .TRUE. 277: COMPR = .TRUE. 278: ELSE 279: ISIDE = -1 280: END IF 281: * 282: INFO = 0 283: IF( ISIDE.LT.0 ) THEN 284: INFO = -1 285: ELSE IF( IHWMNY.LT.0 ) THEN 286: INFO = -2 287: ELSE IF( N.LT.0 ) THEN 288: INFO = -4 289: ELSE IF( LDS.LT.MAX( 1, N ) ) THEN 290: INFO = -6 291: ELSE IF( LDP.LT.MAX( 1, N ) ) THEN 292: INFO = -8 293: END IF 294: IF( INFO.NE.0 ) THEN 295: CALL XERBLA( 'DTGEVC', -INFO ) 296: RETURN 297: END IF 298: * 299: * Count the number of eigenvectors to be computed 300: * 301: IF( .NOT.ILALL ) THEN 302: IM = 0 303: ILCPLX = .FALSE. 304: DO 10 J = 1, N 305: IF( ILCPLX ) THEN 306: ILCPLX = .FALSE. 307: GO TO 10 308: END IF 309: IF( J.LT.N ) THEN 310: IF( S( J+1, J ).NE.ZERO ) 311: $ ILCPLX = .TRUE. 312: END IF 313: IF( ILCPLX ) THEN 314: IF( SELECT( J ) .OR. SELECT( J+1 ) ) 315: $ IM = IM + 2 316: ELSE 317: IF( SELECT( J ) ) 318: $ IM = IM + 1 319: END IF 320: 10 CONTINUE 321: ELSE 322: IM = N 323: END IF 324: * 325: * Check 2-by-2 diagonal blocks of A, B 326: * 327: ILABAD = .FALSE. 328: ILBBAD = .FALSE. 329: DO 20 J = 1, N - 1 330: IF( S( J+1, J ).NE.ZERO ) THEN 331: IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR. 332: $ P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE. 333: IF( J.LT.N-1 ) THEN 334: IF( S( J+2, J+1 ).NE.ZERO ) 335: $ ILABAD = .TRUE. 336: END IF 337: END IF 338: 20 CONTINUE 339: * 340: IF( ILABAD ) THEN 341: INFO = -5 342: ELSE IF( ILBBAD ) THEN 343: INFO = -7 344: ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN 345: INFO = -10 346: ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN 347: INFO = -12 348: ELSE IF( MM.LT.IM ) THEN 349: INFO = -13 350: END IF 351: IF( INFO.NE.0 ) THEN 352: CALL XERBLA( 'DTGEVC', -INFO ) 353: RETURN 354: END IF 355: * 356: * Quick return if possible 357: * 358: M = IM 359: IF( N.EQ.0 ) 360: $ RETURN 361: * 362: * Machine Constants 363: * 364: SAFMIN = DLAMCH( 'Safe minimum' ) 365: BIG = ONE / SAFMIN 366: CALL DLABAD( SAFMIN, BIG ) 367: ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 368: SMALL = SAFMIN*N / ULP 369: BIG = ONE / SMALL 370: BIGNUM = ONE / ( SAFMIN*N ) 371: * 372: * Compute the 1-norm of each column of the strictly upper triangular 373: * part (i.e., excluding all elements belonging to the diagonal 374: * blocks) of A and B to check for possible overflow in the 375: * triangular solver. 376: * 377: ANORM = ABS( S( 1, 1 ) ) 378: IF( N.GT.1 ) 379: $ ANORM = ANORM + ABS( S( 2, 1 ) ) 380: BNORM = ABS( P( 1, 1 ) ) 381: WORK( 1 ) = ZERO 382: WORK( N+1 ) = ZERO 383: * 384: DO 50 J = 2, N 385: TEMP = ZERO 386: TEMP2 = ZERO 387: IF( S( J, J-1 ).EQ.ZERO ) THEN 388: IEND = J - 1 389: ELSE 390: IEND = J - 2 391: END IF 392: DO 30 I = 1, IEND 393: TEMP = TEMP + ABS( S( I, J ) ) 394: TEMP2 = TEMP2 + ABS( P( I, J ) ) 395: 30 CONTINUE 396: WORK( J ) = TEMP 397: WORK( N+J ) = TEMP2 398: DO 40 I = IEND + 1, MIN( J+1, N ) 399: TEMP = TEMP + ABS( S( I, J ) ) 400: TEMP2 = TEMP2 + ABS( P( I, J ) ) 401: 40 CONTINUE 402: ANORM = MAX( ANORM, TEMP ) 403: BNORM = MAX( BNORM, TEMP2 ) 404: 50 CONTINUE 405: * 406: ASCALE = ONE / MAX( ANORM, SAFMIN ) 407: BSCALE = ONE / MAX( BNORM, SAFMIN ) 408: * 409: * Left eigenvectors 410: * 411: IF( COMPL ) THEN 412: IEIG = 0 413: * 414: * Main loop over eigenvalues 415: * 416: ILCPLX = .FALSE. 417: DO 220 JE = 1, N 418: * 419: * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 420: * (b) this would be the second of a complex pair. 421: * Check for complex eigenvalue, so as to be sure of which 422: * entry(-ies) of SELECT to look at. 423: * 424: IF( ILCPLX ) THEN 425: ILCPLX = .FALSE. 426: GO TO 220 427: END IF 428: NW = 1 429: IF( JE.LT.N ) THEN 430: IF( S( JE+1, JE ).NE.ZERO ) THEN 431: ILCPLX = .TRUE. 432: NW = 2 433: END IF 434: END IF 435: IF( ILALL ) THEN 436: ILCOMP = .TRUE. 437: ELSE IF( ILCPLX ) THEN 438: ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 ) 439: ELSE 440: ILCOMP = SELECT( JE ) 441: END IF 442: IF( .NOT.ILCOMP ) 443: $ GO TO 220 444: * 445: * Decide if (a) singular pencil, (b) real eigenvalue, or 446: * (c) complex eigenvalue. 447: * 448: IF( .NOT.ILCPLX ) THEN 449: IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND. 450: $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN 451: * 452: * Singular matrix pencil -- return unit eigenvector 453: * 454: IEIG = IEIG + 1 455: DO 60 JR = 1, N 456: VL( JR, IEIG ) = ZERO 457: 60 CONTINUE 458: VL( IEIG, IEIG ) = ONE 459: GO TO 220 460: END IF 461: END IF 462: * 463: * Clear vector 464: * 465: DO 70 JR = 1, NW*N 466: WORK( 2*N+JR ) = ZERO 467: 70 CONTINUE 468: * T 469: * Compute coefficients in ( a A - b B ) y = 0 470: * a is ACOEF 471: * b is BCOEFR + i*BCOEFI 472: * 473: IF( .NOT.ILCPLX ) THEN 474: * 475: * Real eigenvalue 476: * 477: TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE, 478: $ ABS( P( JE, JE ) )*BSCALE, SAFMIN ) 479: SALFAR = ( TEMP*S( JE, JE ) )*ASCALE 480: SBETA = ( TEMP*P( JE, JE ) )*BSCALE 481: ACOEF = SBETA*ASCALE 482: BCOEFR = SALFAR*BSCALE 483: BCOEFI = ZERO 484: * 485: * Scale to avoid underflow 486: * 487: SCALE = ONE 488: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL 489: LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. 490: $ SMALL 491: IF( LSA ) 492: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 493: IF( LSB ) 494: $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* 495: $ MIN( BNORM, BIG ) ) 496: IF( LSA .OR. LSB ) THEN 497: SCALE = MIN( SCALE, ONE / 498: $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), 499: $ ABS( BCOEFR ) ) ) ) 500: IF( LSA ) THEN 501: ACOEF = ASCALE*( SCALE*SBETA ) 502: ELSE 503: ACOEF = SCALE*ACOEF 504: END IF 505: IF( LSB ) THEN 506: BCOEFR = BSCALE*( SCALE*SALFAR ) 507: ELSE 508: BCOEFR = SCALE*BCOEFR 509: END IF 510: END IF 511: ACOEFA = ABS( ACOEF ) 512: BCOEFA = ABS( BCOEFR ) 513: * 514: * First component is 1 515: * 516: WORK( 2*N+JE ) = ONE 517: XMAX = ONE 518: ELSE 519: * 520: * Complex eigenvalue 521: * 522: CALL DLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP, 523: $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 524: $ BCOEFI ) 525: BCOEFI = -BCOEFI 526: IF( BCOEFI.EQ.ZERO ) THEN 527: INFO = JE 528: RETURN 529: END IF 530: * 531: * Scale to avoid over/underflow 532: * 533: ACOEFA = ABS( ACOEF ) 534: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 535: SCALE = ONE 536: IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) 537: $ SCALE = ( SAFMIN / ULP ) / ACOEFA 538: IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) 539: $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) 540: IF( SAFMIN*ACOEFA.GT.ASCALE ) 541: $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) 542: IF( SAFMIN*BCOEFA.GT.BSCALE ) 543: $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) 544: IF( SCALE.NE.ONE ) THEN 545: ACOEF = SCALE*ACOEF 546: ACOEFA = ABS( ACOEF ) 547: BCOEFR = SCALE*BCOEFR 548: BCOEFI = SCALE*BCOEFI 549: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 550: END IF 551: * 552: * Compute first two components of eigenvector 553: * 554: TEMP = ACOEF*S( JE+1, JE ) 555: TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE ) 556: TEMP2I = -BCOEFI*P( JE, JE ) 557: IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 558: WORK( 2*N+JE ) = ONE 559: WORK( 3*N+JE ) = ZERO 560: WORK( 2*N+JE+1 ) = -TEMP2R / TEMP 561: WORK( 3*N+JE+1 ) = -TEMP2I / TEMP 562: ELSE 563: WORK( 2*N+JE+1 ) = ONE 564: WORK( 3*N+JE+1 ) = ZERO 565: TEMP = ACOEF*S( JE, JE+1 ) 566: WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF* 567: $ S( JE+1, JE+1 ) ) / TEMP 568: WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP 569: END IF 570: XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 571: $ ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) ) 572: END IF 573: * 574: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 575: * 576: * T 577: * Triangular solve of (a A - b B) y = 0 578: * 579: * T 580: * (rowwise in (a A - b B) , or columnwise in (a A - b B) ) 581: * 582: IL2BY2 = .FALSE. 583: * 584: DO 160 J = JE + NW, N 585: IF( IL2BY2 ) THEN 586: IL2BY2 = .FALSE. 587: GO TO 160 588: END IF 589: * 590: NA = 1 591: BDIAG( 1 ) = P( J, J ) 592: IF( J.LT.N ) THEN 593: IF( S( J+1, J ).NE.ZERO ) THEN 594: IL2BY2 = .TRUE. 595: BDIAG( 2 ) = P( J+1, J+1 ) 596: NA = 2 597: END IF 598: END IF 599: * 600: * Check whether scaling is necessary for dot products 601: * 602: XSCALE = ONE / MAX( ONE, XMAX ) 603: TEMP = MAX( WORK( J ), WORK( N+J ), 604: $ ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) ) 605: IF( IL2BY2 ) 606: $ TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ), 607: $ ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) ) 608: IF( TEMP.GT.BIGNUM*XSCALE ) THEN 609: DO 90 JW = 0, NW - 1 610: DO 80 JR = JE, J - 1 611: WORK( ( JW+2 )*N+JR ) = XSCALE* 612: $ WORK( ( JW+2 )*N+JR ) 613: 80 CONTINUE 614: 90 CONTINUE 615: XMAX = XMAX*XSCALE 616: END IF 617: * 618: * Compute dot products 619: * 620: * j-1 621: * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k) 622: * k=je 623: * 624: * To reduce the op count, this is done as 625: * 626: * _ j-1 _ j-1 627: * a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) ) 628: * k=je k=je 629: * 630: * which may cause underflow problems if A or B are close 631: * to underflow. (E.g., less than SMALL.) 632: * 633: * 634: * A series of compiler directives to defeat vectorization 635: * for the next loop 636: * 637: *$PL$ CMCHAR=' ' 638: CDIR$ NEXTSCALAR 639: C$DIR SCALAR 640: CDIR$ NEXT SCALAR 641: CVD$L NOVECTOR 642: CDEC$ NOVECTOR 643: CVD$ NOVECTOR 644: *VDIR NOVECTOR 645: *VOCL LOOP,SCALAR 646: CIBM PREFER SCALAR 647: *$PL$ CMCHAR='*' 648: * 649: DO 120 JW = 1, NW 650: * 651: *$PL$ CMCHAR=' ' 652: CDIR$ NEXTSCALAR 653: C$DIR SCALAR 654: CDIR$ NEXT SCALAR 655: CVD$L NOVECTOR 656: CDEC$ NOVECTOR 657: CVD$ NOVECTOR 658: *VDIR NOVECTOR 659: *VOCL LOOP,SCALAR 660: CIBM PREFER SCALAR 661: *$PL$ CMCHAR='*' 662: * 663: DO 110 JA = 1, NA 664: SUMS( JA, JW ) = ZERO 665: SUMP( JA, JW ) = ZERO 666: * 667: DO 100 JR = JE, J - 1 668: SUMS( JA, JW ) = SUMS( JA, JW ) + 669: $ S( JR, J+JA-1 )* 670: $ WORK( ( JW+1 )*N+JR ) 671: SUMP( JA, JW ) = SUMP( JA, JW ) + 672: $ P( JR, J+JA-1 )* 673: $ WORK( ( JW+1 )*N+JR ) 674: 100 CONTINUE 675: 110 CONTINUE 676: 120 CONTINUE 677: * 678: *$PL$ CMCHAR=' ' 679: CDIR$ NEXTSCALAR 680: C$DIR SCALAR 681: CDIR$ NEXT SCALAR 682: CVD$L NOVECTOR 683: CDEC$ NOVECTOR 684: CVD$ NOVECTOR 685: *VDIR NOVECTOR 686: *VOCL LOOP,SCALAR 687: CIBM PREFER SCALAR 688: *$PL$ CMCHAR='*' 689: * 690: DO 130 JA = 1, NA 691: IF( ILCPLX ) THEN 692: SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 693: $ BCOEFR*SUMP( JA, 1 ) - 694: $ BCOEFI*SUMP( JA, 2 ) 695: SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) + 696: $ BCOEFR*SUMP( JA, 2 ) + 697: $ BCOEFI*SUMP( JA, 1 ) 698: ELSE 699: SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 700: $ BCOEFR*SUMP( JA, 1 ) 701: END IF 702: 130 CONTINUE 703: * 704: * T 705: * Solve ( a A - b B ) y = SUM(,) 706: * with scaling and perturbation of the denominator 707: * 708: CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS, 709: $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR, 710: $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP, 711: $ IINFO ) 712: IF( SCALE.LT.ONE ) THEN 713: DO 150 JW = 0, NW - 1 714: DO 140 JR = JE, J - 1 715: WORK( ( JW+2 )*N+JR ) = SCALE* 716: $ WORK( ( JW+2 )*N+JR ) 717: 140 CONTINUE 718: 150 CONTINUE 719: XMAX = SCALE*XMAX 720: END IF 721: XMAX = MAX( XMAX, TEMP ) 722: 160 CONTINUE 723: * 724: * Copy eigenvector to VL, back transforming if 725: * HOWMNY='B'. 726: * 727: IEIG = IEIG + 1 728: IF( ILBACK ) THEN 729: DO 170 JW = 0, NW - 1 730: CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL, 731: $ WORK( ( JW+2 )*N+JE ), 1, ZERO, 732: $ WORK( ( JW+4 )*N+1 ), 1 ) 733: 170 CONTINUE 734: CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ), 735: $ LDVL ) 736: IBEG = 1 737: ELSE 738: CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ), 739: $ LDVL ) 740: IBEG = JE 741: END IF 742: * 743: * Scale eigenvector 744: * 745: XMAX = ZERO 746: IF( ILCPLX ) THEN 747: DO 180 J = IBEG, N 748: XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+ 749: $ ABS( VL( J, IEIG+1 ) ) ) 750: 180 CONTINUE 751: ELSE 752: DO 190 J = IBEG, N 753: XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) ) 754: 190 CONTINUE 755: END IF 756: * 757: IF( XMAX.GT.SAFMIN ) THEN 758: XSCALE = ONE / XMAX 759: * 760: DO 210 JW = 0, NW - 1 761: DO 200 JR = IBEG, N 762: VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW ) 763: 200 CONTINUE 764: 210 CONTINUE 765: END IF 766: IEIG = IEIG + NW - 1 767: * 768: 220 CONTINUE 769: END IF 770: * 771: * Right eigenvectors 772: * 773: IF( COMPR ) THEN 774: IEIG = IM + 1 775: * 776: * Main loop over eigenvalues 777: * 778: ILCPLX = .FALSE. 779: DO 500 JE = N, 1, -1 780: * 781: * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 782: * (b) this would be the second of a complex pair. 783: * Check for complex eigenvalue, so as to be sure of which 784: * entry(-ies) of SELECT to look at -- if complex, SELECT(JE) 785: * or SELECT(JE-1). 786: * If this is a complex pair, the 2-by-2 diagonal block 787: * corresponding to the eigenvalue is in rows/columns JE-1:JE 788: * 789: IF( ILCPLX ) THEN 790: ILCPLX = .FALSE. 791: GO TO 500 792: END IF 793: NW = 1 794: IF( JE.GT.1 ) THEN 795: IF( S( JE, JE-1 ).NE.ZERO ) THEN 796: ILCPLX = .TRUE. 797: NW = 2 798: END IF 799: END IF 800: IF( ILALL ) THEN 801: ILCOMP = .TRUE. 802: ELSE IF( ILCPLX ) THEN 803: ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 ) 804: ELSE 805: ILCOMP = SELECT( JE ) 806: END IF 807: IF( .NOT.ILCOMP ) 808: $ GO TO 500 809: * 810: * Decide if (a) singular pencil, (b) real eigenvalue, or 811: * (c) complex eigenvalue. 812: * 813: IF( .NOT.ILCPLX ) THEN 814: IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND. 815: $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN 816: * 817: * Singular matrix pencil -- unit eigenvector 818: * 819: IEIG = IEIG - 1 820: DO 230 JR = 1, N 821: VR( JR, IEIG ) = ZERO 822: 230 CONTINUE 823: VR( IEIG, IEIG ) = ONE 824: GO TO 500 825: END IF 826: END IF 827: * 828: * Clear vector 829: * 830: DO 250 JW = 0, NW - 1 831: DO 240 JR = 1, N 832: WORK( ( JW+2 )*N+JR ) = ZERO 833: 240 CONTINUE 834: 250 CONTINUE 835: * 836: * Compute coefficients in ( a A - b B ) x = 0 837: * a is ACOEF 838: * b is BCOEFR + i*BCOEFI 839: * 840: IF( .NOT.ILCPLX ) THEN 841: * 842: * Real eigenvalue 843: * 844: TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE, 845: $ ABS( P( JE, JE ) )*BSCALE, SAFMIN ) 846: SALFAR = ( TEMP*S( JE, JE ) )*ASCALE 847: SBETA = ( TEMP*P( JE, JE ) )*BSCALE 848: ACOEF = SBETA*ASCALE 849: BCOEFR = SALFAR*BSCALE 850: BCOEFI = ZERO 851: * 852: * Scale to avoid underflow 853: * 854: SCALE = ONE 855: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL 856: LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. 857: $ SMALL 858: IF( LSA ) 859: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 860: IF( LSB ) 861: $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* 862: $ MIN( BNORM, BIG ) ) 863: IF( LSA .OR. LSB ) THEN 864: SCALE = MIN( SCALE, ONE / 865: $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), 866: $ ABS( BCOEFR ) ) ) ) 867: IF( LSA ) THEN 868: ACOEF = ASCALE*( SCALE*SBETA ) 869: ELSE 870: ACOEF = SCALE*ACOEF 871: END IF 872: IF( LSB ) THEN 873: BCOEFR = BSCALE*( SCALE*SALFAR ) 874: ELSE 875: BCOEFR = SCALE*BCOEFR 876: END IF 877: END IF 878: ACOEFA = ABS( ACOEF ) 879: BCOEFA = ABS( BCOEFR ) 880: * 881: * First component is 1 882: * 883: WORK( 2*N+JE ) = ONE 884: XMAX = ONE 885: * 886: * Compute contribution from column JE of A and B to sum 887: * (See "Further Details", above.) 888: * 889: DO 260 JR = 1, JE - 1 890: WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) - 891: $ ACOEF*S( JR, JE ) 892: 260 CONTINUE 893: ELSE 894: * 895: * Complex eigenvalue 896: * 897: CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP, 898: $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 899: $ BCOEFI ) 900: IF( BCOEFI.EQ.ZERO ) THEN 901: INFO = JE - 1 902: RETURN 903: END IF 904: * 905: * Scale to avoid over/underflow 906: * 907: ACOEFA = ABS( ACOEF ) 908: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 909: SCALE = ONE 910: IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) 911: $ SCALE = ( SAFMIN / ULP ) / ACOEFA 912: IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) 913: $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) 914: IF( SAFMIN*ACOEFA.GT.ASCALE ) 915: $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) 916: IF( SAFMIN*BCOEFA.GT.BSCALE ) 917: $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) 918: IF( SCALE.NE.ONE ) THEN 919: ACOEF = SCALE*ACOEF 920: ACOEFA = ABS( ACOEF ) 921: BCOEFR = SCALE*BCOEFR 922: BCOEFI = SCALE*BCOEFI 923: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 924: END IF 925: * 926: * Compute first two components of eigenvector 927: * and contribution to sums 928: * 929: TEMP = ACOEF*S( JE, JE-1 ) 930: TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE ) 931: TEMP2I = -BCOEFI*P( JE, JE ) 932: IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 933: WORK( 2*N+JE ) = ONE 934: WORK( 3*N+JE ) = ZERO 935: WORK( 2*N+JE-1 ) = -TEMP2R / TEMP 936: WORK( 3*N+JE-1 ) = -TEMP2I / TEMP 937: ELSE 938: WORK( 2*N+JE-1 ) = ONE 939: WORK( 3*N+JE-1 ) = ZERO 940: TEMP = ACOEF*S( JE-1, JE ) 941: WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF* 942: $ S( JE-1, JE-1 ) ) / TEMP 943: WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP 944: END IF 945: * 946: XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 947: $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) ) 948: * 949: * Compute contribution from columns JE and JE-1 950: * of A and B to the sums. 951: * 952: CREALA = ACOEF*WORK( 2*N+JE-1 ) 953: CIMAGA = ACOEF*WORK( 3*N+JE-1 ) 954: CREALB = BCOEFR*WORK( 2*N+JE-1 ) - 955: $ BCOEFI*WORK( 3*N+JE-1 ) 956: CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) + 957: $ BCOEFR*WORK( 3*N+JE-1 ) 958: CRE2A = ACOEF*WORK( 2*N+JE ) 959: CIM2A = ACOEF*WORK( 3*N+JE ) 960: CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE ) 961: CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE ) 962: DO 270 JR = 1, JE - 2 963: WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) + 964: $ CREALB*P( JR, JE-1 ) - 965: $ CRE2A*S( JR, JE ) + CRE2B*P( JR, JE ) 966: WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) + 967: $ CIMAGB*P( JR, JE-1 ) - 968: $ CIM2A*S( JR, JE ) + CIM2B*P( JR, JE ) 969: 270 CONTINUE 970: END IF 971: * 972: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 973: * 974: * Columnwise triangular solve of (a A - b B) x = 0 975: * 976: IL2BY2 = .FALSE. 977: DO 370 J = JE - NW, 1, -1 978: * 979: * If a 2-by-2 block, is in position j-1:j, wait until 980: * next iteration to process it (when it will be j:j+1) 981: * 982: IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN 983: IF( S( J, J-1 ).NE.ZERO ) THEN 984: IL2BY2 = .TRUE. 985: GO TO 370 986: END IF 987: END IF 988: BDIAG( 1 ) = P( J, J ) 989: IF( IL2BY2 ) THEN 990: NA = 2 991: BDIAG( 2 ) = P( J+1, J+1 ) 992: ELSE 993: NA = 1 994: END IF 995: * 996: * Compute x(j) (and x(j+1), if 2-by-2 block) 997: * 998: CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ), 999: $ LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ), 1000: $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP, 1001: $ IINFO ) 1002: IF( SCALE.LT.ONE ) THEN 1003: * 1004: DO 290 JW = 0, NW - 1 1005: DO 280 JR = 1, JE 1006: WORK( ( JW+2 )*N+JR ) = SCALE* 1007: $ WORK( ( JW+2 )*N+JR ) 1008: 280 CONTINUE 1009: 290 CONTINUE 1010: END IF 1011: XMAX = MAX( SCALE*XMAX, TEMP ) 1012: * 1013: DO 310 JW = 1, NW 1014: DO 300 JA = 1, NA 1015: WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW ) 1016: 300 CONTINUE 1017: 310 CONTINUE 1018: * 1019: * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling 1020: * 1021: IF( J.GT.1 ) THEN 1022: * 1023: * Check whether scaling is necessary for sum. 1024: * 1025: XSCALE = ONE / MAX( ONE, XMAX ) 1026: TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J ) 1027: IF( IL2BY2 ) 1028: $ TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA* 1029: $ WORK( N+J+1 ) ) 1030: TEMP = MAX( TEMP, ACOEFA, BCOEFA ) 1031: IF( TEMP.GT.BIGNUM*XSCALE ) THEN 1032: * 1033: DO 330 JW = 0, NW - 1 1034: DO 320 JR = 1, JE 1035: WORK( ( JW+2 )*N+JR ) = XSCALE* 1036: $ WORK( ( JW+2 )*N+JR ) 1037: 320 CONTINUE 1038: 330 CONTINUE 1039: XMAX = XMAX*XSCALE 1040: END IF 1041: * 1042: * Compute the contributions of the off-diagonals of 1043: * column j (and j+1, if 2-by-2 block) of A and B to the 1044: * sums. 1045: * 1046: * 1047: DO 360 JA = 1, NA 1048: IF( ILCPLX ) THEN 1049: CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 1050: CIMAGA = ACOEF*WORK( 3*N+J+JA-1 ) 1051: CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) - 1052: $ BCOEFI*WORK( 3*N+J+JA-1 ) 1053: CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) + 1054: $ BCOEFR*WORK( 3*N+J+JA-1 ) 1055: DO 340 JR = 1, J - 1 1056: WORK( 2*N+JR ) = WORK( 2*N+JR ) - 1057: $ CREALA*S( JR, J+JA-1 ) + 1058: $ CREALB*P( JR, J+JA-1 ) 1059: WORK( 3*N+JR ) = WORK( 3*N+JR ) - 1060: $ CIMAGA*S( JR, J+JA-1 ) + 1061: $ CIMAGB*P( JR, J+JA-1 ) 1062: 340 CONTINUE 1063: ELSE 1064: CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 1065: CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) 1066: DO 350 JR = 1, J - 1 1067: WORK( 2*N+JR ) = WORK( 2*N+JR ) - 1068: $ CREALA*S( JR, J+JA-1 ) + 1069: $ CREALB*P( JR, J+JA-1 ) 1070: 350 CONTINUE 1071: END IF 1072: 360 CONTINUE 1073: END IF 1074: * 1075: IL2BY2 = .FALSE. 1076: 370 CONTINUE 1077: * 1078: * Copy eigenvector to VR, back transforming if 1079: * HOWMNY='B'. 1080: * 1081: IEIG = IEIG - NW 1082: IF( ILBACK ) THEN 1083: * 1084: DO 410 JW = 0, NW - 1 1085: DO 380 JR = 1, N 1086: WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )* 1087: $ VR( JR, 1 ) 1088: 380 CONTINUE 1089: * 1090: * A series of compiler directives to defeat 1091: * vectorization for the next loop 1092: * 1093: * 1094: DO 400 JC = 2, JE 1095: DO 390 JR = 1, N 1096: WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) + 1097: $ WORK( ( JW+2 )*N+JC )*VR( JR, JC ) 1098: 390 CONTINUE 1099: 400 CONTINUE 1100: 410 CONTINUE 1101: * 1102: DO 430 JW = 0, NW - 1 1103: DO 420 JR = 1, N 1104: VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR ) 1105: 420 CONTINUE 1106: 430 CONTINUE 1107: * 1108: IEND = N 1109: ELSE 1110: DO 450 JW = 0, NW - 1 1111: DO 440 JR = 1, N 1112: VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR ) 1113: 440 CONTINUE 1114: 450 CONTINUE 1115: * 1116: IEND = JE 1117: END IF 1118: * 1119: * Scale eigenvector 1120: * 1121: XMAX = ZERO 1122: IF( ILCPLX ) THEN 1123: DO 460 J = 1, IEND 1124: XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+ 1125: $ ABS( VR( J, IEIG+1 ) ) ) 1126: 460 CONTINUE 1127: ELSE 1128: DO 470 J = 1, IEND 1129: XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) ) 1130: 470 CONTINUE 1131: END IF 1132: * 1133: IF( XMAX.GT.SAFMIN ) THEN 1134: XSCALE = ONE / XMAX 1135: DO 490 JW = 0, NW - 1 1136: DO 480 JR = 1, IEND 1137: VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW ) 1138: 480 CONTINUE 1139: 490 CONTINUE 1140: END IF 1141: 500 CONTINUE 1142: END IF 1143: * 1144: RETURN 1145: * 1146: * End of DTGEVC 1147: * 1148: END