      subroutine simplx (f,n,step,ns,ips,maxnfe,cmode,x,fx,
     *                   nfe,s,fs,iflag,lower,upper,count_lim,dfmin_lim)
      use subplex_global
      use global, only: var_vect
      implicit none
!      INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
!      INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
!      INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision
c
      integer n,ns,maxnfe,nfe,iflag,index
      integer ips(ns)
      real(DP) step(n),x(n),fx,s(ns,ns+3),fs(ns+1)
      real(DP) lower(n),upper(n),dfmin,junk
      logical cmode
c
c                                         Coded by Tom Rowan
c                            Department of Computer Sciences
c                              University of Texas at Austin
c
c simplx uses the Nelder-Mead simplex method to minimize the
c function f on a subspace.
c
c input
c
c   f      - function to be minimized, declared external in
c            calling routine
c
c   n      - problem dimension
c
c   step   - stepsizes for corresponding components of x
c
c   ns     - subspace dimension
c
c   ips    - permutation vector
c
c   maxnfe - maximum number of function evaluations
c
c   cmode  - logical switch
c            = .true.  : continuation of previous call
c            = .false. : first call
c
c   x      - starting guess for minimum
c
c   fx     - value of f at x
c
c   nfe    - number of function evaluations
c
c   s      - double precision work array of dimension .ge.
c            ns*(ns+3) used to store simplex
c
c   fs     - double precision work array of dimension .ge.
c            ns+1 used to store function values of simplex
c            vertices
c
c output
c
c   x      - computed minimum
c
c   fx     - value of f at x
c
c   nfe    - incremented number of function evaluations
c
c   iflag  - error flag
c            = -1 : maxnfe exceeded
c            =  0 : simplex reduced by factor of psi
c            =  1 : limit of machine precision
c            =  2 : reached fstop
c            =  3 : difference between minimum values is less than 1%
c
c common
c
      integer nsmin,nsmax,irepl,ifxsw,nfstop,nfxe
      real(DP) alpha,beta,gamma,delta,psi,omega,
     *     bonus,fstop,fxstat,ftest
      logical minf,initx,newx
c
      common /usubc/ alpha,beta,gamma,delta,psi,omega,nsmin,
     *               nsmax,irepl,ifxsw,bonus,fstop,nfstop,
     *               nfxe,fxstat(4),ftest,minf,initx,newx
c
      real(DP) fbonus,sfstop,sfbest
      logical new
c
      common /isubc/ fbonus,sfstop,sfbest,new
c
c local variables
c
      integer i,icent,ih,il,inew,is,itemp,j,npts,count_lim
      real(DP) dist,dum(1),fc,fe,fr,tol,dfmin_lim
      logical small,updatc
c
      save
c
c subroutines and functions
c
      external f,calcc,dist,evalf,newpt,order,start
c   blas
      external dcopy
c   fortran
      intrinsic min
c
c-----------------------------------------------------------
c
      
      if (cmode) go to 50
      npts = ns+1
      icent = ns+2
      itemp = ns+3
      updatc = .false.
      call start (n,x,step,ns,ips,s,small)
      if (small) then
        iflag = 1
        return
      end if
      if (irepl .gt. 0) then
        new = .false.
        call evalf (f,ns,ips,s(1,1),n,x,fs(1),nfe,lower,upper)
      else
        fs(1) = fx
      end if
      new = .true.
      do 10 j = 2,npts
        call evalf (f,ns,ips,s(1,j),n,x,fs(j),nfe,lower,upper)
   10 continue
      il = 1
      call order (npts,fs,il,is,ih)
      tol = psi*dist(ns,s(1,ih),s(1,il))
c
c     main loop
c
   20 continue
        
        call calcc (ns,s,ih,inew,updatc,s(1,icent))
        updatc = .true.
        inew = ih
c
c       reflect
c
        call newpt (ns,alpha,s(1,icent),s(1,ih),.true.,
     *              s(1,itemp),small)
        if (small) go to 40
        call evalf (f,ns,ips,s(1,itemp),n,x,fr,nfe,lower,upper)
        if (fr .lt. fs(il)) then
c
c         expand
c
          call newpt (ns,-gamma,s(1,icent),s(1,itemp),
     *                .true.,s(1,ih),small)
          if (small) go to 40
          call evalf (f,ns,ips,s(1,ih),n,x,fe,nfe,lower,upper)
          if (fe .lt. fr) then
            fs(ih) = fe
          else
            call dcopy (ns,s(1,itemp),1,s(1,ih),1)
            fs(ih) = fr
          end if
        else if (fr .lt. fs(is)) then
c
c         accept reflected point
c
          call dcopy (ns,s(1,itemp),1,s(1,ih),1)
          fs(ih) = fr
        else
