MB03RY

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

Purpose

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

  This routine is intended to be called only by SLICOT Library
  routine MB03RD. For efficiency purposes, the computations are
  aborted when the infinity norm of an elementary submatrix of X is
  greater than a given value PMAX.

      SUBROUTINE MB03RY( 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 ..
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*)

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 infinity norm of an elementary
          submatrix of X (see METHOD).

  A       (input) DOUBLE PRECISION array, dimension (LDA,M)
          The leading M-by-M part of this array must contain the
          matrix A of the Sylvester equation, in real Schur form.
          The elements below the real Schur form are not referenced.

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

  B       (input) DOUBLE PRECISION array, dimension (LDB,N)
          The leading N-by-N part of this array must contain the
          matrix B of the Sylvester equation, in real Schur form.
          The elements below the real Schur form are not referenced.

  LDB     INTEGER
          The leading dimension of 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 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 elementary submatrix of X (see METHOD)
          has the infinity norm less than or equal to PMAX.
          On exit, if INFO = 1, the solution matrix X has not been
          computed completely, because an elementary submatrix of X
          had the infinity norm 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).

  INFO    INTEGER
          = 0:  successful exit;
          = 1:  an elementary submatrix of X had the infinity norm
                greater than the given value PMAX.

  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
                              p            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:q, and k = p:-1:1, where A  , B  , C  , and X  , are
                                      kk   ll   kl       kl
  block submatrices defined by the partitioning induced by the Schur
  form of A and B, and p and q are the numbers of the diagonal
  blocks of A and B, respectively. So, the elementary submatrices of
  X are found block column by block column, starting from the
  bottom. If any such elementary submatrix has the infinity norm
  greater than the given value PMAX, the calculations are ended.

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

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

  Let

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

  Then

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

  hence Y is an non-orthogonal 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.

Program Text

  None

  None

  None