LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slahqr.f
Go to the documentation of this file.
1*> \brief \b SLAHQR 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 SLAHQR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slahqr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slahqr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slahqr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
22* ILOZ, 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* REAL H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SLAHQR is an auxiliary routine called by SHSEQR to update the
39*> eigenvalues and Schur decomposition already computed by SHSEQR, 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 quasi-triangular in
76*> rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
77*> ILO = 1). SLAHQR works primarily with the Hessenberg
78*> submatrix in rows and columns ILO to IHI, but applies
79*> transformations to all of 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 REAL array, dimension (LDH,N)
86*> On entry, the upper Hessenberg matrix H.
87*> On exit, if INFO is zero and if WANTT is .TRUE., H is upper
88*> quasi-triangular in rows and columns ILO:IHI, with any
89*> 2-by-2 diagonal blocks in standard form. If INFO is zero
90*> and WANTT is .FALSE., the contents of H are unspecified on
91*> exit. The output state of H if INFO is nonzero is given
92*> below under the description of INFO.
93*> \endverbatim
94*>
95*> \param[in] LDH
96*> \verbatim
97*> LDH is INTEGER
98*> The leading dimension of the array H. LDH >= max(1,N).
99*> \endverbatim
100*>
101*> \param[out] WR
102*> \verbatim
103*> WR is REAL array, dimension (N)
104*> \endverbatim
105*>
106*> \param[out] WI
107*> \verbatim
108*> WI is REAL array, dimension (N)
109*> The real and imaginary parts, respectively, of the computed
110*> eigenvalues ILO to IHI are stored in the corresponding
111*> elements of WR and WI. If two eigenvalues are computed as a
112*> complex conjugate pair, they are stored in consecutive
113*> elements of WR and WI, say the i-th and (i+1)th, with
114*> WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
115*> eigenvalues are stored in the same order as on the diagonal
116*> of the Schur form returned in H, with WR(i) = H(i,i), and, if
117*> H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
118*> WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
119*> \endverbatim
120*>
121*> \param[in] ILOZ
122*> \verbatim
123*> ILOZ is INTEGER
124*> \endverbatim
125*>
126*> \param[in] IHIZ
127*> \verbatim
128*> IHIZ is INTEGER
129*> Specify the rows of Z to which transformations must be
130*> applied if WANTZ is .TRUE..
131*> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
132*> \endverbatim
133*>
134*> \param[in,out] Z
135*> \verbatim
136*> Z is REAL array, dimension (LDZ,N)
137*> If WANTZ is .TRUE., on entry Z must contain the current
138*> matrix Z of transformations accumulated by SHSEQR, and on
139*> exit Z has been updated; transformations are applied only to
140*> the submatrix Z(ILOZ:IHIZ,ILO:IHI).
141*> If WANTZ is .FALSE., Z is not referenced.
142*> \endverbatim
143*>
144*> \param[in] LDZ
145*> \verbatim
146*> LDZ is INTEGER
147*> The leading dimension of the array Z. LDZ >= max(1,N).
148*> \endverbatim
149*>
150*> \param[out] INFO
151*> \verbatim
152*> INFO is INTEGER
153*> = 0: successful exit
154*> > 0: If INFO = i, SLAHQR failed to compute all the
155*> eigenvalues ILO to IHI in a total of 30 iterations
156*> per eigenvalue; elements i+1:ihi of WR and WI
157*> contain those eigenvalues which have been
158*> successfully computed.
159*>
160*> If INFO > 0 and WANTT is .FALSE., then on exit,
161*> the remaining unconverged eigenvalues are the
162*> eigenvalues of the upper Hessenberg matrix rows
163*> and columns ILO through INFO of the final, output
164*> value of H.
165*>
166*> If INFO > 0 and WANTT is .TRUE., then on exit
167*> (*) (initial value of H)*U = U*(final value of H)
168*> where U is an orthogonal matrix. The final
169*> value of H is upper Hessenberg and triangular in
170*> rows and columns INFO+1 through IHI.
171*>
172*> If INFO > 0 and WANTZ is .TRUE., then on exit
173*> (final value of Z) = (initial value of Z)*U
174*> where U is the orthogonal matrix in (*)
175*> (regardless of the value of WANTT.)
176*> \endverbatim
177*
178* Authors:
179* ========
180*
181*> \author Univ. of Tennessee
182*> \author Univ. of California Berkeley
183*> \author Univ. of Colorado Denver
184*> \author NAG Ltd.
185*
186*> \ingroup lahqr
187*
188*> \par Further Details:
189* =====================
190*>
191*> \verbatim
192*>
193*> 02-96 Based on modifications by
194*> David Day, Sandia National Laboratory, USA
195*>
196*> 12-04 Further modifications by
197*> Ralph Byers, University of Kansas, USA
198*> This is a modified version of SLAHQR from LAPACK version 3.0.
199*> It is (1) more robust against overflow and underflow and
200*> (2) adopts the more conservative Ahues & Tisseur stopping
201*> criterion (LAWN 122, 1997).
202*> \endverbatim
203*>
204* =====================================================================
205 SUBROUTINE slahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
206 $ ILOZ, IHIZ, Z, LDZ, INFO )
207 IMPLICIT NONE
208*
209* -- LAPACK auxiliary routine --
210* -- LAPACK is a software package provided by Univ. of Tennessee, --
211* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
212*
213* .. Scalar Arguments ..
214 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
215 LOGICAL WANTT, WANTZ
216* ..
217* .. Array Arguments ..
218 REAL H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
219* ..
220*
221* =========================================================
222*
223* .. Parameters ..
224 REAL ZERO, ONE, TWO
225 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
226 REAL DAT1, DAT2
227 parameter( dat1 = 3.0e0 / 4.0e0, dat2 = -0.4375e0 )
228 INTEGER KEXSH
229 parameter( kexsh = 10 )
230* ..
231* .. Local Scalars ..
232 REAL AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
233 $ h22, rt1i, rt1r, rt2i, rt2r, rtdisc, s, safmax,
234 $ safmin, smlnum, sn, sum, t1, t2, t3, tr, tst,
235 $ ulp, v2, v3
236 INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
237 $ kdefl
238* ..
239* .. Local Arrays ..
240 REAL V( 3 )
241* ..
242* .. External Functions ..
243 REAL SLAMCH
244 EXTERNAL slamch
245* ..
246* .. External Subroutines ..
247 EXTERNAL scopy, slanv2, slarfg, srot
248* ..
249* .. Intrinsic Functions ..
250 INTRINSIC abs, max, min, real, sqrt
251* ..
252* .. Executable Statements ..
253*
254 info = 0
255*
256* Quick return if possible
257*
258 IF( n.EQ.0 )
259 $ RETURN
260 IF( ilo.EQ.ihi ) THEN
261 wr( ilo ) = h( ilo, ilo )
262 wi( ilo ) = zero
263 RETURN
264 END IF
265*
266* ==== clear out the trash ====
267 DO 10 j = ilo, ihi - 3
268 h( j+2, j ) = zero
269 h( j+3, j ) = zero
270 10 CONTINUE
271 IF( ilo.LE.ihi-2 )
272 $ h( ihi, ihi-2 ) = zero
273*
274 nh = ihi - ilo + 1
275 nz = ihiz - iloz + 1
276*
277* Set machine-dependent constants for the stopping criterion.
278*
279 safmin = slamch( 'SAFE MINIMUM' )
280 safmax = one / safmin
281 ulp = slamch( 'PRECISION' )
282 smlnum = safmin*( real( nh ) / ulp )
283*
284* I1 and I2 are the indices of the first row and last column of H
285* to which transformations must be applied. If eigenvalues only are
286* being computed, I1 and I2 are set inside the main loop.
287*
288 IF( wantt ) THEN
289 i1 = 1
290 i2 = n
291 END IF
292*
293* ITMAX is the total number of QR iterations allowed.
294*
295 itmax = 30 * max( 10, nh )
296*
297* KDEFL counts the number of iterations since a deflation
298*
299 kdefl = 0
300*
301* The main loop begins here. I is the loop index and decreases from
302* IHI to ILO in steps of 1 or 2. Each iteration of the loop works
303* with the active submatrix in rows and columns L to I.
304* Eigenvalues I+1 to IHI have already converged. Either L = ILO or
305* H(L,L-1) is negligible so that the matrix splits.
306*
307 i = ihi
308 20 CONTINUE
309 l = ilo
310 IF( i.LT.ilo )
311 $ GO TO 160
312*
313* Perform QR iterations on rows and columns ILO to I until a
314* submatrix of order 1 or 2 splits off at the bottom because a
315* subdiagonal element has become negligible.
316*
317 DO 140 its = 0, itmax
318*
319* Look for a single small subdiagonal element.
320*
321 DO 30 k = i, l + 1, -1
322 IF( abs( h( k, k-1 ) ).LE.smlnum )
323 $ GO TO 40
324 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
325 IF( tst.EQ.zero ) THEN
326 IF( k-2.GE.ilo )
327 $ tst = tst + abs( h( k-1, k-2 ) )
328 IF( k+1.LE.ihi )
329 $ tst = tst + abs( h( k+1, k ) )
330 END IF
331* ==== The following is a conservative small subdiagonal
332* . deflation criterion due to Ahues & Tisseur (LAWN 122,
333* . 1997). It has better mathematical foundation and
334* . improves accuracy in some cases. ====
335 IF( abs( h( k, k-1 ) ).LE.ulp*tst ) THEN
336 ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
337 ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
338 aa = max( abs( h( k, k ) ),
339 $ abs( h( k-1, k-1 )-h( k, k ) ) )
340 bb = min( abs( h( k, k ) ),
341 $ abs( h( k-1, k-1 )-h( k, k ) ) )
342 s = aa + ab
343 IF( ba*( ab / s ).LE.max( smlnum,
344 $ ulp*( bb*( aa / s ) ) ) )GO TO 40
345 END IF
346 30 CONTINUE
347 40 CONTINUE
348 l = k
349 IF( l.GT.ilo ) THEN
350*
351* H(L,L-1) is negligible
352*
353 h( l, l-1 ) = zero
354 END IF
355*
356* Exit from loop if a submatrix of order 1 or 2 has split off.
357*
358 IF( l.GE.i-1 )
359 $ GO TO 150
360 kdefl = kdefl + 1
361*
362* Now the active submatrix is in rows and columns L to I. If
363* eigenvalues only are being computed, only the active submatrix
364* need be transformed.
365*
366 IF( .NOT.wantt ) THEN
367 i1 = l
368 i2 = i
369 END IF
370*
371 IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
372*
373* Exceptional shift.
374*
375 s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
376 h11 = dat1*s + h( i, i )
377 h12 = dat2*s
378 h21 = s
379 h22 = h11
380 ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
381*
382* Exceptional shift.
383*
384 s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
385 h11 = dat1*s + h( l, l )
386 h12 = dat2*s
387 h21 = s
388 h22 = h11
389 ELSE
390*
391* Prepare to use Francis' double shift
392* (i.e. 2nd degree generalized Rayleigh quotient)
393*
394 h11 = h( i-1, i-1 )
395 h21 = h( i, i-1 )
396 h12 = h( i-1, i )
397 h22 = h( i, i )
398 END IF
399 s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
400 IF( s.EQ.zero ) THEN
401 rt1r = zero
402 rt1i = zero
403 rt2r = zero
404 rt2i = zero
405 ELSE
406 h11 = h11 / s
407 h21 = h21 / s
408 h12 = h12 / s
409 h22 = h22 / s
410 tr = ( h11+h22 ) / two
411 det = ( h11-tr )*( h22-tr ) - h12*h21
412 rtdisc = sqrt( abs( det ) )
413 IF( det.GE.zero ) THEN
414*
415* ==== complex conjugate shifts ====
416*
417 rt1r = tr*s
418 rt2r = rt1r
419 rt1i = rtdisc*s
420 rt2i = -rt1i
421 ELSE
422*
423* ==== real shifts (use only one of them) ====
424*
425 rt1r = tr + rtdisc
426 rt2r = tr - rtdisc
427 IF( abs( rt1r-h22 ).LE.abs( rt2r-h22 ) ) THEN
428 rt1r = rt1r*s
429 rt2r = rt1r
430 ELSE
431 rt2r = rt2r*s
432 rt1r = rt2r
433 END IF
434 rt1i = zero
435 rt2i = zero
436 END IF
437 END IF
438*
439* Look for two consecutive small subdiagonal elements.
440*
441 DO 50 m = i - 2, l, -1
442* Determine the effect of starting the double-shift QR
443* iteration at row M, and see if this would make H(M,M-1)
444* negligible. (The following uses scaling to avoid
445* overflows and most underflows.)
446*
447 h21s = h( m+1, m )
448 s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
449 h21s = h( m+1, m ) / s
450 v( 1 ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*
451 $ ( ( h( m, m )-rt2r ) / s ) - rt1i*( rt2i / s )
452 v( 2 ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
453 v( 3 ) = h21s*h( m+2, m+1 )
454 s = abs( v( 1 ) ) + abs( v( 2 ) ) + abs( v( 3 ) )
455 v( 1 ) = v( 1 ) / s
456 v( 2 ) = v( 2 ) / s
457 v( 3 ) = v( 3 ) / s
458 IF( m.EQ.l )
459 $ GO TO 60
460 IF( abs( h( m, m-1 ) )*( abs( v( 2 ) )+abs( v( 3 ) ) ).LE.
461 $ ulp*abs( v( 1 ) )*( abs( h( m-1, m-1 ) )+abs( h( m,
462 $ m ) )+abs( h( m+1, m+1 ) ) ) )GO TO 60
463 50 CONTINUE
464 60 CONTINUE
465*
466* Double-shift QR step
467*
468 DO 130 k = m, i - 1
469*
470* The first iteration of this loop determines a reflection G
471* from the vector V and applies it from left and right to H,
472* thus creating a nonzero bulge below the subdiagonal.
473*
474* Each subsequent iteration determines a reflection G to
475* restore the Hessenberg form in the (K-1)th column, and thus
476* chases the bulge one step toward the bottom of the active
477* submatrix. NR is the order of G.
478*
479 nr = min( 3, i-k+1 )
480 IF( k.GT.m )
481 $ CALL scopy( nr, h( k, k-1 ), 1, v, 1 )
482 CALL slarfg( nr, v( 1 ), v( 2 ), 1, t1 )
483 IF( k.GT.m ) THEN
484 h( k, k-1 ) = v( 1 )
485 h( k+1, k-1 ) = zero
486 IF( k.LT.i-1 )
487 $ h( k+2, k-1 ) = zero
488 ELSE IF( m.GT.l ) THEN
489* ==== Use the following instead of
490* . H( K, K-1 ) = -H( K, K-1 ) to
491* . avoid a bug when v(2) and v(3)
492* . underflow. ====
493 h( k, k-1 ) = h( k, k-1 )*( one-t1 )
494 END IF
495 v2 = v( 2 )
496 t2 = t1*v2
497 IF( nr.EQ.3 ) THEN
498 v3 = v( 3 )
499 t3 = t1*v3
500*
501* Apply G from the left to transform the rows of the matrix
502* in columns K to I2.
503*
504 DO 70 j = k, i2
505 sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
506 h( k, j ) = h( k, j ) - sum*t1
507 h( k+1, j ) = h( k+1, j ) - sum*t2
508 h( k+2, j ) = h( k+2, j ) - sum*t3
509 70 CONTINUE
510*
511* Apply G from the right to transform the columns of the
512* matrix in rows I1 to min(K+3,I).
513*
514 DO 80 j = i1, min( k+3, i )
515 sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
516 h( j, k ) = h( j, k ) - sum*t1
517 h( j, k+1 ) = h( j, k+1 ) - sum*t2
518 h( j, k+2 ) = h( j, k+2 ) - sum*t3
519 80 CONTINUE
520*
521 IF( wantz ) THEN
522*
523* Accumulate transformations in the matrix Z
524*
525 DO 90 j = iloz, ihiz
526 sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
527 z( j, k ) = z( j, k ) - sum*t1
528 z( j, k+1 ) = z( j, k+1 ) - sum*t2
529 z( j, k+2 ) = z( j, k+2 ) - sum*t3
530 90 CONTINUE
531 END IF
532 ELSE IF( nr.EQ.2 ) THEN
533*
534* Apply G from the left to transform the rows of the matrix
535* in columns K to I2.
536*
537 DO 100 j = k, i2
538 sum = h( k, j ) + v2*h( k+1, j )
539 h( k, j ) = h( k, j ) - sum*t1
540 h( k+1, j ) = h( k+1, j ) - sum*t2
541 100 CONTINUE
542*
543* Apply G from the right to transform the columns of the
544* matrix in rows I1 to min(K+3,I).
545*
546 DO 110 j = i1, i
547 sum = h( j, k ) + v2*h( j, k+1 )
548 h( j, k ) = h( j, k ) - sum*t1
549 h( j, k+1 ) = h( j, k+1 ) - sum*t2
550 110 CONTINUE
551*
552 IF( wantz ) THEN
553*
554* Accumulate transformations in the matrix Z
555*
556 DO 120 j = iloz, ihiz
557 sum = z( j, k ) + v2*z( j, k+1 )
558 z( j, k ) = z( j, k ) - sum*t1
559 z( j, k+1 ) = z( j, k+1 ) - sum*t2
560 120 CONTINUE
561 END IF
562 END IF
563 130 CONTINUE
564*
565 140 CONTINUE
566*
567* Failure to converge in remaining number of iterations
568*
569 info = i
570 RETURN
571*
572 150 CONTINUE
573*
574 IF( l.EQ.i ) THEN
575*
576* H(I,I-1) is negligible: one eigenvalue has converged.
577*
578 wr( i ) = h( i, i )
579 wi( i ) = zero
580 ELSE IF( l.EQ.i-1 ) THEN
581*
582* H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
583*
584* Transform the 2-by-2 submatrix to standard Schur form,
585* and compute and store the eigenvalues.
586*
587 CALL slanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
588 $ h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
589 $ cs, sn )
590*
591 IF( wantt ) THEN
592*
593* Apply the transformation to the rest of H.
594*
595 IF( i2.GT.i )
596 $ CALL srot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
597 $ cs, sn )
598 CALL srot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
599 END IF
600 IF( wantz ) THEN
601*
602* Apply the transformation to Z.
603*
604 CALL srot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs, sn )
605 END IF
606 END IF
607* reset deflation counter
608 kdefl = 0
609*
610* return to start of the main loop with new value of I.
611*
612 i = l - 1
613 GO TO 20
614*
615 160 CONTINUE
616 RETURN
617*
618* End of SLAHQR
619*
620 END
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine slahqr(wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, iloz, ihiz, z, ldz, info)
SLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition slahqr.f:207
subroutine slanv2(a, b, c, d, rt1r, rt1i, rt2r, rt2i, cs, sn)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition slanv2.f:127
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92