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