001:       SUBROUTINE STGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
002:      $                   LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
003:      $                   IWORK, INFO )
004: *
005: *  -- LAPACK routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          TRANS
011:       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
012:      $                   LWORK, M, N
013:       REAL               DIF, SCALE
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            IWORK( * )
017:       REAL               A( LDA, * ), B( LDB, * ), C( LDC, * ),
018:      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
019:      $                   WORK( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  STGSYL solves the generalized Sylvester equation:
026: *
027: *              A * R - L * B = scale * C                 (1)
028: *              D * R - L * E = scale * F
029: *
030: *  where R and L are unknown m-by-n matrices, (A, D), (B, E) and
031: *  (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
032: *  respectively, with real entries. (A, D) and (B, E) must be in
033: *  generalized (real) Schur canonical form, i.e. A, B are upper quasi
034: *  triangular and D, E are upper triangular.
035: *
036: *  The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
037: *  scaling factor chosen to avoid overflow.
038: *
039: *  In matrix notation (1) is equivalent to solve  Zx = scale b, where
040: *  Z is defined as
041: *
042: *             Z = [ kron(In, A)  -kron(B', Im) ]         (2)
043: *                 [ kron(In, D)  -kron(E', Im) ].
044: *
045: *  Here Ik is the identity matrix of size k and X' is the transpose of
046: *  X. kron(X, Y) is the Kronecker product between the matrices X and Y.
047: *
048: *  If TRANS = 'T', STGSYL solves the transposed system Z'*y = scale*b,
049: *  which is equivalent to solve for R and L in
050: *
051: *              A' * R  + D' * L   = scale *  C           (3)
052: *              R  * B' + L  * E'  = scale * (-F)
053: *
054: *  This case (TRANS = 'T') is used to compute an one-norm-based estimate
055: *  of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
056: *  and (B,E), using SLACON.
057: *
058: *  If IJOB >= 1, STGSYL computes a Frobenius norm-based estimate
059: *  of Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
060: *  reciprocal of the smallest singular value of Z. See [1-2] for more
061: *  information.
062: *
063: *  This is a level 3 BLAS algorithm.
064: *
065: *  Arguments
066: *  =========
067: *
068: *  TRANS   (input) CHARACTER*1
069: *          = 'N', solve the generalized Sylvester equation (1).
070: *          = 'T', solve the 'transposed' system (3).
071: *
072: *  IJOB    (input) INTEGER
073: *          Specifies what kind of functionality to be performed.
074: *           =0: solve (1) only.
075: *           =1: The functionality of 0 and 3.
076: *           =2: The functionality of 0 and 4.
077: *           =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
078: *               (look ahead strategy IJOB  = 1 is used).
079: *           =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
080: *               ( SGECON on sub-systems is used ).
081: *          Not referenced if TRANS = 'T'.
082: *
083: *  M       (input) INTEGER
084: *          The order of the matrices A and D, and the row dimension of
085: *          the matrices C, F, R and L.
086: *
087: *  N       (input) INTEGER
088: *          The order of the matrices B and E, and the column dimension
089: *          of the matrices C, F, R and L.
090: *
091: *  A       (input) REAL array, dimension (LDA, M)
092: *          The upper quasi triangular matrix A.
093: *
094: *  LDA     (input) INTEGER
095: *          The leading dimension of the array A. LDA >= max(1, M).
096: *
097: *  B       (input) REAL array, dimension (LDB, N)
098: *          The upper quasi triangular matrix B.
099: *
100: *  LDB     (input) INTEGER
101: *          The leading dimension of the array B. LDB >= max(1, N).
102: *
103: *  C       (input/output) REAL array, dimension (LDC, N)
104: *          On entry, C contains the right-hand-side of the first matrix
105: *          equation in (1) or (3).
106: *          On exit, if IJOB = 0, 1 or 2, C has been overwritten by
107: *          the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R,
108: *          the solution achieved during the computation of the
109: *          Dif-estimate.
110: *
111: *  LDC     (input) INTEGER
112: *          The leading dimension of the array C. LDC >= max(1, M).
113: *
114: *  D       (input) REAL array, dimension (LDD, M)
115: *          The upper triangular matrix D.
116: *
117: *  LDD     (input) INTEGER
118: *          The leading dimension of the array D. LDD >= max(1, M).
119: *
120: *  E       (input) REAL array, dimension (LDE, N)
121: *          The upper triangular matrix E.
122: *
123: *  LDE     (input) INTEGER
124: *          The leading dimension of the array E. LDE >= max(1, N).
125: *
126: *  F       (input/output) REAL array, dimension (LDF, N)
127: *          On entry, F contains the right-hand-side of the second matrix
128: *          equation in (1) or (3).
129: *          On exit, if IJOB = 0, 1 or 2, F has been overwritten by
130: *          the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L,
131: *          the solution achieved during the computation of the
132: *          Dif-estimate.
133: *
134: *  LDF     (input) INTEGER
135: *          The leading dimension of the array F. LDF >= max(1, M).
136: *
137: *  DIF     (output) REAL
138: *          On exit DIF is the reciprocal of a lower bound of the
139: *          reciprocal of the Dif-function, i.e. DIF is an upper bound of
140: *          Dif[(A,D), (B,E)] = sigma_min(Z), where Z as in (2).
141: *          IF IJOB = 0 or TRANS = 'T', DIF is not touched.
142: *
143: *  SCALE   (output) REAL
144: *          On exit SCALE is the scaling factor in (1) or (3).
145: *          If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
146: *          to a slightly perturbed system but the input matrices A, B, D
147: *          and E have not been changed. If SCALE = 0, C and F hold the
148: *          solutions R and L, respectively, to the homogeneous system
149: *          with C = F = 0. Normally, SCALE = 1.
150: *
151: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
152: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
153: *
154: *  LWORK   (input) INTEGER
155: *          The dimension of the array WORK. LWORK > = 1.
156: *          If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
157: *
158: *          If LWORK = -1, then a workspace query is assumed; the routine
159: *          only calculates the optimal size of the WORK array, returns
160: *          this value as the first entry of the WORK array, and no error
161: *          message related to LWORK is issued by XERBLA.
162: *
163: *  IWORK   (workspace) INTEGER array, dimension (M+N+6)
164: *
165: *  INFO    (output) INTEGER
166: *            =0: successful exit
167: *            <0: If INFO = -i, the i-th argument had an illegal value.
168: *            >0: (A, D) and (B, E) have common or close eigenvalues.
169: *
170: *  Further Details
171: *  ===============
172: *
173: *  Based on contributions by
174: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
175: *     Umea University, S-901 87 Umea, Sweden.
176: *
177: *  [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
178: *      for Solving the Generalized Sylvester Equation and Estimating the
179: *      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
180: *      Department of Computing Science, Umea University, S-901 87 Umea,
181: *      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
182: *      Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
183: *      No 1, 1996.
184: *
185: *  [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
186: *      Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
187: *      Appl., 15(4):1045-1060, 1994
188: *
189: *  [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
190: *      Condition Estimators for Solving the Generalized Sylvester
191: *      Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
192: *      July 1989, pp 745-751.
193: *
194: *  =====================================================================
195: *  Replaced various illegal calls to SCOPY by calls to SLASET.
196: *  Sven Hammarling, 1/5/02.
197: *
198: *     .. Parameters ..
199:       REAL               ZERO, ONE
200:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
201: *     ..
202: *     .. Local Scalars ..
203:       LOGICAL            LQUERY, NOTRAN
204:       INTEGER            I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
205:      $                   LINFO, LWMIN, MB, NB, P, PPQQ, PQ, Q
206:       REAL               DSCALE, DSUM, SCALE2, SCALOC
207: *     ..
208: *     .. External Functions ..
209:       LOGICAL            LSAME
210:       INTEGER            ILAENV
211:       EXTERNAL           LSAME, ILAENV
212: *     ..
213: *     .. External Subroutines ..
214:       EXTERNAL           SGEMM, SLACPY, SLASET, SSCAL, STGSY2, XERBLA
215: *     ..
216: *     .. Intrinsic Functions ..
217:       INTRINSIC          MAX, REAL, SQRT
218: *     ..
219: *     .. Executable Statements ..
220: *
221: *     Decode and test input parameters
222: *
223:       INFO = 0
224:       NOTRAN = LSAME( TRANS, 'N' )
225:       LQUERY = ( LWORK.EQ.-1 )
226: *
227:       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
228:          INFO = -1
229:       ELSE IF( NOTRAN ) THEN
230:          IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN
231:             INFO = -2
232:          END IF
233:       END IF
234:       IF( INFO.EQ.0 ) THEN
235:          IF( M.LE.0 ) THEN
236:             INFO = -3
237:          ELSE IF( N.LE.0 ) THEN
238:             INFO = -4
239:          ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
240:             INFO = -6
241:          ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
242:             INFO = -8
243:          ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
244:             INFO = -10
245:          ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
246:             INFO = -12
247:          ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
248:             INFO = -14
249:          ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
250:             INFO = -16
251:          END IF
252:       END IF
253: *
254:       IF( INFO.EQ.0 ) THEN
255:          IF( NOTRAN ) THEN
256:             IF( IJOB.EQ.1 .OR. IJOB.EQ.2 ) THEN
257:                LWMIN = MAX( 1, 2*M*N )
258:             ELSE
259:                LWMIN = 1
260:             END IF
261:          ELSE
262:             LWMIN = 1
263:          END IF
264:          WORK( 1 ) = LWMIN
265: *
266:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
267:             INFO = -20
268:          END IF
269:       END IF
270: *
271:       IF( INFO.NE.0 ) THEN
272:          CALL XERBLA( 'STGSYL', -INFO )
273:          RETURN
274:       ELSE IF( LQUERY ) THEN
275:          RETURN
276:       END IF
277: *
278: *     Quick return if possible
279: *
280:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
281:          SCALE = 1
282:          IF( NOTRAN ) THEN
283:             IF( IJOB.NE.0 ) THEN
284:                DIF = 0
285:             END IF
286:          END IF
287:          RETURN
288:       END IF
289: *
290: *     Determine optimal block sizes MB and NB
291: *
292:       MB = ILAENV( 2, 'STGSYL', TRANS, M, N, -1, -1 )
293:       NB = ILAENV( 5, 'STGSYL', TRANS, M, N, -1, -1 )
294: *
295:       ISOLVE = 1
296:       IFUNC = 0
297:       IF( NOTRAN ) THEN
298:          IF( IJOB.GE.3 ) THEN
299:             IFUNC = IJOB - 2
300:             CALL SLASET( 'F', M, N, ZERO, ZERO, C, LDC )
301:             CALL SLASET( 'F', M, N, ZERO, ZERO, F, LDF )
302:          ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN
303:             ISOLVE = 2
304:          END IF
305:       END IF
306: *
307:       IF( ( MB.LE.1 .AND. NB.LE.1 ) .OR. ( MB.GE.M .AND. NB.GE.N ) )
308:      $     THEN
309: *
310:          DO 30 IROUND = 1, ISOLVE
311: *
312: *           Use unblocked Level 2 solver
313: *
314:             DSCALE = ZERO
315:             DSUM = ONE
316:             PQ = 0
317:             CALL STGSY2( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D,
318:      $                   LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE,
319:      $                   IWORK, PQ, INFO )
320:             IF( DSCALE.NE.ZERO ) THEN
321:                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
322:                   DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
323:                ELSE
324:                   DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
325:                END IF
326:             END IF
327: *
328:             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
329:                IF( NOTRAN ) THEN
330:                   IFUNC = IJOB
331:                END IF
332:                SCALE2 = SCALE
333:                CALL SLACPY( 'F', M, N, C, LDC, WORK, M )
334:                CALL SLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
335:                CALL SLASET( 'F', M, N, ZERO, ZERO, C, LDC )
336:                CALL SLASET( 'F', M, N, ZERO, ZERO, F, LDF )
337:             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
338:                CALL SLACPY( 'F', M, N, WORK, M, C, LDC )
339:                CALL SLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
340:                SCALE = SCALE2
341:             END IF
342:    30    CONTINUE
343: *
344:          RETURN
345:       END IF
346: *
347: *     Determine block structure of A
348: *
349:       P = 0
350:       I = 1
351:    40 CONTINUE
352:       IF( I.GT.M )
353:      $   GO TO 50
354:       P = P + 1
355:       IWORK( P ) = I
356:       I = I + MB
357:       IF( I.GE.M )
358:      $   GO TO 50
359:       IF( A( I, I-1 ).NE.ZERO )
360:      $   I = I + 1
361:       GO TO 40
362:    50 CONTINUE
363: *
364:       IWORK( P+1 ) = M + 1
365:       IF( IWORK( P ).EQ.IWORK( P+1 ) )
366:      $   P = P - 1
367: *
368: *     Determine block structure of B
369: *
370:       Q = P + 1
371:       J = 1
372:    60 CONTINUE
373:       IF( J.GT.N )
374:      $   GO TO 70
375:       Q = Q + 1
376:       IWORK( Q ) = J
377:       J = J + NB
378:       IF( J.GE.N )
379:      $   GO TO 70
380:       IF( B( J, J-1 ).NE.ZERO )
381:      $   J = J + 1
382:       GO TO 60
383:    70 CONTINUE
384: *
385:       IWORK( Q+1 ) = N + 1
386:       IF( IWORK( Q ).EQ.IWORK( Q+1 ) )
387:      $   Q = Q - 1
388: *
389:       IF( NOTRAN ) THEN
390: *
391:          DO 150 IROUND = 1, ISOLVE
392: *
393: *           Solve (I, J)-subsystem
394: *               A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
395: *               D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
396: *           for I = P, P - 1,..., 1; J = 1, 2,..., Q
397: *
398:             DSCALE = ZERO
399:             DSUM = ONE
400:             PQ = 0
401:             SCALE = ONE
402:             DO 130 J = P + 2, Q
403:                JS = IWORK( J )
404:                JE = IWORK( J+1 ) - 1
405:                NB = JE - JS + 1
406:                DO 120 I = P, 1, -1
407:                   IS = IWORK( I )
408:                   IE = IWORK( I+1 ) - 1
409:                   MB = IE - IS + 1
410:                   PPQQ = 0
411:                   CALL STGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
412:      $                         B( JS, JS ), LDB, C( IS, JS ), LDC,
413:      $                         D( IS, IS ), LDD, E( JS, JS ), LDE,
414:      $                         F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
415:      $                         IWORK( Q+2 ), PPQQ, LINFO )
416:                   IF( LINFO.GT.0 )
417:      $               INFO = LINFO
418: *
419:                   PQ = PQ + PPQQ
420:                   IF( SCALOC.NE.ONE ) THEN
421:                      DO 80 K = 1, JS - 1
422:                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
423:                         CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
424:    80                CONTINUE
425:                      DO 90 K = JS, JE
426:                         CALL SSCAL( IS-1, SCALOC, C( 1, K ), 1 )
427:                         CALL SSCAL( IS-1, SCALOC, F( 1, K ), 1 )
428:    90                CONTINUE
429:                      DO 100 K = JS, JE
430:                         CALL SSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
431:                         CALL SSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
432:   100                CONTINUE
433:                      DO 110 K = JE + 1, N
434:                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
435:                         CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
436:   110                CONTINUE
437:                      SCALE = SCALE*SCALOC
438:                   END IF
439: *
440: *                 Substitute R(I, J) and L(I, J) into remaining
441: *                 equation.
442: *
443:                   IF( I.GT.1 ) THEN
444:                      CALL SGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
445:      $                           A( 1, IS ), LDA, C( IS, JS ), LDC, ONE,
446:      $                           C( 1, JS ), LDC )
447:                      CALL SGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
448:      $                           D( 1, IS ), LDD, C( IS, JS ), LDC, ONE,
449:      $                           F( 1, JS ), LDF )
450:                   END IF
451:                   IF( J.LT.Q ) THEN
452:                      CALL SGEMM( 'N', 'N', MB, N-JE, NB, ONE,
453:      $                           F( IS, JS ), LDF, B( JS, JE+1 ), LDB,
454:      $                           ONE, C( IS, JE+1 ), LDC )
455:                      CALL SGEMM( 'N', 'N', MB, N-JE, NB, ONE,
456:      $                           F( IS, JS ), LDF, E( JS, JE+1 ), LDE,
457:      $                           ONE, F( IS, JE+1 ), LDF )
458:                   END IF
459:   120          CONTINUE
460:   130       CONTINUE
461:             IF( DSCALE.NE.ZERO ) THEN
462:                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
463:                   DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
464:                ELSE
465:                   DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
466:                END IF
467:             END IF
468:             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
469:                IF( NOTRAN ) THEN
470:                   IFUNC = IJOB
471:                END IF
472:                SCALE2 = SCALE
473:                CALL SLACPY( 'F', M, N, C, LDC, WORK, M )
474:                CALL SLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
475:                CALL SLASET( 'F', M, N, ZERO, ZERO, C, LDC )
476:                CALL SLASET( 'F', M, N, ZERO, ZERO, F, LDF )
477:             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
478:                CALL SLACPY( 'F', M, N, WORK, M, C, LDC )
479:                CALL SLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
480:                SCALE = SCALE2
481:             END IF
482:   150    CONTINUE
483: *
484:       ELSE
485: *
486: *        Solve transposed (I, J)-subsystem
487: *             A(I, I)' * R(I, J)  + D(I, I)' * L(I, J)  =  C(I, J)
488: *             R(I, J)  * B(J, J)' + L(I, J)  * E(J, J)' = -F(I, J)
489: *        for I = 1,2,..., P; J = Q, Q-1,..., 1
490: *
491:          SCALE = ONE
492:          DO 210 I = 1, P
493:             IS = IWORK( I )
494:             IE = IWORK( I+1 ) - 1
495:             MB = IE - IS + 1
496:             DO 200 J = Q, P + 2, -1
497:                JS = IWORK( J )
498:                JE = IWORK( J+1 ) - 1
499:                NB = JE - JS + 1
500:                CALL STGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
501:      $                      B( JS, JS ), LDB, C( IS, JS ), LDC,
502:      $                      D( IS, IS ), LDD, E( JS, JS ), LDE,
503:      $                      F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
504:      $                      IWORK( Q+2 ), PPQQ, LINFO )
505:                IF( LINFO.GT.0 )
506:      $            INFO = LINFO
507:                IF( SCALOC.NE.ONE ) THEN
508:                   DO 160 K = 1, JS - 1
509:                      CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
510:                      CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
511:   160             CONTINUE
512:                   DO 170 K = JS, JE
513:                      CALL SSCAL( IS-1, SCALOC, C( 1, K ), 1 )
514:                      CALL SSCAL( IS-1, SCALOC, F( 1, K ), 1 )
515:   170             CONTINUE
516:                   DO 180 K = JS, JE
517:                      CALL SSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
518:                      CALL SSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
519:   180             CONTINUE
520:                   DO 190 K = JE + 1, N
521:                      CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
522:                      CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
523:   190             CONTINUE
524:                   SCALE = SCALE*SCALOC
525:                END IF
526: *
527: *              Substitute R(I, J) and L(I, J) into remaining equation.
528: *
529:                IF( J.GT.P+2 ) THEN
530:                   CALL SGEMM( 'N', 'T', MB, JS-1, NB, ONE, C( IS, JS ),
531:      $                        LDC, B( 1, JS ), LDB, ONE, F( IS, 1 ),
532:      $                        LDF )
533:                   CALL SGEMM( 'N', 'T', MB, JS-1, NB, ONE, F( IS, JS ),
534:      $                        LDF, E( 1, JS ), LDE, ONE, F( IS, 1 ),
535:      $                        LDF )
536:                END IF
537:                IF( I.LT.P ) THEN
538:                   CALL SGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
539:      $                        A( IS, IE+1 ), LDA, C( IS, JS ), LDC, ONE,
540:      $                        C( IE+1, JS ), LDC )
541:                   CALL SGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
542:      $                        D( IS, IE+1 ), LDD, F( IS, JS ), LDF, ONE,
543:      $                        C( IE+1, JS ), LDC )
544:                END IF
545:   200       CONTINUE
546:   210    CONTINUE
547: *
548:       END IF
549: *
550:       WORK( 1 ) = LWMIN
551: *
552:       RETURN
553: *
554: *     End of STGSYL
555: *
556:       END
557: