LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlaqr0 ( 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 
)

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

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

Purpose:
    DLAQR0 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 DLAQR0 does a workspace query.
           In this case, DLAQR0 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, DLAQR0 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 an orthogonal matrix.  The final
                value of H is upper Hessenberg and quasi-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.
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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 258 of file dlaqr0.f.

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