LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlahqr()

subroutine dlahqr ( logical  WANTT,
logical  WANTZ,
integer  N,
integer  ILO,
integer  IHI,
double precision, dimension( ldh, * )  H,
integer  LDH,
double precision, dimension( * )  WR,
double precision, dimension( * )  WI,
integer  ILOZ,
integer  IHIZ,
double precision, dimension( ldz, * )  Z,
integer  LDZ,
integer  INFO 
)

DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.

Download DLAHQR + dependencies [TGZ] [ZIP] [TXT]

Purpose:
    DLAHQR is an auxiliary routine called by DHSEQR to update the
    eigenvalues and Schur decomposition already computed by DHSEQR, by
    dealing with the Hessenberg submatrix in rows and columns ILO to
    IHI.
Parameters
[in]WANTT
          WANTT is LOGICAL
          = .TRUE. : the full Schur form T is required;
          = .FALSE.: only eigenvalues are required.
[in]WANTZ
          WANTZ is LOGICAL
          = .TRUE. : the matrix of Schur vectors Z is required;
          = .FALSE.: Schur vectors are not required.
[in]N
          N is INTEGER
          The order of the matrix H.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER
          It is assumed that H is already upper quasi-triangular in
          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
          ILO = 1). DLAHQR works primarily with the Hessenberg
          submatrix in rows and columns ILO to IHI, but applies
          transformations to all of H if WANTT is .TRUE..
          1 <= ILO <= max(1,IHI); IHI <= N.
[in,out]H
          H is DOUBLE PRECISION array, dimension (LDH,N)
          On entry, the upper Hessenberg matrix H.
          On exit, if INFO is zero and if WANTT is .TRUE., H is upper
          quasi-triangular in rows and columns ILO:IHI, with any
          2-by-2 diagonal blocks in standard form. If INFO is zero
          and WANTT is .FALSE., the contents of H are unspecified on
          exit.  The output state of H if INFO is nonzero is given
          below under the description of INFO.
[in]LDH
          LDH is INTEGER
          The leading dimension of the array H. LDH >= max(1,N).
[out]WR
          WR is DOUBLE PRECISION array, dimension (N)
[out]WI
          WI is DOUBLE PRECISION array, dimension (N)
          The real and imaginary parts, respectively, of the computed
          eigenvalues ILO to IHI are stored in the corresponding
          elements of WR and WI. If two eigenvalues are computed as a
          complex conjugate pair, they are stored in consecutive
          elements of WR and WI, say the i-th and (i+1)th, with
          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
          eigenvalues are stored in the same order as on the diagonal
          of the Schur form returned in H, with WR(i) = H(i,i), and, if
          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
[in]ILOZ
          ILOZ is INTEGER
[in]IHIZ
          IHIZ is INTEGER
          Specify the rows of Z to which transformations must be
          applied if WANTZ is .TRUE..
          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
[in,out]Z
          Z is DOUBLE PRECISION array, dimension (LDZ,N)
          If WANTZ is .TRUE., on entry Z must contain the current
          matrix Z of transformations accumulated by DHSEQR, and on
          exit Z has been updated; transformations are applied only to
          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
          If WANTZ is .FALSE., Z is not referenced.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z. LDZ >= max(1,N).
[out]INFO
          INFO is INTEGER
           = 0:  successful exit
           > 0:  If INFO = i, DLAHQR failed to compute all the
                  eigenvalues ILO to IHI in a total of 30 iterations
                  per eigenvalue; elements i+1:ihi of WR and WI
                  contain those eigenvalues which have been
                  successfully computed.

                  If INFO > 0 and WANTT is .FALSE., then on exit,
                  the remaining unconverged eigenvalues are the
                  eigenvalues of the upper Hessenberg matrix rows
                  and columns ILO through INFO of the final, output
                  value of H.

                  If INFO > 0 and WANTT is .TRUE., then on exit
          (*)       (initial value of H)*U  = U*(final value of H)
                  where U is an orthogonal matrix.    The final
                  value of H is upper Hessenberg and triangular in
                  rows and columns INFO+1 through IHI.

                  If INFO > 0 and WANTZ is .TRUE., then on exit
                      (final value of Z)  = (initial value of Z)*U
                  where U is the orthogonal matrix in (*)
                  (regardless of the value of WANTT.)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
     02-96 Based on modifications by
     David Day, Sandia National Laboratory, USA

     12-04 Further modifications by
     Ralph Byers, University of Kansas, USA
     This is a modified version of DLAHQR from LAPACK version 3.0.
     It is (1) more robust against overflow and underflow and
     (2) adopts the more conservative Ahues & Tisseur stopping
     criterion (LAWN 122, 1997).

Definition at line 205 of file dlahqr.f.

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 DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
219* ..
220*
221* =========================================================
222*
223* .. Parameters ..
224 DOUBLE PRECISION ZERO, ONE, TWO
225 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
226 DOUBLE PRECISION DAT1, DAT2
227 parameter( dat1 = 3.0d0 / 4.0d0, dat2 = -0.4375d0 )
228 INTEGER KEXSH
229 parameter( kexsh = 10 )
230* ..
231* .. Local Scalars ..
232 DOUBLE PRECISION 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 DOUBLE PRECISION V( 3 )
241* ..
242* .. External Functions ..
243 DOUBLE PRECISION DLAMCH
244 EXTERNAL dlamch
245* ..
246* .. External Subroutines ..
247 EXTERNAL dcopy, dlabad, dlanv2, dlarfg, drot
248* ..
249* .. Intrinsic Functions ..
250 INTRINSIC abs, dble, max, min, 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 = dlamch( 'SAFE MINIMUM' )
280 safmax = one / safmin
281 CALL dlabad( safmin, safmax )
282 ulp = dlamch( 'PRECISION' )
283 smlnum = safmin*( dble( nh ) / ulp )
284*
285* I1 and I2 are the indices of the first row and last column of H
286* to which transformations must be applied. If eigenvalues only are
287* being computed, I1 and I2 are set inside the main loop.
288*
289 IF( wantt ) THEN
290 i1 = 1
291 i2 = n
292 END IF
293*
294* ITMAX is the total number of QR iterations allowed.
295*
296 itmax = 30 * max( 10, nh )
297*
298* KDEFL counts the number of iterations since a deflation
299*
300 kdefl = 0
301*
302* The main loop begins here. I is the loop index and decreases from
303* IHI to ILO in steps of 1 or 2. Each iteration of the loop works
304* with the active submatrix in rows and columns L to I.
305* Eigenvalues I+1 to IHI have already converged. Either L = ILO or
306* H(L,L-1) is negligible so that the matrix splits.
307*
308 i = ihi
309 20 CONTINUE
310 l = ilo
311 IF( i.LT.ilo )
312 $ GO TO 160
313*
314* Perform QR iterations on rows and columns ILO to I until a
315* submatrix of order 1 or 2 splits off at the bottom because a
316* subdiagonal element has become negligible.
317*
318 DO 140 its = 0, itmax
319*
320* Look for a single small subdiagonal element.
321*
322 DO 30 k = i, l + 1, -1
323 IF( abs( h( k, k-1 ) ).LE.smlnum )
324 $ GO TO 40
325 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
326 IF( tst.EQ.zero ) THEN
327 IF( k-2.GE.ilo )
328 $ tst = tst + abs( h( k-1, k-2 ) )
329 IF( k+1.LE.ihi )
330 $ tst = tst + abs( h( k+1, k ) )
331 END IF
332* ==== The following is a conservative small subdiagonal
333* . deflation criterion due to Ahues & Tisseur (LAWN 122,
334* . 1997). It has better mathematical foundation and
335* . improves accuracy in some cases. ====
336 IF( abs( h( k, k-1 ) ).LE.ulp*tst ) THEN
337 ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
338 ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
339 aa = max( abs( h( k, k ) ),
340 $ abs( h( k-1, k-1 )-h( k, k ) ) )
341 bb = min( abs( h( k, k ) ),
342 $ abs( h( k-1, k-1 )-h( k, k ) ) )
343 s = aa + ab
344 IF( ba*( ab / s ).LE.max( smlnum,
345 $ ulp*( bb*( aa / s ) ) ) )GO TO 40
346 END IF
347 30 CONTINUE
348 40 CONTINUE
349 l = k
350 IF( l.GT.ilo ) THEN
351*
352* H(L,L-1) is negligible
353*
354 h( l, l-1 ) = zero
355 END IF
356*
357* Exit from loop if a submatrix of order 1 or 2 has split off.
358*
359 IF( l.GE.i-1 )
360 $ GO TO 150
361 kdefl = kdefl + 1
362*
363* Now the active submatrix is in rows and columns L to I. If
364* eigenvalues only are being computed, only the active submatrix
365* need be transformed.
366*
367 IF( .NOT.wantt ) THEN
368 i1 = l
369 i2 = i
370 END IF
371*
372 IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
373*
374* Exceptional shift.
375*
376 s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
377 h11 = dat1*s + h( i, i )
378 h12 = dat2*s
379 h21 = s
380 h22 = h11
381 ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
382*
383* Exceptional shift.
384*
385 s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
386 h11 = dat1*s + h( l, l )
387 h12 = dat2*s
388 h21 = s
389 h22 = h11
390 ELSE
391*
392* Prepare to use Francis' double shift
393* (i.e. 2nd degree generalized Rayleigh quotient)
394*
395 h11 = h( i-1, i-1 )
396 h21 = h( i, i-1 )
397 h12 = h( i-1, i )
398 h22 = h( i, i )
399 END IF
400 s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
401 IF( s.EQ.zero ) THEN
402 rt1r = zero
403 rt1i = zero
404 rt2r = zero
405 rt2i = zero
406 ELSE
407 h11 = h11 / s
408 h21 = h21 / s
409 h12 = h12 / s
410 h22 = h22 / s
411 tr = ( h11+h22 ) / two
412 det = ( h11-tr )*( h22-tr ) - h12*h21
413 rtdisc = sqrt( abs( det ) )
414 IF( det.GE.zero ) THEN
415*
416* ==== complex conjugate shifts ====
417*
418 rt1r = tr*s
419 rt2r = rt1r
420 rt1i = rtdisc*s
421 rt2i = -rt1i
422 ELSE
423*
424* ==== real shifts (use only one of them) ====
425*
426 rt1r = tr + rtdisc
427 rt2r = tr - rtdisc
428 IF( abs( rt1r-h22 ).LE.abs( rt2r-h22 ) ) THEN
429 rt1r = rt1r*s
430 rt2r = rt1r
431 ELSE
432 rt2r = rt2r*s
433 rt1r = rt2r
434 END IF
435 rt1i = zero
436 rt2i = zero
437 END IF
438 END IF
439*
440* Look for two consecutive small subdiagonal elements.
441*
442 DO 50 m = i - 2, l, -1
443* Determine the effect of starting the double-shift QR
444* iteration at row M, and see if this would make H(M,M-1)
445* negligible. (The following uses scaling to avoid
446* overflows and most underflows.)
447*
448 h21s = h( m+1, m )
449 s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
450 h21s = h( m+1, m ) / s
451 v( 1 ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*
452 $ ( ( h( m, m )-rt2r ) / s ) - rt1i*( rt2i / s )
453 v( 2 ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
454 v( 3 ) = h21s*h( m+2, m+1 )
455 s = abs( v( 1 ) ) + abs( v( 2 ) ) + abs( v( 3 ) )
456 v( 1 ) = v( 1 ) / s
457 v( 2 ) = v( 2 ) / s
458 v( 3 ) = v( 3 ) / s
459 IF( m.EQ.l )
460 $ GO TO 60
461 IF( abs( h( m, m-1 ) )*( abs( v( 2 ) )+abs( v( 3 ) ) ).LE.
462 $ ulp*abs( v( 1 ) )*( abs( h( m-1, m-1 ) )+abs( h( m,
463 $ m ) )+abs( h( m+1, m+1 ) ) ) )GO TO 60
464 50 CONTINUE
465 60 CONTINUE
466*
467* Double-shift QR step
468*
469 DO 130 k = m, i - 1
470*
471* The first iteration of this loop determines a reflection G
472* from the vector V and applies it from left and right to H,
473* thus creating a nonzero bulge below the subdiagonal.
474*
475* Each subsequent iteration determines a reflection G to
476* restore the Hessenberg form in the (K-1)th column, and thus
477* chases the bulge one step toward the bottom of the active
478* submatrix. NR is the order of G.
479*
480 nr = min( 3, i-k+1 )
481 IF( k.GT.m )
482 $ CALL dcopy( nr, h( k, k-1 ), 1, v, 1 )
483 CALL dlarfg( nr, v( 1 ), v( 2 ), 1, t1 )
484 IF( k.GT.m ) THEN
485 h( k, k-1 ) = v( 1 )
486 h( k+1, k-1 ) = zero
487 IF( k.LT.i-1 )
488 $ h( k+2, k-1 ) = zero
489 ELSE IF( m.GT.l ) THEN
490* ==== Use the following instead of
491* . H( K, K-1 ) = -H( K, K-1 ) to
492* . avoid a bug when v(2) and v(3)
493* . underflow. ====
494 h( k, k-1 ) = h( k, k-1 )*( one-t1 )
495 END IF
496 v2 = v( 2 )
497 t2 = t1*v2
498 IF( nr.EQ.3 ) THEN
499 v3 = v( 3 )
500 t3 = t1*v3
501*
502* Apply G from the left to transform the rows of the matrix
503* in columns K to I2.
504*
505 DO 70 j = k, i2
506 sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
507 h( k, j ) = h( k, j ) - sum*t1
508 h( k+1, j ) = h( k+1, j ) - sum*t2
509 h( k+2, j ) = h( k+2, j ) - sum*t3
510 70 CONTINUE
511*
512* Apply G from the right to transform the columns of the
513* matrix in rows I1 to min(K+3,I).
514*
515 DO 80 j = i1, min( k+3, i )
516 sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
517 h( j, k ) = h( j, k ) - sum*t1
518 h( j, k+1 ) = h( j, k+1 ) - sum*t2
519 h( j, k+2 ) = h( j, k+2 ) - sum*t3
520 80 CONTINUE
521*
522 IF( wantz ) THEN
523*
524* Accumulate transformations in the matrix Z
525*
526 DO 90 j = iloz, ihiz
527 sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
528 z( j, k ) = z( j, k ) - sum*t1
529 z( j, k+1 ) = z( j, k+1 ) - sum*t2
530 z( j, k+2 ) = z( j, k+2 ) - sum*t3
531 90 CONTINUE
532 END IF
533 ELSE IF( nr.EQ.2 ) THEN
534*
535* Apply G from the left to transform the rows of the matrix
536* in columns K to I2.
537*
538 DO 100 j = k, i2
539 sum = h( k, j ) + v2*h( k+1, j )
540 h( k, j ) = h( k, j ) - sum*t1
541 h( k+1, j ) = h( k+1, j ) - sum*t2
542 100 CONTINUE
543*
544* Apply G from the right to transform the columns of the
545* matrix in rows I1 to min(K+3,I).
546*
547 DO 110 j = i1, i
548 sum = h( j, k ) + v2*h( j, k+1 )
549 h( j, k ) = h( j, k ) - sum*t1
550 h( j, k+1 ) = h( j, k+1 ) - sum*t2
551 110 CONTINUE
552*
553 IF( wantz ) THEN
554*
555* Accumulate transformations in the matrix Z
556*
557 DO 120 j = iloz, ihiz
558 sum = z( j, k ) + v2*z( j, k+1 )
559 z( j, k ) = z( j, k ) - sum*t1
560 z( j, k+1 ) = z( j, k+1 ) - sum*t2
561 120 CONTINUE
562 END IF
563 END IF
564 130 CONTINUE
565*
566 140 CONTINUE
567*
568* Failure to converge in remaining number of iterations
569*
570 info = i
571 RETURN
572*
573 150 CONTINUE
574*
575 IF( l.EQ.i ) THEN
576*
577* H(I,I-1) is negligible: one eigenvalue has converged.
578*
579 wr( i ) = h( i, i )
580 wi( i ) = zero
581 ELSE IF( l.EQ.i-1 ) THEN
582*
583* H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
584*
585* Transform the 2-by-2 submatrix to standard Schur form,
586* and compute and store the eigenvalues.
587*
588 CALL dlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
589 $ h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
590 $ cs, sn )
591*
592 IF( wantt ) THEN
593*
594* Apply the transformation to the rest of H.
595*
596 IF( i2.GT.i )
597 $ CALL drot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
598 $ cs, sn )
599 CALL drot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
600 END IF
601 IF( wantz ) THEN
602*
603* Apply the transformation to Z.
604*
605 CALL drot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs, sn )
606 END IF
607 END IF
608* reset deflation counter
609 kdefl = 0
610*
611* return to start of the main loop with new value of I.
612*
613 i = l - 1
614 GO TO 20
615*
616 160 CONTINUE
617 RETURN
618*
619* End of DLAHQR
620*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dlanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition: dlanv2.f:127
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:106
Here is the call graph for this function:
Here is the caller graph for this function: