C      ALGORITHM 667, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 14, NO. 4, PP. 366-388.
      SUBROUTINE SIGMA1(N,X0,NSUC,IPRINT,XMIN,FMIN,NFEV,IOUT)
C
C  SIGMA1  IS A 'DRIVER' SUBROUTINE WHICH SIMPLY CALLS THE PRINCIPAL
C  SUBROUTINE  SIGMA  AFTER HAVING ASSIGNED DEFAULT VALUES TO A NUMBER
C  OF INPUT PARAMETERS OF  SIGMA , AND HAS THEREFORE A CONSIDERABLY
C  LOWER NUMBER OF INPUT PARAMETERS.
C  IT CAN BE USED AS A SIMPLE EXAMPLE OF HOW TO CALL  SIGMA , BUT
C  ALSO AS AN EASY-TO-USE DRIVER FOR THE AVERAGE USER, WHICH MAY FIND
C  IT EASIER TO CALL  SIGMA1  INSTEAD OF  SIGMA , THUS AVOIDING THE
C  TROUBLE OF ASSIGNING A VALUE TO ALL THE INPUT PARAMETERS OS  SIGMA .
C  ALL THE PARAMETERS IN THE DEFINITION OF  SIGMA1  HAVE THE SAME MEA-
C  NING AS IN  SIGMA .
C
C  THE USER OF  SIGMA1  MUST ONLY GIVE VALUES TO THE INPUT PARAMETERS
C      N, X0, NSUC, IPRINT
C  AND OBTAINS ON OUTPUT THE SAME OUTPUT PARAMETERS OF  SIGMA
C      XMIN, FMIN, NFEV, IOUT
C
C  WE RECALL HERE THE MEANING OF THE ABOVE PARAMETERS
C
C  N       IS THE PROBLEM DIMENSION (NUMBER OF VARIABLES)
C  X0      IS AN  N-VECTOR CONTAINING THE INITIAL VALUES OF THE
C          X-VARIABLES
C  NSUC    IS THE REQUESTED NUMBER OF SUCCESSFUL TRIALS
C          AFTER WHICH THE COMPUTATION IS STOPPED.
C  IPRINT  IS AN INDEX USED TO CONTROL THE AMOUNT OF PRINTED OUTPUT
C          BY CONTROLLING THE CALLS TO THE USER-SUPPLIED OUTPUT SUB-
C          ROUTINES  PTSEG  (END-OF-SEGMENT OUTPUT),  PTRIAL  (END-
C          OF-TRIAL OUTPUT), AND  PTKSUC (END-OF-TRIAL OUTPUT RELATED
C          TO THE COUNT OF SUCCESSFUL TRIALS), WHICH ARE DESCRIBED
C          BELOW.
C          IPRINT.LT.0   NO CALL TO THE OUTPUT SUBROUTINES
C          IPRINT.EQ.0   CALL ONLY  PTRIAL  AND  PTKSUC
C          IPRINT.GT.0   CALL ALL OUTPUT SUBROUTINES.
C  XMIN    IS AN N-VECTOR CONTAINING THE COORDINATES OF THE POINT
C          (OR POSSIBLY ONE OF THE POINTS) WHERE THE FINAL VALUE  FMIN
C          OF  FOPT  WAS FOUND.
C  FMIN    IS THE FINAL VALUE OF THE BEST CURRENT MINIMUM FUNCTION
C          VALUE  FOPT.
C  NFEV    IS THE TOTAL NUMBER OF FUNCTION EVALUATION (INCLUDING
C          THOSE USED FOR THE COMPUTATION OF DERIVATIVES, AND FOR
C          THE REJECTED TIME-INTEGRATION STEPS).
C  IOUT    IS THE INDICATOR OF THE STOPPING CONDITIONS, AS FOLLOWS
C          IF  IOUT = -99  A FATAL ERROR WAS DETECTED WHEN PERFOR-
C             MING SOME PRELIMINARY CHECKING OF THE INPUT DATA, AND
C             THE ALGORITHM WAS NOT EVEN STARTED
C          OTHERWISE THE ALGORITHM WAS STARTED, AND THE MEANING OF
C          IOUT  IS AS FOLLOWS
C          IF  IOUT.EQ.0  NO TRIAL HAD A UNIFORM STOP.
C          IF  IOUT.NE.0  THE SIGN AND THE MAGNITUDE OF  IOUT  HAVE
C          THE FOLLOWING MEANING
C          A) THE SIGN OF  IOUT  INDICATES IF THE BEST UNIFORM STOP
C          (I.E. THE ONE IN WHICH  FTFOPT  WAS OBTAINED)
C          IS TO BE CONSIDERED SUCCESSFUL (IOUT.GT.0) OR UNSUCCESSFUL
C          B) THE MAGNITUDE OF  IOUT  INDICATES WHICH ONE OF THE
C          NUMERICAL EQUALITY CRITERIA WAS SATISFIED IN THE BEST
C          UNIFORM STOP
C          IABS(IOUT).EQ.1  RELATIVE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.2  ABSOLUTE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.3  BOTH CRITERIA
C          (IOUT  IS THE FINAL VALUE OF THE INTERNAL PARAMETER  ISTOPT
C          (AN OUTPUT INDICATOR OF THE USER-SUPPLIED SUBROUTINE
C          PTRIAL )).
C          SUCCESS IS CLAIMED BY THE ALGORITHM IF  IOUT .GT. 0,
C          I.E. IF AT LEAST ONE OF THE TRIALS STOPPED SUCCESSFULLY
C          (I.E. WITH A POSITIVE VALUE OF THE TRIAL STOP INDICATOR
C          ISTOP ), AND THE SUCCESSFUL TRIALS (AS COUNTED BY  KSUC )
C          ARE CURRENTLY VALID.
C
      DOUBLE PRECISION X0,XMIN,FMIN
      DOUBLE PRECISION DX,EPS,H,TOLABS,TOLREL
      DOUBLE PRECISION VRMAX,VRMIN,XRMAX,XRMIN
      DIMENSION X0(N),XMIN(N)
      DIMENSION XRMIN(100),XRMAX(100)
      DATA VRMIN,VRMAX /-1.D4,1.D4/
      DATA NTRILD/50/
C
      H = 1.D-10
      EPS = 1.D0
      DX = 1.D-9
      IRAND = 0
      NTRAJ = 0
      ISEGBR = 0
      INKPBR = 0
      KPBR0 = 0
      NPMIN = 10
      NPMAX0 = 100
      INPMAX = 50
      NTRIAL = MAX0(NTRILD,5*NSUC)
      TOLREL = 1.D-3
      TOLABS = 1.D-6
      KPASCA = 10
      IF(N.GT.5)KPASCA = 300
      INHP = 1
      DO 1 IX = 1,N
         XRMIN(IX)=VRMIN
         XRMAX(IX)=VRMAX
 1    CONTINUE
C
      CALL SIGMA ( N, X0, H, EPS, DX,
     1             NTRAJ, ISEGBR, KPBR0, INKPBR,
     2             NPMIN, NPMAX0, INPMAX,
     3             NSUC, NTRIAL, TOLREL, TOLABS, XRMIN, XRMAX,
     4             KPASCA, IRAND, INHP, IPRINT,
     5             XMIN, FMIN, NFEV, IOUT )
C
      RETURN
      END
      SUBROUTINE SIGMA ( N, X0, H, EPS, DX,
     1                   NTRAJ, ISEGBR, KPBR0, INKPBR,
     2                   NPMIN, NPMAX0, INPMAX,
     3                   NSUC, NTRIAL, TOLREL, TOLABS, XRMIN, XRMAX,
     4                   KPASCA, IRAND, INHP, IPRINT,
     5                   XMIN, FMIN, NFEV, IOUT )
C
C  THE SUBROUTINE  SIGMA  IS THE PRINCIPAL SUBROUTINE OF THE PACKAGE
C  SIGMA, WHICH ATTEMPTS TO FIND A GLOBAL MINIMIZER OF A REAL VALUED
C  FUNCTION  F(X) = F(X1,...,XN)  OF  N  REAL VARIABLES  X1,...,XN.
C  THE ALGORITHM AND THE PACKAGE ARE DESCRIBED IN DETAIL IN THE TWO
C  PAPERS PUBLISHED IN THE SAME ISSUE OF THE A.C.M. TRANSACTIONS ON
C  MATHEMATICAL SOFTWARE, BOTH BY
C     F. ALUFFI-PENTINI, V. PARISI, F. ZIRILLI.
C  (1)  A GLOBAL MINIMIZATION ALGORITHM USING STOCHASTIC DIFFERENTIAL
C       EQUATIONS
C  (2)  ALGORITHM SIGMA, A STOCHASTIC-INTEGRATION GLOBAL MINIMIZATION
C       ALGORITHM.
C  THE SOFTWARE IMPLEMENTATION AND ITS USAGE ARE DESCRIBED IN (2).
C
C  METHOD
C
C  A GLOBAL MINIMIZER OF  F(X)  IS SOUGHT BY MONITORING THE VALUES OF
C  F  ALONG TRAJECTORIES GENERATED BY A SUITABLE (STOCHASTIC) DISCRE-
C  TIZATION OF A FIRST-ORDER STOCHASTIC DIFFERENTIAL EQUATION INSPIRED
C  BY STATISTICAL MECHANICS. STARTING FROM AN INITIAL POINT  X0 ,
C  X  IS UPDATED BY THE (STOCHASTIC) DISCRETIZATION STEP
C       X = X + DX1 + DX2
C  WHERE  DX1 = - H * GAM           (FIRST HALF-STEP)
C         DX2 = EPS * SQRT(H) * U   (SECOND HALF-STEP)
C  AND  H  IS THE TIME-INTEGRATION STEPLENGTH,
C  GAM/N  IS COMPUTED AS A FINITE-DIFFERENCE APPROXIMATION TO THE
C     DIRECTIONAL DERIVATIVE OF  F  ALONG AN ISOTROPICALLY RANDOM
C     DIRECTION,
C  EPS  IS A POSITIVE 'NOISE' COEFFICIENT, AND
C  U  IS A RANDOM SAMPLE FROM AN N-DIMENSIONAL GAUSSIAN DISTRIBUTION.
C  WE CONSIDER THE SIMULTANEOUS EVOLUTION OF A GIVEN FIXED NUMBER
C  NTRAJ  OF TRAJECTORIES DURING AN OBSERVATION PERIOD IN WHICH FOR
C  EACH TRAJECTORY  EPS  IS FIXED WHILE  H  AND THE SPATIAL DISCRETI-
C  ZATION INCREMENT  DX  FOR COMPUTING  GAM  ARE AUTOMATICALLY
C  ADJUSTED BY THE ALGORITHM.
C  AFTER EVERY OBSERVATION PERIOD ONE OF THE TRAJECTORIES IS DISCARDED,
C  ALL OTHER TRAJECTORIES CONTINUE UNPERTURBED, AND ONE OF THEM IS SE-
C  LECTED FOR BRANCHING, I.E. GENERATING A SECOND PERTURBED CONTI-
C  NUATION, WITH DIFFERENT STARTING  EPS  AND  DX  (AND THE SAME
C  'PAST HISTORY' OF THE FIRST).
C  THE SET OF SIMULTANEOUS TRAJECTORIES IS CONSIDERED A SINGLE TRIAL,
C  AND THE COMPLETE ALGORITHM IS A SET OF REPEATED TRIALS.
C  A TRIAL IS STOPPED, AT THE END OF AN OBSERVATION PERIOD, AND AFTER
C  HAVING DISCARDED THE WORST TRAJECTORY, IF ALL THE FINAL VALUES OF
C  F  FOR THE REMAINING TRAJECTORIES ARE EQUAL (WITHIN NUMERICAL TO-
C  LERANCES, AND POSSIBLY AT DIFFERENT POINTS  X) TO THEIR MINIMUM
C  VALUE  FTFMIN ('UNIFORM STOP AT THE LEVEL  FTFMIN ').
C  A UNIFORM STOP IS CONSIDERED SUCCESSFUL ONLY IF THE FINAL VALUE
C  FTFMIN  IS (NUMERICALLY) EQUAL TO THE CURRENT BEST MINIMUM  FOPT
C  FOUND SO FAR FROM ALGORITHM START.
C  A TRIAL IS ALSO ANYWAY STOPPED (UNSUCCESSFULLY) IF A GIVEN MAXIMUM
C  NUMBER  NPMAX  OF OBSERVATION PERIODS HAS ELAPSED.
C  TRIALS ARE REPEATED WITH DIFFERENT OPERATING CONDITIONS (INITIAL
C  POINT, MAX TRIAL LENGTH  NPMAX , SEED OF NOISE GENERATOR, POLICY
C  FOR CHOOSING THE STARTING  EPS  FOR THE PERTURBED CONTINUATION,
C  AND TRIAL-START VALUE OF  EPS ).
C  THE OUTCOMES OF ALL THE COMPLETED TRIALS ARE SUMMARIZED BY
C  FTFOPT (THE BEST CURRENT MIMNIMUM VALUE FOUND SO FAR FOR
C  FTFMIN ) AND BY THE CURRENT COUNT  KSUC  OF SUCCESSFUL TRIALS
C  (WHICH IS ZERO AT ALGORITHM START, AND IS UPDATED AT EVERY
C  SUCCESSFULLY-STOPPING TRIALS).
C  THE ALGORITHM IS STOPPED, AT THE END OF A TRIAL, IF A REQUESTED
C  COUNT  NSUC  OF SUCCESSFUL TRIALS HAS BEEN REACHED,
C  OR ANYWAY IF A GIVEN MAXIMUM NUMBER  NTRIAL  OF TRIALS HAS BEEN
C  REACHED.
C  SUCCESS IS CLAIMED IF THE CURRENT COUNT  KSUC  OF SUCCESSFUL TRIALS
C  IS AT LEAST ONE, AND IF SUCH TRIALS ARE CURRENTLY VALID.
C
C  CALL STATEMENT
C
C  THE CALL STATEMENT IS
C     CALL SIGMA ( N, X0, H, EPS, DX,
C                  NTRAJ, ISEGBR, KPBR0, INKPBR,
C                  NPMIN, NPMAX0, INPMAX,
C                  NSUC, NTRIAL, TOLREL, TOLABS, XRMIN, XRMAX,
C                  KPASCA, IRAND, INHP, IPRINT,
C                  XMIN, FMIN, NFEV, IOUT )
C
C  CALL PARAMETERS
C
C  INPUT PARAMETERS ARE THOSE IN LINES 1,3,4,5 OF THE CALL STATEMENT,
C  INPUT-OUTPUT PARAMETERS ARE THOSE IN LINE 2,
C  OUTPUT PARAMETERS ARE THOSE IN LINE 6.
C  NOTE THAT A NUMBER OF OTHER (INTERNAL) PARAMETERS CAN BE OBTAINED
C  BY MEANS OF THE USER-SUPPLIED OUTPUT SUBROUTINES  PTSEG, PTRIAL,
C  AND  PTKSUC, WHICH ARE DESCRIBED BELOW.
C
C  DESCRIPTION OF THE CALL PARAMETERS
C
C  N        IS THE PROBLEM DIMENSION (NUMBER OF VARIABLES)
C  X0       IS AN  N-VECTOR CONTAINING THE INITIAL VALUES OF THE
C          X-VARIABLES
C  H       IS THE INITIAL VALUE OF THE TIME-INTEGRATION STEPLENGTH.
C  EPS     IS THE INITIAL VALUE OF THE NOISE COEFFICIENT.
C  DX      IS THE INITIAL VALUE OF THE MAGNITUDE OF THE DISCRETIZATION
C          INCREMENT FOR COMPUTING THE FINITE-DIFFERENCE DERIVATIVES.
C  NTRAJ   IS THE NUMBER OF SIMULTANEOUS TRAJECTORIES.
C          (NOTE HOWEVER THAT IF THE INPUT VALUE IS ZERO,  NTRAJ  IS
C          SET TO A DEFAULT VALUE (NTRAJ = 7), AND IF THE INPUT VALUE
C          IS OTHERWISE OUTSIDE THE INTERVAL (3,20) NTRAJ IS SET TO
C          THE NEAREST EXTREME VALUE).
C  ISEGBR, KPBR0, INKPBR  DETERMINE, AT THE END OF AN OBSERVATION
C          PERIOD, WHICH ONE OF THE SIMULTANEOUS TRAJECTORIES
C          IS TO BE BRANCHED, AS FOLLOWS.
C          BRANCHING IS NORMALLY PERFORMED ON THE TRAJECTORY WHICH
C          OCCUPIES THE PLACE  ISEGBR  IN THE TRAJECTORY SELECTION OR-
C          DERING, EXCEPT AT (THE END OF) EXCEPTIONAL OBSERVATION
C          PERIODS, WHERE THE FIRST TRAJECTORY IN THE ORDERING IS
C          BRANCHED. EXCEPTIONAL BRANCHING OCCURS AT THE OBSERVATION
C          PERIODS NUMBERED  KP = KPBR0 + J*INKPBR, (J = 1,2,3,...).
C          THEREFORE  ISEGBR  SELECTS THE LEVEL (IN THE ORDERING) AT
C          WHICH NORMAL BRANCHING OCCURS, WHILE  KPBR0  AND  INKPBR
C          SELECT THE FIRST OCCURRENCE AND THE REPETITION FREQUENCY
C          OF THE EXCEPTIONAL OBSERVATION PERIODS.
C          (NOTE HOWEVER THAT IF ONE OF THE INPUT VALUES IS ZERO,
C          THE CORRESPONDING VARIABLE IS SET TO A DEFAULT VALUE
C          ISEGBR = INT((1+NTRAJ)/2), INKPBR = 10, KPBR0 = 3.
C          IF THE INPUT VALUE FOR  ISEGBR  IS OTHERWISE OUTSIDE THE
C          INTERVAL  (1,NTRAJ),  ISEGBR  IS SET TO THE NEAREST
C          EXTREME VALUE, AND IF KPBR0  HAS A VALUE NOT INSIDE THE
C          INTERVAL  (1,INKPBR), IT IS ASSIGNED THE SAME VALUE
C          MODULO  INKPBR).
C  NPMIN   IS THE MINIMUM DURATION OF A TRIAL, I.E. THE MINIMUM
C          NUMBER OF OBSERVATION PERIODS THAT SHOULD ELAPSE BEFORE
C          STARTING TO CHECK THE TRIAL STOPPING CRITERIA.
C  NPMAX0  IS THE MAXIMUM DURATION OF THE FIRST TRIAL, I.E. THE
C          VALUE, FOR THE FIRST TRIAL, OF MAXIMUM ACCEPTABLE
C          NUMBER  NPMAX  OF OBSERVATION PERIODS IN A TRIAL.
C  INPMAX  IS THE INCREMENT FOR  NPMAX , WHEN  NPMAX  IS VARIED
C          FROM ONE TRIAL TO THE FOLLOWING ONE.
C  NSUC    IS THE REQUESTED NUMBER OF SUCCESSFUL TRIALS
C          AFTER WHICH THE COMPUTATION IS STOPPED.
C  TOLREL  AND  TOLABS  ARE THE RELATIVE AND ABSOLUTE TOLERANCES
C          FOR STOPPING A SINGLE TRIAL.
C  XRMIN,  XRMAX  ARE N-VECTORS DEFINING THE ADMISIBLE REGION FOR
C          THE X-VALUES, WITHIN WHICH THE FUNCTION VALUES CAN BE
C          SAFELY COMPUTED.
C  KPASCA  IS THE MINIMUM NUMBER OF TRAJECTORY SEGMENTS THAT SHOULD
C          ELAPSE BEFORE THE RESCALING PROCEDURES ARE ACTIVATED.
C  IRAND   IS A CONTROL INDEX FOR THE INITIALIZATION OF THE RANDOM
C          NUMBER GENERATOR.
C          IRAND.GT.0    THE GENERATOR IS INITIALIZED, BEFORE STAR-
C                        TING THE TRIAL  KT, WITH SEED  IRAND+KT-1
C          IRAND.LE.0    THE GENERATOR IS INITIALIZED (WITH SEED 0)
C                        ONLY AT THE FIRST CALL OF SIGMA
C  INHP    IS A CONTROL INDEX FOR SELECTING THE NUMBER  NHP  OF TIME-
C          INTEGRATION STEPS FOR OBSERVATION PERIOD  KP (DURATION OF
C          TRIAL  KP) AS FOLLOWS (LOG IS BASE 2)
C          INHP=1  NHP = 1 + INT(LOG(KP))      ('SHORT' DURATION)
C          INHP=2  NHP = INT(SQRT(KP))         ('MEDIUM' DURATION)
C          INHP=3  NHP = KP                    ('LONG' DURATION)
C  IPRINT  IS AN INDEX USED TO CONTROL THE AMOUNT OF PRINTED OUTPUT
C          BY CONTROLLING THE CALLS TO THE USER-SUPPLIED OUTPUT SUB-
C          ROUTINES  PTSEG  (END-OF-SEGMENT OUTPUT),  PTRIAL  (END-
C          OF-TRIAL OUTPUT), AND  PTKSUC (END-OF-TRIAL OUTPUT RELATED
C          TO THE COUNT OF SUCCESSFUL TRIALS), WHICH ARE DESCRIBED
C          BELOW.
C          IPRINT.LT.0   NO CALL TO THE OUTPUT SUBROUTINES
C          IPRINT.EQ.0   CALL ONLY  PTRIAL  AND  PTKSUC
C          IPRINT.GT.0   CALL ALL OUTPUT SUBROUTINES.
C  XMIN    IS AN N-VECTOR CONTAINING THE COORDINATES OF THE POINT
C          (OR POSSIBLY ONE OF THE POINTS) WHERE THE FINAL VALUE  FMIN
C          OF  FOPT  WAS FOUND.
C  FMIN    IS THE FINAL VALUE OF THE BEST CURRENT MINIMUM FUNCTION
C          VALUE  FOPT.
C  NFEV    IS THE TOTAL NUMBER OF FUNCTION EVALUATION (INCLUDING
C          THOSE USED FOR THE COMPUTATION OF DERIVATIVES, AND FOR
C          THE REJECTED TIME-INTEGRATION STEPS).
C  IOUT    IS THE INDICATOR OF THE STOPPING CONDITIONS, AS FOLLOWS
C          IF  IOUT = -99  A FATAL ERROR WAS DETECTED WHEN PERFOR-
C             MING SOME PRELIMINARY CHECKING OF THE INPUT DATA, AND
C             THE ALGORITHM WAS NOT EVEN STARTED
C          OTHERWISE THE ALGORITHM WAS STARTED, AND THE MEANING OF
C          IOUT  IS AS FOLLOWS
C          IF  IOUT.EQ.0  NO TRIAL HAD A UNIFORM STOP.
C          IF  IOUT.NE.0  THE SIGN AND THE MAGNITUDE OF  IOUT  HAVE
C          THE FOLLOWING MEANING
C          A) THE SIGN OF  IOUT  INDICATES IF THE BEST UNIFORM STOP
C          (I.E. THE ONE IN WHICH  FTFOPT  WAS OBTAINED)
C          IS TO BE CONSIDERED SUCCESSFUL (IOUT.GT.0) OR UNSUCCESSFUL
C          B) THE MAGNITUDE OF  IOUT  INDICATES WHICH ONE OF THE
C          NUMERICAL EQUALITY CRITERIA WAS SATISFIED IN THE BEST
C          UNIFORM STOP
C          IABS(IOUT).EQ.1  RELATIVE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.2  ABSOLUTE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.3  BOTH CRITERIA
C          (IOUT  IS THE FINAL VALUE OF THE INTERNAL PARAMETER  ISTOPT
C          (AN OUTPUT INDICATOR OF THE USER-SUPPLIED SUBROUTINE
C          PTRIAL )).
C          SUCCESS IS CLAIMED BY THE ALGORITHM IF  IOUT .GT. 0,
C          I.E. IF AT LEAST ONE OF THE TRIALS STOPPED SUCCESSFULLY
C          (I.E. WITH A POSITIVE VALUE OF THE TRIAL STOP INDICATOR
C          ISTOP ), AND THE SUCCESSFUL TRIALS (AS COUNTED BY  KSUC )
C          ARE CURRENTLY VALID.
C
C
C  USER-SUPPLIED SUPROGRAMS
C
C  THE USER MUST PROVIDE THE FUNCTION  FUNCT  TO COMPUTE  F(X),
C  AND THE THREE OUTPUT SUBROUTINE  PTSEG, PTRIAL, PTKSUC .
C  THE CALLS TO THE OUTPUT SUBROUTINES ARE CONTROLLED BY  IPRINT
C  (INPUT PARAMETER TO  SIGMA).
C  A USER NOT INTERESTED IN USING ANY ONE OF THE OUTPUT SUBROUTINES
C  CAN SIMPLY PUT  IPRINT = -1.
C  SAMPLE VERSIONS OF THE OUTPUT SUBROUTINES ARE PROVIDED
C  WITH THE PACKAGE FOR THE BENEFIT OF THE AVERAGE USER.
C  IN THE FOLLOWING DESCRIPTION ALL NON-INTEGER ARGUMENTS ARE
C  DOUBLE PRECISION (INTEGER ARGUMENTS ARE INDICATED BY MEANS OF THE
C  FORTRAN IMPLICIT TYPE DEFINITION CONVENTION).
C
C     THE FUNCTION  FUNCT
C
C     FUNCT  MUST RETURN AS ITS VALUE THE VALUE AT  X  OF THE FUNCTION
C     TO BE MINIMIZED
C     THE DEFINITION STATEMENT IS
C        DOUBLE PRECISION FUNCTION  FUNCT (N, X)
C     WHERE
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM.
C     X   IS THE (INPUT) N-VECTOR CONTAINING THE COORDINATES OF THE
C         POINT  X  WHERE THE FUNCTION IS TO BE COMPUTED.
C
C     THE SUBROUTINE  PTSEG
C
C     PTSEG  IS CALLED (IF  IPRINT.GT.0) AT THE END OF EVERY OBSER-
C     VATION PERIOD.
C     THE DEFINITION STATEMENT IS
C         SUBROUTINE  PTSEG ( N, XPFMIN, FPFMIN, FPFMAX,
C                                         KP, NFEV, IPRINT )
C     WHERE
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM
C     FPFMIN, FPFMAX  ARE RESPECTIVELY THE MINIMUM AND THE MAXIMUM
C        AMONG THE VALUES OF  F(X) OBTAINED AT THE FINAL POINTS OF
C        THE TRAJECTORY SEGMENTS OF THE (JUST ELAPSED) OBSERVATION
C        PERIOD  KP.
C     XPFMIN   IS AN N-VECTOR CONTAINING THE COORDINATES OF THE
C        (FINAL) POINT (OR POSSIBLY ONE OF THE POINTS) WHERE
C        FPFMIN WAS OBTAINED.
C     KP   IS THE TOTAL NUMBER OF ELAPSED OBSERVATION PERIODS IN
C        THE CURRENT TRIAL.
C     NFEV   IS THE TOTAL NUMBER OF FUNCTION EVALUATIONS PERFORMED
C        FROM ALGORITHM START.
C
C     THE SUBROUTINE  PTRIAL
C
C     PTRIAL  IS CALLED (IF IPRINT.GE.0) AT THE END OF EVERY TRIAL.
C     THE DEFINITION STATEMENT IS
C         SUBROUTINE PTRIAL ( N, XOPT, FOPT,
C                             FTFMIN, FTFMAX, FTFOPT,
C                             ISTOP, ISTOPT, NFEV, KP, IPRINT )
C     WHERE
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM.
C     XOPT   IS AN N-VECTOR CONTINING THE COORDINATES OF THE
C        POINT (OR POSSIBLY ONE OF THE POINTS) WHERE THE CURRENT
C        MINIMUM  FOPT  WAS OBTAINED.
C     FOPT   IS THE CURRENT BEST MINIMUM VALUE FOUND FOR  F  FROM
C        ALGORITHM START ( FOPT  IS UPDATED WHENEVER A FUNCTION
C        VALUE IS COMPUTED).
C     FTFMIN, FTFMAX   ARE RESPECTIVELY THE MINIMUM AND THE MAXIMUM
C        AMONG THE VALUES OF  F(X)  OBTAINED AT THE FINAL POINTS OF
C        THE LAST TRAJECTORY SEGMENTS OF THE CURRENT TRIAL.
C     FTFOPT   IS THE CURRENT MINIMUM VALUE OF  FTFMIN  AMONG THE
C        TRIALS WHICH DID NOT STOP FOR REACHING THE MAXIMUM ALLOWED
C        NUMBER OF SEGMENTS (STOPPING INDICATOR  ISTOP = 0, SEE
C        BELOW).  FTFOPT  IS USED BY  SIGMA  TO COMPUTE  KSUC  (INPUT
C        PARAMETER TO THE OUTPUT SUBROUTINE  PTKSUK , SEE BELOW).
C     KP   IS THE TOTAL NUMBER OF ELAPSED OBSERVATION PERIODS IN
C        THE CURRENT TRIAL.
C     NFEV   IS THE TOTAL NUMBER OF FUNCTION EVALUATIONS PERFORMED
C        FROM ALGORITHM START.
C     ISTOP   IS THE INDICATOR OF THE STOPPING CONDITION OF THE TRIAL,
C        AS FOLLOWS
C        ISTOP = 0
C           THE MAXIMUM NUMBER  NPMAX  OF OBSERVATION PERIODS HAS
C           BEEN REACHED.
C        ISTOP.NE.0
C           ALL THE END-OF-SEGMENT VALUES OF  F(X) , (EXCEPT FOR THE
C           JUST DISCARDED SEGMENT) ARE CLOSE ENOUGH TO THEIR COMMON
C           MINIMUM VALUE  FPFMIN , WITH RESPECT TO AN ABSOLUTE OR
C           RELATIVE DIFFERENCE CRITERION, TO BE CONSIDERED NUMERI-
C           CALLY EQUAL.
C           THE ABSOLUTE VALUE AND THE SIGN OF  ISTOP  HAVE THE
C           FOLLOWING MEANING.
C              THE ABSOLUTE VALUE INDICATES WHICH DIFFERENCE
C                 CRITERION WAS SATISFIED
C                 1   RELATIVE DIFFERENCE CRITERION SATISFIED
C                 2   ABSOLUTE DIFFERENCE CRITERION SATISFIED
C                 3   BOTH CRITERIA SATISFIED
C              THE SIGN OF  ISTOP  INDICATES THE RELATIONSHIP BETWEEN
C                 THE END-OF-TRIAL VALUE  FPFMIN  AND THE CURRENT
C                 BEST MINIMUM VALUE  FOPT  (WHICH IS UPDATED WHEN-
C                 EVER A FUNCTION  VALUE IS COMPUTED
C                 ISTOP.GT.0
C                    FPFMIN  IS NUMERICALLY EQUAL (W.R.T. AT LEAST
C                    ONE OF THE ABOVE  DIFFERENCE CRITERIA) TO  FOPT
C                 ISTOP.LT.0
C                    FPFMIN  IS NOT EVEN NUMERICALLY EQUAL TO  FOPT
C                    (AND THEREFORE CANNOT BE CONSIDERED AS AN
C                    ACCEPTABLE GLOBAL MINIMUM).
C     ISTOPT   IS THE VALUE OF THE TRIAL STOPPING INDICATOR  ISTOP
C        CORRESPONDING TO THE (CURRENT OR PAST) TRIAL WHERE  FTFOPT
C        WAS OBTAINED, WITH THE SIGN WHICH IS UPDATED ACCORDING TO
C        THE COMPARISON BETWEEN  FTFOPT  AND THE PRESENT VALUE OF
C        FOPT , AS DESCRIBED ABOVE.
C        THE FINAL VALUE OF  ISTOP IS RETURNED BY  SIGMA  AS THE VALUE
C        OF THE OUTPUT INDICATOR IOUT OF THE ALGORITHM STOPPING CON-
C        DITIONS (WHENEVER THE ALGORITHM WAS STARTED, IOUT.NE.-99,
C        SEE ABOVE).
C
C     THE SUBROUTINE  PTKSUC
C
C     PTKSUC  IS CALLED ONLY AT THE END OF EVERY SUCCESSFUL TRIAL
C     SUCH THAT AN INCREMENT OCCURRED IN THE VALUE OF THE
C     CURRENT MAXIMUM  MSUC  AMONG ALL THE VALUES OF  KSUC  FROM
C     ALGORITHM START. A CALL TO  PTKSUC  THEREFORE PROVIDES
C     THE USER WITH THE OPERATIONALLY INTERESTING INFORMATION THAT
C     ALGORITHM STOP WOULD HAVE TAKEN PLACE, IF  NSUC  (INPUT
C     PARAMETER TO  SIGMA ) HAD BEEN GIVEN A LOWER VALUE, EQUAL TO
C     THE CURRENT  KSUC .
C     PTKSUC  IS CALLED ONLY IF  IPRINT.GE.0  AND  KSUC.LT.NSUC  .
C     THE DEFINITION STATEMENT IS
C         SUBROUTINE  PTKSUC ( KSUC )
C     WHERE  KSUC  IS THE CURRENT COUNT OF SUCCESSFUL TRIALS
C     (1 .LE. KSUC .LE. NSUC)
C
      DOUBLE PRECISION X0,H,EPS,DX,TOLREL,TOLABS
      DOUBLE PRECISION XRMIN,XRMAX,XMIN,FMIN
      DOUBLE PRECISION EPSAG,EPSAP,EPSC
      DOUBLE PRECISION EPSMAX,EPSR,F,FOPT,FTFMAX
      DOUBLE PRECISION FTFMIN,FTFOPT
C
      DOUBLE PRECISION X,HC,DXC,VMVT,EPSCO,VMCOR,VCOR
      DOUBLE PRECISION XRMIC,XRMAC,XOPT,FOPTC
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
C
      DIMENSION X0(N),XMIN(N),XRMIN(N),XRMAX(N)
C
      COMMON /DINCOM/ X(100,20),HC(20),DXC(20),VMVT(20,19),EPSCO(20),
     1    VMCOR(20),VCOR(20),XRMIC(100),XRMAC(100),XOPT(100),FOPTC,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJC,NTRAJR,
     3    ISEGBC,INKPBC,KPBR0C,NCF,IFEPC,INHPC
C
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
C  DATA FOR THE VARIATION OF NOISE COEFFICIENT
C
      DATA EPSR/1.D4/,EPSAP/10.D0/,EPSAG/1.D3/
      DATA EPSMAX/1.D15/
C
      IFEP = 1
C
C  INITIALIZE COMMON AREA  /DINCOM/
C
      CALL INIT(N,X0,H,EPS,DX,IRAND,F,
     1    NTRAJ,ISEGBR,INKPBR,KPBR0,INHP,IFEP,XRMIN,XRMAX,IOUT)
C
C  CHECK PARAMETER VALUES
C
      IF(NPMIN.LE.0.OR.NPMAX0.LT.0.OR.INPMAX.LE.0.OR.
     1    NSUC.LE.0.OR.NTRIAL.LE.0)IOUT = -99
      IF (IOUT.EQ.(-99))RETURN
C
C  INITIALIZE VARIABLES
C
      EPSC = EPS
      NPMAX = NPMAX0
      ISTOPT = 0
      ISTOP = 0
      ICCOM = (NTRIAL+NTRIAL+4)/5
      NFEV = 0
      NTES = NSUC
      ICTS = NSUC-1
      FTFOPT = F
C
C  START SERIES OF TRIALS
C
      DO 30 IC = 1,NTRIAL
C
C  SET INITIALIZATION INDEX FOR NOISE GENERATOR
C
         IS = 0
         IF(IRAND.GT.0)IS = IRAND+IC-1
C
C  INITIALIZE TRIAL
C
         IF(IC.GT.1.AND.IC.LE.ICCOM)CALL REINIT(N,X0,EPSC,IS,F,IFEP)
         IF(IC.GT.ICCOM)CALL REINIT(N,XMIN,EPSC,IS,FOPT,IFEP)
C
         FTFMIN = F
         FTFMAX = F
         NFEV = NFEV+1
C
C  PRINT INITIAL CONDITIONS OF TRIAL
C
         IF(IPRINT.GT.0)CALL PTSEG(N,X0,FTFMIN,FTFMAX,0,NFEV)
C
C
C  DEACTIVATE SCALING
C
         IF(KPASCA.GT.NPMAX.OR.N.LE.1)CALL NOSCA
C
C  INITIALIZE COMMON AREA  /SCALE/
C
         IF(KPASCA.LE.NPMAX.AND.N.GT.1)CALL INISCA(N,NTRAJ)
C
C  PERFORM A TRIAL
C
         CALL TRIAL(N,NPMIN,NPMAX,KPASCA,TOLREL,TOLABS,
     1       IPRINT,XMIN,FTFMIN,
     1       FTFMAX,NFEV,KP,ISTOP)
C
C
C     EVALUATE PAST TRIAL AND PREPARE NEXT TRIAL
C
C
C  SET TRIAL DURATION
C
         IF(ISTOP.EQ.0)NPMAX = NPMAX+INPMAX
C
C  RESET CURRENT NUMBER OF SUCCESSES REQUIRED BEFORE STOPPING
C  COMPUTE INDICATOR OF TRIAL STOPPING CONDITIONS
C  UPDATE BEST CURRENT VALUES OF TRIAL STOPPING INDICATOR AND
C  OF FUNCTION  F(X)  AT TRIAL STOP
C
         IF((FTFMIN.GT.FTFOPT.OR.ISTOP.EQ.0).AND.ISTOPT.NE.0)GO TO 10
         IF(ITOLCH(FTFOPT,FTFMIN,TOLREL,TOLABS).EQ.0)NTES = NSUC
C
         FTFOPT = FTFMIN
         ISTOPT = ISTOP
  10     CONTINUE
         ISTOPT = IABS(ISTOPT)
         CALL RCLOPT(N,XMIN,FOPT)
         IF(ITOLCH(FTFOPT,FOPT,TOLREL,TOLABS).EQ.0)ISTOPT = -ISTOPT
         IF(ITOLCH(FTFMIN,FOPT,TOLREL,TOLABS).EQ.0)ISTOP = -ISTOP
C
C  END-OF-TRIAL PRINT OUT
C
         IF(IPRINT.GE.0)
     1      CALL PTRIAL ( N, XOPT, FOPT,
     2                    FTFMIN, FTFMAX, FTFOPT,
     3                    ISTOP, ISTOPT, NFEV, KP, IPRINT )
C
C  UPDATE INITIAL VALUE OF NOISE COEFFICIENT FOR NEXT TRIAL
C
         IF (ISTOP.EQ.0)EPSC = EPSC/EPSR
         IF(ISTOP.GT.0)EPSC = EPSC*EPSAG
         IF(ISTOP.LT.0.AND.IC.LE.ICCOM)EPSC = EPSC*EPSAP
         IF(ISTOP.LT.0.AND.IC.GT.ICCOM)EPSC = EPSC/EPSR
C
C  UPDATE OPERATING CONDITIONS FOR SELECTING (IN THE NEXT TRIAL)
C  THE STARTING VALUE OF THE NOISE COEFFICIENT OF THE NEW TRAJECTORY
C  AFTER BRANCHING
C
         IF (ISTOP.EQ.0)IFEP = 1
         IF (ISTOP.NE.0)IFEP = 2
C
C  UPDATE, PRINT, AND CHECK TRIAL STOPPING CONDITIONS
C
         IF(ISTOP.GT.0.AND.ISTOPT.GT.0)NTES = NTES-1
         IF(NTES.LE.0.OR.ISTOPT.LE.0.OR.ISTOP.LE.0.OR.NTES.GT.ICTS.OR.
     1       IPRINT.LT.0) GO TO 20
         KSUC = NSUC-ICTS
         CALL PTKSUC(KSUC)
         ICTS = NTES-1
  20     CONTINUE
         IF(NTES.LE.0.AND.ISTOPT.GT.0.AND.ISTOP.GT.0)GO TO 40
C
C  CONSTRAIN NOISE COEFFICIENT WITHIN BOUNDS
C
         EPSC = DMIN1(EPSC,EPSMAX)
         IF(EPSC.LE.0.D0)EPSC = 1.D0
  30     CONTINUE
C  END OF SERIES OF TRIALS
C
  40  CONTINUE
C
C  INDICATOR OF STOPPING CONDITIONS
C
      IOUT = ISTOPT
      FMIN = FOPT
C
      RETURN
C
      END
      SUBROUTINE INIT(NX,X0,H0,EPS0,DX0,IRAND,F,N1,N2,N3,N4
     1    ,INH,IFE,XRI,XRA,IT)
C
C  INIT  PERFORMS THE INITIALIZATION OF THE FIRST TRIAL.
C  THE PART OF THE INITIALIZATION WHICH IS COMMON ALSO TO THE TRIALS
C  FOLLOWING THE FIRST ONE IS PERFORMED BY CALLING THE SUBROUTINE
C  REINIT.
C
      DOUBLE PRECISION X0,H0,EPS0,DX0,F,XRI,XRA
      DOUBLE PRECISION FUNCT
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION X0(NX),XRI(NX),XRA(NX)
      DATA NPMAX/100/
      DATA NTRAJM/20/
      DATA NTRAJ0/7/
      DATA INKPB0/10/
      DATA KPBR00/3/
C
C CHECK PARAMETER VALUES
C
      IT = 0
      IF (NX.GT.NPMAX.OR.NX.LT.1.OR.H0.LE.0.D0.OR.EPS0.LE.0.D0
     1    .OR.DX0.LE.0.D0) IT = -99
      IF (IT.EQ.(-99)) RETURN
C
C INITIALIZE SOME VARIABLES
C
      INHP = INH
      DO 10 IX = 1,NX
         XRMIN(IX) = XRI(IX)
         XRMAX(IX) = XRA(IX)
  10     CONTINUE
      CALL NOSCA
      NTRAJ = N1
      IF(NTRAJ.EQ.0)NTRAJ = NTRAJ0
      NTRAJ = MIN0(NTRAJ,NTRAJM)
      NTRAJ = MAX0(NTRAJ,3)
      N1 = NTRAJ
      ISEGBR = N2
      IF(ISEGBR.EQ.0)ISEGBR = (1+NTRAJ)/2
      ISEGBR = MIN0(ISEGBR,NTRAJ)
      ISEGBR = MAX0(ISEGBR,1)
      N2 = ISEGBR
      INKPBR = N3
      IF(INKPBR.EQ.0)INKPBR = INKPB0
      N3 = INKPBR
      KPBR0 = N4
      IF(KPBR0.EQ.0)KPBR0 = KPBR00
      KPBR0 = MOD(KPBR0,INKPBR)
      N4 = KPBR0
      NDIM = NX
      NTRAJR = NTRAJ-1
      NCF = 1
      F = FUNCT(NX,X0)
      CALL STOOPT(NX,X0,F)
      DO 20 ID = 1,NTRAJ
         H(ID) = H0
         DX(ID) = DX0
  20     CONTINUE
C
C INITIALIZE REMAINING VARIABLES
C
      CALL REINIT(NX,X0,EPS0,IRAND,F,IFE)
C
      RETURN
      END
      SUBROUTINE REINIT(NX,X0,EPS0,IRAND,F,IFE)
C
C  REINIT  PERFORMS THE INITIALIZATION OF ALL TRIALS FOLLOWING THE
C  FIRST TRIAL, AND PART OF THE INITIALIZATION OF THE FIRST TRIAL.
C
      DOUBLE PRECISION X0,EPS0,F
      DOUBLE PRECISION EPSV,G
      DOUBLE PRECISION CHAOS
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION X0(NX)
      DATA EPSV/1.D0/
C
C INITIALIZE RANDOM NOISE GENERATOR
C
      G = CHAOS(IRAND)
C
      IFEP = IFE
      KGEN = 0
      KTIM = 0
      DO 30 ID = 1,NTRAJ
         DO 10 IC = 1,NDIM
            X(IC,ID) = X0(IC)
  10        CONTINUE
         IE(ID) = 0
         DO 20 IT = 1,NTRAJR
            ISVT(ID,IT) = 1
            VMVT(ID,IT) = F
  20        CONTINUE
         EPS(ID) = EPS0*EPSV**(ISEGBR-ID)
         VMCOR(ID) = F
         VCOR(ID) = F
  30     CONTINUE
      RETURN
      END
      SUBROUTINE TRIAL(N,NPMIN,NPMAX,KPASCA,
     1    TOLREL,TOLABS,IPRINT,XMIN,FMIN,
     2    FMAX,NFEV,NR,ISTOP)
C
C  THE SUBROUTINE  TRIAL  GENERATES A TRIAL, I.E. A SET OF COMPLETE
C  SIMULTANEOUS TRAJECTORIES BY REPEATEDLY PERFORMING
C     A CALL TO THE SUBROUTINE  GENEVA  WHICH GENERATES THE SET OF
C        SIMULTANEOUS TRAJECTORY SEGMENTS CORRESPONDING TO THE CURRENT
C        OBSERVATION PERIOD, AND PERFORMS THE TRAJECTORY SELECTION
C     A (POSSIBLE) CALL TO THE SUBROUTINE  PTSEG  WHICH PERFORMS
C        END-OF-SEGMENT OUTPUT
C     A CHECK OF THE TRIAL STOPPING CRITERIA (WITH THE AID OF THE
C        FUNCTION  ITOLCH
C     A DECISION ABOUT ACTIVATING OR DEACTIVATING THE SCALING OF
C        THE VARIABLES (ACTIONS PERFORMED BY CALLING THE SUBROUTINES
C        ACTSCA  AND  NOSCA).
C
      DOUBLE PRECISION TOLREL,TOLABS,XMIN,FMIN,FMAX
      DIMENSION XMIN(N)
      DATA IRNF/7/
C
      DO 20 IR = 1,NPMAX
C
C ACTIVATE SCALING
C
         IF(IR.GE.KPASCA.AND.IR.GT.N*IRNF)CALL ACTSCA
         NR = IR
C
C GENERATE AND EVALUATE THE SIMULTANEOUS TRAJECTORY SEGMENTS
C PERIOD
C
         CALL GENEVA(N,XMIN,FMIN,FMAX,NFEV)
C
C PRINT RESULTS OF OBSERVATION PERIOD
C
         IF(IPRINT.GT.0)CALL PTSEG(N,XMIN,FMIN,FMAX,IR,NFEV)
C
C CHECK TRIAL STOPPING CONDITIONS
C
         IF(IR.LT.NPMIN)GO TO 10
         ISTOP = ITOLCH(FMAX,FMIN,TOLREL,TOLABS)
         IF(ISTOP.NE.0)GO TO 30
C
  10     CONTINUE
  20     CONTINUE
  30  CONTINUE
      CALL NOSCA
C
      RETURN
      END
      SUBROUTINE GENEVA(NX,XMIN,FMIN,FMAX,NCEF)
C
C  GENEVA  PERFORMS THE GENERATION AND THE FINAL PROCESSING AND
C  EVALUATION OF THE SET OF TRAJECTORY SEGMENTS CORRESPONDING TO
C  THE CURRENT OBSERVATION PERIOD.
C  GENEVA  UPDATES THE SCALING ARRAYS  DIST  AND  BIAS  BY CALLING
C             THE SUBROUTINES  SEGSCA  AND  UPDSCA
C          GENERATES THE TRAJECTORY SEGMENTS BY CALLING THE SUB-
C             ROUTINE  PERIOD
C          DETERMINES SOME END-OF-SEGMENT RESULTS  (FPFMIN, FPFMAX,
C             XPFMIN)  USING THE RESCALING CAPABILITIES OF THE SUB-
C             ROUTINES  SEGSCA  AND  VARSCA.
C
      DOUBLE PRECISION XMIN,FMIN,FMAX
      DOUBLE PRECISION FM
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XMIN(NX)
C
C UPDATE SCALING DATA
C
      DO 10 ID = 1,NTRAJ
         CALL SEGSCA(ID)
         CALL UPDSCA(NX,X(1,ID))
  10     CONTINUE
C
C GENERATE THE SIMULTANEOUS TRAJECTORY SEGMENTS
C
      CALL PERIOD
C
C EXTRACT AND RESCALE SOME FINAL VALUES
C
      FM = VCOR(1)
      FMAX = VCOR(1)
      IFM = 1
      DO 30 ID = 2,NTRAJ
         FMAX = DMAX1(FMAX,VCOR(ID))
         IF(VCOR(ID).GE.FM)GO TO 20
         FM = VCOR(ID)
         IFM = ID
  20     CONTINUE
  30     CONTINUE
      DO 40 IC = 1,NDIM
         XMIN(IC) = X(IC,IFM)
  40     CONTINUE
      FMIN = FM
      NCEF = NCF
      CALL SEGSCA(IFM)
      CALL VARSCA(NX,XMIN)
C
      RETURN
      END
      SUBROUTINE PERIOD
C
C  PERIOD  IS CALLED BY SUBROUTINE  GENEVA  TO PERFORM THE GENERATION
C  OF THE TRAJECTORY SEGMENTS.
C  PERIOD  - COMPUTES THE DURATION OF THE OBSERVATION PERIOD, I.E. THE
C              NUMBER  NHP  OF ACCEPTED INTEGRATION STEPS IN A PERIOD
C          - COMPUTES ALL THE SEGMENT STEPS BY REPEATEDLY CALLING THE
C              SUBROUTINE  STEP
C          - PERFORMS THE SEGMENT SELECTION BY CALLING THE SUBROUTINE
C               BRASI
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
C DETERMINE DURATION OF OBSERVATION PERIOD
C (NUMBER OF TIME INTEGRATION STEPS)
C
      KGEN = KGEN+1
      NKGEN = 1
      IF(INHP.EQ.1)NKGEN=
     1    INT(SNGL(DLOG(DBLE(FLOAT(KGEN))+.5D0)/DLOG(2.D0)))+1
      IF(INHP.EQ.2)NKGEN = INT(SNGL(DSQRT(DBLE(FLOAT(KGEN))+.5D0)))
      IF(INHP.EQ.3)NKGEN = KGEN
C
C PERFORM THE REQUIRED NUMBER OF INTEGRATION STEPS
C (FOR ALL TRAJECTORIES)
C
      DO 10 IK = 1,NKGEN
C
C PERFORM A SINGLE INTEGRATION STEP
C (FOR ALL TRAJECTORIES)
C
         CALL STEP(IK)
  10     CONTINUE
C
C MANAGE THE TRAJECTORY BRANCHING PROCESS
C
      CALL BRASI
      RETURN
      END
      SUBROUTINE BRASI
C
C  BRASI  PERFORMS THE SELECTION PROCESS FOR THE TRAJECTORIES
C  BRASI  - UPDATES THE DATA ABOUT THE PAST TRAJECTORIES
C         - ASKS FOR THE TRAJECTORY-SELECTION ORDERING BY CALLING THE
C              SUBROUTINE  ORDER
C         - DISCARDS ONE OF THE TRAJECTORIES
C         - PERFORMS BRANCHING ON ONE OF THE REMAINING TRAJECTORIES
C         - MOVES THE DATA OF THE UNPERTURBED CONTINUATION
C              TO THE POSITION OF THE PERTURBED CONTINUATION
C         - CALLS THE SUBROUTINE  COMPAS  TO EXAMINE DATA ABOUT PAST
C              HISTORY OF THE TRAJECTORIES AND DISCARD IRRILEVANT DATA
C
      DOUBLE PRECISION DFAC,DLEPMX,DLFACL
      DOUBLE PRECISION EFAC,EPSMAX,FACG
      DOUBLE PRECISION CHAOS
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      DATA FACG/10.D0/
      DATA EFAC/.5D0/
      DATA DFAC/1.D3/
      DATA EPSMAX/1.D15/
      DATA DLFACL /0.301029995663981194D0/
      DATA DLEPMX /25.D0/
C
C UPDATE PAST HISTORY DATA
C
      DO 10 ID = 1,NTRAJ
         VMVT(ID,NTRAJR) = VMCOR(ID)
         ISVT(ID,NTRAJR) = ID
  10     CONTINUE
C
C OBTAIN TRAJECTORY-SELECTION ORDERING
C
      CALL ORDER(IP,IM,IU)
C
C DECIDE WHICH TRAJECTORY IS TO BE BRANCHED
C
      IF(MOD(KGEN,INKPBR).EQ.KPBR0)IM = IP
C
C PERFORM BRANCHING
C
      DO 20 IC = 1,NDIM
         X(IC,IU) = X(IC,IM)
  20     CONTINUE
      H(IU) = H(IM)
      IE(IU) = IE(IM)
      DO 30 IT = 1,NTRAJR
         ISVT(IU,IT) = ISVT(IM,IT)
         VMVT(IU,IT) = VMVT(IM,IT)
  30     CONTINUE
      EPS(IU) = EPS(IM)
      DX(IU) = DX(IM)
      VCOR(IU) = VCOR(IM)
      DO 40 ID = 1,NTRAJ
         VMCOR(ID) = VCOR(ID)
  40     CONTINUE
C
C UPDATE PAST-HISTORY-DATA MATRICES
C
      CALL COMPAS
C
C UPDATE SCALING DATA
C
      CALL MOVSCA(IU,IM)
C
C UPDATE NOISE COEFFICIENT VALUES
C
      IF (EPS(IU).LE.0.D0) GO TO 50
      IF(IFEP.EQ.2)
     1    EPS(IU) = EPS(IU)*FACG**(CHAOS(0)-EFAC)
      IF(IFEP.EQ.1)
     1    EPS(IU) = FACG**(DMIN1(DLEPMX,
     1    DLOG10(EPS(IU))+(CHAOS(-1)-EFAC)*DLFACL))
      EPS(IU) = DMIN1(EPS(IU),EPSMAX)
  50  CONTINUE
C
C UPDATE MAGNITUDE OF SPATIAL DISCRETIZATION INCREMENT
C
      DX(IU) = DX(IU)*DFAC**CHAOS(0)
      RETURN
      END
      SUBROUTINE ORDER(IP,IM,IU)
C
C  ORDER  COMPARES THE TRAJECTORIES FROM THE POINT OV VIEW OF PAST
C  HISTORY (BY CALLING THE FUNCTION  IPREC ) AND FROM THE POINT OF VIEW
C  OF THEIR CURRENT VALUE OF  EPS  (BY CALLING THE FUNCTION  IPRECE)
C  AND PROVIDES THE TRAJECTORY ORDERING NEEDED FOR SELECTING THE TRA-
C  JECTORY WHICH IS TO BE BRANCHED.
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION IORD(20)
      DATA KP/0/
      IR = 0
C
C ASSIGN INITIAL ORDERING
C
  10  CONTINUE
      DO 20 I = 1,NTRAJ
         IORD(I) = I
  20     CONTINUE
C
C SORT TRAJECTORIES ...
C
  30  CONTINUE
      IR = IR+1
C
      DO 60 I = 1,NTRAJR
         I1 = I+1
         DO 50 J = I1,NTRAJ
            K1 = IORD(I)
            K2 = IORD(J)
C
C ... ACCORDING TO PAST HISTORY ...
C
            IF(IR.NE.2)KP = IPREC(K1,K2)
            IF(KP.EQ.0.AND.IR.EQ.1)GO TO 10
C
C ... OR ACCORDING TO NOISE LEVEL
C
            IF(IR.EQ.2)KP = IPRECE(K1,K2)
            IF(KP.EQ.0)GO TO 40
            KM = K1+K2-KP
            IORD(I) = KP
            IORD(J) = KM
  40        CONTINUE
  50        CONTINUE
  60     CONTINUE
      IF(IR.EQ.2)GO TO 30
C
C RETURN THE INDICES OF THE SEGMENTS WHICH IN THE ORDERING
C OCCUPY THE FIRST, THE LAST, AND A SUITABLE MEDIUM LEVEL POSITION
C
      IP = IORD(1)
      IU = IORD(NTRAJ)
      IM = IORD(ISEGBR)
C
      RETURN
      END
      INTEGER FUNCTION IPREC(ID1,ID2)
C
C  IPREC  DETERMINES THE PRECEDENCE RELATION BETWEEN TWO TRAJECTORIES
C  BASED ON THE PAST HISTORY DATA
C
      DOUBLE PRECISION VM1,VM2
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      VM1 = VMVT(ID1,NTRAJR)
      VM2 = VMVT(ID2,NTRAJR)
      DO 10 IIT = 1,NTRAJR
         IT = 1+NTRAJR-IIT
         IF(ISVT(ID1,IT).EQ.ISVT(ID2,IT))GO TO 20
         VM1 = DMIN1(VM1,VMVT(ID1,IT))
         VM2 = DMIN1(VM2,VMVT(ID2,IT))
  10     CONTINUE
  20  CONTINUE
      IPREC = 0
      IF(VM2.LT.VM1)IPREC = ID2
      IF(VM2.GT.VM1)IPREC = ID1
C
      RETURN
      END
      INTEGER FUNCTION IPRECE(ID1,ID2)
C
C  IPRECE  DETERMINES THE PRECEDENCE RELATION BETWEEN TWO TRAJECTORIES
C  BASED ON THEIR CURRENT VALUE OF  EPS
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      IPRECE = 0
      IF(KGEN.GT.ISEGBR*INKPBR)GO TO 10
      IF(EPS(ID2).LT.EPS(ID1))IPRECE = ID1
      IF(EPS(ID2).GT.EPS(ID1))IPRECE = ID2
      RETURN
C
  10  CONTINUE
      IF(EPS(ID2).LT.EPS(ID1))IPRECE = ID2
      IF(EPS(ID2).GT.EPS(ID1))IPRECE = ID1
C
      RETURN
      END
      SUBROUTINE COMPAS
C
C  COMPAS  TAKES CARE OF THE STORAGE OF PAST HISTORY DATA, DISCARDING
C  ALL DATA NOT NEEDED BY THE ONLY USER OF SUCH DATA, THE FUNCTION
C  IPREC .
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      IT1 = 1
      DO 30 IT = 2,NTRAJR
         IT1 = IT-1
         DO 10 ID = 1,NTRAJ
            IF(ISVT(ID,IT).NE.ISVT(ID,IT1))GO TO 20
  10        CONTINUE
         GO TO 120
  20     CONTINUE
  30     CONTINUE
      IT1 = 1
      NCV = 0
      DO 50 ID = 1,NTRAJ
         DO 40 IDD = 1,NTRAJ
            IF(ISVT(ID,IT1).EQ.ISVT(IDD,IT1))NCV = NCV+1
  40        CONTINUE
  50     CONTINUE
      DO 80 IT = 2,NTRAJR
         IT1 = IT-1
         NCN = 0
         DO 70 ID = 1,NTRAJ
            DO 60 IDD = 1,NTRAJ
               IF(ISVT(ID,IT).EQ.ISVT(IDD,IT))NCN = NCN+1
  60           CONTINUE
  70        CONTINUE
         IF(NCN.EQ.NCV)GO TO 110
         NCV = NCN
  80     CONTINUE
      DO 90 ID = 1,NTRAJ
         IF(ISVT(1,1).NE.ISVT(ID,1))GO TO 100
  90     CONTINUE
      IT = 2
      GO TO 140
 100  CONTINUE
 110  CONTINUE
 120  CONTINUE
      IT = IT1+1
      DO 130 ID = 1,NTRAJ
         VMVT(ID,IT) = DMIN1(VMVT(ID,IT),VMVT(ID,IT1))
 130     CONTINUE
 140  CONTINUE
      DO 160 ITC = IT,NTRAJR
         ITM = ITC-1
         DO 150 ID = 1,NTRAJ
            VMVT(ID,ITM) = VMVT(ID,ITC)
            ISVT(ID,ITM) = ISVT(ID,ITC)
 150        CONTINUE
 160     CONTINUE
C
      RETURN
      END
      SUBROUTINE STEP(IK)
C
C  STEP  PERFORMS A SINGLE TIME-INTEGRATION STEP FOR EACH ONE OF THE
C  SIMULTANEOUS TRAJECTORIES BY REPEATEDLY CALLING THE SUBROUTINE  SSTEP
C
      DOUBLE PRECISION F
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      DOUBLE PRECISION XID,HID,EPSID,DXID
      DIMENSION XID(100)
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      KTIM = KTIM+1
C
C LOOP ON THE SIMULTANEOUS TRAJECTORY SEGMENTS
      DO 30 ID = 1,NTRAJ
C
C INFORM THE SCALING SUBPROGRAMS (VIA THE COMMON AREA /SCALE/)
C THAT SCALING OPERATIONS ARE TO BE PERFORMED ON SEGMENT  ID .
         CALL SEGSCA(ID)
         F = VCOR(ID)
C
C PERFORM A TIME-INTEGRATION STEP ON SEGMENT  ID .
         KA = KTIM
         NX = NDIM
         DO 10 IX = 1,NX
            XID(IX) = X(IX,ID)
  10        CONTINUE
         HID = H(ID)
         EPSID = EPS(ID)
         DXID = DX(ID)
         IEID = IE(ID)
C
         CALL SSTEP(KA,NX,XID,HID,EPSID,DXID,IEID,F)
C
         DO 20 IX = 1,NX
            X(IX,ID) = XID(IX)
  20        CONTINUE
         H(ID) = HID
         EPS(ID) = EPSID
         DX(ID) = DXID
         IE(ID) = IEID
         VCOR(ID) = F
         VMCOR(ID) = DMIN1(VMCOR(ID),F)
         IF(IK.EQ.1)VMCOR(ID) = F
         IF(KTIM.EQ.1)IE(ID) = 0
  30     CONTINUE
C
      RETURN
      END
      SUBROUTINE SSTEP(KTIM,NDIM,X,H,EPS,DX,IE,F)
C
C  THE BASIC TIME-INTEGRATION STEP FOR A GIVEN TRAJECTORY IS PERFORMED
C  BY THE SUBROUTINE  STEP  WHICH
C   - CALLS THE FUNCTION  FUNCT0  TO COMPUTE THE VALUE OF  F
C   - CALLS THE SUBROUTINE  UNITRV  TO COMPUTE THE RANDOM DIRECTION
C        ALONG WHICH THE DIRECTIONAL DERIVATIVE  GAM/N  IS TO BE
C        COMPUTED
C   - CALLS THE SUBROUTINE  DERFOR (OR DERCEN ) TO COMPUTE THE FORWARD-
C        (OR CENTRAL-) DIFFERENCES DIRECTIONAL DERIVATIVE  GAM/N
C   - CALLS THE SUBROUTINE  NEWH  TO ACCEPT OR REJECT THE FIRST HALF-
C        STEP AND OBTAIN AN UPDATED VALUE FOR  H
C   - CALLS THE SUBROUTINE  CUMSCA  TO UPDATE THE CUMULATED SCALING DATA
C   - UPDATES THE SPATIAL DISCRETIZATION INCREMENT  DX  BASED ON THE
C        RESULTS OF CALLING THE FUNCTION  ITOLCH
C   - CALLS THE SUBROUTINE  GAUSRV  TO PERFORM THE SECOND HALF-STEP.
C
      DOUBLE PRECISION X,H,EPS,DX,F
      DOUBLE PRECISION DFDX,DFDXV,DXMAX,DXMIN
      DOUBLE PRECISION EPSR,FS,FV,FVS,HMINS,HR,HS
      DOUBLE PRECISION RCD,RDX,STF,TOLABS,TOLRA,TOLRI
      DOUBLE PRECISION W,XP
      DOUBLE PRECISION FUNCT0
      DIMENSION X(NDIM)
      DIMENSION W(100),XP(100)
      DATA RDX/1.D-4/
      DATA DXMIN/1.D-35/
      DATA DXMAX/1.D3/
      DATA HR/1.D-1/
      DATA HMINS/1.D-30/
      DATA STF /100.D0/
      DATA RCD/2.D0/
      DATA TOLRI/1.D-5/
      DATA TOLRA/1.D-11/
      DATA TOLABS/0.D0/
C
      IEC = 0
      FV = F
  10  CONTINUE
  20  CONTINUE
C
C TAKE A RANDOM DIRECTION FOR THE DIRECTIONAL DERIVATIVE
      CALL UNITRV(NDIM,W)
C
C COMPUTE FORWARD-DIFFERENCE DERIVATIVE
      CALL DERFOR(NDIM,X,FV,DX,W,DFDX)
C
C TRY THE FIRST HALF-STEP
      DO 30 IC = 1,NDIM
         XP(IC) = X(IC)-H*W(IC)*DFDX*DBLE(FLOAT(NDIM))
  30     CONTINUE
      HS = H
      F = FUNCT0(NDIM,XP)
      FVS = FV+DX*DABS(DFDX)
C
C UPDATE STEPLENGTH  H  AND ACCEPT OR REJECT THE HALF-STEP
      CALL NEWH(KTIM,FVS,F,H,IE,IEC)
      IF(IEC.LE.0)GO TO 50
      IE = IE-1
      IEC = IEC-1
      H = HS
      DFDXV = DFDX
C
C COMPUTE CENTRAL-DIFFERENCES DERIVATIVE
      CALL DERCEN(NDIM,X,FV,DX,W,DFDX)
C
C TRY AGAIN THE FIRST HALF-STEP
      DO 40 IC = 1,NDIM
         XP(IC) = X(IC)-H*W(IC)*DFDX*DBLE(FLOAT(NDIM))
  40     CONTINUE
      F = FUNCT0(NDIM,XP)
      FVS = FV+DX*DABS(DFDXV-DFDX)
C
C
C UPDATE STEPLENGTH  H  AND ACCEPT OR REJECT THE HALF-STEP
      CALL NEWH(KTIM,FVS,F,H,IE,IEC)
C
C UPDATE CUMULATED PAST SCALING DATA
      IF(IEC.GE.1)CALL CUMSCA(NDIM,W,DFDX)
      IF(IEC.GE.1)GO TO 10
      DX = DX*RDX
C
C FIRST HALF-STEP ACCEPTED
  50  CONTINUE
C
C UPDATE CUMULATED PAST SCALING DATA
      CALL CUMSCA(NDIM,W,DFDX)
      FS = FV+DX*DABS(DFDX)
      IF(ITOLCH(FS,FV,TOLRI,TOLABS).EQ.0)DX = DX/RCD
      IF(ITOLCH(FS,FV,TOLRA,TOLABS).GT.0)DX = DX*RCD
      EPSR = DSQRT(H)*EPS
C
C TAKE A SAMPLE INCREMENT OF THE WIENER PROCESS
      CALL GAUSRV(NDIM,W)
C
C TRY THE SECOND HALF-STEP
      DO 60 IC = 1,NDIM
         XP(IC) = XP(IC)+EPSR*W(IC)
  60     CONTINUE
      F = FUNCT0(NDIM,XP)
C
C ACCEPT OR REJECT THE COMPLETE STEP
      IF (F-FV.LE.EPS*EPS*STF) GO TO 70
      H = H*HR
      IE = IE+1
      IF (H.GT.HMINS) GO TO 20
C
C STEP ACCEPTED
  70  CONTINUE
      DO 80 IC = 1,NDIM
         X(IC) = XP(IC)
  80     CONTINUE
      DX = DMIN1(DX,DXMAX)
      DX = DMAX1(DX,DXMIN)
C
      RETURN
      END
      SUBROUTINE NEWH(K,FV,F,H,IE,IEC)
C
C  NEWH  IS CALLED BY THE SUBROUTINE  SSTEP  TO DECIDE WHETHER TO ACCEPT
C  OR REJECT THE FIRST HALF-STEP, AND TO PROVIDE AN UPDATED VALUE FOR  H
C
      DOUBLE PRECISION FV,F,H
      DOUBLE PRECISION FA,FR,HMAX,HMIN,R
      DIMENSION FR(3),FA(4)
      DATA FR/1.05D0,2.D0,10.D0/
      DATA FA/1.D0,1.1D0,2.D0,10.D0/
      DATA HMIN/1.D-30/
      DATA HMAX/1.D25/
      DATA IECMAX/50/
C
      IF(FV.LT.F)GO TO 20
C
C STEP ACCEPTED, POSSIBLY INCREASE THE STEPLENGTH  H
C
      R = FA(1)
      IF(IE*2.LT.K)R = FA(2)
      IF(IE*3.LT.K)R = FA(3)
      IF(IE.EQ.0.AND.K.GT.1)R = FA(4)
      H = H*R
  10  CONTINUE
      IEC = 0
      GO TO 30
C
  20  CONTINUE
      IE = IE+1
      IEC = IEC+1
      IF(IEC.GT.IECMAX)GO TO 10
C
C STEP REJECTED, DECREASE  H
C
      IC = MIN0(3,IEC)
      R = FR(IC)
      H = H/R
  30  CONTINUE
      H = DMIN1(H,HMAX)
      H = DMAX1(H,HMIN)
C
      RETURN
      END
      SUBROUTINE DERFOR(NDIM,X,F,DX,W,DFDX)
C
C  DERFOR  COMPUTED THE FORWARD FINITE-DIFFERNCES DIRECTIONAL
C  DERIVATIVES (CALLING  FUNCT0 )
C
      DOUBLE PRECISION X,F,DX,W,DFDX
      DOUBLE PRECISION DXF,DXFF,DXMAX,FN,S,XX
      DOUBLE PRECISION FUNCT0
      DIMENSION X(NDIM),W(NDIM)
      DIMENSION XX(100)
      DATA DXFF/1.D6/,DXF/1.D1/
      DATA DXMAX/1.D6/
C
  10  CONTINUE
  20  CONTINUE
      S = 0.D0
      DO 30 IC = 1,NDIM
         XX(IC) = X(IC)+W(IC)*DX
         S = S+(XX(IC)-X(IC))**2
  30     CONTINUE
      IF(S.GT.0.D0)GO TO 40
      DX = DX*DXFF
      GO TO 20
  40  CONTINUE
      FN = FUNCT0(NDIM,XX)
      DFDX = (FN-F)/DX
      IF(DX.GT.DXMAX)RETURN
      IF(DABS(DFDX).GT.1.D0)RETURN
      IF(DFDX**2.GT.0.D0)GO TO 50
      DX = DX*DXF
      GO TO 10
  50  CONTINUE
C
      RETURN
      END
      SUBROUTINE DERCEN(NDIM,X,F,DX,W,DFDX)
C
C  DERFOR  COMPUTED THE CENTRAL FINITE-DIFFERNCES DIRECTIONAL
C  DERIVATIVES (CALLING  FUNCT0 )
C
      DOUBLE PRECISION X,F,DX,W,DFDX
      DOUBLE PRECISION FN,FR,XX
      DOUBLE PRECISION FUNCT0
      DIMENSION X(NDIM),W(NDIM)
      DIMENSION XX(100)
C
      DO 10 IC = 1,NDIM
         XX(IC) = X(IC)-W(IC)*DX
  10     CONTINUE
      FR = FUNCT0(NDIM,XX)
      FN = F+DFDX*DX
      DFDX = (FN-FR)/(2.D0*DX)
C
      RETURN
      END
      SUBROUTINE RCLOPT(N,XO,FO)
C
C  RCLOPT  RECALLS THE CURRENT BEST MINIMUM VALUE  FOPT  FOUND SO FAR
C  FROM ALGORITHM START, AND THE POINT  XOPT  (OR POSSIBLY ONE OF THE
C  POINTS) WHERE  FOPT  WAS OBTAINED
C
      DOUBLE PRECISION XO,FO
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XO(N)
C
      DO 10 I = 1,N
         XO(I) = XOPT(I)
  10     CONTINUE
      FO = FOPT
C
      RETURN
      END
      SUBROUTINE STOOPT(N,XO,FO)
C
C  STOOPT  STORES THE CURRENT BEST MINIMUM VALUE  FOPT  FOUND SO FAR
C  FROM ALGORITHM START, AND THE POINT  XOPT  (OR POSSIBLY ONE OF THE
C  POINTS) WHERE  FOPT  WAS OBTAINED
C
      DOUBLE PRECISION XO,FO
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XO(N)
C
      DO 10 I = 1,N
         XOPT(I) = XO(I)
  10     CONTINUE
      FOPT = FO
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION FUNCT0(N,XX)
C
C  FUNCT0  IS CALLED  WHENEVER THE VALUE OF THE FUNCLION  F  IS REQUIRED
C  IN THE NUMERICAL INTEGRATION PROCESS.
C  THE FUNCTION  FUNCT0
C   - RESCALES THE VARIABLES BY CALLING THE SUBROUTINE  VARSCA
C   - CALLS THE SUBROUTINE  RANGE  TO TAKE CARE OF THE CASES WHERE THE
C      CURRENT POINT  X  IS OUTSIDE A GIVEN ADMISSIBLE REGION
C   - CALLS THE USER-SUPPLIED FUNCTION  FUNCT  TO ACTUALLY COMPUTE THE
C      VALUE OF  F
C   - POSSIBLY CALLS THE SUBROUTINE  STOOPT  TO UPDATE THE CURRENT
C        BEST MINIMUM FUNCTION VALUE  FOPT  AND THE CORRESPONDING
C        MINIMIZER  XOPT
C
      DOUBLE PRECISION XX
      DOUBLE PRECISION F,R,XS
      DOUBLE PRECISION FUNCT
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XX(N)
      DIMENSION XS(100)
C
      DO 10 IX = 1,N
         XS(IX) = XX(IX)
  10     CONTINUE
C
C DESCALE  X-VARIABLES
      CALL VARSCA(N,XS)
C
C CONSTRAIN THE  X-VARIABLES WITHIN BOUNDS
      CALL RANGE(N,XS,XRMIN,XRMAX,R)
C
C COMPUTE THE FUNCTION VALUE...
      F = FUNCT(N,XS)+R
C
C ... AND POSSIBLY UPDATE THE BEST CURRENT MINIMUM
      IF(F.LT.FOPT)CALL STOOPT(N,XS,F)
      FUNCT0 = F
      NCF = NCF+1
C
      RETURN
      END
      SUBROUTINE RANGE (N,XS,XRMIN,XRMAX,R)
C
C  RANGE  IS CALLED BY THE FUNCTION  FUNCT0  TO TAKE CARE OF THE CASES
C  WHERE THE CURRENT POINT  X  IS OUTSIDE A GIVEN ADMISSIBLE REGION
C
      DOUBLE PRECISION XS,XRMIN,XRMAX,R
      DOUBLE PRECISION A,B,C,D,DLRMAX,RMAX,RR,XC
      DIMENSION XS(N),XRMIN(N),XRMAX(N)
      DATA RMAX /1.D35/
      DATA DLRMAX/80.5904782547915990D0/
C
      R = 0.D0
      DO 40 I = 1,N
         A = XRMAX(I)
         C = XRMIN(I)
         XC = XS(I)
         IF (XC.LE.A) GO TO 10
         B = A+A-C
         RR = RMAX
         IF (XC.LT.B) RR = DEXP((XC-A)*DLRMAX/(B-A))-1.D0
         R = R+RR
         XS(I) = XRMAX(I)
         GO TO 30
  10     CONTINUE
         IF(XC.GE.C)GO TO 20
         D = C+C-A
         RR = RMAX
         IF (XC.GT.D) RR = DEXP((C-XC)*DLRMAX/(C-D))-1.D0
         R = R+RR
         XS(I) = XRMIN(I)
  20     CONTINUE
  30     CONTINUE
  40     CONTINUE
C
      RETURN
      END
      INTEGER FUNCTION ITOLCH(FMAX,FMIN,TOLREL,TOLABS)
C
C  ITOLCH  DETERMINES  WHETHER TWO QUANTITIES ARE TO BE CONSIDERED NU-
C  MERICALLY EQUAL WITH RESPECT TO AN ABSOLUTE (OR RELATIVE) DIFFERENCE
C  CRITERION, WITHIN GIVEN TOLERANCES
C
      DOUBLE PRECISION FMAX,FMIN,TOLREL,TOLABS
C
      ISTOP = 0
C
C CHECK RELATIVE DIFFERENCE AGAINST  TOLREL
      IF(DABS(FMAX-FMIN).LE.TOLREL*(DABS(FMIN)+DABS(FMAX))/2.D0)
     1    ISTOP=ISTOP+1
C
C CHECK ABSOLUTE DIFFERENCE AGAINST  TOLABS
      IF(FMAX-FMIN.LE.TOLABS)ISTOP = ISTOP+2
      ITOLCH = ISTOP
C
      RETURN
      END
      SUBROUTINE INISCA(N,ND)
C
C  INISCA  INITIALIZES THE COOMON AREA  /SCALE/  FOR THE SCALING DATA
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DATA NXMSCA/10/
C
      LSCA = -1
      IF(N.GT.NXMSCA.OR.N.EQ.1)RETURN
      LSCA = 0
      NX = N
      NORD = ND
      IDSCA = 1
      DO 1 ID = 1,NORD
         DO 2 IX = 1,NX
            DO 3 IY = 1,NX
               DIST(IX,IY,ID) = 0.D0
               GRAGRA(IX,IY,ID) = 0.D0
   3           CONTINUE
            DIST(IX,IX,ID) = 1.D0
            BIAS(IX,ID) = 0.D0
            GRA(IX,ID) = 0.D0
   2        CONTINUE
         NGRA(ID) = 0
   1     CONTINUE
C
      RETURN
      END
      SUBROUTINE NOSCA
C
C  NOSCA  DEACTIVATES THE SCALING
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      LSCA = -1
C
      RETURN
      END
      SUBROUTINE SEGSCA(ID)
C
C  SEGSCA  SELECTS THE TRAJECTORY WHICH MUST BE RESCALED
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      IDSCA = ID
C
      RETURN
      END
      SUBROUTINE VARSCA(N,X)
C
C  VARSCA  COMPUTES THE RESCALED VARIABLE  AX + B
C
      DOUBLE PRECISION X
      DOUBLE PRECISION XB
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION X(N),XB(10)
C
      IF(LSCA.LE.0)RETURN
      DO 1 I = 1,N
         XB(I) = BIAS(I,IDSCA)
         DO  2 J = 1,N
            XB(I) = XB(I)+DIST(I,J,IDSCA)*X(J)
   2        CONTINUE
   1     CONTINUE
      DO 3 I = 1,N
         X(I) = XB(I)
   3     CONTINUE
C
      RETURN
      END
      SUBROUTINE CUMSCA(N,W,DFDX)
C
C  CUMSCA  STORES CUMULATED STATISTICAL DATA ON THE ILL-CONDITIONING OF
C  F(AX+B) W.R.T. X
C
      DOUBLE PRECISION W,DFDX
      DOUBLE PRECISION DFDXMA
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION W(N)
      DATA DFDXMA/1.D8/
C
      IF(LSCA.LE.0)RETURN
      IF(DABS(DFDX).GT.DFDXMA)RETURN
      DO 1 I = 1,N
         DO 2 J = 1,N
            GRAGRA(I,J,IDSCA) = GRAGRA(I,J,IDSCA)+DFDX*W(I)*DFDX*W(J)
   2        CONTINUE
         GRA(I,IDSCA) = GRA(I,IDSCA)+DFDX*W(I)
   1     CONTINUE
      NGRA(IDSCA) = NGRA(IDSCA)+1
C
      RETURN
      END
      SUBROUTINE ACTSCA
C
C  ACTSCA  ACTIVATES THE RESCALING
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      IF(LSCA.EQ.0)LSCA = 1
C
      RETURN
      END
      SUBROUTINE MOVSCA(IU,IM)
C
C  MOVSCA  MOVES THE SCALING DATA OF THE UNPERTURBED CONTINUATION TO THE
C  POSITION OF THE PERTURBED CONTINUATION
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      IF(LSCA.LT.0)RETURN
      DO 1 I = 1,NX
         DO 2 J = 1,NX
            DIST(I,J,IU) = DIST(I,J,IM)
            GRAGRA(I,J,IU) = GRAGRA(I,J,IM)
   2        CONTINUE
         BIAS(I,IU) = BIAS(I,IM)
         GRA(I,IU) = GRA(I,IM)
   1     CONTINUE
      NGRA(IU) = NGRA(IM)
C
      RETURN
      END
      SUBROUTINE UPDSCA(N,X)
C
C  UPDSCA  UPDATES THE SCALING MATRIX  A  AND THE BIAS VECTOR  B BY
C  CALLING  EIGSCA  AND  VARSCA
C
      DOUBLE PRECISION X
      DOUBLE PRECISION AGRAI,ALA1,ALPHA,AMCOR,BIAST
      DOUBLE PRECISION CD,COR,DISTT,SN
      DOUBLE PRECISION EIGSCA
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION X(N)
      DIMENSION DISTT(10,10),BIAST(10),COR(10,10)
      DATA ALPHA /0.3D0/
C
      IF(LSCA.LE.0)RETURN
      ID = IDSCA
      IF(NGRA(ID).LT.2*NX*NX)GO TO 2
      AGRAI = 1.D0/DBLE(FLOAT(NGRA(ID)))
      AMCOR = 0.D0
      DO 3 I = 1,NX
         DO 4 J = 1,NX
         COR(I,J) = AGRAI*GRAGRA(I,J,ID)-AGRAI*GRA(I,ID)*AGRAI*GRA(J,ID)
            AMCOR = DMAX1(AMCOR,DABS(COR(I,J)))
   4        CONTINUE
   3     CONTINUE
      IF(AMCOR.LE.0.D0)GO TO 12
      DO 1 I = 1,NX
         DO 11 J = 1,NX
            COR(I,J) = COR(I,J)/AMCOR
  11        CONTINUE
   1     CONTINUE
      ALA1 = EIGSCA(COR)
      CD = ALA1*(1.D0+ALPHA)
      DO 5 I = 1,NX
         COR(I,I) = COR(I,I)-CD
         BIAST(I) = X(I)
   5     CONTINUE
      CALL VARSCA(NX,BIAST)
      SN = 0.D0
      DO 6 I = 1,NX
         DO 7 J = 1,NX
            DISTT(I,J) = 0.D0
            DO 8 K = 1,NX
               DISTT(I,J) = DISTT(I,J)-DIST(I,K,ID)*COR(K,J)
   8           CONTINUE
            SN = SN+DISTT(I,J)**2
   7        CONTINUE
   6     CONTINUE
      SN = 1.D0/DSQRT(SN/DBLE(FLOAT(NX)))
      DO 9 I = 1,NX
         BIAS(I,ID) = BIAST(I)
         DO 10 J = 1,N
            DIST(I,J,ID) = SN*DISTT(I,J)
            BIAS(I,ID) = BIAS(I,ID)-DIST(I,J,ID)*X(J)
            GRAGRA(I,J,ID) = 0.D0
  10        CONTINUE
         GRA(I,ID) = 0.D0
   9     CONTINUE
      NGRA(ID) = 0
   2  CONTINUE
  12  CONTINUE
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION EIGSCA(COR)
C
C  EIGSCA COMPUTES THE LARGEST EIGENVALUE OF A MATRIX USED FOR RESCA-
C  LING, STARTING FROM A RANDOMLY-CHOSEN ESTIMATE (OBTAINED BY CALLING
C  THE SUBROUTINE  UNITRV ) OF THE CORRESPONDING EIGENVECTOR
C
      DOUBLE PRECISION COR
      DOUBLE PRECISION ALA1,ALA1I,ALA1O
      DOUBLE PRECISION PREC,SWW,W,WW
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION COR(10,10)
      DIMENSION W(10),WW(10)
      DATA PREC /1.D-3/
      DATA NRMIN /10/
      DATA NRMAX /100/
C
      CALL UNITRV(NX,W)
      ALA1 = 0.D0
      DO 1 IR = 1,NRMAX
         ALA1O = ALA1
         SWW = 0.D0
         DO 2 IX = 1,NX
            WW(IX) = 0.D0
            DO 3 JX = 1,NX
               WW(IX) = WW(IX)+COR(IX,JX)*W(JX)
   3           CONTINUE
            SWW = SWW+WW(IX)**2
   2        CONTINUE
         ALA1 = DSQRT(SWW)
         ALA1I = 1.D0/ALA1
         IF(IR.GE.NRMIN.AND.ALA1*PREC.GT.DABS(ALA1-ALA1O))GO TO 4
         DO 5 IX = 1,NX
            W(IX) = WW(IX)*ALA1I
   5        CONTINUE
   1     CONTINUE
   4  CONTINUE
      EIGSCA = ALA1
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION CHAOS(INIZ)
C
C  CHAOS  GENERATES A RANDOM SAMPLE FROM ONE OUT OF FOUR POSSIBLE
C PROBABILITY DISTRIBUTIONS USING RANDOM NUMBERS UNIFORMLY
C DISTRIBUTED IN (0,1) GENERATED BY THE FUNCTION  UNIFRD
C
C  INIZ  IS AN INPUT PARAMETER AS FOLLOWS
C
C      INIZ>0          INITIALIZATION WITH SEED  INIZ.
C      INIZ=0          STANDARD GAUSSIAN DISTRIBUTION.
C      INIZ=-1         CAUCHY DISTRIBUTION.
C      INIZ=-2         UNIFORM DISTRIBUTION IN (-1,+1).
C      OTHERWISE       UNIFORM DISTRIBUTION IN ( 0,+1).
C
      DOUBLE PRECISION UNIFRD,PAI,A,B
C
      DATA PAI/3.14159265358979324D0/
C
      IF(INIZ.LE.0) GO TO 10
C
C INITIALIZATION.
C
      CHAOS = UNIFRD(INIZ)
      RETURN
C
  10  CONTINUE
      A = UNIFRD(0)
      IF(INIZ.NE.0) GO TO 20
      B = UNIFRD(0)
C
C GAUSSIAN RANDOM NUMBER BY POLAR METHOD
C
      CHAOS = DSQRT(-2.D0*DLOG(A))*DCOS(PAI*B)
C
      RETURN
  20  CONTINUE
C
C UNIFORM RANDOM NUMBER IN (0,1)
C
      CHAOS = A
C
C CAUCHY RANDOM NUMBER BY INVERSE TRANSFORMATION
C
      IF(INIZ.EQ.(-1)) CHAOS = DSIN(PAI*A)/DCOS(PAI*A)
C
C UNIFORM RANDOM NUMBER IN (-1,+1)
C
      IF(INIZ.EQ.(-2)) CHAOS = 2.D0*A-1.D0
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION UNIFRD(INIZ)
C
C UNIFRD GENERATES THE RANDOM NUMBERS UNIFORMLY DISTRIBUTED IN (0,1)
C EXPLOITING THOSE GENERATED BY  ALKNUT  WITH A FURTHER RANDOMIZATION
C IF THE INPUT PARAMETER INIZ  IS NOT 0
C THE RANDOM NUMBER GENERATOR IS INITIALIZED
C
      DOUBLE PRECISION A,B,C,X
      DOUBLE PRECISION X0,P,P0,P1,P2,R1,R2
      DOUBLE PRECISION FINV
C
      DIMENSION X(61)
C
      DATA NREM/61/,NRIP/100/
      DATA A,B,C/-1.5D0,5.5D0,-2.0D0/
      DATA FINV/3.5D-5/
      DATA IREM/0/
      DATA P0/3.D0/,P1,P2/1.D0,3.D0/,R1,R2/0.25D0,0.75D0/
C
      IF(INIZ.NE.0.OR.IREM.EQ.0) GO TO 10
C
      I0 = IREM
      X0 = X(I0)
C
C NONLINEARIZATION OF X0 TO AVOID LONG-DISTANCE LINEAR RELATIONSHIP.
C
      IF(X0.GE.FINV)X0 = DMOD(1.D0/X0,1.D0)
C
C UPDATE A COMPONENT OF THE VECTOR X ...
C
      CALL ALKNUT(NREM,X,IREM)
C
C ... AND FURTHER RANDOMIZE
C
      UNIFRD = DMOD(X0+X(I0),1.D0)
C
      RETURN
C
C INITIALIZATION OF THE RANDOM NUMBER GENERATOR
C
  10  CONTINUE
C
      P = P0-1.D0/DBLE(FLOAT(IABS(INIZ))+100.0)
      DO 20 K = 1,NREM
C
         P = C+P*(B+P*A)
C
C
         X(K) = R1+(R2-R1)*(P-P1)/(P2-P1)
  20     CONTINUE
C
      IREM = 0
C
      DO 30 K = 1,NRIP
C
         CALL ALKNUT(NREM,X,IREM)
C
  30     CONTINUE
C
      UNIFRD = X(1)
C
      RETURN
      END
      SUBROUTINE ALKNUT(NREM,X,IREM)
C
C UPDATES THE COMPONENT IREM OF THE NREM-VECTOR X WITH A RANDOM NUMBER
C UNIFORMLY DISTRIBUTED IN (0,1) BY MEANS OF THE ALGORITHM
C OF MITCHELL-MOORE, MODIFIED AS SUGGESTED BY BRENT, QUOTED IN
C D.E.KNUTH, THE ART OF COMPUTER PROGRAMMING, SECOND EDITION,
C SECOND VOLUME, SEMINUMERICAL ALGORITHMS, ADDISON-WESLEY
C PUB. CO., READING (1981), PP. 26-28.
C
      DOUBLE PRECISION X
C
      DIMENSION X(NREM)
C
      DATA N1,N2/24,55/
C
      IF(IREM.NE.0) GO TO 10
C
      IREM = NREM
      I1 = NREM-N1
      I2 = NREM-N2
C
  10  CONTINUE
C
      X(IREM) = DMOD(X(I1)+X(I2),1.D0)
C
      IREM = 1+MOD(IREM,NREM)
      I1 = 1+MOD(I1,NREM)
      I2 = 1+MOD(I2,NREM)
C
      RETURN
      END
      SUBROUTINE GAUSRV(N,W)
C
C GENERATES A RANDOM VECTOR SAMPLE FROM AN N-DIMENSIONAL
C NORMAL DISTRIBUTION
C
      DOUBLE PRECISION W,X,Y,R
      DOUBLE PRECISION CHAOS
C
      DIMENSION W(N)
C
      NN = (N+1)/2
C
      DO 20 I = 1,NN
         II = 1+N-I
C
  10     CONTINUE
         X = CHAOS(-2)
         Y = CHAOS(-2)
         R = X*X+Y*Y
C
         IF(R.GT.1.D0) GO TO 10
C
         R = DSQRT(-2.D0*DLOG(R)/R)
C
         W(I) = X*R
         W(II) = Y*R
C
  20     CONTINUE
C
      RETURN
      END
      SUBROUTINE UNITRV(N,W)
C
C GENERATES A RANDOM VECTOR UNIFORMLY DISTRIBUTED
C ON THE UNIT SPHERE.
C
      DOUBLE PRECISION W,WW
C
      DIMENSION W(N)
C
      CALL GAUSRV(N,W)
      WW = 0.D0
C
      DO 10 I = 1,N
         WW = WW+W(I)**2
  10     CONTINUE
      WW = 1.D0/DSQRT(WW)
C
      DO 20 I = 1,N
         W(I) = WW*W(I)
  20     CONTINUE
C
      RETURN
      END
C
C  MAIN PROGRAM (SAMPLE VERSION)
C  CALLS  SIGMA  VIA THE DRIVER SUBROUTINE  SIGMA1
C
      DOUBLE PRECISION X0,XMIN,FMIN
C
      DIMENSION X0(2),XMIN(2)
C
C TEST PROBLEM DATA
C
C PROBLEM DIMENSION
      N = 2
C
C INITIAL POINT
      X0(1) = 0.D0
      X0(2) = 0.D0
C
C SET INPUT PARAMETERS
      NSUC = 3
      IPRINT = 0
C
C CALL DRIVER SUBROUTINE  SIGMA1
      CALL SIGMA1(N,X0,NSUC,IPRINT,XMIN,FMIN,NFEV,IOUT)
C
      STOP
      END
      DOUBLE PRECISION FUNCTION FUNCT (N,X)
C
C COMPUTES THE VALUE AT  X  OF THE SIX-HUMP CAMEL FUNCTION
C
      DOUBLE PRECISION X,XX,YY
      DIMENSION X(N)
      XX = X(1)*X(1)
      YY = X(2)*X(2)
      FUNCT = ((XX/3.D0-2.1D0)*XX+4.D0)*XX+X(1)*X(2)
     1    +4.D0*(YY-1.D0)*YY
      RETURN
      END
      SUBROUTINE SIGMA1(N,X0,NSUC,IPRINT,XMIN,FMIN,NFEV,IOUT)
C
C  SIGMA1  IS A 'DRIVER' SUBROUTINE WHICH SIMPLY CALLS THE PRINCIPAL
C  SUBROUTINE  SIGMA  AFTER HAVING ASSIGNED DEFAULT VALUES TO A NUMBER
C  OF INPUT PARAMETERS OF  SIGMA , AND HAS THEREFORE A CONSIDERABLY
C  LOWER NUMBER OF INPUT PARAMETERS.
C  IT CAN BE USED AS A SIMPLE EXAMPLE OF HOW TO CALL  SIGMA , BUT
C  ALSO AS AN EASY-TO-USE DRIVER FOR THE AVERAGE USER, WHICH MAY FIND
C  IT EASIER TO CALL  SIGMA1  INSTEAD OF  SIGMA , THUS AVOIDING THE
C  TROUBLE OF ASSIGNING A VALUE TO ALL THE INPUT PARAMETERS OS  SIGMA .
C  ALL THE PARAMETERS IN THE DEFINITION OF  SIGMA1  HAVE THE SAME MEA-
C  NING AS IN  SIGMA .
C
C  THE USER OF  SIGMA1  MUST ONLY GIVE VALUES TO THE INPUT PARAMETERS
C      N, X0, NSUC, IPRINT
C  AND OBTAINS ON OUTPUT THE SAME OUTPUT PARAMETERS OF  SIGMA
C      XMIN, FMIN, NFEV, IOUT
C
C  WE RECALL HERE THE MEANING OF THE ABOVE PARAMETERS
C
C  N       IS THE PROBLEM DIMENSION (NUMBER OF VARIABLES)
C  X0      IS AN  N-VECTOR CONTAINING THE INITIAL VALUES OF THE
C          X-VARIABLES
C  NSUC    IS THE REQUESTED NUMBER OF SUCCESSFUL TRIALS
C          AFTER WHICH THE COMPUTATION IS STOPPED.
C  IPRINT  IS AN INDEX USED TO CONTROL THE AMOUNT OF PRINTED OUTPUT
C          BY CONTROLLING THE CALLS TO THE USER-SUPPLIED OUTPUT SUB-
C          ROUTINES  PTSEG  (END-OF-SEGMENT OUTPUT),  PTRIAL  (END-
C          OF-TRIAL OUTPUT), AND  PTKSUC (END-OF-TRIAL OUTPUT RELATED
C          TO THE COUNT OF SUCCESSFUL TRIALS), WHICH ARE DESCRIBED
C          BELOW.
C          IPRINT.LT.0   NO CALL TO THE OUTPUT SUBROUTINES
C          IPRINT.EQ.0   CALL ONLY  PTRIAL  AND  PTKSUC
C          IPRINT.GT.0   CALL ALL OUTPUT SUBROUTINES.
C  XMIN    IS AN N-VECTOR CONTAINING THE COORDINATES OF THE POINT
C          (OR POSSIBLY ONE OF THE POINTS) WHERE THE FINAL VALUE  FMIN
C          OF  FOPT  WAS FOUND.
C  FMIN    IS THE FINAL VALUE OF THE BEST CURRENT MINIMUM FUNCTION
C          VALUE  FOPT.
C  NFEV    IS THE TOTAL NUMBER OF FUNCTION EVALUATION (INCLUDING
C          THOSE USED FOR THE COMPUTATION OF DERIVATIVES, AND FOR
C          THE REJECTED TIME-INTEGRATION STEPS).
C  IOUT    IS THE INDICATOR OF THE STOPPING CONDITIONS, AS FOLLOWS
C          IF  IOUT = -99  A FATAL ERROR WAS DETECTED WHEN PERFOR-
C             MING SOME PRELIMINARY CHECKING OF THE INPUT DATA, AND
C             THE ALGORITHM WAS NOT EVEN STARTED
C          OTHERWISE THE ALGORITHM WAS STARTED, AND THE MEANING OF
C          IOUT  IS AS FOLLOWS
C          IF  IOUT.EQ.0  NO TRIAL HAD A UNIFORM STOP.
C          IF  IOUT.NE.0  THE SIGN AND THE MAGNITUDE OF  IOUT  HAVE
C          THE FOLLOWING MEANING
C          A) THE SIGN OF  IOUT  INDICATES IF THE BEST UNIFORM STOP
C          (I.E. THE ONE IN WHICH  FTFOPT  WAS OBTAINED)
C          IS TO BE CONSIDERED SUCCESSFUL (IOUT.GT.0) OR UNSUCCESSFUL
C          B) THE MAGNITUDE OF  IOUT  INDICATES WHICH ONE OF THE
C          NUMERICAL EQUALITY CRITERIA WAS SATISFIED IN THE BEST
C          UNIFORM STOP
C          IABS(IOUT).EQ.1  RELATIVE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.2  ABSOLUTE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.3  BOTH CRITERIA
C          (IOUT  IS THE FINAL VALUE OF THE INTERNAL PARAMETER  ISTOPT
C          (AN OUTPUT INDICATOR OF THE USER-SUPPLIED SUBROUTINE
C          PTRIAL )).
C          SUCCESS IS CLAIMED BY THE ALGORITHM IF  IOUT .GT. 0,
C          I.E. IF AT LEAST ONE OF THE TRIALS STOPPED SUCCESSFULLY
C          (I.E. WITH A POSITIVE VALUE OF THE TRIAL STOP INDICATOR
C          ISTOP ), AND THE SUCCESSFUL TRIALS (AS COUNTED BY  KSUC )
C          ARE CURRENTLY VALID.
C
      DOUBLE PRECISION X0,XMIN,FMIN
      DOUBLE PRECISION DX,EPS,H,TOLABS,TOLREL
      DOUBLE PRECISION VRMAX,VRMIN,XRMAX,XRMIN
      DIMENSION X0(N),XMIN(N)
      DIMENSION XRMIN(100),XRMAX(100)
      DATA VRMIN,VRMAX /-1.D4,1.D4/
      DATA NTRILD/50/
