## MB03XP

### Computing periodic Schur decomposition and eigenvalues of a matrix product A B, with A upper Hessenberg and B upper triangular

[Specification] [Arguments] [Method] [References] [Comments] [Example]

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 .

```
References
```   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.

 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.

 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.

```
```  None
```
Example

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
```