![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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.2) -- 6: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 8: * November 2006 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', Im) ] (2) 44: * [ kron(In, D) -kron(E', Im) ]. 45: * 46: * Here Ik is the identity matrix of size k and X' 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'*y = scale*b, 50: * which is equivalent to solve for R and L in 51: * 52: * A' * R + D' * L = scale * C (3) 53: * R * B' + L * E' = 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)' * 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 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