DGHUST

Computing the eigenvalues of a real skew-Hamiltonian/skew-Hamiltonian pencil (unfactored version)

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

     To compute the eigenvalues of a real N-by-N skew-Hamiltonian/
     skew-Hamiltonian pencil aS - bT with

           (  A  D  )         (  B  F  )
       S = (        ) and T = (        ).                           (1)
           (  E  A' )         (  G  B' )

     Optionally, if JOB = 'T', the pencil aS - bT will be transformed
     to the structured Schur form: an orthogonal transformation matrix
     Q is computed such that

                     (  Aout  Dout  )
       J Q' J' S Q = (              ), and
                     (   0    Aout' )
                                                                    (2)
                     (  Bout  Fout  )            (  0  I  )
       J Q' J' T Q = (              ), where J = (        ),
                     (   0    Bout' )            ( -I  0  )

     Aout is upper triangular, and Bout is upper quasi-triangular. The
     notation M' denotes the transpose of the matrix M.
     Optionally, if COMPQ = 'I' or COMPQ = 'U', the orthogonal
     transformation matrix Q will be computed.  
Specification
      SUBROUTINE DGHUST( JOB, COMPQ, N, A, LDA, DE, LDDE, B, LDB,
     $                   FG, LDFG, Q, LDQ, ALPHAR, ALPHAI, BETA, DWORK,
     $                   LDWORK, INFO )
C
C     .. Scalar Arguments ..
      CHARACTER          COMPQ, JOB
      INTEGER            INFO, LDA, LDB, LDDE, LDFG, LDQ, LDWORK, N
C
C     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
     $                   B( LDB, * ), BETA( * ), DE( LDDE, * ),
     $                   DWORK( * ), FG( LDFG, * ), Q( LDQ, * )
Arguments

Mode Parameters

     JOB     CHARACTER*1
             Specifies the computation to be performed, as follows:
             = 'E':  compute the eigenvalues only; S and T will not
                     necessarily be put into skew-Hamiltonian
                     triangular form (2);
             = 'T':  put S and T into skew-Hamiltonian triangular form
                     (2), and return the eigenvalues in ALPHAR, ALPHAI
                     and BETA.

     COMPQ   CHARACTER*1
             Specifies whether to compute the orthogonal transformation
             matrix Q as follows:
             = 'N':  Q is not computed;
             = 'I':  the array Q is initialized internally to the unit
                     matrix, and the orthogonal matrix Q is returned;
             = 'U':  the array Q contains an orthogonal matrix Q0 on
                     entry, and the product Q0*Q is returned, where Q
                     is the product of the orthogonal transformations
                     that are applied to the pencil aS - bT to reduce
                     S and T to the forms in (2), for COMPQ = 'I'.             
Input/Output Parameters
     N       (input) INTEGER
             The order of the pencil aS - bT.  N >= 0, even.

     A       (input/output) DOUBLE PRECISION array, dimension
                            (LDA, N/2)
             On entry, the leading N/2-by-N/2 part of this array must
             contain the matrix A.
             On exit, if JOB = 'T', the leading N/2-by-N/2 part of this
             array contains the matrix Aout; otherwise, it contains
             meaningless elements.

     LDA     INTEGER
             The leading dimension of the array A.  LDA >= MAX(1, N/2).

     DE      (input/output) DOUBLE PRECISION array, dimension
                            (LDDE, N/2+1)
             On entry, the leading N/2-by-N/2 strictly lower triangular
             part of this array must contain the strictly lower
             triangular part of the skew-symmetric matrix E, and the
             N/2-by-N/2 strictly upper triangular part of the submatrix
             in the columns 2 to N/2+1 of this array must contain the
             strictly upper triangular part of the skew-symmetric
             matrix D.
             The entries on the diagonal and the first superdiagonal of
             this array are not referenced, but are assumed to be zero.
             On exit, if JOB = 'T', the leading N/2-by-N/2 strictly
             upper triangular part of the submatrix in the columns
             2 to N/2+1 of this array contains the strictly upper
             triangular part of the skew-symmetric matrix Dout.
             If JOB = 'E', the leading N/2-by-N/2 strictly upper
             triangular part of the submatrix in the columns 2 to N/2+1
             of this array contains the strictly upper triangular part
             of the skew-symmetric matrix D just before the application
             of the QZ algorithm. The remaining entries are
             meaningless.

     LDDE    INTEGER
             The leading dimension of the array DE.
             LDDE >= MAX(1, N/2).

     B       (input/output) DOUBLE PRECISION array, dimension
                            (LDB, N/2)
             On entry, the leading N/2-by-N/2 part of this array must
             contain the matrix B.
             On exit, if JOB = 'T', the leading N/2-by-N/2 part of this
             array contains the matrix Bout; otherwise, it contains
             meaningless elements.

     LDB     INTEGER
             The leading dimension of the array B.  LDB >= MAX(1, N/2).

     FG      (input/output) DOUBLE PRECISION array, dimension
                            (LDFG, N/2+1)
             On entry, the leading N/2-by-N/2 strictly lower triangular
             part of this array must contain the strictly lower
             triangular part of the skew-symmetric matrix G, and the
             N/2-by-N/2 strictly upper triangular part of the submatrix
             in the columns 2 to N/2+1 of this array must contain the
             strictly upper triangular part of the skew-symmetric
             matrix F.
             The entries on the diagonal and the first superdiagonal of
             this array are not referenced, but are assumed to be zero.
             On exit, if JOB = 'T', the leading N/2-by-N/2 strictly
             upper triangular part of the submatrix in the columns
             2 to N/2+1 of this array contains the strictly upper
             triangular part of the skew-symmetric matrix Fout.
             If JOB = 'E', the leading N/2-by-N/2 strictly upper
             triangular part of the submatrix in the columns 2 to N/2+1
             of this array contains the strictly upper triangular part
             of the skew-symmetric matrix F just before the application
             of the QZ algorithm. The remaining entries are
             meaningless.

     LDFG    INTEGER
             The leading dimension of the array FG.
             LDFG >= MAX(1, N/2).

     Q       (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
             On entry, if COMPQ = 'U', then the leading N-by-N part of
             this array must contain a given matrix Q0, and on exit,
             the leading N-by-N part of this array contains the product
             of the input matrix Q0 and the transformation matrix Q
             used to transform the matrices S and T.
             On exit, if COMPQ = 'I', then the leading N-by-N part of
             this array contains the orthogonal transformation matrix
             Q.
             If COMPQ = 'N', this array is not referenced.

     LDQ     INTEGER
             The leading dimension of the array Q.
             LDQ >= 1,         if COMPQ = 'N';
             LDQ >= MAX(1, N), if COMPQ = 'I' or COMPQ = 'U'.

     ALPHAR  (output) DOUBLE PRECISION array, dimension (N/2)
             The real parts of each scalar alpha defining an eigenvalue
             of the pencil aS - bT.

     ALPHAI  (output) DOUBLE PRECISION array, dimension (N/2)
             The imaginary parts of each scalar alpha defining an
             eigenvalue of the pencil aS - bT.
             If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
             positive, then the j-th and (j+1)-st eigenvalues are a
             complex conjugate pair.

     BETA    (output) DOUBLE PRECISION array, dimension (N/2)
             The scalars beta that define the eigenvalues of the pencil
             aS - bT.
             Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
             beta = BETA(j) represent the j-th eigenvalue of the pencil
             aS - bT, in the form lambda = alpha/beta. Since lambda may
             overflow, the ratios should not, in general, be computed.
             Due to the skew-Hamiltonian/skew-Hamiltonian structure of
             the pencil, every eigenvalue occurs twice and thus it has
             only to be saved once in ALPHAR, ALPHAI and BETA.
Workspace
     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
             On exit, if INFO = 0, DWORK(1) returns the optimal LDWORK.
             On exit, if INFO = -18, DWORK(1) returns the minimum
             value of LDWORK.

     LDWORK  INTEGER
             The dimension of the array DWORK.
             LDWORK >= MAX(1,N/2),        if JOB = 'E' and COMPQ = 'N';
             LDWORK >= MAX(1,N**2/4+N/2), if JOB = 'T' and COMPQ = 'N';
             LDWORK >= MAX(1,3*N**2/4),   if               COMPQ<> 'N'.
             For good performance LDWORK should generally be larger.

             If LDWORK = -1, then a workspace query is assumed;
             the routine only calculates the optimal size of the
             DWORK array, returns this value as the first entry of
             the DWORK array, and no error message related to LDWORK
             is issued by XERBLA.
Error Indicator
     INFO    INTEGER
             = 0: succesful exit;
             < 0: if INFO = -i, the i-th argument had an illegal value;
             = 1: QZ iteration failed in the LAPACK Library routine
                  DHGEQZ. (QZ iteration did not converge or computation
                  of the shifts failed.)
Method
     The algorithm uses Givens rotations and Householder reflections to
     annihilate elements in S and T such that S is in skew-Hamiltonian
     triangular form and T is in skew-Hamiltonian Hessenberg form:

         (  A1  D1  )      (  B1  F1  )
     S = (          ), T = (          ),
         (   0  A1' )      (   0  B1' )

     where A1 is upper triangular and B1 is upper Hessenberg.
     Subsequently, the QZ algorithm is applied to the pencil aA1 - bB1
     to determine orthogonal matrices Q1 and Q2 such that
     Q2' A1 Q1 is upper triangular and Q2' B1 Q1 is upper quasi-
     triangular.
     See also page 40 in [1] for more details.
References
     [1] Benner, P., Byers, R., Mehrmann, V. and Xu, H.
         Numerical Computation of Deflating Subspaces of Embedded
         Hamiltonian Pencils.
         Tech. Rep. SFB393/99-15, Technical University Chemnitz,
         Germany, June 1999.
Numerical Aspects
     
                                                               3
     The algorithm is numerically backward stable and needs O(N )
     real floating point operations.
Further Comments
   
     None.
Example

Program Text

*     DGHUST EXAMPLE PROGRAM TEXT
*
*     .. Parameters ..
      INTEGER            NIN, NOUT
      PARAMETER          ( NIN = 5, NOUT = 6 )
      INTEGER            NMAX
      PARAMETER          ( NMAX = 50 )
      INTEGER            LDA, LDB, LDDE, LDFG, LDQ, LDWORK
      PARAMETER          ( LDA  = NMAX/2, LDB = NMAX/2, LDDE = NMAX/2,
     $                     LDFG = NMAX/2, LDQ = NMAX, LDWORK =
     $                     3*NMAX*NMAX/4 )
*
*     .. Local Scalars ..
      CHARACTER          COMPQ, JOB
      INTEGER            I, INFO, J, N
*
*     .. Local Arrays ..
      DOUBLE PRECISION   A( LDA, NMAX/2 ),  ALPHAI( NMAX/2 ),
     $                   ALPHAR( NMAX/2 ),  B( LDB, NMAX/2 ),
     $                   BETA( NMAX/2 ),  DE( LDDE, NMAX/2+1 ),
     $                   DWORK( LDWORK ), FG( LDFG, NMAX/2+1 ),
     $                   Q( LDQ, NMAX )
*
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*
*     .. External Subroutines ..
      EXTERNAL           DGHUST
*
*     .. Intrinsic Functions ..
      INTRINSIC          MOD
*
*     .. Executable statements ..
*
      WRITE( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read in the data.
      READ( NIN, FMT = * )
      READ( NIN, FMT = * ) JOB, COMPQ, N
      READ( NIN, FMT = * ) ( (  A( I, J ), J = 1, N/2 ),   I = 1, N/2 )
      READ( NIN, FMT = * ) ( ( DE( I, J ), J = 1, N/2+1 ), I = 1, N/2 )
      READ( NIN, FMT = * ) ( (  B( I, J ), J = 1, N/2 ),   I = 1, N/2 )
      READ( NIN, FMT = * ) ( ( FG( I, J ), J = 1, N/2+1 ), I = 1, N/2 )
      IF( LSAME( COMPQ, 'U' ) )
     $   READ( NIN, FMT = * ) ( ( Q( I, J ), J = 1, N ), I = 1, N )
      IF( N.LT.0 .OR. N.GT.NMAX .OR. MOD( N, 2 ).NE.0 ) THEN
         WRITE( NOUT, FMT = 99998 ) N
      ELSE
*        Compute the eigenvalues of a real skew-Hamiltonian/
*        skew-Hamiltonian pencil (unfactored version).
         CALL DGHUST( JOB, COMPQ, N, A, LDA, DE, LDDE, B, LDB, FG, LDFG,
     $                Q, LDQ, ALPHAR, ALPHAI, BETA, DWORK, LDWORK,
     $                INFO )
         IF( INFO.NE.0 ) THEN
            WRITE( NOUT, FMT = 99997 ) INFO
         ELSE
            IF( LSAME( JOB, 'T' ) ) THEN
               WRITE( NOUT, FMT = 99996 )
               DO 10 I = 1, N/2
                  WRITE( NOUT, FMT = 99995 ) ( A( I, J ), J = 1, N/2 )
   10          CONTINUE
            END IF
            WRITE( NOUT, FMT = 99994 )
            DO 20 I = 1, N/2
               WRITE( NOUT, FMT = 99995 ) ( DE( I, J ), J = 1, N/2+1 )
   20       CONTINUE
            IF( LSAME( JOB, 'T' ) ) THEN
               WRITE( NOUT, FMT = 99993 )
               DO 30 I = 1, N/2
                  WRITE( NOUT, FMT = 99995 ) ( B( I, J ), J = 1, N/2 )
   30          CONTINUE
            END IF
            WRITE( NOUT, FMT = 99992 )
            DO 40 I = 1, N/2
               WRITE( NOUT, FMT = 99995 ) ( FG( I, J ), J = 1, N/2+1 )
   40       CONTINUE
            IF( .NOT.LSAME( COMPQ, 'N' ) ) THEN
               WRITE( NOUT, FMT = 99991 )
               DO 50 I = 1, N
                  WRITE( NOUT, FMT = 99995 ) ( Q( I, J ), J = 1, N )
   50          CONTINUE
            END IF
            WRITE( NOUT, FMT = 99990 )
            WRITE( NOUT, FMT = 99995 ) ( ALPHAR( I ), I = 1, N/2 )
            WRITE( NOUT, FMT = 99989 )
            WRITE( NOUT, FMT = 99995 ) ( ALPHAI( I ), I = 1, N/2 )
            WRITE( NOUT, FMT = 99988 )
            WRITE( NOUT, FMT = 99995 ) (   BETA( I ), I = 1, N/2 )
         END IF
      END IF
      STOP
*
99999 FORMAT( 'DGHUST EXAMPLE PROGRAM RESULTS', 1X )
99998 FORMAT( 'N is out of range.', /, 'N = ', I5 )
99997 FORMAT( 'INFO on exit from DGHUST = ', I2 )
99996 FORMAT( 'The transformed matrix A is ' )
99995 FORMAT( 51( 1X, F8.4 ) )
99994 FORMAT( 'The transformed matrix DE is ' )
99993 FORMAT( 'The transformed matrix B is ' )
99992 FORMAT( 'The transformed matrix FG is ' )
99991 FORMAT( 'The transformed matrix Q is ' )
99990 FORMAT( 'The vector ALPHAR is ' )
99989 FORMAT( 'The vector ALPHAI is ' )
99988 FORMAT( 'The vector BETA is ' )
      END
Program Data
DGHUST EXAMPLE PROGRAM DATA
    T   I   8
    0.8147    0.6323    0.9575    0.9571
    0.9057    0.0975    0.9648    0.4853
    0.1269    0.2784    0.1576    0.8002
    0.9133    0.5468    0.9705    0.1418
    0.4217    0.6557    0.6787    0.6554    0.2769
    0.9157    0.0357    0.7577    0.1711    0.0461
    0.7922    0.8491    0.7431    0.7060    0.0971
    0.9594    0.9339    0.3922    0.0318    0.8234
    0.6948    0.4387    0.1868    0.7093
    0.3170    0.3815    0.4897    0.7546
    0.9502    0.7655    0.4455    0.2760
    0.0344    0.7951    0.6463    0.6797
    0.6550    0.9597    0.7512    0.8909    0.1492
    0.1626    0.3403    0.2550    0.9592    0.2575
    0.1189    0.5852    0.5059    0.5472    0.8407
    0.4983    0.2238    0.6990    0.1386    0.2542
Program Results
DGHUST EXAMPLE PROGRAM RESULTS
The transformed matrix A is 
   0.0550  -0.3064   0.0026  -0.2663
   0.0000   1.2189  -1.0643   0.9107
   0.0000   0.0000   1.7796  -1.7196
   0.0000   0.0000   0.0000   1.0940
The transformed matrix DE is 
   0.4217   0.6557   0.4277  -0.1704   0.4606
  -1.5448   0.0357   0.7577  -0.2528   0.4803
   0.3220   1.1004   0.7431   0.7060  -1.1102
   0.3899   0.3463   0.4843   0.0318   0.8234
The transformed matrix B is 
   0.8482  -0.4425   0.7842  -0.4605
   0.0000  -0.5919  -0.3853   0.6993
   0.0000   0.0000   1.3990  -0.8832
   0.0000   0.0000   0.0000   1.5590
The transformed matrix FG is 
   0.6550   0.9597   0.3266  -0.4615   1.0389
  -0.3151   0.3403   0.2550   0.0161   0.5637
   0.1987  -0.0141   0.5059   0.5472  -0.4172
   0.7914   0.0371   0.0000   0.1386   0.2542
The transformed matrix Q is 
   0.0762   0.7824  -0.1002   0.6100   0.0000   0.0000   0.0000   0.0000
   0.1786  -0.2450  -0.6597   0.1836   0.1843  -0.5857  -0.1791  -0.1733
  -0.2452   0.1425  -0.4346  -0.2236  -0.1834  -0.0856   0.7473   0.2851
  -0.2123  -0.3870   0.3710   0.5839   0.3061  -0.1791   0.4481   0.0027
  -0.4085   0.0154  -0.1068   0.0137   0.3908   0.0160  -0.3745   0.7267
   0.6418  -0.2868  -0.2150   0.2523  -0.0286   0.4768   0.1132   0.3938
  -0.0077  -0.1340   0.1904   0.2041  -0.7708  -0.3910  -0.1748   0.3557
  -0.5275  -0.2394  -0.3666   0.3127  -0.3018   0.4869  -0.1581  -0.2812
The vector ALPHAR is 
   0.8482  -0.5919   1.3990   1.5590
The vector ALPHAI is 
   0.0000   0.0000   0.0000   0.0000
The vector BETA is 
   0.0550   1.2189   1.7796   1.0940

Return to index