![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 2: $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, 3: $ IWORK, PQ, INFO ) 4: * 5: * -- LAPACK auxiliary 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: * January 2007 9: * 10: * .. Scalar Arguments .. 11: CHARACTER TRANS 12: INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N, 13: $ PQ 14: DOUBLE PRECISION RDSCAL, RDSUM, 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: * .. 21: * 22: * Purpose 23: * ======= 24: * 25: * DTGSY2 solves the generalized Sylvester equation: 26: * 27: * A * R - L * B = scale * C (1) 28: * D * R - L * E = scale * F, 29: * 30: * using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices, 31: * (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M, 32: * N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E) 33: * must be in generalized Schur canonical form, i.e. A, B are upper 34: * quasi triangular and D, E are upper triangular. The solution (R, L) 35: * overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor 36: * chosen to avoid overflow. 37: * 38: * In matrix notation solving equation (1) corresponds to solve 39: * Z*x = scale*b, where Z is defined as 40: * 41: * Z = [ kron(In, A) -kron(B', Im) ] (2) 42: * [ kron(In, D) -kron(E', Im) ], 43: * 44: * Ik is the identity matrix of size k and X' is the transpose of X. 45: * kron(X, Y) is the Kronecker product between the matrices X and Y. 46: * In the process of solving (1), we solve a number of such systems 47: * where Dim(In), Dim(In) = 1 or 2. 48: * 49: * If TRANS = 'T', solve the transposed system Z'*y = scale*b for y, 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 is used to compute an estimate of Dif[(A, D), (B, E)] = 56: * sigma_min(Z) using reverse communicaton with DLACON. 57: * 58: * DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL 59: * of an upper bound on the separation between to matrix pairs. Then 60: * the input (A, D), (B, E) are sub-pencils of the matrix pair in 61: * DTGSYL. See DTGSYL for details. 62: * 63: * Arguments 64: * ========= 65: * 66: * TRANS (input) CHARACTER*1 67: * = 'N', solve the generalized Sylvester equation (1). 68: * = 'T': solve the 'transposed' system (3). 69: * 70: * IJOB (input) INTEGER 71: * Specifies what kind of functionality to be performed. 72: * = 0: solve (1) only. 73: * = 1: A contribution from this subsystem to a Frobenius 74: * norm-based estimate of the separation between two matrix 75: * pairs is computed. (look ahead strategy is used). 76: * = 2: A contribution from this subsystem to a Frobenius 77: * norm-based estimate of the separation between two matrix 78: * pairs is computed. (DGECON on sub-systems is used.) 79: * Not referenced if TRANS = 'T'. 80: * 81: * M (input) INTEGER 82: * On entry, M specifies the order of A and D, and the row 83: * dimension of C, F, R and L. 84: * 85: * N (input) INTEGER 86: * On entry, N specifies the order of B and E, and the column 87: * dimension of C, F, R and L. 88: * 89: * A (input) DOUBLE PRECISION array, dimension (LDA, M) 90: * On entry, A contains an upper quasi triangular matrix. 91: * 92: * LDA (input) INTEGER 93: * The leading dimension of the matrix A. LDA >= max(1, M). 94: * 95: * B (input) DOUBLE PRECISION array, dimension (LDB, N) 96: * On entry, B contains an upper quasi triangular matrix. 97: * 98: * LDB (input) INTEGER 99: * The leading dimension of the matrix B. LDB >= max(1, N). 100: * 101: * C (input/output) DOUBLE PRECISION array, dimension (LDC, N) 102: * On entry, C contains the right-hand-side of the first matrix 103: * equation in (1). 104: * On exit, if IJOB = 0, C has been overwritten by the 105: * solution R. 106: * 107: * LDC (input) INTEGER 108: * The leading dimension of the matrix C. LDC >= max(1, M). 109: * 110: * D (input) DOUBLE PRECISION array, dimension (LDD, M) 111: * On entry, D contains an upper triangular matrix. 112: * 113: * LDD (input) INTEGER 114: * The leading dimension of the matrix D. LDD >= max(1, M). 115: * 116: * E (input) DOUBLE PRECISION array, dimension (LDE, N) 117: * On entry, E contains an upper triangular matrix. 118: * 119: * LDE (input) INTEGER 120: * The leading dimension of the matrix E. LDE >= max(1, N). 121: * 122: * F (input/output) DOUBLE PRECISION array, dimension (LDF, N) 123: * On entry, F contains the right-hand-side of the second matrix 124: * equation in (1). 125: * On exit, if IJOB = 0, F has been overwritten by the 126: * solution L. 127: * 128: * LDF (input) INTEGER 129: * The leading dimension of the matrix F. LDF >= max(1, M). 130: * 131: * SCALE (output) DOUBLE PRECISION 132: * On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions 133: * R and L (C and F on entry) will hold the solutions to a 134: * slightly perturbed system but the input matrices A, B, D and 135: * E have not been changed. If SCALE = 0, R and L will hold the 136: * solutions to the homogeneous system with C = F = 0. Normally, 137: * SCALE = 1. 138: * 139: * RDSUM (input/output) DOUBLE PRECISION 140: * On entry, the sum of squares of computed contributions to 141: * the Dif-estimate under computation by DTGSYL, where the 142: * scaling factor RDSCAL (see below) has been factored out. 143: * On exit, the corresponding sum of squares updated with the 144: * contributions from the current sub-system. 145: * If TRANS = 'T' RDSUM is not touched. 146: * NOTE: RDSUM only makes sense when DTGSY2 is called by DTGSYL. 147: * 148: * RDSCAL (input/output) DOUBLE PRECISION 149: * On entry, scaling factor used to prevent overflow in RDSUM. 150: * On exit, RDSCAL is updated w.r.t. the current contributions 151: * in RDSUM. 152: * If TRANS = 'T', RDSCAL is not touched. 153: * NOTE: RDSCAL only makes sense when DTGSY2 is called by 154: * DTGSYL. 155: * 156: * IWORK (workspace) INTEGER array, dimension (M+N+2) 157: * 158: * PQ (output) INTEGER 159: * On exit, the number of subsystems (of size 2-by-2, 4-by-4 and 160: * 8-by-8) solved by this routine. 161: * 162: * INFO (output) INTEGER 163: * On exit, if INFO is set to 164: * =0: Successful exit 165: * <0: If INFO = -i, the i-th argument had an illegal value. 166: * >0: The matrix pairs (A, D) and (B, E) have common or very 167: * close eigenvalues. 168: * 169: * Further Details 170: * =============== 171: * 172: * Based on contributions by 173: * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 174: * Umea University, S-901 87 Umea, Sweden. 175: * 176: * ===================================================================== 177: * Replaced various illegal calls to DCOPY by calls to DLASET. 178: * Sven Hammarling, 27/5/02. 179: * 180: * .. Parameters .. 181: INTEGER LDZ 182: PARAMETER ( LDZ = 8 ) 183: DOUBLE PRECISION ZERO, ONE 184: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 185: * .. 186: * .. Local Scalars .. 187: LOGICAL NOTRAN 188: INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1, 189: $ K, MB, NB, P, Q, ZDIM 190: DOUBLE PRECISION ALPHA, SCALOC 191: * .. 192: * .. Local Arrays .. 193: INTEGER IPIV( LDZ ), JPIV( LDZ ) 194: DOUBLE PRECISION RHS( LDZ ), Z( LDZ, LDZ ) 195: * .. 196: * .. External Functions .. 197: LOGICAL LSAME 198: EXTERNAL LSAME 199: * .. 200: * .. External Subroutines .. 201: EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DGER, DGESC2, 202: $ DGETC2, DLASET, DLATDF, DSCAL, XERBLA 203: * .. 204: * .. Intrinsic Functions .. 205: INTRINSIC MAX 206: * .. 207: * .. Executable Statements .. 208: * 209: * Decode and test input parameters 210: * 211: INFO = 0 212: IERR = 0 213: NOTRAN = LSAME( TRANS, 'N' ) 214: IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN 215: INFO = -1 216: ELSE IF( NOTRAN ) THEN 217: IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN 218: INFO = -2 219: END IF 220: END IF 221: IF( INFO.EQ.0 ) THEN 222: IF( M.LE.0 ) THEN 223: INFO = -3 224: ELSE IF( N.LE.0 ) THEN 225: INFO = -4 226: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 227: INFO = -5 228: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 229: INFO = -8 230: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 231: INFO = -10 232: ELSE IF( LDD.LT.MAX( 1, M ) ) THEN 233: INFO = -12 234: ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 235: INFO = -14 236: ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 237: INFO = -16 238: END IF 239: END IF 240: IF( INFO.NE.0 ) THEN 241: CALL XERBLA( 'DTGSY2', -INFO ) 242: RETURN 243: END IF 244: * 245: * Determine block structure of A 246: * 247: PQ = 0 248: P = 0 249: I = 1 250: 10 CONTINUE 251: IF( I.GT.M ) 252: $ GO TO 20 253: P = P + 1 254: IWORK( P ) = I 255: IF( I.EQ.M ) 256: $ GO TO 20 257: IF( A( I+1, I ).NE.ZERO ) THEN 258: I = I + 2 259: ELSE 260: I = I + 1 261: END IF 262: GO TO 10 263: 20 CONTINUE 264: IWORK( P+1 ) = M + 1 265: * 266: * Determine block structure of B 267: * 268: Q = P + 1 269: J = 1 270: 30 CONTINUE 271: IF( J.GT.N ) 272: $ GO TO 40 273: Q = Q + 1 274: IWORK( Q ) = J 275: IF( J.EQ.N ) 276: $ GO TO 40 277: IF( B( J+1, J ).NE.ZERO ) THEN 278: J = J + 2 279: ELSE 280: J = J + 1 281: END IF 282: GO TO 30 283: 40 CONTINUE 284: IWORK( Q+1 ) = N + 1 285: PQ = P*( Q-P-1 ) 286: * 287: IF( NOTRAN ) THEN 288: * 289: * Solve (I, J) - subsystem 290: * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) 291: * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) 292: * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q 293: * 294: SCALE = ONE 295: SCALOC = ONE 296: DO 120 J = P + 2, Q 297: JS = IWORK( J ) 298: JSP1 = JS + 1 299: JE = IWORK( J+1 ) - 1 300: NB = JE - JS + 1 301: DO 110 I = P, 1, -1 302: * 303: IS = IWORK( I ) 304: ISP1 = IS + 1 305: IE = IWORK( I+1 ) - 1 306: MB = IE - IS + 1 307: ZDIM = MB*NB*2 308: * 309: IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN 310: * 311: * Build a 2-by-2 system Z * x = RHS 312: * 313: Z( 1, 1 ) = A( IS, IS ) 314: Z( 2, 1 ) = D( IS, IS ) 315: Z( 1, 2 ) = -B( JS, JS ) 316: Z( 2, 2 ) = -E( JS, JS ) 317: * 318: * Set up right hand side(s) 319: * 320: RHS( 1 ) = C( IS, JS ) 321: RHS( 2 ) = F( IS, JS ) 322: * 323: * Solve Z * x = RHS 324: * 325: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 326: IF( IERR.GT.0 ) 327: $ INFO = IERR 328: * 329: IF( IJOB.EQ.0 ) THEN 330: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 331: $ SCALOC ) 332: IF( SCALOC.NE.ONE ) THEN 333: DO 50 K = 1, N 334: CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 335: CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 336: 50 CONTINUE 337: SCALE = SCALE*SCALOC 338: END IF 339: ELSE 340: CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 341: $ RDSCAL, IPIV, JPIV ) 342: END IF 343: * 344: * Unpack solution vector(s) 345: * 346: C( IS, JS ) = RHS( 1 ) 347: F( IS, JS ) = RHS( 2 ) 348: * 349: * Substitute R(I, J) and L(I, J) into remaining 350: * equation. 351: * 352: IF( I.GT.1 ) THEN 353: ALPHA = -RHS( 1 ) 354: CALL DAXPY( IS-1, ALPHA, A( 1, IS ), 1, C( 1, JS ), 355: $ 1 ) 356: CALL DAXPY( IS-1, ALPHA, D( 1, IS ), 1, F( 1, JS ), 357: $ 1 ) 358: END IF 359: IF( J.LT.Q ) THEN 360: CALL DAXPY( N-JE, RHS( 2 ), B( JS, JE+1 ), LDB, 361: $ C( IS, JE+1 ), LDC ) 362: CALL DAXPY( N-JE, RHS( 2 ), E( JS, JE+1 ), LDE, 363: $ F( IS, JE+1 ), LDF ) 364: END IF 365: * 366: ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN 367: * 368: * Build a 4-by-4 system Z * x = RHS 369: * 370: Z( 1, 1 ) = A( IS, IS ) 371: Z( 2, 1 ) = ZERO 372: Z( 3, 1 ) = D( IS, IS ) 373: Z( 4, 1 ) = ZERO 374: * 375: Z( 1, 2 ) = ZERO 376: Z( 2, 2 ) = A( IS, IS ) 377: Z( 3, 2 ) = ZERO 378: Z( 4, 2 ) = D( IS, IS ) 379: * 380: Z( 1, 3 ) = -B( JS, JS ) 381: Z( 2, 3 ) = -B( JS, JSP1 ) 382: Z( 3, 3 ) = -E( JS, JS ) 383: Z( 4, 3 ) = -E( JS, JSP1 ) 384: * 385: Z( 1, 4 ) = -B( JSP1, JS ) 386: Z( 2, 4 ) = -B( JSP1, JSP1 ) 387: Z( 3, 4 ) = ZERO 388: Z( 4, 4 ) = -E( JSP1, JSP1 ) 389: * 390: * Set up right hand side(s) 391: * 392: RHS( 1 ) = C( IS, JS ) 393: RHS( 2 ) = C( IS, JSP1 ) 394: RHS( 3 ) = F( IS, JS ) 395: RHS( 4 ) = F( IS, JSP1 ) 396: * 397: * Solve Z * x = RHS 398: * 399: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 400: IF( IERR.GT.0 ) 401: $ INFO = IERR 402: * 403: IF( IJOB.EQ.0 ) THEN 404: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 405: $ SCALOC ) 406: IF( SCALOC.NE.ONE ) THEN 407: DO 60 K = 1, N 408: CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 409: CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 410: 60 CONTINUE 411: SCALE = SCALE*SCALOC 412: END IF 413: ELSE 414: CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 415: $ RDSCAL, IPIV, JPIV ) 416: END IF 417: * 418: * Unpack solution vector(s) 419: * 420: C( IS, JS ) = RHS( 1 ) 421: C( IS, JSP1 ) = RHS( 2 ) 422: F( IS, JS ) = RHS( 3 ) 423: F( IS, JSP1 ) = RHS( 4 ) 424: * 425: * Substitute R(I, J) and L(I, J) into remaining 426: * equation. 427: * 428: IF( I.GT.1 ) THEN 429: CALL DGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ), 430: $ 1, C( 1, JS ), LDC ) 431: CALL DGER( IS-1, NB, -ONE, D( 1, IS ), 1, RHS( 1 ), 432: $ 1, F( 1, JS ), LDF ) 433: END IF 434: IF( J.LT.Q ) THEN 435: CALL DAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB, 436: $ C( IS, JE+1 ), LDC ) 437: CALL DAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE, 438: $ F( IS, JE+1 ), LDF ) 439: CALL DAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB, 440: $ C( IS, JE+1 ), LDC ) 441: CALL DAXPY( N-JE, RHS( 4 ), E( JSP1, JE+1 ), LDE, 442: $ F( IS, JE+1 ), LDF ) 443: END IF 444: * 445: ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN 446: * 447: * Build a 4-by-4 system Z * x = RHS 448: * 449: Z( 1, 1 ) = A( IS, IS ) 450: Z( 2, 1 ) = A( ISP1, IS ) 451: Z( 3, 1 ) = D( IS, IS ) 452: Z( 4, 1 ) = ZERO 453: * 454: Z( 1, 2 ) = A( IS, ISP1 ) 455: Z( 2, 2 ) = A( ISP1, ISP1 ) 456: Z( 3, 2 ) = D( IS, ISP1 ) 457: Z( 4, 2 ) = D( ISP1, ISP1 ) 458: * 459: Z( 1, 3 ) = -B( JS, JS ) 460: Z( 2, 3 ) = ZERO 461: Z( 3, 3 ) = -E( JS, JS ) 462: Z( 4, 3 ) = ZERO 463: * 464: Z( 1, 4 ) = ZERO 465: Z( 2, 4 ) = -B( JS, JS ) 466: Z( 3, 4 ) = ZERO 467: Z( 4, 4 ) = -E( JS, JS ) 468: * 469: * Set up right hand side(s) 470: * 471: RHS( 1 ) = C( IS, JS ) 472: RHS( 2 ) = C( ISP1, JS ) 473: RHS( 3 ) = F( IS, JS ) 474: RHS( 4 ) = F( ISP1, JS ) 475: * 476: * Solve Z * x = RHS 477: * 478: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 479: IF( IERR.GT.0 ) 480: $ INFO = IERR 481: IF( IJOB.EQ.0 ) THEN 482: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 483: $ SCALOC ) 484: IF( SCALOC.NE.ONE ) THEN 485: DO 70 K = 1, N 486: CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 487: CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 488: 70 CONTINUE 489: SCALE = SCALE*SCALOC 490: END IF 491: ELSE 492: CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 493: $ RDSCAL, IPIV, JPIV ) 494: END IF 495: * 496: * Unpack solution vector(s) 497: * 498: C( IS, JS ) = RHS( 1 ) 499: C( ISP1, JS ) = RHS( 2 ) 500: F( IS, JS ) = RHS( 3 ) 501: F( ISP1, JS ) = RHS( 4 ) 502: * 503: * Substitute R(I, J) and L(I, J) into remaining 504: * equation. 505: * 506: IF( I.GT.1 ) THEN 507: CALL DGEMV( 'N', IS-1, MB, -ONE, A( 1, IS ), LDA, 508: $ RHS( 1 ), 1, ONE, C( 1, JS ), 1 ) 509: CALL DGEMV( 'N', IS-1, MB, -ONE, D( 1, IS ), LDD, 510: $ RHS( 1 ), 1, ONE, F( 1, JS ), 1 ) 511: END IF 512: IF( J.LT.Q ) THEN 513: CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1, 514: $ B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC ) 515: CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1, 516: $ E( JS, JE+1 ), LDE, F( IS, JE+1 ), LDF ) 517: END IF 518: * 519: ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN 520: * 521: * Build an 8-by-8 system Z * x = RHS 522: * 523: CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ ) 524: * 525: Z( 1, 1 ) = A( IS, IS ) 526: Z( 2, 1 ) = A( ISP1, IS ) 527: Z( 5, 1 ) = D( IS, IS ) 528: * 529: Z( 1, 2 ) = A( IS, ISP1 ) 530: Z( 2, 2 ) = A( ISP1, ISP1 ) 531: Z( 5, 2 ) = D( IS, ISP1 ) 532: Z( 6, 2 ) = D( ISP1, ISP1 ) 533: * 534: Z( 3, 3 ) = A( IS, IS ) 535: Z( 4, 3 ) = A( ISP1, IS ) 536: Z( 7, 3 ) = D( IS, IS ) 537: * 538: Z( 3, 4 ) = A( IS, ISP1 ) 539: Z( 4, 4 ) = A( ISP1, ISP1 ) 540: Z( 7, 4 ) = D( IS, ISP1 ) 541: Z( 8, 4 ) = D( ISP1, ISP1 ) 542: * 543: Z( 1, 5 ) = -B( JS, JS ) 544: Z( 3, 5 ) = -B( JS, JSP1 ) 545: Z( 5, 5 ) = -E( JS, JS ) 546: Z( 7, 5 ) = -E( JS, JSP1 ) 547: * 548: Z( 2, 6 ) = -B( JS, JS ) 549: Z( 4, 6 ) = -B( JS, JSP1 ) 550: Z( 6, 6 ) = -E( JS, JS ) 551: Z( 8, 6 ) = -E( JS, JSP1 ) 552: * 553: Z( 1, 7 ) = -B( JSP1, JS ) 554: Z( 3, 7 ) = -B( JSP1, JSP1 ) 555: Z( 7, 7 ) = -E( JSP1, JSP1 ) 556: * 557: Z( 2, 8 ) = -B( JSP1, JS ) 558: Z( 4, 8 ) = -B( JSP1, JSP1 ) 559: Z( 8, 8 ) = -E( JSP1, JSP1 ) 560: * 561: * Set up right hand side(s) 562: * 563: K = 1 564: II = MB*NB + 1 565: DO 80 JJ = 0, NB - 1 566: CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 ) 567: CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 ) 568: K = K + MB 569: II = II + MB 570: 80 CONTINUE 571: * 572: * Solve Z * x = RHS 573: * 574: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 575: IF( IERR.GT.0 ) 576: $ INFO = IERR 577: IF( IJOB.EQ.0 ) THEN 578: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 579: $ SCALOC ) 580: IF( SCALOC.NE.ONE ) THEN 581: DO 90 K = 1, N 582: CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 583: CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 584: 90 CONTINUE 585: SCALE = SCALE*SCALOC 586: END IF 587: ELSE 588: CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 589: $ RDSCAL, IPIV, JPIV ) 590: END IF 591: * 592: * Unpack solution vector(s) 593: * 594: K = 1 595: II = MB*NB + 1 596: DO 100 JJ = 0, NB - 1 597: CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 ) 598: CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 ) 599: K = K + MB 600: II = II + MB 601: 100 CONTINUE 602: * 603: * Substitute R(I, J) and L(I, J) into remaining 604: * equation. 605: * 606: IF( I.GT.1 ) THEN 607: CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE, 608: $ A( 1, IS ), LDA, RHS( 1 ), MB, ONE, 609: $ C( 1, JS ), LDC ) 610: CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE, 611: $ D( 1, IS ), LDD, RHS( 1 ), MB, ONE, 612: $ F( 1, JS ), LDF ) 613: END IF 614: IF( J.LT.Q ) THEN 615: K = MB*NB + 1 616: CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ), 617: $ MB, B( JS, JE+1 ), LDB, ONE, 618: $ C( IS, JE+1 ), LDC ) 619: CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ), 620: $ MB, E( JS, JE+1 ), LDE, ONE, 621: $ F( IS, JE+1 ), LDF ) 622: END IF 623: * 624: END IF 625: * 626: 110 CONTINUE 627: 120 CONTINUE 628: ELSE 629: * 630: * Solve (I, J) - subsystem 631: * A(I, I)' * R(I, J) + D(I, I)' * L(J, J) = C(I, J) 632: * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) 633: * for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1 634: * 635: SCALE = ONE 636: SCALOC = ONE 637: DO 200 I = 1, P 638: * 639: IS = IWORK( I ) 640: ISP1 = IS + 1 641: IE = IWORK ( I+1 ) - 1 642: MB = IE - IS + 1 643: DO 190 J = Q, P + 2, -1 644: * 645: JS = IWORK( J ) 646: JSP1 = JS + 1 647: JE = IWORK( J+1 ) - 1 648: NB = JE - JS + 1 649: ZDIM = MB*NB*2 650: IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN 651: * 652: * Build a 2-by-2 system Z' * x = RHS 653: * 654: Z( 1, 1 ) = A( IS, IS ) 655: Z( 2, 1 ) = -B( JS, JS ) 656: Z( 1, 2 ) = D( IS, IS ) 657: Z( 2, 2 ) = -E( JS, JS ) 658: * 659: * Set up right hand side(s) 660: * 661: RHS( 1 ) = C( IS, JS ) 662: RHS( 2 ) = F( IS, JS ) 663: * 664: * Solve Z' * x = RHS 665: * 666: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 667: IF( IERR.GT.0 ) 668: $ INFO = IERR 669: * 670: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 671: IF( SCALOC.NE.ONE ) THEN 672: DO 130 K = 1, N 673: CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 674: CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 675: 130 CONTINUE 676: SCALE = SCALE*SCALOC 677: END IF 678: * 679: * Unpack solution vector(s) 680: * 681: C( IS, JS ) = RHS( 1 ) 682: F( IS, JS ) = RHS( 2 ) 683: * 684: * Substitute R(I, J) and L(I, J) into remaining 685: * equation. 686: * 687: IF( J.GT.P+2 ) THEN 688: ALPHA = RHS( 1 ) 689: CALL DAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ), 690: $ LDF ) 691: ALPHA = RHS( 2 ) 692: CALL DAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ), 693: $ LDF ) 694: END IF 695: IF( I.LT.P ) THEN 696: ALPHA = -RHS( 1 ) 697: CALL DAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA, 698: $ C( IE+1, JS ), 1 ) 699: ALPHA = -RHS( 2 ) 700: CALL DAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD, 701: $ C( IE+1, JS ), 1 ) 702: END IF 703: * 704: ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN 705: * 706: * Build a 4-by-4 system Z' * x = RHS 707: * 708: Z( 1, 1 ) = A( IS, IS ) 709: Z( 2, 1 ) = ZERO 710: Z( 3, 1 ) = -B( JS, JS ) 711: Z( 4, 1 ) = -B( JSP1, JS ) 712: * 713: Z( 1, 2 ) = ZERO 714: Z( 2, 2 ) = A( IS, IS ) 715: Z( 3, 2 ) = -B( JS, JSP1 ) 716: Z( 4, 2 ) = -B( JSP1, JSP1 ) 717: * 718: Z( 1, 3 ) = D( IS, IS ) 719: Z( 2, 3 ) = ZERO 720: Z( 3, 3 ) = -E( JS, JS ) 721: Z( 4, 3 ) = ZERO 722: * 723: Z( 1, 4 ) = ZERO 724: Z( 2, 4 ) = D( IS, IS ) 725: Z( 3, 4 ) = -E( JS, JSP1 ) 726: Z( 4, 4 ) = -E( JSP1, JSP1 ) 727: * 728: * Set up right hand side(s) 729: * 730: RHS( 1 ) = C( IS, JS ) 731: RHS( 2 ) = C( IS, JSP1 ) 732: RHS( 3 ) = F( IS, JS ) 733: RHS( 4 ) = F( IS, JSP1 ) 734: * 735: * Solve Z' * x = RHS 736: * 737: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 738: IF( IERR.GT.0 ) 739: $ INFO = IERR 740: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 741: IF( SCALOC.NE.ONE ) THEN 742: DO 140 K = 1, N 743: CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 744: CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 745: 140 CONTINUE 746: SCALE = SCALE*SCALOC 747: END IF 748: * 749: * Unpack solution vector(s) 750: * 751: C( IS, JS ) = RHS( 1 ) 752: C( IS, JSP1 ) = RHS( 2 ) 753: F( IS, JS ) = RHS( 3 ) 754: F( IS, JSP1 ) = RHS( 4 ) 755: * 756: * Substitute R(I, J) and L(I, J) into remaining 757: * equation. 758: * 759: IF( J.GT.P+2 ) THEN 760: CALL DAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1, 761: $ F( IS, 1 ), LDF ) 762: CALL DAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1, 763: $ F( IS, 1 ), LDF ) 764: CALL DAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1, 765: $ F( IS, 1 ), LDF ) 766: CALL DAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1, 767: $ F( IS, 1 ), LDF ) 768: END IF 769: IF( I.LT.P ) THEN 770: CALL DGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA, 771: $ RHS( 1 ), 1, C( IE+1, JS ), LDC ) 772: CALL DGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD, 773: $ RHS( 3 ), 1, C( IE+1, JS ), LDC ) 774: END IF 775: * 776: ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN 777: * 778: * Build a 4-by-4 system Z' * x = RHS 779: * 780: Z( 1, 1 ) = A( IS, IS ) 781: Z( 2, 1 ) = A( IS, ISP1 ) 782: Z( 3, 1 ) = -B( JS, JS ) 783: Z( 4, 1 ) = ZERO 784: * 785: Z( 1, 2 ) = A( ISP1, IS ) 786: Z( 2, 2 ) = A( ISP1, ISP1 ) 787: Z( 3, 2 ) = ZERO 788: Z( 4, 2 ) = -B( JS, JS ) 789: * 790: Z( 1, 3 ) = D( IS, IS ) 791: Z( 2, 3 ) = D( IS, ISP1 ) 792: Z( 3, 3 ) = -E( JS, JS ) 793: Z( 4, 3 ) = ZERO 794: * 795: Z( 1, 4 ) = ZERO 796: Z( 2, 4 ) = D( ISP1, ISP1 ) 797: Z( 3, 4 ) = ZERO 798: Z( 4, 4 ) = -E( JS, JS ) 799: * 800: * Set up right hand side(s) 801: * 802: RHS( 1 ) = C( IS, JS ) 803: RHS( 2 ) = C( ISP1, JS ) 804: RHS( 3 ) = F( IS, JS ) 805: RHS( 4 ) = F( ISP1, JS ) 806: * 807: * Solve Z' * x = RHS 808: * 809: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 810: IF( IERR.GT.0 ) 811: $ INFO = IERR 812: * 813: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 814: IF( SCALOC.NE.ONE ) THEN 815: DO 150 K = 1, N 816: CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 817: CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 818: 150 CONTINUE 819: SCALE = SCALE*SCALOC 820: END IF 821: * 822: * Unpack solution vector(s) 823: * 824: C( IS, JS ) = RHS( 1 ) 825: C( ISP1, JS ) = RHS( 2 ) 826: F( IS, JS ) = RHS( 3 ) 827: F( ISP1, JS ) = RHS( 4 ) 828: * 829: * Substitute R(I, J) and L(I, J) into remaining 830: * equation. 831: * 832: IF( J.GT.P+2 ) THEN 833: CALL DGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ), 834: $ 1, F( IS, 1 ), LDF ) 835: CALL DGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ), 836: $ 1, F( IS, 1 ), LDF ) 837: END IF 838: IF( I.LT.P ) THEN 839: CALL DGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ), 840: $ LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ), 841: $ 1 ) 842: CALL DGEMV( 'T', MB, M-IE, -ONE, D( IS, IE+1 ), 843: $ LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ), 844: $ 1 ) 845: END IF 846: * 847: ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN 848: * 849: * Build an 8-by-8 system Z' * x = RHS 850: * 851: CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ ) 852: * 853: Z( 1, 1 ) = A( IS, IS ) 854: Z( 2, 1 ) = A( IS, ISP1 ) 855: Z( 5, 1 ) = -B( JS, JS ) 856: Z( 7, 1 ) = -B( JSP1, JS ) 857: * 858: Z( 1, 2 ) = A( ISP1, IS ) 859: Z( 2, 2 ) = A( ISP1, ISP1 ) 860: Z( 6, 2 ) = -B( JS, JS ) 861: Z( 8, 2 ) = -B( JSP1, JS ) 862: * 863: Z( 3, 3 ) = A( IS, IS ) 864: Z( 4, 3 ) = A( IS, ISP1 ) 865: Z( 5, 3 ) = -B( JS, JSP1 ) 866: Z( 7, 3 ) = -B( JSP1, JSP1 ) 867: * 868: Z( 3, 4 ) = A( ISP1, IS ) 869: Z( 4, 4 ) = A( ISP1, ISP1 ) 870: Z( 6, 4 ) = -B( JS, JSP1 ) 871: Z( 8, 4 ) = -B( JSP1, JSP1 ) 872: * 873: Z( 1, 5 ) = D( IS, IS ) 874: Z( 2, 5 ) = D( IS, ISP1 ) 875: Z( 5, 5 ) = -E( JS, JS ) 876: * 877: Z( 2, 6 ) = D( ISP1, ISP1 ) 878: Z( 6, 6 ) = -E( JS, JS ) 879: * 880: Z( 3, 7 ) = D( IS, IS ) 881: Z( 4, 7 ) = D( IS, ISP1 ) 882: Z( 5, 7 ) = -E( JS, JSP1 ) 883: Z( 7, 7 ) = -E( JSP1, JSP1 ) 884: * 885: Z( 4, 8 ) = D( ISP1, ISP1 ) 886: Z( 6, 8 ) = -E( JS, JSP1 ) 887: Z( 8, 8 ) = -E( JSP1, JSP1 ) 888: * 889: * Set up right hand side(s) 890: * 891: K = 1 892: II = MB*NB + 1 893: DO 160 JJ = 0, NB - 1 894: CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 ) 895: CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 ) 896: K = K + MB 897: II = II + MB 898: 160 CONTINUE 899: * 900: * 901: * Solve Z' * x = RHS 902: * 903: CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 904: IF( IERR.GT.0 ) 905: $ INFO = IERR 906: * 907: CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 908: IF( SCALOC.NE.ONE ) THEN 909: DO 170 K = 1, N 910: CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 911: CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 912: 170 CONTINUE 913: SCALE = SCALE*SCALOC 914: END IF 915: * 916: * Unpack solution vector(s) 917: * 918: K = 1 919: II = MB*NB + 1 920: DO 180 JJ = 0, NB - 1 921: CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 ) 922: CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 ) 923: K = K + MB 924: II = II + MB 925: 180 CONTINUE 926: * 927: * Substitute R(I, J) and L(I, J) into remaining 928: * equation. 929: * 930: IF( J.GT.P+2 ) THEN 931: CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE, 932: $ C( IS, JS ), LDC, B( 1, JS ), LDB, ONE, 933: $ F( IS, 1 ), LDF ) 934: CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE, 935: $ F( IS, JS ), LDF, E( 1, JS ), LDE, ONE, 936: $ F( IS, 1 ), LDF ) 937: END IF 938: IF( I.LT.P ) THEN 939: CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE, 940: $ A( IS, IE+1 ), LDA, C( IS, JS ), LDC, 941: $ ONE, C( IE+1, JS ), LDC ) 942: CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE, 943: $ D( IS, IE+1 ), LDD, F( IS, JS ), LDF, 944: $ ONE, C( IE+1, JS ), LDC ) 945: END IF 946: * 947: END IF 948: * 949: 190 CONTINUE 950: 200 CONTINUE 951: * 952: END IF 953: RETURN 954: * 955: * End of DTGSY2 956: * 957: END