C**********************************************************
C
C     Powell's Method.
C
C**********************************************************
C     Neculai Andrei, 
C     Research Institute for Informatics
C     Advanced Modeling and Optimization Laboratory     
C**********************************************************
C
	SUBROUTINE POWEL(N,NN,X,X1,X2,X3,S,XD,ITECH)
C
	DOUBLE PRECISION X(N),X1(N),X2(N),X3(N),S(N),XD(NN)
	DOUBLE PRECISION BL,EP,EPS,F1,F2,F3,FS,T,V,DELTA,DEL
	DOUBLE PRECISION FA,FB
C
	CHARACTER*1 FD
C
	COMMON /PREC/EPS
	COMMON /TYPE/NO
	COMMON /WRIT/NX
	COMMON /L3A/ICOD
	COMMON /NOF/NBL
C
	WRITE(NX,150)
150	FORMAT(10X,'Start Optimization Process by Powell Method')
C
	GO TO(151,152,153,154) ITECH
151	WRITE(NX,161)
161	FORMAT(10X,'Using Golden Searching Technique')
	GO TO 170
152	WRITE(NX,162)
162	FORMAT(10X,'Using Fibonacci Searching Technique')
	GO TO 170
153	WRITE(NX,163)
163	FORMAT(10X,'Using Quadratic Interpolation of Powell')
	GO TO 170
154	WRITE(NX,164)
164	FORMAT(10X,'Using Dichotomous Searching Technique')
C
170	CONTINUE
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
	NFUNC=0
	NBL=0
	JTYPE=0
	EP=1.0D-3
C
	ITER=0
	CALL FUNC(N,X,F1)
	NFUNC=NFUNC+1
C
	WRITE(NO,50)
50	FORMAT(10X,'***** Optimization Method of Powell *****')
	IF(ITECH.EQ.1) WRITE(NO,161)
	IF(ITECH.EQ.2) WRITE(NO,162)
	IF(ITECH.EQ.3) WRITE(NO,163)
	IF(ITECH.EQ.4) WRITE(NO,164)
C
	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)F1
54	FORMAT(10X,'Function Value in Initial Point = ',F18.6)
C--------------------------------------------------------------
C
	DO 3 I=1,N
	DO 2 J=1,N
	L=J+(I-1)*N
2	XD(L)=0.0
	L=I+(I-1)*N
3	XD(L)=1.0
C
4	CONTINUE
	ITER=ITER+1
C
	IF(ITER.GT.1000) THEN
	WRITE(NX,56)
56	FORMAT(10X,'Too many iterations.',/,
     +         10X,'This function is too complicated for the power',/,
     +         10X,'of the Powell 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) WRITE(NO,300) ITER
300           FORMAT(10X,'ITER = ',I4)
C
      if(jtype.eq.1) then
              write(no,360)
360           format(4x,'The Current Point:')
	        WRITE(NO,302)(I1,X(I1),I1=1,N)
*
              WRITE(NO,301)
301           FORMAT(10X,'Directions:')
              WRITE(NO,302)(I1,XD(I1),I1=1,NN)
302           FORMAT(10X,I4,4X,F18.6)
      endif
C
	I=1
	DO 5 J=1,N
5	S(J)=XD(J)
C--------------------------- 
	DO 25 J=1,N
25	X1(J)=X(J)+EP*S(J)
C
	CALL FUNC(N,X1,FA)
	NFUNC=NFUNC+1
C--------------------------- 
	DO 26 J=1,N
26	X1(J)=X(J)-EP*S(J)
C
	CALL FUNC(N,X1,FB)
	NFUNC=NFUNC+1
C
	IF(FA.GT.FB) THEN
	DO 27 J=1,N
27	S(J)=-S(J)
	ENDIF
C
C---------------------------------------------------
	GO TO(171,172,173,174) ITECH
171	CALL L1(N,X,S,BL)
	GO TO 180
172	CALL L2(N,X,S,BL)
	GO TO 180
173	CALL L3(N,X,S,BL)
	IF(ICOD.EQ.1) RETURN
	GO TO 180
174	CALL L4(N,X,S,BL)
C
180	CONTINUE
C-----------------------------------------------------
C
	NFUNC=NFUNC+NBL
C
	DO 6 J=1,N
6	X2(J)=X(J)+BL*S(J)
C
	CALL FUNC(N,X2,F2)
	NFUNC=NFUNC+1
C
              IF(JTYPE.EQ.1) WRITE(NO,304)I,BL
              IF(JTYPE.EQ.1) WRITE(NO,302)(I1,X2(I1),I1=1,N)
C
	DELTA=(F1-F2)
	M=1
C
7	I=I+1
	IF(I.GT.N) GO TO 11
C
	DO 8 J=1,N
	L=J+(I-1)*N
8	S(J)=XD(L)
C
	DO 30 J=1,N
30	X1(J)=X(J)+EP*S(J)
C
	CALL FUNC(N,X1,FA)
	NFUNC=NFUNC+1
C
	DO 31 J=1,N
31	X1(J)=X(J)-EP*S(J)
C
	CALL FUNC(N,X1,FB)
	NFUNC=NFUNC+1
C
	IF(FA.GT.FB) THEN
	DO 32 J=1,N
32	S(J)=-S(J)
	ENDIF
C
C-----------------------------------------------------
	GO TO(181,182,183,184) ITECH
181	CALL L1(N,X2,S,BL)
	GO TO 190
182	CALL L2(N,X2,S,BL)
	GO TO 190
183	CALL L3(N,X2,S,BL)
	IF(ICOD.EQ.1) RETURN
	GO TO 190
184	CALL L4(N,X2,S,BL)
C
190	CONTINUE
C----------------------------------------------------
C
	NFUNC=NFUNC+NBL
C
	DO 9 J=1,N
9	X3(J)=X2(J)+BL*S(J)
C
	CALL FUNC(N,X3,F3)
	NFUNC=NFUNC+1
C
	        IF(JTYPE.EQ.1) WRITE(NO,304)I,BL
304           FORMAT(10X,'I = ',I2,4X,'BL=',F18.6,'  New Point')
	        IF(JTYPE.EQ.1) WRITE(NO,302)(I1,X3(I1),I1=1,N)
C
	DEL=(F2-F3)
	IF(DEL.GT.DELTA) THEN
	DELTA=DEL
	M=I
	ENDIF
C
              IF(JTYPE.EQ.1) WRITE(NO,305)M,DELTA
305           FORMAT(10X,'M = ',I3,4X,'DELTA =',F18.8)
C
	DO 10 J=1,N
10	X2(J)=X3(J)
	F2=F3
	GO TO 7
C
11	CONTINUE
C             
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,313)
313           FORMAT(10X,'End of a cycle of N Searchings')
              WRITE(NO,314)
314           FORMAT(10X,'New Directions:')
              WRITE(NO,302)(I1,XD(I1),I1=1,NN)
      ENDIF        
C
	IS=0
	DO 12 I=1,N
	IF(DABS(X3(I)-X(I)).LE.EPS) GO TO 12
	IS=1
C
              IF(JTYPE.EQ.1) WRITE(NO,307)X3(I)-X(I)
307           FORMAT(2X,'X3(I)-X(I) =',F18.8)
C
	GO TO 121
12	CONTINUE
C
121	CONTINUE
C
	IF(IS.EQ.0) GO TO 100	
C
	DO 13 I=1,N
13	S(I)=2.*X3(I)-X(I)
C
	CALL FUNC(N,S,FS)
	NFUNC=NFUNC+1
C          
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,320)
320           FORMAT(/,10X,'Point :  2.*X3(I)-X(I)')
              WRITE(NO,302)(I1,S(I1),I1=1,N)
              WRITE(NO,321)FS
321           FORMAT(/,10X,'Function Value in 2*X3(I)-X(I)=',F18.6)
      ENDIF
C  
	IF(FS.GE.F1) THEN
	DO 14 I=1,N
14	X(I)=X3(I)
	F1=F3
	GO TO 4
	ENDIF
C
	T=(F1-2.*F3+FS)*(F1-F3-DELTA)**2
	V=0.5*DELTA*(F1-FS)**2
C     
      if(jtype.eq.1) then
              write(no,3221)f1,f3,delta
3221          format(10x,'f1   =',f18.6,/,
     1               10x,'f2   =',f18.6,/,
     1               10x,'delta=',f18.6)	
              WRITE(NO,322)T,V
322           FORMAT(/,10X,'T=(f1-2f2+f3)*(f1-f2-delta)**2 = ',F18.6,/,
     +                 10X,'V=0.5*(f1-f3)**2               = ',F18.6)
      endif
C
	IF(T.GE.V) THEN
	        if(jtype.eq.1) write(no,355)
355           format(/,4x,'Here T > V: Consider the old directions')
C-----------------------------------------------------------
	    DO 15 I=1,N
15        X(I)=X3(I)
	    F1=F3
	    GO TO 4
	ENDIF
C
	DO 16 I=1,N
16	S(I)=X3(I)-X(I)
C
              if(jtype.eq.1) write(no,350)
350           format(4x,'Here T < V: Compute the direction X3-X')
C--------------------------------------------------------
	GO TO(191,192,193,194) ITECH
191	CALL L1(N,X3,S,BL)
	GO TO 200	
192	CALL L2(N,X3,S,BL)
	GO TO 200
193	CALL L3(N,X3,S,BL)
	IF(ICOD.EQ.1) RETURN
	GO TO 200
194	CALL L4(N,X3,S,BL)
C
200	CONTINUE
C---------------------------------------------
C
	NFUNC=NFUNC+NBL
C             
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,323)
323           FORMAT(10X,'Direction : X3(I)-X(I)')
              WRITE(NO,302)(I1,S(I1),I1=1,N)
              WRITE(NO,324)BL
324           FORMAT(10X,'-------- BL in Direction X3-X = ',F18.6)
      ENDIF
C
	DO 17 I=1,N
17	X(I)=X3(I)+BL*S(I)
C
	CALL FUNC(N,X,F1)
	NFUNC=NFUNC+1
C     
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,325)
325           FORMAT(/,10X,'New point X=X3+BL*(X3-X)')
              WRITE(NO,302)(I1,X(I1),I1=1,N)
              WRITE(NO,326)F1
326           FORMAT(/,10X,'Function Value in X = ',F18.6)
      ENDIF
C
	IS=0
	DO 18 I=1,N
	IF(DABS(BL*S(I)).LE.EPS) GO TO 18
	IS=1 
C
              IF(JTYPE.EQ.1) WRITE(NO,308)BL*S(I)
308           FORMAT(2X,'BL*S(I) =',F18.8)
C
	GO TO 118
18	CONTINUE
C
118	CONTINUE
C
	IF(IS.EQ.0) GO TO 100
C
	IF(M.EQ.N) THEN
	DO 19 J=1,N
	L=J+(N-1)*N
19	XD(L)=S(J)
	GO TO 4
	ENDIF
C
	DO 20 I=M+1,N
	DO 20 J=1,N
	J1=J+(I-2)*N
	J2=J+(I-1)*N
20	XD(J1)=XD(J2)
C
	DO 21 J=1,N
	L=J+(N-1)*N
21	XD(L)=S(J)
	GO TO 4
C
100	CONTINUE
C
	WRITE(NO,79)
79	FORMAT(/5X,'------------------------------------------')
	WRITE(NO,80)ITER,NFUNC
80	FORMAT(10X,'Iterations = ',I4,/,
     +       10X,'No. of Function Evaluations = ',I5)
	WRITE(NO,81) EPS
81	FORMAT(10X,'Accuracy   = ',F12.9)
	WRITE(NO,82)
82	FORMAT(5X,'********** Minimum Point **********')
	WRITE(NO,83)(X(I),I=1,N)
83	FORMAT(10X,F18.6)
	WRITE(NO,84) F1
84	FORMAT(10X,'Function Value in Minimum Point =',F18.6)
C
C
	RETURN
	END
C
C*********************************Last card of POWEL.FOR