Purpose
To reorder the diagonal blocks of the formal matrix product
T22_K^S(K) * T22_K-1^S(K-1) * ... * T22_1^S(1), (1)
of length K, in the generalized periodic Schur form
[ T11_k T12_k T13_k ]
T_k = [ 0 T22_k T23_k ], k = 1, ..., K, (2)
[ 0 0 T33_k ]
where
- the submatrices T11_k are NI(k+1)-by-NI(k), if S(k) = 1, or
NI(k)-by-NI(k+1), if S(k) = -1, and contain dimension-induced
infinite eigenvalues,
- the submatrices T22_k are NC-by-NC and contain core eigenvalues,
which are generically neither zero nor infinite,
- the submatrices T33_k contain dimension-induced zero
eigenvalues,
such that the block with starting row index IFST in (1) is moved
to row index ILST. The indices refer to the T22_k submatrices.
Optionally, the transformation matrices Q_1,...,Q_K from the
reduction into generalized periodic Schur form are updated with
respect to the performed reordering.
Specification
SUBROUTINE MB03KA( COMPQ, WHICHQ, WS, K, NC, KSCHUR, IFST, ILST,
$ N, NI, S, T, LDT, IXT, Q, LDQ, IXQ, TOL, IWORK,
$ DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER COMPQ
LOGICAL WS
INTEGER IFST, ILST, INFO, K, KSCHUR, LDWORK, NC
C .. Array Arguments ..
INTEGER IWORK( * ), IXQ( * ), IXT( * ), LDQ( * ),
$ LDT( * ), N( * ), NI( * ), S( * ), WHICHQ( * )
DOUBLE PRECISION DWORK( * ), Q( * ), T( * ), TOL( * )
Arguments
Mode Parameters
COMPQ CHARACTER*1
= 'N': do not compute any of the matrices Q_k;
= 'U': each coefficient of Q must contain an orthogonal
matrix Q1_k on entry, and the products Q1_k*Q_k are
returned, where Q_k, k = 1, ..., K, performed the
reordering;
= 'W': the computation of each Q_k is specified
individually in the array WHICHQ.
WHICHQ INTEGER array, dimension (K)
If COMPQ = 'W', WHICHQ(k) specifies the computation of Q_k
as follows:
= 0: do not compute Q_k;
> 0: the kth coefficient of Q must contain an orthogonal
matrix Q1_k on entry, and the product Q1_k*Q_k is
returned.
This array is not referenced if COMPQ <> 'W'.
WS LOGICAL
= .FALSE. : do not perform the strong stability tests;
= .TRUE. : perform the strong stability tests; often,
this is not needed, and omitting them can save
some computations.
Input/Output Parameters
K (input) INTEGER
The period of the periodic matrix sequences T and Q (the
number of factors in the matrix product). K >= 2.
(For K = 1, a standard eigenvalue reordering problem is
obtained.)
NC (input) INTEGER
The number of core eigenvalues. 0 <= NC <= min(N).
KSCHUR (input) INTEGER
The index for which the matrix T22_kschur is upper quasi-
triangular. All other T22 matrices are upper triangular.
IFST (input/output) INTEGER
ILST (input/output) INTEGER
Specify the reordering of the diagonal blocks, as follows:
The block with starting row index IFST in (1) is moved to
row index ILST by a sequence of direct swaps between adjacent
blocks in the product.
On exit, if IFST pointed on entry to the second row of a
2-by-2 block in the product, it is changed to point to the
first row; ILST always points to the first row of the block
in its final position in the product (which may differ from
its input value by +1 or -1).
1 <= IFST <= NC, 1 <= ILST <= NC.
N (input) INTEGER array, dimension (K)
The leading K elements of this array must contain the
dimensions of the factors of the formal matrix product T,
such that the k-th coefficient T_k is an N(k+1)-by-N(k)
matrix, if S(k) = 1, or an N(k)-by-N(k+1) matrix,
if S(k) = -1, k = 1, ..., K, where N(K+1) = N(1).
NI (input) INTEGER array, dimension (K)
The leading K elements of this array must contain the
dimensions of the factors of the matrix sequence T11_k.
N(k) >= NI(k) + NC >= 0.
S (input) INTEGER array, dimension (K)
The leading K elements of this array must contain the
signatures (exponents) of the factors in the K-periodic
matrix sequence. Each entry in S must be either 1 or -1;
the value S(k) = -1 corresponds to using the inverse of
the factor T_k.
T (input/output) DOUBLE PRECISION array, dimension (*)
On entry, this array must contain at position IXT(k) the
matrix T_k, which is at least N(k+1)-by-N(k), if S(k) = 1,
or at least N(k)-by-N(k+1), if S(k) = -1, in periodic
Schur form.
On exit, the matrices T_k are overwritten by the reordered
periodic Schur form.
LDT INTEGER array, dimension (K)
The leading dimensions of the matrices T_k in the one-
dimensional array T.
LDT(k) >= max(1,N(k+1)), if S(k) = 1,
LDT(k) >= max(1,N(k)), if S(k) = -1.
IXT INTEGER array, dimension (K)
Start indices of the matrices T_k in the one-dimensional
array T.
Q (input/output) DOUBLE PRECISION array, dimension (*)
On entry, this array must contain at position IXQ(k) a
matrix Q_k of size at least N(k)-by-N(k), provided that
COMPQ = 'U', or COMPQ = 'W' and WHICHQ(k) > 0.
On exit, if COMPQ = 'U', or COMPQ = 'W' and WHICHQ(k) > 0,
Q_k is post-multiplied with the orthogonal matrix that
performed the reordering.
This array is not referenced if COMPQ = 'N'.
LDQ INTEGER array, dimension (K)
The leading dimensions of the matrices Q_k in the one-
dimensional array Q.
LDQ(k) >= max(1,N(k)), if COMPQ = 'U', or COMPQ = 'W' and
WHICHQ(k) > 0;
This array is not referenced if COMPQ = 'N'.
IXQ INTEGER array, dimension (K)
Start indices of the matrices Q_k in the one-dimensional
array Q.
This array is not referenced if COMPQ = 'N'.
Tolerances
TOL DOUBLE PRECISION array, dimension (3)
This array contains tolerance parameters. The weak and
strong stability tests use a threshold computed by the
formula MAX( c*EPS*NRM, SMLNUM ), where c is a constant,
NRM is the Frobenius norm of the current matrix formed by
concatenating K pairs of adjacent diagonal blocks of sizes
1 and/or 2 in the T22_k submatrices from (2), which are
swapped, and EPS and SMLNUM are the machine precision and
safe minimum divided by EPS, respectively (see LAPACK
Library routine DLAMCH). The norm NRM is computed by this
routine; the other values are stored in the array TOL.
TOL(1), TOL(2), and TOL(3) contain c, EPS, and SMLNUM,
respectively. TOL(1) should normally be at least 10.
Workspace
IWORK INTEGER array, dimension (4*K)
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal LDWORK.
LDWORK INTEGER
The dimension of the array DWORK.
LDWORK >= 10*K + MN, if all blocks between IFST and ILST
have order 1;
LDWORK >= 25*K + MN, if there is at least a block of
order 2, but no adjacent blocks of
order 2 can appear between IFST and
ILST during reordering;
LDWORK >= MAX(42*K + MN, 80*K - 48), if at least a pair of
adjacent blocks of order 2 can appear
between IFST and ILST during
reordering;
where MN = MXN, if MXN > 10, and MN = 0, otherwise, with
MXN = MAX(N(k),k=1,...,K).
If LDWORK = -1 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 is issued by XERBLA.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -21, the LDWORK argument was too small;
= 1: the reordering of T failed because some eigenvalues
are too close to separate (the problem is very ill-
conditioned); T may have been partially reordered.
The returned value of ILST is the index where this
was detected.
Method
An adaptation of the LAPACK Library routine DTGEXC is used.Numerical Aspects
The implemented method is numerically backward stable.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None