**Purpose**

To compute an orthogonal matrix Q and an orthogonal symplectic matrix U for a real regular 2-by-2 or 4-by-4 skew-Hamiltonian/ Hamiltonian pencil a J B' J' B - b D with ( B11 B12 ) ( D11 D12 ) ( 0 I ) B = ( ), D = ( ), J = ( ), ( 0 B22 ) ( 0 -D11' ) ( -I 0 ) such that J Q' J' D Q and U' B Q keep block triangular form, but the eigenvalues are reordered. The notation M' denotes the transpose of the matrix M.

SUBROUTINE MB03GD( N, B, LDB, D, LDD, MACPAR, Q, LDQ, U, LDU, $ DWORK, LDWORK, INFO ) C .. Scalar Arguments .. INTEGER INFO, LDB, LDD, LDQ, LDU, LDWORK, N C .. Array Arguments .. DOUBLE PRECISION B( LDB, * ), D( LDD, * ), DWORK( * ), $ MACPAR( * ), Q( LDQ, * ), U( LDU, * )

**Input/Output Parameters**

N (input) INTEGER The order of the pencil a J B' J' B - b D. N = 2 or N = 4. B (input) DOUBLE PRECISION array, dimension (LDB, N) The leading N-by-N part of this array must contain the non-trivial factor of the decomposition of the skew-Hamiltonian input matrix J B' J' B. The (2,1) block is not referenced. LDB INTEGER The leading dimension of the array B. LDB >= N. D (input) DOUBLE PRECISION array, dimension (LDD, N) The leading N/2-by-N part of this array must contain the first block row of the second matrix of a J B' J' B - b D. The matrix D has to be Hamiltonian. The strict lower triangle of the (1,2) block is not referenced. LDD INTEGER The leading dimension of the array D. LDD >= N/2. MACPAR (input) DOUBLE PRECISION array, dimension (2) Machine parameters: MACPAR(1) (machine precision)*base, DLAMCH( 'P' ); MACPAR(2) safe minimum, DLAMCH( 'S' ). This argument is not used for N = 2. Q (output) DOUBLE PRECISION array, dimension (LDQ, N) The leading N-by-N part of this array contains the orthogonal transformation matrix Q. LDQ INTEGER The leading dimension of the array Q. LDQ >= N. U (output) DOUBLE PRECISION array, dimension (LDU, N) The leading N-by-N part of this array contains the orthogonal symplectic transformation matrix U. LDU INTEGER The leading dimension of the array U. LDU >= N.

DWORK DOUBLE PRECISION array, dimension (LDWORK) If N = 2 then DWORK is not referenced. LDWORK INTEGER The length of the array DWORK. If N = 2 then LDWORK >= 0; if N = 4 then LDWORK >= 12.

INFO INTEGER = 0: succesful exit; = 1: B11 or B22 is a (numerically) singular matrix.

The algorithm uses orthogonal transformations as described on page 22 in [1], but with an improved implementation.

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

The algorithm is numerically backward stable.

None

**Program Text**

None

None

None