001:SUBROUTINECHESVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, 002: $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, 003: $ RWORK, INFO ) 004:*005:* -- LAPACK driver 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, UPLO 012: INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 013: REAL RCOND 014:* ..015:* .. Array Arguments ..016: INTEGERIPIV( * ) 017: REALBERR( * ),FERR( * ),RWORK( * ) 018: COMPLEXA( LDA, * ),AF( LDAF, * ),B( LDB, * ), 019: $WORK( * ),X( LDX, * ) 020:* ..021:*022:* Purpose023:* =======024:*025:* CHESVX uses the diagonal pivoting factorization to compute the026:* solution to a complex system of linear equations A * X = B,027:* where A is an N-by-N Hermitian matrix 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 diagonal pivoting method is used to factor A.039:* The form of the factorization is040:* A = U * D * U**H, if UPLO = 'U', or041:* A = L * D * L**H, if UPLO = 'L',042:* where U (or L) is a product of permutation and unit upper (lower)043:* triangular matrices, and D is Hermitian and block diagonal with044:* 1-by-1 and 2-by-2 diagonal blocks.045:*046:* 2. If some D(i,i)=0, so that D is exactly singular, then the routine047:* returns with INFO = i. Otherwise, the factored form of A is used048:* to estimate the condition number of the matrix A. If the049:* reciprocal of the condition number is less than machine precision,050:* INFO = N+1 is returned as a warning, but the routine still goes on051:* to solve for X and compute error bounds as described below.052:*053:* 3. The system of equations is solved for X using the factored form054:* of A.055:*056:* 4. Iterative refinement is applied to improve the computed solution057:* matrix and calculate error bounds and backward error estimates058:* for it.059:*060:* Arguments061:* =========062:*063:* FACT (input) CHARACTER*1064:* Specifies whether or not the factored form of A has been065:* supplied on entry.066:* = 'F': On entry, AF and IPIV contain the factored form067:* of A. A, AF and IPIV will not be modified.068:* = 'N': The matrix A will be copied to AF and factored.069:*070:* UPLO (input) CHARACTER*1071:* = 'U': Upper triangle of A is stored;072:* = 'L': Lower triangle of A is stored.073:*074:* N (input) INTEGER075:* The number of linear equations, i.e., the order of the076:* matrix A. N >= 0.077:*078:* NRHS (input) INTEGER079:* The number of right hand sides, i.e., the number of columns080:* of the matrices B and X. NRHS >= 0.081:*082:* A (input) COMPLEX array, dimension (LDA,N)083:* The Hermitian matrix A. If UPLO = 'U', the leading N-by-N084:* upper triangular part of A contains the upper triangular part085:* of the matrix A, and the strictly lower triangular part of A086:* is not referenced. If UPLO = 'L', the leading N-by-N lower087:* triangular part of A contains the lower triangular part of088:* the matrix A, and the strictly upper triangular part of A is089:* not referenced.090:*091:* LDA (input) INTEGER092:* The leading dimension of the array A. LDA >= max(1,N).093:*094:* AF (input or output) COMPLEX array, dimension (LDAF,N)095:* If FACT = 'F', then AF is an input argument and on entry096:* contains the block diagonal matrix D and the multipliers used097:* to obtain the factor U or L from the factorization098:* A = U*D*U**H or A = L*D*L**H as computed by CHETRF.099:*100:* If FACT = 'N', then AF is an output argument and on exit101:* returns the block diagonal matrix D and the multipliers used102:* to obtain the factor U or L from the factorization103:* A = U*D*U**H or A = L*D*L**H.104:*105:* LDAF (input) INTEGER106:* The leading dimension of the array AF. LDAF >= max(1,N).107:*108:* IPIV (input or output) INTEGER array, dimension (N)109:* If FACT = 'F', then IPIV is an input argument and on entry110:* contains details of the interchanges and the block structure111:* of D, as determined by CHETRF.112:* If IPIV(k) > 0, then rows and columns k and IPIV(k) were113:* interchanged and D(k,k) is a 1-by-1 diagonal block.114:* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and115:* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)116:* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =117:* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were118:* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.119:*120:* If FACT = 'N', then IPIV is an output argument and on exit121:* contains details of the interchanges and the block structure122:* of D, as determined by CHETRF.123:*124:* B (input) COMPLEX array, dimension (LDB,NRHS)125:* The N-by-NRHS right hand side matrix B.126:*127:* LDB (input) INTEGER128:* The leading dimension of the array B. LDB >= max(1,N).129:*130:* X (output) COMPLEX array, dimension (LDX,NRHS)131:* If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.132:*133:* LDX (input) INTEGER134:* The leading dimension of the array X. LDX >= max(1,N).135:*136:* RCOND (output) REAL137:* The estimate of the reciprocal condition number of the matrix138:* A. If RCOND is less than the machine precision (in139:* particular, if RCOND = 0), the matrix is singular to working140:* precision. This condition is indicated by a return code of141:* INFO > 0.142:*143:* FERR (output) REAL array, dimension (NRHS)144:* The estimated forward error bound for each solution vector145:* X(j) (the j-th column of the solution matrix X).146:* If XTRUE is the true solution corresponding to X(j), FERR(j)147:* is an estimated upper bound for the magnitude of the largest148:* element in (X(j) - XTRUE) divided by the magnitude of the149:* largest element in X(j). The estimate is as reliable as150:* the estimate for RCOND, and is almost always a slight151:* overestimate of the true error.152:*153:* BERR (output) REAL array, dimension (NRHS)154:* The componentwise relative backward error of each solution155:* vector X(j) (i.e., the smallest relative change in156:* any element of A or B that makes X(j) an exact solution).157:*158:* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))159:* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.160:*161:* LWORK (input) INTEGER162:* The length of WORK. LWORK >= max(1,2*N), and for best163:* performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where164:* NB is the optimal blocksize for CHETRF.165:*166:* If LWORK = -1, then a workspace query is assumed; the routine167:* only calculates the optimal size of the WORK array, returns168:* this value as the first entry of the WORK array, and no error169:* message related to LWORK is issued by XERBLA.170:*171:* RWORK (workspace) REAL array, dimension (N)172:*173:* INFO (output) INTEGER174:* = 0: successful exit175:* < 0: if INFO = -i, the i-th argument had an illegal value176:* > 0: if INFO = i, and i is177:* <= N: D(i,i) is exactly zero. The factorization178:* has been completed but the factor D is exactly179:* singular, so the solution and error bounds could180:* not be computed. RCOND = 0 is returned.181:* = N+1: D is nonsingular, but RCOND is less than machine182:* precision, meaning that the matrix is singular183:* to working precision. Nevertheless, the184:* solution and error bounds are computed because185:* there are a number of situations where the186:* computed solution can be more accurate than the187:* value of RCOND would suggest.188:*189:* =====================================================================190:*191:* .. Parameters ..192: REAL ZERO 193:PARAMETER( ZERO = 0.0E+0 ) 194:* ..195:* .. Local Scalars ..196:LOGICALLQUERY, NOFACT 197: INTEGER LWKOPT, NB 198: REAL ANORM 199:* ..200:* .. External Functions ..201:LOGICALLSAME 202: INTEGER ILAENV 203: REAL CLANHE, SLAMCH 204:EXTERNALILAENV, LSAME, CLANHE, SLAMCH 205:* ..206:* .. External Subroutines ..207:EXTERNALCHECON, CHERFS, CHETRF, CHETRS, CLACPY, XERBLA 208:* ..209:* .. Intrinsic Functions ..210:INTRINSICMAX 211:* ..212:* .. Executable Statements ..213:*214:* Test the input parameters.215:*216: INFO = 0 217: NOFACT =LSAME( FACT, 'N' ) 218: LQUERY = ( LWORK.EQ.-1 ) 219:IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) )THEN220: INFO = -1 221:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 222: $THEN223: INFO = -2 224:ELSEIF( N.LT.0 )THEN225: INFO = -3 226:ELSEIF( NRHS.LT.0 )THEN227: INFO = -4 228:ELSEIF( LDA.LT.MAX( 1, N ) )THEN229: INFO = -6 230:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN231: INFO = -8 232:ELSEIF( LDB.LT.MAX( 1, N ) )THEN233: INFO = -11 234:ELSEIF( LDX.LT.MAX( 1, N ) )THEN235: INFO = -13 236:ELSEIF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY )THEN237: INFO = -18 238:ENDIF239:*240:IF( INFO.EQ.0 )THEN241: LWKOPT =MAX( 1, 2*N ) 242:IF( NOFACT )THEN243: NB =ILAENV( 1, 'CHETRF', UPLO, N, -1, -1, -1 ) 244: LWKOPT =MAX( LWKOPT, N*NB ) 245:ENDIF246:WORK( 1 ) = LWKOPT 247:ENDIF248:*249:IF( INFO.NE.0 )THEN250:CALLXERBLA( 'CHESVX', -INFO ) 251:RETURN252:ELSEIF( LQUERY )THEN253:RETURN254:ENDIF255:*256:IF( NOFACT )THEN257:*258:* Compute the factorization A = U*D*U' or A = L*D*L'.259:*260:CALLCLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 261:CALLCHETRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO ) 262:*263:* Return if INFO is non-zero.264:*265:IF( INFO.GT.0 )THEN266: RCOND = ZERO 267:RETURN268:ENDIF269:ENDIF270:*271:* Compute the norm of the matrix A.272:*273: ANORM =CLANHE( 'I', UPLO, N, A, LDA, RWORK ) 274:*275:* Compute the reciprocal of the condition number of A.276:*277:CALLCHECON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, INFO ) 278:*279:* Compute the solution vectors X.280:*281:CALLCLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 282:CALLCHETRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 283:*284:* Use iterative refinement to improve the computed solutions and285:* compute error bounds and backward error estimates for them.286:*287:CALLCHERFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 288: $ LDX, FERR, BERR, WORK, RWORK, INFO ) 289:*290:* Set INFO = N+1 if the matrix is singular to working precision.291:*292:IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 293: $ INFO = N + 1 294:*295:WORK( 1 ) = LWKOPT 296:*297:RETURN298:*299:* End of CHESVX300:*301:END302: