001:SUBROUTINESPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, 002: $ RCOND, FERR, BERR, WORK, 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: REALB( LDB, * ),BERR( * ),D( * ),DF( * ), 015: $E( * ),EF( * ),FERR( * ),WORK( * ), 016: $X( LDX, * ) 017:* ..018:*019:* Purpose020:* =======021:*022:* SPTSVX uses the factorization A = L*D*L**T to compute the solution023:* to a real system of linear equations A*X = B, where A is an N-by-N024:* symmetric positive definite tridiagonal matrix and X and B are025:* N-by-NRHS matrices.026:*027:* Error bounds on the solution and a condition estimate are also028:* provided.029:*030:* Description031:* ===========032:*033:* The following steps are performed:034:*035:* 1. If FACT = 'N', the matrix A is factored as A = L*D*L**T, where L036:* is a unit lower bidiagonal matrix and D is diagonal. The037:* factorization can also be regarded as having the form038:* A = U**T*D*U.039:*040:* 2. If the leading i-by-i principal minor is not positive definite,041:* then the routine returns with INFO = i. Otherwise, the factored042:* form of A is used to estimate the condition number of the matrix043:* A. If the reciprocal of the condition number is less than machine044:* precision, INFO = N+1 is returned as a warning, but the routine045:* still goes on to solve for X and compute error bounds as046:* described below.047:*048:* 3. The system of equations is solved for X using the factored form049:* of A.050:*051:* 4. Iterative refinement is applied to improve the computed solution052:* matrix and calculate error bounds and backward error estimates053:* for it.054:*055:* Arguments056:* =========057:*058:* FACT (input) CHARACTER*1059:* Specifies whether or not the factored form of A has been060:* supplied on entry.061:* = 'F': On entry, DF and EF contain the factored form of A.062:* D, E, DF, and EF will not be modified.063:* = 'N': The matrix A will be copied to DF and EF and064:* factored.065:*066:* N (input) INTEGER067:* The order of the matrix A. N >= 0.068:*069:* NRHS (input) INTEGER070:* The number of right hand sides, i.e., the number of columns071:* of the matrices B and X. NRHS >= 0.072:*073:* D (input) REAL array, dimension (N)074:* The n diagonal elements of the tridiagonal matrix A.075:*076:* E (input) REAL array, dimension (N-1)077:* The (n-1) subdiagonal elements of the tridiagonal matrix A.078:*079:* DF (input or output) REAL array, dimension (N)080:* If FACT = 'F', then DF is an input argument and on entry081:* contains the n diagonal elements of the diagonal matrix D082:* from the L*D*L**T factorization of A.083:* If FACT = 'N', then DF is an output argument and on exit084:* contains the n diagonal elements of the diagonal matrix D085:* from the L*D*L**T factorization of A.086:*087:* EF (input or output) REAL array, dimension (N-1)088:* If FACT = 'F', then EF is an input argument and on entry089:* contains the (n-1) subdiagonal elements of the unit090:* bidiagonal factor L from the L*D*L**T factorization of A.091:* If FACT = 'N', then EF is an output argument and on exit092:* contains the (n-1) subdiagonal elements of the unit093:* bidiagonal factor L from the L*D*L**T factorization of A.094:*095:* B (input) REAL array, dimension (LDB,NRHS)096:* The N-by-NRHS right hand side matrix B.097:*098:* LDB (input) INTEGER099:* The leading dimension of the array B. LDB >= max(1,N).100:*101:* X (output) REAL array, dimension (LDX,NRHS)102:* If INFO = 0 of INFO = N+1, the N-by-NRHS solution matrix X.103:*104:* LDX (input) INTEGER105:* The leading dimension of the array X. LDX >= max(1,N).106:*107:* RCOND (output) REAL108:* The reciprocal condition number of the matrix A. If RCOND109:* is less than the machine precision (in particular, if110:* RCOND = 0), the matrix is singular to working precision.111:* This condition is indicated by a return code of INFO > 0.112:*113:* FERR (output) REAL array, dimension (NRHS)114:* The forward error bound for each solution vector115:* X(j) (the j-th column of the solution matrix X).116:* If XTRUE is the true solution corresponding to X(j), FERR(j)117:* is an estimated upper bound for the magnitude of the largest118:* element in (X(j) - XTRUE) divided by the magnitude of the119:* largest element in X(j).120:*121:* BERR (output) REAL array, dimension (NRHS)122:* The componentwise relative backward error of each solution123:* vector X(j) (i.e., the smallest relative change in any124:* element of A or B that makes X(j) an exact solution).125:*126:* WORK (workspace) REAL array, dimension (2*N)127:*128:* INFO (output) INTEGER129:* = 0: successful exit130:* < 0: if INFO = -i, the i-th argument had an illegal value131:* > 0: if INFO = i, and i is132:* <= N: the leading minor of order i of A is133:* not positive definite, so the factorization134:* could not be completed, and the solution has not135:* been computed. RCOND = 0 is returned.136:* = N+1: U is nonsingular, but RCOND is less than machine137:* precision, meaning that the matrix is singular138:* to working precision. Nevertheless, the139:* solution and error bounds are computed because140:* there are a number of situations where the141:* computed solution can be more accurate than the142:* value of RCOND would suggest.143:*144:* =====================================================================145:*146:* .. Parameters ..147: REAL ZERO 148:PARAMETER( ZERO = 0.0E+0 ) 149:* ..150:* .. Local Scalars ..151:LOGICALNOFACT 152: REAL ANORM 153:* ..154:* .. External Functions ..155:LOGICALLSAME 156: REAL SLAMCH, SLANST 157:EXTERNALLSAME, SLAMCH, SLANST 158:* ..159:* .. External Subroutines ..160:EXTERNALSCOPY, SLACPY, SPTCON, SPTRFS, SPTTRF, SPTTRS, 161: $ XERBLA 162:* ..163:* .. Intrinsic Functions ..164:INTRINSICMAX 165:* ..166:* .. Executable Statements ..167:*168:* Test the input parameters.169:*170: INFO = 0 171: NOFACT =LSAME( FACT, 'N' ) 172:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN173: INFO = -1 174:ELSEIF( N.LT.0 )THEN175: INFO = -2 176:ELSEIF( NRHS.LT.0 )THEN177: INFO = -3 178:ELSEIF( LDB.LT.MAX( 1, N ) )THEN179: INFO = -9 180:ELSEIF( LDX.LT.MAX( 1, N ) )THEN181: INFO = -11 182:ENDIF183:IF( INFO.NE.0 )THEN184:CALLXERBLA( 'SPTSVX', -INFO ) 185:RETURN186:ENDIF187:*188:IF( NOFACT )THEN189:*190:* Compute the L*D*L' (or U'*D*U) factorization of A.191:*192:CALLSCOPY( N, D, 1, DF, 1 ) 193:IF( N.GT.1 ) 194: $CALLSCOPY( N-1, E, 1, EF, 1 ) 195:CALLSPTTRF( N, DF, EF, INFO ) 196:*197:* Return if INFO is non-zero.198:*199:IF( INFO.GT.0 )THEN200: RCOND = ZERO 201:RETURN202:ENDIF203:ENDIF204:*205:* Compute the norm of the matrix A.206:*207: ANORM =SLANST( '1', N, D, E ) 208:*209:* Compute the reciprocal of the condition number of A.210:*211:CALLSPTCON( N, DF, EF, ANORM, RCOND, WORK, INFO ) 212:*213:* Compute the solution vectors X.214:*215:CALLSLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 216:CALLSPTTRS( N, NRHS, DF, EF, X, LDX, INFO ) 217:*218:* Use iterative refinement to improve the computed solutions and219:* compute error bounds and backward error estimates for them.220:*221:CALLSPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR, 222: $ WORK, INFO ) 223:*224:* Set INFO = N+1 if the matrix is singular to working precision.225:*226:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 227: $ INFO = N + 1 228:*229:RETURN230:*231:* End of SPTSVX232:*233:END234: