001:SUBROUTINECHESVXX( 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, RWORK, 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: INTEGERIPIV( * ) 024: COMPLEXA( LDA, * ),AF( LDAF, * ),B( LDB, * ), 025: $WORK( * ),X( LDX, * ) 026: REALS( * ),PARAMS( * ),BERR( * ),RWORK( * ), 027: $ERR_BNDS_NORM( NRHS, * ), 028: $ERR_BNDS_COMP( NRHS, * ) 029:* ..030:*031:* Purpose032:* =======033:*034:* CHESVXX uses the diagonal pivoting factorization to compute the035:* solution to a complex system of linear equations A * X = B, where036:* A is an N-by-N symmetric matrix and X and B are N-by-NRHS037:* matrices.038:*039:* If requested, both normwise and maximum componentwise error bounds040:* are returned. CHESVXX 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:* CHESVXX 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:* CHESVXX 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 CHESVXX 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 LU decomposition is used to factor069:* the matrix A (after equilibration if FACT = 'E') as070:*071:* A = U * D * U**T, if UPLO = 'U', or072:* A = L * D * L**T, if UPLO = 'L',073:*074:* where U (or L) is a product of permutation and unit upper (lower)075:* triangular matrices, and D is symmetric and block diagonal with076:* 1-by-1 and 2-by-2 diagonal blocks.077:*078:* 3. If some D(i,i)=0, so that D is exactly singular, then the079:* routine returns with INFO = i. Otherwise, the factored form of A080:* is used to estimate the condition number of the matrix A (see081:* argument RCOND). If the reciprocal of the condition number is082:* less than machine precision, the routine still goes on to solve083:* for X and compute error bounds as described below.084:*085:* 4. The system of equations is solved for X using the factored form086:* of A.087:*088:* 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),089:* the routine will use iterative refinement to try to get a small090:* error and error bounds. Refinement calculates the residual to at091:* least twice the working precision.092:*093:* 6. If equilibration was used, the matrix X is premultiplied by094:* diag(R) so that it solves the original system before095:* equilibration.096:*097:* Arguments098:* =========099:*100:* Some optional parameters are bundled in the PARAMS array. These101:* settings determine how refinement is performed, but often the102:* defaults are acceptable. If the defaults are acceptable, users103:* can pass NPARAMS = 0 which prevents the source code from accessing104:* the PARAMS argument.105:*106:* FACT (input) CHARACTER*1107:* Specifies whether or not the factored form of the matrix A is108:* supplied on entry, and if not, whether the matrix A should be109:* equilibrated before it is factored.110:* = 'F': On entry, AF and IPIV contain the factored form of A.111:* If EQUED is not 'N', the matrix A has been112:* equilibrated with scaling factors given by S.113:* A, AF, and IPIV are not modified.114:* = 'N': The matrix A will be copied to AF and factored.115:* = 'E': The matrix A will be equilibrated if necessary, then116:* copied to AF and factored.117:*118:* N (input) INTEGER119:* The number of linear equations, i.e., the order of the120:* matrix A. N >= 0.121:*122:* NRHS (input) INTEGER123:* The number of right hand sides, i.e., the number of columns124:* of the matrices B and X. NRHS >= 0.125:*126:* A (input/output) COMPLEX array, dimension (LDA,N)127:* The symmetric matrix A. If UPLO = 'U', the leading N-by-N128:* upper triangular part of A contains the upper triangular129:* part of the matrix A, and the strictly lower triangular130:* part of A is not referenced. If UPLO = 'L', the leading131:* N-by-N lower triangular part of A contains the lower132:* triangular part of the matrix A, and the strictly upper133:* triangular part of A is not referenced.134:*135:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by136:* diag(S)*A*diag(S).137:*138:* LDA (input) INTEGER139:* The leading dimension of the array A. LDA >= max(1,N).140:*141:* AF (input or output) COMPLEX array, dimension (LDAF,N)142:* If FACT = 'F', then AF is an input argument and on entry143:* contains the block diagonal matrix D and the multipliers144:* used to obtain the factor U or L from the factorization A =145:* U*D*U**T or A = L*D*L**T as computed by SSYTRF.146:*147:* If FACT = 'N', then AF is an output argument and on exit148:* returns the block diagonal matrix D and the multipliers149:* used to obtain the factor U or L from the factorization A =150:* U*D*U**T or A = L*D*L**T.151:*152:* LDAF (input) INTEGER153:* The leading dimension of the array AF. LDAF >= max(1,N).154:*155:* IPIV (input or output) INTEGER array, dimension (N)156:* If FACT = 'F', then IPIV is an input argument and on entry157:* contains details of the interchanges and the block158:* structure of D, as determined by SSYTRF. If IPIV(k) > 0,159:* then rows and columns k and IPIV(k) were interchanged and160:* D(k,k) is a 1-by-1 diagonal block. If UPLO = 'U' and161:* IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and162:* -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2163:* diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0,164:* then rows and columns k+1 and -IPIV(k) were interchanged165:* and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.166:*167:* If FACT = 'N', then IPIV is an output argument and on exit168:* contains details of the interchanges and the block169:* structure of D, as determined by SSYTRF.170:*171:* EQUED (input or output) CHARACTER*1172:* Specifies the form of equilibration that was done.173:* = 'N': No equilibration (always true if FACT = 'N').174:* = 'Y': Both row and column equilibration, i.e., A has been175:* replaced by diag(S) * A * diag(S).176:* EQUED is an input argument if FACT = 'F'; otherwise, it is an177:* output argument.178:*179:* S (input or output) REAL array, dimension (N)180:* The scale factors for A. If EQUED = 'Y', A is multiplied on181:* the left and right by diag(S). S is an input argument if FACT =182:* 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED183:* = 'Y', each element of S must be positive. If S is output, each184:* element of S is a power of the radix. If S is input, each element185:* of S should be a power of the radix to ensure a reliable solution186:* and error estimates. Scaling by powers of the radix does not cause187:* rounding errors unless the result underflows or overflows.188:* Rounding errors during scaling lead to refining with a matrix that189:* is not equivalent to the input matrix, producing error estimates190:* that may not be reliable.191:*192:* B (input/output) COMPLEX array, dimension (LDB,NRHS)193:* On entry, the N-by-NRHS right hand side matrix B.194:* On exit,195:* if EQUED = 'N', B is not modified;196:* if EQUED = 'Y', B is overwritten by diag(S)*B;197:*198:* LDB (input) INTEGER199:* The leading dimension of the array B. LDB >= max(1,N).200:*201:* X (output) COMPLEX array, dimension (LDX,NRHS)202:* If INFO = 0, the N-by-NRHS solution matrix X to the original203:* system of equations. Note that A and B are modified on exit if204:* EQUED .ne. 'N', and the solution to the equilibrated system is205:* inv(diag(S))*X.206:*207:* LDX (input) INTEGER208:* The leading dimension of the array X. LDX >= max(1,N).209:*210:* RCOND (output) REAL211:* Reciprocal scaled condition number. This is an estimate of the212:* reciprocal Skeel condition number of the matrix A after213:* equilibration (if done). If this is less than the machine214:* precision (in particular, if it is zero), the matrix is singular215:* to working precision. Note that the error may still be small even216:* if this number is very small and the matrix appears ill-217:* conditioned.218:*219:* RPVGRW (output) REAL220:* Reciprocal pivot growth. On exit, this contains the reciprocal221:* pivot growth factor norm(A)/norm(U). The "max absolute element"222:* norm is used. If this is much less than 1, then the stability of223:* the LU factorization of the (equilibrated) matrix A could be poor.224:* This also means that the solution X, estimated condition numbers,225:* and error bounds could be unreliable. If factorization fails with226:* 0<INFO<=N, then this contains the reciprocal pivot growth factor227:* for the leading INFO columns of A.228:*229:* BERR (output) REAL array, dimension (NRHS)230:* Componentwise relative backward error. This is the231:* componentwise relative backward error of each solution vector X(j)232:* (i.e., the smallest relative change in any element of A or B that233:* makes X(j) an exact solution).234:*235:* N_ERR_BNDS (input) INTEGER236:* Number of error bounds to return for each right hand side237:* and each type (normwise or componentwise). See ERR_BNDS_NORM and238:* ERR_BNDS_COMP below.239:*240:* ERR_BNDS_NORM (output) REAL array, dimension (NRHS, N_ERR_BNDS)241:* For each right-hand side, this array contains information about242:* various error bounds and condition numbers corresponding to the243:* normwise relative error, which is defined as follows:244:*245:* Normwise relative error in the ith solution vector:246:* max_j (abs(XTRUE(j,i) - X(j,i)))247:* ------------------------------248:* max_j abs(X(j,i))249:*250:* The array is indexed by the type of error information as described251:* below. There currently are up to three pieces of information252:* returned.253:*254:* The first index in ERR_BNDS_NORM(i,:) corresponds to the ith255:* right-hand side.256:*257:* The second index in ERR_BNDS_NORM(:,err) contains the following258:* three fields:259:* err = 1 "Trust/don't trust" boolean. Trust the answer if the260:* reciprocal condition number is less than the threshold261:* sqrt(n) * slamch('Epsilon').262:*263:* err = 2 "Guaranteed" error bound: The estimated forward error,264:* almost certainly within a factor of 10 of the true error265:* so long as the next entry is greater than the threshold266:* sqrt(n) * slamch('Epsilon'). This error bound should only267:* be trusted if the previous boolean is true.268:*269:* err = 3 Reciprocal condition number: Estimated normwise270:* reciprocal condition number. Compared with the threshold271:* sqrt(n) * slamch('Epsilon') to determine if the error272:* estimate is "guaranteed". These reciprocal condition273:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some274:* appropriately scaled matrix Z.275:* Let Z = S*A, where S scales each row by a power of the276:* radix so all absolute row sums of Z are approximately 1.277:*278:* See Lapack Working Note 165 for further details and extra279:* cautions.280:*281:* ERR_BNDS_COMP (output) REAL array, dimension (NRHS, N_ERR_BNDS)282:* For each right-hand side, this array contains information about283:* various error bounds and condition numbers corresponding to the284:* componentwise relative error, which is defined as follows:285:*286:* Componentwise relative error in the ith solution vector:287:* abs(XTRUE(j,i) - X(j,i))288:* max_j ----------------------289:* abs(X(j,i))290:*291:* The array is indexed by the right-hand side i (on which the292:* componentwise relative error depends), and the type of error293:* information as described below. There currently are up to three294:* pieces of information returned for each right-hand side. If295:* componentwise accuracy is not requested (PARAMS(3) = 0.0), then296:* ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most297:* the first (:,N_ERR_BNDS) entries are returned.298:*299:* The first index in ERR_BNDS_COMP(i,:) corresponds to the ith300:* right-hand side.301:*302:* The second index in ERR_BNDS_COMP(:,err) contains the following303:* three fields:304:* err = 1 "Trust/don't trust" boolean. Trust the answer if the305:* reciprocal condition number is less than the threshold306:* sqrt(n) * slamch('Epsilon').307:*308:* err = 2 "Guaranteed" error bound: The estimated forward error,309:* almost certainly within a factor of 10 of the true error310:* so long as the next entry is greater than the threshold311:* sqrt(n) * slamch('Epsilon'). This error bound should only312:* be trusted if the previous boolean is true.313:*314:* err = 3 Reciprocal condition number: Estimated componentwise315:* reciprocal condition number. Compared with the threshold316:* sqrt(n) * slamch('Epsilon') to determine if the error317:* estimate is "guaranteed". These reciprocal condition318:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some319:* appropriately scaled matrix Z.320:* Let Z = S*(A*diag(x)), where x is the solution for the321:* current right-hand side and S scales each row of322:* A*diag(x) by a power of the radix so all absolute row323:* sums of Z are approximately 1.324:*325:* See Lapack Working Note 165 for further details and extra326:* cautions.327:*328:* NPARAMS (input) INTEGER329:* Specifies the number of parameters set in PARAMS. If .LE. 0, the330:* PARAMS array is never referenced and default values are used.331:*332:* PARAMS (input / output) REAL array, dimension NPARAMS333:* Specifies algorithm parameters. If an entry is .LT. 0.0, then334:* that entry will be filled with default value used for that335:* parameter. Only positions up to NPARAMS are accessed; defaults336:* are used for higher-numbered parameters.337:*338:* PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative339:* refinement or not.340:* Default: 1.0341:* = 0.0 : No refinement is performed, and no error bounds are342:* computed.343:* = 1.0 : Use the double-precision refinement algorithm,344:* possibly with doubled-single computations if the345:* compilation environment does not support DOUBLE346:* PRECISION.347:* (other values are reserved for future use)348:*349:* PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual350:* computations allowed for refinement.351:* Default: 10352:* Aggressive: Set to 100 to permit convergence using approximate353:* factorizations or factorizations other than LU. If354:* the factorization uses a technique other than355:* Gaussian elimination, the guarantees in356:* err_bnds_norm and err_bnds_comp may no longer be357:* trustworthy.358:*359:* PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code360:* will attempt to find a solution with small componentwise361:* relative error in the double-precision algorithm. Positive362:* is true, 0.0 is false.363:* Default: 1.0 (attempt componentwise convergence)364:*365:* WORK (workspace) COMPLEX array, dimension (2*N)366:*367:* RWORK (workspace) REAL array, dimension (3*N)368:*369:* INFO (output) INTEGER370:* = 0: Successful exit. The solution to every right-hand side is371:* guaranteed.372:* < 0: If INFO = -i, the i-th argument had an illegal value373:* > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization374:* has been completed, but the factor U is exactly singular, so375:* the solution and error bounds could not be computed. RCOND = 0376:* is returned.377:* = N+J: The solution corresponding to the Jth right-hand side is378:* not guaranteed. The solutions corresponding to other right-379:* hand sides K with K > J may not be guaranteed as well, but380:* only the first such right-hand side is reported. If a small381:* componentwise error is not requested (PARAMS(3) = 0.0) then382:* the Jth right-hand side is the first with a normwise error383:* bound that is not guaranteed (the smallest J such384:* that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)385:* the Jth right-hand side is the first with either a normwise or386:* componentwise error bound that is not guaranteed (the smallest387:* J such that either ERR_BNDS_NORM(J,1) = 0.0 or388:* ERR_BNDS_COMP(J,1) = 0.0). See the definition of389:* ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information390:* about all of the right-hand sides check ERR_BNDS_NORM or391:* ERR_BNDS_COMP.392:*393:* ==================================================================394:*395:* .. Parameters ..396: REAL ZERO, ONE 397:PARAMETER( ZERO = 0.0E+0, ONE = 1.0E+0 ) 398: INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 399: INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 400: INTEGER CMP_ERR_I, PIV_GROWTH_I 401:PARAMETER( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 402: $ BERR_I = 3 ) 403:PARAMETER( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 404:PARAMETER( CMP_RCOND_I = 7, CMP_ERR_I = 8, 405: $ PIV_GROWTH_I = 9 ) 406:* ..407:* .. Local Scalars ..408:LOGICALEQUIL, NOFACT, RCEQU 409: INTEGER INFEQU, J 410: REAL AMAX, BIGNUM, SMIN, SMAX, SCOND, SMLNUM 411:* ..412:* .. External Functions ..413:EXTERNALLSAME, SLAMCH, CLA_HERPVGRW 414:LOGICALLSAME 415: REAL SLAMCH, CLA_HERPVGRW 416:* ..417:* .. External Subroutines ..418:EXTERNALCHECON, CHEEQUB, CHETRF, CHETRS, CLACPY, 419: $ CLAQHE, XERBLA, CLASCL2, CHERFSX 420:* ..421:* .. Intrinsic Functions ..422:INTRINSICMAX, MIN 423:* ..424:* .. Executable Statements ..425:*426: INFO = 0 427: NOFACT =LSAME( FACT, 'N' ) 428: EQUIL =LSAME( FACT, 'E' ) 429: SMLNUM =SLAMCH( 'Safe minimum' ) 430: BIGNUM = ONE / SMLNUM 431:IF( NOFACT .OR. EQUIL )THEN432: EQUED = 'N' 433: RCEQU = .FALSE. 434:ELSE435: RCEQU =LSAME( EQUED, 'Y' ) 436:ENDIF437:*438:* Default is failure. If an input parameter is wrong or439:* factorization fails, make everything look horrible. Only the440:* pivot growth is set here, the rest is initialized in CHERFSX.441:*442: RPVGRW = ZERO 443:*444:* Test the input parameters. PARAMS is not tested until CHERFSX.445:*446:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT. 447: $LSAME( FACT, 'F' ) )THEN448: INFO = -1 449:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. 450: $ .NOT.LSAME( UPLO, 'L' ) )THEN451: INFO = -2 452:ELSEIF( N.LT.0 )THEN453: INFO = -3 454:ELSEIF( NRHS.LT.0 )THEN455: INFO = -4 456:ELSEIF( LDA.LT.MAX( 1, N ) )THEN457: INFO = -6 458:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN459: INFO = -8 460:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 461: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN462: INFO = -9 463:ELSE464:IF( RCEQU )THEN465: SMIN = BIGNUM 466: SMAX = ZERO 467:DO10 J = 1, N 468: SMIN =MIN( SMIN,S( J ) ) 469: SMAX =MAX( SMAX,S( J ) ) 470: 10CONTINUE471:IF( SMIN.LE.ZERO )THEN472: INFO = -10 473:ELSEIF( N.GT.0 )THEN474: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 475:ELSE476: SCOND = ONE 477:ENDIF478:ENDIF479:IF( INFO.EQ.0 )THEN480:IF( LDB.LT.MAX( 1, N ) )THEN481: INFO = -12 482:ELSEIF( LDX.LT.MAX( 1, N ) )THEN483: INFO = -14 484:ENDIF485:ENDIF486:ENDIF487:*488:IF( INFO.NE.0 )THEN489:CALLXERBLA( 'CHESVXX', -INFO ) 490:RETURN491:ENDIF492:*493:IF( EQUIL )THEN494:*495:* Compute row and column scalings to equilibrate the matrix A.496:*497:CALLCHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFEQU ) 498:IF( INFEQU.EQ.0 )THEN499:*500:* Equilibrate the matrix.501:*502:CALLCLAQHE( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 503: RCEQU =LSAME( EQUED, 'Y' ) 504:ENDIF505:ENDIF506:*507:* Scale the right-hand side.508:*509:IF( RCEQU )CALLCLASCL2( N, NRHS, S, B, LDB ) 510:*511:IF( NOFACT .OR. EQUIL )THEN512:*513:* Compute the LU factorization of A.514:*515:CALLCLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 516:CALLCHETRF( UPLO, N, AF, LDAF, IPIV, WORK, 5*MAX(1,N), INFO ) 517:*518:* Return if INFO is non-zero.519:*520:IF( INFO.GT.0 )THEN521:*522:* Pivot in column INFO is exactly 0523:* Compute the reciprocal pivot growth factor of the524:* leading rank-deficient INFO columns of A.525:*526:IF( N.GT.0 ) 527: $ RPVGRW =CLA_HERPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF, 528: $ IPIV, WORK ) 529:RETURN530:ENDIF531:ENDIF532:*533:* Compute the reciprocal pivot growth factor RPVGRW.534:*535:IF( N.GT.0 ) 536: $ RPVGRW =CLA_HERPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF, IPIV, 537: $ WORK ) 538:*539:* Compute the solution matrix X.540:*541:CALLCLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 542:CALLCHETRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 543:*544:* Use iterative refinement to improve the computed solution and545:* compute error bounds and backward error estimates for it.546:*547:CALLCHERFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, 548: $ S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, 549: $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO ) 550:*551:* Scale solutions.552:*553:IF( RCEQU )THEN554:CALLCLASCL2( N, NRHS, S, X, LDX ) 555:ENDIF556:*557:RETURN558:*559:* End of CHESVXX560:*561:END562: