001:       SUBROUTINE ZLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S,
002:      $                   H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV,
003:      $                   WV, LDWV, NH, WH, LDWH )
004: *
005: *  -- LAPACK auxiliary routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
011:      $                   LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
012:       LOGICAL            WANTT, WANTZ
013: *     ..
014: *     .. Array Arguments ..
015:       COMPLEX*16         H( LDH, * ), S( * ), U( LDU, * ), V( LDV, * ),
016:      $                   WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
017: *     ..
018: *
019: *     This auxiliary subroutine called by ZLAQR0 performs a
020: *     single small-bulge multi-shift QR sweep.
021: *
022: *      WANTT  (input) logical scalar
023: *             WANTT = .true. if the triangular Schur factor
024: *             is being computed.  WANTT is set to .false. otherwise.
025: *
026: *      WANTZ  (input) logical scalar
027: *             WANTZ = .true. if the unitary Schur factor is being
028: *             computed.  WANTZ is set to .false. otherwise.
029: *
030: *      KACC22 (input) integer with value 0, 1, or 2.
031: *             Specifies the computation mode of far-from-diagonal
032: *             orthogonal updates.
033: *        = 0: ZLAQR5 does not accumulate reflections and does not
034: *             use matrix-matrix multiply to update far-from-diagonal
035: *             matrix entries.
036: *        = 1: ZLAQR5 accumulates reflections and uses matrix-matrix
037: *             multiply to update the far-from-diagonal matrix entries.
038: *        = 2: ZLAQR5 accumulates reflections, uses matrix-matrix
039: *             multiply to update the far-from-diagonal matrix entries,
040: *             and takes advantage of 2-by-2 block structure during
041: *             matrix multiplies.
042: *
043: *      N      (input) integer scalar
044: *             N is the order of the Hessenberg matrix H upon which this
045: *             subroutine operates.
046: *
047: *      KTOP   (input) integer scalar
048: *      KBOT   (input) integer scalar
049: *             These are the first and last rows and columns of an
050: *             isolated diagonal block upon which the QR sweep is to be
051: *             applied. It is assumed without a check that
052: *                       either KTOP = 1  or   H(KTOP,KTOP-1) = 0
053: *             and
054: *                       either KBOT = N  or   H(KBOT+1,KBOT) = 0.
055: *
056: *      NSHFTS (input) integer scalar
057: *             NSHFTS gives the number of simultaneous shifts.  NSHFTS
058: *             must be positive and even.
059: *
060: *      S      (input/output) COMPLEX*16 array of size (NSHFTS)
061: *             S contains the shifts of origin that define the multi-
062: *             shift QR sweep.  On output S may be reordered.
063: *
064: *      H      (input/output) COMPLEX*16 array of size (LDH,N)
065: *             On input H contains a Hessenberg matrix.  On output a
066: *             multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
067: *             to the isolated diagonal block in rows and columns KTOP
068: *             through KBOT.
069: *
070: *      LDH    (input) integer scalar
071: *             LDH is the leading dimension of H just as declared in the
072: *             calling procedure.  LDH.GE.MAX(1,N).
073: *
074: *      ILOZ   (input) INTEGER
075: *      IHIZ   (input) INTEGER
076: *             Specify the rows of Z to which transformations must be
077: *             applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
078: *
079: *      Z      (input/output) COMPLEX*16 array of size (LDZ,IHI)
080: *             If WANTZ = .TRUE., then the QR Sweep unitary
081: *             similarity transformation is accumulated into
082: *             Z(ILOZ:IHIZ,ILO:IHI) from the right.
083: *             If WANTZ = .FALSE., then Z is unreferenced.
084: *
085: *      LDZ    (input) integer scalar
086: *             LDA is the leading dimension of Z just as declared in
087: *             the calling procedure. LDZ.GE.N.
088: *
089: *      V      (workspace) COMPLEX*16 array of size (LDV,NSHFTS/2)
090: *
091: *      LDV    (input) integer scalar
092: *             LDV is the leading dimension of V as declared in the
093: *             calling procedure.  LDV.GE.3.
094: *
095: *      U      (workspace) COMPLEX*16 array of size
096: *             (LDU,3*NSHFTS-3)
097: *
098: *      LDU    (input) integer scalar
099: *             LDU is the leading dimension of U just as declared in the
100: *             in the calling subroutine.  LDU.GE.3*NSHFTS-3.
101: *
102: *      NH     (input) integer scalar
103: *             NH is the number of columns in array WH available for
104: *             workspace. NH.GE.1.
105: *
106: *      WH     (workspace) COMPLEX*16 array of size (LDWH,NH)
107: *
108: *      LDWH   (input) integer scalar
109: *             Leading dimension of WH just as declared in the
110: *             calling procedure.  LDWH.GE.3*NSHFTS-3.
111: *
112: *      NV     (input) integer scalar
113: *             NV is the number of rows in WV agailable for workspace.
114: *             NV.GE.1.
115: *
116: *      WV     (workspace) COMPLEX*16 array of size
117: *             (LDWV,3*NSHFTS-3)
118: *
119: *      LDWV   (input) integer scalar
120: *             LDWV is the leading dimension of WV as declared in the
121: *             in the calling subroutine.  LDWV.GE.NV.
122: *
123: *     ================================================================
124: *     Based on contributions by
125: *        Karen Braman and Ralph Byers, Department of Mathematics,
126: *        University of Kansas, USA
127: *
128: *     ================================================================
129: *     Reference:
130: *
131: *     K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
132: *     Algorithm Part I: Maintaining Well Focused Shifts, and
133: *     Level 3 Performance, SIAM Journal of Matrix Analysis,
134: *     volume 23, pages 929--947, 2002.
135: *
136: *     ================================================================
137: *     .. Parameters ..
138:       COMPLEX*16         ZERO, ONE
139:       PARAMETER          ( ZERO = ( 0.0d0, 0.0d0 ),
140:      $                   ONE = ( 1.0d0, 0.0d0 ) )
141:       DOUBLE PRECISION   RZERO, RONE
142:       PARAMETER          ( RZERO = 0.0d0, RONE = 1.0d0 )
143: *     ..
144: *     .. Local Scalars ..
145:       COMPLEX*16         ALPHA, BETA, CDUM, REFSUM
146:       DOUBLE PRECISION   H11, H12, H21, H22, SAFMAX, SAFMIN, SCL,
147:      $                   SMLNUM, TST1, TST2, ULP
148:       INTEGER            I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN,
149:      $                   JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
150:      $                   M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL,
151:      $                   NS, NU
152:       LOGICAL            ACCUM, BLK22, BMP22
153: *     ..
154: *     .. External Functions ..
155:       DOUBLE PRECISION   DLAMCH
156:       EXTERNAL           DLAMCH
157: *     ..
158: *     .. Intrinsic Functions ..
159: *
160:       INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN, MOD
161: *     ..
162: *     .. Local Arrays ..
163:       COMPLEX*16         VT( 3 )
164: *     ..
165: *     .. External Subroutines ..
166:       EXTERNAL           DLABAD, ZGEMM, ZLACPY, ZLAQR1, ZLARFG, ZLASET,
167:      $                   ZTRMM
168: *     ..
169: *     .. Statement Functions ..
170:       DOUBLE PRECISION   CABS1
171: *     ..
172: *     .. Statement Function definitions ..
173:       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
174: *     ..
175: *     .. Executable Statements ..
176: *
177: *     ==== If there are no shifts, then there is nothing to do. ====
178: *
179:       IF( NSHFTS.LT.2 )
180:      $   RETURN
181: *
182: *     ==== If the active block is empty or 1-by-1, then there
183: *     .    is nothing to do. ====
184: *
185:       IF( KTOP.GE.KBOT )
186:      $   RETURN
187: *
188: *     ==== NSHFTS is supposed to be even, but if it is odd,
189: *     .    then simply reduce it by one.  ====
190: *
191:       NS = NSHFTS - MOD( NSHFTS, 2 )
192: *
193: *     ==== Machine constants for deflation ====
194: *
195:       SAFMIN = DLAMCH( 'SAFE MINIMUM' )
196:       SAFMAX = RONE / SAFMIN
197:       CALL DLABAD( SAFMIN, SAFMAX )
198:       ULP = DLAMCH( 'PRECISION' )
199:       SMLNUM = SAFMIN*( DBLE( N ) / ULP )
200: *
201: *     ==== Use accumulated reflections to update far-from-diagonal
202: *     .    entries ? ====
203: *
204:       ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )
205: *
206: *     ==== If so, exploit the 2-by-2 block structure? ====
207: *
208:       BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 )
209: *
210: *     ==== clear trash ====
211: *
212:       IF( KTOP+2.LE.KBOT )
213:      $   H( KTOP+2, KTOP ) = ZERO
214: *
215: *     ==== NBMPS = number of 2-shift bulges in the chain ====
216: *
217:       NBMPS = NS / 2
218: *
219: *     ==== KDU = width of slab ====
220: *
221:       KDU = 6*NBMPS - 3
222: *
223: *     ==== Create and chase chains of NBMPS bulges ====
224: *
225:       DO 210 INCOL = 3*( 1-NBMPS ) + KTOP - 1, KBOT - 2, 3*NBMPS - 2
226:          NDCOL = INCOL + KDU
227:          IF( ACCUM )
228:      $      CALL ZLASET( 'ALL', KDU, KDU, ZERO, ONE, U, LDU )
229: *
230: *        ==== Near-the-diagonal bulge chase.  The following loop
231: *        .    performs the near-the-diagonal part of a small bulge
232: *        .    multi-shift QR sweep.  Each 6*NBMPS-2 column diagonal
233: *        .    chunk extends from column INCOL to column NDCOL
234: *        .    (including both column INCOL and column NDCOL). The
235: *        .    following loop chases a 3*NBMPS column long chain of
236: *        .    NBMPS bulges 3*NBMPS-2 columns to the right.  (INCOL
237: *        .    may be less than KTOP and and NDCOL may be greater than
238: *        .    KBOT indicating phantom columns from which to chase
239: *        .    bulges before they are actually introduced or to which
240: *        .    to chase bulges beyond column KBOT.)  ====
241: *
242:          DO 140 KRCOL = INCOL, MIN( INCOL+3*NBMPS-3, KBOT-2 )
243: *
244: *           ==== Bulges number MTOP to MBOT are active double implicit
245: *           .    shift bulges.  There may or may not also be small
246: *           .    2-by-2 bulge, if there is room.  The inactive bulges
247: *           .    (if any) must wait until the active bulges have moved
248: *           .    down the diagonal to make room.  The phantom matrix
249: *           .    paradigm described above helps keep track.  ====
250: *
251:             MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 )
252:             MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 )
253:             M22 = MBOT + 1
254:             BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+3*( M22-1 ) ).EQ.
255:      $              ( KBOT-2 )
256: *
257: *           ==== Generate reflections to chase the chain right
258: *           .    one column.  (The minimum value of K is KTOP-1.) ====
259: *
260:             DO 10 M = MTOP, MBOT
261:                K = KRCOL + 3*( M-1 )
262:                IF( K.EQ.KTOP-1 ) THEN
263:                   CALL ZLAQR1( 3, H( KTOP, KTOP ), LDH, S( 2*M-1 ),
264:      $                         S( 2*M ), V( 1, M ) )
265:                   ALPHA = V( 1, M )
266:                   CALL ZLARFG( 3, ALPHA, V( 2, M ), 1, V( 1, M ) )
267:                ELSE
268:                   BETA = H( K+1, K )
269:                   V( 2, M ) = H( K+2, K )
270:                   V( 3, M ) = H( K+3, K )
271:                   CALL ZLARFG( 3, BETA, V( 2, M ), 1, V( 1, M ) )
272: *
273: *                 ==== A Bulge may collapse because of vigilant
274: *                 .    deflation or destructive underflow.  In the
275: *                 .    underflow case, try the two-small-subdiagonals
276: *                 .    trick to try to reinflate the bulge.  ====
277: *
278:                   IF( H( K+3, K ).NE.ZERO .OR. H( K+3, K+1 ).NE.
279:      $                ZERO .OR. H( K+3, K+2 ).EQ.ZERO ) THEN
280: *
281: *                    ==== Typical case: not collapsed (yet). ====
282: *
283:                      H( K+1, K ) = BETA
284:                      H( K+2, K ) = ZERO
285:                      H( K+3, K ) = ZERO
286:                   ELSE
287: *
288: *                    ==== Atypical case: collapsed.  Attempt to
289: *                    .    reintroduce ignoring H(K+1,K) and H(K+2,K).
290: *                    .    If the fill resulting from the new
291: *                    .    reflector is too large, then abandon it.
292: *                    .    Otherwise, use the new one. ====
293: *
294:                      CALL ZLAQR1( 3, H( K+1, K+1 ), LDH, S( 2*M-1 ),
295:      $                            S( 2*M ), VT )
296:                      ALPHA = VT( 1 )
297:                      CALL ZLARFG( 3, ALPHA, VT( 2 ), 1, VT( 1 ) )
298:                      REFSUM = DCONJG( VT( 1 ) )*
299:      $                        ( H( K+1, K )+DCONJG( VT( 2 ) )*
300:      $                        H( K+2, K ) )
301: *
302:                      IF( CABS1( H( K+2, K )-REFSUM*VT( 2 ) )+
303:      $                   CABS1( REFSUM*VT( 3 ) ).GT.ULP*
304:      $                   ( CABS1( H( K, K ) )+CABS1( H( K+1,
305:      $                   K+1 ) )+CABS1( H( K+2, K+2 ) ) ) ) THEN
306: *
307: *                       ==== Starting a new bulge here would
308: *                       .    create non-negligible fill.  Use
309: *                       .    the old one with trepidation. ====
310: *
311:                         H( K+1, K ) = BETA
312:                         H( K+2, K ) = ZERO
313:                         H( K+3, K ) = ZERO
314:                      ELSE
315: *
316: *                       ==== Stating a new bulge here would
317: *                       .    create only negligible fill.
318: *                       .    Replace the old reflector with
319: *                       .    the new one. ====
320: *
321:                         H( K+1, K ) = H( K+1, K ) - REFSUM
322:                         H( K+2, K ) = ZERO
323:                         H( K+3, K ) = ZERO
324:                         V( 1, M ) = VT( 1 )
325:                         V( 2, M ) = VT( 2 )
326:                         V( 3, M ) = VT( 3 )
327:                      END IF
328:                   END IF
329:                END IF
330:    10       CONTINUE
331: *
332: *           ==== Generate a 2-by-2 reflection, if needed. ====
333: *
334:             K = KRCOL + 3*( M22-1 )
335:             IF( BMP22 ) THEN
336:                IF( K.EQ.KTOP-1 ) THEN
337:                   CALL ZLAQR1( 2, H( K+1, K+1 ), LDH, S( 2*M22-1 ),
338:      $                         S( 2*M22 ), V( 1, M22 ) )
339:                   BETA = V( 1, M22 )
340:                   CALL ZLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
341:                ELSE
342:                   BETA = H( K+1, K )
343:                   V( 2, M22 ) = H( K+2, K )
344:                   CALL ZLARFG( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
345:                   H( K+1, K ) = BETA
346:                   H( K+2, K ) = ZERO
347:                END IF
348:             END IF
349: *
350: *           ==== Multiply H by reflections from the left ====
351: *
352:             IF( ACCUM ) THEN
353:                JBOT = MIN( NDCOL, KBOT )
354:             ELSE IF( WANTT ) THEN
355:                JBOT = N
356:             ELSE
357:                JBOT = KBOT
358:             END IF
359:             DO 30 J = MAX( KTOP, KRCOL ), JBOT
360:                MEND = MIN( MBOT, ( J-KRCOL+2 ) / 3 )
361:                DO 20 M = MTOP, MEND
362:                   K = KRCOL + 3*( M-1 )
363:                   REFSUM = DCONJG( V( 1, M ) )*
364:      $                     ( H( K+1, J )+DCONJG( V( 2, M ) )*
365:      $                     H( K+2, J )+DCONJG( V( 3, M ) )*H( K+3, J ) )
366:                   H( K+1, J ) = H( K+1, J ) - REFSUM
367:                   H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M )
368:                   H( K+3, J ) = H( K+3, J ) - REFSUM*V( 3, M )
369:    20          CONTINUE
370:    30       CONTINUE
371:             IF( BMP22 ) THEN
372:                K = KRCOL + 3*( M22-1 )
373:                DO 40 J = MAX( K+1, KTOP ), JBOT
374:                   REFSUM = DCONJG( V( 1, M22 ) )*
375:      $                     ( H( K+1, J )+DCONJG( V( 2, M22 ) )*
376:      $                     H( K+2, J ) )
377:                   H( K+1, J ) = H( K+1, J ) - REFSUM
378:                   H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 )
379:    40          CONTINUE
380:             END IF
381: *
382: *           ==== Multiply H by reflections from the right.
383: *           .    Delay filling in the last row until the
384: *           .    vigilant deflation check is complete. ====
385: *
386:             IF( ACCUM ) THEN
387:                JTOP = MAX( KTOP, INCOL )
388:             ELSE IF( WANTT ) THEN
389:                JTOP = 1
390:             ELSE
391:                JTOP = KTOP
392:             END IF
393:             DO 80 M = MTOP, MBOT
394:                IF( V( 1, M ).NE.ZERO ) THEN
395:                   K = KRCOL + 3*( M-1 )
396:                   DO 50 J = JTOP, MIN( KBOT, K+3 )
397:                      REFSUM = V( 1, M )*( H( J, K+1 )+V( 2, M )*
398:      $                        H( J, K+2 )+V( 3, M )*H( J, K+3 ) )
399:                      H( J, K+1 ) = H( J, K+1 ) - REFSUM
400:                      H( J, K+2 ) = H( J, K+2 ) -
401:      $                             REFSUM*DCONJG( V( 2, M ) )
402:                      H( J, K+3 ) = H( J, K+3 ) -
403:      $                             REFSUM*DCONJG( V( 3, M ) )
404:    50             CONTINUE
405: *
406:                   IF( ACCUM ) THEN
407: *
408: *                    ==== Accumulate U. (If necessary, update Z later
409: *                    .    with with an efficient matrix-matrix
410: *                    .    multiply.) ====
411: *
412:                      KMS = K - INCOL
413:                      DO 60 J = MAX( 1, KTOP-INCOL ), KDU
414:                         REFSUM = V( 1, M )*( U( J, KMS+1 )+V( 2, M )*
415:      $                           U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) )
416:                         U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
417:                         U( J, KMS+2 ) = U( J, KMS+2 ) -
418:      $                                  REFSUM*DCONJG( V( 2, M ) )
419:                         U( J, KMS+3 ) = U( J, KMS+3 ) -
420:      $                                  REFSUM*DCONJG( V( 3, M ) )
421:    60                CONTINUE
422:                   ELSE IF( WANTZ ) THEN
423: *
424: *                    ==== U is not accumulated, so update Z
425: *                    .    now by multiplying by reflections
426: *                    .    from the right. ====
427: *
428:                      DO 70 J = ILOZ, IHIZ
429:                         REFSUM = V( 1, M )*( Z( J, K+1 )+V( 2, M )*
430:      $                           Z( J, K+2 )+V( 3, M )*Z( J, K+3 ) )
431:                         Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
432:                         Z( J, K+2 ) = Z( J, K+2 ) -
433:      $                                REFSUM*DCONJG( V( 2, M ) )
434:                         Z( J, K+3 ) = Z( J, K+3 ) -
435:      $                                REFSUM*DCONJG( V( 3, M ) )
436:    70                CONTINUE
437:                   END IF
438:                END IF
439:    80       CONTINUE
440: *
441: *           ==== Special case: 2-by-2 reflection (if needed) ====
442: *
443:             K = KRCOL + 3*( M22-1 )
444:             IF( BMP22 .AND. ( V( 1, M22 ).NE.ZERO ) ) THEN
445:                DO 90 J = JTOP, MIN( KBOT, K+3 )
446:                   REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )*
447:      $                     H( J, K+2 ) )
448:                   H( J, K+1 ) = H( J, K+1 ) - REFSUM
449:                   H( J, K+2 ) = H( J, K+2 ) -
450:      $                          REFSUM*DCONJG( V( 2, M22 ) )
451:    90          CONTINUE
452: *
453:                IF( ACCUM ) THEN
454:                   KMS = K - INCOL
455:                   DO 100 J = MAX( 1, KTOP-INCOL ), KDU
456:                      REFSUM = V( 1, M22 )*( U( J, KMS+1 )+V( 2, M22 )*
457:      $                        U( J, KMS+2 ) )
458:                      U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
459:                      U( J, KMS+2 ) = U( J, KMS+2 ) -
460:      $                               REFSUM*DCONJG( V( 2, M22 ) )
461:   100             CONTINUE
462:                ELSE IF( WANTZ ) THEN
463:                   DO 110 J = ILOZ, IHIZ
464:                      REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )*
465:      $                        Z( J, K+2 ) )
466:                      Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
467:                      Z( J, K+2 ) = Z( J, K+2 ) -
468:      $                             REFSUM*DCONJG( V( 2, M22 ) )
469:   110             CONTINUE
470:                END IF
471:             END IF
472: *
473: *           ==== Vigilant deflation check ====
474: *
475:             MSTART = MTOP
476:             IF( KRCOL+3*( MSTART-1 ).LT.KTOP )
477:      $         MSTART = MSTART + 1
478:             MEND = MBOT
479:             IF( BMP22 )
480:      $         MEND = MEND + 1
481:             IF( KRCOL.EQ.KBOT-2 )
482:      $         MEND = MEND + 1
483:             DO 120 M = MSTART, MEND
484:                K = MIN( KBOT-1, KRCOL+3*( M-1 ) )
485: *
486: *              ==== The following convergence test requires that
487: *              .    the tradition small-compared-to-nearby-diagonals
488: *              .    criterion and the Ahues & Tisseur (LAWN 122, 1997)
489: *              .    criteria both be satisfied.  The latter improves
490: *              .    accuracy in some examples. Falling back on an
491: *              .    alternate convergence criterion when TST1 or TST2
492: *              .    is zero (as done here) is traditional but probably
493: *              .    unnecessary. ====
494: *
495:                IF( H( K+1, K ).NE.ZERO ) THEN
496:                   TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) )
497:                   IF( TST1.EQ.RZERO ) THEN
498:                      IF( K.GE.KTOP+1 )
499:      $                  TST1 = TST1 + CABS1( H( K, K-1 ) )
500:                      IF( K.GE.KTOP+2 )
501:      $                  TST1 = TST1 + CABS1( H( K, K-2 ) )
502:                      IF( K.GE.KTOP+3 )
503:      $                  TST1 = TST1 + CABS1( H( K, K-3 ) )
504:                      IF( K.LE.KBOT-2 )
505:      $                  TST1 = TST1 + CABS1( H( K+2, K+1 ) )
506:                      IF( K.LE.KBOT-3 )
507:      $                  TST1 = TST1 + CABS1( H( K+3, K+1 ) )
508:                      IF( K.LE.KBOT-4 )
509:      $                  TST1 = TST1 + CABS1( H( K+4, K+1 ) )
510:                   END IF
511:                   IF( CABS1( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) )
512:      $                 THEN
513:                      H12 = MAX( CABS1( H( K+1, K ) ),
514:      $                     CABS1( H( K, K+1 ) ) )
515:                      H21 = MIN( CABS1( H( K+1, K ) ),
516:      $                     CABS1( H( K, K+1 ) ) )
517:                      H11 = MAX( CABS1( H( K+1, K+1 ) ),
518:      $                     CABS1( H( K, K )-H( K+1, K+1 ) ) )
519:                      H22 = MIN( CABS1( H( K+1, K+1 ) ),
520:      $                     CABS1( H( K, K )-H( K+1, K+1 ) ) )
521:                      SCL = H11 + H12
522:                      TST2 = H22*( H11 / SCL )
523: *
524:                      IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE.
525:      $                   MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO
526:                   END IF
527:                END IF
528:   120       CONTINUE
529: *
530: *           ==== Fill in the last row of each bulge. ====
531: *
532:             MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )
533:             DO 130 M = MTOP, MEND
534:                K = KRCOL + 3*( M-1 )
535:                REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )
536:                H( K+4, K+1 ) = -REFSUM
537:                H( K+4, K+2 ) = -REFSUM*DCONJG( V( 2, M ) )
538:                H( K+4, K+3 ) = H( K+4, K+3 ) -
539:      $                         REFSUM*DCONJG( V( 3, M ) )
540:   130       CONTINUE
541: *
542: *           ==== End of near-the-diagonal bulge chase. ====
543: *
544:   140    CONTINUE
545: *
546: *        ==== Use U (if accumulated) to update far-from-diagonal
547: *        .    entries in H.  If required, use U to update Z as
548: *        .    well. ====
549: *
550:          IF( ACCUM ) THEN
551:             IF( WANTT ) THEN
552:                JTOP = 1
553:                JBOT = N
554:             ELSE
555:                JTOP = KTOP
556:                JBOT = KBOT
557:             END IF
558:             IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR.
559:      $          ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN
560: *
561: *              ==== Updates not exploiting the 2-by-2 block
562: *              .    structure of U.  K1 and NU keep track of
563: *              .    the location and size of U in the special
564: *              .    cases of introducing bulges and chasing
565: *              .    bulges off the bottom.  In these special
566: *              .    cases and in case the number of shifts
567: *              .    is NS = 2, there is no 2-by-2 block
568: *              .    structure to exploit.  ====
569: *
570:                K1 = MAX( 1, KTOP-INCOL )
571:                NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1
572: *
573: *              ==== Horizontal Multiply ====
574: *
575:                DO 150 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
576:                   JLEN = MIN( NH, JBOT-JCOL+1 )
577:                   CALL ZGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ),
578:      $                        LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH,
579:      $                        LDWH )
580:                   CALL ZLACPY( 'ALL', NU, JLEN, WH, LDWH,
581:      $                         H( INCOL+K1, JCOL ), LDH )
582:   150          CONTINUE
583: *
584: *              ==== Vertical multiply ====
585: *
586:                DO 160 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV
587:                   JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )
588:                   CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE,
589:      $                        H( JROW, INCOL+K1 ), LDH, U( K1, K1 ),
590:      $                        LDU, ZERO, WV, LDWV )
591:                   CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV,
592:      $                         H( JROW, INCOL+K1 ), LDH )
593:   160          CONTINUE
594: *
595: *              ==== Z multiply (also vertical) ====
596: *
597:                IF( WANTZ ) THEN
598:                   DO 170 JROW = ILOZ, IHIZ, NV
599:                      JLEN = MIN( NV, IHIZ-JROW+1 )
600:                      CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE,
601:      $                           Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ),
602:      $                           LDU, ZERO, WV, LDWV )
603:                      CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV,
604:      $                            Z( JROW, INCOL+K1 ), LDZ )
605:   170             CONTINUE
606:                END IF
607:             ELSE
608: *
609: *              ==== Updates exploiting U's 2-by-2 block structure.
610: *              .    (I2, I4, J2, J4 are the last rows and columns
611: *              .    of the blocks.) ====
612: *
613:                I2 = ( KDU+1 ) / 2
614:                I4 = KDU
615:                J2 = I4 - I2
616:                J4 = KDU
617: *
618: *              ==== KZS and KNZ deal with the band of zeros
619: *              .    along the diagonal of one of the triangular
620: *              .    blocks. ====
621: *
622:                KZS = ( J4-J2 ) - ( NS+1 )
623:                KNZ = NS + 1
624: *
625: *              ==== Horizontal multiply ====
626: *
627:                DO 180 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
628:                   JLEN = MIN( NH, JBOT-JCOL+1 )
629: *
630: *                 ==== Copy bottom of H to top+KZS of scratch ====
631: *                  (The first KZS rows get multiplied by zero.) ====
632: *
633:                   CALL ZLACPY( 'ALL', KNZ, JLEN, H( INCOL+1+J2, JCOL ),
634:      $                         LDH, WH( KZS+1, 1 ), LDWH )
635: *
636: *                 ==== Multiply by U21' ====
637: *
638:                   CALL ZLASET( 'ALL', KZS, JLEN, ZERO, ZERO, WH, LDWH )
639:                   CALL ZTRMM( 'L', 'U', 'C', 'N', KNZ, JLEN, ONE,
640:      $                        U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ),
641:      $                        LDWH )
642: *
643: *                 ==== Multiply top of H by U11' ====
644: *
645:                   CALL ZGEMM( 'C', 'N', I2, JLEN, J2, ONE, U, LDU,
646:      $                        H( INCOL+1, JCOL ), LDH, ONE, WH, LDWH )
647: *
648: *                 ==== Copy top of H to bottom of WH ====
649: *
650:                   CALL ZLACPY( 'ALL', J2, JLEN, H( INCOL+1, JCOL ), LDH,
651:      $                         WH( I2+1, 1 ), LDWH )
652: *
653: *                 ==== Multiply by U21' ====
654: *
655:                   CALL ZTRMM( 'L', 'L', 'C', 'N', J2, JLEN, ONE,
656:      $                        U( 1, I2+1 ), LDU, WH( I2+1, 1 ), LDWH )
657: *
658: *                 ==== Multiply by U22 ====
659: *
660:                   CALL ZGEMM( 'C', 'N', I4-I2, JLEN, J4-J2, ONE,
661:      $                        U( J2+1, I2+1 ), LDU,
662:      $                        H( INCOL+1+J2, JCOL ), LDH, ONE,
663:      $                        WH( I2+1, 1 ), LDWH )
664: *
665: *                 ==== Copy it back ====
666: *
667:                   CALL ZLACPY( 'ALL', KDU, JLEN, WH, LDWH,
668:      $                         H( INCOL+1, JCOL ), LDH )
669:   180          CONTINUE
670: *
671: *              ==== Vertical multiply ====
672: *
673:                DO 190 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV
674:                   JLEN = MIN( NV, MAX( INCOL, KTOP )-JROW )
675: *
676: *                 ==== Copy right of H to scratch (the first KZS
677: *                 .    columns get multiplied by zero) ====
678: *
679:                   CALL ZLACPY( 'ALL', JLEN, KNZ, H( JROW, INCOL+1+J2 ),
680:      $                         LDH, WV( 1, 1+KZS ), LDWV )
681: *
682: *                 ==== Multiply by U21 ====
683: *
684:                   CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, LDWV )
685:                   CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
686:      $                        U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
687:      $                        LDWV )
688: *
689: *                 ==== Multiply by U11 ====
690: *
691:                   CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE,
692:      $                        H( JROW, INCOL+1 ), LDH, U, LDU, ONE, WV,
693:      $                        LDWV )
694: *
695: *                 ==== Copy left of H to right of scratch ====
696: *
697:                   CALL ZLACPY( 'ALL', JLEN, J2, H( JROW, INCOL+1 ), LDH,
698:      $                         WV( 1, 1+I2 ), LDWV )
699: *
700: *                 ==== Multiply by U21 ====
701: *
702:                   CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
703:      $                        U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV )
704: *
705: *                 ==== Multiply by U22 ====
706: *
707:                   CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
708:      $                        H( JROW, INCOL+1+J2 ), LDH,
709:      $                        U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ),
710:      $                        LDWV )
711: *
712: *                 ==== Copy it back ====
713: *
714:                   CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV,
715:      $                         H( JROW, INCOL+1 ), LDH )
716:   190          CONTINUE
717: *
718: *              ==== Multiply Z (also vertical) ====
719: *
720:                IF( WANTZ ) THEN
721:                   DO 200 JROW = ILOZ, IHIZ, NV
722:                      JLEN = MIN( NV, IHIZ-JROW+1 )
723: *
724: *                    ==== Copy right of Z to left of scratch (first
725: *                    .     KZS columns get multiplied by zero) ====
726: *
727:                      CALL ZLACPY( 'ALL', JLEN, KNZ,
728:      $                            Z( JROW, INCOL+1+J2 ), LDZ,
729:      $                            WV( 1, 1+KZS ), LDWV )
730: *
731: *                    ==== Multiply by U12 ====
732: *
733:                      CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV,
734:      $                            LDWV )
735:                      CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
736:      $                           U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
737:      $                           LDWV )
738: *
739: *                    ==== Multiply by U11 ====
740: *
741:                      CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE,
742:      $                           Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE,
743:      $                           WV, LDWV )
744: *
745: *                    ==== Copy left of Z to right of scratch ====
746: *
747:                      CALL ZLACPY( 'ALL', JLEN, J2, Z( JROW, INCOL+1 ),
748:      $                            LDZ, WV( 1, 1+I2 ), LDWV )
749: *
750: *                    ==== Multiply by U21 ====
751: *
752:                      CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
753:      $                           U( 1, I2+1 ), LDU, WV( 1, 1+I2 ),
754:      $                           LDWV )
755: *
756: *                    ==== Multiply by U22 ====
757: *
758:                      CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
759:      $                           Z( JROW, INCOL+1+J2 ), LDZ,
760:      $                           U( J2+1, I2+1 ), LDU, ONE,
761:      $                           WV( 1, 1+I2 ), LDWV )
762: *
763: *                    ==== Copy the result back to Z ====
764: *
765:                      CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV,
766:      $                            Z( JROW, INCOL+1 ), LDZ )
767:   200             CONTINUE
768:                END IF
769:             END IF
770:          END IF
771:   210 CONTINUE
772: *
773: *     ==== End of ZLAQR5 ====
774: *
775:       END
776: