## MB04SU

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

[Specification] [Arguments] [Method] [References] [Comments] [Example]

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.

```
Specification
```      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(*)

```
Arguments

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.

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

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

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

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

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

```
Numerical Aspects
```  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
```
Example

Program Text

```  None
```
Program Data
```  None
```
Program Results
```  None
```