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