![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO ) 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 INFO, JOB, N 10: DOUBLE PRECISION TOL 11: * .. 12: * .. Array Arguments .. 13: INTEGER IN( * ) 14: DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DLAGTS may be used to solve one of the systems of equations 21: * 22: * (T - lambda*I)*x = y or (T - lambda*I)'*x = y, 23: * 24: * where T is an n by n tridiagonal matrix, for x, following the 25: * factorization of (T - lambda*I) as 26: * 27: * (T - lambda*I) = P*L*U , 28: * 29: * by routine DLAGTF. The choice of equation to be solved is 30: * controlled by the argument JOB, and in each case there is an option 31: * to perturb zero or very small diagonal elements of U, this option 32: * being intended for use in applications such as inverse iteration. 33: * 34: * Arguments 35: * ========= 36: * 37: * JOB (input) INTEGER 38: * Specifies the job to be performed by DLAGTS as follows: 39: * = 1: The equations (T - lambda*I)x = y are to be solved, 40: * but diagonal elements of U are not to be perturbed. 41: * = -1: The equations (T - lambda*I)x = y are to be solved 42: * and, if overflow would otherwise occur, the diagonal 43: * elements of U are to be perturbed. See argument TOL 44: * below. 45: * = 2: The equations (T - lambda*I)'x = y are to be solved, 46: * but diagonal elements of U are not to be perturbed. 47: * = -2: The equations (T - lambda*I)'x = y are to be solved 48: * and, if overflow would otherwise occur, the diagonal 49: * elements of U are to be perturbed. See argument TOL 50: * below. 51: * 52: * N (input) INTEGER 53: * The order of the matrix T. 54: * 55: * A (input) DOUBLE PRECISION array, dimension (N) 56: * On entry, A must contain the diagonal elements of U as 57: * returned from DLAGTF. 58: * 59: * B (input) DOUBLE PRECISION array, dimension (N-1) 60: * On entry, B must contain the first super-diagonal elements of 61: * U as returned from DLAGTF. 62: * 63: * C (input) DOUBLE PRECISION array, dimension (N-1) 64: * On entry, C must contain the sub-diagonal elements of L as 65: * returned from DLAGTF. 66: * 67: * D (input) DOUBLE PRECISION array, dimension (N-2) 68: * On entry, D must contain the second super-diagonal elements 69: * of U as returned from DLAGTF. 70: * 71: * IN (input) INTEGER array, dimension (N) 72: * On entry, IN must contain details of the matrix P as returned 73: * from DLAGTF. 74: * 75: * Y (input/output) DOUBLE PRECISION array, dimension (N) 76: * On entry, the right hand side vector y. 77: * On exit, Y is overwritten by the solution vector x. 78: * 79: * TOL (input/output) DOUBLE PRECISION 80: * On entry, with JOB .lt. 0, TOL should be the minimum 81: * perturbation to be made to very small diagonal elements of U. 82: * TOL should normally be chosen as about eps*norm(U), where eps 83: * is the relative machine precision, but if TOL is supplied as 84: * non-positive, then it is reset to eps*max( abs( u(i,j) ) ). 85: * If JOB .gt. 0 then TOL is not referenced. 86: * 87: * On exit, TOL is changed as described above, only if TOL is 88: * non-positive on entry. Otherwise TOL is unchanged. 89: * 90: * INFO (output) INTEGER 91: * = 0 : successful exit 92: * .lt. 0: if INFO = -i, the i-th argument had an illegal value 93: * .gt. 0: overflow would occur when computing the INFO(th) 94: * element of the solution vector x. This can only occur 95: * when JOB is supplied as positive and either means 96: * that a diagonal element of U is very small, or that 97: * the elements of the right-hand side vector y are very 98: * large. 99: * 100: * ===================================================================== 101: * 102: * .. Parameters .. 103: DOUBLE PRECISION ONE, ZERO 104: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 105: * .. 106: * .. Local Scalars .. 107: INTEGER K 108: DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP 109: * .. 110: * .. Intrinsic Functions .. 111: INTRINSIC ABS, MAX, SIGN 112: * .. 113: * .. External Functions .. 114: DOUBLE PRECISION DLAMCH 115: EXTERNAL DLAMCH 116: * .. 117: * .. External Subroutines .. 118: EXTERNAL XERBLA 119: * .. 120: * .. Executable Statements .. 121: * 122: INFO = 0 123: IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN 124: INFO = -1 125: ELSE IF( N.LT.0 ) THEN 126: INFO = -2 127: END IF 128: IF( INFO.NE.0 ) THEN 129: CALL XERBLA( 'DLAGTS', -INFO ) 130: RETURN 131: END IF 132: * 133: IF( N.EQ.0 ) 134: $ RETURN 135: * 136: EPS = DLAMCH( 'Epsilon' ) 137: SFMIN = DLAMCH( 'Safe minimum' ) 138: BIGNUM = ONE / SFMIN 139: * 140: IF( JOB.LT.0 ) THEN 141: IF( TOL.LE.ZERO ) THEN 142: TOL = ABS( A( 1 ) ) 143: IF( N.GT.1 ) 144: $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) ) 145: DO 10 K = 3, N 146: TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ), 147: $ ABS( D( K-2 ) ) ) 148: 10 CONTINUE 149: TOL = TOL*EPS 150: IF( TOL.EQ.ZERO ) 151: $ TOL = EPS 152: END IF 153: END IF 154: * 155: IF( ABS( JOB ).EQ.1 ) THEN 156: DO 20 K = 2, N 157: IF( IN( K-1 ).EQ.0 ) THEN 158: Y( K ) = Y( K ) - C( K-1 )*Y( K-1 ) 159: ELSE 160: TEMP = Y( K-1 ) 161: Y( K-1 ) = Y( K ) 162: Y( K ) = TEMP - C( K-1 )*Y( K ) 163: END IF 164: 20 CONTINUE 165: IF( JOB.EQ.1 ) THEN 166: DO 30 K = N, 1, -1 167: IF( K.LE.N-2 ) THEN 168: TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) 169: ELSE IF( K.EQ.N-1 ) THEN 170: TEMP = Y( K ) - B( K )*Y( K+1 ) 171: ELSE 172: TEMP = Y( K ) 173: END IF 174: AK = A( K ) 175: ABSAK = ABS( AK ) 176: IF( ABSAK.LT.ONE ) THEN 177: IF( ABSAK.LT.SFMIN ) THEN 178: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 179: $ THEN 180: INFO = K 181: RETURN 182: ELSE 183: TEMP = TEMP*BIGNUM 184: AK = AK*BIGNUM 185: END IF 186: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 187: INFO = K 188: RETURN 189: END IF 190: END IF 191: Y( K ) = TEMP / AK 192: 30 CONTINUE 193: ELSE 194: DO 50 K = N, 1, -1 195: IF( K.LE.N-2 ) THEN 196: TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) 197: ELSE IF( K.EQ.N-1 ) THEN 198: TEMP = Y( K ) - B( K )*Y( K+1 ) 199: ELSE 200: TEMP = Y( K ) 201: END IF 202: AK = A( K ) 203: PERT = SIGN( TOL, AK ) 204: 40 CONTINUE 205: ABSAK = ABS( AK ) 206: IF( ABSAK.LT.ONE ) THEN 207: IF( ABSAK.LT.SFMIN ) THEN 208: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 209: $ THEN 210: AK = AK + PERT 211: PERT = 2*PERT 212: GO TO 40 213: ELSE 214: TEMP = TEMP*BIGNUM 215: AK = AK*BIGNUM 216: END IF 217: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 218: AK = AK + PERT 219: PERT = 2*PERT 220: GO TO 40 221: END IF 222: END IF 223: Y( K ) = TEMP / AK 224: 50 CONTINUE 225: END IF 226: ELSE 227: * 228: * Come to here if JOB = 2 or -2 229: * 230: IF( JOB.EQ.2 ) THEN 231: DO 60 K = 1, N 232: IF( K.GE.3 ) THEN 233: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) 234: ELSE IF( K.EQ.2 ) THEN 235: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) 236: ELSE 237: TEMP = Y( K ) 238: END IF 239: AK = A( K ) 240: ABSAK = ABS( AK ) 241: IF( ABSAK.LT.ONE ) THEN 242: IF( ABSAK.LT.SFMIN ) THEN 243: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 244: $ THEN 245: INFO = K 246: RETURN 247: ELSE 248: TEMP = TEMP*BIGNUM 249: AK = AK*BIGNUM 250: END IF 251: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 252: INFO = K 253: RETURN 254: END IF 255: END IF 256: Y( K ) = TEMP / AK 257: 60 CONTINUE 258: ELSE 259: DO 80 K = 1, N 260: IF( K.GE.3 ) THEN 261: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) 262: ELSE IF( K.EQ.2 ) THEN 263: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) 264: ELSE 265: TEMP = Y( K ) 266: END IF 267: AK = A( K ) 268: PERT = SIGN( TOL, AK ) 269: 70 CONTINUE 270: ABSAK = ABS( AK ) 271: IF( ABSAK.LT.ONE ) THEN 272: IF( ABSAK.LT.SFMIN ) THEN 273: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 274: $ THEN 275: AK = AK + PERT 276: PERT = 2*PERT 277: GO TO 70 278: ELSE 279: TEMP = TEMP*BIGNUM 280: AK = AK*BIGNUM 281: END IF 282: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 283: AK = AK + PERT 284: PERT = 2*PERT 285: GO TO 70 286: END IF 287: END IF 288: Y( K ) = TEMP / AK 289: 80 CONTINUE 290: END IF 291: * 292: DO 90 K = N, 2, -1 293: IF( IN( K-1 ).EQ.0 ) THEN 294: Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K ) 295: ELSE 296: TEMP = Y( K-1 ) 297: Y( K-1 ) = Y( K ) 298: Y( K ) = TEMP - C( K-1 )*Y( K ) 299: END IF 300: 90 CONTINUE 301: END IF 302: * 303: * End of DLAGTS 304: * 305: END