LAPACK: SRC/ctgsyl.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
00002      $                   LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
00003      $                   IWORK, INFO )
00004 *
00005 *  -- LAPACK routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     January 2007
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          TRANS
00012       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
00013      $                   LWORK, M, N
00014       REAL               DIF, SCALE
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            IWORK( * )
00018       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
00019      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
00020      $                   WORK( * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  CTGSYL solves the generalized Sylvester equation:
00027 *
00028 *              A * R - L * B = scale * C            (1)
00029 *              D * R - L * E = scale * F
00030 *
00031 *  where R and L are unknown m-by-n matrices, (A, D), (B, E) and
00032 *  (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
00033 *  respectively, with complex entries. A, B, D and E are upper
00034 *  triangular (i.e., (A,D) and (B,E) in generalized Schur form).
00035 *
00036 *  The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1
00037 *  is an output scaling factor chosen to avoid overflow.
00038 *
00039 *  In matrix notation (1) is equivalent to solve Zx = scale*b, where Z
00040 *  is defined as
00041 *
00042 *         Z = [ kron(In, A)  -kron(B', Im) ]        (2)
00043 *             [ kron(In, D)  -kron(E', Im) ],
00044 *
00045 *  Here Ix is the identity matrix of size x and X' is the conjugate
00046 *  transpose of X. Kron(X, Y) is the Kronecker product between the
00047 *  matrices X and Y.
00048 *
00049 *  If TRANS = 'C', y in the conjugate transposed system Z'*y = scale*b
00050 *  is solved for, which is equivalent to solve for R and L in
00051 *
00052 *              A' * R + D' * L = scale * C           (3)
00053 *              R * B' + L * E' = scale * -F
00054 *
00055 *  This case (TRANS = 'C') is used to compute an one-norm-based estimate
00056 *  of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
00057 *  and (B,E), using CLACON.
00058 *
00059 *  If IJOB >= 1, CTGSYL computes a Frobenius norm-based estimate of
00060 *  Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
00061 *  reciprocal of the smallest singular value of Z.
00062 *
00063 *  This is a level-3 BLAS algorithm.
00064 *
00065 *  Arguments
00066 *  =========
00067 *
00068 *  TRANS   (input) CHARACTER*1
00069 *          = 'N': solve the generalized sylvester equation (1).
00070 *          = 'C': solve the "conjugate transposed" system (3).
00071 *
00072 *  IJOB    (input) INTEGER
00073 *          Specifies what kind of functionality to be performed.
00074 *          =0: solve (1) only.
00075 *          =1: The functionality of 0 and 3.
00076 *          =2: The functionality of 0 and 4.
00077 *          =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
00078 *              (look ahead strategy is used).
00079 *          =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
00080 *              (CGECON on sub-systems is used).
00081 *          Not referenced if TRANS = 'C'.
00082 *
00083 *  M       (input) INTEGER
00084 *          The order of the matrices A and D, and the row dimension of
00085 *          the matrices C, F, R and L.
00086 *
00087 *  N       (input) INTEGER
00088 *          The order of the matrices B and E, and the column dimension
00089 *          of the matrices C, F, R and L.
00090 *
00091 *  A       (input) COMPLEX array, dimension (LDA, M)
00092 *          The upper triangular matrix A.
00093 *
00094 *  LDA     (input) INTEGER
00095 *          The leading dimension of the array A. LDA >= max(1, M).
00096 *
00097 *  B       (input) COMPLEX array, dimension (LDB, N)
00098 *          The upper triangular matrix B.
00099 *
00100 *  LDB     (input) INTEGER
00101 *          The leading dimension of the array B. LDB >= max(1, N).
00102 *
00103 *  C       (input/output) COMPLEX array, dimension (LDC, N)
00104 *          On entry, C contains the right-hand-side of the first matrix
00105 *          equation in (1) or (3).
00106 *          On exit, if IJOB = 0, 1 or 2, C has been overwritten by
00107 *          the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R,
00108 *          the solution achieved during the computation of the
00109 *          Dif-estimate.
00110 *
00111 *  LDC     (input) INTEGER
00112 *          The leading dimension of the array C. LDC >= max(1, M).
00113 *
00114 *  D       (input) COMPLEX array, dimension (LDD, M)
00115 *          The upper triangular matrix D.
00116 *
00117 *  LDD     (input) INTEGER
00118 *          The leading dimension of the array D. LDD >= max(1, M).
00119 *
00120 *  E       (input) COMPLEX array, dimension (LDE, N)
00121 *          The upper triangular matrix E.
00122 *
00123 *  LDE     (input) INTEGER
00124 *          The leading dimension of the array E. LDE >= max(1, N).
00125 *
00126 *  F       (input/output) COMPLEX array, dimension (LDF, N)
00127 *          On entry, F contains the right-hand-side of the second matrix
00128 *          equation in (1) or (3).
00129 *          On exit, if IJOB = 0, 1 or 2, F has been overwritten by
00130 *          the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L,
00131 *          the solution achieved during the computation of the
00132 *          Dif-estimate.
00133 *
00134 *  LDF     (input) INTEGER
00135 *          The leading dimension of the array F. LDF >= max(1, M).
00136 *
00137 *  DIF     (output) REAL
00138 *          On exit DIF is the reciprocal of a lower bound of the
00139 *          reciprocal of the Dif-function, i.e. DIF is an upper bound of
00140 *          Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2).
00141 *          IF IJOB = 0 or TRANS = 'C', DIF is not referenced.
00142 *
00143 *  SCALE   (output) REAL
00144 *          On exit SCALE is the scaling factor in (1) or (3).
00145 *          If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
00146 *          to a slightly perturbed system but the input matrices A, B,
00147 *          D and E have not been changed. If SCALE = 0, R and L will
00148 *          hold the solutions to the homogenious system with C = F = 0.
00149 *
00150 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00151 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00152 *
00153 *  LWORK   (input) INTEGER
00154 *          The dimension of the array WORK. LWORK > = 1.
00155 *          If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
00156 *
00157 *          If LWORK = -1, then a workspace query is assumed; the routine
00158 *          only calculates the optimal size of the WORK array, returns
00159 *          this value as the first entry of the WORK array, and no error
00160 *          message related to LWORK is issued by XERBLA.
00161 *
00162 *  IWORK   (workspace) INTEGER array, dimension (M+N+2)
00163 *
00164 *  INFO    (output) INTEGER
00165 *            =0: successful exit
00166 *            <0: If INFO = -i, the i-th argument had an illegal value.
00167 *            >0: (A, D) and (B, E) have common or very close
00168 *                eigenvalues.
00169 *
00170 *  Further Details
00171 *  ===============
00172 *
00173 *  Based on contributions by
00174 *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00175 *     Umea University, S-901 87 Umea, Sweden.
00176 *
00177 *  [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
00178 *      for Solving the Generalized Sylvester Equation and Estimating the
00179 *      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
00180 *      Department of Computing Science, Umea University, S-901 87 Umea,
00181 *      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
00182 *      Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
00183 *      No 1, 1996.
00184 *
00185 *  [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
00186 *      Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
00187 *      Appl., 15(4):1045-1060, 1994.
00188 *
00189 *  [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
00190 *      Condition Estimators for Solving the Generalized Sylvester
00191 *      Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
00192 *      July 1989, pp 745-751.
00193 *
00194 *  =====================================================================
00195 *  Replaced various illegal calls to CCOPY by calls to CLASET.
00196 *  Sven Hammarling, 1/5/02.
00197 *
00198 *     .. Parameters ..
00199       REAL               ZERO, ONE
00200       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00201       COMPLEX            CZERO
00202       PARAMETER          ( CZERO = (0.0E+0, 0.0E+0) )
00203 *     ..
00204 *     .. Local Scalars ..
00205       LOGICAL            LQUERY, NOTRAN
00206       INTEGER            I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
00207      $                   LINFO, LWMIN, MB, NB, P, PQ, Q
00208       REAL               DSCALE, DSUM, SCALE2, SCALOC
00209 *     ..
00210 *     .. External Functions ..
00211       LOGICAL            LSAME
00212       INTEGER            ILAENV
00213       EXTERNAL           LSAME, ILAENV
00214 *     ..
00215 *     .. External Subroutines ..
00216       EXTERNAL           CGEMM, CLACPY, CLASET, CSCAL, CTGSY2, XERBLA
00217 *     ..
00218 *     .. Intrinsic Functions ..
00219       INTRINSIC          CMPLX, MAX, REAL, SQRT
00220 *     ..
00221 *     .. Executable Statements ..
00222 *
00223 *     Decode and test input parameters
00224 *
00225       INFO = 0
00226       NOTRAN = LSAME( TRANS, 'N' )
00227       LQUERY = ( LWORK.EQ.-1 )
00228 *
00229       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
00230          INFO = -1
00231       ELSE IF( NOTRAN ) THEN
00232          IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN
00233             INFO = -2
00234          END IF
00235       END IF
00236       IF( INFO.EQ.0 ) THEN
00237          IF( M.LE.0 ) THEN
00238             INFO = -3
00239          ELSE IF( N.LE.0 ) THEN
00240             INFO = -4
00241          ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00242             INFO = -6
00243          ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00244             INFO = -8
00245          ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
00246             INFO = -10
00247          ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
00248             INFO = -12
00249          ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
00250             INFO = -14
00251          ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
00252             INFO = -16
00253          END IF
00254       END IF
00255 *
00256       IF( INFO.EQ.0 ) THEN
00257          IF( NOTRAN ) THEN
00258             IF( IJOB.EQ.1 .OR. IJOB.EQ.2 ) THEN
00259                LWMIN = MAX( 1, 2*M*N )
00260             ELSE
00261                LWMIN = 1
00262             END IF
00263          ELSE
00264             LWMIN = 1
00265          END IF
00266          WORK( 1 ) = LWMIN
00267 *
00268          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00269             INFO = -20
00270          END IF
00271       END IF
00272 *
00273       IF( INFO.NE.0 ) THEN
00274          CALL XERBLA( 'CTGSYL', -INFO )
00275          RETURN
00276       ELSE IF( LQUERY ) THEN
00277          RETURN
00278       END IF
00279 *
00280 *     Quick return if possible
00281 *
00282       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00283          SCALE = 1
00284          IF( NOTRAN ) THEN
00285             IF( IJOB.NE.0 ) THEN
00286                DIF = 0
00287             END IF
00288          END IF
00289          RETURN
00290       END IF
00291 *
00292 *     Determine  optimal block sizes MB and NB
00293 *
00294       MB = ILAENV( 2, 'CTGSYL', TRANS, M, N, -1, -1 )
00295       NB = ILAENV( 5, 'CTGSYL', TRANS, M, N, -1, -1 )
00296 *
00297       ISOLVE = 1
00298       IFUNC = 0
00299       IF( NOTRAN ) THEN
00300          IF( IJOB.GE.3 ) THEN
00301             IFUNC = IJOB - 2
00302             CALL CLASET( 'F', M, N, CZERO, CZERO, C, LDC )
00303             CALL CLASET( 'F', M, N, CZERO, CZERO, F, LDF )
00304          ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN
00305             ISOLVE = 2
00306          END IF
00307       END IF
00308 *
00309       IF( ( MB.LE.1 .AND. NB.LE.1 ) .OR. ( MB.GE.M .AND. NB.GE.N ) )
00310      $     THEN
00311 *
00312 *        Use unblocked Level 2 solver
00313 *
00314          DO 30 IROUND = 1, ISOLVE
00315 *
00316             SCALE = ONE
00317             DSCALE = ZERO
00318             DSUM = ONE
00319             PQ = M*N
00320             CALL CTGSY2( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D,
00321      $                   LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE,
00322      $                   INFO )
00323             IF( DSCALE.NE.ZERO ) THEN
00324                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
00325                   DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
00326                ELSE
00327                   DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
00328                END IF
00329             END IF
00330             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
00331                IF( NOTRAN ) THEN
00332                   IFUNC = IJOB
00333                END IF
00334                SCALE2 = SCALE
00335                CALL CLACPY( 'F', M, N, C, LDC, WORK, M )
00336                CALL CLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
00337                CALL CLASET( 'F', M, N, CZERO, CZERO, C, LDC )
00338                CALL CLASET( 'F', M, N, CZERO, CZERO, F, LDF )
00339             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
00340                CALL CLACPY( 'F', M, N, WORK, M, C, LDC )
00341                CALL CLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
00342                SCALE = SCALE2
00343             END IF
00344    30    CONTINUE
00345 *
00346          RETURN
00347 *
00348       END IF
00349 *
00350 *     Determine block structure of A
00351 *
00352       P = 0
00353       I = 1
00354    40 CONTINUE
00355       IF( I.GT.M )
00356      $   GO TO 50
00357       P = P + 1
00358       IWORK( P ) = I
00359       I = I + MB
00360       IF( I.GE.M )
00361      $   GO TO 50
00362       GO TO 40
00363    50 CONTINUE
00364       IWORK( P+1 ) = M + 1
00365       IF( IWORK( P ).EQ.IWORK( P+1 ) )
00366      $   P = P - 1
00367 *
00368 *     Determine block structure of B
00369 *
00370       Q = P + 1
00371       J = 1
00372    60 CONTINUE
00373       IF( J.GT.N )
00374      $   GO TO 70
00375 *
00376       Q = Q + 1
00377       IWORK( Q ) = J
00378       J = J + NB
00379       IF( J.GE.N )
00380      $   GO TO 70
00381       GO TO 60
00382 *
00383    70 CONTINUE
00384       IWORK( Q+1 ) = N + 1
00385       IF( IWORK( Q ).EQ.IWORK( Q+1 ) )
00386      $   Q = Q - 1
00387 *
00388       IF( NOTRAN ) THEN
00389          DO 150 IROUND = 1, ISOLVE
00390 *
00391 *           Solve (I, J) - subsystem
00392 *               A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
00393 *               D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
00394 *           for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
00395 *
00396             PQ = 0
00397             SCALE = ONE
00398             DSCALE = ZERO
00399             DSUM = ONE
00400             DO 130 J = P + 2, Q
00401                JS = IWORK( J )
00402                JE = IWORK( J+1 ) - 1
00403                NB = JE - JS + 1
00404                DO 120 I = P, 1, -1
00405                   IS = IWORK( I )
00406                   IE = IWORK( I+1 ) - 1
00407                   MB = IE - IS + 1
00408                   CALL CTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
00409      $                         B( JS, JS ), LDB, C( IS, JS ), LDC,
00410      $                         D( IS, IS ), LDD, E( JS, JS ), LDE,
00411      $                         F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
00412      $                         LINFO )
00413                   IF( LINFO.GT.0 )
00414      $               INFO = LINFO
00415                   PQ = PQ + MB*NB
00416                   IF( SCALOC.NE.ONE ) THEN
00417                      DO 80 K = 1, JS - 1
00418                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
00419      $                              1 )
00420                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
00421      $                              1 )
00422    80                CONTINUE
00423                      DO 90 K = JS, JE
00424                         CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ),
00425      $                              C( 1, K ), 1 )
00426                         CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ),
00427      $                              F( 1, K ), 1 )
00428    90                CONTINUE
00429                      DO 100 K = JS, JE
00430                         CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ),
00431      $                              C( IE+1, K ), 1 )
00432                         CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ),
00433      $                              F( IE+1, K ), 1 )
00434   100                CONTINUE
00435                      DO 110 K = JE + 1, N
00436                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
00437      $                              1 )
00438                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
00439      $                              1 )
00440   110                CONTINUE
00441                      SCALE = SCALE*SCALOC
00442                   END IF
00443 *
00444 *                 Substitute R(I,J) and L(I,J) into remaining equation.
00445 *
00446                   IF( I.GT.1 ) THEN
00447                      CALL CGEMM( 'N', 'N', IS-1, NB, MB,
00448      $                           CMPLX( -ONE, ZERO ), A( 1, IS ), LDA,
00449      $                           C( IS, JS ), LDC, CMPLX( ONE, ZERO ),
00450      $                           C( 1, JS ), LDC )
00451                      CALL CGEMM( 'N', 'N', IS-1, NB, MB,
00452      $                           CMPLX( -ONE, ZERO ), D( 1, IS ), LDD,
00453      $                           C( IS, JS ), LDC, CMPLX( ONE, ZERO ),
00454      $                           F( 1, JS ), LDF )
00455                   END IF
00456                   IF( J.LT.Q ) THEN
00457                      CALL CGEMM( 'N', 'N', MB, N-JE, NB,
00458      $                           CMPLX( ONE, ZERO ), F( IS, JS ), LDF,
00459      $                           B( JS, JE+1 ), LDB, CMPLX( ONE, ZERO ),
00460      $                           C( IS, JE+1 ), LDC )
00461                      CALL CGEMM( 'N', 'N', MB, N-JE, NB,
00462      $                           CMPLX( ONE, ZERO ), F( IS, JS ), LDF,
00463      $                           E( JS, JE+1 ), LDE, CMPLX( ONE, ZERO ),
00464      $                           F( IS, JE+1 ), LDF )
00465                   END IF
00466   120          CONTINUE
00467   130       CONTINUE
00468             IF( DSCALE.NE.ZERO ) THEN
00469                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
00470                   DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
00471                ELSE
00472                   DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
00473                END IF
00474             END IF
00475             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
00476                IF( NOTRAN ) THEN
00477                   IFUNC = IJOB
00478                END IF
00479                SCALE2 = SCALE
00480                CALL CLACPY( 'F', M, N, C, LDC, WORK, M )
00481                CALL CLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
00482                CALL CLASET( 'F', M, N, CZERO, CZERO, C, LDC )
00483                CALL CLASET( 'F', M, N, CZERO, CZERO, F, LDF )
00484             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
00485                CALL CLACPY( 'F', M, N, WORK, M, C, LDC )
00486                CALL CLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
00487                SCALE = SCALE2
00488             END IF
00489   150    CONTINUE
00490       ELSE
00491 *
00492 *        Solve transposed (I, J)-subsystem
00493 *            A(I, I)' * R(I, J) + D(I, I)' * L(I, J) = C(I, J)
00494 *            R(I, J) * B(J, J)  + L(I, J) * E(J, J) = -F(I, J)
00495 *        for I = 1,2,..., P; J = Q, Q-1,..., 1
00496 *
00497          SCALE = ONE
00498          DO 210 I = 1, P
00499             IS = IWORK( I )
00500             IE = IWORK( I+1 ) - 1
00501             MB = IE - IS + 1
00502             DO 200 J = Q, P + 2, -1
00503                JS = IWORK( J )
00504                JE = IWORK( J+1 ) - 1
00505                NB = JE - JS + 1
00506                CALL CTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
00507      $                      B( JS, JS ), LDB, C( IS, JS ), LDC,
00508      $                      D( IS, IS ), LDD, E( JS, JS ), LDE,
00509      $                      F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
00510      $                      LINFO )
00511                IF( LINFO.GT.0 )
00512      $            INFO = LINFO
00513                IF( SCALOC.NE.ONE ) THEN
00514                   DO 160 K = 1, JS - 1
00515                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
00516      $                           1 )
00517                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
00518      $                           1 )
00519   160             CONTINUE
00520                   DO 170 K = JS, JE
00521                      CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ), C( 1, K ),
00522      $                           1 )
00523                      CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ), F( 1, K ),
00524      $                           1 )
00525   170             CONTINUE
00526                   DO 180 K = JS, JE
00527                      CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ),
00528      $                           C( IE+1, K ), 1 )
00529                      CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ),
00530      $                           F( IE+1, K ), 1 )
00531   180             CONTINUE
00532                   DO 190 K = JE + 1, N
00533                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
00534      $                           1 )
00535                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
00536      $                           1 )
00537   190             CONTINUE
00538                   SCALE = SCALE*SCALOC
00539                END IF
00540 *
00541 *              Substitute R(I,J) and L(I,J) into remaining equation.
00542 *
00543                IF( J.GT.P+2 ) THEN
00544                   CALL CGEMM( 'N', 'C', MB, JS-1, NB,
00545      $                        CMPLX( ONE, ZERO ), C( IS, JS ), LDC,
00546      $                        B( 1, JS ), LDB, CMPLX( ONE, ZERO ),
00547      $                        F( IS, 1 ), LDF )
00548                   CALL CGEMM( 'N', 'C', MB, JS-1, NB,
00549      $                        CMPLX( ONE, ZERO ), F( IS, JS ), LDF,
00550      $                        E( 1, JS ), LDE, CMPLX( ONE, ZERO ),
00551      $                        F( IS, 1 ), LDF )
00552                END IF
00553                IF( I.LT.P ) THEN
00554                   CALL CGEMM( 'C', 'N', M-IE, NB, MB,
00555      $                        CMPLX( -ONE, ZERO ), A( IS, IE+1 ), LDA,
00556      $                        C( IS, JS ), LDC, CMPLX( ONE, ZERO ),
00557      $                        C( IE+1, JS ), LDC )
00558                   CALL CGEMM( 'C', 'N', M-IE, NB, MB,
00559      $                        CMPLX( -ONE, ZERO ), D( IS, IE+1 ), LDD,
00560      $                        F( IS, JS ), LDF, CMPLX( ONE, ZERO ),
00561      $                        C( IE+1, JS ), LDC )
00562                END IF
00563   200       CONTINUE
00564   210    CONTINUE
00565       END IF
00566 *
00567       WORK( 1 ) = LWMIN
00568 *
00569       RETURN
00570 *
00571 *     End of CTGSYL
00572 *
00573       END