C***************************************************************
C
C                              FIBO
C                           ==========
C
C     Subrutina de calcul a minimului unei functii neliniare
C     de o variabila reala pe un interval dat:
C
C     min(f(x)), x apartine intervalului inchis [a,b].
C
C     Subrutina implementeaza metoda de cautare directa Fibonacci.
C     Nu necesita evaluarea derivatelor.
C
C     Secventa de apel:
C     call fibo(a,b,eps)
C     unde:
C     a   = marginea inferioara a intervalului de cautare
C     b   = marginea superioara a intervalului de cautare
C     eps = acuratetea de calcul a solutiei, lungimea intervalului
C           in care se gaseste solutia.
C
C     Acuratetea de calcul a solutiei determina numarul de 
C     numere Fibonacci utilizate in procesul de calcul.
C
C     Subrutine apelate: 
C     Subroutine func(x,v)
C       x punctul in care se calculeaza valoarea functiei
C       v = f(x)  valoarea functiei in punctul x.                                                
C
C                                                 Neculai Andrei
C                                          1980, ICI - Bucuresti
C***************************************************************
C               
      subroutine fibo(a,b,eps)
      real*8 a,b,eps

      integer n,j,k,ii,nm2
      real*8 fib(100)
      real*8 fn,vfn,vi,rvi,vb,va,v,w,t,u
      real*8 dif,eps1,dep,w1,pre
      
C--- Output file

      open(unit=1,file='fibo.sol',status='unknown')      
* 
      if(b .lt. a) then
        write(1,108)      
        return
      end if  
        
      dif=b-a

C--- Calculul numerelor Fibonacci. 1,1,2,3,5,8,13,....
      
      fib(1)=1.d0
      fib(2)=1.d0
1       fn=1.d0/eps      
        if(fn - 2.d0) 2,2,3
2       go to 6
3       continue
        n=2
4       n=n+1      
      fib(n)=fib(n-1) + fib(n-2)
        vfn = fib(n)
        if(vfn-fn) 5,7,7
5       go to 4    
6       write(1,100)
      eps = 0.01d0
      go to 1          

C--- Prima iteratie      

7     j=n
      k=n-2      
      vi=b-a
      rvi=fib(k)*vi/fib(j)
      w=a+rvi
      v=b-rvi
      
      call func(w,t)
      call func(v,u)
      
      ii=1

      write(1,101)
      write(1,102) ii,vi,a,b,rvi,v,u

C--- Continuarea iteratiilor.

      k=k-1
      j=j-1
      nm2=n-2
      
      do 14 i=1,nm2
        if(u-t) 8,8,11
8         a=a+rvi
          vi=b-a
          w=v
          t=u
          rvi=fib(k)*vi/fib(j)
          v=b-rvi
          
          call func(v,u)
          
          ii=i+1
          k=k-1
          j=j-1
          if(k-1) 9,10,10
9         k=1
10        continue

          write(1,102) ii,vi,a,b,rvi,w,t
          write(1,103) v,u
          
          go to 14
          
11        b=b-rvi
          vi=b-a
          v=w
          u=t
          rvi=fib(k)*vi/fib(j)
          w=a+rvi
          
          call func(w,t)
          
          ii=i+1
          k=k-1
          j=j-1
          if(k-1) 12,13,13
12        k=1
13        continue
          write(1,102) ii,vi,a,b,rvi,v,u
          write(1,103) w,t                                                
14    continue

C--- Ultima iteratie. Calculul intervalului final.
      
      eps1=0.001*w
      dep=w+eps1
      
      call func(dep,w1)
      
      if(w1-t) 15,15,16
      
15      call func(b,vb)

        write(1,104) w,b
        write(1,105) t,vb
        go to 17

16      call func(a,va)

        write(1,104) w,a
        write(1,105) t,va

C--- Calculul preciziei.

17    pre=(w-a)/dif

      write(1,106) pre
      write(1,107) eps
      
C--- Formate

100   format(4x,'Acuratetea de calcul eps insuficienta.'/,
     *       4x,'Se considera eps=0.01')              

101   format(2x,1hk,5x,2hvi,10x,2hAK,11x,2hBK,9x,3hRVI,
     *       11x,1hX,10x,4hF(x))

102   format(1x,i2,2x,e10.3,2x,e10.3,2x,e10.3,2x,e10.3,
     *             2x,e10.3,2x,e10.3)
     
103   format(53x,e10.3,2x,e10.3)

104   format(10x,'Intervalul final',4x,6h  X1 =,e10.3,
     *        3x,6h  X2 =,e10.3)

105   format(10x,'Valoarea functiei',3x,6hf(X1)=,e10.3,
     *        3x,6hf(X2)=,e10.3)

106   format(10x,'Acuratetea de calcul obtinuta:',7x,e10.3)

107   format(10x,'Acuratetea ceruta:            ',7x,e10.3)

108   format(10x,'Marginea superioara mai mica decat cea inferioara')

      return
      end
      
C**************************************************************


      subroutine func(x,y)
      real*8 x,y

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

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

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

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

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

           
      return
      end
C**************************************************************




C**************************************************************
C
C    Main program
      
      real*8 a,b,eps

* Exemplul 1
c1       a=0.d0
c1       b=10.d0
c1       eps=0.0000003d0
       
* Exemplul 2      
c2       a=0.75d0
c2       b=1.25d0
c2       eps=0.0000003d0

* Exemplul 3
c3       a=4.d0
c3       b=10.d0
c3       eps=0.0000003d0

* Exemplul 4
c4       a=-2.d0
c4       b= 2.d0
c4       eps=0.0000003d0

* Exemplul 5
c5       a=-3.d0
c5       b= 0.d0
c5       eps=0.0000003d0

* Exemplul 6
        a=1.d0
        b=3.d0
        eps=0.0000003d0
             
      call fibo(a,b,eps)
      
      stop
      end
C*****************************************************Last line                                                                       
      
      