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