001:SUBROUTINEDSYSVXX( 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, IWORK, INFO ) 005:*006:* -- LAPACK driver routine (version 3.2) --007:* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --008:* -- Jason Riedy of Univ. of California Berkeley. --009:* -- November 2008 --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( * ),IWORK( * ) 024: DOUBLE PRECISIONA( LDA, * ),AF( LDAF, * ),B( LDB, * ), 025: $X( LDX, * ),WORK( * ) 026: DOUBLE PRECISIONS( * ),PARAMS( * ),BERR( * ), 027: $ERR_BNDS_NORM( NRHS, * ), 028: $ERR_BNDS_COMP( NRHS, * ) 029:* ..030:*031:* Purpose032:* =======033:*034:* DSYSVXX uses the diagonal pivoting factorization to compute the035:* solution to a double precision system of linear equations A * X = B, where A036:* is an N-by-N symmetric matrix and X and B are N-by-NRHS matrices.037:*038:* If requested, both normwise and maximum componentwise error bounds039:* are returned. DSYSVXX will return a solution with a tiny040:* guaranteed error (O(eps) where eps is the working machine041:* precision) unless the matrix is very ill-conditioned, in which042:* case a warning is returned. Relevant condition numbers also are043:* calculated and returned.044:*045:* DSYSVXX accepts user-provided factorizations and equilibration046:* factors; see the definitions of the FACT and EQUED options.047:* Solving with refinement and using a factorization from a previous048:* DSYSVXX call will also produce a solution with either O(eps)049:* errors or warnings, but we cannot make that claim for general050:* user-provided factorizations and equilibration factors if they051:* differ from what DSYSVXX would itself produce.052:*053:* Description054:* ===========055:*056:* The following steps are performed:057:*058:* 1. If FACT = 'E', double precision scaling factors are computed to equilibrate059:* the system:060:*061:* diag(S)*A*diag(S) *inv(diag(S))*X = diag(S)*B062:*063:* Whether or not the system will be equilibrated depends on the064:* scaling of the matrix A, but if equilibration is used, A is065:* overwritten by diag(S)*A*diag(S) and B by diag(S)*B.066:*067:* 2. If FACT = 'N' or 'E', the LU decomposition is used to factor068:* the matrix A (after equilibration if FACT = 'E') as069:*070:* A = U * D * U**T, if UPLO = 'U', or071:* A = L * D * L**T, if UPLO = 'L',072:*073:* where U (or L) is a product of permutation and unit upper (lower)074:* triangular matrices, and D is symmetric and block diagonal with075:* 1-by-1 and 2-by-2 diagonal blocks.076:*077:* 3. If some D(i,i)=0, so that D is exactly singular, then the078:* routine returns with INFO = i. Otherwise, the factored form of A079:* is used to estimate the condition number of the matrix A (see080:* argument RCOND). If the reciprocal of the condition number is081:* less than machine precision, the routine still goes on to solve082:* for X and compute error bounds as described below.083:*084:* 4. The system of equations is solved for X using the factored form085:* of A.086:*087:* 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),088:* the routine will use iterative refinement to try to get a small089:* error and error bounds. Refinement calculates the residual to at090:* least twice the working precision.091:*092:* 6. If equilibration was used, the matrix X is premultiplied by093:* diag(R) so that it solves the original system before094:* equilibration.095:*096:* Arguments097:* =========098:*099:* Some optional parameters are bundled in the PARAMS array. These100:* settings determine how refinement is performed, but often the101:* defaults are acceptable. If the defaults are acceptable, users102:* can pass NPARAMS = 0 which prevents the source code from accessing103:* the PARAMS argument.104:*105:* FACT (input) CHARACTER*1106:* Specifies whether or not the factored form of the matrix A is107:* supplied on entry, and if not, whether the matrix A should be108:* equilibrated before it is factored.109:* = 'F': On entry, AF and IPIV contain the factored form of A.110:* If EQUED is not 'N', the matrix A has been111:* equilibrated with scaling factors given by S.112:* A, AF, and IPIV are not modified.113:* = 'N': The matrix A will be copied to AF and factored.114:* = 'E': The matrix A will be equilibrated if necessary, then115:* copied to AF and factored.116:*117:* N (input) INTEGER118:* The number of linear equations, i.e., the order of the119:* matrix A. N >= 0.120:*121:* NRHS (input) INTEGER122:* The number of right hand sides, i.e., the number of columns123:* of the matrices B and X. NRHS >= 0.124:*125:* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)126:* The symmetric matrix A. If UPLO = 'U', the leading N-by-N127:* upper triangular part of A contains the upper triangular128:* part of the matrix A, and the strictly lower triangular129:* part of A is not referenced. If UPLO = 'L', the leading130:* N-by-N lower triangular part of A contains the lower131:* triangular part of the matrix A, and the strictly upper132:* triangular part of A is not referenced.133:*134:* On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by135:* diag(S)*A*diag(S).136:*137:* LDA (input) INTEGER138:* The leading dimension of the array A. LDA >= max(1,N).139:*140:* AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N)141:* If FACT = 'F', then AF is an input argument and on entry142:* contains the block diagonal matrix D and the multipliers143:* used to obtain the factor U or L from the factorization A =144:* U*D*U**T or A = L*D*L**T as computed by DSYTRF.145:*146:* If FACT = 'N', then AF is an output argument and on exit147:* returns the block diagonal matrix D and the multipliers148:* used to obtain the factor U or L from the factorization A =149:* U*D*U**T or A = L*D*L**T.150:*151:* LDAF (input) INTEGER152:* The leading dimension of the array AF. LDAF >= max(1,N).153:*154:* IPIV (input or output) INTEGER array, dimension (N)155:* If FACT = 'F', then IPIV is an input argument and on entry156:* contains details of the interchanges and the block157:* structure of D, as determined by DSYTRF. If IPIV(k) > 0,158:* then rows and columns k and IPIV(k) were interchanged and159:* D(k,k) is a 1-by-1 diagonal block. If UPLO = 'U' and160:* IPIV(k) = IPIV(k-1) < 0, then rows and columns k-1 and161:* -IPIV(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2162:* diagonal block. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0,163:* then rows and columns k+1 and -IPIV(k) were interchanged164:* and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.165:*166:* If FACT = 'N', then IPIV is an output argument and on exit167:* contains details of the interchanges and the block168:* structure of D, as determined by DSYTRF.169:*170:* EQUED (input or output) CHARACTER*1171:* Specifies the form of equilibration that was done.172:* = 'N': No equilibration (always true if FACT = 'N').173:* = 'Y': Both row and column equilibration, i.e., A has been174:* replaced by diag(S) * A * diag(S).175:* EQUED is an input argument if FACT = 'F'; otherwise, it is an176:* output argument.177:*178:* S (input or output) DOUBLE PRECISION array, dimension (N)179:* The scale factors for A. If EQUED = 'Y', A is multiplied on180:* the left and right by diag(S). S is an input argument if FACT =181:* 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED182:* = 'Y', each element of S must be positive. If S is output, each183:* element of S is a power of the radix. If S is input, each element184:* of S should be a power of the radix to ensure a reliable solution185:* and error estimates. Scaling by powers of the radix does not cause186:* rounding errors unless the result underflows or overflows.187:* Rounding errors during scaling lead to refining with a matrix that188:* is not equivalent to the input matrix, producing error estimates189:* that may not be reliable.190:*191:* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)192:* On entry, the N-by-NRHS right hand side matrix B.193:* On exit,194:* if EQUED = 'N', B is not modified;195:* if EQUED = 'Y', B is overwritten by diag(S)*B;196:*197:* LDB (input) INTEGER198:* The leading dimension of the array B. LDB >= max(1,N).199:*200:* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)201:* If INFO = 0, the N-by-NRHS solution matrix X to the original202:* system of equations. Note that A and B are modified on exit if203:* EQUED .ne. 'N', and the solution to the equilibrated system is204:* inv(diag(S))*X.205:*206:* LDX (input) INTEGER207:* The leading dimension of the array X. LDX >= max(1,N).208:*209:* RCOND (output) DOUBLE PRECISION210:* Reciprocal scaled condition number. This is an estimate of the211:* reciprocal Skeel condition number of the matrix A after212:* equilibration (if done). If this is less than the machine213:* precision (in particular, if it is zero), the matrix is singular214:* to working precision. Note that the error may still be small even215:* if this number is very small and the matrix appears ill-216:* conditioned.217:*218:* RPVGRW (output) DOUBLE PRECISION219:* Reciprocal pivot growth. On exit, this contains the reciprocal220:* pivot growth factor norm(A)/norm(U). The "max absolute element"221:* norm is used. If this is much less than 1, then the stability of222:* the LU factorization of the (equilibrated) matrix A could be poor.223:* This also means that the solution X, estimated condition numbers,224:* and error bounds could be unreliable. If factorization fails with225:* 0<INFO<=N, then this contains the reciprocal pivot growth factor226:* for the leading INFO columns of A.227:*228:* BERR (output) DOUBLE PRECISION array, dimension (NRHS)229:* Componentwise relative backward error. This is the230:* componentwise relative backward error of each solution vector X(j)231:* (i.e., the smallest relative change in any element of A or B that232:* makes X(j) an exact solution).233:*234:* N_ERR_BNDS (input) INTEGER235:* Number of error bounds to return for each right hand side236:* and each type (normwise or componentwise). See ERR_BNDS_NORM and237:* ERR_BNDS_COMP below.238:*239:* ERR_BNDS_NORM (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)240:* For each right-hand side, this array contains information about241:* various error bounds and condition numbers corresponding to the242:* normwise relative error, which is defined as follows:243:*244:* Normwise relative error in the ith solution vector:245:* max_j (abs(XTRUE(j,i) - X(j,i)))246:* ------------------------------247:* max_j abs(X(j,i))248:*249:* The array is indexed by the type of error information as described250:* below. There currently are up to three pieces of information251:* returned.252:*253:* The first index in ERR_BNDS_NORM(i,:) corresponds to the ith254:* right-hand side.255:*256:* The second index in ERR_BNDS_NORM(:,err) contains the following257:* three fields:258:* err = 1 "Trust/don't trust" boolean. Trust the answer if the259:* reciprocal condition number is less than the threshold260:* sqrt(n) * dlamch('Epsilon').261:*262:* err = 2 "Guaranteed" error bound: The estimated forward error,263:* almost certainly within a factor of 10 of the true error264:* so long as the next entry is greater than the threshold265:* sqrt(n) * dlamch('Epsilon'). This error bound should only266:* be trusted if the previous boolean is true.267:*268:* err = 3 Reciprocal condition number: Estimated normwise269:* reciprocal condition number. Compared with the threshold270:* sqrt(n) * dlamch('Epsilon') to determine if the error271:* estimate is "guaranteed". These reciprocal condition272:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some273:* appropriately scaled matrix Z.274:* Let Z = S*A, where S scales each row by a power of the275:* radix so all absolute row sums of Z are approximately 1.276:*277:* See Lapack Working Note 165 for further details and extra278:* cautions.279:*280:* ERR_BNDS_COMP (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)281:* For each right-hand side, this array contains information about282:* various error bounds and condition numbers corresponding to the283:* componentwise relative error, which is defined as follows:284:*285:* Componentwise relative error in the ith solution vector:286:* abs(XTRUE(j,i) - X(j,i))287:* max_j ----------------------288:* abs(X(j,i))289:*290:* The array is indexed by the right-hand side i (on which the291:* componentwise relative error depends), and the type of error292:* information as described below. There currently are up to three293:* pieces of information returned for each right-hand side. If294:* componentwise accuracy is not requested (PARAMS(3) = 0.0), then295:* ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most296:* the first (:,N_ERR_BNDS) entries are returned.297:*298:* The first index in ERR_BNDS_COMP(i,:) corresponds to the ith299:* right-hand side.300:*301:* The second index in ERR_BNDS_COMP(:,err) contains the following302:* three fields:303:* err = 1 "Trust/don't trust" boolean. Trust the answer if the304:* reciprocal condition number is less than the threshold305:* sqrt(n) * dlamch('Epsilon').306:*307:* err = 2 "Guaranteed" error bound: The estimated forward error,308:* almost certainly within a factor of 10 of the true error309:* so long as the next entry is greater than the threshold310:* sqrt(n) * dlamch('Epsilon'). This error bound should only311:* be trusted if the previous boolean is true.312:*313:* err = 3 Reciprocal condition number: Estimated componentwise314:* reciprocal condition number. Compared with the threshold315:* sqrt(n) * dlamch('Epsilon') to determine if the error316:* estimate is "guaranteed". These reciprocal condition317:* numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some318:* appropriately scaled matrix Z.319:* Let Z = S*(A*diag(x)), where x is the solution for the320:* current right-hand side and S scales each row of321:* A*diag(x) by a power of the radix so all absolute row322:* sums of Z are approximately 1.323:*324:* See Lapack Working Note 165 for further details and extra325:* cautions.326:*327:* NPARAMS (input) INTEGER328:* Specifies the number of parameters set in PARAMS. If .LE. 0, the329:* PARAMS array is never referenced and default values are used.330:*331:* PARAMS (input / output) DOUBLE PRECISION array, dimension NPARAMS332:* Specifies algorithm parameters. If an entry is .LT. 0.0, then333:* that entry will be filled with default value used for that334:* parameter. Only positions up to NPARAMS are accessed; defaults335:* are used for higher-numbered parameters.336:*337:* PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative338:* refinement or not.339:* Default: 1.0D+0340:* = 0.0 : No refinement is performed, and no error bounds are341:* computed.342:* = 1.0 : Use the extra-precise refinement algorithm.343:* (other values are reserved for future use)344:*345:* PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual346:* computations allowed for refinement.347:* Default: 10348:* Aggressive: Set to 100 to permit convergence using approximate349:* factorizations or factorizations other than LU. If350:* the factorization uses a technique other than351:* Gaussian elimination, the guarantees in352:* err_bnds_norm and err_bnds_comp may no longer be353:* trustworthy.354:*355:* PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code356:* will attempt to find a solution with small componentwise357:* relative error in the double-precision algorithm. Positive358:* is true, 0.0 is false.359:* Default: 1.0 (attempt componentwise convergence)360:*361:* WORK (workspace) DOUBLE PRECISION array, dimension (4*N)362:*363:* IWORK (workspace) INTEGER array, dimension (N)364:*365:* INFO (output) INTEGER366:* = 0: Successful exit. The solution to every right-hand side is367:* guaranteed.368:* < 0: If INFO = -i, the i-th argument had an illegal value369:* > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization370:* has been completed, but the factor U is exactly singular, so371:* the solution and error bounds could not be computed. RCOND = 0372:* is returned.373:* = N+J: The solution corresponding to the Jth right-hand side is374:* not guaranteed. The solutions corresponding to other right-375:* hand sides K with K > J may not be guaranteed as well, but376:* only the first such right-hand side is reported. If a small377:* componentwise error is not requested (PARAMS(3) = 0.0) then378:* the Jth right-hand side is the first with a normwise error379:* bound that is not guaranteed (the smallest J such380:* that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)381:* the Jth right-hand side is the first with either a normwise or382:* componentwise error bound that is not guaranteed (the smallest383:* J such that either ERR_BNDS_NORM(J,1) = 0.0 or384:* ERR_BNDS_COMP(J,1) = 0.0). See the definition of385:* ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information386:* about all of the right-hand sides check ERR_BNDS_NORM or387:* ERR_BNDS_COMP.388:*389:* ==================================================================390:*391:* .. Parameters ..392: DOUBLE PRECISION ZERO, ONE 393:PARAMETER( ZERO = 0.0D+0, ONE = 1.0D+0 ) 394: INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I 395: INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I 396: INTEGER CMP_ERR_I, PIV_GROWTH_I 397:PARAMETER( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2, 398: $ BERR_I = 3 ) 399:PARAMETER( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 ) 400:PARAMETER( CMP_RCOND_I = 7, CMP_ERR_I = 8, 401: $ PIV_GROWTH_I = 9 ) 402:* ..403:* .. Local Scalars ..404:LOGICALEQUIL, NOFACT, RCEQU 405: INTEGER INFEQU, J 406: DOUBLE PRECISION AMAX, BIGNUM, SMIN, SMAX, SCOND, SMLNUM 407:* ..408:* .. External Functions ..409:EXTERNALLSAME, DLAMCH, DLA_SYRPVGRW 410:LOGICALLSAME 411: DOUBLE PRECISION DLAMCH, DLA_SYRPVGRW 412:* ..413:* .. External Subroutines ..414:EXTERNALDSYCON, DSYEQUB, DSYTRF, DSYTRS, 415: $ DLACPY, DLAQSY, XERBLA, DLASCL2, DSYRFSX 416:* ..417:* .. Intrinsic Functions ..418:INTRINSICMAX, MIN 419:* ..420:* .. Executable Statements ..421:*422: INFO = 0 423: NOFACT =LSAME( FACT, 'N' ) 424: EQUIL =LSAME( FACT, 'E' ) 425: SMLNUM =DLAMCH( 'Safe minimum' ) 426: BIGNUM = ONE / SMLNUM 427:IF( NOFACT .OR. EQUIL )THEN428: EQUED = 'N' 429: RCEQU = .FALSE. 430:ELSE431: RCEQU =LSAME( EQUED, 'Y' ) 432:ENDIF433:*434:* Default is failure. If an input parameter is wrong or435:* factorization fails, make everything look horrible. Only the436:* pivot growth is set here, the rest is initialized in DSYRFSX.437:*438: RPVGRW = ZERO 439:*440:* Test the input parameters. PARAMS is not tested until DSYRFSX.441:*442:IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT. 443: $LSAME( FACT, 'F' ) )THEN444: INFO = -1 445:ELSEIF( .NOT.LSAME(UPLO, 'U') .AND. 446: $ .NOT.LSAME(UPLO, 'L') )THEN447: INFO = -2 448:ELSEIF( N.LT.0 )THEN449: INFO = -3 450:ELSEIF( NRHS.LT.0 )THEN451: INFO = -4 452:ELSEIF( LDA.LT.MAX( 1, N ) )THEN453: INFO = -6 454:ELSEIF( LDAF.LT.MAX( 1, N ) )THEN455: INFO = -8 456:ELSEIF(LSAME( FACT, 'F' ) .AND. .NOT. 457: $ ( RCEQU .OR.LSAME( EQUED, 'N' ) ) )THEN458: INFO = -9 459:ELSE460:IF( RCEQU )THEN461: SMIN = BIGNUM 462: SMAX = ZERO 463:DO10 J = 1, N 464: SMIN =MIN( SMIN,S( J ) ) 465: SMAX =MAX( SMAX,S( J ) ) 466: 10CONTINUE467:IF( SMIN.LE.ZERO )THEN468: INFO = -10 469:ELSEIF( N.GT.0 )THEN470: SCOND =MAX( SMIN, SMLNUM ) /MIN( SMAX, BIGNUM ) 471:ELSE472: SCOND = ONE 473:ENDIF474:ENDIF475:IF( INFO.EQ.0 )THEN476:IF( LDB.LT.MAX( 1, N ) )THEN477: INFO = -12 478:ELSEIF( LDX.LT.MAX( 1, N ) )THEN479: INFO = -14 480:ENDIF481:ENDIF482:ENDIF483:*484:IF( INFO.NE.0 )THEN485:CALLXERBLA( 'DSYSVXX', -INFO ) 486:RETURN487:ENDIF488:*489:IF( EQUIL )THEN490:*491:* Compute row and column scalings to equilibrate the matrix A.492:*493:CALLDSYEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFEQU ) 494:IF( INFEQU.EQ.0 )THEN495:*496:* Equilibrate the matrix.497:*498:CALLDLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED ) 499: RCEQU =LSAME( EQUED, 'Y' ) 500:ENDIF501:ENDIF502:*503:* Scale the right-hand side.504:*505:IF( RCEQU )CALLDLASCL2( N, NRHS, S, B, LDB ) 506:*507:IF( NOFACT .OR. EQUIL )THEN508:*509:* Compute the LU factorization of A.510:*511:CALLDLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 512:CALLDSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, 5*MAX(1,N), INFO ) 513:*514:* Return if INFO is non-zero.515:*516:IF( INFO.GT.0 )THEN517:*518:* Pivot in column INFO is exactly 0519:* Compute the reciprocal pivot growth factor of the520:* leading rank-deficient INFO columns of A.521:*522:IF( N.GT.0 ) 523: $ RPVGRW =DLA_SYRPVGRW(UPLO, N, INFO, A, LDA, AF, 524: $ LDAF, IPIV, WORK ) 525:RETURN526:ENDIF527:ENDIF528:*529:* Compute the reciprocal pivot growth factor RPVGRW.530:*531:IF( N.GT.0 ) 532: $ RPVGRW =DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF, LDAF, 533: $ IPIV, WORK ) 534:*535:* Compute the solution matrix X.536:*537:CALLDLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 538:CALLDSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 539:*540:* Use iterative refinement to improve the computed solution and541:* compute error bounds and backward error estimates for it.542:*543:CALLDSYRFSX( UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, 544: $ S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, 545: $ ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO ) 546:*547:* Scale solutions.548:*549:IF( RCEQU )THEN550:CALLDLASCL2( N, NRHS, S, X, LDX ) 551:ENDIF552:*553:RETURN554:*555:* End of DSYSVXX556:*557:END558: