![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DTBRFS 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DTBRFS + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtbrfs.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtbrfs.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtbrfs.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B, 00022 * LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER DIAG, TRANS, UPLO 00026 * INTEGER INFO, KD, LDAB, LDB, LDX, N, NRHS 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * DOUBLE PRECISION AB( LDAB, * ), B( LDB, * ), BERR( * ), 00031 * $ FERR( * ), WORK( * ), X( LDX, * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> DTBRFS provides error bounds and backward error estimates for the 00041 *> solution to a system of linear equations with a triangular band 00042 *> coefficient matrix. 00043 *> 00044 *> The solution matrix X must be computed by DTBTRS or some other 00045 *> means before entering this routine. DTBRFS does not do iterative 00046 *> refinement because doing so cannot improve the backward error. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] UPLO 00053 *> \verbatim 00054 *> UPLO is CHARACTER*1 00055 *> = 'U': A is upper triangular; 00056 *> = 'L': A is lower triangular. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] TRANS 00060 *> \verbatim 00061 *> TRANS is CHARACTER*1 00062 *> Specifies the form of the system of equations: 00063 *> = 'N': A * X = B (No transpose) 00064 *> = 'T': A**T * X = B (Transpose) 00065 *> = 'C': A**H * X = B (Conjugate transpose = Transpose) 00066 *> \endverbatim 00067 *> 00068 *> \param[in] DIAG 00069 *> \verbatim 00070 *> DIAG is CHARACTER*1 00071 *> = 'N': A is non-unit triangular; 00072 *> = 'U': A is unit triangular. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] N 00076 *> \verbatim 00077 *> N is INTEGER 00078 *> The order of the matrix A. N >= 0. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] KD 00082 *> \verbatim 00083 *> KD is INTEGER 00084 *> The number of superdiagonals or subdiagonals of the 00085 *> triangular band matrix A. KD >= 0. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] NRHS 00089 *> \verbatim 00090 *> NRHS is INTEGER 00091 *> The number of right hand sides, i.e., the number of columns 00092 *> of the matrices B and X. NRHS >= 0. 00093 *> \endverbatim 00094 *> 00095 *> \param[in] AB 00096 *> \verbatim 00097 *> AB is DOUBLE PRECISION array, dimension (LDAB,N) 00098 *> The upper or lower triangular band matrix A, stored in the 00099 *> first kd+1 rows of the array. The j-th column of A is stored 00100 *> in the j-th column of the array AB as follows: 00101 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00102 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00103 *> If DIAG = 'U', the diagonal elements of A are not referenced 00104 *> and are assumed to be 1. 00105 *> \endverbatim 00106 *> 00107 *> \param[in] LDAB 00108 *> \verbatim 00109 *> LDAB is INTEGER 00110 *> The leading dimension of the array AB. LDAB >= KD+1. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] B 00114 *> \verbatim 00115 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00116 *> The right hand side matrix B. 00117 *> \endverbatim 00118 *> 00119 *> \param[in] LDB 00120 *> \verbatim 00121 *> LDB is INTEGER 00122 *> The leading dimension of the array B. LDB >= max(1,N). 00123 *> \endverbatim 00124 *> 00125 *> \param[in] X 00126 *> \verbatim 00127 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 00128 *> The solution matrix X. 00129 *> \endverbatim 00130 *> 00131 *> \param[in] LDX 00132 *> \verbatim 00133 *> LDX is INTEGER 00134 *> The leading dimension of the array X. LDX >= max(1,N). 00135 *> \endverbatim 00136 *> 00137 *> \param[out] FERR 00138 *> \verbatim 00139 *> FERR is DOUBLE PRECISION array, dimension (NRHS) 00140 *> The estimated forward error bound for each solution vector 00141 *> X(j) (the j-th column of the solution matrix X). 00142 *> If XTRUE is the true solution corresponding to X(j), FERR(j) 00143 *> is an estimated upper bound for the magnitude of the largest 00144 *> element in (X(j) - XTRUE) divided by the magnitude of the 00145 *> largest element in X(j). The estimate is as reliable as 00146 *> the estimate for RCOND, and is almost always a slight 00147 *> overestimate of the true error. 00148 *> \endverbatim 00149 *> 00150 *> \param[out] BERR 00151 *> \verbatim 00152 *> BERR is DOUBLE PRECISION array, dimension (NRHS) 00153 *> The componentwise relative backward error of each solution 00154 *> vector X(j) (i.e., the smallest relative change in 00155 *> any element of A or B that makes X(j) an exact solution). 00156 *> \endverbatim 00157 *> 00158 *> \param[out] WORK 00159 *> \verbatim 00160 *> WORK is DOUBLE PRECISION array, dimension (3*N) 00161 *> \endverbatim 00162 *> 00163 *> \param[out] IWORK 00164 *> \verbatim 00165 *> IWORK is INTEGER array, dimension (N) 00166 *> \endverbatim 00167 *> 00168 *> \param[out] INFO 00169 *> \verbatim 00170 *> INFO is INTEGER 00171 *> = 0: successful exit 00172 *> < 0: if INFO = -i, the i-th argument had an illegal value 00173 *> \endverbatim 00174 * 00175 * Authors: 00176 * ======== 00177 * 00178 *> \author Univ. of Tennessee 00179 *> \author Univ. of California Berkeley 00180 *> \author Univ. of Colorado Denver 00181 *> \author NAG Ltd. 00182 * 00183 *> \date November 2011 00184 * 00185 *> \ingroup doubleOTHERcomputational 00186 * 00187 * ===================================================================== 00188 SUBROUTINE DTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B, 00189 $ LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO ) 00190 * 00191 * -- LAPACK computational routine (version 3.4.0) -- 00192 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00193 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00194 * November 2011 00195 * 00196 * .. Scalar Arguments .. 00197 CHARACTER DIAG, TRANS, UPLO 00198 INTEGER INFO, KD, LDAB, LDB, LDX, N, NRHS 00199 * .. 00200 * .. Array Arguments .. 00201 INTEGER IWORK( * ) 00202 DOUBLE PRECISION AB( LDAB, * ), B( LDB, * ), BERR( * ), 00203 $ FERR( * ), WORK( * ), X( LDX, * ) 00204 * .. 00205 * 00206 * ===================================================================== 00207 * 00208 * .. Parameters .. 00209 DOUBLE PRECISION ZERO 00210 PARAMETER ( ZERO = 0.0D+0 ) 00211 DOUBLE PRECISION ONE 00212 PARAMETER ( ONE = 1.0D+0 ) 00213 * .. 00214 * .. Local Scalars .. 00215 LOGICAL NOTRAN, NOUNIT, UPPER 00216 CHARACTER TRANST 00217 INTEGER I, J, K, KASE, NZ 00218 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK 00219 * .. 00220 * .. Local Arrays .. 00221 INTEGER ISAVE( 3 ) 00222 * .. 00223 * .. External Subroutines .. 00224 EXTERNAL DAXPY, DCOPY, DLACN2, DTBMV, DTBSV, XERBLA 00225 * .. 00226 * .. Intrinsic Functions .. 00227 INTRINSIC ABS, MAX, MIN 00228 * .. 00229 * .. External Functions .. 00230 LOGICAL LSAME 00231 DOUBLE PRECISION DLAMCH 00232 EXTERNAL LSAME, DLAMCH 00233 * .. 00234 * .. Executable Statements .. 00235 * 00236 * Test the input parameters. 00237 * 00238 INFO = 0 00239 UPPER = LSAME( UPLO, 'U' ) 00240 NOTRAN = LSAME( TRANS, 'N' ) 00241 NOUNIT = LSAME( DIAG, 'N' ) 00242 * 00243 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00244 INFO = -1 00245 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 00246 $ LSAME( TRANS, 'C' ) ) THEN 00247 INFO = -2 00248 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 00249 INFO = -3 00250 ELSE IF( N.LT.0 ) THEN 00251 INFO = -4 00252 ELSE IF( KD.LT.0 ) THEN 00253 INFO = -5 00254 ELSE IF( NRHS.LT.0 ) THEN 00255 INFO = -6 00256 ELSE IF( LDAB.LT.KD+1 ) THEN 00257 INFO = -8 00258 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00259 INFO = -10 00260 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00261 INFO = -12 00262 END IF 00263 IF( INFO.NE.0 ) THEN 00264 CALL XERBLA( 'DTBRFS', -INFO ) 00265 RETURN 00266 END IF 00267 * 00268 * Quick return if possible 00269 * 00270 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 00271 DO 10 J = 1, NRHS 00272 FERR( J ) = ZERO 00273 BERR( J ) = ZERO 00274 10 CONTINUE 00275 RETURN 00276 END IF 00277 * 00278 IF( NOTRAN ) THEN 00279 TRANST = 'T' 00280 ELSE 00281 TRANST = 'N' 00282 END IF 00283 * 00284 * NZ = maximum number of nonzero elements in each row of A, plus 1 00285 * 00286 NZ = KD + 2 00287 EPS = DLAMCH( 'Epsilon' ) 00288 SAFMIN = DLAMCH( 'Safe minimum' ) 00289 SAFE1 = NZ*SAFMIN 00290 SAFE2 = SAFE1 / EPS 00291 * 00292 * Do for each right hand side 00293 * 00294 DO 250 J = 1, NRHS 00295 * 00296 * Compute residual R = B - op(A) * X, 00297 * where op(A) = A or A**T, depending on TRANS. 00298 * 00299 CALL DCOPY( N, X( 1, J ), 1, WORK( N+1 ), 1 ) 00300 CALL DTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK( N+1 ), 00301 $ 1 ) 00302 CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK( N+1 ), 1 ) 00303 * 00304 * Compute componentwise relative backward error from formula 00305 * 00306 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) ) 00307 * 00308 * where abs(Z) is the componentwise absolute value of the matrix 00309 * or vector Z. If the i-th component of the denominator is less 00310 * than SAFE2, then SAFE1 is added to the i-th components of the 00311 * numerator and denominator before dividing. 00312 * 00313 DO 20 I = 1, N 00314 WORK( I ) = ABS( B( I, J ) ) 00315 20 CONTINUE 00316 * 00317 IF( NOTRAN ) THEN 00318 * 00319 * Compute abs(A)*abs(X) + abs(B). 00320 * 00321 IF( UPPER ) THEN 00322 IF( NOUNIT ) THEN 00323 DO 40 K = 1, N 00324 XK = ABS( X( K, J ) ) 00325 DO 30 I = MAX( 1, K-KD ), K 00326 WORK( I ) = WORK( I ) + 00327 $ ABS( AB( KD+1+I-K, K ) )*XK 00328 30 CONTINUE 00329 40 CONTINUE 00330 ELSE 00331 DO 60 K = 1, N 00332 XK = ABS( X( K, J ) ) 00333 DO 50 I = MAX( 1, K-KD ), K - 1 00334 WORK( I ) = WORK( I ) + 00335 $ ABS( AB( KD+1+I-K, K ) )*XK 00336 50 CONTINUE 00337 WORK( K ) = WORK( K ) + XK 00338 60 CONTINUE 00339 END IF 00340 ELSE 00341 IF( NOUNIT ) THEN 00342 DO 80 K = 1, N 00343 XK = ABS( X( K, J ) ) 00344 DO 70 I = K, MIN( N, K+KD ) 00345 WORK( I ) = WORK( I ) + ABS( AB( 1+I-K, K ) )*XK 00346 70 CONTINUE 00347 80 CONTINUE 00348 ELSE 00349 DO 100 K = 1, N 00350 XK = ABS( X( K, J ) ) 00351 DO 90 I = K + 1, MIN( N, K+KD ) 00352 WORK( I ) = WORK( I ) + ABS( AB( 1+I-K, K ) )*XK 00353 90 CONTINUE 00354 WORK( K ) = WORK( K ) + XK 00355 100 CONTINUE 00356 END IF 00357 END IF 00358 ELSE 00359 * 00360 * Compute abs(A**T)*abs(X) + abs(B). 00361 * 00362 IF( UPPER ) THEN 00363 IF( NOUNIT ) THEN 00364 DO 120 K = 1, N 00365 S = ZERO 00366 DO 110 I = MAX( 1, K-KD ), K 00367 S = S + ABS( AB( KD+1+I-K, K ) )* 00368 $ ABS( X( I, J ) ) 00369 110 CONTINUE 00370 WORK( K ) = WORK( K ) + S 00371 120 CONTINUE 00372 ELSE 00373 DO 140 K = 1, N 00374 S = ABS( X( K, J ) ) 00375 DO 130 I = MAX( 1, K-KD ), K - 1 00376 S = S + ABS( AB( KD+1+I-K, K ) )* 00377 $ ABS( X( I, J ) ) 00378 130 CONTINUE 00379 WORK( K ) = WORK( K ) + S 00380 140 CONTINUE 00381 END IF 00382 ELSE 00383 IF( NOUNIT ) THEN 00384 DO 160 K = 1, N 00385 S = ZERO 00386 DO 150 I = K, MIN( N, K+KD ) 00387 S = S + ABS( AB( 1+I-K, K ) )*ABS( X( I, J ) ) 00388 150 CONTINUE 00389 WORK( K ) = WORK( K ) + S 00390 160 CONTINUE 00391 ELSE 00392 DO 180 K = 1, N 00393 S = ABS( X( K, J ) ) 00394 DO 170 I = K + 1, MIN( N, K+KD ) 00395 S = S + ABS( AB( 1+I-K, K ) )*ABS( X( I, J ) ) 00396 170 CONTINUE 00397 WORK( K ) = WORK( K ) + S 00398 180 CONTINUE 00399 END IF 00400 END IF 00401 END IF 00402 S = ZERO 00403 DO 190 I = 1, N 00404 IF( WORK( I ).GT.SAFE2 ) THEN 00405 S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) ) 00406 ELSE 00407 S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) / 00408 $ ( WORK( I )+SAFE1 ) ) 00409 END IF 00410 190 CONTINUE 00411 BERR( J ) = S 00412 * 00413 * Bound error from formula 00414 * 00415 * norm(X - XTRUE) / norm(X) .le. FERR = 00416 * norm( abs(inv(op(A)))* 00417 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X) 00418 * 00419 * where 00420 * norm(Z) is the magnitude of the largest component of Z 00421 * inv(op(A)) is the inverse of op(A) 00422 * abs(Z) is the componentwise absolute value of the matrix or 00423 * vector Z 00424 * NZ is the maximum number of nonzeros in any row of A, plus 1 00425 * EPS is machine epsilon 00426 * 00427 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B)) 00428 * is incremented by SAFE1 if the i-th component of 00429 * abs(op(A))*abs(X) + abs(B) is less than SAFE2. 00430 * 00431 * Use DLACN2 to estimate the infinity-norm of the matrix 00432 * inv(op(A)) * diag(W), 00433 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) 00434 * 00435 DO 200 I = 1, N 00436 IF( WORK( I ).GT.SAFE2 ) THEN 00437 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) 00438 ELSE 00439 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1 00440 END IF 00441 200 CONTINUE 00442 * 00443 KASE = 0 00444 210 CONTINUE 00445 CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ), 00446 $ KASE, ISAVE ) 00447 IF( KASE.NE.0 ) THEN 00448 IF( KASE.EQ.1 ) THEN 00449 * 00450 * Multiply by diag(W)*inv(op(A)**T). 00451 * 00452 CALL DTBSV( UPLO, TRANST, DIAG, N, KD, AB, LDAB, 00453 $ WORK( N+1 ), 1 ) 00454 DO 220 I = 1, N 00455 WORK( N+I ) = WORK( I )*WORK( N+I ) 00456 220 CONTINUE 00457 ELSE 00458 * 00459 * Multiply by inv(op(A))*diag(W). 00460 * 00461 DO 230 I = 1, N 00462 WORK( N+I ) = WORK( I )*WORK( N+I ) 00463 230 CONTINUE 00464 CALL DTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, 00465 $ WORK( N+1 ), 1 ) 00466 END IF 00467 GO TO 210 00468 END IF 00469 * 00470 * Normalize error. 00471 * 00472 LSTRES = ZERO 00473 DO 240 I = 1, N 00474 LSTRES = MAX( LSTRES, ABS( X( I, J ) ) ) 00475 240 CONTINUE 00476 IF( LSTRES.NE.ZERO ) 00477 $ FERR( J ) = FERR( J ) / LSTRES 00478 * 00479 250 CONTINUE 00480 * 00481 RETURN 00482 * 00483 * End of DTBRFS 00484 * 00485 END