SB04OW

Solving a periodic Sylvester equation with matrices in periodic Schur form

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

Purpose

  To solve a periodic Sylvester equation

           A * R - L * B = scale * C                           (1)
           D * L - R * E = scale * F,

  using Level 1 and 2 BLAS, where R and L are unknown M-by-N
  matrices, (A, D), (B, E) and (C, F) are given matrix pairs of
  size M-by-M, N-by-N and M-by-N, respectively, with real entries.
  (A, D) and (B, E) must be in periodic Schur form, i.e. A, B are
  upper quasi triangular and D, E are upper triangular. The solution
  (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output scaling
  factor chosen to avoid overflow.

  This routine is largely based on the LAPACK routine DTGSY2
  developed by Bo Kagstrom and Peter Poromaa.

Specification
      SUBROUTINE SB04OW( M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE,
     $                   F, LDF, SCALE, IWORK, INFO )
C     .. Scalar Arguments ..
      INTEGER           INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
      DOUBLE PRECISION  SCALE
C     .. Array Arguments ..
      INTEGER           IWORK(*)
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
     $                  E(LDE,*), F(LDF,*)

Arguments

Input/Output Parameters

  M       (input) INTEGER
          The order of A and D, and the row dimension of C, F, R
          and L.  M >= 0.

  N       (input) INTEGER
          The order of B and E, and the column dimension of C, F, R
          and L.  N >= 0.

  A       (input) DOUBLE PRECISION array, dimension (LDA,M)
          On entry, the leading M-by-M part of this array must
          contain the upper quasi triangular matrix A.

  LDA     INTEGER
          The leading dimension of the array A.  LDA >= MAX(1,M).

  B       (input) DOUBLE PRECISION array, dimension (LDB,N)
          On entry, the leading N-by-N part of this array must
          contain the upper quasi triangular matrix B.

  LDB     INTEGER
          The leading dimension of the array B.  LDB >= MAX(1,N).

  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
          On entry, the leading M-by-N part of this array must
          contain the right-hand-side of the first matrix equation
          in (1).
          On exit, the leading M-by-N part of this array contains
          the solution R.

  LDC     INTEGER
          The leading dimension of the array C.  LDC >= MAX(1,M).

  D       (input) DOUBLE PRECISION array, dimension (LDD,M)
          On entry, the leading M-by-M part of this array must
          contain the upper triangular matrix D.

  LDD     INTEGER
          The leading dimension of the array D.  LDD >= MAX(1,M).

  E       (input) DOUBLE PRECISION array, dimension (LDE,N)
          On entry, the leading N-by-N part of this array must
          contain the upper triangular matrix E.

  LDE     INTEGER
          The leading dimension of the array E.  LDE >= MAX(1,N).

  F       (input/output) DOUBLE PRECISION array, dimension (LDF,N)
          On entry, the leading M-by-N part of this array must
          contain the right-hand-side of the second matrix equation
          in (1).
          On exit, the leading M-by-N part of this array contains
          the solution L.

  LDF     INTEGER
          The leading dimension of the array F.  LDF >= MAX(1,M).

  SCALE   (output) DOUBLE PRECISION
          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the arrays
          C and F will hold the solutions R and L, respectively, to
          a slightly perturbed system but the input matrices A, B, D
          and E have not been changed. If SCALE = 0, C and F will
          hold solutions to the homogeneous system with C = F = 0.
          Normally, SCALE = 1.

Workspace
  IWORK   INTEGER array, dimension (M+N+2)

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          > 0:  the matrix products A*D and B*E have common or very
                close eigenvalues.

Method
  In matrix notation solving equation (1) corresponds to solving
  Z*x = scale*b, where Z is defined as

      Z = [  kron(In, A)  -kron(B', Im) ]            (2)
          [ -kron(E', Im)  kron(In, D)  ],

  Ik is the identity matrix of size k and X' is the transpose of X.
  kron(X, Y) is the Kronecker product between the matrices X and Y.
  In the process of solving (1), we solve a number of such systems
  where Dim(Im), Dim(In) = 1 or 2.

References
  [1] Kagstrom, B.
      A Direct Method for Reordering Eigenvalues in the Generalized
      Real Schur Form of a Regular Matrix Pair (A,B). M.S. Moonen
      et al (eds.), Linear Algebra for Large Scale and Real-Time
      Applications, Kluwer Academic Publ., pp. 195-218, 1993.

  [2] Sreedhar, J. and Van Dooren, P.
      A Schur approach for solving some periodic matrix equations.
      U. Helmke et al (eds.), Systems and Networks: Mathematical
      Theory and Applications, Akademie Verlag, Berlin, vol. 77,
      pp. 339-362, 1994.

Further Comments
  None
Example

Program Text

  None
Program Data
  None
Program Results
  None

Return to Supporting Routines index