Purpose
To compute the periodic Schur decomposition and the eigenvalues of
a product of matrices, H = A*B, with A upper Hessenberg and B
upper triangular without evaluating any part of the product.
Specifically, the matrices Q and Z are computed, so that
Q' * A * Z = S, Z' * B * Q = T
where S is in real Schur form, and T is upper triangular.
Specification
SUBROUTINE MB03XP( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB,
$ Q, LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA, DWORK,
$ LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER COMPQ, COMPZ, JOB
INTEGER IHI, ILO, 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
JOB CHARACTER*1
Indicates whether the user wishes to compute the full
Schur form or the eigenvalues only, as follows:
= 'E': Compute the eigenvalues only;
= 'S': compute the factors S and T of the full
Schur form.
COMPQ CHARACTER*1
Indicates whether or not the user wishes to accumulate
the matrix Q as follows:
= 'N': The matrix Q is not required;
= 'I': Q is initialized to the unit matrix and the
orthogonal transformation matrix Q is returned;
= 'V': Q must contain an orthogonal matrix U on entry,
and the product U*Q is returned.
COMPZ CHARACTER*1
Indicates whether or not the user wishes to accumulate
the matrix Z as follows:
= 'N': The matrix Z is not required;
= 'I': Z is initialized to the unit matrix and the
orthogonal transformation matrix Z is returned;
= 'V': Z must contain an orthogonal matrix U on entry,
and the product U*Z is returned.
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 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
JOB = 'S'.
1 <= ILO <= max(1,N+1); min(ILO,N) <= IHI <= N.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array A must
contain the upper Hessenberg matrix A.
On exit, if JOB = 'S', the leading N-by-N part of this
array is upper quasi-triangular with any 2-by-2 diagonal
blocks corresponding to a pair of complex conjugated
eigenvalues.
If JOB = 'E', 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 B must
contain the upper triangular matrix B.
On exit, if JOB = 'S', 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 JOB = 'E', 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 COMPQ = 'V', then the leading N-by-N part of
this array must contain a matrix Q which is assumed to be
equal to the unit matrix except for the submatrix
Q(ILO:IHI,ILO:IHI).
If COMPQ = 'I', Q need not be set on entry.
On exit, if COMPQ = 'V' or COMPQ = 'I' the leading N-by-N
part of this array contains the transformation matrix
which produced the Schur form.
If COMPQ = 'N', Q is not referenced.
LDQ INTEGER
The leading dimension of the array Q. LDQ >= 1.
If COMPQ <> 'N', LDQ >= MAX(1,N).
Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
On entry, if COMPZ = 'V', then the leading N-by-N part of
this array must contain a matrix Z which is assumed to be
equal to the unit matrix except for the submatrix
Z(ILO:IHI,ILO:IHI).
If COMPZ = 'I', Z need not be set on entry.
On exit, if COMPZ = 'V' or COMPZ = 'I' the leading N-by-N
part of this array contains the transformation matrix
which produced the Schur form.
If COMPZ = 'N', Z is not referenced.
LDZ INTEGER
The leading dimension of the array Z. LDZ >= 1.
If COMPZ <> 'N', 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 (1 <= i <= N) 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 JOB = 'S', the eigenvalues are stored in the
same order as on the diagonales of the Schur forms of A
and B.
Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal
value of 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 MB03XP failed to compute the Schur
form in a total of 30*(IHI-ILO+1) iterations;
elements 1:ilo-1 and i+1:n of ALPHAR, ALPHAI and
BETA contain successfully computed eigenvalues.
Method
The implemented algorithm is a multi-shift version of the periodic QR algorithm described in [1,3] with some minor modifications proposed in [2].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
* MB03XP EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX
PARAMETER ( NMAX = 200 )
INTEGER LDA, LDB, LDQ, LDRES, LDZ, LDWORK
PARAMETER ( LDA = NMAX, LDB = NMAX, LDQ = NMAX,
$ LDRES = NMAX, LDWORK = NMAX, LDZ = NMAX )
* .. Local Scalars ..
INTEGER I, IHI, ILO, INFO, J, N
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), ALPHAI(NMAX), ALPHAR(NMAX),
$ B(LDA,NMAX), BETA(NMAX), DWORK(LDWORK),
$ Q(LDQ,NMAX), RES(LDRES,3*NMAX), Z(LDZ,NMAX)
* .. External Functions ..
DOUBLE PRECISION DLANGE
EXTERNAL DLANGE
* .. External Subroutines ..
EXTERNAL DGEMM, MB03XP
* .. Executable Statements ..
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, ILO, IHI
IF( N.LE.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99990 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
CALL DLACPY( 'All', N, N, A, LDA, RES(1,N+1), LDRES )
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,N )
CALL DLACPY( 'All', N, N, B, LDB, RES(1,2*N+1), LDRES )
CALL MB03XP( 'S', 'I', 'I', N, ILO, IHI, A, LDA, B, LDB, Q,
$ LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA, DWORK, LDWORK,
$ INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99996 )
DO 10 I = 1, N
WRITE (NOUT, FMT = 99991) ( A(I,J), J = 1,N )
10 CONTINUE
CALL DGEMM( 'No Transpose', 'No Transpose', N, N, N, ONE,
$ RES(1,N+1), LDRES, Z, LDZ, ZERO, RES, LDRES )
CALL DGEMM( 'No Transpose', 'No Transpose', N, N, N, -ONE,
$ Q, LDQ, A, LDA, ONE, RES, LDRES )
WRITE ( NOUT, FMT = 99989 ) DLANGE( 'Frobenius', N, N, RES,
$ LDRES, DWORK )
WRITE ( NOUT, FMT = 99995 )
DO 20 I = 1, N
WRITE (NOUT, FMT = 99991) ( B(I,J), J = 1,N )
20 CONTINUE
CALL DGEMM( 'No Transpose', 'No Transpose', N, N, N, ONE,
$ RES(1,2*N+1), LDRES, Q, LDQ, ZERO, RES, LDRES )
CALL DGEMM( 'No Transpose', 'No Transpose', N, N, N, -ONE,
$ Z, LDZ, B, LDB, ONE, RES, LDRES )
WRITE ( NOUT, FMT = 99988 ) DLANGE( 'Frobenius', N, N, RES,
$ LDRES, DWORK )
WRITE ( NOUT, FMT = 99994 )
DO 30 I = 1, N
WRITE (NOUT, FMT = 99991) ( Q(I,J), J = 1,N )
30 CONTINUE
CALL DGEMM( 'Transpose', 'No Transpose', N, N, N, ONE, Q,
$ LDQ, Q, LDQ, ONE, RES, LDRES )
DO 40 I = 1, N
RES(I,I) = RES(I,I) - ONE
40 CONTINUE
WRITE ( NOUT, FMT = 99987 ) DLANGE( 'Frobenius', N, N, RES,
$ LDRES, DWORK )
WRITE ( NOUT, FMT = 99993 )
DO 50 I = 1, N
WRITE (NOUT, FMT = 99991) ( Z(I,J), J = 1,N )
50 CONTINUE
CALL DGEMM( 'Transpose', 'No Transpose', N, N, N, ONE, Z,
$ LDZ, Z, LDZ, ONE, RES, LDRES )
DO 60 I = 1, N
RES(I,I) = RES(I,I) - ONE
60 CONTINUE
WRITE ( NOUT, FMT = 99986 ) DLANGE( 'Frobenius', N, N, RES,
$ LDRES, DWORK )
WRITE ( NOUT, FMT = 99992 )
DO 70 I = 1, N
WRITE ( NOUT, FMT = 99991 )
$ ALPHAR(I), ALPHAI(I), BETA(I)
70 CONTINUE
END IF
END IF
*
STOP
*
99999 FORMAT (' MB03XP EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB03XP = ',I2)
99996 FORMAT (' The reduced matrix A is ')
99995 FORMAT (/' The reduced matrix B is ')
99994 FORMAT (/' The orthogonal factor Q is ')
99993 FORMAT (/' The orthogonal factor Z is ')
99992 FORMAT (/4X,'ALPHAR',4X,'ALPHAI',4X,'BETA')
99991 FORMAT (1000(1X,F9.4))
99990 FORMAT (/' N is out of range.',/' N = ',I5)
99989 FORMAT (/' Residual: || A*Z - Q*S ||_F = ',G7.2)
99988 FORMAT (/' Residual: || B*Q - Z*T ||_F = ',G7.2)
99987 FORMAT (/' Orthogonality of Q: || Q''*Q - I ||_F = ',G7.2)
99986 FORMAT (/' Orthogonality of Z: || Z''*Z - I ||_F = ',G7.2)
END
Program Data
MB03XP EXAMPLE PROGRAM DATA
8 1 8
0.9708 -1.1156 -0.0884 -0.2684 0.2152 0.0402 0.0333 0.5141
-1.6142 2.8635 1.0420 -0.2295 -0.3560 0.4885 0.1026 -0.0164
0 1.1138 0.3509 -0.0963 0.0875 0.2158 0.2444 -0.2838
0 0 -0.5975 0.1021 -0.1026 -0.0062 -0.2646 -0.0745
0 0 0 0.6181 0.1986 0.3612 -0.1750 0.3332
0 0 0 0 -0.7387 -0.5201 0.0713 0.0501
0 0 0 0 0 -0.2677 -0.4918 -0.2838
0 0 0 0 0 0 0.3011 0.3389
0.9084 0.1739 0.5915 0.8729 0.8188 0.1911 0.4122 0.5527
0 0.1708 0.1197 0.2379 0.4302 0.4225 0.9016 0.4001
0 0 0.0381 0.6458 0.8903 0.8560 0.0056 0.1988
0 0 0 0.9669 0.7349 0.4902 0.2974 0.6252
0 0 0 0 0.6873 0.8159 0.0492 0.7334
0 0 0 0 0 0.4608 0.6932 0.3759
0 0 0 0 0 0 0.6501 0.0099
0 0 0 0 0 0 0 0.4199
Program Results
MB03XP EXAMPLE PROGRAM RESULTS
The reduced matrix A is
-0.6290 -0.1397 -0.0509 0.1603 -0.3248 0.2381 0.0694 0.0103
1.5112 -3.4273 -0.4485 -0.4357 -0.3456 0.4619 0.5998 0.5654
0.0000 0.0000 0.0547 -0.4360 0.1714 -0.2103 -0.0900 -0.4011
0.0000 0.0000 0.6623 0.2038 0.2796 -0.2629 0.3837 0.2382
0.0000 0.0000 0.0000 0.0000 -0.6315 0.2071 -0.0174 -0.3538
0.0000 0.0000 0.0000 0.0000 0.0000 -0.5850 -0.1813 0.2435
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.7884 0.1535
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.2832
Residual: || A*Z - Q*S ||_F = .60E-14
The reduced matrix B is
-0.9231 0.0000 -0.9834 0.1805 0.4428 0.3655 -0.4300 0.8498
0.0000 -0.1837 -0.1873 0.0681 0.8412 -0.0556 0.0538 0.6113
0.0000 0.0000 -1.8997 0.0000 0.5651 -0.2785 0.2882 1.0458
0.0000 0.0000 0.0000 -0.2602 0.3527 -0.0020 -0.3396 0.2739
0.0000 0.0000 0.0000 0.0000 0.8521 -0.0164 0.2115 0.5446
0.0000 0.0000 0.0000 0.0000 0.0000 0.0283 -0.5128 0.0153
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.4153 0.4587
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.5894
Residual: || B*Q - Z*T ||_F = .55E-14
The orthogonal factor Q is
-0.5333 0.3661 -0.1179 0.0264 0.0026 0.7527 0.0018 0.0189
0.0583 -0.8833 -0.0666 -0.0007 0.0017 0.4603 0.0050 0.0092
-0.8414 -0.2927 0.0347 0.0452 -0.0005 -0.4498 -0.0269 0.0001
0.0077 0.0046 -0.5687 -0.4810 0.0227 -0.0708 -0.6500 0.1312
0.0598 0.0059 -0.6128 0.7656 0.1348 -0.0863 0.0038 0.0954
-0.0242 -0.0016 -0.4295 -0.4163 0.3871 -0.0709 0.6964 -0.0417
0.0027 0.0001 0.3109 0.0620 0.8615 0.0378 -0.2267 0.3231
0.0012 0.0000 0.0188 -0.0514 -0.2987 -0.0172 0.2010 0.9312
Orthogonality of Q: || Q'*Q - I ||_F = .63E-14
The orthogonal factor Z is
0.9957 -0.0786 0.0397 -0.0032 0.0006 0.0227 0.0104 0.0123
0.0764 0.9956 0.0200 0.0073 -0.0009 0.0389 0.0263 0.0193
-0.0062 0.0235 0.6714 -0.0229 0.0271 -0.4461 -0.5354 -0.2486
-0.0445 -0.0437 0.6098 0.4197 -0.0656 0.6125 0.1248 0.2302
-0.0242 -0.0148 0.4049 -0.6041 0.2808 -0.1328 0.5972 0.1311
0.0096 0.0037 -0.0183 0.6539 0.5114 -0.4136 0.3620 -0.0913
-0.0019 -0.0004 -0.1055 -0.1544 0.7891 0.2944 -0.4436 0.2426
-0.0005 0.0000 -0.0039 0.0826 -0.1786 -0.3853 -0.1119 0.8946
Orthogonality of Z: || Z'*Z - I ||_F = .78E-14
ALPHAR ALPHAI BETA
0.4723 0.1464 1.2811
0.4723 -0.1464 1.2811
-0.0318 0.1527 2.4691
-0.0318 -0.1527 2.4691
-0.6315 0.0000 0.8521
-0.5850 0.0000 0.0283
-0.7884 0.0000 0.4153
0.2832 0.0000 0.5894