![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB ) 2: * 3: * -- LAPACK auxiliary routine (version 3.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * November 2006 7: * 8: * .. Scalar Arguments .. 9: INTEGER ITRANS, LDB, N, NRHS 10: * .. 11: * .. Array Arguments .. 12: INTEGER IPIV( * ) 13: DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * DGTTS2 solves one of the systems of equations 20: * A*X = B or A'*X = B, 21: * with a tridiagonal matrix A using the LU factorization computed 22: * by DGTTRF. 23: * 24: * Arguments 25: * ========= 26: * 27: * ITRANS (input) INTEGER 28: * Specifies the form of the system of equations. 29: * = 0: A * X = B (No transpose) 30: * = 1: A'* X = B (Transpose) 31: * = 2: A'* X = B (Conjugate transpose = Transpose) 32: * 33: * N (input) INTEGER 34: * The order of the matrix A. 35: * 36: * NRHS (input) INTEGER 37: * The number of right hand sides, i.e., the number of columns 38: * of the matrix B. NRHS >= 0. 39: * 40: * DL (input) DOUBLE PRECISION array, dimension (N-1) 41: * The (n-1) multipliers that define the matrix L from the 42: * LU factorization of A. 43: * 44: * D (input) DOUBLE PRECISION array, dimension (N) 45: * The n diagonal elements of the upper triangular matrix U from 46: * the LU factorization of A. 47: * 48: * DU (input) DOUBLE PRECISION array, dimension (N-1) 49: * The (n-1) elements of the first super-diagonal of U. 50: * 51: * DU2 (input) DOUBLE PRECISION array, dimension (N-2) 52: * The (n-2) elements of the second super-diagonal of U. 53: * 54: * IPIV (input) INTEGER array, dimension (N) 55: * The pivot indices; for 1 <= i <= n, row i of the matrix was 56: * interchanged with row IPIV(i). IPIV(i) will always be either 57: * i or i+1; IPIV(i) = i indicates a row interchange was not 58: * required. 59: * 60: * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 61: * On entry, the matrix of right hand side vectors B. 62: * On exit, B is overwritten by the solution vectors X. 63: * 64: * LDB (input) INTEGER 65: * The leading dimension of the array B. LDB >= max(1,N). 66: * 67: * ===================================================================== 68: * 69: * .. Local Scalars .. 70: INTEGER I, IP, J 71: DOUBLE PRECISION TEMP 72: * .. 73: * .. Executable Statements .. 74: * 75: * Quick return if possible 76: * 77: IF( N.EQ.0 .OR. NRHS.EQ.0 ) 78: $ RETURN 79: * 80: IF( ITRANS.EQ.0 ) THEN 81: * 82: * Solve A*X = B using the LU factorization of A, 83: * overwriting each right hand side vector with its solution. 84: * 85: IF( NRHS.LE.1 ) THEN 86: J = 1 87: 10 CONTINUE 88: * 89: * Solve L*x = b. 90: * 91: DO 20 I = 1, N - 1 92: IP = IPIV( I ) 93: TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J ) 94: B( I, J ) = B( IP, J ) 95: B( I+1, J ) = TEMP 96: 20 CONTINUE 97: * 98: * Solve U*x = b. 99: * 100: B( N, J ) = B( N, J ) / D( N ) 101: IF( N.GT.1 ) 102: $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / 103: $ D( N-1 ) 104: DO 30 I = N - 2, 1, -1 105: B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* 106: $ B( I+2, J ) ) / D( I ) 107: 30 CONTINUE 108: IF( J.LT.NRHS ) THEN 109: J = J + 1 110: GO TO 10 111: END IF 112: ELSE 113: DO 60 J = 1, NRHS 114: * 115: * Solve L*x = b. 116: * 117: DO 40 I = 1, N - 1 118: IF( IPIV( I ).EQ.I ) THEN 119: B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J ) 120: ELSE 121: TEMP = B( I, J ) 122: B( I, J ) = B( I+1, J ) 123: B( I+1, J ) = TEMP - DL( I )*B( I, J ) 124: END IF 125: 40 CONTINUE 126: * 127: * Solve U*x = b. 128: * 129: B( N, J ) = B( N, J ) / D( N ) 130: IF( N.GT.1 ) 131: $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / 132: $ D( N-1 ) 133: DO 50 I = N - 2, 1, -1 134: B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* 135: $ B( I+2, J ) ) / D( I ) 136: 50 CONTINUE 137: 60 CONTINUE 138: END IF 139: ELSE 140: * 141: * Solve A' * X = B. 142: * 143: IF( NRHS.LE.1 ) THEN 144: * 145: * Solve U'*x = b. 146: * 147: J = 1 148: 70 CONTINUE 149: B( 1, J ) = B( 1, J ) / D( 1 ) 150: IF( N.GT.1 ) 151: $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) 152: DO 80 I = 3, N 153: B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )* 154: $ B( I-2, J ) ) / D( I ) 155: 80 CONTINUE 156: * 157: * Solve L'*x = b. 158: * 159: DO 90 I = N - 1, 1, -1 160: IP = IPIV( I ) 161: TEMP = B( I, J ) - DL( I )*B( I+1, J ) 162: B( I, J ) = B( IP, J ) 163: B( IP, J ) = TEMP 164: 90 CONTINUE 165: IF( J.LT.NRHS ) THEN 166: J = J + 1 167: GO TO 70 168: END IF 169: * 170: ELSE 171: DO 120 J = 1, NRHS 172: * 173: * Solve U'*x = b. 174: * 175: B( 1, J ) = B( 1, J ) / D( 1 ) 176: IF( N.GT.1 ) 177: $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) 178: DO 100 I = 3, N 179: B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )- 180: $ DU2( I-2 )*B( I-2, J ) ) / D( I ) 181: 100 CONTINUE 182: DO 110 I = N - 1, 1, -1 183: IF( IPIV( I ).EQ.I ) THEN 184: B( I, J ) = B( I, J ) - DL( I )*B( I+1, J ) 185: ELSE 186: TEMP = B( I+1, J ) 187: B( I+1, J ) = B( I, J ) - DL( I )*TEMP 188: B( I, J ) = TEMP 189: END IF 190: 110 CONTINUE 191: 120 CONTINUE 192: END IF 193: END IF 194: * 195: * End of DGTTS2 196: * 197: END