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 =