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