![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DTPRFS( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, X, LDX, 2: $ FERR, BERR, WORK, IWORK, INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH. 10: * 11: * .. Scalar Arguments .. 12: CHARACTER DIAG, TRANS, UPLO 13: INTEGER INFO, LDB, LDX, N, NRHS 14: * .. 15: * .. Array Arguments .. 16: INTEGER IWORK( * ) 17: DOUBLE PRECISION AP( * ), B( LDB, * ), BERR( * ), FERR( * ), 18: $ WORK( * ), X( LDX, * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * DTPRFS provides error bounds and backward error estimates for the 25: * solution to a system of linear equations with a triangular packed 26: * coefficient matrix. 27: * 28: * The solution matrix X must be computed by DTPTRS or some other 29: * means before entering this routine. DTPRFS does not do iterative 30: * refinement because doing so cannot improve the backward error. 31: * 32: * Arguments 33: * ========= 34: * 35: * UPLO (input) CHARACTER*1 36: * = 'U': A is upper triangular; 37: * = 'L': A is lower triangular. 38: * 39: * TRANS (input) CHARACTER*1 40: * Specifies the form of the system of equations: 41: * = 'N': A * X = B (No transpose) 42: * = 'T': A**T * X = B (Transpose) 43: * = 'C': A**H * X = B (Conjugate transpose = Transpose) 44: * 45: * DIAG (input) CHARACTER*1 46: * = 'N': A is non-unit triangular; 47: * = 'U': A is unit triangular. 48: * 49: * N (input) INTEGER 50: * The order of the matrix A. N >= 0. 51: * 52: * NRHS (input) INTEGER 53: * The number of right hand sides, i.e., the number of columns 54: * of the matrices B and X. NRHS >= 0. 55: * 56: * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2) 57: * The upper or lower triangular matrix A, packed columnwise in 58: * a linear array. The j-th column of A is stored in the array 59: * AP as follows: 60: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 61: * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 62: * If DIAG = 'U', the diagonal elements of A are not referenced 63: * and are assumed to be 1. 64: * 65: * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) 66: * The right hand side matrix B. 67: * 68: * LDB (input) INTEGER 69: * The leading dimension of the array B. LDB >= max(1,N). 70: * 71: * X (input) DOUBLE PRECISION array, dimension (LDX,NRHS) 72: * The solution matrix X. 73: * 74: * LDX (input) INTEGER 75: * The leading dimension of the array X. LDX >= max(1,N). 76: * 77: * FERR (output) DOUBLE PRECISION array, dimension (NRHS) 78: * The estimated forward error bound for each solution vector 79: * X(j) (the j-th column of the solution matrix X). 80: * If XTRUE is the true solution corresponding to X(j), FERR(j) 81: * is an estimated upper bound for the magnitude of the largest 82: * element in (X(j) - XTRUE) divided by the magnitude of the 83: * largest element in X(j). The estimate is as reliable as 84: * the estimate for RCOND, and is almost always a slight 85: * overestimate of the true error. 86: * 87: * BERR (output) DOUBLE PRECISION array, dimension (NRHS) 88: * The componentwise relative backward error of each solution 89: * vector X(j) (i.e., the smallest relative change in 90: * any element of A or B that makes X(j) an exact solution). 91: * 92: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) 93: * 94: * IWORK (workspace) INTEGER array, dimension (N) 95: * 96: * INFO (output) INTEGER 97: * = 0: successful exit 98: * < 0: if INFO = -i, the i-th argument had an illegal value 99: * 100: * ===================================================================== 101: * 102: * .. Parameters .. 103: DOUBLE PRECISION ZERO 104: PARAMETER ( ZERO = 0.0D+0 ) 105: DOUBLE PRECISION ONE 106: PARAMETER ( ONE = 1.0D+0 ) 107: * .. 108: * .. Local Scalars .. 109: LOGICAL NOTRAN, NOUNIT, UPPER 110: CHARACTER TRANST 111: INTEGER I, J, K, KASE, KC, NZ 112: DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK 113: * .. 114: * .. Local Arrays .. 115: INTEGER ISAVE( 3 ) 116: * .. 117: * .. External Subroutines .. 118: EXTERNAL DAXPY, DCOPY, DLACN2, DTPMV, DTPSV, XERBLA 119: * .. 120: * .. Intrinsic Functions .. 121: INTRINSIC ABS, MAX 122: * .. 123: * .. External Functions .. 124: LOGICAL LSAME 125: DOUBLE PRECISION DLAMCH 126: EXTERNAL LSAME, DLAMCH 127: * .. 128: * .. Executable Statements .. 129: * 130: * Test the input parameters. 131: * 132: INFO = 0 133: UPPER = LSAME( UPLO, 'U' ) 134: NOTRAN = LSAME( TRANS, 'N' ) 135: NOUNIT = LSAME( DIAG, 'N' ) 136: * 137: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 138: INFO = -1 139: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 140: $ LSAME( TRANS, 'C' ) ) THEN 141: INFO = -2 142: ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 143: INFO = -3 144: ELSE IF( N.LT.0 ) THEN 145: INFO = -4 146: ELSE IF( NRHS.LT.0 ) THEN 147: INFO = -5 148: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 149: INFO = -8 150: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 151: INFO = -10 152: END IF 153: IF( INFO.NE.0 ) THEN 154: CALL XERBLA( 'DTPRFS', -INFO ) 155: RETURN 156: END IF 157: * 158: * Quick return if possible 159: * 160: IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 161: DO 10 J = 1, NRHS 162: FERR( J ) = ZERO 163: BERR( J ) = ZERO 164: 10 CONTINUE 165: RETURN 166: END IF 167: * 168: IF( NOTRAN ) THEN 169: TRANST = 'T' 170: ELSE 171: TRANST = 'N' 172: END IF 173: * 174: * NZ = maximum number of nonzero elements in each row of A, plus 1 175: * 176: NZ = N + 1 177: EPS = DLAMCH( 'Epsilon' ) 178: SAFMIN = DLAMCH( 'Safe minimum' ) 179: SAFE1 = NZ*SAFMIN 180: SAFE2 = SAFE1 / EPS 181: * 182: * Do for each right hand side 183: * 184: DO 250 J = 1, NRHS 185: * 186: * Compute residual R = B - op(A) * X, 187: * where op(A) = A or A', depending on TRANS. 188: * 189: CALL DCOPY( N, X( 1, J ), 1, WORK( N+1 ), 1 ) 190: CALL DTPMV( UPLO, TRANS, DIAG, N, AP, WORK( N+1 ), 1 ) 191: CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK( N+1 ), 1 ) 192: * 193: * Compute componentwise relative backward error from formula 194: * 195: * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) ) 196: * 197: * where abs(Z) is the componentwise absolute value of the matrix 198: * or vector Z. If the i-th component of the denominator is less 199: * than SAFE2, then SAFE1 is added to the i-th components of the 200: * numerator and denominator before dividing. 201: * 202: DO 20 I = 1, N 203: WORK( I ) = ABS( B( I, J ) ) 204: 20 CONTINUE 205: * 206: IF( NOTRAN ) THEN 207: * 208: * Compute abs(A)*abs(X) + abs(B). 209: * 210: IF( UPPER ) THEN 211: KC = 1 212: IF( NOUNIT ) THEN 213: DO 40 K = 1, N 214: XK = ABS( X( K, J ) ) 215: DO 30 I = 1, K 216: WORK( I ) = WORK( I ) + ABS( AP( KC+I-1 ) )*XK 217: 30 CONTINUE 218: KC = KC + K 219: 40 CONTINUE 220: ELSE 221: DO 60 K = 1, N 222: XK = ABS( X( K, J ) ) 223: DO 50 I = 1, K - 1 224: WORK( I ) = WORK( I ) + ABS( AP( KC+I-1 ) )*XK 225: 50 CONTINUE 226: WORK( K ) = WORK( K ) + XK 227: KC = KC + K 228: 60 CONTINUE 229: END IF 230: ELSE 231: KC = 1 232: IF( NOUNIT ) THEN 233: DO 80 K = 1, N 234: XK = ABS( X( K, J ) ) 235: DO 70 I = K, N 236: WORK( I ) = WORK( I ) + ABS( AP( KC+I-K ) )*XK 237: 70 CONTINUE 238: KC = KC + N - K + 1 239: 80 CONTINUE 240: ELSE 241: DO 100 K = 1, N 242: XK = ABS( X( K, J ) ) 243: DO 90 I = K + 1, N 244: WORK( I ) = WORK( I ) + ABS( AP( KC+I-K ) )*XK 245: 90 CONTINUE 246: WORK( K ) = WORK( K ) + XK 247: KC = KC + N - K + 1 248: 100 CONTINUE 249: END IF 250: END IF 251: ELSE 252: * 253: * Compute abs(A')*abs(X) + abs(B). 254: * 255: IF( UPPER ) THEN 256: KC = 1 257: IF( NOUNIT ) THEN 258: DO 120 K = 1, N 259: S = ZERO 260: DO 110 I = 1, K 261: S = S + ABS( AP( KC+I-1 ) )*ABS( X( I, J ) ) 262: 110 CONTINUE 263: WORK( K ) = WORK( K ) + S 264: KC = KC + K 265: 120 CONTINUE 266: ELSE 267: DO 140 K = 1, N 268: S = ABS( X( K, J ) ) 269: DO 130 I = 1, K - 1 270: S = S + ABS( AP( KC+I-1 ) )*ABS( X( I, J ) ) 271: 130 CONTINUE 272: WORK( K ) = WORK( K ) + S 273: KC = KC + K 274: 140 CONTINUE 275: END IF 276: ELSE 277: KC = 1 278: IF( NOUNIT ) THEN 279: DO 160 K = 1, N 280: S = ZERO 281: DO 150 I = K, N 282: S = S + ABS( AP( KC+I-K ) )*ABS( X( I, J ) ) 283: 150 CONTINUE 284: WORK( K ) = WORK( K ) + S 285: KC = KC + N - K + 1 286: 160 CONTINUE 287: ELSE 288: DO 180 K = 1, N 289: S = ABS( X( K, J ) ) 290: DO 170 I = K + 1, N 291: S = S + ABS( AP( KC+I-K ) )*ABS( X( I, J ) ) 292: 170 CONTINUE 293: WORK( K ) = WORK( K ) + S 294: KC = KC + N - K + 1 295: 180 CONTINUE 296: END IF 297: END IF 298: END IF 299: S = ZERO 300: DO 190 I = 1, N 301: IF( WORK( I ).GT.SAFE2 ) THEN 302: S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) ) 303: ELSE 304: S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) / 305: $ ( WORK( I )+SAFE1 ) ) 306: END IF 307: 190 CONTINUE 308: BERR( J ) = S 309: * 310: * Bound error from formula 311: * 312: * norm(X - XTRUE) / norm(X) .le. FERR = 313: * norm( abs(inv(op(A)))* 314: * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X) 315: * 316: * where 317: * norm(Z) is the magnitude of the largest component of Z 318: * inv(op(A)) is the inverse of op(A) 319: * abs(Z) is the componentwise absolute value of the matrix or 320: * vector Z 321: * NZ is the maximum number of nonzeros in any row of A, plus 1 322: * EPS is machine epsilon 323: * 324: * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B)) 325: * is incremented by SAFE1 if the i-th component of 326: * abs(op(A))*abs(X) + abs(B) is less than SAFE2. 327: * 328: * Use DLACN2 to estimate the infinity-norm of the matrix 329: * inv(op(A)) * diag(W), 330: * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) 331: * 332: DO 200 I = 1, N 333: IF( WORK( I ).GT.SAFE2 ) THEN 334: WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) 335: ELSE 336: WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1 337: END IF 338: 200 CONTINUE 339: * 340: KASE = 0 341: 210 CONTINUE 342: CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ), 343: $ KASE, ISAVE ) 344: IF( KASE.NE.0 ) THEN 345: IF( KASE.EQ.1 ) THEN 346: * 347: * Multiply by diag(W)*inv(op(A)'). 348: * 349: CALL DTPSV( UPLO, TRANST, DIAG, N, AP, WORK( N+1 ), 1 ) 350: DO 220 I = 1, N 351: WORK( N+I ) = WORK( I )*WORK( N+I ) 352: 220 CONTINUE 353: ELSE 354: * 355: * Multiply by inv(op(A))*diag(W). 356: * 357: DO 230 I = 1, N 358: WORK( N+I ) = WORK( I )*WORK( N+I ) 359: 230 CONTINUE 360: CALL DTPSV( UPLO, TRANS, DIAG, N, AP, WORK( N+1 ), 1 ) 361: END IF 362: GO TO 210 363: END IF 364: * 365: * Normalize error. 366: * 367: LSTRES = ZERO 368: DO 240 I = 1, N 369: LSTRES = MAX( LSTRES, ABS( X( I, J ) ) ) 370: 240 CONTINUE 371: IF( LSTRES.NE.ZERO ) 372: $ FERR( J ) = FERR( J ) / LSTRES 373: * 374: 250 CONTINUE 375: * 376: RETURN 377: * 378: * End of DTPRFS 379: * 380: END