LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlaqr4 ( 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,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition.

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

Purpose: 
    DLAQR4 implements one level of recursion for DLAQR0.
    It is a complete implementation of the small bulge multi-shift
    QR algorithm.  It may be called by DLAQR0 and, for large enough
    deflation window size, it may be called by DLAQR3.  This
    subroutine is identical to DLAQR0 except that it calls DLAQR2
    instead of DLAQR3.

    DLAQR4 computes the eigenvalues of a Hessenberg matrix H
    and, optionally, the matrices T and Z from the Schur decomposition
    H = Z T Z**T, where T is an upper quasi-triangular matrix (the
    Schur form), and Z is the orthogonal matrix of Schur vectors.

    Optionally Z may be postmultiplied into an input orthogonal
    matrix Q so that this routine can give the Schur factorization
    of a matrix A which has been reduced to the Hessenberg form H
    by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
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 .GE. 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER
           It is assumed that H is already upper triangular in rows
           and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
           H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
           previous call to DGEBAL, and then passed to DGEHRD when the
           matrix output by DGEBAL is reduced to Hessenberg form.
           Otherwise, ILO and IHI should be set to 1 and N,
           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
           If N = 0, then ILO = 1 and IHI = 0.
[in,out]H
          H is DOUBLE PRECISION array, dimension (LDH,N)
           On entry, the upper Hessenberg matrix H.
           On exit, if INFO = 0 and WANTT is .TRUE., then H contains
           the upper quasi-triangular matrix T from the Schur
           decomposition (the Schur form); 2-by-2 diagonal blocks
           (corresponding to complex conjugate pairs of eigenvalues)
           are returned in standard form, with H(i,i) = H(i+1,i+1)
           and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is
           .FALSE., then the contents of H are unspecified on exit.
           (The output value of H when INFO.GT.0 is given under the
           description of INFO below.)

           This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
           j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
[in]LDH
          LDH is INTEGER
           The leading dimension of the array H. LDH .GE. max(1,N).
[out]WR
          WR is DOUBLE PRECISION array, dimension (IHI)
[out]WI
          WI is DOUBLE PRECISION array, dimension (IHI)
           The real and imaginary parts, respectively, of the computed
           eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
           and WI(ILO:IHI). 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) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then
           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 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
[in,out]Z
          Z is DOUBLE PRECISION array, dimension (LDZ,IHI)
           If WANTZ is .FALSE., then Z is not referenced.
           If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
           replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
           orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
           (The output value of Z when INFO.GT.0 is given under
           the description of INFO below.)
[in]LDZ
          LDZ is INTEGER
           The leading dimension of the array Z.  if WANTZ is .TRUE.
           then LDZ.GE.MAX(1,IHIZ).  Otherwize, LDZ.GE.1.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension LWORK
           On exit, if LWORK = -1, WORK(1) returns an estimate of
           the optimal value for LWORK.
[in]LWORK
          LWORK is INTEGER
           The dimension of the array WORK.  LWORK .GE. max(1,N)
           is sufficient, but LWORK typically as large as 6*N may
           be required for optimal performance.  A workspace query
           to determine the optimal workspace size is recommended.

           If LWORK = -1, then DLAQR4 does a workspace query.
           In this case, DLAQR4 checks the input parameters and
           estimates the optimal workspace size for the given
           values of N, ILO and IHI.  The estimate is returned
           in WORK(1).  No error message related to LWORK is
           issued by XERBLA.  Neither H nor Z are accessed.
[out]INFO
          INFO is INTEGER
             =  0:  successful exit
           .GT. 0:  if INFO = i, DLAQR4 failed to compute all of
                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
                and WI contain those eigenvalues which have been
                successfully computed.  (Failures are rare.)

                If INFO .GT. 0 and WANT is .FALSE., then on exit,
                the remaining unconverged eigenvalues are the eigen-
                values of the upper Hessenberg matrix rows and
                columns ILO through 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 a orthogonal 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(ILO:IHI,ILOZ:IHIZ)
                   =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U

                where U is the orthogonal matrix in (*) (regard-
                less of the value of WANTT.)

                If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
                accessed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA
References:
K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 Performance, SIAM Journal of Matrix Analysis, volume 23, pages 929–947, 2002.
K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part II: Aggressive Early Deflation, SIAM Journal of Matrix Analysis, volume 23, pages 948–973, 2002.

Definition at line 265 of file dlaqr4.f.

265 *
266 * -- LAPACK auxiliary routine (version 3.4.2) --
267 * -- LAPACK is a software package provided by Univ. of Tennessee, --
268 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269 * September 2012
270 *
271 * .. Scalar Arguments ..
272  INTEGER ihi, ihiz, ilo, iloz, info, ldh, ldz, lwork, n
273  LOGICAL wantt, wantz
274 * ..
275 * .. Array Arguments ..
276  DOUBLE PRECISION h( ldh, * ), wi( * ), work( * ), wr( * ),
277  $ z( ldz, * )
278 * ..
279 *
280 * ================================================================
281 * .. Parameters ..
282 *
283 * ==== Matrices of order NTINY or smaller must be processed by
284 * . DLAHQR because of insufficient subdiagonal scratch space.
285 * . (This is a hard limit.) ====
286  INTEGER ntiny
287  parameter ( ntiny = 11 )
288 *
289 * ==== Exceptional deflation windows: try to cure rare
290 * . slow convergence by varying the size of the
291 * . deflation window after KEXNW iterations. ====
292  INTEGER kexnw
293  parameter ( kexnw = 5 )
294 *
295 * ==== Exceptional shifts: try to cure rare slow convergence
296 * . with ad-hoc exceptional shifts every KEXSH iterations.
297 * . ====
298  INTEGER kexsh
299  parameter ( kexsh = 6 )
300 *
301 * ==== The constants WILK1 and WILK2 are used to form the
302 * . exceptional shifts. ====
303  DOUBLE PRECISION wilk1, wilk2
304  parameter ( wilk1 = 0.75d0, wilk2 = -0.4375d0 )
305  DOUBLE PRECISION zero, one
306  parameter ( zero = 0.0d0, one = 1.0d0 )
307 * ..
308 * .. Local Scalars ..
309  DOUBLE PRECISION aa, bb, cc, cs, dd, sn, ss, swap
310  INTEGER i, inf, it, itmax, k, kacc22, kbot, kdu, ks,
311  $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
312  $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
313  $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
314  LOGICAL sorted
315  CHARACTER jbcmpz*2
316 * ..
317 * .. External Functions ..
318  INTEGER ilaenv
319  EXTERNAL ilaenv
320 * ..
321 * .. Local Arrays ..
322  DOUBLE PRECISION zdum( 1, 1 )
323 * ..
324 * .. External Subroutines ..
325  EXTERNAL dlacpy, dlahqr, dlanv2, dlaqr2, dlaqr5
326 * ..
327 * .. Intrinsic Functions ..
328  INTRINSIC abs, dble, int, max, min, mod
329 * ..
330 * .. Executable Statements ..
331  info = 0
332 *
333 * ==== Quick return for N = 0: nothing to do. ====
334 *
335  IF( n.EQ.0 ) THEN
336  work( 1 ) = one
337  RETURN
338  END IF
339 *
340  IF( n.LE.ntiny ) THEN
341 *
342 * ==== Tiny matrices must use DLAHQR. ====
343 *
344  lwkopt = 1
345  IF( lwork.NE.-1 )
346  $ CALL dlahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,
347  $ iloz, ihiz, z, ldz, info )
348  ELSE
349 *
350 * ==== Use small bulge multi-shift QR with aggressive early
351 * . deflation on larger-than-tiny matrices. ====
352 *
353 * ==== Hope for the best. ====
354 *
355  info = 0
356 *
357 * ==== Set up job flags for ILAENV. ====
358 *
359  IF( wantt ) THEN
360  jbcmpz( 1: 1 ) = 'S'
361  ELSE
362  jbcmpz( 1: 1 ) = 'E'
363  END IF
364  IF( wantz ) THEN
365  jbcmpz( 2: 2 ) = 'V'
366  ELSE
367  jbcmpz( 2: 2 ) = 'N'
368  END IF
369 *
370 * ==== NWR = recommended deflation window size. At this
371 * . point, N .GT. NTINY = 11, so there is enough
372 * . subdiagonal workspace for NWR.GE.2 as required.
373 * . (In fact, there is enough subdiagonal space for
374 * . NWR.GE.3.) ====
375 *
376  nwr = ilaenv( 13, 'DLAQR4', jbcmpz, n, ilo, ihi, lwork )
377  nwr = max( 2, nwr )
378  nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
379 *
380 * ==== NSR = recommended number of simultaneous shifts.
381 * . At this point N .GT. NTINY = 11, so there is at
382 * . enough subdiagonal workspace for NSR to be even
383 * . and greater than or equal to two as required. ====
384 *
385  nsr = ilaenv( 15, 'DLAQR4', jbcmpz, n, ilo, ihi, lwork )
386  nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
387  nsr = max( 2, nsr-mod( nsr, 2 ) )
388 *
389 * ==== Estimate optimal workspace ====
390 *
391 * ==== Workspace query call to DLAQR2 ====
392 *
393  CALL dlaqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
394  $ ihiz, z, ldz, ls, ld, wr, wi, h, ldh, n, h, ldh,
395  $ n, h, ldh, work, -1 )
396 *
397 * ==== Optimal workspace = MAX(DLAQR5, DLAQR2) ====
398 *
399  lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
400 *
401 * ==== Quick return in case of workspace query. ====
402 *
403  IF( lwork.EQ.-1 ) THEN
404  work( 1 ) = dble( lwkopt )
405  RETURN
406  END IF
407 *
408 * ==== DLAHQR/DLAQR0 crossover point ====
409 *
410  nmin = ilaenv( 12, 'DLAQR4', jbcmpz, n, ilo, ihi, lwork )
411  nmin = max( ntiny, nmin )
412 *
413 * ==== Nibble crossover point ====
414 *
415  nibble = ilaenv( 14, 'DLAQR4', jbcmpz, n, ilo, ihi, lwork )
416  nibble = max( 0, nibble )
417 *
418 * ==== Accumulate reflections during ttswp? Use block
419 * . 2-by-2 structure during matrix-matrix multiply? ====
420 *
421  kacc22 = ilaenv( 16, 'DLAQR4', jbcmpz, n, ilo, ihi, lwork )
422  kacc22 = max( 0, kacc22 )
423  kacc22 = min( 2, kacc22 )
424 *
425 * ==== NWMAX = the largest possible deflation window for
426 * . which there is sufficient workspace. ====
427 *
428  nwmax = min( ( n-1 ) / 3, lwork / 2 )
429  nw = nwmax
430 *
431 * ==== NSMAX = the Largest number of simultaneous shifts
432 * . for which there is sufficient workspace. ====
433 *
434  nsmax = min( ( n+6 ) / 9, 2*lwork / 3 )
435  nsmax = nsmax - mod( nsmax, 2 )
436 *
437 * ==== NDFL: an iteration count restarted at deflation. ====
438 *
439  ndfl = 1
440 *
441 * ==== ITMAX = iteration limit ====
442 *
443  itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
444 *
445 * ==== Last row and column in the active block ====
446 *
447  kbot = ihi
448 *
449 * ==== Main Loop ====
450 *
451  DO 80 it = 1, itmax
452 *
453 * ==== Done when KBOT falls below ILO ====
454 *
455  IF( kbot.LT.ilo )
456  $ GO TO 90
457 *
458 * ==== Locate active block ====
459 *
460  DO 10 k = kbot, ilo + 1, -1
461  IF( h( k, k-1 ).EQ.zero )
462  $ GO TO 20
463  10 CONTINUE
464  k = ilo
465  20 CONTINUE
466  ktop = k
467 *
468 * ==== Select deflation window size:
469 * . Typical Case:
470 * . If possible and advisable, nibble the entire
471 * . active block. If not, use size MIN(NWR,NWMAX)
472 * . or MIN(NWR+1,NWMAX) depending upon which has
473 * . the smaller corresponding subdiagonal entry
474 * . (a heuristic).
475 * .
476 * . Exceptional Case:
477 * . If there have been no deflations in KEXNW or
478 * . more iterations, then vary the deflation window
479 * . size. At first, because, larger windows are,
480 * . in general, more powerful than smaller ones,
481 * . rapidly increase the window to the maximum possible.
482 * . Then, gradually reduce the window size. ====
483 *
484  nh = kbot - ktop + 1
485  nwupbd = min( nh, nwmax )
486  IF( ndfl.LT.kexnw ) THEN
487  nw = min( nwupbd, nwr )
488  ELSE
489  nw = min( nwupbd, 2*nw )
490  END IF
491  IF( nw.LT.nwmax ) THEN
492  IF( nw.GE.nh-1 ) THEN
493  nw = nh
494  ELSE
495  kwtop = kbot - nw + 1
496  IF( abs( h( kwtop, kwtop-1 ) ).GT.
497  $ abs( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
498  END IF
499  END IF
500  IF( ndfl.LT.kexnw ) THEN
501  ndec = -1
502  ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd ) THEN
503  ndec = ndec + 1
504  IF( nw-ndec.LT.2 )
505  $ ndec = 0
506  nw = nw - ndec
507  END IF
508 *
509 * ==== Aggressive early deflation:
510 * . split workspace under the subdiagonal into
511 * . - an nw-by-nw work array V in the lower
512 * . left-hand-corner,
513 * . - an NW-by-at-least-NW-but-more-is-better
514 * . (NW-by-NHO) horizontal work array along
515 * . the bottom edge,
516 * . - an at-least-NW-but-more-is-better (NHV-by-NW)
517 * . vertical work array along the left-hand-edge.
518 * . ====
519 *
520  kv = n - nw + 1
521  kt = nw + 1
522  nho = ( n-nw-1 ) - kt + 1
523  kwv = nw + 2
524  nve = ( n-nw ) - kwv + 1
525 *
526 * ==== Aggressive early deflation ====
527 *
528  CALL dlaqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
529  $ ihiz, z, ldz, ls, ld, wr, wi, h( kv, 1 ), ldh,
530  $ nho, h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh,
531  $ work, lwork )
532 *
533 * ==== Adjust KBOT accounting for new deflations. ====
534 *
535  kbot = kbot - ld
536 *
537 * ==== KS points to the shifts. ====
538 *
539  ks = kbot - ls + 1
540 *
541 * ==== Skip an expensive QR sweep if there is a (partly
542 * . heuristic) reason to expect that many eigenvalues
543 * . will deflate without it. Here, the QR sweep is
544 * . skipped if many eigenvalues have just been deflated
545 * . or if the remaining active block is small.
546 *
547  IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
548  $ ktop+1.GT.min( nmin, nwmax ) ) ) ) THEN
549 *
550 * ==== NS = nominal number of simultaneous shifts.
551 * . This may be lowered (slightly) if DLAQR2
552 * . did not provide that many shifts. ====
553 *
554  ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
555  ns = ns - mod( ns, 2 )
556 *
557 * ==== If there have been no deflations
558 * . in a multiple of KEXSH iterations,
559 * . then try exceptional shifts.
560 * . Otherwise use shifts provided by
561 * . DLAQR2 above or from the eigenvalues
562 * . of a trailing principal submatrix. ====
563 *
564  IF( mod( ndfl, kexsh ).EQ.0 ) THEN
565  ks = kbot - ns + 1
566  DO 30 i = kbot, max( ks+1, ktop+2 ), -2
567  ss = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
568  aa = wilk1*ss + h( i, i )
569  bb = ss
570  cc = wilk2*ss
571  dd = aa
572  CALL dlanv2( aa, bb, cc, dd, wr( i-1 ), wi( i-1 ),
573  $ wr( i ), wi( i ), cs, sn )
574  30 CONTINUE
575  IF( ks.EQ.ktop ) THEN
576  wr( ks+1 ) = h( ks+1, ks+1 )
577  wi( ks+1 ) = zero
578  wr( ks ) = wr( ks+1 )
579  wi( ks ) = wi( ks+1 )
580  END IF
581  ELSE
582 *
583 * ==== Got NS/2 or fewer shifts? Use DLAHQR
584 * . on a trailing principal submatrix to
585 * . get more. (Since NS.LE.NSMAX.LE.(N+6)/9,
586 * . there is enough space below the subdiagonal
587 * . to fit an NS-by-NS scratch array.) ====
588 *
589  IF( kbot-ks+1.LE.ns / 2 ) THEN
590  ks = kbot - ns + 1
591  kt = n - ns + 1
592  CALL dlacpy( 'A', ns, ns, h( ks, ks ), ldh,
593  $ h( kt, 1 ), ldh )
594  CALL dlahqr( .false., .false., ns, 1, ns,
595  $ h( kt, 1 ), ldh, wr( ks ), wi( ks ),
596  $ 1, 1, zdum, 1, inf )
597  ks = ks + inf
598 *
599 * ==== In case of a rare QR failure use
600 * . eigenvalues of the trailing 2-by-2
601 * . principal submatrix. ====
602 *
603  IF( ks.GE.kbot ) THEN
604  aa = h( kbot-1, kbot-1 )
605  cc = h( kbot, kbot-1 )
606  bb = h( kbot-1, kbot )
607  dd = h( kbot, kbot )
608  CALL dlanv2( aa, bb, cc, dd, wr( kbot-1 ),
609  $ wi( kbot-1 ), wr( kbot ),
610  $ wi( kbot ), cs, sn )
611  ks = kbot - 1
612  END IF
613  END IF
614 *
615  IF( kbot-ks+1.GT.ns ) THEN
616 *
617 * ==== Sort the shifts (Helps a little)
618 * . Bubble sort keeps complex conjugate
619 * . pairs together. ====
620 *
621  sorted = .false.
622  DO 50 k = kbot, ks + 1, -1
623  IF( sorted )
624  $ GO TO 60
625  sorted = .true.
626  DO 40 i = ks, k - 1
627  IF( abs( wr( i ) )+abs( wi( i ) ).LT.
628  $ abs( wr( i+1 ) )+abs( wi( i+1 ) ) ) THEN
629  sorted = .false.
630 *
631  swap = wr( i )
632  wr( i ) = wr( i+1 )
633  wr( i+1 ) = swap
634 *
635  swap = wi( i )
636  wi( i ) = wi( i+1 )
637  wi( i+1 ) = swap
638  END IF
639  40 CONTINUE
640  50 CONTINUE
641  60 CONTINUE
642  END IF
643 *
644 * ==== Shuffle shifts into pairs of real shifts
645 * . and pairs of complex conjugate shifts
646 * . assuming complex conjugate shifts are
647 * . already adjacent to one another. (Yes,
648 * . they are.) ====
649 *
650  DO 70 i = kbot, ks + 2, -2
651  IF( wi( i ).NE.-wi( i-1 ) ) THEN
652 *
653  swap = wr( i )
654  wr( i ) = wr( i-1 )
655  wr( i-1 ) = wr( i-2 )
656  wr( i-2 ) = swap
657 *
658  swap = wi( i )
659  wi( i ) = wi( i-1 )
660  wi( i-1 ) = wi( i-2 )
661  wi( i-2 ) = swap
662  END IF
663  70 CONTINUE
664  END IF
665 *
666 * ==== If there are only two shifts and both are
667 * . real, then use only one. ====
668 *
669  IF( kbot-ks+1.EQ.2 ) THEN
670  IF( wi( kbot ).EQ.zero ) THEN
671  IF( abs( wr( kbot )-h( kbot, kbot ) ).LT.
672  $ abs( wr( kbot-1 )-h( kbot, kbot ) ) ) THEN
673  wr( kbot-1 ) = wr( kbot )
674  ELSE
675  wr( kbot ) = wr( kbot-1 )
676  END IF
677  END IF
678  END IF
679 *
680 * ==== Use up to NS of the the smallest magnatiude
681 * . shifts. If there aren't NS shifts available,
682 * . then use them all, possibly dropping one to
683 * . make the number of shifts even. ====
684 *
685  ns = min( ns, kbot-ks+1 )
686  ns = ns - mod( ns, 2 )
687  ks = kbot - ns + 1
688 *
689 * ==== Small-bulge multi-shift QR sweep:
690 * . split workspace under the subdiagonal into
691 * . - a KDU-by-KDU work array U in the lower
692 * . left-hand-corner,
693 * . - a KDU-by-at-least-KDU-but-more-is-better
694 * . (KDU-by-NHo) horizontal work array WH along
695 * . the bottom edge,
696 * . - and an at-least-KDU-but-more-is-better-by-KDU
697 * . (NVE-by-KDU) vertical work WV arrow along
698 * . the left-hand-edge. ====
699 *
700  kdu = 3*ns - 3
701  ku = n - kdu + 1
702  kwh = kdu + 1
703  nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
704  kwv = kdu + 4
705  nve = n - kdu - kwv + 1
706 *
707 * ==== Small-bulge multi-shift QR sweep ====
708 *
709  CALL dlaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
710  $ wr( ks ), wi( ks ), h, ldh, iloz, ihiz, z,
711  $ ldz, work, 3, h( ku, 1 ), ldh, nve,
712  $ h( kwv, 1 ), ldh, nho, h( ku, kwh ), ldh )
713  END IF
714 *
715 * ==== Note progress (or the lack of it). ====
716 *
717  IF( ld.GT.0 ) THEN
718  ndfl = 1
719  ELSE
720  ndfl = ndfl + 1
721  END IF
722 *
723 * ==== End of main loop ====
724  80 CONTINUE
725 *
726 * ==== Iteration limit exceeded. Set INFO to show where
727 * . the problem occurred and exit. ====
728 *
729  info = kbot
730  90 CONTINUE
731  END IF
732 *
733 * ==== Return the optimal value of LWORK. ====
734 *
735  work( 1 ) = dble( lwkopt )
736 *
737 * ==== End of DLAQR4 ====
738 *
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
dlahqr
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.
Definition: dlahqr.f:209
dlaqr2
subroutine dlaqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
DLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
Definition: dlaqr2.f:280
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
dlaqr5
subroutine dlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
DLAQR5 performs a single small-bulge multi-shift QR sweep.
Definition: dlaqr5.f:261
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: