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

◆ dhseqr()

subroutine dhseqr ( character  job,
character  compz,
integer  n,
integer  ilo,
integer  ihi,
double precision, dimension( ldh, * )  h,
integer  ldh,
double precision, dimension( * )  wr,
double precision, dimension( * )  wi,
double precision, dimension( ldz, * )  z,
integer  ldz,
double precision, dimension( * )  work,
integer  lwork,
integer  info 
)

DHSEQR

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

Purpose:
    DHSEQR 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 DGEBAL, and then passed to ZGEHRD
           when the matrix output by DGEBAL 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 DOUBLE PRECISION 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 DHSEQR, 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 DOUBLE PRECISION array, dimension (N)
[out]WI
          WI is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DORGHR
           after the call to DGEHRD 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 DOUBLE PRECISION 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 DHSEQR does a workspace query.
           In this case, DHSEQR 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, DHSEQR 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,'DHSEQR',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 DLAHQR vs DLAQR0 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
                       DLAHQR 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 dhseqr.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 DOUBLE PRECISION 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* . DLAHQR 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 DLAHQR 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 DOUBLE PRECISION ZERO, ONE
349 parameter( zero = 0.0d0, one = 1.0d0 )
350* ..
351* .. Local Arrays ..
352 DOUBLE PRECISION 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 dlacpy, dlahqr, dlaqr0, dlaset, xerbla
365* ..
366* .. Intrinsic Functions ..
367 INTRINSIC dble, max, min
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 ) = dble( 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( 'DHSEQR', -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 dlaqr0( 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( dble( max( 1, n ) ), work( 1 ) )
420 RETURN
421*
422 ELSE
423*
424* ==== copy eigenvalues isolated by DGEBAL ====
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 dlaset( '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* ==== DLAHQR/DLAQR0 crossover point ====
449*
450 nmin = ilaenv( 12, 'DHSEQR', job( : 1 ) // compz( : 1 ), n,
451 $ ilo, ihi, lwork )
452 nmin = max( ntiny, nmin )
453*
454* ==== DLAQR0 for big matrices; DLAHQR for small ones ====
455*
456 IF( n.GT.nmin ) THEN
457 CALL dlaqr0( 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 dlahqr( 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 DLAHQR failure! DLAQR0 sometimes succeeds
469* . when DLAHQR fails. ====
470*
471 kbot = info
472*
473 IF( n.GE.nl ) THEN
474*
475* ==== Larger matrices have enough subdiagonal scratch
476* . space to call DLAQR0 directly. ====
477*
478 CALL dlaqr0( 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 DLAQR0. Hence,
485* . tiny matrices must be copied into a larger
486* . array before calling DLAQR0. ====
487*
488 CALL dlacpy( 'A', n, n, h, ldh, hl, nl )
489 hl( n+1, n ) = zero
490 CALL dlaset( 'A', nl, nl-n, zero, zero, hl( 1, n+1 ),
491 $ nl )
492 CALL dlaqr0( 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 dlacpy( '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 dlaset( '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( dble( max( 1, n ) ), work( 1 ) )
509 END IF
510*
511* ==== End of DHSEQR ====
512*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
subroutine dlahqr(wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, iloz, ihiz, z, ldz, info)
DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition dlahqr.f:207
subroutine dlaqr0(wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, iloz, ihiz, z, ldz, work, lwork, info)
DLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition dlaqr0.f:256
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: