001:SUBROUTINEZGTSVX( FACT, TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, 002: $ DU2, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, 003: $ WORK, RWORK, 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( * ) 017: DOUBLE PRECISIONBERR( * ),FERR( * ),RWORK( * ) 018: COMPLEX*16B( LDB, * ),D( * ),DF( * ),DL( * ), 019: $DLF( * ),DU( * ),DU2( * ),DUF( * ), 020: $WORK( * ),X( LDX, * ) 021:* ..022:*023:* Purpose024:* =======025:*026:* ZGTSVX uses the LU factorization to compute the solution to a complex027:* system of linear equations A * X = B, A**T * X = B, or A**H * X = B,028:* where A is a tridiagonal matrix of order N and X and B are N-by-NRHS029:* matrices.030:*031:* Error bounds on the solution and a condition estimate are also032:* provided.033:*034:* Description035:* ===========036:*037:* The following steps are performed:038:*039:* 1. If FACT = 'N', the LU decomposition is used to factor the matrix A040:* as A = L * U, where L is a product of permutation and unit lower041:* bidiagonal matrices and U is upper triangular with nonzeros in042:* only the main diagonal and first two superdiagonals.043:*044:* 2. If some U(i,i)=0, so that U is exactly singular, then the routine045:* returns with INFO = i. Otherwise, the factored form of A is used046:* to estimate the condition number of the matrix A. If the047:* reciprocal of the condition number is less than machine precision,048:* INFO = N+1 is returned as a warning, but the routine still goes on049:* to solve for X and compute error bounds as described below.050:*051:* 3. The system of equations is solved for X using the factored form052:* of A.053:*054:* 4. Iterative refinement is applied to improve the computed solution055:* matrix and calculate error bounds and backward error estimates056:* for it.057:*058:* Arguments059:* =========060:*061:* FACT (input) CHARACTER*1062:* Specifies whether or not the factored form of A has been063:* supplied on entry.064:* = 'F': DLF, DF, DUF, DU2, and IPIV contain the factored form065:* of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV will not066:* be modified.067:* = 'N': The matrix will be copied to DLF, DF, and DUF068:* and factored.069:*070:* TRANS (input) CHARACTER*1071:* Specifies the form of the system of equations:072:* = 'N': A * X = B (No transpose)073:* = 'T': A**T * X = B (Transpose)074:* = 'C': A**H * X = B (Conjugate transpose)075:*076:* N (input) INTEGER077:* The order of the matrix A. N >= 0.078:*079:* NRHS (input) INTEGER080:* The number of right hand sides, i.e., the number of columns081:* of the matrix B. NRHS >= 0.082:*083:* DL (input) COMPLEX*16 array, dimension (N-1)084:* The (n-1) subdiagonal elements of A.085:*086:* D (input) COMPLEX*16 array, dimension (N)087:* The n diagonal elements of A.088:*089:* DU (input) COMPLEX*16 array, dimension (N-1)090:* The (n-1) superdiagonal elements of A.091:*092:* DLF (input or output) COMPLEX*16 array, dimension (N-1)093:* If FACT = 'F', then DLF is an input argument and on entry094:* contains the (n-1) multipliers that define the matrix L from095:* the LU factorization of A as computed by ZGTTRF.096:*097:* If FACT = 'N', then DLF is an output argument and on exit098:* contains the (n-1) multipliers that define the matrix L from099:* the LU factorization of A.100:*101:* DF (input or output) COMPLEX*16 array, dimension (N)102:* If FACT = 'F', then DF is an input argument and on entry103:* contains the n diagonal elements of the upper triangular104:* matrix U from the LU factorization of A.105:*106:* If FACT = 'N', then DF is an output argument and on exit107:* contains the n diagonal elements of the upper triangular108:* matrix U from the LU factorization of A.109:*110:* DUF (input or output) COMPLEX*16 array, dimension (N-1)111:* If FACT = 'F', then DUF is an input argument and on entry112:* contains the (n-1) elements of the first superdiagonal of U.113:*114:* If FACT = 'N', then DUF is an output argument and on exit115:* contains the (n-1) elements of the first superdiagonal of U.116:*117:* DU2 (input or output) COMPLEX*16 array, dimension (N-2)118:* If FACT = 'F', then DU2 is an input argument and on entry119:* contains the (n-2) elements of the second superdiagonal of120:* U.121:*122:* If FACT = 'N', then DU2 is an output argument and on exit123:* contains the (n-2) elements of the second superdiagonal of124:* U.125:*126:* IPIV (input or output) INTEGER array, dimension (N)127:* If FACT = 'F', then IPIV is an input argument and on entry128:* contains the pivot indices from the LU factorization of A as129:* computed by ZGTTRF.130:*131:* If FACT = 'N', then IPIV is an output argument and on exit132:* contains the pivot indices from the LU factorization of A;133:* row i of the matrix was interchanged with row IPIV(i).134:* IPIV(i) will always be either i or i+1; IPIV(i) = i indicates135:* a row interchange was not required.136:*137:* B (input) COMPLEX*16 array, dimension (LDB,NRHS)138:* The N-by-NRHS right hand side matrix B.139:*140:* LDB (input) INTEGER141:* The leading dimension of the array B. LDB >= max(1,N).142:*143:* X (output) COMPLEX*16 array, dimension (LDX,NRHS)144:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.145:*146:* LDX (input) INTEGER147:* The leading dimension of the array X. LDX >= max(1,N).148:*149:* RCOND (output) DOUBLE PRECISION150:* The estimate of the reciprocal condition number of the matrix151:* A. If RCOND is less than the machine precision (in152:* particular, if RCOND = 0), the matrix is singular to working153:* precision. This condition is indicated by a return code of154:* INFO > 0.155:*156:* FERR (output) DOUBLE PRECISION array, dimension (NRHS)157:* The estimated forward error bound for each solution vector158:* X(j) (the j-th column of the solution matrix X).159:* If XTRUE is the true solution corresponding to X(j), FERR(j)160:* is an estimated upper bound for the magnitude of the largest161:* element in (X(j) - XTRUE) divided by the magnitude of the162:* largest element in X(j). The estimate is as reliable as163:* the estimate for RCOND, and is almost always a slight164:* overestimate of the true error.165:*166:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)167:* The componentwise relative backward error of each solution168:* vector X(j) (i.e., the smallest relative change in169:* any element of A or B that makes X(j) an exact solution).170:*171:* WORK (workspace) COMPLEX*16 array, dimension (2*N)172:*173:* RWORK (workspace) DOUBLE PRECISION array, dimension (N)174:*175:* INFO (output) INTEGER176:* = 0: successful exit177:* < 0: if INFO = -i, the i-th argument had an illegal value178:* > 0: if INFO = i, and i is179:* <= N: U(i,i) is exactly zero. The factorization180:* has not been completed unless i = N, but the181:* factor U is exactly singular, so the solution182:* and error bounds could not be computed.183:* RCOND = 0 is returned.184:* = N+1: U is nonsingular, but RCOND is less than machine185:* precision, meaning that the matrix is singular186:* to working precision. Nevertheless, the187:* solution and error bounds are computed because188:* there are a number of situations where the189:* computed solution can be more accurate than the190:* value of RCOND would suggest.191:*192:* =====================================================================193:*194:* .. Parameters ..195: DOUBLE PRECISION ZERO 196:PARAMETER( ZERO = 0.0D+0 ) 197:* ..198:* .. Local Scalars ..199:LOGICALNOFACT, NOTRAN 200: CHARACTER NORM 201: DOUBLE PRECISION ANORM 202:* ..203:* .. External Functions ..204:LOGICALLSAME 205: DOUBLE PRECISION DLAMCH, ZLANGT 206:EXTERNALLSAME, DLAMCH, ZLANGT 207:* ..208:* .. External Subroutines ..209:EXTERNALXERBLA, ZCOPY, ZGTCON, ZGTRFS, ZGTTRF, ZGTTRS, 210: $ ZLACPY 211:* ..212:* .. Intrinsic Functions ..213:INTRINSICMAX 214:* ..215:* .. Executable Statements ..216:*217: INFO = 0 218: NOFACT =LSAME( FACT, 'N' ) 219: NOTRAN =LSAME( TRANS, 'N' ) 220:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN221: INFO = -1 222:ELSEIF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. 223: $LSAME( TRANS, 'C' ) )THEN224: INFO = -2 225:ELSEIF( N.LT.0 )THEN226: INFO = -3 227:ELSEIF( NRHS.LT.0 )THEN228: INFO = -4 229:ELSEIF( LDB.LT.MAX( 1, N ) )THEN230: INFO = -14 231:ELSEIF( LDX.LT.MAX( 1, N ) )THEN232: INFO = -16 233:ENDIF234:IF( INFO.NE.0 )THEN235:CALLXERBLA( 'ZGTSVX', -INFO ) 236:RETURN237:ENDIF238:*239:IF( NOFACT )THEN240:*241:* Compute the LU factorization of A.242:*243:CALLZCOPY( N, D, 1, DF, 1 ) 244:IF( N.GT.1 )THEN245:CALLZCOPY( N-1, DL, 1, DLF, 1 ) 246:CALLZCOPY( N-1, DU, 1, DUF, 1 ) 247:ENDIF248:CALLZGTTRF( N, DLF, DF, DUF, DU2, IPIV, INFO ) 249:*250:* Return if INFO is non-zero.251:*252:IF( INFO.GT.0 )THEN253: RCOND = ZERO 254:RETURN255:ENDIF256:ENDIF257:*258:* Compute the norm of the matrix A.259:*260:IF( NOTRAN )THEN261: NORM = '1' 262:ELSE263: NORM = 'I' 264:ENDIF265: ANORM =ZLANGT( NORM, N, DL, D, DU ) 266:*267:* Compute the reciprocal of the condition number of A.268:*269:CALLZGTCON( NORM, N, DLF, DF, DUF, DU2, IPIV, ANORM, RCOND, WORK, 270: $ INFO ) 271:*272:* Compute the solution vectors X.273:*274:CALLZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 275:CALLZGTTRS( TRANS, N, NRHS, DLF, DF, DUF, DU2, IPIV, X, LDX, 276: $ INFO ) 277:*278:* Use iterative refinement to improve the computed solutions and279:* compute error bounds and backward error estimates for them.280:*281:CALLZGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2, IPIV, 282: $ B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO ) 283:*284:* Set INFO = N+1 if the matrix is singular to working precision.285:*286:IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 287: $ INFO = N + 1 288:*289:RETURN290:*291:* End of ZGTSVX292:*293:END294: