Purpose
To deal with small subtasks of the product eigenvalue problem. MB03YD is an auxiliary routine called by SLICOT Library routine MB03XP.Specification
SUBROUTINE MB03YD( WANTT, WANTQ, WANTZ, N, ILO, IHI, ILOQ, IHIQ,
$ A, LDA, B, LDB, Q, LDQ, Z, LDZ, ALPHAR, ALPHAI,
$ BETA, DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
LOGICAL WANTQ, WANTT, WANTZ
INTEGER IHI, IHIQ, ILO, ILOQ, INFO, LDA, LDB, LDQ,
$ LDWORK, LDZ, N
C .. Array Arguments ..
DOUBLE PRECISION A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*),
$ BETA(*), DWORK(*), Q(LDQ,*), Z(LDZ,*)
Arguments
Mode Parameters
WANTT LOGICAL
Indicates whether the user wishes to compute the full
Schur form or the eigenvalues only, as follows:
= .TRUE. : Compute the full Schur form;
= .FALSE.: compute the eigenvalues only.
WANTQ LOGICAL
Indicates whether or not the user wishes to accumulate
the matrix Q as follows:
= .TRUE. : The matrix Q is updated;
= .FALSE.: the matrix Q is not required.
WANTZ LOGICAL
Indicates whether or not the user wishes to accumulate
the matrix Z as follows:
= .TRUE. : The matrix Z is updated;
= .FALSE.: the matrix Z is not required.
Input/Output Parameters
N (input) INTEGER
The order of the matrices A and B. N >= 0.
ILO (input) INTEGER
IHI (input) INTEGER
It is assumed that the matrices A and B are already
(quasi) upper triangular in rows and columns 1:ILO-1 and
IHI+1:N. The routine works primarily with the submatrices
in rows and columns ILO to IHI, but applies the
transformations to all the rows and columns of the
matrices A and B, if WANTT = .TRUE..
1 <= ILO <= max(1,N); min(ILO,N) <= IHI <= N.
ILOQ (input) INTEGER
IHIQ (input) INTEGER
Specify the rows of Q and Z to which transformations
must be applied if WANTQ = .TRUE. and WANTZ = .TRUE.,
respectively.
1 <= ILOQ <= ILO; IHI <= IHIQ <= N.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the upper Hessenberg matrix A.
On exit, if WANTT = .TRUE., the leading N-by-N part of
this array is upper quasi-triangular in rows and columns
ILO:IHI.
If WANTT = .FALSE., the diagonal elements and 2-by-2
diagonal blocks of A will be correct, but the remaining
parts of A are unspecified on exit.
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
On entry, the leading N-by-N part of this array must
contain the upper triangular matrix B.
On exit, if WANTT = .TRUE., the leading N-by-N part of
this array contains the transformed upper triangular
matrix. 2-by-2 blocks in B corresponding to 2-by-2 blocks
in A will be reduced to positive diagonal form. (I.e., if
A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j)
and B(j+1,j+1) will be positive.)
If WANTT = .FALSE., the elements corresponding to diagonal
elements and 2-by-2 diagonal blocks in A will be correct,
but the remaining parts of B are unspecified on exit.
LDB INTEGER
The leading dimension of the array B. LDB >= MAX(1,N).
Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
On entry, if WANTQ = .TRUE., then the leading N-by-N part
of this array must contain the current matrix Q of
transformations accumulated by MB03XP.
On exit, if WANTQ = .TRUE., then the leading N-by-N part
of this array contains the matrix Q updated in the
submatrix Q(ILOQ:IHIQ,ILO:IHI).
If WANTQ = .FALSE., Q is not referenced.
LDQ INTEGER
The leading dimension of the array Q. LDQ >= 1.
If WANTQ = .TRUE., LDQ >= MAX(1,N).
Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
On entry, if WANTZ = .TRUE., then the leading N-by-N part
of this array must contain the current matrix Z of
transformations accumulated by MB03XP.
On exit, if WANTZ = .TRUE., then the leading N-by-N part
of this array contains the matrix Z updated in the
submatrix Z(ILOQ:IHIQ,ILO:IHI).
If WANTZ = .FALSE., Z is not referenced.
LDZ INTEGER
The leading dimension of the array Z. LDZ >= 1.
If WANTZ = .TRUE., LDZ >= MAX(1,N).
ALPHAR (output) DOUBLE PRECISION array, dimension (N)
ALPHAI (output) DOUBLE PRECISION array, dimension (N)
BETA (output) DOUBLE PRECISION array, dimension (N)
The i-th (ILO <= i <= IHI) computed eigenvalue is given
by BETA(I) * ( ALPHAR(I) + sqrt(-1)*ALPHAI(I) ). If two
eigenvalues are computed as a complex conjugate pair,
they are stored in consecutive elements of ALPHAR, ALPHAI
and BETA. If WANTT = .TRUE., the eigenvalues are stored in
the same order as on the diagonals of the Schur forms of
A and B.
Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = -19, DWORK(1) returns the minimum
value of LDWORK.
LDWORK INTEGER
The length of the array DWORK. LDWORK >= MAX(1,N).
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
> 0: if INFO = i, then MB03YD failed to compute the Schur
form in a total of 30*(IHI-ILO+1) iterations;
elements i+1:n of ALPHAR, ALPHAI and BETA contain
successfully computed eigenvalues.
Method
The implemented algorithm is a double-shift version of the periodic QR algorithm described in [1,3] with some minor modifications [2]. The eigenvalues are computed via an implicit complex single shift algorithm.References
[1] Bojanczyk, A.W., Golub, G.H., and Van Dooren, P.
The periodic Schur decomposition: Algorithms and applications.
Proc. of the SPIE Conference (F.T. Luk, Ed.), 1770, pp. 31-42,
1992.
[2] Kressner, D.
An efficient and reliable implementation of the periodic QZ
algorithm. Proc. of the IFAC Workshop on Periodic Control
Systems, pp. 187-192, 2001.
[3] Van Loan, C.
Generalized Singular Values with Algorithms and Applications.
Ph. D. Thesis, University of Michigan, 1973.
Numerical Aspects
The algorithm requires O(N**3) floating point operations and is backward stable.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None