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