![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLA_SYRFSX_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_SYRFSX_EXTENDED + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cla_syrfsx_extended.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cla_syrfsx_extended.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cla_syrfsx_extended.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLA_SYRFSX_EXTENDED( PREC_TYPE, UPLO, N, NRHS, A, LDA, 00022 * AF, LDAF, IPIV, COLEQU, C, B, LDB, 00023 * Y, LDY, BERR_OUT, N_NORMS, 00024 * ERR_BNDS_NORM, ERR_BNDS_COMP, RES, 00025 * AYB, DY, Y_TAIL, RCOND, ITHRESH, 00026 * RTHRESH, DZ_UB, IGNORE_CWISE, 00027 * INFO ) 00028 * 00029 * .. Scalar Arguments .. 00030 * INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE, 00031 * $ N_NORMS, ITHRESH 00032 * CHARACTER UPLO 00033 * LOGICAL COLEQU, IGNORE_CWISE 00034 * REAL RTHRESH, DZ_UB 00035 * .. 00036 * .. Array Arguments .. 00037 * INTEGER IPIV( * ) 00038 * COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00039 * $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * ) 00040 * REAL C( * ), AYB( * ), RCOND, BERR_OUT( * ), 00041 * $ ERR_BNDS_NORM( NRHS, * ), 00042 * $ ERR_BNDS_COMP( NRHS, * ) 00043 * .. 00044 * 00045 * 00046 *> \par Purpose: 00047 * ============= 00048 *> 00049 *> \verbatim 00050 *> 00051 *> CLA_SYRFSX_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 CSYRFSX 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] UPLO 00077 *> \verbatim 00078 *> UPLO is CHARACTER*1 00079 *> = 'U': Upper triangle of A is stored; 00080 *> = 'L': Lower triangle of A is stored. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] N 00084 *> \verbatim 00085 *> N is INTEGER 00086 *> The number of linear equations, i.e., the order of the 00087 *> matrix A. N >= 0. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] NRHS 00091 *> \verbatim 00092 *> NRHS is INTEGER 00093 *> The number of right-hand-sides, i.e., the number of columns of the 00094 *> matrix B. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] A 00098 *> \verbatim 00099 *> A is COMPLEX array, dimension (LDA,N) 00100 *> On entry, the N-by-N matrix A. 00101 *> \endverbatim 00102 *> 00103 *> \param[in] LDA 00104 *> \verbatim 00105 *> LDA is INTEGER 00106 *> The leading dimension of the array A. LDA >= max(1,N). 00107 *> \endverbatim 00108 *> 00109 *> \param[in] AF 00110 *> \verbatim 00111 *> AF is COMPLEX array, dimension (LDAF,N) 00112 *> The block diagonal matrix D and the multipliers used to 00113 *> obtain the factor U or L as computed by CSYTRF. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] LDAF 00117 *> \verbatim 00118 *> LDAF is INTEGER 00119 *> The leading dimension of the array AF. LDAF >= max(1,N). 00120 *> \endverbatim 00121 *> 00122 *> \param[in] IPIV 00123 *> \verbatim 00124 *> IPIV is INTEGER array, dimension (N) 00125 *> Details of the interchanges and the block structure of D 00126 *> as determined by CSYTRF. 00127 *> \endverbatim 00128 *> 00129 *> \param[in] COLEQU 00130 *> \verbatim 00131 *> COLEQU is LOGICAL 00132 *> If .TRUE. then column equilibration was done to A before calling 00133 *> this routine. This is needed to compute the solution and error 00134 *> bounds correctly. 00135 *> \endverbatim 00136 *> 00137 *> \param[in] C 00138 *> \verbatim 00139 *> C is REAL array, dimension (N) 00140 *> The column scale factors for A. If COLEQU = .FALSE., C 00141 *> is not accessed. If C is input, each element of C should be a power 00142 *> of the radix to ensure a reliable solution and error estimates. 00143 *> Scaling by powers of the radix does not cause rounding errors unless 00144 *> the result underflows or overflows. Rounding errors during scaling 00145 *> lead to refining with a matrix that is not equivalent to the 00146 *> input matrix, producing error estimates that may not be 00147 *> reliable. 00148 *> \endverbatim 00149 *> 00150 *> \param[in] B 00151 *> \verbatim 00152 *> B is COMPLEX array, dimension (LDB,NRHS) 00153 *> The right-hand-side matrix B. 00154 *> \endverbatim 00155 *> 00156 *> \param[in] LDB 00157 *> \verbatim 00158 *> LDB is INTEGER 00159 *> The leading dimension of the array B. LDB >= max(1,N). 00160 *> \endverbatim 00161 *> 00162 *> \param[in,out] Y 00163 *> \verbatim 00164 *> Y is COMPLEX array, dimension 00165 *> (LDY,NRHS) 00166 *> On entry, the solution matrix X, as computed by CSYTRS. 00167 *> On exit, the improved solution matrix Y. 00168 *> \endverbatim 00169 *> 00170 *> \param[in] LDY 00171 *> \verbatim 00172 *> LDY is INTEGER 00173 *> The leading dimension of the array Y. LDY >= max(1,N). 00174 *> \endverbatim 00175 *> 00176 *> \param[out] BERR_OUT 00177 *> \verbatim 00178 *> BERR_OUT is REAL array, dimension (NRHS) 00179 *> On exit, BERR_OUT(j) contains the componentwise relative backward 00180 *> error for right-hand-side j from the formula 00181 *> max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) 00182 *> where abs(Z) is the componentwise absolute value of the matrix 00183 *> or vector Z. This is computed by CLA_LIN_BERR. 00184 *> \endverbatim 00185 *> 00186 *> \param[in] N_NORMS 00187 *> \verbatim 00188 *> N_NORMS is INTEGER 00189 *> Determines which error bounds to return (see ERR_BNDS_NORM 00190 *> and ERR_BNDS_COMP). 00191 *> If N_NORMS >= 1 return normwise error bounds. 00192 *> If N_NORMS >= 2 return componentwise error bounds. 00193 *> \endverbatim 00194 *> 00195 *> \param[in,out] ERR_BNDS_NORM 00196 *> \verbatim 00197 *> ERR_BNDS_NORM is REAL array, dimension 00198 *> (NRHS, N_ERR_BNDS) 00199 *> For each right-hand side, this array contains information about 00200 *> various error bounds and condition numbers corresponding to the 00201 *> normwise relative error, which is defined as follows: 00202 *> 00203 *> Normwise relative error in the ith solution vector: 00204 *> max_j (abs(XTRUE(j,i) - X(j,i))) 00205 *> ------------------------------ 00206 *> max_j abs(X(j,i)) 00207 *> 00208 *> The array is indexed by the type of error information as described 00209 *> below. There currently are up to three pieces of information 00210 *> returned. 00211 *> 00212 *> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith 00213 *> right-hand side. 00214 *> 00215 *> The second index in ERR_BNDS_NORM(:,err) contains the following 00216 *> three fields: 00217 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the 00218 *> reciprocal condition number is less than the threshold 00219 *> sqrt(n) * slamch('Epsilon'). 00220 *> 00221 *> err = 2 "Guaranteed" error bound: The estimated forward error, 00222 *> almost certainly within a factor of 10 of the true error 00223 *> so long as the next entry is greater than the threshold 00224 *> sqrt(n) * slamch('Epsilon'). This error bound should only 00225 *> be trusted if the previous boolean is true. 00226 *> 00227 *> err = 3 Reciprocal condition number: Estimated normwise 00228 *> reciprocal condition number. Compared with the threshold 00229 *> sqrt(n) * slamch('Epsilon') to determine if the error 00230 *> estimate is "guaranteed". These reciprocal condition 00231 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 00232 *> appropriately scaled matrix Z. 00233 *> Let Z = S*A, where S scales each row by a power of the 00234 *> radix so all absolute row sums of Z are approximately 1. 00235 *> 00236 *> This subroutine is only responsible for setting the second field 00237 *> above. 00238 *> See Lapack Working Note 165 for further details and extra 00239 *> cautions. 00240 *> \endverbatim 00241 *> 00242 *> \param[in,out] ERR_BNDS_COMP 00243 *> \verbatim 00244 *> ERR_BNDS_COMP is REAL array, dimension 00245 *> (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 *> ERR_BNDS_COMP 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 ERR_BNDS_COMP(i,:) corresponds to the ith 00264 *> right-hand side. 00265 *> 00266 *> The second index in ERR_BNDS_COMP(:,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 *> ERR_BNDS_NORM and ERR_BNDS_COMP 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 CLA_SYRFSX_EXTENDED 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 complexSYcomputational 00391 * 00392 * ===================================================================== 00393 SUBROUTINE CLA_SYRFSX_EXTENDED( PREC_TYPE, UPLO, N, NRHS, A, LDA, 00394 $ AF, LDAF, IPIV, COLEQU, C, B, LDB, 00395 $ Y, LDY, BERR_OUT, N_NORMS, 00396 $ ERR_BNDS_NORM, ERR_BNDS_COMP, RES, 00397 $ AYB, DY, Y_TAIL, RCOND, ITHRESH, 00398 $ RTHRESH, DZ_UB, IGNORE_CWISE, 00399 $ INFO ) 00400 * 00401 * -- LAPACK computational routine (version 3.4.0) -- 00402 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00403 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00404 * November 2011 00405 * 00406 * .. Scalar Arguments .. 00407 INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE, 00408 $ N_NORMS, ITHRESH 00409 CHARACTER UPLO 00410 LOGICAL COLEQU, IGNORE_CWISE 00411 REAL RTHRESH, DZ_UB 00412 * .. 00413 * .. Array Arguments .. 00414 INTEGER IPIV( * ) 00415 COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00416 $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * ) 00417 REAL C( * ), AYB( * ), RCOND, BERR_OUT( * ), 00418 $ ERR_BNDS_NORM( NRHS, * ), 00419 $ ERR_BNDS_COMP( NRHS, * ) 00420 * .. 00421 * 00422 * ===================================================================== 00423 * 00424 * .. Local Scalars .. 00425 INTEGER UPLO2, CNT, I, J, X_STATE, Z_STATE, 00426 $ Y_PREC_STATE 00427 REAL YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT, 00428 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX, 00429 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z, 00430 $ EPS, HUGEVAL, INCR_THRESH 00431 LOGICAL INCR_PREC, UPPER 00432 COMPLEX ZDUM 00433 * .. 00434 * .. Parameters .. 00435 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE, 00436 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL, 00437 $ EXTRA_Y 00438 PARAMETER ( UNSTABLE_STATE = 0, WORKING_STATE = 1, 00439 $ CONV_STATE = 2, NOPROG_STATE = 3 ) 00440 PARAMETER ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1, 00441 $ EXTRA_Y = 2 ) 00442 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 00443 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 00444 INTEGER CMP_ERR_I, PIV_GROWTH_I 00445 PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 00446 $ BERR_I = 3 ) 00447 PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 00448 PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8, 00449 $ PIV_GROWTH_I = 9 ) 00450 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I, 00451 $ LA_LINRX_CWISE_I 00452 PARAMETER ( LA_LINRX_ITREF_I = 1, 00453 $ LA_LINRX_ITHRESH_I = 2 ) 00454 PARAMETER ( LA_LINRX_CWISE_I = 3 ) 00455 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I, 00456 $ LA_LINRX_RCOND_I 00457 PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 ) 00458 PARAMETER ( LA_LINRX_RCOND_I = 3 ) 00459 * .. 00460 * .. External Functions .. 00461 LOGICAL LSAME 00462 EXTERNAL ILAUPLO 00463 INTEGER ILAUPLO 00464 * .. 00465 * .. External Subroutines .. 00466 EXTERNAL CAXPY, CCOPY, CSYTRS, CSYMV, BLAS_CSYMV_X, 00467 $ BLAS_CSYMV2_X, CLA_SYAMV, CLA_WWADDW, 00468 $ CLA_LIN_BERR 00469 REAL SLAMCH 00470 * .. 00471 * .. Intrinsic Functions .. 00472 INTRINSIC ABS, REAL, AIMAG, MAX, MIN 00473 * .. 00474 * .. Statement Functions .. 00475 REAL CABS1 00476 * .. 00477 * .. Statement Function Definitions .. 00478 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 00479 * .. 00480 * .. Executable Statements .. 00481 * 00482 INFO = 0 00483 UPPER = LSAME( UPLO, 'U' ) 00484 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00485 INFO = -2 00486 ELSE IF( N.LT.0 ) THEN 00487 INFO = -3 00488 ELSE IF( NRHS.LT.0 ) THEN 00489 INFO = -4 00490 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00491 INFO = -6 00492 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00493 INFO = -8 00494 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00495 INFO = -13 00496 ELSE IF( LDY.LT.MAX( 1, N ) ) THEN 00497 INFO = -15 00498 END IF 00499 IF( INFO.NE.0 ) THEN 00500 CALL XERBLA( 'CLA_SYRFSX_EXTENDED', -INFO ) 00501 RETURN 00502 END IF 00503 EPS = SLAMCH( 'Epsilon' ) 00504 HUGEVAL = SLAMCH( 'Overflow' ) 00505 * Force HUGEVAL to Inf 00506 HUGEVAL = HUGEVAL * HUGEVAL 00507 * Using HUGEVAL may lead to spurious underflows. 00508 INCR_THRESH = REAL( N ) * EPS 00509 00510 IF ( LSAME ( UPLO, 'L' ) ) THEN 00511 UPLO2 = ILAUPLO( 'L' ) 00512 ELSE 00513 UPLO2 = ILAUPLO( 'U' ) 00514 ENDIF 00515 00516 DO J = 1, NRHS 00517 Y_PREC_STATE = EXTRA_RESIDUAL 00518 IF ( Y_PREC_STATE .EQ. EXTRA_Y ) THEN 00519 DO I = 1, N 00520 Y_TAIL( I ) = 0.0 00521 END DO 00522 END IF 00523 00524 DXRAT = 0.0 00525 DXRATMAX = 0.0 00526 DZRAT = 0.0 00527 DZRATMAX = 0.0 00528 FINAL_DX_X = HUGEVAL 00529 FINAL_DZ_Z = HUGEVAL 00530 PREVNORMDX = HUGEVAL 00531 PREV_DZ_Z = HUGEVAL 00532 DZ_Z = HUGEVAL 00533 DX_X = HUGEVAL 00534 00535 X_STATE = WORKING_STATE 00536 Z_STATE = UNSTABLE_STATE 00537 INCR_PREC = .FALSE. 00538 00539 DO CNT = 1, ITHRESH 00540 * 00541 * Compute residual RES = B_s - op(A_s) * Y, 00542 * op(A) = A, A**T, or A**H depending on TRANS (and type). 00543 * 00544 CALL CCOPY( N, B( 1, J ), 1, RES, 1 ) 00545 IF ( Y_PREC_STATE .EQ. BASE_RESIDUAL ) THEN 00546 CALL CSYMV( UPLO, N, CMPLX(-1.0), A, LDA, Y(1,J), 1, 00547 $ CMPLX(1.0), RES, 1 ) 00548 ELSE IF ( Y_PREC_STATE .EQ. EXTRA_RESIDUAL ) THEN 00549 CALL BLAS_CSYMV_X( UPLO2, N, CMPLX(-1.0), A, LDA, 00550 $ Y( 1, J ), 1, CMPLX(1.0), RES, 1, PREC_TYPE ) 00551 ELSE 00552 CALL BLAS_CSYMV2_X(UPLO2, N, CMPLX(-1.0), A, LDA, 00553 $ Y(1, J), Y_TAIL, 1, CMPLX(1.0), RES, 1, PREC_TYPE) 00554 END IF 00555 00556 ! XXX: RES is no longer needed. 00557 CALL CCOPY( N, RES, 1, DY, 1 ) 00558 CALL CSYTRS( UPLO, N, 1, AF, LDAF, IPIV, DY, N, INFO ) 00559 * 00560 * Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT. 00561 * 00562 NORMX = 0.0 00563 NORMY = 0.0 00564 NORMDX = 0.0 00565 DZ_Z = 0.0 00566 YMIN = HUGEVAL 00567 00568 DO I = 1, N 00569 YK = CABS1( Y( I, J ) ) 00570 DYK = CABS1( DY( I ) ) 00571 00572 IF ( YK .NE. 0.0 ) THEN 00573 DZ_Z = MAX( DZ_Z, DYK / YK ) 00574 ELSE IF ( DYK .NE. 0.0 ) THEN 00575 DZ_Z = HUGEVAL 00576 END IF 00577 00578 YMIN = MIN( YMIN, YK ) 00579 00580 NORMY = MAX( NORMY, YK ) 00581 00582 IF ( COLEQU ) THEN 00583 NORMX = MAX( NORMX, YK * C( I ) ) 00584 NORMDX = MAX( NORMDX, DYK * C( I ) ) 00585 ELSE 00586 NORMX = NORMY 00587 NORMDX = MAX( NORMDX, DYK ) 00588 END IF 00589 END DO 00590 00591 IF ( NORMX .NE. 0.0 ) THEN 00592 DX_X = NORMDX / NORMX 00593 ELSE IF ( NORMDX .EQ. 0.0 ) THEN 00594 DX_X = 0.0 00595 ELSE 00596 DX_X = HUGEVAL 00597 END IF 00598 00599 DXRAT = NORMDX / PREVNORMDX 00600 DZRAT = DZ_Z / PREV_DZ_Z 00601 * 00602 * Check termination criteria. 00603 * 00604 IF ( YMIN*RCOND .LT. INCR_THRESH*NORMY 00605 $ .AND. Y_PREC_STATE .LT. EXTRA_Y ) 00606 $ INCR_PREC = .TRUE. 00607 00608 IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH ) 00609 $ X_STATE = WORKING_STATE 00610 IF ( X_STATE .EQ. WORKING_STATE ) THEN 00611 IF ( DX_X .LE. EPS ) THEN 00612 X_STATE = CONV_STATE 00613 ELSE IF ( DXRAT .GT. RTHRESH ) THEN 00614 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 00615 INCR_PREC = .TRUE. 00616 ELSE 00617 X_STATE = NOPROG_STATE 00618 END IF 00619 ELSE 00620 IF (DXRAT .GT. DXRATMAX) DXRATMAX = DXRAT 00621 END IF 00622 IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X 00623 END IF 00624 00625 IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB ) 00626 $ Z_STATE = WORKING_STATE 00627 IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH ) 00628 $ Z_STATE = WORKING_STATE 00629 IF ( Z_STATE .EQ. WORKING_STATE ) THEN 00630 IF ( DZ_Z .LE. EPS ) THEN 00631 Z_STATE = CONV_STATE 00632 ELSE IF ( DZ_Z .GT. DZ_UB ) THEN 00633 Z_STATE = UNSTABLE_STATE 00634 DZRATMAX = 0.0 00635 FINAL_DZ_Z = HUGEVAL 00636 ELSE IF ( DZRAT .GT. RTHRESH ) THEN 00637 IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN 00638 INCR_PREC = .TRUE. 00639 ELSE 00640 Z_STATE = NOPROG_STATE 00641 END IF 00642 ELSE 00643 IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT 00644 END IF 00645 IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 00646 END IF 00647 00648 IF ( X_STATE.NE.WORKING_STATE.AND. 00649 $ ( IGNORE_CWISE.OR.Z_STATE.NE.WORKING_STATE ) ) 00650 $ GOTO 666 00651 00652 IF ( INCR_PREC ) THEN 00653 INCR_PREC = .FALSE. 00654 Y_PREC_STATE = Y_PREC_STATE + 1 00655 DO I = 1, N 00656 Y_TAIL( I ) = 0.0 00657 END DO 00658 END IF 00659 00660 PREVNORMDX = NORMDX 00661 PREV_DZ_Z = DZ_Z 00662 * 00663 * Update soluton. 00664 * 00665 IF ( Y_PREC_STATE .LT. EXTRA_Y ) THEN 00666 CALL CAXPY( N, CMPLX(1.0), DY, 1, Y(1,J), 1 ) 00667 ELSE 00668 CALL CLA_WWADDW( N, Y(1,J), Y_TAIL, DY ) 00669 END IF 00670 00671 END DO 00672 * Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT. 00673 666 CONTINUE 00674 * 00675 * Set final_* when cnt hits ithresh. 00676 * 00677 IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X 00678 IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z 00679 * 00680 * Compute error bounds. 00681 * 00682 IF ( N_NORMS .GE. 1 ) THEN 00683 ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 00684 $ FINAL_DX_X / (1 - DXRATMAX) 00685 END IF 00686 IF ( N_NORMS .GE. 2 ) THEN 00687 ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 00688 $ FINAL_DZ_Z / (1 - DZRATMAX) 00689 END IF 00690 * 00691 * Compute componentwise relative backward error from formula 00692 * max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) 00693 * where abs(Z) is the componentwise absolute value of the matrix 00694 * or vector Z. 00695 * 00696 * Compute residual RES = B_s - op(A_s) * Y, 00697 * op(A) = A, A**T, or A**H depending on TRANS (and type). 00698 * 00699 CALL CCOPY( N, B( 1, J ), 1, RES, 1 ) 00700 CALL CSYMV( UPLO, N, CMPLX(-1.0), A, LDA, Y(1,J), 1, 00701 $ CMPLX(1.0), RES, 1 ) 00702 00703 DO I = 1, N 00704 AYB( I ) = CABS1( B( I, J ) ) 00705 END DO 00706 * 00707 * Compute abs(op(A_s))*abs(Y) + abs(B_s). 00708 * 00709 CALL CLA_SYAMV ( UPLO2, N, 1.0, 00710 $ A, LDA, Y(1, J), 1, 1.0, AYB, 1 ) 00711 00712 CALL CLA_LIN_BERR ( N, N, 1, RES, AYB, BERR_OUT( J ) ) 00713 * 00714 * End of loop for each RHS. 00715 * 00716 END DO 00717 * 00718 RETURN 00719 END