**Purpose**

To deal with small subtasks of the product eigenvalue problem. MB03YD is an auxiliary routine called by SLICOT Library routine MB03XP.

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,*)

**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.

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.

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).

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.

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.

[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.

The algorithm requires O(N**3) floating point operations and is backward stable.

None

**Program Text**

None

None

None