LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ shseqr()

subroutine shseqr ( character  JOB,
character  COMPZ,
integer  N,
integer  ILO,
integer  IHI,
real, dimension( ldh, * )  H,
integer  LDH,
real, dimension( * )  WR,
real, dimension( * )  WI,
real, dimension( ldz, * )  Z,
integer  LDZ,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

SHSEQR

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

Purpose:
    SHSEQR 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]JOB
          JOB is CHARACTER*1
           = 'E':  compute eigenvalues only;
           = 'S':  compute eigenvalues and the Schur form T.
[in]COMPZ
          COMPZ is CHARACTER*1
           = 'N':  no Schur vectors are computed;
           = 'I':  Z is initialized to the unit matrix and the matrix Z
                   of Schur vectors of H is returned;
           = 'V':  Z must contain an orthogonal matrix Q on entry, and
                   the product Q*Z is returned.
[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. ILO and IHI are normally
           set by a previous call to SGEBAL, and then passed to ZGEHRD
           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 JOB = 'S', 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 JOB = 'E', the
           contents of H are unspecified on exit.  (The output value of
           H when INFO > 0 is given under the description of INFO
           below.)

           Unlike earlier versions of SHSEQR, this subroutine may
           explicitly 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 (N)
[out]WI
          WI is REAL array, dimension (N)

           The real and imaginary parts, respectively, of the computed
           eigenvalues. 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 JOB = 'S', 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,out]Z
          Z is REAL array, dimension (LDZ,N)
           If COMPZ = 'N', Z is not referenced.
           If COMPZ = 'I', on entry Z need not be set and on exit,
           if INFO = 0, Z contains the orthogonal matrix Z of the Schur
           vectors of H.  If COMPZ = 'V', on entry Z must contain an
           N-by-N matrix Q, which is assumed to be equal to the unit
           matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
           if INFO = 0, Z contains Q*Z.
           Normally Q is the orthogonal matrix generated by SORGHR
           after the call to SGEHRD which formed the Hessenberg matrix
           H. (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 COMPZ = 'I' or
           COMPZ = 'V', then LDZ >= MAX(1,N).  Otherwise, LDZ >= 1.
[out]WORK
          WORK is REAL array, dimension (LWORK)
           On exit, if INFO = 0, 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 and delivers very good and sometimes
           optimal performance.  However, LWORK as large as 11*N
           may be required for optimal performance.  A workspace
           query is recommended to determine the optimal workspace
           size.

           If LWORK = -1, then SHSEQR does a workspace query.
           In this case, SHSEQR 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, the i-th argument had an illegal
                    value
             > 0:  if INFO = i, SHSEQR 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 JOB = 'E', 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 JOB   = 'S', 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 COMPZ = 'V', then on exit

                  (final value of Z)  =  (initial value of Z)*U

                where U is the orthogonal matrix in (*) (regard-
                less of the value of JOB.)

                If INFO > 0 and COMPZ = 'I', then on exit
                      (final value of Z)  = U
                where U is the orthogonal matrix in (*) (regard-
                less of the value of JOB.)

                If INFO > 0 and COMPZ = 'N', 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
Further Details:
             Default values supplied by
             ILAENV(ISPEC,'SHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
             It is suggested that these defaults be adjusted in order
             to attain best performance in each particular
             computational environment.

            ISPEC=12: The SLAHQR vs SLAQR0 crossover point.
                      Default: 75. (Must be at least 11.)

            ISPEC=13: Recommended deflation window size.
                      This depends on ILO, IHI and NS.  NS is the
                      number of simultaneous shifts returned
                      by ILAENV(ISPEC=15).  (See ISPEC=15 below.)
                      The default for (IHI-ILO+1) <= 500 is NS.
                      The default for (IHI-ILO+1) >  500 is 3*NS/2.

            ISPEC=14: Nibble crossover point. (See IPARMQ for
                      details.)  Default: 14% of deflation window
                      size.

            ISPEC=15: Number of simultaneous shifts in a multishift
                      QR iteration.

                      If IHI-ILO+1 is ...

                      greater than      ...but less    ... the
                      or equal to ...      than        default is

                           1               30          NS =   2(+)
                          30               60          NS =   4(+)
                          60              150          NS =  10(+)
                         150              590          NS =  **
                         590             3000          NS =  64
                        3000             6000          NS = 128
                        6000             infinity      NS = 256

                  (+)  By default some or all matrices of this order
                       are passed to the implicit double shift routine
                       SLAHQR and this parameter is ignored.  See
                       ISPEC=12 above and comments in IPARMQ for
                       details.

                 (**)  The asterisks (**) indicate an ad-hoc
                       function of N increasing from 10 to 64.

            ISPEC=16: Select structured matrix multiply.
                      If the number of simultaneous shifts (specified
                      by ISPEC=15) is less than 14, then the default
                      for ISPEC=16 is 0.  Otherwise the default for
                      ISPEC=16 is 2.
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 314 of file shseqr.f.

316 *
317 * -- LAPACK computational routine --
318 * -- LAPACK is a software package provided by Univ. of Tennessee, --
319 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
320 *
321 * .. Scalar Arguments ..
322  INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N
323  CHARACTER COMPZ, JOB
324 * ..
325 * .. Array Arguments ..
326  REAL H( LDH, * ), WI( * ), WORK( * ), WR( * ),
327  $ Z( LDZ, * )
328 * ..
329 *
330 * =====================================================================
331 *
332 * .. Parameters ..
333 *
334 * ==== Matrices of order NTINY or smaller must be processed by
335 * . SLAHQR because of insufficient subdiagonal scratch space.
336 * . (This is a hard limit.) ====
337  INTEGER NTINY
338  parameter( ntiny = 15 )
339 *
340 * ==== NL allocates some local workspace to help small matrices
341 * . through a rare SLAHQR failure. NL > NTINY = 15 is
342 * . required and NL <= NMIN = ILAENV(ISPEC=12,...) is recom-
343 * . mended. (The default value of NMIN is 75.) Using NL = 49
344 * . allows up to six simultaneous shifts and a 16-by-16
345 * . deflation window. ====
346  INTEGER NL
347  parameter( nl = 49 )
348  REAL ZERO, ONE
349  parameter( zero = 0.0e0, one = 1.0e0 )
350 * ..
351 * .. Local Arrays ..
352  REAL HL( NL, NL ), WORKL( NL )
353 * ..
354 * .. Local Scalars ..
355  INTEGER I, KBOT, NMIN
356  LOGICAL INITZ, LQUERY, WANTT, WANTZ
357 * ..
358 * .. External Functions ..
359  INTEGER ILAENV
360  LOGICAL LSAME
361  EXTERNAL ilaenv, lsame
362 * ..
363 * .. External Subroutines ..
364  EXTERNAL slacpy, slahqr, slaqr0, slaset, xerbla
365 * ..
366 * .. Intrinsic Functions ..
367  INTRINSIC max, min, real
368 * ..
369 * .. Executable Statements ..
370 *
371 * ==== Decode and check the input parameters. ====
372 *
373  wantt = lsame( job, 'S' )
374  initz = lsame( compz, 'I' )
375  wantz = initz .OR. lsame( compz, 'V' )
376  work( 1 ) = real( max( 1, n ) )
377  lquery = lwork.EQ.-1
378 *
379  info = 0
380  IF( .NOT.lsame( job, 'E' ) .AND. .NOT.wantt ) THEN
381  info = -1
382  ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
383  info = -2
384  ELSE IF( n.LT.0 ) THEN
385  info = -3
386  ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
387  info = -4
388  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
389  info = -5
390  ELSE IF( ldh.LT.max( 1, n ) ) THEN
391  info = -7
392  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.max( 1, n ) ) ) THEN
393  info = -11
394  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
395  info = -13
396  END IF
397 *
398  IF( info.NE.0 ) THEN
399 *
400 * ==== Quick return in case of invalid argument. ====
401 *
402  CALL xerbla( 'SHSEQR', -info )
403  RETURN
404 *
405  ELSE IF( n.EQ.0 ) THEN
406 *
407 * ==== Quick return in case N = 0; nothing to do. ====
408 *
409  RETURN
410 *
411  ELSE IF( lquery ) THEN
412 *
413 * ==== Quick return in case of a workspace query ====
414 *
415  CALL slaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,
416  $ ihi, z, ldz, work, lwork, info )
417 * ==== Ensure reported workspace size is backward-compatible with
418 * . previous LAPACK versions. ====
419  work( 1 ) = max( real( max( 1, n ) ), work( 1 ) )
420  RETURN
421 *
422  ELSE
423 *
424 * ==== copy eigenvalues isolated by SGEBAL ====
425 *
426  DO 10 i = 1, ilo - 1
427  wr( i ) = h( i, i )
428  wi( i ) = zero
429  10 CONTINUE
430  DO 20 i = ihi + 1, n
431  wr( i ) = h( i, i )
432  wi( i ) = zero
433  20 CONTINUE
434 *
435 * ==== Initialize Z, if requested ====
436 *
437  IF( initz )
438  $ CALL slaset( 'A', n, n, zero, one, z, ldz )
439 *
440 * ==== Quick return if possible ====
441 *
442  IF( ilo.EQ.ihi ) THEN
443  wr( ilo ) = h( ilo, ilo )
444  wi( ilo ) = zero
445  RETURN
446  END IF
447 *
448 * ==== SLAHQR/SLAQR0 crossover point ====
449 *
450  nmin = ilaenv( 12, 'SHSEQR', job( : 1 ) // compz( : 1 ), n,
451  $ ilo, ihi, lwork )
452  nmin = max( ntiny, nmin )
453 *
454 * ==== SLAQR0 for big matrices; SLAHQR for small ones ====
455 *
456  IF( n.GT.nmin ) THEN
457  CALL slaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,
458  $ ihi, z, ldz, work, lwork, info )
459  ELSE
460 *
461 * ==== Small matrix ====
462 *
463  CALL slahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,
464  $ ihi, z, ldz, info )
465 *
466  IF( info.GT.0 ) THEN
467 *
468 * ==== A rare SLAHQR failure! SLAQR0 sometimes succeeds
469 * . when SLAHQR fails. ====
470 *
471  kbot = info
472 *
473  IF( n.GE.nl ) THEN
474 *
475 * ==== Larger matrices have enough subdiagonal scratch
476 * . space to call SLAQR0 directly. ====
477 *
478  CALL slaqr0( wantt, wantz, n, ilo, kbot, h, ldh, wr,
479  $ wi, ilo, ihi, z, ldz, work, lwork, info )
480 *
481  ELSE
482 *
483 * ==== Tiny matrices don't have enough subdiagonal
484 * . scratch space to benefit from SLAQR0. Hence,
485 * . tiny matrices must be copied into a larger
486 * . array before calling SLAQR0. ====
487 *
488  CALL slacpy( 'A', n, n, h, ldh, hl, nl )
489  hl( n+1, n ) = zero
490  CALL slaset( 'A', nl, nl-n, zero, zero, hl( 1, n+1 ),
491  $ nl )
492  CALL slaqr0( wantt, wantz, nl, ilo, kbot, hl, nl, wr,
493  $ wi, ilo, ihi, z, ldz, workl, nl, info )
494  IF( wantt .OR. info.NE.0 )
495  $ CALL slacpy( 'A', n, n, hl, nl, h, ldh )
496  END IF
497  END IF
498  END IF
499 *
500 * ==== Clear out the trash, if necessary. ====
501 *
502  IF( ( wantt .OR. info.NE.0 ) .AND. n.GT.2 )
503  $ CALL slaset( 'L', n-2, n-2, zero, zero, h( 3, 1 ), ldh )
504 *
505 * ==== Ensure reported workspace size is backward-compatible with
506 * . previous LAPACK versions. ====
507 *
508  work( 1 ) = max( real( max( 1, n ) ), work( 1 ) )
509  END IF
510 *
511 * ==== End of SHSEQR ====
512 *
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine slaqr0(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
SLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: slaqr0.f:256
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
Here is the call graph for this function:
Here is the caller graph for this function: