## MB04ZD

### Transforming a Hamiltonian matrix into a square-reduced Hamiltonian matrix

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

Purpose

```  To transform a Hamiltonian matrix

( A   G  )
H = (      T )                                           (1)
( Q  -A  )

into a square-reduced Hamiltonian matrix

( A'  G'  )
H' = (       T )                                         (2)
( Q' -A'  )
T
by an orthogonal symplectic similarity transformation H' = U H U,
where
(  U1   U2 )
U = (          ).                                        (3)
( -U2   U1 )
T
The square-reduced Hamiltonian matrix satisfies Q'A' - A' Q' = 0,
and

2       T     2     ( A''   G''  )
H'  :=  (U  H U)   =  (          T ).
( 0     A''  )

In addition, A'' is upper Hessenberg and G'' is skew symmetric.
The square roots of the eigenvalues of A'' = A'*A' + G'*Q' are the
eigenvalues of H.

```
Specification
```      SUBROUTINE MB04ZD( COMPU, N, A, LDA, QG, LDQG, U, LDU, DWORK, INFO
\$                 )
C     .. Scalar Arguments ..
INTEGER           INFO, LDA, LDQG, LDU, N
CHARACTER         COMPU
C     .. Array Arguments ..
DOUBLE PRECISION  A(LDA,*), DWORK(*), QG(LDQG,*), U(LDU,*)

```
Arguments

Mode Parameters

```  COMPU   CHARACTER*1
Indicates whether the orthogonal symplectic similarity
transformation matrix U in (3) is returned or
accumulated into an orthogonal symplectic matrix, or if
the transformation matrix is not required, as follows:
= 'N':         U is not required;
= 'I' or 'F':  on entry, U need not be set;
on exit, U contains the orthogonal
symplectic matrix U from (3);
= 'V' or 'A':  the orthogonal symplectic similarity
transformations are accumulated into U;
on input, U must contain an orthogonal
symplectic matrix S;
on exit, U contains S*U with U from (3).
See the description of U below for details.

```
Input/Output Parameters
```  N       (input) INTEGER
The order of the matrices A, G, and Q.  N >= 0.

A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On input, the leading N-by-N part of this array must
contain the upper left block A of the Hamiltonian matrix H
in (1).
On output, the leading N-by-N part of this array contains
the upper left block A' of the square-reduced Hamiltonian
matrix H' in (2).

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

QG      (input/output) DOUBLE PRECISION array, dimension
(LDQG,N+1)
On input, the leading N-by-N lower triangular part of this
array must contain the lower triangle of the lower left
symmetric block Q of the Hamiltonian matrix H in (1), and
the N-by-N upper triangular part of the submatrix in the
columns 2 to N+1 of this array must contain the upper
triangle of the upper right symmetric block G of H in (1).
So, if i >= j, then Q(i,j) = Q(j,i) is stored in QG(i,j)
and G(i,j) = G(j,i) is stored in QG(j,i+1).
On output, the leading N-by-N lower triangular part of
this array contains the lower triangle of the lower left
symmetric block Q', and the N-by-N upper triangular part
of the submatrix in the columns 2 to N+1 of this array
contains the upper triangle of the upper right symmetric
block G' of the square-reduced Hamiltonian matrix H'
in (2).

LDQG    INTEGER
The leading dimension of the array QG.  LDQG >= MAX(1,N).

U       (input/output) DOUBLE PRECISION array, dimension (LDU,2*N)
If COMPU = 'N', then this array is not referenced.
If COMPU = 'I' or 'F', then the input contents of this
array are not specified.  On output, the leading
N-by-(2*N) part of this array contains the first N rows
of the orthogonal symplectic matrix U in (3).
If COMPU = 'V' or 'A', then, on input, the leading
N-by-(2*N) part of this array must contain the first N
rows of an orthogonal symplectic matrix S. On output, the
leading N-by-(2*N) part of this array contains the first N
rows of the product S*U where U is the orthogonal
symplectic matrix from (3).
The storage scheme implied by (3) is used for orthogonal
symplectic matrices, i.e., only the first N rows are
stored, as they contain all relevant information.

LDU     INTEGER
The leading dimension of the array U.
LDU >= MAX(1,N), if COMPU <> 'N';
LDU >= 1,        if COMPU =  'N'.

```
Workspace
```  DWORK   DOUBLE PRECISION array, dimension (2*N)

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

```
Method
```  The Hamiltonian matrix H is transformed into a square-reduced
Hamiltonian matrix H' using the implicit version of Van Loan's
method as proposed in [1,2,3].

```
References
```   Van Loan, C. F.
A Symplectic Method for Approximating All the Eigenvalues of
a Hamiltonian Matrix.
Linear Algebra and its Applications, 61, pp. 233-251, 1984.

 Byers, R.
Hamiltonian and Symplectic Algorithms for the Algebraic
Riccati Equation.
Ph. D. Thesis, Cornell University, Ithaca, NY, January 1983.

 Benner, P., Byers, R., and Barth, E.
Fortran 77 Subroutines for Computing the Eigenvalues of
Hamiltonian Matrices. I: The Square-Reduced Method.
ACM Trans. Math. Software, 26, 1, pp. 49-77, 2000.

```
Numerical Aspects
```  This algorithm requires approximately 20*N**3 flops for
transforming H into square-reduced form. If the transformations
are required, this adds another 8*N**3 flops. The method is
strongly backward stable in the sense that if H' and U are the
computed square-reduced Hamiltonian and computed orthogonal
symplectic similarity transformation, then there is an orthogonal
symplectic matrix T and a Hamiltonian matrix M such that

H T  =  T M

|| T - U ||   <=  c1 * eps

|| H' - M ||  <=  c2 * eps * || H ||

where c1, c2 are modest constants depending on the dimension N and
eps is the machine precision.

Eigenvalues computed by explicitly forming the upper Hessenberg
matrix  A'' = A'A' + G'Q', with A', G', and Q' as in (2), and
applying the Hessenberg QR iteration to A'' are exactly
eigenvalues of a perturbed Hamiltonian matrix H + E,  where

|| E ||  <=  c3 * sqrt(eps) * || H ||,

and c3 is a modest constant depending on the dimension N and eps
is the machine precision.  Moreover, if the norm of H and an
eigenvalue lambda are of roughly the same magnitude, the computed
eigenvalue is essentially as accurate as the computed eigenvalue
from traditional methods.  See  or .

```
```  None
```
Example

