ctgsyl.f

*> \brief \b CTGSYL
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CTGSYL + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsyl.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsyl.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsyl.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE CTGSYL( 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
*       REAL               DIF, SCALE
*       ..
*       .. Array Arguments ..
*       INTEGER            IWORK( * )
*       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
*      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
*      $                   WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CTGSYL 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 complex entries. A, B, D and E are upper
*> triangular (i.e., (A,D) and (B,E) in generalized Schur form).
*>
*> 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**H, Im) ]        (2)
*>            [ kron(In, D)  -kron(E**H, Im) ],
*>
*> Here Ix is the identity matrix of size x and X**H is the conjugate
*> transpose of X. Kron(X, Y) is the Kronecker product between the
*> matrices X and Y.
*>
*> If TRANS = 'C', y in the conjugate transposed system Z**H *y = scale*b
*> is solved for, which is equivalent to solve for R and L in
*>
*>             A**H * R + D**H * L = scale * C           (3)
*>             R * B**H + L * E**H = scale * -F
*>
*> This case (TRANS = 'C') 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 CLACON.
*>
*> If IJOB >= 1, CTGSYL 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.
*>
*> This is a level-3 BLAS algorithm.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] TRANS
*> \verbatim
*>          TRANS is CHARACTER*1
*>          = 'N': solve the generalized sylvester equation (1).
*>          = 'C': solve the "conjugate 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 is used).
*>          =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
*>              (CGECON on sub-systems is used).
*>          Not referenced if TRANS = 'C'.
*> \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 COMPLEX array, dimension (LDA, M)
*>          The upper 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 COMPLEX array, dimension (LDB, N)
*>          The upper 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 REAL
*>          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 = 'C', DIF is not referenced.
*> \endverbatim
*>
*> \param[out] SCALE
*> \verbatim
*>          SCALE is REAL
*>          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, R and L will
*>          hold the solutions to the homogeneous system with C = F = 0.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX 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+2)
*> \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 very 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:
*  ================
*>
*>  [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.
*> \n
*>  [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.
*> \n
*>  [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.
*>
*  =====================================================================
      SUBROUTINE ctgsyl( 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
      REAL               DIF, SCALE
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      COMPLEX            A( lda, * ), B( ldb, * ), C( ldc, * ),
     $                   d( ldd, * ), e( lde, * ), f( ldf, * ),
     $                   work( * )
*     ..
*
*  =====================================================================
*  Replaced various illegal calls to CCOPY by calls to CLASET.
*  Sven Hammarling, 1/5/02.
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
      COMPLEX            CZERO
      parameter( czero = (0.0e+0, 0.0e+0) )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, NOTRAN
      INTEGER            I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
     $                   linfo, lwmin, mb, nb, p, pq, q
      REAL               DSCALE, DSUM, SCALE2, SCALOC
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      EXTERNAL           lsame, ilaenv
*     ..
*     .. External Subroutines ..
      EXTERNAL           cgemm, clacpy, claset, cscal, ctgsy2, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          cmplx, max, REAL, 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, 'C' ) ) 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( 'CTGSYL', -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, 'CTGSYL', trans, m, n, -1, -1 )
      nb = ilaenv( 5, 'CTGSYL', trans, m, n, -1, -1 )
*
      isolve = 1
      ifunc = 0
      IF( notran ) THEN
         IF( ijob.GE.3 ) THEN
            ifunc = ijob - 2
            CALL claset( 'F', m, n, czero, czero, c, ldc )
            CALL claset( 'F', m, n, czero, czero, f, ldf )
         ELSE IF( ijob.GE.1 .AND. notran ) THEN
            isolve = 2
         END IF
      END IF
*
      IF( ( mb.LE.1 .AND. nb.LE.1 ) .OR. ( mb.GE.m .AND. nb.GE.n ) )
     $     THEN
*
*        Use unblocked Level 2 solver
*
         DO 30 iround = 1, isolve
*
            scale = one
            dscale = zero
            dsum = one
            pq = m*n
            CALL ctgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc, d,
     $                   ldd, e, lde, f, ldf, scale, dsum, dscale,
     $                   info )
            IF( dscale.NE.zero ) THEN
               IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
                  dif = sqrt( REAL( 2*M*N ) ) / ( DSCALE*SQRT( dsum ) )
               ELSE
                  dif = sqrt( REAL( 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 clacpy( 'F', m, n, c, ldc, work, m )
               CALL clacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
               CALL claset( 'F', m, n, czero, czero, c, ldc )
               CALL claset( 'F', m, n, czero, czero, f, ldf )
            ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
               CALL clacpy( 'F', m, n, work, m, c, ldc )
               CALL clacpy( '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
      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
      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
*
            pq = 0
            scale = one
            dscale = zero
            dsum = 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
                  CALL ctgsy2( 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,
     $                         linfo )
                  IF( linfo.GT.0 )
     $               info = linfo
                  pq = pq + mb*nb
                  IF( scaloc.NE.one ) THEN
                     DO 80 k = 1, js - 1
                        CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
     $                              1 )
                        CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
     $                              1 )
   80                CONTINUE
                     DO 90 k = js, je
                        CALL cscal( is-1, cmplx( scaloc, zero ),
     $                              c( 1, k ), 1 )
                        CALL cscal( is-1, cmplx( scaloc, zero ),
     $                              f( 1, k ), 1 )
   90                CONTINUE
                     DO 100 k = js, je
                        CALL cscal( m-ie, cmplx( scaloc, zero ),
     $                              c( ie+1, k ), 1 )
                        CALL cscal( m-ie, cmplx( scaloc, zero ),
     $                              f( ie+1, k ), 1 )
  100                CONTINUE
                     DO 110 k = je + 1, n
                        CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
     $                              1 )
                        CALL cscal( m, cmplx( scaloc, zero ), 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 cgemm( 'N', 'N', is-1, nb, mb,
     $                           cmplx( -one, zero ), a( 1, is ), lda,
     $                           c( is, js ), ldc, cmplx( one, zero ),
     $                           c( 1, js ), ldc )
                     CALL cgemm( 'N', 'N', is-1, nb, mb,
     $                           cmplx( -one, zero ), d( 1, is ), ldd,
     $                           c( is, js ), ldc, cmplx( one, zero ),
     $                           f( 1, js ), ldf )
                  END IF
                  IF( j.LT.q ) THEN
                     CALL cgemm( 'N', 'N', mb, n-je, nb,
     $                           cmplx( one, zero ), f( is, js ), ldf,
     $                           b( js, je+1 ), ldb, cmplx( one, zero ),
     $                           c( is, je+1 ), ldc )
                     CALL cgemm( 'N', 'N', mb, n-je, nb,
     $                           cmplx( one, zero ), f( is, js ), ldf,
     $                           e( js, je+1 ), lde, cmplx( one, zero ),
     $                           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( REAL( 2*M*N ) ) / ( DSCALE*SQRT( dsum ) )
               ELSE
                  dif = sqrt( REAL( 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 clacpy( 'F', m, n, c, ldc, work, m )
               CALL clacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
               CALL claset( 'F', m, n, czero, czero, c, ldc )
               CALL claset( 'F', m, n, czero, czero, f, ldf )
            ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
               CALL clacpy( 'F', m, n, work, m, c, ldc )
               CALL clacpy( '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)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
*            R(I, J) * B(J, J)  + L(I, J) * E(J, J) = -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 ctgsy2( 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,
     $                      linfo )
               IF( linfo.GT.0 )
     $            info = linfo
               IF( scaloc.NE.one ) THEN
                  DO 160 k = 1, js - 1
                     CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
     $                           1 )
                     CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
     $                           1 )
  160             CONTINUE
                  DO 170 k = js, je
                     CALL cscal( is-1, cmplx( scaloc, zero ), c( 1, k ),
     $                           1 )
                     CALL cscal( is-1, cmplx( scaloc, zero ), f( 1, k ),
     $                           1 )
  170             CONTINUE
                  DO 180 k = js, je
                     CALL cscal( m-ie, cmplx( scaloc, zero ),
     $                           c( ie+1, k ), 1 )
                     CALL cscal( m-ie, cmplx( scaloc, zero ),
     $                           f( ie+1, k ), 1 )
  180             CONTINUE
                  DO 190 k = je + 1, n
                     CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
     $                           1 )
                     CALL cscal( m, cmplx( scaloc, zero ), 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 cgemm( 'N', 'C', mb, js-1, nb,
     $                        cmplx( one, zero ), c( is, js ), ldc,
     $                        b( 1, js ), ldb, cmplx( one, zero ),
     $                        f( is, 1 ), ldf )
                  CALL cgemm( 'N', 'C', mb, js-1, nb,
     $                        cmplx( one, zero ), f( is, js ), ldf,
     $                        e( 1, js ), lde, cmplx( one, zero ),
     $                        f( is, 1 ), ldf )
               END IF
               IF( i.LT.p ) THEN
                  CALL cgemm( 'C', 'N', m-ie, nb, mb,
     $                        cmplx( -one, zero ), a( is, ie+1 ), lda,
     $                        c( is, js ), ldc, cmplx( one, zero ),
     $                        c( ie+1, js ), ldc )
                  CALL cgemm( 'C', 'N', m-ie, nb, mb,
     $                        cmplx( -one, zero ), d( is, ie+1 ), ldd,
     $                        f( is, js ), ldf, cmplx( one, zero ),
     $                        c( ie+1, js ), ldc )
               END IF
  200       CONTINUE
  210    CONTINUE
      END IF
*
      work( 1 ) = lwmin
*
      RETURN
*
*     End of CTGSYL
*
      END