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