LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clahqr.f
Go to the documentation of this file.
1*> \brief \b CLAHQR 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 CLAHQR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clahqr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clahqr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clahqr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLAHQR( 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 H( LDH, * ), W( * ), Z( LDZ, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CLAHQR 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*> CLAHQR 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 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 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 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, CLAHQR 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 lahqr
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 CLAHQR 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 clahqr( 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 H( LDH, * ), W( * ), Z( LDZ, * )
207* ..
208*
209* =========================================================
210*
211* .. Parameters ..
212 COMPLEX ZERO, ONE
213 parameter( zero = ( 0.0e0, 0.0e0 ),
214 $ one = ( 1.0e0, 0.0e0 ) )
215 REAL RZERO, RONE, HALF
216 parameter( rzero = 0.0e0, rone = 1.0e0, half = 0.5e0 )
217 REAL DAT1
218 parameter( dat1 = 3.0e0 / 4.0e0 )
219 INTEGER KEXSH
220 parameter( kexsh = 10 )
221* ..
222* .. Local Scalars ..
223 COMPLEX CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
224 $ v2, x, y
225 REAL 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 V( 2 )
232* ..
233* .. External Functions ..
234 COMPLEX CLADIV
235 REAL SLAMCH
236 EXTERNAL cladiv, slamch
237* ..
238* .. External Subroutines ..
239 EXTERNAL ccopy, clarfg, cscal
240* ..
241* .. Statement Functions ..
242 REAL CABS1
243* ..
244* .. Intrinsic Functions ..
245 INTRINSIC abs, aimag, conjg, max, min, real, sqrt
246* ..
247* .. Statement Function definitions ..
248 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( 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( aimag( 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 = conjg( sc ) / abs( sc )
285 h( i, i-1 ) = abs( h( i, i-1 ) )
286 CALL cscal( jhi-i+1, sc, h( i, i ), ldh )
287 CALL cscal( min( jhi, i+1 )-jlo+1, conjg( sc ), h( jlo, i ),
288 $ 1 )
289 IF( wantz )
290 $ CALL cscal( ihiz-iloz+1, conjg( 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 = slamch( 'SAFE MINIMUM' )
300 safmax = rone / safmin
301 ulp = slamch( 'PRECISION' )
302 smlnum = safmin*( real( nh ) / ulp )
303*
304* I1 and I2 are the indices of the first row and last column of H
305* to which transformations must be applied. If eigenvalues only are
306* being computed, I1 and I2 are set inside the main loop.
307*
308 IF( wantt ) THEN
309 i1 = 1
310 i2 = n
311 END IF
312*
313* ITMAX is the total number of QR iterations allowed.
314*
315 itmax = 30 * max( 10, nh )
316*
317* KDEFL counts the number of iterations since a deflation
318*
319 kdefl = 0
320*
321* The main loop begins here. I is the loop index and decreases from
322* IHI to ILO in steps of 1. Each iteration of the loop works
323* with the active submatrix in rows and columns L to I.
324* Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
325* H(L,L-1) is negligible so that the matrix splits.
326*
327 i = ihi
328 30 CONTINUE
329 IF( i.LT.ilo )
330 $ GO TO 150
331*
332* Perform QR iterations on rows and columns ILO to I until a
333* submatrix of order 1 splits off at the bottom because a
334* subdiagonal element has become negligible.
335*
336 l = ilo
337 DO 130 its = 0, itmax
338*
339* Look for a single small subdiagonal element.
340*
341 DO 40 k = i, l + 1, -1
342 IF( cabs1( h( k, k-1 ) ).LE.smlnum )
343 $ GO TO 50
344 tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
345 IF( tst.EQ.zero ) THEN
346 IF( k-2.GE.ilo )
347 $ tst = tst + abs( real( h( k-1, k-2 ) ) )
348 IF( k+1.LE.ihi )
349 $ tst = tst + abs( real( h( k+1, k ) ) )
350 END IF
351* ==== The following is a conservative small subdiagonal
352* . deflation criterion due to Ahues & Tisseur (LAWN 122,
353* . 1997). It has better mathematical foundation and
354* . improves accuracy in some examples. ====
355 IF( abs( real( h( k, k-1 ) ) ).LE.ulp*tst ) THEN
356 ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
357 ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
358 aa = max( cabs1( h( k, k ) ),
359 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
360 bb = min( cabs1( h( k, k ) ),
361 $ cabs1( h( k-1, k-1 )-h( k, k ) ) )
362 s = aa + ab
363 IF( ba*( ab / s ).LE.max( smlnum,
364 $ ulp*( bb*( aa / s ) ) ) )GO TO 50
365 END IF
366 40 CONTINUE
367 50 CONTINUE
368 l = k
369 IF( l.GT.ilo ) THEN
370*
371* H(L,L-1) is negligible
372*
373 h( l, l-1 ) = zero
374 END IF
375*
376* Exit from loop if a submatrix of order 1 has split off.
377*
378 IF( l.GE.i )
379 $ GO TO 140
380 kdefl = kdefl + 1
381*
382* Now the active submatrix is in rows and columns L to I. If
383* eigenvalues only are being computed, only the active submatrix
384* need be transformed.
385*
386 IF( .NOT.wantt ) THEN
387 i1 = l
388 i2 = i
389 END IF
390*
391 IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
392*
393* Exceptional shift.
394*
395 s = dat1*abs( real( h( i, i-1 ) ) )
396 t = s + h( i, i )
397 ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
398*
399* Exceptional shift.
400*
401 s = dat1*abs( real( h( l+1, l ) ) )
402 t = s + h( l, l )
403 ELSE
404*
405* Wilkinson's shift.
406*
407 t = h( i, i )
408 u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
409 s = cabs1( u )
410 IF( s.NE.rzero ) THEN
411 x = half*( h( i-1, i-1 )-t )
412 sx = cabs1( x )
413 s = max( s, cabs1( x ) )
414 y = s*sqrt( ( x / s )**2+( u / s )**2 )
415 IF( sx.GT.rzero ) THEN
416 IF( real( x / sx )*real( y )+aimag( x / sx )*
417 $ aimag( y ).LT.rzero )y = -y
418 END IF
419 t = t - u*cladiv( u, ( x+y ) )
420 END IF
421 END IF
422*
423* Look for two consecutive small subdiagonal elements.
424*
425 DO 60 m = i - 1, l + 1, -1
426*
427* Determine the effect of starting the single-shift QR
428* iteration at row M, and see if this would make H(M,M-1)
429* negligible.
430*
431 h11 = h( m, m )
432 h22 = h( m+1, m+1 )
433 h11s = h11 - t
434 h21 = real( h( m+1, m ) )
435 s = cabs1( h11s ) + abs( h21 )
436 h11s = h11s / s
437 h21 = h21 / s
438 v( 1 ) = h11s
439 v( 2 ) = h21
440 h10 = real( h( m, m-1 ) )
441 IF( abs( h10 )*abs( h21 ).LE.ulp*
442 $ ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
443 $ GO TO 70
444 60 CONTINUE
445 h11 = h( l, l )
446 h22 = h( l+1, l+1 )
447 h11s = h11 - t
448 h21 = real( h( l+1, l ) )
449 s = cabs1( h11s ) + abs( h21 )
450 h11s = h11s / s
451 h21 = h21 / s
452 v( 1 ) = h11s
453 v( 2 ) = h21
454 70 CONTINUE
455*
456* Single-shift QR step
457*
458 DO 120 k = m, i - 1
459*
460* The first iteration of this loop determines a reflection G
461* from the vector V and applies it from left and right to H,
462* thus creating a nonzero bulge below the subdiagonal.
463*
464* Each subsequent iteration determines a reflection G to
465* restore the Hessenberg form in the (K-1)th column, and thus
466* chases the bulge one step toward the bottom of the active
467* submatrix.
468*
469* V(2) is always real before the call to CLARFG, and hence
470* after the call T2 ( = T1*V(2) ) is also real.
471*
472 IF( k.GT.m )
473 $ CALL ccopy( 2, h( k, k-1 ), 1, v, 1 )
474 CALL clarfg( 2, v( 1 ), v( 2 ), 1, t1 )
475 IF( k.GT.m ) THEN
476 h( k, k-1 ) = v( 1 )
477 h( k+1, k-1 ) = zero
478 END IF
479 v2 = v( 2 )
480 t2 = real( t1*v2 )
481*
482* Apply G from the left to transform the rows of the matrix
483* in columns K to I2.
484*
485 DO 80 j = k, i2
486 sum = conjg( t1 )*h( k, j ) + t2*h( k+1, j )
487 h( k, j ) = h( k, j ) - sum
488 h( k+1, j ) = h( k+1, j ) - sum*v2
489 80 CONTINUE
490*
491* Apply G from the right to transform the columns of the
492* matrix in rows I1 to min(K+2,I).
493*
494 DO 90 j = i1, min( k+2, i )
495 sum = t1*h( j, k ) + t2*h( j, k+1 )
496 h( j, k ) = h( j, k ) - sum
497 h( j, k+1 ) = h( j, k+1 ) - sum*conjg( v2 )
498 90 CONTINUE
499*
500 IF( wantz ) THEN
501*
502* Accumulate transformations in the matrix Z
503*
504 DO 100 j = iloz, ihiz
505 sum = t1*z( j, k ) + t2*z( j, k+1 )
506 z( j, k ) = z( j, k ) - sum
507 z( j, k+1 ) = z( j, k+1 ) - sum*conjg( v2 )
508 100 CONTINUE
509 END IF
510*
511 IF( k.EQ.m .AND. m.GT.l ) THEN
512*
513* If the QR step was started at row M > L because two
514* consecutive small subdiagonals were found, then extra
515* scaling must be performed to ensure that H(M,M-1) remains
516* real.
517*
518 temp = one - t1
519 temp = temp / abs( temp )
520 h( m+1, m ) = h( m+1, m )*conjg( temp )
521 IF( m+2.LE.i )
522 $ h( m+2, m+1 ) = h( m+2, m+1 )*temp
523 DO 110 j = m, i
524 IF( j.NE.m+1 ) THEN
525 IF( i2.GT.j )
526 $ CALL cscal( i2-j, temp, h( j, j+1 ), ldh )
527 CALL cscal( j-i1, conjg( temp ), h( i1, j ), 1 )
528 IF( wantz ) THEN
529 CALL cscal( nz, conjg( temp ), z( iloz, j ), 1 )
530 END IF
531 END IF
532 110 CONTINUE
533 END IF
534 120 CONTINUE
535*
536* Ensure that H(I,I-1) is real.
537*
538 temp = h( i, i-1 )
539 IF( aimag( temp ).NE.rzero ) THEN
540 rtemp = abs( temp )
541 h( i, i-1 ) = rtemp
542 temp = temp / rtemp
543 IF( i2.GT.i )
544 $ CALL cscal( i2-i, conjg( temp ), h( i, i+1 ), ldh )
545 CALL cscal( i-i1, temp, h( i1, i ), 1 )
546 IF( wantz ) THEN
547 CALL cscal( nz, temp, z( iloz, i ), 1 )
548 END IF
549 END IF
550*
551 130 CONTINUE
552*
553* Failure to converge in remaining number of iterations
554*
555 info = i
556 RETURN
557*
558 140 CONTINUE
559*
560* H(I,I-1) is negligible: one eigenvalue has converged.
561*
562 w( i ) = h( i, i )
563* reset deflation counter
564 kdefl = 0
565*
566* return to start of the main loop with new value of I.
567*
568 i = l - 1
569 GO TO 30
570*
571 150 CONTINUE
572 RETURN
573*
574* End of CLAHQR
575*
576 END
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine clahqr(wantt, wantz, n, ilo, ihi, h, ldh, w, iloz, ihiz, z, ldz, info)
CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition clahqr.f:195
subroutine clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:106
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78