LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
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
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 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 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
Here is the call graph for this function:
Here is the caller graph for this function: