001:SUBROUTINECPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, 002: $ RCOND, FERR, BERR, WORK, RWORK, INFO ) 003:*004:* -- LAPACK routine (version 3.2) --005:* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..006:* November 2006007:*008:* .. Scalar Arguments ..009: CHARACTER FACT 010: INTEGER INFO, LDB, LDX, N, NRHS 011: REAL RCOND 012:* ..013:* .. Array Arguments ..014: REALBERR( * ),D( * ),DF( * ),FERR( * ), 015: $RWORK( * ) 016: COMPLEXB( LDB, * ),E( * ),EF( * ),WORK( * ), 017: $X( LDX, * ) 018:* ..019:*020:* Purpose021:* =======022:*023:* CPTSVX uses the factorization A = L*D*L**H to compute the solution024:* to a complex system of linear equations A*X = B, where A is an025:* N-by-N Hermitian positive definite tridiagonal matrix and X and B026:* 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 matrix A is factored as A = L*D*L**H, where L037:* is a unit lower bidiagonal matrix and D is diagonal. The038:* factorization can also be regarded as having the form039:* A = U**H*D*U.040:*041:* 2. If the leading i-by-i principal minor is not positive definite,042:* then the routine returns with INFO = i. Otherwise, the factored043:* form of A is used to estimate the condition number of the matrix044:* A. If the reciprocal of the condition number is less than machine045:* precision, INFO = N+1 is returned as a warning, but the routine046:* still goes on to solve for X and compute error bounds as047:* described below.048:*049:* 3. The system of equations is solved for X using the factored form050:* of A.051:*052:* 4. Iterative refinement is applied to improve the computed solution053:* matrix and calculate error bounds and backward error estimates054:* for it.055:*056:* Arguments057:* =========058:*059:* FACT (input) CHARACTER*1060:* Specifies whether or not the factored form of the matrix061:* A is supplied on entry.062:* = 'F': On entry, DF and EF contain the factored form of A.063:* D, E, DF, and EF will not be modified.064:* = 'N': The matrix A will be copied to DF and EF and065:* factored.066:*067:* N (input) INTEGER068:* The order of the matrix A. N >= 0.069:*070:* NRHS (input) INTEGER071:* The number of right hand sides, i.e., the number of columns072:* of the matrices B and X. NRHS >= 0.073:*074:* D (input) REAL array, dimension (N)075:* The n diagonal elements of the tridiagonal matrix A.076:*077:* E (input) COMPLEX array, dimension (N-1)078:* The (n-1) subdiagonal elements of the tridiagonal matrix A.079:*080:* DF (input or output) REAL array, dimension (N)081:* If FACT = 'F', then DF is an input argument and on entry082:* contains the n diagonal elements of the diagonal matrix D083:* from the L*D*L**H factorization of A.084:* If FACT = 'N', then DF is an output argument and on exit085:* contains the n diagonal elements of the diagonal matrix D086:* from the L*D*L**H factorization of A.087:*088:* EF (input or output) COMPLEX array, dimension (N-1)089:* If FACT = 'F', then EF is an input argument and on entry090:* contains the (n-1) subdiagonal elements of the unit091:* bidiagonal factor L from the L*D*L**H factorization of A.092:* If FACT = 'N', then EF is an output argument and on exit093:* contains the (n-1) subdiagonal elements of the unit094:* bidiagonal factor L from the L*D*L**H factorization of A.095:*096:* B (input) COMPLEX array, dimension (LDB,NRHS)097:* The N-by-NRHS right hand side matrix B.098:*099:* LDB (input) INTEGER100:* The leading dimension of the array B. LDB >= max(1,N).101:*102:* X (output) COMPLEX array, dimension (LDX,NRHS)103:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.104:*105:* LDX (input) INTEGER106:* The leading dimension of the array X. LDX >= max(1,N).107:*108:* RCOND (output) REAL109:* The reciprocal condition number of the matrix A. If RCOND110:* is less than the machine precision (in particular, if111:* RCOND = 0), the matrix is singular to working precision.112:* This condition is indicated by a return code of INFO > 0.113:*114:* FERR (output) REAL array, dimension (NRHS)115:* The forward error bound for each solution vector116:* X(j) (the j-th column of the solution matrix X).117:* If XTRUE is the true solution corresponding to X(j), FERR(j)118:* is an estimated upper bound for the magnitude of the largest119:* element in (X(j) - XTRUE) divided by the magnitude of the120:* largest element in X(j).121:*122:* BERR (output) REAL array, dimension (NRHS)123:* The componentwise relative backward error of each solution124:* vector X(j) (i.e., the smallest relative change in any125:* element of A or B that makes X(j) an exact solution).126:*127:* WORK (workspace) COMPLEX array, dimension (N)128:*129:* RWORK (workspace) REAL array, dimension (N)130:*131:* INFO (output) INTEGER132:* = 0: successful exit133:* < 0: if INFO = -i, the i-th argument had an illegal value134:* > 0: if INFO = i, and i is135:* <= N: the leading minor of order i of A is136:* not positive definite, so the factorization137:* could not be completed, and the solution has not138:* been computed. RCOND = 0 is returned.139:* = N+1: U is nonsingular, but RCOND is less than machine140:* precision, meaning that the matrix is singular141:* to working precision. Nevertheless, the142:* solution and error bounds are computed because143:* there are a number of situations where the144:* computed solution can be more accurate than the145:* value of RCOND would suggest.146:*147:* =====================================================================148:*149:* .. Parameters ..150: REAL ZERO 151:PARAMETER( ZERO = 0.0E+0 ) 152:* ..153:* .. Local Scalars ..154:LOGICALNOFACT 155: REAL ANORM 156:* ..157:* .. External Functions ..158:LOGICALLSAME 159: REAL CLANHT, SLAMCH 160:EXTERNALLSAME, CLANHT, SLAMCH 161:* ..162:* .. External Subroutines ..163:EXTERNALCCOPY, CLACPY, CPTCON, CPTRFS, CPTTRF, CPTTRS, 164: $ SCOPY, XERBLA 165:* ..166:* .. Intrinsic Functions ..167:INTRINSICMAX 168:* ..169:* .. Executable Statements ..170:*171:* Test the input parameters.172:*173: INFO = 0 174: NOFACT =LSAME( FACT, 'N' ) 175:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN176: INFO = -1 177:ELSEIF( N.LT.0 )THEN178: INFO = -2 179:ELSEIF( NRHS.LT.0 )THEN180: INFO = -3 181:ELSEIF( LDB.LT.MAX( 1, N ) )THEN182: INFO = -9 183:ELSEIF( LDX.LT.MAX( 1, N ) )THEN184: INFO = -11 185:ENDIF186:IF( INFO.NE.0 )THEN187:CALLXERBLA( 'CPTSVX', -INFO ) 188:RETURN189:ENDIF190:*191:IF( NOFACT )THEN192:*193:* Compute the L*D*L' (or U'*D*U) factorization of A.194:*195:CALLSCOPY( N, D, 1, DF, 1 ) 196:IF( N.GT.1 ) 197: $CALLCCOPY( N-1, E, 1, EF, 1 ) 198:CALLCPTTRF( N, DF, EF, INFO ) 199:*200:* Return if INFO is non-zero.201:*202:IF( INFO.GT.0 )THEN203: RCOND = ZERO 204:RETURN205:ENDIF206:ENDIF207:*208:* Compute the norm of the matrix A.209:*210: ANORM =CLANHT( '1', N, D, E ) 211:*212:* Compute the reciprocal of the condition number of A.213:*214:CALLCPTCON( N, DF, EF, ANORM, RCOND, RWORK, INFO ) 215:*216:* Compute the solution vectors X.217:*218:CALLCLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 219:CALLCPTTRS( 'Lower', N, NRHS, DF, EF, X, LDX, INFO ) 220:*221:* Use iterative refinement to improve the computed solutions and222:* compute error bounds and backward error estimates for them.223:*224:CALLCPTRFS( 'Lower', N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, 225: $ BERR, WORK, RWORK, INFO ) 226:*227:* Set INFO = N+1 if the matrix is singular to working precision.228:*229:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 230: $ INFO = N + 1 231:*232:RETURN233:*234:* End of CPTSVX235:*236:END237: