Purpose
To balance a complex Hamiltonian matrix,
[ A G ]
H = [ H ] ,
[ Q -A ]
where A is an N-by-N matrix and G, Q are N-by-N Hermitian
matrices. This involves, first, permuting H 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. Assuming ILO = 1, let D be a
diagonal matrix of order N with the scaling factors on the
diagonal. The scaled Hamiltonian is defined by
[ D**-1*A*D D**-1*G*D**-1 ]
Hs = [ H ] .
[ D*Q*D -D*A *D**-1 ]
Specification
SUBROUTINE MB04DZ( 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 SCALE(*)
COMPLEX*16 A(LDA,*), QG(LDQG,*)
Arguments
Mode Parameters
JOB CHARACTER*1
Specifies the operations to be performed on H:
= '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) COMPLEX*16 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 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) COMPLEX*16 array, dimension
(LDQG,N+1)
On entry, the leading N-by-N+1 part of this array must
contain the lower triangular part of the matrix Q and
the upper triangular part of the matrix G.
On exit, the leading N-by-N+1 part of this array contains
the lower and upper triangular parts of the matrices Q and
G, respectively, of the balanced Hamiltonian. In
particular, the 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 Hamiltonian matrix.
SCALE (output) DOUBLE PRECISION array of dimension (N)
Details of the permutations and scaling factors applied to
H. 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.
Further Comments
NoneExample
Program Text
* MB04DZ 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 ..
COMPLEX*16 A(LDA, NMAX), QG(LDQG, NMAX+1)
DOUBLE PRECISION DUMMY(1), SCALE(NMAX)
* .. External Functions ..
DOUBLE PRECISION DLAPY2, ZLANTR
EXTERNAL DLAPY2, ZLANTR
* .. External Subroutines ..
EXTERNAL MB04DZ
* .. 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 MB04DZ( 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 10 I = 1, N
WRITE (NOUT, FMT = 99995) ( A(I,J), J = 1,N )
10 CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, N
WRITE (NOUT, FMT = 99995) ( QG(I,J), J = 1,N+1 )
20 CONTINUE
WRITE (NOUT, FMT = 99993) ILO
IF ( ILO.GT.1 ) THEN
WRITE (NOUT, FMT = 99992) DLAPY2( ZLANTR( 'Frobenius',
$ 'Lower', 'No Unit', N-1, ILO-1, A(2,1), LDA,
$ DUMMY ), ZLANTR( 'Frobenius', 'Lower', 'No Unit',
$ N, ILO-1, QG(1,1), LDQG, DUMMY ) )
END IF
END IF
END IF
*
99999 FORMAT (' MB04DZ EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB04DZ = ',I2)
99997 FORMAT (' The balanced matrix A is ')
99996 FORMAT (/' The balanced matrix QG is ')
99995 FORMAT (20(1X,F9.4,SP,F9.4,S,'i '))
99994 FORMAT (/' N is out of range.',/' N = ',I5)
99993 FORMAT (/' ILO = ',I4)
99992 FORMAT (/' Norm of subdiagonal blocks: ',G7.2)
END
Program Data
MB04DZ EXAMPLE PROGRAM DATA
6 B
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(.0994,0) (0,0) (0,0) (0,0) (0,0) (.9696,0)
(.3248,0) (0,0) (0,0) (0,0) (.4372,0) (.8308,0)
(0,0) (0,0) (0,0) (.0717,0) (0,0) (0,0)
(0,0) (0,0) (0,0) (0,0) (0,0) (.1976,0)
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(0,0) (0,0) (0,0) (0,0) (.0651,0) (0,0) (0,0)
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(0,0) (0,0) (.0444,0) (0,0) (0,0) (.1957,0) (0,0)
(.8144,0) (0,0) (0,0) (0,0) (.3652,0) (0,0) (.9121,0)
(.9023,0) (0,0) (0,0) (0,0) (0,0) (0,0) 1.0945,0)
Program Results
MB04DZ EXAMPLE PROGRAM RESULTS
The balanced matrix A is
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.9696 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i -0.8144 +0.0000i -0.9023 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.1093 +0.0000i 0.2077 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0717 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.1976 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i
The balanced matrix QG is
0.0000 +0.0000i 0.0000 +0.0000i 0.0994 +0.0000i 0.0000 +0.0000i 0.0651 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0812 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.1776 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.1957 +0.0000i 0.0000 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.3652 +0.0000i 0.0000 +0.0000i 0.9121 +0.0000i
0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 0.0000 +0.0000i 1.0945 +0.0000i
ILO = 3
Norm of subdiagonal blocks: 0.0