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

◆ slaqr4()

subroutine slaqr4 ( logical  WANTT,
logical  WANTZ,
integer  N,
integer  ILO,
integer  IHI,
real, dimension( ldh, * )  H,
integer  LDH,
real, dimension( * )  WR,
real, dimension( * )  WI,
integer  ILOZ,
integer  IHIZ,
real, dimension( ldz, * )  Z,
integer  LDZ,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

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

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

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

    SLAQR4 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 SGEBAL, and then passed to SGEHRD when the
           matrix output by SGEBAL 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 REAL 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 REAL array, dimension (IHI)
[out]WI
          WI is REAL 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 REAL 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 REAL 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 SLAQR4 does a workspace query.
           In this case, SLAQR4 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
 \verbatim
          INFO is INTEGER
             = 0:  successful exit
             > 0:  if INFO = i, SLAQR4 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 a 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(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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 263 of file slaqr4.f.

265*
266* -- LAPACK auxiliary routine --
267* -- LAPACK is a software package provided by Univ. of Tennessee, --
268* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269*
270* .. Scalar Arguments ..
271 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
272 LOGICAL WANTT, WANTZ
273* ..
274* .. Array Arguments ..
275 REAL H( LDH, * ), WI( * ), WORK( * ), WR( * ),
276 $ Z( LDZ, * )
277* ..
278*
279* ================================================================
280*
281* .. Parameters ..
282*
283* ==== Matrices of order NTINY or smaller must be processed by
284* . SLAHQR because of insufficient subdiagonal scratch space.
285* . (This is a hard limit.) ====
286 INTEGER NTINY
287 parameter( ntiny = 15 )
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 REAL WILK1, WILK2
304 parameter( wilk1 = 0.75e0, wilk2 = -0.4375e0 )
305 REAL ZERO, ONE
306 parameter( zero = 0.0e0, one = 1.0e0 )
307* ..
308* .. Local Scalars ..
309 REAL 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 REAL ZDUM( 1, 1 )
323* ..
324* .. External Subroutines ..
325 EXTERNAL slacpy, slahqr, slanv2, slaqr2, slaqr5
326* ..
327* .. Intrinsic Functions ..
328 INTRINSIC abs, int, max, min, mod, real
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 SLAHQR. ====
343*
344 lwkopt = 1
345 IF( lwork.NE.-1 )
346 $ CALL slahqr( 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 = 15, 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.4.) ====
375*
376 nwr = ilaenv( 13, 'SLAQR4', 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 = 15, 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, 'SLAQR4', jbcmpz, n, ilo, ihi, lwork )
386 nsr = min( nsr, ( n-3 ) / 6, ihi-ilo )
387 nsr = max( 2, nsr-mod( nsr, 2 ) )
388*
389* ==== Estimate optimal workspace ====
390*
391* ==== Workspace query call to SLAQR2 ====
392*
393 CALL slaqr2( 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(SLAQR5, SLAQR2) ====
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 ) = real( lwkopt )
405 RETURN
406 END IF
407*
408* ==== SLAHQR/SLAQR0 crossover point ====
409*
410 nmin = ilaenv( 12, 'SLAQR4', jbcmpz, n, ilo, ihi, lwork )
411 nmin = max( ntiny, nmin )
412*
413* ==== Nibble crossover point ====
414*
415 nibble = ilaenv( 14, 'SLAQR4', 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, 'SLAQR4', 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-3 ) / 6, 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 slaqr2( 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 SLAQR2
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* . SLAQR2 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 slanv2( 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 SLAHQR
584* . on a trailing principal submatrix to
585* . get more. (Since NS.LE.NSMAX.LE.(N-3)/6,
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 slacpy( 'A', ns, ns, h( ks, ks ), ldh,
593 $ h( kt, 1 ), ldh )
594 CALL slahqr( .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 slanv2( 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 magnitude
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 = 2*ns
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 slaqr5( 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 ) = real( lwkopt )
736*
737* ==== End of SLAQR4 ====
738*
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine slanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition: slanv2.f:127
subroutine slaqr2(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)
SLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
Definition: slaqr2.f:278
subroutine slaqr5(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)
SLAQR5 performs a single small-bulge multi-shift QR sweep.
Definition: slaqr5.f:265
subroutine slahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
SLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition: slahqr.f:207
Here is the call graph for this function:
Here is the caller graph for this function: