LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
zlaqr5.f
Go to the documentation of this file.
1 *> \brief \b ZLAQR5 performs a single small-bulge multi-shift QR sweep.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZLAQR5 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqr5.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr5.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr5.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S,
22 * H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV,
23 * WV, LDWV, NH, WH, LDWH )
24 *
25 * .. Scalar Arguments ..
26 * INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
27 * $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
28 * LOGICAL WANTT, WANTZ
29 * ..
30 * .. Array Arguments ..
31 * COMPLEX*16 H( LDH, * ), S( * ), U( LDU, * ), V( LDV, * ),
32 * $ WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> ZLAQR5, called by ZLAQR0, performs a
42 *> single small-bulge multi-shift QR sweep.
43 *> \endverbatim
44 *
45 * Arguments:
46 * ==========
47 *
48 *> \param[in] WANTT
49 *> \verbatim
50 *> WANTT is LOGICAL
51 *> WANTT = .true. if the triangular Schur factor
52 *> is being computed. WANTT is set to .false. otherwise.
53 *> \endverbatim
54 *>
55 *> \param[in] WANTZ
56 *> \verbatim
57 *> WANTZ is LOGICAL
58 *> WANTZ = .true. if the unitary Schur factor is being
59 *> computed. WANTZ is set to .false. otherwise.
60 *> \endverbatim
61 *>
62 *> \param[in] KACC22
63 *> \verbatim
64 *> KACC22 is INTEGER with value 0, 1, or 2.
65 *> Specifies the computation mode of far-from-diagonal
66 *> orthogonal updates.
67 *> = 0: ZLAQR5 does not accumulate reflections and does not
68 *> use matrix-matrix multiply to update far-from-diagonal
69 *> matrix entries.
70 *> = 1: ZLAQR5 accumulates reflections and uses matrix-matrix
71 *> multiply to update the far-from-diagonal matrix entries.
72 *> = 2: ZLAQR5 accumulates reflections, uses matrix-matrix
73 *> multiply to update the far-from-diagonal matrix entries,
74 *> and takes advantage of 2-by-2 block structure during
75 *> matrix multiplies.
76 *> \endverbatim
77 *>
78 *> \param[in] N
79 *> \verbatim
80 *> N is INTEGER
81 *> N is the order of the Hessenberg matrix H upon which this
82 *> subroutine operates.
83 *> \endverbatim
84 *>
85 *> \param[in] KTOP
86 *> \verbatim
87 *> KTOP is INTEGER
88 *> \endverbatim
89 *>
90 *> \param[in] KBOT
91 *> \verbatim
92 *> KBOT is INTEGER
93 *> These are the first and last rows and columns of an
94 *> isolated diagonal block upon which the QR sweep is to be
95 *> applied. It is assumed without a check that
96 *> either KTOP = 1 or H(KTOP,KTOP-1) = 0
97 *> and
98 *> either KBOT = N or H(KBOT+1,KBOT) = 0.
99 *> \endverbatim
100 *>
101 *> \param[in] NSHFTS
102 *> \verbatim
103 *> NSHFTS is INTEGER
104 *> NSHFTS gives the number of simultaneous shifts. NSHFTS
105 *> must be positive and even.
106 *> \endverbatim
107 *>
108 *> \param[in,out] S
109 *> \verbatim
110 *> S is COMPLEX*16 array, dimension (NSHFTS)
111 *> S contains the shifts of origin that define the multi-
112 *> shift QR sweep. On output S may be reordered.
113 *> \endverbatim
114 *>
115 *> \param[in,out] H
116 *> \verbatim
117 *> H is COMPLEX*16 array, dimension (LDH,N)
118 *> On input H contains a Hessenberg matrix. On output a
119 *> multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
120 *> to the isolated diagonal block in rows and columns KTOP
121 *> through KBOT.
122 *> \endverbatim
123 *>
124 *> \param[in] LDH
125 *> \verbatim
126 *> LDH is INTEGER
127 *> LDH is the leading dimension of H just as declared in the
128 *> calling procedure. LDH >= MAX(1,N).
129 *> \endverbatim
130 *>
131 *> \param[in] ILOZ
132 *> \verbatim
133 *> ILOZ is INTEGER
134 *> \endverbatim
135 *>
136 *> \param[in] IHIZ
137 *> \verbatim
138 *> IHIZ is INTEGER
139 *> Specify the rows of Z to which transformations must be
140 *> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N
141 *> \endverbatim
142 *>
143 *> \param[in,out] Z
144 *> \verbatim
145 *> Z is COMPLEX*16 array, dimension (LDZ,IHIZ)
146 *> If WANTZ = .TRUE., then the QR Sweep unitary
147 *> similarity transformation is accumulated into
148 *> Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
149 *> If WANTZ = .FALSE., then Z is unreferenced.
150 *> \endverbatim
151 *>
152 *> \param[in] LDZ
153 *> \verbatim
154 *> LDZ is INTEGER
155 *> LDA is the leading dimension of Z just as declared in
156 *> the calling procedure. LDZ >= N.
157 *> \endverbatim
158 *>
159 *> \param[out] V
160 *> \verbatim
161 *> V is COMPLEX*16 array, dimension (LDV,NSHFTS/2)
162 *> \endverbatim
163 *>
164 *> \param[in] LDV
165 *> \verbatim
166 *> LDV is INTEGER
167 *> LDV is the leading dimension of V as declared in the
168 *> calling procedure. LDV >= 3.
169 *> \endverbatim
170 *>
171 *> \param[out] U
172 *> \verbatim
173 *> U is COMPLEX*16 array, dimension (LDU,3*NSHFTS-3)
174 *> \endverbatim
175 *>
176 *> \param[in] LDU
177 *> \verbatim
178 *> LDU is INTEGER
179 *> LDU is the leading dimension of U just as declared in the
180 *> in the calling subroutine. LDU >= 3*NSHFTS-3.
181 *> \endverbatim
182 *>
183 *> \param[in] NV
184 *> \verbatim
185 *> NV is INTEGER
186 *> NV is the number of rows in WV agailable for workspace.
187 *> NV >= 1.
188 *> \endverbatim
189 *>
190 *> \param[out] WV
191 *> \verbatim
192 *> WV is COMPLEX*16 array, dimension (LDWV,3*NSHFTS-3)
193 *> \endverbatim
194 *>
195 *> \param[in] LDWV
196 *> \verbatim
197 *> LDWV is INTEGER
198 *> LDWV is the leading dimension of WV as declared in the
199 *> in the calling subroutine. LDWV >= NV.
200 *> \endverbatim
201 *
202 *> \param[in] NH
203 *> \verbatim
204 *> NH is INTEGER
205 *> NH is the number of columns in array WH available for
206 *> workspace. NH >= 1.
207 *> \endverbatim
208 *>
209 *> \param[out] WH
210 *> \verbatim
211 *> WH is COMPLEX*16 array, dimension (LDWH,NH)
212 *> \endverbatim
213 *>
214 *> \param[in] LDWH
215 *> \verbatim
216 *> LDWH is INTEGER
217 *> Leading dimension of WH just as declared in the
218 *> calling procedure. LDWH >= 3*NSHFTS-3.
219 *> \endverbatim
220 *>
221 * Authors:
222 * ========
223 *
224 *> \author Univ. of Tennessee
225 *> \author Univ. of California Berkeley
226 *> \author Univ. of Colorado Denver
227 *> \author NAG Ltd.
228 *
229 *> \date June 2016
230 *
231 *> \ingroup complex16OTHERauxiliary
232 *
233 *> \par Contributors:
234 * ==================
235 *>
236 *> Karen Braman and Ralph Byers, Department of Mathematics,
237 *> University of Kansas, USA
238 *
239 *> \par References:
240 * ================
241 *>
242 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
243 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
244 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
245 *> 929--947, 2002.
246 *>
247 * =====================================================================
248  SUBROUTINE zlaqr5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S,
249  $ H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV,
250  $ WV, LDWV, NH, WH, LDWH )
251 *
252 * -- LAPACK auxiliary routine (version 3.7.1) --
253 * -- LAPACK is a software package provided by Univ. of Tennessee, --
254 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
255 * June 2016
256 *
257 * .. Scalar Arguments ..
258  INTEGER IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
259  $ LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
260  LOGICAL WANTT, WANTZ
261 * ..
262 * .. Array Arguments ..
263  COMPLEX*16 H( LDH, * ), S( * ), U( LDU, * ), V( LDV, * ),
264  $ WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
265 * ..
266 *
267 * ================================================================
268 * .. Parameters ..
269  COMPLEX*16 ZERO, ONE
270  PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ),
271  $ one = ( 1.0d0, 0.0d0 ) )
272  DOUBLE PRECISION RZERO, RONE
273  PARAMETER ( RZERO = 0.0d0, rone = 1.0d0 )
274 * ..
275 * .. Local Scalars ..
276  COMPLEX*16 ALPHA, BETA, CDUM, REFSUM
277  DOUBLE PRECISION H11, H12, H21, H22, SAFMAX, SAFMIN, SCL,
278  $ smlnum, tst1, tst2, ulp
279  INTEGER I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN,
280  $ JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
281  $ m, m22, mbot, mend, mstart, mtop, nbmps, ndcol,
282  $ ns, nu
283  LOGICAL ACCUM, BLK22, BMP22
284 * ..
285 * .. External Functions ..
286  DOUBLE PRECISION DLAMCH
287  EXTERNAL DLAMCH
288 * ..
289 * .. Intrinsic Functions ..
290 *
291  INTRINSIC abs, dble, dconjg, dimag, max, min, mod
292 * ..
293 * .. Local Arrays ..
294  COMPLEX*16 VT( 3 )
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL dlabad, zgemm, zlacpy, zlaqr1, zlarfg, zlaset,
298  $ ztrmm
299 * ..
300 * .. Statement Functions ..
301  DOUBLE PRECISION CABS1
302 * ..
303 * .. Statement Function definitions ..
304  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
305 * ..
306 * .. Executable Statements ..
307 *
308 * ==== If there are no shifts, then there is nothing to do. ====
309 *
310  IF( nshfts.LT.2 )
311  $ RETURN
312 *
313 * ==== If the active block is empty or 1-by-1, then there
314 * . is nothing to do. ====
315 *
316  IF( ktop.GE.kbot )
317  $ RETURN
318 *
319 * ==== NSHFTS is supposed to be even, but if it is odd,
320 * . then simply reduce it by one. ====
321 *
322  ns = nshfts - mod( nshfts, 2 )
323 *
324 * ==== Machine constants for deflation ====
325 *
326  safmin = dlamch( 'SAFE MINIMUM' )
327  safmax = rone / safmin
328  CALL dlabad( safmin, safmax )
329  ulp = dlamch( 'PRECISION' )
330  smlnum = safmin*( dble( n ) / ulp )
331 *
332 * ==== Use accumulated reflections to update far-from-diagonal
333 * . entries ? ====
334 *
335  accum = ( kacc22.EQ.1 ) .OR. ( kacc22.EQ.2 )
336 *
337 * ==== If so, exploit the 2-by-2 block structure? ====
338 *
339  blk22 = ( ns.GT.2 ) .AND. ( kacc22.EQ.2 )
340 *
341 * ==== clear trash ====
342 *
343  IF( ktop+2.LE.kbot )
344  $ h( ktop+2, ktop ) = zero
345 *
346 * ==== NBMPS = number of 2-shift bulges in the chain ====
347 *
348  nbmps = ns / 2
349 *
350 * ==== KDU = width of slab ====
351 *
352  kdu = 6*nbmps - 3
353 *
354 * ==== Create and chase chains of NBMPS bulges ====
355 *
356  DO 210 incol = 3*( 1-nbmps ) + ktop - 1, kbot - 2, 3*nbmps - 2
357  ndcol = incol + kdu
358  IF( accum )
359  $ CALL zlaset( 'ALL', kdu, kdu, zero, one, u, ldu )
360 *
361 * ==== Near-the-diagonal bulge chase. The following loop
362 * . performs the near-the-diagonal part of a small bulge
363 * . multi-shift QR sweep. Each 6*NBMPS-2 column diagonal
364 * . chunk extends from column INCOL to column NDCOL
365 * . (including both column INCOL and column NDCOL). The
366 * . following loop chases a 3*NBMPS column long chain of
367 * . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL
368 * . may be less than KTOP and and NDCOL may be greater than
369 * . KBOT indicating phantom columns from which to chase
370 * . bulges before they are actually introduced or to which
371 * . to chase bulges beyond column KBOT.) ====
372 *
373  DO 140 krcol = incol, min( incol+3*nbmps-3, kbot-2 )
374 *
375 * ==== Bulges number MTOP to MBOT are active double implicit
376 * . shift bulges. There may or may not also be small
377 * . 2-by-2 bulge, if there is room. The inactive bulges
378 * . (if any) must wait until the active bulges have moved
379 * . down the diagonal to make room. The phantom matrix
380 * . paradigm described above helps keep track. ====
381 *
382  mtop = max( 1, ( ( ktop-1 )-krcol+2 ) / 3+1 )
383  mbot = min( nbmps, ( kbot-krcol ) / 3 )
384  m22 = mbot + 1
385  bmp22 = ( mbot.LT.nbmps ) .AND. ( krcol+3*( m22-1 ) ).EQ.
386  $ ( kbot-2 )
387 *
388 * ==== Generate reflections to chase the chain right
389 * . one column. (The minimum value of K is KTOP-1.) ====
390 *
391  DO 10 m = mtop, mbot
392  k = krcol + 3*( m-1 )
393  IF( k.EQ.ktop-1 ) THEN
394  CALL zlaqr1( 3, h( ktop, ktop ), ldh, s( 2*m-1 ),
395  $ s( 2*m ), v( 1, m ) )
396  alpha = v( 1, m )
397  CALL zlarfg( 3, alpha, v( 2, m ), 1, v( 1, m ) )
398  ELSE
399  beta = h( k+1, k )
400  v( 2, m ) = h( k+2, k )
401  v( 3, m ) = h( k+3, k )
402  CALL zlarfg( 3, beta, v( 2, m ), 1, v( 1, m ) )
403 *
404 * ==== A Bulge may collapse because of vigilant
405 * . deflation or destructive underflow. In the
406 * . underflow case, try the two-small-subdiagonals
407 * . trick to try to reinflate the bulge. ====
408 *
409  IF( h( k+3, k ).NE.zero .OR. h( k+3, k+1 ).NE.
410  $ zero .OR. h( k+3, k+2 ).EQ.zero ) THEN
411 *
412 * ==== Typical case: not collapsed (yet). ====
413 *
414  h( k+1, k ) = beta
415  h( k+2, k ) = zero
416  h( k+3, k ) = zero
417  ELSE
418 *
419 * ==== Atypical case: collapsed. Attempt to
420 * . reintroduce ignoring H(K+1,K) and H(K+2,K).
421 * . If the fill resulting from the new
422 * . reflector is too large, then abandon it.
423 * . Otherwise, use the new one. ====
424 *
425  CALL zlaqr1( 3, h( k+1, k+1 ), ldh, s( 2*m-1 ),
426  $ s( 2*m ), vt )
427  alpha = vt( 1 )
428  CALL zlarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
429  refsum = dconjg( vt( 1 ) )*
430  $ ( h( k+1, k )+dconjg( vt( 2 ) )*
431  $ h( k+2, k ) )
432 *
433  IF( cabs1( h( k+2, k )-refsum*vt( 2 ) )+
434  $ cabs1( refsum*vt( 3 ) ).GT.ulp*
435  $ ( cabs1( h( k, k ) )+cabs1( h( k+1,
436  $ k+1 ) )+cabs1( h( k+2, k+2 ) ) ) ) THEN
437 *
438 * ==== Starting a new bulge here would
439 * . create non-negligible fill. Use
440 * . the old one with trepidation. ====
441 *
442  h( k+1, k ) = beta
443  h( k+2, k ) = zero
444  h( k+3, k ) = zero
445  ELSE
446 *
447 * ==== Stating a new bulge here would
448 * . create only negligible fill.
449 * . Replace the old reflector with
450 * . the new one. ====
451 *
452  h( k+1, k ) = h( k+1, k ) - refsum
453  h( k+2, k ) = zero
454  h( k+3, k ) = zero
455  v( 1, m ) = vt( 1 )
456  v( 2, m ) = vt( 2 )
457  v( 3, m ) = vt( 3 )
458  END IF
459  END IF
460  END IF
461  10 CONTINUE
462 *
463 * ==== Generate a 2-by-2 reflection, if needed. ====
464 *
465  k = krcol + 3*( m22-1 )
466  IF( bmp22 ) THEN
467  IF( k.EQ.ktop-1 ) THEN
468  CALL zlaqr1( 2, h( k+1, k+1 ), ldh, s( 2*m22-1 ),
469  $ s( 2*m22 ), v( 1, m22 ) )
470  beta = v( 1, m22 )
471  CALL zlarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
472  ELSE
473  beta = h( k+1, k )
474  v( 2, m22 ) = h( k+2, k )
475  CALL zlarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
476  h( k+1, k ) = beta
477  h( k+2, k ) = zero
478  END IF
479  END IF
480 *
481 * ==== Multiply H by reflections from the left ====
482 *
483  IF( accum ) THEN
484  jbot = min( ndcol, kbot )
485  ELSE IF( wantt ) THEN
486  jbot = n
487  ELSE
488  jbot = kbot
489  END IF
490  DO 30 j = max( ktop, krcol ), jbot
491  mend = min( mbot, ( j-krcol+2 ) / 3 )
492  DO 20 m = mtop, mend
493  k = krcol + 3*( m-1 )
494  refsum = dconjg( v( 1, m ) )*
495  $ ( h( k+1, j )+dconjg( v( 2, m ) )*
496  $ h( k+2, j )+dconjg( v( 3, m ) )*h( k+3, j ) )
497  h( k+1, j ) = h( k+1, j ) - refsum
498  h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m )
499  h( k+3, j ) = h( k+3, j ) - refsum*v( 3, m )
500  20 CONTINUE
501  30 CONTINUE
502  IF( bmp22 ) THEN
503  k = krcol + 3*( m22-1 )
504  DO 40 j = max( k+1, ktop ), jbot
505  refsum = dconjg( v( 1, m22 ) )*
506  $ ( h( k+1, j )+dconjg( v( 2, m22 ) )*
507  $ h( k+2, j ) )
508  h( k+1, j ) = h( k+1, j ) - refsum
509  h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m22 )
510  40 CONTINUE
511  END IF
512 *
513 * ==== Multiply H by reflections from the right.
514 * . Delay filling in the last row until the
515 * . vigilant deflation check is complete. ====
516 *
517  IF( accum ) THEN
518  jtop = max( ktop, incol )
519  ELSE IF( wantt ) THEN
520  jtop = 1
521  ELSE
522  jtop = ktop
523  END IF
524  DO 80 m = mtop, mbot
525  IF( v( 1, m ).NE.zero ) THEN
526  k = krcol + 3*( m-1 )
527  DO 50 j = jtop, min( kbot, k+3 )
528  refsum = v( 1, m )*( h( j, k+1 )+v( 2, m )*
529  $ h( j, k+2 )+v( 3, m )*h( j, k+3 ) )
530  h( j, k+1 ) = h( j, k+1 ) - refsum
531  h( j, k+2 ) = h( j, k+2 ) -
532  $ refsum*dconjg( v( 2, m ) )
533  h( j, k+3 ) = h( j, k+3 ) -
534  $ refsum*dconjg( v( 3, m ) )
535  50 CONTINUE
536 *
537  IF( accum ) THEN
538 *
539 * ==== Accumulate U. (If necessary, update Z later
540 * . with with an efficient matrix-matrix
541 * . multiply.) ====
542 *
543  kms = k - incol
544  DO 60 j = max( 1, ktop-incol ), kdu
545  refsum = v( 1, m )*( u( j, kms+1 )+v( 2, m )*
546  $ u( j, kms+2 )+v( 3, m )*u( j, kms+3 ) )
547  u( j, kms+1 ) = u( j, kms+1 ) - refsum
548  u( j, kms+2 ) = u( j, kms+2 ) -
549  $ refsum*dconjg( v( 2, m ) )
550  u( j, kms+3 ) = u( j, kms+3 ) -
551  $ refsum*dconjg( v( 3, m ) )
552  60 CONTINUE
553  ELSE IF( wantz ) THEN
554 *
555 * ==== U is not accumulated, so update Z
556 * . now by multiplying by reflections
557 * . from the right. ====
558 *
559  DO 70 j = iloz, ihiz
560  refsum = v( 1, m )*( z( j, k+1 )+v( 2, m )*
561  $ z( j, k+2 )+v( 3, m )*z( j, k+3 ) )
562  z( j, k+1 ) = z( j, k+1 ) - refsum
563  z( j, k+2 ) = z( j, k+2 ) -
564  $ refsum*dconjg( v( 2, m ) )
565  z( j, k+3 ) = z( j, k+3 ) -
566  $ refsum*dconjg( v( 3, m ) )
567  70 CONTINUE
568  END IF
569  END IF
570  80 CONTINUE
571 *
572 * ==== Special case: 2-by-2 reflection (if needed) ====
573 *
574  k = krcol + 3*( m22-1 )
575  IF( bmp22 ) THEN
576  IF ( v( 1, m22 ).NE.zero ) THEN
577  DO 90 j = jtop, min( kbot, k+3 )
578  refsum = v( 1, m22 )*( h( j, k+1 )+v( 2, m22 )*
579  $ h( j, k+2 ) )
580  h( j, k+1 ) = h( j, k+1 ) - refsum
581  h( j, k+2 ) = h( j, k+2 ) -
582  $ refsum*dconjg( v( 2, m22 ) )
583  90 CONTINUE
584 *
585  IF( accum ) THEN
586  kms = k - incol
587  DO 100 j = max( 1, ktop-incol ), kdu
588  refsum = v( 1, m22 )*( u( j, kms+1 )+
589  $ v( 2, m22 )*u( j, kms+2 ) )
590  u( j, kms+1 ) = u( j, kms+1 ) - refsum
591  u( j, kms+2 ) = u( j, kms+2 ) -
592  $ refsum*dconjg( v( 2, m22 ) )
593  100 CONTINUE
594  ELSE IF( wantz ) THEN
595  DO 110 j = iloz, ihiz
596  refsum = v( 1, m22 )*( z( j, k+1 )+v( 2, m22 )*
597  $ z( j, k+2 ) )
598  z( j, k+1 ) = z( j, k+1 ) - refsum
599  z( j, k+2 ) = z( j, k+2 ) -
600  $ refsum*dconjg( v( 2, m22 ) )
601  110 CONTINUE
602  END IF
603  END IF
604  END IF
605 *
606 * ==== Vigilant deflation check ====
607 *
608  mstart = mtop
609  IF( krcol+3*( mstart-1 ).LT.ktop )
610  $ mstart = mstart + 1
611  mend = mbot
612  IF( bmp22 )
613  $ mend = mend + 1
614  IF( krcol.EQ.kbot-2 )
615  $ mend = mend + 1
616  DO 120 m = mstart, mend
617  k = min( kbot-1, krcol+3*( m-1 ) )
618 *
619 * ==== The following convergence test requires that
620 * . the tradition small-compared-to-nearby-diagonals
621 * . criterion and the Ahues & Tisseur (LAWN 122, 1997)
622 * . criteria both be satisfied. The latter improves
623 * . accuracy in some examples. Falling back on an
624 * . alternate convergence criterion when TST1 or TST2
625 * . is zero (as done here) is traditional but probably
626 * . unnecessary. ====
627 *
628  IF( h( k+1, k ).NE.zero ) THEN
629  tst1 = cabs1( h( k, k ) ) + cabs1( h( k+1, k+1 ) )
630  IF( tst1.EQ.rzero ) THEN
631  IF( k.GE.ktop+1 )
632  $ tst1 = tst1 + cabs1( h( k, k-1 ) )
633  IF( k.GE.ktop+2 )
634  $ tst1 = tst1 + cabs1( h( k, k-2 ) )
635  IF( k.GE.ktop+3 )
636  $ tst1 = tst1 + cabs1( h( k, k-3 ) )
637  IF( k.LE.kbot-2 )
638  $ tst1 = tst1 + cabs1( h( k+2, k+1 ) )
639  IF( k.LE.kbot-3 )
640  $ tst1 = tst1 + cabs1( h( k+3, k+1 ) )
641  IF( k.LE.kbot-4 )
642  $ tst1 = tst1 + cabs1( h( k+4, k+1 ) )
643  END IF
644  IF( cabs1( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
645  $ THEN
646  h12 = max( cabs1( h( k+1, k ) ),
647  $ cabs1( h( k, k+1 ) ) )
648  h21 = min( cabs1( h( k+1, k ) ),
649  $ cabs1( h( k, k+1 ) ) )
650  h11 = max( cabs1( h( k+1, k+1 ) ),
651  $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
652  h22 = min( cabs1( h( k+1, k+1 ) ),
653  $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
654  scl = h11 + h12
655  tst2 = h22*( h11 / scl )
656 *
657  IF( tst2.EQ.rzero .OR. h21*( h12 / scl ).LE.
658  $ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
659  END IF
660  END IF
661  120 CONTINUE
662 *
663 * ==== Fill in the last row of each bulge. ====
664 *
665  mend = min( nbmps, ( kbot-krcol-1 ) / 3 )
666  DO 130 m = mtop, mend
667  k = krcol + 3*( m-1 )
668  refsum = v( 1, m )*v( 3, m )*h( k+4, k+3 )
669  h( k+4, k+1 ) = -refsum
670  h( k+4, k+2 ) = -refsum*dconjg( v( 2, m ) )
671  h( k+4, k+3 ) = h( k+4, k+3 ) -
672  $ refsum*dconjg( v( 3, m ) )
673  130 CONTINUE
674 *
675 * ==== End of near-the-diagonal bulge chase. ====
676 *
677  140 CONTINUE
678 *
679 * ==== Use U (if accumulated) to update far-from-diagonal
680 * . entries in H. If required, use U to update Z as
681 * . well. ====
682 *
683  IF( accum ) THEN
684  IF( wantt ) THEN
685  jtop = 1
686  jbot = n
687  ELSE
688  jtop = ktop
689  jbot = kbot
690  END IF
691  IF( ( .NOT.blk22 ) .OR. ( incol.LT.ktop ) .OR.
692  $ ( ndcol.GT.kbot ) .OR. ( ns.LE.2 ) ) THEN
693 *
694 * ==== Updates not exploiting the 2-by-2 block
695 * . structure of U. K1 and NU keep track of
696 * . the location and size of U in the special
697 * . cases of introducing bulges and chasing
698 * . bulges off the bottom. In these special
699 * . cases and in case the number of shifts
700 * . is NS = 2, there is no 2-by-2 block
701 * . structure to exploit. ====
702 *
703  k1 = max( 1, ktop-incol )
704  nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
705 *
706 * ==== Horizontal Multiply ====
707 *
708  DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
709  jlen = min( nh, jbot-jcol+1 )
710  CALL zgemm( 'C', 'N', nu, jlen, nu, one, u( k1, k1 ),
711  $ ldu, h( incol+k1, jcol ), ldh, zero, wh,
712  $ ldwh )
713  CALL zlacpy( 'ALL', nu, jlen, wh, ldwh,
714  $ h( incol+k1, jcol ), ldh )
715  150 CONTINUE
716 *
717 * ==== Vertical multiply ====
718 *
719  DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
720  jlen = min( nv, max( ktop, incol )-jrow )
721  CALL zgemm( 'N', 'N', jlen, nu, nu, one,
722  $ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
723  $ ldu, zero, wv, ldwv )
724  CALL zlacpy( 'ALL', jlen, nu, wv, ldwv,
725  $ h( jrow, incol+k1 ), ldh )
726  160 CONTINUE
727 *
728 * ==== Z multiply (also vertical) ====
729 *
730  IF( wantz ) THEN
731  DO 170 jrow = iloz, ihiz, nv
732  jlen = min( nv, ihiz-jrow+1 )
733  CALL zgemm( 'N', 'N', jlen, nu, nu, one,
734  $ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
735  $ ldu, zero, wv, ldwv )
736  CALL zlacpy( 'ALL', jlen, nu, wv, ldwv,
737  $ z( jrow, incol+k1 ), ldz )
738  170 CONTINUE
739  END IF
740  ELSE
741 *
742 * ==== Updates exploiting U's 2-by-2 block structure.
743 * . (I2, I4, J2, J4 are the last rows and columns
744 * . of the blocks.) ====
745 *
746  i2 = ( kdu+1 ) / 2
747  i4 = kdu
748  j2 = i4 - i2
749  j4 = kdu
750 *
751 * ==== KZS and KNZ deal with the band of zeros
752 * . along the diagonal of one of the triangular
753 * . blocks. ====
754 *
755  kzs = ( j4-j2 ) - ( ns+1 )
756  knz = ns + 1
757 *
758 * ==== Horizontal multiply ====
759 *
760  DO 180 jcol = min( ndcol, kbot ) + 1, jbot, nh
761  jlen = min( nh, jbot-jcol+1 )
762 *
763 * ==== Copy bottom of H to top+KZS of scratch ====
764 * (The first KZS rows get multiplied by zero.) ====
765 *
766  CALL zlacpy( 'ALL', knz, jlen, h( incol+1+j2, jcol ),
767  $ ldh, wh( kzs+1, 1 ), ldwh )
768 *
769 * ==== Multiply by U21**H ====
770 *
771  CALL zlaset( 'ALL', kzs, jlen, zero, zero, wh, ldwh )
772  CALL ztrmm( 'L', 'U', 'C', 'N', knz, jlen, one,
773  $ u( j2+1, 1+kzs ), ldu, wh( kzs+1, 1 ),
774  $ ldwh )
775 *
776 * ==== Multiply top of H by U11**H ====
777 *
778  CALL zgemm( 'C', 'N', i2, jlen, j2, one, u, ldu,
779  $ h( incol+1, jcol ), ldh, one, wh, ldwh )
780 *
781 * ==== Copy top of H to bottom of WH ====
782 *
783  CALL zlacpy( 'ALL', j2, jlen, h( incol+1, jcol ), ldh,
784  $ wh( i2+1, 1 ), ldwh )
785 *
786 * ==== Multiply by U21**H ====
787 *
788  CALL ztrmm( 'L', 'L', 'C', 'N', j2, jlen, one,
789  $ u( 1, i2+1 ), ldu, wh( i2+1, 1 ), ldwh )
790 *
791 * ==== Multiply by U22 ====
792 *
793  CALL zgemm( 'C', 'N', i4-i2, jlen, j4-j2, one,
794  $ u( j2+1, i2+1 ), ldu,
795  $ h( incol+1+j2, jcol ), ldh, one,
796  $ wh( i2+1, 1 ), ldwh )
797 *
798 * ==== Copy it back ====
799 *
800  CALL zlacpy( 'ALL', kdu, jlen, wh, ldwh,
801  $ h( incol+1, jcol ), ldh )
802  180 CONTINUE
803 *
804 * ==== Vertical multiply ====
805 *
806  DO 190 jrow = jtop, max( incol, ktop ) - 1, nv
807  jlen = min( nv, max( incol, ktop )-jrow )
808 *
809 * ==== Copy right of H to scratch (the first KZS
810 * . columns get multiplied by zero) ====
811 *
812  CALL zlacpy( 'ALL', jlen, knz, h( jrow, incol+1+j2 ),
813  $ ldh, wv( 1, 1+kzs ), ldwv )
814 *
815 * ==== Multiply by U21 ====
816 *
817  CALL zlaset( 'ALL', jlen, kzs, zero, zero, wv, ldwv )
818  CALL ztrmm( 'R', 'U', 'N', 'N', jlen, knz, one,
819  $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
820  $ ldwv )
821 *
822 * ==== Multiply by U11 ====
823 *
824  CALL zgemm( 'N', 'N', jlen, i2, j2, one,
825  $ h( jrow, incol+1 ), ldh, u, ldu, one, wv,
826  $ ldwv )
827 *
828 * ==== Copy left of H to right of scratch ====
829 *
830  CALL zlacpy( 'ALL', jlen, j2, h( jrow, incol+1 ), ldh,
831  $ wv( 1, 1+i2 ), ldwv )
832 *
833 * ==== Multiply by U21 ====
834 *
835  CALL ztrmm( 'R', 'L', 'N', 'N', jlen, i4-i2, one,
836  $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ), ldwv )
837 *
838 * ==== Multiply by U22 ====
839 *
840  CALL zgemm( 'N', 'N', jlen, i4-i2, j4-j2, one,
841  $ h( jrow, incol+1+j2 ), ldh,
842  $ u( j2+1, i2+1 ), ldu, one, wv( 1, 1+i2 ),
843  $ ldwv )
844 *
845 * ==== Copy it back ====
846 *
847  CALL zlacpy( 'ALL', jlen, kdu, wv, ldwv,
848  $ h( jrow, incol+1 ), ldh )
849  190 CONTINUE
850 *
851 * ==== Multiply Z (also vertical) ====
852 *
853  IF( wantz ) THEN
854  DO 200 jrow = iloz, ihiz, nv
855  jlen = min( nv, ihiz-jrow+1 )
856 *
857 * ==== Copy right of Z to left of scratch (first
858 * . KZS columns get multiplied by zero) ====
859 *
860  CALL zlacpy( 'ALL', jlen, knz,
861  $ z( jrow, incol+1+j2 ), ldz,
862  $ wv( 1, 1+kzs ), ldwv )
863 *
864 * ==== Multiply by U12 ====
865 *
866  CALL zlaset( 'ALL', jlen, kzs, zero, zero, wv,
867  $ ldwv )
868  CALL ztrmm( 'R', 'U', 'N', 'N', jlen, knz, one,
869  $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
870  $ ldwv )
871 *
872 * ==== Multiply by U11 ====
873 *
874  CALL zgemm( 'N', 'N', jlen, i2, j2, one,
875  $ z( jrow, incol+1 ), ldz, u, ldu, one,
876  $ wv, ldwv )
877 *
878 * ==== Copy left of Z to right of scratch ====
879 *
880  CALL zlacpy( 'ALL', jlen, j2, z( jrow, incol+1 ),
881  $ ldz, wv( 1, 1+i2 ), ldwv )
882 *
883 * ==== Multiply by U21 ====
884 *
885  CALL ztrmm( 'R', 'L', 'N', 'N', jlen, i4-i2, one,
886  $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ),
887  $ ldwv )
888 *
889 * ==== Multiply by U22 ====
890 *
891  CALL zgemm( 'N', 'N', jlen, i4-i2, j4-j2, one,
892  $ z( jrow, incol+1+j2 ), ldz,
893  $ u( j2+1, i2+1 ), ldu, one,
894  $ wv( 1, 1+i2 ), ldwv )
895 *
896 * ==== Copy the result back to Z ====
897 *
898  CALL zlacpy( 'ALL', jlen, kdu, wv, ldwv,
899  $ z( jrow, incol+1 ), ldz )
900  200 CONTINUE
901  END IF
902  END IF
903  END IF
904  210 CONTINUE
905 *
906 * ==== End of ZLAQR5 ====
907 *
908  END
zlaqr5
subroutine zlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
ZLAQR5 performs a single small-bulge multi-shift QR sweep.
Definition: zlaqr5.f:251
zlaqr1
subroutine zlaqr1(N, H, LDH, S1, S2, V)
ZLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and spe...
Definition: zlaqr1.f:109
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:108
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
ztrmm
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:179
zlarfg
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76