![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaqr5.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zlaqr5.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S,
2: $ H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV,
3: $ WV, LDWV, NH, WH, LDWH )
4: *
5: * -- LAPACK auxiliary routine (version 3.2) --
6: * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
11: $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
12: LOGICAL WANTT, WANTZ
13: * ..
14: * .. Array Arguments ..
15: COMPLEX*16 H( LDH, * ), S( * ), U( LDU, * ), V( LDV, * ),
16: $ WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
17: * ..
18: *
19: * This auxiliary subroutine called by ZLAQR0 performs a
20: * single small-bulge multi-shift QR sweep.
21: *
22: * WANTT (input) logical scalar
23: * WANTT = .true. if the triangular Schur factor
24: * is being computed. WANTT is set to .false. otherwise.
25: *
26: * WANTZ (input) logical scalar
27: * WANTZ = .true. if the unitary Schur factor is being
28: * computed. WANTZ is set to .false. otherwise.
29: *
30: * KACC22 (input) integer with value 0, 1, or 2.
31: * Specifies the computation mode of far-from-diagonal
32: * orthogonal updates.
33: * = 0: ZLAQR5 does not accumulate reflections and does not
34: * use matrix-matrix multiply to update far-from-diagonal
35: * matrix entries.
36: * = 1: ZLAQR5 accumulates reflections and uses matrix-matrix
37: * multiply to update the far-from-diagonal matrix entries.
38: * = 2: ZLAQR5 accumulates reflections, uses matrix-matrix
39: * multiply to update the far-from-diagonal matrix entries,
40: * and takes advantage of 2-by-2 block structure during
41: * matrix multiplies.
42: *
43: * N (input) integer scalar
44: * N is the order of the Hessenberg matrix H upon which this
45: * subroutine operates.
46: *
47: * KTOP (input) integer scalar
48: * KBOT (input) integer scalar
49: * These are the first and last rows and columns of an
50: * isolated diagonal block upon which the QR sweep is to be
51: * applied. It is assumed without a check that
52: * either KTOP = 1 or H(KTOP,KTOP-1) = 0
53: * and
54: * either KBOT = N or H(KBOT+1,KBOT) = 0.
55: *
56: * NSHFTS (input) integer scalar
57: * NSHFTS gives the number of simultaneous shifts. NSHFTS
58: * must be positive and even.
59: *
60: * S (input/output) COMPLEX*16 array of size (NSHFTS)
61: * S contains the shifts of origin that define the multi-
62: * shift QR sweep. On output S may be reordered.
63: *
64: * H (input/output) COMPLEX*16 array of size (LDH,N)
65: * On input H contains a Hessenberg matrix. On output a
66: * multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
67: * to the isolated diagonal block in rows and columns KTOP
68: * through KBOT.
69: *
70: * LDH (input) integer scalar
71: * LDH is the leading dimension of H just as declared in the
72: * calling procedure. LDH.GE.MAX(1,N).
73: *
74: * ILOZ (input) INTEGER
75: * IHIZ (input) INTEGER
76: * Specify the rows of Z to which transformations must be
77: * applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
78: *
79: * Z (input/output) COMPLEX*16 array of size (LDZ,IHI)
80: * If WANTZ = .TRUE., then the QR Sweep unitary
81: * similarity transformation is accumulated into
82: * Z(ILOZ:IHIZ,ILO:IHI) from the right.
83: * If WANTZ = .FALSE., then Z is unreferenced.
84: *
85: * LDZ (input) integer scalar
86: * LDA is the leading dimension of Z just as declared in
87: * the calling procedure. LDZ.GE.N.
88: *
89: * V (workspace) COMPLEX*16 array of size (LDV,NSHFTS/2)
90: *
91: * LDV (input) integer scalar
92: * LDV is the leading dimension of V as declared in the
93: * calling procedure. LDV.GE.3.
94: *
95: * U (workspace) COMPLEX*16 array of size
96: * (LDU,3*NSHFTS-3)
97: *
98: * LDU (input) integer scalar
99: * 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