![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTRRFS 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZTRRFS + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrrfs.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrrfs.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrrfs.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, 00022 * LDX, FERR, BERR, WORK, RWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER DIAG, TRANS, UPLO 00026 * INTEGER INFO, LDA, LDB, LDX, N, NRHS 00027 * .. 00028 * .. Array Arguments .. 00029 * DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ) 00030 * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ), 00031 * $ X( LDX, * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> ZTRRFS 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 ZTRTRS or some other 00045 *> means before entering this routine. ZTRRFS 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) 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (2*N) 00155 *> \endverbatim 00156 *> 00157 *> \param[out] RWORK 00158 *> \verbatim 00159 *> RWORK is DOUBLE PRECISION 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 complex16OTHERcomputational 00180 * 00181 * ===================================================================== 00182 SUBROUTINE ZTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, 00183 $ LDX, FERR, BERR, WORK, RWORK, 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 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ) 00196 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ), 00197 $ X( LDX, * ) 00198 * .. 00199 * 00200 * ===================================================================== 00201 * 00202 * .. Parameters .. 00203 DOUBLE PRECISION ZERO 00204 PARAMETER ( ZERO = 0.0D+0 ) 00205 COMPLEX*16 ONE 00206 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) 00207 * .. 00208 * .. Local Scalars .. 00209 LOGICAL NOTRAN, NOUNIT, UPPER 00210 CHARACTER TRANSN, TRANST 00211 INTEGER I, J, K, KASE, NZ 00212 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK 00213 COMPLEX*16 ZDUM 00214 * .. 00215 * .. Local Arrays .. 00216 INTEGER ISAVE( 3 ) 00217 * .. 00218 * .. External Subroutines .. 00219 EXTERNAL XERBLA, ZAXPY, ZCOPY, ZLACN2, ZTRMV, ZTRSV 00220 * .. 00221 * .. Intrinsic Functions .. 00222 INTRINSIC ABS, DBLE, DIMAG, MAX 00223 * .. 00224 * .. External Functions .. 00225 LOGICAL LSAME 00226 DOUBLE PRECISION DLAMCH 00227 EXTERNAL LSAME, DLAMCH 00228 * .. 00229 * .. Statement Functions .. 00230 DOUBLE PRECISION CABS1 00231 * .. 00232 * .. Statement Function definitions .. 00233 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00234 * .. 00235 * .. Executable Statements .. 00236 * 00237 * Test the input parameters. 00238 * 00239 INFO = 0 00240 UPPER = LSAME( UPLO, 'U' ) 00241 NOTRAN = LSAME( TRANS, 'N' ) 00242 NOUNIT = LSAME( DIAG, 'N' ) 00243 * 00244 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00245 INFO = -1 00246 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 00247 $ LSAME( TRANS, 'C' ) ) THEN 00248 INFO = -2 00249 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 00250 INFO = -3 00251 ELSE IF( N.LT.0 ) THEN 00252 INFO = -4 00253 ELSE IF( NRHS.LT.0 ) THEN 00254 INFO = -5 00255 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00256 INFO = -7 00257 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00258 INFO = -9 00259 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00260 INFO = -11 00261 END IF 00262 IF( INFO.NE.0 ) THEN 00263 CALL XERBLA( 'ZTRRFS', -INFO ) 00264 RETURN 00265 END IF 00266 * 00267 * Quick return if possible 00268 * 00269 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN 00270 DO 10 J = 1, NRHS 00271 FERR( J ) = ZERO 00272 BERR( J ) = ZERO 00273 10 CONTINUE 00274 RETURN 00275 END IF 00276 * 00277 IF( NOTRAN ) THEN 00278 TRANSN = 'N' 00279 TRANST = 'C' 00280 ELSE 00281 TRANSN = 'C' 00282 TRANST = 'N' 00283 END IF 00284 * 00285 * NZ = maximum number of nonzero elements in each row of A, plus 1 00286 * 00287 NZ = N + 1 00288 EPS = DLAMCH( 'Epsilon' ) 00289 SAFMIN = DLAMCH( 'Safe minimum' ) 00290 SAFE1 = NZ*SAFMIN 00291 SAFE2 = SAFE1 / EPS 00292 * 00293 * Do for each right hand side 00294 * 00295 DO 250 J = 1, NRHS 00296 * 00297 * Compute residual R = B - op(A) * X, 00298 * where op(A) = A, A**T, or A**H, depending on TRANS. 00299 * 00300 CALL ZCOPY( N, X( 1, J ), 1, WORK, 1 ) 00301 CALL ZTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 ) 00302 CALL ZAXPY( N, -ONE, B( 1, J ), 1, WORK, 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 RWORK( I ) = CABS1( 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 = CABS1( X( K, J ) ) 00325 DO 30 I = 1, K 00326 RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK 00327 30 CONTINUE 00328 40 CONTINUE 00329 ELSE 00330 DO 60 K = 1, N 00331 XK = CABS1( X( K, J ) ) 00332 DO 50 I = 1, K - 1 00333 RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK 00334 50 CONTINUE 00335 RWORK( K ) = RWORK( K ) + XK 00336 60 CONTINUE 00337 END IF 00338 ELSE 00339 IF( NOUNIT ) THEN 00340 DO 80 K = 1, N 00341 XK = CABS1( X( K, J ) ) 00342 DO 70 I = K, N 00343 RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK 00344 70 CONTINUE 00345 80 CONTINUE 00346 ELSE 00347 DO 100 K = 1, N 00348 XK = CABS1( X( K, J ) ) 00349 DO 90 I = K + 1, N 00350 RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK 00351 90 CONTINUE 00352 RWORK( K ) = RWORK( K ) + XK 00353 100 CONTINUE 00354 END IF 00355 END IF 00356 ELSE 00357 * 00358 * Compute abs(A**H)*abs(X) + abs(B). 00359 * 00360 IF( UPPER ) THEN 00361 IF( NOUNIT ) THEN 00362 DO 120 K = 1, N 00363 S = ZERO 00364 DO 110 I = 1, K 00365 S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) ) 00366 110 CONTINUE 00367 RWORK( K ) = RWORK( K ) + S 00368 120 CONTINUE 00369 ELSE 00370 DO 140 K = 1, N 00371 S = CABS1( X( K, J ) ) 00372 DO 130 I = 1, K - 1 00373 S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) ) 00374 130 CONTINUE 00375 RWORK( K ) = RWORK( K ) + S 00376 140 CONTINUE 00377 END IF 00378 ELSE 00379 IF( NOUNIT ) THEN 00380 DO 160 K = 1, N 00381 S = ZERO 00382 DO 150 I = K, N 00383 S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) ) 00384 150 CONTINUE 00385 RWORK( K ) = RWORK( K ) + S 00386 160 CONTINUE 00387 ELSE 00388 DO 180 K = 1, N 00389 S = CABS1( X( K, J ) ) 00390 DO 170 I = K + 1, N 00391 S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) ) 00392 170 CONTINUE 00393 RWORK( K ) = RWORK( K ) + S 00394 180 CONTINUE 00395 END IF 00396 END IF 00397 END IF 00398 S = ZERO 00399 DO 190 I = 1, N 00400 IF( RWORK( I ).GT.SAFE2 ) THEN 00401 S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) ) 00402 ELSE 00403 S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) / 00404 $ ( RWORK( I )+SAFE1 ) ) 00405 END IF 00406 190 CONTINUE 00407 BERR( J ) = S 00408 * 00409 * Bound error from formula 00410 * 00411 * norm(X - XTRUE) / norm(X) .le. FERR = 00412 * norm( abs(inv(op(A)))* 00413 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X) 00414 * 00415 * where 00416 * norm(Z) is the magnitude of the largest component of Z 00417 * inv(op(A)) is the inverse of op(A) 00418 * abs(Z) is the componentwise absolute value of the matrix or 00419 * vector Z 00420 * NZ is the maximum number of nonzeros in any row of A, plus 1 00421 * EPS is machine epsilon 00422 * 00423 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B)) 00424 * is incremented by SAFE1 if the i-th component of 00425 * abs(op(A))*abs(X) + abs(B) is less than SAFE2. 00426 * 00427 * Use ZLACN2 to estimate the infinity-norm of the matrix 00428 * inv(op(A)) * diag(W), 00429 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) 00430 * 00431 DO 200 I = 1, N 00432 IF( RWORK( I ).GT.SAFE2 ) THEN 00433 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) 00434 ELSE 00435 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) + 00436 $ SAFE1 00437 END IF 00438 200 CONTINUE 00439 * 00440 KASE = 0 00441 210 CONTINUE 00442 CALL ZLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE ) 00443 IF( KASE.NE.0 ) THEN 00444 IF( KASE.EQ.1 ) THEN 00445 * 00446 * Multiply by diag(W)*inv(op(A)**H). 00447 * 00448 CALL ZTRSV( UPLO, TRANST, DIAG, N, A, LDA, WORK, 1 ) 00449 DO 220 I = 1, N 00450 WORK( I ) = RWORK( I )*WORK( I ) 00451 220 CONTINUE 00452 ELSE 00453 * 00454 * Multiply by inv(op(A))*diag(W). 00455 * 00456 DO 230 I = 1, N 00457 WORK( I ) = RWORK( I )*WORK( I ) 00458 230 CONTINUE 00459 CALL ZTRSV( UPLO, TRANSN, DIAG, N, A, LDA, WORK, 1 ) 00460 END IF 00461 GO TO 210 00462 END IF 00463 * 00464 * Normalize error. 00465 * 00466 LSTRES = ZERO 00467 DO 240 I = 1, N 00468 LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) ) 00469 240 CONTINUE 00470 IF( LSTRES.NE.ZERO ) 00471 $ FERR( J ) = FERR( J ) / LSTRES 00472 * 00473 250 CONTINUE 00474 * 00475 RETURN 00476 * 00477 * End of ZTRRFS 00478 * 00479 END