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

Here is the call graph for this function:

Here is the caller graph for this function: