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