C
      H = 1.D-10
      EPS = 1.D0
      DX = 1.D-9
      IRAND = 0
      NTRAJ = 0
      ISEGBR = 0
      INKPBR = 0
      KPBR0 = 0
      NPMIN = 10
      NPMAX0 = 100
      INPMAX = 50
      NTRIAL = MAX0(NTRILD,5*NSUC)
      TOLREL = 1.D-3
      TOLABS = 1.D-6
      KPASCA = 10
      IF(N.GT.5)KPASCA = 300
      INHP = 1
      DO 1 IX = 1,N
         XRMIN(IX)=VRMIN
         XRMAX(IX)=VRMAX
 1    CONTINUE
C
      CALL SIGMA ( N, X0, H, EPS, DX,
     1             NTRAJ, ISEGBR, KPBR0, INKPBR,
     2             NPMIN, NPMAX0, INPMAX,
     3             NSUC, NTRIAL, TOLREL, TOLABS, XRMIN, XRMAX,
     4             KPASCA, IRAND, INHP, IPRINT,
     5             XMIN, FMIN, NFEV, IOUT )
C
      RETURN
      END
      SUBROUTINE SIGMA ( N, X0, H, EPS, DX,
     1                   NTRAJ, ISEGBR, KPBR0, INKPBR,
     2                   NPMIN, NPMAX0, INPMAX,
     3                   NSUC, NTRIAL, TOLREL, TOLABS, XRMIN, XRMAX,
     4                   KPASCA, IRAND, INHP, IPRINT,
     5                   XMIN, FMIN, NFEV, IOUT )
C
C  THE SUBROUTINE  SIGMA  IS THE PRINCIPAL SUBROUTINE OF THE PACKAGE
C  SIGMA, WHICH ATTEMPTS TO FIND A GLOBAL MINIMIZER OF A REAL VALUED
C  FUNCTION  F(X) = F(X1,...,XN)  OF  N  REAL VARIABLES  X1,...,XN.
C  THE ALGORITHM AND THE PACKAGE ARE DESCRIBED IN DETAIL IN THE TWO
C  PAPERS PUBLISHED IN THE SAME ISSUE OF THE A.C.M. TRANSACTIONS ON
C  MATHEMATICAL SOFTWARE, BOTH BY
C     F. ALUFFI-PENTINI, V. PARISI, F. ZIRILLI.
C  (1)  A GLOBAL MINIMIZATION ALGORITHM USING STOCHASTIC DIFFERENTIAL
C       EQUATIONS
C  (2)  ALGORITHM SIGMA, A STOCHASTIC-INTEGRATION GLOBAL MINIMIZATION
C       ALGORITHM.
C  THE SOFTWARE IMPLEMENTATION AND ITS USAGE ARE DESCRIBED IN (2).
C
C  METHOD
C
C  A GLOBAL MINIMIZER OF  F(X)  IS SOUGHT BY MONITORING THE VALUES OF
C  F  ALONG TRAJECTORIES GENERATED BY A SUITABLE (STOCHASTIC) DISCRE-
C  TIZATION OF A FIRST-ORDER STOCHASTIC DIFFERENTIAL EQUATION INSPIRED
C  BY STATISTICAL MECHANICS. STARTING FROM AN INITIAL POINT  X0 ,
C  X  IS UPDATED BY THE (STOCHASTIC) DISCRETIZATION STEP
C       X = X + DX1 + DX2
C  WHERE  DX1 = - H * GAM           (FIRST HALF-STEP)
C         DX2 = EPS * SQRT(H) * U   (SECOND HALF-STEP)
C  AND  H  IS THE TIME-INTEGRATION STEPLENGTH,
C  GAM/N  IS COMPUTED AS A FINITE-DIFFERENCE APPROXIMATION TO THE
C     DIRECTIONAL DERIVATIVE OF  F  ALONG AN ISOTROPICALLY RANDOM
C     DIRECTION,
C  EPS  IS A POSITIVE 'NOISE' COEFFICIENT, AND
C  U  IS A RANDOM SAMPLE FROM AN N-DIMENSIONAL GAUSSIAN DISTRIBUTION.
C  WE CONSIDER THE SIMULTANEOUS EVOLUTION OF A GIVEN FIXED NUMBER
C  NTRAJ  OF TRAJECTORIES DURING AN OBSERVATION PERIOD IN WHICH FOR
C  EACH TRAJECTORY  EPS  IS FIXED WHILE  H  AND THE SPATIAL DISCRETI-
C  ZATION INCREMENT  DX  FOR COMPUTING  GAM  ARE AUTOMATICALLY
C  ADJUSTED BY THE ALGORITHM.
C  AFTER EVERY OBSERVATION PERIOD ONE OF THE TRAJECTORIES IS DISCARDED,
C  ALL OTHER TRAJECTORIES CONTINUE UNPERTURBED, AND ONE OF THEM IS SE-
C  LECTED FOR BRANCHING, I.E. GENERATING A SECOND PERTURBED CONTI-
C  NUATION, WITH DIFFERENT STARTING  EPS  AND  DX  (AND THE SAME
C  'PAST HISTORY' OF THE FIRST).
C  THE SET OF SIMULTANEOUS TRAJECTORIES IS CONSIDERED A SINGLE TRIAL,
C  AND THE COMPLETE ALGORITHM IS A SET OF REPEATED TRIALS.
C  A TRIAL IS STOPPED, AT THE END OF AN OBSERVATION PERIOD, AND AFTER
C  HAVING DISCARDED THE WORST TRAJECTORY, IF ALL THE FINAL VALUES OF
C  F  FOR THE REMAINING TRAJECTORIES ARE EQUAL (WITHIN NUMERICAL TO-
C  LERANCES, AND POSSIBLY AT DIFFERENT POINTS  X) TO THEIR MINIMUM
C  VALUE  FTFMIN ('UNIFORM STOP AT THE LEVEL  FTFMIN ').
C  A UNIFORM STOP IS CONSIDERED SUCCESSFUL ONLY IF THE FINAL VALUE
C  FTFMIN  IS (NUMERICALLY) EQUAL TO THE CURRENT BEST MINIMUM  FOPT
C  FOUND SO FAR FROM ALGORITHM START.
C  A TRIAL IS ALSO ANYWAY STOPPED (UNSUCCESSFULLY) IF A GIVEN MAXIMUM
C  NUMBER  NPMAX  OF OBSERVATION PERIODS HAS ELAPSED.
C  TRIALS ARE REPEATED WITH DIFFERENT OPERATING CONDITIONS (INITIAL
C  POINT, MAX TRIAL LENGTH  NPMAX , SEED OF NOISE GENERATOR, POLICY
C  FOR CHOOSING THE STARTING  EPS  FOR THE PERTURBED CONTINUATION,
C  AND TRIAL-START VALUE OF  EPS ).
C  THE OUTCOMES OF ALL THE COMPLETED TRIALS ARE SUMMARIZED BY
C  FTFOPT (THE BEST CURRENT MIMNIMUM VALUE FOUND SO FAR FOR
C  FTFMIN ) AND BY THE CURRENT COUNT  KSUC  OF SUCCESSFUL TRIALS
C  (WHICH IS ZERO AT ALGORITHM START, AND IS UPDATED AT EVERY
C  SUCCESSFULLY-STOPPING TRIALS).
C  THE ALGORITHM IS STOPPED, AT THE END OF A TRIAL, IF A REQUESTED
C  COUNT  NSUC  OF SUCCESSFUL TRIALS HAS BEEN REACHED,
C  OR ANYWAY IF A GIVEN MAXIMUM NUMBER  NTRIAL  OF TRIALS HAS BEEN
C  REACHED.
C  SUCCESS IS CLAIMED IF THE CURRENT COUNT  KSUC  OF SUCCESSFUL TRIALS
C  IS AT LEAST ONE, AND IF SUCH TRIALS ARE CURRENTLY VALID.
C
C  CALL STATEMENT
C
C  THE CALL STATEMENT IS
C     CALL SIGMA ( N, X0, H, EPS, DX,
C                  NTRAJ, ISEGBR, KPBR0, INKPBR,
C                  NPMIN, NPMAX0, INPMAX,
C                  NSUC, NTRIAL, TOLREL, TOLABS, XRMIN, XRMAX,
C                  KPASCA, IRAND, INHP, IPRINT,
C                  XMIN, FMIN, NFEV, IOUT )
C
C  CALL PARAMETERS
C
C  INPUT PARAMETERS ARE THOSE IN LINES 1,3,4,5 OF THE CALL STATEMENT,
C  INPUT-OUTPUT PARAMETERS ARE THOSE IN LINE 2,
C  OUTPUT PARAMETERS ARE THOSE IN LINE 6.
C  NOTE THAT A NUMBER OF OTHER (INTERNAL) PARAMETERS CAN BE OBTAINED
C  BY MEANS OF THE USER-SUPPLIED OUTPUT SUBROUTINES  PTSEG, PTRIAL,
C  AND  PTKSUC, WHICH ARE DESCRIBED BELOW.
C
C  DESCRIPTION OF THE CALL PARAMETERS
C
C  N        IS THE PROBLEM DIMENSION (NUMBER OF VARIABLES)
C  X0       IS AN  N-VECTOR CONTAINING THE INITIAL VALUES OF THE
C          X-VARIABLES
C  H       IS THE INITIAL VALUE OF THE TIME-INTEGRATION STEPLENGTH.
C  EPS     IS THE INITIAL VALUE OF THE NOISE COEFFICIENT.
C  DX      IS THE INITIAL VALUE OF THE MAGNITUDE OF THE DISCRETIZATION
C          INCREMENT FOR COMPUTING THE FINITE-DIFFERENCE DERIVATIVES.
C  NTRAJ   IS THE NUMBER OF SIMULTANEOUS TRAJECTORIES.
C          (NOTE HOWEVER THAT IF THE INPUT VALUE IS ZERO,  NTRAJ  IS
C          SET TO A DEFAULT VALUE (NTRAJ = 7), AND IF THE INPUT VALUE
C          IS OTHERWISE OUTSIDE THE INTERVAL (3,20) NTRAJ IS SET TO
C          THE NEAREST EXTREME VALUE).
C  ISEGBR, KPBR0, INKPBR  DETERMINE, AT THE END OF AN OBSERVATION
C          PERIOD, WHICH ONE OF THE SIMULTANEOUS TRAJECTORIES
C          IS TO BE BRANCHED, AS FOLLOWS.
C          BRANCHING IS NORMALLY PERFORMED ON THE TRAJECTORY WHICH
C          OCCUPIES THE PLACE  ISEGBR  IN THE TRAJECTORY SELECTION OR-
C          DERING, EXCEPT AT (THE END OF) EXCEPTIONAL OBSERVATION
C          PERIODS, WHERE THE FIRST TRAJECTORY IN THE ORDERING IS
C          BRANCHED. EXCEPTIONAL BRANCHING OCCURS AT THE OBSERVATION
C          PERIODS NUMBERED  KP = KPBR0 + J*INKPBR, (J = 1,2,3,...).
C          THEREFORE  ISEGBR  SELECTS THE LEVEL (IN THE ORDERING) AT
C          WHICH NORMAL BRANCHING OCCURS, WHILE  KPBR0  AND  INKPBR
C          SELECT THE FIRST OCCURRENCE AND THE REPETITION FREQUENCY
C          OF THE EXCEPTIONAL OBSERVATION PERIODS.
C          (NOTE HOWEVER THAT IF ONE OF THE INPUT VALUES IS ZERO,
C          THE CORRESPONDING VARIABLE IS SET TO A DEFAULT VALUE
C          ISEGBR = INT((1+NTRAJ)/2), INKPBR = 10, KPBR0 = 3.
C          IF THE INPUT VALUE FOR  ISEGBR  IS OTHERWISE OUTSIDE THE
C          INTERVAL  (1,NTRAJ),  ISEGBR  IS SET TO THE NEAREST
C          EXTREME VALUE, AND IF KPBR0  HAS A VALUE NOT INSIDE THE
C          INTERVAL  (1,INKPBR), IT IS ASSIGNED THE SAME VALUE
C          MODULO  INKPBR).
C  NPMIN   IS THE MINIMUM DURATION OF A TRIAL, I.E. THE MINIMUM
C          NUMBER OF OBSERVATION PERIODS THAT SHOULD ELAPSE BEFORE
C          STARTING TO CHECK THE TRIAL STOPPING CRITERIA.
C  NPMAX0  IS THE MAXIMUM DURATION OF THE FIRST TRIAL, I.E. THE
C          VALUE, FOR THE FIRST TRIAL, OF MAXIMUM ACCEPTABLE
C          NUMBER  NPMAX  OF OBSERVATION PERIODS IN A TRIAL.
C  INPMAX  IS THE INCREMENT FOR  NPMAX , WHEN  NPMAX  IS VARIED
C          FROM ONE TRIAL TO THE FOLLOWING ONE.
C  NSUC    IS THE REQUESTED NUMBER OF SUCCESSFUL TRIALS
C          AFTER WHICH THE COMPUTATION IS STOPPED.
C  TOLREL  AND  TOLABS  ARE THE RELATIVE AND ABSOLUTE TOLERANCES
C          FOR STOPPING A SINGLE TRIAL.
C  XRMIN,  XRMAX  ARE N-VECTORS DEFINING THE ADMISIBLE REGION FOR
C          THE X-VALUES, WITHIN WHICH THE FUNCTION VALUES CAN BE
C          SAFELY COMPUTED.
C  KPASCA  IS THE MINIMUM NUMBER OF TRAJECTORY SEGMENTS THAT SHOULD
C          ELAPSE BEFORE THE RESCALING PROCEDURES ARE ACTIVATED.
C  IRAND   IS A CONTROL INDEX FOR THE INITIALIZATION OF THE RANDOM
C          NUMBER GENERATOR.
C          IRAND.GT.0    THE GENERATOR IS INITIALIZED, BEFORE STAR-
C                        TING THE TRIAL  KT, WITH SEED  IRAND+KT-1
C          IRAND.LE.0    THE GENERATOR IS INITIALIZED (WITH SEED 0)
C                        ONLY AT THE FIRST CALL OF SIGMA
C  INHP    IS A CONTROL INDEX FOR SELECTING THE NUMBER  NHP  OF TIME-
C          INTEGRATION STEPS FOR OBSERVATION PERIOD  KP (DURATION OF
C          TRIAL  KP) AS FOLLOWS (LOG IS BASE 2)
C          INHP=1  NHP = 1 + INT(LOG(KP))      ('SHORT' DURATION)
C          INHP=2  NHP = INT(SQRT(KP))         ('MEDIUM' DURATION)
C          INHP=3  NHP = KP                    ('LONG' DURATION)
C  IPRINT  IS AN INDEX USED TO CONTROL THE AMOUNT OF PRINTED OUTPUT
C          BY CONTROLLING THE CALLS TO THE USER-SUPPLIED OUTPUT SUB-
C          ROUTINES  PTSEG  (END-OF-SEGMENT OUTPUT),  PTRIAL  (END-
C          OF-TRIAL OUTPUT), AND  PTKSUC (END-OF-TRIAL OUTPUT RELATED
C          TO THE COUNT OF SUCCESSFUL TRIALS), WHICH ARE DESCRIBED
C          BELOW.
C          IPRINT.LT.0   NO CALL TO THE OUTPUT SUBROUTINES
C          IPRINT.EQ.0   CALL ONLY  PTRIAL  AND  PTKSUC
C          IPRINT.GT.0   CALL ALL OUTPUT SUBROUTINES.
C  XMIN    IS AN N-VECTOR CONTAINING THE COORDINATES OF THE POINT
C          (OR POSSIBLY ONE OF THE POINTS) WHERE THE FINAL VALUE  FMIN
C          OF  FOPT  WAS FOUND.
C  FMIN    IS THE FINAL VALUE OF THE BEST CURRENT MINIMUM FUNCTION
C          VALUE  FOPT.
C  NFEV    IS THE TOTAL NUMBER OF FUNCTION EVALUATION (INCLUDING
C          THOSE USED FOR THE COMPUTATION OF DERIVATIVES, AND FOR
C          THE REJECTED TIME-INTEGRATION STEPS).
C  IOUT    IS THE INDICATOR OF THE STOPPING CONDITIONS, AS FOLLOWS
C          IF  IOUT = -99  A FATAL ERROR WAS DETECTED WHEN PERFOR-
C             MING SOME PRELIMINARY CHECKING OF THE INPUT DATA, AND
C             THE ALGORITHM WAS NOT EVEN STARTED
C          OTHERWISE THE ALGORITHM WAS STARTED, AND THE MEANING OF
C          IOUT  IS AS FOLLOWS
C          IF  IOUT.EQ.0  NO TRIAL HAD A UNIFORM STOP.
C          IF  IOUT.NE.0  THE SIGN AND THE MAGNITUDE OF  IOUT  HAVE
C          THE FOLLOWING MEANING
C          A) THE SIGN OF  IOUT  INDICATES IF THE BEST UNIFORM STOP
C          (I.E. THE ONE IN WHICH  FTFOPT  WAS OBTAINED)
C          IS TO BE CONSIDERED SUCCESSFUL (IOUT.GT.0) OR UNSUCCESSFUL
C          B) THE MAGNITUDE OF  IOUT  INDICATES WHICH ONE OF THE
C          NUMERICAL EQUALITY CRITERIA WAS SATISFIED IN THE BEST
C          UNIFORM STOP
C          IABS(IOUT).EQ.1  RELATIVE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.2  ABSOLUTE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.3  BOTH CRITERIA
C          (IOUT  IS THE FINAL VALUE OF THE INTERNAL PARAMETER  ISTOPT
C          (AN OUTPUT INDICATOR OF THE USER-SUPPLIED SUBROUTINE
C          PTRIAL )).
C          SUCCESS IS CLAIMED BY THE ALGORITHM IF  IOUT .GT. 0,
C          I.E. IF AT LEAST ONE OF THE TRIALS STOPPED SUCCESSFULLY
C          (I.E. WITH A POSITIVE VALUE OF THE TRIAL STOP INDICATOR
C          ISTOP ), AND THE SUCCESSFUL TRIALS (AS COUNTED BY  KSUC )
C          ARE CURRENTLY VALID.
C
C
C  USER-SUPPLIED SUPROGRAMS
C
C  THE USER MUST PROVIDE THE FUNCTION  FUNCT  TO COMPUTE  F(X),
C  AND THE THREE OUTPUT SUBROUTINE  PTSEG, PTRIAL, PTKSUC .
C  THE CALLS TO THE OUTPUT SUBROUTINES ARE CONTROLLED BY  IPRINT
C  (INPUT PARAMETER TO  SIGMA).
C  A USER NOT INTERESTED IN USING ANY ONE OF THE OUTPUT SUBROUTINES
C  CAN SIMPLY PUT  IPRINT = -1.
C  SAMPLE VERSIONS OF THE OUTPUT SUBROUTINES ARE PROVIDED
C  WITH THE PACKAGE FOR THE BENEFIT OF THE AVERAGE USER.
C  IN THE FOLLOWING DESCRIPTION ALL NON-INTEGER ARGUMENTS ARE
C  DOUBLE PRECISION (INTEGER ARGUMENTS ARE INDICATED BY MEANS OF THE
C  FORTRAN IMPLICIT TYPE DEFINITION CONVENTION).
C
C     THE FUNCTION  FUNCT
C
C     FUNCT  MUST RETURN AS ITS VALUE THE VALUE AT  X  OF THE FUNCTION
C     TO BE MINIMIZED
C     THE DEFINITION STATEMENT IS
C        DOUBLE PRECISION FUNCTION  FUNCT (N, X)
C     WHERE
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM.
C     X   IS THE (INPUT) N-VECTOR CONTAINING THE COORDINATES OF THE
C         POINT  X  WHERE THE FUNCTION IS TO BE COMPUTED.
C
C     THE SUBROUTINE  PTSEG
C
C     PTSEG  IS CALLED (IF  IPRINT.GT.0) AT THE END OF EVERY OBSER-
C     VATION PERIOD.
C     THE DEFINITION STATEMENT IS
C         SUBROUTINE  PTSEG ( N, XPFMIN, FPFMIN, FPFMAX,
C                                         KP, NFEV, IPRINT )
C     WHERE
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM
C     FPFMIN, FPFMAX  ARE RESPECTIVELY THE MINIMUM AND THE MAXIMUM
C        AMONG THE VALUES OF  F(X) OBTAINED AT THE FINAL POINTS OF
C        THE TRAJECTORY SEGMENTS OF THE (JUST ELAPSED) OBSERVATION
C        PERIOD  KP.
C     XPFMIN   IS AN N-VECTOR CONTAINING THE COORDINATES OF THE
C        (FINAL) POINT (OR POSSIBLY ONE OF THE POINTS) WHERE
C        FPFMIN WAS OBTAINED.
C     KP   IS THE TOTAL NUMBER OF ELAPSED OBSERVATION PERIODS IN
C        THE CURRENT TRIAL.
C     NFEV   IS THE TOTAL NUMBER OF FUNCTION EVALUATIONS PERFORMED
C        FROM ALGORITHM START.
C
C     THE SUBROUTINE  PTRIAL
C
C     PTRIAL  IS CALLED (IF IPRINT.GE.0) AT THE END OF EVERY TRIAL.
C     THE DEFINITION STATEMENT IS
C         SUBROUTINE PTRIAL ( N, XOPT, FOPT,
C                             FTFMIN, FTFMAX, FTFOPT,
C                             ISTOP, ISTOPT, NFEV, KP, IPRINT )
C     WHERE
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM.
C     XOPT   IS AN N-VECTOR CONTINING THE COORDINATES OF THE
C        POINT (OR POSSIBLY ONE OF THE POINTS) WHERE THE CURRENT
C        MINIMUM  FOPT  WAS OBTAINED.
C     FOPT   IS THE CURRENT BEST MINIMUM VALUE FOUND FOR  F  FROM
C        ALGORITHM START ( FOPT  IS UPDATED WHENEVER A FUNCTION
C        VALUE IS COMPUTED).
C     FTFMIN, FTFMAX   ARE RESPECTIVELY THE MINIMUM AND THE MAXIMUM
C        AMONG THE VALUES OF  F(X)  OBTAINED AT THE FINAL POINTS OF
C        THE LAST TRAJECTORY SEGMENTS OF THE CURRENT TRIAL.
C     FTFOPT   IS THE CURRENT MINIMUM VALUE OF  FTFMIN  AMONG THE
C        TRIALS WHICH DID NOT STOP FOR REACHING THE MAXIMUM ALLOWED
C        NUMBER OF SEGMENTS (STOPPING INDICATOR  ISTOP = 0, SEE
C        BELOW).  FTFOPT  IS USED BY  SIGMA  TO COMPUTE  KSUC  (INPUT
C        PARAMETER TO THE OUTPUT SUBROUTINE  PTKSUK , SEE BELOW).
C     KP   IS THE TOTAL NUMBER OF ELAPSED OBSERVATION PERIODS IN
C        THE CURRENT TRIAL.
C     NFEV   IS THE TOTAL NUMBER OF FUNCTION EVALUATIONS PERFORMED
C        FROM ALGORITHM START.
C     ISTOP   IS THE INDICATOR OF THE STOPPING CONDITION OF THE TRIAL,
C        AS FOLLOWS
C        ISTOP = 0
C           THE MAXIMUM NUMBER  NPMAX  OF OBSERVATION PERIODS HAS
C           BEEN REACHED.
C        ISTOP.NE.0
C           ALL THE END-OF-SEGMENT VALUES OF  F(X) , (EXCEPT FOR THE
C           JUST DISCARDED SEGMENT) ARE CLOSE ENOUGH TO THEIR COMMON
C           MINIMUM VALUE  FPFMIN , WITH RESPECT TO AN ABSOLUTE OR
C           RELATIVE DIFFERENCE CRITERION, TO BE CONSIDERED NUMERI-
C           CALLY EQUAL.
C           THE ABSOLUTE VALUE AND THE SIGN OF  ISTOP  HAVE THE
C           FOLLOWING MEANING.
C              THE ABSOLUTE VALUE INDICATES WHICH DIFFERENCE
C                 CRITERION WAS SATISFIED
C                 1   RELATIVE DIFFERENCE CRITERION SATISFIED
C                 2   ABSOLUTE DIFFERENCE CRITERION SATISFIED
C                 3   BOTH CRITERIA SATISFIED
C              THE SIGN OF  ISTOP  INDICATES THE RELATIONSHIP BETWEEN
C                 THE END-OF-TRIAL VALUE  FPFMIN  AND THE CURRENT
C                 BEST MINIMUM VALUE  FOPT  (WHICH IS UPDATED WHEN-
C                 EVER A FUNCTION  VALUE IS COMPUTED
C                 ISTOP.GT.0
C                    FPFMIN  IS NUMERICALLY EQUAL (W.R.T. AT LEAST
C                    ONE OF THE ABOVE  DIFFERENCE CRITERIA) TO  FOPT
C                 ISTOP.LT.0
C                    FPFMIN  IS NOT EVEN NUMERICALLY EQUAL TO  FOPT
C                    (AND THEREFORE CANNOT BE CONSIDERED AS AN
C                    ACCEPTABLE GLOBAL MINIMUM).
C     ISTOPT   IS THE VALUE OF THE TRIAL STOPPING INDICATOR  ISTOP
C        CORRESPONDING TO THE (CURRENT OR PAST) TRIAL WHERE  FTFOPT
C        WAS OBTAINED, WITH THE SIGN WHICH IS UPDATED ACCORDING TO
C        THE COMPARISON BETWEEN  FTFOPT  AND THE PRESENT VALUE OF
C        FOPT , AS DESCRIBED ABOVE.
C        THE FINAL VALUE OF  ISTOP IS RETURNED BY  SIGMA  AS THE VALUE
C        OF THE OUTPUT INDICATOR IOUT OF THE ALGORITHM STOPPING CON-
C        DITIONS (WHENEVER THE ALGORITHM WAS STARTED, IOUT.NE.-99,
C        SEE ABOVE).
C
C     THE SUBROUTINE  PTKSUC
C
C     PTKSUC  IS CALLED ONLY AT THE END OF EVERY SUCCESSFUL TRIAL
C     SUCH THAT AN INCREMENT OCCURRED IN THE VALUE OF THE
C     CURRENT MAXIMUM  MSUC  AMONG ALL THE VALUES OF  KSUC  FROM
C     ALGORITHM START. A CALL TO  PTKSUC  THEREFORE PROVIDES
C     THE USER WITH THE OPERATIONALLY INTERESTING INFORMATION THAT
C     ALGORITHM STOP WOULD HAVE TAKEN PLACE, IF  NSUC  (INPUT
C     PARAMETER TO  SIGMA ) HAD BEEN GIVEN A LOWER VALUE, EQUAL TO
C     THE CURRENT  KSUC .
C     PTKSUC  IS CALLED ONLY IF  IPRINT.GE.0  AND  KSUC.LT.NSUC  .
C     THE DEFINITION STATEMENT IS
C         SUBROUTINE  PTKSUC ( KSUC )
C     WHERE  KSUC  IS THE CURRENT COUNT OF SUCCESSFUL TRIALS
C     (1 .LE. KSUC .LE. NSUC)
C
      DOUBLE PRECISION X0,H,EPS,DX,TOLREL,TOLABS
      DOUBLE PRECISION XRMIN,XRMAX,XMIN,FMIN
      DOUBLE PRECISION EPSAG,EPSAP,EPSC
      DOUBLE PRECISION EPSMAX,EPSR,F,FOPT,FTFMAX
      DOUBLE PRECISION FTFMIN,FTFOPT
C
      DOUBLE PRECISION X,HC,DXC,VMVT,EPSCO,VMCOR,VCOR
      DOUBLE PRECISION XRMIC,XRMAC,XOPT,FOPTC
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
C
      DIMENSION X0(N),XMIN(N),XRMIN(N),XRMAX(N)
C
      COMMON /DINCOM/ X(100,20),HC(20),DXC(20),VMVT(20,19),EPSCO(20),
     1    VMCOR(20),VCOR(20),XRMIC(100),XRMAC(100),XOPT(100),FOPTC,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJC,NTRAJR,
     3    ISEGBC,INKPBC,KPBR0C,NCF,IFEPC,INHPC
C
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
C  DATA FOR THE VARIATION OF NOISE COEFFICIENT
C
      DATA EPSR/1.D4/,EPSAP/10.D0/,EPSAG/1.D3/
      DATA EPSMAX/1.D15/
C
      IFEP = 1
C
C  INITIALIZE COMMON AREA  /DINCOM/
C
      CALL INIT(N,X0,H,EPS,DX,IRAND,F,
     1    NTRAJ,ISEGBR,INKPBR,KPBR0,INHP,IFEP,XRMIN,XRMAX,IOUT)
C
C  CHECK PARAMETER VALUES
C
      IF(NPMIN.LE.0.OR.NPMAX0.LT.0.OR.INPMAX.LE.0.OR.
     1    NSUC.LE.0.OR.NTRIAL.LE.0)IOUT = -99
      IF (IOUT.EQ.(-99))RETURN
C
C  INITIALIZE VARIABLES
C
      EPSC = EPS
      NPMAX = NPMAX0
      ISTOPT = 0
      ISTOP = 0
      ICCOM = (NTRIAL+NTRIAL+4)/5
      NFEV = 0
      NTES = NSUC
      ICTS = NSUC-1
      FTFOPT = F
C
C  START SERIES OF TRIALS
C
      DO 30 IC = 1,NTRIAL
C
C  SET INITIALIZATION INDEX FOR NOISE GENERATOR
C
         IS = 0
         IF(IRAND.GT.0)IS = IRAND+IC-1
C
C  INITIALIZE TRIAL
C
         IF(IC.GT.1.AND.IC.LE.ICCOM)CALL REINIT(N,X0,EPSC,IS,F,IFEP)
         IF(IC.GT.ICCOM)CALL REINIT(N,XMIN,EPSC,IS,FOPT,IFEP)
C
         FTFMIN = F
         FTFMAX = F
         NFEV = NFEV+1
C
C  PRINT INITIAL CONDITIONS OF TRIAL
C
         IF(IPRINT.GT.0)CALL PTSEG(N,X0,FTFMIN,FTFMAX,0,NFEV)
C
C
C  DEACTIVATE SCALING
C
         IF(KPASCA.GT.NPMAX.OR.N.LE.1)CALL NOSCA
C
C  INITIALIZE COMMON AREA  /SCALE/
C
         IF(KPASCA.LE.NPMAX.AND.N.GT.1)CALL INISCA(N,NTRAJ)
C
C  PERFORM A TRIAL
C
         CALL TRIAL(N,NPMIN,NPMAX,KPASCA,TOLREL,TOLABS,
     1       IPRINT,XMIN,FTFMIN,
     1       FTFMAX,NFEV,KP,ISTOP)
C
C
C     EVALUATE PAST TRIAL AND PREPARE NEXT TRIAL
C
C
C  SET TRIAL DURATION
C
         IF(ISTOP.EQ.0)NPMAX = NPMAX+INPMAX
C
C  RESET CURRENT NUMBER OF SUCCESSES REQUIRED BEFORE STOPPING
C  COMPUTE INDICATOR OF TRIAL STOPPING CONDITIONS
C  UPDATE BEST CURRENT VALUES OF TRIAL STOPPING INDICATOR AND
C  OF FUNCTION  F(X)  AT TRIAL STOP
C
         IF((FTFMIN.GT.FTFOPT.OR.ISTOP.EQ.0).AND.ISTOPT.NE.0)GO TO 10
         IF(ITOLCH(FTFOPT,FTFMIN,TOLREL,TOLABS).EQ.0)NTES = NSUC
C
         FTFOPT = FTFMIN
         ISTOPT = ISTOP
  10     CONTINUE
         ISTOPT = IABS(ISTOPT)
         CALL RCLOPT(N,XMIN,FOPT)
         IF(ITOLCH(FTFOPT,FOPT,TOLREL,TOLABS).EQ.0)ISTOPT = -ISTOPT
         IF(ITOLCH(FTFMIN,FOPT,TOLREL,TOLABS).EQ.0)ISTOP = -ISTOP
C
C  END-OF-TRIAL PRINT OUT
C
         IF(IPRINT.GE.0)
     1      CALL PTRIAL ( N, XOPT, FOPT,
     2                    FTFMIN, FTFMAX, FTFOPT,
     3                    ISTOP, ISTOPT, NFEV, KP, IPRINT )
C
C  UPDATE INITIAL VALUE OF NOISE COEFFICIENT FOR NEXT TRIAL
C
         IF (ISTOP.EQ.0)EPSC = EPSC/EPSR
         IF(ISTOP.GT.0)EPSC = EPSC*EPSAG
         IF(ISTOP.LT.0.AND.IC.LE.ICCOM)EPSC = EPSC*EPSAP
         IF(ISTOP.LT.0.AND.IC.GT.ICCOM)EPSC = EPSC/EPSR
C
C  UPDATE OPERATING CONDITIONS FOR SELECTING (IN THE NEXT TRIAL)
C  THE STARTING VALUE OF THE NOISE COEFFICIENT OF THE NEW TRAJECTORY
C  AFTER BRANCHING
C
         IF (ISTOP.EQ.0)IFEP = 1
         IF (ISTOP.NE.0)IFEP = 2
C
C  UPDATE, PRINT, AND CHECK TRIAL STOPPING CONDITIONS
C
         IF(ISTOP.GT.0.AND.ISTOPT.GT.0)NTES = NTES-1
         IF(NTES.LE.0.OR.ISTOPT.LE.0.OR.ISTOP.LE.0.OR.NTES.GT.ICTS.OR.
     1       IPRINT.LT.0) GO TO 20
         KSUC = NSUC-ICTS
         CALL PTKSUC(KSUC)
         ICTS = NTES-1
  20     CONTINUE
         IF(NTES.LE.0.AND.ISTOPT.GT.0.AND.ISTOP.GT.0)GO TO 40
C
C  CONSTRAIN NOISE COEFFICIENT WITHIN BOUNDS
C
         EPSC = DMIN1(EPSC,EPSMAX)
         IF(EPSC.LE.0.D0)EPSC = 1.D0
  30     CONTINUE
C  END OF SERIES OF TRIALS
C
  40  CONTINUE
C
C  INDICATOR OF STOPPING CONDITIONS
C
      IOUT = ISTOPT
      FMIN = FOPT
C
      RETURN
C
      END
      SUBROUTINE INIT(NX,X0,H0,EPS0,DX0,IRAND,F,N1,N2,N3,N4
     1    ,INH,IFE,XRI,XRA,IT)
C
C  INIT  PERFORMS THE INITIALIZATION OF THE FIRST TRIAL.
C  THE PART OF THE INITIALIZATION WHICH IS COMMON ALSO TO THE TRIALS
C  FOLLOWING THE FIRST ONE IS PERFORMED BY CALLING THE SUBROUTINE
C  REINIT.
C
      DOUBLE PRECISION X0,H0,EPS0,DX0,F,XRI,XRA
      DOUBLE PRECISION FUNCT
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION X0(NX),XRI(NX),XRA(NX)
      DATA NPMAX/100/
      DATA NTRAJM/20/
      DATA NTRAJ0/7/
      DATA INKPB0/10/
      DATA KPBR00/3/
C
C CHECK PARAMETER VALUES
C
      IT = 0
      IF (NX.GT.NPMAX.OR.NX.LT.1.OR.H0.LE.0.D0.OR.EPS0.LE.0.D0
     1    .OR.DX0.LE.0.D0) IT = -99
      IF (IT.EQ.(-99)) RETURN
C
C INITIALIZE SOME VARIABLES
C
      INHP = INH
      DO 10 IX = 1,NX
         XRMIN(IX) = XRI(IX)
         XRMAX(IX) = XRA(IX)
  10     CONTINUE
      CALL NOSCA
      NTRAJ = N1
      IF(NTRAJ.EQ.0)NTRAJ = NTRAJ0
      NTRAJ = MIN0(NTRAJ,NTRAJM)
      NTRAJ = MAX0(NTRAJ,3)
      N1 = NTRAJ
      ISEGBR = N2
      IF(ISEGBR.EQ.0)ISEGBR = (1+NTRAJ)/2
      ISEGBR = MIN0(ISEGBR,NTRAJ)
      ISEGBR = MAX0(ISEGBR,1)
      N2 = ISEGBR
      INKPBR = N3
      IF(INKPBR.EQ.0)INKPBR = INKPB0
      N3 = INKPBR
      KPBR0 = N4
      IF(KPBR0.EQ.0)KPBR0 = KPBR00
      KPBR0 = MOD(KPBR0,INKPBR)
      N4 = KPBR0
      NDIM = NX
      NTRAJR = NTRAJ-1
      NCF = 1
      F = FUNCT(NX,X0)
      CALL STOOPT(NX,X0,F)
      DO 20 ID = 1,NTRAJ
         H(ID) = H0
         DX(ID) = DX0
  20     CONTINUE
C
C INITIALIZE REMAINING VARIABLES
C
      CALL REINIT(NX,X0,EPS0,IRAND,F,IFE)
C
      RETURN
      END
      SUBROUTINE REINIT(NX,X0,EPS0,IRAND,F,IFE)
C
C  REINIT  PERFORMS THE INITIALIZATION OF ALL TRIALS FOLLOWING THE
C  FIRST TRIAL, AND PART OF THE INITIALIZATION OF THE FIRST TRIAL.
C
      DOUBLE PRECISION X0,EPS0,F
      DOUBLE PRECISION EPSV,G
      DOUBLE PRECISION CHAOS
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION X0(NX)
      DATA EPSV/1.D0/
C
C INITIALIZE RANDOM NOISE GENERATOR
C
      G = CHAOS(IRAND)
C
      IFEP = IFE
      KGEN = 0
      KTIM = 0
      DO 30 ID = 1,NTRAJ
         DO 10 IC = 1,NDIM
            X(IC,ID) = X0(IC)
  10        CONTINUE
         IE(ID) = 0
         DO 20 IT = 1,NTRAJR
            ISVT(ID,IT) = 1
            VMVT(ID,IT) = F
  20        CONTINUE
         EPS(ID) = EPS0*EPSV**(ISEGBR-ID)
         VMCOR(ID) = F
         VCOR(ID) = F
  30     CONTINUE
      RETURN
      END
      SUBROUTINE TRIAL(N,NPMIN,NPMAX,KPASCA,
     1    TOLREL,TOLABS,IPRINT,XMIN,FMIN,
     2    FMAX,NFEV,NR,ISTOP)
C
C  THE SUBROUTINE  TRIAL  GENERATES A TRIAL, I.E. A SET OF COMPLETE
C  SIMULTANEOUS TRAJECTORIES BY REPEATEDLY PERFORMING
C     A CALL TO THE SUBROUTINE  GENEVA  WHICH GENERATES THE SET OF
C        SIMULTANEOUS TRAJECTORY SEGMENTS CORRESPONDING TO THE CURRENT
C        OBSERVATION PERIOD, AND PERFORMS THE TRAJECTORY SELECTION
C     A (POSSIBLE) CALL TO THE SUBROUTINE  PTSEG  WHICH PERFORMS
C        END-OF-SEGMENT OUTPUT
C     A CHECK OF THE TRIAL STOPPING CRITERIA (WITH THE AID OF THE
C        FUNCTION  ITOLCH
C     A DECISION ABOUT ACTIVATING OR DEACTIVATING THE SCALING OF
C        THE VARIABLES (ACTIONS PERFORMED BY CALLING THE SUBROUTINES
C        ACTSCA  AND  NOSCA).
C
      DOUBLE PRECISION TOLREL,TOLABS,XMIN,FMIN,FMAX
      DIMENSION XMIN(N)
      DATA IRNF/7/
C
      DO 20 IR = 1,NPMAX
C
C ACTIVATE SCALING
C
         IF(IR.GE.KPASCA.AND.IR.GT.N*IRNF)CALL ACTSCA
         NR = IR
C
C GENERATE AND EVALUATE THE SIMULTANEOUS TRAJECTORY SEGMENTS
C PERIOD
C
         CALL GENEVA(N,XMIN,FMIN,FMAX,NFEV)
C
C PRINT RESULTS OF OBSERVATION PERIOD
C
         IF(IPRINT.GT.0)CALL PTSEG(N,XMIN,FMIN,FMAX,IR,NFEV)
C
C CHECK TRIAL STOPPING CONDITIONS
C
         IF(IR.LT.NPMIN)GO TO 10
         ISTOP = ITOLCH(FMAX,FMIN,TOLREL,TOLABS)
         IF(ISTOP.NE.0)GO TO 30
C
  10     CONTINUE
  20     CONTINUE
  30  CONTINUE
      CALL NOSCA
C
      RETURN
      END
      SUBROUTINE GENEVA(NX,XMIN,FMIN,FMAX,NCEF)
C
C  GENEVA  PERFORMS THE GENERATION AND THE FINAL PROCESSING AND
C  EVALUATION OF THE SET OF TRAJECTORY SEGMENTS CORRESPONDING TO
C  THE CURRENT OBSERVATION PERIOD.
C  GENEVA  UPDATES THE SCALING ARRAYS  DIST  AND  BIAS  BY CALLING
C             THE SUBROUTINES  SEGSCA  AND  UPDSCA
C          GENERATES THE TRAJECTORY SEGMENTS BY CALLING THE SUB-
C             ROUTINE  PERIOD
C          DETERMINES SOME END-OF-SEGMENT RESULTS  (FPFMIN, FPFMAX,
C             XPFMIN)  USING THE RESCALING CAPABILITIES OF THE SUB-
C             ROUTINES  SEGSCA  AND  VARSCA.
C
      DOUBLE PRECISION XMIN,FMIN,FMAX
      DOUBLE PRECISION FM
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XMIN(NX)
C
C UPDATE SCALING DATA
C
      DO 10 ID = 1,NTRAJ
         CALL SEGSCA(ID)
         CALL UPDSCA(NX,X(1,ID))
  10     CONTINUE
C
C GENERATE THE SIMULTANEOUS TRAJECTORY SEGMENTS
C
      CALL PERIOD
C
C EXTRACT AND RESCALE SOME FINAL VALUES
C
      FM = VCOR(1)
      FMAX = VCOR(1)
      IFM = 1
      DO 30 ID = 2,NTRAJ
         FMAX = DMAX1(FMAX,VCOR(ID))
         IF(VCOR(ID).GE.FM)GO TO 20
         FM = VCOR(ID)
         IFM = ID
  20     CONTINUE
  30     CONTINUE
      DO 40 IC = 1,NDIM
         XMIN(IC) = X(IC,IFM)
  40     CONTINUE
      FMIN = FM
      NCEF = NCF
      CALL SEGSCA(IFM)
      CALL VARSCA(NX,XMIN)
C
      RETURN
      END
      SUBROUTINE PERIOD
C
C  PERIOD  IS CALLED BY SUBROUTINE  GENEVA  TO PERFORM THE GENERATION
C  OF THE TRAJECTORY SEGMENTS.
C  PERIOD  - COMPUTES THE DURATION OF THE OBSERVATION PERIOD, I.E. THE
C              NUMBER  NHP  OF ACCEPTED INTEGRATION STEPS IN A PERIOD
C          - COMPUTES ALL THE SEGMENT STEPS BY REPEATEDLY CALLING THE
C              SUBROUTINE  STEP
C          - PERFORMS THE SEGMENT SELECTION BY CALLING THE SUBROUTINE
C               BRASI
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
C DETERMINE DURATION OF OBSERVATION PERIOD
C (NUMBER OF TIME INTEGRATION STEPS)
C
      KGEN = KGEN+1
      NKGEN = 1
      IF(INHP.EQ.1)NKGEN=
     1    INT(SNGL(DLOG(DBLE(FLOAT(KGEN))+.5D0)/DLOG(2.D0)))+1
      IF(INHP.EQ.2)NKGEN = INT(SNGL(DSQRT(DBLE(FLOAT(KGEN))+.5D0)))
      IF(INHP.EQ.3)NKGEN = KGEN
C
C PERFORM THE REQUIRED NUMBER OF INTEGRATION STEPS
C (FOR ALL TRAJECTORIES)
C
      DO 10 IK = 1,NKGEN
C
C PERFORM A SINGLE INTEGRATION STEP
C (FOR ALL TRAJECTORIES)
C
         CALL STEP(IK)
  10     CONTINUE
C
C MANAGE THE TRAJECTORY BRANCHING PROCESS
C
      CALL BRASI
      RETURN
      END
      SUBROUTINE BRASI
C
C  BRASI  PERFORMS THE SELECTION PROCESS FOR THE TRAJECTORIES
C  BRASI  - UPDATES THE DATA ABOUT THE PAST TRAJECTORIES
C         - ASKS FOR THE TRAJECTORY-SELECTION ORDERING BY CALLING THE
C              SUBROUTINE  ORDER
C         - DISCARDS ONE OF THE TRAJECTORIES
C         - PERFORMS BRANCHING ON ONE OF THE REMAINING TRAJECTORIES
C         - MOVES THE DATA OF THE UNPERTURBED CONTINUATION
C              TO THE POSITION OF THE PERTURBED CONTINUATION
C         - CALLS THE SUBROUTINE  COMPAS  TO EXAMINE DATA ABOUT PAST
C              HISTORY OF THE TRAJECTORIES AND DISCARD IRRILEVANT DATA
C
      DOUBLE PRECISION DFAC,DLEPMX,DLFACL
      DOUBLE PRECISION EFAC,EPSMAX,FACG
      DOUBLE PRECISION CHAOS
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      DATA FACG/10.D0/
      DATA EFAC/.5D0/
      DATA DFAC/1.D3/
      DATA EPSMAX/1.D15/
      DATA DLFACL /0.301029995663981194D0/
      DATA DLEPMX /25.D0/
C
C UPDATE PAST HISTORY DATA
C
      DO 10 ID = 1,NTRAJ
         VMVT(ID,NTRAJR) = VMCOR(ID)
         ISVT(ID,NTRAJR) = ID
  10     CONTINUE
C
C OBTAIN TRAJECTORY-SELECTION ORDERING
C
      CALL ORDER(IP,IM,IU)
C
C DECIDE WHICH TRAJECTORY IS TO BE BRANCHED
C
      IF(MOD(KGEN,INKPBR).EQ.KPBR0)IM = IP
C
C PERFORM BRANCHING
C
      DO 20 IC = 1,NDIM
         X(IC,IU) = X(IC,IM)
  20     CONTINUE
      H(IU) = H(IM)
      IE(IU) = IE(IM)
      DO 30 IT = 1,NTRAJR
         ISVT(IU,IT) = ISVT(IM,IT)
         VMVT(IU,IT) = VMVT(IM,IT)
  30     CONTINUE
      EPS(IU) = EPS(IM)
      DX(IU) = DX(IM)
      VCOR(IU) = VCOR(IM)
      DO 40 ID = 1,NTRAJ
         VMCOR(ID) = VCOR(ID)
  40     CONTINUE
C
C UPDATE PAST-HISTORY-DATA MATRICES
C
      CALL COMPAS
C
C UPDATE SCALING DATA
C
      CALL MOVSCA(IU,IM)
C
C UPDATE NOISE COEFFICIENT VALUES
C
      IF (EPS(IU).LE.0.D0) GO TO 50
      IF(IFEP.EQ.2)
     1    EPS(IU) = EPS(IU)*FACG**(CHAOS(0)-EFAC)
      IF(IFEP.EQ.1)
     1    EPS(IU) = FACG**(DMIN1(DLEPMX,
     1    DLOG10(EPS(IU))+(CHAOS(-1)-EFAC)*DLFACL))
      EPS(IU) = DMIN1(EPS(IU),EPSMAX)
  50  CONTINUE
C
C UPDATE MAGNITUDE OF SPATIAL DISCRETIZATION INCREMENT
C
      DX(IU) = DX(IU)*DFAC**CHAOS(0)
      RETURN
      END
      SUBROUTINE ORDER(IP,IM,IU)
C
C  ORDER  COMPARES THE TRAJECTORIES FROM THE POINT OV VIEW OF PAST
C  HISTORY (BY CALLING THE FUNCTION  IPREC ) AND FROM THE POINT OF VIEW
C  OF THEIR CURRENT VALUE OF  EPS  (BY CALLING THE FUNCTION  IPRECE)
C  AND PROVIDES THE TRAJECTORY ORDERING NEEDED FOR SELECTING THE TRA-
C  JECTORY WHICH IS TO BE BRANCHED.
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION IORD(20)
      DATA KP/0/
      IR = 0
C
C ASSIGN INITIAL ORDERING
C
  10  CONTINUE
      DO 20 I = 1,NTRAJ
         IORD(I) = I
  20     CONTINUE
C
C SORT TRAJECTORIES ...
C
  30  CONTINUE
      IR = IR+1
C
      DO 60 I = 1,NTRAJR
         I1 = I+1
         DO 50 J = I1,NTRAJ
            K1 = IORD(I)
            K2 = IORD(J)
C
C ... ACCORDING TO PAST HISTORY ...
C
            IF(IR.NE.2)KP = IPREC(K1,K2)
            IF(KP.EQ.0.AND.IR.EQ.1)GO TO 10
C
C ... OR ACCORDING TO NOISE LEVEL
C
            IF(IR.EQ.2)KP = IPRECE(K1,K2)
            IF(KP.EQ.0)GO TO 40
            KM = K1+K2-KP
            IORD(I) = KP
            IORD(J) = KM
  40        CONTINUE
  50        CONTINUE
  60     CONTINUE
      IF(IR.EQ.2)GO TO 30
C
C RETURN THE INDICES OF THE SEGMENTS WHICH IN THE ORDERING
C OCCUPY THE FIRST, THE LAST, AND A SUITABLE MEDIUM LEVEL POSITION
C
      IP = IORD(1)
      IU = IORD(NTRAJ)
      IM = IORD(ISEGBR)
C
      RETURN
      END
      INTEGER FUNCTION IPREC(ID1,ID2)
C
C  IPREC  DETERMINES THE PRECEDENCE RELATION BETWEEN TWO TRAJECTORIES
C  BASED ON THE PAST HISTORY DATA
C
      DOUBLE PRECISION VM1,VM2
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      VM1 = VMVT(ID1,NTRAJR)
      VM2 = VMVT(ID2,NTRAJR)
      DO 10 IIT = 1,NTRAJR
         IT = 1+NTRAJR-IIT
         IF(ISVT(ID1,IT).EQ.ISVT(ID2,IT))GO TO 20
         VM1 = DMIN1(VM1,VMVT(ID1,IT))
         VM2 = DMIN1(VM2,VMVT(ID2,IT))
  10     CONTINUE
  20  CONTINUE
      IPREC = 0
      IF(VM2.LT.VM1)IPREC = ID2
      IF(VM2.GT.VM1)IPREC = ID1
C
      RETURN
      END
      INTEGER FUNCTION IPRECE(ID1,ID2)
C
C  IPRECE  DETERMINES THE PRECEDENCE RELATION BETWEEN TWO TRAJECTORIES
C  BASED ON THEIR CURRENT VALUE OF  EPS
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      IPRECE = 0
      IF(KGEN.GT.ISEGBR*INKPBR)GO TO 10
      IF(EPS(ID2).LT.EPS(ID1))IPRECE = ID1
      IF(EPS(ID2).GT.EPS(ID1))IPRECE = ID2
      RETURN
C
  10  CONTINUE
      IF(EPS(ID2).LT.EPS(ID1))IPRECE = ID2
      IF(EPS(ID2).GT.EPS(ID1))IPRECE = ID1
C
      RETURN
      END
      SUBROUTINE COMPAS
C
C  COMPAS  TAKES CARE OF THE STORAGE OF PAST HISTORY DATA, DISCARDING
C  ALL DATA NOT NEEDED BY THE ONLY USER OF SUCH DATA, THE FUNCTION
C  IPREC .
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      IT1 = 1
      DO 30 IT = 2,NTRAJR
         IT1 = IT-1
         DO 10 ID = 1,NTRAJ
            IF(ISVT(ID,IT).NE.ISVT(ID,IT1))GO TO 20
  10        CONTINUE
         GO TO 120
  20     CONTINUE
  30     CONTINUE
      IT1 = 1
      NCV = 0
      DO 50 ID = 1,NTRAJ
         DO 40 IDD = 1,NTRAJ
            IF(ISVT(ID,IT1).EQ.ISVT(IDD,IT1))NCV = NCV+1
  40        CONTINUE
  50     CONTINUE
      DO 80 IT = 2,NTRAJR
         IT1 = IT-1
         NCN = 0
         DO 70 ID = 1,NTRAJ
            DO 60 IDD = 1,NTRAJ
               IF(ISVT(ID,IT).EQ.ISVT(IDD,IT))NCN = NCN+1
  60           CONTINUE
  70        CONTINUE
         IF(NCN.EQ.NCV)GO TO 110
         NCV = NCN
  80     CONTINUE
      DO 90 ID = 1,NTRAJ
         IF(ISVT(1,1).NE.ISVT(ID,1))GO TO 100
  90     CONTINUE
      IT = 2
      GO TO 140
 100  CONTINUE
 110  CONTINUE
 120  CONTINUE
      IT = IT1+1
      DO 130 ID = 1,NTRAJ
         VMVT(ID,IT) = DMIN1(VMVT(ID,IT),VMVT(ID,IT1))
 130     CONTINUE
 140  CONTINUE
      DO 160 ITC = IT,NTRAJR
         ITM = ITC-1
         DO 150 ID = 1,NTRAJ
            VMVT(ID,ITM) = VMVT(ID,ITC)
            ISVT(ID,ITM) = ISVT(ID,ITC)
 150        CONTINUE
 160     CONTINUE
C
      RETURN
      END
      SUBROUTINE STEP(IK)
C
C  STEP  PERFORMS A SINGLE TIME-INTEGRATION STEP FOR EACH ONE OF THE
C  SIMULTANEOUS TRAJECTORIES BY REPEATEDLY CALLING THE SUBROUTINE  SSTEP
C
      DOUBLE PRECISION F
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      DOUBLE PRECISION XID,HID,EPSID,DXID
      DIMENSION XID(100)
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      KTIM = KTIM+1
C
C LOOP ON THE SIMULTANEOUS TRAJECTORY SEGMENTS
      DO 30 ID = 1,NTRAJ
C
C INFORM THE SCALING SUBPROGRAMS (VIA THE COMMON AREA /SCALE/)
C THAT SCALING OPERATIONS ARE TO BE PERFORMED ON SEGMENT  ID .
         CALL SEGSCA(ID)
         F = VCOR(ID)
C
C PERFORM A TIME-INTEGRATION STEP ON SEGMENT  ID .
         KA = KTIM
         NX = NDIM
         DO 10 IX = 1,NX
            XID(IX) = X(IX,ID)
  10        CONTINUE
         HID = H(ID)
         EPSID = EPS(ID)
         DXID = DX(ID)
         IEID = IE(ID)
C
         CALL SSTEP(KA,NX,XID,HID,EPSID,DXID,IEID,F)
C
         DO 20 IX = 1,NX
            X(IX,ID) = XID(IX)
  20        CONTINUE
         H(ID) = HID
         EPS(ID) = EPSID
         DX(ID) = DXID
         IE(ID) = IEID
         VCOR(ID) = F
         VMCOR(ID) = DMIN1(VMCOR(ID),F)
         IF(IK.EQ.1)VMCOR(ID) = F
         IF(KTIM.EQ.1)IE(ID) = 0
  30     CONTINUE
C
      RETURN
      END
      SUBROUTINE SSTEP(KTIM,NDIM,X,H,EPS,DX,IE,F)
C
C  THE BASIC TIME-INTEGRATION STEP FOR A GIVEN TRAJECTORY IS PERFORMED
C  BY THE SUBROUTINE  STEP  WHICH
C   - CALLS THE FUNCTION  FUNCT0  TO COMPUTE THE VALUE OF  F
C   - CALLS THE SUBROUTINE  UNITRV  TO COMPUTE THE RANDOM DIRECTION
C        ALONG WHICH THE DIRECTIONAL DERIVATIVE  GAM/N  IS TO BE
C        COMPUTED
C   - CALLS THE SUBROUTINE  DERFOR (OR DERCEN ) TO COMPUTE THE FORWARD-
C        (OR CENTRAL-) DIFFERENCES DIRECTIONAL DERIVATIVE  GAM/N
C   - CALLS THE SUBROUTINE  NEWH  TO ACCEPT OR REJECT THE FIRST HALF-
C        STEP AND OBTAIN AN UPDATED VALUE FOR  H
C   - CALLS THE SUBROUTINE  CUMSCA  TO UPDATE THE CUMULATED SCALING DATA
C   - UPDATES THE SPATIAL DISCRETIZATION INCREMENT  DX  BASED ON THE
C        RESULTS OF CALLING THE FUNCTION  ITOLCH
C   - CALLS THE SUBROUTINE  GAUSRV  TO PERFORM THE SECOND HALF-STEP.
C
      DOUBLE PRECISION X,H,EPS,DX,F
      DOUBLE PRECISION DFDX,DFDXV,DXMAX,DXMIN
      DOUBLE PRECISION EPSR,FS,FV,FVS,HMINS,HR,HS
      DOUBLE PRECISION RCD,RDX,STF,TOLABS,TOLRA,TOLRI
      DOUBLE PRECISION W,XP
      DOUBLE PRECISION FUNCT0
      DIMENSION X(NDIM)
      DIMENSION W(100),XP(100)
      DATA RDX/1.D-4/
      DATA DXMIN/1.D-35/
      DATA DXMAX/1.D3/
      DATA HR/1.D-1/
      DATA HMINS/1.D-30/
      DATA STF /100.D0/
      DATA RCD/2.D0/
      DATA TOLRI/1.D-5/
      DATA TOLRA/1.D-11/
      DATA TOLABS/0.D0/
C
      IEC = 0
      FV = F
  10  CONTINUE
  20  CONTINUE
C
C TAKE A RANDOM DIRECTION FOR THE DIRECTIONAL DERIVATIVE
      CALL UNITRV(NDIM,W)
C
C COMPUTE FORWARD-DIFFERENCE DERIVATIVE
      CALL DERFOR(NDIM,X,FV,DX,W,DFDX)
C
C TRY THE FIRST HALF-STEP
      DO 30 IC = 1,NDIM
         XP(IC) = X(IC)-H*W(IC)*DFDX*DBLE(FLOAT(NDIM))
  30     CONTINUE
      HS = H
      F = FUNCT0(NDIM,XP)
      FVS = FV+DX*DABS(DFDX)
C
C UPDATE STEPLENGTH  H  AND ACCEPT OR REJECT THE HALF-STEP
      CALL NEWH(KTIM,FVS,F,H,IE,IEC)
      IF(IEC.LE.0)GO TO 50
      IE = IE-1
      IEC = IEC-1
      H = HS
      DFDXV = DFDX
C
C COMPUTE CENTRAL-DIFFERENCES DERIVATIVE
      CALL DERCEN(NDIM,X,FV,DX,W,DFDX)
C
C TRY AGAIN THE FIRST HALF-STEP
      DO 40 IC = 1,NDIM
         XP(IC) = X(IC)-H*W(IC)*DFDX*DBLE(FLOAT(NDIM))
  40     CONTINUE
      F = FUNCT0(NDIM,XP)
      FVS = FV+DX*DABS(DFDXV-DFDX)
C
C
C UPDATE STEPLENGTH  H  AND ACCEPT OR REJECT THE HALF-STEP
      CALL NEWH(KTIM,FVS,F,H,IE,IEC)
C
C UPDATE CUMULATED PAST SCALING DATA
      IF(IEC.GE.1)CALL CUMSCA(NDIM,W,DFDX)
      IF(IEC.GE.1)GO TO 10
      DX = DX*RDX
C
C FIRST HALF-STEP ACCEPTED
  50  CONTINUE
C
C UPDATE CUMULATED PAST SCALING DATA
      CALL CUMSCA(NDIM,W,DFDX)
      FS = FV+DX*DABS(DFDX)
      IF(ITOLCH(FS,FV,TOLRI,TOLABS).EQ.0)DX = DX/RCD
      IF(ITOLCH(FS,FV,TOLRA,TOLABS).GT.0)DX = DX*RCD
      EPSR = DSQRT(H)*EPS
C
C TAKE A SAMPLE INCREMENT OF THE WIENER PROCESS
      CALL GAUSRV(NDIM,W)
C
C TRY THE SECOND HALF-STEP
      DO 60 IC = 1,NDIM
         XP(IC) = XP(IC)+EPSR*W(IC)
  60     CONTINUE
      F = FUNCT0(NDIM,XP)
C
C ACCEPT OR REJECT THE COMPLETE STEP
      IF (F-FV.LE.EPS*EPS*STF) GO TO 70
      H = H*HR
      IE = IE+1
      IF (H.GT.HMINS) GO TO 20
C
C STEP ACCEPTED
  70  CONTINUE
      DO 80 IC = 1,NDIM
         X(IC) = XP(IC)
  80     CONTINUE
      DX = DMIN1(DX,DXMAX)
      DX = DMAX1(DX,DXMIN)
C
      RETURN
      END
      SUBROUTINE NEWH(K,FV,F,H,IE,IEC)
C
C  NEWH  IS CALLED BY THE SUBROUTINE  SSTEP  TO DECIDE WHETHER TO ACCEPT
C  OR REJECT THE FIRST HALF-STEP, AND TO PROVIDE AN UPDATED VALUE FOR  H
C
      DOUBLE PRECISION FV,F,H
      DOUBLE PRECISION FA,FR,HMAX,HMIN,R
      DIMENSION FR(3),FA(4)
      DATA FR/1.05D0,2.D0,10.D0/
      DATA FA/1.D0,1.1D0,2.D0,10.D0/
      DATA HMIN/1.D-30/
      DATA HMAX/1.D25/
      DATA IECMAX/50/
C
      IF(FV.LT.F)GO TO 20
C
C STEP ACCEPTED, POSSIBLY INCREASE THE STEPLENGTH  H
C
      R = FA(1)
      IF(IE*2.LT.K)R = FA(2)
      IF(IE*3.LT.K)R = FA(3)
      IF(IE.EQ.0.AND.K.GT.1)R = FA(4)
      H = H*R
  10  CONTINUE
      IEC = 0
      GO TO 30
C
  20  CONTINUE
      IE = IE+1
      IEC = IEC+1
      IF(IEC.GT.IECMAX)GO TO 10
C
C STEP REJECTED, DECREASE  H
C
      IC = MIN0(3,IEC)
      R = FR(IC)
      H = H/R
  30  CONTINUE
      H = DMIN1(H,HMAX)
      H = DMAX1(H,HMIN)
C
      RETURN
      END
      SUBROUTINE DERFOR(NDIM,X,F,DX,W,DFDX)
C
C  DERFOR  COMPUTED THE FORWARD FINITE-DIFFERNCES DIRECTIONAL
C  DERIVATIVES (CALLING  FUNCT0 )
C
      DOUBLE PRECISION X,F,DX,W,DFDX
      DOUBLE PRECISION DXF,DXFF,DXMAX,FN,S,XX
      DOUBLE PRECISION FUNCT0
      DIMENSION X(NDIM),W(NDIM)
      DIMENSION XX(100)
      DATA DXFF/1.D6/,DXF/1.D1/
      DATA DXMAX/1.D6/
C
  10  CONTINUE
  20  CONTINUE
      S = 0.D0
      DO 30 IC = 1,NDIM
         XX(IC) = X(IC)+W(IC)*DX
         S = S+(XX(IC)-X(IC))**2
  30     CONTINUE
      IF(S.GT.0.D0)GO TO 40
      DX = DX*DXFF
      GO TO 20
  40  CONTINUE
      FN = FUNCT0(NDIM,XX)
      DFDX = (FN-F)/DX
      IF(DX.GT.DXMAX)RETURN
      IF(DABS(DFDX).GT.1.D0)RETURN
      IF(DFDX**2.GT.0.D0)GO TO 50
      DX = DX*DXF
      GO TO 10
  50  CONTINUE
C
      RETURN
      END
      SUBROUTINE DERCEN(NDIM,X,F,DX,W,DFDX)
C
C  DERFOR  COMPUTED THE CENTRAL FINITE-DIFFERNCES DIRECTIONAL
C  DERIVATIVES (CALLING  FUNCT0 )
C
      DOUBLE PRECISION X,F,DX,W,DFDX
      DOUBLE PRECISION FN,FR,XX
      DOUBLE PRECISION FUNCT0
      DIMENSION X(NDIM),W(NDIM)
      DIMENSION XX(100)
C
      DO 10 IC = 1,NDIM
         XX(IC) = X(IC)-W(IC)*DX
  10     CONTINUE
      FR = FUNCT0(NDIM,XX)
      FN = F+DFDX*DX
      DFDX = (FN-FR)/(2.D0*DX)
C
      RETURN
      END
      SUBROUTINE RCLOPT(N,XO,FO)
C
C  RCLOPT  RECALLS THE CURRENT BEST MINIMUM VALUE  FOPT  FOUND SO FAR
C  FROM ALGORITHM START, AND THE POINT  XOPT  (OR POSSIBLY ONE OF THE
C  POINTS) WHERE  FOPT  WAS OBTAINED
C
      DOUBLE PRECISION XO,FO
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XO(N)
C
      DO 10 I = 1,N
         XO(I) = XOPT(I)
  10     CONTINUE
      FO = FOPT
C
      RETURN
      END
      SUBROUTINE STOOPT(N,XO,FO)
C
C  STOOPT  STORES THE CURRENT BEST MINIMUM VALUE  FOPT  FOUND SO FAR
C  FROM ALGORITHM START, AND THE POINT  XOPT  (OR POSSIBLY ONE OF THE
C  POINTS) WHERE  FOPT  WAS OBTAINED
C
      DOUBLE PRECISION XO,FO
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XO(N)
C
      DO 10 I = 1,N
         XOPT(I) = XO(I)
  10     CONTINUE
      FOPT = FO
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION FUNCT0(N,XX)
C
C  FUNCT0  IS CALLED  WHENEVER THE VALUE OF THE FUNCLION  F  IS REQUIRED
C  IN THE NUMERICAL INTEGRATION PROCESS.
C  THE FUNCTION  FUNCT0
C   - RESCALES THE VARIABLES BY CALLING THE SUBROUTINE  VARSCA
C   - CALLS THE SUBROUTINE  RANGE  TO TAKE CARE OF THE CASES WHERE THE
C      CURRENT POINT  X  IS OUTSIDE A GIVEN ADMISSIBLE REGION
C   - CALLS THE USER-SUPPLIED FUNCTION  FUNCT  TO ACTUALLY COMPUTE THE
C      VALUE OF  F
C   - POSSIBLY CALLS THE SUBROUTINE  STOOPT  TO UPDATE THE CURRENT
C        BEST MINIMUM FUNCTION VALUE  FOPT  AND THE CORRESPONDING
C        MINIMIZER  XOPT
C
      DOUBLE PRECISION XX
      DOUBLE PRECISION F,R,XS
      DOUBLE PRECISION FUNCT
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XX(N)
      DIMENSION XS(100)
C
      DO 10 IX = 1,N
         XS(IX) = XX(IX)
  10     CONTINUE
C
C DESCALE  X-VARIABLES
      CALL VARSCA(N,XS)
C
C CONSTRAIN THE  X-VARIABLES WITHIN BOUNDS
      CALL RANGE(N,XS,XRMIN,XRMAX,R)
C
C COMPUTE THE FUNCTION VALUE...
      F = FUNCT(N,XS)+R
C
C ... AND POSSIBLY UPDATE THE BEST CURRENT MINIMUM
      IF(F.LT.FOPT)CALL STOOPT(N,XS,F)
      FUNCT0 = F
      NCF = NCF+1
C
      RETURN
      END
      SUBROUTINE RANGE (N,XS,XRMIN,XRMAX,R)
C
C  RANGE  IS CALLED BY THE FUNCTION  FUNCT0  TO TAKE CARE OF THE CASES
C  WHERE THE CURRENT POINT  X  IS OUTSIDE A GIVEN ADMISSIBLE REGION
C
      DOUBLE PRECISION XS,XRMIN,XRMAX,R
      DOUBLE PRECISION A,B,C,D,DLRMAX,RMAX,RR,XC
      DIMENSION XS(N),XRMIN(N),XRMAX(N)
      DATA RMAX /1.D35/
      DATA DLRMAX/80.5904782547915990D0/
C
      R = 0.D0
      DO 40 I = 1,N
         A = XRMAX(I)
         C = XRMIN(I)
         XC = XS(I)
         IF (XC.LE.A) GO TO 10
         B = A+A-C
         RR = RMAX
         IF (XC.LT.B) RR = DEXP((XC-A)*DLRMAX/(B-A))-1.D0
         R = R+RR
         XS(I) = XRMAX(I)
         GO TO 30
  10     CONTINUE
         IF(XC.GE.C)GO TO 20
         D = C+C-A
         RR = RMAX
         IF (XC.GT.D) RR = DEXP((C-XC)*DLRMAX/(C-D))-1.D0
         R = R+RR
         XS(I) = XRMIN(I)
  20     CONTINUE
  30     CONTINUE
  40     CONTINUE
C
      RETURN
      END
      INTEGER FUNCTION ITOLCH(FMAX,FMIN,TOLREL,TOLABS)
C
C  ITOLCH  DETERMINES  WHETHER TWO QUANTITIES ARE TO BE CONSIDERED NU-
C  MERICALLY EQUAL WITH RESPECT TO AN ABSOLUTE (OR RELATIVE) DIFFERENCE
C  CRITERION, WITHIN GIVEN TOLERANCES
C
      DOUBLE PRECISION FMAX,FMIN,TOLREL,TOLABS
C
      ISTOP = 0
C
C CHECK RELATIVE DIFFERENCE AGAINST  TOLREL
      IF(DABS(FMAX-FMIN).LE.TOLREL*(DABS(FMIN)+DABS(FMAX))/2.D0)
     1    ISTOP=ISTOP+1
C
C CHECK ABSOLUTE DIFFERENCE AGAINST  TOLABS
      IF(FMAX-FMIN.LE.TOLABS)ISTOP = ISTOP+2
      ITOLCH = ISTOP
C
      RETURN
      END
      SUBROUTINE INISCA(N,ND)
C
C  INISCA  INITIALIZES THE COOMON AREA  /SCALE/  FOR THE SCALING DATA
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DATA NXMSCA/10/
C
      LSCA = -1
      IF(N.GT.NXMSCA.OR.N.EQ.1)RETURN
      LSCA = 0
      NX = N
      NORD = ND
      IDSCA = 1
      DO 1 ID = 1,NORD
         DO 2 IX = 1,NX
            DO 3 IY = 1,NX
               DIST(IX,IY,ID) = 0.D0
               GRAGRA(IX,IY,ID) = 0.D0
   3           CONTINUE
            DIST(IX,IX,ID) = 1.D0
            BIAS(IX,ID) = 0.D0
            GRA(IX,ID) = 0.D0
   2        CONTINUE
         NGRA(ID) = 0
   1     CONTINUE
C
      RETURN
      END
      SUBROUTINE NOSCA
C
C  NOSCA  DEACTIVATES THE SCALING
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      LSCA = -1
C
      RETURN
      END
      SUBROUTINE SEGSCA(ID)
C
C  SEGSCA  SELECTS THE TRAJECTORY WHICH MUST BE RESCALED
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      IDSCA = ID
C
      RETURN
      END
      SUBROUTINE VARSCA(N,X)
C
C  VARSCA  COMPUTES THE RESCALED VARIABLE  AX + B
C
      DOUBLE PRECISION X
      DOUBLE PRECISION XB
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION X(N),XB(10)
C
      IF(LSCA.LE.0)RETURN
      DO 1 I = 1,N
         XB(I) = BIAS(I,IDSCA)
         DO  2 J = 1,N
            XB(I) = XB(I)+DIST(I,J,IDSCA)*X(J)
   2        CONTINUE
   1     CONTINUE
      DO 3 I = 1,N
         X(I) = XB(I)
   3     CONTINUE
C
      RETURN
      END
      SUBROUTINE CUMSCA(N,W,DFDX)
C
C  CUMSCA  STORES CUMULATED STATISTICAL DATA ON THE ILL-CONDITIONING OF
C  F(AX+B) W.R.T. X
C
      DOUBLE PRECISION W,DFDX
      DOUBLE PRECISION DFDXMA
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION W(N)
      DATA DFDXMA/1.D8/
C
      IF(LSCA.LE.0)RETURN
      IF(DABS(DFDX).GT.DFDXMA)RETURN
      DO 1 I = 1,N
         DO 2 J = 1,N
            GRAGRA(I,J,IDSCA) = GRAGRA(I,J,IDSCA)+DFDX*W(I)*DFDX*W(J)
   2        CONTINUE
         GRA(I,IDSCA) = GRA(I,IDSCA)+DFDX*W(I)
   1     CONTINUE
      NGRA(IDSCA) = NGRA(IDSCA)+1
C
      RETURN
      END
      SUBROUTINE ACTSCA
C
C  ACTSCA  ACTIVATES THE RESCALING
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      IF(LSCA.EQ.0)LSCA = 1
C
      RETURN
      END
      SUBROUTINE MOVSCA(IU,IM)
C
C  MOVSCA  MOVES THE SCALING DATA OF THE UNPERTURBED CONTINUATION TO THE
C  POSITION OF THE PERTURBED CONTINUATION
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      IF(LSCA.LT.0)RETURN
      DO 1 I = 1,NX
         DO 2 J = 1,NX
            DIST(I,J,IU) = DIST(I,J,IM)
            GRAGRA(I,J,IU) = GRAGRA(I,J,IM)
   2        CONTINUE
         BIAS(I,IU) = BIAS(I,IM)
         GRA(I,IU) = GRA(I,IM)
   1     CONTINUE
      NGRA(IU) = NGRA(IM)
C
      RETURN
      END
      SUBROUTINE UPDSCA(N,X)
C
C  UPDSCA  UPDATES THE SCALING MATRIX  A  AND THE BIAS VECTOR  B BY
C  CALLING  EIGSCA  AND  VARSCA
C
      DOUBLE PRECISION X
      DOUBLE PRECISION AGRAI,ALA1,ALPHA,AMCOR,BIAST
      DOUBLE PRECISION CD,COR,DISTT,SN
      DOUBLE PRECISION EIGSCA
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION X(N)
      DIMENSION DISTT(10,10),BIAST(10),COR(10,10)
      DATA ALPHA /0.3D0/
C
      IF(LSCA.LE.0)RETURN
      ID = IDSCA
      IF(NGRA(ID).LT.2*NX*NX)GO TO 2
      AGRAI = 1.D0/DBLE(FLOAT(NGRA(ID)))
      AMCOR = 0.D0
      DO 3 I = 1,NX
         DO 4 J = 1,NX
         COR(I,J) = AGRAI*GRAGRA(I,J,ID)-AGRAI*GRA(I,ID)*AGRAI*GRA(J,ID)
            AMCOR = DMAX1(AMCOR,DABS(COR(I,J)))
   4        CONTINUE
   3     CONTINUE
      IF(AMCOR.LE.0.D0)GO TO 12
      DO 1 I = 1,NX
         DO 11 J = 1,NX
            COR(I,J) = COR(I,J)/AMCOR
  11        CONTINUE
   1     CONTINUE
      ALA1 = EIGSCA(COR)
      CD = ALA1*(1.D0+ALPHA)
      DO 5 I = 1,NX
         COR(I,I) = COR(I,I)-CD
         BIAST(I) = X(I)
   5     CONTINUE
      CALL VARSCA(NX,BIAST)
      SN = 0.D0
      DO 6 I = 1,NX
         DO 7 J = 1,NX
            DISTT(I,J) = 0.D0
            DO 8 K = 1,NX
               DISTT(I,J) = DISTT(I,J)-DIST(I,K,ID)*COR(K,J)
   8           CONTINUE
            SN = SN+DISTT(I,J)**2
   7        CONTINUE
   6     CONTINUE
      SN = 1.D0/DSQRT(SN/DBLE(FLOAT(NX)))
      DO 9 I = 1,NX
         BIAS(I,ID) = BIAST(I)
         DO 10 J = 1,N
            DIST(I,J,ID) = SN*DISTT(I,J)
            BIAS(I,ID) = BIAS(I,ID)-DIST(I,J,ID)*X(J)
            GRAGRA(I,J,ID) = 0.D0
  10        CONTINUE
         GRA(I,ID) = 0.D0
   9     CONTINUE
      NGRA(ID) = 0
   2  CONTINUE
  12  CONTINUE
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION EIGSCA(COR)
C
C  EIGSCA COMPUTES THE LARGEST EIGENVALUE OF A MATRIX USED FOR RESCA-
C  LING, STARTING FROM A RANDOMLY-CHOSEN ESTIMATE (OBTAINED BY CALLING
C  THE SUBROUTINE  UNITRV ) OF THE CORRESPONDING EIGENVECTOR
C
      DOUBLE PRECISION COR
      DOUBLE PRECISION ALA1,ALA1I,ALA1O
      DOUBLE PRECISION PREC,SWW,W,WW
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION COR(10,10)
      DIMENSION W(10),WW(10)
      DATA PREC /1.D-3/
      DATA NRMIN /10/
      DATA NRMAX /100/
C
      CALL UNITRV(NX,W)
      ALA1 = 0.D0
      DO 1 IR = 1,NRMAX
         ALA1O = ALA1
         SWW = 0.D0
         DO 2 IX = 1,NX
            WW(IX) = 0.D0
            DO 3 JX = 1,NX
               WW(IX) = WW(IX)+COR(IX,JX)*W(JX)
   3           CONTINUE
            SWW = SWW+WW(IX)**2
   2        CONTINUE
         ALA1 = DSQRT(SWW)
         ALA1I = 1.D0/ALA1
         IF(IR.GE.NRMIN.AND.ALA1*PREC.GT.DABS(ALA1-ALA1O))GO TO 4
         DO 5 IX = 1,NX
            W(IX) = WW(IX)*ALA1I
   5        CONTINUE
   1     CONTINUE
   4  CONTINUE
      EIGSCA = ALA1
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION CHAOS(INIZ)
C
C  CHAOS  GENERATES A RANDOM SAMPLE FROM ONE OUT OF FOUR POSSIBLE
C PROBABILITY DISTRIBUTIONS USING RANDOM NUMBERS UNIFORMLY
C DISTRIBUTED IN (0,1) GENERATED BY THE FUNCTION  UNIFRD
C
C  INIZ  IS AN INPUT PARAMETER AS FOLLOWS
C
C      INIZ>0          INITIALIZATION WITH SEED  INIZ.
C      INIZ=0          STANDARD GAUSSIAN DISTRIBUTION.
C      INIZ=-1         CAUCHY DISTRIBUTION.
C      INIZ=-2         UNIFORM DISTRIBUTION IN (-1,+1).
C      OTHERWISE       UNIFORM DISTRIBUTION IN ( 0,+1).
C
      DOUBLE PRECISION UNIFRD,PAI,A,B
C
      DATA PAI/3.14159265358979324D0/
C
      IF(INIZ.LE.0) GO TO 10
C
C INITIALIZATION.
C
      CHAOS = UNIFRD(INIZ)
      RETURN
C
  10  CONTINUE
      A = UNIFRD(0)
      IF(INIZ.NE.0) GO TO 20
      B = UNIFRD(0)
C
C GAUSSIAN RANDOM NUMBER BY POLAR METHOD
C
      CHAOS = DSQRT(-2.D0*DLOG(A))*DCOS(PAI*B)
C
      RETURN
  20  CONTINUE
C
C UNIFORM RANDOM NUMBER IN (0,1)
C
      CHAOS = A
C
C CAUCHY RANDOM NUMBER BY INVERSE TRANSFORMATION
C
      IF(INIZ.EQ.(-1)) CHAOS = DSIN(PAI*A)/DCOS(PAI*A)
C
C UNIFORM RANDOM NUMBER IN (-1,+1)
C
      IF(INIZ.EQ.(-2)) CHAOS = 2.D0*A-1.D0
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION UNIFRD(INIZ)
C
C UNIFRD GENERATES THE RANDOM NUMBERS UNIFORMLY DISTRIBUTED IN (0,1)
C EXPLOITING THOSE GENERATED BY  ALKNUT  WITH A FURTHER RANDOMIZATION
C IF THE INPUT PARAMETER INIZ  IS NOT 0
C THE RANDOM NUMBER GENERATOR IS INITIALIZED
C
      DOUBLE PRECISION A,B,C,X
      DOUBLE PRECISION X0,P,P0,P1,P2,R1,R2
      DOUBLE PRECISION FINV
C
      DIMENSION X(61)
C
      DATA NREM/61/,NRIP/100/
      DATA A,B,C/-1.5D0,5.5D0,-2.0D0/
      DATA FINV/3.5D-5/
      DATA IREM/0/
      DATA P0/3.D0/,P1,P2/1.D0,3.D0/,R1,R2/0.25D0,0.75D0/
C
      IF(INIZ.NE.0.OR.IREM.EQ.0) GO TO 10
C
      I0 = IREM
      X0 = X(I0)
C
C NONLINEARIZATION OF X0 TO AVOID LONG-DISTANCE LINEAR RELATIONSHIP.
C
      IF(X0.GE.FINV)X0 = DMOD(1.D0/X0,1.D0)
C
C UPDATE A COMPONENT OF THE VECTOR X ...
C
      CALL ALKNUT(NREM,X,IREM)
C
C ... AND FURTHER RANDOMIZE
C
      UNIFRD = DMOD(X0+X(I0),1.D0)
C
      RETURN
C
C INITIALIZATION OF THE RANDOM NUMBER GENERATOR
C
  10  CONTINUE
C
      P = P0-1.D0/DBLE(FLOAT(IABS(INIZ))+100.0)
      DO 20 K = 1,NREM
C
         P = C+P*(B+P*A)
C
C
         X(K) = R1+(R2-R1)*(P-P1)/(P2-P1)
  20     CONTINUE
C
      IREM = 0
C
      DO 30 K = 1,NRIP
C
         CALL ALKNUT(NREM,X,IREM)
C
  30     CONTINUE
C
      UNIFRD = X(1)
C
      RETURN
      END
      SUBROUTINE ALKNUT(NREM,X,IREM)
C
C UPDATES THE COMPONENT IREM OF THE NREM-VECTOR X WITH A RANDOM NUMBER
C UNIFORMLY DISTRIBUTED IN (0,1) BY MEANS OF THE ALGORITHM
C OF MITCHELL-MOORE, MODIFIED AS SUGGESTED BY BRENT, QUOTED IN
C D.E.KNUTH, THE ART OF COMPUTER PROGRAMMING, SECOND EDITION,
C SECOND VOLUME, SEMINUMERICAL ALGORITHMS, ADDISON-WESLEY
C PUB. CO., READING (1981), PP. 26-28.
C
      DOUBLE PRECISION X
C
      DIMENSION X(NREM)
C
      DATA N1,N2/24,55/
C
      IF(IREM.NE.0) GO TO 10
C
      IREM = NREM
      I1 = NREM-N1
      I2 = NREM-N2
C
  10  CONTINUE
C
      X(IREM) = DMOD(X(I1)+X(I2),1.D0)
C
      IREM = 1+MOD(IREM,NREM)
      I1 = 1+MOD(I1,NREM)
      I2 = 1+MOD(I2,NREM)
C
      RETURN
      END
      SUBROUTINE GAUSRV(N,W)
C
C GENERATES A RANDOM VECTOR SAMPLE FROM AN N-DIMENSIONAL
C NORMAL DISTRIBUTION
C
      DOUBLE PRECISION W,X,Y,R
      DOUBLE PRECISION CHAOS
C
      DIMENSION W(N)
C
      NN = (N+1)/2
C
      DO 20 I = 1,NN
         II = 1+N-I
C
  10     CONTINUE
         X = CHAOS(-2)
         Y = CHAOS(-2)
         R = X*X+Y*Y
C
         IF(R.GT.1.D0) GO TO 10
C
         R = DSQRT(-2.D0*DLOG(R)/R)
C
         W(I) = X*R
         W(II) = Y*R
C
  20     CONTINUE
C
      RETURN
      END
      SUBROUTINE UNITRV(N,W)
C
C GENERATES A RANDOM VECTOR UNIFORMLY DISTRIBUTED
C ON THE UNIT SPHERE.
C
      DOUBLE PRECISION W,WW
C
      DIMENSION W(N)
C
      CALL GAUSRV(N,W)
      WW = 0.D0
C
      DO 10 I = 1,N
         WW = WW+W(I)**2
  10     CONTINUE
      WW = 1.D0/DSQRT(WW)
C
      DO 20 I = 1,N
         W(I) = WW*W(I)
  20     CONTINUE
C
      RETURN
      END
      SUBROUTINE PTSEG(N,XPFMIN,FPFMIN,FPFMAX,KP,NFEV)
C
C  SAMPLE VERSION OF THE OUTPUT SUBROUTINE  PTSEG
C  (PERFORMS END-OF-SEGMENT OUTPUT)
C  WHICH MUST BY SUPPLIED BY THE USER AND WHICH IS DESCRIBED IN
C  DETAIL WITHIN THE SUBROUTINE  SIGMA.
C
      DOUBLE PRECISION XPFMIN,FPFMIN,FPFMAX
      DIMENSION XPFMIN(N)
C
      WRITE(6,100)KP,NFEV,FPFMIN,FPFMAX
 100  FORMAT(1X,24HOBSERVATION PERIOD  KP= ,I4,
     1    32H,   FUNCTION EVALUATIONS  NFEV= ,I7/
     2    2X,33HFINAL BEST, WORST FUNCT. VALUES  ,
     3    8HFPFMIN= ,G13.5,11H,  FPFMAX= ,G13.5)
      WRITE(6,200)XPFMIN
 200  FORMAT(2X,24HBEST FINAL POINT  XPFMIN/(1X,6G13.5))
C
C
      RETURN
      END
      SUBROUTINE PTRIAL ( N, XOPT, FOPT,
     1                    FTFMIN, FTFMAX, FTFOPT,
     2                    ISTOP, ISTOPT, NFEV, KP, IPRINT )
C
C  SAMPLE VERSION OF THE OUTPUT SUBROUTINE  PTRIAL
C  (PERFORMS END-OF-TRIAL OUTPUT)
C  WHICH MUST BY SUPPLIED BY THE USER AND WHICH IS DESCRIBED IN
C  DETAIL WITHIN THE SUBROUTINE  SIGMA.
C
      DOUBLE PRECISION XOPT,FTFMIN,FTFMAX,FTFOPT,FOPT
      DIMENSION XOPT(N)
C
      WRITE(6,100)
 100  FORMAT(//16H END OF A TRIAL )
      IF(IPRINT.EQ.0)WRITE(6,200)KP,NFEV,FTFMIN,FTFMAX
 200  FORMAT(2X,24HOBSERVATION PERIOD  KP= ,I4,
     1    32H,   FUNCTION EVALUATIONS  NFEV= ,I7/
     2    2X,33HFINAL BEST, WORST FUNCT. VALUES  ,
     3    8HFTFMIN= ,G13.5,11H,  FTFMAX= ,G13.5)
      WRITE(6,300)ISTOP,ISTOPT,FTFOPT,FOPT
 300  FORMAT(/2X,29HTRIAL STOP INDICATOR  ISTOP= ,I2,
     1    34H,   PAST STOPS INDICATOR  ISTOPT= ,I2/
     2    2X,43HEND-OF-TRIAL BEST FUNCTION VALUE   FTFOPT= ,G14.6/
     3    2X,43HBEST CURRENT MINIMUM FUNCTION VALUE  FOPT= ,G14.6)
      WRITE(6,400)XOPT
 400  FORMAT(2X,29HBEST CURRENT MINIMIZER   XOPT/(1X,6G13.5))
C
C
      RETURN
      END
      SUBROUTINE PTKSUC(KSUC)
C
C  SAMPLE VERSION OF THE OUTPUT SUBROUTINE  PTKSUC
C  (PERFORMS END-OF-TRIAL OUTPUT RELATED TO THE COUNT OF
C    SUCCESSFUL TRIALS)
C  WHICH MUST BY SUPPLIED BY THE USER AND WHICH IS DESCRIBED IN
C  DETAIL WITHIN THE SUBROUTINE  SIGMA.
C
      WRITE(6,100)KSUC,KSUC
 100  FORMAT(///1X,
     1    50HTHE CURRENT COUNT  KSUC  OF SUCCESSFUL TRIALS HAS ,
     2    21HREACHED FOR THE FIRST/2X,15HTIME THE VALUE ,I2,
     3    53H (IF THE REQUESTED COUNT  NSUC  OF SUCCESSFUL TRIALS /
     4    2X,25HHAD BEEN GIVEN THE VALUE ,I2,21H THE ALGORITHM WOULD ,
     5    18HHAVE STOPPED HERE)///)
C
C
      RETURN
      END
C
C (ALGORITHM SIGMA)
C MAIN PROGRAM (TEST VERSION)
C (CALL  SIGMA  VIA THE DRIVER  SIGMA1 )
C
      DOUBLE PRECISION FMIN,X0,XMAXGL,XMIN,XMINGL
C
C COMMON AREA TO PASS TEST-PROBLEM NUMBER  NPROB
C TO THE FUNCTION  FUNCT  WHICH WILL COMPUTE
C THE FUNCTION VALUES OF TEST-PROBLEM  NPROB
C BY CALLING THE TEST-PROBLEM COLLECTION SUBROUTINE  GLOMTF
C
      COMMON /TUN/ NPROB
C
C  X0     INITIAL POINT
C  XMIN   FINAL ESTIMATE OF GLOBAL MINIMUM
C  XMINGL, XMAXGL  MUST BE DIMENSIONED HERE IN ORDER TO CALL
C         THE PRE-EXISTING SUBROUTINE  GLOMIP.
      DIMENSION X0(100),XMIN(100),XMINGL(100),XMAXGL(100)
C
  10  CONTINUE
C
C INPUT PROBLEM NUMBER
      WRITE(6,20)
  20  FORMAT(//////41H INPUT PROBLEM NUMBER (1 TO 37, 0 = STOP))
      READ(5,30)NPROB
  30  FORMAT(I2)
      WRITE(6,40)NPROB
  40  FORMAT(///18H PROBLEM NUMBER = ,I2//////)
C
C TERMINATE OR CONTINUE
      IF(NPROB.EQ.0)GO TO 50
C
C CALL  GLOMIP  TO GET PROBLEM DIMENSION  N  AND INITIAL POINT  X0
C NOTE THAT  GLOMIP  RETURNS ALSO THE BOUNDARIES  XMINGL ,  XMAXGL
C OF THE OBSERVATION REGION (NOT NEEDED HERE)
      CALL GLOMIP(NPROB,N,X0,XMINGL,XMAXGL)
C
C SET  NSUC  SO AS TO HAVE GOOD CHANCES, WITHOUT PROHIBITIVE
C COMPUTATIONAL EFFORT
      NSUC = 5
C
C SET  IPRINT  SO AS TO HAVE A MODERATE OUTPUT
      IPRINT = 0
C
C CALL DRIVER SUBROUTINE  SIGMA1
      CALL SIGMA1(N,X0,NSUC,IPRINT,XMIN,FMIN,NFEV,IOUT)
C
C GO TO THE NEXT PROBLEM
      GO TO 10
C
C END OF TEST PROBLEMS
  50  CONTINUE
      WRITE(6,60)
  60  FORMAT(/22H END OF TEST PROBLEMS /)
C
      STOP
C
      END
      DOUBLE PRECISION FUNCTION FUNCT(N,X)
C
C COMPUTES THE FUNCTION VALUES OF TEST-PROBLEM  NPROB
C BY CALLING THE SUBROUTINE  GLOMTF .
C
      DOUBLE PRECISION X
      DOUBLE PRECISION F
      DIMENSION X(N)
C
      COMMON /TUN/ NPROB
C
      CALL GLOMTF(NPROB,N,X,F)
      FUNCT = F
C
      RETURN
      END
      SUBROUTINE SIGMA1(N,X0,NSUC,IPRINT,XMIN,FMIN,NFEV,IOUT)
C
C  SIGMA1  IS A 'DRIVER' SUBROUTINE WHICH SIMPLY CALLS THE PRINCIPAL
C  SUBROUTINE  SIGMA  AFTER HAVING ASSIGNED DEFAULT VALUES TO A NUMBER
C  OF INPUT PARAMETERS OF  SIGMA , AND HAS THEREFORE A CONSIDERABLY
C  LOWER NUMBER OF INPUT PARAMETERS.
C  IT CAN BE USED AS A SIMPLE EXAMPLE OF HOW TO CALL  SIGMA , BUT
C  ALSO AS AN EASY-TO-USE DRIVER FOR THE AVERAGE USER, WHICH MAY FIND
C  IT EASIER TO CALL  SIGMA1  INSTEAD OF  SIGMA , THUS AVOIDING THE
C  TROUBLE OF ASSIGNING A VALUE TO ALL THE INPUT PARAMETERS OS  SIGMA .
C  ALL THE PARAMETERS IN THE DEFINITION OF  SIGMA1  HAVE THE SAME MEA-
C  NING AS IN  SIGMA .
C
C  THE USER OF  SIGMA1  MUST ONLY GIVE VALUES TO THE INPUT PARAMETERS
C      N, X0, NSUC, IPRINT
C  AND OBTAINS ON OUTPUT THE SAME OUTPUT PARAMETERS OF  SIGMA
C      XMIN, FMIN, NFEV, IOUT
C
C  WE RECALL HERE THE MEANING OF THE ABOVE PARAMETERS
C
C  N       IS THE PROBLEM DIMENSION (NUMBER OF VARIABLES)
C  X0      IS AN  N-VECTOR CONTAINING THE INITIAL VALUES OF THE
C          X-VARIABLES
C  NSUC    IS THE REQUESTED NUMBER OF SUCCESSFUL TRIALS
C          AFTER WHICH THE COMPUTATION IS STOPPED.
C  IPRINT  IS AN INDEX USED TO CONTROL THE AMOUNT OF PRINTED OUTPUT
C          BY CONTROLLING THE CALLS TO THE USER-SUPPLIED OUTPUT SUB-
C          ROUTINES  PTSEG  (END-OF-SEGMENT OUTPUT),  PTRIAL  (END-
C          OF-TRIAL OUTPUT), AND  PTKSUC (END-OF-TRIAL OUTPUT RELATED
C          TO THE COUNT OF SUCCESSFUL TRIALS), WHICH ARE DESCRIBED
C          BELOW.
C          IPRINT.LT.0   NO CALL TO THE OUTPUT SUBROUTINES
C          IPRINT.EQ.0   CALL ONLY  PTRIAL  AND  PTKSUC
C          IPRINT.GT.0   CALL ALL OUTPUT SUBROUTINES.
C  XMIN    IS AN N-VECTOR CONTAINING THE COORDINATES OF THE POINT
C          (OR POSSIBLY ONE OF THE POINTS) WHERE THE FINAL VALUE  FMIN
C          OF  FOPT  WAS FOUND.
C  FMIN    IS THE FINAL VALUE OF THE BEST CURRENT MINIMUM FUNCTION
C          VALUE  FOPT.
C  NFEV    IS THE TOTAL NUMBER OF FUNCTION EVALUATION (INCLUDING
C          THOSE USED FOR THE COMPUTATION OF DERIVATIVES, AND FOR
C          THE REJECTED TIME-INTEGRATION STEPS).
C  IOUT    IS THE INDICATOR OF THE STOPPING CONDITIONS, AS FOLLOWS
C          IF  IOUT = -99  A FATAL ERROR WAS DETECTED WHEN PERFOR-
C             MING SOME PRELIMINARY CHECKING OF THE INPUT DATA, AND
C             THE ALGORITHM WAS NOT EVEN STARTED
C          OTHERWISE THE ALGORITHM WAS STARTED, AND THE MEANING OF
C          IOUT  IS AS FOLLOWS
C          IF  IOUT.EQ.0  NO TRIAL HAD A UNIFORM STOP.
C          IF  IOUT.NE.0  THE SIGN AND THE MAGNITUDE OF  IOUT  HAVE
C          THE FOLLOWING MEANING
C          A) THE SIGN OF  IOUT  INDICATES IF THE BEST UNIFORM STOP
C          (I.E. THE ONE IN WHICH  FTFOPT  WAS OBTAINED)
C          IS TO BE CONSIDERED SUCCESSFUL (IOUT.GT.0) OR UNSUCCESSFUL
C          B) THE MAGNITUDE OF  IOUT  INDICATES WHICH ONE OF THE
C          NUMERICAL EQUALITY CRITERIA WAS SATISFIED IN THE BEST
C          UNIFORM STOP
C          IABS(IOUT).EQ.1  RELATIVE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.2  ABSOLUTE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.3  BOTH CRITERIA
C          (IOUT  IS THE FINAL VALUE OF THE INTERNAL PARAMETER  ISTOPT
C          (AN OUTPUT INDICATOR OF THE USER-SUPPLIED SUBROUTINE
C          PTRIAL )).
C          SUCCESS IS CLAIMED BY THE ALGORITHM IF  IOUT .GT. 0,
C          I.E. IF AT LEAST ONE OF THE TRIALS STOPPED SUCCESSFULLY
C          (I.E. WITH A POSITIVE VALUE OF THE TRIAL STOP INDICATOR
C          ISTOP ), AND THE SUCCESSFUL TRIALS (AS COUNTED BY  KSUC )
C          ARE CURRENTLY VALID.
C
      DOUBLE PRECISION X0,XMIN,FMIN
      DOUBLE PRECISION DX,EPS,H,TOLABS,TOLREL
      DOUBLE PRECISION VRMAX,VRMIN,XRMAX,XRMIN
      DIMENSION X0(N),XMIN(N)
      DIMENSION XRMIN(100),XRMAX(100)
      DATA VRMIN,VRMAX /-1.D4,1.D4/
      DATA NTRILD/50/
C
      H = 1.D-10
      EPS = 1.D0
      DX = 1.D-9
      IRAND = 0
      NTRAJ = 0
      ISEGBR = 0
      INKPBR = 0
      KPBR0 = 0
      NPMIN = 10
      NPMAX0 = 100
      INPMAX = 50
      NTRIAL = MAX0(NTRILD,5*NSUC)
      TOLREL = 1.D-3
      TOLABS = 1.D-6
      KPASCA = 10
      IF(N.GT.5)KPASCA = 300
      INHP = 1
      DO 1 IX = 1,N
         XRMIN(IX)=VRMIN
         XRMAX(IX)=VRMAX
 1    CONTINUE
C
      CALL SIGMA ( N, X0, H, EPS, DX,
     1             NTRAJ, ISEGBR, KPBR0, INKPBR,
     2             NPMIN, NPMAX0, INPMAX,
     3             NSUC, NTRIAL, TOLREL, TOLABS, XRMIN, XRMAX,
     4             KPASCA, IRAND, INHP, IPRINT,
     5             XMIN, FMIN, NFEV, IOUT )
C
      RETURN
      END
      SUBROUTINE SIGMA ( N, X0, H, EPS, DX,
     1                   NTRAJ, ISEGBR, KPBR0, INKPBR,
     2                   NPMIN, NPMAX0, INPMAX,
     3                   NSUC, NTRIAL, TOLREL, TOLABS, XRMIN, XRMAX,
     4                   KPASCA, IRAND, INHP, IPRINT,
     5                   XMIN, FMIN, NFEV, IOUT )
C
C  THE SUBROUTINE  SIGMA  IS THE PRINCIPAL SUBROUTINE OF THE PACKAGE
C  SIGMA, WHICH ATTEMPTS TO FIND A GLOBAL MINIMIZER OF A REAL VALUED
C  FUNCTION  F(X) = F(X1,...,XN)  OF  N  REAL VARIABLES  X1,...,XN.
C  THE ALGORITHM AND THE PACKAGE ARE DESCRIBED IN DETAIL IN THE TWO
C  PAPERS PUBLISHED IN THE SAME ISSUE OF THE A.C.M. TRANSACTIONS ON
C  MATHEMATICAL SOFTWARE, BOTH BY
C     F. ALUFFI-PENTINI, V. PARISI, F. ZIRILLI.
C  (1)  A GLOBAL MINIMIZATION ALGORITHM USING STOCHASTIC DIFFERENTIAL
C       EQUATIONS
C  (2)  ALGORITHM SIGMA, A STOCHASTIC-INTEGRATION GLOBAL MINIMIZATION
C       ALGORITHM.
C  THE SOFTWARE IMPLEMENTATION AND ITS USAGE ARE DESCRIBED IN (2).
C
C  METHOD
C
C  A GLOBAL MINIMIZER OF  F(X)  IS SOUGHT BY MONITORING THE VALUES OF
C  F  ALONG TRAJECTORIES GENERATED BY A SUITABLE (STOCHASTIC) DISCRE-
C  TIZATION OF A FIRST-ORDER STOCHASTIC DIFFERENTIAL EQUATION INSPIRED
C  BY STATISTICAL MECHANICS. STARTING FROM AN INITIAL POINT  X0 ,
C  X  IS UPDATED BY THE (STOCHASTIC) DISCRETIZATION STEP
C       X = X + DX1 + DX2
C  WHERE  DX1 = - H * GAM           (FIRST HALF-STEP)
C         DX2 = EPS * SQRT(H) * U   (SECOND HALF-STEP)
C  AND  H  IS THE TIME-INTEGRATION STEPLENGTH,
C  GAM/N  IS COMPUTED AS A FINITE-DIFFERENCE APPROXIMATION TO THE
C     DIRECTIONAL DERIVATIVE OF  F  ALONG AN ISOTROPICALLY RANDOM
C     DIRECTION,
C  EPS  IS A POSITIVE 'NOISE' COEFFICIENT, AND
C  U  IS A RANDOM SAMPLE FROM AN N-DIMENSIONAL GAUSSIAN DISTRIBUTION.
C  WE CONSIDER THE SIMULTANEOUS EVOLUTION OF A GIVEN FIXED NUMBER
C  NTRAJ  OF TRAJECTORIES DURING AN OBSERVATION PERIOD IN WHICH FOR
C  EACH TRAJECTORY  EPS  IS FIXED WHILE  H  AND THE SPATIAL DISCRETI-
C  ZATION INCREMENT  DX  FOR COMPUTING  GAM  ARE AUTOMATICALLY
C  ADJUSTED BY THE ALGORITHM.
C  AFTER EVERY OBSERVATION PERIOD ONE OF THE TRAJECTORIES IS DISCARDED,
C  ALL OTHER TRAJECTORIES CONTINUE UNPERTURBED, AND ONE OF THEM IS SE-
C  LECTED FOR BRANCHING, I.E. GENERATING A SECOND PERTURBED CONTI-
C  NUATION, WITH DIFFERENT STARTING  EPS  AND  DX  (AND THE SAME
C  'PAST HISTORY' OF THE FIRST).
C  THE SET OF SIMULTANEOUS TRAJECTORIES IS CONSIDERED A SINGLE TRIAL,
C  AND THE COMPLETE ALGORITHM IS A SET OF REPEATED TRIALS.
C  A TRIAL IS STOPPED, AT THE END OF AN OBSERVATION PERIOD, AND AFTER
C  HAVING DISCARDED THE WORST TRAJECTORY, IF ALL THE FINAL VALUES OF
C  F  FOR THE REMAINING TRAJECTORIES ARE EQUAL (WITHIN NUMERICAL TO-
C  LERANCES, AND POSSIBLY AT DIFFERENT POINTS  X) TO THEIR MINIMUM
C  VALUE  FTFMIN ('UNIFORM STOP AT THE LEVEL  FTFMIN ').
C  A UNIFORM STOP IS CONSIDERED SUCCESSFUL ONLY IF THE FINAL VALUE
C  FTFMIN  IS (NUMERICALLY) EQUAL TO THE CURRENT BEST MINIMUM  FOPT
C  FOUND SO FAR FROM ALGORITHM START.
C  A TRIAL IS ALSO ANYWAY STOPPED (UNSUCCESSFULLY) IF A GIVEN MAXIMUM
C  NUMBER  NPMAX  OF OBSERVATION PERIODS HAS ELAPSED.
C  TRIALS ARE REPEATED WITH DIFFERENT OPERATING CONDITIONS (INITIAL
C  POINT, MAX TRIAL LENGTH  NPMAX , SEED OF NOISE GENERATOR, POLICY
C  FOR CHOOSING THE STARTING  EPS  FOR THE PERTURBED CONTINUATION,
C  AND TRIAL-START VALUE OF  EPS ).
C  THE OUTCOMES OF ALL THE COMPLETED TRIALS ARE SUMMARIZED BY
C  FTFOPT (THE BEST CURRENT MIMNIMUM VALUE FOUND SO FAR FOR
C  FTFMIN ) AND BY THE CURRENT COUNT  KSUC  OF SUCCESSFUL TRIALS
C  (WHICH IS ZERO AT ALGORITHM START, AND IS UPDATED AT EVERY
C  SUCCESSFULLY-STOPPING TRIALS).
C  THE ALGORITHM IS STOPPED, AT THE END OF A TRIAL, IF A REQUESTED
C  COUNT  NSUC  OF SUCCESSFUL TRIALS HAS BEEN REACHED,
C  OR ANYWAY IF A GIVEN MAXIMUM NUMBER  NTRIAL  OF TRIALS HAS BEEN
C  REACHED.
C  SUCCESS IS CLAIMED IF THE CURRENT COUNT  KSUC  OF SUCCESSFUL TRIALS
C  IS AT LEAST ONE, AND IF SUCH TRIALS ARE CURRENTLY VALID.
C
C  CALL STATEMENT
C
C  THE CALL STATEMENT IS
C     CALL SIGMA ( N, X0, H, EPS, DX,
C                  NTRAJ, ISEGBR, KPBR0, INKPBR,
C                  NPMIN, NPMAX0, INPMAX,
C                  NSUC, NTRIAL, TOLREL, TOLABS, XRMIN, XRMAX,
C                  KPASCA, IRAND, INHP, IPRINT,
C                  XMIN, FMIN, NFEV, IOUT )
C
C  CALL PARAMETERS
C
C  INPUT PARAMETERS ARE THOSE IN LINES 1,3,4,5 OF THE CALL STATEMENT,
C  INPUT-OUTPUT PARAMETERS ARE THOSE IN LINE 2,
C  OUTPUT PARAMETERS ARE THOSE IN LINE 6.
C  NOTE THAT A NUMBER OF OTHER (INTERNAL) PARAMETERS CAN BE OBTAINED
C  BY MEANS OF THE USER-SUPPLIED OUTPUT SUBROUTINES  PTSEG, PTRIAL,
C  AND  PTKSUC, WHICH ARE DESCRIBED BELOW.
C
C  DESCRIPTION OF THE CALL PARAMETERS
C
C  N        IS THE PROBLEM DIMENSION (NUMBER OF VARIABLES)
C  X0       IS AN  N-VECTOR CONTAINING THE INITIAL VALUES OF THE
C          X-VARIABLES
C  H       IS THE INITIAL VALUE OF THE TIME-INTEGRATION STEPLENGTH.
C  EPS     IS THE INITIAL VALUE OF THE NOISE COEFFICIENT.
C  DX      IS THE INITIAL VALUE OF THE MAGNITUDE OF THE DISCRETIZATION
C          INCREMENT FOR COMPUTING THE FINITE-DIFFERENCE DERIVATIVES.
C  NTRAJ   IS THE NUMBER OF SIMULTANEOUS TRAJECTORIES.
C          (NOTE HOWEVER THAT IF THE INPUT VALUE IS ZERO,  NTRAJ  IS
C          SET TO A DEFAULT VALUE (NTRAJ = 7), AND IF THE INPUT VALUE
C          IS OTHERWISE OUTSIDE THE INTERVAL (3,20) NTRAJ IS SET TO
C          THE NEAREST EXTREME VALUE).
C  ISEGBR, KPBR0, INKPBR  DETERMINE, AT THE END OF AN OBSERVATION
C          PERIOD, WHICH ONE OF THE SIMULTANEOUS TRAJECTORIES
C          IS TO BE BRANCHED, AS FOLLOWS.
C          BRANCHING IS NORMALLY PERFORMED ON THE TRAJECTORY WHICH
C          OCCUPIES THE PLACE  ISEGBR  IN THE TRAJECTORY SELECTION OR-
C          DERING, EXCEPT AT (THE END OF) EXCEPTIONAL OBSERVATION
C          PERIODS, WHERE THE FIRST TRAJECTORY IN THE ORDERING IS
C          BRANCHED. EXCEPTIONAL BRANCHING OCCURS AT THE OBSERVATION
C          PERIODS NUMBERED  KP = KPBR0 + J*INKPBR, (J = 1,2,3,...).
C          THEREFORE  ISEGBR  SELECTS THE LEVEL (IN THE ORDERING) AT
C          WHICH NORMAL BRANCHING OCCURS, WHILE  KPBR0  AND  INKPBR
C          SELECT THE FIRST OCCURRENCE AND THE REPETITION FREQUENCY
C          OF THE EXCEPTIONAL OBSERVATION PERIODS.
C          (NOTE HOWEVER THAT IF ONE OF THE INPUT VALUES IS ZERO,
C          THE CORRESPONDING VARIABLE IS SET TO A DEFAULT VALUE
C          ISEGBR = INT((1+NTRAJ)/2), INKPBR = 10, KPBR0 = 3.
C          IF THE INPUT VALUE FOR  ISEGBR  IS OTHERWISE OUTSIDE THE
C          INTERVAL  (1,NTRAJ),  ISEGBR  IS SET TO THE NEAREST
C          EXTREME VALUE, AND IF KPBR0  HAS A VALUE NOT INSIDE THE
C          INTERVAL  (1,INKPBR), IT IS ASSIGNED THE SAME VALUE
C          MODULO  INKPBR).
C  NPMIN   IS THE MINIMUM DURATION OF A TRIAL, I.E. THE MINIMUM
C          NUMBER OF OBSERVATION PERIODS THAT SHOULD ELAPSE BEFORE
C          STARTING TO CHECK THE TRIAL STOPPING CRITERIA.
C  NPMAX0  IS THE MAXIMUM DURATION OF THE FIRST TRIAL, I.E. THE
C          VALUE, FOR THE FIRST TRIAL, OF MAXIMUM ACCEPTABLE
C          NUMBER  NPMAX  OF OBSERVATION PERIODS IN A TRIAL.
C  INPMAX  IS THE INCREMENT FOR  NPMAX , WHEN  NPMAX  IS VARIED
C          FROM ONE TRIAL TO THE FOLLOWING ONE.
C  NSUC    IS THE REQUESTED NUMBER OF SUCCESSFUL TRIALS
C          AFTER WHICH THE COMPUTATION IS STOPPED.
C  TOLREL  AND  TOLABS  ARE THE RELATIVE AND ABSOLUTE TOLERANCES
C          FOR STOPPING A SINGLE TRIAL.
C  XRMIN,  XRMAX  ARE N-VECTORS DEFINING THE ADMISIBLE REGION FOR
C          THE X-VALUES, WITHIN WHICH THE FUNCTION VALUES CAN BE
C          SAFELY COMPUTED.
C  KPASCA  IS THE MINIMUM NUMBER OF TRAJECTORY SEGMENTS THAT SHOULD
C          ELAPSE BEFORE THE RESCALING PROCEDURES ARE ACTIVATED.
C  IRAND   IS A CONTROL INDEX FOR THE INITIALIZATION OF THE RANDOM
C          NUMBER GENERATOR.
C          IRAND.GT.0    THE GENERATOR IS INITIALIZED, BEFORE STAR-
C                        TING THE TRIAL  KT, WITH SEED  IRAND+KT-1
C          IRAND.LE.0    THE GENERATOR IS INITIALIZED (WITH SEED 0)
C                        ONLY AT THE FIRST CALL OF SIGMA
C  INHP    IS A CONTROL INDEX FOR SELECTING THE NUMBER  NHP  OF TIME-
C          INTEGRATION STEPS FOR OBSERVATION PERIOD  KP (DURATION OF
C          TRIAL  KP) AS FOLLOWS (LOG IS BASE 2)
C          INHP=1  NHP = 1 + INT(LOG(KP))      ('SHORT' DURATION)
C          INHP=2  NHP = INT(SQRT(KP))         ('MEDIUM' DURATION)
C          INHP=3  NHP = KP                    ('LONG' DURATION)
C  IPRINT  IS AN INDEX USED TO CONTROL THE AMOUNT OF PRINTED OUTPUT
C          BY CONTROLLING THE CALLS TO THE USER-SUPPLIED OUTPUT SUB-
C          ROUTINES  PTSEG  (END-OF-SEGMENT OUTPUT),  PTRIAL  (END-
C          OF-TRIAL OUTPUT), AND  PTKSUC (END-OF-TRIAL OUTPUT RELATED
C          TO THE COUNT OF SUCCESSFUL TRIALS), WHICH ARE DESCRIBED
C          BELOW.
C          IPRINT.LT.0   NO CALL TO THE OUTPUT SUBROUTINES
C          IPRINT.EQ.0   CALL ONLY  PTRIAL  AND  PTKSUC
C          IPRINT.GT.0   CALL ALL OUTPUT SUBROUTINES.
C  XMIN    IS AN N-VECTOR CONTAINING THE COORDINATES OF THE POINT
C          (OR POSSIBLY ONE OF THE POINTS) WHERE THE FINAL VALUE  FMIN
C          OF  FOPT  WAS FOUND.
C  FMIN    IS THE FINAL VALUE OF THE BEST CURRENT MINIMUM FUNCTION
C          VALUE  FOPT.
C  NFEV    IS THE TOTAL NUMBER OF FUNCTION EVALUATION (INCLUDING
C          THOSE USED FOR THE COMPUTATION OF DERIVATIVES, AND FOR
C          THE REJECTED TIME-INTEGRATION STEPS).
C  IOUT    IS THE INDICATOR OF THE STOPPING CONDITIONS, AS FOLLOWS
C          IF  IOUT = -99  A FATAL ERROR WAS DETECTED WHEN PERFOR-
C             MING SOME PRELIMINARY CHECKING OF THE INPUT DATA, AND
C             THE ALGORITHM WAS NOT EVEN STARTED
C          OTHERWISE THE ALGORITHM WAS STARTED, AND THE MEANING OF
C          IOUT  IS AS FOLLOWS
C          IF  IOUT.EQ.0  NO TRIAL HAD A UNIFORM STOP.
C          IF  IOUT.NE.0  THE SIGN AND THE MAGNITUDE OF  IOUT  HAVE
C          THE FOLLOWING MEANING
C          A) THE SIGN OF  IOUT  INDICATES IF THE BEST UNIFORM STOP
C          (I.E. THE ONE IN WHICH  FTFOPT  WAS OBTAINED)
C          IS TO BE CONSIDERED SUCCESSFUL (IOUT.GT.0) OR UNSUCCESSFUL
C          B) THE MAGNITUDE OF  IOUT  INDICATES WHICH ONE OF THE
C          NUMERICAL EQUALITY CRITERIA WAS SATISFIED IN THE BEST
C          UNIFORM STOP
C          IABS(IOUT).EQ.1  RELATIVE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.2  ABSOLUTE DIFFERENCE CRITERION
C          IABS(IOUT).EQ.3  BOTH CRITERIA
C          (IOUT  IS THE FINAL VALUE OF THE INTERNAL PARAMETER  ISTOPT
C          (AN OUTPUT INDICATOR OF THE USER-SUPPLIED SUBROUTINE
C          PTRIAL )).
C          SUCCESS IS CLAIMED BY THE ALGORITHM IF  IOUT .GT. 0,
C          I.E. IF AT LEAST ONE OF THE TRIALS STOPPED SUCCESSFULLY
C          (I.E. WITH A POSITIVE VALUE OF THE TRIAL STOP INDICATOR
C          ISTOP ), AND THE SUCCESSFUL TRIALS (AS COUNTED BY  KSUC )
C          ARE CURRENTLY VALID.
C
C
C  USER-SUPPLIED SUPROGRAMS
C
C  THE USER MUST PROVIDE THE FUNCTION  FUNCT  TO COMPUTE  F(X),
C  AND THE THREE OUTPUT SUBROUTINE  PTSEG, PTRIAL, PTKSUC .
C  THE CALLS TO THE OUTPUT SUBROUTINES ARE CONTROLLED BY  IPRINT
C  (INPUT PARAMETER TO  SIGMA).
C  A USER NOT INTERESTED IN USING ANY ONE OF THE OUTPUT SUBROUTINES
C  CAN SIMPLY PUT  IPRINT = -1.
C  SAMPLE VERSIONS OF THE OUTPUT SUBROUTINES ARE PROVIDED
C  WITH THE PACKAGE FOR THE BENEFIT OF THE AVERAGE USER.
C  IN THE FOLLOWING DESCRIPTION ALL NON-INTEGER ARGUMENTS ARE
C  DOUBLE PRECISION (INTEGER ARGUMENTS ARE INDICATED BY MEANS OF THE
C  FORTRAN IMPLICIT TYPE DEFINITION CONVENTION).
C
C     THE FUNCTION  FUNCT
C
C     FUNCT  MUST RETURN AS ITS VALUE THE VALUE AT  X  OF THE FUNCTION
C     TO BE MINIMIZED
C     THE DEFINITION STATEMENT IS
C        DOUBLE PRECISION FUNCTION  FUNCT (N, X)
C     WHERE
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM.
C     X   IS THE (INPUT) N-VECTOR CONTAINING THE COORDINATES OF THE
C         POINT  X  WHERE THE FUNCTION IS TO BE COMPUTED.
C
C     THE SUBROUTINE  PTSEG
C
C     PTSEG  IS CALLED (IF  IPRINT.GT.0) AT THE END OF EVERY OBSER-
C     VATION PERIOD.
C     THE DEFINITION STATEMENT IS
C         SUBROUTINE  PTSEG ( N, XPFMIN, FPFMIN, FPFMAX,
C                                         KP, NFEV, IPRINT )
C     WHERE
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM
C     FPFMIN, FPFMAX  ARE RESPECTIVELY THE MINIMUM AND THE MAXIMUM
C        AMONG THE VALUES OF  F(X) OBTAINED AT THE FINAL POINTS OF
C        THE TRAJECTORY SEGMENTS OF THE (JUST ELAPSED) OBSERVATION
C        PERIOD  KP.
C     XPFMIN   IS AN N-VECTOR CONTAINING THE COORDINATES OF THE
C        (FINAL) POINT (OR POSSIBLY ONE OF THE POINTS) WHERE
C        FPFMIN WAS OBTAINED.
C     KP   IS THE TOTAL NUMBER OF ELAPSED OBSERVATION PERIODS IN
C        THE CURRENT TRIAL.
C     NFEV   IS THE TOTAL NUMBER OF FUNCTION EVALUATIONS PERFORMED
C        FROM ALGORITHM START.
C
C     THE SUBROUTINE  PTRIAL
C
C     PTRIAL  IS CALLED (IF IPRINT.GE.0) AT THE END OF EVERY TRIAL.
C     THE DEFINITION STATEMENT IS
C         SUBROUTINE PTRIAL ( N, XOPT, FOPT,
C                             FTFMIN, FTFMAX, FTFOPT,
C                             ISTOP, ISTOPT, NFEV, KP, IPRINT )
C     WHERE
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM.
C     XOPT   IS AN N-VECTOR CONTINING THE COORDINATES OF THE
C        POINT (OR POSSIBLY ONE OF THE POINTS) WHERE THE CURRENT
C        MINIMUM  FOPT  WAS OBTAINED.
C     FOPT   IS THE CURRENT BEST MINIMUM VALUE FOUND FOR  F  FROM
C        ALGORITHM START ( FOPT  IS UPDATED WHENEVER A FUNCTION
C        VALUE IS COMPUTED).
C     FTFMIN, FTFMAX   ARE RESPECTIVELY THE MINIMUM AND THE MAXIMUM
C        AMONG THE VALUES OF  F(X)  OBTAINED AT THE FINAL POINTS OF
C        THE LAST TRAJECTORY SEGMENTS OF THE CURRENT TRIAL.
C     FTFOPT   IS THE CURRENT MINIMUM VALUE OF  FTFMIN  AMONG THE
C        TRIALS WHICH DID NOT STOP FOR REACHING THE MAXIMUM ALLOWED
C        NUMBER OF SEGMENTS (STOPPING INDICATOR  ISTOP = 0, SEE
C        BELOW).  FTFOPT  IS USED BY  SIGMA  TO COMPUTE  KSUC  (INPUT
C        PARAMETER TO THE OUTPUT SUBROUTINE  PTKSUK , SEE BELOW).
C     KP   IS THE TOTAL NUMBER OF ELAPSED OBSERVATION PERIODS IN
C        THE CURRENT TRIAL.
C     NFEV   IS THE TOTAL NUMBER OF FUNCTION EVALUATIONS PERFORMED
C        FROM ALGORITHM START.
C     ISTOP   IS THE INDICATOR OF THE STOPPING CONDITION OF THE TRIAL,
C        AS FOLLOWS
C        ISTOP = 0
C           THE MAXIMUM NUMBER  NPMAX  OF OBSERVATION PERIODS HAS
C           BEEN REACHED.
C        ISTOP.NE.0
C           ALL THE END-OF-SEGMENT VALUES OF  F(X) , (EXCEPT FOR THE
C           JUST DISCARDED SEGMENT) ARE CLOSE ENOUGH TO THEIR COMMON
C           MINIMUM VALUE  FPFMIN , WITH RESPECT TO AN ABSOLUTE OR
C           RELATIVE DIFFERENCE CRITERION, TO BE CONSIDERED NUMERI-
C           CALLY EQUAL.
C           THE ABSOLUTE VALUE AND THE SIGN OF  ISTOP  HAVE THE
C           FOLLOWING MEANING.
C              THE ABSOLUTE VALUE INDICATES WHICH DIFFERENCE
C                 CRITERION WAS SATISFIED
C                 1   RELATIVE DIFFERENCE CRITERION SATISFIED
C                 2   ABSOLUTE DIFFERENCE CRITERION SATISFIED
C                 3   BOTH CRITERIA SATISFIED
C              THE SIGN OF  ISTOP  INDICATES THE RELATIONSHIP BETWEEN
C                 THE END-OF-TRIAL VALUE  FPFMIN  AND THE CURRENT
C                 BEST MINIMUM VALUE  FOPT  (WHICH IS UPDATED WHEN-
C                 EVER A FUNCTION  VALUE IS COMPUTED
C                 ISTOP.GT.0
C                    FPFMIN  IS NUMERICALLY EQUAL (W.R.T. AT LEAST
C                    ONE OF THE ABOVE  DIFFERENCE CRITERIA) TO  FOPT
C                 ISTOP.LT.0
C                    FPFMIN  IS NOT EVEN NUMERICALLY EQUAL TO  FOPT
C                    (AND THEREFORE CANNOT BE CONSIDERED AS AN
C                    ACCEPTABLE GLOBAL MINIMUM).
C     ISTOPT   IS THE VALUE OF THE TRIAL STOPPING INDICATOR  ISTOP
C        CORRESPONDING TO THE (CURRENT OR PAST) TRIAL WHERE  FTFOPT
C        WAS OBTAINED, WITH THE SIGN WHICH IS UPDATED ACCORDING TO
C        THE COMPARISON BETWEEN  FTFOPT  AND THE PRESENT VALUE OF
C        FOPT , AS DESCRIBED ABOVE.
C        THE FINAL VALUE OF  ISTOP IS RETURNED BY  SIGMA  AS THE VALUE
C        OF THE OUTPUT INDICATOR IOUT OF THE ALGORITHM STOPPING CON-
C        DITIONS (WHENEVER THE ALGORITHM WAS STARTED, IOUT.NE.-99,
C        SEE ABOVE).
C
C     THE SUBROUTINE  PTKSUC
C
C     PTKSUC  IS CALLED ONLY AT THE END OF EVERY SUCCESSFUL TRIAL
C     SUCH THAT AN INCREMENT OCCURRED IN THE VALUE OF THE
C     CURRENT MAXIMUM  MSUC  AMONG ALL THE VALUES OF  KSUC  FROM
C     ALGORITHM START. A CALL TO  PTKSUC  THEREFORE PROVIDES
C     THE USER WITH THE OPERATIONALLY INTERESTING INFORMATION THAT
C     ALGORITHM STOP WOULD HAVE TAKEN PLACE, IF  NSUC  (INPUT
C     PARAMETER TO  SIGMA ) HAD BEEN GIVEN A LOWER VALUE, EQUAL TO
C     THE CURRENT  KSUC .
C     PTKSUC  IS CALLED ONLY IF  IPRINT.GE.0  AND  KSUC.LT.NSUC  .
C     THE DEFINITION STATEMENT IS
C         SUBROUTINE  PTKSUC ( KSUC )
C     WHERE  KSUC  IS THE CURRENT COUNT OF SUCCESSFUL TRIALS
C     (1 .LE. KSUC .LE. NSUC)
C
      DOUBLE PRECISION X0,H,EPS,DX,TOLREL,TOLABS
      DOUBLE PRECISION XRMIN,XRMAX,XMIN,FMIN
      DOUBLE PRECISION EPSAG,EPSAP,EPSC
      DOUBLE PRECISION EPSMAX,EPSR,F,FOPT,FTFMAX
      DOUBLE PRECISION FTFMIN,FTFOPT
C
      DOUBLE PRECISION X,HC,DXC,VMVT,EPSCO,VMCOR,VCOR
      DOUBLE PRECISION XRMIC,XRMAC,XOPT,FOPTC
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
C
      DIMENSION X0(N),XMIN(N),XRMIN(N),XRMAX(N)
C
      COMMON /DINCOM/ X(100,20),HC(20),DXC(20),VMVT(20,19),EPSCO(20),
     1    VMCOR(20),VCOR(20),XRMIC(100),XRMAC(100),XOPT(100),FOPTC,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJC,NTRAJR,
     3    ISEGBC,INKPBC,KPBR0C,NCF,IFEPC,INHPC
C
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
C  DATA FOR THE VARIATION OF NOISE COEFFICIENT
C
      DATA EPSR/1.D4/,EPSAP/10.D0/,EPSAG/1.D3/
      DATA EPSMAX/1.D15/
C
      IFEP = 1
C
C  INITIALIZE COMMON AREA  /DINCOM/
C
      CALL INIT(N,X0,H,EPS,DX,IRAND,F,
     1    NTRAJ,ISEGBR,INKPBR,KPBR0,INHP,IFEP,XRMIN,XRMAX,IOUT)
C
C  CHECK PARAMETER VALUES
C
      IF(NPMIN.LE.0.OR.NPMAX0.LT.0.OR.INPMAX.LE.0.OR.
     1    NSUC.LE.0.OR.NTRIAL.LE.0)IOUT = -99
      IF (IOUT.EQ.(-99))RETURN
C
C  INITIALIZE VARIABLES
C
      EPSC = EPS
      NPMAX = NPMAX0
      ISTOPT = 0
      ISTOP = 0
      ICCOM = (NTRIAL+NTRIAL+4)/5
      NFEV = 0
      NTES = NSUC
      ICTS = NSUC-1
      FTFOPT = F
C
C  START SERIES OF TRIALS
C
      DO 30 IC = 1,NTRIAL
C
C  SET INITIALIZATION INDEX FOR NOISE GENERATOR
C
         IS = 0
         IF(IRAND.GT.0)IS = IRAND+IC-1
C
C  INITIALIZE TRIAL
C
         IF(IC.GT.1.AND.IC.LE.ICCOM)CALL REINIT(N,X0,EPSC,IS,F,IFEP)
         IF(IC.GT.ICCOM)CALL REINIT(N,XMIN,EPSC,IS,FOPT,IFEP)
C
         FTFMIN = F
         FTFMAX = F
         NFEV = NFEV+1
C
C  PRINT INITIAL CONDITIONS OF TRIAL
C
         IF(IPRINT.GT.0)CALL PTSEG(N,X0,FTFMIN,FTFMAX,0,NFEV)
C
C
C  DEACTIVATE SCALING
C
         IF(KPASCA.GT.NPMAX.OR.N.LE.1)CALL NOSCA
C
C  INITIALIZE COMMON AREA  /SCALE/
C
         IF(KPASCA.LE.NPMAX.AND.N.GT.1)CALL INISCA(N,NTRAJ)
C
C  PERFORM A TRIAL
C
         CALL TRIAL(N,NPMIN,NPMAX,KPASCA,TOLREL,TOLABS,
     1       IPRINT,XMIN,FTFMIN,
     1       FTFMAX,NFEV,KP,ISTOP)
C
C
C     EVALUATE PAST TRIAL AND PREPARE NEXT TRIAL
C
C
C  SET TRIAL DURATION
C
         IF(ISTOP.EQ.0)NPMAX = NPMAX+INPMAX
C
C  RESET CURRENT NUMBER OF SUCCESSES REQUIRED BEFORE STOPPING
C  COMPUTE INDICATOR OF TRIAL STOPPING CONDITIONS
C  UPDATE BEST CURRENT VALUES OF TRIAL STOPPING INDICATOR AND
C  OF FUNCTION  F(X)  AT TRIAL STOP
C
         IF((FTFMIN.GT.FTFOPT.OR.ISTOP.EQ.0).AND.ISTOPT.NE.0)GO TO 10
         IF(ITOLCH(FTFOPT,FTFMIN,TOLREL,TOLABS).EQ.0)NTES = NSUC
C
         FTFOPT = FTFMIN
         ISTOPT = ISTOP
  10     CONTINUE
         ISTOPT = IABS(ISTOPT)
         CALL RCLOPT(N,XMIN,FOPT)
         IF(ITOLCH(FTFOPT,FOPT,TOLREL,TOLABS).EQ.0)ISTOPT = -ISTOPT
         IF(ITOLCH(FTFMIN,FOPT,TOLREL,TOLABS).EQ.0)ISTOP = -ISTOP
C
C  END-OF-TRIAL PRINT OUT
C
         IF(IPRINT.GE.0)
     1      CALL PTRIAL ( N, XOPT, FOPT,
     2                    FTFMIN, FTFMAX, FTFOPT,
     3                    ISTOP, ISTOPT, NFEV, KP, IPRINT )
C
C  UPDATE INITIAL VALUE OF NOISE COEFFICIENT FOR NEXT TRIAL
C
         IF (ISTOP.EQ.0)EPSC = EPSC/EPSR
         IF(ISTOP.GT.0)EPSC = EPSC*EPSAG
         IF(ISTOP.LT.0.AND.IC.LE.ICCOM)EPSC = EPSC*EPSAP
         IF(ISTOP.LT.0.AND.IC.GT.ICCOM)EPSC = EPSC/EPSR
C
C  UPDATE OPERATING CONDITIONS FOR SELECTING (IN THE NEXT TRIAL)
C  THE STARTING VALUE OF THE NOISE COEFFICIENT OF THE NEW TRAJECTORY
C  AFTER BRANCHING
C
         IF (ISTOP.EQ.0)IFEP = 1
         IF (ISTOP.NE.0)IFEP = 2
C
C  UPDATE, PRINT, AND CHECK TRIAL STOPPING CONDITIONS
C
         IF(ISTOP.GT.0.AND.ISTOPT.GT.0)NTES = NTES-1
         IF(NTES.LE.0.OR.ISTOPT.LE.0.OR.ISTOP.LE.0.OR.NTES.GT.ICTS.OR.
     1       IPRINT.LT.0) GO TO 20
         KSUC = NSUC-ICTS
         CALL PTKSUC(KSUC)
         ICTS = NTES-1
  20     CONTINUE
         IF(NTES.LE.0.AND.ISTOPT.GT.0.AND.ISTOP.GT.0)GO TO 40
C
C  CONSTRAIN NOISE COEFFICIENT WITHIN BOUNDS
C
         EPSC = DMIN1(EPSC,EPSMAX)
         IF(EPSC.LE.0.D0)EPSC = 1.D0
  30     CONTINUE
C  END OF SERIES OF TRIALS
C
  40  CONTINUE
C
C  INDICATOR OF STOPPING CONDITIONS
C
      IOUT = ISTOPT
      FMIN = FOPT
C
      RETURN
C
      END
      SUBROUTINE INIT(NX,X0,H0,EPS0,DX0,IRAND,F,N1,N2,N3,N4
     1    ,INH,IFE,XRI,XRA,IT)
C
C  INIT  PERFORMS THE INITIALIZATION OF THE FIRST TRIAL.
C  THE PART OF THE INITIALIZATION WHICH IS COMMON ALSO TO THE TRIALS
C  FOLLOWING THE FIRST ONE IS PERFORMED BY CALLING THE SUBROUTINE
C  REINIT.
C
      DOUBLE PRECISION X0,H0,EPS0,DX0,F,XRI,XRA
      DOUBLE PRECISION FUNCT
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION X0(NX),XRI(NX),XRA(NX)
      DATA NPMAX/100/
      DATA NTRAJM/20/
      DATA NTRAJ0/7/
      DATA INKPB0/10/
      DATA KPBR00/3/
C
C CHECK PARAMETER VALUES
C
      IT = 0
      IF (NX.GT.NPMAX.OR.NX.LT.1.OR.H0.LE.0.D0.OR.EPS0.LE.0.D0
     1    .OR.DX0.LE.0.D0) IT = -99
      IF (IT.EQ.(-99)) RETURN
C
C INITIALIZE SOME VARIABLES
C
      INHP = INH
      DO 10 IX = 1,NX
         XRMIN(IX) = XRI(IX)
         XRMAX(IX) = XRA(IX)
  10     CONTINUE
      CALL NOSCA
      NTRAJ = N1
      IF(NTRAJ.EQ.0)NTRAJ = NTRAJ0
      NTRAJ = MIN0(NTRAJ,NTRAJM)
      NTRAJ = MAX0(NTRAJ,3)
      N1 = NTRAJ
      ISEGBR = N2
      IF(ISEGBR.EQ.0)ISEGBR = (1+NTRAJ)/2
      ISEGBR = MIN0(ISEGBR,NTRAJ)
      ISEGBR = MAX0(ISEGBR,1)
      N2 = ISEGBR
      INKPBR = N3
      IF(INKPBR.EQ.0)INKPBR = INKPB0
      N3 = INKPBR
      KPBR0 = N4
      IF(KPBR0.EQ.0)KPBR0 = KPBR00
      KPBR0 = MOD(KPBR0,INKPBR)
      N4 = KPBR0
      NDIM = NX
      NTRAJR = NTRAJ-1
      NCF = 1
      F = FUNCT(NX,X0)
      CALL STOOPT(NX,X0,F)
      DO 20 ID = 1,NTRAJ
         H(ID) = H0
         DX(ID) = DX0
  20     CONTINUE
C
C INITIALIZE REMAINING VARIABLES
C
      CALL REINIT(NX,X0,EPS0,IRAND,F,IFE)
C
      RETURN
      END
      SUBROUTINE REINIT(NX,X0,EPS0,IRAND,F,IFE)
C
C  REINIT  PERFORMS THE INITIALIZATION OF ALL TRIALS FOLLOWING THE
C  FIRST TRIAL, AND PART OF THE INITIALIZATION OF THE FIRST TRIAL.
C
      DOUBLE PRECISION X0,EPS0,F
      DOUBLE PRECISION EPSV,G
      DOUBLE PRECISION CHAOS
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION X0(NX)
      DATA EPSV/1.D0/
C
C INITIALIZE RANDOM NOISE GENERATOR
C
      G = CHAOS(IRAND)
C
      IFEP = IFE
      KGEN = 0
      KTIM = 0
      DO 30 ID = 1,NTRAJ
         DO 10 IC = 1,NDIM
            X(IC,ID) = X0(IC)
  10        CONTINUE
         IE(ID) = 0
         DO 20 IT = 1,NTRAJR
            ISVT(ID,IT) = 1
            VMVT(ID,IT) = F
  20        CONTINUE
         EPS(ID) = EPS0*EPSV**(ISEGBR-ID)
         VMCOR(ID) = F
         VCOR(ID) = F
  30     CONTINUE
      RETURN
      END
      SUBROUTINE TRIAL(N,NPMIN,NPMAX,KPASCA,
     1    TOLREL,TOLABS,IPRINT,XMIN,FMIN,
     2    FMAX,NFEV,NR,ISTOP)
C
C  THE SUBROUTINE  TRIAL  GENERATES A TRIAL, I.E. A SET OF COMPLETE
C  SIMULTANEOUS TRAJECTORIES BY REPEATEDLY PERFORMING
C     A CALL TO THE SUBROUTINE  GENEVA  WHICH GENERATES THE SET OF
C        SIMULTANEOUS TRAJECTORY SEGMENTS CORRESPONDING TO THE CURRENT
C        OBSERVATION PERIOD, AND PERFORMS THE TRAJECTORY SELECTION
C     A (POSSIBLE) CALL TO THE SUBROUTINE  PTSEG  WHICH PERFORMS
C        END-OF-SEGMENT OUTPUT
C     A CHECK OF THE TRIAL STOPPING CRITERIA (WITH THE AID OF THE
C        FUNCTION  ITOLCH
C     A DECISION ABOUT ACTIVATING OR DEACTIVATING THE SCALING OF
C        THE VARIABLES (ACTIONS PERFORMED BY CALLING THE SUBROUTINES
C        ACTSCA  AND  NOSCA).
C
      DOUBLE PRECISION TOLREL,TOLABS,XMIN,FMIN,FMAX
      DIMENSION XMIN(N)
      DATA IRNF/7/
C
      DO 20 IR = 1,NPMAX
C
C ACTIVATE SCALING
C
         IF(IR.GE.KPASCA.AND.IR.GT.N*IRNF)CALL ACTSCA
         NR = IR
C
C GENERATE AND EVALUATE THE SIMULTANEOUS TRAJECTORY SEGMENTS
C PERIOD
C
         CALL GENEVA(N,XMIN,FMIN,FMAX,NFEV)
C
C PRINT RESULTS OF OBSERVATION PERIOD
C
         IF(IPRINT.GT.0)CALL PTSEG(N,XMIN,FMIN,FMAX,IR,NFEV)
C
C CHECK TRIAL STOPPING CONDITIONS
C
         IF(IR.LT.NPMIN)GO TO 10
         ISTOP = ITOLCH(FMAX,FMIN,TOLREL,TOLABS)
         IF(ISTOP.NE.0)GO TO 30
C
  10     CONTINUE
  20     CONTINUE
  30  CONTINUE
      CALL NOSCA
C
      RETURN
      END
      SUBROUTINE GENEVA(NX,XMIN,FMIN,FMAX,NCEF)
C
C  GENEVA  PERFORMS THE GENERATION AND THE FINAL PROCESSING AND
C  EVALUATION OF THE SET OF TRAJECTORY SEGMENTS CORRESPONDING TO
C  THE CURRENT OBSERVATION PERIOD.
C  GENEVA  UPDATES THE SCALING ARRAYS  DIST  AND  BIAS  BY CALLING
C             THE SUBROUTINES  SEGSCA  AND  UPDSCA
C          GENERATES THE TRAJECTORY SEGMENTS BY CALLING THE SUB-
C             ROUTINE  PERIOD
C          DETERMINES SOME END-OF-SEGMENT RESULTS  (FPFMIN, FPFMAX,
C             XPFMIN)  USING THE RESCALING CAPABILITIES OF THE SUB-
C             ROUTINES  SEGSCA  AND  VARSCA.
C
      DOUBLE PRECISION XMIN,FMIN,FMAX
      DOUBLE PRECISION FM
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XMIN(NX)
C
C UPDATE SCALING DATA
C
      DO 10 ID = 1,NTRAJ
         CALL SEGSCA(ID)
         CALL UPDSCA(NX,X(1,ID))
  10     CONTINUE
C
C GENERATE THE SIMULTANEOUS TRAJECTORY SEGMENTS
C
      CALL PERIOD
C
C EXTRACT AND RESCALE SOME FINAL VALUES
C
      FM = VCOR(1)
      FMAX = VCOR(1)
      IFM = 1
      DO 30 ID = 2,NTRAJ
         FMAX = DMAX1(FMAX,VCOR(ID))
         IF(VCOR(ID).GE.FM)GO TO 20
         FM = VCOR(ID)
         IFM = ID
  20     CONTINUE
  30     CONTINUE
      DO 40 IC = 1,NDIM
         XMIN(IC) = X(IC,IFM)
  40     CONTINUE
      FMIN = FM
      NCEF = NCF
      CALL SEGSCA(IFM)
      CALL VARSCA(NX,XMIN)
C
      RETURN
      END
      SUBROUTINE PERIOD
C
C  PERIOD  IS CALLED BY SUBROUTINE  GENEVA  TO PERFORM THE GENERATION
C  OF THE TRAJECTORY SEGMENTS.
C  PERIOD  - COMPUTES THE DURATION OF THE OBSERVATION PERIOD, I.E. THE
C              NUMBER  NHP  OF ACCEPTED INTEGRATION STEPS IN A PERIOD
C          - COMPUTES ALL THE SEGMENT STEPS BY REPEATEDLY CALLING THE
C              SUBROUTINE  STEP
C          - PERFORMS THE SEGMENT SELECTION BY CALLING THE SUBROUTINE
C               BRASI
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
C DETERMINE DURATION OF OBSERVATION PERIOD
C (NUMBER OF TIME INTEGRATION STEPS)
C
      KGEN = KGEN+1
      NKGEN = 1
      IF(INHP.EQ.1)NKGEN=
     1    INT(SNGL(DLOG(DBLE(FLOAT(KGEN))+.5D0)/DLOG(2.D0)))+1
      IF(INHP.EQ.2)NKGEN = INT(SNGL(DSQRT(DBLE(FLOAT(KGEN))+.5D0)))
      IF(INHP.EQ.3)NKGEN = KGEN
C
C PERFORM THE REQUIRED NUMBER OF INTEGRATION STEPS
C (FOR ALL TRAJECTORIES)
C
      DO 10 IK = 1,NKGEN
C
C PERFORM A SINGLE INTEGRATION STEP
C (FOR ALL TRAJECTORIES)
C
         CALL STEP(IK)
  10     CONTINUE
C
C MANAGE THE TRAJECTORY BRANCHING PROCESS
C
      CALL BRASI
      RETURN
      END
      SUBROUTINE BRASI
C
C  BRASI  PERFORMS THE SELECTION PROCESS FOR THE TRAJECTORIES
C  BRASI  - UPDATES THE DATA ABOUT THE PAST TRAJECTORIES
C         - ASKS FOR THE TRAJECTORY-SELECTION ORDERING BY CALLING THE
C              SUBROUTINE  ORDER
C         - DISCARDS ONE OF THE TRAJECTORIES
C         - PERFORMS BRANCHING ON ONE OF THE REMAINING TRAJECTORIES
C         - MOVES THE DATA OF THE UNPERTURBED CONTINUATION
C              TO THE POSITION OF THE PERTURBED CONTINUATION
C         - CALLS THE SUBROUTINE  COMPAS  TO EXAMINE DATA ABOUT PAST
C              HISTORY OF THE TRAJECTORIES AND DISCARD IRRILEVANT DATA
C
      DOUBLE PRECISION DFAC,DLEPMX,DLFACL
      DOUBLE PRECISION EFAC,EPSMAX,FACG
      DOUBLE PRECISION CHAOS
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      DATA FACG/10.D0/
      DATA EFAC/.5D0/
      DATA DFAC/1.D3/
      DATA EPSMAX/1.D15/
      DATA DLFACL /0.301029995663981194D0/
      DATA DLEPMX /25.D0/
C
C UPDATE PAST HISTORY DATA
C
      DO 10 ID = 1,NTRAJ
         VMVT(ID,NTRAJR) = VMCOR(ID)
         ISVT(ID,NTRAJR) = ID
  10     CONTINUE
C
C OBTAIN TRAJECTORY-SELECTION ORDERING
C
      CALL ORDER(IP,IM,IU)
C
C DECIDE WHICH TRAJECTORY IS TO BE BRANCHED
C
      IF(MOD(KGEN,INKPBR).EQ.KPBR0)IM = IP
C
C PERFORM BRANCHING
C
      DO 20 IC = 1,NDIM
         X(IC,IU) = X(IC,IM)
  20     CONTINUE
      H(IU) = H(IM)
      IE(IU) = IE(IM)
      DO 30 IT = 1,NTRAJR
         ISVT(IU,IT) = ISVT(IM,IT)
         VMVT(IU,IT) = VMVT(IM,IT)
  30     CONTINUE
      EPS(IU) = EPS(IM)
      DX(IU) = DX(IM)
      VCOR(IU) = VCOR(IM)
      DO 40 ID = 1,NTRAJ
         VMCOR(ID) = VCOR(ID)
  40     CONTINUE
C
C UPDATE PAST-HISTORY-DATA MATRICES
C
      CALL COMPAS
C
C UPDATE SCALING DATA
C
      CALL MOVSCA(IU,IM)
C
C UPDATE NOISE COEFFICIENT VALUES
C
      IF (EPS(IU).LE.0.D0) GO TO 50
      IF(IFEP.EQ.2)
     1    EPS(IU) = EPS(IU)*FACG**(CHAOS(0)-EFAC)
      IF(IFEP.EQ.1)
     1    EPS(IU) = FACG**(DMIN1(DLEPMX,
     1    DLOG10(EPS(IU))+(CHAOS(-1)-EFAC)*DLFACL))
      EPS(IU) = DMIN1(EPS(IU),EPSMAX)
  50  CONTINUE
C
C UPDATE MAGNITUDE OF SPATIAL DISCRETIZATION INCREMENT
C
      DX(IU) = DX(IU)*DFAC**CHAOS(0)
      RETURN
      END
      SUBROUTINE ORDER(IP,IM,IU)
C
C  ORDER  COMPARES THE TRAJECTORIES FROM THE POINT OV VIEW OF PAST
C  HISTORY (BY CALLING THE FUNCTION  IPREC ) AND FROM THE POINT OF VIEW
C  OF THEIR CURRENT VALUE OF  EPS  (BY CALLING THE FUNCTION  IPRECE)
C  AND PROVIDES THE TRAJECTORY ORDERING NEEDED FOR SELECTING THE TRA-
C  JECTORY WHICH IS TO BE BRANCHED.
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION IORD(20)
      DATA KP/0/
      IR = 0
C
C ASSIGN INITIAL ORDERING
C
  10  CONTINUE
      DO 20 I = 1,NTRAJ
         IORD(I) = I
  20     CONTINUE
C
C SORT TRAJECTORIES ...
C
  30  CONTINUE
      IR = IR+1
C
      DO 60 I = 1,NTRAJR
         I1 = I+1
         DO 50 J = I1,NTRAJ
            K1 = IORD(I)
            K2 = IORD(J)
C
C ... ACCORDING TO PAST HISTORY ...
C
            IF(IR.NE.2)KP = IPREC(K1,K2)
            IF(KP.EQ.0.AND.IR.EQ.1)GO TO 10
C
C ... OR ACCORDING TO NOISE LEVEL
C
            IF(IR.EQ.2)KP = IPRECE(K1,K2)
            IF(KP.EQ.0)GO TO 40
            KM = K1+K2-KP
            IORD(I) = KP
            IORD(J) = KM
  40        CONTINUE
  50        CONTINUE
  60     CONTINUE
      IF(IR.EQ.2)GO TO 30
C
C RETURN THE INDICES OF THE SEGMENTS WHICH IN THE ORDERING
C OCCUPY THE FIRST, THE LAST, AND A SUITABLE MEDIUM LEVEL POSITION
C
      IP = IORD(1)
      IU = IORD(NTRAJ)
      IM = IORD(ISEGBR)
C
      RETURN
      END
      INTEGER FUNCTION IPREC(ID1,ID2)
C
C  IPREC  DETERMINES THE PRECEDENCE RELATION BETWEEN TWO TRAJECTORIES
C  BASED ON THE PAST HISTORY DATA
C
      DOUBLE PRECISION VM1,VM2
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      VM1 = VMVT(ID1,NTRAJR)
      VM2 = VMVT(ID2,NTRAJR)
      DO 10 IIT = 1,NTRAJR
         IT = 1+NTRAJR-IIT
         IF(ISVT(ID1,IT).EQ.ISVT(ID2,IT))GO TO 20
         VM1 = DMIN1(VM1,VMVT(ID1,IT))
         VM2 = DMIN1(VM2,VMVT(ID2,IT))
  10     CONTINUE
  20  CONTINUE
      IPREC = 0
      IF(VM2.LT.VM1)IPREC = ID2
      IF(VM2.GT.VM1)IPREC = ID1
C
      RETURN
      END
      INTEGER FUNCTION IPRECE(ID1,ID2)
C
C  IPRECE  DETERMINES THE PRECEDENCE RELATION BETWEEN TWO TRAJECTORIES
C  BASED ON THEIR CURRENT VALUE OF  EPS
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      IPRECE = 0
      IF(KGEN.GT.ISEGBR*INKPBR)GO TO 10
      IF(EPS(ID2).LT.EPS(ID1))IPRECE = ID1
      IF(EPS(ID2).GT.EPS(ID1))IPRECE = ID2
      RETURN
C
  10  CONTINUE
      IF(EPS(ID2).LT.EPS(ID1))IPRECE = ID2
      IF(EPS(ID2).GT.EPS(ID1))IPRECE = ID1
C
      RETURN
      END
      SUBROUTINE COMPAS
C
C  COMPAS  TAKES CARE OF THE STORAGE OF PAST HISTORY DATA, DISCARDING
C  ALL DATA NOT NEEDED BY THE ONLY USER OF SUCH DATA, THE FUNCTION
C  IPREC .
C
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      IT1 = 1
      DO 30 IT = 2,NTRAJR
         IT1 = IT-1
         DO 10 ID = 1,NTRAJ
            IF(ISVT(ID,IT).NE.ISVT(ID,IT1))GO TO 20
  10        CONTINUE
         GO TO 120
  20     CONTINUE
  30     CONTINUE
      IT1 = 1
      NCV = 0
      DO 50 ID = 1,NTRAJ
         DO 40 IDD = 1,NTRAJ
            IF(ISVT(ID,IT1).EQ.ISVT(IDD,IT1))NCV = NCV+1
  40        CONTINUE
  50     CONTINUE
      DO 80 IT = 2,NTRAJR
         IT1 = IT-1
         NCN = 0
         DO 70 ID = 1,NTRAJ
            DO 60 IDD = 1,NTRAJ
               IF(ISVT(ID,IT).EQ.ISVT(IDD,IT))NCN = NCN+1
  60           CONTINUE
  70        CONTINUE
         IF(NCN.EQ.NCV)GO TO 110
         NCV = NCN
  80     CONTINUE
      DO 90 ID = 1,NTRAJ
         IF(ISVT(1,1).NE.ISVT(ID,1))GO TO 100
  90     CONTINUE
      IT = 2
      GO TO 140
 100  CONTINUE
 110  CONTINUE
 120  CONTINUE
      IT = IT1+1
      DO 130 ID = 1,NTRAJ
         VMVT(ID,IT) = DMIN1(VMVT(ID,IT),VMVT(ID,IT1))
 130     CONTINUE
 140  CONTINUE
      DO 160 ITC = IT,NTRAJR
         ITM = ITC-1
         DO 150 ID = 1,NTRAJ
            VMVT(ID,ITM) = VMVT(ID,ITC)
            ISVT(ID,ITM) = ISVT(ID,ITC)
 150        CONTINUE
 160     CONTINUE
C
      RETURN
      END
      SUBROUTINE STEP(IK)
C
C  STEP  PERFORMS A SINGLE TIME-INTEGRATION STEP FOR EACH ONE OF THE
C  SIMULTANEOUS TRAJECTORIES BY REPEATEDLY CALLING THE SUBROUTINE  SSTEP
C
      DOUBLE PRECISION F
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      DOUBLE PRECISION XID,HID,EPSID,DXID
      DIMENSION XID(100)
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
C
      KTIM = KTIM+1
C
C LOOP ON THE SIMULTANEOUS TRAJECTORY SEGMENTS
      DO 30 ID = 1,NTRAJ
C
C INFORM THE SCALING SUBPROGRAMS (VIA THE COMMON AREA /SCALE/)
C THAT SCALING OPERATIONS ARE TO BE PERFORMED ON SEGMENT  ID .
         CALL SEGSCA(ID)
         F = VCOR(ID)
C
C PERFORM A TIME-INTEGRATION STEP ON SEGMENT  ID .
         KA = KTIM
         NX = NDIM
         DO 10 IX = 1,NX
            XID(IX) = X(IX,ID)
  10        CONTINUE
         HID = H(ID)
         EPSID = EPS(ID)
         DXID = DX(ID)
         IEID = IE(ID)
C
         CALL SSTEP(KA,NX,XID,HID,EPSID,DXID,IEID,F)
C
         DO 20 IX = 1,NX
            X(IX,ID) = XID(IX)
  20        CONTINUE
         H(ID) = HID
         EPS(ID) = EPSID
         DX(ID) = DXID
         IE(ID) = IEID
         VCOR(ID) = F
         VMCOR(ID) = DMIN1(VMCOR(ID),F)
         IF(IK.EQ.1)VMCOR(ID) = F
         IF(KTIM.EQ.1)IE(ID) = 0
  30     CONTINUE
C
      RETURN
      END
      SUBROUTINE SSTEP(KTIM,NDIM,X,H,EPS,DX,IE,F)
C
C  THE BASIC TIME-INTEGRATION STEP FOR A GIVEN TRAJECTORY IS PERFORMED
C  BY THE SUBROUTINE  STEP  WHICH
C   - CALLS THE FUNCTION  FUNCT0  TO COMPUTE THE VALUE OF  F
C   - CALLS THE SUBROUTINE  UNITRV  TO COMPUTE THE RANDOM DIRECTION
C        ALONG WHICH THE DIRECTIONAL DERIVATIVE  GAM/N  IS TO BE
C        COMPUTED
C   - CALLS THE SUBROUTINE  DERFOR (OR DERCEN ) TO COMPUTE THE FORWARD-
C        (OR CENTRAL-) DIFFERENCES DIRECTIONAL DERIVATIVE  GAM/N
C   - CALLS THE SUBROUTINE  NEWH  TO ACCEPT OR REJECT THE FIRST HALF-
C        STEP AND OBTAIN AN UPDATED VALUE FOR  H
C   - CALLS THE SUBROUTINE  CUMSCA  TO UPDATE THE CUMULATED SCALING DATA
C   - UPDATES THE SPATIAL DISCRETIZATION INCREMENT  DX  BASED ON THE
C        RESULTS OF CALLING THE FUNCTION  ITOLCH
C   - CALLS THE SUBROUTINE  GAUSRV  TO PERFORM THE SECOND HALF-STEP.
C
      DOUBLE PRECISION X,H,EPS,DX,F
      DOUBLE PRECISION DFDX,DFDXV,DXMAX,DXMIN
      DOUBLE PRECISION EPSR,FS,FV,FVS,HMINS,HR,HS
      DOUBLE PRECISION RCD,RDX,STF,TOLABS,TOLRA,TOLRI
      DOUBLE PRECISION W,XP
      DOUBLE PRECISION FUNCT0
      DIMENSION X(NDIM)
      DIMENSION W(100),XP(100)
      DATA RDX/1.D-4/
      DATA DXMIN/1.D-35/
      DATA DXMAX/1.D3/
      DATA HR/1.D-1/
      DATA HMINS/1.D-30/
      DATA STF /100.D0/
      DATA RCD/2.D0/
      DATA TOLRI/1.D-5/
      DATA TOLRA/1.D-11/
      DATA TOLABS/0.D0/
C
      IEC = 0
      FV = F
  10  CONTINUE
  20  CONTINUE
C
C TAKE A RANDOM DIRECTION FOR THE DIRECTIONAL DERIVATIVE
      CALL UNITRV(NDIM,W)
C
C COMPUTE FORWARD-DIFFERENCE DERIVATIVE
      CALL DERFOR(NDIM,X,FV,DX,W,DFDX)
C
C TRY THE FIRST HALF-STEP
      DO 30 IC = 1,NDIM
         XP(IC) = X(IC)-H*W(IC)*DFDX*DBLE(FLOAT(NDIM))
  30     CONTINUE
      HS = H
      F = FUNCT0(NDIM,XP)
      FVS = FV+DX*DABS(DFDX)
C
C UPDATE STEPLENGTH  H  AND ACCEPT OR REJECT THE HALF-STEP
      CALL NEWH(KTIM,FVS,F,H,IE,IEC)
      IF(IEC.LE.0)GO TO 50
      IE = IE-1
      IEC = IEC-1
      H = HS
      DFDXV = DFDX
C
C COMPUTE CENTRAL-DIFFERENCES DERIVATIVE
      CALL DERCEN(NDIM,X,FV,DX,W,DFDX)
C
C TRY AGAIN THE FIRST HALF-STEP
      DO 40 IC = 1,NDIM
         XP(IC) = X(IC)-H*W(IC)*DFDX*DBLE(FLOAT(NDIM))
  40     CONTINUE
      F = FUNCT0(NDIM,XP)
      FVS = FV+DX*DABS(DFDXV-DFDX)
C
C
C UPDATE STEPLENGTH  H  AND ACCEPT OR REJECT THE HALF-STEP
      CALL NEWH(KTIM,FVS,F,H,IE,IEC)
C
C UPDATE CUMULATED PAST SCALING DATA
      IF(IEC.GE.1)CALL CUMSCA(NDIM,W,DFDX)
      IF(IEC.GE.1)GO TO 10
      DX = DX*RDX
C
C FIRST HALF-STEP ACCEPTED
  50  CONTINUE
C
C UPDATE CUMULATED PAST SCALING DATA
      CALL CUMSCA(NDIM,W,DFDX)
      FS = FV+DX*DABS(DFDX)
      IF(ITOLCH(FS,FV,TOLRI,TOLABS).EQ.0)DX = DX/RCD
      IF(ITOLCH(FS,FV,TOLRA,TOLABS).GT.0)DX = DX*RCD
      EPSR = DSQRT(H)*EPS
C
C TAKE A SAMPLE INCREMENT OF THE WIENER PROCESS
      CALL GAUSRV(NDIM,W)
C
C TRY THE SECOND HALF-STEP
      DO 60 IC = 1,NDIM
         XP(IC) = XP(IC)+EPSR*W(IC)
  60     CONTINUE
      F = FUNCT0(NDIM,XP)
C
C ACCEPT OR REJECT THE COMPLETE STEP
      IF (F-FV.LE.EPS*EPS*STF) GO TO 70
      H = H*HR
      IE = IE+1
      IF (H.GT.HMINS) GO TO 20
C
C STEP ACCEPTED
  70  CONTINUE
      DO 80 IC = 1,NDIM
         X(IC) = XP(IC)
  80     CONTINUE
      DX = DMIN1(DX,DXMAX)
      DX = DMAX1(DX,DXMIN)
C
      RETURN
      END
      SUBROUTINE NEWH(K,FV,F,H,IE,IEC)
C
C  NEWH  IS CALLED BY THE SUBROUTINE  SSTEP  TO DECIDE WHETHER TO ACCEPT
C  OR REJECT THE FIRST HALF-STEP, AND TO PROVIDE AN UPDATED VALUE FOR  H
C
      DOUBLE PRECISION FV,F,H
      DOUBLE PRECISION FA,FR,HMAX,HMIN,R
      DIMENSION FR(3),FA(4)
      DATA FR/1.05D0,2.D0,10.D0/
      DATA FA/1.D0,1.1D0,2.D0,10.D0/
      DATA HMIN/1.D-30/
      DATA HMAX/1.D25/
      DATA IECMAX/50/
C
      IF(FV.LT.F)GO TO 20
C
C STEP ACCEPTED, POSSIBLY INCREASE THE STEPLENGTH  H
C
      R = FA(1)
      IF(IE*2.LT.K)R = FA(2)
      IF(IE*3.LT.K)R = FA(3)
      IF(IE.EQ.0.AND.K.GT.1)R = FA(4)
      H = H*R
  10  CONTINUE
      IEC = 0
      GO TO 30
C
  20  CONTINUE
      IE = IE+1
      IEC = IEC+1
      IF(IEC.GT.IECMAX)GO TO 10
C
C STEP REJECTED, DECREASE  H
C
      IC = MIN0(3,IEC)
      R = FR(IC)
      H = H/R
  30  CONTINUE
      H = DMIN1(H,HMAX)
      H = DMAX1(H,HMIN)
C
      RETURN
      END
      SUBROUTINE DERFOR(NDIM,X,F,DX,W,DFDX)
C
C  DERFOR  COMPUTED THE FORWARD FINITE-DIFFERNCES DIRECTIONAL
C  DERIVATIVES (CALLING  FUNCT0 )
C
      DOUBLE PRECISION X,F,DX,W,DFDX
      DOUBLE PRECISION DXF,DXFF,DXMAX,FN,S,XX
      DOUBLE PRECISION FUNCT0
      DIMENSION X(NDIM),W(NDIM)
      DIMENSION XX(100)
      DATA DXFF/1.D6/,DXF/1.D1/
      DATA DXMAX/1.D6/
C
  10  CONTINUE
  20  CONTINUE
      S = 0.D0
      DO 30 IC = 1,NDIM
         XX(IC) = X(IC)+W(IC)*DX
         S = S+(XX(IC)-X(IC))**2
  30     CONTINUE
      IF(S.GT.0.D0)GO TO 40
      DX = DX*DXFF
      GO TO 20
  40  CONTINUE
      FN = FUNCT0(NDIM,XX)
      DFDX = (FN-F)/DX
      IF(DX.GT.DXMAX)RETURN
      IF(DABS(DFDX).GT.1.D0)RETURN
      IF(DFDX**2.GT.0.D0)GO TO 50
      DX = DX*DXF
      GO TO 10
  50  CONTINUE
C
      RETURN
      END
      SUBROUTINE DERCEN(NDIM,X,F,DX,W,DFDX)
C
C  DERFOR  COMPUTED THE CENTRAL FINITE-DIFFERNCES DIRECTIONAL
C  DERIVATIVES (CALLING  FUNCT0 )
C
      DOUBLE PRECISION X,F,DX,W,DFDX
      DOUBLE PRECISION FN,FR,XX
      DOUBLE PRECISION FUNCT0
      DIMENSION X(NDIM),W(NDIM)
      DIMENSION XX(100)
C
      DO 10 IC = 1,NDIM
         XX(IC) = X(IC)-W(IC)*DX
  10     CONTINUE
      FR = FUNCT0(NDIM,XX)
      FN = F+DFDX*DX
      DFDX = (FN-FR)/(2.D0*DX)
C
      RETURN
      END
      SUBROUTINE RCLOPT(N,XO,FO)
C
C  RCLOPT  RECALLS THE CURRENT BEST MINIMUM VALUE  FOPT  FOUND SO FAR
C  FROM ALGORITHM START, AND THE POINT  XOPT  (OR POSSIBLY ONE OF THE
C  POINTS) WHERE  FOPT  WAS OBTAINED
C
      DOUBLE PRECISION XO,FO
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XO(N)
C
      DO 10 I = 1,N
         XO(I) = XOPT(I)
  10     CONTINUE
      FO = FOPT
C
      RETURN
      END
      SUBROUTINE STOOPT(N,XO,FO)
C
C  STOOPT  STORES THE CURRENT BEST MINIMUM VALUE  FOPT  FOUND SO FAR
C  FROM ALGORITHM START, AND THE POINT  XOPT  (OR POSSIBLY ONE OF THE
C  POINTS) WHERE  FOPT  WAS OBTAINED
C
      DOUBLE PRECISION XO,FO
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XO(N)
C
      DO 10 I = 1,N
         XOPT(I) = XO(I)
  10     CONTINUE
      FOPT = FO
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION FUNCT0(N,XX)
C
C  FUNCT0  IS CALLED  WHENEVER THE VALUE OF THE FUNCLION  F  IS REQUIRED
C  IN THE NUMERICAL INTEGRATION PROCESS.
C  THE FUNCTION  FUNCT0
C   - RESCALES THE VARIABLES BY CALLING THE SUBROUTINE  VARSCA
C   - CALLS THE SUBROUTINE  RANGE  TO TAKE CARE OF THE CASES WHERE THE
C      CURRENT POINT  X  IS OUTSIDE A GIVEN ADMISSIBLE REGION
C   - CALLS THE USER-SUPPLIED FUNCTION  FUNCT  TO ACTUALLY COMPUTE THE
C      VALUE OF  F
C   - POSSIBLY CALLS THE SUBROUTINE  STOOPT  TO UPDATE THE CURRENT
C        BEST MINIMUM FUNCTION VALUE  FOPT  AND THE CORRESPONDING
C        MINIMIZER  XOPT
C
      DOUBLE PRECISION XX
      DOUBLE PRECISION F,R,XS
      DOUBLE PRECISION FUNCT
      DOUBLE PRECISION X,H,DX,VMVT,EPS,VMCOR,VCOR
      DOUBLE PRECISION XRMIN,XRMAX,XOPT,FOPT
      COMMON /DINCOM/ X(100,20),H(20),DX(20),VMVT(20,19),EPS(20),
     1    VMCOR(20),VCOR(20),XRMIN(100),XRMAX(100),XOPT(100),FOPT,
     2    IE(20),ISVT(20,19),KGEN,KTIM,NDIM,NTRAJ,NTRAJR,
     3    ISEGBR,INKPBR,KPBR0,NCF,IFEP,INHP
      DIMENSION XX(N)
      DIMENSION XS(100)
C
      DO 10 IX = 1,N
         XS(IX) = XX(IX)
  10     CONTINUE
C
C DESCALE  X-VARIABLES
      CALL VARSCA(N,XS)
C
C CONSTRAIN THE  X-VARIABLES WITHIN BOUNDS
      CALL RANGE(N,XS,XRMIN,XRMAX,R)
C
C COMPUTE THE FUNCTION VALUE...
      F = FUNCT(N,XS)+R
C
C ... AND POSSIBLY UPDATE THE BEST CURRENT MINIMUM
      IF(F.LT.FOPT)CALL STOOPT(N,XS,F)
      FUNCT0 = F
      NCF = NCF+1
C
      RETURN
      END
      SUBROUTINE RANGE (N,XS,XRMIN,XRMAX,R)
C
C  RANGE  IS CALLED BY THE FUNCTION  FUNCT0  TO TAKE CARE OF THE CASES
C  WHERE THE CURRENT POINT  X  IS OUTSIDE A GIVEN ADMISSIBLE REGION
C
      DOUBLE PRECISION XS,XRMIN,XRMAX,R
      DOUBLE PRECISION A,B,C,D,DLRMAX,RMAX,RR,XC
      DIMENSION XS(N),XRMIN(N),XRMAX(N)
      DATA RMAX /1.D35/
      DATA DLRMAX/80.5904782547915990D0/
C
      R = 0.D0
      DO 40 I = 1,N
         A = XRMAX(I)
         C = XRMIN(I)
         XC = XS(I)
         IF (XC.LE.A) GO TO 10
         B = A+A-C
         RR = RMAX
         IF (XC.LT.B) RR = DEXP((XC-A)*DLRMAX/(B-A))-1.D0
         R = R+RR
         XS(I) = XRMAX(I)
         GO TO 30
  10     CONTINUE
         IF(XC.GE.C)GO TO 20
         D = C+C-A
         RR = RMAX
         IF (XC.GT.D) RR = DEXP((C-XC)*DLRMAX/(C-D))-1.D0
         R = R+RR
         XS(I) = XRMIN(I)
  20     CONTINUE
  30     CONTINUE
  40     CONTINUE
C
      RETURN
      END
      INTEGER FUNCTION ITOLCH(FMAX,FMIN,TOLREL,TOLABS)
C
C  ITOLCH  DETERMINES  WHETHER TWO QUANTITIES ARE TO BE CONSIDERED NU-
C  MERICALLY EQUAL WITH RESPECT TO AN ABSOLUTE (OR RELATIVE) DIFFERENCE
C  CRITERION, WITHIN GIVEN TOLERANCES
C
      DOUBLE PRECISION FMAX,FMIN,TOLREL,TOLABS
C
      ISTOP = 0
C
C CHECK RELATIVE DIFFERENCE AGAINST  TOLREL
      IF(DABS(FMAX-FMIN).LE.TOLREL*(DABS(FMIN)+DABS(FMAX))/2.D0)
     1    ISTOP=ISTOP+1
C
C CHECK ABSOLUTE DIFFERENCE AGAINST  TOLABS
      IF(FMAX-FMIN.LE.TOLABS)ISTOP = ISTOP+2
      ITOLCH = ISTOP
C
      RETURN
      END
      SUBROUTINE INISCA(N,ND)
C
C  INISCA  INITIALIZES THE COOMON AREA  /SCALE/  FOR THE SCALING DATA
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DATA NXMSCA/10/
C
      LSCA = -1
      IF(N.GT.NXMSCA.OR.N.EQ.1)RETURN
      LSCA = 0
      NX = N
      NORD = ND
      IDSCA = 1
      DO 1 ID = 1,NORD
         DO 2 IX = 1,NX
            DO 3 IY = 1,NX
               DIST(IX,IY,ID) = 0.D0
               GRAGRA(IX,IY,ID) = 0.D0
   3           CONTINUE
            DIST(IX,IX,ID) = 1.D0
            BIAS(IX,ID) = 0.D0
            GRA(IX,ID) = 0.D0
   2        CONTINUE
         NGRA(ID) = 0
   1     CONTINUE
C
      RETURN
      END
      SUBROUTINE NOSCA
C
C  NOSCA  DEACTIVATES THE SCALING
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      LSCA = -1
C
      RETURN
      END
      SUBROUTINE SEGSCA(ID)
C
C  SEGSCA  SELECTS THE TRAJECTORY WHICH MUST BE RESCALED
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      IDSCA = ID
C
      RETURN
      END
      SUBROUTINE VARSCA(N,X)
C
C  VARSCA  COMPUTES THE RESCALED VARIABLE  AX + B
C
      DOUBLE PRECISION X
      DOUBLE PRECISION XB
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION X(N),XB(10)
C
      IF(LSCA.LE.0)RETURN
      DO 1 I = 1,N
         XB(I) = BIAS(I,IDSCA)
         DO  2 J = 1,N
            XB(I) = XB(I)+DIST(I,J,IDSCA)*X(J)
   2        CONTINUE
   1     CONTINUE
      DO 3 I = 1,N
         X(I) = XB(I)
   3     CONTINUE
C
      RETURN
      END
      SUBROUTINE CUMSCA(N,W,DFDX)
C
C  CUMSCA  STORES CUMULATED STATISTICAL DATA ON THE ILL-CONDITIONING OF
C  F(AX+B) W.R.T. X
C
      DOUBLE PRECISION W,DFDX
      DOUBLE PRECISION DFDXMA
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION W(N)
      DATA DFDXMA/1.D8/
C
      IF(LSCA.LE.0)RETURN
      IF(DABS(DFDX).GT.DFDXMA)RETURN
      DO 1 I = 1,N
         DO 2 J = 1,N
            GRAGRA(I,J,IDSCA) = GRAGRA(I,J,IDSCA)+DFDX*W(I)*DFDX*W(J)
   2        CONTINUE
         GRA(I,IDSCA) = GRA(I,IDSCA)+DFDX*W(I)
   1     CONTINUE
      NGRA(IDSCA) = NGRA(IDSCA)+1
C
      RETURN
      END
      SUBROUTINE ACTSCA
C
C  ACTSCA  ACTIVATES THE RESCALING
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      IF(LSCA.EQ.0)LSCA = 1
C
      RETURN
      END
      SUBROUTINE MOVSCA(IU,IM)
C
C  MOVSCA  MOVES THE SCALING DATA OF THE UNPERTURBED CONTINUATION TO THE
C  POSITION OF THE PERTURBED CONTINUATION
C
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
C
      IF(LSCA.LT.0)RETURN
      DO 1 I = 1,NX
         DO 2 J = 1,NX
            DIST(I,J,IU) = DIST(I,J,IM)
            GRAGRA(I,J,IU) = GRAGRA(I,J,IM)
   2        CONTINUE
         BIAS(I,IU) = BIAS(I,IM)
         GRA(I,IU) = GRA(I,IM)
   1     CONTINUE
      NGRA(IU) = NGRA(IM)
C
      RETURN
      END
      SUBROUTINE UPDSCA(N,X)
C
C  UPDSCA  UPDATES THE SCALING MATRIX  A  AND THE BIAS VECTOR  B BY
C  CALLING  EIGSCA  AND  VARSCA
C
      DOUBLE PRECISION X
      DOUBLE PRECISION AGRAI,ALA1,ALPHA,AMCOR,BIAST
      DOUBLE PRECISION CD,COR,DISTT,SN
      DOUBLE PRECISION EIGSCA
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION X(N)
      DIMENSION DISTT(10,10),BIAST(10),COR(10,10)
      DATA ALPHA /0.3D0/
C
      IF(LSCA.LE.0)RETURN
      ID = IDSCA
      IF(NGRA(ID).LT.2*NX*NX)GO TO 2
      AGRAI = 1.D0/DBLE(FLOAT(NGRA(ID)))
      AMCOR = 0.D0
      DO 3 I = 1,NX
         DO 4 J = 1,NX
         COR(I,J) = AGRAI*GRAGRA(I,J,ID)-AGRAI*GRA(I,ID)*AGRAI*GRA(J,ID)
            AMCOR = DMAX1(AMCOR,DABS(COR(I,J)))
   4        CONTINUE
   3     CONTINUE
      IF(AMCOR.LE.0.D0)GO TO 12
      DO 1 I = 1,NX
         DO 11 J = 1,NX
            COR(I,J) = COR(I,J)/AMCOR
  11        CONTINUE
   1     CONTINUE
      ALA1 = EIGSCA(COR)
      CD = ALA1*(1.D0+ALPHA)
      DO 5 I = 1,NX
         COR(I,I) = COR(I,I)-CD
         BIAST(I) = X(I)
   5     CONTINUE
      CALL VARSCA(NX,BIAST)
      SN = 0.D0
      DO 6 I = 1,NX
         DO 7 J = 1,NX
            DISTT(I,J) = 0.D0
            DO 8 K = 1,NX
               DISTT(I,J) = DISTT(I,J)-DIST(I,K,ID)*COR(K,J)
   8           CONTINUE
            SN = SN+DISTT(I,J)**2
   7        CONTINUE
   6     CONTINUE
      SN = 1.D0/DSQRT(SN/DBLE(FLOAT(NX)))
      DO 9 I = 1,NX
         BIAS(I,ID) = BIAST(I)
         DO 10 J = 1,N
            DIST(I,J,ID) = SN*DISTT(I,J)
            BIAS(I,ID) = BIAS(I,ID)-DIST(I,J,ID)*X(J)
            GRAGRA(I,J,ID) = 0.D0
  10        CONTINUE
         GRA(I,ID) = 0.D0
   9     CONTINUE
      NGRA(ID) = 0
   2  CONTINUE
  12  CONTINUE
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION EIGSCA(COR)
C
C  EIGSCA COMPUTES THE LARGEST EIGENVALUE OF A MATRIX USED FOR RESCA-
C  LING, STARTING FROM A RANDOMLY-CHOSEN ESTIMATE (OBTAINED BY CALLING
C  THE SUBROUTINE  UNITRV ) OF THE CORRESPONDING EIGENVECTOR
C
      DOUBLE PRECISION COR
      DOUBLE PRECISION ALA1,ALA1I,ALA1O
      DOUBLE PRECISION PREC,SWW,W,WW
      DOUBLE PRECISION DIST,BIAS,GRAGRA,GRA
      COMMON /SCALE/ DIST(10,10,20),BIAS(10,20),GRAGRA(10,10,20),
     1    GRA(10,20),NGRA(20),LSCA,IDSCA,NX,NORD
      DIMENSION COR(10,10)
      DIMENSION W(10),WW(10)
      DATA PREC /1.D-3/
      DATA NRMIN /10/
      DATA NRMAX /100/
C
      CALL UNITRV(NX,W)
      ALA1 = 0.D0
      DO 1 IR = 1,NRMAX
         ALA1O = ALA1
         SWW = 0.D0
         DO 2 IX = 1,NX
            WW(IX) = 0.D0
            DO 3 JX = 1,NX
               WW(IX) = WW(IX)+COR(IX,JX)*W(JX)
   3           CONTINUE
            SWW = SWW+WW(IX)**2
   2        CONTINUE
         ALA1 = DSQRT(SWW)
         ALA1I = 1.D0/ALA1
         IF(IR.GE.NRMIN.AND.ALA1*PREC.GT.DABS(ALA1-ALA1O))GO TO 4
         DO 5 IX = 1,NX
            W(IX) = WW(IX)*ALA1I
   5        CONTINUE
   1     CONTINUE
   4  CONTINUE
      EIGSCA = ALA1
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION CHAOS(INIZ)
C
C  CHAOS  GENERATES A RANDOM SAMPLE FROM ONE OUT OF FOUR POSSIBLE
C PROBABILITY DISTRIBUTIONS USING RANDOM NUMBERS UNIFORMLY
C DISTRIBUTED IN (0,1) GENERATED BY THE FUNCTION  UNIFRD
C
C  INIZ  IS AN INPUT PARAMETER AS FOLLOWS
C
C      INIZ>0          INITIALIZATION WITH SEED  INIZ.
C      INIZ=0          STANDARD GAUSSIAN DISTRIBUTION.
C      INIZ=-1         CAUCHY DISTRIBUTION.
C      INIZ=-2         UNIFORM DISTRIBUTION IN (-1,+1).
C      OTHERWISE       UNIFORM DISTRIBUTION IN ( 0,+1).
C
      DOUBLE PRECISION UNIFRD,PAI,A,B
C
      DATA PAI/3.14159265358979324D0/
C
      IF(INIZ.LE.0) GO TO 10
C
C INITIALIZATION.
C
      CHAOS = UNIFRD(INIZ)
      RETURN
C
  10  CONTINUE
      A = UNIFRD(0)
      IF(INIZ.NE.0) GO TO 20
      B = UNIFRD(0)
C
C GAUSSIAN RANDOM NUMBER BY POLAR METHOD
C
      CHAOS = DSQRT(-2.D0*DLOG(A))*DCOS(PAI*B)
C
      RETURN
  20  CONTINUE
C
C UNIFORM RANDOM NUMBER IN (0,1)
C
      CHAOS = A
C
C CAUCHY RANDOM NUMBER BY INVERSE TRANSFORMATION
C
      IF(INIZ.EQ.(-1)) CHAOS = DSIN(PAI*A)/DCOS(PAI*A)
C
C UNIFORM RANDOM NUMBER IN (-1,+1)
C
      IF(INIZ.EQ.(-2)) CHAOS = 2.D0*A-1.D0
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION UNIFRD(INIZ)
C
C UNIFRD GENERATES THE RANDOM NUMBERS UNIFORMLY DISTRIBUTED IN (0,1)
C EXPLOITING THOSE GENERATED BY  ALKNUT  WITH A FURTHER RANDOMIZATION
C IF THE INPUT PARAMETER INIZ  IS NOT 0
C THE RANDOM NUMBER GENERATOR IS INITIALIZED
C
      DOUBLE PRECISION A,B,C,X
      DOUBLE PRECISION X0,P,P0,P1,P2,R1,R2
      DOUBLE PRECISION FINV
C
      DIMENSION X(61)
C
      DATA NREM/61/,NRIP/100/
      DATA A,B,C/-1.5D0,5.5D0,-2.0D0/
      DATA FINV/3.5D-5/
      DATA IREM/0/
      DATA P0/3.D0/,P1,P2/1.D0,3.D0/,R1,R2/0.25D0,0.75D0/
C
      IF(INIZ.NE.0.OR.IREM.EQ.0) GO TO 10
C
      I0 = IREM
      X0 = X(I0)
C
C NONLINEARIZATION OF X0 TO AVOID LONG-DISTANCE LINEAR RELATIONSHIP.
C
      IF(X0.GE.FINV)X0 = DMOD(1.D0/X0,1.D0)
C
C UPDATE A COMPONENT OF THE VECTOR X ...
C
      CALL ALKNUT(NREM,X,IREM)
C
C ... AND FURTHER RANDOMIZE
C
      UNIFRD = DMOD(X0+X(I0),1.D0)
C
      RETURN
C
C INITIALIZATION OF THE RANDOM NUMBER GENERATOR
C
  10  CONTINUE
C
      P = P0-1.D0/DBLE(FLOAT(IABS(INIZ))+100.0)
      DO 20 K = 1,NREM
C
         P = C+P*(B+P*A)
C
C
         X(K) = R1+(R2-R1)*(P-P1)/(P2-P1)
  20     CONTINUE
C
      IREM = 0
C
      DO 30 K = 1,NRIP
C
         CALL ALKNUT(NREM,X,IREM)
C
  30     CONTINUE
C
      UNIFRD = X(1)
C
      RETURN
      END
      SUBROUTINE ALKNUT(NREM,X,IREM)
C
C UPDATES THE COMPONENT IREM OF THE NREM-VECTOR X WITH A RANDOM NUMBER
C UNIFORMLY DISTRIBUTED IN (0,1) BY MEANS OF THE ALGORITHM
C OF MITCHELL-MOORE, MODIFIED AS SUGGESTED BY BRENT, QUOTED IN
C D.E.KNUTH, THE ART OF COMPUTER PROGRAMMING, SECOND EDITION,
C SECOND VOLUME, SEMINUMERICAL ALGORITHMS, ADDISON-WESLEY
C PUB. CO., READING (1981), PP. 26-28.
C
      DOUBLE PRECISION X
C
      DIMENSION X(NREM)
C
      DATA N1,N2/24,55/
C
      IF(IREM.NE.0) GO TO 10
C
      IREM = NREM
      I1 = NREM-N1
      I2 = NREM-N2
C
  10  CONTINUE
C
      X(IREM) = DMOD(X(I1)+X(I2),1.D0)
C
      IREM = 1+MOD(IREM,NREM)
      I1 = 1+MOD(I1,NREM)
      I2 = 1+MOD(I2,NREM)
C
      RETURN
      END
      SUBROUTINE GAUSRV(N,W)
C
C GENERATES A RANDOM VECTOR SAMPLE FROM AN N-DIMENSIONAL
C NORMAL DISTRIBUTION
C
      DOUBLE PRECISION W,X,Y,R
      DOUBLE PRECISION CHAOS
C
      DIMENSION W(N)
C
      NN = (N+1)/2
C
      DO 20 I = 1,NN
         II = 1+N-I
C
  10     CONTINUE
         X = CHAOS(-2)
         Y = CHAOS(-2)
         R = X*X+Y*Y
C
         IF(R.GT.1.D0) GO TO 10
C
         R = DSQRT(-2.D0*DLOG(R)/R)
C
         W(I) = X*R
         W(II) = Y*R
C
  20     CONTINUE
C
      RETURN
      END
      SUBROUTINE UNITRV(N,W)
C
C GENERATES A RANDOM VECTOR UNIFORMLY DISTRIBUTED
C ON THE UNIT SPHERE.
C
      DOUBLE PRECISION W,WW
C
      DIMENSION W(N)
C
      CALL GAUSRV(N,W)
      WW = 0.D0
C
      DO 10 I = 1,N
         WW = WW+W(I)**2
  10     CONTINUE
      WW = 1.D0/DSQRT(WW)
C
      DO 20 I = 1,N
         W(I) = WW*W(I)
  20     CONTINUE
C
      RETURN
      END
      SUBROUTINE GLOMTF (NPROB, N, X, FUNZ)
C
C     THE SUBROUTINE  GLOMTF  PROVIDES THE CODING OF 37 REAL-VALUED
C     FUNCTIONS OF  N  REAL VARIABLES, TO BE USED, TOGETHER WITH THE
C     SUBROUTINE  GLOMIP , TO DEFINE 37 TEST PROBLEMS FOR GLOBAL
C     MINIMIZATION SOFTWARE.
C
C     THE SUBROUTINE  GLOMTF  RETURNS IN  FUNZ  THE FUNCTION VALUE
C     AT THE POINT  X = (X1,X2,...,XN)  FOR THE FUNCTION DEFINED BY
C     PROBLEM NUMBER  NPROB .
C
C     CALLING STATEMENT
C
C        CALL GLOMTF (NPROB, N, X, FUNZ)
C
C     DESCRIPTION OF THE CALL PARAMETERS
C        (THE FORTRAN IMPLICIT TYPE DEFINITION FOR INTEGERS IS USED.
C        ALL NON-INTEGER PARAMETERS ARE DOUBLE-PRECISION).
C
C     NPROB   IS THE (INPUT) TEST-PROBLEM NUMBER
C
C     N   IS THE (INPUT) DIMENSION OF THE PROBLEM
C
C     X   IS AN (INPUT) N-VECTOR CONTAINING THE INDEPENDENT VARIABLES
C
C     FUNZ   IS THE (OUTPUT) VALUE AT  X  OF THE FUNCTION DEFINED BY
C        PROBLEM NUMBER  NPROB .
C
C
C
      DOUBLE PRECISION X
      DOUBLE PRECISION FUNZ
C
      DOUBLE PRECISION PI
      DOUBLE PRECISION P1A
      DOUBLE PRECISION P4A
      DOUBLE PRECISION P5A
      DOUBLE PRECISION P6A
      DOUBLE PRECISION P9A,P9B
      DOUBLE PRECISION P17A
      DOUBLE PRECISION P20A,P20B
      DOUBLE PRECISION P21A,P21B,P21C,P21D
      DOUBLE PRECISION P22A,P22B,P22C,P22D
      DOUBLE PRECISION P36A,P36B,P36C,P36D
      DOUBLE PRECISION P37A,P37B,P37C,P37D
      DOUBLE PRECISION Y
      DOUBLE PRECISION XX,DI,YY,BETA,S1,S2,A,B,R2,R4,R8,U,V,UU,VV,S
      DOUBLE PRECISION RG,RP,H,PUND,RANGE
      DOUBLE PRECISION FACT,VAR,PENFUN,RAN
C
      DIMENSION X(N)
C
      DIMENSION P20A(10,4),P20B(10)
      DIMENSION P21A(4,3),P21B(4,3),P21C(4)
      DIMENSION P22A(4,6),P22B(4,6),P22C(4)
      DIMENSION Y(10)
C
      DATA PI /3.14159265358979323846D0/
      DATA P1A /0.1D0/
      DATA P4A /0.1D0/
      DATA P5A /0.1D0/
      DATA P6A /0.1D0/
      DATA P9A,P9B /-1.4251284283197609708D0,-0.80032110047197312466D0/
      DATA P17A /1.275D0/
      DATA P20A /4.D0,1.D0,8.D0,6.D0,3.D0,2.D0,5.D0,8.D0,6.D0,7.D0,
     1           4.D0,1.D0,8.D0,6.D0,7.D0,9.D0,5.D0,1.D0,2.D0,3.6D0,
     2           4.D0,1.D0,8.D0,6.D0,3.D0,2.D0,3.D0,8.D0,6.D0,7.D0,
     3           4.D0,1.D0,8.D0,6.D0,7.D0,9.D0,3.D0,1.D0,2.D0,3.6D0/
      DATA P20B /0.1D0,0.2D0,0.2D0,0.4D0,0.4D0,0.6D0,0.3D0,0.7D0,
     1    0.5D0,0.5D0/
      DATA P21A /3.D0,0.1D0,3.D0,0.1D0,
     1          10.D0,10.D0,10.D0,10.D0,
     2          30.D0,35.D0,30.D0,35.D0/
      DATA P21B /0.3689D0,0.4699D0,0.1091D0,0.03815D0,
     1           0.1170D0,0.4387D0,0.8732D0,0.5743D0,
     2           0.2673D0,0.7470D0,0.5547D0,0.8828D0/
      DATA P21C /1.D0,1.2D0,3.D0,3.2D0/
      DATA P21D /-69.D0/
      DATA P22A /10.D0,0.05D0,3.D0,17.D0,3.D0,10.D0,3.5D0,8.D0,
     1           17.D0,17.D0,1.7D0,0.05D0,3.5D0,0.1D0,10.D0,10.D0,
     2           1.7D0,8.D0,17.D0,0.1D0,8.D0,14.D0,8.D0,14.D0/
      DATA P22B /0.1312D0,0.2329D0,0.2348D0,0.4047D0,
     1           0.1696D0,0.4135D0,0.1451D0,0.8828D0,
     2           0.5569D0,0.8307D0,0.3522D0,0.8732D0,
     3           0.0124D0,0.3736D0,0.2883D0,0.5743D0,
     4           0.8283D0,0.1004D0,0.3047D0,0.1091D0,
     5           0.5886D0,0.9991D0,0.6650D0,0.0381D0/
      DATA P22C /1.D0,1.2D0,3.D0,3.2D0/
      DATA P22D /-69.D0/
      DATA P36A,P36B,P36C,P36D /10.D1,1.D0,10.D0,0.98D0/
      DATA P37A,P37B,P37C,P37D /10.D0,1.D0,10.D0,0.98D0/
C
C     PENALIZATION FUNCTION
C
      PENFUN (VAR,RAN,FACT,IEXP) = FACT*(DABS(VAR)-RAN)**IEXP
C
      GO TO (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,
     1       160,170,180,190,200,210,220,230,240,250,260,270,280,
     2       290,300,310,320,330,340,350,360,370), NPROB
C
C    1  A FOURTH ORDER POLYNOMIAL  (N = 1)
C
   10 CONTINUE
      FUNZ = ((0.25D0*X(1)*X(1)-0.5D0)*X(1)+P1A)*X(1)
      RETURN
C
C    2  GOLDSTEIN SIXTH ORDER POLYNOMIAL  (N = 1)
C
   20 CONTINUE
      XX = X(1)*X(1)
      FUNZ = ((XX-15.D0)*XX+27.D0)*XX+250.D0
      RETURN
C
C    3  ONE-DIMENSIONAL PENALIZED SHUBERT FUNCTION  (N = 1)
C
   30 CONTINUE
      FUNZ = 0.D0
      DO 35 I = 1,5
         DI = DBLE(FLOAT(I))
         FUNZ = FUNZ+DI*DCOS((DI+1.D0)*X(1)+DI)
   35 CONTINUE
      IF (DABS(X(1)).GT.10.D0) FUNZ = FUNZ+PENFUN(X(1),10.D0,100.D0,2)
      RETURN
C
C    4  A FOURTH ORDER POLYNOMIAL IN TWO VARIABLES  (N = 2)
C
   40 CONTINUE
      XX = X(1)*X(1)
      YY = X(2)*X(2)
      FUNZ = 0.25D0*XX*XX-0.5D0*XX+P4A*X(1)+0.5D0*YY
      RETURN
C
C    5  A FUNCTION WITH A SINGLE ROW OF LOCAL MINIMA  (N = 2)
C
   50 CONTINUE
      FUNZ = 0.5D0*(P5A*X(1)*X(1)+1.D0-DCOS(2.D0*X(1)))+X(2)*X(2)
      RETURN
C
C    6  SIX-HUMP CAMEL FUNCTION  (N = 2)
C
   60 CONTINUE
      XX = X(1)*X(1)
      YY = X(2)*X(2)
      FUNZ = ((XX/3.D0-(2.D0+P6A))*XX+4.D0)*XX+X(1)*X(2)
     1    +4.D0*(YY-1.D0)*YY
      RETURN
C
C    7  TWO-DIMENSIONAL PENALIZED SHUBERT FUNCTION, BETA = 0  (N = 2)
C
   70 CONTINUE
      BETA = 0.D0
      GO TO 93
C
C    8  TWO-DIMENSIONAL PENALIZED SHUBERT FUNCTION, BETA = 1/2  (N = 2)
C
   80 CONTINUE
      BETA = 0.5D0
      GO TO 93
C
C    9  TWO-DIMENSIONAL PENALIZED SHUBERT FUNCTION, BETA = 1  (N = 2)
C
   90 CONTINUE
      BETA = 1.D0
C
   93 CONTINUE
      FUNZ = ((X(1)-P9A)**2+(X(2)-P9B)**2)*BETA
      S1 = 0.D0
      S2 = 0.D0
      DO 95 I = 1,5
         DI = DBLE(FLOAT(I))
         S1 = S1+DI*DCOS((DI+1.D0)*X(1)+DI)
         S2 = S2+DI*DCOS((DI+1.D0)*X(2)+DI)
   95 CONTINUE
      FUNZ = FUNZ+S1*S2
      IF (DABS(X(1)).GT.10.D0) FUNZ = FUNZ+PENFUN(X(1),10.D0,100.D0,2)
      IF (DABS(X(2)).GT.10.D0) FUNZ = FUNZ+PENFUN(X(2),10.D0,100.D0,2)
      RETURN
C
C   10  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10  (N = 2)
C
  100 CONTINUE
      A = 1.D1
      B = 1.D-1
      GO TO 155
C
C   11  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**2  (N = 2)
C
  110 CONTINUE
      A = 1.D2
      B = 1.D-2
      GO TO 155
C
C   12  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**3  (N = 2)
C
  120 CONTINUE
      A = 1.D3
      B = 1.D-3
      GO TO 155
C
C   13  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**4  (N = 2)
C
  130 CONTINUE
      A = 1.D4
      B = 1.D-4
      GO TO 155
C
C   14  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**5  (N = 2)
C
  140 CONTINUE
      A = 1.D5
      B = 1.D-5
      GO TO 155
C
C   15  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**6  (N = 2)
C
  150 CONTINUE
      A = 1.D6
      B = 1.D-6
C
  155 CONTINUE
      XX = X(1)*X(1)
      YY = X(2)*X(2)
      R2 = XX+YY
      R4 = R2*R2
      R8 = R4*R4
      FUNZ = A*XX+YY-R4+B*R8
      RETURN
C
C   16  GOLDSTEIN-PRICE FUNCTION  (N = 2)
C
  160 CONTINUE
      U = X(1)+X(2)+1.D0
      V = 2.D0*X(1)-3.D0*X(2)
      UU = U*U
      VV = V*V
      FUNZ = (1.D0+UU*(36.D0-20.D0*U+3.D0*UU))
     1    * (30.D0+VV*(18.D0-16.D0*V+3.D0*VV))
      RETURN
C
C   17  PENALIZED BRANIN FUNCTION  (N = 2)
C
  170 CONTINUE
      FUNZ = (X(2)-P17A*(X(1)/PI)**2+(5.D0/PI)*X(1)-6.D0)**2
     1    +10.D0*(1.D0-1.D0/(8.D0*PI))*DCOS(X(1))+10.D0
      IF (DABS(X(1)-2.5D0).GT.7.5D0)
     1   FUNZ = FUNZ+PENFUN(X(1)-2.5D0,7.5D0,100.D0,2)
      IF (DABS(X(2)-7.5D0).GT.7.5D0)
     1   FUNZ = FUNZ+PENFUN(X(2)-7.5D0,7.5D0,100.D0,2)
      RETURN
C
C   18  PENALIZED SHEKEL FUNCTION, M = 5  (N = 4)
C
  180 CONTINUE
      M = 5
      GO TO 203
C
C   19  PENALIZED SHEKEL FUNCTION, M = 7  (N = 4)
C
  190 CONTINUE
      M = 7
      GO TO 203
C
C   20  PENALIZED SHEKEL FUNCTION, M = 10  (N = 4)
C
  200 CONTINUE
      M = 10
C
  203 CONTINUE
      FUNZ = 0.D0
      DO 207 I = 1,M
         S = P20B(I)
         DO 205 J = 1,N
            S = S+(X(J)-P20A(I,J))**2
  205    CONTINUE
         FUNZ = FUNZ - 1.D0/S
  207 CONTINUE
      DO 209 I = 1,N
         IF (DABS(X(I)-5.D0).GT.5.D0)
     1      FUNZ = FUNZ+PENFUN(X(I)-5.D0,5.D0,100.D0,2)
  209 CONTINUE
      RETURN
C
C   21  PENALIZED THREE-DIMENSIONAL HARTMAN FUNCTION  (N = 3)
C
  210 CONTINUE
      M = 4
      FUNZ = 0.D0
      DO 217 I = 1,M
         S = 0.D0
         DO 213 J = 1,N
            S = S-P21A(I,J)*(X(J)-P21B(I,J))**2
  213    CONTINUE
         IF (S.GE.P21D) FUNZ = FUNZ-P21C(I)*DEXP(S)
  217 CONTINUE
      GO TO 227
C
C   22  PENALIZED SIX-DIMENSIONAL HARTMAN FUNCTION  (N = 6)
C
  220 CONTINUE
      M = 4
      FUNZ = 0.D0
      DO 225 I = 1,M
         S = 0.D0
         DO 223 J = 1,N
            S = S-P22A(I,J)*(X(J)-P22B(I,J))**2
  223    CONTINUE
         IF (S.GE.P22D) FUNZ = FUNZ-P22C(I)*DEXP(S)
  225 CONTINUE
  227 CONTINUE
         DO 229 I = 1,N
         IF (DABS(X(I)-0.5D0).GT.0.5D0)
     1      FUNZ = FUNZ+PENFUN(X(I)-0.5D0,0.5D0,100.D0,2)
  229 CONTINUE
      RETURN
C
C   23  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 1  (N = 2)
C
  230 CONTINUE
C
C   24  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 1  (N = 3)
C
  240 CONTINUE
C
C   25  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 1  (N = 4)
C
  250 CONTINUE
      DO 255 I = 1,N
         Y(I) = 1.D0+0.25D0*(X(I)-1.D0)
  255 CONTINUE
      GO TO 285
C
C   26  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 2  (N = 5)
C
  260 CONTINUE
C
C   27  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 2  (N = 8)
C
  270 CONTINUE
C
C   28  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 2  (N = 10)
C
  280 CONTINUE
      DO 283 I = 1,N
         Y(I) = X(I)
  283 CONTINUE
C
  285 CONTINUE
      FUNZ = 10.D0*DSIN(PI*Y(1))**2+(Y(N)-1.D0)**2
      DO 287 I=2,N
         FUNZ = FUNZ+(Y(I-1)-1.D0)**2*(1.D0+10.D0*DSIN(PI*Y(I))**2)
  287 CONTINUE
      FUNZ = FUNZ*PI/DBLE(FLOAT(N))
      RANGE = 10.D0
      GO TO 347
C
C   29  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 10   (N = 2)
C
  290 CONTINUE
C
C   30  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 10   (N = 3)
C
  300 CONTINUE
C
C   31  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 10   (N = 4)
C
  310 CONTINUE
      RANGE = 10.D0
      GO TO 342
C
C   32  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 5   (N = 5)
C
  320 CONTINUE
C
C   33  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 5   (N = 6)
C
  330 CONTINUE
C
C   34  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 5  (N = 7)
C
  340 CONTINUE
      RANGE = 5.D0
  342 CONTINUE
      FUNZ = DSIN(3.D0*PI*X(1))**2
     1   +(X(N)-1.D0)**2*(1.D0+DSIN(2.D0*PI*X(N))**2)
      DO 345 I=2,N
         FUNZ = FUNZ+(X(I-1)-1.D0)**2*(1.D0+DSIN(3.D0*PI*X(I))**2)
  345 CONTINUE
      FUNZ = FUNZ*0.1D0
  347 CONTINUE
      DO 349 I = 1,N
         IF (DABS(X(I)).GT.RANGE)
     1      FUNZ = FUNZ+PENFUN(X(I),RANGE,100.D0,4)
  349 CONTINUE
      RETURN
C
C   35  A FUNCTION WITH A SINGLE CUSP-SHAPED MINIMUM  (N = 5)
C
  350 CONTINUE
      FUNZ = 0.D0
      DO 355 I = 1,N
         FUNZ = FUNZ+DBLE(FLOAT(I))*X(I)*X(I)
  355 CONTINUE
      FUNZ = DSQRT(DSQRT(FUNZ))
      RETURN
C
C   36  A FUNCTION WITH A SMALL-ATTRACTION-REGION GLOBAL MINIMUM (N = 2)
C
  360 CONTINUE
      RG = P36A
      RP = P36B
      H = P36C
      PUND = P36D
      GOTO 373
C
C   37  A FUNCTION WITH A SMALL-ATTRACTION-REGION GLOBAL MINIMUM (N = 5)
C
  370 CONTINUE
      RG = P37A
      RP = P37B
      H = P37C
      PUND = P37D
  373 CONTINUE
      S = 0.D0
      DO 377 I = 2,N
         S = S+X(I)*X(I)
  377 CONTINUE
      FUNZ = S+X(1)*X(1)
      S = S+(X(1)-RG)**2
      IF (S.LT.RP*RP*PUND) FUNZ = FUNZ-(RG*RG+H)*DEXP(-S/(RP*RP-S))
      RETURN
C
      END
      SUBROUTINE GLOMIP (NPROB, N, X0, XMIN, XMAX)
C
C     THE SUBROUTINE  GLOMIP  PROVIDES THE CODING FOR THE NUMBER
C     OF VARIABLES, THE INITIAL POINT, AND THE OBSERVATION REGION
C     TO BE USED, TOGETHER WITH THE 37 TEST FUNCTIONS GIVEN BY
C     SUBROUTINE  GLOMTF , TO DEFINE 37 TEST PROBLEMS FOR GLOBAL
C     MINIMIZATION SOFTWARE.
C
C     THE SUBROUTINE  GLOMIP  RETURNS IN  N , X0 , AND  XMIN ,
C     XMAX  THE NUMBER OF VARIABLES, THE INITIAL POINT, AND THE
C     BOUNDARIES OF THE OBSERVATION REGION.
C
C     CALLING STATMENT
C        CALL GLOMIP (NPROB, N, X0, XMIN, XMAX)
C
C     DESCRIPTION OF THE CALL PARAMETERS
C        (THE FORTRAN IMPLICIT TYPE DEFINITION FOR INTEGERS IS USED.
C        ALL NON-INTEGER PARAMETERS ARE DOUBLE-PRECISION).
C
C     NPROB   IS THE (INPUT) NUMBER OF THE TEST PROBLEM TO BE
C             CONSIDERED
C     N       IS THE (OUTPUT) NUMBER OF VARIABLES (DIMENSION) OF
C             THE PROBLEM
C     XMIN , XMAX   ARE THE (OUTPUT) N-VECTORS CONTAINING THE LEFT
C             AND RIGHT BOUNDARIES OF THE OBSERVATION REGION
C             DEFINED BY THE POINTS  X = (X1,...,XN)  SUCH THAT
C                XMIN(I) .LE. X(I) .LE. XMAX(I),  I = 1,...,N
C
      DOUBLE PRECISION X0,XMIN,XMAX
      DOUBLE PRECISION V0,VMIN,VMAX
C
      DIMENSION X0(1),XMIN(1),XMAX(1)
C
      GO TO (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,
     1       160,170,180,190,200,210,220,230,240,250,260,270,280,
     2       290,300,310,320,330,340,350,360,370),NPROB
C
C    1  A FOURTH-ORDER POLYNOMIAL  (N = 1)
C
   10 CONTINUE
      N = 1
      V0 = 1.0D0
      VMIN = -10.D0
      VMAX = 10.D0
      GO TO 800
C
C    2  GOLDSTEIN SIXTH ORDER POLYNOMIAL  (N = 1)
C
   20 CONTINUE
      N = 1
      V0 = 0.D0
      VMIN = -4.D0
      VMAX = 4.D0
      GOTO 800
C
C    3  ONE-DIMENSIONAL PENALIZED SHUBERT FUNCTION  (N = 1)
C
   30 CONTINUE
      N = 1
      V0 = 0.D0
      VMIN = -10.D0
      VMAX = 10.D0
      GOTO 800
C
C    4  A FOURTH ORDER POLYNOMIAL IN THO VARIABLES  (N = 2)
C
   40 CONTINUE
      N = 2
      X0(1) = 1.D0
      X0(2) = 0.D0
      VMIN = -10.D0
      VMAX = 10.D0
      GO TO 700
C
C    5  A FUNCTION WITH A SINGLE ROW OF LOCAL MINIMA  (N = 2)
C
   50 CONTINUE
      N = 2
      X0(1) = -3.D0
      X0(2) = 0.D0
      XMIN(1) = -15.D0
      XMAX(1) = 25.D0
      XMIN(2) = -5.D0
      XMAX(2) = 15.D0
      RETURN
C
C    6  SIX-HUMP CAMEL FUNCTION  (N = 2)
C
   60 CONTINUE
      N = 2
      V0 = 0.D0
      XMIN(1) = -3.D0
      XMAX(1) = 3.D0
      XMIN(2) = -2.D0
      XMAX(2) = 2.D0
      GO TO 900
C
C    7  TWO-DIMENSIONAL PENALIZED SHUBERT FUNCTION, BETA = 0  (N = 2)
C
   70 CONTINUE
C
C    8  TWO-DIMENSIONAL PENALIZED SHUBERT FUNCTION, BETA = 1/2  (N = 2)
C
   80 CONTINUE
C
C    9  TWO-DIMENSIONAL PENALIZED SHUBERT FUNCTION, BETA = 1  (N = 2)
C
   90 CONTINUE
      N = 2
      V0 = 0.D0
      VMIN = -10.D0
      VMAX = 10.D0
      GO TO 800
C
C   10  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10  (N = 2)
C
  100 CONTINUE
C
C   11  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**2  (N = 2)
C
  110 CONTINUE
C
C   12  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**3  (N = 2)
C
  120 CONTINUE
C
C   13  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**4  (N = 2)
C
  130 CONTINUE
C
C   14  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**5  (N = 2)
C
  140 CONTINUE
C
C   15  A FUNCTION WITH THREE ILL-CONDITIONED MINIMA, A=10**6  (N = 2)
C
  150 CONTINUE
      N = 2
      X0(1) = 0.D0
      X0(2) = 0.D0
      XMIN(1) = -10.D0
      XMAX(1) = 10.D0
      XMIN(2) = -100.D0
      XMAX(2) = 100.D0
      RETURN
C
C   16  GOLDSTEIN-PRICE FUNCTION  (N = 2)
C
  160 CONTINUE
      N = 2
      V0 = 1.D0
      VMIN = -2.D0
      VMAX = 2.D0
      GO TO 800
C
C   17  PENALIZED BRANIN FUNCTION  (N = 2)
C
  170 CONTINUE
      N = 2
      X0(1) = 2.5D0
      X0(2) = 7.5D0
      XMIN(1) = -5.D0
      XMAX(1) = 10.D0
      XMIN(2) = 0.D0
      XMAX(2) = 15.D0
      RETURN
C
C   18  PENALIZED SHEKEL FUNCTION, M = 5  (N = 4)
C
  180 CONTINUE
C
C   19  PENALIZED SHEKEL FUNCTION, M = 7  (N = 4)
C
  190 CONTINUE
C
C   20  PENALIZED SHEKEL FUNCTION, M = 10  (N = 4)
C
  200 CONTINUE
      N = 4
      V0 = 9.D0
      VMIN = 0.D0
      VMAX = 10.D0
      GO TO 800
C
C   21  PENALIZED THREE-DIMENSIONAL HARTMAN FUNCTION  (N = 3)
C
  210 CONTINUE
      N = 3
      GO TO 225
C
C   22  PENALIZED SIX-DIMENSIONAL HARTMAN FUNCTION  (N = 6)
C
  220 CONTINUE
      N = 6
  225 CONTINUE
      V0 = 0.5D0
      VMIN = 0.D0
      VMAX = 1.D0
      GO TO 800
C
C   23  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 1  (N = 2)
C
  230 CONTINUE
      N = 2
      GO TO 315
C
C   24  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 1  (N = 3)
C
  240 CONTINUE
      N = 3
      GO TO 315
C
C   25  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 1  (N = 4)
C
  250 CONTINUE
      N = 4
      GO TO 315
C
C   26  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 2  (N = 5)
C
  260 CONTINUE
      N = 5
      GO TO 315
C
C   27  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 2  (N = 8)
C
  270 CONTINUE
      N = 8
      GO TO 315
C
C   28  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 2  (N = 10)
C
  280 CONTINUE
      N = 10
      GO TO 315
C
C   29  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 10   (N = 2)
C
  290 CONTINUE
      N = 2
      GO TO 315
C
C   30  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 10   (N = 3)
C
  300 CONTINUE
      N = 3
      GO TO 315
C
C   31  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 10   (N = 4)
C
  310 CONTINUE
      N = 4
C
  315 CONTINUE
      V0 = 0.D0
      VMIN = -10.D0
      VMAX = 10.D0
      GO TO 800
C
C   32  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 5   (N = 5)
C
  320 CONTINUE
      N = 5
      GO TO 345
C
C   33  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 5   (N = 6)
C
  330 CONTINUE
      N = 6
      GO TO 345
C
C   34  PENALIZED LEVY-MONTALVO FUNCTION, TYPE 3, RANGE = 5  (N = 7)
C
  340 CONTINUE
      N = 7
C
  345 CONTINUE
      V0 = 0.D0
      VMIN = -5.D0
      VMAX = 5.D0
      GO TO 800
C
C   35  A FUNCTION WITH A SINGLE CUSP-SHAPED MINIMUM  (N = 5)
C
  350 CONTINUE
      N = 5
      V0 = 1000.D0
      VMIN = -20000.D0
      VMAX = 10000.D0
      GO TO 800
C
C   36  A FUNCTION WITH A SMALL-ATTRACTION-REGION GLOBAL MINIMUM (N = 2)
C
  360 CONTINUE
      N = 2
      X0(1) = 0.D0
      X0(2) = 100.D0
      VMIN = -1000.D0
      VMAX = 1000.D0
      GO TO 700
C
C   37  A FUNCTION WITH A SMALL-ATTRACTION-REGION GLOBAL MINIMUM (N = 5)
C
  370 CONTINUE
      N = 5
      X0(1) = 0.D0
      X0(2) = 0.D0
      X0(3) = 0.D0
      X0(4) = 0.D0
      X0(5) = 10.D0
      VMIN = -100.D0
      VMAX = 100.D0
      GO TO 700
C
  700 CONTINUE
      DO 710 I = 1,N
         XMIN(I) = VMIN
         XMAX(I) = VMAX
  710 CONTINUE
      RETURN
C
  800 CONTINUE
      DO 810 I = 1,N
         XMIN(I) = VMIN
         XMAX(I) = VMAX
  810 CONTINUE
C
  900 CONTINUE
      DO 910 I = 1,N
         X0(I) = V0
  910 CONTINUE
      RETURN
      END
      SUBROUTINE PTSEG(N,XPFMIN,FPFMIN,FPFMAX,KP,NFEV)
C
C  SAMPLE VERSION OF THE OUTPUT SUBROUTINE  PTSEG
C  (PERFORMS END-OF-SEGMENT OUTPUT)
C  WHICH MUST BY SUPPLIED BY THE USER AND WHICH IS DESCRIBED IN
C  DETAIL WITHIN THE SUBROUTINE  SIGMA.
C
      DOUBLE PRECISION XPFMIN,FPFMIN,FPFMAX
      DIMENSION XPFMIN(N)
C
      WRITE(6,100)KP,NFEV,FPFMIN,FPFMAX
 100  FORMAT(1X,24HOBSERVATION PERIOD  KP= ,I4,
     1    32H,   FUNCTION EVALUATIONS  NFEV= ,I7/
     2    2X,33HFINAL BEST, WORST FUNCT. VALUES  ,
     3    8HFPFMIN= ,G13.5,11H,  FPFMAX= ,G13.5)
      WRITE(6,200)XPFMIN
 200  FORMAT(2X,24HBEST FINAL POINT  XPFMIN/(1X,6G13.5))
C
C
      RETURN
      END
      SUBROUTINE PTRIAL ( N, XOPT, FOPT,
     1                    FTFMIN, FTFMAX, FTFOPT,
     2                    ISTOP, ISTOPT, NFEV, KP, IPRINT )
C
C  SAMPLE VERSION OF THE OUTPUT SUBROUTINE  PTRIAL
C  (PERFORMS END-OF-TRIAL OUTPUT)
C  WHICH MUST BY SUPPLIED BY THE USER AND WHICH IS DESCRIBED IN
C  DETAIL WITHIN THE SUBROUTINE  SIGMA.
C
      DOUBLE PRECISION XOPT,FTFMIN,FTFMAX,FTFOPT,FOPT
      DIMENSION XOPT(N)
C
      WRITE(6,100)
 100  FORMAT(//16H END OF A TRIAL )
      IF(IPRINT.EQ.0)WRITE(6,200)KP,NFEV,FTFMIN,FTFMAX
 200  FORMAT(2X,24HOBSERVATION PERIOD  KP= ,I4,
     1    32H,   FUNCTION EVALUATIONS  NFEV= ,I7/
     2    2X,33HFINAL BEST, WORST FUNCT. VALUES  ,
     3    8HFTFMIN= ,G13.5,11H,  FTFMAX= ,G13.5)
      WRITE(6,300)ISTOP,ISTOPT,FTFOPT,FOPT
 300  FORMAT(/2X,29HTRIAL STOP INDICATOR  ISTOP= ,I2,
     1    34H,   PAST STOPS INDICATOR  ISTOPT= ,I2/
     2    2X,43HEND-OF-TRIAL BEST FUNCTION VALUE   FTFOPT= ,G14.6/
     3    2X,43HBEST CURRENT MINIMUM FUNCTION VALUE  FOPT= ,G14.6)
      WRITE(6,400)XOPT
 400  FORMAT(2X,29HBEST CURRENT MINIMIZER   XOPT/(1X,6G13.5))
C
C
      RETURN
      END
      SUBROUTINE PTKSUC(KSUC)
C
C  SAMPLE VERSION OF THE OUTPUT SUBROUTINE  PTKSUC
C  (PERFORMS END-OF-TRIAL OUTPUT RELATED TO THE COUNT OF
C    SUCCESSFUL TRIALS)
C  WHICH MUST BY SUPPLIED BY THE USER AND WHICH IS DESCRIBED IN
C  DETAIL WITHIN THE SUBROUTINE  SIGMA.
C
      WRITE(6,100)KSUC,KSUC
 100  FORMAT(///1X,
     1    50HTHE CURRENT COUNT  KSUC  OF SUCCESSFUL TRIALS HAS ,
     2    21HREACHED FOR THE FIRST/2X,15HTIME THE VALUE ,I2,
     3    53H (IF THE REQUESTED COUNT  NSUC  OF SUCCESSFUL TRIALS /
     4    2X,25HHAD BEEN GIVEN THE VALUE ,I2,21H THE ALGORITHM WOULD ,
     5    18HHAVE STOPPED HERE)///)
C
C
      RETURN
      END