001:SUBROUTINEDGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, 002: $ DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, 003: $ WORK, IWORK, INFO ) 004:*005:* -- LAPACK routine (version 3.2) --006:* -- LAPACK is a software package provided by Univ. of Tennessee, --007:* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--008:* November 2006009:*010:* .. Scalar Arguments ..011: CHARACTER FACT, TRANS 012: INTEGER INFO, LDB, LDX, N, NRHS 013: DOUBLE PRECISION RCOND 014:* ..015:* .. Array Arguments ..016: INTEGERIPIV( * ),IWORK( * ) 017: DOUBLE PRECISIONB( LDB, * ),BERR( * ),D( * ),DF( * ), 018: $DL( * ),DLF( * ),DU( * ),DU2( * ),DUF( * ), 019: $FERR( * ),WORK( * ),X( LDX, * ) 020:* ..021:*022:* Purpose023:* =======024:*025:* DGTSVX uses the LU factorization to compute the solution to a real026:* system of linear equations A * X = B or A**T * X = B,027:* where A is a tridiagonal matrix of order N and X and B are N-by-NRHS028:* matrices.029:*030:* Error bounds on the solution and a condition estimate are also031:* provided.032:*033:* Description034:* ===========035:*036:* The following steps are performed:037:*038:* 1. If FACT = 'N', the LU decomposition is used to factor the matrix A039:* as A = L * U, where L is a product of permutation and unit lower040:* bidiagonal matrices and U is upper triangular with nonzeros in041:* only the main diagonal and first two superdiagonals.042:*043:* 2. If some U(i,i)=0, so that U 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': DLF, DF, DUF, DU2, and IPIV contain the factored064:* form of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV065:* will not be modified.066:* = 'N': The matrix will be copied to DLF, DF, and DUF067:* and factored.068:*069:* TRANS (input) CHARACTER*1070:* Specifies the form of the system of equations:071:* = 'N': A * X = B (No transpose)072:* = 'T': A**T * X = B (Transpose)073:* = 'C': A**H * X = B (Conjugate transpose = Transpose)074:*075:* N (input) INTEGER076:* The order of the matrix A. N >= 0.077:*078:* NRHS (input) INTEGER079:* The number of right hand sides, i.e., the number of columns080:* of the matrix B. NRHS >= 0.081:*082:* DL (input) DOUBLE PRECISION array, dimension (N-1)083:* The (n-1) subdiagonal elements of A.084:*085:* D (input) DOUBLE PRECISION array, dimension (N)086:* The n diagonal elements of A.087:*088:* DU (input) DOUBLE PRECISION array, dimension (N-1)089:* The (n-1) superdiagonal elements of A.090:*091:* DLF (input or output) DOUBLE PRECISION array, dimension (N-1)092:* If FACT = 'F', then DLF is an input argument and on entry093:* contains the (n-1) multipliers that define the matrix L from094:* the LU factorization of A as computed by DGTTRF.095:*096:* If FACT = 'N', then DLF is an output argument and on exit097:* contains the (n-1) multipliers that define the matrix L from098:* the LU factorization of A.099:*100:* DF (input or output) DOUBLE PRECISION array, dimension (N)101:* If FACT = 'F', then DF is an input argument and on entry102:* contains the n diagonal elements of the upper triangular103:* matrix U from the LU factorization of A.104:*105:* If FACT = 'N', then DF is an output argument and on exit106:* contains the n diagonal elements of the upper triangular107:* matrix U from the LU factorization of A.108:*109:* DUF (input or output) DOUBLE PRECISION array, dimension (N-1)110:* If FACT = 'F', then DUF is an input argument and on entry111:* contains the (n-1) elements of the first superdiagonal of U.112:*113:* If FACT = 'N', then DUF is an output argument and on exit114:* contains the (n-1) elements of the first superdiagonal of U.115:*116:* DU2 (input or output) DOUBLE PRECISION array, dimension (N-2)117:* If FACT = 'F', then DU2 is an input argument and on entry118:* contains the (n-2) elements of the second superdiagonal of119:* U.120:*121:* If FACT = 'N', then DU2 is an output argument and on exit122:* contains the (n-2) elements of the second superdiagonal of123:* U.124:*125:* IPIV (input or output) INTEGER array, dimension (N)126:* If FACT = 'F', then IPIV is an input argument and on entry127:* contains the pivot indices from the LU factorization of A as128:* computed by DGTTRF.129:*130:* If FACT = 'N', then IPIV is an output argument and on exit131:* contains the pivot indices from the LU factorization of A;132:* row i of the matrix was interchanged with row IPIV(i).133:* IPIV(i) will always be either i or i+1; IPIV(i) = i indicates134:* a row interchange was not required.135:*136:* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)137:* The N-by-NRHS right hand side matrix B.138:*139:* LDB (input) INTEGER140:* The leading dimension of the array B. LDB >= max(1,N).141:*142:* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)143:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.144:*145:* LDX (input) INTEGER146:* The leading dimension of the array X. LDX >= max(1,N).147:*148:* RCOND (output) DOUBLE PRECISION149:* The estimate of the reciprocal condition number of the matrix150:* A. If RCOND is less than the machine precision (in151:* particular, if RCOND = 0), the matrix is singular to working152:* precision. This condition is indicated by a return code of153:* INFO > 0.154:*155:* FERR (output) DOUBLE PRECISION array, dimension (NRHS)156:* The estimated forward error bound for each solution vector157:* X(j) (the j-th column of the solution matrix X).158:* If XTRUE is the true solution corresponding to X(j), FERR(j)159:* is an estimated upper bound for the magnitude of the largest160:* element in (X(j) - XTRUE) divided by the magnitude of the161:* largest element in X(j). The estimate is as reliable as162:* the estimate for RCOND, and is almost always a slight163:* overestimate of the true error.164:*165:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)166:* The componentwise relative backward error of each solution167:* vector X(j) (i.e., the smallest relative change in168:* any element of A or B that makes X(j) an exact solution).169:*170:* WORK (workspace) DOUBLE PRECISION array, dimension (3*N)171:*172:* IWORK (workspace) INTEGER array, dimension (N)173:*174:* INFO (output) INTEGER175:* = 0: successful exit176:* < 0: if INFO = -i, the i-th argument had an illegal value177:* > 0: if INFO = i, and i is178:* <= N: U(i,i) is exactly zero. The factorization179:* has not been completed unless i = N, but the180:* factor U is exactly singular, so the solution181:* and error bounds could not be computed.182:* RCOND = 0 is returned.183:* = N+1: U is nonsingular, but RCOND is less than machine184:* precision, meaning that the matrix is singular185:* to working precision. Nevertheless, the186:* solution and error bounds are computed because187:* there are a number of situations where the188:* computed solution can be more accurate than the189:* value of RCOND would suggest.190:*191:* =====================================================================192:*193:* .. Parameters ..194: DOUBLE PRECISION ZERO 195:PARAMETER( ZERO = 0.0D+0 ) 196:* ..197:* .. Local Scalars ..198:LOGICALNOFACT, NOTRAN 199: CHARACTER NORM 200: DOUBLE PRECISION ANORM 201:* ..202:* .. External Functions ..203:LOGICALLSAME 204: DOUBLE PRECISION DLAMCH, DLANGT 205:EXTERNALLSAME, DLAMCH, DLANGT 206:* ..207:* .. External Subroutines ..208:EXTERNALDCOPY, DGTCON, DGTRFS, DGTTRF, DGTTRS, DLACPY, 209: $ XERBLA 210:* ..211:* .. Intrinsic Functions ..212:INTRINSICMAX 213:* ..214:* .. Executable Statements ..215:*216: INFO = 0 217: NOFACT =LSAME( FACT, 'N' ) 218: NOTRAN =LSAME( TRANS, 'N' ) 219:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN220: INFO = -1 221:ELSEIF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 222: $LSAME( TRANS, 'C' ) )THEN223: INFO = -2 224:ELSEIF( N.LT.0 )THEN225: INFO = -3 226:ELSEIF( NRHS.LT.0 )THEN227: INFO = -4 228:ELSEIF( LDB.LT.MAX( 1, N ) )THEN229: INFO = -14 230:ELSEIF( LDX.LT.MAX( 1, N ) )THEN231: INFO = -16 232:ENDIF233:IF( INFO.NE.0 )THEN234:CALLXERBLA( 'DGTSVX', -INFO ) 235:RETURN236:ENDIF237:*238:IF( NOFACT )THEN239:*240:* Compute the LU factorization of A.241:*242:CALLDCOPY( N, D, 1, DF, 1 ) 243:IF( N.GT.1 )THEN244:CALLDCOPY( N-1, DL, 1, DLF, 1 ) 245:CALLDCOPY( N-1, DU, 1, DUF, 1 ) 246:ENDIF247:CALLDGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO ) 248:*249:* Return if INFO is non-zero.250:*251:IF( INFO.GT.0 )THEN252: RCOND = ZERO 253:RETURN254:ENDIF255:ENDIF256:*257:* Compute the norm of the matrix A.258:*259:IF( NOTRAN )THEN260: NORM = '1' 261:ELSE262: NORM = 'I' 263:ENDIF264: ANORM =DLANGT( NORM, N, DL, D, DU ) 265:*266:* Compute the reciprocal of the condition number of A.267:*268:CALLDGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK, 269: $ IWORK, INFO ) 270:*271:* Compute the solution vectors X.272:*273:CALLDLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 274:CALLDGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX, 275: $ INFO ) 276:*277:* Use iterative refinement to improve the computed solutions and278:* compute error bounds and backward error estimates for them.279:*280:CALLDGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV, 281: $ B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO ) 282:*283:* Set INFO = N+1 if the matrix is singular to working precision.284:*285:IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 286: $ INFO = N + 1 287:*288:RETURN289:*290:* End of DGTSVX291:*292:END293: