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