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