1*> \brief <b> DPOSVXX computes the solution to system of linear equations A * X = B for PO matrices</b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DPOSVXX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dposvxx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dposvxx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dposvxx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DPOSVXX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 22* S, B, LDB, X, LDX, RCOND, RPVGRW, BERR, 23* N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, 24* NPARAMS, PARAMS, WORK, IWORK, INFO ) 25* 26* .. Scalar Arguments .. 27* CHARACTER EQUED, FACT, UPLO 28* INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS, 29* $ N_ERR_BNDS 30* DOUBLE PRECISION RCOND, RPVGRW 31* .. 32* .. Array Arguments .. 33* INTEGER IWORK( * ) 34* DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 35* $ X( LDX, * ), WORK( * ) 36* DOUBLE PRECISION S( * ), PARAMS( * ), BERR( * ), 37* $ ERR_BNDS_NORM( NRHS, * ), 38* $ ERR_BNDS_COMP( NRHS, * ) 39* .. 40* 41* 42*> \par Purpose: 43* ============= 44*> 45*> \verbatim 46*> 47*> DPOSVXX uses the Cholesky factorization A = U**T*U or A = L*L**T 48*> to compute the solution to a double precision system of linear equations 49*> A * X = B, where A is an N-by-N symmetric positive definite matrix 50*> and X and B are N-by-NRHS matrices. 51*> 52*> If requested, both normwise and maximum componentwise error bounds 53*> are returned. DPOSVXX will return a solution with a tiny 54*> guaranteed error (O(eps) where eps is the working machine 55*> precision) unless the matrix is very ill-conditioned, in which 56*> case a warning is returned. Relevant condition numbers also are 57*> calculated and returned. 58*> 59*> DPOSVXX accepts user-provided factorizations and equilibration 60*> factors; see the definitions of the FACT and EQUED options. 61*> Solving with refinement and using a factorization from a previous 62*> DPOSVXX call will also produce a solution with either O(eps) 63*> errors or warnings, but we cannot make that claim for general 64*> user-provided factorizations and equilibration factors if they 65*> differ from what DPOSVXX would itself produce. 66*> \endverbatim 67* 68*> \par Description: 69* ================= 70*> 71*> \verbatim 72*> 73*> The following steps are performed: 74*> 75*> 1. If FACT = 'E', double precision scaling factors are computed to equilibrate 76*> the system: 77*> 78*> diag(S)*A*diag(S) *inv(diag(S))*X = diag(S)*B 79*> 80*> Whether or not the system will be equilibrated depends on the 81*> scaling of the matrix A, but if equilibration is used, A is 82*> overwritten by diag(S)*A*diag(S) and B by diag(S)*B. 83*> 84*> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to 85*> factor the matrix A (after equilibration if FACT = 'E') as 86*> A = U**T* U, if UPLO = 'U', or 87*> A = L * L**T, if UPLO = 'L', 88*> where U is an upper triangular matrix and L is a lower triangular 89*> matrix. 90*> 91*> 3. If the leading i-by-i principal minor is not positive definite, 92*> then the routine returns with INFO = i. Otherwise, the factored 93*> form of A is used to estimate the condition number of the matrix 94*> A (see argument RCOND). If the reciprocal of the condition number 95*> is less than machine precision, the routine still goes on to solve 96*> for X and compute error bounds as described below. 97*> 98*> 4. The system of equations is solved for X using the factored form 99*> of A. 100*> 101*> 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero), 102*> the routine will use iterative refinement to try to get a small 103*> error and error bounds. Refinement calculates the residual to at 104*> least twice the working precision. 105*> 106*> 6. If equilibration was used, the matrix X is premultiplied by 107*> diag(S) so that it solves the original system before 108*> equilibration. 109*> \endverbatim 110* 111* Arguments: 112* ========== 113* 114*> \verbatim 115*> Some optional parameters are bundled in the PARAMS array. These 116*> settings determine how refinement is performed, but often the 117*> defaults are acceptable. If the defaults are acceptable, users 118*> can pass NPARAMS = 0 which prevents the source code from accessing 119*> the PARAMS argument. 120*> \endverbatim 121*> 122*> \param[in] FACT 123*> \verbatim 124*> FACT is CHARACTER*1 125*> Specifies whether or not the factored form of the matrix A is 126*> supplied on entry, and if not, whether the matrix A should be 127*> equilibrated before it is factored. 128*> = 'F': On entry, AF contains the factored form of A. 129*> If EQUED is not 'N', the matrix A has been 130*> equilibrated with scaling factors given by S. 131*> A and AF are not modified. 132*> = 'N': The matrix A will be copied to AF and factored. 133*> = 'E': The matrix A will be equilibrated if necessary, then 134*> copied to AF and factored. 135*> \endverbatim 136*> 137*> \param[in] UPLO 138*> \verbatim 139*> UPLO is CHARACTER*1 140*> = 'U': Upper triangle of A is stored; 141*> = 'L': Lower triangle of A is stored. 142*> \endverbatim 143*> 144*> \param[in] N 145*> \verbatim 146*> N is INTEGER 147*> The number of linear equations, i.e., the order of the 148*> matrix A. N >= 0. 149*> \endverbatim 150*> 151*> \param[in] NRHS 152*> \verbatim 153*> NRHS is INTEGER 154*> The number of right hand sides, i.e., the number of columns 155*> of the matrices B and X. NRHS >= 0. 156*> \endverbatim 157*> 158*> \param[in,out] A 159*> \verbatim 160*> A is DOUBLE PRECISION array, dimension (LDA,N) 161*> On entry, the symmetric matrix A, except if FACT = 'F' and EQUED = 162*> 'Y', then A must contain the equilibrated matrix 163*> diag(S)*A*diag(S). If UPLO = 'U', the leading N-by-N upper 164*> triangular part of A contains the upper triangular part of the 165*> matrix A, and the strictly lower triangular part of A is not 166*> referenced. If UPLO = 'L', the leading N-by-N lower triangular 167*> part of A contains the lower triangular part of the matrix A, and 168*> the strictly upper triangular part of A is not referenced. A is 169*> not modified if FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 170*> 'N' on exit. 171*> 172*> On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by 173*> diag(S)*A*diag(S). 174*> \endverbatim 175*> 176*> \param[in] LDA 177*> \verbatim 178*> LDA is INTEGER 179*> The leading dimension of the array A. LDA >= max(1,N). 180*> \endverbatim 181*> 182*> \param[in,out] AF 183*> \verbatim 184*> AF is DOUBLE PRECISION array, dimension (LDAF,N) 185*> If FACT = 'F', then AF is an input argument and on entry 186*> contains the triangular factor U or L from the Cholesky 187*> factorization A = U**T*U or A = L*L**T, in the same storage 188*> format as A. If EQUED .ne. 'N', then AF is the factored 189*> form of the equilibrated matrix diag(S)*A*diag(S). 190*> 191*> If FACT = 'N', then AF is an output argument and on exit 192*> returns the triangular factor U or L from the Cholesky 193*> factorization A = U**T*U or A = L*L**T of the original 194*> matrix A. 195*> 196*> If FACT = 'E', then AF is an output argument and on exit 197*> returns the triangular factor U or L from the Cholesky 198*> factorization A = U**T*U or A = L*L**T of the equilibrated 199*> matrix A (see the description of A for the form of the 200*> equilibrated matrix). 201*> \endverbatim 202*> 203*> \param[in] LDAF 204*> \verbatim 205*> LDAF is INTEGER 206*> The leading dimension of the array AF. LDAF >= max(1,N). 207*> \endverbatim 208*> 209*> \param[in,out] EQUED 210*> \verbatim 211*> EQUED is CHARACTER*1 212*> Specifies the form of equilibration that was done. 213*> = 'N': No equilibration (always true if FACT = 'N'). 214*> = 'Y': Both row and column equilibration, i.e., A has been 215*> replaced by diag(S) * A * diag(S). 216*> EQUED is an input argument if FACT = 'F'; otherwise, it is an 217*> output argument. 218*> \endverbatim 219*> 220*> \param[in,out] S 221*> \verbatim 222*> S is DOUBLE PRECISION array, dimension (N) 223*> The row scale factors for A. If EQUED = 'Y', A is multiplied on 224*> the left and right by diag(S). S is an input argument if FACT = 225*> 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED 226*> = 'Y', each element of S must be positive. If S is output, each 227*> element of S is a power of the radix. If S is input, each element 228*> of S should be a power of the radix to ensure a reliable solution 229*> and error estimates. Scaling by powers of the radix does not cause 230*> rounding errors unless the result underflows or overflows. 231*> Rounding errors during scaling lead to refining with a matrix that 232*> is not equivalent to the input matrix, producing error estimates 233*> that may not be reliable. 234*> \endverbatim 235*> 236*> \param[in,out] B 237*> \verbatim 238*> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 239*> On entry, the N-by-NRHS right hand side matrix B. 240*> On exit, 241*> if EQUED = 'N', B is not modified; 242*> if EQUED = 'Y', B is overwritten by diag(S)*B; 243*> \endverbatim 244*> 245*> \param[in] LDB 246*> \verbatim 247*> LDB is INTEGER 248*> The leading dimension of the array B. LDB >= max(1,N). 249*> \endverbatim 250*> 251*> \param[out] X 252*> \verbatim 253*> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 254*> If INFO = 0, the N-by-NRHS solution matrix X to the original 255*> system of equations. Note that A and B are modified on exit if 256*> EQUED .ne. 'N', and the solution to the equilibrated system is 257*> inv(diag(S))*X. 258*> \endverbatim 259*> 260*> \param[in] LDX 261*> \verbatim 262*> LDX is INTEGER 263*> The leading dimension of the array X. LDX >= max(1,N). 264*> \endverbatim 265*> 266*> \param[out] RCOND 267*> \verbatim 268*> RCOND is DOUBLE PRECISION 269*> Reciprocal scaled condition number. This is an estimate of the 270*> reciprocal Skeel condition number of the matrix A after 271*> equilibration (if done). If this is less than the machine 272*> precision (in particular, if it is zero), the matrix is singular 273*> to working precision. Note that the error may still be small even 274*> if this number is very small and the matrix appears ill- 275*> conditioned. 276*> \endverbatim 277*> 278*> \param[out] RPVGRW 279*> \verbatim 280*> RPVGRW is DOUBLE PRECISION 281*> Reciprocal pivot growth. On exit, this contains the reciprocal 282*> pivot growth factor norm(A)/norm(U). The "max absolute element" 283*> norm is used. If this is much less than 1, then the stability of 284*> the LU factorization of the (equilibrated) matrix A could be poor. 285*> This also means that the solution X, estimated condition numbers, 286*> and error bounds could be unreliable. If factorization fails with 287*> 0<INFO<=N, then this contains the reciprocal pivot growth factor 288*> for the leading INFO columns of A. 289*> \endverbatim 290*> 291*> \param[out] BERR 292*> \verbatim 293*> BERR is DOUBLE PRECISION array, dimension (NRHS) 294*> Componentwise relative backward error. This is the 295*> componentwise relative backward error of each solution vector X(j) 296*> (i.e., the smallest relative change in any element of A or B that 297*> makes X(j) an exact solution). 298*> \endverbatim 299*> 300*> \param[in] N_ERR_BNDS 301*> \verbatim 302*> N_ERR_BNDS is INTEGER 303*> Number of error bounds to return for each right hand side 304*> and each type (normwise or componentwise). See ERR_BNDS_NORM and 305*> ERR_BNDS_COMP below. 306*> \endverbatim 307*> 308*> \param[out] ERR_BNDS_NORM 309*> \verbatim 310*> ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS) 311*> For each right-hand side, this array contains information about 312*> various error bounds and condition numbers corresponding to the 313*> normwise relative error, which is defined as follows: 314*> 315*> Normwise relative error in the ith solution vector: 316*> max_j (abs(XTRUE(j,i) - X(j,i))) 317*> ------------------------------ 318*> max_j abs(X(j,i)) 319*> 320*> The array is indexed by the type of error information as described 321*> below. There currently are up to three pieces of information 322*> returned. 323*> 324*> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith 325*> right-hand side. 326*> 327*> The second index in ERR_BNDS_NORM(:,err) contains the following 328*> three fields: 329*> err = 1 "Trust/don't trust" boolean. Trust the answer if the 330*> reciprocal condition number is less than the threshold 331*> sqrt(n) * dlamch('Epsilon'). 332*> 333*> err = 2 "Guaranteed" error bound: The estimated forward error, 334*> almost certainly within a factor of 10 of the true error 335*> so long as the next entry is greater than the threshold 336*> sqrt(n) * dlamch('Epsilon'). This error bound should only 337*> be trusted if the previous boolean is true. 338*> 339*> err = 3 Reciprocal condition number: Estimated normwise 340*> reciprocal condition number. Compared with the threshold 341*> sqrt(n) * dlamch('Epsilon') to determine if the error 342*> estimate is "guaranteed". These reciprocal condition 343*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 344*> appropriately scaled matrix Z. 345*> Let Z = S*A, where S scales each row by a power of the 346*> radix so all absolute row sums of Z are approximately 1. 347*> 348*> See Lapack Working Note 165 for further details and extra 349*> cautions. 350*> \endverbatim 351*> 352*> \param[out] ERR_BNDS_COMP 353*> \verbatim 354*> ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS) 355*> For each right-hand side, this array contains information about 356*> various error bounds and condition numbers corresponding to the 357*> componentwise relative error, which is defined as follows: 358*> 359*> Componentwise relative error in the ith solution vector: 360*> abs(XTRUE(j,i) - X(j,i)) 361*> max_j ---------------------- 362*> abs(X(j,i)) 363*> 364*> The array is indexed by the right-hand side i (on which the 365*> componentwise relative error depends), and the type of error 366*> information as described below. There currently are up to three 367*> pieces of information returned for each right-hand side. If 368*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then 369*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most 370*> the first (:,N_ERR_BNDS) entries are returned. 371*> 372*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith 373*> right-hand side. 374*> 375*> The second index in ERR_BNDS_COMP(:,err) contains the following 376*> three fields: 377*> err = 1 "Trust/don't trust" boolean. Trust the answer if the 378*> reciprocal condition number is less than the threshold 379*> sqrt(n) * dlamch('Epsilon'). 380*> 381*> err = 2 "Guaranteed" error bound: The estimated forward error, 382*> almost certainly within a factor of 10 of the true error 383*> so long as the next entry is greater than the threshold 384*> sqrt(n) * dlamch('Epsilon'). This error bound should only 385*> be trusted if the previous boolean is true. 386*> 387*> err = 3 Reciprocal condition number: Estimated componentwise 388*> reciprocal condition number. Compared with the threshold 389*> sqrt(n) * dlamch('Epsilon') to determine if the error 390*> estimate is "guaranteed". These reciprocal condition 391*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some 392*> appropriately scaled matrix Z. 393*> Let Z = S*(A*diag(x)), where x is the solution for the 394*> current right-hand side and S scales each row of 395*> A*diag(x) by a power of the radix so all absolute row 396*> sums of Z are approximately 1. 397*> 398*> See Lapack Working Note 165 for further details and extra 399*> cautions. 400*> \endverbatim 401*> 402*> \param[in] NPARAMS 403*> \verbatim 404*> NPARAMS is INTEGER 405*> Specifies the number of parameters set in PARAMS. If .LE. 0, the 406*> PARAMS array is never referenced and default values are used. 407*> \endverbatim 408*> 409*> \param[in,out] PARAMS 410*> \verbatim 411*> PARAMS is DOUBLE PRECISION array, dimension NPARAMS 412*> Specifies algorithm parameters. If an entry is .LT. 0.0, then 413*> that entry will be filled with default value used for that 414*> parameter. Only positions up to NPARAMS are accessed; defaults 415*> are used for higher-numbered parameters. 416*> 417*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative 418*> refinement or not. 419*> Default: 1.0D+0 420*> = 0.0 : No refinement is performed, and no error bounds are 421*> computed. 422*> = 1.0 : Use the extra-precise refinement algorithm. 423*> (other values are reserved for future use) 424*> 425*> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual 426*> computations allowed for refinement. 427*> Default: 10 428*> Aggressive: Set to 100 to permit convergence using approximate 429*> factorizations or factorizations other than LU. If 430*> the factorization uses a technique other than 431*> Gaussian elimination, the guarantees in 432*> err_bnds_norm and err_bnds_comp may no longer be 433*> trustworthy. 434*> 435*> PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code 436*> will attempt to find a solution with small componentwise 437*> relative error in the double-precision algorithm. Positive 438*> is true, 0.0 is false. 439*> Default: 1.0 (attempt componentwise convergence) 440*> \endverbatim 441*> 442*> \param[out] WORK 443*> \verbatim 444*> WORK is DOUBLE PRECISION array, dimension (4*N) 445*> \endverbatim 446*> 447*> \param[out] IWORK 448*> \verbatim 449*> IWORK is INTEGER array, dimension (N) 450*> \endverbatim 451*> 452*> \param[out] INFO 453*> \verbatim 454*> INFO is INTEGER 455*> = 0: Successful exit. The solution to every right-hand side is 456*> guaranteed. 457*> < 0: If INFO = -i, the i-th argument had an illegal value 458*> > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization 459*> has been completed, but the factor U is exactly singular, so 460*> the solution and error bounds could not be computed. RCOND = 0 461*> is returned. 462*> = N+J: The solution corresponding to the Jth right-hand side is 463*> not guaranteed. The solutions corresponding to other right- 464*> hand sides K with K > J may not be guaranteed as well, but 465*> only the first such right-hand side is reported. If a small 466*> componentwise error is not requested (PARAMS(3) = 0.0) then 467*> the Jth right-hand side is the first with a normwise error 468*> bound that is not guaranteed (the smallest J such 469*> that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0) 470*> the Jth right-hand side is the first with either a normwise or 471*> componentwise error bound that is not guaranteed (the smallest 472*> J such that either ERR_BNDS_NORM(J,1) = 0.0 or 473*> ERR_BNDS_COMP(J,1) = 0.0). See the definition of 474*> ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information 475*> about all of the right-hand sides check ERR_BNDS_NORM or 476*> ERR_BNDS_COMP. 477*> \endverbatim 478* 479* Authors: 480* ======== 481* 482*> \author Univ. of Tennessee 483*> \author Univ. of California Berkeley 484*> \author Univ. of Colorado Denver 485*> \author NAG Ltd. 486* 487*> \date April 2012 488* 489*> \ingroup doublePOsolve 490* 491* ===================================================================== 492 SUBROUTINE DPOSVXX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 493 $ S, B, LDB, X, LDX, RCOND, RPVGRW, BERR, 494 $ N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, 495 $ NPARAMS, PARAMS, WORK, IWORK, INFO ) 496* 497* -- LAPACK driver routine (version 3.4.1) -- 498* -- LAPACK is a software package provided by Univ. of Tennessee, -- 499* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 500* April 2012 501* 502* .. Scalar Arguments .. 503 CHARACTER EQUED, FACT, UPLO 504 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS, 505 $ N_ERR_BNDS 506 DOUBLE PRECISION RCOND, RPVGRW 507* .. 508* .. Array Arguments .. 509 INTEGER IWORK( * ) 510 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 511 $ X( LDX, * ), WORK( * ) 512 DOUBLE PRECISION S( * ), PARAMS( * ), BERR( * ), 513 $ ERR_BNDS_NORM( NRHS, * ), 514 $ ERR_BNDS_COMP( NRHS, * ) 515* .. 516* 517* ================================================================== 518* 519* .. Parameters .. 520 DOUBLE PRECISION ZERO, ONE 521 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 522 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 523 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 524 INTEGER CMP_ERR_I, PIV_GROWTH_I 525 PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 526 $ BERR_I = 3 ) 527 PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 528 PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8, 529 $ PIV_GROWTH_I = 9 ) 530* .. 531* .. Local Scalars .. 532 LOGICAL EQUIL, NOFACT, RCEQU 533 INTEGER INFEQU, J 534 DOUBLE PRECISION AMAX, BIGNUM, SMIN, SMAX, 535 $ SCOND, SMLNUM 536* .. 537* .. External Functions .. 538 EXTERNAL LSAME, DLAMCH, DLA_PORPVGRW 539 LOGICAL LSAME 540 DOUBLE PRECISION DLAMCH, DLA_PORPVGRW 541* .. 542* .. External Subroutines .. 543 EXTERNAL DPOEQUB, DPOTRF, DPOTRS, DLACPY, DLAQSY, 544 $ XERBLA, DLASCL2, DPORFSX 545* .. 546* .. Intrinsic Functions .. 547 INTRINSIC MAX, MIN 548* .. 549* .. Executable Statements .. 550* 551 INFO = 0 552 NOFACT = LSAME( FACT, 'N' ) 553 EQUIL = LSAME( FACT, 'E' ) 554 SMLNUM = DLAMCH( 'Safe minimum' ) 555 BIGNUM = ONE / SMLNUM 556 IF( NOFACT .OR. EQUIL ) THEN 557 EQUED = 'N' 558 RCEQU = .FALSE. 559 ELSE 560 RCEQU = LSAME( EQUED, 'Y' ) 561 ENDIF 562* 563* Default is failure. If an input parameter is wrong or 564* factorization fails, make everything look horrible. Only the 565* pivot growth is set here, the rest is initialized in DPORFSX. 566* 567 RPVGRW = ZERO 568* 569* Test the input parameters. PARAMS is not tested until DPORFSX. 570* 571 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT. 572 $ LSAME( FACT, 'F' ) ) THEN 573 INFO = -1 574 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. 575 $ .NOT.LSAME( UPLO, 'L' ) ) THEN 576 INFO = -2 577 ELSE IF( N.LT.0 ) THEN 578 INFO = -3 579 ELSE IF( NRHS.LT.0 ) THEN 580 INFO = -4 581 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 582 INFO = -6 583 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 584 INFO = -8 585 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT. 586 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN 587 INFO = -9 588 ELSE 589 IF ( RCEQU ) THEN 590 SMIN = BIGNUM 591 SMAX = ZERO 592 DO 10 J = 1, N 593 SMIN = MIN( SMIN, S( J ) ) 594 SMAX = MAX( SMAX, S( J ) ) 595 10 CONTINUE 596 IF( SMIN.LE.ZERO ) THEN 597 INFO = -10 598 ELSE IF( N.GT.0 ) THEN 599 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) 600 ELSE 601 SCOND = ONE 602 END IF 603 END IF 604 IF( INFO.EQ.0 ) THEN 605 IF( LDB.LT.MAX( 1, N ) ) THEN 606 INFO = -12 607 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 608 INFO = -14 609 END IF 610 END IF 611 END IF 612* 613 IF( INFO.NE.0 ) THEN 614 CALL XERBLA( 'DPOSVXX', -INFO ) 615 RETURN 616 END IF 617* 618 IF( EQUIL ) THEN 619* 620* Compute row and column scalings to equilibrate the matrix A. 621* 622 CALL DPOEQUB( N, A, LDA, S, SCOND, AMAX, INFEQU ) 623 IF( INFEQU.EQ.0 ) THEN 624* 625* Equilibrate the matrix. 626* 627 CALL DLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 628 RCEQU = LSAME( EQUED, 'Y' ) 629 END IF 630 END IF 631* 632* Scale the right-hand side. 633* 634 IF( RCEQU ) CALL DLASCL2( N, NRHS, S, B, LDB ) 635* 636 IF( NOFACT .OR. EQUIL ) THEN 637* 638* Compute the Cholesky factorization of A. 639* 640 CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 641 CALL DPOTRF( UPLO, N, AF, LDAF, INFO ) 642* 643* Return if INFO is non-zero. 644* 645 IF( INFO.NE.0 ) THEN 646* 647* Pivot in column INFO is exactly 0 648* Compute the reciprocal pivot growth factor of the 649* leading rank-deficient INFO columns of A. 650* 651 RPVGRW = DLA_PORPVGRW( UPLO, INFO, A, LDA, AF, LDAF, WORK ) 652 RETURN 653 ENDIF 654 END IF 655* 656* Compute the reciprocal growth factor RPVGRW. 657* 658 RPVGRW = DLA_PORPVGRW( UPLO, N, A, LDA, AF, LDAF, WORK ) 659* 660* Compute the solution matrix X. 661* 662 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 663 CALL DPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO ) 664* 665* Use iterative refinement to improve the computed solution and 666* compute error bounds and backward error estimates for it. 667* 668 CALL DPORFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, 669 $ S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, 670 $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO ) 671 672* 673* Scale solutions. 674* 675 IF ( RCEQU ) THEN 676 CALL DLASCL2 ( N, NRHS, S, X, LDX ) 677 END IF 678* 679 RETURN 680* 681* End of DPOSVXX 682* 683 END 684