001:SUBROUTINESPOSVXX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 002: $ S, B, LDB, X, LDX, RCOND, RPVGRW, BERR, 003: $ N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, 004: $ NPARAMS, PARAMS, WORK, IWORK, INFO ) 005:*006:* -- LAPACK driver routine (version 3.2) --007:* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --008:* -- Jason Riedy of Univ. of California Berkeley. --009:* -- November 2008 --010:*011:* -- LAPACK is a software package provided by Univ. of Tennessee, --012:* -- Univ. of California Berkeley and NAG Ltd. --013:*014:IMPLICITNONE 015:* ..016:* .. Scalar Arguments ..017: CHARACTER EQUED, FACT, UPLO 018: INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS, 019: $ N_ERR_BNDS 020: REAL RCOND, RPVGRW 021:* ..022:* .. Array Arguments ..023: INTEGERIWORK( * ) 024: REALA( LDA, * ),AF( LDAF, * ),B( LDB, * ), 025: $X( LDX, * ),WORK( * ) 026: REALS( * ),PARAMS( * ),BERR( * ), 027: $ERR_BNDS_NORM( NRHS, * ), 028: $ERR_BNDS_COMP( NRHS, * ) 029:* ..030:*031:* Purpose032:* =======033:*034:* SPOSVXX uses the Cholesky factorization A = U**T*U or A = L*L**T035:* to compute the solution to a real system of linear equations036:* A * X = B, where A is an N-by-N symmetric positive definite matrix037:* and X and B are N-by-NRHS matrices.038:*039:* If requested, both normwise and maximum componentwise error bounds040:* are returned. SPOSVXX will return a solution with a tiny041:* guaranteed error (O(eps) where eps is the working machine042:* precision) unless the matrix is very ill-conditioned, in which043:* case a warning is returned. Relevant condition numbers also are044:* calculated and returned.045:*046:* SPOSVXX accepts user-provided factorizations and equilibration047:* factors; see the definitions of the FACT and EQUED options.048:* Solving with refinement and using a factorization from a previous049:* SPOSVXX call will also produce a solution with either O(eps)050:* errors or warnings, but we cannot make that claim for general051:* user-provided factorizations and equilibration factors if they052:* differ from what SPOSVXX would itself produce.053:*054:* Description055:* ===========056:*057:* The following steps are performed:058:*059:* 1. If FACT = 'E', real scaling factors are computed to equilibrate060:* the system:061:*062:* diag(S)*A*diag(S) *inv(diag(S))*X = diag(S)*B063:*064:* Whether or not the system will be equilibrated depends on the065:* scaling of the matrix A, but if equilibration is used, A is066:* overwritten by diag(S)*A*diag(S) and B by diag(S)*B.067:*068:* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to069:* factor the matrix A (after equilibration if FACT = 'E') as070:* A = U**T* U, if UPLO = 'U', or071:* A = L * L**T, if UPLO = 'L',072:* where U is an upper triangular matrix and L is a lower triangular073:* matrix.074:*075:* 3. If the leading i-by-i principal minor is not positive definite,076:* then the routine returns with INFO = i. Otherwise, the factored077:* form of A is used to estimate the condition number of the matrix078:* A (see argument RCOND). If the reciprocal of the condition number079:* is less than machine precision, the routine still goes on to solve080:* for X and compute error bounds as described below.081:*082:* 4. The system of equations is solved for X using the factored form083:* of A.084:*085:* 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),086:* the routine will use iterative refinement to try to get a small087:* error and error bounds. Refinement calculates the residual to at088:* least twice the working precision.089:*090:* 6. If equilibration was used, the matrix X is premultiplied by091:* diag(S) so that it solves the original system before092:* equilibration.093:*094:* Arguments095:* =========096:*097:* Some optional parameters are bundled in the PARAMS array. These098:* settings determine how refinement is performed, but often the099:* defaults are acceptable. If the defaults are acceptable, users100:* can pass NPARAMS = 0 which prevents the source code from accessing101:* the PARAMS argument.102:*103:* FACT (input) CHARACTER*1104:* Specifies whether or not the factored form of the matrix A is105:* supplied on entry, and if not, whether the matrix A should be106:* equilibrated before it is factored.107:* = 'F': On entry, AF contains the factored form of A.108:* If EQUED is not 'N', the matrix A has been109:* equilibrated with scaling factors given by S.110:* A and AF are not modified.111:* = 'N': The matrix A will be copied to AF and factored.112:* = 'E': The matrix A will be equilibrated if necessary, then113:* copied to AF and factored.114:*115:* UPLO (input) CHARACTER*1116:* = 'U': Upper triangle of A is stored;117:* = 'L': Lower triangle of A is stored.118:*119:* N (input) INTEGER120:* The number of linear equations, i.e., the order of the121:* matrix A. N >= 0.122:*123:* NRHS (input) INTEGER124:* The number of right hand sides, i.e., the number of columns125:* of the matrices B and X. NRHS >= 0.126:*127:* A (input/output) REAL array, dimension (LDA,N)128:* On entry, the symmetric matrix A, except if FACT = 'F' and EQUED =129:* 'Y', then A must contain the equilibrated matrix130:* diag(S)*A*diag(S). If UPLO = 'U', the leading N-by-N upper131:* triangular part of A contains the upper triangular part of the132:* matrix A, and the strictly lower triangular part of A is not133:* referenced. If UPLO = 'L', the leading N-by-N lower triangular134:* part of A contains the lower triangular part of the matrix A, and135:* the strictly upper triangular part of A is not referenced. A is136:* not modified if FACT = 'F' or 'N', or if FACT = 'E' and EQUED =137:* 'N' on exit.138:*139:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by140:* diag(S)*A*diag(S).141:*142:* LDA (input) INTEGER143:* The leading dimension of the array A. LDA >= max(1,N).144:*145:* AF (input or output) REAL array, dimension (LDAF,N)146:* If FACT = 'F', then AF is an input argument and on entry147:* contains the triangular factor U or L from the Cholesky148:* factorization A = U**T*U or A = L*L**T, in the same storage149:* format as A. If EQUED .ne. 'N', then AF is the factored150:* form of the equilibrated matrix diag(S)*A*diag(S).151:*152:* If FACT = 'N', then AF is an output argument and on exit153:* returns the triangular factor U or L from the Cholesky154:* factorization A = U**T*U or A = L*L**T of the original155:* matrix A.156:*157:* If FACT = 'E', then AF is an output argument and on exit158:* returns the triangular factor U or L from the Cholesky159:* factorization A = U**T*U or A = L*L**T of the equilibrated160:* matrix A (see the description of A for the form of the161:* equilibrated matrix).162:*163:* LDAF (input) INTEGER164:* The leading dimension of the array AF. LDAF >= max(1,N).165:*166:* EQUED (input or output) CHARACTER*1167:* Specifies the form of equilibration that was done.168:* = 'N': No equilibration (always true if FACT = 'N').169:* = 'Y': Both row and column equilibration, i.e., A has been170:* replaced by diag(S) * A * diag(S).171:* EQUED is an input argument if FACT = 'F'; otherwise, it is an172:* output argument.173:*174:* S (input or output) REAL array, dimension (N)175:* The row scale factors for A. If EQUED = 'Y', A is multiplied on176:* the left and right by diag(S). S is an input argument if FACT =177:* 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED178:* = 'Y', each element of S must be positive. If S is output, each179:* element of S is a power of the radix. If S is input, each element180:* of S should be a power of the radix to ensure a reliable solution181:* and error estimates. Scaling by powers of the radix does not cause182:* rounding errors unless the result underflows or overflows.183:* Rounding errors during scaling lead to refining with a matrix that184:* is not equivalent to the input matrix, producing error estimates185:* that may not be reliable.186:*187:* B (input/output) REAL array, dimension (LDB,NRHS)188:* On entry, the N-by-NRHS right hand side matrix B.189:* On exit,190:* if EQUED = 'N', B is not modified;191:* if EQUED = 'Y', B is overwritten by diag(S)*B;192:*193:* LDB (input) INTEGER194:* The leading dimension of the array B. LDB >= max(1,N).195:*196:* X (output) REAL array, dimension (LDX,NRHS)197:* If INFO = 0, the N-by-NRHS solution matrix X to the original198:* system of equations. Note that A and B are modified on exit if199:* EQUED .ne. 'N', and the solution to the equilibrated system is200:* inv(diag(S))*X.201:*202:* LDX (input) INTEGER203:* The leading dimension of the array X. LDX >= max(1,N).204:*205:* RCOND (output) REAL206:* Reciprocal scaled condition number. This is an estimate of the207:* reciprocal Skeel condition number of the matrix A after208:* equilibration (if done). If this is less than the machine209:* precision (in particular, if it is zero), the matrix is singular210:* to working precision. Note that the error may still be small even211:* if this number is very small and the matrix appears ill-212:* conditioned.213:*214:* RPVGRW (output) REAL215:* Reciprocal pivot growth. On exit, this contains the reciprocal216:* pivot growth factor norm(A)/norm(U). The "max absolute element"217:* norm is used. If this is much less than 1, then the stability of218:* the LU factorization of the (equilibrated) matrix A could be poor.219:* This also means that the solution X, estimated condition numbers,220:* and error bounds could be unreliable. If factorization fails with221:* 0<INFO<=N, then this contains the reciprocal pivot growth factor222:* for the leading INFO columns of A.223:*224:* BERR (output) REAL array, dimension (NRHS)225:* Componentwise relative backward error. This is the226:* componentwise relative backward error of each solution vector X(j)227:* (i.e., the smallest relative change in any element of A or B that228:* makes X(j) an exact solution).229:*230:* N_ERR_BNDS (input) INTEGER231:* Number of error bounds to return for each right hand side232:* and each type (normwise or componentwise). See ERR_BNDS_NORM and233:* ERR_BNDS_COMP below.234:*235:* ERR_BNDS_NORM (output) REAL array, dimension (NRHS, N_ERR_BNDS)236:* For each right-hand side, this array contains information about237:* various error bounds and condition numbers corresponding to the238:* normwise relative error, which is defined as follows:239:*240:* Normwise relative error in the ith solution vector:241:* max_j (abs(XTRUE(j,i) - X(j,i)))242:* ------------------------------243:* max_j abs(X(j,i))244:*245:* The array is indexed by the type of error information as described246:* below. There currently are up to three pieces of information247:* returned.248:*249:* The first index in ERR_BNDS_NORM(i,:) corresponds to the ith250:* right-hand side.251:*252:* The second index in ERR_BNDS_NORM(:,err) contains the following253:* three fields:254:* err = 1 "Trust/don't trust" boolean. Trust the answer if the255:* reciprocal condition number is less than the threshold256:* sqrt(n) * slamch('Epsilon').257:*258:* err = 2 "Guaranteed" error bound: The estimated forward error,259:* almost certainly within a factor of 10 of the true error260:* so long as the next entry is greater than the threshold261:* sqrt(n) * slamch('Epsilon'). This error bound should only262:* be trusted if the previous boolean is true.263:*264:* err = 3 Reciprocal condition number: Estimated normwise265:* reciprocal condition number. Compared with the threshold266:* sqrt(n) * slamch('Epsilon') to determine if the error267:* estimate is "guaranteed". These reciprocal condition268:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some269:* appropriately scaled matrix Z.270:* Let Z = S*A, where S scales each row by a power of the271:* radix so all absolute row sums of Z are approximately 1.272:*273:* See Lapack Working Note 165 for further details and extra274:* cautions.275:*276:* ERR_BNDS_COMP (output) REAL array, dimension (NRHS, N_ERR_BNDS)277:* For each right-hand side, this array contains information about278:* various error bounds and condition numbers corresponding to the279:* componentwise relative error, which is defined as follows:280:*281:* Componentwise relative error in the ith solution vector:282:* abs(XTRUE(j,i) - X(j,i))283:* max_j ----------------------284:* abs(X(j,i))285:*286:* The array is indexed by the right-hand side i (on which the287:* componentwise relative error depends), and the type of error288:* information as described below. There currently are up to three289:* pieces of information returned for each right-hand side. If290:* componentwise accuracy is not requested (PARAMS(3) = 0.0), then291:* ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most292:* the first (:,N_ERR_BNDS) entries are returned.293:*294:* The first index in ERR_BNDS_COMP(i,:) corresponds to the ith295:* right-hand side.296:*297:* The second index in ERR_BNDS_COMP(:,err) contains the following298:* three fields:299:* err = 1 "Trust/don't trust" boolean. Trust the answer if the300:* reciprocal condition number is less than the threshold301:* sqrt(n) * slamch('Epsilon').302:*303:* err = 2 "Guaranteed" error bound: The estimated forward error,304:* almost certainly within a factor of 10 of the true error305:* so long as the next entry is greater than the threshold306:* sqrt(n) * slamch('Epsilon'). This error bound should only307:* be trusted if the previous boolean is true.308:*309:* err = 3 Reciprocal condition number: Estimated componentwise310:* reciprocal condition number. Compared with the threshold311:* sqrt(n) * slamch('Epsilon') to determine if the error312:* estimate is "guaranteed". These reciprocal condition313:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some314:* appropriately scaled matrix Z.315:* Let Z = S*(A*diag(x)), where x is the solution for the316:* current right-hand side and S scales each row of317:* A*diag(x) by a power of the radix so all absolute row318:* sums of Z are approximately 1.319:*320:* See Lapack Working Note 165 for further details and extra321:* cautions.322:*323:* NPARAMS (input) INTEGER324:* Specifies the number of parameters set in PARAMS. If .LE. 0, the325:* PARAMS array is never referenced and default values are used.326:*327:* PARAMS (input / output) REAL array, dimension NPARAMS328:* Specifies algorithm parameters. If an entry is .LT. 0.0, then329:* that entry will be filled with default value used for that330:* parameter. Only positions up to NPARAMS are accessed; defaults331:* are used for higher-numbered parameters.332:*333:* PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative334:* refinement or not.335:* Default: 1.0336:* = 0.0 : No refinement is performed, and no error bounds are337:* computed.338:* = 1.0 : Use the double-precision refinement algorithm,339:* possibly with doubled-single computations if the340:* compilation environment does not support DOUBLE341:* PRECISION.342:* (other values are reserved for future use)343:*344:* PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual345:* computations allowed for refinement.346:* Default: 10347:* Aggressive: Set to 100 to permit convergence using approximate348:* factorizations or factorizations other than LU. If349:* the factorization uses a technique other than350:* Gaussian elimination, the guarantees in351:* err_bnds_norm and err_bnds_comp may no longer be352:* trustworthy.353:*354:* PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code355:* will attempt to find a solution with small componentwise356:* relative error in the double-precision algorithm. Positive357:* is true, 0.0 is false.358:* Default: 1.0 (attempt componentwise convergence)359:*360:* WORK (workspace) REAL array, dimension (4*N)361:*362:* IWORK (workspace) INTEGER array, dimension (N)363:*364:* INFO (output) INTEGER365:* = 0: Successful exit. The solution to every right-hand side is366:* guaranteed.367:* < 0: If INFO = -i, the i-th argument had an illegal value368:* > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization369:* has been completed, but the factor U is exactly singular, so370:* the solution and error bounds could not be computed. RCOND = 0371:* is returned.372:* = N+J: The solution corresponding to the Jth right-hand side is373:* not guaranteed. The solutions corresponding to other right-374:* hand sides K with K > J may not be guaranteed as well, but375:* only the first such right-hand side is reported. If a small376:* componentwise error is not requested (PARAMS(3) = 0.0) then377:* the Jth right-hand side is the first with a normwise error378:* bound that is not guaranteed (the smallest J such379:* that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)380:* the Jth right-hand side is the first with either a normwise or381:* componentwise error bound that is not guaranteed (the smallest382:* J such that either ERR_BNDS_NORM(J,1) = 0.0 or383:* ERR_BNDS_COMP(J,1) = 0.0). See the definition of384:* ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information385:* about all of the right-hand sides check ERR_BNDS_NORM or386:* ERR_BNDS_COMP.387:*388:* ==================================================================389:*390:* .. Parameters ..391: REAL ZERO, ONE 392:PARAMETER( ZERO = 0.0E+0, ONE = 1.0E+0 ) 393: INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 394: INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 395: INTEGER CMP_ERR_I, PIV_GROWTH_I 396:PARAMETER( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 397: $ BERR_I = 3 ) 398:PARAMETER( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 399:PARAMETER( CMP_RCOND_I = 7, CMP_ERR_I = 8, 400: $ PIV_GROWTH_I = 9 ) 401:* ..402:* .. Local Scalars ..403:LOGICALEQUIL, NOFACT, RCEQU 404: INTEGER INFEQU, J 405: REAL AMAX, BIGNUM, SMIN, SMAX, 406: $ SCOND, SMLNUM 407:* ..408:* .. External Functions ..409:EXTERNALLSAME, SLAMCH, SLA_PORPVGRW 410:LOGICALLSAME 411: REAL SLAMCH, SLA_PORPVGRW 412:* ..413:* .. External Subroutines ..414:EXTERNALSPOEQUB, SPOTRF, SPOTRS, SLACPY, SLAQSY, 415: $ XERBLA, SLASCL2, SPORFSX 416:* ..417:* .. Intrinsic Functions ..418:INTRINSICMAX, MIN 419:* ..420:* .. Executable Statements ..421:*422: INFO = 0 423: NOFACT =LSAME( FACT, 'N' ) 424: EQUIL =LSAME( FACT, 'E' ) 425: SMLNUM =SLAMCH( 'Safe minimum' ) 426: BIGNUM = ONE / SMLNUM 427:IF( NOFACT .OR. EQUIL )THEN428: EQUED = 'N' 429: RCEQU = .FALSE. 430:ELSE431: RCEQU =LSAME( EQUED, 'Y' ) 432:ENDIF433:*434:* Default is failure. If an input parameter is wrong or435:* factorization fails, make everything look horrible. Only the436:* pivot growth is set here, the rest is initialized in SPORFSX.437:*438: RPVGRW = ZERO 439:*440:* Test the input parameters. PARAMS is not tested until SPORFSX.441:*442:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT. 443: $LSAME( FACT, 'F' ) )THEN444: INFO = -1 445:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. 446: $ .NOT.LSAME( UPLO, 'L' ) )THEN447: INFO = -2 448:ELSEIF( N.LT.0 )THEN449: INFO = -3 450:ELSEIF( NRHS.LT.0 )THEN451: INFO = -4 452:ELSEIF( LDA.LT.MAX( 1, N ) )THEN453: INFO = -6 454:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN455: INFO = -8 456:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 457: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN458: INFO = -9 459:ELSE460:IF( RCEQU )THEN461: SMIN = BIGNUM 462: SMAX = ZERO 463:DO10 J = 1, N 464: SMIN =MIN( SMIN,S( J ) ) 465: SMAX =MAX( SMAX,S( J ) ) 466: 10CONTINUE467:IF( SMIN.LE.ZERO )THEN468: INFO = -10 469:ELSEIF( N.GT.0 )THEN470: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 471:ELSE472: SCOND = ONE 473:ENDIF474:ENDIF475:IF( INFO.EQ.0 )THEN476:IF( LDB.LT.MAX( 1, N ) )THEN477: INFO = -12 478:ELSEIF( LDX.LT.MAX( 1, N ) )THEN479: INFO = -14 480:ENDIF481:ENDIF482:ENDIF483:*484:IF( INFO.NE.0 )THEN485:CALLXERBLA( 'SPOSVXX', -INFO ) 486:RETURN487:ENDIF488:*489:IF( EQUIL )THEN490:*491:* Compute row and column scalings to equilibrate the matrix A.492:*493:CALLSPOEQUB( N, A, LDA, S, SCOND, AMAX, INFEQU ) 494:IF( INFEQU.EQ.0 )THEN495:*496:* Equilibrate the matrix.497:*498:CALLSLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 499: RCEQU =LSAME( EQUED, 'Y' ) 500:ENDIF501:ENDIF502:*503:* Scale the right-hand side.504:*505:IF( RCEQU )CALLSLASCL2( N, NRHS, S, B, LDB ) 506:*507:IF( NOFACT .OR. EQUIL )THEN508:*509:* Compute the LU factorization of A.510:*511:CALLSLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 512:CALLSPOTRF( UPLO, N, AF, LDAF, INFO ) 513:*514:* Return if INFO is non-zero.515:*516:IF( INFO.NE.0 )THEN517:*518:* Pivot in column INFO is exactly 0519:* Compute the reciprocal pivot growth factor of the520:* leading rank-deficient INFO columns of A.521:*522: RPVGRW =SLA_PORPVGRW( UPLO, INFO, A, LDA, AF, LDAF, WORK ) 523:RETURN524:ENDIF525:ENDIF526:*527:* Compute the reciprocal growth factor RPVGRW.528:*529: RPVGRW =SLA_PORPVGRW( UPLO, N, A, LDA, AF, LDAF, WORK ) 530:*531:* Compute the solution matrix X.532:*533:CALLSLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 534:CALLSPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO ) 535:*536:* Use iterative refinement to improve the computed solution and537:* compute error bounds and backward error estimates for it.538:*539:CALLSPORFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, 540: $ S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, 541: $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO ) 542: 543:*544:* Scale solutions.545:*546:IF( RCEQU )THEN547:CALLSLASCL2( N, NRHS, S, X, LDX ) 548:ENDIF549:*550:RETURN551:*552:* End of SPOSVXX553:*554:END555: