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