## MB03KE

### Solving periodic Sylvester-like equations with matrices of order at most 2

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

Purpose

```  To solve small periodic Sylvester-like equations (PSLE)

op(A(i))*X( i ) + isgn*X(i+1)*op(B(i)) = -scale*C(i), S(i) =  1,
op(A(i))*X(i+1) + isgn*X( i )*op(B(i)) = -scale*C(i), S(i) = -1.

i = 1, ..., K, where op(A) means A or A**T, for the K-periodic
matrix sequence X(i) = X(i+K), where A, B and C are K-periodic
matrix sequences and A and B are in periodic real Schur form. The
matrices A(i) are M-by-M and B(i) are N-by-N, with 1 <= M, N <= 2.

```
Specification
```      SUBROUTINE MB03KE( TRANA, TRANB, ISGN, K, M, N, PREC, SMIN, S, A,
\$                   B, C, SCALE, DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
LOGICAL            TRANA, TRANB
INTEGER            INFO, ISGN, K, LDWORK, M, N
DOUBLE PRECISION   PREC, SCALE, SMIN
C     .. Array Arguments ..
INTEGER            S( * )
DOUBLE PRECISION   A( * ), B( * ), C( * ), DWORK( * )

```
Arguments

Mode Parameters

```  TRANA   LOGICAL
Specifies the form of op(A) to be used, as follows:
= .FALSE.:  op(A) = A,
= .TRUE. :  op(A) = A**T.

TRANB   LOGICAL
Specifies the form of op(B) to be used, as follows:
= .FALSE.:  op(B) = B,
= .TRUE. :  op(B) = B**T.

ISGN    INTEGER
Specifies which sign variant of the equations to solve.
ISGN = 1 or ISGN = -1.

```
Input/Output Parameters
```  K       (input) INTEGER
The period of the periodic matrix sequences A, B, C and X.
K >= 2. (For K = 1, a standard Sylvester equation is
obtained.)

M       (input) INTEGER
The order of the matrices A(i) and the number of rows of
the matrices C(i) and X(i), i = 1, ..., K.  1 <= M <= 2.

N       (input) INTEGER
The order of the matrices B(i) and the number of columns
of the matrices C(i) and X(i), i = 1, ..., K.
1 <= N <= 2.

PREC    (input) DOUBLE PRECISION
The relative machine precision. See the LAPACK Library
routine DLAMCH.

SMIN    (input) DOUBLE PRECISION
The machine safe minimum divided by PREC.

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 sequences for A and B. Each entry in S must be
either 1 or -1. Notice that it is assumed that the same
exponents are tied to both A and B on reduction to the
periodic Schur form.

A       (input) DOUBLE PRECISION array, dimension (M*M*K)
On entry, this array must contain the M-by-M matrices
A(i), for i = 1, ..., K, stored with the leading dimension
M. Matrix A(i) is stored starting at position M*M*(i-1)+1.

B       (input) DOUBLE PRECISION array, dimension (N*N*K)
On entry, this array must contain the N-by-N matrices
B(i), for i = 1, ..., K, stored with the leading dimension
N. Matrix B(i) is stored starting at position N*N*(i-1)+1.

C       (input/output) DOUBLE PRECISION array, dimension (M*N*K)
On entry, this array must contain the M-by-N matrices
C(i), for i = 1, ..., K, stored with the leading dimension
M. Matrix C(i) is stored starting at position M*N*(i-1)+1.
On exit, the matrices C(i) are overwritten by the solution
sequence X(i).

SCALE   (output) DOUBLE PRECISION
The scale factor, scale, set less than or equal to 1 to
avoid overflow in X.

```
Workspace
```  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal LDWORK.
On exit, if INFO = -21, DWORK(1) returns the minimum value
of LDWORK.

LDWORK  INTEGER
The dimension of the array DWORK.
LDWORK >= (4*K-3) * (M*N)**2 + K * M*N.

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, then LDWORK is too small; appropriate
value for LDWORK is returned in DWORK(1); the other
arguments are not tested, for efficiency;
= 1:  the solution would overflow with scale = 1, so
SCALE was set less than 1. This is a warning, not
an error.

```
Method
```  A version of the algorithm described in [1] is used. The routine
uses a sparse Kronecker product representation Z of the PSLE and
solves for X(i) from an associated linear system Z*x = c using
structured (overlapping) variants of QR factorization and backward
substitution.

```
References
```  [1] Granat, R., Kagstrom, B. and Kressner, D.
Computing periodic deflating subspaces associated with a
specified set of eigenvalues.
BIT Numerical Mathematics, vol. 47, 763-791, 2007.

```
Numerical Aspects
```  The implemented method is numerically backward stable.

```
```  None
```
Example

Program Text

```  None
```
Program Data
```  None
```
Program Results
```  None
```