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