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