001:SUBROUTINECSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, 002: $ LDX, RCOND, FERR, BERR, WORK, RWORK, INFO ) 003:*004:* -- LAPACK driver routine (version 3.2) --005:* -- LAPACK is a software package provided by Univ. of Tennessee, --006:* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--007:* November 2006008:*009:* .. Scalar Arguments ..010: CHARACTER FACT, UPLO 011: INTEGER INFO, LDB, LDX, N, NRHS 012: REAL RCOND 013:* ..014:* .. Array Arguments ..015: INTEGERIPIV( * ) 016: REALBERR( * ),FERR( * ),RWORK( * ) 017: COMPLEXAFP( * ),AP( * ),B( LDB, * ),WORK( * ), 018: $X( LDX, * ) 019:* ..020:*021:* Purpose022:* =======023:*024:* CSPSVX uses the diagonal pivoting factorization A = U*D*U**T or025:* A = L*D*L**T to compute the solution to a complex system of linear026:* equations A * X = B, where A is an N-by-N symmetric matrix stored027:* in packed format and X and B 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 = 'N', the diagonal pivoting method is used to factor A as038:* A = U * D * U**T, if UPLO = 'U', or039:* A = L * D * L**T, if UPLO = 'L',040:* where U (or L) is a product of permutation and unit upper (lower)041:* triangular matrices and D is symmetric and block diagonal with042:* 1-by-1 and 2-by-2 diagonal blocks.043:*044:* 2. If some D(i,i)=0, so that D is exactly singular, then the routine045:* returns with INFO = i. Otherwise, the factored form of A is used046:* to estimate the condition number of the matrix A. If the047:* reciprocal of the condition number is less than machine precision,048:* INFO = N+1 is returned as a warning, but the routine still goes on049:* to solve for X and compute error bounds as described below.050:*051:* 3. The system of equations is solved for X using the factored form052:* of A.053:*054:* 4. Iterative refinement is applied to improve the computed solution055:* matrix and calculate error bounds and backward error estimates056:* for it.057:*058:* Arguments059:* =========060:*061:* FACT (input) CHARACTER*1062:* Specifies whether or not the factored form of A has been063:* supplied on entry.064:* = 'F': On entry, AFP and IPIV contain the factored form065:* of A. AP, AFP and IPIV will not be modified.066:* = 'N': The matrix A will be copied to AFP and factored.067:*068:* UPLO (input) CHARACTER*1069:* = 'U': Upper triangle of A is stored;070:* = 'L': Lower triangle of A is stored.071:*072:* N (input) INTEGER073:* The number of linear equations, i.e., the order of the074:* matrix A. N >= 0.075:*076:* NRHS (input) INTEGER077:* The number of right hand sides, i.e., the number of columns078:* of the matrices B and X. NRHS >= 0.079:*080:* AP (input) COMPLEX array, dimension (N*(N+1)/2)081:* The upper or lower triangle of the symmetric matrix A, packed082:* columnwise in a linear array. The j-th column of A is stored083:* in the array AP as follows:084:* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;085:* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.086:* See below for further details.087:*088:* AFP (input or output) COMPLEX array, dimension (N*(N+1)/2)089:* If FACT = 'F', then AFP is an input argument and on entry090:* contains the block diagonal matrix D and the multipliers used091:* to obtain the factor U or L from the factorization092:* A = U*D*U**T or A = L*D*L**T as computed by CSPTRF, stored as093:* a packed triangular matrix in the same storage format as A.094:*095:* If FACT = 'N', then AFP is an output argument and on exit096:* contains the block diagonal matrix D and the multipliers used097:* to obtain the factor U or L from the factorization098:* A = U*D*U**T or A = L*D*L**T as computed by CSPTRF, stored as099:* a packed triangular matrix in the same storage format as A.100:*101:* IPIV (input or output) INTEGER array, dimension (N)102:* If FACT = 'F', then IPIV is an input argument and on entry103:* contains details of the interchanges and the block structure104:* of D, as determined by CSPTRF.105:* If IPIV(k) > 0, then rows and columns k and IPIV(k) were106:* interchanged and D(k,k) is a 1-by-1 diagonal block.107:* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and108:* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)109:* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =110:* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were111:* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.112:*113:* If FACT = 'N', then IPIV is an output argument and on exit114:* contains details of the interchanges and the block structure115:* of D, as determined by CSPTRF.116:*117:* B (input) COMPLEX array, dimension (LDB,NRHS)118:* The N-by-NRHS right hand side matrix B.119:*120:* LDB (input) INTEGER121:* The leading dimension of the array B. LDB >= max(1,N).122:*123:* X (output) COMPLEX array, dimension (LDX,NRHS)124:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.125:*126:* LDX (input) INTEGER127:* The leading dimension of the array X. LDX >= max(1,N).128:*129:* RCOND (output) REAL130:* The estimate of the reciprocal condition number of the matrix131:* A. If RCOND is less than the machine precision (in132:* particular, if RCOND = 0), the matrix is singular to working133:* precision. This condition is indicated by a return code of134:* INFO > 0.135:*136:* FERR (output) REAL array, dimension (NRHS)137:* The estimated forward error bound for each solution vector138:* X(j) (the j-th column of the solution matrix X).139:* If XTRUE is the true solution corresponding to X(j), FERR(j)140:* is an estimated upper bound for the magnitude of the largest141:* element in (X(j) - XTRUE) divided by the magnitude of the142:* largest element in X(j). The estimate is as reliable as143:* the estimate for RCOND, and is almost always a slight144:* overestimate of the true error.145:*146:* BERR (output) REAL array, dimension (NRHS)147:* The componentwise relative backward error of each solution148:* vector X(j) (i.e., the smallest relative change in149:* any element of A or B that makes X(j) an exact solution).150:*151:* WORK (workspace) COMPLEX array, dimension (2*N)152:*153:* RWORK (workspace) REAL array, dimension (N)154:*155:* INFO (output) INTEGER156:* = 0: successful exit157:* < 0: if INFO = -i, the i-th argument had an illegal value158:* > 0: if INFO = i, and i is159:* <= N: D(i,i) is exactly zero. The factorization160:* has been completed but the factor D is exactly161:* singular, so the solution and error bounds could162:* not be computed. RCOND = 0 is returned.163:* = N+1: D is nonsingular, but RCOND is less than machine164:* precision, meaning that the matrix is singular165:* to working precision. Nevertheless, the166:* solution and error bounds are computed because167:* there are a number of situations where the168:* computed solution can be more accurate than the169:* value of RCOND would suggest.170:*171:* Further Details172:* ===============173:*174:* The packed storage scheme is illustrated by the following example175:* when N = 4, UPLO = 'U':176:*177:* Two-dimensional storage of the symmetric matrix A:178:*179:* a11 a12 a13 a14180:* a22 a23 a24181:* a33 a34 (aij = aji)182:* a44183:*184:* Packed storage of the upper triangle of A:185:*186:* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]187:*188:* =====================================================================189:*190:* .. Parameters ..191: REAL ZERO 192:PARAMETER( ZERO = 0.0E+0 ) 193:* ..194:* .. Local Scalars ..195:LOGICALNOFACT 196: REAL ANORM 197:* ..198:* .. External Functions ..199:LOGICALLSAME 200: REAL CLANSP, SLAMCH 201:EXTERNALLSAME, CLANSP, SLAMCH 202:* ..203:* .. External Subroutines ..204:EXTERNALCCOPY, CLACPY, CSPCON, CSPRFS, CSPTRF, CSPTRS, 205: $ XERBLA 206:* ..207:* .. Intrinsic Functions ..208:INTRINSICMAX 209:* ..210:* .. Executable Statements ..211:*212:* Test the input parameters.213:*214: INFO = 0 215: NOFACT =LSAME( FACT, 'N' ) 216:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN217: INFO = -1 218:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 219: $THEN220: INFO = -2 221:ELSEIF( N.LT.0 )THEN222: INFO = -3 223:ELSEIF( NRHS.LT.0 )THEN224: INFO = -4 225:ELSEIF( LDB.LT.MAX( 1, N ) )THEN226: INFO = -9 227:ELSEIF( LDX.LT.MAX( 1, N ) )THEN228: INFO = -11 229:ENDIF230:IF( INFO.NE.0 )THEN231:CALLXERBLA( 'CSPSVX', -INFO ) 232:RETURN233:ENDIF234:*235:IF( NOFACT )THEN236:*237:* Compute the factorization A = U*D*U' or A = L*D*L'.238:*239:CALLCCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 ) 240:CALLCSPTRF( UPLO, N, AFP, IPIV, INFO ) 241:*242:* Return if INFO is non-zero.243:*244:IF( INFO.GT.0 )THEN245: RCOND = ZERO 246:RETURN247:ENDIF248:ENDIF249:*250:* Compute the norm of the matrix A.251:*252: ANORM =CLANSP( 'I', UPLO, N, AP, RWORK ) 253:*254:* Compute the reciprocal of the condition number of A.255:*256:CALLCSPCON( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, INFO ) 257:*258:* Compute the solution vectors X.259:*260:CALLCLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 261:CALLCSPTRS( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO ) 262:*263:* Use iterative refinement to improve the computed solutions and264:* compute error bounds and backward error estimates for them.265:*266:CALLCSPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR, 267: $ BERR, WORK, RWORK, INFO ) 268:*269:* Set INFO = N+1 if the matrix is singular to working precision.270:*271:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 272: $ INFO = N + 1 273:*274:RETURN275:*276:* End of CSPSVX277:*278:END279: