      subroutine nsdtst ( title, option, tol, idlist, flag )

c***********************************************************************
c
cc NSDTST is a test package for nonstiff differential equation solvers.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
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 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 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<INFINITY
c         THUS   FORMING  12X12  PIGEONHOLES.  EACH   INTEGRATION   STEP
c         CONTRIBUTES  A DATA POINT (R1,R2)  WHICH  IS  ENTERED  IN  ONE
c         PIGEONHOLE. THE  COUNTS  OF  THE  NUMBER  OF  ENTRIES  IN EACH
c         PIGEONHOLE ARE EXPRESSED AS integer PERCENTAGES OF  THE  TOTAL
c         NUMBER  OF INTEGRATION STEPS AND PRINTED OUT IN A 12X12 ARRAY,
c         ZERO ENTRIES BEING LEFT BLANK,  AND  POSITIVE  VALUES  BELOW 1
c         BEING SHOWN BY A DOT '.'.
c
c         STEP-LUMPING (SEE [4]) IS DEEMED TO MAKE THIS ANALYSIS USELESS
c         SO  STATISTICS ARE ONLY GATHERED ON UNLUMPED STEPS. IT  IS  AT
c         PRESENT ALSO  NOT  CONSIDERED USEFUL TO PRODUCE SUMMARY TABLES
c         OVER SEVERAL PROBLEMS (AND WOULD BE COSTLY IN ARRAY SPACE).
c
c
c NORMEF  ONE  OF  0   1   OR   2   ,   SELECTS   NORMALIZED  EFFICIENCY
c         STATISTICS.    THESE  TRY  TO  COMPENSATE  FOR  THE  FACT THAT
c         ACHIEVED  ACCURACY  MAY  BE MUCH HIGHER OR LOWER   THAN   THAT
c         REQUESTED  BY  TOL, AND THIS RELATIONSHIP IS VERY PROBLEM- AND
c         METHOD- DEPENDENT.  FOR EACH PROBLEM, A LEAST-SQUARES  FIT  IS
c         MADE OF LOG10(ACTUAL ERROR) VS LOG10(TOL) AND USED TO ESTIMATE
c         WHAT THE VARIOUS COST STATISTICS WOULD BE FOR AN ACTUAL  ERROR
c         OF 10**N.  THIS IS ACHIEVED BY INTERPOLATION, FOR THOSE N SUCH
c         THAT 10**N LIES WITHIN THE RANGE OF ACCURACIES  ACHIEVED  WITH
c         THE USER-SPECIFIED TOLERANCES.
c     0   NO NORMALIZED STATISTICS
c     1   NORMALIZED STATISTICS ARE PRODUCED TAKING THE  'ACTUAL  ERROR'
c         USED IN THE LEAST SQUARES FIT TO BE THE ENDPOINT GLOBAL ERROR.
c     2   NORMALIZED STATISTICS ARE PRODUCED TAKING  'ACTUAL  ERROR'  AS
c         THE  MAXIMUM  GLOBAL ERROR OVER THE RANGE.  N.B.  IN THIS CASE
c         OPT MUST BE AT LEAST 2.
c
c NRMTYP  ONE OF 1, 2  OR 3, SELECTS THE NORM USED IN ASSESSING THE SIZE
c         OF LOCAL AND GLOBAL ERRORS. IT SHOULD BE CHOSEN BY THE USER TO
c         AGREE WITH THE NORM USED IN SOLVER. WE OFFER:
c     1   MAX-NORM.
c     2   2-NORM (EUCLIDEAN NORM).
c     3   R.M.S. NORM, THAT IS (2-NORM OF X)/SQRT(N) FOR AN N-VECTOR X.
c
c TOL     (INPUT) REAL ARRAY, HOLDS LIST OF UP TO 10  TOLERANCES  TO  BE
c         USED,  IN  STRICTLY  DECREASING  ORDER,  WITH 0 AS TERMINATOR.
c         EACH PROBLEM IS INTEGRATED AT EACH TOLERANCE IN TURN.
c         EXAMPLE:  IN CALLING PROGRAM
c                   REAL TOL(11)
c                   DATA TOL/1E-1,1E-3,1E-5,1E-7,7*0D0/
c         REQUESTS THE FOUR TOLERANCES .1, .001, .00001, .0000001.
c
c IDLIST  (INPUT) integer ARRAY, HOLDS LIST OF GROUPS OF  PROBLEMS,  AND
c         SPECIFIES  FOR  EACH  ONE  WHETHER  IT  IS TO BE INTEGRATED IN
c         SCALED OR UNSCALED  FORM  (SEE  GENERAL  NOTES  ABOVE).   EACH
c         PROBLEM  IS SPECIFIED BY A NUMERIC CODE, 11 TO 14 FOR PROBLEMS
c         A1 TO A4, 21 TO 25 FOR B1 TO B5  ETC.   A  ZERO  TERMINATES  A
c         GROUP AND TWO ZEROS TERMINATE THE LIST OF GROUPS.
c         IF THE PROBLEM CODE IS GIVEN A NEGATIVE SIGN,  THE  SYSTEM  IS
c         INTEGRATED  IN  UNSCALED  FORM;  IF A POSITIVE SIGN, IN SCALED
c         FORM.
c         EXAMPLE:  IN CALLING PROGRAM
c                   integer IDLIST(7)
c                   DATA IDLIST/11,22,0,-31,-51,0,0/
c         SPECIFIES GROUP 1 CONSISTING OF PROBLEMS A1,B2 AND GROUP 2  OF
c         PROBLEMS  C1,E1.  THE FIRST TWO ARE TO BE SOLVED IN THE SCALED
c         FORM AND THE LAST TWO  IN  UNSCALED  FORM.
c
c         THE TOTAL LENGTH OF THE LIST INCLUDING ZEROS MUST BE  AT  MOST
c         60 ITEMS.
c
c FLAG
c         (OUTPUT) REAL.  A NONZERO VALUE INDICATES  THAT  THE  CALL  TO
c         NSDTST  WAS  ABORTED BECAUSE OF ARGUMENT ERRORS, IN WHICH CASE
c         THE VALUES OF THE DECIMAL DIGITS OF FLAG INDICATE THE ERROR(S)
c         THAT HAVE OCCURRED, AS FOLLOWS:
c           1:  OPT INVALID.
c           2:  NORMEF INVALID.
c           3:  NORMEF = 2 WAS REQUESTED WITH OPT = 1.
c           4:  A NEGATIVE  TOLERANCE  WAS  SUPPLIED,  OR THE  LIST  OF
c               TOLERANCES WAS NOT IN DECREASING ORDER.
c           5:  THE LIST OF TOLERANCES WAS EMPTY OR NOT TERMINATED BY A
c               ZERO.
c           6:  AN INVALID PROBLEM-ID WAS FOUND IN IDLIST.
c           7:  THE LIST  OF  GROUPS  IN  IDLIST  IS  EMPTY  OR  IS NOT
c               TERMINATED  BY  TWO  ZEROS OR HAS MORE THAN THE MAXIMUM
c               ALLOWED NUMBER (6) OF GROUPS.
c           8:  NRMTYP INVALID.
c         EG.  A VALUE FLAG = 0.245E 03 INDICATES THAT ERRORS 2, 4 AND 5
c         IN  THE  ABOVE  LIST  HAVE  OCCURRED.  ITS VALUE IF NONZERO IS
c         PRINTED BY NSDTST ANYWAY, BUT FLAG IS MEANT TO BE INSPECTED IF
c         FURTHER  ACTION  OF  THE  MAIN PROGRAM DEPENDS ON A SUCCESSFUL
c         CALL TO NSDTST.
c
c
c 3. INTERFACE ROUTINE METHOD
c  
c
c THIS INVOKES THE CODE BEING TESTED, CALL IT SOLVER.  THE SPECIFICATION
c IS
c         SUBROUTINE METHOD(N,X,Y,XEND,TOL,HMAX,HSTART)
c         integer N
c         double precision X,Y(N),XEND,TOL,HMAX,HSTART
c         EXTERNAL FCN
c
c METHOD IS TO BE WRITTEN BY THE USER AS A SIMPLE INTEGRATOR TO  ADVANCE
c THE  SOLUTION OF N DIFFERENTIAL EQUATIONS FROM THE INITIAL VALUES HELD
c IN X,Y UP TO XEND, WITH AN UNWEIGHTED ABSOLUTE ERROR CONTROL  OF  TOL.
c HMAX  IS  A  RECOMMENDED  MAXIMUM STEPSIZE AND HSTART IS A RECOMMENDED
c INITIAL STEPSIZE.  IF SOLVER CAN MAKE USE OF THESE TWO PARAMETERS, THE
c STATISTICS WILL PROBABLY BE MORE FAVORABLE AND RELIABLE, BUT THEIR USE
c IS NOT CRUCIAL.
c
c THE DERIVATIVES OF THE PROBLEM ARE  COMPUTED  BY  PACKAGE ROUTINE FCN.
c THUS FCN WILL BE AN ARGUMENT TO  SOLVER, AND MUST BE DECLARED EXTERNAL
c IN METHOD.
c
c METHOD SHOULD CALL SOLVER IN ONE-STEP MODE  SO  THAT  A  CALL  TO  THE
c PACKAGE  ROUTINE  STATS  CAN  BE  MADE AFTER EACH SUCCESSFUL STEP.  IF
c SOLVER DOES NOT HAVE THIS FACILITY, SOLVER MUST HAVE A CALL  TO  STATS
c INSERTED AT THE APPROPRIATE POINT IN THE CODE.
c
c SOME  CALLS  TO  METHOD  ARE  INTENDED  TO  BE  ABORTED  AFTER  A  FEW
c INTEGRATION  STEPS  BY  THE  STATS CALL SETTING X = XEND.  THUS A TEST
c SHOULD BE MADE AFTER EACH CALL TO STATS, OF THE FORM
c         IF STATS HAS SET X = XEND THEN EXIT.
c
c NB:  IF THE ACTUAL X  ARGUMENT  TO  STATS  IS  DIFFERENT  FROM  THE  X
c ARGUMENT  OF METHOD (WHICH MAY BE NECESSARY WITH SOME SOLVERS), ENSURE
c THAT THE X ARGUMENT OF METHOD IS SET TO XEND  BEFORE  EXIT,  ELSE  THE
c PACKAGE WILL REPORT 'METHOD FAILED TO START'.
c
c THE ALGORITHM FOR METHOD SHOULD THUS BE OF THE FORM:
c - DECLARE ALL ARGUMENTS AND WORKSPACE EXPECTED BY SOLVER
c - SET APPROPRIATE OPTIONS  INCLUDING  ABSOLUTE  ERROR  CONTROL  AND
c    ONE-STEP MODE
c - INITIALIZE EXTRA ARGUMENTS IF REQUIRED
c - FOR EACH SUCCESSFUL STEP DO
c    - CALL SOLVER( ...  ,FCN, ...  )
c      EXIT IF SOLVER IS IN TROUBLE.
c    - SET X,Y TO THE JUST COMPUTED MESHPOINT X AND SOLUTION VECTOR Y
c    - SET ERRBND TO THE BOUND THAT IS  SATISFIED  BY ||ERREST||, AND
c      HENCE IS INTENDED TO BE SATISFIED BY ||LE||, AT THIS STEP.
c    - SET ERREST  TO THE  LOCAL ERROR ESTIMATE  VECTOR Y-Y*  DEFINED
c      ABOVE
c
c      (SEE   [3]   FOR   DISCUSSION  AND  NOTE  THAT X,Y ARE IGNORED
c      UNLESS  OPT.GE.2,  ERRBND   IS  IGNORED  UNLESS  OPT.GE.3, AND
c      ERREST IS IGNORED UNLESS OPT.GE.4.)
c
c    - CALL STATS(X,Y,ERRBND,ERREST)
c    - EXIT IF X .GE.  XEND.
c - ENDLOOP
c
c
c ON NORMAL EXIT X,Y MUST HOLD XEND AND THE SOLUTION AT XEND.   ON  EXIT
c BECAUSE  SOLVER  WAS  IN TROUBLE, X MUST HOLD THE FINAL POINT REACHED.
c ON AN EXIT FORCED BY STATS, X MUST HOLD XEND.
c
c
c 4. CONTROLLING THE DESTINATION OF OUTPUT
c  
c
c THE UNIT NUMBER ON WHICH THE PACKAGE WRITES ITS OUTPUT  IS  SET  BY  A
c CALL  TO ONE OF THE PACKAGE ROUTINES, AND YOU CAN FIND OUT WHAT IT IS,
c BY PUTTING THE STATEMENT
c
c       IOUT = CONST ( 3 )
c
c IN YOUR MAIN PROGRAM.  PROBABLY OUTPUT WILL DEFAULT TO YOUR  TERMINAL,
c WHICH  IS  GOOD  FOR DEBUGGING.  FOR MORE SERIOUS WORK YOU MAY WANT TO
c SEND OUTPUT TO A FILE.  THE STATEMENTS
c
c       IOUT = CONST ( 3 )
c       OPEN(IOUT, FILE=FILENAME, OTHER OPTIONS..  )
c
c WILL DO THIS FOR YOU, ASSUMING YOUR  FORTRAN   I/O  IS CONSISTENT WITH
c THE  1977 STANDARD.
c
c
c
c 5. THE ROUTINE FCN
c  
c
c THE SPECIFICATION OF FCN IS
c         SUBROUTINE FCN(X,Y,YP)
c         double precision X,Y(51),YP(51)
c
c ON ENTRY X HOLDS THE INDEPENDENT VARIABLE AND Y HOLDS  THE  VECTOR  OF
c DEPENDENT  VARIABLES.   ON EXIT YP HOLDS THE VECTOR OF DERIVATIVES FOR
c THE PROBLEM BEING SOLVED (SELECTED BY A SWITCH IN COMMON).
c
c
c 6. THE DUMMY NSDTST FOR DEBUGGING
c  
c
c TO THE USER:
c
c THIS WILL PROBABLY BE IMPLEMENTED  AT  YOUR  SITE  AS  A  SOURCE  FILE
c CONTAINING  CUT-DOWN  VERSIONS  OF NSDTST AND STATS (AND OTHER PACKAGE
c ROUTINES OF NO CONCERN TO THE  USER).   THIS  FILE  MAKES  A  COMPLETE
c PROGRAM WHEN COMBINED WITH THE NSPROB FILE AND THE USER'S MAIN PROGRAM
c AND METHOD (AND OF COURSE SOLVER).  THE  CUT-DOWN  ROUTINES  HAVE  THE
c SAME CALLING SEQUENCE AS THE PROPER ONES.
c
c THE RESULTING PROGRAM USES METHOD TO SOLVE THE FIRST PROBLEM SPECIFIED
c IN IDLIST, AT THE FIRST TOLERANCE SPECIFIED IN TOL.  IT WILL PRINT OUT
c THE VALUES OF THE ARGUMENTS PASSED BY METHOD TO STATS AND ALSO THE  LU
c DECOMP  COUNTER  NLUD,  FOR  5 STEPS, AND THEN SET X = XEND.  THE USER
c SHOULD CHECK THAT THE VALUES OF X, Y, ERREST, ERRBND LOOK RIGHT;  THAT
c X = XEND  FORCES  TERMINATION  AS  IT  SHOULD;  AND THAT NLUD IS BEING
c COUNTED UP CORRECTLY.
c
c FEEL FREE TO MODIFY THESE ROUTINES TO WORK INTERACTIVELY.
c
c TO THE PERSON IMPLEMENTING THE PACKAGE:
c
c PLEASE MODIFY THESE ROUTINES TO MATCH THE USER ENVIRONMENT.
c
c
c
c 7. IMPLEMENTATION NOTES
c  
c
c  7.1.  MACHINE-DEPENDENT CONSTANTS
c
c     THESE  ARE   ISOLATED  IN  THE  ROUTINE   CONST   WHICH   HAS  THE
c     SPECIFICATION   REAL  FUNCTION CONST ( I ).  YOU MUST SET THE ARRAY C
c     AND THE STRING MCNAME IN THE DATA STATEMENT:
c
c     C(1)   APPROXIMATELY  THE  DOUBLE  PRECISION  UNIT  ROUNDOFF, USED
c            IN STATS AND TRUE.
c     C(2)   A NUMBER NEAR THE UNDERFLOW THRESHOLD, USED IN TRUE.
c     C(3)   THE  STANDARD  OUTPUT  UNIT NUMBER IOUT, USED IN NSDTST AND
c            TRUE.   WE SUGGEST OUTPUT BE TO THE TERMINAL BY DEFAULT.
c     C(4)   TSTTIM, USED IN CNTROL (SEE CLOCK ROUTINE).
c     MCNAME TITLING  INFORMATION  FOR PRINTOUT, GIVING THE  NAME OF THE
c            COMPUTER AND OPERATING SYSTEM.
c
c     IN ADDITION, A CALL OF CONST(0) (EXECUTED NEAR THE TOP OF  NSDTST)
c     IS  INTENDED  TO  INVOKE  CALLS  TO  SYSTEM  ROUTINES  TO SUPPRESS
c     UNDERFLOW  REPORTING  (WHICH  MAY  SPOIL  THE  APPEARANCE  OF  THE
c     OUTPUT), ETC.
c
c     IT MAY BE CONVENIENT  TO  ALLOW  IOUT  (C(3)  ABOVE)  TO BE SET BY
c     INTERACTION WITH THE USER AT THIS POINT.
c
c  7.2.  CLOCK ROUTINE
c
c     IF   IT  IS  DECIDED  TO  IMPLEMENT  THE  TIMING  FACILITIES,  THE
c     IMPLEMENTER  SHOULD  PROVIDE  A  TIMING ROUTINE  WHICH  CALLS  THE
c     SYSTEM CLOCK AND HAS THE SPECIFICATION
c          REAL FUNCTION CLOCK(S)
c          REAL(S)
c     IT SHOULD BE SUCH THAT IT IS 'RESET TO ZERO' BY THE STATEMENT
c          S = CLOCK(0.0)
c     AND (AS LONG AS S IS LEFT ALONE) CAN THEN BE 'READ'  AS  OFTEN  AS
c     DESIRED BY STATEMENTS LIKE
c          TIME = CLOCK(S)
c     WHICH SETS TIME TO THE NUMBER OF SECONDS OF PROCESSOR  TIME  SINCE
c     CLOCK WAS 'RESET'.
c
c     THE LARGER IS  THE  VALUE  OF  TSTTIM  (IE.   CONST(4))  THE  MORE
c     ACCURATE,  AND EXPENSIVE, IS THE TIMING PROCESS.  IT SHOULD BE SET
c     TO A VALUE REFLECTING THE SPEED OF THE HARDWARE AND THE RESOLUTION
c     OF  THE  SYSTEM CLOCK.  WE CANNOT GIVE MUCH GUIDANCE HERE, AND OUR
c     EXPERIENCE IS THAT TIMINGS INEVITABLY VARY SIGNIFICANTLY FROM  RUN
c     TO RUN ON A TIME-SHARED COMPUTING SYSTEM.
c
c     IF TIMING IS LEFT UNIMPLEMENTED, GIVE  TSTTIM  THE  VALUE ZERO AND
c     LEAVE THE TIMING DATA IN IVALU AS ALL ZERO  TO CAUSE ALL VALUES OF
c     TIMING STATISTICS TO BE PRINTED OUT AS ZERO.
c
c  7.3.  THE TIMING DATA IN IVALU
c
c     ROUTINE IVALU CONTAINS VALUES OF THE  QUANTITY  FCNTIM  FOR   EACH
c     PROBLEM:  THESE ARE THE COST OF ONE CALL TO  FCN  AS  MEASURED  BY
c     CLOCK, AND ARE USED IN COMPUTING THE "OVERHEAD"   STATISTICS. THEY
c     SHOULD BE  RECOMPUTED  FOR  ANOTHER  MACHINE.  THE UTILITY PROGRAM
c     NSGTIM  PROVIDED  WITH  THE PACKAGE, WHEN SUPPLIED  WITH  A  CLOCK
c     ROUTINE, CAN EITHER PRODUCE  A COMPLETE REVISED IVALU FILE, OR FOR
c     SELECTED PROBLEMS WILL PRODUCE BLOCKS OF OUTPUT OF THE FORM
c
c     C PROBLEM XX
c           fcntim = ...
c
c     SUITABLE FOR INCLUSION IN THE TEXT OF IVALU.
c
c
c  7.4.  ADDING EXTRA PROBLEMS
c
c     SAY YOU WISH TO ADD THREE EXTRA PROBLEMS TO CLASS B  OF  THE  SET.
c     THEY  WILL  THEN  BE  CALLED  B6,  B7  AND B8 (FOR THE SAKE OF THE
c     CHECKING ROUTINE PARCHK THEY  MUST  FOLLOW  CONSECUTIVELY  ON  THE
c     EXISTING  PROBLEM-IDS).   THEIR NUMERIC CODES WHICH YOU SPECIFY IN
c     THE IDLIST ARGUMENT OF NSDTST WILL THEN BE 26, 27, 27.   YOU  NEED
c     TO BE AWARE THAT THE INTERNAL CODE, PUT IN VARIABLE ID AND USED IN
c     FCN, IVALU AND EVALU TO SELECT THE CORRECT SECTION  OF ROUTINE
c     TO EXECUTE, IS 10 LESS THAN THIS, IE.  16, 17 OR 17.
c
c     THE  STEPS  INVOLVED  ARE:
c     A)  CODE THE  DEFINITION  OF   THE   DIFFERENTIAL   EQUATIONS   AT
c         THE APPROPRIATE PLACE IN FCN.  CHANGE THE COMPUTED GOTO AT THE
c         HEAD  OF  FCN SO THAT THE VALUE ID = 16, 17 OR 18 GIVES A JUMP
c         TO  THE CORRECT PLACE.
c     B)  CODE THE  INITIAL  VALUES,  "TRUE"  FINAL   VALUES  AND  OTHER
c         DATA INTO  THE  APPROPRIATE  PLACES  IN IVALU  AND  EVALU IN A
c         SIMILAR  WAY.   THE  TRUE  FINAL  VALUES  SHOULD  PROBABLY  BE
c         COMPUTED   BY   AN  INTEGRATOR   USING   HIGHER   THAN  DOUBLE
c         PRECISION, BUT THE ONLY CONSEQUENCE OF SLIGHT  INACCURACIES IS
c         TO  AFFECT  THE  END  PT   GLB ERR  STATISTIC   AT   STRINGENT
c         TOLERANCES.  AT THIS STAGE IGNORE THE  WEIGHTS  W(I)  AND  THE
c         TIMING DATA FCNTIM.
c     C)  IN THE ARGUMENT-CHECKING  ROUTINE   PARCHK  CHANGE  THE   DATA
c         STATEMENT  WHICH DEFINES ARRAY NSYSTM,  TO INDICATE THAT CLASS
c         B NOW HAS 8 MEMBERS.  (IE.   CHANGE ITS SECOND ELEMENT FROM  5
c         TO 7.)
c     D)  RUN  THE UTILITY  PROGRAM  NSGWT.F ON THE TAPE TO COMPUTE  THE
c         VALUES  OF  THE  WEIGHTS  W(I).   SIMILARLY  RUN  NSGTIM.F  TO
c         DETERMINE FCNTIM FOR YOUR PROBLEMS.
c
c ADDING  AN  ENTIRE  NEW  PROBLEM CLASS IS  NO  MORE  DIFFICULT.   NOTE
c THAT  IT  INVOLVES  INCREASING  THE   VALUE   OF   NCLASS  IN THE DATA
c STATEMENT AND THE LENGTH OF NSYSTM IN THE  DIMENSIONING   STATEMENT IN
c PARCHK; ALSO CHECK THE STRING IDCLAS  IN  NSDTST HAS ENOUGH LETTERS IN
c IT.
c
c  7.5.  OTHER STATISTICS TO PRINT
c
c STATISTICS  WHICH ARE GATHERED BUT DO  NOT  APPEAR   IN   THE   OUTPUT
c TABLES  INCLUDE   NSTART,  NSTL  AND  TRUTIM.  THEY ARE DEFINED IN THE
c DESCRIPTION   OF   COMMON  /NSCOM3/  BELOW.   NSTART   ASSESSES    THE
c EFFICIENCY  OF  THE  STARTING  PHASE  OF  A CODE AND MAY BE OF GENERAL
c INTEREST.  TRUTIM IS OF USE  IF  YOU  ARE  TROUBLED  BY  THE OVERHEADS
c OF CALLS TO TRUE WITH OPT  >=  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. ROUTINES 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 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
      implicit none

      double precision FLAG
      character*80      TITLE

      double precision             TOL(11)
      integer           IDLIST(60), OPTION(10)

      double precision  ERRTOL, XFIN, XTRUE
      double precision  GEMX, GEND, LEMXSC, OVHD, TIME, TRUTIM
      integer           ID, IOUT, IWT, NBAD, NDCV, NFCN, NRMTYP, NSTART,
     *                  NSTL, NSTP, NTRU, OPT, XTRAP

      double precision    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

      double precision 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)

      double precision CONST, RATIO
      EXTERNAL          CONST, RATIO

      EXTERNAL          CNTROL, EFSTAT, LSQFIT, PARCHK, PLOT

      INTRINSIC         CHAR, DBLE, IABS, ISIGN

      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

      data IDCLAS/'ABCDEFGHIJ'/
      data BL/' '/
      data BIG/1.E20/
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 I = 1, 32
         ICH = CONST(-I)
         MCNAME(I:I) = CHAR(ICH)
      end do
      IOUT = CONST(3)

      WRITE (IOUT,FMT=99999) OPT, NORMEF, NRMTYP, MCNAME

      CALL PARCHK(OPT,NORMEF,NRMTYP,TOL,IDLIST,NTOL,NGRP,GRPLST,LGTOL,
     *            FLAG)

      if ( FLAG .ne. 0.0D+00 ) then
        WRITE (IOUT,FMT=99998) FLAG
        RETURN
      end if
