001:       SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
002:      $                   IHIZ, Z, LDZ, INFO )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
010:       LOGICAL            WANTT, WANTZ
011: *     ..
012: *     .. Array Arguments ..
013:       COMPLEX            H( LDH, * ), W( * ), Z( LDZ, * )
014: *     ..
015: *
016: *     Purpose
017: *     =======
018: *
019: *     CLAHQR is an auxiliary routine called by CHSEQR to update the
020: *     eigenvalues and Schur decomposition already computed by CHSEQR, by
021: *     dealing with the Hessenberg submatrix in rows and columns ILO to
022: *     IHI.
023: *
024: *     Arguments
025: *     =========
026: *
027: *     WANTT   (input) LOGICAL
028: *          = .TRUE. : the full Schur form T is required;
029: *          = .FALSE.: only eigenvalues are required.
030: *
031: *     WANTZ   (input) LOGICAL
032: *          = .TRUE. : the matrix of Schur vectors Z is required;
033: *          = .FALSE.: Schur vectors are not required.
034: *
035: *     N       (input) INTEGER
036: *          The order of the matrix H.  N >= 0.
037: *
038: *     ILO     (input) INTEGER
039: *     IHI     (input) INTEGER
040: *          It is assumed that H is already upper triangular in rows and
041: *          columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
042: *          CLAHQR works primarily with the Hessenberg submatrix in rows
043: *          and columns ILO to IHI, but applies transformations to all of
044: *          H if WANTT is .TRUE..
045: *          1 <= ILO <= max(1,IHI); IHI <= N.
046: *
047: *     H       (input/output) COMPLEX array, dimension (LDH,N)
048: *          On entry, the upper Hessenberg matrix H.
049: *          On exit, if INFO is zero and if WANTT is .TRUE., then H
050: *          is upper triangular in rows and columns ILO:IHI.  If INFO
051: *          is zero and if WANTT is .FALSE., then the contents of H
052: *          are unspecified on exit.  The output state of H in case
053: *          INF is positive is below under the description of INFO.
054: *
055: *     LDH     (input) INTEGER
056: *          The leading dimension of the array H. LDH >= max(1,N).
057: *
058: *     W       (output) COMPLEX array, dimension (N)
059: *          The computed eigenvalues ILO to IHI are stored in the
060: *          corresponding elements of W. If WANTT is .TRUE., the
061: *          eigenvalues are stored in the same order as on the diagonal
062: *          of the Schur form returned in H, with W(i) = H(i,i).
063: *
064: *     ILOZ    (input) INTEGER
065: *     IHIZ    (input) INTEGER
066: *          Specify the rows of Z to which transformations must be
067: *          applied if WANTZ is .TRUE..
068: *          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
069: *
070: *     Z       (input/output) COMPLEX array, dimension (LDZ,N)
071: *          If WANTZ is .TRUE., on entry Z must contain the current
072: *          matrix Z of transformations accumulated by CHSEQR, and on
073: *          exit Z has been updated; transformations are applied only to
074: *          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
075: *          If WANTZ is .FALSE., Z is not referenced.
076: *
077: *     LDZ     (input) INTEGER
078: *          The leading dimension of the array Z. LDZ >= max(1,N).
079: *
080: *     INFO    (output) INTEGER
081: *           =   0: successful exit
082: *          .GT. 0: if INFO = i, CLAHQR failed to compute all the
083: *                  eigenvalues ILO to IHI in a total of 30 iterations
084: *                  per eigenvalue; elements i+1:ihi of W contain
085: *                  those eigenvalues which have been successfully
086: *                  computed.
087: *
088: *                  If INFO .GT. 0 and WANTT is .FALSE., then on exit,
089: *                  the remaining unconverged eigenvalues are the
090: *                  eigenvalues of the upper Hessenberg matrix
091: *                  rows and columns ILO thorugh INFO of the final,
092: *                  output value of H.
093: *
094: *                  If INFO .GT. 0 and WANTT is .TRUE., then on exit
095: *          (*)       (initial value of H)*U  = U*(final value of H)
096: *                  where U is an orthognal matrix.    The final
097: *                  value of H is upper Hessenberg and triangular in
098: *                  rows and columns INFO+1 through IHI.
099: *
100: *                  If INFO .GT. 0 and WANTZ is .TRUE., then on exit
101: *                      (final value of Z)  = (initial value of Z)*U
102: *                  where U is the orthogonal matrix in (*)
103: *                  (regardless of the value of WANTT.)
104: *
105: *     Further Details
106: *     ===============
107: *
108: *     02-96 Based on modifications by
109: *     David Day, Sandia National Laboratory, USA
110: *
111: *     12-04 Further modifications by
112: *     Ralph Byers, University of Kansas, USA
113: *     This is a modified version of CLAHQR from LAPACK version 3.0.
114: *     It is (1) more robust against overflow and underflow and
115: *     (2) adopts the more conservative Ahues & Tisseur stopping
116: *     criterion (LAWN 122, 1997).
117: *
118: *     =========================================================
119: *
120: *     .. Parameters ..
121:       INTEGER            ITMAX
122:       PARAMETER          ( ITMAX = 30 )
123:       COMPLEX            ZERO, ONE
124:       PARAMETER          ( ZERO = ( 0.0e0, 0.0e0 ),
125:      $                   ONE = ( 1.0e0, 0.0e0 ) )
126:       REAL               RZERO, RONE, HALF
127:       PARAMETER          ( RZERO = 0.0e0, RONE = 1.0e0, HALF = 0.5e0 )
128:       REAL               DAT1
129:       PARAMETER          ( DAT1 = 3.0e0 / 4.0e0 )
130: *     ..
131: *     .. Local Scalars ..
132:       COMPLEX            CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
133:      $                   V2, X, Y
134:       REAL               AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
135:      $                   SAFMIN, SMLNUM, SX, T2, TST, ULP
136:       INTEGER            I, I1, I2, ITS, J, JHI, JLO, K, L, M, NH, NZ
137: *     ..
138: *     .. Local Arrays ..
139:       COMPLEX            V( 2 )
140: *     ..
141: *     .. External Functions ..
142:       COMPLEX            CLADIV
143:       REAL               SLAMCH
144:       EXTERNAL           CLADIV, SLAMCH
145: *     ..
146: *     .. External Subroutines ..
147:       EXTERNAL           CCOPY, CLARFG, CSCAL, SLABAD
148: *     ..
149: *     .. Statement Functions ..
150:       REAL               CABS1
151: *     ..
152: *     .. Intrinsic Functions ..
153:       INTRINSIC          ABS, AIMAG, CONJG, MAX, MIN, REAL, SQRT
154: *     ..
155: *     .. Statement Function definitions ..
156:       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
157: *     ..
158: *     .. Executable Statements ..
159: *
160:       INFO = 0
161: *
162: *     Quick return if possible
163: *
164:       IF( N.EQ.0 )
165:      $   RETURN
166:       IF( ILO.EQ.IHI ) THEN
167:          W( ILO ) = H( ILO, ILO )
168:          RETURN
169:       END IF
170: *
171: *     ==== clear out the trash ====
172:       DO 10 J = ILO, IHI - 3
173:          H( J+2, J ) = ZERO
174:          H( J+3, J ) = ZERO
175:    10 CONTINUE
176:       IF( ILO.LE.IHI-2 )
177:      $   H( IHI, IHI-2 ) = ZERO
178: *     ==== ensure that subdiagonal entries are real ====
179:       IF( WANTT ) THEN
180:          JLO = 1
181:          JHI = N
182:       ELSE
183:          JLO = ILO
184:          JHI = IHI
185:       END IF
186:       DO 20 I = ILO + 1, IHI
187:          IF( AIMAG( H( I, I-1 ) ).NE.RZERO ) THEN
188: *           ==== The following redundant normalization
189: *           .    avoids problems with both gradual and
190: *           .    sudden underflow in ABS(H(I,I-1)) ====
191:             SC = H( I, I-1 ) / CABS1( H( I, I-1 ) )
192:             SC = CONJG( SC ) / ABS( SC )
193:             H( I, I-1 ) = ABS( H( I, I-1 ) )
194:             CALL CSCAL( JHI-I+1, SC, H( I, I ), LDH )
195:             CALL CSCAL( MIN( JHI, I+1 )-JLO+1, CONJG( SC ), H( JLO, I ),
196:      $                  1 )
197:             IF( WANTZ )
198:      $         CALL CSCAL( IHIZ-ILOZ+1, CONJG( SC ), Z( ILOZ, I ), 1 )
199:          END IF
200:    20 CONTINUE
201: *
202:       NH = IHI - ILO + 1
203:       NZ = IHIZ - ILOZ + 1
204: *
205: *     Set machine-dependent constants for the stopping criterion.
206: *
207:       SAFMIN = SLAMCH( 'SAFE MINIMUM' )
208:       SAFMAX = RONE / SAFMIN
209:       CALL SLABAD( SAFMIN, SAFMAX )
210:       ULP = SLAMCH( 'PRECISION' )
211:       SMLNUM = SAFMIN*( REAL( NH ) / ULP )
212: *
213: *     I1 and I2 are the indices of the first row and last column of H
214: *     to which transformations must be applied. If eigenvalues only are
215: *     being computed, I1 and I2 are set inside the main loop.
216: *
217:       IF( WANTT ) THEN
218:          I1 = 1
219:          I2 = N
220:       END IF
221: *
222: *     The main loop begins here. I is the loop index and decreases from
223: *     IHI to ILO in steps of 1. Each iteration of the loop works
224: *     with the active submatrix in rows and columns L to I.
225: *     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
226: *     H(L,L-1) is negligible so that the matrix splits.
227: *
228:       I = IHI
229:    30 CONTINUE
230:       IF( I.LT.ILO )
231:      $   GO TO 150
232: *
233: *     Perform QR iterations on rows and columns ILO to I until a
234: *     submatrix of order 1 splits off at the bottom because a
235: *     subdiagonal element has become negligible.
236: *
237:       L = ILO
238:       DO 130 ITS = 0, ITMAX
239: *
240: *        Look for a single small subdiagonal element.
241: *
242:          DO 40 K = I, L + 1, -1
243:             IF( CABS1( H( K, K-1 ) ).LE.SMLNUM )
244:      $         GO TO 50
245:             TST = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
246:             IF( TST.EQ.ZERO ) THEN
247:                IF( K-2.GE.ILO )
248:      $            TST = TST + ABS( REAL( H( K-1, K-2 ) ) )
249:                IF( K+1.LE.IHI )
250:      $            TST = TST + ABS( REAL( H( K+1, K ) ) )
251:             END IF
252: *           ==== The following is a conservative small subdiagonal
253: *           .    deflation criterion due to Ahues & Tisseur (LAWN 122,
254: *           .    1997). It has better mathematical foundation and
255: *           .    improves accuracy in some examples.  ====
256:             IF( ABS( REAL( H( K, K-1 ) ) ).LE.ULP*TST ) THEN
257:                AB = MAX( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
258:                BA = MIN( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
259:                AA = MAX( CABS1( H( K, K ) ),
260:      $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
261:                BB = MIN( CABS1( H( K, K ) ),
262:      $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
263:                S = AA + AB
264:                IF( BA*( AB / S ).LE.MAX( SMLNUM,
265:      $             ULP*( BB*( AA / S ) ) ) )GO TO 50
266:             END IF
267:    40    CONTINUE
268:    50    CONTINUE
269:          L = K
270:          IF( L.GT.ILO ) THEN
271: *
272: *           H(L,L-1) is negligible
273: *
274:             H( L, L-1 ) = ZERO
275:          END IF
276: *
277: *        Exit from loop if a submatrix of order 1 has split off.
278: *
279:          IF( L.GE.I )
280:      $      GO TO 140
281: *
282: *        Now the active submatrix is in rows and columns L to I. If
283: *        eigenvalues only are being computed, only the active submatrix
284: *        need be transformed.
285: *
286:          IF( .NOT.WANTT ) THEN
287:             I1 = L
288:             I2 = I
289:          END IF
290: *
291:          IF( ITS.EQ.10 ) THEN
292: *
293: *           Exceptional shift.
294: *
295:             S = DAT1*ABS( REAL( H( L+1, L ) ) )
296:             T = S + H( L, L )
297:          ELSE IF( ITS.EQ.20 ) THEN
298: *
299: *           Exceptional shift.
300: *
301:             S = DAT1*ABS( REAL( H( I, I-1 ) ) )
302:             T = S + H( I, I )
303:          ELSE
304: *
305: *           Wilkinson's shift.
306: *
307:             T = H( I, I )
308:             U = SQRT( H( I-1, I ) )*SQRT( H( I, I-1 ) )
309:             S = CABS1( U )
310:             IF( S.NE.RZERO ) THEN
311:                X = HALF*( H( I-1, I-1 )-T )
312:                SX = CABS1( X )
313:                S = MAX( S, CABS1( X ) )
314:                Y = S*SQRT( ( X / S )**2+( U / S )**2 )
315:                IF( SX.GT.RZERO ) THEN
316:                   IF( REAL( X / SX )*REAL( Y )+AIMAG( X / SX )*
317:      $                AIMAG( Y ).LT.RZERO )Y = -Y
318:                END IF
319:                T = T - U*CLADIV( U, ( X+Y ) )
320:             END IF
321:          END IF
322: *
323: *        Look for two consecutive small subdiagonal elements.
324: *
325:          DO 60 M = I - 1, L + 1, -1
326: *
327: *           Determine the effect of starting the single-shift QR
328: *           iteration at row M, and see if this would make H(M,M-1)
329: *           negligible.
330: *
331:             H11 = H( M, M )
332:             H22 = H( M+1, M+1 )
333:             H11S = H11 - T
334:             H21 = REAL( H( M+1, M ) )
335:             S = CABS1( H11S ) + ABS( H21 )
336:             H11S = H11S / S
337:             H21 = H21 / S
338:             V( 1 ) = H11S
339:             V( 2 ) = H21
340:             H10 = REAL( H( M, M-1 ) )
341:             IF( ABS( H10 )*ABS( H21 ).LE.ULP*
342:      $          ( CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) ) ) )
343:      $          GO TO 70
344:    60    CONTINUE
345:          H11 = H( L, L )
346:          H22 = H( L+1, L+1 )
347:          H11S = H11 - T
348:          H21 = REAL( H( L+1, L ) )
349:          S = CABS1( H11S ) + ABS( H21 )
350:          H11S = H11S / S
351:          H21 = H21 / S
352:          V( 1 ) = H11S
353:          V( 2 ) = H21
354:    70    CONTINUE
355: *
356: *        Single-shift QR step
357: *
358:          DO 120 K = M, I - 1
359: *
360: *           The first iteration of this loop determines a reflection G
361: *           from the vector V and applies it from left and right to H,
362: *           thus creating a nonzero bulge below the subdiagonal.
363: *
364: *           Each subsequent iteration determines a reflection G to
365: *           restore the Hessenberg form in the (K-1)th column, and thus
366: *           chases the bulge one step toward the bottom of the active
367: *           submatrix.
368: *
369: *           V(2) is always real before the call to CLARFG, and hence
370: *           after the call T2 ( = T1*V(2) ) is also real.
371: *
372:             IF( K.GT.M )
373:      $         CALL CCOPY( 2, H( K, K-1 ), 1, V, 1 )
374:             CALL CLARFG( 2, V( 1 ), V( 2 ), 1, T1 )
375:             IF( K.GT.M ) THEN
376:                H( K, K-1 ) = V( 1 )
377:                H( K+1, K-1 ) = ZERO
378:             END IF
379:             V2 = V( 2 )
380:             T2 = REAL( T1*V2 )
381: *
382: *           Apply G from the left to transform the rows of the matrix
383: *           in columns K to I2.
384: *
385:             DO 80 J = K, I2
386:                SUM = CONJG( T1 )*H( K, J ) + T2*H( K+1, J )
387:                H( K, J ) = H( K, J ) - SUM
388:                H( K+1, J ) = H( K+1, J ) - SUM*V2
389:    80       CONTINUE
390: *
391: *           Apply G from the right to transform the columns of the
392: *           matrix in rows I1 to min(K+2,I).
393: *
394:             DO 90 J = I1, MIN( K+2, I )
395:                SUM = T1*H( J, K ) + T2*H( J, K+1 )
396:                H( J, K ) = H( J, K ) - SUM
397:                H( J, K+1 ) = H( J, K+1 ) - SUM*CONJG( V2 )
398:    90       CONTINUE
399: *
400:             IF( WANTZ ) THEN
401: *
402: *              Accumulate transformations in the matrix Z
403: *
404:                DO 100 J = ILOZ, IHIZ
405:                   SUM = T1*Z( J, K ) + T2*Z( J, K+1 )
406:                   Z( J, K ) = Z( J, K ) - SUM
407:                   Z( J, K+1 ) = Z( J, K+1 ) - SUM*CONJG( V2 )
408:   100          CONTINUE
409:             END IF
410: *
411:             IF( K.EQ.M .AND. M.GT.L ) THEN
412: *
413: *              If the QR step was started at row M > L because two
414: *              consecutive small subdiagonals were found, then extra
415: *              scaling must be performed to ensure that H(M,M-1) remains
416: *              real.
417: *
418:                TEMP = ONE - T1
419:                TEMP = TEMP / ABS( TEMP )
420:                H( M+1, M ) = H( M+1, M )*CONJG( TEMP )
421:                IF( M+2.LE.I )
422:      $            H( M+2, M+1 ) = H( M+2, M+1 )*TEMP
423:                DO 110 J = M, I
424:                   IF( J.NE.M+1 ) THEN
425:                      IF( I2.GT.J )
426:      $                  CALL CSCAL( I2-J, TEMP, H( J, J+1 ), LDH )
427:                      CALL CSCAL( J-I1, CONJG( TEMP ), H( I1, J ), 1 )
428:                      IF( WANTZ ) THEN
429:                         CALL CSCAL( NZ, CONJG( TEMP ), Z( ILOZ, J ), 1 )
430:                      END IF
431:                   END IF
432:   110          CONTINUE
433:             END IF
434:   120    CONTINUE
435: *
436: *        Ensure that H(I,I-1) is real.
437: *
438:          TEMP = H( I, I-1 )
439:          IF( AIMAG( TEMP ).NE.RZERO ) THEN
440:             RTEMP = ABS( TEMP )
441:             H( I, I-1 ) = RTEMP
442:             TEMP = TEMP / RTEMP
443:             IF( I2.GT.I )
444:      $         CALL CSCAL( I2-I, CONJG( TEMP ), H( I, I+1 ), LDH )
445:             CALL CSCAL( I-I1, TEMP, H( I1, I ), 1 )
446:             IF( WANTZ ) THEN
447:                CALL CSCAL( NZ, TEMP, Z( ILOZ, I ), 1 )
448:             END IF
449:          END IF
450: *
451:   130 CONTINUE
452: *
453: *     Failure to converge in remaining number of iterations
454: *
455:       INFO = I
456:       RETURN
457: *
458:   140 CONTINUE
459: *
460: *     H(I,I-1) is negligible: one eigenvalue has converged.
461: *
462:       W( I ) = H( I, I )
463: *
464: *     return to start of the main loop with new value of I.
465: *
466:       I = L - 1
467:       GO TO 30
468: *
469:   150 CONTINUE
470:       RETURN
471: *
472: *     End of CLAHQR
473: *
474:       END
475: