**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.

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,*)

**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.

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.

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

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.

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

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.

[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. SIAM, Philadelphia, 1995. [2] Benner, P. Symplectic balancing of Hamiltonian matrices. SIAM J. Sci. Comput., 22 (5), pp. 1885-1904, 2001.

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

**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

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

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