LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zlahqr.f
Go to the documentation of this file.
1*> \brief \b ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZLAHQR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahqr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahqr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahqr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
22* IHIZ, Z, LDZ, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
26* LOGICAL WANTT, WANTZ
27* ..
28* .. Array Arguments ..
29* COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> ZLAHQR is an auxiliary routine called by CHSEQR to update the
39*> eigenvalues and Schur decomposition already computed by CHSEQR, by
40*> dealing with the Hessenberg submatrix in rows and columns ILO to
41*> IHI.
42*> \endverbatim
43*
44* Arguments:
45* ==========
46*
47*> \param[in] WANTT
48*> \verbatim
49*> WANTT is LOGICAL
50*> = .TRUE. : the full Schur form T is required;
51*> = .FALSE.: only eigenvalues are required.
52*> \endverbatim
53*>
54*> \param[in] WANTZ
55*> \verbatim
56*> WANTZ is LOGICAL
57*> = .TRUE. : the matrix of Schur vectors Z is required;
58*> = .FALSE.: Schur vectors are not required.
59*> \endverbatim
60*>
61*> \param[in] N
62*> \verbatim
63*> N is INTEGER
64*> The order of the matrix H. N >= 0.
65*> \endverbatim
66*>
67*> \param[in] ILO
68*> \verbatim
69*> ILO is INTEGER
70*> \endverbatim
71*>
72*> \param[in] IHI
73*> \verbatim
74*> IHI is INTEGER
75*> It is assumed that H is already upper triangular in rows and
76*> columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
77*> ZLAHQR works primarily with the Hessenberg submatrix in rows
78*> and columns ILO to IHI, but applies transformations to all of
79*> H if WANTT is .TRUE..
80*> 1 <= ILO <= max(1,IHI); IHI <= N.
81*> \endverbatim
82*>
83*> \param[in,out] H
84*> \verbatim
85*> H is COMPLEX*16 array, dimension (LDH,N)
86*> On entry, the upper Hessenberg matrix H.
87*> On exit, if INFO is zero and if WANTT is .TRUE., then H
88*> is upper triangular in rows and columns ILO:IHI. If INFO
89*> is zero and if WANTT is .FALSE., then the contents of H
90*> are unspecified on exit. The output state of H in case
91*> INF is positive is below under the description of INFO.
92*> \endverbatim
93*>
94*> \param[in] LDH
95*> \verbatim
96*> LDH is INTEGER
97*> The leading dimension of the array H. LDH >= max(1,N).
98*> \endverbatim
99*>
100*> \param[out] W
101*> \verbatim
102*> W is COMPLEX*16 array, dimension (N)
103*> The computed eigenvalues ILO to IHI are stored in the
104*> corresponding elements of W. If WANTT is .TRUE., the
105*> eigenvalues are stored in the same order as on the diagonal
106*> of the Schur form returned in H, with W(i) = H(i,i).
107*> \endverbatim
108*>
109*> \param[in] ILOZ
110*> \verbatim
111*> ILOZ is INTEGER
112*> \endverbatim
113*>
114*> \param[in] IHIZ
115*> \verbatim
116*> IHIZ is INTEGER
117*> Specify the rows of Z to which transformations must be
118*> applied if WANTZ is .TRUE..
119*> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
120*> \endverbatim
121*>
122*> \param[in,out] Z
123*> \verbatim
124*> Z is COMPLEX*16 array, dimension (LDZ,N)
125*> If WANTZ is .TRUE., on entry Z must contain the current
126*> matrix Z of transformations accumulated by CHSEQR, and on
127*> exit Z has been updated; transformations are applied only to
128*> the submatrix Z(ILOZ:IHIZ,ILO:IHI).
129*> If WANTZ is .FALSE., Z is not referenced.
130*> \endverbatim
131*>
132*> \param[in] LDZ
133*> \verbatim
134*> LDZ is INTEGER
135*> The leading dimension of the array Z. LDZ >= max(1,N).
136*> \endverbatim
137*>
138*> \param[out] INFO
139*> \verbatim
140*> INFO is INTEGER
141*> = 0: successful exit
142*> > 0: if INFO = i, ZLAHQR failed to compute all the
143*> eigenvalues ILO to IHI in a total of 30 iterations
144*> per eigenvalue; elements i+1:ihi of W contain
145*> those eigenvalues which have been successfully
146*> computed.
147*>
148*> If INFO > 0 and WANTT is .FALSE., then on exit,
149*> the remaining unconverged eigenvalues are the
150*> eigenvalues of the upper Hessenberg matrix
151*> rows and columns ILO through INFO of the final,
152*> output value of H.
153*>
154*> If INFO > 0 and WANTT is .TRUE., then on exit
155*> (*) (initial value of H)*U = U*(final value of H)
156*> where U is an orthogonal matrix. The final
157*> value of H is upper Hessenberg and triangular in
158*> rows and columns INFO+1 through IHI.
159*>
160*> If INFO > 0 and WANTZ is .TRUE., then on exit
161*> (final value of Z) = (initial value of Z)*U
162*> where U is the orthogonal matrix in (*)
163*> (regardless of the value of WANTT.)
164*> \endverbatim
165*
166* Authors:
167* ========
168*
169*> \author Univ. of Tennessee
170*> \author Univ. of California Berkeley
171*> \author Univ. of Colorado Denver
172*> \author NAG Ltd.
173*
174*> \ingroup complex16OTHERauxiliary
175*
176*> \par Contributors:
177* ==================
178*>
179*> \verbatim
180*>
181*> 02-96 Based on modifications by
182*> David Day, Sandia National Laboratory, USA
183*>
184*> 12-04 Further modifications by
185*> Ralph Byers, University of Kansas, USA
186*> This is a modified version of ZLAHQR from LAPACK version 3.0.
187*> It is (1) more robust against overflow and underflow and
188*> (2) adopts the more conservative Ahues & Tisseur stopping
189*> criterion (LAWN 122, 1997).
190*> \endverbatim
191*>
192* =====================================================================
193 SUBROUTINE zlahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
194 $ IHIZ, Z, LDZ, INFO )
195 IMPLICIT NONE
196*
197* -- LAPACK auxiliary routine --
198* -- LAPACK is a software package provided by Univ. of Tennessee, --
199* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
200*
201* .. Scalar Arguments ..
202 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
203 LOGICAL WANTT, WANTZ
204* ..
205* .. Array Arguments ..
206 COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * )
207* ..
208*
209* =========================================================
210*
211* .. Parameters ..
212 COMPLEX*16 ZERO, ONE
213 parameter( zero = ( 0.0d0, 0.0d0 ),
214 $ one = ( 1.0d0, 0.0d0 ) )
215 DOUBLE PRECISION RZERO, RONE, HALF
216 parameter( rzero = 0.0d0, rone = 1.0d0, half = 0.5d0 )
217 DOUBLE PRECISION DAT1
218 parameter( dat1 = 3.0d0 / 4.0d0 )
219 INTEGER KEXSH
220 parameter( kexsh = 10 )
221* ..
222* .. Local Scalars ..
223 COMPLEX*16 CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
224 $ v2, x, y
225 DOUBLE PRECISION AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
226 $ safmin, smlnum, sx, t2, tst, ulp
227 INTEGER I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
228 $ nh, nz, kdefl
229* ..
230* .. Local Arrays ..
231 COMPLEX*16 V( 2 )
232* ..
233* .. External Functions ..
234 COMPLEX*16 ZLADIV
235 DOUBLE PRECISION DLAMCH
236 EXTERNAL zladiv, dlamch
237* ..
238* .. External Subroutines ..
239 EXTERNAL dlabad, zcopy, zlarfg, zscal
240* ..
241* .. Statement Functions ..
242 DOUBLE PRECISION CABS1
243* ..
244* .. Intrinsic Functions ..
245 INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
246* ..
247* .. Statement Function definitions ..
248 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
249* ..
250* .. Executable Statements ..
251*
252 info = 0
253*
254* Quick return if possible
255*
256 IF( n.EQ.0 )
257 $ RETURN
258 IF( ilo.EQ.ihi ) THEN
259 w( ilo ) = h( ilo, ilo )
260 RETURN
261 END IF
262*
263* ==== clear out the trash ====
264 DO 10 j = ilo, ihi - 3
265 h( j+2, j ) = zero
266 h( j+3, j ) = zero
267 10 CONTINUE
268 IF( ilo.LE.ihi-2 )
269 $ h( ihi, ihi-2 ) = zero
270* ==== ensure that subdiagonal entries are real ====
271 IF( wantt ) THEN
272 jlo = 1
273 jhi = n
274 ELSE
275 jlo = ilo
276 jhi = ihi
277 END IF
278 DO 20 i = ilo + 1, ihi
279 IF( dimag( h( i, i-1 ) ).NE.rzero ) THEN
280* ==== The following redundant normalization
281* . avoids problems with both gradual and
282* . sudden underflow in ABS(H(I,I-1)) ====
283 sc = h( i, i-1 ) / cabs1( h( i, i-1 ) )
284 sc = dconjg( sc ) / abs( sc )
285 h( i, i-1 ) = abs( h( i, i-1 ) )
286 CALL zscal( jhi-i+1, sc, h( i, i ), ldh )
287 CALL zscal( min( jhi, i+1 )-jlo+1, dconjg( sc ),
288 $ h( jlo, i ), 1 )
289 IF( wantz )
290 $ CALL zscal( ihiz-iloz+1, dconjg( sc ), z( iloz, i ), 1 )
291 END IF
292 20 CONTINUE
293*
294 nh = ihi - ilo + 1
295 nz = ihiz - iloz + 1
296*
297* Set machine-dependent constants for the stopping criterion.
298*
299 safmin = dlamch( 'SAFE MINIMUM' )
300 safmax = rone / safmin
301 CALL dlabad( safmin, safmax )
302 ulp = dlamch( 'PRECISION' )
303 smlnum = safmin*( dble( nh ) / ulp )
304*
305* I1 and I2 are the indices of the first row and last column of H
306* to which transformations must be applied. If eigenvalues only are
307* being computed, I1 and I2 are set inside the main loop.
308*
309 IF( wantt ) THEN
310 i1 = 1
311 i2 = n
312 END IF
313*
314* ITMAX is the total number of QR iterations allowed.
315*
316 itmax = 30 * max( 10, nh )
317*
318* KDEFL counts the number of iterations since a deflation
319*
320 kdefl = 0
321*
322* The main loop begins here. I is the loop index and decreases from
323* IHI to ILO in steps of 1. Each iteration of the loop works
324* with the active submatrix in rows and columns L to I.
325* Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
326* H(L,L-1) is negligible so that the matrix splits.
327*
328 i = ihi
329 30 CONTINUE
330 IF( i.LT.ilo )
331 $ GO TO 150
332*
333* Perform QR iterations on rows and columns ILO to I until a
334* submatrix of order 1 splits off at the bottom because a
335* subdiagonal element has become negligible.
336*
337 l = ilo
338 DO 130 its = 0, itmax
339*
340* Look for a single small subdiagonal element.
341*
342 DO 40 k = i, l + 1, -1
343 IF( cabs1( h( k, k-1 ) ).LE.smlnum )
344 $ GO TO 50
345 tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
346 IF( tst.EQ.zero ) THEN
347 IF( k-2.GE.ilo )
348 $ tst = tst + abs( dble( h( k-1, k-2 ) ) )
349 IF( k+1.LE.ihi )
350 $ tst = tst + abs( dble( h( k+1, k ) ) )
351 END IF
352* ==== The following is a conservative small subdiagonal
353* . deflation criterion due to Ahues & Tisseur (LAWN 122,
354* . 1997). It has better mathematical foundation and
355* . improves accuracy in some examples. ====
356 IF( abs( dble( h( k, k-1 ) ) ).LE.ulp*tst ) THEN
357 ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
358 ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
359 aa = max( cabs1( h( k, k ) ),
360 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
361 bb = min( cabs1( h( k, k ) ),
362 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
363 s = aa + ab
364 IF( ba*( ab / s ).LE.max( smlnum,
365 $ ulp*( bb*( aa / s ) ) ) )GO TO 50
366 END IF
367 40 CONTINUE
368 50 CONTINUE
369 l = k
370 IF( l.GT.ilo ) THEN
371*
372* H(L,L-1) is negligible
373*
374 h( l, l-1 ) = zero
375 END IF
376*
377* Exit from loop if a submatrix of order 1 has split off.
378*
379 IF( l.GE.i )
380 $ GO TO 140
381 kdefl = kdefl + 1
382*
383* Now the active submatrix is in rows and columns L to I. If
384* eigenvalues only are being computed, only the active submatrix
385* need be transformed.
386*
387 IF( .NOT.wantt ) THEN
388 i1 = l
389 i2 = i
390 END IF
391*
392 IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
393*
394* Exceptional shift.
395*
396 s = dat1*abs( dble( h( i, i-1 ) ) )
397 t = s + h( i, i )
398 ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
399*
400* Exceptional shift.
401*
402 s = dat1*abs( dble( h( l+1, l ) ) )
403 t = s + h( l, l )
404 ELSE
405*
406* Wilkinson's shift.
407*
408 t = h( i, i )
409 u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
410 s = cabs1( u )
411 IF( s.NE.rzero ) THEN
412 x = half*( h( i-1, i-1 )-t )
413 sx = cabs1( x )
414 s = max( s, cabs1( x ) )
415 y = s*sqrt( ( x / s )**2+( u / s )**2 )
416 IF( sx.GT.rzero ) THEN
417 IF( dble( x / sx )*dble( y )+dimag( x / sx )*
418 $ dimag( y ).LT.rzero )y = -y
419 END IF
420 t = t - u*zladiv( u, ( x+y ) )
421 END IF
422 END IF
423*
424* Look for two consecutive small subdiagonal elements.
425*
426 DO 60 m = i - 1, l + 1, -1
427*
428* Determine the effect of starting the single-shift QR
429* iteration at row M, and see if this would make H(M,M-1)
430* negligible.
431*
432 h11 = h( m, m )
433 h22 = h( m+1, m+1 )
434 h11s = h11 - t
435 h21 = dble( h( m+1, m ) )
436 s = cabs1( h11s ) + abs( h21 )
437 h11s = h11s / s
438 h21 = h21 / s
439 v( 1 ) = h11s
440 v( 2 ) = h21
441 h10 = dble( h( m, m-1 ) )
442 IF( abs( h10 )*abs( h21 ).LE.ulp*
443 $ ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
444 $ GO TO 70
445 60 CONTINUE
446 h11 = h( l, l )
447 h22 = h( l+1, l+1 )
448 h11s = h11 - t
449 h21 = dble( h( l+1, l ) )
450 s = cabs1( h11s ) + abs( h21 )
451 h11s = h11s / s
452 h21 = h21 / s
453 v( 1 ) = h11s
454 v( 2 ) = h21
455 70 CONTINUE
456*
457* Single-shift QR step
458*
459 DO 120 k = m, i - 1
460*
461* The first iteration of this loop determines a reflection G
462* from the vector V and applies it from left and right to H,
463* thus creating a nonzero bulge below the subdiagonal.
464*
465* Each subsequent iteration determines a reflection G to
466* restore the Hessenberg form in the (K-1)th column, and thus
467* chases the bulge one step toward the bottom of the active
468* submatrix.
469*
470* V(2) is always real before the call to ZLARFG, and hence
471* after the call T2 ( = T1*V(2) ) is also real.
472*
473 IF( k.GT.m )
474 $ CALL zcopy( 2, h( k, k-1 ), 1, v, 1 )
475 CALL zlarfg( 2, v( 1 ), v( 2 ), 1, t1 )
476 IF( k.GT.m ) THEN
477 h( k, k-1 ) = v( 1 )
478 h( k+1, k-1 ) = zero
479 END IF
480 v2 = v( 2 )
481 t2 = dble( t1*v2 )
482*
483* Apply G from the left to transform the rows of the matrix
484* in columns K to I2.
485*
486 DO 80 j = k, i2
487 sum = dconjg( t1 )*h( k, j ) + t2*h( k+1, j )
488 h( k, j ) = h( k, j ) - sum
489 h( k+1, j ) = h( k+1, j ) - sum*v2
490 80 CONTINUE
491*
492* Apply G from the right to transform the columns of the
493* matrix in rows I1 to min(K+2,I).
494*
495 DO 90 j = i1, min( k+2, i )
496 sum = t1*h( j, k ) + t2*h( j, k+1 )
497 h( j, k ) = h( j, k ) - sum
498 h( j, k+1 ) = h( j, k+1 ) - sum*dconjg( v2 )
499 90 CONTINUE
500*
501 IF( wantz ) THEN
502*
503* Accumulate transformations in the matrix Z
504*
505 DO 100 j = iloz, ihiz
506 sum = t1*z( j, k ) + t2*z( j, k+1 )
507 z( j, k ) = z( j, k ) - sum
508 z( j, k+1 ) = z( j, k+1 ) - sum*dconjg( v2 )
509 100 CONTINUE
510 END IF
511*
512 IF( k.EQ.m .AND. m.GT.l ) THEN
513*
514* If the QR step was started at row M > L because two
515* consecutive small subdiagonals were found, then extra
516* scaling must be performed to ensure that H(M,M-1) remains
517* real.
518*
519 temp = one - t1
520 temp = temp / abs( temp )
521 h( m+1, m ) = h( m+1, m )*dconjg( temp )
522 IF( m+2.LE.i )
523 $ h( m+2, m+1 ) = h( m+2, m+1 )*temp
524 DO 110 j = m, i
525 IF( j.NE.m+1 ) THEN
526 IF( i2.GT.j )
527 $ CALL zscal( i2-j, temp, h( j, j+1 ), ldh )
528 CALL zscal( j-i1, dconjg( temp ), h( i1, j ), 1 )
529 IF( wantz ) THEN
530 CALL zscal( nz, dconjg( temp ), z( iloz, j ),
531 $ 1 )
532 END IF
533 END IF
534 110 CONTINUE
535 END IF
536 120 CONTINUE
537*
538* Ensure that H(I,I-1) is real.
539*
540 temp = h( i, i-1 )
541 IF( dimag( temp ).NE.rzero ) THEN
542 rtemp = abs( temp )
543 h( i, i-1 ) = rtemp
544 temp = temp / rtemp
545 IF( i2.GT.i )
546 $ CALL zscal( i2-i, dconjg( temp ), h( i, i+1 ), ldh )
547 CALL zscal( i-i1, temp, h( i1, i ), 1 )
548 IF( wantz ) THEN
549 CALL zscal( nz, temp, z( iloz, i ), 1 )
550 END IF
551 END IF
552*
553 130 CONTINUE
554*
555* Failure to converge in remaining number of iterations
556*
557 info = i
558 RETURN
559*
560 140 CONTINUE
561*
562* H(I,I-1) is negligible: one eigenvalue has converged.
563*
564 w( i ) = h( i, i )
565* reset deflation counter
566 kdefl = 0
567*
568* return to start of the main loop with new value of I.
569*
570 i = l - 1
571 GO TO 30
572*
573 150 CONTINUE
574 RETURN
575*
576* End of ZLAHQR
577*
578 END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
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,...
Definition: zlahqr.f:195
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:106