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

  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.
          No computations are performed if N = 0 or M = 0, but SCALE
          and DIF (if JOB <> 'N') are set to 1.

  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 computed generalized (real)
          Schur canonical form matrices P'*A*Q, U'*B*V, P'*D*Q and
          U'*E*V (see METHOD), or input matrices A, B, D, and E (if
          already reduced to these forms), 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,11*MN,10*MN+23,2*M*N) if REDUCE = 'R';
             LDWORK = MAX(1,11*M, 10*M+23, 2*M*N) if REDUCE = 'A';
             LDWORK = MAX(1,11*N, 10*N+23, 2*M*N) if REDUCE = 'B';
             LDWORK = MAX(1,2*M*N)                if REDUCE = 'N',
          where MN = max(M,N).
          Otherwise, the term 2*M*N above should be omitted.
          For optimum performance LDWORK should be larger.

          If LDWORK = -1  a workspace query is assumed; the
          routine only calculates the optimal size of the DWORK
          array, returns this value as the first entry of the DWORK
          array, and no error message is issued by XERBLA.

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-2017 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(11*MAX(MMAX,NMAX),
     $                                10*MAX(MMAX,NMAX)+23,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

Return to index