1.
program forktest
! use the forking method to find the minimum
! of a function
  interface
    real function f(x)
      real, intent(in) :: x
    end function
  end interface
  real x0
  x0 = fork(f, 5.0, 10.0)
  write(*,*) x0, f(x0)
contains
real function fork(f, x0, x1)
  interface
    real function f(x)
      real, intent(in) :: x
    end function
  end interface
  real, intent(in) :: x0, x1
  real, parameter :: limit=0.001
  real :: q, a, b, c, fa, fb, fc, x, fx
  q=(3-sqrt(5.0))/2
  a=x0; fa=f(a)
  c=x1; fc=f(c)
  b=a+q*(c-a); fb=f(b)

  if (fa < fb .or. fc < fb) then
     write(*,'("erroneous initial values",3F10.4)') fa,fb,fc
     fork=0.0
     return
  end if
  do while (abs(a-c)>limit)
  if (c-b > b-a) then
    x=b+q*(c-b); fx=f(x)
    if (fx < fb) then
       a=b; fa=fb
       b=x; fb=fx
    else
       c=x; fc=fx
    end if
  else
    x=b-q*(b-a); fx=f(x)
    if (fx < fb) then
       c=b; fc=fb
       b=x; fb=fx
    else
       a=x; fa=fb
    end if
  end if   
  write(*,*) a,b,c
  end do
  fork=x
end function
end program

real function f(x)
  real, intent(in) :: x
  f = sqrt(x)+sin(x-sqrt(x))
end function

 5.00000000 6.90983009 8.09016991
 6.18033981 6.90983009 8.09016991
 6.90983009 7.36067963 8.09016991
  ...
 7.15514851 7.15601397 7.15654850
 7.15514851 7.15568352 7.15601397  18. iteration
 7.15568352 1.70173728

2.



program fork2
! find the minimum of a function by forking
! use also the derivative  of the function
  interface
    real function f(x)
      real, intent(in) :: x
    end function
    real function df(x)
      real, intent(in) :: x
    end function
  end interface
  real x0
  x0 = fork(f, df, 5.0, 10.0)
  write(*,*) x0, f(x0)
contains
real function fork(f, df, x0, x1)
  interface
    real function f(x)
      real, intent(in) :: x
    end function
    real function df(x)
      real, intent(in) :: x
    end function
  end interface
  real, intent(in) :: x0, x1
  real, parameter :: limit=0.001
  real :: q, a, b, c, fa, fb, fc, dfa, dfb, dfc, x, fx, dfx
  q=(3-sqrt(5.0))/2
  a=x0; fa=f(a) ; dfa=df(a) 
  c=x1; fc=f(c) ; dfc=df(c)
  b=a+q*(c-a); fb=f(b) ; dfb=df(b)
  if (fa < fb .or. fc < fb) then
     write(*,'("erroneous initial values",3F10.4)') fa,fb,fc
     fork=0.0
     return
  end if
  do while (abs(a-c)>limit)
  if (dfb < 0) then
    x=b+q*(c-b); fx=f(x); dfx=df(x)
    a=b; fa=fb; dfa=dfb
    b=x; fb=fx; dfb=dfx
  else
    x=b-q*(b-a); fx=f(x); dfx=df(x)
    c=b; fc=fb; dfc=dfb
    b=x; fb=fx; dfb=dfx
  end if   
  write(*,*) a,b,c
  end do
  fork=x
end function
end program

real function f(x)
  real, intent(in) :: x
  f = sqrt(x)+sin(x-sqrt(x))
end function
real function df(x)
  real, intent(in) :: x
  real :: t
  t = 1.0/(2*sqrt(x))
  df = t+cos(x-sqrt(x))*(1-t)
end function

 6.90983009 8.09016991 10.0000000
 6.90983009 7.63932037 8.09016991
 6.90983009 7.36068010 7.63932037
  ...
 7.15514898 7.15601444 7.15741444
 7.15514898 7.15568399 7.15601444   13. iteration
 7.15568399 1.70173717


3.
program amotest
  use amoeba
  implicit none
  interface
    real function funk(p, ndim)
      integer, intent(in) :: ndim
      real, dimension(ndim) :: p
    end function
  end interface
  integer, parameter :: ndim = 2
  real p(ndim+1, ndim), y(ndim+1)

  ! vertices of the initial simplex
  p(1,:) = (/ 0.0, 0.0 /)
  p(2,:) = (/ 0.0, 2.0 /)
  p(3,:) = (/ 2.0, 0.0 /)

  call simplex(funk, 2, p, y)

  ! write the vertices and function values, which should now
  ! be almost equal since the simplex has shrunk to a small
  ! volume around the optimum
  write(*,*) p(1,:), y(1)
  write(*,*) p(2,:), y(2)
  write(*,*) p(3,:), y(3)

end program

real function funk(p, ndim)
! Rosenbrock's banana function
! amoeba calculates the relative error of the minimum =>
! the function must be positive => add a constant
  implicit none
  integer, intent(in) :: ndim
  real, dimension(ndim) :: p
  funk=100*(p(2)-p(1)**2)**2+(1-p(1))**2+1.0
end function

rtol=   0.00005
nfunk=   86
  0.99709773      0.99482119       1.0000465    
  0.99611020      0.99200487       1.0000205    
   1.0055230       1.0111275       1.0000308    

4.
program amolin
  use amoeba
  implicit none
  interface
    real function funk(p, ndim)
      integer, intent(in) :: ndim
      real, dimension(ndim) :: p
    end function
  end interface

  integer, parameter :: ndim = 2     ! dimension of the space
  real p(ndim+1, ndim), &
       y(ndim+1)

  ! vertices of the initial simplex
  p(1,:) = (/ 0.0, 0.0 /)
  p(2,:) = (/ 0.0, 3.0 /)
  p(3,:) = (/ 3.0, 0.0 /)

  call simplex(funk, 2, p, y)

  write(*,*) p(1,:), y(1)
  write(*,*) p(2,:), y(2)
  write(*,*) p(3,:), y(3)

end program

real function funk(p, ndim)
! L1 residual to be optimised
  implicit none
  integer, intent(in) :: ndim
  real, dimension(ndim) :: p
  integer, parameter :: ndata=6
  real, dimension(ndata) :: &
        xx = (/ 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 /), &
        yy = (/ 0.0, 1.0, 0.0, 2.0, 4.0, 5.0 /)
  integer i
  real s

  s=0.0
  do i=1,ndata
     s=s+abs(yy(i)-p(2)-p(1)*xx(i))
  end do
  funk=s
end function


rtol=   0.00009
nfunk=   84
 0.999904156 -5.40369656E-05 3.00058722
 0.999935985 4.93799598E-05 3.00031996
 1.00009656 -4.58408496E-04 3.00048280