C*********************************************************
C
C     Nelder - Mead Optimization method.
C     Simplex Method.
C
C*********************************************************   
C     Neculai Andrei, 
C     Research Institute for Informatics
C     Advanced Modeling and Optimization Laboratory  
C*********************************************************
C
	SUBROUTINE NELMED(N,N1,NN,X,X0,XR,XE,XD,VF)
C
	DOUBLE PRECISION X(N),X0(N),XR(N),XE(N)
	DOUBLE PRECISION XD(NN),VF(N1),VF1(100)
C
	DOUBLE PRECISION ALFA,BETA,GAMA,EPS
	DOUBLE PRECISION FX,F0,FR,FE,VFI,VFM,S,T
C
	CHARACTER*1 FD
C
	COMMON /PREC/EPS
	COMMON /TYPE/NO
	COMMON /WRIT/NX
C
	NFUNC=0
	ITER=0
	ALFA=1.0
	BETA=0.5
	GAMA=2.0
	JTYPE=0
C
	WRITE(NX,150)
150	FORMAT(10X,'Start Optimization Process by Nelder-Mead',
     +              ' Method')
C
110	WRITE(*,111)'   Printing [F=File/D=Display] ? '
111	FORMAT(A,$)
	READ(5,112,ERR=113) FD
112	FORMAT(A1)
	IF(FD.NE.'F'.AND.FD.NE.'D') THEN
	WRITE(NX,114)
114	FORMAT(10X,'Error.   Please try again')
	GO TO 110
	ENDIF
	IF(FD.EQ.'F') NO=6
	IF(FD.EQ.'D') NO=5
	GO TO 115
113	WRITE(NX,114)
	GO TO 110
115	CONTINUE
C
	IF(NO.EQ.6) THEN
	OPEN(6,FILE='PRINT.DAT',ACCESS='SEQUENTIAL',
     +         FORM='FORMATTED')
	ENDIF
C
	CALL FUNC(N,X,FX)
	NFUNC=NFUNC+1
C
	WRITE(NO,50)
50	FORMAT(10X,'***** Optimization Method of Nelder-Mead *****')
	WRITE(NO,51)
51	FORMAT(10X,'Initial Point : ')
	WRITE(NO,52)(X(I),I=1,N)
52	FORMAT(10X,F18.6)
	WRITE(NO,53)EPS
53	FORMAT(10X,'Accuracy EPS = ',g30.15)
	WRITE(NO,54)FX
54	FORMAT(10X,'Val. Func. in Initial Point = ',F18.6)
C
C------------- Build up the initial simplex
C              using the initial point X.
C
	N2=N*N
	DO 1 J=1,N
	L=J+N2
1	XD(L)=X(J)
C
	DO 3 I=1,N
	DO 2 J=1,N
	L=J+(I-1)*N
2	XD(L)=X(J)
	L=I+(I-1)*N
3       XD(L)=XD(L)+0.5
C     
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,201)
201           FORMAT(10X,'Initial simplex')
              WRITE(NO,202)(I1,XD(I1),I1=1,NN)
202           FORMAT(10X,I4,5X,F18.6)
      ENDIF
C
C-----------------Compute the values of the objective 
C                 function in vertices of the simplex
C
	DO 5 I=1,N1
	DO 4 J=1,N
	L=J+(I-1)*N
4	X(J)=XD(L)
	CALL FUNC(N,X,FX)
	NFUNC=NFUNC+1
	VF(I)=FX
5	CONTINUE
C
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,203)
203           FORMAT(10X,'Val. Func in points of initial simplex')
              WRITE(NO,202)(I1,VF(I1),I1=1,N1)
      ENDIF        
C
6	ITER=ITER+1
C
	IF(ITER.GT.1001)THEN
	IF(EPS.EQ.0.0) THEN
	WRITE(NX,41)
41	FORMAT(10X,'Too many iterations.',/,
     +         10X,'The accuracy is too high, EPS=0.0',//,
     +         10X,'Please consider another finite, small,',/,
     +         10X,'positive value for accuracy EPS.')
	RETURN
	ENDIF
C
	WRITE(NX,40)
40	FORMAT(10X,'Too many iterations.',/,
     +         10X,'This function is too complicated for the',/,
     +         10X,'power of the Nelder-Mead Optimization Method.',/,
     +         10X,'May be this function is a ridge or a valley.',//,
     +         10X,'Please consider another initial point, or',/,
     +         10X,'try another optimization method')
	RETURN
	ENDIF
C     
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,204)ITER
204           FORMAT(5X,'ITERATION = ',I4)
              write(no,220)
220           format(3x,'Current simplex')
              WRITE(NO,202)(I1,XD(I1),I1=1,NN)
              WRITE(NO,223)
223           FORMAT(10X,'Func. value in points of the current simplex')
      ENDIF
*
	do i=1,n1
	   do j=1,n
 	      l=j+(i-1)*n
	      x(j)=xd(l)
	   end do
	   call func(n,x,fx)
	   vf1(i)=fx
	end do	                
C
              IF(JTYPE.EQ.1) WRITE(NO,202)(I1,VF1(I1),I1=1,N1)
C
	VFI=VF1(1)
	VFM=VFI
	IMI=1
	IMA=1
	DO 7 I=2,N1
	IF(VF1(I).LT.VFI) THEN
	  VFI=VF1(I)
	  IMI=I
	ENDIF
	IF(VF1(I).GT.VFM) THEN
	  VFM=VF1(I)
	  IMA=I
	ENDIF
7	CONTINUE
C
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,205)IMI,VFI
205           FORMAT(10X,'IMI = ',I4,4X,'VFI = ',F18.6)
              WRITE(NO,206)IMA,VFM
206           FORMAT(10X,'IMA = ',I4,4X,'VFM = ',F18.6)
      ENDIF
C
C
C-----------------Compute the centroid of the simplex.
C
	DO 9 J=1,N
	X0(J)=0.0
	DO 8 I=1,N1
	IF(I.EQ.IMA) GO TO 8
	L=J+(I-1)*N
	X0(J)=X0(J)+XD(L)
8	CONTINUE
9	CONTINUE
C
	DO 10 J=1,N
10	X0(J)=X0(J)/FLOAT(N)
	CALL FUNC(N,X0,F0)
	NFUNC=NFUNC+1
C     
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,208)
208           FORMAT(10X,'Centroid Point')
              WRITE(NO,202)(I1,X0(I1),I1=1,N)
              WRITE(NO,211)F0
211           FORMAT(10X,'Val. Func. in Centr. Point = ',F18.6)
      ENDIF
C
C----------------Compute the reflection point.
C
	DO 11 J=1,N
	L=J+(IMA-1)*N
11	XR(J)=(1.+ALFA)*X0(J)-ALFA*XD(L)
	CALL FUNC(N,XR,FR)
	NFUNC=NFUNC+1
C     
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,209)
209           FORMAT(10X,'Reflection Point')
              WRITE(NO,202)(I1,XR(I1),I1=1,N)
              WRITE(NO,210) FR
210           FORMAT(10X,'Func. Val. in Ref. Point = ',F18.6)
      ENDIF
C
C----------------Decision on the way.
C
	IF(FR.LT.VFI) THEN
C
	DO 12 J=1,N
12	XE(J)=GAMA*XR(J)+(1.-GAMA)*X0(J)
	CALL FUNC(N,XE,FE)
	NFUNC=NFUNC+1
C     
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,212)
212           FORMAT(10X,'FR.LT.VFI -- Expandation point')
              WRITE(NO,202)(I1,XE(I1),I1=1,N)
              WRITE(NO,213)FE
213           FORMAT(10X,'Func. Val. in Expan. Point = ',F18.6)
      ENDIF
C
	IF(FE.LT.VFI) THEN
	DO 13 J=1,N
	L=J+(IMA-1)*N
13	XD(L)=XE(J)
	VF1(IMA)=FE
	GO TO 23
	ELSE
	DO 16 J=1,N
	L=J+(IMA-1)*N
16	XD(L)=XR(J)
	VF1(IMA)=FR
	GO TO 23
	ENDIF
C
	ELSE
C
	IS=1
	DO 14 I=1,N1
	IF(I.EQ.IMA) GO TO 14
	IF(FR.GT.VF1(I)) GO TO 14
	IS=0
	GO TO 15
14	CONTINUE
15	IF(IS.EQ.1) GO TO 18
	DO 17 J=1,N
	L=J+(IMA-1)*N
17	XD(L)=XR(J)
	VF1(IMA)=FR
	GO TO 23
C
	ENDIF
C
C----------------Compute the contraction point.
C
18	CONTINUE
	IF(FR.GT.VFM) THEN
	DO 19 J=1,N
	L=J+(IMA-1)*N
19	X(J)=BETA*XD(L)+(1.-BETA)*X0(J)
	ELSE
	DO 20 J=1,N
20	X(J)=BETA*XR(J)+(1.-BETA)*X0(J)
	ENDIF
C
	CALL FUNC(N,X,FX)
	NFUNC=NFUNC+1
C
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,214)
214           FORMAT(10X,'Contraction Point')
              WRITE(NO,202)(I1,X(I1),I1=1,N)
              WRITE(NO,215)FX
215           FORMAT(10X,'Func. Val. in Contr. Point = ',F18.6)
      ENDIF
C
C----------------Reduce the simplex.
C
	IF(FX.GT.VFM) THEN
	DO 21 I=1,N1
	DO 21 J=1,N
	L=J+(I-1)*N
	K=J+(IMI-1)*N
21	XD(L)=(XD(L)+XD(K))/2.
C
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,216)
216           FORMAT(10X,'New Reduced Simplex ')
              WRITE(NO,202)(I1,XD(I1),I1=1,NN)
      ENDIF
C
	GO TO 23
C
	ELSE
C
	DO 22 J=1,N
	L=J+(IMA-1)*N
22	XD(L)=X(J)
	GO TO 23
C
	ENDIF
C
C--------------Test of continuation (accuracy).
C
23	CONTINUE
	S=0.0
	DO 24 I=1,N1
24	S=S+(VF1(I)-F0)**2
	S=S/FLOAT(N+1)
	T=SQRT(S)
C
              IF(JTYPE.EQ.1) WRITE(NO,218)T
218           FORMAT(10X,'*** Accuracy test T = ',F18.6)
C
	IF(T.LE.EPS) GO TO 30
	GO TO 6
C
30	CONTINUE
	DO 31 J=1,N
	L=J+(IMI-1)*N
31	X(J)=XD(L)
C
	WRITE(NO,32)
32	FORMAT(/5X,'------------------------------------------')
	WRITE(NO,33)ITER,NFUNC
33	FORMAT(10X,'Iterations = ',I4,/,
     +         10X,'No. Function Evaluations = ',I5)
	WRITE(NO,34)EPS
34	FORMAT(10X,'Accuracy =   ',F12.9)
	WRITE(NO,35)
35	FORMAT(5X,'********** Minimum Point **********')
	WRITE(NO,36)(X(J),J=1,N)
36	FORMAT(10X,F18.6)
	WRITE(NO,37) VFI
37	FORMAT(10X,'Val. Func. in Minimum Point =',F18.6)
C
	RETURN
	END
C
C*****************************Last card of NELMED.FOR