c
c  INITIALIZE OVERALL- AND GROUPS-SUMMARY STATISTICS.
c  
      OSTIME = 0.0D+00
      OSOVHD = 0.0D+00
      OSNFCN = 0
      OSNSTP = 0
      OSNTRU = 0
      OSLEMX = 0.0D+00
      OSNDCV = 0
      OSNBAD = 0
      DO I = 1, NTOL
         GSTIME(I) = 0.0D+00
         GSOVHD(I) = 0.0D+00
         GSNFCN(I) = 0
         GSNSTP(I) = 0
         GSNTRU(I) = 0
         GSLEMX(I) = 0.0D+00
         GSNDCV(I) = 0
         GSNBAD(I) = 0
      end do
c  
c  LOOP OVER GROUPS OF PROBLEMS
c  
      DO 300 KGRP = 1, NGRP
c  
c  OUTPUT HEADING, ON NEW PAGE FOR GROUPS AFTER FIRST.
c  SELECT GROUP OF DIFFERENTIAL EQUATIONS.
c  GET number OF SYSTEMS IN THIS GROUP, and 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

         NSYST = GRPLST(1,KGRP)
         IDSUB = GRPLST(2,KGRP)

         DO I = 1, NTOL
            PSTIME(I) = 0.0D+00
            PSOVHD(I) = 0.0D+00
            PSNFCN(I) = 0
            PSNSTP(I) = 0
            PSNTRU(I) = 0
            PSLEMX(I) = 0.0D+00
            PSNDCV(I) = 0
            PSNBAD(I) = 0
            PSGEMX(I) = 0.0D+00
            PSGEND(I) = 0.0D+00
         end do
c
c  LOOP OVER PROBLEMS WITHIN A GROUP
c  
         DO 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  LOOP OVER TOLERANCES FOR ONE PROBLEM
c  
            DO 220 KTOL = 1, NTOL
c  
c  CALL PLOT TO INITIALIZE LOCAL-ERROR SCATTER DIAGRAM IF OPT=4.
c  CALL CNTROL TO ORGANIZE THE COLLECTION OF 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 LOCAL OR GLOBAL 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
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 IF OPT=4
c  
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 
               TOLK = TOL(KTOL)
               ERRTOL = DBLE(TOLK)
               if (OPT.EQ.4) CALL PLOT(0.0D+00,0.0D+00,0)

               CALL CNTROL(CMPLET,INDG1,INDL1)

               ERFLGE(KTOL) = CMPLET .LT. 0 .AND. CMPLET .GT. -4
               ERFLG1(KTOL) = CMPLET .LT. 1 .AND. CMPLET .GT. -4
               GENDSC = BIG

               if ( .not. ERFLGE(KTOL) ) then
                 GENDSC = GEND/TOLK
                 LGGEND(KTOL) = DLOG10(DMAX1(GEND,0.01D+00*TOLK))
               end if

               GEMXSC = GEMX/TOLK
               FDECEV = RATIO(NDCV,NTRU)
               FBADEC = RATIO(NBAD,NTRU)

               if ( CMPLET .ne. -3 ) then

                 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

               end if

               if (CMPLET.EQ.-4) WRITE (IOUT,FMT=99989)
               if (CMPLET.EQ.-3) WRITE (IOUT,FMT=99988) LGTOL(KTOL)

               if (CMPLET.EQ.-2) WRITE (IOUT,FMT=99987) XTRUE, INDG1,
     *             INDL1, XFIN

               if (CMPLET.EQ.-1) WRITE (IOUT,FMT=99986) XFIN

               if (CMPLET.EQ.0) WRITE (IOUT,FMT=99985) XTRUE, INDG1,
     *             INDL1

               if (OPT.EQ.4) THEN

                  WRITE (IOUT,FMT=99984) XTRAP

                  CALL PLOT(0.0D+00,0.0D+00,2)
               END IF
c
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
               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) = DMAX1(PSGEND(KTOL),GENDSC)

               if ( 2 .le. OPT ) then
                 PSGEMX(KTOL) = DMAX1(PSGEMX(KTOL),GEMXSC)
                 LGGEMX(KTOL) = DLOG10(DMAX1(GEMX,0.01D+00*TOLK))
               end if

               if (3 .le. OPT ) then
                 PSNTRU(KTOL) = PSNTRU(KTOL) + NTRU
                 PSLEMX(KTOL) = DMAX1(PSLEMX(KTOL),LEMXSC)
                 PSNDCV(KTOL) = PSNDCV(KTOL) + NDCV
                 PSNBAD(KTOL) = PSNBAD(KTOL) + NBAD
               end if

  180          CONTINUE

               if (NORMEF .ne. 0) then
                 NSTIME(KTOL) = TIME
                 NSOVHD(KTOL) = OVHD
                 NSNFCN(KTOL) = NFCN
                 NSNSTP(KTOL) = NSTP
               end if
c  
c  END OF LOOP OVER TOLERANCES FOR ONE PROBLEM
c  
  220       CONTINUE
c
c  SMOOTHNESS AND NORMALIZED EFFICIENCY CALCULATIONS BEGIN
c  
            WRITE (IOUT,FMT=99983)

            WRITE (IOUT,FMT=99982)

            CALL LSQFIT(LGTOL,LGGEND,ERFLGE,NTOL,NOK,C,E,RES)

            CTEN = 10.0D+00**C
            if (NOK.LE.2) WRITE (IOUT,FMT=99981) NOK

            if (NOK.GT.2) WRITE (IOUT,FMT=99980) CTEN, E, RES, NOK

            if ( 2 .le. OPT ) then

              WRITE (IOUT,FMT=99979)

              CALL LSQFIT(LGTOL,LGGEMX,ERFLG1,NTOL,NOK1,C1,E1,RES1)

              CTEN1 = 10.0D+00**C1
              if (NOK1.LE.2) WRITE (IOUT,FMT=99981) NOK1
              if (NOK1.GT.2) WRITE (IOUT,FMT=99980) CTEN1, E1, RES1, NOK1

            end if

            if (NORMEF.EQ.1) then
              CALL EFSTAT(C,E,LGTOL,NTOL,NOK,ERFLGE,
     *                                   'ENDPOINT',IOUT,NSTIME,NSOVHD,
     *                                   NSNFCN,NSNSTP)
            end if

            if (NORMEF.EQ.2) then
              CALL EFSTAT(C1,E1,LGTOL,NTOL,NOK1,ERFLG1,
     *                                   'MAXIMUM ',IOUT,NSTIME,NSOVHD,
     *                                   NSNFCN,NSNSTP)
            end if
c  
c  SMOOTHNESS AND NORMALIZED EFFICIENCY CALCULATIONS END
c
c  END OF LOOP OVER PROBLEMS IN A GROUP.
c  
         end do
c
c  OUTPUT PROBLEMS-SUMMARY STATISTICS
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))

            if (OPT.EQ.1) WRITE (IOUT,FMT=99991) LGTOL(KTOL),
     *          PSTIME(KTOL), PSOVHD(KTOL), PSNFCN(KTOL), PSNSTP(KTOL),
     *          PSGEND(KTOL)

            if (OPT.EQ.2) WRITE (IOUT,FMT=99991) LGTOL(KTOL),
     *          PSTIME(KTOL), PSOVHD(KTOL), PSNFCN(KTOL), PSNSTP(KTOL),
     *          PSGEND(KTOL), PSGEMX(KTOL)

            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

            if (OPT.GE.3 .AND. PSNSTP(KTOL).NE.PSNTRU(KTOL))
     *          WRITE (IOUT,FMT=99990) PSNTRU(KTOL)
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)

            if ( 3 .le. OPT ) then
              GSNTRU(KTOL) = GSNTRU(KTOL) + PSNTRU(KTOL)
              GSLEMX(KTOL) = DMAX1(GSLEMX(KTOL),PSLEMX(KTOL))
              GSNDCV(KTOL) = GSNDCV(KTOL) + PSNDCV(KTOL)
              GSNBAD(KTOL) = GSNBAD(KTOL) + PSNBAD(KTOL)
            end if

  280   continue
c
c  END OF LOOP OVER GROUPS
c  
  300 CONTINUE
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 I = 1, NTOL
         WRITE (IOUT,FMT=99975) LGTOL(I), GSTIME(I), GSOVHD(I),
     *     GSNFCN(I), GSNSTP(I)
      end do

      GO TO 380
  340 DO 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

         if (GSNSTP(I).NE.GSNTRU(I)) WRITE (IOUT,FMT=99990) GSNTRU(I)
      end do

  380 CONTINUE
c
c  COMPUTE OVERALL-SUMMARY STATISTICS.
c  
      DO I = 1, NTOL

         OSTIME = OSTIME + GSTIME(I)
         OSOVHD = OSOVHD + GSOVHD(I)
         OSNFCN = OSNFCN + GSNFCN(I)
         OSNSTP = OSNSTP + GSNSTP(I)

         if ( 3 .le. OPT ) then
           OSNTRU = OSNTRU + GSNTRU(I)
           OSNDCV = OSNDCV + GSNDCV(I)
           OSNBAD = OSNBAD + GSNBAD(I)
           OSLEMX = DMAX1(OSLEMX,GSLEMX(I))
         end if

      end do

      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

      if (OPT.GE.3) WRITE (IOUT,FMT=99974) OSTIME, OSOVHD, OSNFCN,
     *    OSNSTP, OSLEMX, FDECEV, FBADEC

      RETURN
c
99999 FORMAT ('NONSTIFF DETEST PACKAGE    OPTION=',I2,', NORMEF=',I2,
     *       ', NRMTYP=',I2,19X,'ON ',A,//)
99998 FORMAT ('PARAMETER ERRORS AS SHOWN BY FLAG=',E15.8,/' ',49('*')
     *       ,//)
99997 FORMAT (' ')
99996 FORMAT ('GROUP',I3,18X,A)
99995 FORMAT (/A3,I1,'   (SCALED)',/)
99994 FORMAT (/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 (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 (20X,
     *      '***** UNEXPECTED FAILURE OF METHOD WHILE BEING TIMED *****'
     *       ,/)
99988 FORMAT (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 (/17X,'SMOOTHNESS FIT OF LOG10(ERROR) VS LOG10(TOL)')
99982 FORMAT (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 (17X,'MAXIMUM  GLOBAL ERROR')
99978 FORMAT (/'SUMMARY OVER GROUP',I3)
99977 FORMAT ('SUMMARY 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 (6X,F6.2,2X,2F7.3,1X,2I8,1X,3F11.3)
99974 FORMAT (5X,'OVERALL',/6X,'SUMMARY',2X,2F7.3,1X,2I8,1X,3F11.3)
      end
      function clock ( s )

c***********************************************************************
c
cc CLOCK is meant to return the elapsed CPU seconds.
c
c  Modified:
c
c    15 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
c    Input, double precision S, the previous clock reading.
c
c    Output, double precision CLOCK, the elapsed CPU time, in seconds,
c    since S.
c
      implicit none

      double precision clock
      double precision s

      clock = 0.0D+00

      return
      end
      subroutine cntrol ( cmplet, indg1, indl1 )

c***********************************************************************
c
cc CNTROL calls METHOD to get statistics for one problem and tolerance.
c
c  Discussion:
c
c    This routine 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  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
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  COMMON AREAS
c
      implicit none

      integer           CMPLET, INDG1, INDL1

      double precision  ERRTOL, HSTART, XEND, XFIN, XTRUE
      double precision  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

      double precision  WT(51)

      double precision  DUMMY, HINIT, HMAX, X, XSTART
      double precision              FCNTIM, S, TIMCUM, TSTTIM
      integer           COUNT, I
      LOGICAL           NOSTRT, OKMETH, TIMERR

      double precision  Y(51), YEND(51), YSTART(51)

      double precision              CLOCK, CONST, DIFNRM
      EXTERNAL          CLOCK, CONST, DIFNRM

      EXTERNAL          EVALU, IVALU, METHOD, STATS

      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
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 ROUTINE 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  
      CALL IVALU(N,XSTART,XEND,HINIT,HMAX,YSTART,FCNTIM,WT,IWT,ID)

      N1 = N
      ID1 = ID
      IWT1 = IWT
      X = XSTART
      DO I = 1, N
        y(I) = YSTART(I)
      end do

      IFLAG = 1
      CALL STATS(X,Y,DUMMY,Y)

      NFCN1 = 0
      NSTP = 0
      IFLAG = 2

      CALL METHOD(N,X,Y,XEND,ERRTOL,HMAX,HINIT)

      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.0D+00
      NDCV = 0
      NBAD = 0
      GEMX = 0.0D+00
      TRUTIM = 0.0D+00
      NTRU = 0

      if (NOSTRT) GO TO 180

      X = XSTART
      DO I = 1, N
        y(I) = YSTART(I)
      end do
      IFLAG = 3
      S = CLOCK(0.0D+00)

      CALL METHOD(N,X,Y,XEND,ERRTOL,HMAX,HSTART)

      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.5D+00*TSTTIM) GO TO 120
      COUNT = 0
      IFLAG = 0
      S = CLOCK(0.0D+00)
c  
c  LOOP TILL 'TSTTIM' TIME UNITS HAVE ELAPSED:
c  
   60 CONTINUE

      X = XSTART
      DO I = 1, N
        y(I) = YSTART(I)
      end do
      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 ) then
        GO TO 60
      end if

  100 if (COUNT.GE.1) TIME = TIMCUM/dble(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 - dble(NFCN)*FCNTIM
      CALL EVALU(YEND,N,WT,IWT,ID)
      GEND = DIFNRM(YEND,Y,N)

      if (TIMERR) GO TO 200
c  
c  SET THE OUTPUT VALUE OF CMPLET, INDG1 AND INDL1.
c  
      CMPLET = 1
      if (INDL.LT.0 .OR. INDG.LT.0) then
        CMPLET = 0
      end if

  140 INDG1 = INDG
      INDL1 = INDL
      RETURN
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 = 1.0D+20
      OVHD = 1.0D+20
      GEND = 1.0D+20
      indg1 = indg
      indl1 = indl
      return
c
c  METHOD FAILED TO START
c  
  180 CMPLET = -3
      NFCN = 0
      TIME = 1.0D+20
      OVHD = 1.0D+20
      GEND = 1.0D+20
      indg1 = indg
      indl1 = indl
      return
c  
c  INTEGRATION FAILED IN TIMING LOOP
c  
  200 CMPLET = -4
      indg1 = indg
      indl1 = indl
      return
      end
      function const ( i )

c***********************************************************************
c
cc CONST returns machine dependent constants.
c
c  Modified:
c
c    15 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
c    Input, integer I, indicates the value to be returned.  If
c    -32 <= I <= -1, then character (-I) of a 32 character string
c    defining the machine name is to be returned.  Otherwise,
c    1, unit roundoff;
c    2, underflow threshold.  setting this to the smallest representable
c       real number seems to have undesirable effects!  So we just make
c       if a "small" value.
c    3, standard output unit;
c    4, value of TSTTIM.
c
c    Output, double precision CONST, the requested value.
c
      implicit none

      double precision const
      integer i
      character ( len = 32 ) :: mcname = 'Dec Alpha'

      if ( i .eq. 1 ) then
        const = 2.220446049250313D-16
      else if ( i .eq. 2 ) then
c       const = 1.225073858507201D-308
        const = 1.0D-30
      else if ( i .eq. 3 ) then
        const = dble ( 6 )
      else if ( i .eq. 4 ) then
        const = 0.0D+00
      else if ( i .lt. 0 ) then
        const = dble ( ichar ( mcname(-i:-i) ) )
      else
        const = 0.0D+00
      end if

      return
      end
      function difnrm ( a, b, n )

c***********************************************************************
c
cc DIFNRM computes the norm of the difference of two vectors.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
c     NRMTYP=1,2,3 CHOOSES MAX-NORM, 2-NORM, R.M.S.-NORM.
c  
      implicit none

      integer n

      double precision A(N)
      double precision B(N)
      double precision difnrm
      double precision     ERRTOL
      integer              ID, IOUT, IWT, NRMTYP, OPT, XTRAP

      integer              I

      INTRINSIC            DABS, SQRT

      COMMON               /NSCOM1/ERRTOL, OPT, NRMTYP, XTRAP, ID, IWT,
     *                     IOUT

      if (NRMTYP.EQ.1) THEN
         DIFNRM = 0.0D+00
         DO I = 1, N
            DIFNRM = DMAX1(DIFNRM,dble(DABS(A(I)-B(I))))
         end do
      else
         DIFNRM = 0.0D+00
         DO I = 1, N
            DIFNRM = DIFNRM + dble(DABS(A(I)-B(I)))**2
         end do

         if (NRMTYP.EQ.2) DIFNRM = SQRT(DIFNRM)
         if (NRMTYP.EQ.3) DIFNRM = SQRT(DIFNRM/N)
      END IF

      return
      end
      subroutine efstat ( C, E, LGTOL, NTOL, NOK, ERFLG, TITLE,
     &  iout, w1, w2, w3, w4 )

c***********************************************************************
c
cc EFSTAT computes and printes normalized efficiency statistics.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
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
      implicit none

      double precision              C, E
      integer           IOUT, NOK, NTOL
      character*8       TITLE
      double precision              LGTOL(NTOL), W1(NTOL), W2(NTOL)
      integer           W3(NTOL), W4(NTOL)
      LOGICAL           ERFLG(NTOL)
      double precision              EQVTOL, S0, THETA, W1INT, W2INT, X
      integer           I, MSINT, NHI, NLO, SHI, SINT, SLO, W3INT, W4INT
      double precision              S(10)
      INTRINSIC          INT
      integer           FLOOR
c
c  STATEMENT FUNCTION
c  FLOOR FUNCTION VALID IF ARGUMENT X.GE.-100 WHICH IS OK HERE.
c
      FLOOR(X) = INT(X+100.0D+00) - 100

      if (NOK.LE.2) GO TO 200
c
c  TRANSFORM THE LOG10(TOL)'S TO NORMALIZED-EFFICIENCY VARIABLE:
c
      DO I = 1, NTOL
         S(I) = -(C+E*LGTOL(I))
      end do
c
c  FIND SET OF CONSECUTIVE TOL'S FOR WHICH INTEGRATION SUCCEEDED:
c
      DO NLO = 1, NTOL
         if ( .NOT. ERFLG(NLO)) GO TO 60
      end do
c
c  ELSE ALL INTEGRATIONS FOR THIS PROBLEM FAILED:
c
      GO TO 200
   60 CONTINUE
      NHI = NLO - 1
      DO I = NLO, NTOL
         if (ERFLG(I)) GO TO 100
         NHI = I
      end do

  100 CONTINUE

      if (NHI.LE.NLO) GO TO 200
      if (E.LE.0.0D+00) GO TO 220
c
c  FORM RANGE OF integer POWERS OF 10 FOR WHICH NORMALIZED STATISTICS
c  ARE TO BE PRINTED:
c
      SLO = -FLOOR(-S(NLO)+0.1D+00)
      SHI = FLOOR(S(NHI)+0.1D+00)
 
      if (SHI.LT.SLO) GO TO 240

      WRITE (IOUT,FMT=99999) TITLE
c
c  START OF LOOP TO PRINT A LINE OF STATISTICS FOR EACH POWER OF 10:
c
      I = NLO + 1
c
c  WHICH IS KNOWN TO BE .LE. NHI
c
      DO SINT = SLO, SHI

         S0 = dble(SINT)
c
c  MOVE INTERVAL S(I-1)..S(I) TO RIGHT WHILE S(I).LT.SINT:
c
  120    if ( S(I) .GE. S0 .OR. I .GE. NHI ) GO TO 140
         I = I + 1
         GO TO 120
  140    CONTINUE
c
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:
c
         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))

         MSINT = -SINT
         EQVTOL = -(C+S0)/E
         WRITE (IOUT,FMT=99998) MSINT, EQVTOL, W1INT, W2INT, W3INT,
     *     W4INT

      end do

  180 RETURN

  200 WRITE (IOUT,FMT=99997)
      return

  220 WRITE (IOUT,FMT=99996)
      return

  240 WRITE (IOUT,FMT=99995)
      return

99999 FORMAT (/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 (6X,'10**',I3,F8.2,F9.3,F7.3,1X,2I8)
99997 FORMAT (10X,'NOT ENOUGH SUCCESSFUL INTEGRATIONS TO FORM',1X,
     *       'NORMALIZED STATISTICS')
99996 FORMAT (10X,'DEPENDENCE OF ACCURACY ON TOLERANCE IS TOO',1X,
     *       'UNRELIABLE TO FORM NORMALIZED STATISTICS')
99995 FORMAT (10X,'NO POWERS OF TEN WITHIN RANGE OF TOLERANCES',1X,
     *       'USED: NO NORMALIZED STATISTICS')
      end
      subroutine evalu ( y, n, w, iwt, id )

c***********************************************************************
c
cc EVALU provides "true" values of the ODE solution at the endpoint.
c
c  Discussion:
c
c    This routine provides the "true" solution of the differential 
c    equation, evaluated at the endpoing of the integration interval.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
c    Output, double precision Y(N), the true solution evaluated at the endpoint.
c
c    Input, integer N, the problem dimension.
c
c    Input, double precision W(N), the optional weight vector.
c
c    Input, integer IWT, is positive if the weight vector is to be used.
c
c    Input, integer ID, the ID number of the problem being solved.
c
      implicit none

      integer n_max
      parameter ( n_max = 51 )

      integer i
      integer id
      integer iwt
      integer n
      double precision w(n_max)
      double precision y(n_max)
c
c  Problem A1
c
      if ( id .eq. 1 ) then

        y(1) = 2.061153353012535D-09
c
c  Problem A2
c
      else if ( id .eq. 2 ) then

        y(1) = 2.182178902359887D-01
c
c  Problem A3
c
      else if ( id .eq. 3 ) then

        y(1) = 2.491650271850414D+00
c
c  Problem A4
c
      else if ( id .eq. 4 ) then

        y(1) = 1.773016648131483D+01
c
c  Problem A5
c
      else if ( id .eq. 5 ) then

        y(1) = -7.887826688964196D-01
c
c  Problem B1
c
      else if ( id .eq. 11 ) then

        y(1) = 6.761876008576667D-01
        y(2) = 1.860816099640036D-01
c
c  Problem B2
c
      else if ( id .eq. 12 ) then

        y(1) = 1.000000001030576D+00
        y(2) = 1.000000000000000D+00
        y(3) = 9.999999989694235D-01
c
c  Problem B3
c
      else if ( id .eq. 13 ) then

        y(1) = 2.061153488557776D-09
        y(2) = 5.257228022048349D-02
        y(3) = 9.474277177183630D-01
c
c  Problem B4
c
      else if ( id .eq. 14 ) then

        y(1) = 9.826950928005993D-01
        y(2) = 2.198447081694832D+00
        y(3) = 9.129452507276399D-01
c
c  Problem B5
c
      else if ( id .eq. 15 ) then

        y(1) = -9.396570798729192D-01
        y(2) = -3.421177754000779D-01
        y(3) = 7.414126596199957D-01
c
c  Problem C1
c
      else if ( id .eq. 21 ) then

        y(1) = 2.061153622240064D-09
        y(2) = 4.122307244619555D-08
        y(3) = 4.122307244716968D-07
        y(4) = 2.748204829855288D-06
        y(5) = 1.374102414941961D-05
        y(6) = 5.496409659803266D-05
        y(7) = 1.832136553274552D-04
        y(8) = 5.234675866508716D-04
        y(9) = 1.308668966628220D-03
        y(10) = 9.979127409508656D-01
c
c  Problem C2
c
      else if ( id .eq. 22 ) then

        y(1) = 2.061153577984930D-09
        y(2) = 2.061153573736588D-09
        y(3) = 2.061153569488245D-09
        y(4) = 2.061153565239902D-09
        y(5) = 2.061153560991560D-09
        y(6) = 2.061153556743217D-09
        y(7) = 2.061153552494874D-09
        y(8) = 2.061153548246532D-09
        y(9) = 2.061153543998189D-09
        y(10) = 9.999999814496180D-01
c
c  Problem C3
c
      else if ( id .eq. 23 ) then

        y(1) = 2.948119211022058D-03
        y(2) = 5.635380154844266D-03
        y(3) = 7.829072515926013D-03
        y(4) = 9.348257908594937D-03
        y(5) = 1.007943610301970D-02
        y(6) = 9.982674171429909D-03
        y(7) = 9.088693332766085D-03
        y(8) = 7.489115195185912D-03
        y(9) = 5.322964130953349D-03
        y(10) = 2.762434379029886D-03
c
c  Problem C4
c
      else if ( id .eq. 24 ) then

        y(1) = 3.124111453721466D-03
        y(2) = 6.015416842150318D-03
        y(3) = 8.470021834842650D-03
        y(4) = 1.033682931733337D-02
        y(5) = 1.153249572873923D-02
        y(6) = 1.204549525737964D-02
        y(7) = 1.192957068015293D-02
        y(8) = 1.128883207111195D-02
        y(9) = 1.025804501391024D-02
        y(10) = 8.982017581934167D-03
        y(11) = 7.597500902492453D-03
        y(12) = 6.219920556824985D-03
        y(13) = 4.935916341009131D-03
        y(14) = 3.801432544256119D-03
        y(15) = 2.844213677587894D-03
        y(16) = 2.069123394222672D-03
        y(17) = 1.464687282843915D-03
        y(18) = 1.009545263941126D-03
        y(19) = 6.779354330227017D-04
        y(20) = 4.437815269118510D-04
        y(21) = 2.833264542938954D-04
        y(22) = 1.765005798796805D-04
        y(23) = 1.073342592697238D-04
        y(24) = 6.374497601777217D-05
        y(25) = 3.698645309704183D-05
        y(26) = 2.097466832643746D-05
        y(27) = 1.162956710412555D-05
        y(28) = 6.306710405783322D-06
        y(29) = 3.346286430868515D-06
        y(30) = 1.737760074184334D-06
        y(31) = 8.835366904275847D-07
        y(32) = 4.399520411127637D-07
        y(33) = 2.146181897152360D-07
        y(34) = 1.025981211654928D-07
        y(35) = 4.807864068784215D-08
        y(36) = 2.209175152474847D-08
        y(37) = 9.956251263138180D-09
        y(38) = 4.402193653748924D-09
        y(39) = 1.910149382204028D-09
        y(40) = 8.135892921473050D-10
        y(41) = 3.402477118549235D-10
        y(42) = 1.397485617545782D-10
        y(43) = 5.638575303049199D-11
        y(44) = 2.235459707956947D-11
        y(45) = 8.710498036398032D-12
        y(46) = 3.336554275346643D-12
        y(47) = 1.256679567784939D-12
        y(48) = 4.654359053128788D-13
        y(49) = 1.693559145599857D-13
        y(50) = 5.996593816663054D-14
        y(51) = 1.891330702629865D-14
c
c  Problem C5
c
      else if ( id .eq. 25 ) then

        y(1) = -4.792730224323733D+00
        y(2) = -2.420550725448973D+00
        y(3) = -9.212509306014886D-01
        y(4) = -4.217310404035213D+00
        y(5) = 7.356202947498970D+00
        y(6) = 3.223785985421212D+00
        y(7) = 4.035559443262270D+00
        y(8) = 1.719865528670555D+01
        y(9) = 7.478910794233703D+00
        y(10) = -2.998759326324844D+01
        y(11) = -4.107310937550929D+00
        y(12) = -9.277008321754407D-01
        y(13) = -2.442125302518482D+01
        y(14) = 2.381459045746554D+01
        y(15) = 1.492096306951359D+01
        y(16) = 3.499208963063806D-01
        y(17) = -5.748487687912825D-01
        y(18) = -2.551694020879149D-01
        y(19) = -5.237040978903326D-01
        y(20) = -2.493000463579661D-01
        y(21) = -8.045341642044464D-02
        y(22) = -3.875289237334110D-01
        y(23) = 5.648603288767891D-02
        y(24) = 3.023606472143342D-02
        y(25) = 4.133856546712445D-02
        y(26) = -2.862393029841379D-01
        y(27) = -1.183032405136207D-01
        y(28) = -1.511986457359206D-01
        y(29) = -2.460068894318766D-01
        y(30) = -3.189687411323877D-02
c
c  Problem D1
c
      else if ( id .eq. 31 ) then

        y(1) = 2.198835352008397D-01
        y(2) = 9.427076846341813D-01
        y(3) = -9.787659841058176D-01
        y(4) = 3.287977990962036D-01
c
c  Problem D2
c
      else if ( id .eq. 32 ) then

        y(1) = -1.777027357140412D-01
        y(2) = 9.467784719905892D-01
        y(3) = -1.030294163192969D+00
        y(4) = 1.211074890053952D-01
c
c  Problem D3
c
      else if ( id .eq. 33 ) then

        y(1) = -5.780432953035361D-01
        y(2) = 8.633840009194193D-01
        y(3) = -9.595083730380727D-01
        y(4) = -6.504915126712089D-02
c
c  Problem D4
c
      else if ( id .eq. 34 ) then

        y(1) = -9.538990293416394D-01
        y(2) = 6.907409024219432D-01
        y(3) = -8.212674270877433D-01
        y(4) = -1.539574259125825D-01
c
c  Problem D5
c
      else if ( id .eq. 35 ) then

        y(1) = -1.295266250987574D+00
        y(2) = 4.003938963792321D-01
        y(3) = -6.775390924707566D-01
        y(4) = -1.270838154278686D-01
c
c  Problem E1
c
      else if ( id .eq. 41 ) then

        y(1) = 1.456723600728308D-01
        y(2) = -9.883500195574063D-02
c
c  Problem E2
c
      else if ( id .eq. 42 ) then

        y(1) = 2.008149762174948D+00
        y(2) = -4.250887527320057D-02
c
c  Problem E3
c
      else if ( id .eq. 43 ) then

        y(1) = -1.004178858647128D-01
        y(2) = 2.411400132095954D-01
c
c  Problem E4
c
      else if ( id .eq. 44 ) then

        y(1) = 3.395091444646555D+01
        y(2) = 2.767822659672869D-01
c
c  Problem E5
c
      else if ( id .eq. 45 ) then

        y(1) = 1.411797390542629D+01
        y(2) = 2.400000000000002D+00
c
c  Problem F1
c
      else if ( id .eq. 51 ) then

        y(1) = -1.294460621213470D+01
        y(2) = -2.208575158908672D-15
c
c  Problem F2
c
      else if ( id .eq. 52 ) then

        y(1) = 70.03731057008607D+00
c
c  Problem F3
c
      else if ( id .eq. 53 ) then

        y(1) = -3.726957553088175D-01
        y(2) = -6.230137949234190D-01
c
c  Problem F4
c
      else if ( id .eq. 54 ) then

        y(1) = 9.815017249707434D-11
c
c  Problem F5
c
      else if ( id .eq. 55 ) then

        y(1) = 1.0D+00

      else

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'EVALU - Fatal error!'
        write ( *, '(a,i6)' ) '  Illegal problem ID = ', id
        stop

      end if

      if ( 0 .le. iwt ) then
        do i = 1, n
          y(i) = y(i) / w(i)
        end do
      end if

      return
      end
      subroutine fcn ( x, y, yp )

c***********************************************************************
c
cc FCN evaluates the derivative F(X,Y) for the ODE test problem.
c
c  Discussion:
c
c    This routine evaluates the derivative F(X,Y) corresponding to
c    the differential equation:
c      DY/DX = F(X,Y).
c
c    The routine stores the vector of derivatives in YP.
c
c    The equation being integrated is indicated by the value
c    of ID, which is passed through a common block.
c
c    The differential equation is scaled by the weight vector W
c    if this option has been selected by the flag IWT.
c
c  Modified:
c
c    17 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
c    Input, double precision X, the value of the independent variable.
c
c    Input, double precision Y(N), the value of the dependent variables.
c
c    Output, double precision YP(N), the value of the right hand side of
c    the differential equations associated with the dependent
c    variables.
c
      implicit none

      integer n_max
      parameter ( n_max = 51 )

      double precision c1
      double precision c2
      double precision d
      double precision ex
      integer i
      integer i3
      integer i3m2
      integer id
      integer itemp
      integer iwt
      integer j
      double precision k2
      integer l
      integer ll
      double precision m(5)
      integer mm
      double precision m0
      integer n
      integer nfcn
      double precision p
      double precision q(5,5)
      double precision r(5)
      double precision t
      double precision temp
      double precision w(n_max)
      double precision x
      double precision y(n_max)
      double precision yp(n_max)
      double precision ytemp(n_max)

      common /nscom5/ w, iwt, n, id
      common /nscom6/ nfcn
c
c  The following data is for problem C5, and defines the masses
c  of the five outer planets in solar units.
c  K2 is the gravitational constant.
c
      data m0 / 1.00000597682D+00 /
      data m / 0.954786104043D-03, 0.285583733151D-03, 
     &  0.437273164546D-04, 0.517759138449D-04, 0.277777777778D-05 /
      data k2 /2.95912208286D+00 /
c
c  This data is for problem F1.
c  C1 is pi**2 + 0.1**2.
c
      data c1 / 9.879604401089358D+00 /
c
c  This data iS for problem F5.
c  C2 iS sum i**(4/3) for i=1 to 19.
c
      data ex / 0.33333333333333333D+00 /
      data c2 / 438.4461015267790D+00 /

      nfcn = nfcn + 1

      iF ( 0 .le. iwt ) then
        do i = 1, n
          ytemp(i) = y(i)
          y(i) = y(i) * w(i)
        end do
      end if
c
c  Problem A1
c
      if ( id .eq. 1 ) then

        yp(1) = -y(1)
c
c  Problem A2
c
      else if ( id .eq. 2 ) then

        yp(1) = - 0.5D+00 * y(1) * y(1) * y(1)
c
c  Problem A3
c
      else if ( id .eq. 3 ) then

        yp(1) = y(1) * dcos ( x )
c
c  Problem A4
c
      else if ( id .eq. 4 ) then

        yp(1) = ( 1.0D+00 - y(1) / 20.0D+00 ) * y(1) / 4.0D+00
c
c  Problem A5
c
      else if ( id .eq. 5 ) then

        yp(1) = ( y(1) - x ) / ( y(1) + x )
c
c  Problem B1
c
      else if ( id .eq. 11 ) then

        d = y(1) - y(1) * y(2)
        yp(1) = d + d
        yp(2) = -( y(2) - y(1) * y(2) )
c
c  Problem B2
c
      else if ( id .eq. 12 ) then

        yp(1) = -y(1) + y(2)
        yp(3) = y(2) - y(3)
        yp(2) = -yp(1) - yp(3)
c
c  Problem B3
c
      else if ( id .eq. 13 ) then

        d = y(2) * y(2)
        yp(1) = -y(1)
        yp(2) = y(1) - d
        yp(3) = d
c
c  Problem B4
c
      else if ( id .eq. 14 ) then

        d = dsqrt ( y(1) * y(1) + y(2) * y(2) )
        yp(1) = -y(2) - y(1) * y(3) / d
        yp(2) = y(1) - y(2) * y(3) / d
        yp(3) = y(1) / d
c
c  Problem B5
c
      else if ( id .eq. 15 ) then

        yp(1) = y(2) * y(3)
        yp(2) = -y(1) * y(3)
        yp(3) = -0.51D+00 * y(1) * y(2)
c
c  Problem C1
c
      else if ( id .eq. 21 ) then

        yp(1) = -y(1)
        do i = 2, 9
          yp(i) = y(i-1) - y(i)
        end do
        yp(10) = y(9)
c
c  Problem C2
c
      else if ( id .eq. 22 ) then

        yp(1) = -y(1)
        do i = 2, 9
          yp(i) = dble ( i - 1 ) * y(i-1) - dble ( i ) * y(i)
        end do
        yp(10) = 9.0D+00 * y(9)
c
c  Problem C3
c
      else if ( id .eq. 23 ) then

        yp(1) = -2.0D+00 * y(1) + y(2)
        do i = 2, 9
          yp(i) = y(i-1) - 2.0D+00 * y(i) + y(i+1)
        end do
        yp(10) = y(9) - 2.0D+00 * y(10)
c
c  Problem C4
c
      else if ( id .eq. 24 ) then

        yp(1) = -2.0D+00 * y(1) + y(2)
        do i = 2, 50
          yp(i) = y(i-1) - 2.0D+00 * y(i) + y(i+1)
        end do
        yp(51) = y(50) - 2.0D+00 * y(51)
c
c  Problem C5
c
      else if ( id .eq. 25 ) then

        i = 0
        do l = 3, 15, 3
          i = i + 1
          p = y(l-2)**2 + y(l-1)**2 + y(l)**2
          r(i) = 1.0D+00 / ( p * dsqrt ( p ) )
          j = 0
          do ll = 3, 15, 3
            j = j + 1
            if ( ll .ne. l ) then
              p = ( y(l-2) - y(ll-2) )**2 + ( y(l-1) - y(ll-1) )**2 
     &          + ( y(l) - y(ll) )**2
              q(i,j) = 1.0D+00 / ( p * dsqrt ( p ) )
              q(j,i) = q(i,j)
            end if
          end do
        end do

        i3 = 0
        do i = 1, 5
          i3 = i3 + 3
          i3m2 = i3 - 2
          do ll = i3m2, i3
            mm = ll - i3
            yp(ll) = y(ll+15)
            p = 0.0D+00
            do j = 1, 5
              mm = mm + 3
              if ( j .ne. i ) then
                p = p + m(j) 
     &             * ( y(mm) * ( q(i,j) - r(j) ) - y(ll) * q(i,j) )
              end if
            end do
            yp(ll+15) = k2 * ( - ( m0 + m(i) ) * y(ll) * r(i) + p )
          end do
        end do
c
c  Problems D1, D2, D3, D4, D5
c
      else if ( id .eq. 31 .or.
     &  id .eq. 32 .or.
     &  id .eq. 33 .or.
     &  id .eq. 34 .or.
     &  id .eq. 35 ) then

        yp(1) = y(3)
        yp(2) = y(4)
        d = y(1) * y(1) + y(2) * y(2)
        d = dsqrt ( d * d * d )
        yp(3) = -y(1) / d
        yp(4) = -y(2) / d
c
c  Problem E1
c
      else if ( id .eq. 41 ) then

        yp(1) = y(2)
        yp(2) = - ( y(2) / ( x + 1.0D+00 ) 
     &    + ( 1.0D+00 - 0.25D+00 / ( x + 1.0D+00 )**2 ) * y(1) )
c
c  Problem E2
c
      else if ( id .eq. 42 ) then

        yp(1) = y(2)
        yp(2) = ( 1.0D+00 - y(1) * y(1) ) * y(2) - y(1)
c
c  Problem E3
c
      else if ( id .eq. 43 ) then

        yp(1) = y(2)
        yp(2) = y(1)**3 / 6.0D+00 
     &    - y(1) + 2.0D+00 * dsin ( 2.78535D+00 * x )
c
c  Problem E4
c
      else if ( id .eq. 44 ) then

        yp(1) = y(2)
        yp(2) = 0.032D+00 - 0.4D+00 * y(2) * y(2)
c
c  Problem E5
c
      else if ( id .eq. 45 ) then

        yp(1) = y(2)
        yp(2) = dsqrt ( 1.0D+00 + y(2) * y(2) ) / ( 25.0D+00 - x )
c
c  Problem F1
c
      else if ( id .eq. 51 ) then

        yp(1) = y(2)
        yp(2) = 0.2D+00 * y(2) - c1 * y(1)
        itemp = idint ( x )
        if ( mod ( itemp, 2 ) .eq. 1 ) then
          yp(2) = yp(2) - 1.0D+00
        else
          yp(2) = yp(2) + 1.0D+00
        end if
c
c  Problem F2
c
      else if ( id .eq. 52 ) then

        itemp = idint ( x ) 
        if ( mod ( itemp, 2 ) .eq. 1 ) then
          yp(1) = 55.0D+00 - 0.5D+00 * y(1)
        else
          yp(1) = 55.0D+00 - 1.5D+00 * y(1)
        end if
c
c  Problem F3
c
      else if ( id .eq. 53 ) then

        yp(1) = y(2)
        yp(2) = 0.01D+00 * y(2) * ( 1.0D+00 - y(1)**2 ) - y(1) -
     *    dabs ( dsin ( 3.1415926535897932D+00 * x ) )
c
c  Problem F4
c
      else if ( id .eq. 54 ) then

        if ( x .le. 10.0D+00 ) then
          temp = x - 5.0D+00
          yp(1) = -2.0D+00 / 21.0D+00 
     &      - 120.0D+00 * temp / ( 1.0D+00 + 4.0D+00 * temp**2 )**16
        else
          yp(1) = -2.0D+00 * y(1)
        end if
c
c  Problem F5
c
      else if ( id .eq. 55 ) then

        temp = 0.0D+00

        do i = 1, 19
          t = x - dble ( i )
          temp = temp + dsign ( dabs ( t )**ex, t )
        end do

        yp(1) = y(1) * ( 4.0D+00 / ( 3.0D+00 * c2 ) ) * temp

      else

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'FCN - Fatal error!'
        write ( *, '(a,i6)' ) '  illegal problem iD = ', id
        stop

      end if

      if ( 0 .le. iwt ) then
        do i = 1, n
          yp(i) = yp(i) / w(i)
          y(i) = ytemp(i)
        end do
      end if

      return
      end
      subroutine fcn2 ( n, x, y, yp )

c***********************************************************************
c
cc FCN2 transfers function values from FCN to TRUE.
c
c  Discussion:
c
c    The METHOD package integrates an ODE whose right hand side is
c    evaluated by a routine whose name is an arbitrary external,
c    but which must have the form
c
c      subroutine FCN2 ( N, X, Y, YP )
c
c    while the ODE test package has a generic interface for the
c    function whose name and form are
c
c      subroutine FCN ( X, Y, YP )
c
c    This routine serves as an interface.  METHOD is told to call
c    FCN2 (the name is specified as an external function argument to 
c    METHOD) and FCN2, in turn, calls FCN.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
c    Input, integer N, the number of dependent variables.
c
c    Input, double precision X, the value of the independent variable.
c
c    Input, double precision Y(N), the value of the dependent variables.
c
c    Output, double precision YP(N), the value of the right hand side of
c    the differential equations associated with the dependent
c    variables.
c
      implicit none

      integer n

      double precision x
      double precision y(n)
      double precision yp(n)

      call fcn ( x, y, yp )

      return
      end
      subroutine ivalu ( n, xstart, xend, hbegin, hmax, y, fcntim, 
     &  w, iwt, id )

c***********************************************************************
c
cc IVALU provides initial values to the ODE test problem.
c
c  Discussion:
c
c    This routine provides the initial values that specify the mathematical
c    problem, as well as various algorithmic parameters.
c
c  Modified:
c
c    17 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
c    Output, integer N, the dimension of the problem.
c
c    Output, double precision XSTART, the initial value of the
c    independent variable.
c
c    Output, double precision XEND, the final value of the
c    independent variable.
c
c    Output, double precision HBEGIN, a suggested initial stepsize.
c
c    Output, double precision HMAX, the maximum stepsize allowed.
c
c    Output, double precision Y(N), the initial condition for the
c    dependent variables.
c
c    Output, double precision FCNTIM, the average computer time
c    required for a derivative evaluation.  NOT USED, SET TO 0.
c
c    Output, double precision WT(N), weights used to scale the problem,
c    if the scaling option is selected.
c
c    Input, integer IWT, is nonnegative if weighting is to be used.
c
c    Input, integer ID, the identifier for the problem.
c
      implicit none

      double precision e
      double precision fcntim
      double precision hb
      double precision hbegin
      double precision hm
      double precision hmax
      integer i
      integer id
      integer iwt
      integer n
      double precision w(51)
      double precision xe
      double precision xend
      double precision xs
      double precision xstart
      double precision y(51)

      data hb / 1.0D+00 /
      data hm / 20.0D+00 /
      data xe / 20.0D+00 /
      data xs / 0.0D+00 /

      hmax = hm
      hbegin = hb
      xstart = xs
      xend = xe
c
c  Problem A1
c
      if ( id .eq. 1 ) then

        fcntim = 0.0D+00
        n = 1
        w(1) = 0.100D+01
        y(1) = 1.0D+00
c
c  Problem A2
c
      else if ( id .eq. 2 ) then

        fcntim = 0.0D+00
        n = 1
        w(1) = 0.100D+01
        y(1) = 1.0D+00
c
c  Problem A3
c
      else if ( id .eq. 3 ) then

        fcntim = 0.0D+00
        n = 1
        w(1) = 0.271D+01
        y(1) = 1.0D+00
c
c  Problem A4
c
      else if ( id .eq. 4 ) then

        fcntim = 0.0D+00
        n = 1
        w(1) = 0.177D+02
        y(1) = 1.0D+00
c
c  Problem A5
c
      else if ( id .eq. 5 ) then

        fcntim = 0.0D+00
        n = 1
        w(1) = 0.620D+01
        y(1) = 4.0D+00
c
c  Problem B1
c
      else if ( id .eq. 11 ) then

        fcntim = 0.0D+00
        n = 2
        w(1) = 0.425D+01
        w(2) = 0.300D+01
        y(1) = 1.0D+00
        y(2) = 3.0D+00
c
c  Problem B2
c
      else if ( id .eq. 12 ) then

        fcntim = 0.0D+00
        n = 3
        w(1) = 0.200D+01
        w(2) = 0.100D+01
        w(3) = 0.100D+01
        y(1) = 2.0D+00
        y(2) = 0.0D+00
        y(3) = 1.0D+00
c
c  Problem B3
c
      else if ( id .eq. 13 ) then

        fcntim = 0.0D+00
        n = 3
        w(1) = 0.100D+01
        w(2) = 0.519D+00
        w(3) = 0.947D+00
        y(1) = 1.0D+00
        y(2) = 0.0D+00
        y(3) = 0.0D+00
c
c  Problem B4
c
      else if ( id .eq. 14 ) then

        fcntim = 0.0D+00
        n = 3
        w(1) = 0.300D+01
        w(2) = 0.220D+01
        w(3) = 0.100D+01
        y(1) = 3.0D+00
        y(2) = 0.0D+00
        y(3) = 0.0D+00
c
c  Problem B5
c
      else if ( id .eq. 15 ) then

        fcntim = 0.0D+00
        n = 3
        w(1) = 0.100D+01
        w(2) = 0.100D+01
        w(3) = 0.100D+01
        y(1) = 0.0D+00
        y(2) = 1.0D+00
        y(3) = 1.0D+00
c
c  Problem C1
c
      else if ( id .eq. 21 ) then

        fcntim = 0.0D+00
        n = 10
        w(1) = 0.100D+01
        w(2) = 0.368D+00
        w(3) = 0.271D+00
        w(4) = 0.224D+00
        w(5) = 0.195D+00
        w(6) = 0.175D+00
        w(7) = 0.161D+00
        w(8) = 0.149D+00
        w(9) = 0.139D+00
        w(10) = 0.998D+00
        y(1) = 1.0D+00
        do i = 2, n
          y(i) = 0.0D+00
        end do
c
c  Problem C2
c
      else if ( id .eq. 22 ) then

        fcntim = 0.0D+00
        n = 10
        w(1) = 0.100D+01
        w(2) = 0.250D+00
        w(3) = 0.148D+00
        w(4) = 0.105D+00
        w(5) = 0.818D-01
        w(6) = 0.669D-01
        w(7) = 0.566D-01
        w(8) = 0.491D-01
        w(9) = 0.433D-01
        w(10) = 0.100D+01
        y(1) = 1.0D+00
        do i = 2, n
          y(i) = 0.0D+00
        end do
c
c  Problem C3
c
      else if ( id .eq. 23 ) then

        fcntim = 0.0D+00
        n = 10
        w(1) = 0.100D+01
        w(2) = 0.204D+00
        w(3) = 0.955D-01
        w(4) = 0.553D-01
        w(5) = 0.359D-01
        w(6) = 0.252D-01
        w(7) = 0.184D-01
        w(8) = 0.133D-01
        w(9) = 0.874D-02
        w(10) = 0.435D-02
        y(1) = 1.0D+00
        do i = 2, n
          y(i) = 0.0D+00
        end do
c
c  Problem C4
c
      else if ( id .eq. 24 ) then

        fcntim = 0.0D+00
        n = 51
        w(1) = 0.100D+01
        w(2) = 0.204D+00
        w(3) = 0.955D-01
        w(4) = 0.553D-01
        w(5) = 0.359D-01
        w(6) = 0.252D-01
        w(7) = 0.186D-01
        w(8) = 0.143D-01
        w(9) = 0.113D-01
        w(10) = 0.918D-02
        w(11) = 0.760D-02
        w(12) = 0.622D-02
        w(13) = 0.494D-02
        w(14) = 0.380D-02
        w(15) = 0.284D-02
        w(16) = 0.207D-02
        w(17) = 0.146D-02
        w(18) = 0.101D-02
        w(19) = 0.678D-03
        w(20) = 0.444D-03
        w(21) = 0.283D-03
        w(22) = 0.177D-03
        w(23) = 0.107D-03
        w(24) = 0.637D-04
        w(25) = 0.370D-04
        w(26) = 0.210D-04
        w(27) = 0.116D-04
        w(28) = 0.631D-05
        w(29) = 0.335D-05
        w(30) = 0.174D-05
        w(31) = 0.884D-06
        w(32) = 0.440D-06
        w(33) = 0.215D-06
        w(34) = 0.103D-06
        w(35) = 0.481D-07
        w(36) = 0.221D-07
        w(37) = 0.996D-08
        w(38) = 0.440D-08
        w(39) = 0.191D-08
        w(40) = 0.814D-09
        w(41) = 0.340D-09
        w(42) = 0.140D-09
        w(43) = 0.564D-10
        w(44) = 0.224D-10
        w(45) = 0.871D-11
        w(46) = 0.334D-11
        w(47) = 0.126D-11
        w(48) = 0.465D-12
        w(49) = 0.169D-12
        w(50) = 0.600D-13
        w(51) = 0.189D-13
        y(1) = 1.0D+00
        do i = 2, n
          y(i) = 0.0D+00
        end do
c
c  Problem C5
c
      else if ( id .eq. 25 ) then

        fcntim = 0.0D+00
        n = 30
        w(1) = 0.545D+01
        w(2) = 0.471D+01
        w(3) = 0.203D+01
        w(4) = 0.664D+01
        w(5) = 0.834D+01
        w(6) = 0.346D+01
        w(7) = 0.113D+02
        w(8) = 0.172D+02
        w(9) = 0.748D+01
        w(10) = 0.302D+02
        w(11) = 0.411D+01
        w(12) = 0.144D+01
        w(13) = 0.244D+02
        w(14) = 0.284D+02
        w(15) = 0.154D+02
        w(16) = 0.764D+00
        w(17) = 0.661D+00
        w(18) = 0.284D+00
        w(19) = 0.588D+00
        w(20) = 0.366D+00
        w(21) = 0.169D+00
        w(22) = 0.388D+00
        w(23) = 0.190D+00
        w(24) = 0.877D-01
        w(25) = 0.413D-01
        w(26) = 0.289D+00
        w(27) = 0.119D+00
        w(28) = 0.177D+00
        w(29) = 0.246D+00
        w(30) = 0.319D-01
        y(1) = 3.42947415189D+00
        y(2) = 3.35386959711D+00
        y(3) = 1.35494901715D+00
        y(4) = 6.64145542550D+00
        y(5) = 5.97156957878D+00
        y(6) = 2.18231499728D+00
        y(7) = 11.2630437207D+00
        y(8) = 14.6952576794D+00
        y(9) = 6.27960525067D+00
        y(10) = -30.1552268759D+00
        y(11) = 1.65699966404D+00
        y(12) = 1.43785752721D+00
        y(13) = -21.1238353380D+00
        y(14) = 28.4465098142D+00
        y(15) = 15.3882659679D+00
        y(16) = -0.557160570446D+00
        y(17) = 0.505696783289D+00
        y(18) = 0.230578543901D+00
        y(19) = -0.415570776342D+00
        y(20) = 0.365682722812D+00
        y(21) = 0.169143213293D+00
        y(22) = -0.325325669158D+00
        y(23) = 0.189706021964D+00
        y(24) = 0.0877265322780D+00
        y(25) = -0.0240476254170D+00
        y(26) = -0.287659532608D+00
        y(27) = -0.117219543175D+00
        y(28) = -0.176860753121D+00
        y(29) = -0.216393453025D+00
        y(30) = -0.0148647893090D+00
c
c  Problem D31.
c
      else if ( id .eq. 31 ) then

        fcntim = 0.0D+00
        e = 0.1D+00
        n = 4
        w(1) = 0.110D+01
        w(2) = 0.995D+00
        w(3) = 0.101D+01
        w(4) = 0.111D+01
        y(1) = 1.0D+00 - e
        y(2) = 0.0D+00
        y(3) = 0.0D+00
        y(4) = dsqrt ( ( 1.0D+00 + e ) / ( 1.0D+00 - e ) )
c
c  Problem D32.
c
      else if ( id .eq. 32 ) then

        fcntim = 0.0D+00
        e = 0.3D+00
        n = 4
        w(1) = 0.130D+01
        w(2) = 0.954D+00
        w(3) = 0.105D+01
        w(4) = 0.136D+01
        y(1) = 1.0D+00 - e
        y(2) = 0.0D+00
        y(3) = 0.0D+00
        y(4) = dsqrt ( ( 1.0D+00 + e ) / ( 1.0D+00 - e ) )
c
c  Problem D33.
c
      else if ( id .eq. 33 ) then

        fcntim = 0.0D+00
        e = 0.5D+00
        n = 4
        w(1) = 0.150D+01
        w(2) = 0.866D+00
        w(3) = 0.115D+01
        w(4) = 0.173D+01
        y(1) = 1.0D+00 - e
        y(2) = 0.0D+00
        y(3) = 0.0D+00
        y(4) = dsqrt ( ( 1.0D+00 + e ) / ( 1.0D+00 - e ) )
c
c  Problem D34.
c
      else if ( id .eq. 34 ) then

        fcntim = 0.0D+00
        e = 0.7D+00
        n = 4
        w(1) = 0.170D+01
        w(2) = 0.714D+00
        w(3) = 0.140D+01
        w(4) = 0.238D+01
        y(1) = 1.0D+00 - e
        y(2) = 0.0D+00
        y(3) = 0.0D+00
        y(4) = dsqrt ( ( 1.0D+00 + e ) / ( 1.0D+00 - e ) )
c
c  Problem D35.
c
      else if ( id .eq. 35 ) then

        fcntim = 0.0D+00
        e = 0.9D+00
        n = 4
        w(1) = 0.190D+01
        w(2) = 0.436D+00
        w(3) = 0.229D+01
        w(4) = 0.436D+01
        y(1) = 1.0D+00 - e
        y(2) = 0.0D+00
        y(3) = 0.0D+00
        y(4) = dsqrt ( ( 1.0D+00 + e ) / ( 1.0D+00 - e ) )
c
c  Problem E1
c
      else if ( id .eq. 41 ) then

        fcntim = 0.0D+00
        n = 2
        w(1) = 0.679D+00
        w(2) = 0.478D+00
        e = 0.79788456080286536D+00
        y(1) = e * 0.84147098480789651D+00
        y(2) = e * 0.11956681346419146D+00
c
c  Problem E2
c
      else if ( id .eq. 42 ) then

        fcntim = 0.0D+00
        n = 2
        w(1) = 0.201D+01
        w(2) = 0.268D+01
        y(1) = 2.0D+00
        y(2) = 0.0D+00
c
c  Problem E3
c
      else if ( id .eq. 43 ) then

        fcntim = 0.0D+00
        n = 2
        w(1) = 0.116D+01
        w(2) = 0.128D+01
        y(1) = 0.0D+00
        y(2) = 0.0D+00
c
c  Problem E4
c
      else if ( id .eq. 44 ) then

        fcntim = 0.0D+00
        n = 2
        w(1) = 0.340D+02
        w(2) = 0.277D+00
        y(1) = 3.0D+01
        y(2) = 0.0D+00
c
c  Problem E5
c
      else if ( id .eq. 45 ) then

        fcntim = 0.0D+00
        n = 2
        w(1) = 0.141D+02
        w(2) = 0.240D+01
        y(1) = 0.0D+00
        y(2) = 0.0D+00
c
c  Problem F1
c
      else if ( id .eq. 51 ) then

        fcntim = 0.0D+00
        n = 2
        w(1) = 0.129D+02
        w(2) = 0.384D+02
        y(1) = 0.0D+00
        y(2) = 0.0D+00
        HMAX = 1.0D+00
c
c  Problem F2
c
      else if ( id .eq. 52 ) then

        fcntim = 0.0D+00
        n = 1
        w(1) = 0.110D+03
        y(1) = 110.D0
        hmax = 1.0D+00
c
c  Problem F3
c
      else if ( id .eq. 53 ) then

        fcntim = 0.0D+00
        n = 2
        w(1) = 1.31D+00
        w(2) = 0.737D+00
        y(1) = 0.0D+00
        y(2) = 0.0D+00
        hmax = 1.0D+00
        hbegin = 0.9D+00
c
c  Problem F4
c
      else if ( id .eq. 54 ) then

        fcntim = 0.0D+00
        n = 1
        w(1) = 1.52D+00
        y(1) = 1.0D+00
        hmax = 10.0D+00
c
c  Problem F5
c
      else if ( id .eq. 55 ) then

        fcntim = 0.0D+00
        n = 1
        w(1) = 1.0D+00
        y(1) = 1.0D+00
        hmax = 20.0D+00

      else

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'IVALU - Fatal error!'
        write ( *, '(a,i6)' ) '  Illegal problem ID = ', id
        stop

      end if

      if ( 0 .le. iwt ) then
        do i = 1, n
          y(i) = y(i) / w(i)
        end do
      end if

      return
      end
      subroutine lsqfit ( x, y, miss, n, nn, c0, c1, res )

c***********************************************************************
c
cc LSQFIT fits a linear model to data.
c
c  Discussion:
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  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
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
      implicit none

      integer n

      double precision C0
      double precision C1
      double precision RES
      integer nn

      double precision X(N), Y(N)
      LOGICAL           MISS(N)

      double precision SX, SXX, SXY, SY, XNN
      integer           I

      NN = 0
      SX = 0.0D+00
      SY = 0.0D+00

      DO I = 1, N
        if ( .not. MISS(I) ) then
          NN = NN + 1
          SX = SX + X(I)
          SY = SY + Y(I)
        end if
      end do

      if ( NN .LE. 1 ) then
        return
      end if

      XNN = NN
      SX = SX / XNN
      SY = SY / XNN
      SXX = 0.0D+00
      SXY = 0.0D+00
      DO I = 1, N
         if ( .not. MISS(I)) then
           SXX = SXX + (X(I)-SX)**2
           SXY = SXY + (X(I)-SX)*(Y(I)-SY)
         end if
      end do

      C1 = SXY/SXX
      C0 = SY - C1*SX
      RES = 0.0D+00
      DO I = 1, N
        if ( .NOT. MISS(I)) then
          RES = RES + (Y(I)-SY-C1*(X(I)-SX))**2
        end if
      end do

      RES = SQRT ( RES / XNN )

      return
      end
      subroutine parchk ( OPT, NORMEF, NRMTYP, TOL, IDLIST, ntol, 
     &  ngrp, grplst, lgtol, flag )

c***********************************************************************
c
cc PARCHK checks the parameters.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
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
c          ...  (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  
      implicit none

      double precision FLAG
      integer           NGRP, NORMEF, NRMTYP, NTOL, OPT
      double precision LGTOL(10), TOL(11)
      integer           GRPLST(2,6), IDLIST(60)
      double precision BIG, TOLPRV
      integer           ENDLST, I, ID, IID, ISAV, KCLASS, LENIDS,
     *                  LENTOL, MAXGRP, NCLASS
      integer           NSYSTM(6)

      INTRINSIC         IABS

      DATA              ENDLST/-1/, BIG/1E20/
      DATA              NCLASS/6/, NSYSTM/5, 5, 5, 5, 5, 5/, MAXGRP/6/,
     *                  LENTOL/11/, LENIDS/60/

      FLAG = 0.0D+00
      if (OPT.LT.1 .OR. OPT.GT.4) FLAG = 1.0D+00
      if (NORMEF.LT.0 .OR. NORMEF.GT.2) FLAG = 10.*FLAG + 2.0D+00
      if (OPT.EQ.1 .AND. NORMEF.EQ.2) FLAG = 10.*FLAG + 3.0D+00
      if (NRMTYP.LT.1 .OR. NRMTYP.GT.3) FLAG = 10.*FLAG + 8.0D+00
c
c  TOLERANCES:
c
      NTOL = 0
      TOLPRV = BIG

      DO I = 1, LENTOL
         if (TOL(I).LT.0.0D+00 .OR. TOL(I).GE.TOLPRV) then
           FLAG = 10.*FLAG + 4.0D+00
         end if
         if (TOL(I).EQ.0.0D+00) GO TO 40
         NTOL = NTOL + 1
         TOLPRV = TOL(I)
      end do
c
c  NO TERMINATING 0 IN TOLERANCE LIST:
c
      FLAG = 10.0D+00 * FLAG + 5.0D+00
c
c  CHECK FOR EMPTY TOLERANCE LIST:
c
   40 if (NTOL.EQ.0) FLAG = 10.*FLAG + 5.0D+00
c
c  LIST OF GROUPS OF PROBLEMS:
c
      NGRP = 0
      I = 0
c
c  WHILE NEXT ID IN LIST ISNT 0 OR END OF LIST:
c
   60 I = I + 1
      ID = ENDLST
      if (I.LE.LENIDS) ID = IDLIST(I)

      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:
c
   80 if (ID.EQ.0) GO TO 140
      if (ID.EQ.ENDLST) GO TO 180
c
c  TRANSLATE ID INTO CLASS and NUMBER WITHIN CLASS,
c  IGNORING SIGN (WHICH SELECTS SCALED/UNSCALED OPTION):
c
      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.0D+00

  120 CONTINUE
c
c  GET NEXT ID AS ABOVE:
c
      I = I + 1
      ID = ENDLST
      if (I.LE.LENIDS) ID = IDLIST(I)
      GO TO 80
c
c  NEW GROUP FORMED:
c
  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:
c
  160 if (NGRP.LE.0) GO TO 180
      GO TO 200

  180 FLAG = 10.0D+00 * FLAG + 7.0D+00
c
c  IF ALL OK, COMPUTE LOGS OF TOLERANCES:
c
  200 continue

      if ( flag .eq. 0.0D+00 ) then
        do i = 1, ntol
          lgtol(i) = dlog10 ( tol(i) )
        end do
      end if

      return
      end
      subroutine plot ( x, y, iflag )

c***********************************************************************
c
cc PLOT plots local error information.
c
c  Discussion:
c
c    ROUTINE TO FORM PLOTS OF LOCAL ERROR INFORMATION FOR DETEST, USING
c    AN ARRAY K WHICH IS IN 'SAVE' STORAGE.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
c  IF IFLAG<=0, IT RESETS ARRAY K TO ZERO.
c
c  IF IFLAG=1, THE ROUTINE ENTERS (X,Y) ON THE SCATTER-DIAGRAM
c  REPRESENTED BY K.  HERE X,Y ARE >= 0, AND THE RANGE 0 TO INFINITY IS
c  SPLIT INTO CLASS-INTERVALS NUMBERED I = NLO .. NHI, THE I-TH INTERVAL
c  BEING 2**(I-1) <= X < 2**I EXCEPT THAT THE NLO-TH ONE INCLUDES ALL
c  X BELOW 2**NLO AND THE NHI-TH INCLUDES ALL X >=2**(NHI-1).
c
c  IF IFLAG=2, THE SCATTER DIAGRAM IS PRINTED OUT.
c
      implicit none

      integer         NLO, NHI
      double precision  ALOG2
      integer         NMIN, LINLEN
      double precision            XYMIN
      PARAMETER       (NLO=-7,NHI=4,ALOG2=0.69314718D+00,NMIN=NLO-1,
     *                LINLEN=3*(NHI-NLO+1)+1,XYMIN=2.0D+00**NMIN)

      double precision            X, Y
      integer         IFLAG

      double precision            C, P, T
      integer         I, IOUT, J, JL, KMAX, KTOT
      character*(LINLEN) LINE
      character*(LINLEN) LINE1
      character*(LINLEN) LINE2
      character*(LINLEN) LINE3
      character*(LINLEN) LINE4

      integer         K(NLO:NHI,NLO:NHI)

      double precision            CONST
      character*3     STR3
      EXTERNAL        CONST, STR3

      INTRINSIC       MAX, MIN, NINT

      integer         ICLAS, ICLAS0

      SAVE            K, KTOT, KMAX, IOUT

      DATA            LINE1/'+--+--+--+--+--+--+--+--+--+--+--+--+'/,
     *                LINE2/'+                                   +'/,
     *                LINE3/'|                                   |'/,
     *                LINE4/'  2  2  2  2  2  2  2  2  2  2  2    '/
c
c     .. Statement Function definitions ..
      ICLAS0(T) = NMIN + NINT(DLOG(MAX(1.0D+00,T/XYMIN))/ALOG2)
      ICLAS(T) = MIN(MAX(ICLAS0(T),NLO),NHI)

      if (IFLAG.LE.0) THEN

         IOUT = CONST(3)
         KTOT = 0
         KMAX = 0

         DO I = NLO, NHI
            DO J = NLO, NHI
               K(I,J) = 0
            end do
         end do

      ELSE if (IFLAG.EQ.1) THEN
         if (X.LT.0.0D+00 .OR. Y.LT.0.0D+00) THEN
            WRITE (IOUT,FMT=*)
     *        ' ERROR IN ARGUMENTS TO DETEST PLOT ROUTINE', X, Y
            STOP
         END IF
         I = ICLAS(X)
         J = ICLAS(Y)
         K(I,J) = K(I,J) + 1
         KTOT = KTOT + 1
         KMAX = MAX(KMAX,K(I,J))
      ELSE

         C = KTOT
         DO I = NHI, NLO, -1
            LINE = LINE3
            DO J = NLO, NHI
               JL = J - NLO
               P = K(J,I)/C
               LINE(3*JL+1:3*JL+3) = STR3(P)
            end do
            if (LINE(1:1).EQ.' ') LINE(1:1) = '|'
            if (I.EQ.NHI) THEN
               WRITE (IOUT,FMT='(1X,15X,''INFINITY '',A)') LINE1
               WRITE (IOUT,FMT='(1X,20X,''    '',A)') LINE
            ELSE
               WRITE (IOUT,FMT='(1X,15X,I8,1X,A)') I, LINE2
               WRITE (IOUT,FMT='(1X,20X,''2   '',A)') LINE
            END IF
         end do

         WRITE (IOUT,FMT='(1X,24X,A)') LINE1
         WRITE (IOUT,FMT='(/1X,25X,30I3)') (J,J=NLO,NHI-1)
         WRITE (IOUT,FMT='(1X,24X,A)') LINE4
      END IF

      return
      end
      function ratio ( m, n )

c***********************************************************************
c
cc RATIO computes the ratio of two integers.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
      implicit none

      integer m
      integer n
      double precision ratio

      if (N.NE.0) then
        RATIO = dble(M)/dble(N)
      else
        RATIO = 1.0D+20
      end if

      return
      end
      subroutine stats ( x, y, errbnd, errest )

c***********************************************************************
c
cc STATS instruments the ODE-solver being tested.
c
c  Discussion:
c
c     STATS 'INSTRUMENTS' THE ODE-SOLVER BEING TESTED, BY COMPUTING
c     THE DEVIATION OF THE SOLUTION COMPUTED IN ROUTINE METHOD FROM
c     THE 'TRUE' GLOBAL AND LOCAL SOLUTIONS IF REQUESTED, AND BY
c     ACCUMULATING VARIOUS ASSOCIATED STATISTICS. IT ALSO PERFORMS
c     VARIOUS INITIALIZATION DUTIES, DEPENDING ON THE VALUE OF IFLAG
c     ON ENTRY.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c     ON ENTRY
c     X,Y   - MUST HOLD 'SOLVER' SOLUTION AT CURRENT STEP
c     ERREST- MUST HOLD ESTIMATED LOCAL ERROR VECTOR AT THIS STEP
c             DEFINED AS (COMPUTED Y) - (TRUE LOCAL SOLUTION AT NEW X).
c             SINCE ABSOLUTE ERROR-CONTROL IS SPECIFIED, THIS IS THE
c             VECTOR WHOSE NORM IS MAINTAINED BELOW ERRBND BY 'METHOD'.
c             IT IS ASSUMED THAT 'METHOD' USES ONE OF THE 3 NORMS
c             OFFERED BY THE PACKAGE, AND NRMTYP MUST BE SET SUITABLY.
c     ERRBND- MUST HOLD TOLERANCE BELOW WHICH THE NORM OF ERREST IS
c             BEING HELD AT THIS STEP. USUALLY SAME AS ERRTOL BUT WILL
c             BE DIFFERENT AND VARY WITH STEPSIZE if (EG) A PER-UNIT-
c             STEP ERROR CRITERION IS USED.
c
c     STORAGE FOR VARIOUS SOLUTIONS:
c     X,Y      - CURRENT SOLUTION COMPUTED BY METHOD, PASSED IN
c                VIA ARGUMENT LIST.
c     XOLD,YOLD- VALUES OF X,Y AT AN OLD MESHPOINT OF METHOD,
c                USUALLY THE LAST ONE BUT OLDER IF A LUMPED
c                STEP IS BEING FORMED (SEE BELOW).
c                IF IFLAG = 0, NEITHER XOLD NOR YOLD IS USED.
c                YOLD IS NOT USED UNLESS STATISTICS ON LOCAL ERROR
c                ARE BEING COMPILED (IFLAG=3 AND OPT=3).
c                THE 'TRUE' LOCAL SOLUTION IS OBTAINED BY INTEG-
c                RATING FROM XOLD,YOLD TO THE CURRENT X.
c                XOLD,YOLD ARE USED AS THE ACTUAL ARGUMENTS IN THIS
c                INTEGRATION, AND ARE THEN UPDATED TO HOLD X,Y IN
c                PREPARATION FOR NEXT CALL TO STATS.
c     XT       - LAST MESHPOINT OF METHOD.
c     XOLDG    - INDEP VAR FOR 'TRUE' GLOBAL SOLUTION, IN COMMON.
c     YOLDG    - 'TRUE' GLOBAL SOLUTION AT XOLDG, HELD IN COMMON.
c                UPDATED BY CALLING TRUE AT EACH CALL TO STATS IF
c                DETAILED STATISTICS ARE BEING COMPILED (IFLAG = 3)
c     YSTAR    - ONLY USED IF OPT.EQ.4.  IF SOLVER DOES NOT DO LOCAL
c                EXTRAPOLATION, WE FORM THE LOCALLY EXTRAPOLATED
c                SOLUTION IN YSTAR.
c
      implicit none
  
      double precision CG(24)
      double precision ERRBND
      double precision ERREST(51)
      double precision WG(51,9)
      double precision X
      double precision Y(51)
      double precision YOLD(51)
      double precision YOLDG(51)

      double precision ERLUMP, ERRTOL, HSTART, PRECIS, XEND, XFIN, XOLD,
     *                 XOLD1, XOLDG, XT, XTRUE
      double precision   GEMX, GEND, LEMXSC, OVHD, TIME, TRUTIM
      integer          ID, IFLAG, INDG, INDL, IOUT, IWT, N, NBAD, NDCV,
     *                 NFCN, NFCN1, NRMTYP, NSTART, NSTL, NSTP, NTRU,
     *                 OPT, XTRAP

      double precision HLUMP, HMIN, YNORM
      double precision ESTSC, LEERSC, LESC, TRUT0
      integer          I, NDIM, NNFCN

      double precision CL(24), WL(51,9), YSTAR(51), ZERO(51)

      double precision             CLOCK, CONST, DIFNRM
      EXTERNAL         CLOCK, CONST, DIFNRM

      EXTERNAL         FCN2, PLOT, TRUE

      INTRINSIC        DABS, DMAX1

      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           /NSCOM4/XOLD1, XOLD, YOLD, XOLDG, YOLDG, CG, WG,
     *                 XT, PRECIS, ERLUMP
      COMMON           /NSCOM6/NFCN1

      DATA             NDIM/51/, ZERO/51*0.0D+00/
c
c  IF IFLAG = 0 METHOD IS BEING TIMED.
c  
      if ( IFLAG .EQ. 0 ) then
        return
      end if
c
c  IF IFLAG = 1 INITIALIZE VARIABLES TO DO WITH FINDING FIRST STEP-
c  SIZE, ASSESSING LUMPED STEPS AND COMPUTING TRUE GLOBAL SOLUTION.
c  RESET INDL, OTHERWISE A LOCAL FAILURE (INDL<0) ON A PREVIOUS
c  INTEGRATION WILL BE DEEMED A FAILURE ON THIS ONE.
c  1ST 9 ELEMENTS OF CG MUST BE INITIALIZED; WE INITIALIZE
c  MORE TO AID DIAGNOSTICS.
c  
      if (IFLAG.eq.1) then

        PRECIS = 1000.0D+00*CONST(1)
        ERLUMP = 0.0D+00
        XOLD1 = X
        XOLD = X
        XOLDG = X
        XT = X

        DO I = 1, N
          YOLD(I) = Y(I)
          YOLDG(I) = Y(I)
        end do

        DO I = 1, 24
          CG(I) = 0.0D+00
        end do
        CG(1) = 1.0D+00
        CG(7) = 200.0D+00
        INDG = 2
        INDL = 2
        return

      end if
c  
c  IF IFLAG = 2   DETERMINE THE INITIAL STEPSIZE FOR
c  THE INTEGRATION PROPER.  WE CHOOSE THE SECOND STEP
c  TAKEN AND TERMINATE THE INTEGRATION BY SETTING X
c  EQUAL TO XEND. HSTART THEN HOLDS THE CURRENT STEPSIZE.
c
      if ( IFLAG .eq. 2 ) then
        NSTP = NSTP + 1
        HSTART = X - XOLD1
        XOLD1 = X
        if ( 2 .le. NSTP ) then
          X = XEND
        end if
        return
      end if
c  
c  IF IFLAG = 3   COMPILE STATISTICS.
c  
c  IF THE STEPSIZE AND, HENCE, THE ERROR REQUIREMENT WAS
c  TOO SMALL TO PERMIT AN EFFECTIVE ASSESSMENT AT THIS
c  PRECISION, CONTINUE THE INTEGRATION.  A LUMPED ERROR
c  ESTIMATE IS FORMED IN ERLUMP AND SEVERAL SMALL STEPS
c  ASSESSED AS ONE.
c  THE TEST FOR THE SIZE OF A LUMPED STEP IS MATCHED TO THE
c  MINIMUM STEPSIZE TEST IN 'TRUE' AND IS INTENDED TO ENSURE
c  (VERY CONSERVATIVELY) THAT ROUNDOFF EFFECTS ARE NEGLIGIBLE.
c  MAX-NORM IS USED IRRESPECTIVE OF THE VALUE OF NRMTYP IN /NSCOM1/.
c  THE LUMPED LOCAL ERROR IS TAKEN SIMPLY AS THE SUM OF THE
c  INDIVIDUAL LOCAL ERRORS.
c
      NSTP = NSTP + 1
      HLUMP = X - XOLD
      ERLUMP = ERLUMP + ERRBND
      XT = X
      YNORM = 0.0D+00
      DO I = 1, N
        YNORM = DMAX1(YNORM,DABS(YOLDG(I)),DABS(Y(I)))
      end do

      if ( HLUMP * ERRTOL .lt. YNORM * PRECIS ) then
        return
      end if
c
c  A SUFFICIENTLY LARGE LUMPED STEP HAS BEEN FORMED.
c  INCREMENT THE LUMPED STEP COUNT.
c  
  120 CONTINUE

      NSTL = NSTL + 1
c  
c  GLOBAL ASSESSMENT
c  SAVE COUNTERS THAT WILL BE AFFECTED BY 'TRUE' CALLS. SET MAX
c  STEPSIZE FOR GLOBAL SOLUTION TO X-XOLDG (DEFAULT VALUE IN TRUE IS
c  SIMPLY 2.)
c  CONTINUE TRUE GLOBAL SOLUTION TO CURRENT MESHPOINT AND
c  UPDATE MAX GLOBAL ERROR GEMX.
c  IF FAILURE OCCURS, RECORD POSITION IN XTRUE AND SKIP LOCAL
c  ASSESSMENT ALSO.
c  
      if (OPT.LT.2 .OR. INDG.LT.0) GO TO 240

      NNFCN = NFCN1
      HMIN = 10.0D+00*DMAX1(1.0D-30,CONST(1)*DABS(X))
      CG(3) = HMIN
      CG(6) = 1.1D+00*(X-XOLDG)
      TRUT0 = CLOCK(0.0D+00)

      CALL TRUE(N,FCN2,XOLDG,YOLDG,X,1.0D-02*ERRTOL,INDG,CG,NDIM,WG)

      TRUTIM = TRUTIM + CLOCK(TRUT0)
      CG(7) = CG(24) + 200.0D+00

      if ( INDG .lt. 0 ) then
        XTRUE = XOLDG
        GO TO 220
      end if

  140 GEMX = DMAX1(GEMX,DIFNRM(Y,YOLDG,N))
c  
c  LOCAL ASSESSMENT
c  OBTAIN THE LOCAL SOLUTION THROUGH THE PREVIOUS COMPUTED
c  MESH VALUE TO HIGHER ACCURACY THAN METHOD, PROVIDED NO
c  FAILURES HAVE OCCURRED IN PREVIOUS CALLS TO TRUE  (INDL.GE.0).
c  THE STARTING STEP FOR TRUE IS TAKEN AS .8 * THE LAST RECOMM-
c  ENDED STEPSIZE OF THE GLOBAL SOLUTION.
c  CHECK FOR A FAILURE THIS TIME AFTER THE
c  CALL TO TRUE.  COMPILE THE RELIABILITY STATISTICS.
c  
      if (OPT.LT.3 .OR. INDL.LT.0) GO TO 220

      DO I = 1, 9
        CL(I) = 0.0D+00
      end do
      INDL = 2
      CL(1) = 1.0D+00
      CL(3) = HMIN
      CL(4) = 0.8D+00*CG(14)
      CL(6) = 1.1D+00*(X-XOLD)
      CL(7) = 200.0D+00
      TRUT0 = CLOCK(0.0D+00)

      CALL TRUE(N,FCN2,XOLD,YOLD,X,1.0D-02*ERLUMP,INDL,CL,NDIM,WL)

      TRUTIM = TRUTIM + CLOCK(TRUT0)
      XTRUE = XOLD
c  
c  UPDATE STATISTICS
c  LESC RECORDS THE RATIO OF THE MAGNITUDE OF THE TRUE
c  LOCAL ERROR TO THE ASSUMED LOCAL ERROR BOUND.
c  LEMXSC RECORDS ITS MAXIMUM OVER THE RANGE.
c  NTRU COUNTS THE NO. OF LUMPED STEPS OF METHOD ON WHICH
c  LOCAL ASSESSMENT SUCCEEDED, SO AS TO ALLOW SUMMARY OF PARTIAL
c  RESULTS IF TRUE FAILS AT SOME POINT.
c
c  IF OPT=4, DO THE ANALYSIS OF THE LOCAL ERROR ESTIMATE VECTOR,
c  ERREST, BY FORMING THE SCALED ||ERROR|| IN ERREST.  IF LOCAL
c  EXTRAPOLATION IS DONE THIS IS LESC=||ERREST||/ERLUMP. IF NOT,
c  FORM YSTAR=LOCALLY EXTRAPOLATED SOLUTION AND IT IS THEN
c  ||YSTAR-YOLD||/ERLUMP. FORM A POINT ON THE SCATTER DIAGRAM
c  OF ERROR IN ERREST (VERT AXIS) VS. ERREST (HORIZ AXIS)
c  AND ENTER IT BY A CALL TO 'PLOT'.
c
c  FOR EVALUATING PERFORMANCE OF 'TRUE':
c  CALL TRUCHK(3,INFL)
c
      if ( 0 .le. INDL ) then

        LESC = DIFNRM ( Y, YOLD, N ) / ERLUMP
        LEMXSC = DMAX1 ( LEMXSC, LESC )

        if ( LESC .GT. 1.0D+00 ) then
          NDCV = NDCV + 1
        end if

        if ( LESC .GT. 5.0D+00 ) then
          NBAD = NBAD + 1
        end if
c
c  XTRAP=1 OR 0 ACCORDING AS THE USER HAS TOLD THE PACKAGE THAT
c  LOCAL EXTRAPOLATION IS OR IS NOT BEING DONE BY SOLVER:
c
        if ( OPT .EQ. 4 ) then
          if ( XTRAP .EQ. 0 ) then
            DO I = 1, N
              YSTAR(I) = Y(I) - ERREST(I)
            end do
            LEERSC = DIFNRM(YSTAR,YOLD,N) / ERLUMP
          ELSE
            LEERSC = LESC
          END IF
          ESTSC = DIFNRM(ERREST,ZERO,N) / ERLUMP
          CALL PLOT(ESTSC,LEERSC,1)
        END IF

        NTRU = NTRU + 1
c  
c  UPDATE MEMORY OF LAST COMPUTED VALUES.
c
        do I = 1, N
          YOLD(I) = Y(I)
        end do

      end if
c
c  RESTORE THE COUNTS AFFECTED BY 'TRUE' CALLS.
c
  220 NFCN1 = NNFCN
c
c  RE-INITIALIZE THE DATA PERTAINING TO A LUMPED STEP.
c
  240 ERLUMP = 0.0D+00
      XOLD = X
c
c  RETURN TO METHOD TO CONTINUE THE INTEGRATION.
c
      return
      end
      function str3 ( p )

c***********************************************************************
c
cc STR3 converts P to a 3 character integer percentage.
c
c  Discussion:
c
c    CONVERTS P (MEANT TO BE IN RANGE 0 TO 1) TO A 3 character
c    integer PERCENTAGE. P=0 BECOMES '   ', 0<P<1 BECOMES '  .',
c    P OUTSIDE RANGE BECOMES '***'
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c  Parameters:
c
      implicit none

      integer i
      integer j
      double precision p
      character * 3 str3

      character DIG(0:9)

      DATA DIG / '0', '1', '2', '3', '4', '5', '6',
     *                          '7', '8', '9'/

      DIG(0) = ' '

      if ( P .LT. 0.0D+00 .OR. P .GT. 1.0D+00 ) then
        STR3 = '***'
      else if (P .EQ. 0.0D+00 ) then
        STR3 = '   '
      else if ( P .LT. 0.01D+00 ) then
        STR3 = '  .'
      else

         DO j = 1, 3
            i = int ( p )
            p = p - dble ( i )
            STR3(J:J) = DIG(I)
            if (I.GT.0) DIG(0) = '0'
            P = 10.0D+00 * P
         end do

      end if

      return
      end
      subroutine true ( n, fcn, x, y, xend, tol, ind, c, nw, w )

c***********************************************************************
c
cc TRUE uses a Runge-Kutta scheme to estimate the true solution.
c
c  Modified:
c
c    16 March 2006
c
c  Reference:
c
c    W Enright and J Pryce,
c    Algorithm 648, NSDTST and STDTST,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 1, 1987, pages 28-34.
c
c THIS IS THE SAME AS THE DVERK CODE (DOCUMENTED BELOW) EXCEPT THAT
c THE SETTING OF THE MACHINE-DEPENDENT PARAMETERS C(10),C(11) HAS BEEN
c DONE BY CALLS TO THE 'CONST' FUNCTION OF THE DETEST PACKAGE.
c                             J D PRYCE AUG 1983
c
c  
c  
c     PURPOSE - THIS IS A RUNGE-KUTTA ROUTINE  BASED  ON  VERNER'S
c FIFTH AND SIXTH ORDER PAIR OF FORMULAS FOR FINDING APPROXIMATIONS TO
c THE SOLUTION OF  A  SYSTEM  OF  FIRST  ORDER  ORDINARY  DIFFERENTIAL
c EQUATIONS  WITH  INITIAL  CONDITIONS. IT ATTEMPTS TO KEEP THE GLOBAL
c ERROR PROPORTIONAL TO  A  TOLERANCE  SPECIFIED  BY  THE  USER.  (THE
c PROPORTIONALITY  DEPENDS  ON THE KIND OF ERROR CONTROL THAT IS USED,
c AS WELL AS THE DIFFERENTIAL EQUATION AND THE RANGE OF INTEGRATION.) 
c  
c     VARIOUS OPTIONS ARE AVAILABLE TO THE USER,  INCLUDING  DIFFERENT
c KINDS  OF  ERROR CONTROL, RESTRICTIONS ON STEP SIZES, AND INTERRUPTS
c WHICH PERMIT THE USER TO EXAMINE THE STATE OF THE  CALCULATION  (AND
c PERHAPS MAKE MODIFICATIONS) DURING INTERMEDIATE STAGES.              
c                                                                      
c     THE PROGRAM IS EFFICIENT FOR NON-STIFF SYSTEMS.  HOWEVER, A GOOD 
c VARIABLE-ORDER-ADAMS  METHOD  WILL PROBABLY BE MORE EFFICIENT IF THE 
c FUNCTION EVALUATIONS ARE VERY COSTLY.  SUCH A METHOD WOULD  ALSO  BE 
c MORE SUITABLE IF ONE WANTED TO OBTAIN A LARGE NUMBER OF INTERMEDIATE 
c SOLUTION VALUES BY INTERPOLATION, AS MIGHT BE THE CASE  FOR  EXAMPLE 
c WITH GRAPHICAL OUTPUT.                                               
c                                                                      
c                                    HULL-ENRIGHT-JACKSON   1/10/76    
c                                                                      
c                                                                        
c     USE - THE USER MUST SPECIFY EACH OF THE FOLLOWING                
c                                                                      
c     N  NUMBER OF EQUATIONS                                           
c                                                                      
c   FCN  NAME OF ROUTINE FOR EVALUATING FUNCTIONS - THE ROUTINE 
c           ITSELF MUST ALSO BE PROVIDED BY THE USER - IT SHOULD BE OF 
c           THE FOLLOWING FORM                                         
c              SUBROUTINE FCN(N, X, Y, YPRIME)                         
c              integer N                                               
c              double precision X, Y(N), YPRIME(N)                     
c                      *** ETC ***                                     
c           AND IT SHOULD EVALUATE YPRIME, GIVEN N, X AND Y            
c                                                                      
c     X  INDEPENDENT VARIABLE - INITIAL VALUE SUPPLIED BY USER         
c                                                                      
c     Y  DEPENDENT VARIABLE - INITIAL VALUES OF COMPONENTS Y(1), Y(2), 
c          ..., Y(N) SUPPLIED BY USER                                 
c                                                                      
c  XEND  VALUE OF X TO WHICH INTEGRATION IS TO BE CARRIED OUT - IT MAY 
c           BE LESS THAN THE INITIAL VALUE OF X                        
c                                                                      
c   TOL  TOLERANCE - THE ROUTINE ATTEMPTS TO CONTROL A NORM OF  THE 
c           LOCAL  ERROR  IN  SUCH  A  WAY  THAT  THE  GLOBAL ERROR IS 
c           PROPORTIONAL TO TOL. IN SOME PROBLEMS THERE WILL BE ENOUGH 
c           DAMPING  OF  ERRORS, AS WELL AS SOME CANCELLATION, SO THAT 
c           THE GLOBAL ERROR WILL BE LESS THAN TOL. ALTERNATIVELY, THE 
c           CONTROL   CAN   BE  VIEWED  AS  ATTEMPTING  TO  PROVIDE  A 
c           CALCULATED VALUE OF Y AT XEND WHICH IS THE EXACT  SOLUTION 
c           TO  THE  PROBLEM Y' = F(X,Y) + E(X) WHERE THE NORM OF E(X) 
c           IS PROPORTIONAL TO TOL.  (THE NORM  IS  A  MAX  NORM  WITH 
c           WEIGHTS  THAT  DEPEND ON THE ERROR CONTROL STRATEGY CHOSEN 
c           BY THE USER.  THE DEFAULT WEIGHT FOR THE K-TH COMPONENT IS 
c           1/MAX(1,ABS(Y(K))),  WHICH THEREFORE PROVIDES A MIXTURE OF 
c           ABSOLUTE AND RELATIVE ERROR CONTROL.)                      
c                                                                      
c   IND  INDICATOR - ON INITIAL ENTRY IND MUST BE SET EQUAL TO  EITHER 
c           1  OR  2. IF THE USER DOES NOT WISH TO USE ANY OPTIONS, HE 
c           SHOULD SET IND TO 1 - ALL THAT REMAINS FOR THE USER TO  DO 
c           THEN  IS  TO  DECLARE C AND W, AND TO SPECIFY NW. THE USER 
c           MAY ALSO  SELECT  VARIOUS  OPTIONS  ON  INITIAL  ENTRY  BY 
c           SETTING IND = 2 AND INITIALIZING THE FIRST 9 COMPONENTS OF 
c           C AS DESCRIBED IN THE NEXT SECTION.  HE MAY ALSO  RE-ENTER 
c           THE ROUTINE  WITH IND = 3 AS MENTIONED AGAIN BELOW. IN 
c           ANY EVENT, THE ROUTINE RETURNS WITH IND EQUAL TO        
c              3 AFTER A NORMAL RETURN                                 
c              4, 5, OR 6 AFTER AN INTERRUPT (SEE OPTIONS C(8), C(9))  
c              -1, -2, OR -3 AFTER AN ERROR CONDITION (SEE BELOW)      
c                                                                      
c     C  COMMUNICATIONS VECTOR - THE DIMENSION MUST BE GREATER THAN OR 
c           EQUAL TO 24, UNLESS OPTION C(1) = 4 OR 5 IS USED, IN WHICH 
c           CASE THE DIMENSION MUST BE GREATER THAN OR EQUAL TO N+30   
c                                                                      
c    NW  FIRST DIMENSION OF WORKSPACE W -  MUST  BE  GREATER  THAN  OR 
c           EQUAL TO N                                                 
c                                                                      
c     W  WORKSPACE MATRIX - FIRST DIMENSION MUST BE NW AND SECOND MUST 
c           BE GREATER THAN OR EQUAL TO 9                              
c                                                                      
c     THE ROUTINE  WILL  NORMALLY  RETURN  WITH  IND  =  3,  HAVING 
c REPLACED THE INITIAL VALUES OF X AND Y WITH, RESPECTIVELY, THE VALUE 
c OF XEND AND AN APPROXIMATION TO Y AT XEND.  THE ROUTINE  CAN  BE 
c CALLED  REPEATEDLY  WITH NEW VALUES OF XEND WITHOUT HAVING TO CHANGE 
c ANY OTHER ARGUMENT.  HOWEVER, CHANGES IN TOL, OR ANY OF THE  OPTIONS 
c DESCRIBED BELOW, MAY ALSO BE MADE ON SUCH A RE-ENTRY IF DESIRED.     
c                                                                      
c     THREE ERROR RETURNS ARE ALSO POSSIBLE, IN WHICH  CASE  X  AND  Y 
c WILL BE THE MOST RECENTLY ACCEPTED VALUES -                          
c     WITH IND = -3 THE ROUTINE WAS UNABLE  TO  SATISFY  THE  ERROR 
c        REQUIREMENT  WITH A PARTICULAR STEP-SIZE THAT IS LESS THAN OR 
c        EQUAL TO HMIN, WHICH MAY MEAN THAT TOL IS TOO SMALL           
c     WITH IND = -2 THE VALUE OF HMIN  IS  GREATER  THAN  HMAX,  WHICH 
c        PROBABLY  MEANS  THAT THE REQUESTED TOL (WHICH IS USED IN THE 
c        CALCULATION OF HMIN) IS TOO SMALL                             
c     WITH IND = -1 THE ALLOWED MAXIMUM NUMBER OF FCN EVALUATIONS  HAS 
c        BEEN  EXCEEDED,  BUT  THIS  CAN ONLY OCCUR IF OPTION C(7), AS 
c        DESCRIBED IN THE NEXT SECTION, HAS BEEN USED                  
c                                                                      
c     THERE ARE SEVERAL CIRCUMSTANCES THAT WILL CAUSE THE CALCULATIONS 
c TO  BE  TERMINATED,  ALONG WITH OUTPUT OF INFORMATION THAT WILL HELP 
c THE USER DETERMINE THE CAUSE OF  THE  TROUBLE.  THESE  CIRCUMSTANCES 
c INVOLVE  ENTRY WITH ILLEGAL OR INCONSISTENT VALUES OF THE ARGUMENTS, 
c SUCH AS ATTEMPTING A NORMAL  RE-ENTRY  WITHOUT  FIRST  CHANGING  THE 
c VALUE OF XEND, OR ATTEMPTING TO RE-ENTER WITH IND LESS THAN ZERO.    
c                                                                      
c     OPTIONS - IF THE ROUTINE IS ENTERED WITH IND = 1, THE FIRST 9 
c COMPONENTS OF THE COMMUNICATIONS VECTOR ARE INITIALIZED TO ZERO, AND 
c THE ROUTINE USES ONLY DEFAULT VALUES  FOR  EACH  OPTION.  IF  THE 
c ROUTINE  IS  ENTERED  WITH IND = 2, THE USER MUST SPECIFY EACH OF 
c THESE 9 COMPONENTS - NORMALLY HE WOULD FIRST SET THEM ALL  TO  ZERO, 
c AND  THEN  MAKE  NON-ZERO  THOSE  THAT  CORRESPOND TO THE PARTICULAR 
c OPTIONS HE WISHES TO SELECT. IN ANY EVENT, OPTIONS MAY BE CHANGED ON 
c RE-ENTRY  TO  THE ROUTINE  -  BUT IF THE USER CHANGES ANY OF THE 
c OPTIONS, OR TOL, IN THE COURSE OF A CALCULATION HE SHOULD BE CAREFUL 
c ABOUT  HOW  SUCH CHANGES AFFECT THE ROUTINE - IT MAY BE BETTER TO 
c RESTART WITH IND = 1 OR 2. (COMPONENTS 10 TO 24 OF C ARE USED BY THE 
c PROGRAM  -  THE INFORMATION IS AVAILABLE TO THE USER, BUT SHOULD NOT 
c NORMALLY BE CHANGED BY HIM.)                                         
c                                                                      
c  C(1)  ERROR CONTROL INDICATOR - THE NORM OF THE LOCAL ERROR IS  THE 
c           MAX  NORM  OF  THE  WEIGHTED  ERROR  ESTIMATE  VECTOR, THE 
c           WEIGHTS BEING DETERMINED ACCORDING TO THE VALUE OF C(1) -  
c              IF C(1)=1 THE WEIGHTS ARE 1 (ABSOLUTE ERROR CONTROL)    
c              IF C(1)=2 THE WEIGHTS ARE 1/ABS(Y(K))  (RELATIVE  ERROR 
c                 CONTROL)                                             
c              IF C(1)=3 THE  WEIGHTS  ARE  1/MAX(ABS(C(2)),ABS(Y(K))) 
c                 (RELATIVE  ERROR  CONTROL,  UNLESS ABS(Y(K)) IS LESS 
c                 THAN THE FLOOR VALUE, ABS(C(2)) )                    
c              IF C(1)=4 THE WEIGHTS ARE 1/MAX(ABS(C(K+30)),ABS(Y(K))) 
c                 (HERE INDIVIDUAL FLOOR VALUES ARE USED)              
c              IF C(1)=5 THE WEIGHTS ARE 1/ABS(C(K+30))                
c              FOR ALL OTHER VALUES OF C(1), INCLUDING  C(1) = 0,  THE 
c                 DEFAULT  VALUES  OF  THE  WEIGHTS  ARE  TAKEN  TO BE 
c                 1/MAX(1,ABS(Y(K))), AS MENTIONED EARLIER             
c           (IN THE TWO CASES C(1) = 4 OR 5 THE USER MUST DECLARE  THE 
c           DIMENSION OF C TO BE AT LEAST N+30 AND MUST INITIALIZE THE 
c           COMPONENTS C(31), C(32), ..., C(N+30).)                    
c                                                                      
c  C(2)  FLOOR VALUE - USED WHEN THE INDICATOR C(1) HAS THE VALUE 3    
c                                                                      
c  C(3)  HMIN SPECIFICATION - IF NOT ZERO, THE ROUTINE CHOOSES HMIN 
c           TO BE ABS(C(3)) - OTHERWISE IT USES THE DEFAULT VALUE      
c              10*MAX(DWARF,RREB*MAX(WEIGHTED NORM Y/TOL,ABS(X))),     
c           WHERE DWARF IS A VERY SMALL POSITIVE  MACHINE  NUMBER  AND 
c           RREB IS THE RELATIVE ROUNDOFF ERROR BOUND                  
c                                                                      
c  C(4)  HSTART SPECIFICATION - IF NOT ZERO, THE ROUTINE  WILL  USE 
c           AN  INITIAL  HMAG EQUAL TO ABS(C(4)), EXCEPT OF COURSE FOR 
c           THE RESTRICTIONS IMPOSED BY HMIN AND HMAX  -  OTHERWISE IT 
c           USES THE DEFAULT VALUE OF HMAX*(TOL)**(1/6)                
c                                                                      
c  C(5)  SCALE SPECIFICATION - THIS IS INTENDED TO BE A MEASURE OF THE 
c           SCALE OF THE PROBLEM - LARGER VALUES OF SCALE TEND TO MAKE 
c           THE METHOD MORE RELIABLE, FIRST  BY  POSSIBLY  RESTRICTING 
c           HMAX  (AS  DESCRIBED  BELOW) AND SECOND, BY TIGHTENING THE 
c           ACCEPTANCE REQUIREMENT - IF C(5) IS ZERO, A DEFAULT  VALUE 
c           OF  1  IS  USED.  FOR  LINEAR  HOMOGENEOUS  PROBLEMS  WITH 
c           CONSTANT COEFFICIENTS, AN APPROPRIATE VALUE FOR SCALE IS A 
c           NORM  OF  THE  ASSOCIATED  MATRIX.  FOR OTHER PROBLEMS, AN 
c           APPROXIMATION TO  AN  AVERAGE  VALUE  OF  A  NORM  OF  THE 
c           JACOBIAN ALONG THE TRAJECTORY MAY BE APPROPRIATE           
c                                                                      
c  C(6)  HMAX SPECIFICATION - FOUR CASES ARE POSSIBLE                  
c           IF C(6).NE.0 AND C(5).NE.0, HMAX IS TAKEN TO BE            
c              MIN(ABS(C(6)),2/ABS(C(5)))                              
c           IF C(6).NE.0 AND C(5).EQ.0, HMAX IS TAKEN TO BE  ABS(C(6)) 
c           IF C(6).EQ.0 AND C(5).NE.0, HMAX IS TAKEN TO BE            
c              2/ABS(C(5))                                             
c           IF C(6).EQ.0 AND C(5).EQ.0, HMAX IS GIVEN A DEFAULT  VALUE 
c              OF 2                                                    
c                                                                      
c  C(7)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS  -  IF  NOT  ZERO,  AN 
c           ERROR  RETURN WITH IND = -1 WILL BE CAUSED WHEN THE NUMBER 
c           OF FUNCTION EVALUATIONS EXCEEDS ABS(C(7))                  
c                                                                      
c  C(8)  INTERRUPT NUMBER  1  -  IF  NOT  ZERO,  THE ROUTINE  WILL 
c           INTERRUPT   THE  CALCULATIONS  AFTER  IT  HAS  CHOSEN  ITS 
c           PRELIMINARY VALUE OF HMAG, AND JUST BEFORE CHOOSING HTRIAL 
c           AND  XTRIAL  IN  PREPARATION FOR TAKING A STEP (HTRIAL MAY 
c           DIFFER FROM HMAG IN SIGN, AND MAY  REQUIRE  ADJUSTMENT  IF 
c           XEND  IS  NEAR) - THE ROUTINE RETURNS WITH IND = 4, AND 
c           WILL RESUME CALCULATION AT THE POINT  OF  INTERRUPTION  IF 
c           RE-ENTERED WITH IND = 4                                    
c                                                                      
c  C(9)  INTERRUPT NUMBER  2  -  IF  NOT  ZERO,  THE ROUTINE  WILL 
c           INTERRUPT   THE  CALCULATIONS  IMMEDIATELY  AFTER  IT  HAS 
c           DECIDED WHETHER OR NOT TO ACCEPT THE RESULT  OF  THE  MOST 
c           RECENT  TRIAL STEP, WITH IND = 5 IF IT PLANS TO ACCEPT, OR 
c           IND = 6 IF IT PLANS TO REJECT -  Y(*)  IS  THE  PREVIOUSLY 
c           ACCEPTED  RESULT, WHILE W(*,9) IS THE NEWLY COMPUTED TRIAL 
c           VALUE, AND W(*,2) IS THE UNWEIGHTED ERROR ESTIMATE VECTOR. 
c           THE ROUTINE  WILL  RESUME CALCULATIONS AT THE POINT OF 
c           INTERRUPTION ON RE-ENTRY WITH IND = 5 OR 6. (THE USER  MAY 
c           CHANGE IND IN THIS CASE IF HE WISHES, FOR EXAMPLE TO FORCE 
c           ACCEPTANCE OF A STEP THAT WOULD OTHERWISE BE REJECTED,  OR 
c           VICE VERSA. HE CAN ALSO RESTART WITH IND = 1 OR 2.)        
c
c  SUMMARY OF THE COMPONENTS OF THE COMMUNICATIONS VECTOR              
c                                                                      
c     PRESCRIBED AT THE OPTION       DETERMINED BY THE PROGRAM         
c           OF THE USER                                                
c                                                                      
c                                    C(10) RREB(REL ROUNDOFF ERR BND)  
c     C(1) ERROR CONTROL INDICATOR   C(11) DWARF (VERY SMALL MACH NO)  
c     C(2) FLOOR VALUE               C(12) WEIGHTED NORM Y             
c     C(3) HMIN SPECIFICATION        C(13) HMIN                        
c     C(4) HSTART SPECIFICATION      C(14) HMAG                        
c     C(5) SCALE SPECIFICATION       C(15) SCALE                       
c     C(6) HMAX SPECIFICATION        C(16) HMAX                        
c     C(7) MAX NO OF FCN EVALS       C(17) XTRIAL                      
c     C(8) INTERRUPT NO 1            C(18) HTRIAL                      
c     C(9) INTERRUPT NO 2            C(19) EST                         
c                                    C(20) PREVIOUS XEND               
c                                    C(21) FLAG FOR XEND               
c                                    C(22) NO OF SUCCESSFUL STEPS      
c                                    C(23) NO OF SUCCESSIVE FAILURES   
c                                    C(24) NO OF FCN EVALS             
c                                                                      
c  IF C(1) = 4 OR 5, C(31), C(32), ... C(N+30) ARE FLOOR VALUES        
c                                                                      
c
c                                                                      
c  AN OVERVIEW OF THE PROGRAM                                          
c                                                                      
c     BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES   
c ......ABORT IF IND OUT OF RANGE 1 TO 6                              
c  .     CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES  
c  .     CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2)                      
c  V........ABORT IF N.GT.NW OR TOL.LE.0                               
c  .        IF INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1)              
c  .           SET C(1) TO C(9) EQUAL TO ZERO                          
c  .        ELSE INITIAL ENTRY WITH OPTIONS (IND .EQ. 2)               
c  .           MAKE C(1) TO C(9) NON-NEGATIVE                          
c  .           MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED   
c  .        END IF                                                     
c  .        INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS            
c  .     CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3)                         
c .........ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT    
c  .        RE-INITIALIZE FLAG                                         
c  .     CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6)    
c  V        TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT.......  
c  .     END CASES                                                  .  
c  .  END INITIALIZATION, ETC.                                      .  
c  .                                                                V  
c  .  LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL STEP .  
c  .     STAGE 1 - PREPARE                                          .  
c***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT .  
c  .        CALC SLOPE (ADDING 1 TO NO OF FCN EVALS) IF IND .NE. 6  .  
c  .        CALC HMIN, SCALE, HMAX                                  .  
c***********ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX            .  
c  .        CALC PRELIMINARY HMAG                                   .  
c***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED.......RE-ENTRY.V  
c  .        CALC HMAG, XTRIAL AND HTRIAL                            .  
c  .     END STAGE 1                                                .  
c  V     STAGE 2 - CALC YTRIAL (ADDING 7 TO NO OF FCN EVALS)        .  
c  .     STAGE 3 - CALC THE ERROR ESTIMATE                          .  
c  .     STAGE 4 - MAKE DECISIONS                                   .  
c  .        SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6            .  
c***********INTERRUPT NO 2 IF REQUESTED....................RE-ENTRY.V  
c  .        IF STEP ACCEPTED (IND .EQ. 5)                              
c  .           UPDATE X, Y FROM XTRIAL, YTRIAL                         
c  .           ADD 1 TO NO OF SUCCESSFUL STEPS                         
c  .           SET NO OF SUCCESSIVE FAILURES TO ZERO                   
c**************RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND 
c  .        ELSE STEP NOT ACCEPTED (IND .EQ. 6)                        
c  .           ADD 1 TO NO OF SUCCESSIVE FAILURES                      
c**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN            
c  .        END IF                                                     
c  .     END STAGE 4                                                  
c  .  END LOOP                                                         
c  .                                                                   
c  BEGIN ABORT ACTION                                                  
c     OUTPUT APPROPRIATE  MESSAGE  ABOUT  STOPPING  THE  CALCULATIONS, 
c        ALONG WITH VALUES OF IND, N, NW, TOL, HMIN,  HMAX,  X,  XEND, 
c        PREVIOUS XEND,  NO OF  SUCCESSFUL  STEPS,  NO  OF  SUCCESSIVE 
c        FAILURES, NO OF FCN EVALS, AND THE COMPONENTS OF Y            
c     STOP                                                             
c  END ABORT ACTION                                                    
c                                                                      
      implicit none

      double precision TOL
      double precision X
      double precision XEND
      integer         IND, N, NW
      double precision C(24), W(NW,9), Y(N)
      EXTERNAL        FCN
      double precision TEMP
      integer         IOUT, K
      double precision            CONST
      EXTERNAL        CONST
      INTRINSIC       DABS, DMAX1, DMIN1, DSIGN
c
c  BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES *
c
c  ABORT IF IND OUT OF RANGE 1 TO 6.
c
      if (IND.LT.1 .OR. IND.GT.6) GO TO 1200
c
c  CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES
c
      GO TO (20,20,180,600,1080,1080) IND
c
c  CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2)
c  ABORT IF N.GT.NW OR TOL.LE.0
c
   20 if ( N .GT. NW .OR. TOL .LE. 0.0D+00 ) GO TO 1200
c
c  INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1)
c  SET C(1) TO C(9) EQUAL TO 0
c
      if ( IND .EQ. 2 ) GO TO 60

      DO K = 1, 9
        C(K) = 0.0D+00
      end do

      GO TO 140

   60 CONTINUE
c
c  INITIAL ENTRY WITH OPTIONS (IND .EQ. 2)
c  MAKE C(1) TO C(9) NON-NEGATIVE
c
      DO K = 1, 9
        C(K) = DABS(C(K))
      end do
c
c  MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED
c
      if ( C(1) .eq. 4.0D+00 .or. C(1) .eq. 5.0D+00 ) then

        DO k = 1, n
          c(k+30) = dabs ( c(k+30) )
        end do

      end if

  120 CONTINUE

  140 CONTINUE
c
c  INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS
c
      C(10) = CONST ( 1 )
      C(11) = CONST ( 2 )
c
c  SET PREVIOUS XEND INITIALLY TO INITIAL VALUE OF X
c
      C(20) = X
      DO k = 21, 24
        c(k) = 0.0D+00
      end do

      GO TO 200
c
c  CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3)
c  ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT
c
  180 continue

      if (C(21).NE.0.0D+00 .AND. (X.NE.C(20) .OR. XEND.EQ.C(20)))
     *    GO TO 1200
c
c  RE-INITIALIZE FLAG
c
      C(21) = 0.0D+00
      GO TO 200
c
c  CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6)
c  TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT
c  THIS HAS ALREADY BEEN HANDLED BY THE COMPUTED GO TO
c  END CASES 
c
  200 CONTINUE
c
c  END INITIALIZATION, ETC.
c
c  LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL  STEP 
c  UNTIL THE OCCURRENCE OF ONE OF THE FOLLOWING                   
c  (A) THE NORMAL RETURN (WITH IND .EQ. 3) ON REACHING XEND IN STAGE 4
c  (B) AN ERROR RETURN (WITH IND .LT. 0) IN STAGE 1 OR STAGE 4 
c  (C) AN INTERRUPT RETURN (WITH IND  .EQ.  4,  5  OR  6),  IF 
c      REQUESTED, IN STAGE 1 OR STAGE 4                        
c
  220 CONTINUE
c
c  STAGE 1 - PREPARE - DO CALCULATIONS OF  HMIN,  HMAX,  ETC., 
c  AND SOME PARAMETER  CHECKING,  AND  END  UP  WITH  SUITABLE 
c  VALUES OF HMAG, XTRIAL AND HTRIAL IN PREPARATION FOR TAKING 
c  AN INTEGRATION STEP.                                        
c
c  ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT.
c
      if ( C(7) .ne. 0.0D+00 .and. C(24) .ge. C(7) ) then
        IND = -1
        RETURN
      end if

  240 CONTINUE
c
c  CALCULATE SLOPE (ADDING 1 TO NO OF FCN EVALS) IF IND .NE. 6
c
      if ( IND .ne. 6 ) then
        CALL FCN(N,X,Y,W(1,1))
        C(24) = C(24) + 1.0D+00
      end if
