DGHFST

Computing the eigenvalues of a real skew-Hamiltonian/skew-Hamiltonian pencil (factored 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

                             (  B  F  )            (  0  I  )
       S = J Z' J' Z and T = (        ), where J = (        ).      (1)
                             (  G  B' )            ( -I  0  )

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

                (  Z11  Z12  )
       U' Z Q = (            ) = Zout, and
                (   0   Z22  )
                                                                    (2)
                     (  Bout  Fout  )
       J Q' J' T Q = (              ),
                     (   0    Bout' )

     where Z11 and Z22' are upper triangular and Bout is upper quasi-
     triangular. The notation M' denotes the transpose of the matrix M.
     Optionally, if COMPQ = 'I', the orthogonal transformation matrix Q
     will be computed.
     Optionally, if COMPU = 'I' or COMPU = 'U', the orthogonal
     symplectic transformation matrix

           (  U1  U2  )
       U = (          )
           ( -U2  U1  )

     will be computed.
Specification
      SUBROUTINE DGHFST( JOB, COMPQ, COMPU, N, Z, LDZ, B, LDB, FG, LDFG,
     $                   Q, LDQ, U1, LDU1, U2, LDU2, ALPHAR, ALPHAI,
     $                   BETA, IWORK, LIWORK, DWORK, LDWORK, INFO )
C
C     .. Scalar Arguments ..
      CHARACTER          COMPQ, COMPU, JOB
      INTEGER            INFO, LDB, LDFG, LDQ, LDU1, LDU2, LDWORK, LDZ,
     $                   LIWORK, N
C
C     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   ALPHAI( * ), ALPHAR( * ), B( LDB, * ),
     $                   BETA( * ), DWORK( * ), FG( LDFG, * ),
     $                   Q( LDQ, * ), U1( LDU1, * ), U2( LDU2, * ),
     $                   Z( LDZ, * )
Arguments

Mode Parameters

     JOB     CHARACTER*1
             Specifies the computation to be performed, as follows:
             = 'E':  compute the eigenvalues only; Z and T will not
                     necessarily be put into the forms in (2);
             = 'T':  put Z and T into the forms in (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.

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

     Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
             On entry, the leading N-by-N part of this array must
             contain the matrix Z.
             On exit, if JOB = 'T', the leading N-by-N part of this
             array contains the matrix Zout; otherwise, it contains the
             matrix Z just before the application of the periodic QZ
             algorithm. The entries in the rows N/2+1 to N and the
             first N/2 columns are unchanged.

     LDZ     INTEGER
             The leading dimension of the array Z.  LDZ >= MAX(1, N).

     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 the
             matrix B just before the application of the periodic QZ
             algorithm.

     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.
             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 entries on the diagonal and the first superdiagonal of
             this array are not referenced, but are assumed to be zero.
             Moreover, the diagonal and the first subdiagonal of this
             array on exit coincide to the corresponding diagonals of
             this array on entry.

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

     Q       (output) DOUBLE PRECISION array, dimension (LDQ, N)
             On exit, if COMPQ = 'I', the leading N-by-N part of this
             array contains the orthogonal transformation matrix Q.
             On exit, if COMPQ = 'N', the leading N-by-N part of this
             array contains the orthogonal matrix Q1, such that

                      (  Z11  Z12  )
               Z*Q1 = (            ),
                      (   0   Z22  )

             where Z11 and Z22' are upper triangular (the first step
             of the algorithm).

     LDQ     INTEGER
             The leading dimension of the array Q.  LDQ >= MAX(1, N).

     U1      (input/output) DOUBLE PRECISION array, dimension
                            (LDU1, N/2)
             On entry, if COMPU = 'U', then the leading N/2-by-N/2 part
             of this array must contain the upper left block of a
             given matrix U0, and on exit, the leading N/2-by-N/2 part
             of this array contains the updated upper left block U1 of
             the product of the input matrix U0 and the transformation
             matrix U used to transform the matrices Z and T.
             On exit, if COMPU = 'I', then the leading N/2-by-N/2 part
             of this array contains the upper left block U1 of the
             orthogonal symplectic transformation matrix U.
             If COMPU = 'N' this array is not referenced.

     LDU1    INTEGER
             The leading dimension of the array U1.
             LDU1 >= 1,           if COMPU = 'N';
             LDU1 >= MAX(1, N/2), if COMPU = 'I' or COMPU = 'U'.

     U2      (input/output) DOUBLE PRECISION array, dimension
                            (LDU2, N/2)
             On entry, if COMPU = 'U', then the leading N/2-by-N/2 part
             of this array must contain the upper right block of a
             given matrix U0, and on exit, the leading N/2-by-N/2 part
             of this array contains the updated upper right block U2 of
             the product of the input matrix U0 and the transformation
             matrix U used to transform the matrices Z and T.
             On exit, if COMPU = 'I', then the leading N/2-by-N/2 part
             of this array contains the upper right block U2 of the
             orthogonal symplectic transformation matrix U.
             If COMPU = 'N' this array is not referenced.

     LDU2    INTEGER
             The leading dimension of the array U2.
             LDU2 >= 1,           if COMPU = 'N';
             LDU2 >= MAX(1, N/2), if COMPU = 'I' or COMPU = '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, with ALPHAI(j+1) = -ALPHAI(j).

     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
     IWORK   INTEGER array, dimension (LIWORK)

     LIWORK  INTEGER
             The dimension of the array IWORK.  LIWORK >= N+9.

     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
             On exit, if INFO = 0, DWORK(1) returns the optimal LDWORK.
             On exit, if INFO = -23, DWORK(1) returns the minimum
             value of LDWORK.

     LDWORK  INTEGER
             The dimension of the array DWORK.
             If JOB = 'E' and COMPQ = 'N' and COMPU = 'N',
                   LDWORK >= 3/4*N**2+MAX(N, 24)+3;
             else, LDWORK >= 3/2*N**2+MAX(N, 24)+3.
             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: problem during computation of the eigenvalues;
             = 2: periodic QZ algorithm did not converge in the SLICOT
                  Library routine MB03BD.
Method
     The algorithm uses Givens rotations and Householder reflections to
     annihilate elements in Z and T such that Z is in a special block
     triangular form and T is in skew-Hamiltonian Hessenberg form:

         (  Z11  Z12  )      (  B1  F1  )                  
     Z = (            ), T = (          ),
         (   0   Z22  )      (   0  B1' )

     with Z11 and Z22' upper triangular and B1 upper Hessenberg.
     Subsequently, the periodic QZ algorithm is applied to the pencil
     aZ22' Z11 - bB1 to determine orthogonal matrices Q1, Q2 and U such
     that U' Z11 Q1, Q2' Z22' U are upper triangular and Q2' B1 Q1 is
     upper quasi-triangular. See also page 35 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

*     DGHFST EXAMPLE PROGRAM TEXT
*
*     .. Parameters ..
      INTEGER            NIN, NOUT
      PARAMETER          ( NIN = 5, NOUT = 6 )
      INTEGER            NMAX
      PARAMETER          ( NMAX = 60 )
      INTEGER            LDB, LDFG, LDQ, LDU1, LDU2, LDWORK, LDZ,
     $                   LIWORK
      PARAMETER          ( LDB = NMAX/2, LDFG = NMAX/2,
     $                     LDQ = NMAX,   LDU1 = NMAX/2, LDU2 = NMAX/2,
     $                     LDWORK = 3*NMAX**2/2 + MAX( NMAX, 24 ),
     $                     LDZ = NMAX, LIWORK = NMAX/2 + 9 )
*
*     .. Local Scalars ..
      CHARACTER          COMPQ, COMPU, JOB
      INTEGER            I, INFO, J, N
*
*     .. Local Arrays ..
      INTEGER            IWORK( LIWORK )
      DOUBLE PRECISION   ALPHAI( NMAX/2 ),   ALPHAR( NMAX/2 ),
     $                   B( LDB, NMAX/2 ),     BETA( NMAX/2 ),
     $                   DWORK( LDWORK ),  FG( LDFG, NMAX/2+1 ),
     $                   Q( LDQ, NMAX ),   U1( LDU1, NMAX/2 ),
     $                   U2( LDU2, NMAX/2 ), Z( LDZ, NMAX )
*
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*
*     .. External Subroutines ..
      EXTERNAL           DGHFST
*
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, 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, COMPU, N
      READ( NIN, FMT = * ) ( (  Z( I, J ), J = 1, N ),     I = 1, N )
      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( COMPU, 'U' ) ) THEN
         READ( NIN, FMT = * ) ( ( U1( I, J ), J = 1, N/2 ), I = 1, N/2 )
         READ( NIN, FMT = * ) ( ( U2( I, J ), J = 1, N/2 ), I = 1, N/2 )
      END IF
      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 (factored version).
         CALL DGHFST( JOB, COMPQ, COMPU, N, Z, LDZ, B, LDB, FG, LDFG, Q,
     $                LDQ, U1, LDU1, U2, LDU2, ALPHAR, ALPHAI, BETA,
     $                IWORK, LIWORK, DWORK, LDWORK, INFO )
         IF( INFO.NE.0 ) THEN
            WRITE( NOUT, FMT = 99997 ) INFO
         ELSE
            WRITE( NOUT, FMT = 99996 )
            DO 10 I = 1, N
               WRITE( NOUT, FMT = 99995 ) ( Z( I, J ), J = 1, N )
   10       CONTINUE
            WRITE( NOUT, FMT = 99994 )
            DO 20 I = 1, N/2
               WRITE( NOUT, FMT = 99995 ) ( B( I, J ), J = 1, N/2 )
   20       CONTINUE
            WRITE( NOUT, FMT = 99993 )
            DO 30 I = 1, N/2
               WRITE( NOUT, FMT = 99995 ) ( FG( I, J ), J = 1, N/2+1 )
   30       CONTINUE
            WRITE( NOUT, FMT = 99992 )
            WRITE( NOUT, FMT = 99995 ) ( ALPHAR( I ), I = 1, N/2 )
            WRITE( NOUT, FMT = 99991 )
            WRITE( NOUT, FMT = 99995 ) ( ALPHAI( I ), I = 1, N/2 )
            WRITE( NOUT, FMT = 99990 )
            WRITE( NOUT, FMT = 99995 ) (   BETA( I ), I = 1, N/2 )
            WRITE( NOUT, FMT = 99989 )
            DO 40 I = 1, N
               WRITE( NOUT, FMT = 99995 ) ( Q( I, J ), J = 1, N )
   40       CONTINUE
            IF ( .NOT.LSAME( COMPU, 'N' ) ) THEN
               WRITE( NOUT, FMT = 99988 )
               DO 50 I = 1, N/2
                  WRITE( NOUT, FMT = 99995 ) ( U1( I, J ), J = 1, N/2 )
   50          CONTINUE
               WRITE( NOUT, FMT = 99987 )
               DO 60 I = 1, N/2
                  WRITE( NOUT, FMT = 99995 ) ( U2( I, J ), J = 1, N/2 )
   60          CONTINUE
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT( 'DGHFST EXAMPLE PROGRAM RESULTS', 1X )
99998 FORMAT( 'N is out of range.', /, 'N = ', I5 )
99997 FORMAT( 'INFO on exit from DGHFST = ', I2 )
99996 FORMAT( 'The transformed matrix Z is ' )
99995 FORMAT( 60( 1X, F8.4 ) )
99994 FORMAT( 'The transformed matrix B is ' )
99993 FORMAT( 'The transformed matrix FG is ' )
99992 FORMAT( 'The vector ALPHAR is ' )
99991 FORMAT( 'The vector ALPHAI is ' )
99990 FORMAT( 'The vector BETA is ' )
99989 FORMAT( 'The matrix Q is ' )
99988 FORMAT( 'The upper left block of the matrix U is ' )
99987 FORMAT( 'The upper right block of the matrix U is ' )
      END
Program Data
DGHFST EXAMPLE PROGRAM DATA
   T   I   I   8
   0.0949   3.3613  -4.7663  -0.5534   0.6408  -3.2793   3.4253   2.9654
   0.1138  -1.5903   2.1837  -4.1648  -4.3775  -1.7454   0.1744   2.3262
   2.7505   4.4048   4.4183   3.0478   2.7728   2.3048  -0.6451  -1.2045
   3.6091  -4.1716   3.4461   3.6880  -0.0985   3.8458   0.2528  -1.3859
   0.4352  -3.2829   3.7246   0.4794  -0.3690  -1.5562  -3.4817  -2.2902
   1.3080  -3.9881  -3.5497   3.5020   2.2582   4.4764  -4.4080  -1.6818
   1.1308  -1.5087   2.4730   2.1553  -1.7129  -4.8669  -2.4102   4.2274
   4.7933  -4.3671  -0.0473  -2.0092   1.2439  -4.7385   3.4242  -0.2764
   2.0936   1.5510   4.5974   2.5127
   2.5469  -3.3739  -1.5961  -2.4490
  -2.2397  -3.8100   0.8527   0.0596
   1.7970  -0.0164  -2.7619   1.9908
   1.0000   2.0000  -4.0500   1.3353   0.2899
  -0.4318   2.0000   2.0000  -2.9860  -0.0160
   1.0241   0.9469   2.0000   2.0000   1.3303
   0.0946  -0.1272  -4.4003   2.0000   2.0000
Program Results
DGHFST EXAMPLE PROGRAM RESULTS
The transformed matrix Z is 
  -2.5678  -2.9798   0.4885  -2.8719   2.7331   1.2619   1.7893   2.8246
   0.0000  -3.9711  -6.1190   6.1945  -3.0512  -5.3985  -1.4398   4.0102
   0.0000   0.0000   4.3224   4.5967   0.5414  -4.6471   3.6510  -0.5421
   0.0000   0.0000   0.0000   7.0155   1.5557   2.0499   3.7184  -2.3864
   0.4352  -3.2829   3.7246   0.4794  -5.3205   0.0000   0.0000   0.0000
   1.3080  -3.9881  -3.5497   3.5020   2.2836   6.9817   0.0000   0.0000
   1.1308  -1.5087   2.4730   2.1553  -1.6709  -0.8715   8.1253   0.0000
   4.7933  -4.3671  -0.0473  -2.0092  -3.9547   0.2398   1.0447   5.5977
The transformed matrix B is 
   3.8629  -1.3239   0.1512   2.1882
   0.0000   3.7484  -3.5999  -2.4361
   0.0000  -3.6088  -2.5273  -0.9003
   0.0000   0.0000   0.0000  -6.7384
The transformed matrix FG is 
   1.0000   2.0000  -0.7130  -1.2545  -1.3653
   0.0158   2.0000   2.0000   3.4030   3.2438
  -1.1665   2.5791   2.0000   2.0000  -0.3270
   3.3823  -1.2344   3.9016   2.0000   2.0000
The vector ALPHAR is 
   1.1310  -0.0697  -0.0697  -0.6864
The vector ALPHAI is 
   0.0000   0.6035  -0.6035   0.0000
The vector BETA is 
   4.0000   4.0000   4.0000   4.0000
The matrix Q is 
  -0.6042  -0.4230  -0.4661   0.1400  -0.2947   0.3486  -0.0892   0.0534
  -0.3706   0.1454   0.4415  -0.1381  -0.1210   0.2727   0.7320  -0.0524
   0.1325  -0.2745  -0.0461  -0.5084  -0.3163  -0.2896   0.1565   0.6619
   0.2373   0.5415  -0.5094   0.3373  -0.3852  -0.0092   0.3328   0.1339
   0.4777  -0.4463   0.2826   0.5172  -0.0775   0.3845   0.1186   0.2395
  -0.0116  -0.4407  -0.1758   0.2474   0.1236  -0.6172   0.4616  -0.3228
   0.1237  -0.0394  -0.4293  -0.2090   0.7209   0.3333   0.2984   0.1883
   0.4245  -0.1906  -0.1766  -0.4655  -0.3304   0.2832   0.0695  -0.5843
The upper left block of the matrix U is 
   0.0154  -0.4941  -0.5382   0.6826
   0.4829  -0.1455  -0.2953  -0.3491
   0.4981   0.1510   0.1052   0.1810
  -0.0188   0.6384  -0.5122   0.0587
The upper right block of the matrix U is 
   0.0000   0.0000   0.0000   0.0000
   0.3179  -0.4272  -0.1896  -0.4659
   0.5644   0.0775   0.4498   0.3979
  -0.3137  -0.3404   0.3339   0.0239

Return to index