![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, 2: $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, 3: $ LWORK, INFO ) 4: * 5: * -- LAPACK routine (version 3.2.1) -- 6: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 8: * -- April 2009 -- 9: * 10: * .. Scalar Arguments .. 11: CHARACTER COMPQ, COMPZ, JOB 12: INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N 13: * .. 14: * .. Array Arguments .. 15: DOUBLE PRECISION ALPHAI( * ), ALPHAR( * ), BETA( * ), 16: $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ), 17: $ WORK( * ), Z( LDZ, * ) 18: * .. 19: * 20: * Purpose 21: * ======= 22: * 23: * DHGEQZ computes the eigenvalues of a real matrix pair (H,T), 24: * where H is an upper Hessenberg matrix and T is upper triangular, 25: * using the double-shift QZ method. 26: * Matrix pairs of this type are produced by the reduction to 27: * generalized upper Hessenberg form of a real matrix pair (A,B): 28: * 29: * A = Q1*H*Z1**T, B = Q1*T*Z1**T, 30: * 31: * as computed by DGGHRD. 32: * 33: * If JOB='S', then the Hessenberg-triangular pair (H,T) is 34: * also reduced to generalized Schur form, 35: * 36: * H = Q*S*Z**T, T = Q*P*Z**T, 37: * 38: * where Q and Z are orthogonal matrices, P is an upper triangular 39: * matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2 40: * diagonal blocks. 41: * 42: * The 1-by-1 blocks correspond to real eigenvalues of the matrix pair 43: * (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of 44: * eigenvalues. 45: * 46: * Additionally, the 2-by-2 upper triangular diagonal blocks of P 47: * corresponding to 2-by-2 blocks of S are reduced to positive diagonal 48: * form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0, 49: * P(j,j) > 0, and P(j+1,j+1) > 0. 50: * 51: * Optionally, the orthogonal matrix Q from the generalized Schur 52: * factorization may be postmultiplied into an input matrix Q1, and the 53: * orthogonal matrix Z may be postmultiplied into an input matrix Z1. 54: * If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced 55: * the matrix pair (A,B) to generalized upper Hessenberg form, then the 56: * output matrices Q1*Q and Z1*Z are the orthogonal factors from the 57: * generalized Schur factorization of (A,B): 58: * 59: * A = (Q1*Q)*S*(Z1*Z)**T, B = (Q1*Q)*P*(Z1*Z)**T. 60: * 61: * To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently, 62: * of (A,B)) are computed as a pair of values (alpha,beta), where alpha is 63: * complex and beta real. 64: * If beta is nonzero, lambda = alpha / beta is an eigenvalue of the 65: * generalized nonsymmetric eigenvalue problem (GNEP) 66: * A*x = lambda*B*x 67: * and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the 68: * alternate form of the GNEP 69: * mu*A*y = B*y. 70: * Real eigenvalues can be read directly from the generalized Schur 71: * form: 72: * alpha = S(i,i), beta = P(i,i). 73: * 74: * Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix 75: * Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), 76: * pp. 241--256. 77: * 78: * Arguments 79: * ========= 80: * 81: * JOB (input) CHARACTER*1 82: * = 'E': Compute eigenvalues only; 83: * = 'S': Compute eigenvalues and the Schur form. 84: * 85: * COMPQ (input) CHARACTER*1 86: * = 'N': Left Schur vectors (Q) are not computed; 87: * = 'I': Q is initialized to the unit matrix and the matrix Q 88: * of left Schur vectors of (H,T) is returned; 89: * = 'V': Q must contain an orthogonal matrix Q1 on entry and 90: * the product Q1*Q is returned. 91: * 92: * COMPZ (input) CHARACTER*1 93: * = 'N': Right Schur vectors (Z) are not computed; 94: * = 'I': Z is initialized to the unit matrix and the matrix Z 95: * of right Schur vectors of (H,T) is returned; 96: * = 'V': Z must contain an orthogonal matrix Z1 on entry and 97: * the product Z1*Z is returned. 98: * 99: * N (input) INTEGER 100: * The order of the matrices H, T, Q, and Z. N >= 0. 101: * 102: * ILO (input) INTEGER 103: * IHI (input) INTEGER 104: * ILO and IHI mark the rows and columns of H which are in 105: * Hessenberg form. It is assumed that A is already upper 106: * triangular in rows and columns 1:ILO-1 and IHI+1:N. 107: * If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. 108: * 109: * H (input/output) DOUBLE PRECISION array, dimension (LDH, N) 110: * On entry, the N-by-N upper Hessenberg matrix H. 111: * On exit, if JOB = 'S', H contains the upper quasi-triangular 112: * matrix S from the generalized Schur factorization; 113: * 2-by-2 diagonal blocks (corresponding to complex conjugate 114: * pairs of eigenvalues) are returned in standard form, with 115: * H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. 116: * If JOB = 'E', the diagonal blocks of H match those of S, but 117: * the rest of H is unspecified. 118: * 119: * LDH (input) INTEGER 120: * The leading dimension of the array H. LDH >= max( 1, N ). 121: * 122: * T (input/output) DOUBLE PRECISION array, dimension (LDT, N) 123: * On entry, the N-by-N upper triangular matrix T. 124: * On exit, if JOB = 'S', T contains the upper triangular 125: * matrix P from the generalized Schur factorization; 126: * 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S 127: * are reduced to positive diagonal form, i.e., if H(j+1,j) is 128: * non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and 129: * T(j+1,j+1) > 0. 130: * If JOB = 'E', the diagonal blocks of T match those of P, but 131: * the rest of T is unspecified. 132: * 133: * LDT (input) INTEGER 134: * The leading dimension of the array T. LDT >= max( 1, N ). 135: * 136: * ALPHAR (output) DOUBLE PRECISION array, dimension (N) 137: * The real parts of each scalar alpha defining an eigenvalue 138: * of GNEP. 139: * 140: * ALPHAI (output) DOUBLE PRECISION array, dimension (N) 141: * The imaginary parts of each scalar alpha defining an 142: * eigenvalue of GNEP. 143: * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if 144: * positive, then the j-th and (j+1)-st eigenvalues are a 145: * complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j). 146: * 147: * BETA (output) DOUBLE PRECISION array, dimension (N) 148: * The scalars beta that define the eigenvalues of GNEP. 149: * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and 150: * beta = BETA(j) represent the j-th eigenvalue of the matrix 151: * pair (A,B), in one of the forms lambda = alpha/beta or 152: * mu = beta/alpha. Since either lambda or mu may overflow, 153: * they should not, in general, be computed. 154: * 155: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) 156: * On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in 157: * the reduction of (A,B) to generalized Hessenberg form. 158: * On exit, if COMPZ = 'I', the orthogonal matrix of left Schur 159: * vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix 160: * of left Schur vectors of (A,B). 161: * Not referenced if COMPZ = 'N'. 162: * 163: * LDQ (input) INTEGER 164: * The leading dimension of the array Q. LDQ >= 1. 165: * If COMPQ='V' or 'I', then LDQ >= N. 166: * 167: * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) 168: * On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in 169: * the reduction of (A,B) to generalized Hessenberg form. 170: * On exit, if COMPZ = 'I', the orthogonal matrix of 171: * right Schur vectors of (H,T), and if COMPZ = 'V', the 172: * orthogonal matrix of right Schur vectors of (A,B). 173: * Not referenced if COMPZ = 'N'. 174: * 175: * LDZ (input) INTEGER 176: * The leading dimension of the array Z. LDZ >= 1. 177: * If COMPZ='V' or 'I', then LDZ >= N. 178: * 179: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 180: * On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. 181: * 182: * LWORK (input) INTEGER 183: * The dimension of the array WORK. LWORK >= max(1,N). 184: * 185: * If LWORK = -1, then a workspace query is assumed; the routine 186: * only calculates the optimal size of the WORK array, returns 187: * this value as the first entry of the WORK array, and no error 188: * message related to LWORK is issued by XERBLA. 189: * 190: * INFO (output) INTEGER 191: * = 0: successful exit 192: * < 0: if INFO = -i, the i-th argument had an illegal value 193: * = 1,...,N: the QZ iteration did not converge. (H,T) is not 194: * in Schur form, but ALPHAR(i), ALPHAI(i), and 195: * BETA(i), i=INFO+1,...,N should be correct. 196: * = N+1,...,2*N: the shift calculation failed. (H,T) is not 197: * in Schur form, but ALPHAR(i), ALPHAI(i), and 198: * BETA(i), i=INFO-N+1,...,N should be correct. 199: * 200: * Further Details 201: * =============== 202: * 203: * Iteration counters: 204: * 205: * JITER -- counts iterations. 206: * IITER -- counts iterations run since ILAST was last 207: * changed. This is therefore reset only when a 1-by-1 or 208: * 2-by-2 block deflates off the bottom. 209: * 210: * ===================================================================== 211: * 212: * .. Parameters .. 213: * $ SAFETY = 1.0E+0 ) 214: DOUBLE PRECISION HALF, ZERO, ONE, SAFETY 215: PARAMETER ( HALF = 0.5D+0, ZERO = 0.0D+0, ONE = 1.0D+0, 216: $ SAFETY = 1.0D+2 ) 217: * .. 218: * .. Local Scalars .. 219: LOGICAL ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ, 220: $ LQUERY 221: INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST, 222: $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER, 223: $ JR, MAXIT 224: DOUBLE PRECISION A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11, 225: $ AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L, 226: $ AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I, 227: $ B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE, 228: $ BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL, 229: $ CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX, 230: $ SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1, 231: $ TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L, 232: $ U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR, 233: $ WR2 234: * .. 235: * .. Local Arrays .. 236: DOUBLE PRECISION V( 3 ) 237: * .. 238: * .. External Functions .. 239: LOGICAL LSAME 240: DOUBLE PRECISION DLAMCH, DLANHS, DLAPY2, DLAPY3 241: EXTERNAL LSAME, DLAMCH, DLANHS, DLAPY2, DLAPY3 242: * .. 243: * .. External Subroutines .. 244: EXTERNAL DLAG2, DLARFG, DLARTG, DLASET, DLASV2, DROT, 245: $ XERBLA 246: * .. 247: * .. Intrinsic Functions .. 248: INTRINSIC ABS, DBLE, MAX, MIN, SQRT 249: * .. 250: * .. Executable Statements .. 251: * 252: * Decode JOB, COMPQ, COMPZ 253: * 254: IF( LSAME( JOB, 'E' ) ) THEN 255: ILSCHR = .FALSE. 256: ISCHUR = 1 257: ELSE IF( LSAME( JOB, 'S' ) ) THEN 258: ILSCHR = .TRUE. 259: ISCHUR = 2 260: ELSE 261: ISCHUR = 0 262: END IF 263: * 264: IF( LSAME( COMPQ, 'N' ) ) THEN 265: ILQ = .FALSE. 266: ICOMPQ = 1 267: ELSE IF( LSAME( COMPQ, 'V' ) ) THEN 268: ILQ = .TRUE. 269: ICOMPQ = 2 270: ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 271: ILQ = .TRUE. 272: ICOMPQ = 3 273: ELSE 274: ICOMPQ = 0 275: END IF 276: * 277: IF( LSAME( COMPZ, 'N' ) ) THEN 278: ILZ = .FALSE. 279: ICOMPZ = 1 280: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 281: ILZ = .TRUE. 282: ICOMPZ = 2 283: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 284: ILZ = .TRUE. 285: ICOMPZ = 3 286: ELSE 287: ICOMPZ = 0 288: END IF 289: * 290: * Check Argument Values 291: * 292: INFO = 0 293: WORK( 1 ) = MAX( 1, N ) 294: LQUERY = ( LWORK.EQ.-1 ) 295: IF( ISCHUR.EQ.0 ) THEN 296: INFO = -1 297: ELSE IF( ICOMPQ.EQ.0 ) THEN 298: INFO = -2 299: ELSE IF( ICOMPZ.EQ.0 ) THEN 300: INFO = -3 301: ELSE IF( N.LT.0 ) THEN 302: INFO = -4 303: ELSE IF( ILO.LT.1 ) THEN 304: INFO = -5 305: ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 306: INFO = -6 307: ELSE IF( LDH.LT.N ) THEN 308: INFO = -8 309: ELSE IF( LDT.LT.N ) THEN 310: INFO = -10 311: ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN 312: INFO = -15 313: ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN 314: INFO = -17 315: ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 316: INFO = -19 317: END IF 318: IF( INFO.NE.0 ) THEN 319: CALL XERBLA( 'DHGEQZ', -INFO ) 320: RETURN 321: ELSE IF( LQUERY ) THEN 322: RETURN 323: END IF 324: * 325: * Quick return if possible 326: * 327: IF( N.LE.0 ) THEN 328: WORK( 1 ) = DBLE( 1 ) 329: RETURN 330: END IF 331: * 332: * Initialize Q and Z 333: * 334: IF( ICOMPQ.EQ.3 ) 335: $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 336: IF( ICOMPZ.EQ.3 ) 337: $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) 338: * 339: * Machine Constants 340: * 341: IN = IHI + 1 - ILO 342: SAFMIN = DLAMCH( 'S' ) 343: SAFMAX = ONE / SAFMIN 344: ULP = DLAMCH( 'E' )*DLAMCH( 'B' ) 345: ANORM = DLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK ) 346: BNORM = DLANHS( 'F', IN, T( ILO, ILO ), LDT, WORK ) 347: ATOL = MAX( SAFMIN, ULP*ANORM ) 348: BTOL = MAX( SAFMIN, ULP*BNORM ) 349: ASCALE = ONE / MAX( SAFMIN, ANORM ) 350: BSCALE = ONE / MAX( SAFMIN, BNORM ) 351: * 352: * Set Eigenvalues IHI+1:N 353: * 354: DO 30 J = IHI + 1, N 355: IF( T( J, J ).LT.ZERO ) THEN 356: IF( ILSCHR ) THEN 357: DO 10 JR = 1, J 358: H( JR, J ) = -H( JR, J ) 359: T( JR, J ) = -T( JR, J ) 360: 10 CONTINUE 361: ELSE 362: H( J, J ) = -H( J, J ) 363: T( J, J ) = -T( J, J ) 364: END IF 365: IF( ILZ ) THEN 366: DO 20 JR = 1, N 367: Z( JR, J ) = -Z( JR, J ) 368: 20 CONTINUE 369: END IF 370: END IF 371: ALPHAR( J ) = H( J, J ) 372: ALPHAI( J ) = ZERO 373: BETA( J ) = T( J, J ) 374: 30 CONTINUE 375: * 376: * If IHI < ILO, skip QZ steps 377: * 378: IF( IHI.LT.ILO ) 379: $ GO TO 380 380: * 381: * MAIN QZ ITERATION LOOP 382: * 383: * Initialize dynamic indices 384: * 385: * Eigenvalues ILAST+1:N have been found. 386: * Column operations modify rows IFRSTM:whatever. 387: * Row operations modify columns whatever:ILASTM. 388: * 389: * If only eigenvalues are being computed, then 390: * IFRSTM is the row of the last splitting row above row ILAST; 391: * this is always at least ILO. 392: * IITER counts iterations since the last eigenvalue was found, 393: * to tell when to use an extraordinary shift. 394: * MAXIT is the maximum number of QZ sweeps allowed. 395: * 396: ILAST = IHI 397: IF( ILSCHR ) THEN 398: IFRSTM = 1 399: ILASTM = N 400: ELSE 401: IFRSTM = ILO 402: ILASTM = IHI 403: END IF 404: IITER = 0 405: ESHIFT = ZERO 406: MAXIT = 30*( IHI-ILO+1 ) 407: * 408: DO 360 JITER = 1, MAXIT 409: * 410: * Split the matrix if possible. 411: * 412: * Two tests: 413: * 1: H(j,j-1)=0 or j=ILO 414: * 2: T(j,j)=0 415: * 416: IF( ILAST.EQ.ILO ) THEN 417: * 418: * Special case: j=ILAST 419: * 420: GO TO 80 421: ELSE 422: IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN 423: H( ILAST, ILAST-1 ) = ZERO 424: GO TO 80 425: END IF 426: END IF 427: * 428: IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN 429: T( ILAST, ILAST ) = ZERO 430: GO TO 70 431: END IF 432: * 433: * General case: j<ILAST 434: * 435: DO 60 J = ILAST - 1, ILO, -1 436: * 437: * Test 1: for H(j,j-1)=0 or j=ILO 438: * 439: IF( J.EQ.ILO ) THEN 440: ILAZRO = .TRUE. 441: ELSE 442: IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN 443: H( J, J-1 ) = ZERO 444: ILAZRO = .TRUE. 445: ELSE 446: ILAZRO = .FALSE. 447: END IF 448: END IF 449: * 450: * Test 2: for T(j,j)=0 451: * 452: IF( ABS( T( J, J ) ).LT.BTOL ) THEN 453: T( J, J ) = ZERO 454: * 455: * Test 1a: Check for 2 consecutive small subdiagonals in A 456: * 457: ILAZR2 = .FALSE. 458: IF( .NOT.ILAZRO ) THEN 459: TEMP = ABS( H( J, J-1 ) ) 460: TEMP2 = ABS( H( J, J ) ) 461: TEMPR = MAX( TEMP, TEMP2 ) 462: IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 463: TEMP = TEMP / TEMPR 464: TEMP2 = TEMP2 / TEMPR 465: END IF 466: IF( TEMP*( ASCALE*ABS( H( J+1, J ) ) ).LE.TEMP2* 467: $ ( ASCALE*ATOL ) )ILAZR2 = .TRUE. 468: END IF 469: * 470: * If both tests pass (1 & 2), i.e., the leading diagonal 471: * element of B in the block is zero, split a 1x1 block off 472: * at the top. (I.e., at the J-th row/column) The leading 473: * diagonal element of the remainder can also be zero, so 474: * this may have to be done repeatedly. 475: * 476: IF( ILAZRO .OR. ILAZR2 ) THEN 477: DO 40 JCH = J, ILAST - 1 478: TEMP = H( JCH, JCH ) 479: CALL DLARTG( TEMP, H( JCH+1, JCH ), C, S, 480: $ H( JCH, JCH ) ) 481: H( JCH+1, JCH ) = ZERO 482: CALL DROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH, 483: $ H( JCH+1, JCH+1 ), LDH, C, S ) 484: CALL DROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT, 485: $ T( JCH+1, JCH+1 ), LDT, C, S ) 486: IF( ILQ ) 487: $ CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 488: $ C, S ) 489: IF( ILAZR2 ) 490: $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C 491: ILAZR2 = .FALSE. 492: IF( ABS( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN 493: IF( JCH+1.GE.ILAST ) THEN 494: GO TO 80 495: ELSE 496: IFIRST = JCH + 1 497: GO TO 110 498: END IF 499: END IF 500: T( JCH+1, JCH+1 ) = ZERO 501: 40 CONTINUE 502: GO TO 70 503: ELSE 504: * 505: * Only test 2 passed -- chase the zero to T(ILAST,ILAST) 506: * Then process as in the case T(ILAST,ILAST)=0 507: * 508: DO 50 JCH = J, ILAST - 1 509: TEMP = T( JCH, JCH+1 ) 510: CALL DLARTG( TEMP, T( JCH+1, JCH+1 ), C, S, 511: $ T( JCH, JCH+1 ) ) 512: T( JCH+1, JCH+1 ) = ZERO 513: IF( JCH.LT.ILASTM-1 ) 514: $ CALL DROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT, 515: $ T( JCH+1, JCH+2 ), LDT, C, S ) 516: CALL DROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH, 517: $ H( JCH+1, JCH-1 ), LDH, C, S ) 518: IF( ILQ ) 519: $ CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 520: $ C, S ) 521: TEMP = H( JCH+1, JCH ) 522: CALL DLARTG( TEMP, H( JCH+1, JCH-1 ), C, S, 523: $ H( JCH+1, JCH ) ) 524: H( JCH+1, JCH-1 ) = ZERO 525: CALL DROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1, 526: $ H( IFRSTM, JCH-1 ), 1, C, S ) 527: CALL DROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1, 528: $ T( IFRSTM, JCH-1 ), 1, C, S ) 529: IF( ILZ ) 530: $ CALL DROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1, 531: $ C, S ) 532: 50 CONTINUE 533: GO TO 70 534: END IF 535: ELSE IF( ILAZRO ) THEN 536: * 537: * Only test 1 passed -- work on J:ILAST 538: * 539: IFIRST = J 540: GO TO 110 541: END IF 542: * 543: * Neither test passed -- try next J 544: * 545: 60 CONTINUE 546: * 547: * (Drop-through is "impossible") 548: * 549: INFO = N + 1 550: GO TO 420 551: * 552: * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a 553: * 1x1 block. 554: * 555: 70 CONTINUE 556: TEMP = H( ILAST, ILAST ) 557: CALL DLARTG( TEMP, H( ILAST, ILAST-1 ), C, S, 558: $ H( ILAST, ILAST ) ) 559: H( ILAST, ILAST-1 ) = ZERO 560: CALL DROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1, 561: $ H( IFRSTM, ILAST-1 ), 1, C, S ) 562: CALL DROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1, 563: $ T( IFRSTM, ILAST-1 ), 1, C, S ) 564: IF( ILZ ) 565: $ CALL DROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S ) 566: * 567: * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI, 568: * and BETA 569: * 570: 80 CONTINUE 571: IF( T( ILAST, ILAST ).LT.ZERO ) THEN 572: IF( ILSCHR ) THEN 573: DO 90 J = IFRSTM, ILAST 574: H( J, ILAST ) = -H( J, ILAST ) 575: T( J, ILAST ) = -T( J, ILAST ) 576: 90 CONTINUE 577: ELSE 578: H( ILAST, ILAST ) = -H( ILAST, ILAST ) 579: T( ILAST, ILAST ) = -T( ILAST, ILAST ) 580: END IF 581: IF( ILZ ) THEN 582: DO 100 J = 1, N 583: Z( J, ILAST ) = -Z( J, ILAST ) 584: 100 CONTINUE 585: END IF 586: END IF 587: ALPHAR( ILAST ) = H( ILAST, ILAST ) 588: ALPHAI( ILAST ) = ZERO 589: BETA( ILAST ) = T( ILAST, ILAST ) 590: * 591: * Go to next block -- exit if finished. 592: * 593: ILAST = ILAST - 1 594: IF( ILAST.LT.ILO ) 595: $ GO TO 380 596: * 597: * Reset counters 598: * 599: IITER = 0 600: ESHIFT = ZERO 601: IF( .NOT.ILSCHR ) THEN 602: ILASTM = ILAST 603: IF( IFRSTM.GT.ILAST ) 604: $ IFRSTM = ILO 605: END IF 606: GO TO 350 607: * 608: * QZ step 609: * 610: * This iteration only involves rows/columns IFIRST:ILAST. We 611: * assume IFIRST < ILAST, and that the diagonal of B is non-zero. 612: * 613: 110 CONTINUE 614: IITER = IITER + 1 615: IF( .NOT.ILSCHR ) THEN 616: IFRSTM = IFIRST 617: END IF 618: * 619: * Compute single shifts. 620: * 621: * At this point, IFIRST < ILAST, and the diagonal elements of 622: * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in 623: * magnitude) 624: * 625: IF( ( IITER / 10 )*10.EQ.IITER ) THEN 626: * 627: * Exceptional shift. Chosen for no particularly good reason. 628: * (Single shift only.) 629: * 630: IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT. 631: $ ABS( T( ILAST-1, ILAST-1 ) ) ) THEN 632: ESHIFT = ESHIFT + H( ILAST-1, ILAST ) / 633: $ T( ILAST-1, ILAST-1 ) 634: ELSE 635: ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) ) 636: END IF 637: S1 = ONE 638: WR = ESHIFT 639: * 640: ELSE 641: * 642: * Shifts based on the generalized eigenvalues of the 643: * bottom-right 2x2 block of A and B. The first eigenvalue 644: * returned by DLAG2 is the Wilkinson shift (AEP p.512), 645: * 646: CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH, 647: $ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1, 648: $ S2, WR, WR2, WI ) 649: * 650: TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) ) 651: IF( WI.NE.ZERO ) 652: $ GO TO 200 653: END IF 654: * 655: * Fiddle with shift to avoid overflow 656: * 657: TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX ) 658: IF( S1.GT.TEMP ) THEN 659: SCALE = TEMP / S1 660: ELSE 661: SCALE = ONE 662: END IF 663: * 664: TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX ) 665: IF( ABS( WR ).GT.TEMP ) 666: $ SCALE = MIN( SCALE, TEMP / ABS( WR ) ) 667: S1 = SCALE*S1 668: WR = SCALE*WR 669: * 670: * Now check for two consecutive small subdiagonals. 671: * 672: DO 120 J = ILAST - 1, IFIRST + 1, -1 673: ISTART = J 674: TEMP = ABS( S1*H( J, J-1 ) ) 675: TEMP2 = ABS( S1*H( J, J )-WR*T( J, J ) ) 676: TEMPR = MAX( TEMP, TEMP2 ) 677: IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 678: TEMP = TEMP / TEMPR 679: TEMP2 = TEMP2 / TEMPR 680: END IF 681: IF( ABS( ( ASCALE*H( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )* 682: $ TEMP2 )GO TO 130 683: 120 CONTINUE 684: * 685: ISTART = IFIRST 686: 130 CONTINUE 687: * 688: * Do an implicit single-shift QZ sweep. 689: * 690: * Initial Q 691: * 692: TEMP = S1*H( ISTART, ISTART ) - WR*T( ISTART, ISTART ) 693: TEMP2 = S1*H( ISTART+1, ISTART ) 694: CALL DLARTG( TEMP, TEMP2, C, S, TEMPR ) 695: * 696: * Sweep 697: * 698: DO 190 J = ISTART, ILAST - 1 699: IF( J.GT.ISTART ) THEN 700: TEMP = H( J, J-1 ) 701: CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) 702: H( J+1, J-1 ) = ZERO 703: END IF 704: * 705: DO 140 JC = J, ILASTM 706: TEMP = C*H( J, JC ) + S*H( J+1, JC ) 707: H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC ) 708: H( J, JC ) = TEMP 709: TEMP2 = C*T( J, JC ) + S*T( J+1, JC ) 710: T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC ) 711: T( J, JC ) = TEMP2 712: 140 CONTINUE 713: IF( ILQ ) THEN 714: DO 150 JR = 1, N 715: TEMP = C*Q( JR, J ) + S*Q( JR, J+1 ) 716: Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 717: Q( JR, J ) = TEMP 718: 150 CONTINUE 719: END IF 720: * 721: TEMP = T( J+1, J+1 ) 722: CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) 723: T( J+1, J ) = ZERO 724: * 725: DO 160 JR = IFRSTM, MIN( J+2, ILAST ) 726: TEMP = C*H( JR, J+1 ) + S*H( JR, J ) 727: H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J ) 728: H( JR, J+1 ) = TEMP 729: 160 CONTINUE 730: DO 170 JR = IFRSTM, J 731: TEMP = C*T( JR, J+1 ) + S*T( JR, J ) 732: T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J ) 733: T( JR, J+1 ) = TEMP 734: 170 CONTINUE 735: IF( ILZ ) THEN 736: DO 180 JR = 1, N 737: TEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 738: Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J ) 739: Z( JR, J+1 ) = TEMP 740: 180 CONTINUE 741: END IF 742: 190 CONTINUE 743: * 744: GO TO 350 745: * 746: * Use Francis double-shift 747: * 748: * Note: the Francis double-shift should work with real shifts, 749: * but only if the block is at least 3x3. 750: * This code may break if this point is reached with 751: * a 2x2 block with real eigenvalues. 752: * 753: 200 CONTINUE 754: IF( IFIRST+1.EQ.ILAST ) THEN 755: * 756: * Special case -- 2x2 block with complex eigenvectors 757: * 758: * Step 1: Standardize, that is, rotate so that 759: * 760: * ( B11 0 ) 761: * B = ( ) with B11 non-negative. 762: * ( 0 B22 ) 763: * 764: CALL DLASV2( T( ILAST-1, ILAST-1 ), T( ILAST-1, ILAST ), 765: $ T( ILAST, ILAST ), B22, B11, SR, CR, SL, CL ) 766: * 767: IF( B11.LT.ZERO ) THEN 768: CR = -CR 769: SR = -SR 770: B11 = -B11 771: B22 = -B22 772: END IF 773: * 774: CALL DROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH, 775: $ H( ILAST, ILAST-1 ), LDH, CL, SL ) 776: CALL DROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1, 777: $ H( IFRSTM, ILAST ), 1, CR, SR ) 778: * 779: IF( ILAST.LT.ILASTM ) 780: $ CALL DROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT, 781: $ T( ILAST, ILAST+1 ), LDT, CL, SL ) 782: IF( IFRSTM.LT.ILAST-1 ) 783: $ CALL DROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1, 784: $ T( IFRSTM, ILAST ), 1, CR, SR ) 785: * 786: IF( ILQ ) 787: $ CALL DROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL, 788: $ SL ) 789: IF( ILZ ) 790: $ CALL DROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR, 791: $ SR ) 792: * 793: T( ILAST-1, ILAST-1 ) = B11 794: T( ILAST-1, ILAST ) = ZERO 795: T( ILAST, ILAST-1 ) = ZERO 796: T( ILAST, ILAST ) = B22 797: * 798: * If B22 is negative, negate column ILAST 799: * 800: IF( B22.LT.ZERO ) THEN 801: DO 210 J = IFRSTM, ILAST 802: H( J, ILAST ) = -H( J, ILAST ) 803: T( J, ILAST ) = -T( J, ILAST ) 804: 210 CONTINUE 805: * 806: IF( ILZ ) THEN 807: DO 220 J = 1, N 808: Z( J, ILAST ) = -Z( J, ILAST ) 809: 220 CONTINUE 810: END IF 811: END IF 812: * 813: * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.) 814: * 815: * Recompute shift 816: * 817: CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH, 818: $ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1, 819: $ TEMP, WR, TEMP2, WI ) 820: * 821: * If standardization has perturbed the shift onto real line, 822: * do another (real single-shift) QR step. 823: * 824: IF( WI.EQ.ZERO ) 825: $ GO TO 350 826: S1INV = ONE / S1 827: * 828: * Do EISPACK (QZVAL) computation of alpha and beta 829: * 830: A11 = H( ILAST-1, ILAST-1 ) 831: A21 = H( ILAST, ILAST-1 ) 832: A12 = H( ILAST-1, ILAST ) 833: A22 = H( ILAST, ILAST ) 834: * 835: * Compute complex Givens rotation on right 836: * (Assume some element of C = (sA - wB) > unfl ) 837: * __ 838: * (sA - wB) ( CZ -SZ ) 839: * ( SZ CZ ) 840: * 841: C11R = S1*A11 - WR*B11 842: C11I = -WI*B11 843: C12 = S1*A12 844: C21 = S1*A21 845: C22R = S1*A22 - WR*B22 846: C22I = -WI*B22 847: * 848: IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+ 849: $ ABS( C22R )+ABS( C22I ) ) THEN 850: T1 = DLAPY3( C12, C11R, C11I ) 851: CZ = C12 / T1 852: SZR = -C11R / T1 853: SZI = -C11I / T1 854: ELSE 855: CZ = DLAPY2( C22R, C22I ) 856: IF( CZ.LE.SAFMIN ) THEN 857: CZ = ZERO 858: SZR = ONE 859: SZI = ZERO 860: ELSE 861: TEMPR = C22R / CZ 862: TEMPI = C22I / CZ 863: T1 = DLAPY2( CZ, C21 ) 864: CZ = CZ / T1 865: SZR = -C21*TEMPR / T1 866: SZI = C21*TEMPI / T1 867: END IF 868: END IF 869: * 870: * Compute Givens rotation on left 871: * 872: * ( CQ SQ ) 873: * ( __ ) A or B 874: * ( -SQ CQ ) 875: * 876: AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 ) 877: BN = ABS( B11 ) + ABS( B22 ) 878: WABS = ABS( WR ) + ABS( WI ) 879: IF( S1*AN.GT.WABS*BN ) THEN 880: CQ = CZ*B11 881: SQR = SZR*B22 882: SQI = -SZI*B22 883: ELSE 884: A1R = CZ*A11 + SZR*A12 885: A1I = SZI*A12 886: A2R = CZ*A21 + SZR*A22 887: A2I = SZI*A22 888: CQ = DLAPY2( A1R, A1I ) 889: IF( CQ.LE.SAFMIN ) THEN 890: CQ = ZERO 891: SQR = ONE 892: SQI = ZERO 893: ELSE 894: TEMPR = A1R / CQ 895: TEMPI = A1I / CQ 896: SQR = TEMPR*A2R + TEMPI*A2I 897: SQI = TEMPI*A2R - TEMPR*A2I 898: END IF 899: END IF 900: T1 = DLAPY3( CQ, SQR, SQI ) 901: CQ = CQ / T1 902: SQR = SQR / T1 903: SQI = SQI / T1 904: * 905: * Compute diagonal elements of QBZ 906: * 907: TEMPR = SQR*SZR - SQI*SZI 908: TEMPI = SQR*SZI + SQI*SZR 909: B1R = CQ*CZ*B11 + TEMPR*B22 910: B1I = TEMPI*B22 911: B1A = DLAPY2( B1R, B1I ) 912: B2R = CQ*CZ*B22 + TEMPR*B11 913: B2I = -TEMPI*B11 914: B2A = DLAPY2( B2R, B2I ) 915: * 916: * Normalize so beta > 0, and Im( alpha1 ) > 0 917: * 918: BETA( ILAST-1 ) = B1A 919: BETA( ILAST ) = B2A 920: ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV 921: ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV 922: ALPHAR( ILAST ) = ( WR*B2A )*S1INV 923: ALPHAI( ILAST ) = -( WI*B2A )*S1INV 924: * 925: * Step 3: Go to next block -- exit if finished. 926: * 927: ILAST = IFIRST - 1 928: IF( ILAST.LT.ILO ) 929: $ GO TO 380 930: * 931: * Reset counters 932: * 933: IITER = 0 934: ESHIFT = ZERO 935: IF( .NOT.ILSCHR ) THEN 936: ILASTM = ILAST 937: IF( IFRSTM.GT.ILAST ) 938: $ IFRSTM = ILO 939: END IF 940: GO TO 350 941: ELSE 942: * 943: * Usual case: 3x3 or larger block, using Francis implicit 944: * double-shift 945: * 946: * 2 947: * Eigenvalue equation is w - c w + d = 0, 948: * 949: * -1 2 -1 950: * so compute 1st column of (A B ) - c A B + d 951: * using the formula in QZIT (from EISPACK) 952: * 953: * We assume that the block is at least 3x3 954: * 955: AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) / 956: $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 957: AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) / 958: $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 959: AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) / 960: $ ( BSCALE*T( ILAST, ILAST ) ) 961: AD22 = ( ASCALE*H( ILAST, ILAST ) ) / 962: $ ( BSCALE*T( ILAST, ILAST ) ) 963: U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST ) 964: AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) / 965: $ ( BSCALE*T( IFIRST, IFIRST ) ) 966: AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) / 967: $ ( BSCALE*T( IFIRST, IFIRST ) ) 968: AD12L = ( ASCALE*H( IFIRST, IFIRST+1 ) ) / 969: $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) ) 970: AD22L = ( ASCALE*H( IFIRST+1, IFIRST+1 ) ) / 971: $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) ) 972: AD32L = ( ASCALE*H( IFIRST+2, IFIRST+1 ) ) / 973: $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) ) 974: U12L = T( IFIRST, IFIRST+1 ) / T( IFIRST+1, IFIRST+1 ) 975: * 976: V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 + 977: $ AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L 978: V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )- 979: $ ( AD22-AD11L )+AD21*U12 )*AD21L 980: V( 3 ) = AD32L*AD21L 981: * 982: ISTART = IFIRST 983: * 984: CALL DLARFG( 3, V( 1 ), V( 2 ), 1, TAU ) 985: V( 1 ) = ONE 986: * 987: * Sweep 988: * 989: DO 290 J = ISTART, ILAST - 2 990: * 991: * All but last elements: use 3x3 Householder transforms. 992: * 993: * Zero (j-1)st column of A 994: * 995: IF( J.GT.ISTART ) THEN 996: V( 1 ) = H( J, J-1 ) 997: V( 2 ) = H( J+1, J-1 ) 998: V( 3 ) = H( J+2, J-1 ) 999: * 1000: CALL DLARFG( 3, H( J, J-1 ), V( 2 ), 1, TAU ) 1001: V( 1 ) = ONE 1002: H( J+1, J-1 ) = ZERO 1003: H( J+2, J-1 ) = ZERO 1004: END IF 1005: * 1006: DO 230 JC = J, ILASTM 1007: TEMP = TAU*( H( J, JC )+V( 2 )*H( J+1, JC )+V( 3 )* 1008: $ H( J+2, JC ) ) 1009: H( J, JC ) = H( J, JC ) - TEMP 1010: H( J+1, JC ) = H( J+1, JC ) - TEMP*V( 2 ) 1011: H( J+2, JC ) = H( J+2, JC ) - TEMP*V( 3 ) 1012: TEMP2 = TAU*( T( J, JC )+V( 2 )*T( J+1, JC )+V( 3 )* 1013: $ T( J+2, JC ) ) 1014: T( J, JC ) = T( J, JC ) - TEMP2 1015: T( J+1, JC ) = T( J+1, JC ) - TEMP2*V( 2 ) 1016: T( J+2, JC ) = T( J+2, JC ) - TEMP2*V( 3 ) 1017: 230 CONTINUE 1018: IF( ILQ ) THEN 1019: DO 240 JR = 1, N 1020: TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )* 1021: $ Q( JR, J+2 ) ) 1022: Q( JR, J ) = Q( JR, J ) - TEMP 1023: Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 ) 1024: Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 ) 1025: 240 CONTINUE 1026: END IF 1027: * 1028: * Zero j-th column of B (see DLAGBC for details) 1029: * 1030: * Swap rows to pivot 1031: * 1032: ILPIVT = .FALSE. 1033: TEMP = MAX( ABS( T( J+1, J+1 ) ), ABS( T( J+1, J+2 ) ) ) 1034: TEMP2 = MAX( ABS( T( J+2, J+1 ) ), ABS( T( J+2, J+2 ) ) ) 1035: IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN 1036: SCALE = ZERO 1037: U1 = ONE 1038: U2 = ZERO 1039: GO TO 250 1040: ELSE IF( TEMP.GE.TEMP2 ) THEN 1041: W11 = T( J+1, J+1 ) 1042: W21 = T( J+2, J+1 ) 1043: W12 = T( J+1, J+2 ) 1044: W22 = T( J+2, J+2 ) 1045: U1 = T( J+1, J ) 1046: U2 = T( J+2, J ) 1047: ELSE 1048: W21 = T( J+1, J+1 ) 1049: W11 = T( J+2, J+1 ) 1050: W22 = T( J+1, J+2 ) 1051: W12 = T( J+2, J+2 ) 1052: U2 = T( J+1, J ) 1053: U1 = T( J+2, J ) 1054: END IF 1055: * 1056: * Swap columns if nec. 1057: * 1058: IF( ABS( W12 ).GT.ABS( W11 ) ) THEN 1059: ILPIVT = .TRUE. 1060: TEMP = W12 1061: TEMP2 = W22 1062: W12 = W11 1063: W22 = W21 1064: W11 = TEMP 1065: W21 = TEMP2 1066: END IF 1067: * 1068: * LU-factor 1069: * 1070: TEMP = W21 / W11 1071: U2 = U2 - TEMP*U1 1072: W22 = W22 - TEMP*W12 1073: W21 = ZERO 1074: * 1075: * Compute SCALE 1076: * 1077: SCALE = ONE 1078: IF( ABS( W22 ).LT.SAFMIN ) THEN 1079: SCALE = ZERO 1080: U2 = ONE 1081: U1 = -W12 / W11 1082: GO TO 250 1083: END IF 1084: IF( ABS( W22 ).LT.ABS( U2 ) ) 1085: $ SCALE = ABS( W22 / U2 ) 1086: IF( ABS( W11 ).LT.ABS( U1 ) ) 1087: $ SCALE = MIN( SCALE, ABS( W11 / U1 ) ) 1088: * 1089: * Solve 1090: * 1091: U2 = ( SCALE*U2 ) / W22 1092: U1 = ( SCALE*U1-W12*U2 ) / W11 1093: * 1094: 250 CONTINUE 1095: IF( ILPIVT ) THEN 1096: TEMP = U2 1097: U2 = U1 1098: U1 = TEMP 1099: END IF 1100: * 1101: * Compute Householder Vector 1102: * 1103: T1 = SQRT( SCALE**2+U1**2+U2**2 ) 1104: TAU = ONE + SCALE / T1 1105: VS = -ONE / ( SCALE+T1 ) 1106: V( 1 ) = ONE 1107: V( 2 ) = VS*U1 1108: V( 3 ) = VS*U2 1109: * 1110: * Apply transformations from the right. 1111: * 1112: DO 260 JR = IFRSTM, MIN( J+3, ILAST ) 1113: TEMP = TAU*( H( JR, J )+V( 2 )*H( JR, J+1 )+V( 3 )* 1114: $ H( JR, J+2 ) ) 1115: H( JR, J ) = H( JR, J ) - TEMP 1116: H( JR, J+1 ) = H( JR, J+1 ) - TEMP*V( 2 ) 1117: H( JR, J+2 ) = H( JR, J+2 ) - TEMP*V( 3 ) 1118: 260 CONTINUE 1119: DO 270 JR = IFRSTM, J + 2 1120: TEMP = TAU*( T( JR, J )+V( 2 )*T( JR, J+1 )+V( 3 )* 1121: $ T( JR, J+2 ) ) 1122: T( JR, J ) = T( JR, J ) - TEMP 1123: T( JR, J+1 ) = T( JR, J+1 ) - TEMP*V( 2 ) 1124: T( JR, J+2 ) = T( JR, J+2 ) - TEMP*V( 3 ) 1125: 270 CONTINUE 1126: IF( ILZ ) THEN 1127: DO 280 JR = 1, N 1128: TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )* 1129: $ Z( JR, J+2 ) ) 1130: Z( JR, J ) = Z( JR, J ) - TEMP 1131: Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 ) 1132: Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 ) 1133: 280 CONTINUE 1134: END IF 1135: T( J+1, J ) = ZERO 1136: T( J+2, J ) = ZERO 1137: 290 CONTINUE 1138: * 1139: * Last elements: Use Givens rotations 1140: * 1141: * Rotations from the left 1142: * 1143: J = ILAST - 1 1144: TEMP = H( J, J-1 ) 1145: CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) 1146: H( J+1, J-1 ) = ZERO 1147: * 1148: DO 300 JC = J, ILASTM 1149: TEMP = C*H( J, JC ) + S*H( J+1, JC ) 1150: H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC ) 1151: H( J, JC ) = TEMP 1152: TEMP2 = C*T( J, JC ) + S*T( J+1, JC ) 1153: T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC ) 1154: T( J, JC ) = TEMP2 1155: 300 CONTINUE 1156: IF( ILQ ) THEN 1157: DO 310 JR = 1, N 1158: TEMP = C*Q( JR, J ) + S*Q( JR, J+1 ) 1159: Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 1160: Q( JR, J ) = TEMP 1161: 310 CONTINUE 1162: END IF 1163: * 1164: * Rotations from the right. 1165: * 1166: TEMP = T( J+1, J+1 ) 1167: CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) 1168: T( J+1, J ) = ZERO 1169: * 1170: DO 320 JR = IFRSTM, ILAST 1171: TEMP = C*H( JR, J+1 ) + S*H( JR, J ) 1172: H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J ) 1173: H( JR, J+1 ) = TEMP 1174: 320 CONTINUE 1175: DO 330 JR = IFRSTM, ILAST - 1 1176: TEMP = C*T( JR, J+1 ) + S*T( JR, J ) 1177: T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J ) 1178: T( JR, J+1 ) = TEMP 1179: 330 CONTINUE 1180: IF( ILZ ) THEN 1181: DO 340 JR = 1, N 1182: TEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 1183: Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J ) 1184: Z( JR, J+1 ) = TEMP 1185: 340 CONTINUE 1186: END IF 1187: * 1188: * End of Double-Shift code 1189: * 1190: END IF 1191: * 1192: GO TO 350 1193: * 1194: * End of iteration loop 1195: * 1196: 350 CONTINUE 1197: 360 CONTINUE 1198: * 1199: * Drop-through = non-convergence 1200: * 1201: INFO = ILAST 1202: GO TO 420 1203: * 1204: * Successful completion of all QZ steps 1205: * 1206: 380 CONTINUE 1207: * 1208: * Set Eigenvalues 1:ILO-1 1209: * 1210: DO 410 J = 1, ILO - 1 1211: IF( T( J, J ).LT.ZERO ) THEN 1212: IF( ILSCHR ) THEN 1213: DO 390 JR = 1, J 1214: H( JR, J ) = -H( JR, J ) 1215: T( JR, J ) = -T( JR, J ) 1216: 390 CONTINUE 1217: ELSE 1218: H( J, J ) = -H( J, J ) 1219: T( J, J ) = -T( J, J ) 1220: END IF 1221: IF( ILZ ) THEN 1222: DO 400 JR = 1, N 1223: Z( JR, J ) = -Z( JR, J ) 1224: 400 CONTINUE 1225: END IF 1226: END IF 1227: ALPHAR( J ) = H( J, J ) 1228: ALPHAI( J ) = ZERO 1229: BETA( J ) = T( J, J ) 1230: 410 CONTINUE 1231: * 1232: * Normal Termination 1233: * 1234: INFO = 0 1235: * 1236: * Exit (other than argument error) -- return optimal workspace size 1237: * 1238: 420 CONTINUE 1239: WORK( 1 ) = DBLE( N ) 1240: RETURN 1241: * 1242: * End of DHGEQZ 1243: * 1244: END