c
c  CALCULATE HMIN - USE DEFAULT UNLESS VALUE PRESCRIBED
c
      C(13) = C(3)

      if (C(3).NE.0.0D+00 ) GO TO 500
c
c  CALCULATE DEFAULT VALUE OF HMIN
c  FIRST CALCULATE WEIGHTED NORM Y - C(12) - AS SPECIFIED
c  BY THE ERROR CONTROL INDICATOR C(1)
c
      TEMP = 0.0D+00
c
c  1: ABSOLUTE ERROR CONTROL - WEIGHTS ARE 1
c  2: RELATIVE ERROR CONTROL - WEIGHTS ARE 1/DABS(Y(K)) SO WEIGHTED NORM Y IS 1
c  3: WEIGHTS ARE 1/MAX(C(2),ABS(Y(K)))
c  4: WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K)))
c  5: WEIGHTS ARE 1/C(K+30)
c  DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K)))
c
      if ( C(1) .eq. 1.0D+00 ) then

        DO K = 1, N
          TEMP = DMAX1(TEMP,DABS(Y(K)))
        end do
        C(12) = TEMP

      else if ( C(1) .eq. 2.0D+00 ) then

        C(12) = 1.0D+00

      else if ( C(1) .eq. 3.0D+00 ) then

        DO K = 1, N
          TEMP = DMAX1(TEMP,DABS(Y(K)) / C(2) )
        end do
        C(12) = DMIN1(TEMP,1.0D+00)

      else if ( C(1) .eq. 4.0D+00 ) then

        DO K = 1, N
          TEMP = DMAX1(TEMP,DABS(Y(K))/C(K+30))
        end do
        C(12) = DMIN1(TEMP,1.0D+00)

      else if ( C(1) .eq. 5.0D+00 ) then

        DO K = 1, N
          TEMP = DMAX1(TEMP,DABS(Y(K))/C(K+30))
        end do
        C(12) = TEMP

      else

        DO K = 1, N
          TEMP = DMAX1 ( TEMP, DABS ( Y(K) ) )
        end do

        C(12) = DMIN1 ( TEMP, 1.0D+00 )

      end if

      C(13) = 10.0D+00 
     &  * DMAX1 ( C(11), C(10) * DMAX1 ( C(12) / TOL, DABS ( X ) ) )

  500 CONTINUE
c
c  CALCULATE SCALE - USE DEFAULT UNLESS VALUE PRESCRIBED
c
      if ( C(5) .EQ. 0.0D+00 ) then
        C(15) = 1.0D+00
      else
        C(15) = C(5)
      end if
c
c  CALCULATE HMAX - CONSIDER 4 CASES
c  CASE 1 BOTH HMAX AND SCALE PRESCRIBED
c  CASE 2 - HMAX PRESCRIBED, BUT SCALE NOT
c  CASE 3 - HMAX NOT PRESCRIBED, BUT SCALE IS
c  CASE 4 - NEITHER HMAX NOR SCALE IS PROVIDED
c
      if (C(6).NE.0.0D+00 .AND. C(5).NE.0.0D+00) then
        C(16) = DMIN1(C(6),2.D0/C(5))
      end if

      if (C(6).NE.0.0D+00 .AND. C(5).EQ.0.0D+00) then
        C(16) = C(6)
      end if

      if ( C(6) .EQ. 0.0D+00 .AND. C(5).NE.0.0D+00) then
        C(16) = 2.0D+00 / C(5)
      end if

      if (C(6).EQ.0.0D+00 .AND. C(5).EQ.0.0D+00) then
        C(16) = 2.0D+00
      end if
c
c  ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX.
c
      if ( C(13) .gt. C(16) ) then
        IND = -2
        RETURN
      end if

  520 CONTINUE
