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