**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.

SUBROUTINE MB03ID( 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 .. Scalar Arguments .. CHARACTER COMPQ, COMPU INTEGER INFO, LDA, LDB, LDC, LDD, LDF, LDQ, LDU1, LDU2, $ LDWORK, LIWORK, N, NEIG C .. Array Arguments .. INTEGER IWORK( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ), $ D( LDD, * ), DWORK( * ), F( LDF, * ), $ Q( LDQ, * ), U1( LDU1, * ), U2( LDU2, * )

**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.

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.

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).

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 MB03CD; = 3: an error occured during the execution of MB03GD.

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 MB03CD and MB03GD are called for the corresponding matrix structures.

[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.

3 The algorithm is numerically backward stable and needs O(N ) real floating point operations.

None

**Program Text**

None

None

None