SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, $ ILOZ, IHIZ, Z, LDZ, INFO ) * * -- LAPACK auxiliary routine (instrum. to count ops. version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. LOGICAL WANTT, WANTZ INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N * .. * .. Array Arguments .. DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) * .. * Common block to return operation count. * .. Common blocks .. COMMON / LATIME / OPS, ITCNT * .. * .. Scalars in Common .. DOUBLE PRECISION ITCNT, OPS * .. * * Purpose * ======= * * DLAHQR is an auxiliary routine called by DHSEQR to update the * eigenvalues and Schur decomposition already computed by DHSEQR, by * dealing with the Hessenberg submatrix in rows and columns ILO to IHI. * * Arguments * ========= * * WANTT (input) LOGICAL * = .TRUE. : the full Schur form T is required; * = .FALSE.: only eigenvalues are required. * * WANTZ (input) LOGICAL * = .TRUE. : the matrix of Schur vectors Z is required; * = .FALSE.: Schur vectors are not required. * * N (input) INTEGER * The order of the matrix H. N >= 0. * * ILO (input) INTEGER * IHI (input) INTEGER * It is assumed that H is already upper quasi-triangular in * rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless * ILO = 1). DLAHQR works primarily with the Hessenberg * submatrix in rows and columns ILO to IHI, but applies * transformations to all of H if WANTT is .TRUE.. * 1 <= ILO <= max(1,IHI); IHI <= N. * * H (input/output) DOUBLE PRECISION array, dimension (LDH,N) * On entry, the upper Hessenberg matrix H. * On exit, if WANTT is .TRUE., H is upper quasi-triangular in * rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in * standard form. If WANTT is .FALSE., the contents of H are * unspecified on exit. * * LDH (input) INTEGER * The leading dimension of the array H. LDH >= max(1,N). * * WR (output) DOUBLE PRECISION array, dimension (N) * WI (output) DOUBLE PRECISION array, dimension (N) * The real and imaginary parts, respectively, of the computed * eigenvalues ILO to IHI are stored in the corresponding * elements of WR and WI. If two eigenvalues are computed as a * complex conjugate pair, they are stored in consecutive * elements of WR and WI, say the i-th and (i+1)th, with * WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the * eigenvalues are stored in the same order as on the diagonal * of the Schur form returned in H, with WR(i) = H(i,i), and, if * H(i:i+1,i:i+1) is a 2-by-2 diagonal block, * WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). * * ILOZ (input) INTEGER * IHIZ (input) INTEGER * Specify the rows of Z to which transformations must be * applied if WANTZ is .TRUE.. * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. * * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) * If WANTZ is .TRUE., on entry Z must contain the current * matrix Z of transformations accumulated by DHSEQR, and on * exit Z has been updated; transformations are applied only to * the submatrix Z(ILOZ:IHIZ,ILO:IHI). * If WANTZ is .FALSE., Z is not referenced. * * LDZ (input) INTEGER * The leading dimension of the array Z. LDZ >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * > 0: DLAHQR failed to compute all the eigenvalues ILO to IHI * in a total of 30*(IHI-ILO+1) iterations; if INFO = i, * elements i+1:ihi of WR and WI contain those eigenvalues * which have been successfully computed. * * Further Details * =============== * * 2-96 Based on modifications by * David Day, Sandia National Laboratory, USA * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, HALF PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D0 ) DOUBLE PRECISION DAT1, DAT2 PARAMETER ( DAT1 = 0.75D+0, DAT2 = -0.4375D+0 ) * .. * .. Local Scalars .. INTEGER I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ DOUBLE PRECISION AVE, CS, DISC, H00, H10, H11, H12, H21, H22, $ H33, H33S, H43H34, H44, H44S, OPST, OVFL, S, $ SMLNUM, SN, SUM, T1, T2, T3, TST1, ULP, UNFL, $ V1, V2, V3 * .. * .. Local Arrays .. DOUBLE PRECISION V( 3 ), WORK( 1 ) * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DLANHS EXTERNAL DLAMCH, DLANHS * .. * .. External Subroutines .. EXTERNAL DCOPY, DLABAD, DLANV2, DLARFG, DROT * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SIGN, SQRT * .. * .. Executable Statements .. * INFO = 0 *** * Initialize OPST = 0 *** * * Quick return if possible * IF( N.EQ.0 ) $ RETURN IF( ILO.EQ.IHI ) THEN WR( ILO ) = H( ILO, ILO ) WI( ILO ) = ZERO RETURN END IF * NH = IHI - ILO + 1 NZ = IHIZ - ILOZ + 1 * * Set machine-dependent constants for the stopping criterion. * If norm(H) <= sqrt(OVFL), overflow should not occur. * UNFL = DLAMCH( 'Safe minimum' ) OVFL = ONE / UNFL CALL DLABAD( UNFL, OVFL ) ULP = DLAMCH( 'Precision' ) SMLNUM = UNFL*( NH / ULP ) * * I1 and I2 are the indices of the first row and last column of H * to which transformations must be applied. If eigenvalues only are * being computed, I1 and I2 are set inside the main loop. * IF( WANTT ) THEN I1 = 1 I2 = N END IF * * ITN is the total number of QR iterations allowed. * ITN = 30*NH * * The main loop begins here. I is the loop index and decreases from * IHI to ILO in steps of 1 or 2. Each iteration of the loop works * with the active submatrix in rows and columns L to I. * Eigenvalues I+1 to IHI have already converged. Either L = ILO or * H(L,L-1) is negligible so that the matrix splits. * I = IHI 10 CONTINUE L = ILO IF( I.LT.ILO ) $ GO TO 150 * * Perform QR iterations on rows and columns ILO to I until a * submatrix of order 1 or 2 splits off at the bottom because a * subdiagonal element has become negligible. * DO 130 ITS = 0, ITN * * Look for a single small subdiagonal element. * DO 20 K = I, L + 1, -1 TST1 = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) ) IF( TST1.EQ.ZERO ) THEN TST1 = DLANHS( '1', I-L+1, H( L, L ), LDH, WORK ) *** * Increment op count OPS = OPS + ( I-L+1 )*( I-L+2 ) / 2 *** END IF IF( ABS( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) ) $ GO TO 30 20 CONTINUE 30 CONTINUE L = K *** * Increment op count OPST = OPST + 3*( I-L+1 ) *** IF( L.GT.ILO ) THEN * * H(L,L-1) is negligible * H( L, L-1 ) = ZERO END IF * * Exit from loop if a submatrix of order 1 or 2 has split off. * IF( L.GE.I-1 ) $ GO TO 140 * * Now the active submatrix is in rows and columns L to I. If * eigenvalues only are being computed, only the active submatrix * need be transformed. * IF( .NOT.WANTT ) THEN I1 = L I2 = I END IF * IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN * * Exceptional shift. * S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) ) H44 = DAT1*S + H( I, I ) H33 = H44 H43H34 = DAT2*S*S *** * Increment op count OPST = OPST + 5 *** ELSE * * Prepare to use Francis' double shift * (i.e. 2nd degree generalized Rayleigh quotient) * H44 = H( I, I ) H33 = H( I-1, I-1 ) H43H34 = H( I, I-1 )*H( I-1, I ) S = H( I-1, I-2 )*H( I-1, I-2 ) DISC = ( H33-H44 )*HALF DISC = DISC*DISC + H43H34 *** * Increment op count OPST = OPST + 6 *** IF( DISC.GT.ZERO ) THEN * * Real roots: use Wilkinson's shift twice * DISC = SQRT( DISC ) AVE = HALF*( H33+H44 ) *** * Increment op count OPST = OPST + 2 *** IF( ABS( H33 )-ABS( H44 ).GT.ZERO ) THEN H33 = H33*H44 - H43H34 H44 = H33 / ( SIGN( DISC, AVE )+AVE ) *** * Increment op count OPST = OPST + 4 *** ELSE H44 = SIGN( DISC, AVE ) + AVE *** * Increment op count OPST = OPST + 1 *** END IF H33 = H44 H43H34 = ZERO END IF END IF * * Look for two consecutive small subdiagonal elements. * DO 40 M = I - 2, L, -1 * * Determine the effect of starting the double-shift QR * iteration at row M, and see if this would make H(M,M-1) * negligible. * H11 = H( M, M ) H22 = H( M+1, M+1 ) H21 = H( M+1, M ) H12 = H( M, M+1 ) H44S = H44 - H11 H33S = H33 - H11 V1 = ( H33S*H44S-H43H34 ) / H21 + H12 V2 = H22 - H11 - H33S - H44S V3 = H( M+2, M+1 ) S = ABS( V1 ) + ABS( V2 ) + ABS( V3 ) V1 = V1 / S V2 = V2 / S V3 = V3 / S V( 1 ) = V1 V( 2 ) = V2 V( 3 ) = V3 IF( M.EQ.L ) $ GO TO 50 H00 = H( M-1, M-1 ) H10 = H( M, M-1 ) TST1 = ABS( V1 )*( ABS( H00 )+ABS( H11 )+ABS( H22 ) ) IF( ABS( H10 )*( ABS( V2 )+ABS( V3 ) ).LE.ULP*TST1 ) $ GO TO 50 40 CONTINUE 50 CONTINUE *** * Increment op count OPST = OPST + 20*( I-M-1 ) *** * * Double-shift QR step * DO 120 K = M, I - 1 * * The first iteration of this loop determines a reflection G * from the vector V and applies it from left and right to H, * thus creating a nonzero bulge below the subdiagonal. * * Each subsequent iteration determines a reflection G to * restore the Hessenberg form in the (K-1)th column, and thus * chases the bulge one step toward the bottom of the active * submatrix. NR is the order of G. * NR = MIN( 3, I-K+1 ) IF( K.GT.M ) $ CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 ) CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 ) *** * Increment op count OPST = OPST + 3*NR + 9 *** IF( K.GT.M ) THEN H( K, K-1 ) = V( 1 ) H( K+1, K-1 ) = ZERO IF( K.LT.I-1 ) $ H( K+2, K-1 ) = ZERO ELSE IF( M.GT.L ) THEN H( K, K-1 ) = -H( K, K-1 ) END IF V2 = V( 2 ) T2 = T1*V2 IF( NR.EQ.3 ) THEN V3 = V( 3 ) T3 = T1*V3 * * Apply G from the left to transform the rows of the matrix * in columns K to I2. * DO 60 J = K, I2 SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J ) H( K, J ) = H( K, J ) - SUM*T1 H( K+1, J ) = H( K+1, J ) - SUM*T2 H( K+2, J ) = H( K+2, J ) - SUM*T3 60 CONTINUE * * Apply G from the right to transform the columns of the * matrix in rows I1 to min(K+3,I). * DO 70 J = I1, MIN( K+3, I ) SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 ) H( J, K ) = H( J, K ) - SUM*T1 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 H( J, K+2 ) = H( J, K+2 ) - SUM*T3 70 CONTINUE *** * Increment op count OPS = OPS + 10*( I2-I1+2+MIN( 3, I-K ) ) *** * IF( WANTZ ) THEN * * Accumulate transformations in the matrix Z * DO 80 J = ILOZ, IHIZ SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 ) Z( J, K ) = Z( J, K ) - SUM*T1 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3 80 CONTINUE *** * Increment op count OPS = OPS + 10*NZ *** END IF ELSE IF( NR.EQ.2 ) THEN * * Apply G from the left to transform the rows of the matrix * in columns K to I2. * DO 90 J = K, I2 SUM = H( K, J ) + V2*H( K+1, J ) H( K, J ) = H( K, J ) - SUM*T1 H( K+1, J ) = H( K+1, J ) - SUM*T2 90 CONTINUE * * Apply G from the right to transform the columns of the * matrix in rows I1 to min(K+3,I). * DO 100 J = I1, I SUM = H( J, K ) + V2*H( J, K+1 ) H( J, K ) = H( J, K ) - SUM*T1 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 100 CONTINUE *** * Increment op count OPS = OPS + 6*( I2-I1+3 ) *** * IF( WANTZ ) THEN * * Accumulate transformations in the matrix Z * DO 110 J = ILOZ, IHIZ SUM = Z( J, K ) + V2*Z( J, K+1 ) Z( J, K ) = Z( J, K ) - SUM*T1 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 110 CONTINUE *** * Increment op count OPS = OPS + 6*NZ *** END IF END IF 120 CONTINUE * 130 CONTINUE * * Failure to converge in remaining number of iterations * INFO = I RETURN * 140 CONTINUE * IF( L.EQ.I ) THEN * * H(I,I-1) is negligible: one eigenvalue has converged. * WR( I ) = H( I, I ) WI( I ) = ZERO ELSE IF( L.EQ.I-1 ) THEN * * H(I-1,I-2) is negligible: a pair of eigenvalues have converged. * * Transform the 2-by-2 submatrix to standard Schur form, * and compute and store the eigenvalues. * CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ), $ H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ), $ CS, SN ) * IF( WANTT ) THEN * * Apply the transformation to the rest of H. * IF( I2.GT.I ) $ CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH, $ CS, SN ) CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN ) *** * Increment op count OPS = OPS + 6*( I2-I1-1 ) *** END IF IF( WANTZ ) THEN * * Apply the transformation to Z. * CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN ) *** * Increment op count OPS = OPS + 6*NZ *** END IF END IF * * Decrement number of remaining iterations, and return to start of * the main loop with new value of I. * ITN = ITN - ITS I = L - 1 GO TO 10 * 150 CONTINUE *** * Compute final op count OPS = OPS + OPST *** RETURN * * End of DLAHQR * END