001: SUBROUTINE CGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, 002: $ IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, 003: $ INFO ) 004: * 005: * -- LAPACK routine (version 3.2) -- 006: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 007: * November 2006 008: * 009: * Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH. 010: * 011: * .. Scalar Arguments .. 012: CHARACTER TRANS 013: INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS 014: * .. 015: * .. Array Arguments .. 016: INTEGER IPIV( * ) 017: REAL BERR( * ), FERR( * ), RWORK( * ) 018: COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 019: $ WORK( * ), X( LDX, * ) 020: * .. 021: * 022: * Purpose 023: * ======= 024: * 025: * CGBRFS improves the computed solution to a system of linear 026: * equations when the coefficient matrix is banded, and provides 027: * error bounds and backward error estimates for the solution. 028: * 029: * Arguments 030: * ========= 031: * 032: * TRANS (input) CHARACTER*1 033: * Specifies the form of the system of equations: 034: * = 'N': A * X = B (No transpose) 035: * = 'T': A**T * X = B (Transpose) 036: * = 'C': A**H * X = B (Conjugate transpose) 037: * 038: * N (input) INTEGER 039: * The order of the matrix A. N >= 0. 040: * 041: * KL (input) INTEGER 042: * The number of subdiagonals within the band of A. KL >= 0. 043: * 044: * KU (input) INTEGER 045: * The number of superdiagonals within the band of A. KU >= 0. 046: * 047: * NRHS (input) INTEGER 048: * The number of right hand sides, i.e., the number of columns 049: * of the matrices B and X. NRHS >= 0. 050: * 051: * AB (input) COMPLEX array, dimension (LDAB,N) 052: * The original band matrix A, stored in rows 1 to KL+KU+1. 053: * The j-th column of A is stored in the j-th column of the 054: * array AB as follows: 055: * AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl). 056: * 057: * LDAB (input) INTEGER 058: * The leading dimension of the array AB. LDAB >= KL+KU+1. 059: * 060: * AFB (input) COMPLEX array, dimension (LDAFB,N) 061: * Details of the LU factorization of the band matrix A, as 062: * computed by CGBTRF. U is stored as an upper triangular band 063: * matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and 064: * the multipliers used during the factorization are stored in 065: * rows KL+KU+2 to 2*KL+KU+1. 066: * 067: * LDAFB (input) INTEGER 068: * The leading dimension of the array AFB. LDAFB >= 2*KL*KU+1. 069: * 070: * IPIV (input) INTEGER array, dimension (N) 071: * The pivot indices from CGBTRF; for 1<=i<=N, row i of the 072: * matrix was interchanged with row IPIV(i). 073: * 074: * B (input) COMPLEX array, dimension (LDB,NRHS) 075: * The right hand side matrix B. 076: * 077: * LDB (input) INTEGER 078: * The leading dimension of the array B. LDB >= max(1,N). 079: * 080: * X (input/output) COMPLEX array, dimension (LDX,NRHS) 081: * On entry, the solution matrix X, as computed by CGBTRS. 082: * On exit, the improved solution matrix X. 083: * 084: * LDX (input) INTEGER 085: * The leading dimension of the array X. LDX >= max(1,N). 086: * 087: * FERR (output) REAL array, dimension (NRHS) 088: * The estimated forward error bound for each solution vector 089: * X(j) (the j-th column of the solution matrix X). 090: * If XTRUE is the true solution corresponding to X(j), FERR(j) 091: * is an estimated upper bound for the magnitude of the largest 092: * element in (X(j) - XTRUE) divided by the magnitude of the 093: * largest element in X(j). The estimate is as reliable as 094: * the estimate for RCOND, and is almost always a slight 095: * overestimate of the true error. 096: * 097: * BERR (output) REAL array, dimension (NRHS) 098: * The componentwise relative backward error of each solution 099: * vector X(j) (i.e., the smallest relative change in 100: * any element of A or B that makes X(j) an exact solution). 101: * 102: * WORK (workspace) COMPLEX array, dimension (2*N) 103: * 104: * RWORK (workspace) REAL array, dimension (N) 105: * 106: * INFO (output) INTEGER 107: * = 0: successful exit 108: * < 0: if INFO = -i, the i-th argument had an illegal value 109: * 110: * Internal Parameters 111: * =================== 112: * 113: * ITMAX is the maximum number of steps of iterative refinement. 114: * 115: * ===================================================================== 116: * 117: * .. Parameters .. 118: INTEGER ITMAX 119: PARAMETER ( ITMAX = 5 ) 120: REAL ZERO 121: PARAMETER ( ZERO = 0.0E+0 ) 122: COMPLEX CONE 123: PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 124: REAL TWO 125: PARAMETER ( TWO = 2.0E+0 ) 126: REAL THREE 127: PARAMETER ( THREE = 3.0E+0 ) 128: * .. 129: * .. Local Scalars .. 130: LOGICAL NOTRAN 131: CHARACTER TRANSN, TRANST 132: INTEGER COUNT, I, J, K, KASE, KK, NZ 133: REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK 134: COMPLEX ZDUM 135: * .. 136: * .. Local Arrays .. 137: INTEGER ISAVE( 3 ) 138: * .. 139: * .. External Subroutines .. 140: EXTERNAL CAXPY, CCOPY, CGBMV, CGBTRS, CLACN2, XERBLA 141: * .. 142: * .. Intrinsic Functions .. 143: INTRINSIC ABS, AIMAG, MAX, MIN, REAL 144: * .. 145: * .. External Functions .. 146: LOGICAL LSAME 147: REAL SLAMCH 148: EXTERNAL LSAME, SLAMCH 149: * .. 150: * .. Statement Functions .. 151: REAL CABS1 152: * .. 153: * .. Statement Function definitions .. 154: CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 155: * .. 156: * .. Executable Statements .. 157: * 158: * Test the input parameters. 159: * 160: INFO = 0 161: NOTRAN = LSAME( TRANS, 'N' ) 162: IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 163: $ LSAME( TRANS, 'C' ) ) THEN 164: INFO = -1 165: ELSE IF( N.LT.0 ) THEN 166: INFO = -2 167: ELSE IF( KL.LT.0 ) THEN 168: INFO = -3 169: ELSE IF( KU.LT.0 ) THEN 170: INFO = -4 171: ELSE IF( NRHS.LT.0 ) THEN 172: INFO = -5 173: ELSE IF( LDAB.LT.KL+KU+1 ) THEN 174: INFO = -7 175: ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 176: INFO = -9 177: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 178: INFO = -12 179: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 180: INFO = -14 181: END IF 182: IF( INFO.NE.0 ) THEN 183: CALL XERBLA( 'CGBRFS', -INFO ) 184: RETURN 185: END IF 186: * 187: * Quick return if possible 188: * 189: IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 190: DO 10 J = 1, NRHS 191: FERR( J ) = ZERO 192: BERR( J ) = ZERO 193: 10 CONTINUE 194: RETURN 195: END IF 196: * 197: IF( NOTRAN ) THEN 198: TRANSN = 'N' 199: TRANST = 'C' 200: ELSE 201: TRANSN = 'C' 202: TRANST = 'N' 203: END IF 204: * 205: * NZ = maximum number of nonzero elements in each row of A, plus 1 206: * 207: NZ = MIN( KL+KU+2, N+1 ) 208: EPS = SLAMCH( 'Epsilon' ) 209: SAFMIN = SLAMCH( 'Safe minimum' ) 210: SAFE1 = NZ*SAFMIN 211: SAFE2 = SAFE1 / EPS 212: * 213: * Do for each right hand side 214: * 215: DO 140 J = 1, NRHS 216: * 217: COUNT = 1 218: LSTRES = THREE 219: 20 CONTINUE 220: * 221: * Loop until stopping criterion is satisfied. 222: * 223: * Compute residual R = B - op(A) * X, 224: * where op(A) = A, A**T, or A**H, depending on TRANS. 225: * 226: CALL CCOPY( N, B( 1, J ), 1, WORK, 1 ) 227: CALL CGBMV( TRANS, N, N, KL, KU, -CONE, AB, LDAB, X( 1, J ), 1, 228: $ CONE, WORK, 1 ) 229: * 230: * Compute componentwise relative backward error from formula 231: * 232: * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) ) 233: * 234: * where abs(Z) is the componentwise absolute value of the matrix 235: * or vector Z. If the i-th component of the denominator is less 236: * than SAFE2, then SAFE1 is added to the i-th components of the 237: * numerator and denominator before dividing. 238: * 239: DO 30 I = 1, N 240: RWORK( I ) = CABS1( B( I, J ) ) 241: 30 CONTINUE 242: * 243: * Compute abs(op(A))*abs(X) + abs(B). 244: * 245: IF( NOTRAN ) THEN 246: DO 50 K = 1, N 247: KK = KU + 1 - K 248: XK = CABS1( X( K, J ) ) 249: DO 40 I = MAX( 1, K-KU ), MIN( N, K+KL ) 250: RWORK( I ) = RWORK( I ) + CABS1( AB( KK+I, K ) )*XK 251: 40 CONTINUE 252: 50 CONTINUE 253: ELSE 254: DO 70 K = 1, N 255: S = ZERO 256: KK = KU + 1 - K 257: DO 60 I = MAX( 1, K-KU ), MIN( N, K+KL ) 258: S = S + CABS1( AB( KK+I, K ) )*CABS1( X( I, J ) ) 259: 60 CONTINUE 260: RWORK( K ) = RWORK( K ) + S 261: 70 CONTINUE 262: END IF 263: S = ZERO 264: DO 80 I = 1, N 265: IF( RWORK( I ).GT.SAFE2 ) THEN 266: S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) ) 267: ELSE 268: S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) / 269: $ ( RWORK( I )+SAFE1 ) ) 270: END IF 271: 80 CONTINUE 272: BERR( J ) = S 273: * 274: * Test stopping criterion. Continue iterating if 275: * 1) The residual BERR(J) is larger than machine epsilon, and 276: * 2) BERR(J) decreased by at least a factor of 2 during the 277: * last iteration, and 278: * 3) At most ITMAX iterations tried. 279: * 280: IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND. 281: $ COUNT.LE.ITMAX ) THEN 282: * 283: * Update solution and try again. 284: * 285: CALL CGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV, WORK, N, 286: $ INFO ) 287: CALL CAXPY( N, CONE, WORK, 1, X( 1, J ), 1 ) 288: LSTRES = BERR( J ) 289: COUNT = COUNT + 1 290: GO TO 20 291: END IF 292: * 293: * Bound error from formula 294: * 295: * norm(X - XTRUE) / norm(X) .le. FERR = 296: * norm( abs(inv(op(A)))* 297: * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X) 298: * 299: * where 300: * norm(Z) is the magnitude of the largest component of Z 301: * inv(op(A)) is the inverse of op(A) 302: * abs(Z) is the componentwise absolute value of the matrix or 303: * vector Z 304: * NZ is the maximum number of nonzeros in any row of A, plus 1 305: * EPS is machine epsilon 306: * 307: * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B)) 308: * is incremented by SAFE1 if the i-th component of 309: * abs(op(A))*abs(X) + abs(B) is less than SAFE2. 310: * 311: * Use CLACN2 to estimate the infinity-norm of the matrix 312: * inv(op(A)) * diag(W), 313: * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) 314: * 315: DO 90 I = 1, N 316: IF( RWORK( I ).GT.SAFE2 ) THEN 317: RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) 318: ELSE 319: RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) + 320: $ SAFE1 321: END IF 322: 90 CONTINUE 323: * 324: KASE = 0 325: 100 CONTINUE 326: CALL CLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE ) 327: IF( KASE.NE.0 ) THEN 328: IF( KASE.EQ.1 ) THEN 329: * 330: * Multiply by diag(W)*inv(op(A)**H). 331: * 332: CALL CGBTRS( TRANST, N, KL, KU, 1, AFB, LDAFB, IPIV, 333: $ WORK, N, INFO ) 334: DO 110 I = 1, N 335: WORK( I ) = RWORK( I )*WORK( I ) 336: 110 CONTINUE 337: ELSE 338: * 339: * Multiply by inv(op(A))*diag(W). 340: * 341: DO 120 I = 1, N 342: WORK( I ) = RWORK( I )*WORK( I ) 343: 120 CONTINUE 344: CALL CGBTRS( TRANSN, N, KL, KU, 1, AFB, LDAFB, IPIV, 345: $ WORK, N, INFO ) 346: END IF 347: GO TO 100 348: END IF 349: * 350: * Normalize error. 351: * 352: LSTRES = ZERO 353: DO 130 I = 1, N 354: LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) ) 355: 130 CONTINUE 356: IF( LSTRES.NE.ZERO ) 357: $ FERR( J ) = FERR( J ) / LSTRES 358: * 359: 140 CONTINUE 360: * 361: RETURN 362: * 363: * End of CGBRFS 364: * 365: END 366: