## 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
```  [1] 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.

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

[3] 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 [1] or [2].

```
```  None
```
Example

Program Text

```*     MB04ZD EXAMPLE PROGRAM TEXT.
*     Copyright (c) 2002-2012 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
```

Click here to get a compressed (gzip) tar file containing the source code of the routine, the example program, data, documentation, and related files.