LAPACK: SRC/dtgsyl.f Source File

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