001:       SUBROUTINE ZTGSYL( 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: *     January 2007
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          TRANS
011:       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
012:      $                   LWORK, M, N
013:       DOUBLE PRECISION   DIF, SCALE
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            IWORK( * )
017:       COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * ),
018:      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
019:      $                   WORK( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  ZTGSYL 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 complex entries. A, B, D and E are upper
033: *  triangular (i.e., (A,D) and (B,E) in generalized Schur form).
034: *
035: *  The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1
036: *  is an output scaling factor chosen to avoid overflow.
037: *
038: *  In matrix notation (1) is equivalent to solve Zx = scale*b, where Z
039: *  is defined as
040: *
041: *         Z = [ kron(In, A)  -kron(B', Im) ]        (2)
042: *             [ kron(In, D)  -kron(E', Im) ],
043: *
044: *  Here Ix is the identity matrix of size x and X' is the conjugate
045: *  transpose of X. Kron(X, Y) is the Kronecker product between the
046: *  matrices X and Y.
047: *
048: *  If TRANS = 'C', y in the conjugate transposed system Z'*y = scale*b
049: *  is solved for, 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 = 'C') 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 ZLACON.
057: *
058: *  If IJOB >= 1, ZTGSYL computes a Frobenius norm-based estimate of
059: *  Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
060: *  reciprocal of the smallest singular value of Z.
061: *
062: *  This is a level-3 BLAS algorithm.
063: *
064: *  Arguments
065: *  =========
066: *
067: *  TRANS   (input) CHARACTER*1
068: *          = 'N': solve the generalized sylvester equation (1).
069: *          = 'C': solve the "conjugate transposed" system (3).
070: *
071: *  IJOB    (input) INTEGER
072: *          Specifies what kind of functionality to be performed.
073: *          =0: solve (1) only.
074: *          =1: The functionality of 0 and 3.
075: *          =2: The functionality of 0 and 4.
076: *          =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
077: *              (look ahead strategy is used).
078: *          =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
079: *              (ZGECON on sub-systems is used).
080: *          Not referenced if TRANS = 'C'.
081: *
082: *  M       (input) INTEGER
083: *          The order of the matrices A and D, and the row dimension of
084: *          the matrices C, F, R and L.
085: *
086: *  N       (input) INTEGER
087: *          The order of the matrices B and E, and the column dimension
088: *          of the matrices C, F, R and L.
089: *
090: *  A       (input) COMPLEX*16 array, dimension (LDA, M)
091: *          The upper triangular matrix A.
092: *
093: *  LDA     (input) INTEGER
094: *          The leading dimension of the array A. LDA >= max(1, M).
095: *
096: *  B       (input) COMPLEX*16 array, dimension (LDB, N)
097: *          The upper triangular matrix B.
098: *
099: *  LDB     (input) INTEGER
100: *          The leading dimension of the array B. LDB >= max(1, N).
101: *
102: *  C       (input/output) COMPLEX*16 array, dimension (LDC, N)
103: *          On entry, C contains the right-hand-side of the first matrix
104: *          equation in (1) or (3).
105: *          On exit, if IJOB = 0, 1 or 2, C has been overwritten by
106: *          the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R,
107: *          the solution achieved during the computation of the
108: *          Dif-estimate.
109: *
110: *  LDC     (input) INTEGER
111: *          The leading dimension of the array C. LDC >= max(1, M).
112: *
113: *  D       (input) COMPLEX*16 array, dimension (LDD, M)
114: *          The upper triangular matrix D.
115: *
116: *  LDD     (input) INTEGER
117: *          The leading dimension of the array D. LDD >= max(1, M).
118: *
119: *  E       (input) COMPLEX*16 array, dimension (LDE, N)
120: *          The upper triangular matrix E.
121: *
122: *  LDE     (input) INTEGER
123: *          The leading dimension of the array E. LDE >= max(1, N).
124: *
125: *  F       (input/output) COMPLEX*16 array, dimension (LDF, N)
126: *          On entry, F contains the right-hand-side of the second matrix
127: *          equation in (1) or (3).
128: *          On exit, if IJOB = 0, 1 or 2, F has been overwritten by
129: *          the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L,
130: *          the solution achieved during the computation of the
131: *          Dif-estimate.
132: *
133: *  LDF     (input) INTEGER
134: *          The leading dimension of the array F. LDF >= max(1, M).
135: *
136: *  DIF     (output) DOUBLE PRECISION
137: *          On exit DIF is the reciprocal of a lower bound of the
138: *          reciprocal of the Dif-function, i.e. DIF is an upper bound of
139: *          Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2).
140: *          IF IJOB = 0 or TRANS = 'C', DIF is not referenced.
141: *
142: *  SCALE   (output) DOUBLE PRECISION
143: *          On exit SCALE is the scaling factor in (1) or (3).
144: *          If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
145: *          to a slightly perturbed system but the input matrices A, B,
146: *          D and E have not been changed. If SCALE = 0, R and L will
147: *          hold the solutions to the homogenious system with C = F = 0.
148: *
149: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
150: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
151: *
152: *  LWORK   (input) INTEGER
153: *          The dimension of the array WORK. LWORK > = 1.
154: *          If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
155: *
156: *          If LWORK = -1, then a workspace query is assumed; the routine
157: *          only calculates the optimal size of the WORK array, returns
158: *          this value as the first entry of the WORK array, and no error
159: *          message related to LWORK is issued by XERBLA.
160: *
161: *  IWORK   (workspace) INTEGER array, dimension (M+N+2)
162: *
163: *  INFO    (output) INTEGER
164: *            =0: successful exit
165: *            <0: If INFO = -i, the i-th argument had an illegal value.
166: *            >0: (A, D) and (B, E) have common or very close
167: *                eigenvalues.
168: *
169: *  Further Details
170: *  ===============
171: *
172: *  Based on contributions by
173: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
174: *     Umea University, S-901 87 Umea, Sweden.
175: *
176: *  [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
177: *      for Solving the Generalized Sylvester Equation and Estimating the
178: *      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
179: *      Department of Computing Science, Umea University, S-901 87 Umea,
180: *      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
181: *      Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
182: *      No 1, 1996.
183: *
184: *  [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
185: *      Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
186: *      Appl., 15(4):1045-1060, 1994.
187: *
188: *  [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
189: *      Condition Estimators for Solving the Generalized Sylvester
190: *      Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
191: *      July 1989, pp 745-751.
192: *
193: *  =====================================================================
194: *  Replaced various illegal calls to CCOPY by calls to CLASET.
195: *  Sven Hammarling, 1/5/02.
196: *
197: *     .. Parameters ..
198:       DOUBLE PRECISION   ZERO, ONE
199:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
200:       COMPLEX*16         CZERO
201:       PARAMETER          ( CZERO = (0.0D+0, 0.0D+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, PQ, Q
207:       DOUBLE PRECISION   DSCALE, DSUM, SCALE2, SCALOC
208: *     ..
209: *     .. External Functions ..
210:       LOGICAL            LSAME
211:       INTEGER            ILAENV
212:       EXTERNAL           LSAME, ILAENV
213: *     ..
214: *     .. External Subroutines ..
215:       EXTERNAL           XERBLA, ZGEMM, ZLACPY, ZLASET, ZSCAL, ZTGSY2
216: *     ..
217: *     .. Intrinsic Functions ..
218:       INTRINSIC          DBLE, DCMPLX, MAX, 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, 'C' ) ) 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( 'ZTGSYL', -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, 'ZTGSYL', TRANS, M, N, -1, -1 )
294:       NB = ILAENV( 5, 'ZTGSYL', 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 ZLASET( 'F', M, N, CZERO, CZERO, C, LDC )
302:             CALL ZLASET( 'F', M, N, CZERO, CZERO, 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: *        Use unblocked Level 2 solver
312: *
313:          DO 30 IROUND = 1, ISOLVE
314: *
315:             SCALE = ONE
316:             DSCALE = ZERO
317:             DSUM = ONE
318:             PQ = M*N
319:             CALL ZTGSY2( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D,
320:      $                   LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE,
321:      $                   INFO )
322:             IF( DSCALE.NE.ZERO ) THEN
323:                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
324:                   DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
325:                ELSE
326:                   DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
327:                END IF
328:             END IF
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 ZLACPY( 'F', M, N, C, LDC, WORK, M )
335:                CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
336:                CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC )
337:                CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF )
338:             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
339:                CALL ZLACPY( 'F', M, N, WORK, M, C, LDC )
340:                CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
341:                SCALE = SCALE2
342:             END IF
343:    30    CONTINUE
344: *
345:          RETURN
346: *
347:       END IF
348: *
349: *     Determine block structure of A
350: *
351:       P = 0
352:       I = 1
353:    40 CONTINUE
354:       IF( I.GT.M )
355:      $   GO TO 50
356:       P = P + 1
357:       IWORK( P ) = I
358:       I = I + MB
359:       IF( I.GE.M )
360:      $   GO TO 50
361:       GO TO 40
362:    50 CONTINUE
363:       IWORK( P+1 ) = M + 1
364:       IF( IWORK( P ).EQ.IWORK( P+1 ) )
365:      $   P = P - 1
366: *
367: *     Determine block structure of B
368: *
369:       Q = P + 1
370:       J = 1
371:    60 CONTINUE
372:       IF( J.GT.N )
373:      $   GO TO 70
374: *
375:       Q = Q + 1
376:       IWORK( Q ) = J
377:       J = J + NB
378:       IF( J.GE.N )
379:      $   GO TO 70
380:       GO TO 60
381: *
382:    70 CONTINUE
383:       IWORK( Q+1 ) = N + 1
384:       IF( IWORK( Q ).EQ.IWORK( Q+1 ) )
385:      $   Q = Q - 1
386: *
387:       IF( NOTRAN ) THEN
388:          DO 150 IROUND = 1, ISOLVE
389: *
390: *           Solve (I, J) - subsystem
391: *               A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
392: *               D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
393: *           for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
394: *
395:             PQ = 0
396:             SCALE = ONE
397:             DSCALE = ZERO
398:             DSUM = ONE
399:             DO 130 J = P + 2, Q
400:                JS = IWORK( J )
401:                JE = IWORK( J+1 ) - 1
402:                NB = JE - JS + 1
403:                DO 120 I = P, 1, -1
404:                   IS = IWORK( I )
405:                   IE = IWORK( I+1 ) - 1
406:                   MB = IE - IS + 1
407:                   CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
408:      $                         B( JS, JS ), LDB, C( IS, JS ), LDC,
409:      $                         D( IS, IS ), LDD, E( JS, JS ), LDE,
410:      $                         F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
411:      $                         LINFO )
412:                   IF( LINFO.GT.0 )
413:      $               INFO = LINFO
414:                   PQ = PQ + MB*NB
415:                   IF( SCALOC.NE.ONE ) THEN
416:                      DO 80 K = 1, JS - 1
417:                         CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
418:      $                              C( 1, K ), 1 )
419:                         CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
420:      $                              F( 1, K ), 1 )
421:    80                CONTINUE
422:                      DO 90 K = JS, JE
423:                         CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
424:      $                              C( 1, K ), 1 )
425:                         CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
426:      $                              F( 1, K ), 1 )
427:    90                CONTINUE
428:                      DO 100 K = JS, JE
429:                         CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
430:      $                              C( IE+1, K ), 1 )
431:                         CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
432:      $                              F( IE+1, K ), 1 )
433:   100                CONTINUE
434:                      DO 110 K = JE + 1, N
435:                         CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
436:      $                              C( 1, K ), 1 )
437:                         CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
438:      $                              F( 1, K ), 1 )
439:   110                CONTINUE
440:                      SCALE = SCALE*SCALOC
441:                   END IF
442: *
443: *                 Substitute R(I,J) and L(I,J) into remaining equation.
444: *
445:                   IF( I.GT.1 ) THEN
446:                      CALL ZGEMM( 'N', 'N', IS-1, NB, MB,
447:      $                           DCMPLX( -ONE, ZERO ), A( 1, IS ), LDA,
448:      $                           C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
449:      $                           C( 1, JS ), LDC )
450:                      CALL ZGEMM( 'N', 'N', IS-1, NB, MB,
451:      $                           DCMPLX( -ONE, ZERO ), D( 1, IS ), LDD,
452:      $                           C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
453:      $                           F( 1, JS ), LDF )
454:                   END IF
455:                   IF( J.LT.Q ) THEN
456:                      CALL ZGEMM( 'N', 'N', MB, N-JE, NB,
457:      $                           DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
458:      $                           B( JS, JE+1 ), LDB,
459:      $                           DCMPLX( ONE, ZERO ), C( IS, JE+1 ),
460:      $                           LDC )
461:                      CALL ZGEMM( 'N', 'N', MB, N-JE, NB,
462:      $                           DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
463:      $                           E( JS, JE+1 ), LDE,
464:      $                           DCMPLX( ONE, ZERO ), F( IS, JE+1 ),
465:      $                           LDF )
466:                   END IF
467:   120          CONTINUE
468:   130       CONTINUE
469:             IF( DSCALE.NE.ZERO ) THEN
470:                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
471:                   DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
472:                ELSE
473:                   DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
474:                END IF
475:             END IF
476:             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
477:                IF( NOTRAN ) THEN
478:                   IFUNC = IJOB
479:                END IF
480:                SCALE2 = SCALE
481:                CALL ZLACPY( 'F', M, N, C, LDC, WORK, M )
482:                CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
483:                CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC )
484:                CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF )
485:             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
486:                CALL ZLACPY( 'F', M, N, WORK, M, C, LDC )
487:                CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
488:                SCALE = SCALE2
489:             END IF
490:   150    CONTINUE
491:       ELSE
492: *
493: *        Solve transposed (I, J)-subsystem
494: *            A(I, I)' * R(I, J) + D(I, I)' * L(I, J) = C(I, J)
495: *            R(I, J) * B(J, J)  + L(I, J) * E(J, J) = -F(I, J)
496: *        for I = 1,2,..., P; J = Q, Q-1,..., 1
497: *
498:          SCALE = ONE
499:          DO 210 I = 1, P
500:             IS = IWORK( I )
501:             IE = IWORK( I+1 ) - 1
502:             MB = IE - IS + 1
503:             DO 200 J = Q, P + 2, -1
504:                JS = IWORK( J )
505:                JE = IWORK( J+1 ) - 1
506:                NB = JE - JS + 1
507:                CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
508:      $                      B( JS, JS ), LDB, C( IS, JS ), LDC,
509:      $                      D( IS, IS ), LDD, E( JS, JS ), LDE,
510:      $                      F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
511:      $                      LINFO )
512:                IF( LINFO.GT.0 )
513:      $            INFO = LINFO
514:                IF( SCALOC.NE.ONE ) THEN
515:                   DO 160 K = 1, JS - 1
516:                      CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
517:      $                           1 )
518:                      CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
519:      $                           1 )
520:   160             CONTINUE
521:                   DO 170 K = JS, JE
522:                      CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
523:      $                           C( 1, K ), 1 )
524:                      CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
525:      $                           F( 1, K ), 1 )
526:   170             CONTINUE
527:                   DO 180 K = JS, JE
528:                      CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
529:      $                           C( IE+1, K ), 1 )
530:                      CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
531:      $                           F( IE+1, K ), 1 )
532:   180             CONTINUE
533:                   DO 190 K = JE + 1, N
534:                      CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
535:      $                           1 )
536:                      CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
537:      $                           1 )
538:   190             CONTINUE
539:                   SCALE = SCALE*SCALOC
540:                END IF
541: *
542: *              Substitute R(I,J) and L(I,J) into remaining equation.
543: *
544:                IF( J.GT.P+2 ) THEN
545:                   CALL ZGEMM( 'N', 'C', MB, JS-1, NB,
546:      $                        DCMPLX( ONE, ZERO ), C( IS, JS ), LDC,
547:      $                        B( 1, JS ), LDB, DCMPLX( ONE, ZERO ),
548:      $                        F( IS, 1 ), LDF )
549:                   CALL ZGEMM( 'N', 'C', MB, JS-1, NB,
550:      $                        DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
551:      $                        E( 1, JS ), LDE, DCMPLX( ONE, ZERO ),
552:      $                        F( IS, 1 ), LDF )
553:                END IF
554:                IF( I.LT.P ) THEN
555:                   CALL ZGEMM( 'C', 'N', M-IE, NB, MB,
556:      $                        DCMPLX( -ONE, ZERO ), A( IS, IE+1 ), LDA,
557:      $                        C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
558:      $                        C( IE+1, JS ), LDC )
559:                   CALL ZGEMM( 'C', 'N', M-IE, NB, MB,
560:      $                        DCMPLX( -ONE, ZERO ), D( IS, IE+1 ), LDD,
561:      $                        F( IS, JS ), LDF, DCMPLX( ONE, ZERO ),
562:      $                        C( IE+1, JS ), LDC )
563:                END IF
564:   200       CONTINUE
565:   210    CONTINUE
566:       END IF
567: *
568:       WORK( 1 ) = LWMIN
569: *
570:       RETURN
571: *
572: *     End of ZTGSYL
573: *
574:       END
575: