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