MB04RW

Solution of a generalized complex Sylvester equation with matrix pairs in generalized complex Schur form (blocked version)

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

Purpose

  To solve the generalized complex Sylvester equation

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

  using Level 3 BLAS, where R and L are unknown M-by-N matrices, and
  (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
  N-by-N and M-by-N, respectively. A, B, D and E are complex upper
  triangular (i.e., (A,D) and (B,E) are in generalized Schur form).

  The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an
  output scaling factor chosen to avoid overflow.

  This routine is intended to be called only by SLICOT Library
  routine MB04RZ. For efficiency purposes, the computations are
  aborted when the absolute value of an element of R or L is greater
  than a given value PMAX.

Specification
      SUBROUTINE MB04RW( M, N, PMAX, 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   PMAX, SCALE
C     .. Array Arguments ..
      INTEGER            IWORK( * )
      COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * ),
     $                   D( LDD, * ), E( LDE, * ), F( LDF, * )

Arguments

Input/Output Parameters

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

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

  PMAX    (input) DOUBLE PRECISION
          An upper bound for the "absolute value" of the elements of
          the solution (R, L). (See FURTHER COMMENTS.)
          PMAX >= 1.0D0.

  A       (input) COMPLEX*16 array, dimension (LDA, M)
          On entry, the leading M-by-M upper triangular part of this
          array must contain the matrix A in the generalized complex
          Schur form, as returned by LAPACK routine ZGGES.

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

  B       (input) COMPLEX*16 array, dimension (LDB, N)
          On entry, the leading N-by-N upper triangular part of this
          array must contain the matrix B in the generalized complex
          Schur form.

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

  C       (input/output) COMPLEX*16 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, if INFO = 0, 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) COMPLEX*16 array, dimension (LDD, M)
          On entry, the leading M-by-M upper triangular part of this
          array must contain the matrix D in the generalized complex
          Schur form. The diagonal elements are non-negative real.

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

  E       (input) COMPLEX*16 array, dimension (LDE, N)
          On entry, the leading N-by-N upper triangular part of this
          array must contain the matrix E in the generalized complex
          Schur form. The diagonal elements are non-negative real.

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

  F       (input/output) COMPLEX*16 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, if INFO = 0, 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 solutions
          R and L (C and F on entry) will hold the solutions to a
          slightly perturbed system but the input matrices A, B, D
          and E have not been changed. If SCALE = 0, R and L will
          hold the solutions to the homogeneous system with C = 0
          and F = 0. Normally, SCALE = 1.

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

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          = 1:  an element of R or L had the absolute value greater
                than the given value PMAX.
          = 2:  the matrix pairs (A, D) and (B, E) have common or
                very close eigenvalues. The matrix Z in section
                METHOD is (almost) singular.

Method
  The routine uses an adaptation of the method for solving
  generalized Sylvester equations [1], which controls the magnitude
  of the individual elements of the computed solution [2].

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

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

  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.

References
  [1] Kagstrom, B. and Westin, L.
      Generalized Schur Methods with Condition Estimators for
      Solving the Generalized Sylvester Equation.
      IEEE Trans. Auto. Contr., 34, pp. 745-751, 1989.
  [2] Kagstrom, B. and Westin, L.
      GSYLV - Fortran Routines for the Generalized Schur Method with
      Dif Estimators for Solving the Generalized Sylvester Equation.
      Report UMINF-132.86, Institute of Information Processing,
      Univ. of Umea, Sweden, July 1987.
  [3] Kagstrom, B.
      A Perturbation Analysis of the Generalized Sylvester Equation
      (AR - LB, DR - LE ) = (C, F).
      SIAM J. Matrix Anal. Appl., 15(4), pp. 1045-1060, 1994.
  [4] Kagstrom, B. and Poromaa, P.
      LAPACK-Style Algorithms and Software for Solving the
      Generalized Sylvester Equation and Estimating the Separation
      between Regular Matrix Pairs.
      ACM Trans. on Math. Software, 22(1), pp. 78103, 1996.

Further Comments
  For efficiency reasons, the "absolute value" of a complex number x
  is computed as |real(x)| + |imag(x)|.

Example

Program Text

  None
Program Data
  None
Program Results
  None

Return to index