Purpose
To reorder the diagonal blocks of a principal submatrix of an upper quasi-triangular matrix A together with their eigenvalues by constructing an orthogonal similarity transformation UT. After reordering, the leading block of the selected submatrix of A has eigenvalues in a suitably defined domain of interest, usually related to stability/instability in a continuous- or discrete-time sense.Specification
SUBROUTINE MB03QD( DICO, STDOM, JOBU, N, NLOW, NSUP, ALPHA,
$ A, LDA, U, LDU, NDIM, DWORK, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, JOBU, STDOM
INTEGER INFO, LDA, LDU, N, NDIM, NLOW, NSUP
DOUBLE PRECISION ALPHA
C .. Array Arguments ..
DOUBLE PRECISION A(LDA,*), DWORK(*), U(LDU,*)
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies the type of the spectrum separation to be
performed as follows:
= 'C': continuous-time sense;
= 'D': discrete-time sense.
STDOM CHARACTER*1
Specifies whether the domain of interest is of stability
type (left part of complex plane or inside of a circle)
or of instability type (right part of complex plane or
outside of a circle) as follows:
= 'S': stability type domain;
= 'U': instability type domain.
JOBU CHARACTER*1
Indicates how the performed orthogonal transformations UT
are accumulated, as follows:
= 'I': U is initialized to the unit matrix and the matrix
UT is returned in U;
= 'U': the given matrix U is updated and the matrix U*UT
is returned in U.
Input/Output Parameters
N (input) INTEGER
The order of the matrices A and U. N >= 1.
NLOW, (input) INTEGER
NSUP NLOW and NSUP specify the boundary indices for the rows
and columns of the principal submatrix of A whose diagonal
blocks are to be reordered. 1 <= NLOW <= NSUP <= N.
ALPHA (input) DOUBLE PRECISION
The boundary of the domain of interest for the eigenvalues
of A. If DICO = 'C', ALPHA is the boundary value for the
real parts of eigenvalues, while for DICO = 'D',
ALPHA >= 0 represents the boundary value for the moduli of
eigenvalues.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain a matrix in a real Schur form whose 1-by-1 and
2-by-2 diagonal blocks between positions NLOW and NSUP
are to be reordered.
On exit, the leading N-by-N part contains the ordered
real Schur matrix UT' * A * UT with the elements below the
first subdiagonal set to zero.
The leading NDIM-by-NDIM part of the principal submatrix
D = A(NLOW:NSUP,NLOW:NSUP) has eigenvalues in the domain
of interest and the trailing part of this submatrix has
eigenvalues outside the domain of interest.
The domain of interest for lambda(D), the eigenvalues of
D, is defined by the parameters ALPHA, DICO and STDOM as
follows:
For DICO = 'C':
Real(lambda(D)) < ALPHA if STDOM = 'S';
Real(lambda(D)) > ALPHA if STDOM = 'U'.
For DICO = 'D':
Abs(lambda(D)) < ALPHA if STDOM = 'S';
Abs(lambda(D)) > ALPHA if STDOM = 'U'.
LDA INTEGER
The leading dimension of array A. LDA >= N.
U (input/output) DOUBLE PRECISION array, dimension (LDU,N)
On entry with JOBU = 'U', the leading N-by-N part of this
array must contain a transformation matrix (e.g. from a
previous call to this routine).
On exit, if JOBU = 'U', the leading N-by-N part of this
array contains the product of the input matrix U and the
orthogonal matrix UT used to reorder the diagonal blocks
of A.
On exit, if JOBU = 'I', the leading N-by-N part of this
array contains the matrix UT of the performed orthogonal
transformations.
Array U need not be set on entry if JOBU = 'I'.
LDU INTEGER
The leading dimension of array U. LDU >= N.
NDIM (output) INTEGER
The number of eigenvalues of the selected principal
submatrix lying inside the domain of interest.
If NLOW = 1, NDIM is also the dimension of the invariant
subspace corresponding to the eigenvalues of the leading
NDIM-by-NDIM submatrix. In this case, if U is the
orthogonal transformation matrix used to compute and
reorder the real Schur form of A, its first NDIM columns
form an orthonormal basis for the above invariant
subspace.
Workspace
DWORK DOUBLE PRECISION array, dimension (N)Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: A(NLOW,NLOW-1) is nonzero, i.e. A(NLOW,NLOW) is not
the leading element of a 1-by-1 or 2-by-2 diagonal
block of A, or A(NSUP+1,NSUP) is nonzero, i.e.
A(NSUP,NSUP) is not the bottom element of a 1-by-1
or 2-by-2 diagonal block of A;
= 2: two adjacent blocks are too close to swap (the
problem is very ill-conditioned).
Method
Given an upper quasi-triangular matrix A with 1-by-1 or 2-by-2 diagonal blocks, the routine reorders its diagonal blocks along with its eigenvalues by performing an orthogonal similarity transformation UT' * A * UT. The column transformation UT is also performed on the given (initial) transformation U (resulted from a possible previous step or initialized as the identity matrix). After reordering, the eigenvalues inside the region specified by the parameters ALPHA, DICO and STDOM appear at the top of the selected diagonal block between positions NLOW and NSUP. In other words, lambda(A(NLOW:NSUP,NLOW:NSUP)) are ordered such that lambda(A(NLOW:NLOW+NDIM-1,NLOW:NLOW+NDIM-1)) are inside and lambda(A(NLOW+NDIM:NSUP,NLOW+NDIM:NSUP)) are outside the domain of interest. If NLOW = 1, the first NDIM columns of U*UT span the corresponding invariant subspace of A.References
[1] Stewart, G.W.
HQR3 and EXCHQZ: FORTRAN subroutines for calculating and
ordering the eigenvalues of a real upper Hessenberg matrix.
ACM TOMS, 2, pp. 275-280, 1976.
Numerical Aspects
3 The algorithm requires less than 4*N operations.Further Comments
NoneExample
Program Text
* MB03QD 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, LDU
PARAMETER ( LDA = NMAX, LDU = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = 3*NMAX )
* .. Local Scalars ..
CHARACTER*1 DICO, JOBU, STDOM
INTEGER I, INFO, J, N, NDIM, NLOW, NSUP
DOUBLE PRECISION ALPHA
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), DWORK(LDWORK), U(LDU,NMAX),
$ WI(NMAX), WR(NMAX)
LOGICAL BWORK(NMAX)
* .. External Functions ..
LOGICAL SELECT
* .. External Subroutines ..
EXTERNAL DGEES, MB03QD
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, NLOW, NSUP, ALPHA, DICO, STDOM, JOBU
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99992 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
* Compute Schur form, eigenvalues and Schur vectors.
CALL DGEES( 'Vectors', 'Not sorted', SELECT, N, A, LDA, NDIM,
$ WR, WI, U, LDU, DWORK, LDWORK, BWORK, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
* Block reordering.
CALL MB03QD( DICO, STDOM, JOBU, N, NLOW, NSUP, ALPHA,
$ A, LDA, U, LDU, NDIM, DWORK, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99997 ) INFO
ELSE
WRITE ( NOUT, FMT = 99996 ) NDIM
WRITE ( NOUT, FMT = 99994 )
DO 10 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N )
10 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( U(I,J), J = 1,N )
20 CONTINUE
END IF
END IF
END IF
*
STOP
*
99999 FORMAT (' MB03QD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from DGEES = ',I2)
99997 FORMAT (' INFO on exit from MB03QD = ',I2)
99996 FORMAT (' The number of eigenvalues in the domain is ',I5)
99995 FORMAT (8X,20(1X,F8.4))
99994 FORMAT (/' The ordered Schur form matrix is ')
99993 FORMAT (/' The transformation matrix is ')
99992 FORMAT (/' N is out of range.',/' N = ',I5)
END
Program Data
MB03QD EXAMPLE PROGRAM DATA 4 1 4 0.0 C S U -1.0 37.0 -12.0 -12.0 -1.0 -10.0 0.0 4.0 2.0 -4.0 7.0 -6.0 2.0 2.0 7.0 -9.0Program Results
MB03QD EXAMPLE PROGRAM RESULTS
The number of eigenvalues in the domain is 4
The ordered Schur form matrix is
-3.1300 -26.5066 27.2262 -16.2009
0.9070 -3.1300 13.6254 8.9206
0.0000 0.0000 -3.3700 0.3419
0.0000 0.0000 -1.7879 -3.3700
The transformation matrix is
0.9611 0.1784 0.2064 -0.0440
-0.1468 -0.2704 0.8116 -0.4965
-0.2224 0.7675 0.4555 0.3924
-0.0733 0.5531 -0.3018 -0.7730