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

◆ slaqr3()

subroutine slaqr3 ( logical  WANTT,
logical  WANTZ,
integer  N,
integer  KTOP,
integer  KBOT,
integer  NW,
real, dimension( ldh, * )  H,
integer  LDH,
integer  ILOZ,
integer  IHIZ,
real, dimension( ldz, * )  Z,
integer  LDZ,
integer  NS,
integer  ND,
real, dimension( * )  SR,
real, dimension( * )  SI,
real, dimension( ldv, * )  V,
integer  LDV,
integer  NH,
real, dimension( ldt, * )  T,
integer  LDT,
integer  NV,
real, dimension( ldwv, * )  WV,
integer  LDWV,
real, dimension( * )  WORK,
integer  LWORK 
)

SLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).

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

Purpose:
    Aggressive early deflation:

    SLAQR3 accepts as input an upper Hessenberg matrix
    H and performs an orthogonal similarity transformation
    designed to detect and deflate fully converged eigenvalues from
    a trailing principal submatrix.  On output H has been over-
    written by a new Hessenberg matrix that is a perturbation of
    an orthogonal similarity transformation of H.  It is to be
    hoped that the final version of H has many zero subdiagonal
    entries.
Parameters
[in]WANTT
          WANTT is LOGICAL
          If .TRUE., then the Hessenberg matrix H is fully updated
          so that the quasi-triangular Schur factor may be
          computed (in cooperation with the calling subroutine).
          If .FALSE., then only enough of H is updated to preserve
          the eigenvalues.
[in]WANTZ
          WANTZ is LOGICAL
          If .TRUE., then the orthogonal matrix Z is updated so
          so that the orthogonal Schur factor may be computed
          (in cooperation with the calling subroutine).
          If .FALSE., then Z is not referenced.
[in]N
          N is INTEGER
          The order of the matrix H and (if WANTZ is .TRUE.) the
          order of the orthogonal matrix Z.
[in]KTOP
          KTOP is INTEGER
          It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
          KBOT and KTOP together determine an isolated block
          along the diagonal of the Hessenberg matrix.
[in]KBOT
          KBOT is INTEGER
          It is assumed without a check that either
          KBOT = N or H(KBOT+1,KBOT)=0.  KBOT and KTOP together
          determine an isolated block along the diagonal of the
          Hessenberg matrix.
[in]NW
          NW is INTEGER
          Deflation window size.  1 <= NW <= (KBOT-KTOP+1).
[in,out]H
          H is REAL array, dimension (LDH,N)
          On input the initial N-by-N section of H stores the
          Hessenberg matrix undergoing aggressive early deflation.
          On output H has been transformed by an orthogonal
          similarity transformation, perturbed, and the returned
          to Hessenberg form that (it is to be hoped) has some
          zero subdiagonal entries.
[in]LDH
          LDH is INTEGER
          Leading dimension of H just as declared in the calling
          subroutine.  N <= LDH
