MB03GD

Exchanging eigenvalues of a real 2-by-2 or 4-by-4 block upper triangular skew-Hamiltonian/Hamiltonian pencil (factored version)

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