|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, 00002 $ 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 DIAG, TRANS, UPLO 00013 INTEGER INFO, LDA, LDB, LDX, N, NRHS 00014 * .. 00015 * .. Array Arguments .. 00016 REAL BERR( * ), FERR( * ), RWORK( * ) 00017 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ), 00018 $ X( LDX, * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * CTRRFS provides error bounds and backward error estimates for the 00025 * solution to a system of linear equations with a triangular 00026 * coefficient matrix. 00027 * 00028 * The solution matrix X must be computed by CTRTRS or some other 00029 * means before entering this routine. CTRRFS does not do iterative 00030 * refinement because doing so cannot improve the backward error. 00031 * 00032 * Arguments 00033 * ========= 00034 * 00035 * UPLO (input) CHARACTER*1 00036 * = 'U': A is upper triangular; 00037 * = 'L': A is lower triangular. 00038 * 00039 * TRANS (input) CHARACTER*1 00040 * Specifies the form of the system of equations: 00041 * = 'N': A * X = B (No transpose) 00042 * = 'T': A**T * X = B (Transpose) 00043 * = 'C': A**H * X = B (Conjugate transpose) 00044 * 00045 * DIAG (input) CHARACTER*1 00046 * = 'N': A is non-unit triangular; 00047 * = 'U': A is unit triangular. 00048 * 00049 * N (input) INTEGER 00050 * The order of the matrix A. N >= 0. 00051 * 00052 * NRHS (input) INTEGER 00053 * The number of right hand sides, i.e., the number of columns 00054 * of the matrices B and X. NRHS >= 0. 00055 * 00056 * A (input) COMPLEX array, dimension (LDA,N) 00057 * The triangular matrix A. If UPLO = 'U', the leading N-by-N 00058 * upper triangular part of the array A contains the upper 00059 * triangular matrix, and the strictly lower triangular part of 00060 * A is not referenced. If UPLO = 'L', the leading N-by-N lower 00061 * triangular part of the array A contains the lower triangular 00062 * matrix, and the strictly upper triangular part of A is not 00063 * referenced. If DIAG = 'U', the diagonal elements of A are 00064 * also not referenced and are assumed to be 1. 00065 * 00066 * LDA (input) INTEGER 00067 * The leading dimension of the array A. LDA >= max(1,N). 00068 * 00069 * B (input) COMPLEX array, dimension (LDB,NRHS) 00070 * The right hand side matrix B. 00071 * 00072 * LDB (input) INTEGER 00073 * The leading dimension of the array B. LDB >= max(1,N). 00074 * 00075 * X (input) COMPLEX array, dimension (LDX,NRHS) 00076 * The solution matrix X. 00077 * 00078 * LDX (input) INTEGER 00079 * The leading dimension of the array X. LDX >= max(1,N). 00080 * 00081 * FERR (output) REAL array, dimension (NRHS) 00082 * The estimated forward error bound for each solution vector 00083 * X(j) (the j-th column of the solution matrix X). 00084 * If XTRUE is the true solution corresponding to X(j), FERR(j) 00085 * is an estimated upper bound for the magnitude of the largest 00086 * element in (X(j) - XTRUE) divided by the magnitude of the 00087 * largest element in X(j). The estimate is as reliable as 00088 * the estimate for RCOND, and is almost always a slight 00089 * overestimate of the true error. 00090 * 00091 * BERR (output) REAL array, dimension (NRHS) 00092 * The componentwise relative backward error of each solution 00093 * vector X(j) (i.e., the smallest relative change in 00094 * any element of A or B that makes X(j) an exact solution). 00095 * 00096 * WORK (workspace) COMPLEX array, dimension (2*N) 00097 * 00098 * RWORK (workspace) REAL array, dimension (N) 00099 * 00100 * INFO (output) INTEGER 00101 * = 0: successful exit 00102 * < 0: if INFO = -i, the i-th argument had an illegal value 00103 * 00104 * ===================================================================== 00105 * 00106 * .. Parameters .. 00107 REAL ZERO 00108 PARAMETER ( ZERO = 0.0E+0 ) 00109 COMPLEX ONE 00110 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) 00111 * .. 00112 * .. Local Scalars .. 00113 LOGICAL NOTRAN, NOUNIT, UPPER 00114 CHARACTER TRANSN, TRANST 00115 INTEGER I, J, K, KASE, NZ 00116 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK 00117 COMPLEX ZDUM 00118 * .. 00119 * .. Local Arrays .. 00120 INTEGER ISAVE( 3 ) 00121 * .. 00122 * .. External Subroutines .. 00123 EXTERNAL CAXPY, CCOPY, CLACN2, CTRMV, CTRSV, XERBLA 00124 * .. 00125 * .. Intrinsic Functions .. 00126 INTRINSIC ABS, AIMAG, MAX, REAL 00127 * .. 00128 * .. External Functions .. 00129 LOGICAL LSAME 00130 REAL SLAMCH 00131 EXTERNAL LSAME, SLAMCH 00132 * .. 00133 * .. Statement Functions .. 00134 REAL CABS1 00135 * .. 00136 * .. Statement Function definitions .. 00137 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 00138 * .. 00139 * .. Executable Statements .. 00140 * 00141 * Test the input parameters. 00142 * 00143 INFO = 0 00144 UPPER = LSAME( UPLO, 'U' ) 00145 NOTRAN = LSAME( TRANS, 'N' ) 00146 NOUNIT = LSAME( DIAG, 'N' ) 00147 * 00148 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00149 INFO = -1 00150 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 00151 $ LSAME( TRANS, 'C' ) ) THEN 00152 INFO = -2 00153 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 00154 INFO = -3 00155 ELSE IF( N.LT.0 ) THEN 00156 INFO = -4 00157 ELSE IF( NRHS.LT.0 ) THEN 00158 INFO = -5 00159 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00160 INFO = -7 00161 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00162 INFO = -9 00163 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00164 INFO = -11 00165 END IF 00166 IF( INFO.NE.0 ) THEN 00167 CALL XERBLA( 'CTRRFS', -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 250 J = 1, NRHS 00200 * 00201 * Compute residual R = B - op(A) * X, 00202 * where op(A) = A, A**T, or A**H, depending on TRANS. 00203 * 00204 CALL CCOPY( N, X( 1, J ), 1, WORK, 1 ) 00205 CALL CTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 ) 00206 CALL CAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 ) 00207 * 00208 * Compute componentwise relative backward error from formula 00209 * 00210 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) ) 00211 * 00212 * where abs(Z) is the componentwise absolute value of the matrix 00213 * or vector Z. If the i-th component of the denominator is less 00214 * than SAFE2, then SAFE1 is added to the i-th components of the 00215 * numerator and denominator before dividing. 00216 * 00217 DO 20 I = 1, N 00218 RWORK( I ) = CABS1( B( I, J ) ) 00219 20 CONTINUE 00220 * 00221 IF( NOTRAN ) THEN 00222 * 00223 * Compute abs(A)*abs(X) + abs(B). 00224 * 00225 IF( UPPER ) THEN 00226 IF( NOUNIT ) THEN 00227 DO 40 K = 1, N 00228 XK = CABS1( X( K, J ) ) 00229 DO 30 I = 1, K 00230 RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK 00231 30 CONTINUE 00232 40 CONTINUE 00233 ELSE 00234 DO 60 K = 1, N 00235 XK = CABS1( X( K, J ) ) 00236 DO 50 I = 1, K - 1 00237 RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK 00238 50 CONTINUE 00239 RWORK( K ) = RWORK( K ) + XK 00240 60 CONTINUE 00241 END IF 00242 ELSE 00243 IF( NOUNIT ) THEN 00244 DO 80 K = 1, N 00245 XK = CABS1( X( K, J ) ) 00246 DO 70 I = K, N 00247 RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK 00248 70 CONTINUE 00249 80 CONTINUE 00250 ELSE 00251 DO 100 K = 1, N 00252 XK = CABS1( X( K, J ) ) 00253 DO 90 I = K + 1, N 00254 RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK 00255 90 CONTINUE 00256 RWORK( K ) = RWORK( K ) + XK 00257 100 CONTINUE 00258 END IF 00259 END IF 00260 ELSE 00261 * 00262 * Compute abs(A**H)*abs(X) + abs(B). 00263 * 00264 IF( UPPER ) THEN 00265 IF( NOUNIT ) THEN 00266 DO 120 K = 1, N 00267 S = ZERO 00268 DO 110 I = 1, K 00269 S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) ) 00270 110 CONTINUE 00271 RWORK( K ) = RWORK( K ) + S 00272 120 CONTINUE 00273 ELSE 00274 DO 140 K = 1, N 00275 S = CABS1( X( K, J ) ) 00276 DO 130 I = 1, K - 1 00277 S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) ) 00278 130 CONTINUE 00279 RWORK( K ) = RWORK( K ) + S 00280 140 CONTINUE 00281 END IF 00282 ELSE 00283 IF( NOUNIT ) THEN 00284 DO 160 K = 1, N 00285 S = ZERO 00286 DO 150 I = K, N 00287 S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) ) 00288 150 CONTINUE 00289 RWORK( K ) = RWORK( K ) + S 00290 160 CONTINUE 00291 ELSE 00292 DO 180 K = 1, N 00293 S = CABS1( X( K, J ) ) 00294 DO 170 I = K + 1, N 00295 S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) ) 00296 170 CONTINUE 00297 RWORK( K ) = RWORK( K ) + S 00298 180 CONTINUE 00299 END IF 00300 END IF 00301 END IF 00302 S = ZERO 00303 DO 190 I = 1, N 00304 IF( RWORK( I ).GT.SAFE2 ) THEN 00305 S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) ) 00306 ELSE 00307 S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) / 00308 $ ( RWORK( I )+SAFE1 ) ) 00309 END IF 00310 190 CONTINUE 00311 BERR( J ) = S 00312 * 00313 * Bound error from formula 00314 * 00315 * norm(X - XTRUE) / norm(X) .le. FERR = 00316 * norm( abs(inv(op(A)))* 00317 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X) 00318 * 00319 * where 00320 * norm(Z) is the magnitude of the largest component of Z 00321 * inv(op(A)) is the inverse of op(A) 00322 * abs(Z) is the componentwise absolute value of the matrix or 00323 * vector Z 00324 * NZ is the maximum number of nonzeros in any row of A, plus 1 00325 * EPS is machine epsilon 00326 * 00327 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B)) 00328 * is incremented by SAFE1 if the i-th component of 00329 * abs(op(A))*abs(X) + abs(B) is less than SAFE2. 00330 * 00331 * Use CLACN2 to estimate the infinity-norm of the matrix 00332 * inv(op(A)) * diag(W), 00333 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) 00334 * 00335 DO 200 I = 1, N 00336 IF( RWORK( I ).GT.SAFE2 ) THEN 00337 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) 00338 ELSE 00339 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) + 00340 $ SAFE1 00341 END IF 00342 200 CONTINUE 00343 * 00344 KASE = 0 00345 210 CONTINUE 00346 CALL CLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE ) 00347 IF( KASE.NE.0 ) THEN 00348 IF( KASE.EQ.1 ) THEN 00349 * 00350 * Multiply by diag(W)*inv(op(A)**H). 00351 * 00352 CALL CTRSV( UPLO, TRANST, DIAG, N, A, LDA, WORK, 1 ) 00353 DO 220 I = 1, N 00354 WORK( I ) = RWORK( I )*WORK( I ) 00355 220 CONTINUE 00356 ELSE 00357 * 00358 * Multiply by inv(op(A))*diag(W). 00359 * 00360 DO 230 I = 1, N 00361 WORK( I ) = RWORK( I )*WORK( I ) 00362 230 CONTINUE 00363 CALL CTRSV( UPLO, TRANSN, DIAG, N, A, LDA, WORK, 1 ) 00364 END IF 00365 GO TO 210 00366 END IF 00367 * 00368 * Normalize error. 00369 * 00370 LSTRES = ZERO 00371 DO 240 I = 1, N 00372 LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) ) 00373 240 CONTINUE 00374 IF( LSTRES.NE.ZERO ) 00375 $ FERR( J ) = FERR( J ) / LSTRES 00376 * 00377 250 CONTINUE 00378 * 00379 RETURN 00380 * 00381 * End of CTRRFS 00382 * 00383 END
1.7.3 
