001:SUBROUTINESSYSVXX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, 002: $ EQUED, 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 routine (version 3.2.1) --007:* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --008:* -- Jason Riedy of Univ. of California Berkeley. --009:* -- April 2009 --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: INTEGERIPIV( * ),IWORK( * ) 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:* SSYSVXX uses the diagonal pivoting factorization to compute the035:* solution to a real system of linear equations A * X = B, where A036:* is an N-by-N symmetric matrix and X and B are N-by-NRHS matrices.037:*038:* If requested, both normwise and maximum componentwise error bounds039:* are returned. SSYSVXX will return a solution with a tiny040:* guaranteed error (O(eps) where eps is the working machine041:* precision) unless the matrix is very ill-conditioned, in which042:* case a warning is returned. Relevant condition numbers also are043:* calculated and returned.044:*045:* SSYSVXX accepts user-provided factorizations and equilibration046:* factors; see the definitions of the FACT and EQUED options.047:* Solving with refinement and using a factorization from a previous048:* SSYSVXX call will also produce a solution with either O(eps)049:* errors or warnings, but we cannot make that claim for general050:* user-provided factorizations and equilibration factors if they051:* differ from what SSYSVXX would itself produce.052:*053:* Description054:* ===========055:*056:* The following steps are performed:057:*058:* 1. If FACT = 'E', real scaling factors are computed to equilibrate059:* the system:060:*061:* diag(S)*A*diag(S) *inv(diag(S))*X = diag(S)*B062:*063:* Whether or not the system will be equilibrated depends on the064:* scaling of the matrix A, but if equilibration is used, A is065:* overwritten by diag(S)*A*diag(S) and B by diag(S)*B.066:*067:* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor068:* the matrix A (after equilibration if FACT = 'E') as069:*070:* A = U * D * U**T, if UPLO = 'U', or071:* A = L * D * L**T, if UPLO = 'L',072:*073:* where U (or L) is a product of permutation and unit upper (lower)074:* triangular matrices, and D is symmetric and block diagonal with075:* 1-by-1 and 2-by-2 diagonal blocks.076:*077:* 3. If some D(i,i)=0, so that D is exactly singular, then the078:* routine returns with INFO = i. Otherwise, the factored form of A079:* is used to estimate the condition number of the matrix A (see080:* argument RCOND). If the reciprocal of the condition number is081:* less than machine precision, the routine still goes on to solve082:* for X and compute error bounds as described below.083:*084:* 4. The system of equations is solved for X using the factored form085:* of A.086:*087:* 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),088:* the routine will use iterative refinement to try to get a small089:* error and error bounds. Refinement calculates the residual to at090:* least twice the working precision.091:*092:* 6. If equilibration was used, the matrix X is premultiplied by093:* diag(R) so that it solves the original system before094:* equilibration.095:*096:* Arguments097:* =========098:*099:* Some optional parameters are bundled in the PARAMS array. These100:* settings determine how refinement is performed, but often the101:* defaults are acceptable. If the defaults are acceptable, users102:* can pass NPARAMS = 0 which prevents the source code from accessing103:* the PARAMS argument.104:*105:* FACT (input) CHARACTER*1106:* Specifies whether or not the factored form of the matrix A is107:* supplied on entry, and if not, whether the matrix A should be108:* equilibrated before it is factored.109:* = 'F': On entry, AF and IPIV contain the factored form of A.110:* If EQUED is not 'N', the matrix A has been111:* equilibrated with scaling factors given by S.112:* A, AF, and IPIV are not modified.113:* = 'N': The matrix A will be copied to AF and factored.114:* = 'E': The matrix A will be equilibrated if necessary, then115:* copied to AF and factored.116:*117:* UPLO (input) CHARACTER*1118:* = 'U': Upper triangle of A is stored;119:* = 'L': Lower triangle of A is stored.120:*121:* N (input) INTEGER122:* The number of linear equations, i.e., the order of the123:* matrix A. N >= 0.124:*125:* NRHS (input) INTEGER126:* The number of right hand sides, i.e., the number of columns127:* of the matrices B and X. NRHS >= 0.128:*129:* A (input/output) REAL array, dimension (LDA,N)130:* The symmetric matrix A. If UPLO = 'U', the leading N-by-N131:* upper triangular part of A contains the upper triangular132:* part of the matrix A, and the strictly lower triangular133:* part of A is not referenced. If UPLO = 'L', the leading134:* N-by-N lower triangular part of A contains the lower135:* triangular part of the matrix A, and the strictly upper136:* triangular part of A is not referenced.137:*138:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by139:* diag(S)*A*diag(S).140:*141:* LDA (input) INTEGER142:* The leading dimension of the array A. LDA >= max(1,N).143:*144:* AF (input or output) REAL array, dimension (LDAF,N)145:* If FACT = 'F', then AF is an input argument and on entry146:* contains the block diagonal matrix D and the multipliers147:* used to obtain the factor U or L from the factorization A =148:* U*D*U**T or A = L*D*L**T as computed by SSYTRF.149:*150:* If FACT = 'N', then AF is an output argument and on exit151:* returns the block diagonal matrix D and the multipliers152:* used to obtain the factor U or L from the factorization A =153:* U*D*U**T or A = L*D*L**T.154:*155:* LDAF (input) INTEGER156:* The leading dimension of the array AF. LDAF >= max(1,N).157:*158:* IPIV (input or output) INTEGER array, dimension (N)159:* If FACT = 'F', then IPIV is an input argument and on entry160:* contains details of the interchanges and the block161:* structure of D, as determined by SSYTRF. If IPIV(k) > 0,162:* then rows and columns k and IPIV(k) were interchanged and163:* D(k,k) is a 1-by-1 diagonal block. If UPLO = 'U' and164:* IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and165:* -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2166:* diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0,167:* then rows and columns k+1 and -IPIV(k) were interchanged168:* and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.169:*170:* If FACT = 'N', then IPIV is an output argument and on exit171:* contains details of the interchanges and the block172:* structure of D, as determined by SSYTRF.173:*174:* EQUED (input or output) CHARACTER*1175:* Specifies the form of equilibration that was done.176:* = 'N': No equilibration (always true if FACT = 'N').177:* = 'Y': Both row and column equilibration, i.e., A has been178:* replaced by diag(S) * A * diag(S).179:* EQUED is an input argument if FACT = 'F'; otherwise, it is an180:* output argument.181:*182:* S (input or output) REAL array, dimension (N)183:* The scale factors for A. If EQUED = 'Y', A is multiplied on184:* the left and right by diag(S). S is an input argument if FACT =185:* 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED186:* = 'Y', each element of S must be positive. If S is output, each187:* element of S is a power of the radix. If S is input, each element188:* of S should be a power of the radix to ensure a reliable solution189:* and error estimates. Scaling by powers of the radix does not cause190:* rounding errors unless the result underflows or overflows.191:* Rounding errors during scaling lead to refining with a matrix that192:* is not equivalent to the input matrix, producing error estimates193:* that may not be reliable.194:*195:* B (input/output) REAL array, dimension (LDB,NRHS)196:* On entry, the N-by-NRHS right hand side matrix B.197:* On exit,198:* if EQUED = 'N', B is not modified;199:* if EQUED = 'Y', B is overwritten by diag(S)*B;200:*201:* LDB (input) INTEGER202:* The leading dimension of the array B. LDB >= max(1,N).203:*204:* X (output) REAL array, dimension (LDX,NRHS)205:* If INFO = 0, the N-by-NRHS solution matrix X to the original206:* system of equations. Note that A and B are modified on exit if207:* EQUED .ne. 'N', and the solution to the equilibrated system is208:* inv(diag(S))*X.209:*210:* LDX (input) INTEGER211:* The leading dimension of the array X. LDX >= max(1,N).212:*213:* RCOND (output) REAL214:* Reciprocal scaled condition number. This is an estimate of the215:* reciprocal Skeel condition number of the matrix A after216:* equilibration (if done). If this is less than the machine217:* precision (in particular, if it is zero), the matrix is singular218:* to working precision. Note that the error may still be small even219:* if this number is very small and the matrix appears ill-220:* conditioned.221:*222:* RPVGRW (output) REAL223:* Reciprocal pivot growth. On exit, this contains the reciprocal224:* pivot growth factor norm(A)/norm(U). The "max absolute element"225:* norm is used. If this is much less than 1, then the stability of226:* the LU factorization of the (equilibrated) matrix A could be poor.227:* This also means that the solution X, estimated condition numbers,228:* and error bounds could be unreliable. If factorization fails with229:* 0<INFO<=N, then this contains the reciprocal pivot growth factor230:* for the leading INFO columns of A.231:*232:* BERR (output) REAL array, dimension (NRHS)233:* Componentwise relative backward error. This is the234:* componentwise relative backward error of each solution vector X(j)235:* (i.e., the smallest relative change in any element of A or B that236:* makes X(j) an exact solution).237:*238:* N_ERR_BNDS (input) INTEGER239:* Number of error bounds to return for each right hand side240:* and each type (normwise or componentwise). See ERR_BNDS_NORM and241:* ERR_BNDS_COMP below.242:*243:* ERR_BNDS_NORM (output) REAL array, dimension (NRHS, N_ERR_BNDS)244:* For each right-hand side, this array contains information about245:* various error bounds and condition numbers corresponding to the246:* normwise relative error, which is defined as follows:247:*248:* Normwise relative error in the ith solution vector:249:* max_j (abs(XTRUE(j,i) - X(j,i)))250:* ------------------------------251:* max_j abs(X(j,i))252:*253:* The array is indexed by the type of error information as described254:* below. There currently are up to three pieces of information255:* returned.256:*257:* The first index in ERR_BNDS_NORM(i,:) corresponds to the ith258:* right-hand side.259:*260:* The second index in ERR_BNDS_NORM(:,err) contains the following261:* three fields:262:* err = 1 "Trust/don't trust" boolean. Trust the answer if the263:* reciprocal condition number is less than the threshold264:* sqrt(n) * slamch('Epsilon').265:*266:* err = 2 "Guaranteed" error bound: The estimated forward error,267:* almost certainly within a factor of 10 of the true error268:* so long as the next entry is greater than the threshold269:* sqrt(n) * slamch('Epsilon'). This error bound should only270:* be trusted if the previous boolean is true.271:*272:* err = 3 Reciprocal condition number: Estimated normwise273:* reciprocal condition number. Compared with the threshold274:* sqrt(n) * slamch('Epsilon') to determine if the error275:* estimate is "guaranteed". These reciprocal condition276:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some277:* appropriately scaled matrix Z.278:* Let Z = S*A, where S scales each row by a power of the279:* radix so all absolute row sums of Z are approximately 1.280:*281:* See Lapack Working Note 165 for further details and extra282:* cautions.283:*284:* ERR_BNDS_COMP (output) REAL array, dimension (NRHS, N_ERR_BNDS)285:* For each right-hand side, this array contains information about286:* various error bounds and condition numbers corresponding to the287:* componentwise relative error, which is defined as follows:288:*289:* Componentwise relative error in the ith solution vector:290:* abs(XTRUE(j,i) - X(j,i))291:* max_j ----------------------292:* abs(X(j,i))293:*294:* The array is indexed by the right-hand side i (on which the295:* componentwise relative error depends), and the type of error296:* information as described below. There currently are up to three297:* pieces of information returned for each right-hand side. If298:* componentwise accuracy is not requested (PARAMS(3) = 0.0), then299:* ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most300:* the first (:,N_ERR_BNDS) entries are returned.301:*302:* The first index in ERR_BNDS_COMP(i,:) corresponds to the ith303:* right-hand side.304:*305:* The second index in ERR_BNDS_COMP(:,err) contains the following306:* three fields:307:* err = 1 "Trust/don't trust" boolean. Trust the answer if the308:* reciprocal condition number is less than the threshold309:* sqrt(n) * slamch('Epsilon').310:*311:* err = 2 "Guaranteed" error bound: The estimated forward error,312:* almost certainly within a factor of 10 of the true error313:* so long as the next entry is greater than the threshold314:* sqrt(n) * slamch('Epsilon'). This error bound should only315:* be trusted if the previous boolean is true.316:*317:* err = 3 Reciprocal condition number: Estimated componentwise318:* reciprocal condition number. Compared with the threshold319:* sqrt(n) * slamch('Epsilon') to determine if the error320:* estimate is "guaranteed". These reciprocal condition321:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some322:* appropriately scaled matrix Z.323:* Let Z = S*(A*diag(x)), where x is the solution for the324:* current right-hand side and S scales each row of325:* A*diag(x) by a power of the radix so all absolute row326:* sums of Z are approximately 1.327:*328:* See Lapack Working Note 165 for further details and extra329:* cautions.330:*331:* NPARAMS (input) INTEGER332:* Specifies the number of parameters set in PARAMS. If .LE. 0, the333:* PARAMS array is never referenced and default values are used.334:*335:* PARAMS (input / output) REAL array, dimension NPARAMS336:* Specifies algorithm parameters. If an entry is .LT. 0.0, then337:* that entry will be filled with default value used for that338:* parameter. Only positions up to NPARAMS are accessed; defaults339:* are used for higher-numbered parameters.340:*341:* PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative342:* refinement or not.343:* Default: 1.0344:* = 0.0 : No refinement is performed, and no error bounds are345:* computed.346:* = 1.0 : Use the double-precision refinement algorithm,347:* possibly with doubled-single computations if the348:* compilation environment does not support DOUBLE349:* PRECISION.350:* (other values are reserved for future use)351:*352:* PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual353:* computations allowed for refinement.354:* Default: 10355:* Aggressive: Set to 100 to permit convergence using approximate356:* factorizations or factorizations other than LU. If357:* the factorization uses a technique other than358:* Gaussian elimination, the guarantees in359:* err_bnds_norm and err_bnds_comp may no longer be360:* trustworthy.361:*362:* PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code363:* will attempt to find a solution with small componentwise364:* relative error in the double-precision algorithm. Positive365:* is true, 0.0 is false.366:* Default: 1.0 (attempt componentwise convergence)367:*368:* WORK (workspace) REAL array, dimension (4*N)369:*370:* IWORK (workspace) INTEGER array, dimension (N)371:*372:* INFO (output) INTEGER373:* = 0: Successful exit. The solution to every right-hand side is374:* guaranteed.375:* < 0: If INFO = -i, the i-th argument had an illegal value376:* > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization377:* has been completed, but the factor U is exactly singular, so378:* the solution and error bounds could not be computed. RCOND = 0379:* is returned.380:* = N+J: The solution corresponding to the Jth right-hand side is381:* not guaranteed. The solutions corresponding to other right-382:* hand sides K with K > J may not be guaranteed as well, but383:* only the first such right-hand side is reported. If a small384:* componentwise error is not requested (PARAMS(3) = 0.0) then385:* the Jth right-hand side is the first with a normwise error386:* bound that is not guaranteed (the smallest J such387:* that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)388:* the Jth right-hand side is the first with either a normwise or389:* componentwise error bound that is not guaranteed (the smallest390:* J such that either ERR_BNDS_NORM(J,1) = 0.0 or391:* ERR_BNDS_COMP(J,1) = 0.0). See the definition of392:* ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information393:* about all of the right-hand sides check ERR_BNDS_NORM or394:* ERR_BNDS_COMP.395:*396:* ==================================================================397:*398:* .. Parameters ..399: REAL ZERO, ONE 400:PARAMETER( ZERO = 0.0E+0, ONE = 1.0E+0 ) 401: INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 402: INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 403: INTEGER CMP_ERR_I, PIV_GROWTH_I 404:PARAMETER( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 405: $ BERR_I = 3 ) 406:PARAMETER( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 407:PARAMETER( CMP_RCOND_I = 7, CMP_ERR_I = 8, 408: $ PIV_GROWTH_I = 9 ) 409:* ..410:* .. Local Scalars ..411:LOGICALEQUIL, NOFACT, RCEQU 412: INTEGER INFEQU, J 413: REAL AMAX, BIGNUM, SMIN, SMAX, SCOND, SMLNUM 414:* ..415:* .. External Functions ..416:EXTERNALLSAME, SLAMCH, SLA_SYRPVGRW 417:LOGICALLSAME 418: REAL SLAMCH, SLA_SYRPVGRW 419:* ..420:* .. External Subroutines ..421:EXTERNALSSYCON, SSYEQUB, SSYTRF, SSYTRS, 422: $ SLACPY, SLAQSY, XERBLA, SLASCL2, SSYRFSX 423:* ..424:* .. Intrinsic Functions ..425:INTRINSICMAX, MIN 426:* ..427:* .. Executable Statements ..428:*429: INFO = 0 430: NOFACT =LSAME( FACT, 'N' ) 431: EQUIL =LSAME( FACT, 'E' ) 432: SMLNUM =SLAMCH( 'Safe minimum' ) 433: BIGNUM = ONE / SMLNUM 434:IF( NOFACT .OR. EQUIL )THEN435: EQUED = 'N' 436: RCEQU = .FALSE. 437:ELSE438: RCEQU =LSAME( EQUED, 'Y' ) 439:ENDIF440:*441:* Default is failure. If an input parameter is wrong or442:* factorization fails, make everything look horrible. Only the443:* pivot growth is set here, the rest is initialized in SSYRFSX.444:*445: RPVGRW = ZERO 446:*447:* Test the input parameters. PARAMS is not tested until SSYRFSX.448:*449:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT. 450: $LSAME( FACT, 'F' ) )THEN451: INFO = -1 452:ELSEIF( .NOT.LSAME(UPLO, 'U') .AND. 453: $ .NOT.LSAME(UPLO, 'L') )THEN454: INFO = -2 455:ELSEIF( N.LT.0 )THEN456: INFO = -3 457:ELSEIF( NRHS.LT.0 )THEN458: INFO = -4 459:ELSEIF( LDA.LT.MAX( 1, N ) )THEN460: INFO = -6 461:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN462: INFO = -8 463:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 464: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN465: INFO = -9 466:ELSE467:IF( RCEQU )THEN468: SMIN = BIGNUM 469: SMAX = ZERO 470:DO10 J = 1, N 471: SMIN =MIN( SMIN,S( J ) ) 472: SMAX =MAX( SMAX,S( J ) ) 473: 10CONTINUE474:IF( SMIN.LE.ZERO )THEN475: INFO = -10 476:ELSEIF( N.GT.0 )THEN477: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 478:ELSE479: SCOND = ONE 480:ENDIF481:ENDIF482:IF( INFO.EQ.0 )THEN483:IF( LDB.LT.MAX( 1, N ) )THEN484: INFO = -12 485:ELSEIF( LDX.LT.MAX( 1, N ) )THEN486: INFO = -14 487:ENDIF488:ENDIF489:ENDIF490:*491:IF( INFO.NE.0 )THEN492:CALLXERBLA( 'SSYSVXX', -INFO ) 493:RETURN494:ENDIF495:*496:IF( EQUIL )THEN497:*498:* Compute row and column scalings to equilibrate the matrix A.499:*500:CALLSSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFEQU ) 501:IF( INFEQU.EQ.0 )THEN502:*503:* Equilibrate the matrix.504:*505:CALLSLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 506: RCEQU =LSAME( EQUED, 'Y' ) 507:ENDIF508:ENDIF509:*510:* Scale the right-hand side.511:*512:IF( RCEQU )CALLSLASCL2( N, NRHS, S, B, LDB ) 513:*514:IF( NOFACT .OR. EQUIL )THEN515:*516:* Compute the LU factorization of A.517:*518:CALLSLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 519:CALLSSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, 5*MAX(1,N), INFO ) 520:*521:* Return if INFO is non-zero.522:*523:IF( INFO.GT.0 )THEN524:*525:* Pivot in column INFO is exactly 0526:* Compute the reciprocal pivot growth factor of the527:* leading rank-deficient INFO columns of A.528:*529:IF( N.GT.0 ) 530: $ RPVGRW =SLA_SYRPVGRW(UPLO, N, INFO, A, LDA, AF, 531: $ LDAF, IPIV, WORK ) 532:RETURN533:ENDIF534:ENDIF535:*536:* Compute the reciprocal pivot growth factor RPVGRW.537:*538:IF( N.GT.0 ) 539: $ RPVGRW =SLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF, 540: $ IPIV, WORK ) 541:*542:* Compute the solution matrix X.543:*544:CALLSLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 545:CALLSSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 546:*547:* Use iterative refinement to improve the computed solution and548:* compute error bounds and backward error estimates for it.549:*550:CALLSSYRFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, 551: $ S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, 552: $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO ) 553:*554:* Scale solutions.555:*556:IF( RCEQU )THEN557:CALLSLASCL2( N, NRHS, S, X, LDX ) 558:ENDIF559:*560:RETURN561:*562:* End of SSYSVXX563:*564:END565: