C ALGORITHM 648, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 13, NO. 1, P. 28 Brief information for DETEST 1986 version (sufficient for preliminary tests of the code without the timing facility). Implementation Steps -------------------- 1. In routine CONST in CONCLK.F insert a data statement suitable for your machine. (See implementation notes in STDOC.T or NSDOC.T. Sample settings are included in CONST in comments.) 2. Insert system calls to suppress underflow exceptions, and any other needed initialization, in the I=0 section of CONST. Again, sample code is included. 3. The program consisting of NSTTRU.F, NSDTST.F, NSTRUE.F, NSPROB.F and CONCLK.F should now compile and run, giving output similar to NSTTRU.O. 4. Similarly with STTTRU.F, STDTST.F, STTRUE.F, STPROB.F and CONCLK.F. 1. Purpose ------- Two packages are provided for assessing the performance of initial value solvers. The first package, whose main routine is NSDTST, is designed for assessing the performance of solvers suitable for non-stiff systems, while the second package, whose main routine is STDTST, is designed for assessing the performance of solvers suitable for stiff systems. Each package consists of a number of routines but the user need only be aware of the main routine and the routines FCN (and for the stiff package, PDERV), and STATS, whose role is explained below. In this document we will describe the use of STDTST (double precision version). The requirements and calling sequence of NSDTST are almost identical. In [4] the design of the testing package is discussed and guidance is given on the interpretation of the results it produces. A set of test problems, described in detail in [2,3], is incorporated in the stiff package. The code being tested is run on a selection of these problems at various tolerances. The user selects the problems and the tolerances, and also organizes the problems into groups for statistical reporting purposes, at his discretion. To test a code a user must write an interface routine called METHOD, described below, and then call STDTST with the desired options. Note that STDTST comes in a 'single' and a 'double' precision version. It is best if the version used matches the precision of the SOLVER under test. If this is not possible then great care must be exercised when constructing METHOD. The arguments of STDTST are, in any event, always single precision but those of METHOD are of single or double precision according to the version used. The package divides naturally into five parts: STDTST,CNTROL and various service routines organize the assembling, computation and reporting of statistics. STATS is the routine which 'instruments' the code being tested and passes statistics via COMMON to CNTROL and STDTST. FCN, PDERV, IVALU, EVALU describe the set of test problems. FCN gives the r.h.s. f(y) of the ODE system and PDERV gives the Jacobian matrix df/dy. (At present all the problems are posed in autonomous form). IVALU gives the initial conditions, scaling weights and other data about each problem. EVALU gives accurately computed values at the endpoint. DDCOMP and DSOLVE are standard (double precision) LU decomposition and backsolve routines for full matrices, compatible with the layout of the Jacobian produced by PDERV. They are used by TRUE but are available for use by the code being tested if desired. TRUE and its subordinate routines (alias the Addison-Enright code SECDER) form a reliable stiff solver for computing the 'true' global and local solutions when required. There is also a 'dummy' STDTST and STATS to help the user debug his METHOD routine (described below); a utility GENTIM which can be used on each new machine to generate timing data embedded in the code; and a utility GENWT which can be used if ever a user wishes to add further test problems to the set. Main Lines of Calling Hierarchy (user-supplied routines are in boxes) +--------+ | User's |---STDTST---CNTROL-----IVALU |Program | | +--------+ +--------+ | +------+ |'SOLVER'| |---|METHOD|----|(Code |->-+ | +------+ | being | | | | | tested)| | | | +--------+ |---FCN,PDERV | | | | STATS---TRUE--->--+ | +----EVALU We acknowledge valuable recommendations in Shampine's paper [5]. In particular the package will, by default, integrate each system in scaled form, scaling each solution component by its maximum observed value over the range of integration. That is, the change of variable -1 z = D y is done where D = diag(w(1), .., w(n)) and w(i) =max |i-th component of y| over the range. The problem -1 solved is then z' = D f(x,Dz). The weights w(i) were found by an accurate integration of each problem and are embedded in IVALU. Note that this scaling affects the norms which are used in measuring all errors, and thus can have a considerable effect on the accuracy in some of the problems. If the problem code in IDLIST (see below) is given a negative sign the system is solved in its 'natural' scaling, as was done in the 1975 version of DETEST. 2. Arguments to STDTST: --------- -- ------- TITLE (input) Character of length 80, holds name of method being tested. OPTION (input) Integer array of length 10, only elements 1 to 3 are used and are referred to henceforth as OPT, NORMEF and NRMTYP. (OPTION(4) is also used when OPT=4) OPT one of 1, 2, 3 or 4. OPT selects level of analysis required: 1 gives a report of the following at each tolerance used: - Total time per integration - Overhead time excluding function and Jacobian calls and matrix factorizations. - Number of function calls, Jacobian calls, matrix factorizations and successful steps over range - Global error at endpoint XEND, divided by TOL, ie. ||(computed y) - (true y)||/TOL at x=XEND The norm used throughout the package is that chosen by NRMTYP. 2 reports (in addition to the above statistics): - Maximum global error over range. The 'true' solution over the range is obtained by a reliable integrator at a more stringent tolerance. 3 reports (in addition to the above): - Maximum local error over range, ie. max over all meshpoints of LENRM = ||(computed y) - yloc||/ERRBND where yloc is the true local solution through the previous meshpoint, and ERRBND, the assumed error bound, is explained below. - Fraction of steps where LENRM exceeded 1. - Fraction of steps where LENRM exceeded 5. 4 reports (in addition to the above): - An analysis of the local error estimates used by SOLVER as the basis for its error control. Under development and described more fully in the actual code. NORMEF one of 0 1 or 2 , selects normalized efficiency statistics. These try to compensate for the fact that achieved accuracy may be much higher or lower than that requested by TOL, and this relationship is very problem- and method- dependent. For each problem, a least-squares fit is made of log10(actual error) vs log10(TOL) and used to estimate what the various cost statistics would be for an actual error of 10**n. This is achieved by interpolation, for those n such that 10**n lies within the range of accuracies achieved with the user-specified tolerances. 0 No normalized statistics 1 Normalized statistics are produced taking the 'actual error' used in the least squares fit to be the endpoint global error. 2 Normalized statistics are produced taking 'actual error' as the maximum global error over the range. N.B. In this case OPT must be at least 2. NRMTYP one of 1, 2 or 3, selects the norm used in assessing the size of local and global errors. It should be chosen by the user to agree with the norm used in SOLVER. We offer: 1 Max-norm. 2 2-norm (Euclidean norm). 3 r.m.s. norm, that is (2-norm of x)/sqrt(n) for an n-vector x. TOL (input) Real array, holds list of up to 10 tolerances to be used, in strictly decreasing order, with 0 as terminator. Each Problem is integrated at each tolerance in turn. Example: in calling program REAL TOL(11) DATA TOL/1E-1,1E-3,1E-5,1E-7,7*0E0/ requests the four tolerances .1, .001, .00001, .0000001. IDLIST (input) Integer array, holds list of groups of problems, and specifies for each one whether it is to be integrated in scaled or unscaled form (see General Notes above). Each problem is specified by a numeric code, 11 to 14 for problems A1 to A4, 21 to 25 for B1 to B5 etc. A zero terminates a group and two zeros terminate the list of groups. If the problem code is given a negative sign, the system is integrated in unscaled form; if a positive sign, in scaled form. Example: in calling program INTEGER IDLIST(7) DATA IDLIST/11,22,0,-31,-51,0,0/ specifies Group 1 consisting of Problems A1,B2 and Group 2 of Problems C1,E1. The first two are to be solved in the scaled form and the last two in unscaled form. The total length of the list including zeros must be at most 60 items. FLAG (output) Real. A nonzero value indicates that the call to STDTST was aborted because of argument errors, in which case the values of the decimal digits of FLAG indicate the error(s) that have occurred, as follows: 1: OPT invalid. 2: NORMEF invalid. 3: NORMEF = 2 was requested with OPT = 1. 4: A negative tolerance was supplied, or the list of tolerances was not in decreasing order. 5: The list of tolerances was empty or not terminated by a zero. 6: An invalid Problem-Id was found in IDLIST. 7: The list of groups in IDLIST is empty or is not terminated by two zeros or has more than the maximum allowed number (6) of groups. 8: NRMTYP invalid. Eg. a value FLAG = 0.245E 03 indicates that errors 2, 4 and 5 in the above list have occurred. Its value if nonzero is printed by STDTST anyway, but FLAG is meant to be inspected if further action of the main program depends on a successful call to STDTST. 3. Interface routine METHOD --------- ------- ------ This invokes the code being tested, call it SOLVER. The specification is SUBROUTINE METHOD(N,X,Y,XEND,TOL,HMAX,HSTART) INTEGER N DOUBLE PRECISION X,Y(N),XEND,TOL,HMAX,HSTART EXTERNAL FCN, PDERV METHOD is to be written by the user as a simple integrator to advance the solution of N differential equations from the initial values held in X,Y up to XEND, with an unweighted absolute error control of TOL. HMAX is a recommended maximum stepsize and HSTART is a recommended initial stepsize. If SOLVER can make use of these two parameters, the statistics will probably be more favorable and reliable, but their use is not crucial. The derivatives, and the analytical Jacobian matrix, of the problem are computed by package routines FCN and PDERV respectively. Thus certainly FCN, and in most cases PDERV, must be arguments to SOLVER, and they must be declared EXTERNAL in METHOD. METHOD should call SOLVER in one-step mode so that a call to the package routine STATS can be made after each successful step. If SOLVER does not have this facility, SOLVER must have a call to STATS inserted at the appropriate point in the code. Some calls to METHOD are intended to be aborted after a few integration steps by the STATS call setting X = XEND. Thus a test should be made after each call to STATS, of the form if STATS has set X = XEND then EXIT. NB: If the actual X argument to STATS is different from the X argument of METHOD (which may be necessary with some SOLVERs), ensure that the X argument of METHOD is set to XEND before exit, else the package will report 'METHOD failed to start'. The algorithm for METHOD should thus be of the form: - Declare all arguments and workspace expected by SOLVER - Set appropriate options including absolute error control and one-step mode - Initialize extra arguments if required - FOR each successful step DO - Call SOLVER( ... ,FCN,PDERV, ... ) EXIT if SOLVER is in trouble. - Set X,Y to the just computed meshpoint x and solution vector y - Set ERRBND to the bound that is satisfied by ||ERREST||, and hence is intended to be satisfied by ||LE||, at this step. - Set ERREST to the local error estimate vector (OPT=4 only) (See [4] for discussion and note that X,Y are ignored unless OPT.GE.2, ERRBND is ignored unless OPT.GE.3, and ERREST is ignored unless OPT.GE.4.) - Call STATS(X,Y,ERRBND,ERREST) - EXIT if X .ge. XEND. - ENDLOOP On normal exit X,Y must hold XEND and the solution at XEND. On exit because SOLVER was in trouble, X must hold the final point reached. On an exit forced by STATS, X must hold XEND. 4. Controlling the destination of output ----------- --- ----------- -- ------ The unit number on which the package writes its output is set by a call to one of the package routines, and you can find out what it is, by putting the statement IOUT = CONST(3) in your main program. Probably output will default to your terminal, which is good for debugging. For more serious work you may want to send output to a file. The statements IOUT = CONST(3) OPEN(IOUT, FILE=filename, other options.. ) will do this for you, assuming your Fortran I/O is consistent with the 1977 standard. 5. The routines FCN, PDERV --- -------- ---- ----- The specification of FCN is SUBROUTINE FCN(X,Y,YP) DOUBLE PRECISION X,Y(20),YP(20) On entry X holds the independent variable and Y holds the vector of dependent variables. On exit YP holds the vector of derivatives for the problem being solved (selected by a switch in COMMON). The specification for PDERV is SUBROUTINE PDERV(X,Y,DY) DOUBLE PRECISION X,Y(20),DY(400) where X and Y are as for FCN. The entries of the Jacobian matrix are stored in the first N**2 elements of DY with df(i)/dy(j) being stored in element i+(j-1)*N. Thus DY may be treated as if it were declared DIMENSION DY(N,N) 6. Function, Jacobian and LU Decomposition counts --------- -------- --- -- ------------- ------ These are maintained in three COMMON variables: COMMON/STCOM6/NFCN,NJAC,NLUD Each call to FCN, PDERV and DDCOMP increments NFCN, NJAC and NLUD by 1 respectively. If SOLVER uses its own linear algebra routines it is the user's responsibility to insert the above COMMON at an appropriate place in his code and set NLUD correctly. This may be done by incrementing it at each LU decomposition call, or by setting it equal to an independently maintained count before exit from METHOD. Similar comments apply to NJAC if SOLVER does its own Jacobian evaluation (eg. by numerical differencing). If a method does not use Jacobians, NJAC and NLUD may be used for gathering some other statistics. 7. Sample Program ------ ------- The following driver is the program used to generate the results of Fig. 4 of [4]. C SAMPLE DRIVER FOR STDTST, WITH ONE GROUP CONSISTING ONLY C OF PROBLEM E3 SOLVED IN SCALED FORM, AT FOUR TOLERANCES. C IN THIS CASE THE ARRAYS IDLIST, TOL NEED NOT BE SO LONG. C CHARACTER TITLE*80 INTEGER OPTION(10),IDLIST(60) REAL TOL(11) DATA TITLE/'SECDER, ADDISON-ENRIGHT SECOND DERIVATIVE METHOD'/ * , OPTION/2, 2, 1, 0, 6*0/ * , TOL/1E-2, 1E-4, 1E-6, 1E-8, 7*0E0/ * , IDLIST/53, 0, 58*0/ CALL STDTST(TITLE, OPTION, TOL, IDLIST, FLAG) STOP END C C SUBROUTINE METHOD(N,X,Y,XEND,TOL,HMAX,HSTART) C C DRIVER FOR THE SECDER CODE WHICH IS PART OF THE PACKAGE. C IT IS SOMEWHAT LENGTHY BECAUSE ITS INTERRUPT MECHANISM DOES C NOT ALLOW INTERRUPT IMMEDIATELY AFTER ACCEPTING A STEP. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION X,Y(N),XEND,TOL,HMAX,HSTART EXTERNAL FCN,PDERV DOUBLE PRECISION C(20),YP(20,11),W(400),PD(400),WK(20,12) INTEGER INF(40) C COMMON/STCOM6/NFCN,NJAC,NLUD C DATA NDIM/20/ C IND=2 DO 20 I=1,5 INF(I)=0 C(I)=0.D0 20 CONTINUE C C SET ABS ERROR CONTROL: INF(1); INTERRUPT NO. 2: INF(5); C MIN,MAX & STARTING STEPSIZE: C(2),C(4),C(5). INF(1)=1 INF(5)=1 C(2)=1D-12 C(4)=HMAX C(5)=HSTART C 50 CALL TRUE(FCN,PDERV,NDIM,N,X,Y,XEND,TOL,IND,C,INF,YP,W,PD,WK) IF(IND.EQ.6)GOTO 50 C WRITE(5,999)X,Y,C(13),(WK(I,1),I=1,N) C999 FORMAT(20X,10F10.6) IF(IND.NE.5) GOTO 60 TEMP=C(13) C C(13),WK(*,1) ARE THE ABOUT-TO-BE-ACCEPTED X,Y. C WK(*,12) IS THE ERROR-ESTIMATE VECTOR, DELIVERED C BY A SMALL CHANGE IN 'TRUE'. CALL STATS(C(13),WK(1,1),TOL,WK(1,12)) IF(C(13).NE.TEMP) GOTO 70 GOTO 50 C 60 IF(IND.NE.3) GOTO 70 X = XEND GOTO 80 C C FAILURE EXIT OF SOME KIND: 70 X=C(13) C WRITE(IOUT,110)IND,(INF(I),I=9,15) C110 FORMAT(1H ,'IND,INF(9)..INF(15)=',8I10) 80 CONTINUE NLUD=INF(15) RETURN END * * * * * References ----------- * [1] W H Enright, 'Using a testing package for the automatic assessment of numerical methods for ODEs', in Performance Evaluation of Numerical Software, (Fosdick, ed), IFIP, North Holland Publ Co (1979) 199-213. * * [2] W H Enright and T E Hull, 'Comparing numerical methods for the solution of stiff systems of ODEs arising in chemistry', in Numerical Methods for Differential Systems (Lapidus and Schiesser, eds), Academic Press, New York (1976) 45-65. * [3] W H Enright, T E Hull and B Lindberg, 'Comparing numerical methods for stiff systems of ordinary differential equations', BIT 15(1975) 10-48. * [4] W H Enright and J D Pryce, 'A pair of packages for assessing initial value methods', University of Toronto Technical Report no. 167/83. * [5] L F Shampine 'Evaluation of a test set for stiff ODE solvers', TOMS 7(1981)409-420. * * * REAL FUNCTION CONST(I) C C********+*********+*********+*********+*********+*********+*********+** C .. Scalar Arguments .. INTEGER I C .. Local Scalars .. CHARACTER*32 MCNAME C .. Local Arrays .. REAL C(4) C .. Intrinsic Functions .. INTRINSIC ICHAR C .. Data statements .. C C CONST AND CLOCK ENCLOSE (WE HOPE) ALL THE MACHINE-DEPENDENT PARTS C OF THE STIFF AND NONSTIFF DETEST PACKAGES, EXCEPT THE TIMING C DATA IN THE IVALU ROUTINE OF EACH PACKAGE. C C CALLS WITH VALUES I=1 TO 4 RETURN THE FOLLOWING VALUES IN 'CONST': C I=1 UNIT ROUNDOFF APPROXIMATELY, IN THE PRECISION USED BY THE C ODE-SOLVING PART OF THE PACKAGE. C I=2 NUMBER NEAR UNDERFLOW THRESHOLD C I=3 STANDARD OUTPUT UNIT NUMBER C I=4 VALUE OF TSTTIM (USED IN CNTROL) C C CALLS WITH VALUES -1 TO -32 RETURN THE 'ICHAR' VALUE OF SUCCESSIVE C CHARACTERS OF THE NAME OF THE COMPUTER WE ARE RUNNING ON. (CLUMSY BUT C INTENDED TO ISOLATE MACHINE-DEPENDENCIES HERE) C C A CALL WITH I OUTSIDE THESE RANGES (DONE WITH I=0 NEAR THE HEAD OF TH C MAIN NSDTST & STDTST ROUTINES) RETURNS CONST=0 AND C 1. IS TO BE USED FOR MACHINE-DEPENDENT INITIALIZATIONS SUCH AS THE C SUPPRESSION OF UNDERFLOW MESSAGES. C C****** VALUES FOR IBM3033 MODEL N12 ****** CIBM DATA C/2.25E-16,1E-50,6.0,0.5/ C****** VALUES FOR DEC10 MODEL KL10 ****** CDEC DATA C /2.17E-19,1E-38,5.0,0.5/ DATA C/4*0.0/, MCNAME/ * '..PUT NAME OF COMPUTER HERE..'/ C .. Executable Statements .. C IF (I.GE.1 .AND. I.LE.4) THEN CONST = C(I) ELSE IF (I.LT.0) THEN CONST = ICHAR(MCNAME(-I:-I)) ELSE C SUPPRESS UNDERFLOW REPORTING (DEC): CDEC CALL ERRSET(0,6) C SUPPRESS UNDERFLOW REPORTING (IBM): CIBM CALL ERRSET(208,256,-1,0,0,208) C CONST = 0 END IF RETURN END C C********+*********+*********+*********+*********+*********+*********+** C REAL FUNCTION CLOCK(S) C .. Scalar Arguments .. REAL S C .. Executable Statements .. C C********+*********+*********+*********+*********+*********+*********+** C CLOCK IS 'RESET' TO 0 BY A CALL C S = CLOCK(0.0) C AND THEN (WITH S SET AS ABOVE) DELIVERS THE ELAPSED CPU SECONDS C SINCE LAST RESET, BY CALLS OF FORM C TIME = CLOCK(S) C C THIS WORKS ON AN IBM: CIBM CLOCK = UTTIMR(1) - S C THIS WORKS ON A DEC10 UNDER TOPS10 IN CONJUNCTION WITH THE C COMMAND "SET TIME LLL" AT MONITOR LEVEL WHERE LLL IS A C SUITABLE TIME LIMIT. CDEC CLOCK = -TIM2GO(S) -S C WHEN FIRST MOUNTING THE PACKAGE, LET IT RETURN ZERO AS BELOW: CLOCK = 0.0 RETURN END REAL FUNCTION CONST(I) C C********+*********+*********+*********+*********+*********+*********+** C .. Scalar Arguments .. INTEGER I C .. Local Scalars .. CHARACTER*32 MCNAME C .. Local Arrays .. REAL C(4) C .. Intrinsic Functions .. INTRINSIC ICHAR C .. Data statements .. C C CONST AND CLOCK ENCLOSE (WE HOPE) ALL THE MACHINE-DEPENDENT PARTS C OF THE STIFF AND NONSTIFF DETEST PACKAGES, EXCEPT THE TIMING C DATA IN THE IVALU ROUTINE OF EACH PACKAGE. C C CALLS WITH VALUES I=1 TO 4 RETURN THE FOLLOWING VALUES IN 'CONST': C I=1 UNIT ROUNDOFF APPROXIMATELY, IN THE PRECISION USED BY THE C ODE-SOLVING PART OF THE PACKAGE. C I=2 NUMBER NEAR UNDERFLOW THRESHOLD C I=3 STANDARD OUTPUT UNIT NUMBER C I=4 VALUE OF TSTTIM (USED IN CNTROL) C C CALLS WITH VALUES -1 TO -32 RETURN THE 'ICHAR' VALUE OF SUCCESSIVE C CHARACTERS OF THE NAME OF THE COMPUTER WE ARE RUNNING ON. (CLUMSY BUT C INTENDED TO ISOLATE MACHINE-DEPENDENCIES HERE) C C A CALL WITH I OUTSIDE THESE RANGES (DONE WITH I=0 NEAR THE HEAD OF TH C MAIN NSDTST & STDTST ROUTINES) RETURNS CONST=0 AND C 1. IS TO BE USED FOR MACHINE-DEPENDENT INITIALIZATIONS SUCH AS THE C SUPPRESSION OF UNDERFLOW MESSAGES. C C****** VALUES FOR IBM3033 MODEL N12 ****** CIBM DATA C/2.25E-16,1E-50,6.0,0.5/ C****** VALUES FOR DEC10 MODEL KL10 ****** CDEC DATA C /2.17E-19,1E-38,5.0,0.5/ DATA C/4*0.0/, MCNAME/ * '..PUT NAME OF COMPUTER HERE..'/ C .. Executable Statements .. C IF (I.GE.1 .AND. I.LE.4) THEN CONST = C(I) ELSE IF (I.LT.0) THEN CONST = ICHAR(MCNAME(-I:-I)) ELSE C SUPPRESS UNDERFLOW REPORTING (DEC): CDEC CALL ERRSET(0,6) C SUPPRESS UNDERFLOW REPORTING (IBM): CIBM CALL ERRSET(208,256,-1,0,0,208) C CONST = 0 END IF RETURN END C C********+*********+*********+*********+*********+*********+*********+** C REAL FUNCTION CLOCK(S) C .. Scalar Arguments .. REAL S C .. Executable Statements .. C C********+*********+*********+*********+*********+*********+*********+** C CLOCK IS 'RESET' TO 0 BY A CALL C S = CLOCK(0.0) C AND THEN (WITH S SET AS ABOVE) DELIVERS THE ELAPSED CPU SECONDS C SINCE LAST RESET, BY CALLS OF FORM C TIME = CLOCK(S) C C THIS WORKS ON AN IBM: CIBM CLOCK = UTTIMR(1) - S C THIS WORKS ON A DEC10 UNDER TOPS10 IN CONJUNCTION WITH THE C COMMAND "SET TIME LLL" AT MONITOR LEVEL WHERE LLL IS A C SUITABLE TIME LIMIT. CDEC CLOCK = -TIM2GO(S) -S C WHEN FIRST MOUNTING THE PACKAGE, LET IT RETURN ZERO AS BELOW: CLOCK = 0.0 RETURN END * * * Nonstiff DETEST 1986 version ----- ------ ---- ------- by W H Enright, and J D Pryce, Dept of Computer Science, School of Mathematics University of Toronto, University Walk Toronto M5S 1A4 Bristol BS8 1TW Canada England Tel (416) 978-6025 Tel (272) 303335 * Please inform the authors of any errors in code or documentation. * 1. General Notes ------- ----- * Nonstiff DETEST is a package to test the performance of initial-value codes for nonstiff differential systems. This code is a revision of the 1971 version, used to produce the results reported on in [2,4]. * A set of test problems, described in detail in [2], is incorporated in the package. The code being tested is run on a selection of these problems at various tolerances. The user selects the problems and the tolerances, and also organizes the problems into groups for statistical reporting purposes, at his discretion. * To test a code a user must write an interface routine called METHOD, described below, and then call NSDTST with the desired options. Note that NSDTST comes in a 'single' and a 'double' precision version for use according as the software under test is written in single or double precision. The arguments of NSDTST are single precision but METHOD must be implemented in the appropriate precision. * The package divides naturally into four parts: * NSDTST,CNTROL and various service routines organize the assembling, computation and reporting of statistics. * STATS is the routine which 'instruments' the code being tested and passes statistics via COMMON to CNTROL and NSDTST. * FCN, IVALU, EVALU describe the set of test problems. FCN gives the r.h.s. f(x,y) of the ODE system. IVALU gives the initial conditions, scaling weights and other data about each problem. EVALU gives accurately computed values at the endpoint. * TRUE and its subordinate routines (alias the Hull-Enright-Jackson code DVERK based on Verner's Runge-Kutta formulas) form a reliable nonstiff solver for computing the 'true' global and local solutions when required. * There is also a 'dummy' NSDTST and STATS to help the user debug his METHOD routine (described below); a utility NSGTIM which can be used on each new machine to generate timing data embedded in the code; and a utility NSGWT can be used if ever a user wishes to add further test problems to the set. * Main Lines of Calling Hierarchy (user-supplied routines are in boxes) * * * +--------+ | User's |---NSDTST---CNTROL-----IVALU |Program | | +--------+ +--------+ | +------+ |'SOLVER'| |---|METHOD|----|(Code |->-+ | +------+ | being | | | | | tested)| | | | +--------+ |---FCN | | | | STATS---TRUE--->--+ | +----EVALU * We acknowledge valuable recommendations in Shampine's paper [5]. In particular the package will, by default, integrate each system in scaled form, scaling each solution component by its maximum observed value over the range of integration. That is, the change of variable -1 z = D y is done where D = diag(w(1), .., w(n)) * and w(i) =max |i-th component of y| over the range. The problem -1 solved is then z' = D f(x,Dz). The weights w(i) were found by an accurate integration of each problem and are embedded in IVALU. Note that this scaling affects the norms which are used in measuring all errors, and thus can have a considerable effect on the accuracy in some of the problems. * If the problem code in IDLIST (see below) is given a negative sign the system is solved in its 'natural' scaling, as was done in the 1975 version of DETEST. * * References ----------- * [1] W H Enright, 'Using a testing package for the automatic assessment of numerical methods for ODEs', in Performance Evaluation of Numerical Software, (Fosdick, ed), IFIP, North Holland Publ Co (1979) 199-213. * [2] T E Hull, W H Enright, B M Fellen and A E Sedgwick, 'Comparing numerical methods for ordinary differential equations', SIAM J. Numer. Anal. 9(1972)603-637. * [3] W H Enright and J D Pryce, 'A pair of packages for assessing initial value methods', University of Toronto Technical Report no. 167/83. * [4] W H Enright and T E Hull, 'Test results on initial value methods for nonstiff ordinary differential equations', SIAM J. Numer. Anal. 13(1976)944-961. * [5] L F Shampine 'Evaluation of a test set for stiff ODE solvers', TOMS 7(1981)409-420. * * * * * * * * 2. Arguments to NSDTST: --------- -- ------- * TITLE (input) Character of length 80, holds name of method being tested. * OPTION (input) Integer array of length 10, only elements 1 to 3 are used and are referred to henceforth as OPT, NORMEF and NRMTYP. (OPTION(4) is also used when OPT=4) * OPT one of 1, 2, 3 or 4. OPT selects level of analysis required: 1 gives a report of the following at each tolerance used: - Total time per integration - Overhead time excluding function calls. - Number of function calls and successful steps over range. - Global error at endpoint XEND, divided by TOL, ie. ||(computed y) - (true y)||/TOL at x=XEND The norm used throughout the package is that chosen by NRMTYP. * 2 reports (in addition to the above statistics): - Maximum global error over range. The 'true' solution over the range is obtained by a reliable integrator at a more stringent tolerance. * 3 reports (in addition to the above): - Maximum local error over range, ie. max over all meshpoints of LENRM = ||(computed y) - yloc||/ERRBND where yloc is the true local solution through the previous meshpoint, and ERRBND, the assumed error bound, is explained below. - Fraction of steps where LENRM exceeded 1. - Fraction of steps where LENRM exceeded 5. * 4 reports (in addition to the above): - An analysis of the local error estimates used by SOLVER as the basis for its error control. At this level three assumptions are made. First, that at each step SOLVER forms two approximations, y and y*, to the local solution yloc at the new meshpoint, such that asymptotically as TOL->0, y* is 'more accurate' than y. Second, that the approximation which is taken as the computed solution at the new meshpoint is either always y* (in which case one says local extrapolation is used) or always y (in which case it is not used). The vector LE = y - yloc is the true local error in the 'less accurate' solution y, and ERREST = y - y* is an estimate of LE. It is assumed finally that the error control consists in keeping ||ERREST||, in an appropriate norm, below ERRBND at each step. * Note that some methods, such as Merson's method, cannot be regarded in this way. * At this level DETEST analyses how accurately ERREST approximates to LE, by forming a scatter plot of the values of r1 = ||ERREST - LE||/ERRBND (vertical axis) against r2 = ||ERREST||/ERRBND (horizontal) at each step. Note ERREST - LE = -(y* - yloc) = -LE*, say, so that LENRM defined above is r1 if local extrapolation is being done. For an 'ideal' error control strategy, we expect the plotted points to cluster near (1,0) on the graph, whether or not local extrapolation is used. * To use this level of analysis the user must: a) Ensure that the STATS call in METHOD delivers ERREST as defined above (with the correct sign!). b) Set OPTION(4) as follows. =0 Argument Y to STATS is y above (no local extrapolation). =1 Y is y* above (local extrapolation). * For each integration, a scatter plot is produced. Each of the ratios r1, r2 is put into one of 12 class-intervals -7 -7 -6 2 3 3 0<=r<2 , 2 <=r<2 , ..., 2 <=r<2 , 2 <=r= 2, and have a possibly more efficient code to put in its place. NSTL is relevant if you are interested in the algorithms used by the package, specifically the step-lumping process which takes place in STATS at stringent tolerances. * * * * 8. Subroutines in the Package ----------- -- --- ------- * In order of appearance in the files. The list also shows, for each routine, the other package routines and COMMON areas which it uses. A name in parentheses, like (FCN) denotes a routine which is called at one remove (eg. METHOD calls SOLVER which must call FCN) or which is passed as an argument rather than being an external reference (eg. FCN in TRUE). * In CONCLK file CONST calls: none CLOCK calls: none * In NSDTST file NSDTST calls: PARCHK LSQFIT RATIO EFSTAT CNTROL CONST ; NSCOM1 NSCOM3 PARCHK calls: none LSQFIT calls: none RATIO calls: none EFSTAT calls: none CNTROL calls: DIFNRM STATS CONST CLOCK IVALU EVALU METHOD PLOT ; NSCOM1 NSCOM2 NSCOM3 NSCOM5 NSCOM6 DIFNRM calls: none STATS calls: DIFNRM CONST TRUE FCN PLOT ; NSCOM1 NSCOM2 NSCOM3 NSCOM4 NSCOM6 PLOT calls: none * In NSTRUE file TRUE calls: CONST (FCN2 ) FCN2 calls: FCN * In NSPROB file IVALU calls: none EVALU calls: none FCN calls: ; NSCOM5 NSCOM6 * User-supplied METHOD calls: STATS (FCN ) * * 9. Definition of Common Areas and Dictionary of Data-flow ---------- -- ------ ----- --- ---------- -- --------- * The flow of information between those routines which use COMMON is indicated for each variable by the codes S: the variable is assigned a value (Set) in this routine, possibly by a call to another routine to which the variable is passed as an argument. A: the value is used (Accessed) in this routine. * For counters and similar variables, these codes are used instead of code S: I: the variable is Initialized in this routine. U: the variable is Updated in this routine. * * COMMON /NSCOM1/ passes information from NSDTST to CNTROL and STATS. * NSDTST | CNTROL | | STATS | | | DIFNRM | | | | S A A - ERRTOL DOUBLE. Copy of current error tolerance. S A A - OPT INTEGER. Copy of OPTION(1) argument of NSDTST. S - - A NRMTYP INTEGER. Copy of OPTION(3) argument of NSDTST. S - A - XTRAP INTEGER. Copy of OPTION(4) argument of NSDTST. S A - - ID INTEGER. Internal code of current problem, 1 for A1, ..., 13 for B3, etc. S A - - IWT INTEGER. Flag for scaling (+1: Scaled. -1: Unscaled) S - - - IOUT INTEGER. Standard output unit number. * * * * COMMON /NSCOM2/ communicates between CNTROL and STATS. * CNTROL | STATS | | S A XEND DOUBLE. End of integration range of current problem. A S HSTART DOUBLE. Initial stepsize passed to METHOD for integration proper. S A N INTEGER. No. of equations in current problem. S A IFLAG INTEGER. Set by CNTROL to inform STATS what it is to do: =0 METHOD is being timed. =1 Initializing call of STATS from CNTROL to set up NSCOM4. =2 Preliminary integration to determine HSTART, aborted after 2 steps. =3 Integration proper, compiling statistics. * * A SA INDL,INDG Error flags for the local and global 'true solutions' obtained by calls to routine TRUE. * * * * * COMMON /NSCOM3/ outputs statistics from CNTROL and STATS. * NSDTST | CNTROL | | STATS | | | A S - XFIN DOUBLE. Point of failure of METHOD if it doesn't reach XEND. A - S XTRUE DOUBLE. Point of failure of TRUE if any. If both local and global fail, point of global failure is returned. A S - TIME REAL. CPU time for one integration as measured by CLOCK function. A S - OVHD REAL. Equals TIME less estimated cost of FCN calls. A I U TRUTIM REAL. The time spent in calls to TRUE. Not relevant to performance of METHOD but measures the overhead incurred by the testing package when OPT = 2, 3 or 4. Not printed but available. A S - GEND REAL. Norm of global error of METHOD at XEND. * * A I U GEMX REAL. Maximum of global error over all lumped step meshpoints, ie. usually over all meshpoints of METHOD, except when ERRTOL is very small. A I U LEMXSC REAL. Maximum local error in units of ERRBND, over all lumped step meshpoints. A S - NFCN INTEGER. Copy of NFCN1, see /NSCOM6/. /NSCOM6/ A I U NSTP INTEGER. Counts (unlumped) steps taken by METHOD in current integration. - I U NSTL INTEGER. Counts lumped steps formed in current integration (see STATS). Not printed but available. A I U NDCV,NBAD INTEGER. Count lumped steps on which SOLVER's local error control was deceived, resp. badly deceived. A I U NTRU INTEGER. Counts lumped steps on which true local solution was successfully computed, hence valid local error statistics obtained. Used in computing 'fraction deceived' information. Reported if different from NSTP. Note NTRU <= NSTL <= NSTP. - S - NSTART INTEGER. No. of FCN calls needed by METHOD to start, ie. to do preliminary integration (2 steps). Not printed out but available. * * COMMON /NSCOM4/ is used only by STATS, to preserve information from one call of STATS to another. All variables are set and/or updated in STATS. * XOLD1 DOUBLE. Similar to XOLD but used in preliminary integration. XOLD,YOLD DOUBLE and DOUBLE array. Copy of METHOD's computed solution at end of previous lumped step. Used as actual arguments of TRUE local solution call. XOLDG,YOLDG DOUBLE and DOUBLE array. Hold 'true' global solution updated to end of previous lumped step. Used as actual arguments of TRUE global solution call. CG,PDG,WKG,WG,YPG,INFG Workspace for 'true' global solution. XT DOUBLE. Holds last METHOD meshpoint between calls to STATS. PRECIS DOUBLE. Holds 1000 * (unit roundoff) approx. ERLUMP DOUBLE. Accumulates METHOD's local error estimates to form an estimate over a lumped step. * * COMMON /NSCOM5/ passes information between CNTROL and FCN, (or any replacement a user may provide for FCN). * CNTROL | FCN | | * S A WT DOUBLE. Array of weights used to implement the 'scaled' integration option. S A IWT1,N1,ID1 INTEGER. Copies of IWT,N,ID in /NSCOM1/ or /NSCOM2/. * * COMMON /NSCOM6/ holds a counter. It is initialized in CNTROL, saved-and-restored in STATS, and eventually copied by CNTROL to the corresponding variable in /NSCOM3/. * CNTROL | STATS | | FCN | | | * IA AS U - - NFCN1 INTEGER. Counts calls to FCN. * * There is also a COMMON/NSCOM7/ used by the dummy (debugging) versions of NSDTST and STATS for communication. * SUBROUTINE NSDTST(TITLE,OPTION,TOL,IDLIST,FLAG) C C********+*********+*********+*********+*********+*********+*********+** C G E N E R A L D O C U M E N T A T I O N C--------+---------+---------+---------+---------+---------+---------+-- C C C C NONSTIFF DETEST 1986 VERSION C ----- ------ ---- ------- C BY W H ENRIGHT, AND J D PRYCE, C DEPT OF COMPUTER SCIENCE, SCHOOL OF MATHEMATICS C UNIVERSITY OF TORONTO, UNIVERSITY WALK C TORONTO M5S 1A4 BRISTOL BS8 1TW C CANADA ENGLAND C TEL (416) 978-6025 TEL (272) 303335 C C PLEASE INFORM THE AUTHORS OF ANY ERRORS IN CODE OR C DOCUMENTATION. C C 1. GENERAL NOTES C ------- ----- C C NONSTIFF DETEST IS A PACKAGE TO TEST THE PERFORMANCE OF INITIAL-VALUE C CODES FOR NONSTIFF DIFFERENTIAL SYSTEMS. THIS CODE IS A REVISION OF C THE 1971 VERSION, USED TO PRODUCE THE RESULTS REPORTED ON IN [2,4]. C C A SET OF TEST PROBLEMS, DESCRIBED IN DETAIL IN [2], IS INCORPORATED C IN THE PACKAGE. THE CODE BEING TESTED IS RUN ON A SELECTION OF THESE C PROBLEMS AT VARIOUS TOLERANCES. THE USER SELECTS THE PROBLEMS AND C THE TOLERANCES, AND ALSO ORGANIZES THE PROBLEMS INTO GROUPS FOR C STATISTICAL REPORTING PURPOSES, AT HIS DISCRETION. C C TO TEST A CODE A USER MUST WRITE AN INTERFACE ROUTINE CALLED METHOD, C DESCRIBED BELOW, AND THEN CALL NSDTST WITH THE DESIRED OPTIONS. NOTE C THAT NSDTST COMES IN A 'SINGLE' AND A 'DOUBLE' PRECISION VERSION FOR C USE ACCORDING AS THE SOFTWARE UNDER TEST IS WRITTEN IN SINGLE OR C DOUBLE PRECISION. THE ARGUMENTS OF NSDTST ARE SINGLE PRECISION BUT C METHOD MUST BE IMPLEMENTED IN THE APPROPRIATE PRECISION. C C THE PACKAGE DIVIDES NATURALLY INTO FOUR PARTS: C C NSDTST,CNTROL AND VARIOUS SERVICE ROUTINES C ORGANIZE THE ASSEMBLING, COMPUTATION AND REPORTING OF C STATISTICS. C C STATS C IS THE ROUTINE WHICH 'INSTRUMENTS' THE CODE BEING TESTED AND C PASSES STATISTICS VIA COMMON TO CNTROL AND NSDTST. C C FCN, IVALU, EVALU C DESCRIBE THE SET OF TEST PROBLEMS. FCN GIVES THE R.H.S. C F(X,Y) OF THE ODE SYSTEM. IVALU GIVES THE INITIAL CONDITIONS, C SCALING WEIGHTS AND OTHER DATA ABOUT EACH PROBLEM. EVALU C GIVES ACCURATELY COMPUTED VALUES AT THE ENDPOINT. C C TRUE AND ITS SUBORDINATE ROUTINES C (ALIAS THE HULL-ENRIGHT-JACKSON CODE DVERK BASED ON VERNER'S C RUNGE-KUTTA FORMULAS) FORM A RELIABLE NONSTIFF SOLVER FOR C COMPUTING THE 'TRUE' GLOBAL AND LOCAL SOLUTIONS WHEN REQUIRED. C C THERE IS ALSO A 'DUMMY' NSDTST AND STATS TO HELP THE USER DEBUG HIS C METHOD ROUTINE (DESCRIBED BELOW); A UTILITY NSGTIM WHICH CAN BE USED C ON EACH NEW MACHINE TO GENERATE TIMING DATA EMBEDDED IN THE CODE; AND C A UTILITY NSGWT CAN BE USED IF EVER A USER WISHES TO ADD FURTHER TEST C PROBLEMS TO THE SET. C C MAIN LINES OF CALLING HIERARCHY (USER-SUPPLIED ROUTINES ARE IN BOXES) C C C C +--------+ C | USER'S |---NSDTST---CNTROL-----IVALU C |PROGRAM | | +--------+ C +--------+ | +------+ |'SOLVER'| C |---|METHOD|----|(CODE |->-+ C | +------+ | BEING | | C | | | TESTED)| | C | | +--------+ |---FCN C | | | C | STATS---TRUE--->--+ C | C +----EVALU C C WE ACKNOWLEDGE VALUABLE RECOMMENDATIONS IN SHAMPINE'S PAPER [5]. IN C PARTICULAR THE PACKAGE WILL, BY DEFAULT, INTEGRATE EACH SYSTEM IN C SCALED FORM, SCALING EACH SOLUTION COMPONENT BY ITS MAXIMUM OBSERVED C VALUE OVER THE RANGE OF INTEGRATION. THAT IS, THE CHANGE OF VARIABLE C -1 C Z = D Y IS DONE WHERE C D = DIAG(W(1), .., W(N)) C C AND W(I) =MAX |I-TH COMPONENT OF Y| OVER THE RANGE. THE PROBLEM C -1 C SOLVED IS THEN Z' = D F(X,DZ). THE WEIGHTS W(I) WERE FOUND BY AN C ACCURATE INTEGRATION OF EACH PROBLEM AND ARE EMBEDDED IN IVALU. C NOTE THAT THIS SCALING AFFECTS THE NORMS WHICH ARE USED IN C MEASURING ALL ERRORS, AND THUS CAN HAVE A CONSIDERABLE EFFECT ON THE C ACCURACY IN SOME OF THE PROBLEMS. C C IF THE PROBLEM CODE IN IDLIST (SEE BELOW) IS GIVEN A NEGATIVE SIGN THE C SYSTEM IS SOLVED IN ITS 'NATURAL' SCALING, AS WAS DONE IN THE 1975 C VERSION OF DETEST. C C C REFERENCES C ----------- C C [1] W H ENRIGHT, 'USING A TESTING PACKAGE FOR THE AUTOMATIC C ASSESSMENT OF NUMERICAL METHODS FOR ODES', IN PERFORMANCE C EVALUATION OF NUMERICAL SOFTWARE, (FOSDICK, ED), IFIP, NORTH C HOLLAND PUBL CO (1979) 199-213. C C [2] T E HULL, W H ENRIGHT, B M FELLEN AND A E SEDGWICK, 'COMPARING C NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS', SIAM J. C NUMER. ANAL. 9(1972)603-637. C C [3] W H ENRIGHT AND J D PRYCE, 'A PAIR OF PACKAGES FOR ASSESSING C INITIAL VALUE METHODS', UNIVERSITY OF TORONTO TECHNICAL REPORT C NO. 167/83. C C [4] W H ENRIGHT AND T E HULL, 'TEST RESULTS ON INITIAL VALUE METHODS C FOR NONSTIFF ORDINARY DIFFERENTIAL EQUATIONS', SIAM J. NUMER. C ANAL. 13(1976)944-961. C C [5] L F SHAMPINE 'EVALUATION OF A TEST SET FOR STIFF ODE SOLVERS', C TOMS 7(1981)409-420. C C C C C C C C C 2. ARGUMENTS TO NSDTST: C --------- -- ------- C C TITLE (INPUT) CHARACTER OF LENGTH 80, HOLDS NAME OF METHOD BEING C TESTED. C C OPTION (INPUT) INTEGER ARRAY OF LENGTH 10, ONLY ELEMENTS 1 TO 3 ARE C USED AND ARE REFERRED TO HENCEFORTH AS OPT, NORMEF AND NRMTYP. C (OPTION(4) IS ALSO USED WHEN OPT=4) C C OPT ONE OF 1, 2, 3 OR 4. OPT SELECTS LEVEL OF ANALYSIS REQUIRED: C 1 GIVES A REPORT OF THE FOLLOWING AT EACH TOLERANCE USED: C - TOTAL TIME PER INTEGRATION C - OVERHEAD TIME EXCLUDING FUNCTION CALLS. C - NUMBER OF FUNCTION CALLS AND SUCCESSFUL STEPS OVER RANGE. C - GLOBAL ERROR AT ENDPOINT XEND, DIVIDED BY TOL, IE. C ||(COMPUTED Y) - (TRUE Y)||/TOL AT X=XEND C THE NORM USED THROUGHOUT THE PACKAGE IS THAT CHOSEN BY NRMTYP. C C 2 REPORTS (IN ADDITION TO THE ABOVE STATISTICS): C - MAXIMUM GLOBAL ERROR OVER RANGE. THE 'TRUE' SOLUTION OVER C THE RANGE IS OBTAINED BY A RELIABLE INTEGRATOR AT A MORE C STRINGENT TOLERANCE. C C 3 REPORTS (IN ADDITION TO THE ABOVE): C - MAXIMUM LOCAL ERROR OVER RANGE, IE. MAX OVER ALL MESHPOINTS C OF C LENRM = ||(COMPUTED Y) - YLOC||/ERRBND C WHERE YLOC IS THE TRUE LOCAL SOLUTION THROUGH THE PREVIOUS C MESHPOINT, AND ERRBND, THE ASSUMED ERROR BOUND, IS EXPLAINED C BELOW. C - FRACTION OF STEPS WHERE LENRM EXCEEDED 1. C - FRACTION OF STEPS WHERE LENRM EXCEEDED 5. C C 4 REPORTS (IN ADDITION TO THE ABOVE): C - AN ANALYSIS OF THE LOCAL ERROR ESTIMATES USED BY SOLVER AS THE C BASIS FOR ITS ERROR CONTROL. AT THIS LEVEL THREE ASSUMPTIONS C ARE MADE. FIRST, THAT AT EACH STEP SOLVER FORMS TWO C APPROXIMATIONS, Y AND Y*, TO THE LOCAL SOLUTION YLOC AT THE C NEW MESHPOINT, SUCH THAT ASYMPTOTICALLY AS TOL->0, Y* IS 'MORE C ACCURATE' THAN Y. SECOND, THAT THE APPROXIMATION WHICH IS C TAKEN AS THE COMPUTED SOLUTION AT THE NEW MESHPOINT IS EITHER C ALWAYS Y* (IN WHICH CASE ONE SAYS LOCAL EXTRAPOLATION IS USED) C OR ALWAYS Y (IN WHICH CASE IT IS NOT USED). THE VECTOR C LE = Y - YLOC C IS THE TRUE LOCAL ERROR IN THE 'LESS ACCURATE' SOLUTION Y, C AND C ERREST = Y - Y* C IS AN ESTIMATE OF LE. IT IS ASSUMED FINALLY THAT THE ERROR C CONTROL CONSISTS IN KEEPING ||ERREST||, IN AN APPROPRIATE C NORM, BELOW ERRBND AT EACH STEP. C C NOTE THAT SOME METHODS, SUCH AS MERSON'S METHOD, CANNOT BE C REGARDED IN THIS WAY. C C AT THIS LEVEL DETEST ANALYSES HOW ACCURATELY ERREST C APPROXIMATES TO LE, BY FORMING A SCATTER PLOT OF THE VALUES OF C R1 = ||ERREST - LE||/ERRBND (VERTICAL AXIS) AGAINST R2 = C ||ERREST||/ERRBND (HORIZONTAL) AT EACH STEP. NOTE ERREST - C LE = -(Y* - YLOC) = -LE*, SAY, SO THAT LENRM DEFINED ABOVE IS C R1 IF LOCAL EXTRAPOLATION IS BEING DONE. FOR AN 'IDEAL' ERROR C CONTROL STRATEGY, WE EXPECT THE PLOTTED POINTS TO CLUSTER NEAR C (1,0) ON THE GRAPH, WHETHER OR NOT LOCAL EXTRAPOLATION IS C USED. C C TO USE THIS LEVEL OF ANALYSIS THE USER MUST: C A) ENSURE THAT THE STATS CALL IN METHOD DELIVERS ERREST AS C DEFINED ABOVE (WITH THE CORRECT SIGN!). C B) SET OPTION(4) AS FOLLOWS. C =0 ARGUMENT Y TO STATS IS Y ABOVE (NO LOCAL EXTRAPOLATION). C =1 Y IS Y* ABOVE (LOCAL EXTRAPOLATION). C C FOR EACH INTEGRATION, A SCATTER PLOT IS PRODUCED. EACH OF THE C RATIOS R1, R2 IS PUT INTO ONE OF 12 CLASS-INTERVALS C -7 -7 -6 2 3 3 C 0<=R<2 , 2 <=R<2 , ..., 2 <=R<2 , 2 <=R= 2, AND HAVE A POSSIBLY MORE EFFICIENT C CODE TO PUT IN ITS PLACE. NSTL IS RELEVANT IF YOU ARE C INTERESTED IN THE ALGORITHMS USED BY THE PACKAGE, SPECIFICALLY THE C STEP-LUMPING PROCESS WHICH TAKES PLACE IN STATS AT STRINGENT C TOLERANCES. C C C C C 8. SUBROUTINES IN THE PACKAGE C ----------- -- --- ------- C C IN ORDER OF APPEARANCE IN THE FILES. THE LIST ALSO SHOWS, FOR EACH C ROUTINE, THE OTHER PACKAGE ROUTINES AND COMMON AREAS WHICH IT USES. A C NAME IN PARENTHESES, LIKE (FCN) DENOTES A ROUTINE WHICH IS CALLED AT C ONE REMOVE (EG. METHOD CALLS SOLVER WHICH MUST CALL FCN) OR WHICH IS C PASSED AS AN ARGUMENT RATHER THAN BEING AN EXTERNAL REFERENCE (EG. C FCN IN TRUE). C C IN CONCLK FILE C CONST CALLS: NONE C CLOCK CALLS: NONE C C IN NSDTST FILE C NSDTST CALLS: PARCHK LSQFIT RATIO EFSTAT CNTROL CONST ; NSCOM1 C NSCOM3 C PARCHK CALLS: NONE C LSQFIT CALLS: NONE C RATIO CALLS: NONE C EFSTAT CALLS: NONE C CNTROL CALLS: DIFNRM STATS CONST CLOCK IVALU EVALU METHOD PLOT ; C NSCOM1 NSCOM2 NSCOM3 NSCOM5 NSCOM6 C DIFNRM CALLS: NONE C STATS CALLS: DIFNRM CONST TRUE FCN PLOT ; NSCOM1 NSCOM2 NSCOM3 C NSCOM4 NSCOM6 C PLOT CALLS: NONE C C IN NSTRUE FILE C TRUE CALLS: CONST (FCN2 ) C FCN2 CALLS: FCN C C IN NSPROB FILE C IVALU CALLS: NONE C EVALU CALLS: NONE C FCN CALLS: ; NSCOM5 NSCOM6 C C USER-SUPPLIED C METHOD CALLS: STATS (FCN ) C C C 9. DEFINITION OF COMMON AREAS AND DICTIONARY OF DATA-FLOW C ---------- -- ------ ----- --- ---------- -- --------- C C THE FLOW OF INFORMATION BETWEEN THOSE ROUTINES WHICH USE COMMON IS C INDICATED FOR EACH VARIABLE BY THE CODES C S: THE VARIABLE IS ASSIGNED A VALUE (SET) IN THIS ROUTINE, POSSIBLY C BY A CALL TO ANOTHER ROUTINE TO WHICH THE VARIABLE IS PASSED AS C AN ARGUMENT. C A: THE VALUE IS USED (ACCESSED) IN THIS ROUTINE. C C FOR COUNTERS AND SIMILAR VARIABLES, THESE CODES ARE USED INSTEAD OF C CODE S: C I: THE VARIABLE IS INITIALIZED IN THIS ROUTINE. C U: THE VARIABLE IS UPDATED IN THIS ROUTINE. C C C COMMON /NSCOM1/ PASSES INFORMATION FROM NSDTST TO CNTROL AND STATS. C C NSDTST C | CNTROL C | | STATS C | | | DIFNRM C | | | | C S A A - ERRTOL DOUBLE. COPY OF CURRENT ERROR TOLERANCE. C S A A - OPT INTEGER. COPY OF OPTION(1) ARGUMENT OF NSDTST. C S - - A NRMTYP INTEGER. COPY OF OPTION(3) ARGUMENT OF NSDTST. C S - A - XTRAP INTEGER. COPY OF OPTION(4) ARGUMENT OF NSDTST. C S A - - ID INTEGER. INTERNAL CODE OF CURRENT PROBLEM, 1 FOR A1, C ..., 13 FOR B3, ETC. C S A - - IWT INTEGER. FLAG FOR SCALING (+1: SCALED. -1: C UNSCALED) C S - - - IOUT INTEGER. STANDARD OUTPUT UNIT NUMBER. C C C C C COMMON /NSCOM2/ COMMUNICATES BETWEEN CNTROL AND STATS. C C CNTROL C | STATS C | | C S A XEND DOUBLE. END OF INTEGRATION RANGE OF CURRENT PROBLEM. C A S HSTART DOUBLE. INITIAL STEPSIZE PASSED TO METHOD FOR C INTEGRATION PROPER. C S A N INTEGER. NO. OF EQUATIONS IN CURRENT PROBLEM. C S A IFLAG INTEGER. SET BY CNTROL TO INFORM STATS WHAT IT IS TO C DO: C =0 METHOD IS BEING TIMED. C =1 INITIALIZING CALL OF STATS FROM CNTROL TO SET UP C NSCOM4. C =2 PRELIMINARY INTEGRATION TO DETERMINE HSTART, ABORTED C AFTER 2 STEPS. C =3 INTEGRATION PROPER, COMPILING STATISTICS. C C C A SA INDL,INDG C ERROR FLAGS FOR THE LOCAL AND GLOBAL 'TRUE SOLUTIONS' C OBTAINED BY CALLS TO ROUTINE TRUE. C C C C C C COMMON /NSCOM3/ OUTPUTS STATISTICS FROM CNTROL AND STATS. C C NSDTST C | CNTROL C | | STATS C | | | C A S - XFIN DOUBLE. POINT OF FAILURE OF METHOD IF IT DOESN'T REACH C XEND. C A - S XTRUE DOUBLE. POINT OF FAILURE OF TRUE IF ANY. IF BOTH C LOCAL AND GLOBAL FAIL, POINT OF GLOBAL FAILURE IS C RETURNED. C A S - TIME REAL. CPU TIME FOR ONE INTEGRATION AS MEASURED BY C CLOCK FUNCTION. C A S - OVHD REAL. EQUALS TIME LESS ESTIMATED COST OF FCN CALLS. C A I U TRUTIM REAL. THE TIME SPENT IN CALLS TO TRUE. NOT RELEVANT C TO PERFORMANCE OF METHOD BUT MEASURES THE OVERHEAD C INCURRED BY THE TESTING PACKAGE WHEN OPT = 2, 3 OR 4. C NOT PRINTED BUT AVAILABLE. C A S - GEND REAL. NORM OF GLOBAL ERROR OF METHOD AT XEND. C C C A I U GEMX REAL. MAXIMUM OF GLOBAL ERROR OVER ALL LUMPED STEP C MESHPOINTS, IE. USUALLY OVER ALL MESHPOINTS OF METHOD, C EXCEPT WHEN ERRTOL IS VERY SMALL. C A I U LEMXSC REAL. MAXIMUM LOCAL ERROR IN UNITS OF ERRBND, OVER ALL C LUMPED STEP MESHPOINTS. C A S - NFCN INTEGER. COPY OF NFCN1, SEE /NSCOM6/. C /NSCOM6/ C A I U NSTP INTEGER. COUNTS (UNLUMPED) STEPS TAKEN BY METHOD IN C CURRENT INTEGRATION. C - I U NSTL INTEGER. COUNTS LUMPED STEPS FORMED IN CURRENT C INTEGRATION (SEE STATS). NOT PRINTED BUT AVAILABLE. C A I U NDCV,NBAD C INTEGER. COUNT LUMPED STEPS ON WHICH SOLVER'S LOCAL C ERROR CONTROL WAS DECEIVED, RESP. BADLY DECEIVED. C A I U NTRU INTEGER. COUNTS LUMPED STEPS ON WHICH TRUE LOCAL C SOLUTION WAS SUCCESSFULLY COMPUTED, HENCE VALID LOCAL C ERROR STATISTICS OBTAINED. USED IN COMPUTING 'FRACTION C DECEIVED' INFORMATION. REPORTED IF DIFFERENT FROM C NSTP. NOTE NTRU <= NSTL <= NSTP. C - S - NSTART INTEGER. NO. OF FCN CALLS NEEDED BY METHOD TO START, C IE. TO DO PRELIMINARY INTEGRATION (2 STEPS). NOT C PRINTED OUT BUT AVAILABLE. C C C COMMON /NSCOM4/ IS USED ONLY BY STATS, TO PRESERVE INFORMATION FROM C ONE CALL OF STATS TO ANOTHER. ALL VARIABLES ARE SET AND/OR UPDATED IN C STATS. C C XOLD1 DOUBLE. SIMILAR TO XOLD BUT USED IN PRELIMINARY C INTEGRATION. C XOLD,YOLD C DOUBLE AND DOUBLE ARRAY. COPY OF METHOD'S COMPUTED C SOLUTION AT END OF PREVIOUS LUMPED STEP. USED AS C ACTUAL ARGUMENTS OF TRUE LOCAL SOLUTION CALL. C XOLDG,YOLDG C DOUBLE AND DOUBLE ARRAY. HOLD 'TRUE' GLOBAL SOLUTION C UPDATED TO END OF PREVIOUS LUMPED STEP. USED AS ACTUAL C ARGUMENTS OF TRUE GLOBAL SOLUTION CALL. C CG,PDG,WKG,WG,YPG,INFG C WORKSPACE FOR 'TRUE' GLOBAL SOLUTION. C XT DOUBLE. HOLDS LAST METHOD MESHPOINT BETWEEN CALLS TO C STATS. C PRECIS DOUBLE. HOLDS 1000 * (UNIT ROUNDOFF) APPROX. C ERLUMP DOUBLE. ACCUMULATES METHOD'S LOCAL ERROR ESTIMATES TO C FORM AN ESTIMATE OVER A LUMPED STEP. C C C COMMON /NSCOM5/ PASSES INFORMATION BETWEEN CNTROL AND FCN, (OR ANY C REPLACEMENT A USER MAY PROVIDE FOR FCN). C C CNTROL C | FCN C | | C C S A WT DOUBLE. ARRAY OF WEIGHTS USED TO IMPLEMENT THE C 'SCALED' INTEGRATION OPTION. C S A IWT1,N1,ID1 C INTEGER. COPIES OF IWT,N,ID IN /NSCOM1/ OR /NSCOM2/. C C C COMMON /NSCOM6/ HOLDS A COUNTER. IT IS INITIALIZED IN CNTROL, C SAVED-AND-RESTORED IN STATS, AND EVENTUALLY COPIED BY CNTROL TO THE C CORRESPONDING VARIABLE IN /NSCOM3/. C C CNTROL C | STATS C | | FCN C | | | C C IA AS U - - NFCN1 INTEGER. COUNTS CALLS TO FCN. C C C THERE IS ALSO A COMMON/NSCOM7/ USED BY THE DUMMY (DEBUGGING) VERSIONS C OF NSDTST AND STATS FOR COMMUNICATION. C C--------+---------+---------+---------+---------+---------+---------+-- C E N D O F G E N E R A L D O C U M E N T A T I O N C********+*********+*********+*********+*********+*********+*********+** C C DESCRIPTION OF NSDTST C ----------- -- ------ C C ROUTINE NSDTST INTERPRETS THE LIST OF TOLERANCES AND LIST OF C GROUPS OF PROBLEMS SPECIFIED IN THE ARGUMENTS. USING CNTROL C TO GATHER INDIVIDUAL STATISTICS FOR ONE PROBLEM AT ONE C TOLERANCE, IT ORGANIZES THE FORMATION AND OUTPUT OF SUMMARY C STATISTICS. C INDIVIDUAL STATISTICS ARE INDEXED OVER TOLERANCES, PROBLEMS C AND GROUPS. C 'PROBLEMS-SUMMARY' MEANS SUM OF THESE OVER PROBLEMS IN A GROUP. C 'GROUPS-SUMMARY' MEANS SUM OF PROBLEMS-SUMMARY OVER ALL GROUPS. C 'OVERALL-SUMMARY' MEANS SUM OF GROUPS-SUMMARIES OVER ALL C TOLERANCES. C (READ 'MAX' FOR 'SUM' IN CASE OF SOME OF THE STATISTICS.) C C LOCAL VARIABLES C PSNFCN,PSNSTP,... HOLD THE SUMMARY OVER PROBLEMS IN A GROUP C OF NFCN,NSTP,... (SEE DESCRIPTION OF /NSCOM3/) AT ALL THE C TOLERANCES USED. C GSNFCN,... HOLD SUMMARY OVER GROUPS OF PSNFCN,... C OSNFCN,... HOLD OVERALL SUMMARY (OVER TOLERANCES) OF GSNFCN,... C C LGTOL HOLDS LOGARITHMS TO BASE 10 OF ELEMENTS OF ARRAY TOL, C AND LGGEMX,LGGEND HOLD LOGARITHMS OF CORRESPONDING GEMX C AND GEND VALUES, USED IN SMOOTHNESS CALCULATIONS. C NSNFCN,... STORE NFCN,... FOR ONE PROBLEM AT ALL TOLERANCES C USED, FOR USE IN NORMALIZED EFFICIENCY CALCULATIONS. C ERFLGE,ERFLG1 FLAG 'MISSING VALUES' IN SMOOTHNESS AND NORMALIZED C EFFICIENCY CALCULATIONS. C C C--------+---------+---------+---------+---------+---------+---------+-- C COMMON AREAS C--------+---------+---------+---------+---------+---------+---------+-- C1 C3 C .. Scalar Arguments .. REAL FLAG CHARACTER*80 TITLE C .. Array Arguments .. REAL TOL(11) INTEGER IDLIST(60), OPTION(10) C .. Scalars in Common .. DOUBLE PRECISION ERRTOL, XFIN, XTRUE REAL GEMX, GEND, LEMXSC, OVHD, TIME, TRUTIM INTEGER ID, IOUT, IWT, NBAD, NDCV, NFCN, NRMTYP, NSTART, * NSTL, NSTP, NTRU, OPT, XTRAP C .. Local Scalars .. REAL BIG, C, C1, CTEN, CTEN1, DUM, E, E1, FBADEC, * FDECEV, GEMXSC, GENDSC, OSLEMX, OSOVHD, OSTIME, * RES, RES1, TOLK INTEGER CMPLET, I, ICH, IDSUB, IID, INDG1, INDL1, * KCLASS, KGRP, KSYST, KTOL, NGRP, NOK, NOK1, * NORMEF, NSYST, NTOL, OSNBAD, OSNDCV, OSNFCN, * OSNSTP, OSNTRU CHARACTER BL CHARACTER*10 IDCLAS CHARACTER*32 MCNAME C .. Local Arrays .. REAL GSLEMX(10), GSOVHD(10), GSTIME(10), LGGEMX(10), * LGGEND(10), LGTOL(10), NSOVHD(10), NSTIME(10), * PSGEMX(10), PSGEND(10), PSLEMX(10), PSOVHD(10), * PSTIME(10) INTEGER GRPLST(2,6), GSNBAD(10), GSNDCV(10), GSNFCN(10), * GSNSTP(10), GSNTRU(10), NSNFCN(10), NSNSTP(10), * PSNBAD(10), PSNDCV(10), PSNFCN(10), PSNSTP(10), * PSNTRU(10) LOGICAL ERFLG1(10), ERFLGE(10) C .. External Functions .. REAL CONST, RATIO EXTERNAL CONST, RATIO C .. External Subroutines .. EXTERNAL CNTROL, EFSTAT, LSQFIT, PARCHK, PLOT C .. Intrinsic Functions .. INTRINSIC ALOG10, AMAX1, CHAR, DBLE, IABS, ISIGN C .. Common blocks .. COMMON /NSCOM1/ERRTOL, OPT, NRMTYP, XTRAP, ID, IWT, * IOUT COMMON /NSCOM3/XFIN, XTRUE, TIME, OVHD, TRUTIM, GEND, * GEMX, LEMXSC, NFCN, NSTP, NSTL, NDCV, NBAD, * NTRU, NSTART C .. Data statements .. CE C DATA IDCLAS/'ABCDEFGHIJ'/, BL/' '/, BIG/1.E20/ C .. Executable Statements .. C C--------+---------+---------+---------+---------+---------+---------+-- C COPY THE ENTRIES IN ARRAY 'OPTION'. C DO DUMMY CALL TO CONST TO INVOKE MACHINE-DEPENDENT INITIALIZ- C ATIONS. SET MACHINE NAME. SET OUTPUT UNIT NUMBER. C WRITE OUTPUT-HEADING. CALL ARGUMENT-CHECKING ROUTINE. C--------+---------+---------+---------+---------+---------+---------+-- OPT = OPTION(1) NORMEF = OPTION(2) NRMTYP = OPTION(3) XTRAP = OPTION(4) DUM = CONST(0) DO 20 I = 1, 32 ICH = CONST(-I) MCNAME(I:I) = CHAR(ICH) 20 CONTINUE IOUT = CONST(3) C WRITE (IOUT,FMT=99999) OPT, NORMEF, NRMTYP, MCNAME C CALL PARCHK(OPT,NORMEF,NRMTYP,TOL,IDLIST,NTOL,NGRP,GRPLST,LGTOL, * FLAG) IF (FLAG.EQ.0.) GO TO 40 WRITE (IOUT,FMT=99998) FLAG RETURN C C--------+---------+---------+---------+---------+---------+---------+-- C INITIALIZE OVERALL- AND GROUPS-SUMMARY STATISTICS. C--------+---------+---------+---------+---------+---------+---------+-- 40 OSTIME = 0. OSOVHD = 0. OSNFCN = 0 OSNSTP = 0 OSNTRU = 0 OSLEMX = 0. OSNDCV = 0 OSNBAD = 0 DO 60 I = 1, NTOL GSTIME(I) = 0. GSOVHD(I) = 0. GSNFCN(I) = 0 GSNSTP(I) = 0 GSNTRU(I) = 0 GSLEMX(I) = 0. GSNDCV(I) = 0 GSNBAD(I) = 0 60 CONTINUE C C--------+---------+---------+---------+---------+---------+---------+-- C LOOP OVER GROUPS OF PROBLEMS C--------+---------+---------+---------+---------+---------+---------+-- C DO 300 KGRP = 1, NGRP C C--------+---------+---------+---------+---------+---------+---------+-- C OUTPUT HEADING, ON NEW PAGE FOR GROUPS AFTER FIRST. C SELECT GROUP OF DIFFERENTIAL EQUATIONS. C GET NO. OF SYSTEMS IN THIS GROUP, & OFFSET FOR C POSITION OF ITEM IN GROUP WITHIN IDLIST. C INITIALIZE PROBLEM SUMMARY STATISTICS. C--------+---------+---------+---------+---------+---------+---------+-- IF (KGRP.GT.1) WRITE (IOUT,FMT=99997) WRITE (IOUT,FMT=99996) KGRP, TITLE C NSYST = GRPLST(1,KGRP) IDSUB = GRPLST(2,KGRP) C DO 80 I = 1, NTOL PSTIME(I) = 0. PSOVHD(I) = 0. PSNFCN(I) = 0 PSNSTP(I) = 0 PSNTRU(I) = 0 PSLEMX(I) = 0. PSNDCV(I) = 0 PSNBAD(I) = 0 PSGEMX(I) = 0. PSGEND(I) = 0. 80 CONTINUE C C--------+---------+---------+---------+---------+---------+---------+-- C LOOP OVER PROBLEMS WITHIN A GROUP C--------+---------+---------+---------+---------+---------+---------+-- DO 260 KSYST = 1, NSYST C--------+---------+---------+---------+---------+---------+---------+-- C GET NEXT PROBLEM-ID: C EXTRACT THE WEIGHTING OPTION (IWT=1 OR -1). C UNPACK ID INTO CLASSNAME + INDEX WITHIN CLASS AND TRANSLATE C INTO NSDTST INTERNAL ID BY SUBTRACTING 10: C--------+---------+---------+---------+---------+---------+---------+-- IDSUB = IDSUB + 1 ID = IDLIST(IDSUB) IWT = ISIGN(1,ID) ID = IABS(ID) KCLASS = (ID-1)/10 IID = ID - 10*KCLASS ID = ID - 10 IF (IWT.GT.0) WRITE (IOUT,FMT=99995) IDCLAS(KCLASS:KCLASS), * IID IF (IWT.LE.0) WRITE (IOUT,FMT=99994) IDCLAS(KCLASS:KCLASS), * IID WRITE (IOUT,FMT=99993) (BL,I=1,OPT) WRITE (IOUT,FMT=99992) (BL,I=1,OPT) C C--------+---------+---------+---------+---------+---------+---------+-- C LOOP OVER TOLERANCES FOR ONE PROBLEM C--------+---------+---------+---------+---------+---------+---------+-- DO 220 KTOL = 1, NTOL C--------+---------+---------+---------+---------+---------+---------+-- C CALL PLOT TO INITIALIZE LOCAL-ERROR SCATTER DIAGRAM C IF OPT=4. C CALL CNTROL TO ORGANIZE THE COLLECTION OF C STATISTICS. C ON EXIT FROM CNTROL THE VALUE OF CMPLET WILL C INDICATE WHETHER A FAILURE OCCURRED. C C CMPLET = 1 NO FAILURES. C CMPLET = 0 DETEST FAILED TO OBTAIN TRUE C LOCAL OR GLOBAL SOLUTION. C CMPLET = -1 METHOD FAILED TO REACH THE END C OF RANGE. C CMPLET = -2 DETEST FAILED AND SUBSEQUENTLY C METHOD FAILED. C CMPLET = -3 METHOD COULD NOT START THE C INTEGRATION. C CMPLET = -4 METHOD COMPLETED THE STATISTICS C GATHERING BUT FAILED IN TIMING LOOP. C C ON EXIT INDG1,INDL1 HOLD EXIT-FLAGS OF 'TRUE' C GLOBAL AND LOCAL SOLUTIONS RESPECTIVELY. C C ERFLGE(KTOL) IS TRUE IF METHOD FAILED TO REACH XEND. C ERFLG1(KTOL) IS TRUE IF EITHER METHOD OR C TRUE-SOLUTION FAILED TO REACH XEND (THUS INVALIDATING C GEMX AS DATA FOR SMOOTHNESS CALC WHEN NORMEF=2 ). C C IF CMPLET IS -4,-2,-1,0 OR 1 PRINT A LINE OF STATISTICS: C IF CMPLET ISNT 1, PRINT AN ERROR MESSAGE. C CALL PLOT TO PRINT LOCAL-ERROR SCATTER DIAGRAM C IF OPT=4 C NOTE IF METHOD FAILED TO REACH XEND, ANY STATISTICS FOR C THIS PROBLEM ARE PRINTED BUT DO NOT CONTRIBUTE TO THE C SUMMARY STATISTICS. CONVERSELY IF METHOD REACHED XEND, C ALL STATISTICS CONTRIBUTE TO THE SUMMARIES THOUGH GEMX, C LEMXSC,NDCV,NBAD,NTRU ONLY APPLY TO PART OF THE RANGE C IF 'TRUE' FAILED. C--------+---------+---------+---------+---------+---------+---------+-- C TOLK = TOL(KTOL) ERRTOL = DBLE(TOLK) IF (OPT.EQ.4) CALL PLOT(0.,0.,0) C CALL CNTROL(CMPLET,INDG1,INDL1) C ERFLGE(KTOL) = CMPLET .LT. 0 .AND. CMPLET .GT. -4 ERFLG1(KTOL) = CMPLET .LT. 1 .AND. CMPLET .GT. -4 GENDSC = BIG IF (ERFLGE(KTOL)) GO TO 100 GENDSC = GEND/TOLK LGGEND(KTOL) = ALOG10(AMAX1(GEND,.01*TOLK)) 100 CONTINUE GEMXSC = GEMX/TOLK FDECEV = RATIO(NDCV,NTRU) FBADEC = RATIO(NBAD,NTRU) C IF (CMPLET.EQ.-3) GO TO 120 IF (OPT.EQ.1) WRITE (IOUT,FMT=99991) LGTOL(KTOL), TIME, * OVHD, NFCN, NSTP, GENDSC IF (OPT.EQ.2) WRITE (IOUT,FMT=99991) LGTOL(KTOL), TIME, * OVHD, NFCN, NSTP, GENDSC, GEMXSC IF (OPT.GE.3) WRITE (IOUT,FMT=99991) LGTOL(KTOL), TIME, * OVHD, NFCN, NSTP, GENDSC, GEMXSC, LEMXSC, FDECEV, * FBADEC IF (OPT.GE.3 .AND. NSTP.NE.NTRU) WRITE (IOUT,FMT=99990) * NTRU 120 CONTINUE C C IF (CMPLET.EQ.-4) WRITE (IOUT,FMT=99989) IF (CMPLET.EQ.-3) WRITE (IOUT,FMT=99988) LGTOL(KTOL) C IF (CMPLET.EQ.-2) WRITE (IOUT,FMT=99987) XTRUE, INDG1, * INDL1, XFIN C IF (CMPLET.EQ.-1) WRITE (IOUT,FMT=99986) XFIN C IF (CMPLET.EQ.0) WRITE (IOUT,FMT=99985) XTRUE, INDG1, * INDL1 C IF (OPT.EQ.4) THEN C WRITE (IOUT,FMT=99984) XTRAP C CALL PLOT(0.,0.,2) END IF C FOR EVALUATING PERFORMANCE OF 'TRUE': C CALL TRUCHK(4,IDUM) C C--------+---------+---------+---------+---------+---------+---------+-- C UPDATE PROBLEMS-SUMMARY STATS IF METHOD REACHED XEND. C (IF IT DIDN'T, DON'T UPDATE THE LOCAL-ASSESSMENT INFO: C NTRU,LEMXSC,NDCV,NBAD. THIS IS AN ARBITRARY CHOICE, IT C MAKES IT SIMPLER TO EXPLAIN TO THE USER. C STORE NORMEF STATISTICS: C--------+---------+---------+---------+---------+---------+---------+-- C IF (ERFLGE(KTOL)) GO TO 180 PSTIME(KTOL) = PSTIME(KTOL) + TIME PSOVHD(KTOL) = PSOVHD(KTOL) + OVHD PSNFCN(KTOL) = PSNFCN(KTOL) + NFCN PSNSTP(KTOL) = PSNSTP(KTOL) + NSTP PSGEND(KTOL) = AMAX1(PSGEND(KTOL),GENDSC) C IF (OPT.LT.2) GO TO 140 PSGEMX(KTOL) = AMAX1(PSGEMX(KTOL),GEMXSC) LGGEMX(KTOL) = ALOG10(AMAX1(GEMX,.01*TOLK)) C 140 IF (OPT.LT.3) GO TO 160 PSNTRU(KTOL) = PSNTRU(KTOL) + NTRU PSLEMX(KTOL) = AMAX1(PSLEMX(KTOL),LEMXSC) PSNDCV(KTOL) = PSNDCV(KTOL) + NDCV PSNBAD(KTOL) = PSNBAD(KTOL) + NBAD 160 CONTINUE 180 CONTINUE C IF (NORMEF.EQ.0) GO TO 200 NSTIME(KTOL) = TIME NSOVHD(KTOL) = OVHD NSNFCN(KTOL) = NFCN NSNSTP(KTOL) = NSTP 200 CONTINUE C--------+---------+---------+---------+---------+---------+---------+-- C END OF LOOP OVER TOLERANCES FOR ONE PROBLEM C--------+---------+---------+---------+---------+---------+---------+-- 220 CONTINUE C C--------+---------+---------+---------+---------+---------+---------+-- C SMOOTHNESS AND NORMALIZED EFFICIENCY CALCULATIONS BEGIN C--------+---------+---------+---------+---------+---------+---------+-- WRITE (IOUT,FMT=99983) C WRITE (IOUT,FMT=99982) C CALL LSQFIT(LGTOL,LGGEND,ERFLGE,NTOL,NOK,C,E,RES) C CTEN = 10.**C IF (NOK.LE.2) WRITE (IOUT,FMT=99981) NOK C IF (NOK.GT.2) WRITE (IOUT,FMT=99980) CTEN, E, RES, NOK C IF (OPT.LT.2) GO TO 240 WRITE (IOUT,FMT=99979) C CALL LSQFIT(LGTOL,LGGEMX,ERFLG1,NTOL,NOK1,C1,E1,RES1) C CTEN1 = 10.**C1 IF (NOK1.LE.2) WRITE (IOUT,FMT=99981) NOK1 IF (NOK1.GT.2) WRITE (IOUT,FMT=99980) CTEN1, E1, RES1, NOK1 240 CONTINUE C IF (NORMEF.EQ.1) CALL EFSTAT(C,E,LGTOL,NTOL,NOK,ERFLGE, * 'ENDPOINT',IOUT,NSTIME,NSOVHD, * NSNFCN,NSNSTP) C IF (NORMEF.EQ.2) CALL EFSTAT(C1,E1,LGTOL,NTOL,NOK1,ERFLG1, * 'MAXIMUM ',IOUT,NSTIME,NSOVHD, * NSNFCN,NSNSTP) C C--------+---------+---------+---------+---------+---------+---------+-- C SMOOTHNESS AND NORMALIZED EFFICIENCY CALCULATIONS END C--------+---------+---------+---------+---------+---------+---------+-- C C--------+---------+---------+---------+---------+---------+---------+-- C END OF LOOP OVER PROBLEMS IN A GROUP. C--------+---------+---------+---------+---------+---------+---------+-- 260 CONTINUE C C--------+---------+---------+---------+---------+---------+---------+-- C OUTPUT PROBLEMS-SUMMARY STATISTICS C--------+---------+---------+---------+---------+---------+---------+-- C WRITE (IOUT,FMT=99978) KGRP WRITE (IOUT,FMT=99993) (BL,I=1,OPT) WRITE (IOUT,FMT=99992) (BL,I=1,OPT) DO 280 KTOL = 1, NTOL FDECEV = RATIO(PSNDCV(KTOL),PSNTRU(KTOL)) FBADEC = RATIO(PSNBAD(KTOL),PSNTRU(KTOL)) C IF (OPT.EQ.1) WRITE (IOUT,FMT=99991) LGTOL(KTOL), * PSTIME(KTOL), PSOVHD(KTOL), PSNFCN(KTOL), PSNSTP(KTOL), * PSGEND(KTOL) C IF (OPT.EQ.2) WRITE (IOUT,FMT=99991) LGTOL(KTOL), * PSTIME(KTOL), PSOVHD(KTOL), PSNFCN(KTOL), PSNSTP(KTOL), * PSGEND(KTOL), PSGEMX(KTOL) C IF (OPT.GE.3) WRITE (IOUT,FMT=99991) LGTOL(KTOL), * PSTIME(KTOL), PSOVHD(KTOL), PSNFCN(KTOL), PSNSTP(KTOL), * PSGEND(KTOL), PSGEMX(KTOL), PSLEMX(KTOL), FDECEV, FBADEC C IF (OPT.GE.3 .AND. PSNSTP(KTOL).NE.PSNTRU(KTOL)) * WRITE (IOUT,FMT=99990) PSNTRU(KTOL) C C--------+---------+---------+---------+---------+---------+---------+-- C UPDATE GROUPS-SUMMARY STATISTICS C--------+---------+---------+---------+---------+---------+---------+-- GSTIME(KTOL) = GSTIME(KTOL) + PSTIME(KTOL) GSOVHD(KTOL) = GSOVHD(KTOL) + PSOVHD(KTOL) GSNFCN(KTOL) = GSNFCN(KTOL) + PSNFCN(KTOL) GSNSTP(KTOL) = GSNSTP(KTOL) + PSNSTP(KTOL) C IF (OPT.LT.3) GO TO 280 GSNTRU(KTOL) = GSNTRU(KTOL) + PSNTRU(KTOL) GSLEMX(KTOL) = AMAX1(GSLEMX(KTOL),PSLEMX(KTOL)) GSNDCV(KTOL) = GSNDCV(KTOL) + PSNDCV(KTOL) GSNBAD(KTOL) = GSNBAD(KTOL) + PSNBAD(KTOL) 280 CONTINUE C C--------+---------+---------+---------+---------+---------+---------+-- C END OF LOOP OVER GROUPS C--------+---------+---------+---------+---------+---------+---------+-- 300 CONTINUE C C C--------+---------+---------+---------+---------+---------+---------+-- C OUTPUT HEADINGS FOR GROUPS- AND OVERALL-SUMMARY STATISTICS. C--------+---------+---------+---------+---------+---------+---------+-- WRITE (IOUT,FMT=99977) TITLE, (BL,I=1,OPT) WRITE (IOUT,FMT=99976) (BL,I=1,OPT) C--------+---------+---------+---------+---------+---------+---------+-- C OUTPUT GROUPS-SUMMARY STATISTICS C--------+---------+---------+---------+---------+---------+---------+-- IF (OPT.GE.3) GO TO 340 DO 320 I = 1, NTOL WRITE (IOUT,FMT=99975) LGTOL(I), GSTIME(I), GSOVHD(I), * GSNFCN(I), GSNSTP(I) 320 CONTINUE GO TO 380 340 DO 360 I = 1, NTOL FDECEV = RATIO(GSNDCV(I),GSNTRU(I)) FBADEC = RATIO(GSNBAD(I),GSNTRU(I)) WRITE (IOUT,FMT=99975) LGTOL(I), GSTIME(I), GSOVHD(I), * GSNFCN(I), GSNSTP(I), GSLEMX(I), FDECEV, FBADEC C IF (GSNSTP(I).NE.GSNTRU(I)) WRITE (IOUT,FMT=99990) GSNTRU(I) 360 CONTINUE 380 CONTINUE C C--------+---------+---------+---------+---------+---------+---------+-- C COMPUTE OVERALL-SUMMARY STATISTICS. C--------+---------+---------+---------+---------+---------+---------+-- DO 400 I = 1, NTOL OSTIME = OSTIME + GSTIME(I) OSOVHD = OSOVHD + GSOVHD(I) OSNFCN = OSNFCN + GSNFCN(I) OSNSTP = OSNSTP + GSNSTP(I) C IF (OPT.LT.3) GO TO 400 OSNTRU = OSNTRU + GSNTRU(I) OSNDCV = OSNDCV + GSNDCV(I) OSNBAD = OSNBAD + GSNBAD(I) OSLEMX = AMAX1(OSLEMX,GSLEMX(I)) 400 CONTINUE FDECEV = RATIO(OSNDCV,OSNTRU) FBADEC = RATIO(OSNBAD,OSNTRU) C--------+---------+---------+---------+---------+---------+---------+-- C OUTPUT OVERALL-SUMMARY STATISTICS C--------+---------+---------+---------+---------+---------+---------+-- IF (OPT.LT.3) WRITE (IOUT,FMT=99974) OSTIME, OSOVHD, OSNFCN, * OSNSTP C IF (OPT.GE.3) WRITE (IOUT,FMT=99974) OSTIME, OSOVHD, OSNFCN, * OSNSTP, OSLEMX, FDECEV, FBADEC C C RETURN C 99999 FORMAT ('0NONSTIFF DETEST PACKAGE OPTION=',I2,', NORMEF=',I2, * ', NRMTYP=',I2,19X,'ON ',A,//) 99998 FORMAT ('0PARAMETER ERRORS AS SHOWN BY FLAG=',E15.8,/' ',49('*') * ,//) 99997 FORMAT ('1') 99996 FORMAT ('0GROUP',I3,18X,A) 99995 FORMAT (/'0',A3,I1,' (SCALED)',/) 99994 FORMAT (/'0',A3,I1,' (UNSCALED)',/) 99993 FORMAT (' ',A1,6X,'LOG10',5X,'TIME',3X,'OVHD',5X,'FCN',4X,'NO OF', * 3X,'END PNT',A1,2X,'MAXIMUM',A1,2X,'MAXIMUM',3X,'FRACTION', * 3X,'FRACTION',A1) 99992 FORMAT (' ',A1,7X,'TOL',21X,'CALLS',3X,'STEPS',3X,'GLB ERR',A1,2X, * 'GLB ERR',A1,2X,'LOC ERR',3X,'DECEIVED',3X,'BAD DECV',A1) 99991 FORMAT ('0',6X,F6.2,2X,2F7.3,1X,2I8,2X,F8.2,1X,F9.2,1X,F9.3,1X, * F9.3,1X,F10.3,1X,F10.3) 99990 FORMAT (114X,'(LOC ASSESS ON',I4,')') 99989 FORMAT ('0',20X, * '***** UNEXPECTED FAILURE OF METHOD WHILE BEING TIMED *****' * ,/) 99988 FORMAT ('0',6X,F6.2,' *** METHOD FAILED TO START ***') 99987 FORMAT (15X,'TRUE-SOLUTION OF TEST PACKAGE FAILED AT X = ',1P, * E12.5,', ERROR FLAG (GLOBAL) ',I3,', (LOCAL) ',I3,/21X, * 'AND SUBSEQUENTLY METHOD FAILED AT X = ',1P,E12.5) 99986 FORMAT (21X,'METHOD FAILED AT X = ',1P,E12.5) 99985 FORMAT (21X,'TRUE-SOLUTION OF TEST PACKAGE FAILED AT X = ',1P, * E12.5,', ERROR FLAG (GLOBAL) ',I3,', (LOCAL) ',I3) 99984 FORMAT (/6X,'ERROR ESTIMATE ANALYSIS',10X, * 'EXTRAPOLATION (0=NO 1=YES):',I2,/11X, * 'HORIZONTAL AXIS: R1=||ERREST|| / ERRBND',/11X, * 'VERTICAL AXIS: R2 = ||ERROR IN ERREST|| / ERRBND',/11X, * 'PLOT SHOWS % STEPS WHERE (R1,R2) LAY',1X, * 'IN INDICATED PIGEONHOLE, A DOT MEANS UNDER 1%',/) 99983 FORMAT (/'0',17X,'SMOOTHNESS FIT OF LOG10(ERROR) VS LOG10(TOL)') 99982 FORMAT ('0',17X,'ENDPOINT GLOBAL ERROR') 99981 FORMAT (39X,I2,' VALUES, TOO FEW TO GET STATISTICS') 99980 FORMAT (39X,'=',1P,G10.3,' *(TOL**',0P,F6.3,') APPROX,',6X, * 'R.M.S. RESIDUAL=',1P,E8.1,' OVER',I3,' VALUES') 99979 FORMAT ('0',17X,'MAXIMUM GLOBAL ERROR') 99978 FORMAT (/'0SUMMARY OVER GROUP',I3) 99977 FORMAT ('1SUMMARY OVER ALL GROUPS',6X,A,//' ',A1,6X,'LOG10',5X, * 'TIME',3X,'OVHD',5X,'FCN',4X,'NO OF',2A1,'MAXIMUM',3X, * 'FRACTION',3X,'FRACTION',A1) 99976 FORMAT (' ',A1,7X,'TOL',21X,'CALLS',3X,'STEPS',2A1,'LOC ERR',3X, * 'DECEIVED',3X,'BAD DECV',A1) 99975 FORMAT ('0',6X,F6.2,2X,2F7.3,1X,2I8,1X,3F11.3) 99974 FORMAT ('0',5X,'OVERALL',/6X,'SUMMARY',2X,2F7.3,1X,2I8,1X,3F11.3) END C C C********+*********+*********+*********+*********+*********+*********+** C SUBROUTINE PARCHK(OPT,NORMEF,NRMTYP,TOL,IDLIST,NTOL,NGRP,GRPLST, * LGTOL,FLAG) C C********+*********+*********+*********+*********+*********+*********+** C ROUTINE TO DO PARAMETER CHECKS FOR REVISED NSDTST INTERFACE. C C INPUT: OPT,NORMEF,NRMTYP,TOL,IDLIST C VALID INPUT IS: C OPTION = 1 2 3 OR 4 C NORMEF = 0 1 OR 2 C NRMTYP = 1 2 OR 3 C TOL = LIST OF UP TO 10 POSITIVE REAL'S FOLLOWED BY A 0., C IN STRICTLY DECREASING ORDER C IDLIST = LIST OF GROUPS OF PROBLEM-IDS SEPARATED BY ZEROS C WITH 2 ZEROS AFTER LAST GROUP, AT MOST 60 ITEMS TOTAL. C EACH ID MAY HAVE A MINUS SIGN TO SELECT THE 'UNSCALED' C ERROR CONTROL OPTION. C VALID PROBLEM-IDS ARE IN RANGES C 11-15 21-25 31-35 41-45 51-55 61-65 C FOR PROBLEM CLASSES A1-A5 B1-B5 ETC. C OUTPUT: NTOL = NO. OF TOLERANCES IN TOL LIST C NGRP = NO. OF GROUPS IN IDLIST LIST C GRPLST(1,I) = SIZE OF I-TH GROUP OF PROBLEMS CC ... (2,I) = POINTER TO (START OF I-TH GROUP)-1 IN IDLIST C LGTOL(I) = LOG10(TOL(I)) C FLAG IS ERROR FLAG, 0.0 IF ALL OK, ELSE ITS DECIMAL DIGITS C INDICATE WHICH PARAMETER ERRORS WERE FOUND: C 1: OPT INVALID C 2: NORMEF INVALID C 3: NORMEF = 2 REQUESTED WITH OPT = 1 C 4: TOL(I) < 0, OR LIST NOT IN DECREASING ORDER C 5: TOL LIST EMPTY OR NOT TERMINATED BY ZERO C 6: INVALID PROBLEM-ID FOUND C 7: LIST OF GROUPS IN IDLIST EMPTY,NOT TERMINATED BY C 2 ZEROS OR HAS MORE THAN MAXGRP GROUPS C 8: NRMTYP INVALID C--------+---------+---------+---------+---------+---------+---------+-- C C .. Scalar Arguments .. REAL FLAG INTEGER NGRP, NORMEF, NRMTYP, NTOL, OPT C .. Array Arguments .. REAL LGTOL(10), TOL(11) INTEGER GRPLST(2,6), IDLIST(60) C .. Local Scalars .. REAL BIG, TOLPRV INTEGER ENDLST, I, ID, IID, ISAV, KCLASS, LENIDS, * LENTOL, MAXGRP, NCLASS C .. Local Arrays .. INTEGER NSYSTM(6) C .. Intrinsic Functions .. INTRINSIC ALOG10, IABS C .. Data statements .. DATA ENDLST/-1/, BIG/1E20/ DATA NCLASS/6/, NSYSTM/5, 5, 5, 5, 5, 5/, MAXGRP/6/, * LENTOL/11/, LENIDS/60/ C .. Executable Statements .. C FLAG = 0. IF (OPT.LT.1 .OR. OPT.GT.4) FLAG = 1. IF (NORMEF.LT.0 .OR. NORMEF.GT.2) FLAG = 10.*FLAG + 2. IF (OPT.EQ.1 .AND. NORMEF.EQ.2) FLAG = 10.*FLAG + 3. IF (NRMTYP.LT.1 .OR. NRMTYP.GT.3) FLAG = 10.*FLAG + 8. C C TOLERANCES: NTOL = 0 TOLPRV = BIG DO 20 I = 1, LENTOL IF (TOL(I).LT.0. .OR. TOL(I).GE.TOLPRV) FLAG = 10.*FLAG + 4. IF (TOL(I).EQ.0.) GO TO 40 NTOL = NTOL + 1 TOLPRV = TOL(I) 20 CONTINUE C C NO TERMINATING 0 IN TOLERANCE LIST: FLAG = 10.*FLAG + 5. C C CHECK FOR EMPTY TOLERANCE LIST: 40 IF (NTOL.EQ.0) FLAG = 10.*FLAG + 5. C C LIST OF GROUPS OF PROBLEMS: NGRP = 0 I = 0 C C WHILE NEXT ID IN LIST ISNT 0 OR END OF LIST: 60 I = I + 1 ID = ENDLST IF (I.LE.LENIDS) ID = IDLIST(I) C IF (ID.EQ.0) GO TO 160 IF (NGRP.GE.MAXGRP) GO TO 180 ISAV = I - 1 C C WHILE ID ISNT 0, GET ONE GROUP: 80 IF (ID.EQ.0) GO TO 140 IF (ID.EQ.ENDLST) GO TO 180 C TRANSLATE ID INTO CLASS & NUMBER WITHIN CLASS, C IGNORING SIGN (WHICH SELECTS SCALED/UNSCALED OPTION): ID = IABS(ID) KCLASS = (ID-1)/10 IID = ID - 10*KCLASS IF ( .NOT. (KCLASS.GE.1 .AND. KCLASS.LE.NCLASS)) GO TO 100 IF (IID.LE.NSYSTM(KCLASS)) GO TO 120 100 FLAG = 10.*FLAG + 6. 120 CONTINUE C GET NEXT ID AS ABOVE: I = I + 1 ID = ENDLST IF (I.LE.LENIDS) ID = IDLIST(I) GO TO 80 C C NEW GROUP FORMED: 140 NGRP = NGRP + 1 GRPLST(1,NGRP) = I - ISAV - 1 GRPLST(2,NGRP) = ISAV GO TO 60 C C CHECK IF NO GROUPS WERE SPECIFIED: 160 IF (NGRP.LE.0) GO TO 180 GO TO 200 C 180 FLAG = 10.*FLAG + 7. C C IF ALL OK, COMPUTE LOGS OF TOLERANCES: C 200 IF (FLAG.NE.0.) GO TO 240 DO 220 I = 1, NTOL LGTOL(I) = ALOG10(TOL(I)) 220 CONTINUE 240 RETURN END C C********+*********+*********+*********+*********+*********+*********+** C SUBROUTINE LSQFIT(X,Y,MISS,N,NN,C0,C1,RES) C .. Scalar Arguments .. REAL C0, C1, RES INTEGER N, NN C .. Array Arguments .. REAL X(N), Y(N) LOGICAL MISS(N) C .. Local Scalars .. REAL SX, SXX, SXY, SY, XNN INTEGER I C .. Intrinsic Functions .. INTRINSIC SQRT C .. Executable Statements .. C C********+*********+*********+*********+*********+*********+*********+** C FITS MODEL Y = C0 + C1*X TO DATA X(I),Y(I),I = 1..N WHERE DATA C FOR WHICH MISS(I) IS .TRUE. IS REGARDED AS MISSING. C C ON EXIT C X,Y,MISS,N ARE UNCHANGED. C NN = NO. OF NONMISSING VALUES C C0,C1 = FITTED COEFFICIENTS C RES = ROOT MEAN SQUARE RESIDUAL C C EXCEPT THAT IF NN.LE.1 NO COMPUTATION OF THE COEFFICIENTS IS DONE. C--------+---------+---------+---------+---------+---------+---------+-- C NN = 0 SX = 0. SY = 0. DO 20 I = 1, N IF (MISS(I)) GO TO 20 NN = NN + 1 SX = SX + X(I) SY = SY + Y(I) 20 CONTINUE IF (NN.LE.1) GO TO 80 XNN = NN SX = SX/XNN SY = SY/XNN SXX = 0. SXY = 0. DO 40 I = 1, N IF (MISS(I)) GO TO 40 SXX = SXX + (X(I)-SX)**2 SXY = SXY + (X(I)-SX)*(Y(I)-SY) 40 CONTINUE C1 = SXY/SXX C0 = SY - C1*SX RES = 0. DO 60 I = 1, N IF ( .NOT. MISS(I)) RES = RES + (Y(I)-SY-C1*(X(I)-SX))**2 60 CONTINUE C RES = SQRT(RES/XNN) C 80 RETURN END C C********+*********+*********+*********+*********+*********+*********+** C REAL FUNCTION RATIO(M,N) C C********+*********+*********+*********+*********+*********+*********+** C .. Scalar Arguments .. INTEGER M, N C .. Intrinsic Functions .. INTRINSIC FLOAT C .. Executable Statements .. RATIO = 1E20 IF (N.NE.0) RATIO = FLOAT(M)/FLOAT(N) RETURN END C C********+*********+*********+*********+*********+*********+*********+** C SUBROUTINE EFSTAT(C,E,LGTOL,NTOL,NOK,ERFLG,TITLE,IOUT,W1,W2,W3,W4) C C********+*********+*********+*********+*********+*********+*********+** C ROUTINE TO COMPUTE AND PRINT NORMALIZED EFFICIENCY STATISTICS. C C PARAMETERS (ALL INPUT): C C,E - COEFFICIENTS IN LEAST-SQUARES FIT OF ACHIEVED ACCURACY C (EITHER AT ENDPOINT OR MAX-OVER-RANGE) TO TOLERANCE. C LGTOL - LIST OF LOGS TO BASE 10 OF TOLERANCES C NTOL - NO. OF TOLERANCES. C NOK - NO. OF .FALSE. ENTRIES IN ERFLG (FROM LSQFIT CALL) C ERFLG - LOGICAL VECTOR INDICATING FOR WHICH TOLERANCES DATA C IS TO BE REGARDED AS MISSING. C TITLE C - IDENTIFYING CHARACTER STRING. C IOUT - OUTPUT UNIT NUMBER. C W1,...,W6 C - VECTORS OF STATISTICS, INDEXED OVER TOLERANCES, FOR C WHICH NORMALIZED STATISTICS ARE TO BE PRODUCED. C (NOTE SOME ARE REAL, SOME INTEGER: REFER TO ACTUAL CALL C IN NSDTST.) C IT IS ASSUMED THAT NTOL.LE.10, OTHERWISE ARRAY S MUST BE LONGER. C--------+---------+---------+---------+---------+---------+---------+-- C C LOCAL VARIABLES C .. Scalar Arguments .. REAL C, E INTEGER IOUT, NOK, NTOL CHARACTER*8 TITLE C .. Array Arguments .. REAL LGTOL(NTOL), W1(NTOL), W2(NTOL) INTEGER W3(NTOL), W4(NTOL) LOGICAL ERFLG(NTOL) C .. Local Scalars .. REAL EQVTOL, S0, THETA, W1INT, W2INT, X INTEGER I, MSINT, NHI, NLO, SHI, SINT, SLO, W3INT, W4INT C .. Local Arrays .. REAL S(10) C .. Intrinsic Functions .. INTRINSIC FLOAT, INT C .. Statement Functions .. INTEGER FLOOR C .. Statement Function definitions .. C C STATEMENT FUNCTION C FLOOR FUNCTION VALID IF ARGUMENT X.GE.-100 WHICH IS OK HERE. FLOOR(X) = INT(X+100.) - 100 C .. Executable Statements .. C IF (NOK.LE.2) GO TO 200 C C TRANSFORM THE LOG10(TOL)'S TO NORMALIZED-EFFICIENCY VARIABLE: DO 20 I = 1, NTOL S(I) = -(C+E*LGTOL(I)) 20 CONTINUE C C FIND SET OF CONSECUTIVE TOL'S FOR WHICH INTEGRATION SUCCEEDED: DO 40 NLO = 1, NTOL IF ( .NOT. ERFLG(NLO)) GO TO 60 40 CONTINUE C ELSE ALL INTEGRATIONS FOR THIS PROBLEM FAILED: GO TO 200 60 CONTINUE NHI = NLO - 1 DO 80 I = NLO, NTOL IF (ERFLG(I)) GO TO 100 NHI = I 80 CONTINUE 100 CONTINUE C IF (NHI.LE.NLO) GO TO 200 IF (E.LE.0.) GO TO 220 C C FORM RANGE OF INTEGER POWERS OF 10 FOR WHICH NORMALIZED STATISTICS C ARE TO BE PRINTED: SLO = -FLOOR(-S(NLO)+0.1) SHI = FLOOR(S(NHI)+0.1) IF (SHI.LT.SLO) GO TO 240 C WRITE (IOUT,FMT=99999) TITLE C C START OF LOOP TO PRINT A LINE OF STATISTICS FOR EACH POWER OF 10: I = NLO + 1 CC ... WHICH IS KNOWN TO BE .LE. NHI C DO 160 SINT = SLO, SHI S0 = FLOAT(SINT) C C MOVE INTERVAL S(I-1)..S(I) TO RIGHT WHILE S(I).LT.SINT: 120 IF (S(I).GE.S0 .OR. I.GE.NHI) GO TO 140 I = I + 1 GO TO 120 140 CONTINUE C NECESSARILY NOW NLO + 1 .LE. I .LE. NHI C C NOW DO INTERPOLATION (POSSIBLY EXTRAPOLATION A SHORT DISTANCE) C USING DATA FOR I AND I + 1: THETA = (S0-S(I-1))/(S(I)-S(I-1)) W1INT = W1(I-1) + THETA*(W1(I)-W1(I-1)) W2INT = W2(I-1) + THETA*(W2(I)-W2(I-1)) W3INT = W3(I-1) + THETA*(W3(I)-W3(I-1)) W4INT = W4(I-1) + THETA*(W4(I)-W4(I-1)) C MSINT = -SINT EQVTOL = -(C+S0)/E WRITE (IOUT,FMT=99998) MSINT, EQVTOL, W1INT, W2INT, W3INT, * W4INT C 160 CONTINUE C 180 RETURN C 200 WRITE (IOUT,FMT=99997) GO TO 180 C 220 WRITE (IOUT,FMT=99996) GO TO 180 C 240 WRITE (IOUT,FMT=99995) GO TO 180 C 99999 FORMAT (/'0',6X,'NORMALIZED EFFICIENCY - ',A8,' GLOBAL ERROR', * //7X,'EXPECTED',3X,'EQUIV',4X,'TIME',3X,'OVHD',5X,'FCN',4X, * 'NO OF',/7X,'ACCURACY',1X,'LOG10 TOL',17X,'CALLS',3X, * 'STEPS') 99998 FORMAT ('0',6X,'10**',I3,F8.2,F9.3,F7.3,1X,2I8) 99997 FORMAT ('0',10X,'NOT ENOUGH SUCCESSFUL INTEGRATIONS TO FORM',1X, * 'NORMALIZED STATISTICS') 99996 FORMAT ('0',10X,'DEPENDENCE OF ACCURACY ON TOLERANCE IS TOO',1X, * 'UNRELIABLE TO FORM NORMALIZED STATISTICS') 99995 FORMAT ('0',10X,'NO POWERS OF TEN WITHIN RANGE OF TOLERANCES',1X, * 'USED: NO NORMALIZED STATISTICS') END C C C********+*********+*********+*********+*********+*********+*********+** C SUBROUTINE CNTROL(CMPLET,INDG1,INDL1) C C********+*********+*********+*********+*********+*********+*********+** C CNTROL ORGANIZES THE CALLS TO METHOD NEEDED TO GATHER C STATISTICS FOR ONE PROBLEM AND ONE TOLERANCE AT THE LEVEL OF C DETAIL SPECIFIED BY OPT, WITH SCALING TURNED ON OR OFF BY IWT. C C ON EXIT FROM CNTROL C CMPLET INDICATES WHETHER A FAILURE OCCURRED: C CMPLET = 1 NO FAILURES. C CMPLET = 0 DETEST FAILED TO OBTAIN TRUE LOCAL OR GLOBAL C SOLUTION. C CMPLET = -1 METHOD FAILED TO REACH THE END OF RANGE. C CMPLET = -2 DETEST FAILED AND SUBSEQUENTLY METHOD FAILED C CMPLET = -3 METHOD COULD NOT START THE INTEGRATION. C CMPLET = -4 METHOD COMPLETED THE STATISTICS GATHERING CALL C BUT (UNEXPECTEDLY) FAILED IN THE TIMING LOOP. C C INDG1, INDL1 RETURN THE ERROR FLAGS OF THE 'TRUE' GLOBAL C AND LOCAL SOLUTIONS RESPECTIVELY. C C THE MAIN OUTPUT FROM CNTROL CONSISTS OF THE STATISTICS HELD C IN COMMON /NSCOM3/ C--------+---------+---------+---------+---------+---------+---------+-- C--------+---------+---------+---------+---------+---------+---------+-- C COMMON AREAS C--------+---------+---------+---------+---------+---------+---------+-- C1 C2 C3 C5 C6 C .. Scalar Arguments .. INTEGER CMPLET, INDG1, INDL1 C .. Scalars in Common .. DOUBLE PRECISION ERRTOL, HSTART, XEND, XFIN, XTRUE REAL GEMX, GEND, LEMXSC, OVHD, TIME, TRUTIM INTEGER ID, ID1, IFLAG, INDG, INDL, IOUT, IWT, IWT1, N, * N1, NBAD, NDCV, NFCN, NFCN1, NRMTYP, NSTART, * NSTL, NSTP, NTRU, OPT, XTRAP C .. Arrays in Common .. DOUBLE PRECISION WT(51) C .. Local Scalars .. DOUBLE PRECISION DUMMY, HINIT, HMAX, X, XSTART REAL FCNTIM, S, TIMCUM, TSTTIM INTEGER COUNT, I LOGICAL NOSTRT, OKMETH, TIMERR C .. Local Arrays .. DOUBLE PRECISION Y(51), YEND(51), YSTART(51) C .. External Functions .. REAL CLOCK, CONST, DIFNRM EXTERNAL CLOCK, CONST, DIFNRM C .. External Subroutines .. EXTERNAL EVALU, IVALU, METHOD, STATS C .. Intrinsic Functions .. INTRINSIC FLOAT C .. Common blocks .. COMMON /NSCOM1/ERRTOL, OPT, NRMTYP, XTRAP, ID, IWT, * IOUT COMMON /NSCOM2/XEND, HSTART, N, IFLAG, INDL, INDG COMMON /NSCOM3/XFIN, XTRUE, TIME, OVHD, TRUTIM, GEND, * GEMX, LEMXSC, NFCN, NSTP, NSTL, NDCV, NBAD, * NTRU, NSTART COMMON /NSCOM5/WT, IWT1, N1, ID1 COMMON /NSCOM6/NFCN1 C .. Executable Statements .. CE C C--------+---------+---------+---------+---------+---------+---------+-- C NOTE ON INDL, INDG IN /NSCOM2/: C THESE ARE ERROR INDICATORS FOR THE 'TRUE' LOCAL AND C GLOBAL SOLUTION RESPECTIVELY. THEY ARE SET INSIDE STATS C WHICH IS CALLED BY METHOD. C ON RETURN FROM METHOD, INDL IS: C 2 IF NO CALL TO TRUE TO COMPUTE LOCAL SOLUTION HAS C YET BEEN MADE (SET BY INITIALIZING CALL TO STATS). C .GT.0 IF ALL CALLS TO TRUE FOR CALCULATION OF LOCAL C SOLUTION WERE SUCCESSFUL. C .LT.0 IF AN UNSUCCESSFUL CALL TO TRUE FOR THE LOCAL C SOLUTION WAS MADE. C THE VALUE ON EXIT IF NOT 0 IS THE VALUE RETURNED IN THE C FLAG 'IND' OF SUBROUTINE TRUE. C INDG IS THE SAME, BUT FOR THE GLOBAL SOLUTION. C C INDL,INDG ARE USED ON RE-ENTRY TO STATS TO TEST IF A C FAILURE OF THE TRUE SOLUTIONS OCCURRED ON A PREVIOUS STEP C AND SHOULD THUS BE LEFT ALONE BETWEEN STEPS. C--------+---------+---------+---------+---------+---------+---------+-- C C ACTION OF THE ROUTINE: C CALL IVALU TO SET INTEGRATION PARAMETERS. C COPY N,ID,IWT INTO /NSCOM5/ FOR USE BY FCN. C SET IFLAG = 1 AND CALL STATS TO INITIALIZE ITS COMMON AREAS. C (THE ARGUMENTS FOR THIS CALL ARE DUMMIES.) C SET X,Y,NSTP,NFCN FOR USE IN STATS. SET IFLAG = 2 SO THAT C THE CALL TO METHOD WILL SET THE FIRST STEP SIZE (HSTART) C AND RETURN. C SET NSTART = NO. OF FCN CALLS NEEDED BY METHOD TO START. C--------+---------+---------+---------+---------+---------+---------+-- C CALL IVALU(N,XSTART,XEND,HINIT,HMAX,YSTART,FCNTIM,WT,IWT,ID) C N1 = N ID1 = ID IWT1 = IWT X = XSTART DO 20 I = 1, N Y(I) = YSTART(I) 20 CONTINUE C IFLAG = 1 CALL STATS(X,Y,DUMMY,Y) C NFCN1 = 0 NSTP = 0 IFLAG = 2 C CALL METHOD(N,X,Y,XEND,ERRTOL,HMAX,HINIT) C NOSTRT = X .LT. XEND NSTART = NFCN1 C--------+---------+---------+---------+---------+---------+---------+-- C INITIALIZE THE COUNTERS ETC. IN /NSCOM3/,/NSCOM6/. C IF METHOD FAILED TO START, SET FLAGS AND EXIT. C SET IFLAG = 3 SO THAT THE CALL TO METHOD WILL DO A COMPLETE C INTEGRATION, COMPILING STATISTICS ON EACH STEP. C START THE CLOCK. C--------+---------+---------+---------+---------+---------+---------+-- NFCN1 = 0 NSTP = 0 NSTL = 0 LEMXSC = 0. NDCV = 0 NBAD = 0 GEMX = 0. TRUTIM = 0. NTRU = 0 C IF (NOSTRT) GO TO 180 C X = XSTART DO 40 I = 1, N Y(I) = YSTART(I) 40 CONTINUE IFLAG = 3 S = CLOCK(0.0) C CALL METHOD(N,X,Y,XEND,ERRTOL,HMAX,HSTART) C TIME = CLOCK(S) OKMETH = X .GE. XEND XFIN = X NFCN = NFCN1 IF ( .NOT. OKMETH) GO TO 160 C--------+---------+---------+---------+---------+---------+---------+-- C IF OPT.GT.1, OR IF OPT = 1 BUT THE TIMING ESTIMATE ALREADY C OBTAINED WAS TOO SMALL TO BE RELIABLE, DO A TIMING COMPUTATION C PROVIDED THAT METHOD REACHED THE ENDPOINT IN THE PREVIOUS CALL. C SET IFLAG = 0, START THE CLOCK, AND CALL C METHOD SUFFICIENTLY MANY TIMES FOR THE SOLUTION TIME TO C BE OBTAINED ACCURATELY. COMPUTE THE OVERHEAD AS THE C TOTAL TIME EXCLUSIVE OF FUNCTION EVALUATIONS C--------+---------+---------+---------+---------+---------+---------+-- TSTTIM = CONST(4) TIMERR = .FALSE. IF (TSTTIM.LE.0) GO TO 120 IF (OPT.EQ.1 .AND. TIME.GE.0.5*TSTTIM) GO TO 120 COUNT = 0 IFLAG = 0 S = CLOCK(0.0) C--------+---------+---------+---------+---------+---------+---------+-- C LOOP TILL 'TSTTIM' TIME UNITS HAVE ELAPSED: C--------+---------+---------+---------+---------+---------+---------+-- 60 CONTINUE X = XSTART DO 80 I = 1, N Y(I) = YSTART(I) 80 CONTINUE CALL METHOD(N,X,Y,XEND,ERRTOL,HMAX,HSTART) TIMERR = X .LT. XEND IF (TIMERR) GO TO 100 TIMCUM = CLOCK(S) COUNT = COUNT + 1 IF (TIMCUM.LT.TSTTIM .AND. COUNT.LT.10) GO TO 60 C 100 IF (COUNT.GE.1) TIME = TIMCUM/FLOAT(COUNT) 120 CONTINUE C--------+---------+---------+---------+---------+---------+---------+-- C WE NOW HAVE A VALUE FOR TIME: THE ONE OBTAINED BEFORE THE C TIMING LOOP IF WE SKIPPED THE LATTER OR IN THE UNLIKELY C EVENT OF AN ERROR IN THE 1ST TIMING ITERATION; OTHERWISE C THE ONE FROM THE TIMING LOOP. C COMPUTE OVERHEAD AND ENDPOINT GLOBAL ERROR. C--------+---------+---------+---------+---------+---------+---------+-- OVHD = TIME - FLOAT(NFCN)*FCNTIM CALL EVALU(YEND,N,WT,IWT,ID) GEND = DIFNRM(YEND,Y,N) C IF (TIMERR) GO TO 200 C C--------+---------+---------+---------+---------+---------+---------+-- C SET THE OUTPUT VALUE OF CMPLET, INDG1 AND INDL1. C--------+---------+---------+---------+---------+---------+---------+-- CMPLET = 1 IF (INDL.LT.0 .OR. INDG.LT.0) CMPLET = 0 140 INDG1 = INDG INDL1 = INDL RETURN C C--------+---------+---------+---------+---------+---------+---------+-- C *********** ERROR EXITS *********** C--------+---------+---------+---------+---------+---------+---------+-- C METHOD FAILED TO REACH XEND C--------+---------+---------+---------+---------+---------+---------+-- 160 CMPLET = -1 IF (INDL.LT.0 .OR. INDG.LT.0) CMPLET = -2 TIME = 1E20 OVHD = 1E20 GEND = 1E20 GO TO 140 C C--------+---------+---------+---------+---------+---------+---------+-- C METHOD FAILED TO START C--------+---------+---------+---------+---------+---------+---------+-- 180 CMPLET = -3 NFCN = 0 TIME = 1E20 OVHD = 1E20 GEND = 1E20 GO TO 140 C--------+---------+---------+---------+