MB03KE

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

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.

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

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.

  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.

  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.

  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.

  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.

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

  The implemented method is numerically backward stable.

  None

Program Text

  None

  None

  None