c
c  CALCULATE PRELIMINARY HMAG - CONSIDER 3 CASES
c
      if (IND.GT.2) GO TO 540
c
c  CASE 1 - INITIAL ENTRY - USE PRESCRIBED VALUE OF HSTART, IF
c  ANY, ELSE DEFAULT
c
      C(14) = C(4)
      if (C(4).EQ.0.0D+00) C(14) = C(16)*TOL**(1.0D+00/6.0D+00)
      GO TO 580

  540 continue

      if (C(23).GT.1.0D+00) GO TO 560
c
c  CASE 2 - AFTER A SUCCESSFUL STEP, OR AT MOST  ONE  FAILURE,
c  USE MIN(2, .9*(TOL/EST)**(1/6))*HMAG, BUT AVOID POSSIBLE
c  OVERFLOW. THEN AVOID REDUCTION BY MORE THAN HALF.
c
      TEMP = 2.0D+00 * C(14)

      if ( TOL .LT. ( 2.0D0 / 0.9D0 )**6 * C(19) ) then
        TEMP = 0.9D0 * ( TOL / C(19) )**(1.0/6.0) * C(14)
      end if

      C(14) = DMAX1 ( TEMP, 0.5D+00 * C(14) )
      GO TO 580
c
c  CASE 3 - AFTER TWO OR MORE SUCCESSIVE FAILURES
c
  560 CONTINUE

      C(14) = 0.5D+00 * C(14)

  580 CONTINUE
c
c  CHECK AGAINST HMAX
c
      C(14) = DMIN1(C(14),C(16))
c
c  CHECK AGAINST HMIN
c
      C(14) = DMAX1(C(14),C(13))
c
c  INTERRUPT NO 1 (WITH IND=4) IF REQUESTED
c
      if ( C(8) .ne. 0.0D+00 ) then
        IND = 4
        RETURN
      end if
c
c  RESUME HERE ON RE-ENTRY WITH IND .EQ. 4   ........RE-ENTRY..
c
  600 CONTINUE
c
c  CALCULATE HMAG, XTRIAL - DEPENDING ON PRELIMINARY HMAG, XEND.
c
c  DO NOT STEP MORE THAN HALF WAY TO XEND
c
      if ( C(14) .lt. DABS(XEND-X) ) then
        C(14) = DMIN1 ( C(14), 0.5D+00*DABS(XEND-X) )
        C(17) = X + DSIGN(C(14),XEND-X)
      else
        C(14) = DABS ( XEND - X )
        C(17) = XEND
      end if
c
c  CALCULATE HTRIAL
c
      C(18) = C(17) - X
c
c  END STAGE 1
c
c  STAGE 2 - CALCULATE YTRIAL (ADDING 7 TO NO OF  FCN  EVALS). 
c  W(*,2), ... W(*,8)  HOLD  INTERMEDIATE  RESULTS  NEEDED  IN 
c  STAGE 3. W(*,9) IS TEMPORARY STORAGE UNTIL FINALLY IT HOLDS 
c        * YTRIAL.
c
      TEMP = C(18) / 1398169080000.0D+00

      DO K = 1, N
         w(K,9) = Y(K) + TEMP*W(K,1)*233028180000.0D+00
      end do

      CALL FCN(N,X+C(18)/6.0D+00,W(1,9),W(1,2))

      DO K = 1, N
         w(K,9) = Y(K) + TEMP * (
     &      W(K,1) * 74569017600.0D+00
     &    + W(K,2) * 298276070400.0D+00 )
      end do

      CALL FCN(N,X+C(18)*(4.0D+00/15.0D+00),W(1,9),W(1,3))

      DO K = 1, N
         w(K,9) = Y(K) + TEMP * ( 
     &      W(K,1) * 1165140900000.D0
     &    - W(K,2) * 3728450880000.D0
     &    + W(K,3) * 3495422700000.D0 )
      end do

      CALL FCN(N,X+C(18)*(2.0D+00/3.0D+00),W(1,9),W(1,4))

      DO K = 1, N
         w(K,9) = Y(K) + TEMP * (
     &    - W(K,1) * 3604654659375.D0
     &    + W(K,2) * 12816549900000.D0
     &    - W(K,3) * 9284716546875.D0
     &    + W(K,4) * 1237962206250.D0 )
      end do

      CALL FCN(N,X+C(18)*(5.0D+00/6.0D+00),W(1,9),W(1,5))

      DO K = 1, N
         w(K,9) = Y(K) + TEMP*(W(K,1)*3355605792000.D0-W(K,2)
     *            *11185352640000.D0+W(K,3)*9172628850000.D0-W(K,4)
     *            *427218330000.D0+W(K,5)*482505408000.D0)
      end do

      CALL FCN(N,X+C(18),W(1,9),W(1,6))

      DO K = 1, N
         w(K,9) = Y(K) + TEMP*(-W(K,1)*770204740536.D0+W(K,2)
     *            *2311639545600.D0-W(K,3)*1322092233000.D0-W(K,4)
     *            *453006781920.D0+W(K,5)*326875481856.D0)
      end do

      CALL FCN(N,X+C(18)/15.0D+00,W(1,9),W(1,7))

      DO K = 1, N
         w(K,9) = Y(K) + TEMP*(W(K,1)*2845924389000.D0-W(K,2)
     *            *9754668000000.D0+W(K,3)*7897110375000.D0-W(K,4)
     *            *192082660000.D0+W(K,5)*400298976000.D0+W(K,7)
     *            *201586000000.D0)
      end do

      CALL FCN(N,X+C(18),W(1,9),W(1,8))
c
c  CALCULATE YTRIAL, THE EXTRAPOLATED APPROXIMATION AND STORE IN W(*,9)
c
      DO K = 1, N
         w(K,9) = Y(K) + TEMP*(W(K,1)*104862681000.D0+W(K,3)
     *            *545186250000.D0+W(K,4)*446637345000.D0+W(K,5)
     *            *188806464000.D0+W(K,7)*15076875000.D0+W(K,8)
     *            *97599465000.D0)
      end do
c
c  ADD 7 TO THE NO OF FCN EVALS
c
      C(24) = C(24) + 7.0D+00
c
c  END STAGE 2
c
c  STAGE 3 - CALCULATE THE ERROR ESTIMATE EST. FIRST CALCULATE 
c  THE  UNWEIGHTED  ABSOLUTE  ERROR  ESTIMATE VECTOR (PER UNIT 
c  STEP) FOR THE UNEXTRAPOLATED APPROXIMATION AND STORE IT  IN 
c  W(*,2).  THEN  CALCULATE THE WEIGHTED MAX NORM OF W(*,2) AS 
c  SPECIFIED BY THE ERROR  CONTROL  INDICATOR  C(1).  FINALLY, 
c  MODIFY  THIS RESULT TO PRODUCE EST, THE ERROR ESTIMATE (PER 
c  UNIT STEP) FOR THE EXTRAPOLATED APPROXIMATION YTRIAL. 
c
c  CALCULATE THE UNWEIGHTED ABSOLUTE ERROR ESTIMATE VECTOR
c
      DO K = 1, N
        w(K,2) = ( 
     &      w(K,1) * 8738556750.0D+00
     &    + W(K,3) * 9735468750.0D+00
     &    - W(K,4) * 9709507500.0D+00
     &    + W(K,5) * 8582112000.0D+00
     &    + W(K,6) * 95329710000.0D+00
     &    - W(K,7) * 15076875000.0D+00
     &    - W(K,8) * 97599465000.0D+00 )
     &    / 1398169080000.0D+00
      end do
c
c  CALCULATE THE WEIGHTED MAX NORM OF W(*,2) AS SPECIFIED BY
c  THE ERROR CONTROL INDICATOR C(1)
c
      TEMP = 0.0D+00
c
c  ABSOLUTE ERROR CONTROL
c
      if ( C(1) .eq. 1.0D+00 ) then

        do k = 1, n
          TEMP = DMAX1(TEMP,DABS(W(K,2)))
        end do
c
c  RELATIVE ERROR CONTROL
c
      else if ( C(1) .eq. 2.0D+00 ) then

        do k = 1, n
          TEMP = DMAX1(TEMP,DABS(W(K,2)/Y(K)))
        end do
c
c  WEIGHTS ARE 1/MAX(C(2),ABS(Y(K)))
c
      else if ( C(1) .eq. 3.0D+00 ) then

        do k = 1, n
          TEMP = DMAX1(TEMP,DABS(W(K,2))/DMAX1(C(2),DABS(Y(K))))
        end do
c
c  WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K)))
c
      else if ( C(1) .eq. 4.0D+00 ) then

        do k = 1, n
          TEMP = DMAX1(TEMP,DABS(W(K,2))/DMAX1(C(K+30),DABS(Y(K))))
        end do
c
c  WEIGHTS ARE 1/C(K+30)
c
      else if ( C(1) .eq. 5.0D+00 ) then

        do k = 1, n
          TEMP = DMAX1(TEMP,DABS(W(K,2)/C(K+30)))
        end do
c
c  DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K)))
c
      else

        do k = 1, n
          TEMP = DMAX1(TEMP,DABS(W(K,2))/DMAX1(1.0D+00,DABS(Y(K))))
        end do

      end if
c
c  CALCULATE EST - (THE WEIGHTED MAX NORM OF W(*,2))*HMAG*SCALE
c  - EST IS INTENDED TO BE A MEASURE OF THE ERROR  PER  UNIT
c  STEP IN YTRIAL
c
      C(19) = TEMP * C(14) * C(15)
c
c  END STAGE 3
c
c  STAGE 4 - MAKE DECISIONS.
c
c  SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6.
c
      IND = 5
      if (C(19).GT.TOL) IND = 6
c
c  INTERRUPT NO 2 IF REQUESTED
c
      if ( C(9) .ne. 0.0D+00 ) then
        return
      end if
c
c  RESUME HERE ON RE-ENTRY WITH IND .EQ. 5 OR 6   ...RE-ENTRY..
c
 1080 CONTINUE

      if (IND.EQ.6) GO TO 1140
c
c  STEP ACCEPTED (IND .EQ. 5), SO UPDATE X, Y FROM XTRIAL,
c  YTRIAL, ADD 1 TO THE NO OF SUCCESSFUL STEPS, AND SET
c  THE NO OF SUCCESSIVE FAILURES TO ZERO.
c
      X = C(17)
      do k = 1, n
        y(k) = w(k,9)
      end do
      C(22) = C(22) + 1.0D+00
      C(23) = 0.0D+00
c
c  RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND
c
      if ( x .eq. xend ) then
        ind = 3
        c(20) = xend
        c(21) = 1.0D+00
        return
      end if

      GO TO 1180

 1140 CONTINUE
c
c  STEP NOT ACCEPTED (IND .EQ. 6), SO ADD 1 TO THE NO OF
c  SUCCESSIVE FAILURES.
c
      C(23) = C(23) + 1.0D+00
c
c  ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN
c
      if ( c(14) .le. c(13) ) then
        ind = -3
        return
      end if

 1180 CONTINUE
c
c  END STAGE 4
c
      GO TO 220
c
c  END LOOP
c
c  BEGIN ABORT ACTION
c
 1200 CONTINUE

      IOUT = CONST(3)
      WRITE (IOUT,FMT=99999) IND, TOL, X, N, C(13), XEND, NW, C(16),
     *  C(20), C(22), C(23), C(24), (Y(K),K=1,N)

      STOP
c
c  END ABORT ACTION
c
99999 FORMAT (///
     *      'COMPUTATION STOPPED IN TRUE WITH THE FOLLOWING VALUES - '
     *       ,/'IND =',I4,5X,'TOL  =',1P,D13.6,5X,'X         =',1P,
     *       D22.15,/' N   =',I4,5X,'HMIN =',1P,D13.6,5X,'XEND      =',
     *       1P,D22.15,/' NW  =',I4,5X,'HMAX =',1P,D13.6,5X,
     *       'PREV XEND =',1P,D22.15,/ 14X,
     *       'NO OF SUCCESSFUL STEPS    =',0P,F8.0,/15X,
     *       'NO OF SUCCESSIVE FAILURES =',0P,F8.0,/15X,
     *       'NO OF FUNCTION EVALS      =',0P,F8.0,
     *       /'THE COMPONENTS OF Y ARE',//(' ',1P,5D24.15))
      end
