Purpose
To compute the eigenvalues of a real N-by-N skew-Hamiltonian/
skew-Hamiltonian pencil aS - bT with
( A D ) ( B F )
S = ( ) and T = ( ). (1)
( E A' ) ( G B' )
Optionally, if JOB = 'T', the pencil aS - bT will be transformed
to the structured Schur form: an orthogonal transformation matrix
Q is computed such that
( Aout Dout )
J Q' J' S Q = ( ), and
( 0 Aout' )
(2)
( Bout Fout ) ( 0 I )
J Q' J' T Q = ( ), where J = ( ),
( 0 Bout' ) ( -I 0 )
Aout is upper triangular, and Bout is upper quasi-triangular. The
notation M' denotes the transpose of the matrix M.
Optionally, if COMPQ = 'I' or COMPQ = 'U', the orthogonal
transformation matrix Q will be computed.
Specification
SUBROUTINE MB04FD( JOB, COMPQ, N, A, LDA, DE, LDDE, B, LDB, FG,
$ LDFG, Q, LDQ, ALPHAR, ALPHAI, BETA, IWORK,
$ DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER COMPQ, JOB
INTEGER INFO, LDA, LDB, LDDE, LDFG, LDQ, LDWORK, N
C .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
$ B( LDB, * ), BETA( * ), DE( LDDE, * ),
$ DWORK( * ), FG( LDFG, * ), Q( LDQ, * )
INTEGER IWORK( * )
Arguments
Mode Parameters
JOB CHARACTER*1
Specifies the computation to be performed, as follows:
= 'E': compute the eigenvalues only; S and T will not
necessarily be put into skew-Hamiltonian
triangular form (2);
= 'T': put S and T into skew-Hamiltonian triangular form
(2), and return the eigenvalues in ALPHAR, ALPHAI
and BETA.
COMPQ CHARACTER*1
Specifies whether to compute the orthogonal transformation
matrix Q as follows:
= 'N': Q is not computed;
= 'I': the array Q is initialized internally to the unit
matrix, and the orthogonal matrix Q is returned;
= 'U': the array Q contains an orthogonal matrix Q0 on
entry, and the product Q0*Q is returned, where Q
is the product of the orthogonal transformations
that are applied to the pencil aS - bT to reduce
S and T to the forms in (2), for COMPQ = 'I'.
Input/Output Parameters
N (input) INTEGER
The order of the pencil aS - bT. N >= 0, even.
A (input/output) DOUBLE PRECISION array, dimension
(LDA, N/2)
On entry, the leading N/2-by-N/2 part of this array must
contain the matrix A.
On exit, if JOB = 'T', the leading N/2-by-N/2 part of this
array contains the matrix Aout; otherwise, it contains
meaningless elements, except for the diagonal blocks,
which are correctly set.
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1, N/2).
DE (input/output) DOUBLE PRECISION array, dimension
(LDDE, N/2+1)
On entry, the leading N/2-by-N/2 strictly lower triangular
part of this array must contain the strictly lower
triangular part of the skew-symmetric matrix E, and the
N/2-by-N/2 strictly upper triangular part of the submatrix
in the columns 2 to N/2+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 are not referenced, but are assumed to be zero.
On exit, if JOB = 'T', the leading N/2-by-N/2 strictly
upper triangular part of the submatrix in the columns
2 to N/2+1 of this array contains the strictly upper
triangular part of the skew-symmetric matrix Dout.
If JOB = 'E', the leading N/2-by-N/2 strictly upper
triangular part of the submatrix in the columns 2 to N/2+1
of this array contains the strictly upper triangular part
of the skew-symmetric matrix D just before the application
of the QZ algorithm. The remaining entries are
meaningless.
LDDE INTEGER
The leading dimension of the array DE.
LDDE >= MAX(1, N/2).
B (input/output) DOUBLE PRECISION array, dimension
(LDB, N/2)
On entry, the leading N/2-by-N/2 part of this array must
contain the matrix B.
On exit, if JOB = 'T', the leading N/2-by-N/2 part of this
array contains the matrix Bout; otherwise, it contains
meaningless elements, except for the diagonal 1-by-1 and
2-by-2 blocks, which are correctly set.
LDB INTEGER
The leading dimension of the array B. LDB >= MAX(1, N/2).
FG (input/output) DOUBLE PRECISION array, dimension
(LDFG, N/2+1)
On entry, the leading N/2-by-N/2 strictly lower triangular
part of this array must contain the strictly lower
triangular part of the skew-symmetric matrix G, and the
N/2-by-N/2 strictly upper triangular part of the submatrix
in the columns 2 to N/2+1 of this array must contain the
strictly upper triangular part of the skew-symmetric
matrix F.
The entries on the diagonal and the first superdiagonal of
this array are not referenced, but are assumed to be zero.
On exit, if JOB = 'T', the leading N/2-by-N/2 strictly
upper triangular part of the submatrix in the columns
2 to N/2+1 of this array contains the strictly upper
triangular part of the skew-symmetric matrix Fout.
If JOB = 'E', the leading N/2-by-N/2 strictly upper
triangular part of the submatrix in the columns 2 to N/2+1
of this array contains the strictly upper triangular part
of the skew-symmetric matrix F just before the application
of the QZ algorithm. The remaining entries are
meaningless.
LDFG INTEGER
The leading dimension of the array FG.
LDFG >= MAX(1, N/2).
Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
On entry, if COMPQ = 'U', then the leading N-by-N part of
this array must contain a given matrix Q0, and on exit,
the leading N-by-N part of this array contains the product
of the input matrix Q0 and the transformation matrix Q
used to transform the matrices S and T.
On exit, if COMPQ = 'I', then the leading N-by-N part of
this array contains the orthogonal transformation matrix
Q.
If COMPQ = 'N', this array is not referenced.
LDQ INTEGER
The leading dimension of the array Q.
LDQ >= 1, if COMPQ = 'N';
LDQ >= MAX(1, N), if COMPQ = 'I' or COMPQ = 'U'.
ALPHAR (output) DOUBLE PRECISION array, dimension (N/2)
The real parts of each scalar alpha defining an eigenvalue
of the pencil aS - bT.
ALPHAI (output) DOUBLE PRECISION array, dimension (N/2)
The imaginary parts of each scalar alpha defining an
eigenvalue of the pencil aS - bT.
If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
positive, then the j-th and (j+1)-st eigenvalues are a
complex conjugate pair.
BETA (output) DOUBLE PRECISION array, dimension (N/2)
The scalars beta that define the eigenvalues of the pencil
aS - bT.
Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
beta = BETA(j) represent the j-th eigenvalue of the pencil
aS - bT, in the form lambda = alpha/beta. Since lambda may
overflow, the ratios should not, in general, be computed.
Due to the skew-Hamiltonian/skew-Hamiltonian structure of
the pencil, every eigenvalue occurs twice and thus it has
only to be saved once in ALPHAR, ALPHAI and BETA.
Workspace
IWORK INTEGER array, dimension (N/2+1)
On exit, IWORK(1) contains the number of (pairs of)
possibly inaccurate eigenvalues, q <= N/2, and the
absolute values in IWORK(2), ..., IWORK(q+1) are their
indices, as well as of the corresponding diagonal blocks.
Specifically, a positive value is an index of a real
eigenvalue, corresponding to a 1-by-1 block pair, while
the absolute value of a negative entry in IWORK is an
index to the first eigenvalue in a pair of consecutively
stored eigenvalues, corresponding to a 2-by-2 block pair.
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal LDWORK;
DWORK(2) and DWORK(3) contain the Frobenius norms of the
matrices S and T on entry. These norms are used in the
tests to decide that some eigenvalues are considered as
unreliable.
On exit, if INFO = -19, DWORK(1) returns the minimum
value of LDWORK.
LDWORK INTEGER
The dimension of the array DWORK.
LDWORK >= MAX(3,N/2), if JOB = 'E' and COMPQ = 'N';
LDWORK >= MAX(3,N**2/4+N/2), if JOB = 'T' and COMPQ = 'N';
LDWORK >= MAX(1,3*N**2/4), if COMPQ<> 'N'.
For good performance LDWORK should generally be larger.
If LDWORK = -1, then a workspace query is assumed;
the routine only calculates the optimal size of the
DWORK array, returns this value as the first entry of
the DWORK array, and no error message related to LDWORK
is issued by XERBLA.
Error Indicator
INFO INTEGER
= 0: succesful exit;
< 0: if INFO = -i, the i-th argument had an illegal value;
= 1: QZ iteration failed in the LAPACK Library routine
DHGEQZ. (QZ iteration did not converge or computation
of the shifts failed.)
= 2: warning: the pencil is numerically singular.
Method
The algorithm uses Givens rotations and Householder reflections to
annihilate elements in S and T such that S is in skew-Hamiltonian
triangular form and T is in skew-Hamiltonian Hessenberg form:
( A1 D1 ) ( B1 F1 )
S = ( ), T = ( ),
( 0 A1' ) ( 0 B1' )
where A1 is upper triangular and B1 is upper Hessenberg.
Subsequently, the QZ algorithm is applied to the pencil aA1 - bB1
to determine orthogonal matrices Q1 and Q2 such that
Q2' A1 Q1 is upper triangular and Q2' B1 Q1 is upper quasi-
triangular.
See also page 40 in [1] for more details.
References
[1] Benner, P., Byers, R., Mehrmann, V. and Xu, H.
Numerical Computation of Deflating Subspaces of Embedded
Hamiltonian Pencils.
Tech. Rep. SFB393/99-15, Technical University Chemnitz,
Germany, June 1999.
Numerical Aspects
3 The algorithm is numerically backward stable and needs O(N ) real floating point operations.Further Comments
NoneExample
Program Text
* MB04FD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX
PARAMETER ( NMAX = 50 )
INTEGER LDA, LDB, LDDE, LDFG, LDQ, LDWORK
PARAMETER ( LDA = NMAX/2, LDB = NMAX/2, LDDE = NMAX/2,
$ LDFG = NMAX/2, LDQ = NMAX, LDWORK =
$ 3*NMAX*NMAX/4 )
*
* .. Local Scalars ..
CHARACTER COMPQ, JOB
INTEGER I, INFO, J, N
*
* .. Local Arrays ..
DOUBLE PRECISION A( LDA, NMAX/2 ), ALPHAI( NMAX/2 ),
$ ALPHAR( NMAX/2 ), B( LDB, NMAX/2 ),
$ BETA( NMAX/2 ), DE( LDDE, NMAX/2+1 ),
$ DWORK( LDWORK ), FG( LDFG, NMAX/2+1 ),
$ Q( LDQ, NMAX )
INTEGER IWORK( NMAX/2+1 )
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
*
* .. External Subroutines ..
EXTERNAL MB04FD
*
* .. Intrinsic Functions ..
INTRINSIC MOD
*
* .. Executable statements ..
*
WRITE( NOUT, FMT = 99999 )
*
* Skip first line in data file.
*
READ( NIN, FMT = * )
READ( NIN, FMT = * ) JOB, COMPQ, N
READ( NIN, FMT = * ) ( ( A( I, J ), J = 1, N/2 ), I = 1, N/2 )
READ( NIN, FMT = * ) ( ( DE( I, J ), J = 1, N/2+1 ), I = 1, N/2 )
READ( NIN, FMT = * ) ( ( B( I, J ), J = 1, N/2 ), I = 1, N/2 )
READ( NIN, FMT = * ) ( ( FG( I, J ), J = 1, N/2+1 ), I = 1, N/2 )
IF( LSAME( COMPQ, 'U' ) )
$ READ( NIN, FMT = * ) ( ( Q( I, J ), J = 1, N ), I = 1, N )
IF( N.LT.0 .OR. N.GT.NMAX .OR. MOD( N, 2 ).NE.0 ) THEN
WRITE( NOUT, FMT = 99998 ) N
ELSE
*
* Test of MB04FD.
*
CALL MB04FD( JOB, COMPQ, N, A, LDA, DE, LDDE, B, LDB, FG, LDFG,
$ Q, LDQ, ALPHAR, ALPHAI, BETA, IWORK, DWORK,
$ LDWORK, INFO )
IF( INFO.NE.0 ) THEN
WRITE( NOUT, FMT = 99997 ) INFO
ELSE
IF( LSAME( JOB, 'T' ) ) THEN
WRITE( NOUT, FMT = 99996 )
DO 10 I = 1, N/2
WRITE( NOUT, FMT = 99995 ) ( A( I, J ), J = 1, N/2 )
10 CONTINUE
END IF
WRITE( NOUT, FMT = 99994 )
DO 20 I = 1, N/2
WRITE( NOUT, FMT = 99995 ) ( DE( I, J ), J = 1, N/2+1 )
20 CONTINUE
IF( LSAME( JOB, 'T' ) ) THEN
WRITE( NOUT, FMT = 99993 )
DO 30 I = 1, N/2
WRITE( NOUT, FMT = 99995 ) ( B( I, J ), J = 1, N/2 )
30 CONTINUE
END IF
WRITE( NOUT, FMT = 99992 )
DO 40 I = 1, N/2
WRITE( NOUT, FMT = 99995 ) ( FG( I, J ), J = 1, N/2+1 )
40 CONTINUE
IF( .NOT.LSAME( COMPQ, 'N' ) ) THEN
WRITE( NOUT, FMT = 99991 )
DO 50 I = 1, N
WRITE( NOUT, FMT = 99995 ) ( Q( I, J ), J = 1, N )
50 CONTINUE
END IF
WRITE( NOUT, FMT = 99990 )
WRITE( NOUT, FMT = 99995 ) ( ALPHAR(I), I = 1, N/2 )
WRITE( NOUT, FMT = 99995 ) ( ALPHAI(I), I = 1, N/2 )
WRITE( NOUT, FMT = 99995 ) ( BETA(I), I = 1, N/2 )
END IF
END IF
STOP
99999 FORMAT ( 'MB04FD EXAMPLE PROGRAM RESULTS', 1X )
99998 FORMAT ( 'N is out of range.', /, 'N = ', I5 )
99997 FORMAT ( 'INFO on exit from MB04FD = ', I2 )
99996 FORMAT (/' The transformed matrix A is' )
99995 FORMAT ( 51( 1X, F8.4 ) )
99994 FORMAT (/' The transformed matrix DE is' )
99993 FORMAT (/' The transformed matrix B is' )
99992 FORMAT (/' The transformed matrix FG is' )
99991 FORMAT (/' The transformed matrix Q is ' )
99990 FORMAT (/' The real, imaginary, and beta parts of eigenvalues are'
$ )
END
Program Data
MB04FD EXAMPLE PROGRAM DATA
T I 8
0.8147 0.6323 0.9575 0.9571
0.9057 0.0975 0.9648 0.4853
0.1269 0.2784 0.1576 0.8002
0.9133 0.5468 0.9705 0.1418
0.4217 0.6557 0.6787 0.6554 0.2769
0.9157 0.0357 0.7577 0.1711 0.0461
0.7922 0.8491 0.7431 0.7060 0.0971
0.9594 0.9339 0.3922 0.0318 0.8234
0.6948 0.4387 0.1868 0.7093
0.3170 0.3815 0.4897 0.7546
0.9502 0.7655 0.4455 0.2760
0.0344 0.7951 0.6463 0.6797
0.6550 0.9597 0.7512 0.8909 0.1492
0.1626 0.3403 0.2550 0.9592 0.2575
0.1189 0.5852 0.5059 0.5472 0.8407
0.4983 0.2238 0.6990 0.1386 0.2542
Program Results
MB04FD EXAMPLE PROGRAM RESULTS The transformed matrix A is 0.0550 -0.3064 0.1543 0.2170 0.0000 1.2189 0.3267 -1.3622 0.0000 0.0000 0.7734 -0.6215 0.0000 0.0000 0.0000 2.5172 The transformed matrix DE is 0.4217 0.6557 0.4277 -0.2877 -0.3980 -1.5448 0.0357 0.7577 -0.2582 -0.4775 0.3220 1.1004 0.7431 0.7060 -1.1102 0.3899 0.3463 0.4843 0.0318 0.8234 The transformed matrix B is 0.8482 -0.4425 -0.3643 0.8333 0.0000 -0.5919 -0.0987 -0.7923 0.0000 0.0000 1.1021 0.1926 0.0000 0.0000 0.0000 1.9788 The transformed matrix FG is 0.6550 0.9597 0.3266 -0.6059 -0.9618 -0.3151 0.3403 0.2550 -0.4771 -0.3005 0.1987 -0.0141 0.5059 0.5472 -0.4172 0.7914 0.0371 0.0000 0.1386 0.2542 The transformed matrix Q is 0.0762 0.7824 -0.2771 -0.5526 0.0000 0.0000 0.0000 0.0000 0.1786 -0.2450 0.4261 -0.5360 0.1843 -0.5857 0.2437 -0.0522 -0.2452 0.1425 0.4831 -0.0742 -0.1834 -0.0856 -0.6529 0.4620 -0.2123 -0.3870 -0.6431 -0.2548 0.3061 -0.1791 -0.2518 0.3707 -0.4085 0.0154 0.0784 -0.0738 0.3908 0.0160 -0.3950 -0.7157 0.6418 -0.2868 0.0258 -0.3305 -0.0286 0.4768 -0.3901 -0.1253 -0.0077 -0.1340 -0.2740 -0.0534 -0.7708 -0.3910 -0.1981 -0.3432 -0.5275 -0.2394 0.1131 -0.4684 -0.3018 0.4869 0.3216 0.0252 The real, imaginary, and beta parts of eigenvalues are 0.8482 -0.5919 1.1021 1.9788 0.0000 0.0000 0.0000 0.0000 0.0550 1.2189 0.7734 2.5172