001:SUBROUTINEZHESVXX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, 002: $ EQUED, S, B, LDB, X, LDX, RCOND, RPVGRW, BERR, 003: $ N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, 004: $ NPARAMS, PARAMS, WORK, RWORK, INFO ) 005:*006:* -- LAPACK driver routine (version 3.2.1) --007:* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --008:* -- Jason Riedy of Univ. of California Berkeley. --009:* -- April 2009 --010:*011:* -- LAPACK is a software package provided by Univ. of Tennessee, --012:* -- Univ. of California Berkeley and NAG Ltd. --013:*014:IMPLICITNONE 015:* ..016:* .. Scalar Arguments ..017: CHARACTER EQUED, FACT, UPLO 018: INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS, 019: $ N_ERR_BNDS 020: DOUBLE PRECISION RCOND, RPVGRW 021:* ..022:* .. Array Arguments ..023: INTEGERIPIV( * ) 024: COMPLEX*16A( LDA, * ),AF( LDAF, * ),B( LDB, * ), 025: $WORK( * ),X( LDX, * ) 026: DOUBLE PRECISIONS( * ),PARAMS( * ),BERR( * ),RWORK( * ), 027: $ERR_BNDS_NORM( NRHS, * ), 028: $ERR_BNDS_COMP( NRHS, * ) 029:* ..030:*031:* Purpose032:* =======033:*034:* ZHESVXX uses the diagonal pivoting factorization to compute the035:* solution to a complex*16 system of linear equations A * X = B, where036:* A is an N-by-N symmetric matrix and X and B are N-by-NRHS037:* matrices.038:*039:* If requested, both normwise and maximum componentwise error bounds040:* are returned. ZHESVXX will return a solution with a tiny041:* guaranteed error (O(eps) where eps is the working machine042:* precision) unless the matrix is very ill-conditioned, in which043:* case a warning is returned. Relevant condition numbers also are044:* calculated and returned.045:*046:* ZHESVXX accepts user-provided factorizations and equilibration047:* factors; see the definitions of the FACT and EQUED options.048:* Solving with refinement and using a factorization from a previous049:* ZHESVXX call will also produce a solution with either O(eps)050:* errors or warnings, but we cannot make that claim for general051:* user-provided factorizations and equilibration factors if they052:* differ from what ZHESVXX would itself produce.053:*054:* Description055:* ===========056:*057:* The following steps are performed:058:*059:* 1. If FACT = 'E', double precision scaling factors are computed to equilibrate060:* the system:061:*062:* diag(S)*A*diag(S) *inv(diag(S))*X = diag(S)*B063:*064:* Whether or not the system will be equilibrated depends on the065:* scaling of the matrix A, but if equilibration is used, A is066:* overwritten by diag(S)*A*diag(S) and B by diag(S)*B.067:*068:* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor069:* the matrix A (after equilibration if FACT = 'E') as070:*071:* A = U * D * U**T, if UPLO = 'U', or072:* A = L * D * L**T, if UPLO = 'L',073:*074:* where U (or L) is a product of permutation and unit upper (lower)075:* triangular matrices, and D is symmetric and block diagonal with076:* 1-by-1 and 2-by-2 diagonal blocks.077:*078:* 3. If some D(i,i)=0, so that D is exactly singular, then the079:* routine returns with INFO = i. Otherwise, the factored form of A080:* is used to estimate the condition number of the matrix A (see081:* argument RCOND). If the reciprocal of the condition number is082:* less than machine precision, the routine still goes on to solve083:* for X and compute error bounds as described below.084:*085:* 4. The system of equations is solved for X using the factored form086:* of A.087:*088:* 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),089:* the routine will use iterative refinement to try to get a small090:* error and error bounds. Refinement calculates the residual to at091:* least twice the working precision.092:*093:* 6. If equilibration was used, the matrix X is premultiplied by094:* diag(R) so that it solves the original system before095:* equilibration.096:*097:* Arguments098:* =========099:*100:* Some optional parameters are bundled in the PARAMS array. These101:* settings determine how refinement is performed, but often the102:* defaults are acceptable. If the defaults are acceptable, users103:* can pass NPARAMS = 0 which prevents the source code from accessing104:* the PARAMS argument.105:*106:* FACT (input) CHARACTER*1107:* Specifies whether or not the factored form of the matrix A is108:* supplied on entry, and if not, whether the matrix A should be109:* equilibrated before it is factored.110:* = 'F': On entry, AF and IPIV contain the factored form of A.111:* If EQUED is not 'N', the matrix A has been112:* equilibrated with scaling factors given by S.113:* A, AF, and IPIV are not modified.114:* = 'N': The matrix A will be copied to AF and factored.115:* = 'E': The matrix A will be equilibrated if necessary, then116:* copied to AF and factored.117:*118:* N (input) INTEGER119:* The number of linear equations, i.e., the order of the120:* matrix A. N >= 0.121:*122:* NRHS (input) INTEGER123:* The number of right hand sides, i.e., the number of columns124:* of the matrices B and X. NRHS >= 0.125:*126:* A (input/output) COMPLEX*16 array, dimension (LDA,N)127:* The symmetric matrix A. If UPLO = 'U', the leading N-by-N128:* upper triangular part of A contains the upper triangular129:* part of the matrix A, and the strictly lower triangular130:* part of A is not referenced. If UPLO = 'L', the leading131:* N-by-N lower triangular part of A contains the lower132:* triangular part of the matrix A, and the strictly upper133:* triangular part of A is not referenced.134:*135:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by136:* diag(S)*A*diag(S).137:*138:* LDA (input) INTEGER139:* The leading dimension of the array A. LDA >= max(1,N).140:*141:* AF (input or output) COMPLEX*16 array, dimension (LDAF,N)142:* If FACT = 'F', then AF is an input argument and on entry143:* contains the block diagonal matrix D and the multipliers144:* used to obtain the factor U or L from the factorization A =145:* U*D*U**T or A = L*D*L**T as computed by DSYTRF.146:*147:* If FACT = 'N', then AF is an output argument and on exit148:* returns the block diagonal matrix D and the multipliers149:* used to obtain the factor U or L from the factorization A =150:* U*D*U**T or A = L*D*L**T.151:*152:* LDAF (input) INTEGER153:* The leading dimension of the array AF. LDAF >= max(1,N).154:*155:* IPIV (input or output) INTEGER array, dimension (N)156:* If FACT = 'F', then IPIV is an input argument and on entry157:* contains details of the interchanges and the block158:* structure of D, as determined by ZHETRF. If IPIV(k) > 0,159:* then rows and columns k and IPIV(k) were interchanged and160:* D(k,k) is a 1-by-1 diagonal block. If UPLO = 'U' and161:* IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and162:* -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2163:* diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0,164:* then rows and columns k+1 and -IPIV(k) were interchanged165:* and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.166:*167:* If FACT = 'N', then IPIV is an output argument and on exit168:* contains details of the interchanges and the block169:* structure of D, as determined by ZHETRF.170:*171:* EQUED (input or output) CHARACTER*1172:* Specifies the form of equilibration that was done.173:* = 'N': No equilibration (always true if FACT = 'N').174:* = 'Y': Both row and column equilibration, i.e., A has been175:* replaced by diag(S) * A * diag(S).176:* EQUED is an input argument if FACT = 'F'; otherwise, it is an177:* output argument.178:*179:* S (input or output) DOUBLE PRECISION array, dimension (N)180:* The scale factors for A. If EQUED = 'Y', A is multiplied on181:* the left and right by diag(S). S is an input argument if FACT =182:* 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED183:* = 'Y', each element of S must be positive. If S is output, each184:* element of S is a power of the radix. If S is input, each element185:* of S should be a power of the radix to ensure a reliable solution186:* and error estimates. Scaling by powers of the radix does not cause187:* rounding errors unless the result underflows or overflows.188:* Rounding errors during scaling lead to refining with a matrix that189:* is not equivalent to the input matrix, producing error estimates190:* that may not be reliable.191:*192:* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)193:* On entry, the N-by-NRHS right hand side matrix B.194:* On exit,195:* if EQUED = 'N', B is not modified;196:* if EQUED = 'Y', B is overwritten by diag(S)*B;197:*198:* LDB (input) INTEGER199:* The leading dimension of the array B. LDB >= max(1,N).200:*201:* X (output) COMPLEX*16 array, dimension (LDX,NRHS)202:* If INFO = 0, the N-by-NRHS solution matrix X to the original203:* system of equations. Note that A and B are modified on exit if204:* EQUED .ne. 'N', and the solution to the equilibrated system is205:* inv(diag(S))*X.206:*207:* LDX (input) INTEGER208:* The leading dimension of the array X. LDX >= max(1,N).209:*210:* RCOND (output) DOUBLE PRECISION211:* Reciprocal scaled condition number. This is an estimate of the212:* reciprocal Skeel condition number of the matrix A after213:* equilibration (if done). If this is less than the machine214:* precision (in particular, if it is zero), the matrix is singular215:* to working precision. Note that the error may still be small even216:* if this number is very small and the matrix appears ill-217:* conditioned.218:*219:* RPVGRW (output) DOUBLE PRECISION220:* Reciprocal pivot growth. On exit, this contains the reciprocal221:* pivot growth factor norm(A)/norm(U). The "max absolute element"222:* norm is used. If this is much less than 1, then the stability of223:* the LU factorization of the (equilibrated) matrix A could be poor.224:* This also means that the solution X, estimated condition numbers,225:* and error bounds could be unreliable. If factorization fails with226:* 0<INFO<=N, then this contains the reciprocal pivot growth factor227:* for the leading INFO columns of A.228:*229:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)230:* Componentwise relative backward error. This is the231:* componentwise relative backward error of each solution vector X(j)232:* (i.e., the smallest relative change in any element of A or B that233:* makes X(j) an exact solution).234:*235:* N_ERR_BNDS (input) INTEGER236:* Number of error bounds to return for each right hand side237:* and each type (normwise or componentwise). See ERR_BNDS_NORM and238:* ERR_BNDS_COMP below.239:*240:* ERR_BNDS_NORM (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)241:* For each right-hand side, this array contains information about242:* various error bounds and condition numbers corresponding to the243:* normwise relative error, which is defined as follows:244:*245:* Normwise relative error in the ith solution vector:246:* max_j (abs(XTRUE(j,i) - X(j,i)))247:* ------------------------------248:* max_j abs(X(j,i))249:*250:* The array is indexed by the type of error information as described251:* below. There currently are up to three pieces of information252:* returned.253:*254:* The first index in ERR_BNDS_NORM(i,:) corresponds to the ith255:* right-hand side.256:*257:* The second index in ERR_BNDS_NORM(:,err) contains the following258:* three fields:259:* err = 1 "Trust/don't trust" boolean. Trust the answer if the260:* reciprocal condition number is less than the threshold261:* sqrt(n) * dlamch('Epsilon').262:*263:* err = 2 "Guaranteed" error bound: The estimated forward error,264:* almost certainly within a factor of 10 of the true error265:* so long as the next entry is greater than the threshold266:* sqrt(n) * dlamch('Epsilon'). This error bound should only267:* be trusted if the previous boolean is true.268:*269:* err = 3 Reciprocal condition number: Estimated normwise270:* reciprocal condition number. Compared with the threshold271:* sqrt(n) * dlamch('Epsilon') to determine if the error272:* estimate is "guaranteed". These reciprocal condition273:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some274:* appropriately scaled matrix Z.275:* Let Z = S*A, where S scales each row by a power of the276:* radix so all absolute row sums of Z are approximately 1.277:*278:* See Lapack Working Note 165 for further details and extra279:* cautions.280:*281:* ERR_BNDS_COMP (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)282:* For each right-hand side, this array contains information about283:* various error bounds and condition numbers corresponding to the284:* componentwise relative error, which is defined as follows:285:*286:* Componentwise relative error in the ith solution vector:287:* abs(XTRUE(j,i) - X(j,i))288:* max_j ----------------------289:* abs(X(j,i))290:*291:* The array is indexed by the right-hand side i (on which the292:* componentwise relative error depends), and the type of error293:* information as described below. There currently are up to three294:* pieces of information returned for each right-hand side. If295:* componentwise accuracy is not requested (PARAMS(3) = 0.0), then296:* ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most297:* the first (:,N_ERR_BNDS) entries are returned.298:*299:* The first index in ERR_BNDS_COMP(i,:) corresponds to the ith300:* right-hand side.301:*302:* The second index in ERR_BNDS_COMP(:,err) contains the following303:* three fields:304:* err = 1 "Trust/don't trust" boolean. Trust the answer if the305:* reciprocal condition number is less than the threshold306:* sqrt(n) * dlamch('Epsilon').307:*308:* err = 2 "Guaranteed" error bound: The estimated forward error,309:* almost certainly within a factor of 10 of the true error310:* so long as the next entry is greater than the threshold311:* sqrt(n) * dlamch('Epsilon'). This error bound should only312:* be trusted if the previous boolean is true.313:*314:* err = 3 Reciprocal condition number: Estimated componentwise315:* reciprocal condition number. Compared with the threshold316:* sqrt(n) * dlamch('Epsilon') to determine if the error317:* estimate is "guaranteed". These reciprocal condition318:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some319:* appropriately scaled matrix Z.320:* Let Z = S*(A*diag(x)), where x is the solution for the321:* current right-hand side and S scales each row of322:* A*diag(x) by a power of the radix so all absolute row323:* sums of Z are approximately 1.324:*325:* See Lapack Working Note 165 for further details and extra326:* cautions.327:*328:* NPARAMS (input) INTEGER329:* Specifies the number of parameters set in PARAMS. If .LE. 0, the330:* PARAMS array is never referenced and default values are used.331:*332:* PARAMS (input / output) DOUBLE PRECISION array, dimension NPARAMS333:* Specifies algorithm parameters. If an entry is .LT. 0.0, then334:* that entry will be filled with default value used for that335:* parameter. Only positions up to NPARAMS are accessed; defaults336:* are used for higher-numbered parameters.337:*338:* PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative339:* refinement or not.340:* Default: 1.0D+0341:* = 0.0 : No refinement is performed, and no error bounds are342:* computed.343:* = 1.0 : Use the extra-precise refinement algorithm.344:* (other values are reserved for future use)345:*346:* PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual347:* computations allowed for refinement.348:* Default: 10349:* Aggressive: Set to 100 to permit convergence using approximate350:* factorizations or factorizations other than LU. If351:* the factorization uses a technique other than352:* Gaussian elimination, the guarantees in353:* err_bnds_norm and err_bnds_comp may no longer be354:* trustworthy.355:*356:* PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code357:* will attempt to find a solution with small componentwise358:* relative error in the double-precision algorithm. Positive359:* is true, 0.0 is false.360:* Default: 1.0 (attempt componentwise convergence)361:*362:* WORK (workspace) COMPLEX*16 array, dimension (2*N)363:*364:* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)365:*366:* INFO (output) INTEGER367:* = 0: Successful exit. The solution to every right-hand side is368:* guaranteed.369:* < 0: If INFO = -i, the i-th argument had an illegal value370:* > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization371:* has been completed, but the factor U is exactly singular, so372:* the solution and error bounds could not be computed. RCOND = 0373:* is returned.374:* = N+J: The solution corresponding to the Jth right-hand side is375:* not guaranteed. The solutions corresponding to other right-376:* hand sides K with K > J may not be guaranteed as well, but377:* only the first such right-hand side is reported. If a small378:* componentwise error is not requested (PARAMS(3) = 0.0) then379:* the Jth right-hand side is the first with a normwise error380:* bound that is not guaranteed (the smallest J such381:* that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)382:* the Jth right-hand side is the first with either a normwise or383:* componentwise error bound that is not guaranteed (the smallest384:* J such that either ERR_BNDS_NORM(J,1) = 0.0 or385:* ERR_BNDS_COMP(J,1) = 0.0). See the definition of386:* ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information387:* about all of the right-hand sides check ERR_BNDS_NORM or388:* ERR_BNDS_COMP.389:*390:* ==================================================================391:*392:* .. Parameters ..393: DOUBLE PRECISION ZERO, ONE 394:PARAMETER( ZERO = 0.0D+0, ONE = 1.0D+0 ) 395: INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 396: INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 397: INTEGER CMP_ERR_I, PIV_GROWTH_I 398:PARAMETER( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 399: $ BERR_I = 3 ) 400:PARAMETER( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 401:PARAMETER( CMP_RCOND_I = 7, CMP_ERR_I = 8, 402: $ PIV_GROWTH_I = 9 ) 403:* ..404:* .. Local Scalars ..405:LOGICALEQUIL, NOFACT, RCEQU 406: INTEGER INFEQU, J 407: DOUBLE PRECISION AMAX, BIGNUM, SMIN, SMAX, SCOND, SMLNUM 408:* ..409:* .. External Functions ..410:EXTERNALLSAME, DLAMCH, ZLA_HERPVGRW 411:LOGICALLSAME 412: DOUBLE PRECISION DLAMCH, ZLA_HERPVGRW 413:* ..414:* .. External Subroutines ..415:EXTERNALZHECON, ZHEEQUB, ZHETRF, ZHETRS, ZLACPY, 416: $ ZLAQHE, XERBLA, ZLASCL2, ZHERFSX 417:* ..418:* .. Intrinsic Functions ..419:INTRINSICMAX, MIN 420:* ..421:* .. Executable Statements ..422:*423: INFO = 0 424: NOFACT =LSAME( FACT, 'N' ) 425: EQUIL =LSAME( FACT, 'E' ) 426: SMLNUM =DLAMCH( 'Safe minimum' ) 427: BIGNUM = ONE / SMLNUM 428:IF( NOFACT .OR. EQUIL )THEN429: EQUED = 'N' 430: RCEQU = .FALSE. 431:ELSE432: RCEQU =LSAME( EQUED, 'Y' ) 433:ENDIF434:*435:* Default is failure. If an input parameter is wrong or436:* factorization fails, make everything look horrible. Only the437:* pivot growth is set here, the rest is initialized in ZHERFSX.438:*439: RPVGRW = ZERO 440:*441:* Test the input parameters. PARAMS is not tested until ZHERFSX.442:*443:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT. 444: $LSAME( FACT, 'F' ) )THEN445: INFO = -1 446:ELSEIF( .NOT.LSAME( UPLO, 'U' ) .AND. 447: $ .NOT.LSAME( UPLO, 'L' ) )THEN448: INFO = -2 449:ELSEIF( N.LT.0 )THEN450: INFO = -3 451:ELSEIF( NRHS.LT.0 )THEN452: INFO = -4 453:ELSEIF( LDA.LT.MAX( 1, N ) )THEN454: INFO = -6 455:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN456: INFO = -8 457:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 458: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN459: INFO = -9 460:ELSE461:IF( RCEQU )THEN462: SMIN = BIGNUM 463: SMAX = ZERO 464:DO10 J = 1, N 465: SMIN =MIN( SMIN,S( J ) ) 466: SMAX =MAX( SMAX,S( J ) ) 467: 10CONTINUE468:IF( SMIN.LE.ZERO )THEN469: INFO = -10 470:ELSEIF( N.GT.0 )THEN471: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 472:ELSE473: SCOND = ONE 474:ENDIF475:ENDIF476:IF( INFO.EQ.0 )THEN477:IF( LDB.LT.MAX( 1, N ) )THEN478: INFO = -12 479:ELSEIF( LDX.LT.MAX( 1, N ) )THEN480: INFO = -14 481:ENDIF482:ENDIF483:ENDIF484:*485:IF( INFO.NE.0 )THEN486:CALLXERBLA( 'ZHESVXX', -INFO ) 487:RETURN488:ENDIF489:*490:IF( EQUIL )THEN491:*492:* Compute row and column scalings to equilibrate the matrix A.493:*494:CALLZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFEQU ) 495:IF( INFEQU.EQ.0 )THEN496:*497:* Equilibrate the matrix.498:*499:CALLZLAQHE( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 500: RCEQU =LSAME( EQUED, 'Y' ) 501:ENDIF502:ENDIF503:*504:* Scale the right-hand side.505:*506:IF( RCEQU )CALLZLASCL2( N, NRHS, S, B, LDB ) 507:*508:IF( NOFACT .OR. EQUIL )THEN509:*510:* Compute the LU factorization of A.511:*512:CALLZLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 513:CALLZHETRF( UPLO, N, AF, LDAF, IPIV, WORK, 5*MAX(1,N), INFO ) 514:*515:* Return if INFO is non-zero.516:*517:IF( INFO.GT.0 )THEN518:*519:* Pivot in column INFO is exactly 0520:* Compute the reciprocal pivot growth factor of the521:* leading rank-deficient INFO columns of A.522:*523:IF( N.GT.0 ) 524: $ RPVGRW =ZLA_HERPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF, 525: $ IPIV, RWORK ) 526:RETURN527:ENDIF528:ENDIF529:*530:* Compute the reciprocal pivot growth factor RPVGRW.531:*532:IF( N.GT.0 ) 533: $ RPVGRW =ZLA_HERPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF, IPIV, 534: $ RWORK ) 535:*536:* Compute the solution matrix X.537:*538:CALLZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 539:CALLZHETRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 540:*541:* Use iterative refinement to improve the computed solution and542:* compute error bounds and backward error estimates for it.543:*544:CALLZHERFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, 545: $ S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, 546: $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO ) 547:*548:* Scale solutions.549:*550:IF( RCEQU )THEN551:CALLZLASCL2( N, NRHS, S, X, LDX ) 552:ENDIF553:*554:RETURN555:*556:* End of ZHESVXX557:*558:END559: