001:SUBROUTINEZSPSVX( 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:* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..006:* November 2006007:*008:* .. Scalar Arguments ..009: CHARACTER FACT, UPLO 010: INTEGER INFO, LDB, LDX, N, NRHS 011: DOUBLE PRECISION RCOND 012:* ..013:* .. Array Arguments ..014: INTEGERIPIV( * ) 015: DOUBLE PRECISIONBERR( * ),FERR( * ),RWORK( * ) 016: COMPLEX*16AFP( * ),AP( * ),B( LDB, * ),WORK( * ), 017: $X( LDX, * ) 018:* ..019:*020:* Purpose021:* =======022:*023:* ZSPSVX uses the diagonal pivoting factorization A = U*D*U**T or024:* A = L*D*L**T to compute the solution to a complex system of linear025:* equations A * X = B, where A is an N-by-N symmetric matrix stored026:* in 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 = 'N', the diagonal pivoting method is used to factor A as037:* A = U * D * U**T, if UPLO = 'U', or038:* A = L * D * L**T, if UPLO = 'L',039:* where U (or L) is a product of permutation and unit upper (lower)040:* triangular matrices and D is symmetric and block diagonal with041:* 1-by-1 and 2-by-2 diagonal blocks.042:*043:* 2. If some D(i,i)=0, so that D is exactly singular, then the routine044:* returns with INFO = i. Otherwise, the factored form of A is used045:* to estimate the condition number of the matrix A. If the046:* reciprocal of the condition number is less than machine precision,047:* INFO = N+1 is returned as a warning, but the routine still goes on048:* to solve for X and compute error bounds as described below.049:*050:* 3. The system of equations is solved for X using the factored form051:* of A.052:*053:* 4. Iterative refinement is applied to improve the computed solution054:* matrix and calculate error bounds and backward error estimates055:* for it.056:*057:* Arguments058:* =========059:*060:* FACT (input) CHARACTER*1061:* Specifies whether or not the factored form of A has been062:* supplied on entry.063:* = 'F': On entry, AFP and IPIV contain the factored form064:* of A. AP, AFP and IPIV will not be modified.065:* = 'N': The matrix A will be copied to AFP and factored.066:*067:* UPLO (input) CHARACTER*1068:* = 'U': Upper triangle of A is stored;069:* = 'L': Lower triangle of A is stored.070:*071:* N (input) INTEGER072:* The number of linear equations, i.e., the order of the073:* matrix A. N >= 0.074:*075:* NRHS (input) INTEGER076:* The number of right hand sides, i.e., the number of columns077:* of the matrices B and X. NRHS >= 0.078:*079:* AP (input) COMPLEX*16 array, dimension (N*(N+1)/2)080:* The upper or lower triangle of the symmetric matrix A, packed081:* columnwise in a linear array. The j-th column of A is stored082:* in the array AP as follows:083:* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;084:* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.085:* See below for further details.086:*087:* AFP (input or output) COMPLEX*16 array, dimension (N*(N+1)/2)088:* If FACT = 'F', then AFP is an input argument and on entry089:* contains the block diagonal matrix D and the multipliers used090:* to obtain the factor U or L from the factorization091:* A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as092:* a packed triangular matrix in the same storage format as A.093:*094:* If FACT = 'N', then AFP is an output argument and on exit095:* contains the block diagonal matrix D and the multipliers used096:* to obtain the factor U or L from the factorization097:* A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as098:* a packed triangular matrix in the same storage format as A.099:*100:* IPIV (input or output) INTEGER array, dimension (N)101:* If FACT = 'F', then IPIV is an input argument and on entry102:* contains details of the interchanges and the block structure103:* of D, as determined by ZSPTRF.104:* If IPIV(k) > 0, then rows and columns k and IPIV(k) were105:* interchanged and D(k,k) is a 1-by-1 diagonal block.106:* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and107:* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)108:* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =109:* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were110:* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.111:*112:* If FACT = 'N', then IPIV is an output argument and on exit113:* contains details of the interchanges and the block structure114:* of D, as determined by ZSPTRF.115:*116:* B (input) COMPLEX*16 array, dimension (LDB,NRHS)117:* The N-by-NRHS right hand side matrix B.118:*119:* LDB (input) INTEGER120:* The leading dimension of the array B. LDB >= max(1,N).121:*122:* X (output) COMPLEX*16 array, dimension (LDX,NRHS)123:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.124:*125:* LDX (input) INTEGER126:* The leading dimension of the array X. LDX >= max(1,N).127:*128:* RCOND (output) DOUBLE PRECISION129:* The estimate of the reciprocal condition number of the matrix130:* A. If RCOND is less than the machine precision (in131:* particular, if RCOND = 0), the matrix is singular to working132:* precision. This condition is indicated by a return code of133:* INFO > 0.134:*135:* FERR (output) DOUBLE PRECISION array, dimension (NRHS)136:* The estimated forward error bound for each solution vector137:* X(j) (the j-th column of the solution matrix X).138:* If XTRUE is the true solution corresponding to X(j), FERR(j)139:* is an estimated upper bound for the magnitude of the largest140:* element in (X(j) - XTRUE) divided by the magnitude of the141:* largest element in X(j). The estimate is as reliable as142:* the estimate for RCOND, and is almost always a slight143:* overestimate of the true error.144:*145:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)146:* The componentwise relative backward error of each solution147:* vector X(j) (i.e., the smallest relative change in148:* any element of A or B that makes X(j) an exact solution).149:*150:* WORK (workspace) COMPLEX*16 array, dimension (2*N)151:*152:* RWORK (workspace) DOUBLE PRECISION array, dimension (N)153:*154:* INFO (output) INTEGER155:* = 0: successful exit156:* < 0: if INFO = -i, the i-th argument had an illegal value157:* > 0: if INFO = i, and i is158:* <= N: D(i,i) is exactly zero. The factorization159:* has been completed but the factor D is exactly160:* singular, so the solution and error bounds could161:* not be computed. RCOND = 0 is returned.162:* = N+1: D is nonsingular, but RCOND is less than machine163:* precision, meaning that the matrix is singular164:* to working precision. Nevertheless, the165:* solution and error bounds are computed because166:* there are a number of situations where the167:* computed solution can be more accurate than the168:* value of RCOND would suggest.169:*170:* Further Details171:* ===============172:*173:* The packed storage scheme is illustrated by the following example174:* when N = 4, UPLO = 'U':175:*176:* Two-dimensional storage of the symmetric matrix A:177:*178:* a11 a12 a13 a14179:* a22 a23 a24180:* a33 a34 (aij = aji)181:* a44182:*183:* Packed storage of the upper triangle of A:184:*185:* AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]186:*187:* =====================================================================188:*189:* .. Parameters ..190: DOUBLE PRECISION ZERO 191:PARAMETER( ZERO = 0.0D+0 ) 192:* ..193:* .. Local Scalars ..194:LOGICALNOFACT 195: DOUBLE PRECISION ANORM 196:* ..197:* .. External Functions ..198:LOGICALLSAME 199: DOUBLE PRECISION DLAMCH, ZLANSP 200:EXTERNALLSAME, DLAMCH, ZLANSP 201:* ..202:* .. External Subroutines ..203:EXTERNALXERBLA, ZCOPY, ZLACPY, ZSPCON, ZSPRFS, ZSPTRF, 204: $ ZSPTRS 205:* ..206:* .. Intrinsic Functions ..207:INTRINSICMAX 208:* ..209:* .. Executable Statements ..210:*211:* Test the input parameters.212:*213: INFO = 0 214: NOFACT =LSAME( FACT, 'N' ) 215:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN216: INFO = -1 217:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 218: $THEN219: INFO = -2 220:ELSEIF( N.LT.0 )THEN221: INFO = -3 222:ELSEIF( NRHS.LT.0 )THEN223: INFO = -4 224:ELSEIF( LDB.LT.MAX( 1, N ) )THEN225: INFO = -9 226:ELSEIF( LDX.LT.MAX( 1, N ) )THEN227: INFO = -11 228:ENDIF229:IF( INFO.NE.0 )THEN230:CALLXERBLA( 'ZSPSVX', -INFO ) 231:RETURN232:ENDIF233:*234:IF( NOFACT )THEN235:*236:* Compute the factorization A = U*D*U' or A = L*D*L'.237:*238:CALLZCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 ) 239:CALLZSPTRF( UPLO, N, AFP, IPIV, INFO ) 240:*241:* Return if INFO is non-zero.242:*243:IF( INFO.GT.0 )THEN244: RCOND = ZERO 245:RETURN246:ENDIF247:ENDIF248:*249:* Compute the norm of the matrix A.250:*251: ANORM =ZLANSP( 'I', UPLO, N, AP, RWORK ) 252:*253:* Compute the reciprocal of the condition number of A.254:*255:CALLZSPCON( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, INFO ) 256:*257:* Compute the solution vectors X.258:*259:CALLZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 260:CALLZSPTRS( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO ) 261:*262:* Use iterative refinement to improve the computed solutions and263:* compute error bounds and backward error estimates for them.264:*265:CALLZSPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR, 266: $ BERR, WORK, RWORK, INFO ) 267:*268:* Set INFO = N+1 if the matrix is singular to working precision.269:*270:IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 271: $ INFO = N + 1 272:*273:RETURN274:*275:* End of ZSPSVX276:*277:END278: