LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
dlahqr.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dlahqr (WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
 DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. More...
 

Function/Subroutine Documentation

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
          .GT. 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 .GT. 0 and WANTT is .FALSE., then on exit,
                  the remaining unconverged eigenvalues are the
                  eigenvalues of the upper Hessenberg matrix rows
                  and columns ILO thorugh INFO of the final, output
                  value of H.

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

                  If INFO .GT. 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.
Date
September 2012
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 209 of file dlahqr.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: