SB04OD

Solution of generalized Sylvester equations with separation estimation

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

Purpose

  To solve for R and L one of the generalized Sylvester equations

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

  or

     A' * R + D' * L = scale * C    )
                                    )                            (2)
     R * B' + L * E' = scale * (-F) )

  where A and D are M-by-M matrices, B and E are N-by-N matrices and
  C, F, R and L are M-by-N matrices.

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

  The routine also optionally computes a Dif estimate, which
  measures the separation of the spectrum of the matrix pair (A,D)
  from the spectrum of the matrix pair (B,E), Dif[(A,D),(B,E)].

Specification
      SUBROUTINE SB04OD( REDUCE, TRANS, JOBD, M, N, A, LDA, B, LDB, C,
     $                   LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, P,
     $                   LDP, Q, LDQ, U, LDU, V, LDV, IWORK, DWORK,
     $                   LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER         JOBD, REDUCE, TRANS
      INTEGER           INFO, LDA, LDB, LDC, LDD, LDE, LDF, LDP, LDQ,
     $                  LDU, LDV, LDWORK, M, N
      DOUBLE PRECISION  DIF, SCALE
C     .. Array Arguments ..
      INTEGER           IWORK(*)
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
     $                  DWORK(*), E(LDE,*), F(LDF,*), P(LDP,*),
     $                  Q(LDQ,*), U(LDU,*), V(LDV,*)

Arguments

Mode Parameters

  REDUCE  CHARACTER*1
          Indicates whether the matrix pairs (A,D) and/or (B,E) are
          to be reduced to generalized Schur form as follows:
          = 'R':  The matrix pairs (A,D) and (B,E) are to be reduced
                  to generalized (real) Schur canonical form;
          = 'A':  The matrix pair (A,D) only is to be reduced
                  to generalized (real) Schur canonical form,
                  and the matrix pair (B,E) already is in this form;
          = 'B':  The matrix pair (B,E) only is to be reduced
                  to generalized (real) Schur canonical form,
                  and the matrix pair (A,D) already is in this form;
          = 'N':  The matrix pairs (A,D) and (B,E) are already in
                  generalized (real) Schur canonical form, as
                  produced by LAPACK routine DGEES.

  TRANS   CHARACTER*1
          Indicates which of the equations, (1) or (2), is to be
          solved as follows:
          = 'N':  The generalized Sylvester equation (1) is to be
                  solved;
          = 'T':  The "transposed" generalized Sylvester equation
                  (2) is to be solved.

  JOBD    CHARACTER*1
          Indicates whether the Dif estimator is to be computed as
          follows:
          = '1':  Only the one-norm-based Dif estimate is computed
                  and stored in DIF;
          = '2':  Only the Frobenius norm-based Dif estimate is
                  computed and stored in DIF;
          = 'D':  The equation (1) is solved and the one-norm-based
                  Dif estimate is computed and stored in DIF;
          = 'F':  The equation (1) is solved and the Frobenius norm-
                  based Dif estimate is computed and stored in DIF;
          = 'N':  The Dif estimator is not required and hence DIF is
                  not referenced. (Solve either (1) or (2) only.)
          JOBD is not referenced if TRANS = 'T'.

Input/Output Parameters
  M       (input) INTEGER
          The order of the matrices A and D and the number of rows
          of the matrices C, F, R and L.  M >= 0.

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

  A       (input/output) DOUBLE PRECISION array, dimension (LDA,M)
          On entry, the leading M-by-M part of this array must
          contain the coefficient matrix A of the equation; A must
          be in upper quasi-triangular form if REDUCE = 'B' or 'N'.
          On exit, the leading M-by-M part of this array contains
          the upper quasi-triangular form of A.

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

  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
          On entry, the leading N-by-N part of this array must
          contain the coefficient matrix B of the equation; B must
          be in upper quasi-triangular form if REDUCE = 'A' or 'N'.
          On exit, the leading N-by-N part of this array contains
          the upper quasi-triangular form of B.

  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 right-hand side matrix C of the first equation
          in (1) or (2).
          On exit, if JOBD = 'N', 'D' or 'F', the leading M-by-N
          part of this array contains the solution matrix R of the
          problem; if JOBD = '1' or '2' and TRANS = 'N', the leading
          M-by-N part of this array contains the solution matrix R
          achieved during the computation of the Dif estimate.

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

  D       (input/output) DOUBLE PRECISION array, dimension (LDD,M)
          On entry, the leading M-by-M part of this array must
          contain the coefficient matrix D of the equation; D must
          be in upper triangular form if REDUCE = 'B' or 'N'.
          On exit, the leading M-by-M part of this array contains
          the upper triangular form of D.

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

  E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
          On entry, the leading N-by-N part of this array must
          contain the coefficient matrix E of the equation; E must
          be in upper triangular form if REDUCE = 'A' or 'N'.
          On exit, the leading N-by-N part of this array contains
          the upper triangular form of E.

  LDE     INTEGER
          The leading dimension of 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 matrix F of the second
          equation in (1) or (2).
          On exit, if JOBD = 'N', 'D' or 'F', the leading M-by-N
          part of this array contains the solution matrix L of the
          problem; if JOBD = '1' or '2' and TRANS = 'N', the leading
          M-by-N part of this array contains the solution matrix L
          achieved during the computation of the Dif estimate.

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

  SCALE   (output) DOUBLE PRECISION
          The scaling factor in (1) or (2). If 0 < SCALE < 1, C and
          F hold the solutions R and L, respectively, to a slightly
          perturbed system (but the input or computed generalized
          (real) Schur canonical form matrices A, B, D, and E
          have not been changed). If SCALE = 0, C and F hold the
          solutions R and L, respectively, to the homogeneous system
          with C = F = 0. Normally, SCALE = 1.

  DIF     (output) DOUBLE PRECISION
          If TRANS = 'N' and JOBD <> 'N', then DIF contains the
          value of the Dif estimator, which is an upper bound of
                                                 -1
          Dif[(A,D),(B,E)] = sigma_min(Z) = 1/||Z  ||, in either the
          one-norm, or Frobenius norm, respectively (see METHOD).
          Otherwise, DIF is not referenced.

  P       (output) DOUBLE PRECISION array, dimension (LDP,*)
          If REDUCE = 'R' or 'A', then the leading M-by-M part of
          this array contains the (left) transformation matrix used
          to reduce (A,D) to generalized Schur form.
          Otherwise, P is not referenced and can be supplied as a
          dummy array (i.e. set parameter LDP = 1 and declare this
          array to be P(1,1) in the calling program).

  LDP     INTEGER
          The leading dimension of array P.
          LDP >= MAX(1,M) if REDUCE = 'R' or 'A',
          LDP >= 1        if REDUCE = 'B' or 'N'.

  Q       (output) DOUBLE PRECISION array, dimension (LDQ,*)
          If REDUCE = 'R' or 'A', then the leading M-by-M part of
          this array contains the (right) transformation matrix used
          to reduce (A,D) to generalized Schur form.
          Otherwise, Q is not referenced and can be supplied as a
          dummy array (i.e. set parameter LDQ = 1 and declare this
          array to be Q(1,1) in the calling program).

  LDQ     INTEGER
          The leading dimension of array Q.
          LDQ >= MAX(1,M) if REDUCE = 'R' or 'A',
          LDQ >= 1        if REDUCE = 'B' or 'N'.

  U       (output) DOUBLE PRECISION array, dimension (LDU,*)
          If REDUCE = 'R' or 'B', then the leading N-by-N part of
          this array contains the (left) transformation matrix used
          to reduce (B,E) to generalized Schur form.
          Otherwise, U is not referenced and can be supplied as a
          dummy array (i.e. set parameter LDU = 1 and declare this
          array to be U(1,1) in the calling program).

  LDU     INTEGER
          The leading dimension of array U.
          LDU >= MAX(1,N) if REDUCE = 'R' or 'B',
          LDU >= 1        if REDUCE = 'A' or 'N'.

  V       (output) DOUBLE PRECISION array, dimension (LDV,*)
          If REDUCE = 'R' or 'B', then the leading N-by-N part of
          this array contains the (right) transformation matrix used
          to reduce (B,E) to generalized Schur form.
          Otherwise, V is not referenced and can be supplied as a
          dummy array (i.e. set parameter LDV = 1 and declare this
          array to be V(1,1) in the calling program).

  LDV     INTEGER
          The leading dimension of array V.
          LDV >= MAX(1,N) if REDUCE = 'R' or 'B',
          LDV >= 1        if REDUCE = 'A' or 'N'.

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

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, DWORK(1) returns the optimal value
          of LDWORK.

  LDWORK  INTEGER
          The length of the array DWORK.
          If TRANS = 'N' and JOBD = 'D' or 'F', then
             LDWORK = MAX(1,7*M,7*N,2*M*N) if REDUCE = 'R';
             LDWORK = MAX(1,7*M,2*M*N)     if REDUCE = 'A';
             LDWORK = MAX(1,7*N,2*M*N)     if REDUCE = 'B';
             LDWORK = MAX(1,2*M*N)         if REDUCE = 'N'.
          Otherwise, the term 2*M*N above should be omitted.
          For optimum performance LDWORK should be larger.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  if REDUCE <> 'N' and either (A,D) and/or (B,E)
                cannot be reduced to generalized Schur form;
          = 2:  if REDUCE = 'N' and either A or B is not in
                upper quasi-triangular form;
          = 3:  if a singular matrix was encountered during the
                computation of the solution matrices R and L, that
                is (A,D) and (B,E) have common or close eigenvalues.

Method
  For the case TRANS = 'N', and REDUCE = 'R' or 'N', the algorithm
  used by the routine consists of four steps (see [1] and [2]) as
  follows:

     (a) if REDUCE = 'R', then the matrix pairs (A,D) and (B,E) are
         transformed to generalized Schur form, i.e. orthogonal
         matrices P, Q, U and V are computed such that P' * A * Q
         and U' * B * V are in upper quasi-triangular form and
         P' * D * Q and U' * E * V are in upper triangular form;
     (b) if REDUCE = 'R', then the matrices C and F are transformed
         to give P' * C * V and P' * F * V respectively;
     (c) if REDUCE = 'R', then the transformed system

         P' * A * Q * R1 - L1 * U' * B * V = scale * P' * C * V
         P' * D * Q * R1 - L1 * U' * E * V = scale * P' * F * V

         is solved to give R1 and L1; otherwise, equation (1) is
         solved to give R and L directly. The Dif estimator
         is also computed if JOBD <> 'N'.
     (d) if REDUCE = 'R', then the solution is transformed back
         to give R = Q * R1 * V' and L = P * L1 * U'.

  By using Kronecker products, equation (1) can also be written as
  the system of linear equations Z * x = scale*y (see [1]), where

         | I*A    I*D  |
     Z = |             |.
         |-B'*I  -E'*I |

                                           -1
  If JOBD <> 'N', then a lower bound on ||Z  ||, in either the one-
  norm or Frobenius norm, is computed, which in most cases is
  a reliable estimate of the true value. Notice that since Z is a
  matrix of order 2 * M * N, the exact value of Dif (i.e., in the
  Frobenius norm case, the smallest singular value of Z) may be very
  expensive to compute.

  The case TRANS = 'N', and REDUCE = 'A' or 'B', is similar, but
  only one of the matrix pairs should be reduced and the
  calculations simplify.

  For the case TRANS = 'T', and REDUCE = 'R' or 'N', the algorithm
  is similar, but the steps (b), (c), and (d) are as follows:

     (b) if REDUCE = 'R', then the matrices C and F are transformed
         to give Q' * C * V and P' * F * U respectively;
     (c) if REDUCE = 'R', then the transformed system

         Q' * A' * P * R1 + Q' * D' * P * L1 =  scale * Q' * C * V
         R1 * V' * B' * U + L1 * V' * E' * U = -scale * P' * F * U

         is solved to give R1 and L1; otherwise, equation (2) is
         solved to give R and L directly.
     (d) if REDUCE = 'R', then the solution is transformed back
         to give R = P * R1 * V' and L = P * L1 * V'.

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] Golub, G.H., Nash, S. and Van Loan, C.F.
      A Hessenberg-Schur Method for the Problem AX + XB = C.
      IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979.
  [4] Kagstrom, B. and Van Dooren, P.
      Additive Decomposition of a Transfer Function with respect to
      a Specified Region.
      In: "Signal Processing, Scattering and Operator Theory, and
      Numerical Methods" (Eds. M.A. Kaashoek et al.).
      Proceedings of MTNS-89, Vol. 3, pp. 469-477, Birkhauser Boston
      Inc., 1990.
  [5] Kagstrom, B. and Van Dooren, P.
      A Generalized State-space Approach for the Additive
      Decomposition of a Transfer Matrix.
      Report UMINF-91.12, Institute of Information Processing, Univ.
      of Umea, Sweden, April 1991.

Numerical Aspects
  The algorithm is backward stable. A reliable estimate for the
  condition number of Z in the Frobenius norm, is (see [1])

     K(Z) = SQRT(  ||A||**2 + ||B||**2 + ||C||**2 + ||D||**2 )/DIF.

  If mu is an upper bound on the relative error of the elements of
  the matrices A, B, C, D, E and F, then the relative error in the
  actual solution is approximately mu * K(Z).

  The relative error in the computed solution (due to rounding
  errors) is approximately EPS * K(Z), where EPS is the machine
  precision (see LAPACK Library routine DLAMCH).

Further Comments
  For applications of the generalized Sylvester equation in control
  theory, see [4] and [5].

Example

Program Text

*     SB04OD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2010 NICONET e.V.
*
*     .. Parameters ..
      DOUBLE PRECISION ONE
      PARAMETER        ( ONE = 1.0D0 )
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          MMAX, NMAX
      PARAMETER        ( MMAX = 10, NMAX = 10 )
      INTEGER          LDA, LDB, LDC, LDD, LDE, LDF, LDP, LDQ, LDU, LDV
      PARAMETER        ( LDA = MMAX, LDB = NMAX, LDC = MMAX, LDD = MMAX,
     $                   LDE = NMAX, LDF = MMAX, LDP = MMAX, LDQ = MMAX,
     $                   LDU = NMAX, LDV = NMAX )
      INTEGER          LDWORK, LIWORK
      PARAMETER        ( LDWORK = MAX(7*MAX(MMAX,NMAX),2*MMAX*NMAX),
     $                   LIWORK = MMAX+NMAX+6 )
*     .. Local Scalars ..
      DOUBLE PRECISION DIF, SCALE
      INTEGER          I, INFO, J, M, N
      CHARACTER*1      JOBD, REDUCE, TRANS
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,MMAX), B(LDB,NMAX), C(LDC,NMAX),
     $                 D(LDD,MMAX), DWORK(LDWORK), E(LDE,NMAX),
     $                 F(LDF,NMAX), P(LDP,MMAX), Q(LDQ,MMAX),
     $                 U(LDU,NMAX), V(LDV,NMAX)
      INTEGER          IWORK(LIWORK)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         SB04OD
*     .. Intrinsic Functions ..
      INTRINSIC        MAX
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) M, N, REDUCE, TRANS, JOBD
      IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
         WRITE ( NOUT, FMT = 99989 ) M
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,M ), I = 1,M )
         IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
            WRITE ( NOUT, FMT = 99988 ) N
         ELSE
            READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,N )
            READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,M )
            READ ( NIN, FMT = * ) ( ( D(I,J), J = 1,M ), I = 1,M )
            READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,N )
            READ ( NIN, FMT = * ) ( ( F(I,J), J = 1,N ), I = 1,M )
*           Find the solution matrices L and R.
            CALL SB04OD( REDUCE, TRANS, JOBD, M, N, A, LDA, B, LDB, C,
     $                   LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, P,
     $                   LDP, Q, LDQ, U, LDU, V, LDV, IWORK, DWORK,
     $                   LDWORK, INFO )
*
            IF ( INFO.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99998 ) INFO
            ELSE
               WRITE ( NOUT, FMT = 99997 )
               DO 20 I = 1, M
                  WRITE ( NOUT, FMT = 99991 ) ( F(I,J), J = 1,N )
   20          CONTINUE
               WRITE ( NOUT, FMT = 99996 )
               DO 40 I = 1, M
                  WRITE ( NOUT, FMT = 99991 ) ( C(I,J), J = 1,N )
   40          CONTINUE
               IF ( LSAME( REDUCE, 'R' ).OR.LSAME( REDUCE, 'A' ) ) THEN
                  WRITE ( NOUT, FMT = 99995 )
                  DO 60 I = 1, M
                     WRITE ( NOUT, FMT = 99991 ) ( P(I,J), J = 1,M )
   60             CONTINUE
                  WRITE ( NOUT, FMT = 99994 )
                  DO 80 I = 1, M
                     WRITE ( NOUT, FMT = 99991 ) ( Q(I,J), J = 1,M )
   80             CONTINUE
               END IF
               IF ( LSAME( REDUCE, 'R' ).OR.LSAME( REDUCE, 'B' ) ) THEN
                  WRITE ( NOUT, FMT = 99993 )
                  DO 100 I = 1, N
                     WRITE ( NOUT, FMT = 99991 ) ( U(I,J), J = 1,N )
  100             CONTINUE
                  WRITE ( NOUT, FMT = 99992 )
                  DO 120 I = 1, N
                     WRITE ( NOUT, FMT = 99991 ) ( V(I,J), J = 1,N )
  120             CONTINUE
               END IF
               IF ( SCALE.NE.ONE ) WRITE ( NOUT, FMT = 99987 ) SCALE
               IF ( .NOT.LSAME( JOBD, 'N' ) )
     $            WRITE ( NOUT, FMT = 99990 ) DIF
            END IF
         END IF
      END IF
*
      STOP
*
99999 FORMAT (' SB04OD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB04OD = ',I2)
99997 FORMAT (' The solution matrix L is ')
99996 FORMAT (/' The solution matrix R is ')
99995 FORMAT (/' The left transformation matrix P is ')
99994 FORMAT (/' The right transformation matrix Q is ')
99993 FORMAT (/' The left transformation matrix U is ')
99992 FORMAT (/' The right transformation matrix V is ')
99991 FORMAT (20(1X,F8.4))
99990 FORMAT (/' DIF = ',F8.4)
99989 FORMAT (/' M is out of range.',/' M = ',I5)
99988 FORMAT (/' N is out of range.',/' N = ',I5)
99987 FORMAT (/' SCALE = ',F8.4)
      END
Program Data
 SB04OD EXAMPLE PROGRAM DATA
   3     2     R     N     D
    1.6   -3.1    1.9
   -3.8    4.2    2.4
    0.5    2.2   -4.5
    1.1    0.1
   -1.3   -3.1
   -2.0   28.9
   -5.7  -11.8
   12.9  -31.7
    2.5    0.1    1.7
   -2.5    0.0    0.9
    0.1    5.1   -7.3
    6.0    2.4
   -3.6    2.5
    0.5   23.8
  -11.0  -10.4
   39.5  -74.8
Program Results
 SB04OD EXAMPLE PROGRAM RESULTS

 The solution matrix L is 
  -0.7538  -1.6210
   2.1778   1.7005
  -3.5029   2.7961

 The solution matrix R is 
   1.3064   2.7989
   0.3698  -5.3376
  -0.8767   6.7500

 The left transformation matrix P is 
  -0.3093  -0.9502   0.0383
   0.9366  -0.2974   0.1851
  -0.1645   0.0932   0.9820

 The right transformation matrix Q is 
  -0.6097  -0.7920  -0.0314
   0.6310  -0.5090   0.5854
   0.4796  -0.3371  -0.8102

 The left transformation matrix U is 
  -0.8121   0.5835
   0.5835   0.8121

 The right transformation matrix V is 
  -0.9861   0.1660
   0.1660   0.9861

 DIF =   0.1147

Click here to get a compressed (gzip) tar file containing the source code of the routine, the example program, data, documentation, and related files.

Return to index