PROGRAM PTML
************************************************************************
*  PTML v1.0,  Mar 2003
*
*  This is a parallel implementation of the method:
*
*   "Topographical Multi Level Single Linkage"
*
*  due to A.H.G. Rinnooy Kan and G.T. Timmer
*  described in the following two articles:
*
*  1) A.H.G. Rinnooy Kan and G.T.Timmer,
*     Stochastic global optimization methods
*     Part I: Clustering methods
*     Mathematical Programming, Vol 39, (1987) pp 27-56.
*
*  2) A.H.G. Rinnooy Kan and G.T.Timmer,
*     Stochastic global optimization methods
*     Part II: Multi level methods
*     Mathematical Programming, Vol 39, (1987) pp 57-78.
*
*  with modifications suggested by Aimo Torn and Sami Viitanen
*  described in:
*
*     A. Torn and S. Viitanen,
*     Topographical global optimization using pre-sampled points,
*     Journal of Global Optimization, Vol 5, (1994) pp 267-276.
*
*  and modifications due to M.M. Ali and C. Storey
*  described in:
*
*     M.M. Ali and C. Storey,
*     Topographical multilevel single linkage,
*     Journal of Global Optimization, Vol 5, (1994) pp 349-358.
*
*  Newer versions of this code are available from:
*  http://nrt.cs.uoi.gr/merlin
*
************************************************************************
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	INCLUDE 'mpif.h'
	INCLUDE 'ptml.h'
C  ---------------------------------------------------------------------
C  Workspaces:
C  Total set of points, set of minima.
	DIMENSION S(MXV,MXNT), XMINZ(MXV,MXNM)
	DIMENSION SS(MXV)
C  Corresponding Values of the objective function for the above sets.
	DIMENSION FVAL(MXNT), FMINZ(MXNM)
C  Indices (pointers) for the Graph-Minima and for the local starters.
	DIMENSION IGMIN(MXNT), ISTART(MXNT)
C  Indices for the nearest neighbors, Euclidean distances.
	DIMENSION IREGIO(MXNNN), DIS(MXNT)
	DIMENSION POINT(MXV), XLL(MXV), XRL(MXV), IXAT(MXV), GRAD(MXV) 
	DIMENSION TP(MXV)
	DIMENSION ICODE(4)
	DIMENSION P(MXV,MXPROC), V(MXPROC), NCALL(4,MXPROC)
	DIMENSION IRP(MXPROC), IRV(MXPROC), IRC(MXPROC), IXPROC(MXPROC)
	COMMON /PTMLCN/ NMPERF, NGPERF, NFPROC, NGPROC
C  ---------------------------------------------------------------------
C  Local variables:
	CHARACTER*20 TMPOUT
	CHARACTER*(*) TMPIN
	PARAMETER ( TMPIN = 'tmpin' )
	CHARACTER*10 CIT, CEST, CPROC
	CHARACTER CMD
	DIMENSION ISTAT(MPI_STATUS_SIZE)
	PARAMETER ( ITAG = 0 )
	LOGICAL CEXI
	DATA ICODE  / 1, 1, 1, 1 /
	DATA KPG / 0 /
C  ---------------------------------------------------------------------
C  Default values for the input parameters.
	DATA NOD    / 2 /
	DATA NOF    / 0 /
	DATA ITTHR  / 10 /
	DATA HEAL   / 5 /
	DATA NEARN  / 2 /
	DATA NSAMPL / 100 /
	DATA SIGMA  / 4. /
	DATA COVTHR / 0.995 /
	DATA XDEPS  / 1.D-2 /
	DATA FDEPS  / 1.D-4 /
	DATA GDEPS  / 1.D-3 /
	DATA IREST  / 0 /
	DATA NDUMP  / 100 /
	DATA FINP   / 'in.dat' /
	DATA FOUT   / 'out.dat' /
	DATA FMIN   / 'MINIMA' /
	DATA POIFIL / 'POIMAR.DAT' /
	DATA INDIR  / 'input' /
C  Default values for bounds and fix statuses.
	DATA XLL  / MXV * -1.D0 /
	DATA XRL  / MXV *  1.D0 /
	DATA IXAT / MXV *  1    /
C  ---------------------------------------------------------------------
C
C  Initialize the MPI environment.
	CALL MPI_INIT(IERR)
	IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_INIT',IERR)
	CALL MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)
	IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_COMM_RANK',IERR)
	CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NPROC,IERR)
	IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_COMM_RANK',IERR)
C
C  Check number of processors.
	IF (NPROC.GT.MXPROC) THEN
		IF (IP.EQ.0) WRITE (*,*)
     &	             'PTML: Maximum nuber of processors exceeded.'
		CALL MPI_FINALIZE(IERR)
		STOP
	END IF
	IF (NPROC.LT.2) THEN
		IF (IP.EQ.0) WRITE (*,*)
     &	             'PTML: You need at least two processors.'
		CALL MPI_FINALIZE(IERR)
		STOP
	END IF
C
C  Initialize partial counters (for this processor only).
	NMPERF = 0
	NGPERF = 0
	NFPROC = 0
	NGPROC = 0
C
	IF (IP.NE.0) CALL WORKER(IP,NPROC,POINT,XLL,XRL,IXAT,ICODE,
     &                         XMINZ,FMINZ,S,FVAL,DIS,IREGIO)
C  ---------------------------------------------------------------------
C  Only processor 0 continues below this point.
C  ---------------------------------------------------------------------
C
C  Read the input file.
	CALL FROMFL(NPROC,POINT,XLL,XRL,IXAT)
C
C  Initialize Function, Gradient, Hessian and Jacobian call-counters
C  (for all processors).
	NFTOT = 0
	NGTOT = 0
	NHTOT = 0
	NJTOT = 0
	MXSP = 0
	TNOS = 0
	ITC = 0
	NOM = 0
	LASTP = NPROC-1
C
C  If restore was specified, read the dump file.
	IF (IREST.EQ.1)
     &    CALL RESTOR(ITC,NOM,NFTOT,NGTOT,NHTOT,NJTOT,TNOS,XMINZ,
     &                FMINZ,XLL,XRL,IXAT,MXSP,NLOC)
C
	NOE = NSAMPL
	NNB = NEARN
C
C  Informative printout.
	WRITE (*,*)
	WRITE (*,*) 'Parallel TML running with ...'
	WRITE (*,*) 'Iteration threshold: ', ITTHR
	WRITE (*,*) 'Healing parameter:   ', HEAL
	WRITE (*,*) 'Number of nearest neighbors: ', NEARN
	WRITE (*,*) 'Sample size: ', NSAMPL
	WRITE (*,*) 'Sigma parameter: ', SIGMA
	WRITE (*,*) 'Coverage threshold: ', COVTHR
	WRITE (*,*) 'X distance criterion: ', XDEPS
	WRITE (*,*) 'F distance criterion: ', FDEPS
	WRITE (*,*) 'Gradient convergence criterion: ', GDEPS
	WRITE (*,*) 'Saving every ', NDUMP, ' iterations'
	WRITE (*,*) 'Merlin input file:  ', FINP 
	WRITE (*,*) 'Merlin output file: ', FOUT 
	WRITE (*,*) 'Directory containing input files: ', INDIR
	WRITE (*,*)
C
C  Now every worker has to go to its own directory.
	CMD = 'T'
	DO 10,I=1,LASTP
		CALL MPI_SEND(CMD,1,MPI_CHARACTER,I,ITAG,MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_SEND',IERR)
10	CONTINUE
C
C  Send necessary variables to all worker processors.
	CALL DSEND(NPROC,NOD,NOF,POINT,XLL,XRL,IXAT,ICODE,FINP,FOUT,
     &           NSAMPL,NEARN)
C
C  If continuing a previous run, send the current collection of minima
C  to the rest of the processors.
	IF (IREST.EQ.1) THEN
		DO 20,I=1,NOM
			CALL MSEND(NPROC,NOD,XMINZ(1,I),FMINZ(I))
20		CONTINUE
	END IF
C  ---------------------------------------------------------------------
C  At least one Merlin invocation is needed in order to initialize
C  its internal data structures. 
C
C  Prepare the OPTIMA input file.
	OPEN (1,FILE=TMPIN)
	WRITE (1,*) 'EVALUATE'
	WRITE (1,*) 'RETURN'
	CLOSE (1)
	DO 50,I=1,NOD
		TP(I) = (XLL(I)+XRL(I))/2.
50	CONTINUE
	TMPOUT = FOUT
	CALL OPTIMA(NOD,NOF,TP,VAL,XLL,XRL,IXAT,ICODE,TMPIN,TMPOUT,
     &            GRMS,NF,NG,NH,NJ)
C
C  Update call counters.
	NFTOT = NFTOT + NF
	NGTOT = NGTOT + NG
	NHTOT = NHTOT + NH
	NJTOT = NJTOT + NJ
	NFPROC = NFPROC + NF
	NGPROC = NGPROC + NG
C
C  Remove the OPTIMA input file.
	OPEN (1,FILE=TMPIN)
	CLOSE (1,STATUS='DELETE')
C  ---------------------------------------------------------------------
	NOFV = 0
	ALVOL = 0.
	DO 1,I=1,NOD
		IF (IXAT(I) .EQ. 1) THEN 
			NOFV = NOFV + 1
			ALVOL = ALVOL + LOG(XRL(I)-XLL(I))
		END IF
1	CONTINUE
	CALL CALCBD(XLL,XRL,IXAT,NOD,BOXDI)
	IF (KPG.EQ.0) KPG = LASTP
C
C /\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\
C   Main loop starts below the line: ( 30 CONTINUE )
C /\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\
C
30	CONTINUE
	ITC = ITC + 1
C
C  Calculate the critical radious RC.
	CALL CRITIC(ITC,NOE,NOFV,ALVOL,SIGMA,RC) 
	IF (LP) WRITE (*,*) 'Calculated critical radious: ', RC
C
C  Sample NOE points in array S(,).
	DO 2,I=1,NOE
		DO 3,J=1,NOD
			SS(J) = RANM()*(XRL(J)-XLL(J))+XLL(J)
			IF (IXAT(J).EQ.0) SS(J) = POINT(J)
			S(J,I) = SS(J)
3		CONTINUE
		FVAL(I) = ACSQ(SS,NOD)
		NFTOT = NFTOT+1
		NFPROC = NFPROC+1
2	CONTINUE
C
C  Send the random sample to all worker processors.
	CALL RSEND(NPROC,NOD,S,FVAL,MXV,NOE)
C
C  Add the NOM minima to the set S(,).
	DO 4,I=1,NOM
		NI = NOE+I
		DO 5,J=1,NOD
			S(J,NI) = XMINZ(J,I)
5		CONTINUE
		FVAL(NI) = FMINZ(I)
4	CONTINUE
	NOV = NOE + NOM
C
C  Calculate the Graph minima, with respect to NNB nearest neighbours. 
C  IGMIN[I] = 0 means that the I_th point in S, is not a graph minimum
C  IGMIN[I] = 1 means that the I_th point in S, is a graph minimum
C
	CALL GRAPHM(NOV,S,NNB,IGMIN,NOGM,FVAL,IREGIO,DIS,NPROC)
	IF (LP) WRITE (*,*) 'Calculated graph minima.'
C
C  Determine start points.
C  ISTART(I) is a Flag-Array.
C  If the I_th point in S  is a Start Point then ISTART(I) = 1
C  Otherwise ISTART(I) is set to 0 (ZERO). 
C  Start Points are chosen among the Graph minima.
C  A previously found minimum cannot be a Start Point.
C
	NOSP = 0
	DO 6,I = 1,NOE
		ISTART(I) = IGMIN(I)
C  Consider only points that correspond to graph minima.
		IF (ISTART(I) .EQ. 0) GOTO 70
C  Calculate the distances dis(j) = |r(i)-r(j)| for j=1,NOV 
C  (r(i) = S(.,i))
		CALL DIST(S,NOD,NOV,I,DIS,MXV)
C
C  Check if there is a point within the critical radious RC,
C  that corresponds to a lower function value. 
C  In that case r(i) cannot be a start point.
C
		DO 7,J=1,NOV
			IF (DIS(J).LT.RC) THEN
				IF (FVAL(J) .LT. FVAL(I)) THEN
					ISTART(I) = 0
					GOTO 70
				END IF
			END IF
7		CONTINUE
70		CONTINUE
		NOSP = NOSP + ISTART(I)
6	CONTINUE
C  Start points are determined. (A total of NOSP start points are found)
	IF (LP) WRITE (*,*) 'Start points: ', NOSP
C
	EXX = EXP(-ITC/HEAL)
	TNOS = TNOS + NOSP + (NOGM-NOSP)*(1-EXX)/(1+EXX)
C
	KP = KPG
	NOMST = NOM
	NSE = 0
	DO 200,I=1,NOE
		IF (NOM .EQ. MXNM) THEN
			WRITE (*,*)
     &	            'PTML: MAX number of allowed MINIMA reached.'
			CALL STOPER(NPROC)
		END IF
		IF (ISTART(I).EQ.0) GOTO 200
C
C  Determine the next processor to be used for minimization.
		KP = KP+1
		IF (KP.GT.LASTP) KP = 1
C
C  Send the starting point to a worker processor for further minimization.
		CMD = 'M'
		CALL MPI_SEND(CMD,1,MPI_CHARACTER,KP,ITAG,MPI_COMM_WORLD,
     &	              IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_SEND',IERR)
		CALL MPI_SEND(S(1,I),NOD,MPI_DOUBLE_PRECISION,KP,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_SEND',IERR)
		NSE = NSE+1
C
C  Post non-blocking receives for the results (point, value and
C  call counters).
		CALL MPI_IRECV(P(1,NSE),NOD,MPI_DOUBLE_PRECISION,KP,ITAG,
     &	               MPI_COMM_WORLD,IRP(NSE),IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_IRECV',IERR)
		CALL MPI_IRECV(V(NSE),1,MPI_DOUBLE_PRECISION,KP,ITAG,
     &	               MPI_COMM_WORLD,IRV(NSE),IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_IRECV',IERR)
		CALL MPI_IRECV(NCALL(1,NSE),4,MPI_INTEGER,KP,ITAG,
     &	               MPI_COMM_WORLD,IRC(NSE),IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_IRECV',IERR)
		IXPROC(NSE) = NSE
C
C  A full round over all processors has been completed. Now get the
C  results.
		IF (KP.EQ.KPG) THEN
			NTP = NSE
			DO 250,J=1,NSE
				CALL MPI_WAITANY(NTP,IRP,IDX,ISTAT,IERR)
				IF (IERR.NE.MPI_SUCCESS)
     &			    CALL ERRSTP('MAIN','MPI_WAITANY',IERR)
				CALL MPI_WAIT(IRV(IDX),ISTAT,IERR)
				IF (IERR.NE.MPI_SUCCESS)
     &			    CALL ERRSTP('MAIN','MPI_WAIT',IERR)
				CALL MPI_WAIT(IRC(IDX),ISTAT,IERR)
				IF (IERR.NE.MPI_SUCCESS)
     &			    CALL ERRSTP('MAIN','MPI_WAIT',IERR)
				IXX = IXPROC(IDX)
				NFTOT = NFTOT + NCALL(1,IXX)
				NGTOT = NGTOT + NCALL(2,IXX)
				NHTOT = NHTOT + NCALL(3,IXX)
				NJTOT = NJTOT + NCALL(4,IXX)
				CALL ADDMIN(LP,XDEPS,FDEPS,GDEPS,P(1,IXX),NOD,V(IXX),
     &		            XLL,XRL,IXAT,DIS,GRAD,NOM,FMIN,XMINZ,FMINZ,
     &                        MXV,MXNM,MXNT,NFTOT,NGTOT,NPROC)
				IF (IDX.NE.NTP) THEN
					IRP(IDX) = IRP(NTP)
					IRV(IDX) = IRV(NTP)
					IRC(IDX) = IRC(NTP)
					IXPROC(IDX) = IXPROC(NTP)
					NTP = NTP-1
				END IF
250			CONTINUE
			NSE = 0
		END IF
200	CONTINUE
	NTP = NSE
	DO 260,J=1,NSE
		CALL MPI_WAITANY(NTP,IRP,IDX,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_WAITANY',IERR)
		CALL MPI_WAIT(IRV(IDX),ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_WAIT',IERR)
		CALL MPI_WAIT(IRC(IDX),ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MAIN','MPI_WAIT',IERR)
		IXX = IXPROC(IDX)
		NFTOT = NFTOT + NCALL(1,IXX)
		NGTOT = NGTOT + NCALL(2,IXX)
		NHTOT = NHTOT + NCALL(3,IXX)
		NJTOT = NJTOT + NCALL(4,IXX)
		CALL ADDMIN(LP,XDEPS,FDEPS,GDEPS,P(1,IXX),NOD,V(IXX),
     &	            XLL,XRL,IXAT,DIS,GRAD,NOM,FMIN,XMINZ,FMINZ,
     &                  MXV,MXNM,MXNT,NFTOT,NGTOT,NPROC)
		IF (IDX.NE.NTP) THEN
			IRP(IDX) = IRP(NTP)
			IRV(IDX) = IRV(NTP)
			IRC(IDX) = IRC(NTP)
			IXPROC(IDX) = IXPROC(NTP)
			NTP = NTP-1
		END IF
260	CONTINUE
	KPG = KP
C
C  Send all minima found during this iteration to the rest of the
C  processors.
	DO 210,I=NOMST+1,NOM
		CALL MSEND(NPROC,NOD,XMINZ(1,I),FMINZ(I))
210	CONTINUE
	MXSP = MAX(MXSP,NOSP)
C
	ESTNOM =  NOM*(TNOS-1)/(TNOS-NOM-2) 
	IF (ESTNOM.LT.0.) ESTNOM = NOM+1
	COVER = (TNOS-NOM-1)*(TNOS+NOM)/TNOS/(TNOS-1) 
	IF (COVER.LT.0.) COVER = 0.
C
C  Informative printout.
	CALL I2STR(ITC,CIT,LECIT)
	WRITE (CEST,'(F10.1)') ESTNOM
	DO 40,IEST=1,10
		IF (CEST(IEST:IEST).NE.' ') GOTO 41
40	CONTINUE
41	CONTINUE
	WRITE (*,130) CIT(1:LECIT), COVER, CEST(IEST:)
C
C  Stop the program if the command file exists.
	INQUIRE (FILE=CMDFIL,EXIST=CEXI)
	IF (CEXI) THEN
		CALL DUMPER(ITC,NOM,NFTOT,NGTOT,NHTOT,NJTOT,
     &                  TNOS,XMINZ,FMINZ,XLL,XRL,IXAT,MXSP,NLOC)
		OPEN (10,FILE=CMDFIL)
		CLOSE (10,STATUS='DELETE')
		CALL STOPER(NPROC)
	END IF
C
C  Save the current state every NDUMP iterations.
	IF (NDUMP.NE.0) THEN
		IF (MOD(ITC,NDUMP).EQ.0)
     &	    CALL DUMPER(ITC,NOM,NFTOT,NGTOT,NHTOT,NJTOT,
     &	                TNOS,XMINZ,FMINZ,XLL,XRL,IXAT,MXSP,NLOC)
	END IF
C
	IF (ITC .LT. ITTHR) GOTO 30
	IF ( TNOS .LE. NOM + 2) GOTO 30
	IF ( ESTNOM  .LT.  NOM + 0.5  .AND. COVER .GT. COVTHR ) THEN
		WRITE (*,*)
		WRITE (*,*) 'PTML run completed.'
		IF (NOM.EQ.0) THEN
			WRITE (*,*) 'No local minima were found.'
		ELSE
			CLOSE (1)
			WRITE (*,*) 'Number of local minimizers found: ', NOM
			WRITE (*,*) 'The minimizers are disposed to file: ', FMIN
		END IF
		WRITE (*,*) 'Total number of Function calls: ', NFTOT
		WRITE (*,*) 'Total number of Gradient calls: ', NGTOT
		WRITE (*,*) 'Total number of Jacobian calls: ', NJTOT
		WRITE (*,*) 'Total number of Hessian calls : ', NHTOT
		WRITE (*,*) 'Total number of local optimizations : ', NLOC
		WRITE (*,*) 'Maximum number of starting points: ', MXSP
		CALL DUMPER(ITC,NOM,NFTOT,NGTOT,NHTOT,NJTOT,
     &                  TNOS,XMINZ,FMINZ,XLL,XRL,IXAT,MXSP,NLOC)
		CALL STOPER(NPROC)
	END IF
	GOTO 30
C /\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\
C   End of main loop.
C /\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\
C
120   FORMAT (F8.5,3X,2(F4.1,3X),I1)
130	FORMAT ('Iterations: ',A,5X,'Coverage: ',F9.7,5X,
     &        'Est. minima: ',A)
C
	END
C  ---------------------------------------------------------------------
	SUBROUTINE CHEMIN ( GN, EPS, POINT, NOD, VAL, XLL, XRL, GRAD,
     &                    NOC, NGR )
C  ---------------------------------------------------------------------
C  Description:
C    Calculates the root mean square gradient, for the
C    components that do not "coincide" with one of their margins.
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  ---------------------------------------------------------------------
C  Arguments:
	DIMENSION POINT(NOD), GRAD(NOD), XLL(NOD), XRL(NOD)
C  ---------------------------------------------------------------------
C
	CALL DERMER(POINT,VAL,GRAD,NOC,NGR,1)
	GN = 0.
	DO 10,I = 1,NOD
		RGRAD = GRAD(I)
		YL = ABS( XLL(I) - POINT(I) )
		YR = ABS( XRL(I) - POINT(I) )
		YLX = MAX( ABS( XLL(I) ), ABS( POINT(I) ) )*EPS
		YRX = MAX( ABS( XRL(I) ), ABS( POINT(I) ) )*EPS
		IF (YL .LT. YLX) RGRAD = MIN(0.D0,GRAD(I))
		IF (YR .LT. YRX) RGRAD = MAX(0.D0,GRAD(I))
		GN = GN + RGRAD**2
10	CONTINUE
	GN = SQRT(GN/NOD)
	END
C  ---------------------------------------------------------------------
	SUBROUTINE CRITIC ( IT, NOE, NOFV, ALVOL, SIGMA, RC ) 
C  ---------------------------------------------------------------------
C  Description:
C    Calculates the critical distance RC for the MLSL algorithm.
C    Uses logarithms to avoid intermediate overflows.
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  ---------------------------------------------------------------------
C
      PI = ACOS(-1.D0)
C
C  The following IF structure, calculates the
C  Logarithm of the Gamma function:  
C  algam = Gamma(1+nofv/2)
C  for two cases: NOFV is even, and NOFV is odd.
C
	IF ( (NOFV/2)*2 .EQ. NOFV ) THEN
		ALGAM = 0.
		DO 10,I=2,NOFV/2
			ALGAM = ALGAM + LOG(DBLE(I))
10		CONTINUE
	ELSE
		ALGAM = LOG(SQRT(PI)) - (NOFV+1)/2*LOG(2.D0)
		DO 20,I=1,NOFV,2
			ALGAM = ALGAM + LOG(DBLE(I))
20		CONTINUE
	END IF
C 
	DIN = DBLE(IT*NOE)
	ALRC = -LOG(SQRT(PI))+LOG( SIGMA*LOG(DIN)/DIN )/NOFV
	ALRC = ALRC + ALVOL/NOFV + ALGAM/NOFV
	RC = EXP(ALRC)
	END
C  ---------------------------------------------------------------------
	SUBROUTINE MINEL ( DIS, NOE, JN, VMIN )
C  ---------------------------------------------------------------------
C  Description:
C    Given an array (DIS) with NOE elements, subroutine MINEL 
C    calculates the minimum element (VMIN) as well as its index (JN).
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  ---------------------------------------------------------------------
C  Arguments:
	DIMENSION DIS(*)
C  ---------------------------------------------------------------------
C
	JN = 1
	VMIN = DIS(JN)
	DO 10,I=2,NOE
		IF (DIS(I) .LT. VMIN) THEN
			VMIN = DIS(I)
			JN = I
		END IF
10	CONTINUE
	END
C  ---------------------------------------------------------------------
	SUBROUTINE MAXEL ( DIS, NOE, JX, VMAX )
C  ---------------------------------------------------------------------
C  Description:
C    Given an array (DIS) with NOE elements, subroutine MAXEL 
C    calculates the maximum element (VMAX) as well as its index (JX).
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  ---------------------------------------------------------------------
C  Arguments:
	DIMENSION DIS(*)
C  ---------------------------------------------------------------------
C
	JX = 1
	VMAX = DIS(JX)
	DO 10,I=2,NOE
		IF (DIS(I) .GT. VMAX) THEN
			VMAX = DIS(I)
			JX = I
		END IF
10	CONTINUE
	END
C  ---------------------------------------------------------------------
	SUBROUTINE DIST ( S, NOD, NOV, I, DIS, MXV )
C  ---------------------------------------------------------------------
C
C  Description:
C    Given:
C    1) An array S containing NOV vectors, with NOD components each,
C    2) An Integer I, in [1,NOV],
C    It calculates the Euclidean distances
C    between the Ith and the Jth vector in S for all J=1,2,...,NOV
C    and returns them in array DIS
C
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  ---------------------------------------------------------------------
C  Arguments:
	DIMENSION S(MXV,*), DIS(*)
C  ---------------------------------------------------------------------
C
	DO 10,J=1,NOV
		D = 0.D0
		DO 20,K=1,NOD
			D = D + (S(K,J)-S(K,I))**2
20		CONTINUE
		DIS(J) = SQRT(D)
10	CONTINUE
	END
C  ---------------------------------------------------------------------
	SUBROUTINE NEARNE ( DIS, NOV, I, IREGIO, NNB, BOXDI )
C  ---------------------------------------------------------------------
C
C  Description:
C    Calculates the NNB nearest neighbours of point with index I
C    and returns their indices in IREGIO(1:NNB)
C    DIS(1:NOV) contains the distances |X(I)-X(J)| for J=1,2,...,NOV
C    The array DIS is destroyed.
C
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  ---------------------------------------------------------------------
C  Arguments:
	DIMENSION DIS(*), IREGIO(*)
C  ---------------------------------------------------------------------
C
	DIS(I) = BOXDI
	DO 10,J=1,NNB
		CALL MINEL(DIS,NOV,JN,VN)
		IREGIO(J) = JN
		DIS(JN) = BOXDI
10	CONTINUE
	END
C  ---------------------------------------------------------------------
	SUBROUTINE WORKER ( IP, NPROC, POINT, XLL, XRL, IXAT, ICODE,
     &                    XMINZ, FMINZ, S, FVAL, DIS, IREGIO )
C  ---------------------------------------------------------------------
C  Description:
C    This is the main subroutine executed by all worker processors.
C    It accepts commands from the main processor and acts accordingly.
C    The following commands are recognized:
C      T  Change directory.
C      D  Receive initial data.
C      S  Stop this worker process.
C      M  Mimimize a point.
C      A  Add to our collection of minima.
C      R  Receive random sample.
C      E  Enhance the random sample S with the set of minima.
C      G  Calculate graph minimum.
C      
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	INCLUDE 'mpif.h'
	INCLUDE 'ptml.h'
	DIMENSION POINT(*), XLL(*), XRL(*), IXAT(*), ICODE(*)
	DIMENSION ICOUNT(4), ISTAT(MPI_STATUS_SIZE), IREQ(3)
	DIMENSION ISTAT2(MPI_STATUS_SIZE,3)
	DIMENSION XMINZ(MXV,MXNM), FMINZ(MXNM)
	DIMENSION S(MXV,MXNT), FVAL(MXNT)
	DIMENSION IREGIO(MXNNN), DIS(MXNT)
	CHARACTER*10 STR
	CHARACTER CMD
	PARAMETER ( ITAG = 0 )
	PARAMETER ( IPROC = 0 )
	INTEGER CHDIR
	COMMON /PTMLCN/ NMPERF, NGPERF, NFPROC, NGPROC
C
	NOM = 0
C
100	CONTINUE
	CALL MPI_RECV(CMD,1,MPI_CHARACTER,IPROC,ITAG,MPI_COMM_WORLD,
     &              ISTAT,IERR)
	IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
C  ---------------------------------------------------------------------
	IF (CMD.EQ.'T') THEN
C  ---------------------------------------------------------------------
C  Change directory.
		CALL I2STR(IP,STR,LES)
		II = CHDIR('tmp.'//STR(1:LES))
		IF (II.NE.0) THEN
			WRITE (*,*) 'PTML was unable to change directory for ',
     &		            'processor ', IP
			STOP 'ERR'
		END IF
C  ---------------------------------------------------------------------
	ELSE IF (CMD.EQ.'D') THEN
C  ---------------------------------------------------------------------
C  Receive initial data.
		CALL MPI_RECV(NOD,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL MPI_RECV(NOF,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL MPI_RECV(XLL,NOD,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL MPI_RECV(XRL,NOD,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL MPI_RECV(IXAT,NOD,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL MPI_RECV(FINP,LEN(FINP),MPI_CHARACTER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL MPI_RECV(FOUT,LEN(FOUT),MPI_CHARACTER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL MPI_RECV(NOE,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL MPI_RECV(NNB,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL CALCBD(XLL,XRL,IXAT,NOD,BOXDI)
C  ---------------------------------------------------------------------
	ELSE IF (CMD.EQ.'S') THEN
C  ---------------------------------------------------------------------
C  Stop this worker process.
		CALL MPI_SEND(NMPERF,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_SEND',IERR)
		CALL MPI_SEND(NGPERF,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_SEND',IERR)
		CALL MPI_SEND(NFPROC,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_SEND',IERR)
		CALL MPI_SEND(NGPROC,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_SEND',IERR)
		CALL MPI_FINALIZE(IERR)
		STOP
C  ---------------------------------------------------------------------
	ELSE IF (CMD.EQ.'M') THEN
C  ---------------------------------------------------------------------
C  Mimimize a point.
		CALL MPI_RECV(POINT,NOD,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL OPTIMA(NOD,NOF,POINT,VAL,XLL,XRL,IXAT,ICODE,FINP,FOUT,
     &	            GRMS,NF,NG,NH,NJ)
		NMPERF = NMPERF+1
		NFPROC = NFPROC+NF
		NGPROC = NGPROC+NG
		CALL MPI_ISEND(POINT,NOD,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	               MPI_COMM_WORLD,IREQ(1),IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_ISEND',IERR)
		CALL MPI_ISEND(VAL,1,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	               MPI_COMM_WORLD,IREQ(2),IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_ISEND',IERR)
		ICOUNT(1) = NF
		ICOUNT(2) = NG
		ICOUNT(3) = NH
		ICOUNT(4) = NJ
		CALL MPI_ISEND(ICOUNT,4,MPI_INTEGER,IPROC,ITAG,
     &	               MPI_COMM_WORLD,IREQ(3),IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_ISEND',IERR)
		CALL MPI_WAITALL(3,IREQ,ISTAT2,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_WAITALL',IERR)
C  ---------------------------------------------------------------------
	ELSE IF (CMD.EQ.'A') THEN
C  ---------------------------------------------------------------------
C  Add to our collection of minima.
		CALL MPI_RECV(POINT,NOD,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL MPI_RECV(VAL,1,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		NOM = NOM+1
		DO 10,J=1,NOD
			XMINZ(J,NOM) = POINT(J)
10		CONTINUE
		FMINZ(NOM) = VAL
C  ---------------------------------------------------------------------
	ELSE IF (CMD.EQ.'R') THEN
C  ---------------------------------------------------------------------
C  Receive random sample.
		DO 20,J=1,NOE
			CALL MPI_RECV(S(1,J),NOD,MPI_DOUBLE_PRECISION,IPROC,
     &		              ITAG,MPI_COMM_WORLD,ISTAT,IERR)
			IF (IERR.NE.MPI_SUCCESS)
     &		    CALL ERRSTP('WORKER','MPI_RECV',IERR)
			CALL MPI_RECV(FVAL(J),1,MPI_DOUBLE_PRECISION,IPROC,
     &		              ITAG,MPI_COMM_WORLD,ISTAT,IERR)
			IF (IERR.NE.MPI_SUCCESS)
     &		    CALL ERRSTP('WORKER','MPI_RECV',IERR)
20		CONTINUE
C  ---------------------------------------------------------------------
	ELSE IF (CMD.EQ.'E') THEN
C  ---------------------------------------------------------------------
C  Enhance the random sample S with the set of minima.
		DO 4,I=1,NOM
			NI = NOE+I
			DO 5,J=1,NOD
				S(J,NI) = XMINZ(J,I)
5			CONTINUE
			FVAL(NI) = FMINZ(I)
4		CONTINUE
		NOV = NOE + NOM
C  ---------------------------------------------------------------------
	ELSE IF (CMD.EQ.'G') THEN
C  ---------------------------------------------------------------------
C  Calculate graph minimum.
		NGPERF = NGPERF+1
		CALL MPI_RECV(IM,1,MPI_INTEGER,IPROC,ITAG,MPI_COMM_WORLD,
     &		        ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_RECV',IERR)
		CALL GMIN1(IM,S,NOD,NOV,NNB,FVAL,IREGIO,DIS,MXV,BOXDI,IGM)
		CALL MPI_SEND(IGM,1,MPI_INTEGER,IPROC,ITAG,MPI_COMM_WORLD,
     &	              IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('WORKER','MPI_SEND',IERR)
C  ---------------------------------------------------------------------
	ELSE
C  ---------------------------------------------------------------------
		WRITE (*,*) 'PTML: Processor ', IP, ' received an ',
     &	            'incorrect command:', CMD
		STOP 'ERR'
	END IF
	GOTO 100
C
	END
C  ---------------------------------------------------------------------
	SUBROUTINE DSEND ( NPROC, NOD, NOF, POINT, XLL, XRL, IXAT, ICODE,
     &                   FINP, FOUT, NOE, NNB )
C  ---------------------------------------------------------------------
C  Description:
C    Sends initial data to all worker processors.
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	INCLUDE 'mpif.h'
	DIMENSION POINT(NOD), XLL(NOD), XRL(NOD), IXAT(NOD), ICODE(NOD)
	CHARACTER*(*) FINP, FOUT
	PARAMETER ( ITAG = 0 )
	CHARACTER CMD
C
	CMD = 'D'
	LASTP = NPROC-1
	DO 10,IPROC=1,LASTP
		CALL MPI_SEND(CMD,1,MPI_CHARACTER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
		CALL MPI_SEND(NOD,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
		CALL MPI_SEND(NOF,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
		CALL MPI_SEND(XLL,NOD,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
		CALL MPI_SEND(XRL,NOD,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
		CALL MPI_SEND(IXAT,NOD,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
		CALL MPI_SEND(FINP,LEN(FINP),MPI_CHARACTER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
		CALL MPI_SEND(FOUT,LEN(FOUT),MPI_CHARACTER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
		CALL MPI_SEND(NOE,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
		CALL MPI_SEND(NNB,1,MPI_INTEGER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('DSEND','MPI_SEND',IERR)
10	CONTINUE
	END
C  ---------------------------------------------------------------------
	SUBROUTINE ADDMIN ( LP, XDEPS, FDEPS, GDEPS, POINT, NOD, VAL, XLL,
     &                    XRL, IXAT, DIS, GRAD, NOM, FMIN, XMINZ, FMINZ,
     &                    MXV, MXNM, MXNT, NFTOT, NGTOT, NPROC )
C  ---------------------------------------------------------------------
C  Description:
C    Adds a newly found minimum to our collection.
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	DIMENSION POINT(*), XLL(*), XRL(*), GRAD(*), IXAT(*)
	DIMENSION XMINZ(MXV,MXNM), FMINZ(MXNM), DIS(MXNT)
	LOGICAL LP
	CHARACTER*(*) FMIN
	COMMON /PTMLCN/ NMPERF, NGPERF, NFPROC, NGPROC
c
C  Calculate the effective root mean square of the gradient.
C  Components "touching" the margins are not taken into account.
	CALL CHEMIN(GN,XDEPS,POINT,NOD,VAL,XLL,XRL,GRAD,NFC,NGC)
	NFPROC = NFPROC+NFC
	NGPROC = NGPROC+NGC
	NFTOT = NFTOT+NFC
	NGTOT = NGTOT+NGC
	IF (LP) WRITE (*,*) 'GN = ',GN
	IF (GN .GT. GDEPS) return
	IF (NOM .EQ. 0) THEN
		OPEN (1,FILE=FMIN)
		CLOSE (1,STATUS='DELETE')
		OPEN (1,FILE=FMIN)
		WRITE (1,*) NOD
		NOM = 1
C  Always accept the first minimum.
		DO 12,J=1,NOD
			XMINZ(J,NOM) = POINT(J)
12		CONTINUE
		FMINZ(NOM) = VAL
		CALL WRITER(POINT,VAL,1,1,0)
		WRITE (*,*) '... The first minimum is found ...',
     &            '   F: ', NFTOT, '   G: ', NGTOT
		RETURN
	END IF
C Check if this is a new minimum. Follow the steps: 
C  I) At first accept it as a new minimum. 
C II) Then check if it is true and if not discard it.
	NOM = NOM + 1
	DO 13,J=1,NOD
		XMINZ(J,NOM) = POINT(J)
13	CONTINUE
	FMINZ(NOM) = VAL
	CALL DIST(XMINZ,NOD,NOM,NOM,DIS,MXV)
	CALL MINEL(DIS,NOM-1,JN,VMIN)
	IF ( ABS(FMINZ(JN)-VAL) .LE. 
     &     FDEPS*MAX(1.,ABS(VAL),ABS(FMINZ(JN))) .AND.
     &     VMIN .LE. XDEPS ) THEN
		NOM = NOM - 1
		IF (LP) WRITE (*,*) 'NOT a new MINIMUM'
		RETURN
	END IF
	CALL WRITER(POINT,VAL,1,1,0)
	WRITE (*,*) '... Number of minima found: ',NOM,
     &            '   F: ', NFTOT, '   G: ', NGTOT 
C
	END
C  ---------------------------------------------------------------------
	SUBROUTINE STOPER ( NPROC )
C  ---------------------------------------------------------------------
C  Description:
C    Stops the program by sending the "S" (stop) command to all worker
C    processors.
C  ---------------------------------------------------------------------
	INCLUDE 'mpif.h'
C  ---------------------------------------------------------------------
C  Local variables:
	CHARACTER CMD
	PARAMETER ( ITAG = 0 )
	COMMON /PTMLCN/ NMPERF, NGPERF, NFPROC, NGPROC
	DIMENSION ISTAT(MPI_STATUS_SIZE)
C  ---------------------------------------------------------------------
C
	LASTP = NPROC-1
	CMD = 'S'
C
C  Send the "stop" command to all workers.
	WRITE (*,*)
	WRITE (*,*) 'Processor Utilization'
	WRITE (*,30) 
	WRITE (*,20) 0, NMPERF, NGPERF, NFPROC, NGPROC
	DO 10,I=1,LASTP
		CALL MPI_SEND(CMD,1,MPI_CHARACTER,I,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('STOPER','MPI_SEND',IERR)
		CALL MPI_RECV(NM1,1,MPI_INTEGER,I,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('STOPER','MPI_RECV',IERR)
		CALL MPI_RECV(NG1,1,MPI_INTEGER,I,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('STOPER','MPI_RECV',IERR)
		CALL MPI_RECV(NFP1,1,MPI_INTEGER,I,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('STOPER','MPI_RECV',IERR)
		CALL MPI_RECV(NGP1,1,MPI_INTEGER,I,ITAG,
     &	              MPI_COMM_WORLD,ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('STOPER','MPI_RECV',IERR)
		WRITE (*,20) I, NM1, NG1, NFP1, NGP1
10	CONTINUE
20	FORMAT (1X,I4,4(2X,I10))
30	FORMAT (1X,'Proc','       Minim','       Graph','      Fcalls',
     &        '      Gcalls')
C
C  End the program.
	CALL MPI_FINALIZE(IERR)
	STOP
C
	END
C  ---------------------------------------------------------------------
	SUBROUTINE MSEND ( NPROC, NOD, POINT, VAL )
C  ---------------------------------------------------------------------
C  Description:
C    Sends a minimum to all processors.
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	INCLUDE 'mpif.h'
	DIMENSION POINT(NOD)
	PARAMETER ( ITAG = 0 )
	CHARACTER CMD
C
	CMD = 'A'
	LASTP = NPROC-1
	DO 10,IPROC=1,LASTP
		CALL MPI_SEND(CMD,1,MPI_CHARACTER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MSEND','MPI_SEND',IERR)
		CALL MPI_SEND(POINT,NOD,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MSEND','MPI_SEND',IERR)
		CALL MPI_SEND(VAL,1,MPI_DOUBLE_PRECISION,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('MSEND','MPI_SEND',IERR)
10	CONTINUE
	END
C  ---------------------------------------------------------------------
	SUBROUTINE RSEND ( NPROC, NOD, S, FVAL, MXV, NOE )
C  ---------------------------------------------------------------------
C  Description:
C    Sends the random sample to all processors.
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	INCLUDE 'mpif.h'
	DIMENSION S(MXV,*), FVAL(*)
	PARAMETER ( ITAG = 0 )
	CHARACTER CMD
C
	CMD = 'R'
	LASTP = NPROC-1
	DO 10,IPROC=1,LASTP
		CALL MPI_SEND(CMD,1,MPI_CHARACTER,IPROC,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('RSEND','MPI_SEND',IERR)
		DO 20,J=1,NOE
			CALL MPI_SEND(S(1,J),NOD,MPI_DOUBLE_PRECISION,
     &		              IPROC,ITAG,MPI_COMM_WORLD,IERR)
			IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('RSEND','MPI_SEND',IERR)
			CALL MPI_SEND(FVAL(J),1,MPI_DOUBLE_PRECISION,
     &		              IPROC,ITAG,MPI_COMM_WORLD,IERR)
			IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('RSEND','MPI_SEND',IERR)
20		CONTINUE
10	CONTINUE
	END
C  ---------------------------------------------------------------------
	SUBROUTINE GRAPHM ( NOV, S, NNB, IGMIN, NOGM, FVAL, IREGIO, DIS,
     &                    NPROC )
C  ---------------------------------------------------------------------
C
C  Description:
C    Calculates the graph minima, for NNB number of nearest neighbours.
C    The graph is the set of points contained in array S(MXV,*)
C    It returns the array  IGMIN(*)
C    IGMIN(I) = 0 means that the I_th point in S, is not a graph minimum
C    IGMIN(I) = 1 means that the I_th point in S, is a graph minimum
C    NOGM is output. Corresponds to the number of graph minima
C    FVAL(I) is the value corresponding to the I_th point in S
C
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	INCLUDE 'mpif.h'
	INCLUDE 'ptml.h'
C  ---------------------------------------------------------------------
C  Arguments:
      DIMENSION S(MXV,*), IGMIN(*), IREGIO(*), FVAL(*), DIS(*)
C  ---------------------------------------------------------------------
	PARAMETER ( ITAG = 0 )
	CHARACTER CMD
	DIMENSION IV(MXPROC)
	DIMENSION ISTAT(MPI_STATUS_SIZE)
	COMMON /PTMLCN/ NMPERF, NGPERF, NFPROC, NGPROC
	DATA KPG / 0 /
	SAVE KPG
C
	LASTP = NPROC-1
	IF (KPG.EQ.0) KPG = LASTP
C
C  Add the minima to the random sample.
	CMD = 'E'
	DO 10,I=1,LASTP
		CALL MPI_SEND(CMD,1,MPI_CHARACTER,I,ITAG,MPI_COMM_WORLD,
     &	              IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('GRAPHM','MPI_SEND',IERR)
10	CONTINUE
C
	NOGM = 0
	KP = KPG
	NSE = 0
	CMD = 'G'
	DO 20,I=1,NOV
		KP = KP+1
		IF (KP.GT.LASTP) KP = 1
C
		CALL MPI_SEND(CMD,1,MPI_CHARACTER,KP,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('GRAPHM','MPI_SEND',IERR)
		CALL MPI_SEND(I,1,MPI_INTEGER,KP,ITAG,
     &	              MPI_COMM_WORLD,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('GRAPHM','MPI_SEND',IERR)
		NSE = NSE+1
		IV(NSE) = I
C
		IF (KP.EQ.KPG) THEN
			KR = KPG
			DO 100,J=1,NSE
				KR = KR+1
				IF (KR.GT.LASTP) KR = 1
				CALL MPI_RECV(IGM,1,MPI_INTEGER,KR,ITAG,
     &			              MPI_COMM_WORLD,ISTAT,IERR)
				IF (IERR.NE.MPI_SUCCESS)
     &			    CALL ERRSTP('GRAPHM','MPI_SEND',IERR)
				IGMIN(IV(J)) = IGM
				NOGM = NOGM + IGM
100			CONTINUE
			NSE = 0
		END IF
20	CONTINUE
C
	KR = KPG
	DO 110,J=1,NSE
		KR = KR+1
		IF (KR.GT.LASTP) KR = 1
		CALL MPI_RECV(IGM,1,MPI_INTEGER,KR,ITAG,MPI_COMM_WORLD,
     &	              ISTAT,IERR)
		IF (IERR.NE.MPI_SUCCESS) CALL ERRSTP('GRAPHM','MPI_RECV',IERR)
		IGMIN(IV(J)) = IGM
		NOGM = NOGM + IGM
110	CONTINUE
	KPG = KP
      END
C  ---------------------------------------------------------------------
	SUBROUTINE CALCBD ( XLL, XRL, IXAT, NOD, BOXDI )
C  ---------------------------------------------------------------------
C  Description:
C    Calculates the size of the box as specified by the margins.
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  Arguments:
      DIMENSION XLL(NOD), XRL(NOD), IXAT(NOD)
C  ---------------------------------------------------------------------
C
	BOXDI = 0.
	DO 1,I=1,NOD
		IF (IXAT(I) .EQ. 1)
     &	    BOXDI = BOXDI + (XRL(I)-XLL(I))**2
1	CONTINUE
	BOXDI = 2*SQRT(BOXDI)
	END
C  ---------------------------------------------------------------------
	SUBROUTINE GMIN1 ( IM, S, NOD, NOV, NNB, FVAL, IREGIO,
     &                   DIS, MXV, BOXDI, IGM )
C  ---------------------------------------------------------------------
C  Description:
C    Decides whether a specified point is a graph minimum.
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  ---------------------------------------------------------------------
C  Arguments:
      DIMENSION S(MXV,*), IREGIO(*), FVAL(*), DIS(*)
C  ---------------------------------------------------------------------
C
	CALL DIST(S,NOD,NOV,IM,DIS,MXV)
	CALL NEARNE(DIS,NOV,IM,IREGIO,NNB,BOXDI)
	DO 10,J=1,NNB
		DIS(J) = FVAL( IREGIO(J) )
10	CONTINUE
	CALL MINEL(DIS,NNB,JM,FM)
	IF ( FVAL(IM) .LE. FM ) THEN
		IGM = 1
	ELSE
		IGM = 0
	END IF
	END
C  ---------------------------------------------------------------------
	SUBROUTINE FROMFL ( NPROC, POINT, XLL, XRL, IXAT )
C  ---------------------------------------------------------------------
C  Description:
C    Reads the input file and checks whether the supplied values make
C    sense.
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	INCLUDE 'ptml.h'
	DIMENSION POINT(MXV), XLL(MXV), XRL(MXV), IXAT(MXV)
	LOGICAL ERR
C
	LOGICAL EXI1, EXI2, MEXI
C
C  Check if the input file exists. 
C  If it does, read the parameters from there.
C  Otherwise use the defaults & create a sample input file,
C  containing the default values, as an aid to the user.
C
	INQUIRE (FILE=INFIL,EXIST=EXI1)
	IF (EXI1) THEN
C  Open and read the data file.
		OPEN (1,FILE=INFIL)
		READ (1,*) NOD
		READ (1,*) NOF
		READ (1,*) ITTHR
		READ (1,*) HEAL
		READ (1,*) NEARN
		READ (1,*) NSAMPL
		READ (1,*) SIGMA
		READ (1,*) COVTHR
		READ (1,*) XDEPS
		READ (1,*) FDEPS
		READ (1,*) GDEPS
		READ (1,*) IREST
		READ (1,*) NDUMP
		READ (1,'(A)') FINP 
		READ (1,'(A)') FOUT 
		READ (1,'(A)') FMIN
		READ (1,'(A)') POIFIL
		READ (1,'(A)') INDIR
		CLOSE (1)
	ELSE
C  Create a sample input file with default values.
		OPEN (1,FILE=INFIL,STATUS='NEW')
		WRITE (1,'(1X,I8,T25,A)')     NOD,    'NOD'
		WRITE (1,'(1X,I8,T25,A)')     NOF,    'NOF'
		WRITE (1,'(1X,I8,T25,A)')     ITTHR,  'ITTHR'
		WRITE (1,'(1X,G21.14,T25,A)') HEAL,   'HEAL'
		WRITE (1,'(1X,I8,T25,A)')     NEARN,  'NEARN'
		WRITE (1,'(1X,I8,T25,A)')     NSAMPL, 'NSAMPL'
		WRITE (1,'(1X,G21.14,T25,A)') SIGMA,  'SIGMA'
		WRITE (1,'(1X,G21.14,T25,A)') COVTHR, 'COVTHR'
		WRITE (1,'(1X,G21.14,T25,A)') XDEPS,  'XDEPS'
		WRITE (1,'(1X,G21.14,T25,A)') FDEPS,  'FDEPS'
		WRITE (1,'(1X,G21.14,T25,A)') GDEPS,  'GDEPS'
		WRITE (1,'(1X,I8,T25,A)')     IREST,  'IREST'
		WRITE (1,'(1X,I8,T25,A)')     NDUMP,  'NDUMP'
		WRITE (1,'(A,T25,A)')         FINP ,  'FINP'
		WRITE (1,'(A,T25,A)')         FOUT ,  'FOUT'
		WRITE (1,'(A,T25,A)')         FMIN ,  'FMIN'
		WRITE (1,'(A,T25,A)')         POIFIL, 'POIFIL'
		WRITE (1,'(A,T25,A)')         INDIR,  'INDIR'
		CLOSE (1)
		WRITE (*,*) 'PTML: File ',INFIL,' is created with ',
     &	            'default values.'
	END IF
C
C  Make sure NOD is correct.
	IF (NOD.LT.1) THEN
		WRITE (*,*) 'PTML: Incorrect number of parameters ',
     &	            '(must be positive).'
		CALL STOPER(NPROC)
	END IF
	IF (NOD.GT.MXV) THEN
		WRITE (*,*) 'PTML: The number of parameters is set to: ',
     &	            NOD,' > ', MXV,' which is the maximum allowed.'
		CALL STOPER(NPROC)
	END IF
C
C  Perform initialization of required arrays. (POINT,XLL,XRL & IXAT)
C 
	INQUIRE (FILE=POIFIL,EXIST=EXI2)
	IF (EXI2) THEN
C  Read parameter bounds.
		OPEN (1,FILE=POIFIL,STATUS='OLD')
		DO 20,I=1,NOD
			READ (1,*) POINT(I), XLL(I), XRL(I), IXAT(I)
20		CONTINUE
		CLOSE (1)
	ELSE
C  Create a sample parameter file with default values.
		OPEN (1,FILE=POIFIL,STATUS='NEW')
		DO 9,I=1,NOD
			POINT(I) = 2*RANM()-1.
			WRITE (1,120) POINT(I), XLL(I), XRL(I), IXAT(I)
9		CONTINUE
		CLOSE (1)
		WRITE (*,*) 'PTML: File ',POIFIL,' is created with ',
     &	            'default values.'
	END IF
C
C Check if the input is meaningful.
	ERR = .FALSE.
	IF (NOF.LT.0) THEN
		WRITE (*,*) 'PTML: Incorrect number of function terms ',
     &	            '(must be zero or positive).'
		CALL STOPER(NPROC)
	END IF
C
	IF (ITTHR.LE.0) THEN
		WRITE (*,*) 'PTML: Incorrect iteration threshold ',ITTHR,
     &	            ' (must be positive).'
		ERR = .TRUE.
	END IF
C
	IF (NEARN.LE.0) THEN
		WRITE (*,*) 'PTML: Incorrect NEARN ', NEARN,
     &	            ' (must be positive).'
		ERR = .TRUE.
	END IF
	IF (NEARN.GT.MXNNN) THEN
		WRITE (*,*) 'PTML: NEARN is set to: ',NEARN,' > ',
     &	            MXNNN,' which is the maximum allowed.'
		ERR = .TRUE.
	END IF
C
	IF (NSAMPL.LE.0) THEN
		WRITE (*,*) 'PTML: Incorrect NSAMPL ', NSAMPL,
     &	            ' (must be positive).'
		ERR = .TRUE.
	END IF
	IF (NSAMPL.GT.MXNS) THEN
		WRITE (*,*) 'PTML: NSAMPL is set to: ', NSAMPL,' > ',
     &	            MXNS,' which is the maximum allowed.'
		ERR = .TRUE.
	END IF
C
	IF (SIGMA.LE.0.) THEN
		WRITE (*,*) 'PTML: Incorrect SIGMA ', SIGMA
		ERR = .TRUE.
	END IF
C
	IF (COVTHR.LE.0. .OR. COVTHR.GT.1.) THEN
		WRITE (*,*) 'PTML: Incorrect COVTHR ', COVTHR
		ERR = .TRUE.
	END IF
C
	IF (XDEPS.LE.0. .OR. XDEPS.GE.1.) THEN
		WRITE (*,*) 'PTML: Incorrect XDEPS ', XDEPS
		ERR = .TRUE.
	END IF
C
	IF (FDEPS.LE.0. .OR. FDEPS.GE.1.) THEN
		WRITE (*,*) 'PTML: Incorrect FDEPS ', FDEPS
		ERR = .TRUE.
	END IF
C
	IF (GDEPS.LE.0. .OR. GDEPS.GE.1.) THEN
		WRITE (*,*) 'PTML: Incorrect GDEPS ', GDEPS
		ERR = .TRUE.
	END IF
C
	IF (IREST.NE.0 .AND. IREST.NE.1) THEN
		WRITE (*,*) 'PTML: Incorrect IREST ', IREST
		ERR = .TRUE.
	END IF
C
	IF (NDUMP.LT.0) THEN
		WRITE (*,*) 'PTML: Incorrect NDUMP ', NDUMP
		ERR = .TRUE.
	END IF
C
	IF (FINP.EQ.' ') THEN
		WRITE (*,*) 'PTML: Merlin input file not specified.'
		ERR = .TRUE.
	END IF
C
	IF (FOUT.EQ.' ') THEN
		WRITE (*,*) 'PTML: Merlin output file not specified.'
		ERR = .TRUE.
	END IF
C
	INQUIRE (FILE=FINP,EXIST=MEXI)
	IF (.NOT.MEXI) THEN
		LEFI = LENGTH(FINP)
		WRITE (*,*) 'PTML: The Merlin input file ',FINP(1:LEFI),
     &	            ' does not exist.'
	END IF
C
C  Missing input files ?
	IF (.NOT. (EXI1.AND.EXI2)) THEN
		WRITE (*,*) 'PTML: Please edit the iput files ',
     &	            'to provide the desired input.'
		CALL STOPER(NPROC)
	END IF
C
C  Check if parameters, bounds and fix statuses are meaningful.
	DO 100,I=1,NOD
		IF (XLL(I).GE.XRL(I)) THEN
			WRITE (*,*) 'PTML: Incorrect bounds for parameter ', I
			ERR = .TRUE.
		ELSE IF (POINT(I).LT.XLL(I) .OR. POINT(I).GT.XRL(I)) THEN
			WRITE (*,*) 'PTML: Out of bounds parameter ', I
			ERR = .TRUE.
		END IF
		IF (IXAT(I).NE.0 .AND. IXAT(I).NE.1) THEN
			WRITE (*,*) 'PTML: Incorrect fix status for variable ', I
			ERR = .TRUE.
		END IF
100	CONTINUE
	IF (ERR) CALL STOPER(NPROC)
C
120   FORMAT (F8.5,3X,2(F4.1,3X),I1)
C
	END
C  ---------------------------------------------------------------------
	SUBROUTINE DUMPER ( ITC, NOM, NFTOT, NGTOT, NHTOT, NJTOT, TNOS,
     &                    XMINZ, FMINZ, XLL, XRL, IXAT, MXSP, NLOC )
C  ---------------------------------------------------------------------
C  Description:
C    Writes all necessary information to a file (dump file), so that
C    the search can be restart at a later time.
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	INCLUDE 'ptml.h'
	DIMENSION XMINZ(MXV,MXNM), FMINZ(MXNM)
	DIMENSION XLL(NOD), XRL(NOD), IXAT(NOD)
C
	DIMENSION IRLX(25)
C
	OPEN (10,FILE='DUMP',FORM='UNFORMATTED')
	WRITE (10) NOD
	WRITE (10) NOF
	WRITE (10) ITTHR
	WRITE (10) HEAL
	WRITE (10) NEARN
	WRITE (10) NSAMPL
	WRITE (10) SIGMA
	WRITE (10) COVTHR
	WRITE (10) XDEPS
	WRITE (10) FDEPS
	WRITE (10) GDEPS
	WRITE (10) NDUMP
	WRITE (10) FINP 
	WRITE (10) FOUT 
	WRITE (10) FMIN
	WRITE (10) POIFIL
	WRITE (10) INDIR
C
	WRITE (10) ITC
	WRITE (10) NOM
	WRITE (10) NFTOT
	WRITE (10) NGTOT
	WRITE (10) NHTOT
	WRITE (10) NJTOT
	WRITE (10) TNOS
	WRITE (10) MXSP
	WRITE (10) NLOC
C
	DO 10,J=1,NOM
		DO 20,I=1,NOD
			WRITE (10) XMINZ(I,J)
20		CONTINUE
10	CONTINUE
	DO 30,I=1,NOM
		WRITE (10) FMINZ(I)
30	CONTINUE
	DO 40,I=1,NOD
		WRITE (10) XLL(I)
40	CONTINUE
	DO 50,I=1,NOD
		WRITE (10) XRL(I)
50	CONTINUE
	DO 60,I=1,NOD
		WRITE (10) IXAT(I)
60	CONTINUE
C
	CALL RLUXUT(IRLX)
	DO 70,I=1,25
		WRITE (10) IRLX(I)
70	CONTINUE
	CLOSE (10)
	WRITE (*,*) 'Current state has been saved.'
	END
C  ---------------------------------------------------------------------
	SUBROUTINE RESTOR ( ITC, NOM, NFTOT, NGTOT, NHTOT, NJTOT, TNOS,
     &                    XMINZ, FMINZ, XLL, XRL, IXAT, MXSP, NLOC )
C  ---------------------------------------------------------------------
C  Description:
C    Reads a dump file and restarts the search.
C  ---------------------------------------------------------------------
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	INCLUDE 'ptml.h'
	DIMENSION XMINZ(MXV,MXNM), FMINZ(MXNM)
	DIMENSION XLL(MXV), XRL(MXV), IXAT(MXV)
C
	DIMENSION IRLX(25)
C
	OPEN (10,FILE='DUMP',FORM='UNFORMATTED')
	READ (10) NOD
	READ (10) NOF
	READ (10) ITTHR
	READ (10) HEAL
	READ (10) NEARN
	READ (10) NSAMPL
	READ (10) SIGMA
	READ (10) COVTHR
	READ (10) XDEPS
	READ (10) FDEPS
	READ (10) GDEPS
	READ (10) NDUMP
	READ (10) FINP 
	READ (10) FOUT 
	READ (10) FMIN
	READ (10) POIFIL
	READ (10) INDIR
C
	READ (10) ITC
	READ (10) NOM
	READ (10) NFTOT
	READ (10) NGTOT
	READ (10) NHTOT
	READ (10) NJTOT
	READ (10) TNOS
	READ (10) MXSP
	READ (10) NLOC
C
	DO 10,J=1,NOM
		DO 20,I=1,NOD
			READ (10) XMINZ(I,J)
20		CONTINUE
10	CONTINUE
	DO 30,I=1,NOM
		READ (10) FMINZ(I)
30	CONTINUE
	DO 40,I=1,NOD
		READ (10) XLL(I)
40	CONTINUE
	DO 50,I=1,NOD
		READ (10) XRL(I)
50	CONTINUE
	DO 60,I=1,NOD
		READ (10) IXAT(I)
60	CONTINUE
C
	DO 70,I=1,25
		READ (10) IRLX(I)
70	CONTINUE
	CALL RLUXIN(IRLX)
	CLOSE (10)
	WRITE (*,*) 'Previous state has been restored.'
	END
C  ---------------------------------------------------------------------
	SUBROUTINE ERRSTP ( SUB, MPISUB, IERR )
C  ---------------------------------------------------------------------
C  Description:
C    Prints an informative message and aborts the program.
C    This routine should be called when an MPI related error
C    is encountered.
C  ---------------------------------------------------------------------
	CHARACTER*(*) SUB, MPISUB
	WRITE (*,*) '========================================='
	WRITE (*,*) 'Routine ', SUB, ' aborted'
	WRITE (*,*) 'After a call to ', MPISUB
	WRITE (*,*) 'The return code was: ', IERR
	WRITE (*,*) 'PTML ABNORMAL TERMINATION'
	WRITE (*,*) '========================================='
	STOP
	END