001:SUBROUTINECPPSVX( FACT, UPLO, N, NRHS, AP, AFP, EQUED, S, B, LDB, 002: $ X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO ) 003:*004:* -- LAPACK driver routine (version 3.2) --005:* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..006:* November 2006007:*008:* .. Scalar Arguments ..009: CHARACTER EQUED, FACT, UPLO 010: INTEGER INFO, LDB, LDX, N, NRHS 011: REAL RCOND 012:* ..013:* .. Array Arguments ..014: REALBERR( * ),FERR( * ),RWORK( * ),S( * ) 015: COMPLEXAFP( * ),AP( * ),B( LDB, * ),WORK( * ), 016: $X( LDX, * ) 017:* ..018:*019:* Purpose020:* =======021:*022:* CPPSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to023:* compute the solution to a complex system of linear equations024:* A * X = B,025:* where A is an N-by-N Hermitian positive definite matrix stored in026:* packed format and X and B are N-by-NRHS matrices.027:*028:* Error bounds on the solution and a condition estimate are also029:* provided.030:*031:* Description032:* ===========033:*034:* The following steps are performed:035:*036:* 1. If FACT = 'E', real scaling factors are computed to equilibrate037:* the system:038:* diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B039:* Whether or not the system will be equilibrated depends on the040:* scaling of the matrix A, but if equilibration is used, A is041:* overwritten by diag(S)*A*diag(S) and B by diag(S)*B.042:*043:* 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to044:* factor the matrix A (after equilibration if FACT = 'E') as045:* A = U'* U , if UPLO = 'U', or046:* A = L * L', if UPLO = 'L',047:* where U is an upper triangular matrix, L is a lower triangular048:* matrix, and ' indicates conjugate transpose.049:*050:* 3. If the leading i-by-i principal minor is not positive definite,051:* then the routine returns with INFO = i. Otherwise, the factored052:* form of A is used to estimate the condition number of the matrix053:* A. If the reciprocal of the condition number is less than machine054:* precision, INFO = N+1 is returned as a warning, but the routine055:* still goes on to solve for X and compute error bounds as056:* described below.057:*058:* 4. The system of equations is solved for X using the factored form059:* of A.060:*061:* 5. Iterative refinement is applied to improve the computed solution062:* matrix and calculate error bounds and backward error estimates063:* for it.064:*065:* 6. If equilibration was used, the matrix X is premultiplied by066:* diag(S) so that it solves the original system before067:* equilibration.068:*069:* Arguments070:* =========071:*072:* FACT (input) CHARACTER*1073:* Specifies whether or not the factored form of the matrix A is074:* supplied on entry, and if not, whether the matrix A should be075:* equilibrated before it is factored.076:* = 'F': On entry, AFP contains the factored form of A.077:* If EQUED = 'Y', the matrix A has been equilibrated078:* with scaling factors given by S. AP and AFP will not079:* be modified.080:* = 'N': The matrix A will be copied to AFP and factored.081:* = 'E': The matrix A will be equilibrated if necessary, then082:* copied to AFP and factored.083:*084:* UPLO (input) CHARACTER*1085:* = 'U': Upper triangle of A is stored;086:* = 'L': Lower triangle of A is stored.087:*088:* N (input) INTEGER089:* The number of linear equations, i.e., the order of the090:* matrix A. N >= 0.091:*092:* NRHS (input) INTEGER093:* The number of right hand sides, i.e., the number of columns094:* of the matrices B and X. NRHS >= 0.095:*096:* AP (input/output) COMPLEX array, dimension (N*(N+1)/2)097:* On entry, the upper or lower triangle of the Hermitian matrix098:* A, packed columnwise in a linear array, except if FACT = 'F'099:* and EQUED = 'Y', then A must contain the equilibrated matrix100:* diag(S)*A*diag(S). The j-th column of A is stored in the101:* array AP as follows:102:* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;103:* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.104:* See below for further details. A is not modified if105:* FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.106:*107:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by108:* diag(S)*A*diag(S).109:*110:* AFP (input or output) COMPLEX array, dimension (N*(N+1)/2)111:* If FACT = 'F', then AFP is an input argument and on entry112:* contains the triangular factor U or L from the Cholesky113:* factorization A = U**H*U or A = L*L**H, in the same storage114:* format as A. If EQUED .ne. 'N', then AFP is the factored115:* form of the equilibrated matrix A.116:*117:* If FACT = 'N', then AFP is an output argument and on exit118:* returns the triangular factor U or L from the Cholesky119:* factorization A = U**H*U or A = L*L**H of the original120:* matrix A.121:*122:* If FACT = 'E', then AFP 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 equilibrated125:* matrix A (see the description of AP for the form of the126:* equilibrated matrix).127:*128:* EQUED (input or output) CHARACTER*1129:* Specifies the form of equilibration that was done.130:* = 'N': No equilibration (always true if FACT = 'N').131:* = 'Y': Equilibration was done, i.e., A has been replaced by132:* diag(S) * A * diag(S).133:* EQUED is an input argument if FACT = 'F'; otherwise, it is an134:* output argument.135:*136:* S (input or output) REAL array, dimension (N)137:* The scale factors for A; not accessed if EQUED = 'N'. S is138:* an input argument if FACT = 'F'; otherwise, S is an output139:* argument. If FACT = 'F' and EQUED = 'Y', each element of S140:* must be positive.141:*142:* B (input/output) COMPLEX array, dimension (LDB,NRHS)143:* On entry, the N-by-NRHS right hand side matrix B.144:* On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',145:* B is overwritten by diag(S) * B.146:*147:* LDB (input) INTEGER148:* The leading dimension of the array B. LDB >= max(1,N).149:*150:* X (output) COMPLEX array, dimension (LDX,NRHS)151:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to152:* the original system of equations. Note that if EQUED = 'Y',153:* A and B are modified on exit, and the solution to the154:* equilibrated system is inv(diag(S))*X.155:*156:* LDX (input) INTEGER157:* The leading dimension of the array X. LDX >= max(1,N).158:*159:* RCOND (output) REAL160:* The estimate of the reciprocal condition number of the matrix161:* A after equilibration (if done). If RCOND is less than the162:* machine precision (in particular, if RCOND = 0), the matrix163:* is singular to working precision. This condition is164:* indicated by a return code of INFO > 0.165:*166:* FERR (output) REAL array, dimension (NRHS)167:* The estimated forward error bound for each solution vector168:* X(j) (the j-th column of the solution matrix X).169:* If XTRUE is the true solution corresponding to X(j), FERR(j)170:* is an estimated upper bound for the magnitude of the largest171:* element in (X(j) - XTRUE) divided by the magnitude of the172:* largest element in X(j). The estimate is as reliable as173:* the estimate for RCOND, and is almost always a slight174:* overestimate of the true error.175:*176:* BERR (output) REAL array, dimension (NRHS)177:* The componentwise relative backward error of each solution178:* vector X(j) (i.e., the smallest relative change in179:* any element of A or B that makes X(j) an exact solution).180:*181:* WORK (workspace) COMPLEX array, dimension (2*N)182:*183:* RWORK (workspace) REAL array, dimension (N)184:*185:* INFO (output) INTEGER186:* = 0: successful exit187:* < 0: if INFO = -i, the i-th argument had an illegal value188:* > 0: if INFO = i, and i is189:* <= N: the leading minor of order i of A is190:* not positive definite, so the factorization191:* could not be completed, and the solution has not192:* been computed. RCOND = 0 is returned.193:* = N+1: U is nonsingular, but RCOND is less than machine194:* precision, meaning that the matrix is singular195:* to working precision. Nevertheless, the196:* solution and error bounds are computed because197:* there are a number of situations where the198:* computed solution can be more accurate than the199:* value of RCOND would suggest.200:*201:* Further Details202:* ===============203:*204:* The packed storage scheme is illustrated by the following example205:* when N = 4, UPLO = 'U':206:*207:* Two-dimensional storage of the Hermitian matrix A:208:*209:* a11 a12 a13 a14210:* a22 a23 a24211:* a33 a34 (aij = conjg(aji))212:* a44213:*214:* Packed storage of the upper triangle of A:215:*216:* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]217:*218:* =====================================================================219:*220:* .. Parameters ..221: REAL ZERO, ONE 222:PARAMETER( ZERO = 0.0E+0, ONE = 1.0E+0 ) 223:* ..224:* .. Local Scalars ..225:LOGICALEQUIL, NOFACT, RCEQU 226: INTEGER I, INFEQU, J 227: REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM 228:* ..229:* .. External Functions ..230:LOGICALLSAME 231: REAL CLANHP, SLAMCH 232:EXTERNALLSAME, CLANHP, SLAMCH 233:* ..234:* .. External Subroutines ..235:EXTERNALCCOPY, CLACPY, CLAQHP, CPPCON, CPPEQU, CPPRFS, 236: $ CPPTRF, CPPTRS, XERBLA 237:* ..238:* .. Intrinsic Functions ..239:INTRINSICMAX, MIN 240:* ..241:* .. Executable Statements ..242:*243: INFO = 0 244: NOFACT =LSAME( FACT, 'N' ) 245: EQUIL =LSAME( FACT, 'E' ) 246:IF( NOFACT .OR. EQUIL )THEN247: EQUED = 'N' 248: RCEQU = .FALSE. 249:ELSE250: RCEQU =LSAME( EQUED, 'Y' ) 251: SMLNUM =SLAMCH( 'Safe minimum' ) 252: BIGNUM = ONE / SMLNUM 253:ENDIF254:*255:* Test the input parameters.256:*257:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) ) 258: $THEN259: INFO = -1 260:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 261: $THEN262: INFO = -2 263:ELSEIF( N.LT.0 )THEN264: INFO = -3 265:ELSEIF( NRHS.LT.0 )THEN266: INFO = -4 267:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 268: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN269: INFO = -7 270:ELSE271:IF( RCEQU )THEN272: SMIN = BIGNUM 273: SMAX = ZERO 274:DO10 J = 1, N 275: SMIN =MIN( SMIN,S( J ) ) 276: SMAX =MAX( SMAX,S( J ) ) 277: 10CONTINUE278:IF( SMIN.LE.ZERO )THEN279: INFO = -8 280:ELSEIF( N.GT.0 )THEN281: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 282:ELSE283: SCOND = ONE 284:ENDIF285:ENDIF286:IF( INFO.EQ.0 )THEN287:IF( LDB.LT.MAX( 1, N ) )THEN288: INFO = -10 289:ELSEIF( LDX.LT.MAX( 1, N ) )THEN290: INFO = -12 291:ENDIF292:ENDIF293:ENDIF294:*295:IF( INFO.NE.0 )THEN296:CALLXERBLA( 'CPPSVX', -INFO ) 297:RETURN298:ENDIF299:*300:IF( EQUIL )THEN301:*302:* Compute row and column scalings to equilibrate the matrix A.303:*304:CALLCPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFEQU ) 305:IF( INFEQU.EQ.0 )THEN306:*307:* Equilibrate the matrix.308:*309:CALLCLAQHP( UPLO, N, AP, S, SCOND, AMAX, EQUED ) 310: RCEQU =LSAME( EQUED, 'Y' ) 311:ENDIF312:ENDIF313:*314:* Scale the right-hand side.315:*316:IF( RCEQU )THEN317:DO30 J = 1, NRHS 318:DO20 I = 1, N 319:B( I, J ) =S( I )*B( I, J ) 320: 20CONTINUE321: 30CONTINUE322:ENDIF323:*324:IF( NOFACT .OR. EQUIL )THEN325:*326:* Compute the Cholesky factorization A = U'*U or A = L*L'.327:*328:CALLCCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 ) 329:CALLCPPTRF( UPLO, N, AFP, INFO ) 330:*331:* Return if INFO is non-zero.332:*333:IF( INFO.GT.0 )THEN334: RCOND = ZERO 335:RETURN336:ENDIF337:ENDIF338:*339:* Compute the norm of the matrix A.340:*341: ANORM =CLANHP( 'I', UPLO, N, AP, RWORK ) 342:*343:* Compute the reciprocal of the condition number of A.344:*345:CALLCPPCON( UPLO, N, AFP, ANORM, RCOND, WORK, RWORK, INFO ) 346:*347:* Compute the solution matrix X.348:*349:CALLCLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 350:CALLCPPTRS( UPLO, N, NRHS, AFP, X, LDX, INFO ) 351:*352:* Use iterative refinement to improve the computed solution and353:* compute error bounds and backward error estimates for it.354:*355:CALLCPPRFS( UPLO, N, NRHS, AP, AFP, B, LDB, X, LDX, FERR, BERR, 356: $ WORK, RWORK, INFO ) 357:*358:* Transform the solution matrix X to a solution of the original359:* system.360:*361:IF( RCEQU )THEN362:DO50 J = 1, NRHS 363:DO40 I = 1, N 364:X( I, J ) =S( I )*X( I, J ) 365: 40CONTINUE366: 50CONTINUE367:DO60 J = 1, NRHS 368:FERR( J ) =FERR( J ) / SCOND 369: 60CONTINUE370:ENDIF371:*372:* Set INFO = N+1 if the matrix is singular to working precision.373:*374:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 375: $ INFO = N + 1 376:*377:RETURN378:*379:* End of CPPSVX380:*381:END382: