C***************************************************************
C
C                              MAXFUN
C                           ============
C
C  Subrutina de calcul a maximului unei functii neliniare 
C  de o variabila.
C  Variabila nu este supusa nici unei restrictii.
C  Nu se utilizeaza derivata functiei pentru care se cauta
C  punctul de maxim.
C
C                             max(f(x))
C
C  Subrutina implementeaza metoda de cautare directa Powell.
C  
C     Secventa de apel:
C     call maxfun(eps,pas,x1)
C     unde:
C     eps = acuratetea de calcul a solutiei.
C     pas = lungimea pasului de cautare.
C     x1  = punctul initial (de unde incepe cautarea maximului).
C
C     Subrutine apelate: 
C     Subroutine func(x,y)
C       x punctul in care se calculeaza valoarea functiei
C       y = f(x)  valoarea functiei in punctul x.                                                
C
C                                                 Neculai Andrei
C                                          1980, ICI - Bucuresti
C***************************************************************


		subroutine maxfun(eps,pas,x1)
**        
		double precision t1(3), t2(3), vf(3)
        double precision eps,pas
        double precision x1,y1, x2,y2, x3,y3, x4,y4
        double precision a,b,c,d,e,f,p,r,xe,xmax,vfxe,vfmax
        logical iprint

* Open the output file
		open(unit=2,file='maxfun.out',status='unknown')

		iprint=.true.
       
        call func(x1,y1)        
        if(iprint) then
          write(2,100)
100		  format(/,24x,1hX,20x,4hF(X),17x,5hPASUL)          
          write(2,101) x1,y1,pas
101	      format(15x,e20.13,2x,e20.13,2x,e20.13)
		end if

* Incrementarea variabilei
		x2=x1+pas
        ii=0
1		call func(x2,y2)      
		if(iprint) write(2,101) x2,y2,pas  

* Compararea valorilor functiei n cele doua puncte
		if(y2-y1) 3,6,6
3		if(ii-1) 4,4,5
4		ii=ii+1

* Schimbarea directiei de cautare
		pas=-pas
        x2=x1+pas
        go to 1
5		go to 11
6		continue
		x3=x2+pas
        go to 8
7		x3=x4                		
8		call func(x3,y3)
		if(iprint) write(2,101) x3,y3,pas

* Dublarea pasului        
		pas=pas*2.d0
        x4=x3+pas
        call func(x4,y4)
        if(y4-y3) 9,10,10
9		go to 12
10		go to 7

* Evaluarea a trei puncte in jurul maximului
11		t1(1)=x1
        t1(2)=x1+pas/2.d0
        t1(3)=x2
        go to 13
12		t1(1)=x3
        t1(2)=x3+pas/2.d0
        t1(3)=x4
        go to 13

* Evaluarea functiei in aceste puncte
13		do 14 i=1,3
		  call func(t1(i),vf(i))
14		continue

* Interpolarea patratica utilizand aceste puncte
		if(iprint) write(2,102)
102		format(/,24x,2hX*,19x,5hF(X*))

15		continue
		a=t1(2)-t1(3)
        b=t1(3)-t1(1)
        c=t1(1)-t1(2)
        
        d=t1(2)**2-t1(3)**2
        e=t1(3)**2-t1(1)**2
        f=t1(1)**2-t1(2)**2
        
        p=0.5d0*(d*vf(1)+e*vf(2)+f*vf(3))
        r=       a*vf(1)+b*vf(2)+c*vf(3)
        
        xe=p/r
        call func(xe,vfxe)
        if(iprint) write(2,103)xe,vfxe
103		format(15x,e20.13,2x,e20.13)

		do 18 j=1,3
          t2(j)=dabs(t1(j)-xe)
          if(t2(j)-eps) 16,16,17
16			xmax=t1(j)
			go to 25
17		  continue
18		continue

* Determinarea si inlocuirea punctului de valoare maxima
		if(vf(2)-vf(1)) 20,19,19
19		if(vf(3)-vf(1)) 20,21,21
20		if(vf(3)-vf(2)) 23,22,22
21		jk=1
		go to 24
22		jk=2
		go to 24
23		jk=3
24		t1(jk)=xe
		vf(jk)=vfxe
        go to 15
        
25		call func(xmax,vfmax)
		write(2,104) xmax
104	  format(/15x,'Punctul de maxim determinat. xmax=',2x,e20.13)
		write(2,105) vfmax
105	  format(/15x,'Valoarea functiei in punctul xmax=',2x,e20.13)
        write(2,106) eps
106	  format(/15x,'Acuratetea de calcul a solutiei   ',2x,e20.13)

		return
        end
****************************************************************




		subroutine func(x,y)
        double precision x,y

* Exemplul 1
c1       y = -x*x + 6.d0*x - 2.d0
c1       y =-y
       
* Exemplul 2      
c2       y = x**5 - 2.d0*x**3 +10.d0*dsin(5.d0*x)
c2       y =-y

* Exemplul 3
c3       y = 1.d0 - 10.d0*x + 0.01d0*exp(x)
c3       y =-y

* Exemplul 4  
c4       y=-1.d0+2.d0/(x**2+1.d0)

* Exemplul 5
c       y=-x/(x*x+1.d0) 

* Exemplul 6
c6       y=5.d0*x**5-4.d0*x**4+400.d0*x*dsin(4.d0*x-4.d0)
c6       y=-y

* Exemplul 7
c7        y=x**4-12.d0*x**3+15.d0*x**2+56.d0*x-60.d0

* Exemplul 8
 	     y=-x**4+12.d0*x**3-47.d0*x**2+60.d0*x
         
        return
        end
****************************************************************




* Main Program for MAXFUN subroutine

      double precision eps,pas,x1
      eps= 0.00001d0
      pas= 0.001d0
      x1 = -1.5d0
      call maxfun(eps,pas,x1)
      stop
      end
****************************************************** Last line      


        