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