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