SB04QD

Solution of discrete-time Sylvester equations (Hessenberg-Schur method)

Purpose

  To solve for X the discrete-time Sylvester equation

     X + AXB = C,

  where A, B, C and X are general N-by-N, M-by-M, N-by-M and
  N-by-M matrices respectively. A Hessenberg-Schur method, which
  reduces A to upper Hessenberg form, H = U'AU, and B' to real
  Schur form, S = Z'B'Z (with U, Z orthogonal matrices), is used.

      SUBROUTINE SB04QD( N, M, A, LDA, B, LDB, C, LDC, Z, LDZ, IWORK,
     $                   DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
      INTEGER           INFO, LDA, LDB, LDC, LDWORK, LDZ, M, N
C     .. Array Arguments ..
      INTEGER           IWORK(*)
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), Z(LDZ,*)

Input/Output Parameters

  N       (input) INTEGER
          The order of the matrix A.  N >= 0.

  M       (input) INTEGER
          The order of the matrix B.  M >= 0.

  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the leading N-by-N part of this array must
          contain the coefficient matrix A of the equation.
          On exit, the leading N-by-N upper Hessenberg part of this
          array contains the matrix H, and the remainder of the
          leading N-by-N part, together with the elements 2,3,...,N
          of array DWORK, contain the orthogonal transformation
          matrix U (stored in factored form).

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

  B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
          On entry, the leading M-by-M part of this array must
          contain the coefficient matrix B of the equation.
          On exit, the leading M-by-M part of this array contains
          the quasi-triangular Schur factor S of the matrix B'.

  LDB     INTEGER
          The leading dimension of array B.  LDB >= MAX(1,M).

  C       (input/output) DOUBLE PRECISION array, dimension (LDC,M)
          On entry, the leading N-by-M part of this array must
          contain the coefficient matrix C of the equation.
          On exit, the leading N-by-M part of this array contains
          the solution matrix X of the problem.

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

  Z       (output) DOUBLE PRECISION array, dimension (LDZ,M)
          The leading M-by-M part of this array contains the
          orthogonal matrix Z used to transform B' to real upper
          Schur form.

  LDZ     INTEGER
          The leading dimension of array Z.  LDZ >= MAX(1,M).

  IWORK   INTEGER array, dimension (4*N)

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, DWORK(1) returns the optimal value
          of LDWORK, and DWORK(2), DWORK(3),..., DWORK(N) contain
          the scalar factors of the elementary reflectors used to
          reduce A to upper Hessenberg form, as returned by LAPACK
          Library routine DGEHRD.

  LDWORK  INTEGER
          The length of the array DWORK.
          LDWORK = MAX(1, 2*N*N + 9*N, 5*M, N + M).
          For optimum performance LDWORK should be larger.

          If LDWORK = -1, then 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 related to LDWORK is issued by
          XERBLA.

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          > 0:  if INFO = i, 1 <= i <= M, the QR algorithm failed to
                compute all the eigenvalues of B (see LAPACK Library
                routine DGEES);
          > M:  if a singular matrix was encountered whilst solving
                for the (INFO-M)-th column of matrix X.

  The matrix A is transformed to upper Hessenberg form H = U'AU by
  the orthogonal transformation matrix U; matrix B' is transformed
  to real upper Schur form S = Z'B'Z using the orthogonal
  transformation matrix Z. The matrix C is also multiplied by the
  transformations, F = U'CZ, and the solution matrix Y of the
  transformed system

     Y + HYS' = F

  is computed by back substitution. Finally, the matrix Y is then
  multiplied by the orthogonal transformation matrices, X = UYZ', in
  order to obtain the solution matrix X to the original problem.

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

  [2] Sima, V.
      Algorithms for Linear-quadratic Optimization.
      Marcel Dekker, Inc., New York, 1996.

                                      3       3      2         2
  The algorithm requires about (5/3) N  + 10 M  + 5 N M + 2.5 M N
  operations and is backward stable.

  None

Program Text

*     SB04QD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX, MMAX
      PARAMETER        ( NMAX = 20, MMAX = 20 )
      INTEGER          LDA, LDB, LDC, LDZ
      PARAMETER        ( LDA = NMAX, LDB = MMAX, LDC = NMAX,
     $                   LDZ = MMAX )
      INTEGER          LIWORK
      PARAMETER        ( LIWORK = 4*NMAX )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = MAX( 1, 2*NMAX*NMAX+9*NMAX, 5*MMAX,
     $                   NMAX+MMAX ) )
*     .. Local Scalars ..
      INTEGER          I, INFO, J, M, N
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,MMAX),
     $                 DWORK(LDWORK), Z(LDZ,MMAX)
      INTEGER          IWORK(LIWORK)
*     .. External Subroutines ..
      EXTERNAL         SB04QD
*     .. 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 = * ) N, M
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99994 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
         IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
            WRITE ( NOUT, FMT = 99993 ) M
         ELSE
            READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,M )
            READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,M ), I = 1,N )
*           Find the solution matrix X.
            CALL SB04QD( N, M, A, LDA, B, LDB, C, LDC, Z, LDZ, IWORK,
     $                   DWORK, LDWORK, INFO )
*
            IF ( INFO.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99998 ) INFO
            ELSE
               WRITE ( NOUT, FMT = 99997 )
               DO 20 I = 1, N
                  WRITE ( NOUT, FMT = 99996 ) ( C(I,J), J = 1,M )
   20          CONTINUE
               WRITE ( NOUT, FMT = 99995 )
               DO 40 I = 1, M
                  WRITE ( NOUT, FMT = 99996 ) ( Z(I,J), J = 1,M )
   40          CONTINUE
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' SB04QD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB04QD = ',I2)
99997 FORMAT (' The solution matrix X is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' The orthogonal matrix Z is ')
99994 FORMAT (/' N is out of range.',/' N = ',I5)
99993 FORMAT (/' M is out of range.',/' M = ',I5)
      END

 SB04QD EXAMPLE PROGRAM DATA
   3     3
   1.0   2.0   3.0 
   6.0   7.0   8.0 
   9.0   2.0   3.0 
   7.0   2.0   3.0 
   2.0   1.0   2.0 
   3.0   4.0   1.0 
   271.0   135.0   147.0
   923.0   494.0   482.0
   578.0   383.0   287.0

 SB04QD EXAMPLE PROGRAM RESULTS

 The solution matrix X is 
   2.0000   3.0000   6.0000
   4.0000   7.0000   1.0000
   5.0000   3.0000   2.0000

 The orthogonal matrix Z is 
   0.8337   0.5204  -0.1845
   0.3881  -0.7900  -0.4746
   0.3928  -0.3241   0.8606

Return to index