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.3.0) --
6 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
7 * November 2010
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 ) THEN
445 IF ( V( 1, M22 ).NE.ZERO ) THEN
446 DO 90 J = JTOP, MIN( KBOT, K+3 )
447 REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )*
448 $ H( J, K+2 ) )
449 H( J, K+1 ) = H( J, K+1 ) - REFSUM
450 H( J, K+2 ) = H( J, K+2 ) -
451 $ REFSUM*DCONJG( V( 2, M22 ) )
452 90 CONTINUE
453 *
454 IF( ACCUM ) THEN
455 KMS = K - INCOL
456 DO 100 J = MAX( 1, KTOP-INCOL ), KDU
457 REFSUM = V( 1, M22 )*( U( J, KMS+1 )+
458 $ V( 2, M22 )*U( J, KMS+2 ) )
459 U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
460 U( J, KMS+2 ) = U( J, KMS+2 ) -
461 $ REFSUM*DCONJG( V( 2, M22 ) )
462 100 CONTINUE
463 ELSE IF( WANTZ ) THEN
464 DO 110 J = ILOZ, IHIZ
465 REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )*
466 $ Z( J, K+2 ) )
467 Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
468 Z( J, K+2 ) = Z( J, K+2 ) -
469 $ REFSUM*DCONJG( V( 2, M22 ) )
470 110 CONTINUE
471 END IF
472 END IF
473 END IF
474 *
475 * ==== Vigilant deflation check ====
476 *
477 MSTART = MTOP
478 IF( KRCOL+3*( MSTART-1 ).LT.KTOP )
479 $ MSTART = MSTART + 1
480 MEND = MBOT
481 IF( BMP22 )
482 $ MEND = MEND + 1
483 IF( KRCOL.EQ.KBOT-2 )
484 $ MEND = MEND + 1
485 DO 120 M = MSTART, MEND
486 K = MIN( KBOT-1, KRCOL+3*( M-1 ) )
487 *
488 * ==== The following convergence test requires that
489 * . the tradition small-compared-to-nearby-diagonals
490 * . criterion and the Ahues & Tisseur (LAWN 122, 1997)
491 * . criteria both be satisfied. The latter improves
492 * . accuracy in some examples. Falling back on an
493 * . alternate convergence criterion when TST1 or TST2
494 * . is zero (as done here) is traditional but probably
495 * . unnecessary. ====
496 *
497 IF( H( K+1, K ).NE.ZERO ) THEN
498 TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) )
499 IF( TST1.EQ.RZERO ) THEN
500 IF( K.GE.KTOP+1 )
501 $ TST1 = TST1 + CABS1( H( K, K-1 ) )
502 IF( K.GE.KTOP+2 )
503 $ TST1 = TST1 + CABS1( H( K, K-2 ) )
504 IF( K.GE.KTOP+3 )
505 $ TST1 = TST1 + CABS1( H( K, K-3 ) )
506 IF( K.LE.KBOT-2 )
507 $ TST1 = TST1 + CABS1( H( K+2, K+1 ) )
508 IF( K.LE.KBOT-3 )
509 $ TST1 = TST1 + CABS1( H( K+3, K+1 ) )
510 IF( K.LE.KBOT-4 )
511 $ TST1 = TST1 + CABS1( H( K+4, K+1 ) )
512 END IF
513 IF( CABS1( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) )
514 $ THEN
515 H12 = MAX( CABS1( H( K+1, K ) ),
516 $ CABS1( H( K, K+1 ) ) )
517 H21 = MIN( CABS1( H( K+1, K ) ),
518 $ CABS1( H( K, K+1 ) ) )
519 H11 = MAX( CABS1( H( K+1, K+1 ) ),
520 $ CABS1( H( K, K )-H( K+1, K+1 ) ) )
521 H22 = MIN( CABS1( H( K+1, K+1 ) ),
522 $ CABS1( H( K, K )-H( K+1, K+1 ) ) )
523 SCL = H11 + H12
524 TST2 = H22*( H11 / SCL )
525 *
526 IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE.
527 $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO
528 END IF
529 END IF
530 120 CONTINUE
531 *
532 * ==== Fill in the last row of each bulge. ====
533 *
534 MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )
535 DO 130 M = MTOP, MEND
536 K = KRCOL + 3*( M-1 )
537 REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )
538 H( K+4, K+1 ) = -REFSUM
539 H( K+4, K+2 ) = -REFSUM*DCONJG( V( 2, M ) )
540 H( K+4, K+3 ) = H( K+4, K+3 ) -
541 $ REFSUM*DCONJG( V( 3, M ) )
542 130 CONTINUE
543 *
544 * ==== End of near-the-diagonal bulge chase. ====
545 *
546 140 CONTINUE
547 *
548 * ==== Use U (if accumulated) to update far-from-diagonal
549 * . entries in H. If required, use U to update Z as
550 * . well. ====
551 *
552 IF( ACCUM ) THEN
553 IF( WANTT ) THEN
554 JTOP = 1
555 JBOT = N
556 ELSE
557 JTOP = KTOP
558 JBOT = KBOT
559 END IF
560 IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR.
561 $ ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN
562 *
563 * ==== Updates not exploiting the 2-by-2 block
564 * . structure of U. K1 and NU keep track of
565 * . the location and size of U in the special
566 * . cases of introducing bulges and chasing
567 * . bulges off the bottom. In these special
568 * . cases and in case the number of shifts
569 * . is NS = 2, there is no 2-by-2 block
570 * . structure to exploit. ====
571 *
572 K1 = MAX( 1, KTOP-INCOL )
573 NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1
574 *
575 * ==== Horizontal Multiply ====
576 *
577 DO 150 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
578 JLEN = MIN( NH, JBOT-JCOL+1 )
579 CALL ZGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ),
580 $ LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH,
581 $ LDWH )
582 CALL ZLACPY( 'ALL', NU, JLEN, WH, LDWH,
583 $ H( INCOL+K1, JCOL ), LDH )
584 150 CONTINUE
585 *
586 * ==== Vertical multiply ====
587 *
588 DO 160 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV
589 JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )
590 CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE,
591 $ H( JROW, INCOL+K1 ), LDH, U( K1, K1 ),
592 $ LDU, ZERO, WV, LDWV )
593 CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV,
594 $ H( JROW, INCOL+K1 ), LDH )
595 160 CONTINUE
596 *
597 * ==== Z multiply (also vertical) ====
598 *
599 IF( WANTZ ) THEN
600 DO 170 JROW = ILOZ, IHIZ, NV
601 JLEN = MIN( NV, IHIZ-JROW+1 )
602 CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE,
603 $ Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ),
604 $ LDU, ZERO, WV, LDWV )
605 CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV,
606 $ Z( JROW, INCOL+K1 ), LDZ )
607 170 CONTINUE
608 END IF
609 ELSE
610 *
611 * ==== Updates exploiting U's 2-by-2 block structure.
612 * . (I2, I4, J2, J4 are the last rows and columns
613 * . of the blocks.) ====
614 *
615 I2 = ( KDU+1 ) / 2
616 I4 = KDU
617 J2 = I4 - I2
618 J4 = KDU
619 *
620 * ==== KZS and KNZ deal with the band of zeros
621 * . along the diagonal of one of the triangular
622 * . blocks. ====
623 *
624 KZS = ( J4-J2 ) - ( NS+1 )
625 KNZ = NS + 1
626 *
627 * ==== Horizontal multiply ====
628 *
629 DO 180 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
630 JLEN = MIN( NH, JBOT-JCOL+1 )
631 *
632 * ==== Copy bottom of H to top+KZS of scratch ====
633 * (The first KZS rows get multiplied by zero.) ====
634 *
635 CALL ZLACPY( 'ALL', KNZ, JLEN, H( INCOL+1+J2, JCOL ),
636 $ LDH, WH( KZS+1, 1 ), LDWH )
637 *
638 * ==== Multiply by U21**H ====
639 *
640 CALL ZLASET( 'ALL', KZS, JLEN, ZERO, ZERO, WH, LDWH )
641 CALL ZTRMM( 'L', 'U', 'C', 'N', KNZ, JLEN, ONE,
642 $ U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ),
643 $ LDWH )
644 *
645 * ==== Multiply top of H by U11**H ====
646 *
647 CALL ZGEMM( 'C', 'N', I2, JLEN, J2, ONE, U, LDU,
648 $ H( INCOL+1, JCOL ), LDH, ONE, WH, LDWH )
649 *
650 * ==== Copy top of H to bottom of WH ====
651 *
652 CALL ZLACPY( 'ALL', J2, JLEN, H( INCOL+1, JCOL ), LDH,
653 $ WH( I2+1, 1 ), LDWH )
654 *
655 * ==== Multiply by U21**H ====
656 *
657 CALL ZTRMM( 'L', 'L', 'C', 'N', J2, JLEN, ONE,
658 $ U( 1, I2+1 ), LDU, WH( I2+1, 1 ), LDWH )
659 *
660 * ==== Multiply by U22 ====
661 *
662 CALL ZGEMM( 'C', 'N', I4-I2, JLEN, J4-J2, ONE,
663 $ U( J2+1, I2+1 ), LDU,
664 $ H( INCOL+1+J2, JCOL ), LDH, ONE,
665 $ WH( I2+1, 1 ), LDWH )
666 *
667 * ==== Copy it back ====
668 *
669 CALL ZLACPY( 'ALL', KDU, JLEN, WH, LDWH,
670 $ H( INCOL+1, JCOL ), LDH )
671 180 CONTINUE
672 *
673 * ==== Vertical multiply ====
674 *
675 DO 190 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV
676 JLEN = MIN( NV, MAX( INCOL, KTOP )-JROW )
677 *
678 * ==== Copy right of H to scratch (the first KZS
679 * . columns get multiplied by zero) ====
680 *
681 CALL ZLACPY( 'ALL', JLEN, KNZ, H( JROW, INCOL+1+J2 ),
682 $ LDH, WV( 1, 1+KZS ), LDWV )
683 *
684 * ==== Multiply by U21 ====
685 *
686 CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, LDWV )
687 CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
688 $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
689 $ LDWV )
690 *
691 * ==== Multiply by U11 ====
692 *
693 CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE,
694 $ H( JROW, INCOL+1 ), LDH, U, LDU, ONE, WV,
695 $ LDWV )
696 *
697 * ==== Copy left of H to right of scratch ====
698 *
699 CALL ZLACPY( 'ALL', JLEN, J2, H( JROW, INCOL+1 ), LDH,
700 $ WV( 1, 1+I2 ), LDWV )
701 *
702 * ==== Multiply by U21 ====
703 *
704 CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
705 $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV )
706 *
707 * ==== Multiply by U22 ====
708 *
709 CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
710 $ H( JROW, INCOL+1+J2 ), LDH,
711 $ U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ),
712 $ LDWV )
713 *
714 * ==== Copy it back ====
715 *
716 CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV,
717 $ H( JROW, INCOL+1 ), LDH )
718 190 CONTINUE
719 *
720 * ==== Multiply Z (also vertical) ====
721 *
722 IF( WANTZ ) THEN
723 DO 200 JROW = ILOZ, IHIZ, NV
724 JLEN = MIN( NV, IHIZ-JROW+1 )
725 *
726 * ==== Copy right of Z to left of scratch (first
727 * . KZS columns get multiplied by zero) ====
728 *
729 CALL ZLACPY( 'ALL', JLEN, KNZ,
730 $ Z( JROW, INCOL+1+J2 ), LDZ,
731 $ WV( 1, 1+KZS ), LDWV )
732 *
733 * ==== Multiply by U12 ====
734 *
735 CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV,
736 $ LDWV )
737 CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
738 $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
739 $ LDWV )
740 *
741 * ==== Multiply by U11 ====
742 *
743 CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE,
744 $ Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE,
745 $ WV, LDWV )
746 *
747 * ==== Copy left of Z to right of scratch ====
748 *
749 CALL ZLACPY( 'ALL', JLEN, J2, Z( JROW, INCOL+1 ),
750 $ LDZ, WV( 1, 1+I2 ), LDWV )
751 *
752 * ==== Multiply by U21 ====
753 *
754 CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
755 $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ),
756 $ LDWV )
757 *
758 * ==== Multiply by U22 ====
759 *
760 CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
761 $ Z( JROW, INCOL+1+J2 ), LDZ,
762 $ U( J2+1, I2+1 ), LDU, ONE,
763 $ WV( 1, 1+I2 ), LDWV )
764 *
765 * ==== Copy the result back to Z ====
766 *
767 CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV,
768 $ Z( JROW, INCOL+1 ), LDZ )
769 200 CONTINUE
770 END IF
771 END IF
772 END IF
773 210 CONTINUE
774 *
775 * ==== End of ZLAQR5 ====
776 *
777 END
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.3.0) --
6 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
7 * November 2010
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 ) THEN
445 IF ( V( 1, M22 ).NE.ZERO ) THEN
446 DO 90 J = JTOP, MIN( KBOT, K+3 )
447 REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )*
448 $ H( J, K+2 ) )
449 H( J, K+1 ) = H( J, K+1 ) - REFSUM
450 H( J, K+2 ) = H( J, K+2 ) -
451 $ REFSUM*DCONJG( V( 2, M22 ) )
452 90 CONTINUE
453 *
454 IF( ACCUM ) THEN
455 KMS = K - INCOL
456 DO 100 J = MAX( 1, KTOP-INCOL ), KDU
457 REFSUM = V( 1, M22 )*( U( J, KMS+1 )+
458 $ V( 2, M22 )*U( J, KMS+2 ) )
459 U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
460 U( J, KMS+2 ) = U( J, KMS+2 ) -
461 $ REFSUM*DCONJG( V( 2, M22 ) )
462 100 CONTINUE
463 ELSE IF( WANTZ ) THEN
464 DO 110 J = ILOZ, IHIZ
465 REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )*
466 $ Z( J, K+2 ) )
467 Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
468 Z( J, K+2 ) = Z( J, K+2 ) -
469 $ REFSUM*DCONJG( V( 2, M22 ) )
470 110 CONTINUE
471 END IF
472 END IF
473 END IF
474 *
475 * ==== Vigilant deflation check ====
476 *
477 MSTART = MTOP
478 IF( KRCOL+3*( MSTART-1 ).LT.KTOP )
479 $ MSTART = MSTART + 1
480 MEND = MBOT
481 IF( BMP22 )
482 $ MEND = MEND + 1
483 IF( KRCOL.EQ.KBOT-2 )
484 $ MEND = MEND + 1
485 DO 120 M = MSTART, MEND
486 K = MIN( KBOT-1, KRCOL+3*( M-1 ) )
487 *
488 * ==== The following convergence test requires that
489 * . the tradition small-compared-to-nearby-diagonals
490 * . criterion and the Ahues & Tisseur (LAWN 122, 1997)
491 * . criteria both be satisfied. The latter improves
492 * . accuracy in some examples. Falling back on an
493 * . alternate convergence criterion when TST1 or TST2
494 * . is zero (as done here) is traditional but probably
495 * . unnecessary. ====
496 *
497 IF( H( K+1, K ).NE.ZERO ) THEN
498 TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) )
499 IF( TST1.EQ.RZERO ) THEN
500 IF( K.GE.KTOP+1 )
501 $ TST1 = TST1 + CABS1( H( K, K-1 ) )
502 IF( K.GE.KTOP+2 )
503 $ TST1 = TST1 + CABS1( H( K, K-2 ) )
504 IF( K.GE.KTOP+3 )
505 $ TST1 = TST1 + CABS1( H( K, K-3 ) )
506 IF( K.LE.KBOT-2 )
507 $ TST1 = TST1 + CABS1( H( K+2, K+1 ) )
508 IF( K.LE.KBOT-3 )
509 $ TST1 = TST1 + CABS1( H( K+3, K+1 ) )
510 IF( K.LE.KBOT-4 )
511 $ TST1 = TST1 + CABS1( H( K+4, K+1 ) )
512 END IF
513 IF( CABS1( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) )
514 $ THEN
515 H12 = MAX( CABS1( H( K+1, K ) ),
516 $ CABS1( H( K, K+1 ) ) )
517 H21 = MIN( CABS1( H( K+1, K ) ),
518 $ CABS1( H( K, K+1 ) ) )
519 H11 = MAX( CABS1( H( K+1, K+1 ) ),
520 $ CABS1( H( K, K )-H( K+1, K+1 ) ) )
521 H22 = MIN( CABS1( H( K+1, K+1 ) ),
522 $ CABS1( H( K, K )-H( K+1, K+1 ) ) )
523 SCL = H11 + H12
524 TST2 = H22*( H11 / SCL )
525 *
526 IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE.
527 $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO
528 END IF
529 END IF
530 120 CONTINUE
531 *
532 * ==== Fill in the last row of each bulge. ====
533 *
534 MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )
535 DO 130 M = MTOP, MEND
536 K = KRCOL + 3*( M-1 )
537 REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )
538 H( K+4, K+1 ) = -REFSUM
539 H( K+4, K+2 ) = -REFSUM*DCONJG( V( 2, M ) )
540 H( K+4, K+3 ) = H( K+4, K+3 ) -
541 $ REFSUM*DCONJG( V( 3, M ) )
542 130 CONTINUE
543 *
544 * ==== End of near-the-diagonal bulge chase. ====
545 *
546 140 CONTINUE
547 *
548 * ==== Use U (if accumulated) to update far-from-diagonal
549 * . entries in H. If required, use U to update Z as
550 * . well. ====
551 *
552 IF( ACCUM ) THEN
553 IF( WANTT ) THEN
554 JTOP = 1
555 JBOT = N
556 ELSE
557 JTOP = KTOP
558 JBOT = KBOT
559 END IF
560 IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR.
561 $ ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN
562 *
563 * ==== Updates not exploiting the 2-by-2 block
564 * . structure of U. K1 and NU keep track of
565 * . the location and size of U in the special
566 * . cases of introducing bulges and chasing
567 * . bulges off the bottom. In these special
568 * . cases and in case the number of shifts
569 * . is NS = 2, there is no 2-by-2 block
570 * . structure to exploit. ====
571 *
572 K1 = MAX( 1, KTOP-INCOL )
573 NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1
574 *
575 * ==== Horizontal Multiply ====
576 *
577 DO 150 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
578 JLEN = MIN( NH, JBOT-JCOL+1 )
579 CALL ZGEMM( 'C', 'N', NU, JLEN, NU, ONE, U( K1, K1 ),
580 $ LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH,
581 $ LDWH )
582 CALL ZLACPY( 'ALL', NU, JLEN, WH, LDWH,
583 $ H( INCOL+K1, JCOL ), LDH )
584 150 CONTINUE
585 *
586 * ==== Vertical multiply ====
587 *
588 DO 160 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV
589 JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )
590 CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE,
591 $ H( JROW, INCOL+K1 ), LDH, U( K1, K1 ),
592 $ LDU, ZERO, WV, LDWV )
593 CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV,
594 $ H( JROW, INCOL+K1 ), LDH )
595 160 CONTINUE
596 *
597 * ==== Z multiply (also vertical) ====
598 *
599 IF( WANTZ ) THEN
600 DO 170 JROW = ILOZ, IHIZ, NV
601 JLEN = MIN( NV, IHIZ-JROW+1 )
602 CALL ZGEMM( 'N', 'N', JLEN, NU, NU, ONE,
603 $ Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ),
604 $ LDU, ZERO, WV, LDWV )
605 CALL ZLACPY( 'ALL', JLEN, NU, WV, LDWV,
606 $ Z( JROW, INCOL+K1 ), LDZ )
607 170 CONTINUE
608 END IF
609 ELSE
610 *
611 * ==== Updates exploiting U's 2-by-2 block structure.
612 * . (I2, I4, J2, J4 are the last rows and columns
613 * . of the blocks.) ====
614 *
615 I2 = ( KDU+1 ) / 2
616 I4 = KDU
617 J2 = I4 - I2
618 J4 = KDU
619 *
620 * ==== KZS and KNZ deal with the band of zeros
621 * . along the diagonal of one of the triangular
622 * . blocks. ====
623 *
624 KZS = ( J4-J2 ) - ( NS+1 )
625 KNZ = NS + 1
626 *
627 * ==== Horizontal multiply ====
628 *
629 DO 180 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
630 JLEN = MIN( NH, JBOT-JCOL+1 )
631 *
632 * ==== Copy bottom of H to top+KZS of scratch ====
633 * (The first KZS rows get multiplied by zero.) ====
634 *
635 CALL ZLACPY( 'ALL', KNZ, JLEN, H( INCOL+1+J2, JCOL ),
636 $ LDH, WH( KZS+1, 1 ), LDWH )
637 *
638 * ==== Multiply by U21**H ====
639 *
640 CALL ZLASET( 'ALL', KZS, JLEN, ZERO, ZERO, WH, LDWH )
641 CALL ZTRMM( 'L', 'U', 'C', 'N', KNZ, JLEN, ONE,
642 $ U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ),
643 $ LDWH )
644 *
645 * ==== Multiply top of H by U11**H ====
646 *
647 CALL ZGEMM( 'C', 'N', I2, JLEN, J2, ONE, U, LDU,
648 $ H( INCOL+1, JCOL ), LDH, ONE, WH, LDWH )
649 *
650 * ==== Copy top of H to bottom of WH ====
651 *
652 CALL ZLACPY( 'ALL', J2, JLEN, H( INCOL+1, JCOL ), LDH,
653 $ WH( I2+1, 1 ), LDWH )
654 *
655 * ==== Multiply by U21**H ====
656 *
657 CALL ZTRMM( 'L', 'L', 'C', 'N', J2, JLEN, ONE,
658 $ U( 1, I2+1 ), LDU, WH( I2+1, 1 ), LDWH )
659 *
660 * ==== Multiply by U22 ====
661 *
662 CALL ZGEMM( 'C', 'N', I4-I2, JLEN, J4-J2, ONE,
663 $ U( J2+1, I2+1 ), LDU,
664 $ H( INCOL+1+J2, JCOL ), LDH, ONE,
665 $ WH( I2+1, 1 ), LDWH )
666 *
667 * ==== Copy it back ====
668 *
669 CALL ZLACPY( 'ALL', KDU, JLEN, WH, LDWH,
670 $ H( INCOL+1, JCOL ), LDH )
671 180 CONTINUE
672 *
673 * ==== Vertical multiply ====
674 *
675 DO 190 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV
676 JLEN = MIN( NV, MAX( INCOL, KTOP )-JROW )
677 *
678 * ==== Copy right of H to scratch (the first KZS
679 * . columns get multiplied by zero) ====
680 *
681 CALL ZLACPY( 'ALL', JLEN, KNZ, H( JROW, INCOL+1+J2 ),
682 $ LDH, WV( 1, 1+KZS ), LDWV )
683 *
684 * ==== Multiply by U21 ====
685 *
686 CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV, LDWV )
687 CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
688 $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
689 $ LDWV )
690 *
691 * ==== Multiply by U11 ====
692 *
693 CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE,
694 $ H( JROW, INCOL+1 ), LDH, U, LDU, ONE, WV,
695 $ LDWV )
696 *
697 * ==== Copy left of H to right of scratch ====
698 *
699 CALL ZLACPY( 'ALL', JLEN, J2, H( JROW, INCOL+1 ), LDH,
700 $ WV( 1, 1+I2 ), LDWV )
701 *
702 * ==== Multiply by U21 ====
703 *
704 CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
705 $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV )
706 *
707 * ==== Multiply by U22 ====
708 *
709 CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
710 $ H( JROW, INCOL+1+J2 ), LDH,
711 $ U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ),
712 $ LDWV )
713 *
714 * ==== Copy it back ====
715 *
716 CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV,
717 $ H( JROW, INCOL+1 ), LDH )
718 190 CONTINUE
719 *
720 * ==== Multiply Z (also vertical) ====
721 *
722 IF( WANTZ ) THEN
723 DO 200 JROW = ILOZ, IHIZ, NV
724 JLEN = MIN( NV, IHIZ-JROW+1 )
725 *
726 * ==== Copy right of Z to left of scratch (first
727 * . KZS columns get multiplied by zero) ====
728 *
729 CALL ZLACPY( 'ALL', JLEN, KNZ,
730 $ Z( JROW, INCOL+1+J2 ), LDZ,
731 $ WV( 1, 1+KZS ), LDWV )
732 *
733 * ==== Multiply by U12 ====
734 *
735 CALL ZLASET( 'ALL', JLEN, KZS, ZERO, ZERO, WV,
736 $ LDWV )
737 CALL ZTRMM( 'R', 'U', 'N', 'N', JLEN, KNZ, ONE,
738 $ U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
739 $ LDWV )
740 *
741 * ==== Multiply by U11 ====
742 *
743 CALL ZGEMM( 'N', 'N', JLEN, I2, J2, ONE,
744 $ Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE,
745 $ WV, LDWV )
746 *
747 * ==== Copy left of Z to right of scratch ====
748 *
749 CALL ZLACPY( 'ALL', JLEN, J2, Z( JROW, INCOL+1 ),
750 $ LDZ, WV( 1, 1+I2 ), LDWV )
751 *
752 * ==== Multiply by U21 ====
753 *
754 CALL ZTRMM( 'R', 'L', 'N', 'N', JLEN, I4-I2, ONE,
755 $ U( 1, I2+1 ), LDU, WV( 1, 1+I2 ),
756 $ LDWV )
757 *
758 * ==== Multiply by U22 ====
759 *
760 CALL ZGEMM( 'N', 'N', JLEN, I4-I2, J4-J2, ONE,
761 $ Z( JROW, INCOL+1+J2 ), LDZ,
762 $ U( J2+1, I2+1 ), LDU, ONE,
763 $ WV( 1, 1+I2 ), LDWV )
764 *
765 * ==== Copy the result back to Z ====
766 *
767 CALL ZLACPY( 'ALL', JLEN, KDU, WV, LDWV,
768 $ Z( JROW, INCOL+1 ), LDZ )
769 200 CONTINUE
770 END IF
771 END IF
772 END IF
773 210 CONTINUE
774 *
775 * ==== End of ZLAQR5 ====
776 *
777 END