c
c
      subroutine mymain
c
c***********************************************************************
c
c     This is an example of an interactive program which performs
c     a direct search on the objective function "objf". You can do
c     a search and replace to substitute your own objective function
c     and recompile.
c     You can also redirect a file to standard input to avoid user
c     interaction:
c         MySearch < sample.dat
c
c***********************************************************************
c
      external objf
c
      parameter (maxn=10)
      integer type, n, i, maxf
      double precision start(maxn), step1(maxn), step2, aux(maxn+2)
      integer iright, set, istd, iaux
      double precision sigma, alpha, beta, gamma
      logical right, stddev
c
      write(*,*) 'Enter the type of search you would like to perform.'
 5    write(*,*) '1 - Coordinate Search'
      write(*,*) '2 - Compass Search'
      write(*,*) '3 - NLess Search'
      write(*,*) '4 - Hooke and Jeeves'
      write(*,*) '5 - Hooke and Jeeves (edited by E. Dolan)'
      write(*,*) '6 - Simplex based on Spendley, Hext, and Himsworth'
      write(*,*) '7 - Simplex based on Nelder and Meade'
      write(*,*) '8 - Sequential Multidirectional Search'
      read(*,*) type
      if (type.lt.1.or.type.gt.8) then
         write(*,*) 'That was not a valid entry. Choose a number',
     c              ' between 1 and 8'
         go to 5
      endif
      write(*,*) 'Enter the dimension of the problem.'
      write(*,*) 'If you are running the default subroutine objf,',
     c     ' then the dimension must be 2.'
      read(*,*) n
      if (n.gt.maxn) then
         write(*,*) 'Increase the parameter maxn before proceeding'
         stop
      endif
      write(*,*) 'Enter the starting values for x.'
      read(*,*) (start(i), i=1,n)
      write(*,*) 'Enter the maximum number of function calls.',
     c           ' (-1 for no max)'
      read(*,*) maxf
c     Input for pattern searches
      if (type.ge.1.and.type.le.5) then 
         write(*,*) 'Enter the initial step length. Enter -1 to get',
     c              ' the default of 0.25'
         read(*,*) step1(1)
         write(*,*) 'Enter the final step length. Enter -1 to get',
     c              ' the default of 10E-8'
         read(*,*) step2
         call pattrn(type, objf, n, start, step1, step2, maxf, istat, 
     c        aux, iaux)
c     Input for simplex searches
      else
         write(*,*) 'Enter the simplex type: 1-right 2-regular'
         read(*,*) iright
         if (iright.eq.1) then 
            right = .true.
         elseif (iright.eq.2) then
            right = .false.
         else
            write(*,*) 'Input was invalid. Right simplex was chosen.'
            right = .true.
         endif
         write(*,*) 'Do you want to set the edgelengths or use the ',
     c              'default of 2.0? 1-Set 2-Use Default'
         read(*,*) set
         if (set.eq.1) then
            if (right) then
               write(*,*) 'Enter the n edgelengths. (Direct search ',
     c                    'will calculate length n+1.)'
               read(*,*) (step1(i), i=1,n)
            else
               write(*,*) 'Enter one edgelength of the regular simplex.'
               read(*,*) step1(1)
            endif
         elseif (set.eq.2) then
            step1(1) = -1.0
         else
            write(*,*) 'Input was invalid. Default edgelengths used'
            step1(1) = -1.0
         endif
         write(*,*) 'Stopping criteria on standard deviation or delta?',
     c              '1-stddev 2-delta'
         read(*,*) istd
         if (istd.eq.1) then 
            stddev = .true.
         elseif (istd.eq.2) then
            stddev = .false.
         else
            write(*,*) 'Input was invalid. Stopping on stddev.'
            stddev = .true.
         endif
         write(*,*) 'Enter the final step length. Enter -1 to get',
     c              ' the default of 10E-8'
         read(*,*) step2
         if (type.eq.6) then
            call shh(objf, n, start, right, step1, step2, maxf, 
     c               stddev, istat, aux, iaux)
         elseif (type.eq.7) then
            write(*,*) 'Enter the shrinking coefficient, sigma.'
            write(*,*) '(-1 uses a default of 0.5.)'
            read(*,*) sigma
            write(*,*) 'Enter the reflection coefficient, alpha.'
            write(*,*) '(-1 uses a default of 1.0.)'
            read(*,*) alpha
            write(*,*) 'Enter the contraction coefficient, beta.'
            write(*,*) '(-1 uses a default of 0.5.)'
            read(*,*) beta
            write(*,*) 'Enter the expansion coefficient, gamma.'
            write(*,*) '(-1 uses a default of 2.0.)'
            read(*,*) gamma
            call nm(objf, n, start, right, sigma, alpha, beta, gamma, 
     c              step1, step2, maxf, stddev, istat, aux, iaux)
         else
            call smd(objf, n, start, right, step1, step2, maxf, 
     c               stddev, istat, aux, iaux)
         endif
      endif
      
      if (istat.eq.1) then
         write (*,10) aux(1)
         write (*,20) aux(2)
         write (*,30) (aux(i),i=3,n+2)
         write (*,40) iaux
 10      format('The final delta is: ', 1PE11.4)      
 20      format('The minimum function evaluation is: ', F10.6)      
 30      format('The corresponding x is: ', F10.6)      
 40      format('The number of evaluations was ', I6)
      endif
      end
c
      subroutine objf(x, f, n)
c
c***********************************************************************
c
c     This subroutine implements the objective function. Bounds are
c     constructed by making the function evaluation very large when
c     the x values are out of range.
c
c***********************************************************************
c
      double precision x(n), f
      double precision temp1, temp2, temp3
c
      if (x(1).LT.6.d0.OR.x(2).LT.0.d0) then
         f = 10000000.0
      else
         temp1 = 1.d0 / (27.d0 * dsqrt(3.d0))
         temp2 = x(1) - 3.d0
         temp3 = x(2) * x(2) * x(2)
         f  = temp1 * (temp2*temp2 - 9.d0) * temp3
      endif
c
      end