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 ====