DGHFXC

Moving eigenvalues with negative real parts of a real skew-Hamiltonian/Hamiltonian pencil in
structured Schur form to the leading subpencil (factored version)

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

Purpose

     To move the eigenvalues with strictly negative real parts of an
     N-by-N real skew-Hamiltonian/Hamiltonian pencil aS - bH in
     structured Schur form, with

                          (  0  I  )      (  A  D  )      (  B  F  )
       S = J Z' J' Z, J = (        ), Z = (        ), H = (        ),
                          ( -I  0  )      (  0  C  )      (  0 -B' )

     to the leading principal subpencil, while keeping the triangular
     form. Above, A is upper triangular, B is upper quasi-triangular,
     and C is lower triangular.
     The matrices Z and H are transformed by an orthogonal symplectic
     matrix U and an orthogonal matrix Q such that

                       (  Aout  Dout  )
       Zout = U' Z Q = (              ), and
                       (    0   Cout  )
                                                                    (1)
                            (  Bout  Fout  )
       Hout = J Q' J' H Q = (              ),
                            (    0  -Bout' )

     where Aout, Bout and Cout remain in triangular form. The notation
     M' denotes the transpose of the matrix M.
     Optionally, if COMPQ = 'I' or COMPQ = 'U', the orthogonal matrix Q
     that fulfills (1) is computed.
     Optionally, if COMPU = 'I' or COMPU = 'U', the orthogonal
     symplectic matrix

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

     that fulfills (1) is computed.
Specification
      SUBROUTINE DGHFXC( COMPQ, COMPU, N, A, LDA, C, LDC, D, LDD, B,
     $                   LDB, F, LDF, Q, LDQ, U1, LDU1, U2, LDU2, NEIG,
     $                   IWORK, LIWORK, DWORK, LDWORK, INFO )
C
C     .. Scalar Arguments ..
      CHARACTER          COMPQ, COMPU
      INTEGER            INFO, LDA, LDB, LDC, LDD, LDF, LDQ, LDU1, LDU2,
     $                   LDWORK, LIWORK, N, NEIG
C
C     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
     $                   D( LDD, * ), DWORK( *  ), F( LDF, * ),
     $                   Q( LDQ, * ), U1( LDU1, * ), U2( LDU2, * )
Arguments

Mode Parameters

     COMPQ   CHARACTER*1
             Specifies whether or not the orthogonal transformations
             should be accumulated in the array 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 matrix Q0*Q is returned, where Q
                     is the product of the orthogonal transformations
                     that are applied to the pencil aS - bH to reorder
                     the eigenvalues.

     COMPU   CHARACTER*1
             Specifies whether or not the orthogonal symplectic
             transformations should be accumulated in the arrays U1 and
             U2, as follows:
             = 'N':  U1 and U2 are not computed;
             = 'I':  the arrays U1 and U2 are initialized internally,
                     and the submatrices U1 and U2 defining the
                     orthogonal symplectic matrix U are 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 - bH to reorder the eigenvalues.             
Input/Output Parameters
     N       (input) INTEGER
             The order of the pencil aS - bH.  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 upper triangular matrix A. The elements of the
             strictly lower triangular part of this array are not used.
             On exit, the leading  N/2-by-N/2 part of this array
             contains the transformed matrix Aout.

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

     C       (input/output) DOUBLE PRECISION array, dimension
                            (LDC, N/2)
             On entry, the leading N/2-by-N/2 part of this array must
             contain the lower triangular matrix C. The elements of the
             strictly upper triangular part of this array are not used.
             On exit, the leading  N/2-by-N/2 part of this array
             contains the transformed matrix Cout.

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

     D       (input/output) DOUBLE PRECISION array, dimension
                            (LDD, N/2)
             On entry, the leading N/2-by-N/2 part of this array must
             contain the matrix D.
             On exit, the leading  N/2-by-N/2 part of this array
             contains the transformed matrix Dout.

     LDD     INTEGER
             The leading dimension of the array D.  LDD >= 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 upper quasi-triangular matrix B.
             On exit, the leading  N/2-by-N/2 part of this array
             contains the transformed upper quasi-triangular part of
             the matrix Bout.
             The part below the first subdiagonal of this array is
             not referenced.

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

     F       (input/output) DOUBLE PRECISION array, dimension
                            (LDF, N/2)
             On entry, the leading N/2-by-N/2 part of this array must
             contain the upper triangular part of the symmetric matrix
             F.
             On exit, the leading  N/2-by-N/2 part of this array
             contains the transformed upper triangular part of the
             matrix Fout.
             The strictly lower triangular part of this array is not
             referenced, except for the element F(N/2,N/2-1), but its
             initial value is preserved.

     LDF     INTEGER
             The leading dimension of the array F.  LDF >= 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 Z and H.
             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'.

     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 H.
             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 H.
             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'.

     NEIG    (output) INTEGER
             The number of eigenvalues in aS - bH with strictly
             negative real part.
Workspace
     IWORK   INTEGER array, dimension (LIWORK)

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

     DWORK   DOUBLE PRECISION array, dimension (LDWORK)

     LDWORK  INTEGER
             The dimension of the array DWORK.
             If COMPQ = 'N',
                LDWORK >= MAX(2*N+48,171);
             if COMPQ = 'I' or COMPQ = 'U',
                LDWORK >= MAX(4*N+48,171).
Error Indicator
     INFO    INTEGER
             = 0: succesful exit;
             < 0: if INFO = -i, the i-th argument had an illegal value;
             = 1: the periodic QZ algorithm did not converge in SLICOT
                  Library routine MB03BB;
             = 2: an error occured during the execution of DGHFEX;
             = 3: an error occured during the execution of DGHFEY.
Method
     The algorithm reorders the eigenvalues like the following scheme:

     Step 1: Reorder the eigenvalues in the subpencil aC'*A - bB.
          I. Reorder the eigenvalues with negative real parts to the
             top.
         II. Reorder the eigenvalues with positive real parts to the
             bottom.

     Step 2: Reorder the remaining eigenvalues with negative real
             parts in the pencil aS - bH.
          I. Exchange the eigenvalues between the last diagonal block
             in aC'*A - bB and the last diagonal block in aS - bH.
         II. Move the eigenvalues of the R-th block to the (MM+1)-th
             block, where R denotes the number of upper quasi-
             triangular blocks in aC'*A - bB and MM denotes the current
             number of blocks in aC'*A - bB with eigenvalues with
             negative real parts.

     The algorithm uses a sequence of orthogonal transformations as
     described on page 25 in [1]. To achieve those transformations the
     elementary subroutines DGHFEX and DGHFEY are called for the
     corresponding matrix structures.
References
     [1] Benner, P., Byers, R., Losse, P., Mehrmann, V. and Xu, H.
         Numerical Solution of Real Skew-Hamiltonian/Hamiltonian
         Eigenproblems.
         Tech. Rep., Technical University Chemnitz, Germany,
         Nov. 2007.
Numerical Aspects
     
                                                               3
     The algorithm is numerically backward stable and needs O(N ) real
     floating point operations.

Further Comments
   
     None.
Example

Program Text

     None.
Program Data
     None.   
Program Results
     None.

Return to index