![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLA_GBRFSX_EXTENDED 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLA_GBRFSX_EXTENDED + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_gbrfsx_extended.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_gbrfsx_extended.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_gbrfsx_extended.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU, 00022 * NRHS, AB, LDAB, AFB, LDAFB, IPIV, 00023 * COLEQU, C, B, LDB, Y, LDY, 00024 * BERR_OUT, N_NORMS, ERR_BNDS_NORM, 00025 * ERR_BNDS_COMP, RES, AYB, DY, 00026 * Y_TAIL, RCOND, ITHRESH, RTHRESH, 00027 * DZ_UB, IGNORE_CWISE, INFO ) 00028 * 00029 * .. Scalar Arguments .. 00030 * INTEGER INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS, 00031 * $ PREC_TYPE, TRANS_TYPE, N_NORMS, ITHRESH 00032 * LOGICAL COLEQU, IGNORE_CWISE 00033 * DOUBLE PRECISION RTHRESH, DZ_UB 00034 * .. 00035 * .. Array Arguments .. 00036 * INTEGER IPIV( * ) 00037 * DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 00038 * $ Y( LDY, * ), RES(*), DY(*), Y_TAIL(*) 00039 * DOUBLE PRECISION C( * ), AYB(*), RCOND, BERR_OUT(*), 00040 * $ ERR_BNDS_NORM( NRHS, * ), 00041 * $ ERR_BNDS_COMP( NRHS, * ) 00042 * .. 00043 * 00044 * 00045 *> \par Purpose: 00046 * ============= 00047 *> 00048 *> \verbatim 00049 *> 00050 *> 00051 *> DLA_GBRFSX_EXTENDED improves the computed solution to a system of 00052 *> linear equations by performing extra-precise iterative refinement 00053 *> and provides error bounds and backward error estimates for the solution. 00054 *> This subroutine is called by DGBRFSX to perform iterative refinement. 00055 *> In addition to normwise error bound, the code provides maximum 00056 *> componentwise error bound if possible. See comments for ERR_BNDS_NORM 00057 *> and ERR_BNDS_COMP for details of the error bounds. Note that this 00058 *> subroutine is only resonsible for setting the second fields of 00059 *> ERR_BNDS_NORM and ERR_BNDS_COMP. 00060 *> \endverbatim 00061 * 00062 * Arguments: 00063 * ========== 00064 * 00065 *> \param[in] PREC_TYPE 00066 *> \verbatim 00067 *> PREC_TYPE is INTEGER 00068 *> Specifies the intermediate precision to be used in refinement. 00069 *> The value is defined by ILAPREC(P) where P is a CHARACTER and 00070 *> P = 'S': Single 00071 *> = 'D': Double 00072 *> = 'I': Indigenous 00073 *> = 'X', 'E': Extra 00074 *> \endverbatim 00075 *> 00076 *> \param[in] TRANS_TYPE 00077 *> \verbatim 00078 *> TRANS_TYPE is INTEGER 00079 *> Specifies the transposition operation on A. 00080 *> The value is defined by ILATRANS(T) where T is a CHARACTER and 00081 *> T = 'N': No transpose 00082 *> = 'T': Transpose 00083 *> = 'C': Conjugate transpose 00084 *> \endverbatim 00085 *> 00086 *> \param[in] N 00087 *> \verbatim 00088 *> N is INTEGER 00089 *> The number of linear equations, i.e., the order of the 00090 *> matrix A. N >= 0. 00091 *> \endverbatim 00092 *> 00093 *> \param[in] KL 00094 *> \verbatim 00095 *> KL is INTEGER 00096 *> The number of subdiagonals within the band of A. KL >= 0. 00097 *> \endverbatim 00098 *> 00099 *> \param[in] KU 00100 *> \verbatim 00101 *> KU is INTEGER 00102 *> The number of superdiagonals within the band of A. KU >= 0 00103 *> \endverbatim 00104 *> 00105 *> \param[in] NRHS 00106 *> \verbatim 00107 *> NRHS is INTEGER 00108 *> The number of right-hand-sides, i.e., the number of columns of the 00109 *> matrix B. 00110 *> \endverbatim 00111 *> 00112 *> \param[in] AB 00113 *> \verbatim 00114 *> AB is DOUBLE PRECISION array, dimension (LDAB,N) 00115 *> On entry, the N-by-N matrix AB. 00116 *> \endverbatim 00117 *> 00118 *> \param[in] LDAB 00119 *> \verbatim 00120 *> LDAB is INTEGER 00121 *> The leading dimension of the array AB. LDBA >= max(1,N). 00122 *> \endverbatim 00123 *> 00124 *> \param[in] AFB 00125 *> \verbatim 00126 *> AFB is DOUBLE PRECISION array, dimension (LDAFB,N) 00127 *> The factors L and U from the factorization 00128 *> A = P*L*U as computed by DGBTRF. 00129 *> \endverbatim 00130 *> 00131 *> \param[in] LDAFB 00132 *> \verbatim 00133 *> LDAFB is INTEGER 00134 *> The leading dimension of the array AF. LDAFB >= max(1,N). 00135 *> \endverbatim 00136 *> 00137 *> \param[in] IPIV 00138 *> \verbatim 00139 *> IPIV is INTEGER array, dimension (N) 00140 *> The pivot indices from the factorization A = P*L*U 00141 *> as computed by DGBTRF; row i of the matrix was interchanged 00142 *> with row IPIV(i). 00143 *> \endverbatim 00144 *> 00145 *> \param[in] COLEQU 00146 *> \verbatim 00147 *> COLEQU is LOGICAL 00148 *> If .TRUE. then column equilibration was done to A before calling 00149 *> this routine. This is needed to compute the solution and error 00150 *> bounds correctly. 00151 *> \endverbatim 00152 *> 00153 *> \param[in] C 00154 *> \verbatim 00155 *> C is DOUBLE PRECISION array, dimension (N) 00156 *> The column scale factors for A. If COLEQU = .FALSE., C 00157 *> is not accessed. If C is input, each element of C should be a power 00158 *> of the radix to ensure a reliable solution and error estimates. 00159 *> Scaling by powers of the radix does not cause rounding errors unless 00160 *> the result underflows or overflows. Rounding errors during scaling 00161 *> lead to refining with a matrix that is not equivalent to the 00162 *> input matrix, producing error estimates that may not be 00163 *> reliable. 00164 *> \endverbatim 00165 *> 00166 *> \param[in] B 00167 *> \verbatim 00168 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00169 *> The right-hand-side matrix B. 00170 *> \endverbatim 00171 *> 00172 *> \param[in] LDB 00173 *> \verbatim 00174 *> LDB is INTEGER 00175 *> The leading dimension of the array B. LDB >= max(1,N). 00176 *> \endverbatim 00177 *> 00178 *> \param[in,out] Y 00179 *> \verbatim 00180 *> Y is DOUBLE PRECISION array, dimension 00181 *> (LDY,NRHS) 00182 *> On entry, the solution matrix X, as computed by DGBTRS. 00183 *> On exit, the improved solution matrix Y. 00184 *> \endverbatim 00185 *> 00186 *> \param[in] LDY 00187 *> \verbatim 00188 *> LDY is INTEGER 00189 *> The leading dimension of the array Y. LDY >= max(1,N). 00190 *> \endverbatim 00191 *> 00192 *> \param[out] BERR_OUT 00193 *> \verbatim 00194 *> BERR_OUT is DOUBLE PRECISION array, dimension (NRHS) 00195 *> On exit, BERR_OUT(j) contains the componentwise relative backward 00196 *> error for right-hand-side j from the formula 00197 *> max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) 00198 *> where abs(Z) is the componentwise absolute value of the matrix 00199 *> or vector Z. This is computed by DLA_LIN_BERR. 00200 *> \endverbatim 00201 *> 00202 *> \param[in] N_NORMS 00203 *> \verbatim 00204 *> N_NORMS is INTEGER 00205 *> Determines which error bounds to return (see ERR_BNDS_NORM 00206 *> and ERR_BNDS_COMP). 00207 *> If N_NORMS >= 1 return normwise error bounds. 00208 *> If N_NORMS >= 2 return componentwise error bounds. 00209 *> \endverbatim 00210 *> 00211 *> \param[in,out] ERR_BNDS_NORM 00212 *> \verbatim 00213 *> ERR_BNDS_NORM is DOUBLE PRECISION array, dimension 00214 *> (NRHS, N_ERR_BNDS) 00215 *> For each right-hand side, this array contains information about 00216 *> various error bounds and condition numbers corresponding to the 00217 *> normwise relative error, which is defined as follows: 00218 *> 00219 *> Normwise relative error in the ith solution vector: 00220 *> max_j (abs(XTRUE(j,i) - X(j,i))) 00221 *> ------------------------------ 00222 *> max_j abs(X(j,i)) 00223 *> 00224 *> The array is indexed by the type of error information as described 00225 *> below. There currently are up to three pieces of information 00226 *> returned. 00227 *> 00228 *> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith 00229 *> right-hand side. 00230 *> 00231 *> The second index in ERR_BNDS_NORM(:,err) contains the following 00232 *> three fields: 00233 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the 00234 *> reciprocal condition number is less than the threshold 00235 *> sqrt(n) * slamch('Epsilon'). 00236 *> 00237 *> err = 2 "Guaranteed" error bound: The estimated forward error, 00238 *> almost certainly within a factor of 10 of the true error 00239 *> so long as the next entry is greater than the threshold 00240 *> sqrt(n) * slamch('Epsilon'). This error bound should only 00241 *> be trusted if the previous boolean is true. 00242 *> 00243 *> err = 3 Reciprocal condition number: Estimated normwise 00244 *> reciprocal condition number. Compared with the threshold 00245 *> sqrt(n) * slamch('Epsilon') to determine if the error 00246 *> estimate is "guaranteed". These reciprocal condition 00247 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 00248 *> appropriately scaled matrix Z. 00249 *> Let Z = S*A, where S scales each row by a power of the 00250 *> radix so all absolute row sums of Z are approximately 1. 00251 *> 00252 *> This subroutine is only responsible for setting the second field 00253 *> above. 00254 *> See Lapack Working Note 165 for further details and extra 00255 *> cautions. 00256 *> \endverbatim 00257 *> 00258 *> \param[in,out] ERR_BNDS_COMP 00259 *> \verbatim 00260 *> ERR_BNDS_COMP is DOUBLE PRECISION array, dimension 00261 *> (NRHS, N_ERR_BNDS) 00262 *> For each right-hand side, this array contains information about 00263 *> various error bounds and condition numbers corresponding to the 00264 *> componentwise relative error, which is defined as follows: 00265 *> 00266 *> Componentwise relative error in the ith solution vector: 00267 *> abs(XTRUE(j,i) - X(j,i)) 00268 *> max_j ---------------------- 00269 *> abs(X(j,i)) 00270 *> 00271 *> The array is indexed by the right-hand side i (on which the 00272 *> componentwise relative error depends), and the type of error 00273 *> information as described below. There currently are up to three 00274 *> pieces of information returned for each right-hand side. If 00275 *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then 00276 *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most 00277 *> the first (:,N_ERR_BNDS) entries are returned. 00278 *> 00279 *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith 00280 *> right-hand side. 00281 *> 00282 *> The second index in ERR_BNDS_COMP(:,err) contains the following 00283 *> three fields: 00284 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the 00285 *> reciprocal condition number is less than the threshold 00286 *> sqrt(n) * slamch('Epsilon'). 00287 *> 00288 *> err = 2 "Guaranteed" error bound: The estimated forward error, 00289 *> almost certainly within a factor of 10 of the true error 00290 *> so long as the next entry is greater than the threshold 00291 *> sqrt(n) * slamch('Epsilon'). This error bound should only 00292 *> be trusted if the previous boolean is true. 00293 *> 00294 *> err = 3 Reciprocal condition number: Estimated componentwise 00295 *> reciprocal condition number. Compared with the threshold 00296 *> sqrt(n) * slamch('Epsilon') to determine if the error 00297 *> estimate is "guaranteed". These reciprocal condition 00298 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 00299 *> appropriately scaled matrix Z. 00300 *> Let Z = S*(A*diag(x)), where x is the solution for the 00301 *> current right-hand side and S scales each row of 00302 *> A*diag(x) by a power of the radix so all absolute row 00303 *> sums of Z are approximately 1. 00304 *> 00305 *> This subroutine is only responsible for setting the second field 00306 *> above. 00307 *> See Lapack Working Note 165 for further details and extra 00308 *> cautions. 00309 *> \endverbatim 00310 *> 00311 *> \param[in] RES 00312 *> \verbatim 00313 *> RES is DOUBLE PRECISION array, dimension (N) 00314 *> Workspace to hold the intermediate residual. 00315 *> \endverbatim 00316 *> 00317 *> \param[in] AYB 00318 *> \verbatim 00319 *> AYB is DOUBLE PRECISION array, dimension (N) 00320 *> Workspace. This can be the same workspace passed for Y_TAIL. 00321 *> \endverbatim 00322 *> 00323 *> \param[in] DY 00324 *> \verbatim 00325 *> DY is DOUBLE PRECISION array, dimension (N) 00326 *> Workspace to hold the intermediate solution. 00327 *> \endverbatim 00328 *> 00329 *> \param[in] Y_TAIL 00330 *> \verbatim 00331 *> Y_TAIL is DOUBLE PRECISION array, dimension (N) 00332 *> Workspace to hold the trailing bits of the intermediate solution. 00333 *> \endverbatim 00334 *> 00335 *> \param[in] RCOND 00336 *> \verbatim 00337 *> RCOND is DOUBLE PRECISION 00338 *> Reciprocal scaled condition number. This is an estimate of the 00339 *> reciprocal Skeel condition number of the matrix A after 00340 *> equilibration (if done). If this is less than the machine 00341 *> precision (in particular, if it is zero), the matrix is singular 00342 *> to working precision. Note that the error may still be small even 00343 *> if this number is very small and the matrix appears ill- 00344 *> conditioned. 00345 *> \endverbatim 00346 *> 00347 *> \param[in] ITHRESH 00348 *> \verbatim 00349 *> ITHRESH is INTEGER 00350 *> The maximum number of residual computations allowed for 00351 *> refinement. The default is 10. For 'aggressive' set to 100 to 00352 *> permit convergence using approximate factorizations or 00353 *> factorizations other than LU. If the factorization uses a 00354 *> technique other than Gaussian elimination, the guarantees in 00355 *> ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy. 00356 *> \endverbatim 00357 *> 00358 *> \param[in] RTHRESH 00359 *> \verbatim 00360 *> RTHRESH is DOUBLE PRECISION 00361 *> Determines when to stop refinement if the error estimate stops 00362 *> decreasing. Refinement will stop when the next solution no longer 00363 *> satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is 00364 *> the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The 00365 *> default value is 0.5. For 'aggressive' set to 0.9 to permit 00366 *> convergence on extremely ill-conditioned matrices. See LAWN 165 00367 *> for more details. 00368 *> \endverbatim 00369 *> 00370 *> \param[in] DZ_UB 00371 *> \verbatim 00372 *> DZ_UB is DOUBLE PRECISION 00373 *> Determines when to start considering componentwise convergence. 00374 *> Componentwise convergence is only considered after each component 00375 *> of the solution Y is stable, which we definte as the relative 00376 *> change in each component being less than DZ_UB. The default value 00377 *> is 0.25, requiring the first bit to be stable. See LAWN 165 for 00378 *> more details. 00379 *> \endverbatim 00380 *> 00381 *> \param[in] IGNORE_CWISE 00382 *> \verbatim 00383 *> IGNORE_CWISE is LOGICAL 00384 *> If .TRUE. then ignore componentwise convergence. Default value 00385 *> is .FALSE.. 00386 *> \endverbatim 00387 *> 00388 *> \param[out] INFO 00389 *> \verbatim 00390 *> INFO is INTEGER 00391 *> = 0: Successful exit. 00392 *> < 0: if INFO = -i, the ith argument to DGBTRS had an illegal 00393 *> value 00394 *> \endverbatim 00395 * 00396 * Authors: 00397 * ======== 00398 * 00399 *> \author Univ. of Tennessee 00400 *> \author Univ. of California Berkeley 00401 *> \author Univ. of Colorado Denver 00402 *> \author NAG Ltd. 00403 * 00404 *> \date November 2011 00405 * 00406 *> \ingroup doubleGBcomputational 00407 * 00408 * ===================================================================== 00409 SUBROUTINE DLA_GBRFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, KL, KU, 00410 $ NRHS, AB, LDAB, AFB, LDAFB, IPIV, 00411 $ COLEQU, C, B, LDB, Y, LDY, 00412 $ BERR_OUT, N_NORMS, ERR_BNDS_NORM, 00413 $ ERR_BNDS_COMP, RES, AYB, DY, 00414 $ Y_TAIL, RCOND, ITHRESH, RTHRESH, 00415 $ DZ_UB, IGNORE_CWISE, INFO ) 00416 * 00417 * -- LAPACK computational routine (version 3.4.0) -- 00418 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00419 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00420 * November 2011 00421 * 00422 * .. Scalar Arguments .. 00423 INTEGER INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS, 00424 $ PREC_TYPE, TRANS_TYPE, N_NORMS, ITHRESH 00425 LOGICAL COLEQU, IGNORE_CWISE 00426 DOUBLE PRECISION RTHRESH, DZ_UB 00427 * .. 00428 * .. Array Arguments .. 00429 INTEGER IPIV( * ) 00430 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ), 00431 $ Y( LDY, * ), RES(*), DY(*), Y_TAIL(*) 00432 DOUBLE PRECISION C( * ), AYB(*), RCOND, BERR_OUT(*), 00433 $ ERR_BNDS_NORM( NRHS, * ), 00434 $ ERR_BNDS_COMP( NRHS, * ) 00435 * .. 00436 * 00437 * ===================================================================== 00438 * 00439 * .. Local Scalars .. 00440 CHARACTER TRANS 00441 INTEGER CNT, I, J, M, X_STATE, Z_STATE, Y_PREC_STATE 00442 DOUBLE PRECISION YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT, 00443 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX, 00444 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z, 00445 $ EPS, HUGEVAL, INCR_THRESH 00446 LOGICAL INCR_PREC 00447 * .. 00448 * .. Parameters .. 00449 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE, 00450 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL, 00451 $ EXTRA_Y 00452 PARAMETER ( UNSTABLE_STATE = 0, WORKING_STATE = 1, 00453 $ CONV_STATE = 2, NOPROG_STATE = 3 ) 00454 PARAMETER ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1, 00455 $ EXTRA_Y = 2 ) 00456 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 00457 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 00458 INTEGER CMP_ERR_I, PIV_GROWTH_I 00459 PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 00460 $ BERR_I = 3 ) 00461 PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 00462 PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8, 00463 $ PIV_GROWTH_I = 9 ) 00464 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I, 00465 $ LA_LINRX_CWISE_I 00466 PARAMETER ( LA_LINRX_ITREF_I = 1, 00467 $ LA_LINRX_ITHRESH_I = 2 ) 00468 PARAMETER ( LA_LINRX_CWISE_I = 3 ) 00469 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I, 00470 $ LA_LINRX_RCOND_I 00471 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 ) 00472 PARAMETER ( LA_LINRX_RCOND_I = 3 ) 00473 * .. 00474 * .. External Subroutines .. 00475 EXTERNAL DAXPY, DCOPY, DGBTRS, DGBMV, BLAS_DGBMV_X, 00476 $ BLAS_DGBMV2_X, DLA_GBAMV, DLA_WWADDW, DLAMCH, 00477 $ CHLA_TRANSTYPE, DLA_LIN_BERR 00478 DOUBLE PRECISION DLAMCH 00479 CHARACTER CHLA_TRANSTYPE 00480 * .. 00481 * .. Intrinsic Functions .. 00482 INTRINSIC ABS, MAX, MIN 00483 * .. 00484 * .. Executable Statements .. 00485 * 00486 IF (INFO.NE.0) RETURN 00487 TRANS = CHLA_TRANSTYPE(TRANS_TYPE) 00488 EPS = DLAMCH( 'Epsilon' ) 00489 HUGEVAL = DLAMCH( 'Overflow' ) 00490 * Force HUGEVAL to Inf 00491 HUGEVAL = HUGEVAL * HUGEVAL 00492 * Using HUGEVAL may lead to spurious underflows. 00493 INCR_THRESH = DBLE( N ) * EPS 00494 M = KL+KU+1 00495 00496 DO J = 1, NRHS 00497 Y_PREC_STATE = EXTRA_RESIDUAL 00498 IF ( Y_PREC_STATE .EQ. EXTRA_Y ) THEN 00499 DO I = 1, N 00500 Y_TAIL( I ) = 0.0D+0 00501 END DO 00502 END IF 00503 00504 DXRAT = 0.0D+0 00505 DXRATMAX = 0.0D+0 00506 DZRAT = 0.0D+0 00507 DZRATMAX = 0.0D+0 00508 FINAL_DX_X = HUGEVAL 00509 FINAL_DZ_Z = HUGEVAL 00510 PREVNORMDX = HUGEVAL 00511 PREV_DZ_Z = HUGEVAL 00512 DZ_Z = HUGEVAL 00513 DX_X = HUGEVAL 00514 00515 X_STATE = WORKING_STATE 00516 Z_STATE = UNSTABLE_STATE 00517 INCR_PREC = .FALSE. 00518 00519 DO CNT = 1, ITHRESH 00520 * 00521 * Compute residual RES = B_s - op(A_s) * Y, 00522 * op(A) = A, A**T, or A**H depending on TRANS (and type). 00523 * 00524 CALL DCOPY( N, B( 1, J ), 1, RES, 1 ) 00525 IF ( Y_PREC_STATE .EQ. BASE_RESIDUAL ) THEN 00526 CALL DGBMV( TRANS, M, N, KL, KU, -1.0D+0, AB, LDAB, 00527 $ Y( 1, J ), 1, 1.0D+0, RES, 1 ) 00528 ELSE IF ( Y_PREC_STATE .EQ. EXTRA_RESIDUAL ) THEN 00529 CALL BLAS_DGBMV_X( TRANS_TYPE, N, N, KL, KU, 00530 $ -1.0D+0, AB, LDAB, Y( 1, J ), 1, 1.0D+0, RES, 1, 00531 $ PREC_TYPE ) 00532 ELSE 00533 CALL BLAS_DGBMV2_X( TRANS_TYPE, N, N, KL, KU, -1.0D+0, 00534 $ AB, LDAB, Y( 1, J ), Y_TAIL, 1, 1.0D+0, RES, 1, 00535 $ PREC_TYPE ) 00536 END IF 00537 00538 ! XXX: RES is no longer needed. 00539 CALL DCOPY( N, RES, 1, DY, 1 ) 00540 CALL DGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV, DY, N, 00541 $ INFO ) 00542 * 00543 * Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT. 00544 * 00545 NORMX = 0.0D+0 00546 NORMY = 0.0D+0 00547 NORMDX = 0.0D+0 00548 DZ_Z = 0.0D+0 00549 YMIN = HUGEVAL 00550 00551 DO I = 1, N 00552 YK = ABS( Y( I, J ) ) 00553 DYK = ABS( DY( I ) ) 00554 00555 IF ( YK .NE. 0.0D+0 ) THEN 00556 DZ_Z = MAX( DZ_Z, DYK / YK ) 00557 ELSE IF ( DYK .NE. 0.0D+0 ) THEN 00558 DZ_Z = HUGEVAL 00559 END IF 00560 00561 YMIN = MIN( YMIN, YK ) 00562 00563 NORMY = MAX( NORMY, YK ) 00564 00565 IF ( COLEQU ) THEN 00566 NORMX = MAX( NORMX, YK * C( I ) ) 00567 NORMDX = MAX( NORMDX, DYK * C( I ) ) 00568 ELSE 00569 NORMX = NORMY 00570 NORMDX = MAX( NORMDX, DYK ) 00571 END IF 00572 END DO 00573 00574 IF ( NORMX .NE. 0.0D+0 ) THEN 00575 DX_X = NORMDX / NORMX 00576 ELSE IF ( NORMDX .EQ. 0.0D+0 ) THEN 00577 DX_X = 0.0D+0 00578 ELSE 00579 DX_X = HUGEVAL 00580 END IF 00581 00582 DXRAT = NORMDX / PREVNORMDX 00583 DZRAT = DZ_Z / PREV_DZ_Z 00584 * 00585 * Check termination criteria. 00586 * 00587 IF ( .NOT.IGNORE_CWISE 00588 $ .AND. YMIN*RCOND .LT. INCR_THRESH*NORMY 00589 $ .AND. Y_PREC_STATE .LT. EXTRA_Y ) 00590 $ INCR_PREC = .TRUE. 00591 00592 IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH ) 00593 $ X_STATE = WORKING_STATE 00594 IF ( X_STATE .EQ. WORKING_STATE ) THEN 00595 IF ( DX_X .LE. EPS ) THEN 00596 X_STATE = CONV_STATE 00597 ELSE IF ( DXRAT .GT. RTHRESH ) THEN 00598 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 00599 INCR_PREC = .TRUE. 00600 ELSE 00601 X_STATE = NOPROG_STATE 00602 END IF 00603 ELSE 00604 IF ( DXRAT .GT. DXRATMAX ) DXRATMAX = DXRAT 00605 END IF 00606 IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X 00607 END IF 00608 00609 IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB ) 00610 $ Z_STATE = WORKING_STATE 00611 IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH ) 00612 $ Z_STATE = WORKING_STATE 00613 IF ( Z_STATE .EQ. WORKING_STATE ) THEN 00614 IF ( DZ_Z .LE. EPS ) THEN 00615 Z_STATE = CONV_STATE 00616 ELSE IF ( DZ_Z .GT. DZ_UB ) THEN 00617 Z_STATE = UNSTABLE_STATE 00618 DZRATMAX = 0.0D+0 00619 FINAL_DZ_Z = HUGEVAL 00620 ELSE IF ( DZRAT .GT. RTHRESH ) THEN 00621 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 00622 INCR_PREC = .TRUE. 00623 ELSE 00624 Z_STATE = NOPROG_STATE 00625 END IF 00626 ELSE 00627 IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT 00628 END IF 00629 IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 00630 END IF 00631 * 00632 * Exit if both normwise and componentwise stopped working, 00633 * but if componentwise is unstable, let it go at least two 00634 * iterations. 00635 * 00636 IF ( X_STATE.NE.WORKING_STATE ) THEN 00637 IF ( IGNORE_CWISE ) GOTO 666 00638 IF ( Z_STATE.EQ.NOPROG_STATE .OR. Z_STATE.EQ.CONV_STATE ) 00639 $ GOTO 666 00640 IF ( Z_STATE.EQ.UNSTABLE_STATE .AND. CNT.GT.1 ) GOTO 666 00641 END IF 00642 00643 IF ( INCR_PREC ) THEN 00644 INCR_PREC = .FALSE. 00645 Y_PREC_STATE = Y_PREC_STATE + 1 00646 DO I = 1, N 00647 Y_TAIL( I ) = 0.0D+0 00648 END DO 00649 END IF 00650 00651 PREVNORMDX = NORMDX 00652 PREV_DZ_Z = DZ_Z 00653 * 00654 * Update soluton. 00655 * 00656 IF (Y_PREC_STATE .LT. EXTRA_Y) THEN 00657 CALL DAXPY( N, 1.0D+0, DY, 1, Y(1,J), 1 ) 00658 ELSE 00659 CALL DLA_WWADDW( N, Y(1,J), Y_TAIL, DY ) 00660 END IF 00661 00662 END DO 00663 * Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT. 00664 666 CONTINUE 00665 * 00666 * Set final_* when cnt hits ithresh. 00667 * 00668 IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X 00669 IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 00670 * 00671 * Compute error bounds. 00672 * 00673 IF ( N_NORMS .GE. 1 ) THEN 00674 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 00675 $ FINAL_DX_X / (1 - DXRATMAX) 00676 END IF 00677 IF (N_NORMS .GE. 2) THEN 00678 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 00679 $ FINAL_DZ_Z / (1 - DZRATMAX) 00680 END IF 00681 * 00682 * Compute componentwise relative backward error from formula 00683 * max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) 00684 * where abs(Z) is the componentwise absolute value of the matrix 00685 * or vector Z. 00686 * 00687 * Compute residual RES = B_s - op(A_s) * Y, 00688 * op(A) = A, A**T, or A**H depending on TRANS (and type). 00689 * 00690 CALL DCOPY( N, B( 1, J ), 1, RES, 1 ) 00691 CALL DGBMV(TRANS, N, N, KL, KU, -1.0D+0, AB, LDAB, Y(1,J), 00692 $ 1, 1.0D+0, RES, 1 ) 00693 00694 DO I = 1, N 00695 AYB( I ) = ABS( B( I, J ) ) 00696 END DO 00697 * 00698 * Compute abs(op(A_s))*abs(Y) + abs(B_s). 00699 * 00700 CALL DLA_GBAMV( TRANS_TYPE, N, N, KL, KU, 1.0D+0, 00701 $ AB, LDAB, Y(1, J), 1, 1.0D+0, AYB, 1 ) 00702 00703 CALL DLA_LIN_BERR( N, N, 1, RES, AYB, BERR_OUT( J ) ) 00704 * 00705 * End of loop for each RHS 00706 * 00707 END DO 00708 * 00709 RETURN 00710 END