MB03RW

Solution of a Sylvester equation -AX + XB = C, with A and B in complex Schur form, aborting the computations when the norm of X is too large

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

Purpose

  To solve the Sylvester equation -AX + XB = C, where A and B are
  complex M-by-M and N-by-N matrices, respectively, in Schur form.

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

Specification
      SUBROUTINE MB03RW( M, N, PMAX, A, LDA, B, LDB, C, LDC, INFO )
C     .. Scalar Arguments ..
      INTEGER           INFO, LDA, LDB, LDC, M, N
      DOUBLE PRECISION  PMAX
C     .. Array Arguments ..
      COMPLEX*16        A(LDA,*), B(LDB,*), C(LDC,*)

Arguments

Input/Output Parameters

  M       (input) INTEGER
          The order of the matrix A and the number of rows of the
          matrices C and X.  M >= 0.

  N       (input) INTEGER
          The order of the matrix B and the number of columns of the
          matrices C and X.  N >= 0.

  PMAX    (input) DOUBLE PRECISION
          An upper bound for the absolute value of the elements of X
          (see METHOD).

  A       (input) COMPLEX*16 array, dimension (LDA,M)
          The leading M-by-M upper triangular part of this array
          must contain the matrix A of the Sylvester equation.
          The elements below the diagonal are not referenced.

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

  B       (input) COMPLEX*16 array, dimension (LDB,N)
          The leading N-by-N upper triangular part of this array
          must contain the matrix B of the Sylvester equation.
          The elements below the diagonal are not referenced.

  LDB     INTEGER
          The leading dimension of 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 matrix C of the Sylvester equation.
          On exit, if INFO = 0, the leading M-by-N part of this
          array contains the solution matrix X of the Sylvester
          equation, and each element of X (see METHOD) has the
          absolute value less than or equal to PMAX.
          On exit, if INFO = 1, the solution matrix X has not been
          computed completely, because an element of X had the
          absolute value greater than PMAX. Part of the matrix C has
          possibly been overwritten with the corresponding part
          of X.

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

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          = 1:  an element of X had the absolute value greater than
                the given value PMAX.
          = 2:  A and B have common or very close eigenvalues;
                perturbed values were used to solve the equation
                (but the matrices A and B are unchanged). This is a
                warning.

Method
  The routine uses an adaptation of the standard method for solving
  Sylvester equations [1], which controls the magnitude of the
  individual elements of the computed solution [2]. The equation
  -AX + XB = C can be rewritten as
                               m            l-1
    -A  X   + X  B   = C   +  sum  A  X   - sum  X  B
      kk kl    kl ll    kl   i=k+1  ki il   j=1   kj jl

  for l = 1:n, and k = m:-1:1, where A  , B  , C  , and X  , are the
                                      kk   ll   kl       kl
  elements defined by the partitioning induced by the Schur form
  of A and B. So, the elements of X are found column by column,
  starting from the bottom. If any such element has the absolute
  value greater than the given value PMAX, the calculations are
  ended.

References
  [1] Bartels, R.H. and Stewart, G.W.  T
      Solution of the matrix equation A X + XB = C.
      Comm. A.C.M., 15, pp. 820-826, 1972.

  [2] Bavely, C. and Stewart, G.W.
      An Algorithm for Computing Reducing Subspaces by Block
      Diagonalization.
      SIAM J. Numer. Anal., 16, pp. 359-367, 1979.

Numerical Aspects
                            2      2
  The algorithm requires 0(M N + MN ) operations.

Further Comments
  Let

         ( A   C )       ( I   X )
     M = (       ),  Y = (       ).
         ( 0   B )       ( 0   I )

  Then

      -1      ( A   0 )
     Y  M Y = (       ),
              ( 0   B )

  hence Y is a non-unitary transformation matrix which performs the
  reduction of M to a block-diagonal form. Bounding a norm of X is
  equivalent to setting an upper bound to the condition number of
  the transformation matrix Y.

Example

Program Text

  None
Program Data
  None
Program Results
  None

Return to Supporting Routines index