001:SUBROUTINEZPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, 002: $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, 003: $ RWORK, INFO ) 004:*005:* -- LAPACK driver routine (version 3.2) --006:* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..007:* November 2006008:*009:* .. Scalar Arguments ..010: CHARACTER EQUED, FACT, UPLO 011: INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS 012: DOUBLE PRECISION RCOND 013:* ..014:* .. Array Arguments ..015: DOUBLE PRECISIONBERR( * ),FERR( * ),RWORK( * ),S( * ) 016: COMPLEX*16A( LDA, * ),AF( LDAF, * ),B( LDB, * ), 017: $WORK( * ),X( LDX, * ) 018:* ..019:*020:* Purpose021:* =======022:*023:* ZPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to024:* compute the solution to a complex system of linear equations025:* A * X = B,026:* where A is an N-by-N Hermitian positive definite matrix and X and B027:* are N-by-NRHS matrices.028:*029:* Error bounds on the solution and a condition estimate are also030:* provided.031:*032:* Description033:* ===========034:*035:* The following steps are performed:036:*037:* 1. If FACT = 'E', real scaling factors are computed to equilibrate038:* the system:039:* diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B040:* Whether or not the system will be equilibrated depends on the041:* scaling of the matrix A, but if equilibration is used, A is042:* overwritten by diag(S)*A*diag(S) and B by diag(S)*B.043:*044:* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to045:* factor the matrix A (after equilibration if FACT = 'E') as046:* A = U**H* U, if UPLO = 'U', or047:* A = L * L**H, if UPLO = 'L',048:* where U is an upper triangular matrix and L is a lower triangular049:* matrix.050:*051:* 3. If the leading i-by-i principal minor is not positive definite,052:* then the routine returns with INFO = i. Otherwise, the factored053:* form of A is used to estimate the condition number of the matrix054:* A. If the reciprocal of the condition number is less than machine055:* precision, INFO = N+1 is returned as a warning, but the routine056:* still goes on to solve for X and compute error bounds as057:* described below.058:*059:* 4. The system of equations is solved for X using the factored form060:* of A.061:*062:* 5. Iterative refinement is applied to improve the computed solution063:* matrix and calculate error bounds and backward error estimates064:* for it.065:*066:* 6. If equilibration was used, the matrix X is premultiplied by067:* diag(S) so that it solves the original system before068:* equilibration.069:*070:* Arguments071:* =========072:*073:* FACT (input) CHARACTER*1074:* Specifies whether or not the factored form of the matrix A is075:* supplied on entry, and if not, whether the matrix A should be076:* equilibrated before it is factored.077:* = 'F': On entry, AF contains the factored form of A.078:* If EQUED = 'Y', the matrix A has been equilibrated079:* with scaling factors given by S. A and AF will not080:* be modified.081:* = 'N': The matrix A will be copied to AF and factored.082:* = 'E': The matrix A will be equilibrated if necessary, then083:* copied to AF and factored.084:*085:* UPLO (input) CHARACTER*1086:* = 'U': Upper triangle of A is stored;087:* = 'L': Lower triangle of A is stored.088:*089:* N (input) INTEGER090:* The number of linear equations, i.e., the order of the091:* matrix A. N >= 0.092:*093:* NRHS (input) INTEGER094:* The number of right hand sides, i.e., the number of columns095:* of the matrices B and X. NRHS >= 0.096:*097:* A (input/output) COMPLEX*16 array, dimension (LDA,N)098:* On entry, the Hermitian matrix A, except if FACT = 'F' and099:* EQUED = 'Y', then A must contain the equilibrated matrix100:* diag(S)*A*diag(S). If UPLO = 'U', the leading101:* N-by-N upper triangular part of A contains the upper102:* triangular part of the matrix A, and the strictly lower103:* triangular part of A is not referenced. If UPLO = 'L', the104:* leading N-by-N lower triangular part of A contains the lower105:* triangular part of the matrix A, and the strictly upper106:* triangular part of A is not referenced. A is not modified if107:* FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.108:*109:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by110:* diag(S)*A*diag(S).111:*112:* LDA (input) INTEGER113:* The leading dimension of the array A. LDA >= max(1,N).114:*115:* AF (input or output) COMPLEX*16 array, dimension (LDAF,N)116:* If FACT = 'F', then AF is an input argument and on entry117:* contains the triangular factor U or L from the Cholesky118:* factorization A = U**H*U or A = L*L**H, in the same storage119:* format as A. If EQUED .ne. 'N', then AF is the factored form120:* of the equilibrated matrix diag(S)*A*diag(S).121:*122:* If FACT = 'N', then AF is an output argument and on exit123:* returns the triangular factor U or L from the Cholesky124:* factorization A = U**H*U or A = L*L**H of the original125:* matrix A.126:*127:* If FACT = 'E', then AF is an output argument and on exit128:* returns the triangular factor U or L from the Cholesky129:* factorization A = U**H*U or A = L*L**H of the equilibrated130:* matrix A (see the description of A for the form of the131:* equilibrated matrix).132:*133:* LDAF (input) INTEGER134:* The leading dimension of the array AF. LDAF >= max(1,N).135:*136:* EQUED (input or output) CHARACTER*1137:* Specifies the form of equilibration that was done.138:* = 'N': No equilibration (always true if FACT = 'N').139:* = 'Y': Equilibration was done, i.e., A has been replaced by140:* diag(S) * A * diag(S).141:* EQUED is an input argument if FACT = 'F'; otherwise, it is an142:* output argument.143:*144:* S (input or output) DOUBLE PRECISION array, dimension (N)145:* The scale factors for A; not accessed if EQUED = 'N'. S is146:* an input argument if FACT = 'F'; otherwise, S is an output147:* argument. If FACT = 'F' and EQUED = 'Y', each element of S148:* must be positive.149:*150:* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)151:* On entry, the N-by-NRHS righthand side matrix B.152:* On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',153:* B is overwritten by diag(S) * B.154:*155:* LDB (input) INTEGER156:* The leading dimension of the array B. LDB >= max(1,N).157:*158:* X (output) COMPLEX*16 array, dimension (LDX,NRHS)159:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to160:* the original system of equations. Note that if EQUED = 'Y',161:* A and B are modified on exit, and the solution to the162:* equilibrated system is inv(diag(S))*X.163:*164:* LDX (input) INTEGER165:* The leading dimension of the array X. LDX >= max(1,N).166:*167:* RCOND (output) DOUBLE PRECISION168:* The estimate of the reciprocal condition number of the matrix169:* A after equilibration (if done). If RCOND is less than the170:* machine precision (in particular, if RCOND = 0), the matrix171:* is singular to working precision. This condition is172:* indicated by a return code of INFO > 0.173:*174:* FERR (output) DOUBLE PRECISION array, dimension (NRHS)175:* The estimated forward error bound for each solution vector176:* X(j) (the j-th column of the solution matrix X).177:* If XTRUE is the true solution corresponding to X(j), FERR(j)178:* is an estimated upper bound for the magnitude of the largest179:* element in (X(j) - XTRUE) divided by the magnitude of the180:* largest element in X(j). The estimate is as reliable as181:* the estimate for RCOND, and is almost always a slight182:* overestimate of the true error.183:*184:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)185:* The componentwise relative backward error of each solution186:* vector X(j) (i.e., the smallest relative change in187:* any element of A or B that makes X(j) an exact solution).188:*189:* WORK (workspace) COMPLEX*16 array, dimension (2*N)190:*191:* RWORK (workspace) DOUBLE PRECISION array, dimension (N)192:*193:* INFO (output) INTEGER194:* = 0: successful exit195:* < 0: if INFO = -i, the i-th argument had an illegal value196:* > 0: if INFO = i, and i is197:* <= N: the leading minor of order i of A is198:* not positive definite, so the factorization199:* could not be completed, and the solution has not200:* been computed. RCOND = 0 is returned.201:* = N+1: U is nonsingular, but RCOND is less than machine202:* precision, meaning that the matrix is singular203:* to working precision. Nevertheless, the204:* solution and error bounds are computed because205:* there are a number of situations where the206:* computed solution can be more accurate than the207:* value of RCOND would suggest.208:*209:* =====================================================================210:*211:* .. Parameters ..212: DOUBLE PRECISION ZERO, ONE 213:PARAMETER( ZERO = 0.0D+0, ONE = 1.0D+0 ) 214:* ..215:* .. Local Scalars ..216:LOGICALEQUIL, NOFACT, RCEQU 217: INTEGER I, INFEQU, J 218: DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 219:* ..220:* .. External Functions ..221:LOGICALLSAME 222: DOUBLE PRECISION DLAMCH, ZLANHE 223:EXTERNALLSAME, DLAMCH, ZLANHE 224:* ..225:* .. External Subroutines ..226:EXTERNALXERBLA, ZLACPY, ZLAQHE, ZPOCON, ZPOEQU, ZPORFS, 227: $ ZPOTRF, ZPOTRS 228:* ..229:* .. Intrinsic Functions ..230:INTRINSICMAX, MIN 231:* ..232:* .. Executable Statements ..233:*234: INFO = 0 235: NOFACT =LSAME( FACT, 'N' ) 236: EQUIL =LSAME( FACT, 'E' ) 237:IF( NOFACT .OR. EQUIL )THEN238: EQUED = 'N' 239: RCEQU = .FALSE. 240:ELSE241: RCEQU =LSAME( EQUED, 'Y' ) 242: SMLNUM =DLAMCH( 'Safe minimum' ) 243: BIGNUM = ONE / SMLNUM 244:ENDIF245:*246:* Test the input parameters.247:*248:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 249: $THEN250: INFO = -1 251:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 252: $THEN253: INFO = -2 254:ELSEIF( N.LT.0 )THEN255: INFO = -3 256:ELSEIF( NRHS.LT.0 )THEN257: INFO = -4 258:ELSEIF( LDA.LT.MAX( 1, N ) )THEN259: INFO = -6 260:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN261: INFO = -8 262:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 263: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN264: INFO = -9 265:ELSE266:IF( RCEQU )THEN267: SMIN = BIGNUM 268: SMAX = ZERO 269:DO10 J = 1, N 270: SMIN =MIN( SMIN,S( J ) ) 271: SMAX =MAX( SMAX,S( J ) ) 272: 10CONTINUE273:IF( SMIN.LE.ZERO )THEN274: INFO = -10 275:ELSEIF( N.GT.0 )THEN276: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 277:ELSE278: SCOND = ONE 279:ENDIF280:ENDIF281:IF( INFO.EQ.0 )THEN282:IF( LDB.LT.MAX( 1, N ) )THEN283: INFO = -12 284:ELSEIF( LDX.LT.MAX( 1, N ) )THEN285: INFO = -14 286:ENDIF287:ENDIF288:ENDIF289:*290:IF( INFO.NE.0 )THEN291:CALLXERBLA( 'ZPOSVX', -INFO ) 292:RETURN293:ENDIF294:*295:IF( EQUIL )THEN296:*297:* Compute row and column scalings to equilibrate the matrix A.298:*299:CALLZPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU ) 300:IF( INFEQU.EQ.0 )THEN301:*302:* Equilibrate the matrix.303:*304:CALLZLAQHE( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 305: RCEQU =LSAME( EQUED, 'Y' ) 306:ENDIF307:ENDIF308:*309:* Scale the right hand side.310:*311:IF( RCEQU )THEN312:DO30 J = 1, NRHS 313:DO20 I = 1, N 314:B( I, J ) =S( I )*B( I, J ) 315: 20CONTINUE316: 30CONTINUE317:ENDIF318:*319:IF( NOFACT .OR. EQUIL )THEN320:*321:* Compute the Cholesky factorization A = U'*U or A = L*L'.322:*323:CALLZLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 324:CALLZPOTRF( UPLO, N, AF, LDAF, INFO ) 325:*326:* Return if INFO is non-zero.327:*328:IF( INFO.GT.0 )THEN329: RCOND = ZERO 330:RETURN331:ENDIF332:ENDIF333:*334:* Compute the norm of the matrix A.335:*336: ANORM =ZLANHE( '1', UPLO, N, A, LDA, RWORK ) 337:*338:* Compute the reciprocal of the condition number of A.339:*340:CALLZPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO ) 341:*342:* Compute the solution matrix X.343:*344:CALLZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 345:CALLZPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO ) 346:*347:* Use iterative refinement to improve the computed solution and348:* compute error bounds and backward error estimates for it.349:*350:CALLZPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, 351: $ FERR, BERR, WORK, RWORK, INFO ) 352:*353:* Transform the solution matrix X to a solution of the original354:* system.355:*356:IF( RCEQU )THEN357:DO50 J = 1, NRHS 358:DO40 I = 1, N 359:X( I, J ) =S( I )*X( I, J ) 360: 40CONTINUE361: 50CONTINUE362:DO60 J = 1, NRHS 363:FERR( J ) =FERR( J ) / SCOND 364: 60CONTINUE365:ENDIF366:*367:* Set INFO = N+1 if the matrix is singular to working precision.368:*369:IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 370: $ INFO = N + 1 371:*372:RETURN373:*374:* End of ZPOSVX375:*376:END377: