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