Annotation of rpl/lapack/lapack/dtgsyl.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtgsyl.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack