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