[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 <= IHIZ <= N.
[in,out]Z
          Z is REAL array, dimension (LDZ,N)
          IF WANTZ is .TRUE., then on output, the orthogonal
          similarity transformation mentioned above has been
          accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
          If WANTZ is .FALSE., then Z is unreferenced.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of Z just as declared in the
          calling subroutine.  1 <= LDZ.
[out]NS
          NS is INTEGER
          The number of unconverged (ie approximate) eigenvalues
          returned in SR and SI that may be used as shifts by the
          calling subroutine.
[out]ND
          ND is INTEGER
          The number of converged eigenvalues uncovered by this
          subroutine.
[out]SR
          SR is REAL array, dimension (KBOT)
[out]SI
          SI is REAL array, dimension (KBOT)
          On output, the real and imaginary parts of approximate
          eigenvalues that may be used for shifts are stored in
          SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
          SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
          The real and imaginary parts of converged eigenvalues
          are stored in SR(KBOT-ND+1) through SR(KBOT) and
          SI(KBOT-ND+1) through SI(KBOT), respectively.
[out]V
          V is REAL array, dimension (LDV,NW)
          An NW-by-NW work array.
[in]LDV
          LDV is INTEGER
          The leading dimension of V just as declared in the
          calling subroutine.  NW <= LDV
[in]NH
          NH is INTEGER
          The number of columns of T.  NH >= NW.
[out]T
          T is REAL array, dimension (LDT,NW)
[in]LDT
          LDT is INTEGER
          The leading dimension of T just as declared in the
          calling subroutine.  NW <= LDT
[in]NV
          NV is INTEGER
          The number of rows of work array WV available for
          workspace.  NV >= NW.
[out]WV
          WV is REAL array, dimension (LDWV,NW)
[in]LDWV
          LDWV is INTEGER
          The leading dimension of W just as declared in the
          calling subroutine.  NW <= LDV
[out]WORK
          WORK is REAL array, dimension (LWORK)
          On exit, WORK(1) is set to an estimate of the optimal value
          of LWORK for the given values of N, NW, KTOP and KBOT.
[in]LWORK
          LWORK is INTEGER
          The dimension of the work array WORK.  LWORK = 2*NW
          suffices, but greater efficiency may result from larger
          values of LWORK.

          If LWORK = -1, then a workspace query is assumed; SLAQR3
          only estimates the optimal workspace size for the given
          values of N, NW, KTOP and KBOT.  The estimate is returned
          in WORK(1).  No error message related to LWORK is issued
          by XERBLA.  Neither H nor Z are 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

Definition at line 272 of file slaqr3.f.

275*
276* -- LAPACK auxiliary routine --
277* -- LAPACK is a software package provided by Univ. of Tennessee, --
278* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
279*
280* .. Scalar Arguments ..
281 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
282 $ LDZ, LWORK, N, ND, NH, NS, NV, NW
283 LOGICAL WANTT, WANTZ
284* ..
285* .. Array Arguments ..
286 REAL H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
287 $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
288 $ Z( LDZ, * )
289* ..
290*
291* ================================================================
292* .. Parameters ..
293 REAL ZERO, ONE
294 parameter( zero = 0.0e0, one = 1.0e0 )
295* ..
296* .. Local Scalars ..
297 REAL AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S,
298 $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
299 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
300 $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2, LWK3,
301 $ LWKOPT, NMIN
302 LOGICAL BULGE, SORTED
303* ..
304* .. External Functions ..
305 REAL SLAMCH
306 INTEGER ILAENV
307 EXTERNAL slamch, ilaenv
308* ..
309* .. External Subroutines ..
310 EXTERNAL scopy, sgehrd, sgemm, slabad, slacpy, slahqr,
311 $ slanv2, slaqr4, slarf, slarfg, slaset, sormhr,
312 $ strexc
313* ..
314* .. Intrinsic Functions ..
315 INTRINSIC abs, int, max, min, real, sqrt
316* ..
317* .. Executable Statements ..
318*
319* ==== Estimate optimal workspace. ====
320*
321 jw = min( nw, kbot-ktop+1 )
322 IF( jw.LE.2 ) THEN
323 lwkopt = 1
324 ELSE
325*
326* ==== Workspace query call to SGEHRD ====
327*
328 CALL sgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
329 lwk1 = int( work( 1 ) )
330*
331* ==== Workspace query call to SORMHR ====
332*
333 CALL sormhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
334 $ work, -1, info )
335 lwk2 = int( work( 1 ) )
336*
337* ==== Workspace query call to SLAQR4 ====
338*
339 CALL slaqr4( .true., .true., jw, 1, jw, t, ldt, sr, si, 1, jw,
340 $ v, ldv, work, -1, infqr )
341 lwk3 = int( work( 1 ) )
342*
343* ==== Optimal workspace ====
344*
345 lwkopt = max( jw+max( lwk1, lwk2 ), lwk3 )
346 END IF
347*
348* ==== Quick return in case of workspace query. ====
349*
350 IF( lwork.EQ.-1 ) THEN
351 work( 1 ) = real( lwkopt )
352 RETURN
353 END IF
354*
355* ==== Nothing to do ...
356* ... for an empty active block ... ====
357 ns = 0
358 nd = 0
359 work( 1 ) = one
360 IF( ktop.GT.kbot )
361 $ RETURN
362* ... nor for an empty deflation window. ====
363 IF( nw.LT.1 )
364 $ RETURN
365*
366* ==== Machine constants ====
367*
368 safmin = slamch( 'SAFE MINIMUM' )
369 safmax = one / safmin
370 CALL slabad( safmin, safmax )
371 ulp = slamch( 'PRECISION' )
372 smlnum = safmin*( real( n ) / ulp )
373*
374* ==== Setup deflation window ====
375*
376 jw = min( nw, kbot-ktop+1 )
377 kwtop = kbot - jw + 1
378 IF( kwtop.EQ.ktop ) THEN
379 s = zero
380 ELSE
381 s = h( kwtop, kwtop-1 )
382 END IF
383*
384 IF( kbot.EQ.kwtop ) THEN
385*
386* ==== 1-by-1 deflation window: not much to do ====
387*
388 sr( kwtop ) = h( kwtop, kwtop )
389 si( kwtop ) = zero
390 ns = 1
391 nd = 0
392 IF( abs( s ).LE.max( smlnum, ulp*abs( h( kwtop, kwtop ) ) ) )
393 $ THEN
394 ns = 0
395 nd = 1
396 IF( kwtop.GT.ktop )
397 $ h( kwtop, kwtop-1 ) = zero
398 END IF
399 work( 1 ) = one
400 RETURN
401 END IF
402*
403* ==== Convert to spike-triangular form. (In case of a
404* . rare QR failure, this routine continues to do
405* . aggressive early deflation using that part of
406* . the deflation window that converged using INFQR
407* . here and there to keep track.) ====
408*
409 CALL slacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
410 CALL scopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
411*
412 CALL slaset( 'A', jw, jw, zero, one, v, ldv )
413 nmin = ilaenv( 12, 'SLAQR3', 'SV', jw, 1, jw, lwork )
414 IF( jw.GT.nmin ) THEN
415 CALL slaqr4( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
416 $ si( kwtop ), 1, jw, v, ldv, work, lwork, infqr )
417 ELSE
418 CALL slahqr( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
419 $ si( kwtop ), 1, jw, v, ldv, infqr )
420 END IF
421*
422* ==== STREXC needs a clean margin near the diagonal ====
423*
424 DO 10 j = 1, jw - 3
425 t( j+2, j ) = zero
426 t( j+3, j ) = zero
427 10 CONTINUE
428 IF( jw.GT.2 )
429 $ t( jw, jw-2 ) = zero
430*
431* ==== Deflation detection loop ====
432*
433 ns = jw
434 ilst = infqr + 1
435 20 CONTINUE
436 IF( ilst.LE.ns ) THEN
437 IF( ns.EQ.1 ) THEN
438 bulge = .false.
439 ELSE
440 bulge = t( ns, ns-1 ).NE.zero
441 END IF
442*
443* ==== Small spike tip test for deflation ====
444*
445 IF( .NOT. bulge ) THEN
446*
447* ==== Real eigenvalue ====
448*
449 foo = abs( t( ns, ns ) )
450 IF( foo.EQ.zero )
451 $ foo = abs( s )
452 IF( abs( s*v( 1, ns ) ).LE.max( smlnum, ulp*foo ) ) THEN
453*
454* ==== Deflatable ====
455*
456 ns = ns - 1
457 ELSE
458*
459* ==== Undeflatable. Move it up out of the way.
460* . (STREXC can not fail in this case.) ====
461*
462 ifst = ns
463 CALL strexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
464 $ info )
465 ilst = ilst + 1
466 END IF
467 ELSE
468*
469* ==== Complex conjugate pair ====
470*
471 foo = abs( t( ns, ns ) ) + sqrt( abs( t( ns, ns-1 ) ) )*
472 $ sqrt( abs( t( ns-1, ns ) ) )
473 IF( foo.EQ.zero )
474 $ foo = abs( s )
475 IF( max( abs( s*v( 1, ns ) ), abs( s*v( 1, ns-1 ) ) ).LE.
476 $ max( smlnum, ulp*foo ) ) THEN
477*
478* ==== Deflatable ====
479*
480 ns = ns - 2
481 ELSE
482*
483* ==== Undeflatable. Move them up out of the way.
484* . Fortunately, STREXC does the right thing with
485* . ILST in case of a rare exchange failure. ====
486*
487 ifst = ns
488 CALL strexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
489 $ info )
490 ilst = ilst + 2
491 END IF
492 END IF
493*
494* ==== End deflation detection loop ====
495*
496 GO TO 20
497 END IF
498*
499* ==== Return to Hessenberg form ====
500*
501 IF( ns.EQ.0 )
502 $ s = zero
503*
504 IF( ns.LT.jw ) THEN
505*
506* ==== sorting diagonal blocks of T improves accuracy for
507* . graded matrices. Bubble sort deals well with
508* . exchange failures. ====
509*
510 sorted = .false.
511 i = ns + 1
512 30 CONTINUE
513 IF( sorted )
514 $ GO TO 50
515 sorted = .true.
516*
517 kend = i - 1
518 i = infqr + 1
519 IF( i.EQ.ns ) THEN
520 k = i + 1
521 ELSE IF( t( i+1, i ).EQ.zero ) THEN
522 k = i + 1
523 ELSE
524 k = i + 2
525 END IF
526 40 CONTINUE
527 IF( k.LE.kend ) THEN
528 IF( k.EQ.i+1 ) THEN
529 evi = abs( t( i, i ) )
530 ELSE
531 evi = abs( t( i, i ) ) + sqrt( abs( t( i+1, i ) ) )*
532 $ sqrt( abs( t( i, i+1 ) ) )
533 END IF
534*
535 IF( k.EQ.kend ) THEN
536 evk = abs( t( k, k ) )
537 ELSE IF( t( k+1, k ).EQ.zero ) THEN
538 evk = abs( t( k, k ) )
539 ELSE
540 evk = abs( t( k, k ) ) + sqrt( abs( t( k+1, k ) ) )*
541 $ sqrt( abs( t( k, k+1 ) ) )
542 END IF
543*
544 IF( evi.GE.evk ) THEN
545 i = k
546 ELSE
547 sorted = .false.
548 ifst = i
549 ilst = k
550 CALL strexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, work,
551 $ info )
552 IF( info.EQ.0 ) THEN
553 i = ilst
554 ELSE
555 i = k
556 END IF
557 END IF
558 IF( i.EQ.kend ) THEN
559 k = i + 1
560 ELSE IF( t( i+1, i ).EQ.zero ) THEN
561 k = i + 1
562 ELSE
563 k = i + 2
564 END IF
565 GO TO 40
566 END IF
567 GO TO 30
568 50 CONTINUE
569 END IF
570*
571* ==== Restore shift/eigenvalue array from T ====
572*
573 i = jw
574 60 CONTINUE
575 IF( i.GE.infqr+1 ) THEN
576 IF( i.EQ.infqr+1 ) THEN
577 sr( kwtop+i-1 ) = t( i, i )
578 si( kwtop+i-1 ) = zero
579 i = i - 1
580 ELSE IF( t( i, i-1 ).EQ.zero ) THEN
581 sr( kwtop+i-1 ) = t( i, i )
582 si( kwtop+i-1 ) = zero
583 i = i - 1
584 ELSE
585 aa = t( i-1, i-1 )
586 cc = t( i, i-1 )
587 bb = t( i-1, i )
588 dd = t( i, i )
589 CALL slanv2( aa, bb, cc, dd, sr( kwtop+i-2 ),
590 $ si( kwtop+i-2 ), sr( kwtop+i-1 ),
591 $ si( kwtop+i-1 ), cs, sn )
592 i = i - 2
593 END IF
594 GO TO 60
595 END IF
596*
597 IF( ns.LT.jw .OR. s.EQ.zero ) THEN
598 IF( ns.GT.1 .AND. s.NE.zero ) THEN
599*
600* ==== Reflect spike back into lower triangle ====
601*
602 CALL scopy( ns, v, ldv, work, 1 )
603 beta = work( 1 )
604 CALL slarfg( ns, beta, work( 2 ), 1, tau )
605 work( 1 ) = one
606*
607 CALL slaset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
608*
609 CALL slarf( 'L', ns, jw, work, 1, tau, t, ldt,
610 $ work( jw+1 ) )
611 CALL slarf( 'R', ns, ns, work, 1, tau, t, ldt,
612 $ work( jw+1 ) )
613 CALL slarf( 'R', jw, ns, work, 1, tau, v, ldv,
614 $ work( jw+1 ) )
615*
616 CALL sgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
617 $ lwork-jw, info )
618 END IF
619*
620* ==== Copy updated reduced window into place ====
621*
622 IF( kwtop.GT.1 )
623 $ h( kwtop, kwtop-1 ) = s*v( 1, 1 )
624 CALL slacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
625 CALL scopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
626 $ ldh+1 )
627*
628* ==== Accumulate orthogonal matrix in order update
629* . H and Z, if requested. ====
630*
631 IF( ns.GT.1 .AND. s.NE.zero )
632 $ CALL sormhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
633 $ work( jw+1 ), lwork-jw, info )
634*
635* ==== Update vertical slab in H ====
636*
637 IF( wantt ) THEN
638 ltop = 1
639 ELSE
640 ltop = ktop
641 END IF
642 DO 70 krow = ltop, kwtop - 1, nv
643 kln = min( nv, kwtop-krow )
644 CALL sgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
645 $ ldh, v, ldv, zero, wv, ldwv )
646 CALL slacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
647 70 CONTINUE
648*
649* ==== Update horizontal slab in H ====
650*
651 IF( wantt ) THEN
652 DO 80 kcol = kbot + 1, n, nh
653 kln = min( nh, n-kcol+1 )
654 CALL sgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
655 $ h( kwtop, kcol ), ldh, zero, t, ldt )
656 CALL slacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
657 $ ldh )
658 80 CONTINUE
659 END IF
660*
661* ==== Update vertical slab in Z ====
662*
663 IF( wantz ) THEN
664 DO 90 krow = iloz, ihiz, nv
665 kln = min( nv, ihiz-krow+1 )
666 CALL sgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
667 $ ldz, v, ldv, zero, wv, ldwv )
668 CALL slacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
669 $ ldz )
670 90 CONTINUE
671 END IF
672 END IF
673*
674* ==== Return the number of deflations ... ====
675*
676 nd = jw - ns
677*
678* ==== ... and the number of shifts. (Subtracting
679* . INFQR from the spike length takes care
680* . of the case of a rare QR failure while
681* . calculating eigenvalues of the deflation
682* . window.) ====
683*
684 ns = ns - infqr
685*
686* ==== Return optimal workspace. ====
687*
688 work( 1 ) = real( lwkopt )
689*
690* ==== End of SLAQR3 ====
691*
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
slaset
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
slacpy
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
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
sgehrd
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:167
slarfg
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:106
slanv2
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
slarf
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:124
slaqr4
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
slahqr
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
sormhr
subroutine sormhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMHR
Definition: sormhr.f:179
strexc
subroutine strexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
STREXC
Definition: strexc.f:148
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: