LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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,
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 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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
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 sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:167
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:106
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 slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:124
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 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 sormhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMHR
Definition: sormhr.f:179
subroutine strexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
STREXC
Definition: strexc.f:148
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
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: