dtgsyl.f

Go to the documentation of this file.

*> \brief \b DTGSYL
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> Download DTGSYL + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsyl.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsyl.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsyl.f">
*> [TXT]</a>
*
*  Definition:
*  ===========
*
*       SUBROUTINE DTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
*                          LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
*                          IWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          TRANS
*       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
*      $                   LWORK, M, N
*       DOUBLE PRECISION   DIF, SCALE
*       ..
*       .. Array Arguments ..
*       INTEGER            IWORK( * )
*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
*      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
*      $                   WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DTGSYL 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**T, Im) ]         (2)
*>                [ kron(In, D)  -kron(E**T, Im) ].
*>
*> Here Ik is the identity matrix of size k and X**T is the transpose of
*> X. kron(X, Y) is the Kronecker product between the matrices X and Y.
*>
*> If TRANS = 'T', DTGSYL solves the transposed system Z**T*y = scale*b,
*> which is equivalent to solve for R and L in
*>
*>             A**T * R + D**T * L = scale * C           (3)
*>             R * B**T + L * E**T = 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 DLACON.
*>
*> If IJOB >= 1, DTGSYL computes a Frobenius norm-based estimate
*> 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] TRANS
*> \verbatim
*>          TRANS is CHARACTER*1
*>          = 'N': solve the generalized Sylvester equation (1).
*>          = 'T': solve the 'transposed' system (3).
*> \endverbatim
*>
*> \param[in] IJOB
*> \verbatim
*>          IJOB is INTEGER
*>          Specifies what kind of functionality to be performed.
*>          = 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 computed.
*>               (look ahead strategy IJOB  = 1 is used).
*>          = 4: Only an estimate of Dif[(A,D), (B,E)] is computed.
*>               ( DGECON on sub-systems is used ).
*>          Not referenced if TRANS = 'T'.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The order of the matrices A and D, and the row dimension of
*>          the matrices C, F, R and L.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrices B and E, and the column dimension
*>          of the matrices C, F, R and L.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, M)
*>          The upper quasi triangular matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A. LDA >= max(1, M).
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB, N)
*>          The upper quasi triangular matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B. LDB >= max(1, N).
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*>          C is DOUBLE PRECISION array, dimension (LDC, N)
*>          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 computation of the
*>          Dif-estimate.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*>          LDC is INTEGER
*>          The leading dimension of the array C. LDC >= max(1, M).
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (LDD, M)
*>          The upper triangular matrix D.
*> \endverbatim
*>
*> \param[in] LDD
*> \verbatim
*>          LDD is INTEGER
*>          The leading dimension of the array D. LDD >= max(1, M).
*> \endverbatim
*>
*> \param[in] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (LDE, N)
*>          The upper triangular matrix E.
*> \endverbatim
*>
*> \param[in] LDE
*> \verbatim
*>          LDE is INTEGER
*>          The leading dimension of the array E. LDE >= max(1, N).
*> \endverbatim
*>
*> \param[in,out] F
*> \verbatim
*>          F is DOUBLE PRECISION array, dimension (LDF, N)
*>          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 computation of the
*>          Dif-estimate.
*> \endverbatim
*>
*> \param[in] LDF
*> \verbatim
*>          LDF is INTEGER
*>          The leading dimension of the array F. LDF >= max(1, M).
*> \endverbatim
*>
*> \param[out] DIF
*> \verbatim
*>          DIF is DOUBLE PRECISION
*>          On exit DIF is the reciprocal of a lower bound of the
*>          reciprocal of the Dif-function, i.e. DIF is an upper bound of
*>          Dif[(A,D), (B,E)] = sigma_min(Z), where Z as in (2).
*>          IF IJOB = 0 or TRANS = 'T', DIF is not touched.
*> \endverbatim
*>
*> \param[out] SCALE
*> \verbatim
*>          SCALE is DOUBLE PRECISION
*>          On exit SCALE is the scaling factor in (1) or (3).
*>          If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
*>          to a slightly perturbed system but the input 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.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The dimension of the array WORK. LWORK > = 1.
*>          If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,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.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (M+N+6)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>            =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 eigenvalues.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup tgsyl
*
*> \par Contributors:
*  ==================
*>
*>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
*>     Umea University, S-901 87 Umea, Sweden.
*
*> \par References:
*  ================
*>
*> \verbatim
*>
*>  [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 Sylvester
*>      Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
*>      July 1989, pp 745-751.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE dtgsyl( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC,
     $                   D,
     $                   LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
     $                   IWORK, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          TRANS
      INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
     $                   LWORK, M, N
      DOUBLE PRECISION   DIF, SCALE
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
     $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
     $                   work( * )
*     ..
*
*  =====================================================================
*  Replaced various illegal calls to DCOPY by calls to DLASET.
*  Sven Hammarling, 1/5/02.
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, NOTRAN
      INTEGER            I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
     $                   LINFO, LWMIN, MB, NB, P, PPQQ, PQ, Q
      DOUBLE PRECISION   DSCALE, DSUM, SCALE2, SCALOC
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      EXTERNAL           LSAME, ILAENV
*     ..
*     .. External Subroutines ..
      EXTERNAL           dgemm, dlacpy, dlaset, dscal, dtgsy2,
     $                   xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, max, sqrt
*     ..
*     .. Executable Statements ..
*
*     Decode and test input parameters
*
      info = 0
      notran = lsame( trans, 'N' )
      lquery = ( lwork.EQ.-1 )
*
      IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) ) THEN
         info = -1
      ELSE IF( notran ) THEN
         IF( ( ijob.LT.0 ) .OR. ( ijob.GT.4 ) ) THEN
            info = -2
         END IF
      END IF
      IF( info.EQ.0 ) THEN
         IF( m.LE.0 ) THEN
            info = -3
         ELSE IF( n.LE.0 ) THEN
            info = -4
         ELSE IF( lda.LT.max( 1, m ) ) THEN
            info = -6
         ELSE IF( ldb.LT.max( 1, n ) ) THEN
            info = -8
         ELSE IF( ldc.LT.max( 1, m ) ) THEN
            info = -10
         ELSE IF( ldd.LT.max( 1, m ) ) THEN
            info = -12
         ELSE IF( lde.LT.max( 1, n ) ) THEN
            info = -14
         ELSE IF( ldf.LT.max( 1, m ) ) THEN
            info = -16
         END IF
      END IF
*
      IF( info.EQ.0 ) THEN
         IF( notran ) THEN
            IF( ijob.EQ.1 .OR. ijob.EQ.2 ) THEN
               lwmin = max( 1, 2*m*n )
            ELSE
               lwmin = 1
            END IF
         ELSE
            lwmin = 1
         END IF
         work( 1 ) = lwmin
*
         IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
            info = -20
         END IF
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DTGSYL', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( m.EQ.0 .OR. n.EQ.0 ) THEN
         scale = 1
         IF( notran ) THEN
            IF( ijob.NE.0 ) THEN
               dif = 0
            END IF
         END IF
         RETURN
      END IF
*
*     Determine optimal block sizes MB and NB
*
      mb = ilaenv( 2, 'DTGSYL', trans, m, n, -1, -1 )
      nb = ilaenv( 5, 'DTGSYL', trans, m, n, -1, -1 )
*
      isolve = 1
      ifunc = 0
      IF( notran ) THEN
         IF( ijob.GE.3 ) THEN
            ifunc = ijob - 2
            CALL dlaset( 'F', m, n, zero, zero, c, ldc )
            CALL dlaset( 'F', m, n, zero, zero, f, ldf )
         ELSE IF( ijob.GE.1 ) THEN
            isolve = 2
         END IF
      END IF
*
      IF( ( mb.LE.1 .AND. nb.LE.1 ) .OR. ( mb.GE.m .AND. nb.GE.n ) )
     $     THEN
*
         DO 30 iround = 1, isolve
*
*           Use unblocked Level 2 solver
*
            dscale = zero
            dsum = one
            pq = 0
            CALL dtgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc,
     $                   d,
     $                   ldd, e, lde, f, ldf, scale, dsum, dscale,
     $                   iwork, pq, info )
            IF( dscale.NE.zero ) THEN
               IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
                  dif = sqrt( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
               ELSE
                  dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
               END IF
            END IF
*
            IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
               IF( notran ) THEN
                  ifunc = ijob
               END IF
               scale2 = scale
               CALL dlacpy( 'F', m, n, c, ldc, work, m )
               CALL dlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
               CALL dlaset( 'F', m, n, zero, zero, c, ldc )
               CALL dlaset( 'F', m, n, zero, zero, f, ldf )
            ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
               CALL dlacpy( 'F', m, n, work, m, c, ldc )
               CALL dlacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
               scale = scale2
            END IF
   30    CONTINUE
*
         RETURN
      END IF
*
*     Determine block structure of A
*
      p = 0
      i = 1
   40 CONTINUE
      IF( i.GT.m )
     $   GO TO 50
      p = p + 1
      iwork( p ) = i
      i = i + mb
      IF( i.GE.m )
     $   GO TO 50
      IF( a( i, i-1 ).NE.zero )
     $   i = i + 1
      GO TO 40
   50 CONTINUE
*
      iwork( p+1 ) = m + 1
      IF( iwork( p ).EQ.iwork( p+1 ) )
     $   p = p - 1
*
*     Determine block structure of B
*
      q = p + 1
      j = 1
   60 CONTINUE
      IF( j.GT.n )
     $   GO TO 70
      q = q + 1
      iwork( q ) = j
      j = j + nb
      IF( j.GE.n )
     $   GO TO 70
      IF( b( j, j-1 ).NE.zero )
     $   j = j + 1
      GO TO 60
   70 CONTINUE
*
      iwork( q+1 ) = n + 1
      IF( iwork( q ).EQ.iwork( q+1 ) )
     $   q = q - 1
*
      IF( notran ) THEN
*
         DO 150 iround = 1, isolve
*
*           Solve (I, J)-subsystem
*               A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
*               D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
*           for I = P, P - 1,..., 1; J = 1, 2,..., Q
*
            dscale = zero
            dsum = one
            pq = 0
            scale = one
            DO 130 j = p + 2, q
               js = iwork( j )
               je = iwork( j+1 ) - 1
               nb = je - js + 1
               DO 120 i = p, 1, -1
                  is = iwork( i )
                  ie = iwork( i+1 ) - 1
                  mb = ie - is + 1
                  ppqq = 0
                  CALL dtgsy2( trans, ifunc, mb, nb, a( is, is ),
     $                         lda,
     $                         b( js, js ), ldb, c( is, js ), ldc,
     $                         d( is, is ), ldd, e( js, js ), lde,
     $                         f( is, js ), ldf, scaloc, dsum, dscale,
     $                         iwork( q+2 ), ppqq, linfo )
                  IF( linfo.GT.0 )
     $               info = linfo
*
                  pq = pq + ppqq
                  IF( scaloc.NE.one ) THEN
                     DO 80 k = 1, js - 1
                        CALL dscal( m, scaloc, c( 1, k ), 1 )
                        CALL dscal( m, scaloc, f( 1, k ), 1 )
   80                CONTINUE
                     DO 90 k = js, je
                        CALL dscal( is-1, scaloc, c( 1, k ), 1 )
                        CALL dscal( is-1, scaloc, f( 1, k ), 1 )
   90                CONTINUE
                     DO 100 k = js, je
                        CALL dscal( m-ie, scaloc, c( ie+1, k ), 1 )
                        CALL dscal( m-ie, scaloc, f( ie+1, k ), 1 )
  100                CONTINUE
                     DO 110 k = je + 1, n
                        CALL dscal( m, scaloc, c( 1, k ), 1 )
                        CALL dscal( m, scaloc, f( 1, k ), 1 )
  110                CONTINUE
                     scale = scale*scaloc
                  END IF
*
*                 Substitute R(I, J) and L(I, J) into remaining
*                 equation.
*
                  IF( i.GT.1 ) THEN
                     CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
     $                           a( 1, is ), lda, c( is, js ), ldc, one,
     $                           c( 1, js ), ldc )
                     CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
     $                           d( 1, is ), ldd, c( is, js ), ldc, one,
     $                           f( 1, js ), ldf )
                  END IF
                  IF( j.LT.q ) THEN
                     CALL dgemm( 'N', 'N', mb, n-je, nb, one,
     $                           f( is, js ), ldf, b( js, je+1 ), ldb,
     $                           one, c( is, je+1 ), ldc )
                     CALL dgemm( 'N', 'N', mb, n-je, nb, one,
     $                           f( is, js ), ldf, e( js, je+1 ), lde,
     $                           one, f( is, je+1 ), ldf )
                  END IF
  120          CONTINUE
  130       CONTINUE
            IF( dscale.NE.zero ) THEN
               IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
                  dif = sqrt( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
               ELSE
                  dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
               END IF
            END IF
            IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
               IF( notran ) THEN
                  ifunc = ijob
               END IF
               scale2 = scale
               CALL dlacpy( 'F', m, n, c, ldc, work, m )
               CALL dlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
               CALL dlaset( 'F', m, n, zero, zero, c, ldc )
               CALL dlaset( 'F', m, n, zero, zero, f, ldf )
            ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
               CALL dlacpy( 'F', m, n, work, m, c, ldc )
               CALL dlacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
               scale = scale2
            END IF
  150    CONTINUE
*
      ELSE
*
*        Solve transposed (I, J)-subsystem
*             A(I, I)**T * R(I, J)  + D(I, I)**T * L(I, J)  =  C(I, J)
*             R(I, J)  * B(J, J)**T + L(I, J)  * E(J, J)**T = -F(I, J)
*        for I = 1,2,..., P; J = Q, Q-1,..., 1
*
         scale = one
         DO 210 i = 1, p
            is = iwork( i )
            ie = iwork( i+1 ) - 1
            mb = ie - is + 1
            DO 200 j = q, p + 2, -1
               js = iwork( j )
               je = iwork( j+1 ) - 1
               nb = je - js + 1
               CALL dtgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
     $                      b( js, js ), ldb, c( is, js ), ldc,
     $                      d( is, is ), ldd, e( js, js ), lde,
     $                      f( is, js ), ldf, scaloc, dsum, dscale,
     $                      iwork( q+2 ), ppqq, linfo )
               IF( linfo.GT.0 )
     $            info = linfo
               IF( scaloc.NE.one ) THEN
                  DO 160 k = 1, js - 1
                     CALL dscal( m, scaloc, c( 1, k ), 1 )
                     CALL dscal( m, scaloc, f( 1, k ), 1 )
  160             CONTINUE
                  DO 170 k = js, je
                     CALL dscal( is-1, scaloc, c( 1, k ), 1 )
                     CALL dscal( is-1, scaloc, f( 1, k ), 1 )
  170             CONTINUE
                  DO 180 k = js, je
                     CALL dscal( m-ie, scaloc, c( ie+1, k ), 1 )
                     CALL dscal( m-ie, scaloc, f( ie+1, k ), 1 )
  180             CONTINUE
                  DO 190 k = je + 1, n
                     CALL dscal( m, scaloc, c( 1, k ), 1 )
                     CALL dscal( m, scaloc, f( 1, k ), 1 )
  190             CONTINUE
                  scale = scale*scaloc
               END IF
*
*              Substitute R(I, J) and L(I, J) into remaining equation.
*
               IF( j.GT.p+2 ) THEN
                  CALL dgemm( 'N', 'T', mb, js-1, nb, one, c( is,
     $                        js ),
     $                        ldc, b( 1, js ), ldb, one, f( is, 1 ),
     $                        ldf )
                  CALL dgemm( 'N', 'T', mb, js-1, nb, one, f( is,
     $                        js ),
     $                        ldf, e( 1, js ), lde, one, f( is, 1 ),
     $                        ldf )
               END IF
               IF( i.LT.p ) THEN
                  CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
     $                        a( is, ie+1 ), lda, c( is, js ), ldc, one,
     $                        c( ie+1, js ), ldc )
                  CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
     $                        d( is, ie+1 ), ldd, f( is, js ), ldf, one,
     $                        c( ie+1, js ), ldc )
               END IF
  200       CONTINUE
  210    CONTINUE
*
      END IF
*
      work( 1 ) = lwmin
*
      RETURN
*
*     End of DTGSYL
*
      END