1*> \brief \b DGBRFS 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DGBRFS + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbrfs.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbrfs.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbrfs.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, 22* IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, 23* INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER TRANS 27* INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS 28* .. 29* .. Array Arguments .. 30* INTEGER IPIV( * ), IWORK( * ) 31* DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 32* $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> DGBRFS improves the computed solution to a system of linear 42*> equations when the coefficient matrix is banded, and provides 43*> error bounds and backward error estimates for the solution. 44*> \endverbatim 45* 46* Arguments: 47* ========== 48* 49*> \param[in] TRANS 50*> \verbatim 51*> TRANS is CHARACTER*1 52*> Specifies the form of the system of equations: 53*> = 'N': A * X = B (No transpose) 54*> = 'T': A**T * X = B (Transpose) 55*> = 'C': A**H * X = B (Conjugate transpose = Transpose) 56*> \endverbatim 57*> 58*> \param[in] N 59*> \verbatim 60*> N is INTEGER 61*> The order of the matrix A. N >= 0. 62*> \endverbatim 63*> 64*> \param[in] KL 65*> \verbatim 66*> KL is INTEGER 67*> The number of subdiagonals within the band of A. KL >= 0. 68*> \endverbatim 69*> 70*> \param[in] KU 71*> \verbatim 72*> KU is INTEGER 73*> The number of superdiagonals within the band of A. KU >= 0. 74*> \endverbatim 75*> 76*> \param[in] NRHS 77*> \verbatim 78*> NRHS is INTEGER 79*> The number of right hand sides, i.e., the number of columns 80*> of the matrices B and X. NRHS >= 0. 81*> \endverbatim 82*> 83*> \param[in] AB 84*> \verbatim 85*> AB is DOUBLE PRECISION array, dimension (LDAB,N) 86*> The original band matrix A, stored in rows 1 to KL+KU+1. 87*> The j-th column of A is stored in the j-th column of the 88*> array AB as follows: 89*> AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl). 90*> \endverbatim 91*> 92*> \param[in] LDAB 93*> \verbatim 94*> LDAB is INTEGER 95*> The leading dimension of the array AB. LDAB >= KL+KU+1. 96*> \endverbatim 97*> 98*> \param[in] AFB 99*> \verbatim 100*> AFB is DOUBLE PRECISION array, dimension (LDAFB,N) 101*> Details of the LU factorization of the band matrix A, as 102*> computed by DGBTRF. U is stored as an upper triangular band 103*> matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and 104*> the multipliers used during the factorization are stored in 105*> rows KL+KU+2 to 2*KL+KU+1. 106*> \endverbatim 107*> 108*> \param[in] LDAFB 109*> \verbatim 110*> LDAFB is INTEGER 111*> The leading dimension of the array AFB. LDAFB >= 2*KL*KU+1. 112*> \endverbatim 113*> 114*> \param[in] IPIV 115*> \verbatim 116*> IPIV is INTEGER array, dimension (N) 117*> The pivot indices from DGBTRF; for 1<=i<=N, row i of the 118*> matrix was interchanged with row IPIV(i). 119*> \endverbatim 120*> 121*> \param[in] B 122*> \verbatim 123*> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 124*> The right hand side matrix B. 125*> \endverbatim 126*> 127*> \param[in] LDB 128*> \verbatim 129*> LDB is INTEGER 130*> The leading dimension of the array B. LDB >= max(1,N). 131*> \endverbatim 132*> 133*> \param[in,out] X 134*> \verbatim 135*> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 136*> On entry, the solution matrix X, as computed by DGBTRS. 137*> On exit, the improved solution matrix X. 138*> \endverbatim 139*> 140*> \param[in] LDX 141*> \verbatim 142*> LDX is INTEGER 143*> The leading dimension of the array X. LDX >= max(1,N). 144*> \endverbatim 145*> 146*> \param[out] FERR 147*> \verbatim 148*> FERR is DOUBLE PRECISION array, dimension (NRHS) 149*> The estimated forward error bound for each solution vector 150*> X(j) (the j-th column of the solution matrix X). 151*> If XTRUE is the true solution corresponding to X(j), FERR(j) 152*> is an estimated upper bound for the magnitude of the largest 153*> element in (X(j) - XTRUE) divided by the magnitude of the 154*> largest element in X(j). The estimate is as reliable as 155*> the estimate for RCOND, and is almost always a slight 156*> overestimate of the true error. 157*> \endverbatim 158*> 159*> \param[out] BERR 160*> \verbatim 161*> BERR is DOUBLE PRECISION array, dimension (NRHS) 162*> The componentwise relative backward error of each solution 163*> vector X(j) (i.e., the smallest relative change in 164*> any element of A or B that makes X(j) an exact solution). 165*> \endverbatim 166*> 167*> \param[out] WORK 168*> \verbatim 169*> WORK is DOUBLE PRECISION array, dimension (3*N) 170*> \endverbatim 171*> 172*> \param[out] IWORK 173*> \verbatim 174*> IWORK is INTEGER array, dimension (N) 175*> \endverbatim 176*> 177*> \param[out] INFO 178*> \verbatim 179*> INFO is INTEGER 180*> = 0: successful exit 181*> < 0: if INFO = -i, the i-th argument had an illegal value 182*> \endverbatim 183* 184*> \par Internal Parameters: 185* ========================= 186*> 187*> \verbatim 188*> ITMAX is the maximum number of steps of iterative refinement. 189*> \endverbatim 190* 191* Authors: 192* ======== 193* 194*> \author Univ. of Tennessee 195*> \author Univ. of California Berkeley 196*> \author Univ. of Colorado Denver 197*> \author NAG Ltd. 198* 199*> \ingroup doubleGBcomputational 200* 201* ===================================================================== 202 SUBROUTINE DGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, 203 $ IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, 204 $ INFO ) 205* 206* -- LAPACK computational routine -- 207* -- LAPACK is a software package provided by Univ. of Tennessee, -- 208* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 209* 210* .. Scalar Arguments .. 211 CHARACTER TRANS 212 INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS 213* .. 214* .. Array Arguments .. 215 INTEGER IPIV( * ), IWORK( * ) 216 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 217 $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * ) 218* .. 219* 220* ===================================================================== 221* 222* .. Parameters .. 223 INTEGER ITMAX 224 PARAMETER ( ITMAX = 5 ) 225 DOUBLE PRECISION ZERO 226 PARAMETER ( ZERO = 0.0D+0 ) 227 DOUBLE PRECISION ONE 228 PARAMETER ( ONE = 1.0D+0 ) 229 DOUBLE PRECISION TWO 230 PARAMETER ( TWO = 2.0D+0 ) 231 DOUBLE PRECISION THREE 232 PARAMETER ( THREE = 3.0D+0 ) 233* .. 234* .. Local Scalars .. 235 LOGICAL NOTRAN 236 CHARACTER TRANST 237 INTEGER COUNT, I, J, K, KASE, KK, NZ 238 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK 239* .. 240* .. Local Arrays .. 241 INTEGER ISAVE( 3 ) 242* .. 243* .. External Subroutines .. 244 EXTERNAL DAXPY, DCOPY, DGBMV, DGBTRS, DLACN2, XERBLA 245* .. 246* .. Intrinsic Functions .. 247 INTRINSIC ABS, MAX, MIN 248* .. 249* .. External Functions .. 250 LOGICAL LSAME 251 DOUBLE PRECISION DLAMCH 252 EXTERNAL LSAME, DLAMCH 253* .. 254* .. Executable Statements .. 255* 256* Test the input parameters. 257* 258 INFO = 0 259 NOTRAN = LSAME( TRANS, 'N' ) 260 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 261 $ LSAME( TRANS, 'C' ) ) THEN 262 INFO = -1 263 ELSE IF( N.LT.0 ) THEN 264 INFO = -2 265 ELSE IF( KL.LT.0 ) THEN 266 INFO = -3 267 ELSE IF( KU.LT.0 ) THEN 268 INFO = -4 269 ELSE IF( NRHS.LT.0 ) THEN 270 INFO = -5 271 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 272 INFO = -7 273 ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN 274 INFO = -9 275 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 276 INFO = -12 277 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 278 INFO = -14 279 END IF 280 IF( INFO.NE.0 ) THEN 281 CALL XERBLA( 'DGBRFS', -INFO ) 282 RETURN 283 END IF 284* 285* Quick return if possible 286* 287 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 288 DO 10 J = 1, NRHS 289 FERR( J ) = ZERO 290 BERR( J ) = ZERO 291 10 CONTINUE 292 RETURN 293 END IF 294* 295 IF( NOTRAN ) THEN 296 TRANST = 'T' 297 ELSE 298 TRANST = 'N' 299 END IF 300* 301* NZ = maximum number of nonzero elements in each row of A, plus 1 302* 303 NZ = MIN( KL+KU+2, N+1 ) 304 EPS = DLAMCH( 'Epsilon' ) 305 SAFMIN = DLAMCH( 'Safe minimum' ) 306 SAFE1 = NZ*SAFMIN 307 SAFE2 = SAFE1 / EPS 308* 309* Do for each right hand side 310* 311 DO 140 J = 1, NRHS 312* 313 COUNT = 1 314 LSTRES = THREE 315 20 CONTINUE 316* 317* Loop until stopping criterion is satisfied. 318* 319* Compute residual R = B - op(A) * X, 320* where op(A) = A, A**T, or A**H, depending on TRANS. 321* 322 CALL DCOPY( N, B( 1, J ), 1, WORK( N+1 ), 1 ) 323 CALL DGBMV( TRANS, N, N, KL, KU, -ONE, AB, LDAB, X( 1, J ), 1, 324 $ ONE, WORK( N+1 ), 1 ) 325* 326* Compute componentwise relative backward error from formula 327* 328* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) ) 329* 330* where abs(Z) is the componentwise absolute value of the matrix 331* or vector Z. If the i-th component of the denominator is less 332* than SAFE2, then SAFE1 is added to the i-th components of the 333* numerator and denominator before dividing. 334* 335 DO 30 I = 1, N 336 WORK( I ) = ABS( B( I, J ) ) 337 30 CONTINUE 338* 339* Compute abs(op(A))*abs(X) + abs(B). 340* 341 IF( NOTRAN ) THEN 342 DO 50 K = 1, N 343 KK = KU + 1 - K 344 XK = ABS( X( K, J ) ) 345 DO 40 I = MAX( 1, K-KU ), MIN( N, K+KL ) 346 WORK( I ) = WORK( I ) + ABS( AB( KK+I, K ) )*XK 347 40 CONTINUE 348 50 CONTINUE 349 ELSE 350 DO 70 K = 1, N 351 S = ZERO 352 KK = KU + 1 - K 353 DO 60 I = MAX( 1, K-KU ), MIN( N, K+KL ) 354 S = S + ABS( AB( KK+I, K ) )*ABS( X( I, J ) ) 355 60 CONTINUE 356 WORK( K ) = WORK( K ) + S 357 70 CONTINUE 358 END IF 359 S = ZERO 360 DO 80 I = 1, N 361 IF( WORK( I ).GT.SAFE2 ) THEN 362 S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) ) 363 ELSE 364 S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) / 365 $ ( WORK( I )+SAFE1 ) ) 366 END IF 367 80 CONTINUE 368 BERR( J ) = S 369* 370* Test stopping criterion. Continue iterating if 371* 1) The residual BERR(J) is larger than machine epsilon, and 372* 2) BERR(J) decreased by at least a factor of 2 during the 373* last iteration, and 374* 3) At most ITMAX iterations tried. 375* 376 IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND. 377 $ COUNT.LE.ITMAX ) THEN 378* 379* Update solution and try again. 380* 381 CALL DGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV, 382 $ WORK( N+1 ), N, INFO ) 383 CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 ) 384 LSTRES = BERR( J ) 385 COUNT = COUNT + 1 386 GO TO 20 387 END IF 388* 389* Bound error from formula 390* 391* norm(X - XTRUE) / norm(X) .le. FERR = 392* norm( abs(inv(op(A)))* 393* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X) 394* 395* where 396* norm(Z) is the magnitude of the largest component of Z 397* inv(op(A)) is the inverse of op(A) 398* abs(Z) is the componentwise absolute value of the matrix or 399* vector Z 400* NZ is the maximum number of nonzeros in any row of A, plus 1 401* EPS is machine epsilon 402* 403* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B)) 404* is incremented by SAFE1 if the i-th component of 405* abs(op(A))*abs(X) + abs(B) is less than SAFE2. 406* 407* Use DLACN2 to estimate the infinity-norm of the matrix 408* inv(op(A)) * diag(W), 409* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) 410* 411 DO 90 I = 1, N 412 IF( WORK( I ).GT.SAFE2 ) THEN 413 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) 414 ELSE 415 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1 416 END IF 417 90 CONTINUE 418* 419 KASE = 0 420 100 CONTINUE 421 CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ), 422 $ KASE, ISAVE ) 423 IF( KASE.NE.0 ) THEN 424 IF( KASE.EQ.1 ) THEN 425* 426* Multiply by diag(W)*inv(op(A)**T). 427* 428 CALL DGBTRS( TRANST, N, KL, KU, 1, AFB, LDAFB, IPIV, 429 $ WORK( N+1 ), N, INFO ) 430 DO 110 I = 1, N 431 WORK( N+I ) = WORK( N+I )*WORK( I ) 432 110 CONTINUE 433 ELSE 434* 435* Multiply by inv(op(A))*diag(W). 436* 437 DO 120 I = 1, N 438 WORK( N+I ) = WORK( N+I )*WORK( I ) 439 120 CONTINUE 440 CALL DGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV, 441 $ WORK( N+1 ), N, INFO ) 442 END IF 443 GO TO 100 444 END IF 445* 446* Normalize error. 447* 448 LSTRES = ZERO 449 DO 130 I = 1, N 450 LSTRES = MAX( LSTRES, ABS( X( I, J ) ) ) 451 130 CONTINUE 452 IF( LSTRES.NE.ZERO ) 453 $ FERR( J ) = FERR( J ) / LSTRES 454* 455 140 CONTINUE 456* 457 RETURN 458* 459* End of DGBRFS 460* 461 END 462