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