Program Text

```*     MB04ZD EXAMPLE PROGRAM TEXT.
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
INTEGER          NIN, NOUT
PARAMETER        ( NIN = 5, NOUT = 6 )
INTEGER          NMAX
PARAMETER        ( NMAX = 20 )
INTEGER          LDA, LDQG, LDU
PARAMETER        ( LDA = NMAX, LDQG = NMAX, LDU = NMAX )
INTEGER          LDWORK
PARAMETER        ( LDWORK = ( NMAX+NMAX )*( NMAX+NMAX+1 ) )
DOUBLE PRECISION ZERO, ONE
PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0 )
*     .. Local Scalars ..
INTEGER          I, INFO, IJ, J, JI, N, POS, WPOS
CHARACTER*1      COMPU
*     .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), DWORK(LDWORK), QG(LDQG,NMAX+1),
\$                 U(LDU,NMAX)
*     .. External Subroutines ..
EXTERNAL         DCOPY, DGEMM, DSYMV, MB04ZD
*     .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, COMPU
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99998 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J),    J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( QG(J,I+1), I = J,N ), J = 1,N )
READ ( NIN, FMT = * ) ( ( QG(I,J),   I = J,N ), J = 1,N )
*        Square-reduce by symplectic orthogonal similarity.
CALL MB04ZD( COMPU, N, A, LDA, QG, LDQG, U, LDU, DWORK, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99997 ) INFO
ELSE
*           Show the square-reduced Hamiltonian.
WRITE ( NOUT, FMT = 99996 )
DO 10 I = 1, N
WRITE ( NOUT, FMT = 99994 )  ( A(I,J),    J = 1,N ),
\$            ( QG(J,I+1), J = 1,I-1 ), ( QG(I,J+1), J = I,N )
10          CONTINUE
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99994 ) ( QG(I,J), J = 1,I-1 ),
\$               ( QG(J,I), J = I,N ), ( -A(J,I), J = 1,N )
20          CONTINUE
*           Show the square of H.
WRITE ( NOUT, FMT = 99995 )
WPOS = ( NMAX+NMAX )*( NMAX+NMAX )
*                                                    T
*           Compute N11 = A*A + G*Q and set N22 = N11 .
CALL DGEMM( 'N', 'N', N, N, N, ONE, A, LDA, A, LDA, ZERO,
\$                  DWORK, N+N )
DO 30 I = 1, N
CALL DCOPY( N-I+1, QG(I,I), 1, DWORK(WPOS+I), 1 )
CALL DCOPY( I-1, QG(I,1), LDQG, DWORK(WPOS+1), 1 )
CALL DSYMV( 'U', N, ONE, QG(1,2), LDQG, DWORK(WPOS+1), 1,
\$                     ONE, DWORK((I-1)*(N+N)+1), 1 )
POS = N*( N+N ) + N + I
CALL DCOPY( N, DWORK((I-1)*(N+N)+1), 1, DWORK(POS), N+N )
30          CONTINUE
DO 40 I = 1, N
CALL DSYMV( 'U', N, -ONE, QG(1,2), LDQG, A(I,1), LDA,
\$                     ZERO, DWORK((N+I-1)*(N+N)+1), 1 )
CALL DSYMV( 'L', N, ONE, QG, LDQG, A(1,I), 1, ZERO,
\$                     DWORK((I-1)*(N+N)+N+1), 1 )
40          CONTINUE
DO 60 J = 1, N
DO 50 I = J, N
IJ = ( N+J-1 )*( N+N ) + I
JI = ( N+I-1 )*( N+N ) + J
DWORK(IJ) =  DWORK(IJ) - DWORK(JI)
DWORK(JI) = -DWORK(IJ)
IJ = N + I + ( J-1 )*( N+N )
JI = N + J + ( I-1 )*( N+N )
DWORK(IJ) =  DWORK(IJ) - DWORK(JI)
DWORK(JI) = -DWORK(IJ)
50             CONTINUE
60          CONTINUE
DO 70 I = 1, N+N
WRITE ( NOUT, FMT = 99994 )
\$               ( DWORK(I+(J-1)*(N+N) ), J = 1,N+N )
70          CONTINUE
ENDIF
END IF
STOP
*
99999 FORMAT (' MB04ZD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (/' N is out of range.',/' N = ',I5)
99997 FORMAT (' INFO on exit from MB04ZD = ',I2)
99996 FORMAT (/' The square-reduced Hamiltonian is ')
99995 FORMAT (/' The square of the square-reduced Hamiltonian is ')
99994 FORMAT (1X,8(F10.4))
END
```
Program Data
```MB04ZD EXAMPLE PROGRAM DATA
3 N
1.0 2.0 3.0
4.0 5.0 6.0
7.0 8.0 9.0
1.0 1.0 1.0 2.0 2.0 3.0
7.0 6.0 5.0 8.0 4.0 9.0
```
Program Results
``` MB04ZD EXAMPLE PROGRAM RESULTS

The square-reduced Hamiltonian is
1.0000    3.3485    0.3436    1.0000    1.9126   -0.1072
6.7566   11.0750   -0.3014    1.9126    8.4479   -1.0790
2.3478    1.6899   -2.3868   -0.1072   -1.0790   -2.9871
7.0000    8.6275   -0.6352   -1.0000   -6.7566   -2.3478
8.6275   16.2238   -0.1403   -3.3485  -11.0750   -1.6899
-0.6352   -0.1403    1.2371   -0.3436    0.3014    2.3868

The square of the square-reduced Hamiltonian is
48.0000   80.6858   -2.5217    0.0000    1.8590  -10.5824
167.8362  298.4815   -4.0310   -1.8590    0.0000  -33.1160
0.0000    4.5325    2.5185   10.5824   33.1160    0.0000
0.0000    0.0000    0.0000   48.0000  167.8362    0.0000
0.0000    0.0000    0.0000   80.6858  298.4815    4.5325
0.0000    0.0000    0.0000   -2.5217   -4.0310    2.5185
```