1 SUBROUTINE DTGSYL( 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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
19 $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
20 $ WORK( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * DTGSYL 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 real entries. (A, D) and (B, E) must be in
34 * generalized (real) Schur canonical form, i.e. A, B are upper quasi
35 * triangular and D, E are upper triangular.
36 *
37 * The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
38 * scaling factor chosen to avoid overflow.
39 *
40 * In matrix notation (1) is equivalent to solve Zx = scale b, where
41 * Z is defined as
42 *
43 * Z = [ kron(In, A) -kron(B**T, Im) ] (2)
44 * [ kron(In, D) -kron(E**T, Im) ].
45 *
46 * Here Ik is the identity matrix of size k and X**T is the transpose of
47 * X. kron(X, Y) is the Kronecker product between the matrices X and Y.
48 *
49 * If TRANS = 'T', DTGSYL solves the transposed system Z**T*y = scale*b,
50 * which is equivalent to solve for R and L in
51 *
52 * A**T * R + D**T * L = scale * C (3)
53 * R * B**T + L * E**T = scale * -F
54 *
55 * This case (TRANS = 'T') 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 DLACON.
58 *
59 * If IJOB >= 1, DTGSYL computes a Frobenius norm-based estimate
60 * of Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
61 * reciprocal of the smallest singular value of Z. See [1-2] for more
62 * information.
63 *
64 * This is a level 3 BLAS algorithm.
65 *
66 * Arguments
67 * =========
68 *
69 * TRANS (input) CHARACTER*1
70 * = 'N', solve the generalized Sylvester equation (1).
71 * = 'T', solve the 'transposed' system (3).
72 *
73 * IJOB (input) INTEGER
74 * Specifies what kind of functionality to be performed.
75 * =0: solve (1) only.
76 * =1: The functionality of 0 and 3.
77 * =2: The functionality of 0 and 4.
78 * =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
79 * (look ahead strategy IJOB = 1 is used).
80 * =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
81 * ( DGECON on sub-systems is used ).
82 * Not referenced if TRANS = 'T'.
83 *
84 * M (input) INTEGER
85 * The order of the matrices A and D, and the row dimension of
86 * the matrices C, F, R and L.
87 *
88 * N (input) INTEGER
89 * The order of the matrices B and E, and the column dimension
90 * of the matrices C, F, R and L.
91 *
92 * A (input) DOUBLE PRECISION array, dimension (LDA, M)
93 * The upper quasi triangular matrix A.
94 *
95 * LDA (input) INTEGER
96 * The leading dimension of the array A. LDA >= max(1, M).
97 *
98 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
99 * 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION
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) DOUBLE PRECISION
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) DOUBLE PRECISION 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 DCOPY by calls to DLASET.
197 * Sven Hammarling, 1/5/02.
198 *
199 * .. Parameters ..
200 DOUBLE PRECISION ZERO, ONE
201 PARAMETER ( ZERO = 0.0D+0, ONE = 1.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, PPQQ, 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 DGEMM, DLACPY, DLASET, DSCAL, DTGSY2, XERBLA
216 * ..
217 * .. Intrinsic Functions ..
218 INTRINSIC DBLE, 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, '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( 'DTGSYL', -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, 'DTGSYL', TRANS, M, N, -1, -1 )
294 NB = ILAENV( 5, 'DTGSYL', 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 DLASET( 'F', M, N, ZERO, ZERO, C, LDC )
302 CALL DLASET( 'F', M, N, ZERO, ZERO, F, LDF )
303 ELSE IF( IJOB.GE.1 ) 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 DTGSY2( 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( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
324 ELSE
325 DIF = SQRT( DBLE( 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 DLACPY( 'F', M, N, C, LDC, WORK, M )
335 CALL DLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
336 CALL DLASET( 'F', M, N, ZERO, ZERO, C, LDC )
337 CALL DLASET( 'F', M, N, ZERO, ZERO, F, LDF )
338 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
339 CALL DLACPY( 'F', M, N, WORK, M, C, LDC )
340 CALL DLACPY( '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 DTGSY2( 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 DSCAL( M, SCALOC, C( 1, K ), 1 )
424 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
425 80 CONTINUE
426 DO 90 K = JS, JE
427 CALL DSCAL( IS-1, SCALOC, C( 1, K ), 1 )
428 CALL DSCAL( IS-1, SCALOC, F( 1, K ), 1 )
429 90 CONTINUE
430 DO 100 K = JS, JE
431 CALL DSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
432 CALL DSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
433 100 CONTINUE
434 DO 110 K = JE + 1, N
435 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
436 CALL DSCAL( 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 DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
446 $ A( 1, IS ), LDA, C( IS, JS ), LDC, ONE,
447 $ C( 1, JS ), LDC )
448 CALL DGEMM( '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 DGEMM( '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 DGEMM( '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( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
465 ELSE
466 DIF = SQRT( DBLE( 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 DLACPY( 'F', M, N, C, LDC, WORK, M )
475 CALL DLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
476 CALL DLASET( 'F', M, N, ZERO, ZERO, C, LDC )
477 CALL DLASET( 'F', M, N, ZERO, ZERO, F, LDF )
478 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
479 CALL DLACPY( 'F', M, N, WORK, M, C, LDC )
480 CALL DLACPY( '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)**T * R(I, J) + D(I, I)**T * L(I, J) = C(I, J)
489 * R(I, J) * B(J, J)**T + L(I, J) * E(J, J)**T = -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 DTGSY2( 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 DSCAL( M, SCALOC, C( 1, K ), 1 )
511 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
512 160 CONTINUE
513 DO 170 K = JS, JE
514 CALL DSCAL( IS-1, SCALOC, C( 1, K ), 1 )
515 CALL DSCAL( IS-1, SCALOC, F( 1, K ), 1 )
516 170 CONTINUE
517 DO 180 K = JS, JE
518 CALL DSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
519 CALL DSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
520 180 CONTINUE
521 DO 190 K = JE + 1, N
522 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
523 CALL DSCAL( 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 DGEMM( 'N', 'T', MB, JS-1, NB, ONE, C( IS, JS ),
532 $ LDC, B( 1, JS ), LDB, ONE, F( IS, 1 ),
533 $ LDF )
534 CALL DGEMM( '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 DGEMM( '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 DGEMM( '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 DTGSYL
556 *
557 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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
19 $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
20 $ WORK( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * DTGSYL 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 real entries. (A, D) and (B, E) must be in
34 * generalized (real) Schur canonical form, i.e. A, B are upper quasi
35 * triangular and D, E are upper triangular.
36 *
37 * The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
38 * scaling factor chosen to avoid overflow.
39 *
40 * In matrix notation (1) is equivalent to solve Zx = scale b, where
41 * Z is defined as
42 *
43 * Z = [ kron(In, A) -kron(B**T, Im) ] (2)
44 * [ kron(In, D) -kron(E**T, Im) ].
45 *
46 * Here Ik is the identity matrix of size k and X**T is the transpose of
47 * X. kron(X, Y) is the Kronecker product between the matrices X and Y.
48 *
49 * If TRANS = 'T', DTGSYL solves the transposed system Z**T*y = scale*b,
50 * which is equivalent to solve for R and L in
51 *
52 * A**T * R + D**T * L = scale * C (3)
53 * R * B**T + L * E**T = scale * -F
54 *
55 * This case (TRANS = 'T') 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 DLACON.
58 *
59 * If IJOB >= 1, DTGSYL computes a Frobenius norm-based estimate
60 * of Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
61 * reciprocal of the smallest singular value of Z. See [1-2] for more
62 * information.
63 *
64 * This is a level 3 BLAS algorithm.
65 *
66 * Arguments
67 * =========
68 *
69 * TRANS (input) CHARACTER*1
70 * = 'N', solve the generalized Sylvester equation (1).
71 * = 'T', solve the 'transposed' system (3).
72 *
73 * IJOB (input) INTEGER
74 * Specifies what kind of functionality to be performed.
75 * =0: solve (1) only.
76 * =1: The functionality of 0 and 3.
77 * =2: The functionality of 0 and 4.
78 * =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
79 * (look ahead strategy IJOB = 1 is used).
80 * =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
81 * ( DGECON on sub-systems is used ).
82 * Not referenced if TRANS = 'T'.
83 *
84 * M (input) INTEGER
85 * The order of the matrices A and D, and the row dimension of
86 * the matrices C, F, R and L.
87 *
88 * N (input) INTEGER
89 * The order of the matrices B and E, and the column dimension
90 * of the matrices C, F, R and L.
91 *
92 * A (input) DOUBLE PRECISION array, dimension (LDA, M)
93 * The upper quasi triangular matrix A.
94 *
95 * LDA (input) INTEGER
96 * The leading dimension of the array A. LDA >= max(1, M).
97 *
98 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
99 * 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION
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) DOUBLE PRECISION
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) DOUBLE PRECISION 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 DCOPY by calls to DLASET.
197 * Sven Hammarling, 1/5/02.
198 *
199 * .. Parameters ..
200 DOUBLE PRECISION ZERO, ONE
201 PARAMETER ( ZERO = 0.0D+0, ONE = 1.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, PPQQ, 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 DGEMM, DLACPY, DLASET, DSCAL, DTGSY2, XERBLA
216 * ..
217 * .. Intrinsic Functions ..
218 INTRINSIC DBLE, 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, '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( 'DTGSYL', -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, 'DTGSYL', TRANS, M, N, -1, -1 )
294 NB = ILAENV( 5, 'DTGSYL', 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 DLASET( 'F', M, N, ZERO, ZERO, C, LDC )
302 CALL DLASET( 'F', M, N, ZERO, ZERO, F, LDF )
303 ELSE IF( IJOB.GE.1 ) 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 DTGSY2( 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( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
324 ELSE
325 DIF = SQRT( DBLE( 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 DLACPY( 'F', M, N, C, LDC, WORK, M )
335 CALL DLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
336 CALL DLASET( 'F', M, N, ZERO, ZERO, C, LDC )
337 CALL DLASET( 'F', M, N, ZERO, ZERO, F, LDF )
338 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
339 CALL DLACPY( 'F', M, N, WORK, M, C, LDC )
340 CALL DLACPY( '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 DTGSY2( 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 DSCAL( M, SCALOC, C( 1, K ), 1 )
424 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
425 80 CONTINUE
426 DO 90 K = JS, JE
427 CALL DSCAL( IS-1, SCALOC, C( 1, K ), 1 )
428 CALL DSCAL( IS-1, SCALOC, F( 1, K ), 1 )
429 90 CONTINUE
430 DO 100 K = JS, JE
431 CALL DSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
432 CALL DSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
433 100 CONTINUE
434 DO 110 K = JE + 1, N
435 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
436 CALL DSCAL( 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 DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
446 $ A( 1, IS ), LDA, C( IS, JS ), LDC, ONE,
447 $ C( 1, JS ), LDC )
448 CALL DGEMM( '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 DGEMM( '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 DGEMM( '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( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
465 ELSE
466 DIF = SQRT( DBLE( 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 DLACPY( 'F', M, N, C, LDC, WORK, M )
475 CALL DLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
476 CALL DLASET( 'F', M, N, ZERO, ZERO, C, LDC )
477 CALL DLASET( 'F', M, N, ZERO, ZERO, F, LDF )
478 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
479 CALL DLACPY( 'F', M, N, WORK, M, C, LDC )
480 CALL DLACPY( '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)**T * R(I, J) + D(I, I)**T * L(I, J) = C(I, J)
489 * R(I, J) * B(J, J)**T + L(I, J) * E(J, J)**T = -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 DTGSY2( 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 DSCAL( M, SCALOC, C( 1, K ), 1 )
511 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
512 160 CONTINUE
513 DO 170 K = JS, JE
514 CALL DSCAL( IS-1, SCALOC, C( 1, K ), 1 )
515 CALL DSCAL( IS-1, SCALOC, F( 1, K ), 1 )
516 170 CONTINUE
517 DO 180 K = JS, JE
518 CALL DSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
519 CALL DSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
520 180 CONTINUE
521 DO 190 K = JE + 1, N
522 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
523 CALL DSCAL( 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 DGEMM( 'N', 'T', MB, JS-1, NB, ONE, C( IS, JS ),
532 $ LDC, B( 1, JS ), LDB, ONE, F( IS, 1 ),
533 $ LDF )
534 CALL DGEMM( '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 DGEMM( '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 DGEMM( '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 DTGSYL
556 *
557 END