LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zlaqr2.f
Go to the documentation of this file.
1 *> \brief \b ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZLAQR2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqr2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
22 * IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT,
23 * NV, WV, LDWV, WORK, LWORK )
24 *
25 * .. Scalar Arguments ..
26 * INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
27 * $ LDZ, LWORK, N, ND, NH, NS, NV, NW
28 * LOGICAL WANTT, WANTZ
29 * ..
30 * .. Array Arguments ..
31 * COMPLEX*16 H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
32 * $ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> ZLAQR2 is identical to ZLAQR3 except that it avoids
42 *> recursion by calling ZLAHQR instead of ZLAQR4.
43 *>
44 *> Aggressive early deflation:
45 *>
46 *> ZLAQR2 accepts as input an upper Hessenberg matrix
47 *> H and performs an unitary similarity transformation
48 *> designed to detect and deflate fully converged eigenvalues from
49 *> a trailing principal submatrix. On output H has been over-
50 *> written by a new Hessenberg matrix that is a perturbation of
51 *> an unitary similarity transformation of H. It is to be
52 *> hoped that the final version of H has many zero subdiagonal
53 *> entries.
54 *>
55 *> \endverbatim
56 *
57 * Arguments:
58 * ==========
59 *
60 *> \param[in] WANTT
61 *> \verbatim
62 *> WANTT is LOGICAL
63 *> If .TRUE., then the Hessenberg matrix H is fully updated
64 *> so that the triangular Schur factor may be
65 *> computed (in cooperation with the calling subroutine).
66 *> If .FALSE., then only enough of H is updated to preserve
67 *> the eigenvalues.
68 *> \endverbatim
69 *>
70 *> \param[in] WANTZ
71 *> \verbatim
72 *> WANTZ is LOGICAL
73 *> If .TRUE., then the unitary matrix Z is updated so
74 *> so that the unitary Schur factor may be computed
75 *> (in cooperation with the calling subroutine).
76 *> If .FALSE., then Z is not referenced.
77 *> \endverbatim
78 *>
79 *> \param[in] N
80 *> \verbatim
81 *> N is INTEGER
82 *> The order of the matrix H and (if WANTZ is .TRUE.) the
83 *> order of the unitary matrix Z.
84 *> \endverbatim
85 *>
86 *> \param[in] KTOP
87 *> \verbatim
88 *> KTOP is INTEGER
89 *> It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
90 *> KBOT and KTOP together determine an isolated block
91 *> along the diagonal of the Hessenberg matrix.
92 *> \endverbatim
93 *>
94 *> \param[in] KBOT
95 *> \verbatim
96 *> KBOT is INTEGER
97 *> It is assumed without a check that either
98 *> KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together
99 *> determine an isolated block along the diagonal of the
100 *> Hessenberg matrix.
101 *> \endverbatim
102 *>
103 *> \param[in] NW
104 *> \verbatim
105 *> NW is INTEGER
106 *> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).
107 *> \endverbatim
108 *>
109 *> \param[in,out] H
110 *> \verbatim
111 *> H is COMPLEX*16 array, dimension (LDH,N)
112 *> On input the initial N-by-N section of H stores the
113 *> Hessenberg matrix undergoing aggressive early deflation.
114 *> On output H has been transformed by a unitary
115 *> similarity transformation, perturbed, and the returned
116 *> to Hessenberg form that (it is to be hoped) has some
117 *> zero subdiagonal entries.
118 *> \endverbatim
119 *>
120 *> \param[in] LDH
121 *> \verbatim
122 *> LDH is integer
123 *> Leading dimension of H just as declared in the calling
124 *> subroutine. N .LE. LDH
125 *> \endverbatim
126 *>
127 *> \param[in] ILOZ
128 *> \verbatim
129 *> ILOZ is INTEGER
130 *> \endverbatim
131 *>
132 *> \param[in] IHIZ
133 *> \verbatim
134 *> IHIZ is INTEGER
135 *> Specify the rows of Z to which transformations must be
136 *> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
137 *> \endverbatim
138 *>
139 *> \param[in,out] Z
140 *> \verbatim
141 *> Z is COMPLEX*16 array, dimension (LDZ,N)
142 *> IF WANTZ is .TRUE., then on output, the unitary
143 *> similarity transformation mentioned above has been
144 *> accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right.
145 *> If WANTZ is .FALSE., then Z is unreferenced.
146 *> \endverbatim
147 *>
148 *> \param[in] LDZ
149 *> \verbatim
150 *> LDZ is integer
151 *> The leading dimension of Z just as declared in the
152 *> calling subroutine. 1 .LE. LDZ.
153 *> \endverbatim
154 *>
155 *> \param[out] NS
156 *> \verbatim
157 *> NS is integer
158 *> The number of unconverged (ie approximate) eigenvalues
159 *> returned in SR and SI that may be used as shifts by the
160 *> calling subroutine.
161 *> \endverbatim
162 *>
163 *> \param[out] ND
164 *> \verbatim
165 *> ND is integer
166 *> The number of converged eigenvalues uncovered by this
167 *> subroutine.
168 *> \endverbatim
169 *>
170 *> \param[out] SH
171 *> \verbatim
172 *> SH is COMPLEX*16 array, dimension KBOT
173 *> On output, approximate eigenvalues that may
174 *> be used for shifts are stored in SH(KBOT-ND-NS+1)
175 *> through SR(KBOT-ND). Converged eigenvalues are
176 *> stored in SH(KBOT-ND+1) through SH(KBOT).
177 *> \endverbatim
178 *>
179 *> \param[out] V
180 *> \verbatim
181 *> V is COMPLEX*16 array, dimension (LDV,NW)
182 *> An NW-by-NW work array.
183 *> \endverbatim
184 *>
185 *> \param[in] LDV
186 *> \verbatim
187 *> LDV is integer scalar
188 *> The leading dimension of V just as declared in the
189 *> calling subroutine. NW .LE. LDV
190 *> \endverbatim
191 *>
192 *> \param[in] NH
193 *> \verbatim
194 *> NH is integer scalar
195 *> The number of columns of T. NH.GE.NW.
196 *> \endverbatim
197 *>
198 *> \param[out] T
199 *> \verbatim
200 *> T is COMPLEX*16 array, dimension (LDT,NW)
201 *> \endverbatim
202 *>
203 *> \param[in] LDT
204 *> \verbatim
205 *> LDT is integer
206 *> The leading dimension of T just as declared in the
207 *> calling subroutine. NW .LE. LDT
208 *> \endverbatim
209 *>
210 *> \param[in] NV
211 *> \verbatim
212 *> NV is integer
213 *> The number of rows of work array WV available for
214 *> workspace. NV.GE.NW.
215 *> \endverbatim
216 *>
217 *> \param[out] WV
218 *> \verbatim
219 *> WV is COMPLEX*16 array, dimension (LDWV,NW)
220 *> \endverbatim
221 *>
222 *> \param[in] LDWV
223 *> \verbatim
224 *> LDWV is integer
225 *> The leading dimension of W just as declared in the
226 *> calling subroutine. NW .LE. LDV
227 *> \endverbatim
228 *>
229 *> \param[out] WORK
230 *> \verbatim
231 *> WORK is COMPLEX*16 array, dimension LWORK.
232 *> On exit, WORK(1) is set to an estimate of the optimal value
233 *> of LWORK for the given values of N, NW, KTOP and KBOT.
234 *> \endverbatim
235 *>
236 *> \param[in] LWORK
237 *> \verbatim
238 *> LWORK is integer
239 *> The dimension of the work array WORK. LWORK = 2*NW
240 *> suffices, but greater efficiency may result from larger
241 *> values of LWORK.
242 *>
243 *> If LWORK = -1, then a workspace query is assumed; ZLAQR2
244 *> only estimates the optimal workspace size for the given
245 *> values of N, NW, KTOP and KBOT. The estimate is returned
246 *> in WORK(1). No error message related to LWORK is issued
247 *> by XERBLA. Neither H nor Z are accessed.
248 *> \endverbatim
249 *
250 * Authors:
251 * ========
252 *
253 *> \author Univ. of Tennessee
254 *> \author Univ. of California Berkeley
255 *> \author Univ. of Colorado Denver
256 *> \author NAG Ltd.
257 *
258 *> \date September 2012
259 *
260 *> \ingroup complex16OTHERauxiliary
261 *
262 *> \par Contributors:
263 * ==================
264 *>
265 *> Karen Braman and Ralph Byers, Department of Mathematics,
266 *> University of Kansas, USA
267 *>
268 * =====================================================================
269  SUBROUTINE zlaqr2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
270  $ ihiz, z, ldz, ns, nd, sh, v, ldv, nh, t, ldt,
271  $ nv, wv, ldwv, work, lwork )
272 *
273 * -- LAPACK auxiliary routine (version 3.4.2) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * September 2012
277 *
278 * .. Scalar Arguments ..
279  INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
280  $ ldz, lwork, n, nd, nh, ns, nv, nw
281  LOGICAL WANTT, WANTZ
282 * ..
283 * .. Array Arguments ..
284  COMPLEX*16 H( ldh, * ), SH( * ), T( ldt, * ), V( ldv, * ),
285  $ work( * ), wv( ldwv, * ), z( ldz, * )
286 * ..
287 *
288 * ================================================================
289 *
290 * .. Parameters ..
291  COMPLEX*16 ZERO, ONE
292  parameter ( zero = ( 0.0d0, 0.0d0 ),
293  $ one = ( 1.0d0, 0.0d0 ) )
294  DOUBLE PRECISION RZERO, RONE
295  parameter ( rzero = 0.0d0, rone = 1.0d0 )
296 * ..
297 * .. Local Scalars ..
298  COMPLEX*16 BETA, CDUM, S, TAU
299  DOUBLE PRECISION FOO, SAFMAX, SAFMIN, SMLNUM, ULP
300  INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
301  $ knt, krow, kwtop, ltop, lwk1, lwk2, lwkopt
302 * ..
303 * .. External Functions ..
304  DOUBLE PRECISION DLAMCH
305  EXTERNAL dlamch
306 * ..
307 * .. External Subroutines ..
308  EXTERNAL dlabad, zcopy, zgehrd, zgemm, zlacpy, zlahqr,
310 * ..
311 * .. Intrinsic Functions ..
312  INTRINSIC abs, dble, dcmplx, dconjg, dimag, int, max, min
313 * ..
314 * .. Statement Functions ..
315  DOUBLE PRECISION CABS1
316 * ..
317 * .. Statement Function definitions ..
318  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
319 * ..
320 * .. Executable Statements ..
321 *
322 * ==== Estimate optimal workspace. ====
323 *
324  jw = min( nw, kbot-ktop+1 )
325  IF( jw.LE.2 ) THEN
326  lwkopt = 1
327  ELSE
328 *
329 * ==== Workspace query call to ZGEHRD ====
330 *
331  CALL zgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
332  lwk1 = int( work( 1 ) )
333 *
334 * ==== Workspace query call to ZUNMHR ====
335 *
336  CALL zunmhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
337  $ work, -1, info )
338  lwk2 = int( work( 1 ) )
339 *
340 * ==== Optimal workspace ====
341 *
342  lwkopt = jw + max( lwk1, lwk2 )
343  END IF
344 *
345 * ==== Quick return in case of workspace query. ====
346 *
347  IF( lwork.EQ.-1 ) THEN
348  work( 1 ) = dcmplx( lwkopt, 0 )
349  RETURN
350  END IF
351 *
352 * ==== Nothing to do ...
353 * ... for an empty active block ... ====
354  ns = 0
355  nd = 0
356  work( 1 ) = one
357  IF( ktop.GT.kbot )
358  $ RETURN
359 * ... nor for an empty deflation window. ====
360  IF( nw.LT.1 )
361  $ RETURN
362 *
363 * ==== Machine constants ====
364 *
365  safmin = dlamch( 'SAFE MINIMUM' )
366  safmax = rone / safmin
367  CALL dlabad( safmin, safmax )
368  ulp = dlamch( 'PRECISION' )
369  smlnum = safmin*( dble( n ) / ulp )
370 *
371 * ==== Setup deflation window ====
372 *
373  jw = min( nw, kbot-ktop+1 )
374  kwtop = kbot - jw + 1
375  IF( kwtop.EQ.ktop ) THEN
376  s = zero
377  ELSE
378  s = h( kwtop, kwtop-1 )
379  END IF
380 *
381  IF( kbot.EQ.kwtop ) THEN
382 *
383 * ==== 1-by-1 deflation window: not much to do ====
384 *
385  sh( kwtop ) = h( kwtop, kwtop )
386  ns = 1
387  nd = 0
388  IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
389  $ kwtop ) ) ) ) THEN
390  ns = 0
391  nd = 1
392  IF( kwtop.GT.ktop )
393  $ h( kwtop, kwtop-1 ) = zero
394  END IF
395  work( 1 ) = one
396  RETURN
397  END IF
398 *
399 * ==== Convert to spike-triangular form. (In case of a
400 * . rare QR failure, this routine continues to do
401 * . aggressive early deflation using that part of
402 * . the deflation window that converged using INFQR
403 * . here and there to keep track.) ====
404 *
405  CALL zlacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
406  CALL zcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
407 *
408  CALL zlaset( 'A', jw, jw, zero, one, v, ldv )
409  CALL zlahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
410  $ jw, v, ldv, infqr )
411 *
412 * ==== Deflation detection loop ====
413 *
414  ns = jw
415  ilst = infqr + 1
416  DO 10 knt = infqr + 1, jw
417 *
418 * ==== Small spike tip deflation test ====
419 *
420  foo = cabs1( t( ns, ns ) )
421  IF( foo.EQ.rzero )
422  $ foo = cabs1( s )
423  IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
424  $ THEN
425 *
426 * ==== One more converged eigenvalue ====
427 *
428  ns = ns - 1
429  ELSE
430 *
431 * ==== One undeflatable eigenvalue. Move it up out of the
432 * . way. (ZTREXC can not fail in this case.) ====
433 *
434  ifst = ns
435  CALL ztrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
436  ilst = ilst + 1
437  END IF
438  10 CONTINUE
439 *
440 * ==== Return to Hessenberg form ====
441 *
442  IF( ns.EQ.0 )
443  $ s = zero
444 *
445  IF( ns.LT.jw ) THEN
446 *
447 * ==== sorting the diagonal of T improves accuracy for
448 * . graded matrices. ====
449 *
450  DO 30 i = infqr + 1, ns
451  ifst = i
452  DO 20 j = i + 1, ns
453  IF( cabs1( t( j, j ) ).GT.cabs1( t( ifst, ifst ) ) )
454  $ ifst = j
455  20 CONTINUE
456  ilst = i
457  IF( ifst.NE.ilst )
458  $ CALL ztrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
459  30 CONTINUE
460  END IF
461 *
462 * ==== Restore shift/eigenvalue array from T ====
463 *
464  DO 40 i = infqr + 1, jw
465  sh( kwtop+i-1 ) = t( i, i )
466  40 CONTINUE
467 *
468 *
469  IF( ns.LT.jw .OR. s.EQ.zero ) THEN
470  IF( ns.GT.1 .AND. s.NE.zero ) THEN
471 *
472 * ==== Reflect spike back into lower triangle ====
473 *
474  CALL zcopy( ns, v, ldv, work, 1 )
475  DO 50 i = 1, ns
476  work( i ) = dconjg( work( i ) )
477  50 CONTINUE
478  beta = work( 1 )
479  CALL zlarfg( ns, beta, work( 2 ), 1, tau )
480  work( 1 ) = one
481 *
482  CALL zlaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
483 *
484  CALL zlarf( 'L', ns, jw, work, 1, dconjg( tau ), t, ldt,
485  $ work( jw+1 ) )
486  CALL zlarf( 'R', ns, ns, work, 1, tau, t, ldt,
487  $ work( jw+1 ) )
488  CALL zlarf( 'R', jw, ns, work, 1, tau, v, ldv,
489  $ work( jw+1 ) )
490 *
491  CALL zgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
492  $ lwork-jw, info )
493  END IF
494 *
495 * ==== Copy updated reduced window into place ====
496 *
497  IF( kwtop.GT.1 )
498  $ h( kwtop, kwtop-1 ) = s*dconjg( v( 1, 1 ) )
499  CALL zlacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
500  CALL zcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
501  $ ldh+1 )
502 *
503 * ==== Accumulate orthogonal matrix in order update
504 * . H and Z, if requested. ====
505 *
506  IF( ns.GT.1 .AND. s.NE.zero )
507  $ CALL zunmhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
508  $ work( jw+1 ), lwork-jw, info )
509 *
510 * ==== Update vertical slab in H ====
511 *
512  IF( wantt ) THEN
513  ltop = 1
514  ELSE
515  ltop = ktop
516  END IF
517  DO 60 krow = ltop, kwtop - 1, nv
518  kln = min( nv, kwtop-krow )
519  CALL zgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
520  $ ldh, v, ldv, zero, wv, ldwv )
521  CALL zlacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
522  60 CONTINUE
523 *
524 * ==== Update horizontal slab in H ====
525 *
526  IF( wantt ) THEN
527  DO 70 kcol = kbot + 1, n, nh
528  kln = min( nh, n-kcol+1 )
529  CALL zgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
530  $ h( kwtop, kcol ), ldh, zero, t, ldt )
531  CALL zlacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
532  $ ldh )
533  70 CONTINUE
534  END IF
535 *
536 * ==== Update vertical slab in Z ====
537 *
538  IF( wantz ) THEN
539  DO 80 krow = iloz, ihiz, nv
540  kln = min( nv, ihiz-krow+1 )
541  CALL zgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
542  $ ldz, v, ldv, zero, wv, ldwv )
543  CALL zlacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
544  $ ldz )
545  80 CONTINUE
546  END IF
547  END IF
548 *
549 * ==== Return the number of deflations ... ====
550 *
551  nd = jw - ns
552 *
553 * ==== ... and the number of shifts. (Subtracting
554 * . INFQR from the spike length takes care
555 * . of the case of a rare QR failure while
556 * . calculating eigenvalues of the deflation
557 * . window.) ====
558 *
559  ns = ns - infqr
560 *
561 * ==== Return optimal workspace. ====
562 *
563  work( 1 ) = dcmplx( lwkopt, 0 )
564 *
565 * ==== End of ZLAQR2 ====
566 *
567  END
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
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:169
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine zunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMHR
Definition: zunmhr.f:180
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
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine zlaqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
Definition: zlaqr2.f:272
subroutine zlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
Definition: zlahqr.f:197
subroutine ztrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
ZTREXC
Definition: ztrexc.f:126
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130