LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dlaqr0()

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 >= 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 > 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 > 0, then 1 <= ILO <= IHI <= 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) < 0. If INFO = 0 and WANTT is
           .FALSE., then the contents of H are unspecified on exit.
           (The output value of H when INFO > 0 is given under the
           description of INFO below.)

           This subroutine may explicitly set H(i,j) = 0 for i > 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 >= 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) > 0 and WI(i+1) < 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 <= ILOZ <= ILO; IHI <= IHIZ <= 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 > 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 >= MAX(1,IHIZ).  Otherwise, LDZ >= 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 >= 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
             > 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 > 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 > 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 > 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 > 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.

Definition at line 254 of file dlaqr0.f.

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