Contents


NAME

     stgsyl - solve the generalized Sylvester equation

SYNOPSIS

     SUBROUTINE STGSYL(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
           E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)

     CHARACTER * 1 TRANS
     INTEGER IJOB, M, N, LDA, LDB, LDC,  LDD,  LDE,  LDF,  LWORK,
     INFO
     INTEGER IWORK(*)
     REAL SCALE, DIF
     REAL  A(LDA,*),  B(LDB,*),  C(LDC,*),  D(LDD,*),   E(LDE,*),
     F(LDF,*), WORK(*)

     SUBROUTINE STGSYL_64(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
           LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)

     CHARACTER * 1 TRANS
     INTEGER*8 IJOB, M, N, LDA, LDB, LDC, LDD, LDE,  LDF,  LWORK,
     INFO
     INTEGER*8 IWORK(*)
     REAL SCALE, DIF
     REAL  A(LDA,*),  B(LDB,*),  C(LDC,*),  D(LDD,*),   E(LDE,*),
     F(LDF,*), WORK(*)

  F95 INTERFACE
     SUBROUTINE TGSYL(TRANS, IJOB, [M], [N], A, [LDA], B, [LDB], C, [LDC],
            D, [LDD], E, [LDE], F, [LDF], SCALE, DIF, [WORK], [LWORK], [IWORK],
            [INFO])

     CHARACTER(LEN=1) :: TRANS
     INTEGER :: IJOB, M, N, LDA, LDB, LDC, LDD, LDE, LDF,  LWORK,
     INFO
     INTEGER, DIMENSION(:) :: IWORK
     REAL :: SCALE, DIF
     REAL, DIMENSION(:) :: WORK
     REAL, DIMENSION(:,:) :: A, B, C, D, E, F

     SUBROUTINE TGSYL_64(TRANS, IJOB, [M], [N], A, [LDA], B, [LDB], C,
            [LDC], D, [LDD], E, [LDE], F, [LDF], SCALE, DIF, [WORK], [LWORK],
            [IWORK], [INFO])

     CHARACTER(LEN=1) :: TRANS
     INTEGER(8) :: IJOB, M, N, LDA,  LDB,  LDC,  LDD,  LDE,  LDF,
     LWORK, INFO
     INTEGER(8), DIMENSION(:) :: IWORK
     REAL :: SCALE, DIF
     REAL, DIMENSION(:) :: WORK
     REAL, DIMENSION(:,:) :: A, B, C, D, E, F

  C INTERFACE
     #include <sunperf.h>

     void stgsyl(char trans, int ijob, int m, int  n,  float  *a,
               int  lda,  float  *b,  int ldb, float *c, int ldc,
               float *d, int ldd, float *e, int  lde,  float  *f,
               int ldf, float *scale, float *dif, int *info);

     void stgsyl_64(char trans, long ijob, long m, long n,  float
               *a,  long  lda, float *b, long ldb, float *c, long
               ldc, float *d, long ldd, float *e, long lde, float
               *f,  long  ldf,  float  *scale,  float  *dif, long
               *info);

PURPOSE

     stgsyl solves the generalized Sylvester equation:

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

     where R and L are unknown m-by-n matrices, (A,  D),  (B,  E)
     and (C, F) are given matrix pairs of size m-by-m, n-by-n and
     m-by-n, respectively, with real entries. (A, D) and  (B,  E)
     must  be in generalized (real) Schur canonical form, i.e. A,
     B are upper quasi triangular and D, E are upper triangular.

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

     In matrix notation (1) is equivalent to solve  Zx = scale b,
     where Z is defined as

                Z = [ kron(In, A)  -kron(B', Im) ]         (2)
                    [ kron(In, D)  -kron(E', Im) ].

     Here Ik is the identity matrix of size k and X' is the tran-
     spose  of X. kron(X, Y) is the Kronecker product between the
     matrices X and Y.

     If TRANS = 'T', STGSYL solves the transposed system  Z'*y  =
     scale*b, which is equivalent to solve for R and L in

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

     This case (TRANS = 'T') is used to compute an one-norm-based
     estimate  of  Dif[(A,D),  (B,E)], the separation between the
     matrix pairs (A,D) and (B,E), using SLACON.

     If IJOB >= 1, STGSYL computes a Frobenius  norm-based  esti-
     mate of Dif[(A,D),(B,E)]. That is, the reciprocal of a lower
     bound on the reciprocal of the smallest singular value of Z.
     See [1-2] for more information.

     This is a level 3 BLAS algorithm.

ARGUMENTS

     TRANS (input)
               = 'N', solve the  generalized  Sylvester  equation
               (1).  = 'T', solve the 'transposed' system (3).

     IJOB (input)
               Specifies what kind of functionality  to  be  per-
               formed.  =0: solve (1) only.
               =1: The functionality of 0 and 3.
               =2: The functionality of 0 and 4.
               =3: Only an estimate of Dif[(A,D), (B,E)] is  com-
               puted.   (look  ahead strategy IJOB  = 1 is used).
               =4: Only an estimate of Dif[(A,D), (B,E)] is  com-
               puted.   (  SGECON  on sub-systems is used ).  Not
               referenced if TRANS = 'T'.

     M (input) The order of the matrices A and  D,  and  the  row
               dimension of the matrices C, F, R and L.

     N (input) The order of the matrices B and E, and the  column
               dimension of the matrices C, F, R and L.

     A (input) The upper quasi triangular matrix A.

     LDA (input)
               The leading dimension  of  the  array  A.  LDA  >=
               max(1, M).

     B (input) The upper quasi triangular matrix B.

     LDB (input)
               The leading dimension  of  the  array  B.  LDB  >=
               max(1, N).
     C (input/output)
               On entry, C contains the  right-hand-side  of  the
               first  matrix equation in (1) or (3).  On exit, if
               IJOB = 0, 1 or 2, C has been  overwritten  by  the
               solution  R.  If  IJOB = 3 or 4 and TRANS = 'N', C
               holds R, the solution achieved during the computa-
               tion of the Dif-estimate.

     LDC (input)
               The leading dimension  of  the  array  C.  LDC  >=
               max(1, M).

     D (input) The upper triangular matrix D.

     LDD (input)
               The leading dimension  of  the  array  D.  LDD  >=
               max(1, M).

     E (input) The upper triangular matrix E.

     LDE (input)
               The leading dimension  of  the  array  E.  LDE  >=
               max(1, N).

     F (input/output)
               On entry, F contains the  right-hand-side  of  the
               second matrix equation in (1) or (3).  On exit, if
               IJOB = 0, 1 or 2, F has been  overwritten  by  the
               solution  L.  If  IJOB = 3 or 4 and TRANS = 'N', F
               holds L, the solution achieved during the computa-
               tion of the Dif-estimate.

     LDF (input)
               The leading dimension  of  the  array  F.  LDF  >=
               max(1, M).

     DIF (output)
               On exit SCALE is the reciprocal of a  lower  bound
               of  the reciprocal of the Dif-function, i.e. SCALE
               is  an  upper  bound  of   Dif[(A,D),   (B,E)]   =
               sigma_min(Z),  where  Z as in (2).  If IJOB = 0 or
               TRANS = 'T', SCALE is not touched.
     SCALE (output)
               On exit SCALE is the reciprocal of a  lower  bound
               of  the reciprocal of the Dif-function, i.e. SCALE
               is  an  upper  bound  of   Dif[(A,D),   (B,E)]   =
               sigma_min(Z),  where  Z as in (2).  If IJOB = 0 or
               TRANS = 'T', SCALE is not touched.

     WORK (workspace)
               If IJOB = 0, WORK is not  referenced.   Otherwise,
               on  exit, if INFO = 0, WORK(1) returns the optimal
               LWORK.

     LWORK (input)
               The dimension of the array WORK. LWORK > = 1.   If
               IJOB = 1 or 2 and TRANS = 'N', LWORK >= 2*M*N.

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

     IWORK (workspace)
               dimension(M+N+6)

     INFO (output)
               =0: successful exit
               <0: If INFO = -i, the i-th argument had an illegal
               value.
               >0: (A, D) and (B, E) have common or close  eigen-
               values.

FURTHER DETAILS

     Based on contributions by
        Bo Kagstrom and Peter Poromaa,  Department  of  Computing
     Science,
        Umea University, S-901 87 Umea, Sweden.

     [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms  and
     Software
         for  Solving  the  Generalized  Sylvester  Equation  and
     Estimating the
         Separation between Regular Matrix Pairs, Report UMINF  -
     93.23,
         Department of Computing Science, Umea University,  S-901
     87 Umea,
         Sweden, December  1993,  Revised  April  1994,  Also  as
     LAPACK Working
         Note 75.  To appear in ACM Trans. on Math. Software, Vol
     22,
         No 1, 1996.

     [2] B. Kagstrom, A Perturbation Analysis of the  Generalized
     Sylvester
         Equation (AR - LB, DR - LE ) = (C, F),  SIAM  J.  Matrix
     Anal.
         Appl., 15(4):1045-1060, 1994

     [3] B. Kagstrom and L.  Westin,  Generalized  Schur  Methods
     with
         Condition Estimators for Solving the Generalized Sylves-
     ter
         Equation, IEEE Transactions on Automatic  Control,  Vol.
     34, No. 7,
         July 1989, pp 745-751.