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

◆ slaqr2()

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

SLAQR2 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 SLAQR2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!>    SLAQR2 is identical to SLAQR3 except that it avoids
!>    recursion by calling SLAHQR instead of SLAQR4.
!>
!>    Aggressive early deflation:
!>
!>    This subroutine 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; SLAQR2
!>          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 273 of file slaqr2.f.

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