MB04SU

Symplectic QR decomposition of a real 2M-by-N matrix

Purpose

  To compute a symplectic QR decomposition of a real 2M-by-N matrix
  [A; B],

            [ A ]             [ R11  R12 ]
            [   ] = Q * R = Q [          ],
            [ B ]             [ R21  R22 ]

  where Q is a symplectic orthogonal matrix, R11 is upper triangular
  and R21 is strictly upper triangular.
  If [A; B] is symplectic then, theoretically, R21 = 0 and
  R22 = inv(R11)^T. Unblocked version.

      SUBROUTINE MB04SU( M, N, A, LDA, B, LDB, CS, TAU, DWORK, LDWORK,
     $                   INFO )
C     .. Scalar Arguments ..
      INTEGER           INFO, LDA, LDB, LDWORK, M, N
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), CS(*), DWORK(*), TAU(*)

Input/Output Parameters

  M       (input) INTEGER
          The number of rows of A and B. M >= 0.

  N       (input) INTEGER
          The number of columns of A and B. N >= 0.

  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the leading M-by-N part of this array must
          contain the matrix A.
          On exit, the leading M-by-N part of this array contains
          the matrix [ R11  R12 ] and, in the zero parts of R,
          information about the elementary reflectors used to
          compute the symplectic QR decomposition.

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

  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
          On entry, the leading M-by-N part of this array must
          contain the matrix B.
          On exit, the leading M-by-N part of this array contains
          the matrix [ R21  R22 ] and, in the zero parts of B,
          information about the elementary reflectors used to
          compute the symplectic QR decomposition.

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

  CS      (output) DOUBLE PRECISION array, dimension (2 * min(M,N))
          On exit, the first 2*min(M,N) elements of this array
          contain the cosines and sines of the symplectic Givens
          rotations used to compute the symplectic QR decomposition.

  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
          On exit, the first min(M,N) elements of this array
          contain the scalar factors of some of the elementary
          reflectors.

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

  LDWORK  INTEGER
          The length of the array DWORK.  LDWORK >= MAX(1,N).

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value.

  The matrix Q is represented as a product of symplectic reflectors
  and Givens rotations

  Q = diag( H(1),H(1) ) G(1) diag( F(1),F(1) )
      diag( H(2),H(2) ) G(2) diag( F(2),F(2) )
                        ....
      diag( H(k),H(k) ) G(k) diag( F(k),F(k) ),

  where k = min(m,n).

  Each H(i) has the form

        H(i) = I - tau * w * w'

  where tau is a real scalar, and w is a real vector with
  w(1:i-1) = 0 and w(i) = 1; w(i+1:m) is stored on exit in
  B(i+1:m,i), and tau in B(i,i).

  Each F(i) has the form

        F(i) = I - nu * v * v'

  where nu is a real scalar, and v is a real vector with
  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
  A(i+1:m,i), and nu in TAU(i).

  Each G(i) is a Givens rotation acting on rows i of A and B,
  where the cosine is stored in CS(2*i-1) and the sine in
  CS(2*i).

  [1] Bunse-Gerstner, A.
      Matrix factorizations for symplectic QR-like methods.
      Linear Algebra Appl., 83, pp. 49-77, 1986.

  [2] Byers, R.
      Hamiltonian and Symplectic Algorithms for the Algebraic
      Riccati Equation.
      Ph.D. Dissertation, Center for Applied Mathematics,
      Cornell University, Ithaca, NY, 1983.

  The algorithm requires
     8*M*N*N - 8/3*N*N*N +  2*M*N + 6*N*N + 8/3*N,  if M >= N,
     8*M*M*N - 8/3*M*M*M + 14*M*N - 6*M*M + 8/3*N,  if M <= N,
  floating point operations and is numerically backward stable.

  None

Program Text

  None

  None

  None