program simplex
implicit none

integer(4),parameter :: ndim=2 !stevilo neodvisnih spremenljivk
real(8),parameter    :: eps=1.d-28 ! konvergencni kriterij
real(8)    :: p(ndim+1,ndim),y(ndim+1),funk
integer(4) :: iter
external funk

! na zacetku si izberemo eno tocko vec kot je neodvisnih spremenljivk
p(1,1)=.5d0
p(1,2)=.5d0
p(2,1)=30.d0
p(2,2)=30.d0
p(3,1)=-2.d0
p(3,2)=-2.d0

! priredim imena arrayev kot so v podprogramu
y(1)=funk(p(1,:))
y(2)=funk(p(2,:))
y(3)=funk(p(3,:))

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

! klic podprogramov
call amoeba(p,y,ndim+1,ndim,ndim,eps,funk,iter)

!izpis koncnega rezultata

write(*,*)'iter',iter
write(*,*)p(1,1),p(1,2),y(1)
write(*,*)p(2,1),p(2,2),y(2)
write(*,*)p(3,1),p(3,2),y(3)

end


! funkcija; tukaj sem imel nekaj tezav, preden sem nasel "pravo". 

function funk(x)
implicit none
integer(4),parameter :: ndim=2
real(8) :: x(ndim),funk
!funk=(x(1)-1.d0)**2+(x(2)-2.d0)**2
funk=x(1)**2+x(2)**2+1.d0
!funk=sin(x(1))*sin(x(2))
return
end






      SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter)
      implicit none
      integer(4),PARAMETER :: NMAX=20
      integer(4),PARAMETER :: ITMAX=5000
      INTEGER(4) :: iter,mp,ndim,np !,NMAX,ITMAX
      REAL(8) :: ftol,p(mp,np),y(mp),funk
!C     USES amotry,funk
      INTEGER(4) :: i,ihi,ilo,inhi,j,m,n
      REAL(8) :: rtol,suma,swap,ysave,ytry,psum(NMAX),amotry
      EXTERNAL funk

      iter=0

1     do 12 n=1,ndim
        suma=0.
        do 11 m=1,ndim+1
          suma=suma+p(m,n)
11      continue
        psum(n)=suma
12    continue
2     ilo=1
      if (y(1).gt.y(2)) then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 13 i=1,ndim+1
        if(y(i).le.y(ilo)) ilo=i
        if(y(i).gt.y(ihi)) then
          inhi=ihi
          ihi=i
        else if(y(i).gt.y(inhi)) then
          if(i.ne.ihi) inhi=i
        endif
13    continue
      rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))+1.d-20)

      write(*,*)rtol

      if (rtol.lt.ftol) then
        swap=y(1)
        y(1)=y(ilo)
        y(ilo)=swap
        do 14 n=1,ndim
          swap=p(1,n)
          p(1,n)=p(ilo,n)
          p(ilo,n)=swap
14      continue
        return
      endif
      if (iter.ge.ITMAX) pause 'ITMAX exceeded in amoeba'
      iter=iter+2
      ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.0)
      if (ytry.le.y(ilo)) then
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.0)
      else if (ytry.ge.y(inhi)) then
        ysave=y(ihi)
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.5)
        if (ytry.ge.ysave) then
          do 16 i=1,ndim+1
            if(i.ne.ilo)then
              do 15 j=1,ndim
                psum(j)=0.5*(p(i,j)+p(ilo,j))
                p(i,j)=psum(j)
15            continue
              y(i)=funk(psum)
            endif
16        continue
          iter=iter+ndim
          goto 1
        endif
      else
        iter=iter-1
      endif
      goto 2
      END


      FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac)
      implicit none
      integer(4),PARAMETER ::NMAX=20
      INTEGER(4) :: ihi,mp,ndim,np,j
      REAL(8) :: amotry,fac,p(mp,np),psum(np),y(mp),funk
!      PARAMETER (NMAX=20)
      EXTERNAL funk
!CU    USES funk
!      INTEGER j
      REAL(8) ::  fac1,fac2,ytry,ptry(NMAX)

      fac1=(1.-fac)/ndim
      fac2=fac1-fac
      do 11 j=1,ndim
        ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
11    continue
      ytry=funk(ptry)
      if (ytry.lt.y(ihi)) then
        y(ihi)=ytry
        do 12 j=1,ndim
          psum(j)=psum(j)-p(ihi,j)+ptry(j)
          p(ihi,j)=ptry(j)
12      continue
      endif
      amotry=ytry
      return
      END