SUBROUTINE ZHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, $ WORK, LWORK, INFO ) * * -- LAPACK routine (instrumented to count operations, 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 .. CHARACTER COMPZ, JOB INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N * .. * .. Array Arguments .. COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * ) * .. * Common block to return operation count. * .. Common blocks .. COMMON / LATIME / OPS, ITCNT * .. * .. Scalars in Common .. DOUBLE PRECISION ITCNT, OPS * .. * * Purpose * ======= * * ZHSEQR computes the eigenvalues of a complex upper Hessenberg * matrix H, and, optionally, the matrices T and Z from the Schur * decomposition H = Z T Z**H, where T is an upper triangular matrix * (the Schur form), and Z is the unitary matrix of Schur vectors. * * Optionally Z may be postmultiplied into an input unitary matrix Q, * so that this routine can give the Schur factorization of a matrix A * which has been reduced to the Hessenberg form H by the unitary * matrix Q: A = Q*H*Q**H = (QZ)*T*(QZ)**H. * * Arguments * ========= * * JOB (input) CHARACTER*1 * = 'E': compute eigenvalues only; * = 'S': compute eigenvalues and the Schur form T. * * COMPZ (input) CHARACTER*1 * = 'N': no Schur vectors are computed; * = 'I': Z is initialized to the unit matrix and the matrix Z * of Schur vectors of H is returned; * = 'V': Z must contain an unitary matrix Q on entry, and * the product Q*Z is returned. * * 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 triangular in rows * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally * set by a previous call to ZGEBAL, and then passed to CGEHRD * when the matrix output by ZGEBAL is reduced to Hessenberg * form. Otherwise ILO and IHI should be set to 1 and N * respectively. * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. * * H (input/output) COMPLEX*16 array, dimension (LDH,N) * On entry, the upper Hessenberg matrix H. * On exit, if JOB = 'S', H contains the upper triangular matrix * T from the Schur decomposition (the Schur form). If * JOB = 'E', the contents of H are unspecified on exit. * * LDH (input) INTEGER * The leading dimension of the array H. LDH >= max(1,N). * * W (output) COMPLEX*16 array, dimension (N) * The computed eigenvalues. If JOB = 'S', the eigenvalues are * stored in the same order as on the diagonal of the Schur form * returned in H, with W(i) = H(i,i). * * Z (input/output) COMPLEX*16 array, dimension (LDZ,N) * If COMPZ = 'N': Z is not referenced. * If COMPZ = 'I': on entry, Z need not be set, and on exit, Z * contains the unitary matrix Z of the Schur vectors of H. * If COMPZ = 'V': on entry Z must contain an N-by-N matrix Q, * which is assumed to be equal to the unit matrix except for * the submatrix Z(ILO:IHI,ILO:IHI); on exit Z contains Q*Z. * Normally Q is the unitary matrix generated by ZUNGHR after * the call to ZGEHRD which formed the Hessenberg matrix H. * * LDZ (input) INTEGER * The leading dimension of the array Z. * LDZ >= max(1,N) if COMPZ = 'I' or 'V'; LDZ >= 1 otherwise. * * WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= max(1,N). * * If LWORK = -1, then a workspace query is assumed; the routine * only calculates the optimal size of the WORK array, returns * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, ZHSEQR failed to compute all the * eigenvalues in a total of 30*(IHI-ILO+1) iterations; * elements 1:ilo-1 and i+1:n of W contain those * eigenvalues which have been successfully computed. * * ===================================================================== * * .. Parameters .. COMPLEX*16 ZERO, ONE PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), $ ONE = ( 1.0D+0, 0.0D+0 ) ) DOUBLE PRECISION RZERO, RONE, CONST PARAMETER ( RZERO = 0.0D+0, RONE = 1.0D+0, $ CONST = 1.5D+0 ) INTEGER NSMAX, LDS PARAMETER ( NSMAX = 15, LDS = NSMAX ) * .. * .. Local Scalars .. LOGICAL INITZ, LQUERY, WANTT, WANTZ INTEGER I, I1, I2, IERR, II, ITEMP, ITN, ITS, J, K, L, $ MAXB, NH, NR, NS, NV DOUBLE PRECISION OPST, OVFL, RTEMP, SMLNUM, TST1, ULP, UNFL COMPLEX*16 CDUM, TAU, TEMP * .. * .. Local Arrays .. DOUBLE PRECISION RWORK( 1 ) COMPLEX*16 S( LDS, NSMAX ), V( NSMAX+1 ), VV( NSMAX+1 ) * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV, IZAMAX DOUBLE PRECISION DLAMCH, DLAPY2, ZLANHS EXTERNAL LSAME, ILAENV, IZAMAX, DLAMCH, DLAPY2, ZLANHS * .. * .. External Subroutines .. EXTERNAL DLABAD, XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLACPY, $ ZLAHQR, ZLARFG, ZLARFX, ZLASET, ZSCAL * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN * .. * .. Statement Functions .. DOUBLE PRECISION CABS1 * .. * .. Statement Function definitions .. CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) * .. * .. Executable Statements .. * * Decode and test the input parameters * WANTT = LSAME( JOB, 'S' ) INITZ = LSAME( COMPZ, 'I' ) WANTZ = INITZ .OR. LSAME( COMPZ, 'V' ) * INFO = 0 WORK( 1 ) = MAX( 1, N ) LQUERY = ( LWORK.EQ.-1 ) IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN INFO = -1 ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN INFO = -4 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN INFO = -5 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN INFO = -7 ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. LDZ.LT.MAX( 1, N ) ) THEN INFO = -10 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN INFO = -12 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'ZHSEQR', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF *** * Initialize OPST = 0 *** * * Initialize Z, if necessary * IF( INITZ ) $ CALL ZLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) * * Store the eigenvalues isolated by ZGEBAL. * DO 10 I = 1, ILO - 1 W( I ) = H( I, I ) 10 CONTINUE DO 20 I = IHI + 1, N W( I ) = H( I, I ) 20 CONTINUE * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN IF( ILO.EQ.IHI ) THEN W( ILO ) = H( ILO, ILO ) RETURN END IF * * Set rows and columns ILO to IHI to zero below the first * subdiagonal. * DO 40 J = ILO, IHI - 2 DO 30 I = J + 2, N H( I, J ) = ZERO 30 CONTINUE 40 CONTINUE NH = IHI - ILO + 1 * * 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 re-set inside the main loop. * IF( WANTT ) THEN I1 = 1 I2 = N ELSE I1 = ILO I2 = IHI END IF * * Ensure that the subdiagonal elements are real. * DO 50 I = ILO + 1, IHI TEMP = H( I, I-1 ) IF( DIMAG( TEMP ).NE.RZERO ) THEN RTEMP = DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) ) H( I, I-1 ) = RTEMP TEMP = TEMP / RTEMP IF( I2.GT.I ) $ CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH ) CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 ) IF( I.LT.IHI ) $ H( I+1, I ) = TEMP*H( I+1, I ) *** * Increment op count OPST = OPST + 6*( I2-I1+2 ) *** IF( WANTZ ) THEN CALL ZSCAL( NH, TEMP, Z( ILO, I ), 1 ) *** * Increment op count OPST = OPST + 6*NH *** END IF END IF 50 CONTINUE * * Determine the order of the multi-shift QR algorithm to be used. * NS = ILAENV( 4, 'ZHSEQR', JOB // COMPZ, N, ILO, IHI, -1 ) MAXB = ILAENV( 8, 'ZHSEQR', JOB // COMPZ, N, ILO, IHI, -1 ) IF( NS.LE.1 .OR. NS.GT.NH .OR. MAXB.GE.NH ) THEN * * Use the standard double-shift algorithm * CALL ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILO, IHI, Z, $ LDZ, INFO ) RETURN END IF MAXB = MAX( 2, MAXB ) NS = MIN( NS, MAXB, NSMAX ) * * Now 1 < NS <= MAXB < NH. * * Set machine-dependent constants for the stopping criterion. * If norm(H) <= sqrt(OVFL), overflow should not occur. * UNFL = DLAMCH( 'Safe minimum' ) OVFL = RONE / UNFL CALL DLABAD( UNFL, OVFL ) ULP = DLAMCH( 'Precision' ) SMLNUM = UNFL*( NH / ULP ) * * ITN is the total number of multiple-shift 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 at most MAXB. 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 60 CONTINUE IF( I.LT.ILO ) $ GO TO 180 * * Perform multiple-shift QR iterations on rows and columns ILO to I * until a submatrix of order at most MAXB splits off at the bottom * because a subdiagonal element has become negligible. * L = ILO DO 160 ITS = 0, ITN * * Look for a single small subdiagonal element. * DO 70 K = I, L + 1, -1 TST1 = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) ) IF( TST1.EQ.RZERO ) THEN TST1 = ZLANHS( '1', I-L+1, H( L, L ), LDH, RWORK ) *** * Increment op count OPS = OPS + 5*( I-L+1 )*( I-L ) / 2 *** END IF IF( ABS( DBLE( H( K, K-1 ) ) ).LE.MAX( ULP*TST1, SMLNUM ) ) $ GO TO 80 70 CONTINUE 80 CONTINUE L = K *** * Increment op count OPST = OPST + 5*( 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 <= MAXB has split off. * IF( L.GE.I-MAXB+1 ) $ GO TO 170 * * 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.20 .OR. ITS.EQ.30 ) THEN * * Exceptional shifts. * DO 90 II = I - NS + 1, I W( II ) = CONST*( ABS( DBLE( H( II, II-1 ) ) )+ $ ABS( DBLE( H( II, II ) ) ) ) 90 CONTINUE *** * Increment op count OPST = OPST + 2*NS *** ELSE * * Use eigenvalues of trailing submatrix of order NS as shifts. * CALL ZLACPY( 'Full', NS, NS, H( I-NS+1, I-NS+1 ), LDH, S, $ LDS ) CALL ZLAHQR( .FALSE., .FALSE., NS, 1, NS, S, LDS, $ W( I-NS+1 ), 1, NS, Z, LDZ, IERR ) IF( IERR.GT.0 ) THEN * * If ZLAHQR failed to compute all NS eigenvalues, use the * unconverged diagonal elements as the remaining shifts. * DO 100 II = 1, IERR W( I-NS+II ) = S( II, II ) 100 CONTINUE END IF END IF * * Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns)) * where G is the Hessenberg submatrix H(L:I,L:I) and w is * the vector of shifts (stored in W). The result is * stored in the local array V. * V( 1 ) = ONE DO 110 II = 2, NS + 1 V( II ) = ZERO 110 CONTINUE NV = 1 DO 130 J = I - NS + 1, I CALL ZCOPY( NV+1, V, 1, VV, 1 ) CALL ZGEMV( 'No transpose', NV+1, NV, ONE, H( L, L ), LDH, $ VV, 1, -W( J ), V, 1 ) NV = NV + 1 *** * Increment op count OPST = OPST + 8*NV*( N+1 ) + 6*( NV+1 ) *** * * Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero, * reset it to the unit vector. * ITEMP = IZAMAX( NV, V, 1 ) *** * Increment op count OPST = OPST + 2*NV *** RTEMP = CABS1( V( ITEMP ) ) IF( RTEMP.EQ.RZERO ) THEN V( 1 ) = ONE DO 120 II = 2, NV V( II ) = ZERO 120 CONTINUE ELSE RTEMP = MAX( RTEMP, SMLNUM ) CALL ZDSCAL( NV, RONE / RTEMP, V, 1 ) *** * Increment op count OPST = OPST + 2*NV *** END IF 130 CONTINUE * * Multiple-shift QR step * DO 150 K = L, 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( NS+1, I-K+1 ) IF( K.GT.L ) $ CALL ZCOPY( NR, H( K, K-1 ), 1, V, 1 ) CALL ZLARFG( NR, V( 1 ), V( 2 ), 1, TAU ) *** * Increment op count OPST = OPST + 10*NR + 12 *** IF( K.GT.L ) THEN H( K, K-1 ) = V( 1 ) DO 140 II = K + 1, I H( II, K-1 ) = ZERO 140 CONTINUE END IF V( 1 ) = ONE * * Apply G' from the left to transform the rows of the matrix * in columns K to I2. * CALL ZLARFX( 'Left', NR, I2-K+1, V, DCONJG( TAU ), $ H( K, K ), LDH, WORK ) * * Apply G from the right to transform the columns of the * matrix in rows I1 to min(K+NR,I). * CALL ZLARFX( 'Right', MIN( K+NR, I )-I1+1, NR, V, TAU, $ H( I1, K ), LDH, WORK ) *** * Increment op count OPS = OPS + 4*( 4*NR-2 )*( I2-I1+2+MIN( NR, I-K ) ) *** * IF( WANTZ ) THEN * * Accumulate transformations in the matrix Z * CALL ZLARFX( 'Right', NH, NR, V, TAU, Z( ILO, K ), LDZ, $ WORK ) *** * Increment op count OPS = OPS + 4*( 4*NR-2 )*NH *** END IF 150 CONTINUE * * Ensure that H(I,I-1) is real. * TEMP = H( I, I-1 ) IF( DIMAG( TEMP ).NE.RZERO ) THEN RTEMP = DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) ) H( I, I-1 ) = RTEMP TEMP = TEMP / RTEMP IF( I2.GT.I ) $ CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH ) CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 ) *** * Increment op count OPST = OPST + 6*( I2-I1+1 ) *** IF( WANTZ ) THEN CALL ZSCAL( NH, TEMP, Z( ILO, I ), 1 ) *** * Increment op count OPST = OPST + 6*NH *** END IF END IF * 160 CONTINUE * * Failure to converge in remaining number of iterations * INFO = I RETURN * 170 CONTINUE * * A submatrix of order <= MAXB in rows and columns L to I has split * off. Use the double-shift QR algorithm to handle it. * CALL ZLAHQR( WANTT, WANTZ, N, L, I, H, LDH, W, ILO, IHI, Z, LDZ, $ INFO ) IF( INFO.GT.0 ) $ RETURN * * Decrement number of remaining iterations, and return to start of * the main loop with a new value of I. * ITN = ITN - ITS I = L - 1 GO TO 60 * 180 CONTINUE *** * Compute final op count OPS = OPS + OPST *** WORK( 1 ) = MAX( 1, N ) RETURN * * End of ZHSEQR * END