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:* -- 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 011: INTEGER INFO, LDB, LDX, N, NRHS 012: REAL RCOND 013:* ..014:* .. Array Arguments ..015: REALBERR( * ),D( * ),DF( * ),FERR( * ), 016: $RWORK( * ) 017: COMPLEXB( LDB, * ),E( * ),EF( * ),WORK( * ), 018: $X( LDX, * ) 019:* ..020:*021:* Purpose022:* =======023:*024:* CPTSVX uses the factorization A = L*D*L**H to compute the solution025:* to a complex system of linear equations A*X = B, where A is an026:* N-by-N Hermitian positive definite tridiagonal matrix and X and B027:* 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 matrix A is factored as A = L*D*L**H, where L038:* is a unit lower bidiagonal matrix and D is diagonal. The039:* factorization can also be regarded as having the form040:* A = U**H*D*U.041:*042:* 2. If the leading i-by-i principal minor is not positive definite,043:* then the routine returns with INFO = i. Otherwise, the factored044:* form of A is used to estimate the condition number of the matrix045:* A. If the reciprocal of the condition number is less than machine046:* precision, INFO = N+1 is returned as a warning, but the routine047:* still goes on to solve for X and compute error bounds as048:* 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 the matrix062:* A is supplied on entry.063:* = 'F': On entry, DF and EF contain the factored form of A.064:* D, E, DF, and EF will not be modified.065:* = 'N': The matrix A will be copied to DF and EF and066:* factored.067:*068:* N (input) INTEGER069:* The order of the matrix A. N >= 0.070:*071:* NRHS (input) INTEGER072:* The number of right hand sides, i.e., the number of columns073:* of the matrices B and X. NRHS >= 0.074:*075:* D (input) REAL array, dimension (N)076:* The n diagonal elements of the tridiagonal matrix A.077:*078:* E (input) COMPLEX array, dimension (N-1)079:* The (n-1) subdiagonal elements of the tridiagonal matrix A.080:*081:* DF (input or output) REAL array, dimension (N)082:* If FACT = 'F', then DF is an input argument and on entry083:* contains the n diagonal elements of the diagonal matrix D084:* from the L*D*L**H factorization of A.085:* If FACT = 'N', then DF is an output argument and on exit086:* contains the n diagonal elements of the diagonal matrix D087:* from the L*D*L**H factorization of A.088:*089:* EF (input or output) COMPLEX array, dimension (N-1)090:* If FACT = 'F', then EF is an input argument and on entry091:* contains the (n-1) subdiagonal elements of the unit092:* bidiagonal factor L from the L*D*L**H factorization of A.093:* If FACT = 'N', then EF is an output argument and on exit094:* contains the (n-1) subdiagonal elements of the unit095:* bidiagonal factor L from the L*D*L**H factorization of A.096:*097:* B (input) COMPLEX array, dimension (LDB,NRHS)098:* The N-by-NRHS right hand side matrix B.099:*100:* LDB (input) INTEGER101:* The leading dimension of the array B. LDB >= max(1,N).102:*103:* X (output) COMPLEX array, dimension (LDX,NRHS)104:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.105:*106:* LDX (input) INTEGER107:* The leading dimension of the array X. LDX >= max(1,N).108:*109:* RCOND (output) REAL110:* The reciprocal condition number of the matrix A. If RCOND111:* is less than the machine precision (in particular, if112:* RCOND = 0), the matrix is singular to working precision.113:* This condition is indicated by a return code of INFO > 0.114:*115:* FERR (output) REAL array, dimension (NRHS)116:* The forward error bound for each solution vector117:* X(j) (the j-th column of the solution matrix X).118:* If XTRUE is the true solution corresponding to X(j), FERR(j)119:* is an estimated upper bound for the magnitude of the largest120:* element in (X(j) - XTRUE) divided by the magnitude of the121:* largest element in X(j).122:*123:* BERR (output) REAL array, dimension (NRHS)124:* The componentwise relative backward error of each solution125:* vector X(j) (i.e., the smallest relative change in any126:* element of A or B that makes X(j) an exact solution).127:*128:* WORK (workspace) COMPLEX array, dimension (N)129:*130:* RWORK (workspace) REAL array, dimension (N)131:*132:* INFO (output) INTEGER133:* = 0: successful exit134:* < 0: if INFO = -i, the i-th argument had an illegal value135:* > 0: if INFO = i, and i is136:* <= N: the leading minor of order i of A is137:* not positive definite, so the factorization138:* could not be completed, and the solution has not139:* been computed. RCOND = 0 is returned.140:* = N+1: U is nonsingular, but RCOND is less than machine141:* precision, meaning that the matrix is singular142:* to working precision. Nevertheless, the143:* solution and error bounds are computed because144:* there are a number of situations where the145:* computed solution can be more accurate than the146:* value of RCOND would suggest.147:*148:* =====================================================================149:*150:* .. Parameters ..151: REAL ZERO 152:PARAMETER( ZERO = 0.0E+0 ) 153:* ..154:* .. Local Scalars ..155:LOGICALNOFACT 156: REAL ANORM 157:* ..158:* .. External Functions ..159:LOGICALLSAME 160: REAL CLANHT, SLAMCH 161:EXTERNALLSAME, CLANHT, SLAMCH 162:* ..163:* .. External Subroutines ..164:EXTERNALCCOPY, CLACPY, CPTCON, CPTRFS, CPTTRF, CPTTRS, 165: $ SCOPY, XERBLA 166:* ..167:* .. Intrinsic Functions ..168:INTRINSICMAX 169:* ..170:* .. Executable Statements ..171:*172:* Test the input parameters.173:*174: INFO = 0 175: NOFACT =LSAME( FACT, 'N' ) 176:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN177: INFO = -1 178:ELSEIF( N.LT.0 )THEN179: INFO = -2 180:ELSEIF( NRHS.LT.0 )THEN181: INFO = -3 182:ELSEIF( LDB.LT.MAX( 1, N ) )THEN183: INFO = -9 184:ELSEIF( LDX.LT.MAX( 1, N ) )THEN185: INFO = -11 186:ENDIF187:IF( INFO.NE.0 )THEN188:CALLXERBLA( 'CPTSVX', -INFO ) 189:RETURN190:ENDIF191:*192:IF( NOFACT )THEN193:*194:* Compute the L*D*L' (or U'*D*U) factorization of A.195:*196:CALLSCOPY( N, D, 1, DF, 1 ) 197:IF( N.GT.1 ) 198: $CALLCCOPY( N-1, E, 1, EF, 1 ) 199:CALLCPTTRF( N, DF, EF, INFO ) 200:*201:* Return if INFO is non-zero.202:*203:IF( INFO.GT.0 )THEN204: RCOND = ZERO 205:RETURN206:ENDIF207:ENDIF208:*209:* Compute the norm of the matrix A.210:*211: ANORM =CLANHT( '1', N, D, E ) 212:*213:* Compute the reciprocal of the condition number of A.214:*215:CALLCPTCON( N, DF, EF, ANORM, RCOND, RWORK, INFO ) 216:*217:* Compute the solution vectors X.218:*219:CALLCLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 220:CALLCPTTRS( 'Lower', N, NRHS, DF, EF, X, LDX, INFO ) 221:*222:* Use iterative refinement to improve the computed solutions and223:* compute error bounds and backward error estimates for them.224:*225:CALLCPTRFS( 'Lower', N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, 226: $ BERR, WORK, RWORK, INFO ) 227:*228:* Set INFO = N+1 if the matrix is singular to working precision.229:*230:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 231: $ INFO = N + 1 232:*233:RETURN234:*235:* End of CPTSVX236:*237:END238: