1 SUBROUTINE ZTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
2 $ LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
3 $ IWORK, INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER TRANS
12 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
13 $ LWORK, M, N
14 DOUBLE PRECISION DIF, SCALE
15 * ..
16 * .. Array Arguments ..
17 INTEGER IWORK( * )
18 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
19 $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
20 $ WORK( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * ZTGSYL solves the generalized Sylvester equation:
27 *
28 * A * R - L * B = scale * C (1)
29 * D * R - L * E = scale * F
30 *
31 * where R and L are unknown m-by-n matrices, (A, D), (B, E) and
32 * (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
33 * respectively, with complex entries. A, B, D and E are upper
34 * triangular (i.e., (A,D) and (B,E) in generalized Schur form).
35 *
36 * The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1
37 * is an output scaling factor chosen to avoid overflow.
38 *
39 * In matrix notation (1) is equivalent to solve Zx = scale*b, where Z
40 * is defined as
41 *
42 * Z = [ kron(In, A) -kron(B**H, Im) ] (2)
43 * [ kron(In, D) -kron(E**H, Im) ],
44 *
45 * Here Ix is the identity matrix of size x and X**H is the conjugate
46 * transpose of X. Kron(X, Y) is the Kronecker product between the
47 * matrices X and Y.
48 *
49 * If TRANS = 'C', y in the conjugate transposed system Z**H *y = scale*b
50 * is solved for, which is equivalent to solve for R and L in
51 *
52 * A**H * R + D**H * L = scale * C (3)
53 * R * B**H + L * E**H = scale * -F
54 *
55 * This case (TRANS = 'C') is used to compute an one-norm-based estimate
56 * of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
57 * and (B,E), using ZLACON.
58 *
59 * If IJOB >= 1, ZTGSYL computes a Frobenius norm-based estimate of
60 * Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
61 * reciprocal of the smallest singular value of Z.
62 *
63 * This is a level-3 BLAS algorithm.
64 *
65 * Arguments
66 * =========
67 *
68 * TRANS (input) CHARACTER*1
69 * = 'N': solve the generalized sylvester equation (1).
70 * = 'C': solve the "conjugate transposed" system (3).
71 *
72 * IJOB (input) INTEGER
73 * Specifies what kind of functionality to be performed.
74 * =0: solve (1) only.
75 * =1: The functionality of 0 and 3.
76 * =2: The functionality of 0 and 4.
77 * =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
78 * (look ahead strategy is used).
79 * =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
80 * (ZGECON on sub-systems is used).
81 * Not referenced if TRANS = 'C'.
82 *
83 * M (input) INTEGER
84 * The order of the matrices A and D, and the row dimension of
85 * the matrices C, F, R and L.
86 *
87 * N (input) INTEGER
88 * The order of the matrices B and E, and the column dimension
89 * of the matrices C, F, R and L.
90 *
91 * A (input) COMPLEX*16 array, dimension (LDA, M)
92 * The upper triangular matrix A.
93 *
94 * LDA (input) INTEGER
95 * The leading dimension of the array A. LDA >= max(1, M).
96 *
97 * B (input) COMPLEX*16 array, dimension (LDB, N)
98 * The upper triangular matrix B.
99 *
100 * LDB (input) INTEGER
101 * The leading dimension of the array B. LDB >= max(1, N).
102 *
103 * C (input/output) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) DOUBLE PRECISION
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 = 'C', DIF is not referenced.
142 *
143 * SCALE (output) DOUBLE PRECISION
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,
147 * D and E have not been changed. If SCALE = 0, R and L will
148 * hold the solutions to the homogenious system with C = F = 0.
149 *
150 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
151 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
152 *
153 * LWORK (input) INTEGER
154 * The dimension of the array WORK. LWORK > = 1.
155 * If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
156 *
157 * If LWORK = -1, then a workspace query is assumed; the routine
158 * only calculates the optimal size of the WORK array, returns
159 * this value as the first entry of the WORK array, and no error
160 * message related to LWORK is issued by XERBLA.
161 *
162 * IWORK (workspace) INTEGER array, dimension (M+N+2)
163 *
164 * INFO (output) INTEGER
165 * =0: successful exit
166 * <0: If INFO = -i, the i-th argument had an illegal value.
167 * >0: (A, D) and (B, E) have common or very close
168 * 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 CCOPY by calls to CLASET.
196 * Sven Hammarling, 1/5/02.
197 *
198 * .. Parameters ..
199 DOUBLE PRECISION ZERO, ONE
200 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
201 COMPLEX*16 CZERO
202 PARAMETER ( CZERO = (0.0D+0, 0.0D+0) )
203 * ..
204 * .. Local Scalars ..
205 LOGICAL LQUERY, NOTRAN
206 INTEGER I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
207 $ LINFO, LWMIN, MB, NB, P, PQ, Q
208 DOUBLE PRECISION DSCALE, DSUM, SCALE2, SCALOC
209 * ..
210 * .. External Functions ..
211 LOGICAL LSAME
212 INTEGER ILAENV
213 EXTERNAL LSAME, ILAENV
214 * ..
215 * .. External Subroutines ..
216 EXTERNAL XERBLA, ZGEMM, ZLACPY, ZLASET, ZSCAL, ZTGSY2
217 * ..
218 * .. Intrinsic Functions ..
219 INTRINSIC DBLE, DCMPLX, MAX, SQRT
220 * ..
221 * .. Executable Statements ..
222 *
223 * Decode and test input parameters
224 *
225 INFO = 0
226 NOTRAN = LSAME( TRANS, 'N' )
227 LQUERY = ( LWORK.EQ.-1 )
228 *
229 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
230 INFO = -1
231 ELSE IF( NOTRAN ) THEN
232 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN
233 INFO = -2
234 END IF
235 END IF
236 IF( INFO.EQ.0 ) THEN
237 IF( M.LE.0 ) THEN
238 INFO = -3
239 ELSE IF( N.LE.0 ) THEN
240 INFO = -4
241 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
242 INFO = -6
243 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
244 INFO = -8
245 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
246 INFO = -10
247 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
248 INFO = -12
249 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
250 INFO = -14
251 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
252 INFO = -16
253 END IF
254 END IF
255 *
256 IF( INFO.EQ.0 ) THEN
257 IF( NOTRAN ) THEN
258 IF( IJOB.EQ.1 .OR. IJOB.EQ.2 ) THEN
259 LWMIN = MAX( 1, 2*M*N )
260 ELSE
261 LWMIN = 1
262 END IF
263 ELSE
264 LWMIN = 1
265 END IF
266 WORK( 1 ) = LWMIN
267 *
268 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
269 INFO = -20
270 END IF
271 END IF
272 *
273 IF( INFO.NE.0 ) THEN
274 CALL XERBLA( 'ZTGSYL', -INFO )
275 RETURN
276 ELSE IF( LQUERY ) THEN
277 RETURN
278 END IF
279 *
280 * Quick return if possible
281 *
282 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
283 SCALE = 1
284 IF( NOTRAN ) THEN
285 IF( IJOB.NE.0 ) THEN
286 DIF = 0
287 END IF
288 END IF
289 RETURN
290 END IF
291 *
292 * Determine optimal block sizes MB and NB
293 *
294 MB = ILAENV( 2, 'ZTGSYL', TRANS, M, N, -1, -1 )
295 NB = ILAENV( 5, 'ZTGSYL', TRANS, M, N, -1, -1 )
296 *
297 ISOLVE = 1
298 IFUNC = 0
299 IF( NOTRAN ) THEN
300 IF( IJOB.GE.3 ) THEN
301 IFUNC = IJOB - 2
302 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC )
303 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF )
304 ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN
305 ISOLVE = 2
306 END IF
307 END IF
308 *
309 IF( ( MB.LE.1 .AND. NB.LE.1 ) .OR. ( MB.GE.M .AND. NB.GE.N ) )
310 $ THEN
311 *
312 * Use unblocked Level 2 solver
313 *
314 DO 30 IROUND = 1, ISOLVE
315 *
316 SCALE = ONE
317 DSCALE = ZERO
318 DSUM = ONE
319 PQ = M*N
320 CALL ZTGSY2( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D,
321 $ LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE,
322 $ INFO )
323 IF( DSCALE.NE.ZERO ) THEN
324 IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
325 DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
326 ELSE
327 DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
328 END IF
329 END IF
330 IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
331 IF( NOTRAN ) THEN
332 IFUNC = IJOB
333 END IF
334 SCALE2 = SCALE
335 CALL ZLACPY( 'F', M, N, C, LDC, WORK, M )
336 CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
337 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC )
338 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF )
339 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
340 CALL ZLACPY( 'F', M, N, WORK, M, C, LDC )
341 CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
342 SCALE = SCALE2
343 END IF
344 30 CONTINUE
345 *
346 RETURN
347 *
348 END IF
349 *
350 * Determine block structure of A
351 *
352 P = 0
353 I = 1
354 40 CONTINUE
355 IF( I.GT.M )
356 $ GO TO 50
357 P = P + 1
358 IWORK( P ) = I
359 I = I + MB
360 IF( I.GE.M )
361 $ GO TO 50
362 GO TO 40
363 50 CONTINUE
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 *
376 Q = Q + 1
377 IWORK( Q ) = J
378 J = J + NB
379 IF( J.GE.N )
380 $ GO TO 70
381 GO TO 60
382 *
383 70 CONTINUE
384 IWORK( Q+1 ) = N + 1
385 IF( IWORK( Q ).EQ.IWORK( Q+1 ) )
386 $ Q = Q - 1
387 *
388 IF( NOTRAN ) THEN
389 DO 150 IROUND = 1, ISOLVE
390 *
391 * Solve (I, J) - subsystem
392 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
393 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
394 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
395 *
396 PQ = 0
397 SCALE = ONE
398 DSCALE = ZERO
399 DSUM = ONE
400 DO 130 J = P + 2, Q
401 JS = IWORK( J )
402 JE = IWORK( J+1 ) - 1
403 NB = JE - JS + 1
404 DO 120 I = P, 1, -1
405 IS = IWORK( I )
406 IE = IWORK( I+1 ) - 1
407 MB = IE - IS + 1
408 CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
409 $ B( JS, JS ), LDB, C( IS, JS ), LDC,
410 $ D( IS, IS ), LDD, E( JS, JS ), LDE,
411 $ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
412 $ LINFO )
413 IF( LINFO.GT.0 )
414 $ INFO = LINFO
415 PQ = PQ + MB*NB
416 IF( SCALOC.NE.ONE ) THEN
417 DO 80 K = 1, JS - 1
418 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
419 $ C( 1, K ), 1 )
420 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
421 $ F( 1, K ), 1 )
422 80 CONTINUE
423 DO 90 K = JS, JE
424 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
425 $ C( 1, K ), 1 )
426 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
427 $ F( 1, K ), 1 )
428 90 CONTINUE
429 DO 100 K = JS, JE
430 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
431 $ C( IE+1, K ), 1 )
432 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
433 $ F( IE+1, K ), 1 )
434 100 CONTINUE
435 DO 110 K = JE + 1, N
436 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
437 $ C( 1, K ), 1 )
438 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
439 $ F( 1, K ), 1 )
440 110 CONTINUE
441 SCALE = SCALE*SCALOC
442 END IF
443 *
444 * Substitute R(I,J) and L(I,J) into remaining equation.
445 *
446 IF( I.GT.1 ) THEN
447 CALL ZGEMM( 'N', 'N', IS-1, NB, MB,
448 $ DCMPLX( -ONE, ZERO ), A( 1, IS ), LDA,
449 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
450 $ C( 1, JS ), LDC )
451 CALL ZGEMM( 'N', 'N', IS-1, NB, MB,
452 $ DCMPLX( -ONE, ZERO ), D( 1, IS ), LDD,
453 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
454 $ F( 1, JS ), LDF )
455 END IF
456 IF( J.LT.Q ) THEN
457 CALL ZGEMM( 'N', 'N', MB, N-JE, NB,
458 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
459 $ B( JS, JE+1 ), LDB,
460 $ DCMPLX( ONE, ZERO ), C( IS, JE+1 ),
461 $ LDC )
462 CALL ZGEMM( 'N', 'N', MB, N-JE, NB,
463 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
464 $ E( JS, JE+1 ), LDE,
465 $ DCMPLX( ONE, ZERO ), F( IS, JE+1 ),
466 $ LDF )
467 END IF
468 120 CONTINUE
469 130 CONTINUE
470 IF( DSCALE.NE.ZERO ) THEN
471 IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
472 DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
473 ELSE
474 DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
475 END IF
476 END IF
477 IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
478 IF( NOTRAN ) THEN
479 IFUNC = IJOB
480 END IF
481 SCALE2 = SCALE
482 CALL ZLACPY( 'F', M, N, C, LDC, WORK, M )
483 CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
484 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC )
485 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF )
486 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
487 CALL ZLACPY( 'F', M, N, WORK, M, C, LDC )
488 CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
489 SCALE = SCALE2
490 END IF
491 150 CONTINUE
492 ELSE
493 *
494 * Solve transposed (I, J)-subsystem
495 * A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
496 * R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
497 * for I = 1,2,..., P; J = Q, Q-1,..., 1
498 *
499 SCALE = ONE
500 DO 210 I = 1, P
501 IS = IWORK( I )
502 IE = IWORK( I+1 ) - 1
503 MB = IE - IS + 1
504 DO 200 J = Q, P + 2, -1
505 JS = IWORK( J )
506 JE = IWORK( J+1 ) - 1
507 NB = JE - JS + 1
508 CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
509 $ B( JS, JS ), LDB, C( IS, JS ), LDC,
510 $ D( IS, IS ), LDD, E( JS, JS ), LDE,
511 $ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
512 $ LINFO )
513 IF( LINFO.GT.0 )
514 $ INFO = LINFO
515 IF( SCALOC.NE.ONE ) THEN
516 DO 160 K = 1, JS - 1
517 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
518 $ 1 )
519 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
520 $ 1 )
521 160 CONTINUE
522 DO 170 K = JS, JE
523 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
524 $ C( 1, K ), 1 )
525 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
526 $ F( 1, K ), 1 )
527 170 CONTINUE
528 DO 180 K = JS, JE
529 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
530 $ C( IE+1, K ), 1 )
531 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
532 $ F( IE+1, K ), 1 )
533 180 CONTINUE
534 DO 190 K = JE + 1, N
535 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
536 $ 1 )
537 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
538 $ 1 )
539 190 CONTINUE
540 SCALE = SCALE*SCALOC
541 END IF
542 *
543 * Substitute R(I,J) and L(I,J) into remaining equation.
544 *
545 IF( J.GT.P+2 ) THEN
546 CALL ZGEMM( 'N', 'C', MB, JS-1, NB,
547 $ DCMPLX( ONE, ZERO ), C( IS, JS ), LDC,
548 $ B( 1, JS ), LDB, DCMPLX( ONE, ZERO ),
549 $ F( IS, 1 ), LDF )
550 CALL ZGEMM( 'N', 'C', MB, JS-1, NB,
551 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
552 $ E( 1, JS ), LDE, DCMPLX( ONE, ZERO ),
553 $ F( IS, 1 ), LDF )
554 END IF
555 IF( I.LT.P ) THEN
556 CALL ZGEMM( 'C', 'N', M-IE, NB, MB,
557 $ DCMPLX( -ONE, ZERO ), A( IS, IE+1 ), LDA,
558 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
559 $ C( IE+1, JS ), LDC )
560 CALL ZGEMM( 'C', 'N', M-IE, NB, MB,
561 $ DCMPLX( -ONE, ZERO ), D( IS, IE+1 ), LDD,
562 $ F( IS, JS ), LDF, DCMPLX( ONE, ZERO ),
563 $ C( IE+1, JS ), LDC )
564 END IF
565 200 CONTINUE
566 210 CONTINUE
567 END IF
568 *
569 WORK( 1 ) = LWMIN
570 *
571 RETURN
572 *
573 * End of ZTGSYL
574 *
575 END
2 $ LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
3 $ IWORK, INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER TRANS
12 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
13 $ LWORK, M, N
14 DOUBLE PRECISION DIF, SCALE
15 * ..
16 * .. Array Arguments ..
17 INTEGER IWORK( * )
18 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
19 $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
20 $ WORK( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * ZTGSYL solves the generalized Sylvester equation:
27 *
28 * A * R - L * B = scale * C (1)
29 * D * R - L * E = scale * F
30 *
31 * where R and L are unknown m-by-n matrices, (A, D), (B, E) and
32 * (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
33 * respectively, with complex entries. A, B, D and E are upper
34 * triangular (i.e., (A,D) and (B,E) in generalized Schur form).
35 *
36 * The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1
37 * is an output scaling factor chosen to avoid overflow.
38 *
39 * In matrix notation (1) is equivalent to solve Zx = scale*b, where Z
40 * is defined as
41 *
42 * Z = [ kron(In, A) -kron(B**H, Im) ] (2)
43 * [ kron(In, D) -kron(E**H, Im) ],
44 *
45 * Here Ix is the identity matrix of size x and X**H is the conjugate
46 * transpose of X. Kron(X, Y) is the Kronecker product between the
47 * matrices X and Y.
48 *
49 * If TRANS = 'C', y in the conjugate transposed system Z**H *y = scale*b
50 * is solved for, which is equivalent to solve for R and L in
51 *
52 * A**H * R + D**H * L = scale * C (3)
53 * R * B**H + L * E**H = scale * -F
54 *
55 * This case (TRANS = 'C') is used to compute an one-norm-based estimate
56 * of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
57 * and (B,E), using ZLACON.
58 *
59 * If IJOB >= 1, ZTGSYL computes a Frobenius norm-based estimate of
60 * Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
61 * reciprocal of the smallest singular value of Z.
62 *
63 * This is a level-3 BLAS algorithm.
64 *
65 * Arguments
66 * =========
67 *
68 * TRANS (input) CHARACTER*1
69 * = 'N': solve the generalized sylvester equation (1).
70 * = 'C': solve the "conjugate transposed" system (3).
71 *
72 * IJOB (input) INTEGER
73 * Specifies what kind of functionality to be performed.
74 * =0: solve (1) only.
75 * =1: The functionality of 0 and 3.
76 * =2: The functionality of 0 and 4.
77 * =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
78 * (look ahead strategy is used).
79 * =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
80 * (ZGECON on sub-systems is used).
81 * Not referenced if TRANS = 'C'.
82 *
83 * M (input) INTEGER
84 * The order of the matrices A and D, and the row dimension of
85 * the matrices C, F, R and L.
86 *
87 * N (input) INTEGER
88 * The order of the matrices B and E, and the column dimension
89 * of the matrices C, F, R and L.
90 *
91 * A (input) COMPLEX*16 array, dimension (LDA, M)
92 * The upper triangular matrix A.
93 *
94 * LDA (input) INTEGER
95 * The leading dimension of the array A. LDA >= max(1, M).
96 *
97 * B (input) COMPLEX*16 array, dimension (LDB, N)
98 * The upper triangular matrix B.
99 *
100 * LDB (input) INTEGER
101 * The leading dimension of the array B. LDB >= max(1, N).
102 *
103 * C (input/output) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) DOUBLE PRECISION
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 = 'C', DIF is not referenced.
142 *
143 * SCALE (output) DOUBLE PRECISION
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,
147 * D and E have not been changed. If SCALE = 0, R and L will
148 * hold the solutions to the homogenious system with C = F = 0.
149 *
150 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
151 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
152 *
153 * LWORK (input) INTEGER
154 * The dimension of the array WORK. LWORK > = 1.
155 * If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
156 *
157 * If LWORK = -1, then a workspace query is assumed; the routine
158 * only calculates the optimal size of the WORK array, returns
159 * this value as the first entry of the WORK array, and no error
160 * message related to LWORK is issued by XERBLA.
161 *
162 * IWORK (workspace) INTEGER array, dimension (M+N+2)
163 *
164 * INFO (output) INTEGER
165 * =0: successful exit
166 * <0: If INFO = -i, the i-th argument had an illegal value.
167 * >0: (A, D) and (B, E) have common or very close
168 * 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 CCOPY by calls to CLASET.
196 * Sven Hammarling, 1/5/02.
197 *
198 * .. Parameters ..
199 DOUBLE PRECISION ZERO, ONE
200 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
201 COMPLEX*16 CZERO
202 PARAMETER ( CZERO = (0.0D+0, 0.0D+0) )
203 * ..
204 * .. Local Scalars ..
205 LOGICAL LQUERY, NOTRAN
206 INTEGER I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
207 $ LINFO, LWMIN, MB, NB, P, PQ, Q
208 DOUBLE PRECISION DSCALE, DSUM, SCALE2, SCALOC
209 * ..
210 * .. External Functions ..
211 LOGICAL LSAME
212 INTEGER ILAENV
213 EXTERNAL LSAME, ILAENV
214 * ..
215 * .. External Subroutines ..
216 EXTERNAL XERBLA, ZGEMM, ZLACPY, ZLASET, ZSCAL, ZTGSY2
217 * ..
218 * .. Intrinsic Functions ..
219 INTRINSIC DBLE, DCMPLX, MAX, SQRT
220 * ..
221 * .. Executable Statements ..
222 *
223 * Decode and test input parameters
224 *
225 INFO = 0
226 NOTRAN = LSAME( TRANS, 'N' )
227 LQUERY = ( LWORK.EQ.-1 )
228 *
229 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
230 INFO = -1
231 ELSE IF( NOTRAN ) THEN
232 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN
233 INFO = -2
234 END IF
235 END IF
236 IF( INFO.EQ.0 ) THEN
237 IF( M.LE.0 ) THEN
238 INFO = -3
239 ELSE IF( N.LE.0 ) THEN
240 INFO = -4
241 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
242 INFO = -6
243 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
244 INFO = -8
245 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
246 INFO = -10
247 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
248 INFO = -12
249 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
250 INFO = -14
251 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
252 INFO = -16
253 END IF
254 END IF
255 *
256 IF( INFO.EQ.0 ) THEN
257 IF( NOTRAN ) THEN
258 IF( IJOB.EQ.1 .OR. IJOB.EQ.2 ) THEN
259 LWMIN = MAX( 1, 2*M*N )
260 ELSE
261 LWMIN = 1
262 END IF
263 ELSE
264 LWMIN = 1
265 END IF
266 WORK( 1 ) = LWMIN
267 *
268 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
269 INFO = -20
270 END IF
271 END IF
272 *
273 IF( INFO.NE.0 ) THEN
274 CALL XERBLA( 'ZTGSYL', -INFO )
275 RETURN
276 ELSE IF( LQUERY ) THEN
277 RETURN
278 END IF
279 *
280 * Quick return if possible
281 *
282 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
283 SCALE = 1
284 IF( NOTRAN ) THEN
285 IF( IJOB.NE.0 ) THEN
286 DIF = 0
287 END IF
288 END IF
289 RETURN
290 END IF
291 *
292 * Determine optimal block sizes MB and NB
293 *
294 MB = ILAENV( 2, 'ZTGSYL', TRANS, M, N, -1, -1 )
295 NB = ILAENV( 5, 'ZTGSYL', TRANS, M, N, -1, -1 )
296 *
297 ISOLVE = 1
298 IFUNC = 0
299 IF( NOTRAN ) THEN
300 IF( IJOB.GE.3 ) THEN
301 IFUNC = IJOB - 2
302 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC )
303 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF )
304 ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN
305 ISOLVE = 2
306 END IF
307 END IF
308 *
309 IF( ( MB.LE.1 .AND. NB.LE.1 ) .OR. ( MB.GE.M .AND. NB.GE.N ) )
310 $ THEN
311 *
312 * Use unblocked Level 2 solver
313 *
314 DO 30 IROUND = 1, ISOLVE
315 *
316 SCALE = ONE
317 DSCALE = ZERO
318 DSUM = ONE
319 PQ = M*N
320 CALL ZTGSY2( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D,
321 $ LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE,
322 $ INFO )
323 IF( DSCALE.NE.ZERO ) THEN
324 IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
325 DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
326 ELSE
327 DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
328 END IF
329 END IF
330 IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
331 IF( NOTRAN ) THEN
332 IFUNC = IJOB
333 END IF
334 SCALE2 = SCALE
335 CALL ZLACPY( 'F', M, N, C, LDC, WORK, M )
336 CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
337 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC )
338 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF )
339 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
340 CALL ZLACPY( 'F', M, N, WORK, M, C, LDC )
341 CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
342 SCALE = SCALE2
343 END IF
344 30 CONTINUE
345 *
346 RETURN
347 *
348 END IF
349 *
350 * Determine block structure of A
351 *
352 P = 0
353 I = 1
354 40 CONTINUE
355 IF( I.GT.M )
356 $ GO TO 50
357 P = P + 1
358 IWORK( P ) = I
359 I = I + MB
360 IF( I.GE.M )
361 $ GO TO 50
362 GO TO 40
363 50 CONTINUE
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 *
376 Q = Q + 1
377 IWORK( Q ) = J
378 J = J + NB
379 IF( J.GE.N )
380 $ GO TO 70
381 GO TO 60
382 *
383 70 CONTINUE
384 IWORK( Q+1 ) = N + 1
385 IF( IWORK( Q ).EQ.IWORK( Q+1 ) )
386 $ Q = Q - 1
387 *
388 IF( NOTRAN ) THEN
389 DO 150 IROUND = 1, ISOLVE
390 *
391 * Solve (I, J) - subsystem
392 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
393 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
394 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
395 *
396 PQ = 0
397 SCALE = ONE
398 DSCALE = ZERO
399 DSUM = ONE
400 DO 130 J = P + 2, Q
401 JS = IWORK( J )
402 JE = IWORK( J+1 ) - 1
403 NB = JE - JS + 1
404 DO 120 I = P, 1, -1
405 IS = IWORK( I )
406 IE = IWORK( I+1 ) - 1
407 MB = IE - IS + 1
408 CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
409 $ B( JS, JS ), LDB, C( IS, JS ), LDC,
410 $ D( IS, IS ), LDD, E( JS, JS ), LDE,
411 $ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
412 $ LINFO )
413 IF( LINFO.GT.0 )
414 $ INFO = LINFO
415 PQ = PQ + MB*NB
416 IF( SCALOC.NE.ONE ) THEN
417 DO 80 K = 1, JS - 1
418 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
419 $ C( 1, K ), 1 )
420 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
421 $ F( 1, K ), 1 )
422 80 CONTINUE
423 DO 90 K = JS, JE
424 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
425 $ C( 1, K ), 1 )
426 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
427 $ F( 1, K ), 1 )
428 90 CONTINUE
429 DO 100 K = JS, JE
430 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
431 $ C( IE+1, K ), 1 )
432 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
433 $ F( IE+1, K ), 1 )
434 100 CONTINUE
435 DO 110 K = JE + 1, N
436 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
437 $ C( 1, K ), 1 )
438 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
439 $ F( 1, K ), 1 )
440 110 CONTINUE
441 SCALE = SCALE*SCALOC
442 END IF
443 *
444 * Substitute R(I,J) and L(I,J) into remaining equation.
445 *
446 IF( I.GT.1 ) THEN
447 CALL ZGEMM( 'N', 'N', IS-1, NB, MB,
448 $ DCMPLX( -ONE, ZERO ), A( 1, IS ), LDA,
449 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
450 $ C( 1, JS ), LDC )
451 CALL ZGEMM( 'N', 'N', IS-1, NB, MB,
452 $ DCMPLX( -ONE, ZERO ), D( 1, IS ), LDD,
453 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
454 $ F( 1, JS ), LDF )
455 END IF
456 IF( J.LT.Q ) THEN
457 CALL ZGEMM( 'N', 'N', MB, N-JE, NB,
458 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
459 $ B( JS, JE+1 ), LDB,
460 $ DCMPLX( ONE, ZERO ), C( IS, JE+1 ),
461 $ LDC )
462 CALL ZGEMM( 'N', 'N', MB, N-JE, NB,
463 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
464 $ E( JS, JE+1 ), LDE,
465 $ DCMPLX( ONE, ZERO ), F( IS, JE+1 ),
466 $ LDF )
467 END IF
468 120 CONTINUE
469 130 CONTINUE
470 IF( DSCALE.NE.ZERO ) THEN
471 IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
472 DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
473 ELSE
474 DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
475 END IF
476 END IF
477 IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
478 IF( NOTRAN ) THEN
479 IFUNC = IJOB
480 END IF
481 SCALE2 = SCALE
482 CALL ZLACPY( 'F', M, N, C, LDC, WORK, M )
483 CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
484 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC )
485 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF )
486 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
487 CALL ZLACPY( 'F', M, N, WORK, M, C, LDC )
488 CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
489 SCALE = SCALE2
490 END IF
491 150 CONTINUE
492 ELSE
493 *
494 * Solve transposed (I, J)-subsystem
495 * A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
496 * R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
497 * for I = 1,2,..., P; J = Q, Q-1,..., 1
498 *
499 SCALE = ONE
500 DO 210 I = 1, P
501 IS = IWORK( I )
502 IE = IWORK( I+1 ) - 1
503 MB = IE - IS + 1
504 DO 200 J = Q, P + 2, -1
505 JS = IWORK( J )
506 JE = IWORK( J+1 ) - 1
507 NB = JE - JS + 1
508 CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
509 $ B( JS, JS ), LDB, C( IS, JS ), LDC,
510 $ D( IS, IS ), LDD, E( JS, JS ), LDE,
511 $ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
512 $ LINFO )
513 IF( LINFO.GT.0 )
514 $ INFO = LINFO
515 IF( SCALOC.NE.ONE ) THEN
516 DO 160 K = 1, JS - 1
517 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
518 $ 1 )
519 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
520 $ 1 )
521 160 CONTINUE
522 DO 170 K = JS, JE
523 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
524 $ C( 1, K ), 1 )
525 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ),
526 $ F( 1, K ), 1 )
527 170 CONTINUE
528 DO 180 K = JS, JE
529 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
530 $ C( IE+1, K ), 1 )
531 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ),
532 $ F( IE+1, K ), 1 )
533 180 CONTINUE
534 DO 190 K = JE + 1, N
535 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
536 $ 1 )
537 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
538 $ 1 )
539 190 CONTINUE
540 SCALE = SCALE*SCALOC
541 END IF
542 *
543 * Substitute R(I,J) and L(I,J) into remaining equation.
544 *
545 IF( J.GT.P+2 ) THEN
546 CALL ZGEMM( 'N', 'C', MB, JS-1, NB,
547 $ DCMPLX( ONE, ZERO ), C( IS, JS ), LDC,
548 $ B( 1, JS ), LDB, DCMPLX( ONE, ZERO ),
549 $ F( IS, 1 ), LDF )
550 CALL ZGEMM( 'N', 'C', MB, JS-1, NB,
551 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF,
552 $ E( 1, JS ), LDE, DCMPLX( ONE, ZERO ),
553 $ F( IS, 1 ), LDF )
554 END IF
555 IF( I.LT.P ) THEN
556 CALL ZGEMM( 'C', 'N', M-IE, NB, MB,
557 $ DCMPLX( -ONE, ZERO ), A( IS, IE+1 ), LDA,
558 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ),
559 $ C( IE+1, JS ), LDC )
560 CALL ZGEMM( 'C', 'N', M-IE, NB, MB,
561 $ DCMPLX( -ONE, ZERO ), D( IS, IE+1 ), LDD,
562 $ F( IS, JS ), LDF, DCMPLX( ONE, ZERO ),
563 $ C( IE+1, JS ), LDC )
564 END IF
565 200 CONTINUE
566 210 CONTINUE
567 END IF
568 *
569 WORK( 1 ) = LWMIN
570 *
571 RETURN
572 *
573 * End of ZTGSYL
574 *
575 END