## MB04DS

### Balancing a real skew-Hamiltonian matrix

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

Purpose

```  To balance a real skew-Hamiltonian matrix

[  A   G  ]
S =  [       T ] ,
[  Q   A  ]

where A is an N-by-N matrix and G, Q are N-by-N skew-symmetric
matrices. This involves, first, permuting S by a symplectic
similarity transformation to isolate eigenvalues in the first
1:ILO-1 elements on the diagonal of A; and second, applying a
diagonal similarity transformation to rows and columns
ILO:N, N+ILO:2*N to make the rows and columns as close in 1-norm
as possible. Both steps are optional.

```
Specification
```      SUBROUTINE MB04DS( JOB, N, A, LDA, QG, LDQG, ILO, SCALE, INFO )
C     .. Scalar Arguments ..
CHARACTER         JOB
INTEGER           ILO, INFO, LDA, LDQG, N
C     .. Array Arguments ..
DOUBLE PRECISION  A(LDA,*), QG(LDQG,*), SCALE(*)

```
Arguments

Mode Parameters

```  JOB     CHARACTER*1
Specifies the operations to be performed on S:
= 'N':  none, set ILO = 1, SCALE(I) = 1.0, I = 1 .. N;
= 'P':  permute only;
= 'S':  scale only;
= 'B':  both permute and scale.

```
Input/Output Parameters
```  N       (input) INTEGER
The order of the matrix A. N >= 0.

A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the matrix A.
On exit, the leading N-by-N part of this array contains
the matrix A of the balanced skew-Hamiltonian. In
particular, the strictly lower triangular part of the
first ILO-1 columns of A is zero.

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

QG      (input/output) DOUBLE PRECISION array, dimension
(LDQG,N+1)
On entry, the leading N-by-N+1 part of this array must
contain in columns 1:N the strictly lower triangular part
of the matrix Q and in columns 2:N+1 the strictly upper
triangular part of the matrix G. The parts containing the
diagonal and the first supdiagonal of this array are not
referenced.
On exit, the leading N-by-N+1 part of this array contains
the strictly lower and strictly upper triangular parts of
the matrices Q and G, respectively, of the balanced
skew-Hamiltonian. In particular, the strictly lower
triangular part of the first ILO-1 columns of QG is zero.

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

ILO     (output) INTEGER
ILO-1 is the number of deflated eigenvalues in the
balanced skew-Hamiltonian matrix.

SCALE   (output) DOUBLE PRECISION array of dimension (N)
Details of the permutations and scaling factors applied to
S.  For j = 1,...,ILO-1 let P(j) = SCALE(j). If P(j) <= N,
then rows and columns P(j) and P(j)+N are interchanged
with rows and columns j and j+N, respectively. If
P(j) > N, then row and column P(j)-N are interchanged with
row and column j+N by a generalized symplectic
permutation. For j = ILO,...,N the j-th element of SCALE
contains the factor of the scaling applied to row and
column j.

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

```
References
```  [1] Benner, P.
Symplectic balancing of Hamiltonian matrices.
SIAM J. Sci. Comput., 22 (5), pp. 1885-1904, 2001.

```
```  None
```
Example

Program Text

```*     MB04DS EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
INTEGER          NIN, NOUT
PARAMETER        ( NIN = 5, NOUT = 6 )
INTEGER          NMAX
PARAMETER        ( NMAX = 100 )
INTEGER          LDA, LDQG
PARAMETER        ( LDA = NMAX, LDQG = NMAX )
*     .. Local Scalars ..
CHARACTER*1      JOB
INTEGER          I, ILO, INFO, J, N
*     .. Local Arrays ..
DOUBLE PRECISION A(LDA, NMAX), DUMMY(1), QG(LDQG, NMAX+1),
\$                 SCALE(NMAX)
*     .. External Functions ..
DOUBLE PRECISION DLANTR, DLAPY2
EXTERNAL         DLANTR, DLAPY2
*     .. External Subroutines ..
EXTERNAL         MB04DS
*     .. Executable Statements ..
WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * )  N, JOB
IF( N.LE.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99994 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( QG(I,J), J = 1,N+1 ), I = 1,N )
CALL MB04DS( JOB, N, A, LDA, QG, LDQG, ILO, SCALE, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 )
DO 30  I = 1, N
WRITE (NOUT, FMT = 99995) ( A(I,J), J = 1,N )
30          CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 40  I = 1, N
WRITE (NOUT, FMT = 99995) ( QG(I,J), J = 1,N+1 )
40          CONTINUE
WRITE (NOUT, FMT = 99993)  ILO
IF ( ILO.GT.1 ) THEN
WRITE (NOUT, FMT = 99992) DLAPY2( DLANTR( 'Frobenius',
\$                 'Lower', 'No Unit', N-1, ILO-1, A(2,1), LDA,
\$                 DUMMY ), DLANTR( 'Frobenius', 'Lower', 'No Unit',
\$                 N-1, ILO-1, QG(2,1), LDQG, DUMMY ) )
END IF
END IF
END IF
*
99999 FORMAT (' MB04DS EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB04DS = ',I2)
99997 FORMAT (' The balanced matrix A is ')
99996 FORMAT (/' The balanced matrix QG is ')
99995 FORMAT (20(1X,F9.4))
99994 FORMAT (/' N is out of range.',/' N = ',I5)
99993 FORMAT (/' ILO = ',I4)
99992 FORMAT (/' Norm of subdiagonal blocks: ',G7.2)
END
```
Program Data
```MB04DS EXAMPLE PROGRAM DATA
6       B
0.0576         0    0.5208         0    0.7275   -0.7839
0.1901    0.0439    0.1663    0.0928    0.6756   -0.5030
0.5962         0    0.4418         0   -0.5955    0.7176
0.5869         0    0.3939    0.0353    0.6992   -0.0147
0.2222         0   -0.3663         0    0.5548   -0.4608
0         0         0         0         0    0.1338
0         0   -0.9862   -0.4544   -0.4733    0.4435         0
0         0         0   -0.6927    0.6641    0.4453         0
-0.3676         0         0         0    0.0841    0.3533         0
0         0         0         0         0    0.0877         0
0.9561         0    0.4784         0         0         0         0
-0.0164   -0.4514   -0.8289   -0.6831   -0.1536         0         0
```
Program Results
``` MB04DS EXAMPLE PROGRAM RESULTS

The balanced matrix A is
0.1338    0.4514    0.6831    0.8289    0.1536    0.0164
0.0000    0.0439    0.0928    0.1663    0.6756    0.1901
0.0000    0.0000    0.0353    0.3939    0.6992    0.5869
0.0000    0.0000    0.0000    0.4418   -0.5955    0.5962
0.0000    0.0000    0.0000   -0.3663    0.5548    0.2222
0.0000    0.0000    0.0000    0.5208    0.7275    0.0576

The balanced matrix QG is
0.0000    0.0000    0.5030    0.0147   -0.7176    0.4608    0.7839
0.0000    0.0000    0.0000    0.6641   -0.6927    0.4453    0.9862
0.0000    0.0000    0.0000    0.0000   -0.0841    0.0877    0.4733
0.0000    0.0000    0.0000    0.0000    0.0000    0.3533    0.4544
0.0000    0.0000    0.0000    0.4784    0.0000    0.0000   -0.4435
0.0000    0.0000    0.0000    0.3676   -0.9561    0.0000    0.0000

ILO =    4

Norm of subdiagonal blocks: 0.0
```