c
c         contract
c
          if (fr .gt. fs(ih)) then
            call newpt (ns,-beta,s(1,icent),s(1,ih),.true.,
     *                  s(1,itemp),small)
          else
            call newpt (ns,-beta,s(1,icent),s(1,itemp),
     *                  .false.,dum(1),small)
          end if
          if (small) go to 40
          call evalf (f,ns,ips,s(1,itemp),n,x,fc,nfe,lower,upper)
          if (fc .lt. min(fr,fs(ih))) then
            call dcopy (ns,s(1,itemp),1,s(1,ih),1)
            fs(ih) = fc
          else
c
c           shrink simplex
c
            do 30 j = 1,npts
              if (j .ne. il) then
                call newpt (ns,-delta,s(1,il),s(1,j),
     *                      .false.,dum(1),small)
                if (small) go to 40
                call evalf (f,ns,ips,s(1,j),n,x,fs(j),nfe,lower,upper)
              end if
   30       continue
          end if
          updatc = .false.
        end if
        call order (npts,fs,il,is,ih)
c
c       check termination
c
   40   continue
        fmin_subplex(fmin_index,1) = minval(fs)
!        write(*,*) 1
!        read(*,*) junk
        if(fmin_subplex(fmin_index,1) < fmin(1,1))then
            fmin(1,1) = fmin_subplex(fmin_index,1)
            fmin(2,1) = 0  
!            write(*,*) 2
!            read(*,*) junk
        elseif((fmin_subplex(fmin_index,1) == fmin(1,1))
     *                                    .and.(fmin(1,1)<1e8))then
            fmin(2,1) = fmin(2,1) + 1
!            write(*,*) 3
!            read(*,*) junk
        endif
!        write(*,*) 'fmin count   ', fmin(2,1)
!        write(*,*) fmin_index,'     ',fmin_subplex(fmin_index,1)
        fmin_index = fmin_index + 1
        if(fmin_index > 2)then  !
            dfmin = 1 - fmin_subplex(fmin_index - 1,1)/ 
     *              fmin_subplex(fmin_index - 2,1)
!            write(*,*) 4
!            read(*,*) junk    
            if(dfmin == 1.)then
                count = count + 1
            else
                count = 0
            endif
        else
!            write(*,*) 5
!            read(*,*) junk 
            dfmin = 0
            count = 0
        endif
       
        if (irepl .eq. 0) then
          fx = fs(il)
        else
          fx = sfbest
        end if
   50   continue
          if (nfstop .gt. 0 .and. fx .le. sfstop .and.
     *      nfxe .ge. nfstop) then
          iflag = 2
          write(subplex_log,*) 'reached fstop'
!            write(*,*) 6
!            read(*,*) junk 
        else if (nfe .ge. maxnfe) then
          iflag = -1
          write(subplex_log,*) 'maxnfe exceeded'
!          write(*,*) 7
!            read(*,*) junk 
        else if (dist(ns,s(1,ih),s(1,il)) .le. tol .or.
     *           small) then
          iflag = 0
          write(subplex_log,*) 'simplex reduced by factor of psi'
!          write(*,*) 8
!            read(*,*) junk 
        elseif((dfmin < dfmin_lim).and.(dfmin > 0)) then
          iflag = 3
          write(subplex_log,*) 'small difference between minimum values'
!          write(*,*) 9
!            read(*,*) junk 
        elseif((fmin(2,1) > count_lim)) then
          iflag = 3
          write(subplex_log,*) 'repeat minimum value for too many times'
     *                         ,count 
!            write(*,*) 10
!            read(*,*) junk 
        else
            if(count > 2)then
                do index = 2,3
                    if(var_vect(index) == 1)then
                        if(x(index-1)  == upper(index-1))then
                            upper(index - 1) = upper(index - 1)*1.1
                            write(subplex_log,10035) upper(index - 1)
                        elseif(x(index - 1) == lower(index - 1))then
                            lower(index - 1) = lower(index - 1)*0.9
                            write(subplex_log,10036) lower(index - 1)
                        endif
                    elseif(var_vect(index) == 2)then
                        if(x(index-1)  == upper(index-1))then
                            upper(index - 1) = upper(index - 1)*1.1
                            write(subplex_log,10037) upper(index - 1)
                        elseif(x(index - 1) == lower(index - 1))then
                            lower(index - 1) = lower(index - 1)*0.9
                            write(subplex_log,10038) lower(index - 1)
                        endif   
                    endif
                enddo
            endif
          go to 20
        end if
c
c     end main loop, return best point
c
      do 60 i = 1,ns
        x(ips(i)) = s(i,il)
   60 continue
   
10035  FORMAT('Increased upper bound on GHX to ',f12.5)  
10036  FORMAT('Decreased lower bound on GHX to ',f12.5) 
10037  FORMAT('Increased upper bound on CT to ',f12.5)  
10038  FORMAT('Decreased lower bound on CT to ',f12.5) 

     
      return
      end
