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).
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download ZLAQR2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqr2.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr2.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr2.f">
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 )
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
30 * .. Array Arguments ..
31 * COMPLEX*16 H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
32 * $ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
41 *> ZLAQR2 is identical to ZLAQR3 except that it avoids
42 *> recursion by calling ZLAHQR instead of ZLAQR4.
44 *> Aggressive early deflation:
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
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
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.
82 *> The order of the matrix H and (if WANTZ is .TRUE.) the
83 *> order of the unitary matrix Z.
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.
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.
106 *> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).
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.
123 *> Leading dimension of H just as declared in the calling
124 *> subroutine. N .LE. LDH
135 *> Specify the rows of Z to which transformations must be
136 *> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
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.
151 *> The leading dimension of Z just as declared in the
152 *> calling subroutine. 1 .LE. LDZ.
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.
166 *> The number of converged eigenvalues uncovered by this
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).
181 *> V is COMPLEX*16 array, dimension (LDV,NW)
182 *> An NW-by-NW work array.
187 *> LDV is integer scalar
188 *> The leading dimension of V just as declared in the
189 *> calling subroutine. NW .LE. LDV
194 *> NH is integer scalar
195 *> The number of columns of T. NH.GE.NW.
200 *> T is COMPLEX*16 array, dimension (LDT,NW)
206 *> The leading dimension of T just as declared in the
207 *> calling subroutine. NW .LE. LDT
213 *> The number of rows of work array WV available for
214 *> workspace. NV.GE.NW.
219 *> WV is COMPLEX*16 array, dimension (LDWV,NW)
225 *> The leading dimension of W just as declared in the
226 *> calling subroutine. NW .LE. LDV
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.
239 *> The dimension of the work array WORK. LWORK = 2*NW
240 *> suffices, but greater efficiency may result from larger
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.
253 *> \author Univ. of Tennessee
254 *> \author Univ. of California Berkeley
255 *> \author Univ. of Colorado Denver
258 *> \date September 2012
260 *> \ingroup complex16OTHERauxiliary
262 *> \par Contributors:
265 *> Karen Braman and Ralph Byers, Department of Mathematics,
266 *> University of Kansas, USA
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 )
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..--
278 * .. Scalar Arguments ..
279 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
280 $ LDZ, LWORK, N, ND, NH, NS, NV, NW
283 * .. Array Arguments ..
284 COMPLEX*16 H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
285 $ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
288 * ================================================================
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 )
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
303 * .. External Functions ..
304 DOUBLE PRECISION DLAMCH
307 * .. External Subroutines ..
308 EXTERNAL DLABAD, ZCOPY, ZGEHRD, ZGEMM, ZLACPY, ZLAHQR,
309 $ ZLARF, ZLARFG, ZLASET, ZTREXC, ZUNMHR
311 * .. Intrinsic Functions ..
312 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, MAX, MIN
314 * .. Statement Functions ..
315 DOUBLE PRECISION CABS1
317 * .. Statement Function definitions ..
318 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
320 * .. Executable Statements ..
322 * ==== Estimate optimal workspace. ====
324 JW = MIN( NW, KBOT-KTOP+1 )
329 * ==== Workspace query call to ZGEHRD ====
331 CALL ZGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
332 LWK1 = INT( WORK( 1 ) )
334 * ==== Workspace query call to ZUNMHR ====
336 CALL ZUNMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV,
338 LWK2 = INT( WORK( 1 ) )
340 * ==== Optimal workspace ====
342 LWKOPT = JW + MAX( LWK1, LWK2 )
345 * ==== Quick return in case of workspace query. ====
347 IF( LWORK.EQ.-1 ) THEN
348 WORK( 1 ) = DCMPLX( LWKOPT, 0 )
352 * ==== Nothing to do ...
353 * ... for an empty active block ... ====
359 * ... nor for an empty deflation window. ====
363 * ==== Machine constants ====
365 SAFMIN = DLAMCH( 'SAFE MINIMUM' )
366 SAFMAX = RONE / SAFMIN
367 CALL DLABAD( SAFMIN, SAFMAX )
368 ULP = DLAMCH( 'PRECISION' )
369 SMLNUM = SAFMIN*( DBLE( N ) / ULP )
371 * ==== Setup deflation window ====
373 JW = MIN( NW, KBOT-KTOP+1 )
374 KWTOP = KBOT - JW + 1
375 IF( KWTOP.EQ.KTOP ) THEN
378 S = H( KWTOP, KWTOP-1 )
381 IF( KBOT.EQ.KWTOP ) THEN
383 * ==== 1-by-1 deflation window: not much to do ====
385 SH( KWTOP ) = H( KWTOP, KWTOP )
388 IF( CABS1( S ).LE.MAX( SMLNUM, ULP*CABS1( H( KWTOP,
393 $ H( KWTOP, KWTOP-1 ) = ZERO
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.) ====
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 )
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 )
412 * ==== Deflation detection loop ====
416 DO 10 KNT = INFQR + 1, JW
418 * ==== Small spike tip deflation test ====
420 FOO = CABS1( T( NS, NS ) )
423 IF( CABS1( S )*CABS1( V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) )
426 * ==== One more converged eigenvalue ====
431 * ==== One undeflatable eigenvalue. Move it up out of the
432 * . way. (ZTREXC can not fail in this case.) ====
435 CALL ZTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO )
440 * ==== Return to Hessenberg form ====
447 * ==== sorting the diagonal of T improves accuracy for
448 * . graded matrices. ====
450 DO 30 I = INFQR + 1, NS
453 IF( CABS1( T( J, J ) ).GT.CABS1( T( IFST, IFST ) ) )
458 $ CALL ZTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO )
462 * ==== Restore shift/eigenvalue array from T ====
464 DO 40 I = INFQR + 1, JW
465 SH( KWTOP+I-1 ) = T( I, I )
469 IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
470 IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
472 * ==== Reflect spike back into lower triangle ====
474 CALL ZCOPY( NS, V, LDV, WORK, 1 )
476 WORK( I ) = DCONJG( WORK( I ) )
479 CALL ZLARFG( NS, BETA, WORK( 2 ), 1, TAU )
482 CALL ZLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )
484 CALL ZLARF( 'L', NS, JW, WORK, 1, DCONJG( TAU ), T, LDT,
486 CALL ZLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT,
488 CALL ZLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV,
491 CALL ZGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
495 * ==== Copy updated reduced window into place ====
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 ),
503 * ==== Accumulate orthogonal matrix in order update
504 * . H and Z, if requested. ====
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 )
510 * ==== Update vertical slab in H ====
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 )
524 * ==== Update horizontal slab in H ====
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 ),
536 * ==== Update vertical slab in Z ====
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 ),
549 * ==== Return the number of deflations ... ====
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
561 * ==== Return optimal workspace. ====
563 WORK( 1 ) = DCMPLX( LWKOPT, 0 )
565 * ==== End of ZLAQR2 ====