LAPACK 3.12.0
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 275 of file slaqr2.f.

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