LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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: