## MB04DP

### Balancing a real skew-Hamiltonian/Hamiltonian pencil, exploiting the structure

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

Purpose

```  To balance the 2*N-by-2*N skew-Hamiltonian/Hamiltonian pencil
aS - bH, with

(  A  D  )         (  C  V  )
S = (        ) and H = (        ),  A, C N-by-N,             (1)
(  E  A' )         (  W -C' )

where D and E are skew-symmetric, and V and W are symmetric
matrices. This involves, first, permuting aS - bH by a symplectic
equivalence transformation to isolate eigenvalues in the first
1:ILO-1 elements on the diagonal of A and C; and second, applying
a diagonal equivalence transformation to make the pairs of rows
and columns ILO:N and N+ILO:2*N as close in 1-norm as possible.
Both steps are optional. Balancing may reduce the 1-norms of the
matrices S and H.

```
Specification
```      SUBROUTINE MB04DP( JOB, N, THRESH, A, LDA, DE, LDDE, C, LDC, VW,
\$                   LDVW, ILO, LSCALE, RSCALE, DWORK, IWARN, INFO )
C     .. Scalar Arguments ..
CHARACTER          JOB
INTEGER            ILO, INFO, IWARN, LDA, LDC, LDDE, LDVW, N
DOUBLE PRECISION   THRESH
C     .. Array Arguments ..
DOUBLE PRECISION   A(LDA,*), C(LDC,*), DE(LDDE,*), DWORK(*),
\$                   LSCALE(*), RSCALE(*), VW(LDVW,*)

```
Arguments

Mode Parameters

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

```
Input/Output Parameters
```  N       (input) INTEGER
The order of matrices A, D, E, C, V, and W.  N >= 0.

THRESH  (input) DOUBLE PRECISION
If JOB = 'S' or JOB = 'B', and THRESH >= 0, threshold
value for magnitude of the elements to be considered in
the scaling process: elements with magnitude less than or
equal to THRESH*MXNORM are ignored for scaling, where
MXNORM is the maximum of the 1-norms of the original
submatrices S(s,s) and H(s,s), with s = [ILO:N,N+ILO:2*N].
If THRESH < 0, the subroutine finds the scaling factors
for which some conditions, detailed below, are fulfilled.
A sequence of increasing strictly positive threshold
values is used.
If THRESH = -1, the condition is that
max( norm(H(s,s),1)/norm(S(s,s),1),
norm(S(s,s),1)/norm(H(s,s),1) )                (1)
has the smallest value, for the threshold values used,
where S(s,s) and H(s,s) are the scaled submatrices.
If THRESH = -2, the norm ratio reduction (1) is tried, but
the subroutine may return IWARN = 1 and reset the scaling
factors to 1, if this seems suitable. See the description
of the argument IWARN and FURTHER COMMENTS.
If THRESH = -3, the condition is that
norm(H(s,s),1)*norm(S(s,s),1)                       (2)
has the smallest value for the scaled submatrices.
If THRESH = -4, the norm reduction in (2) is tried, but
the subroutine may return IWARN = 1 and reset the scaling
factors to 1, as for THRESH = -2 above.
If THRESH = -VALUE, with VALUE >= 10, the condition
numbers of the left and right scaling transformations will
be bounded by VALUE, i.e., the ratios between the largest
and smallest entries in [LSCALE(ILO:N); RSCALE(ILO:N)]
will be at most VALUE. VALUE should be a power of 10.
If JOB = 'N' or JOB = 'P', the value of THRESH is
irrelevant.

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 matrix S.
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).

DE      (input/output) DOUBLE PRECISION array, dimension
(LDDE, N+1)
On entry, the leading N-by-N strictly lower triangular
part of this array must contain the strictly lower
triangular part of the skew-symmetric matrix E, and the
N-by-N strictly upper triangular part of the submatrix
in the columns 2 to N+1 of this array must contain the
strictly upper triangular part of the skew-symmetric
matrix D.
The entries on the diagonal and the first superdiagonal of
this array need not be set, but are assumed to be zero.
On exit, the leading N-by-N strictly lower triangular
part of this array contains the strictly lower triangular
part of the balanced matrix E, and the N-by-N strictly
upper triangular part of the submatrix in the columns 2 to
N+1 of this array contains the strictly upper triangular
part of the balanced matrix D. In particular, the strictly
lower triangular part of the first ILO-1 columns of DE is
zero.

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

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

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

VW      (input/output) DOUBLE PRECISION array, dimension
(LDVW, N+1)
On entry, the leading N-by-N lower triangular part of
this array must contain the lower triangular part of the
symmetric matrix W, 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 triangular part of the
symmetric matrix V.
On exit, the leading N-by-N lower triangular part of this
array contains the lower triangular part of the balanced
matrix W, and the N-by-N upper triangular part of the
submatrix in the columns 2 to N+1 of this array contains
the upper triangular part of the balanced matrix V. In
particular, the lower triangular part of the first ILO-1
columns of VW is zero.

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

ILO     (output) INTEGER
ILO-1 is the number of deflated eigenvalues in the
balanced skew-Hamiltonian/Hamiltonian matrix pencil.
ILO is set to 1 if JOB = 'N' or JOB = 'S'.

LSCALE  (output) DOUBLE PRECISION array, dimension (N)
Details of the permutations of S and H and scaling applied
to A, D, C, and V from the left. For j = 1,...,ILO-1 let
P(j) = LSCALE(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 LSCALE contains the factor of the scaling applied to
row j of the matrices A, D, C, and V.

RSCALE  (output) DOUBLE PRECISION array, dimension (N)
Details of the permutations of S and H and scaling applied
to A, E, C, and W from the right. For j = 1,...,ILO-1 let
P(j) = RSCALE(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 RSCALE contains the factor of the scaling applied to
column j of the matrices A, E, C, and W.

```
Workspace
```  DWORK   DOUBLE PRECISION array, dimension (LDWORK) where
LDWORK = 0,   if  JOB = 'N' or JOB = 'P', or N = 0;
LDWORK = 6*N, if (JOB = 'S' or JOB = 'B') and THRESH >= 0;
LDWORK = 8*N, if (JOB = 'S' or JOB = 'B') and THRESH <  0.
On exit, if JOB = 'S' or JOB = 'B', DWORK(1) and DWORK(2)
contain the initial 1-norms of S(s,s) and H(s,s), and
DWORK(3) and DWORK(4) contain their final 1-norms,
respectively. Moreover, DWORK(5) contains the THRESH value
used (irrelevant if IWARN = 1 or ILO = N).

```
Warning Indicator
```  IWARN   INTEGER
= 0:  no warning;
= 1:  scaling has been requested, for THRESH = -2 or
THRESH = -4, but it most probably would not improve
the accuracy of the computed solution for a related
eigenproblem (since maximum norm increased
significantly compared to the original pencil
matrices and (very) high and/or small scaling
factors occurred). The returned scaling factors have
been reset to 1, but information about permutations,
if requested, has been preserved.

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

```
Method
```  Balancing consists of applying a (symplectic) equivalence
transformation to isolate eigenvalues and/or to make the 1-norms
of each pair of rows and columns indexed by s of S and H nearly
equal. If THRESH < 0, a search is performed to find those scaling
factors giving the smallest norm ratio or product defined above
(see the description of the parameter THRESH).

Assuming JOB = 'S', let Dl and Dr be diagonal matrices containing
the vectors LSCALE and RSCALE, respectively. The returned matrices
are obtained using the equivalence transformation

( Dl  0 ) ( A  D  ) ( Dr  0 )   ( Dl  0 ) ( C  V  ) ( Dr  0 )
(       ) (       ) (       ),  (       ) (       ) (       ).
( 0  Dr ) ( E  A' ) ( 0  Dl )   ( 0  Dr ) ( W -C' ) ( 0  Dl )

For THRESH = 0, the routine returns essentially the same results
as the LAPACK subroutine DGGBAL [1]. Setting THRESH < 0, usually
gives better results than DGGBAL for badly scaled matrix pencils.

```
References
```  [1] Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J.,
Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A.,
Ostrouchov, S., and Sorensen, D.
LAPACK Users' Guide: Second Edition.

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

```
Numerical Aspects
```  The transformations used preserve the skew-Hamiltonian/Hamiltonian
structure and do not introduce significant rounding errors.
No rounding errors appear if JOB = 'P'. If T is the global
transformation matrix applied to the right, then J'*T*J is the
global transformation matrix applied to the left, where
J = [ 0 I; -I 0 ], with blocks of order N.

```
```  If THRESH = -2, the increase of the maximum norm of the scaled
submatrices, compared to the maximum norm of the initial
submatrices, is bounded by MXGAIN = 100.
If THRESH = -2, or THRESH = -4, the maximum condition number of
the scaling transformations is bounded by MXCOND = 1/SQRT(EPS),
where EPS is the machine precision (see LAPACK Library routine
DLAMCH).

```
Example

Program Text

```*     MB04DP EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
INTEGER          NIN, NOUT
PARAMETER        ( NIN = 5, NOUT = 6 )
INTEGER          NMAX
PARAMETER        ( NMAX = 10 )
INTEGER          LDA, LDC, LDDE, LDVW
PARAMETER        ( LDA  = NMAX, LDC = NMAX, LDDE = NMAX,
\$                   LDVW = NMAX )
*     .. Local Scalars ..
CHARACTER*1      JOB
INTEGER          I, ILO, INFO, IWARN, J, N
DOUBLE PRECISION THRESH
*     .. Local Arrays ..
DOUBLE PRECISION A(LDA, NMAX), DWORK(8*NMAX), C(LDC, NMAX),
\$                 DE(LDDE, NMAX+1), LSCALE(NMAX), RSCALE(NMAX),
\$                 VW(LDVW, NMAX+1)
*     .. External Functions ..
LOGICAL          LSAME
EXTERNAL         LSAME
*     .. External Subroutines ..
EXTERNAL         MB04DP
*     .. 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, THRESH
IF( N.LE.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99985 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J),  J = 1,N ),   I = 1,N )
READ ( NIN, FMT = * ) ( ( DE(I,J), J = 1,N+1 ), I = 1,N )
READ ( NIN, FMT = * ) ( ( C(I,J),  J = 1,N ),   I = 1,N )
READ ( NIN, FMT = * ) ( ( VW(I,J), J = 1,N+1 ), I = 1,N )
CALL MB04DP( JOB, N, THRESH, A, LDA, DE, LDDE, C, LDC, VW,
\$                LDVW, ILO, LSCALE, RSCALE, DWORK, IWARN, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 )
DO 10  I = 1, N
WRITE ( NOUT, FMT = 99993 ) ( A(I,J), J = 1,N )
10          CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 20  I = 1, N
WRITE ( NOUT, FMT = 99993 ) ( DE(I,J), J = 1,N+1 )
20          CONTINUE
WRITE ( NOUT, FMT = 99995 )
DO 30  I = 1, N
WRITE ( NOUT, FMT = 99993 ) ( C(I,J), J = 1,N )
30          CONTINUE
WRITE ( NOUT, FMT = 99994 )
DO 40  I = 1, N
WRITE ( NOUT, FMT = 99993 ) ( VW(I,J), J = 1,N+1 )
40          CONTINUE
WRITE ( NOUT, FMT = 99992 )  ILO
WRITE ( NOUT, FMT = 99991 )
WRITE ( NOUT, FMT = 99993 ) ( LSCALE(I), I = 1,N )
WRITE ( NOUT, FMT = 99990 )
WRITE ( NOUT, FMT = 99993 ) ( RSCALE(I), I = 1,N )
IF ( LSAME( JOB, 'S' ) .OR. LSAME( JOB, 'B' ) ) THEN
IF ( .NOT.( THRESH.EQ.-2 .OR. THRESH.EQ.-4 ) ) THEN
WRITE ( NOUT, FMT = 99989 )
WRITE ( NOUT, FMT = 99993 ) ( DWORK(I), I = 1,2 )
WRITE ( NOUT, FMT = 99988 )
WRITE ( NOUT, FMT = 99993 ) ( DWORK(I), I = 3,4 )
WRITE ( NOUT, FMT = 99987 )
WRITE ( NOUT, FMT = 99993 ) ( DWORK(5) )
ELSE
WRITE ( NOUT, FMT = 99986 ) IWARN
END IF
END IF
END IF
END IF
*
99999 FORMAT (' MB04DP EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB04DP = ',I2)
99997 FORMAT (' The balanced matrix A is ')
99996 FORMAT (/' The balanced matrix DE is ')
99995 FORMAT (' The balanced matrix C is ')
99994 FORMAT (/' The balanced matrix VW is ')
99993 FORMAT (20(1X,G12.4))
99992 FORMAT (/' ILO = ',I4)
99991 FORMAT (/' The permutations and left scaling factors are ')
99990 FORMAT (/' The permutations and right scaling factors are ')
99989 FORMAT (/' The initial 1-norms of the (sub)matrices are ')
99988 FORMAT (/' The final 1-norms of the (sub)matrices are ')
99987 FORMAT (/' The threshold value finally used is ')
99986 FORMAT (/' IWARN on exit from MB04DP = ',I2)
99985 FORMAT (/' N is out of range.',/' N = ',I5)
END
```
Program Data
```MB04DP EXAMPLE PROGRAM DATA
2       B      -3
1         0
0         1
0         0         0
0         0         0
1         0
0        -2
-1  -1.0e-12         0
-1        -1         0

```
Program Results
``` MB04DP EXAMPLE PROGRAM RESULTS

The balanced matrix A is
1.000        0.000
0.000        1.000

The balanced matrix DE is
0.000        0.000        0.000
0.000        0.000        0.000
The balanced matrix C is
2.000        1.000
0.000        1.000

The balanced matrix VW is
0.000        1.000        0.000
0.000       -1.000      -0.1000E-11

ILO =    2

The permutations and left scaling factors are
4.000        1.000

The permutations and right scaling factors are
4.000        1.000

The initial 1-norms of the (sub)matrices are
1.000        2.000

The final 1-norms of the (sub)matrices are
1.000        2.000

The threshold value finally used is
-3.000
```