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

◆ slaqr0()

subroutine slaqr0 ( 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 
)

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

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

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