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