module opt_global
      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   
      INTEGER :: JONES,logfile   
      integer :: exit_early,first_val
      real(DP) :: exit_early_factor,exit_early_mat(5000),emin
end module 

!+-----------------------------------------------------------------------+
!| Program       : main.f                                                |
!| Last modified : 07-16-2001                                            |
!| Written by    : Owen Esslinger, Joerg Gablonsky, Alton Patrick        |
!| The main program to run DIRECT on test problems.                      |
!+-----------------------------------------------------------------------+
      subroutine main_direct(L_hyb_m,C_CCCT,DIRx,DIRf,DIRmaxf,diff,m_or_s,s,l,u,x_in)
	  use opt_global
	  use global, only: double,strat,var_vect,downstream,Qset,conv,guess,select_file,var_vect_pl,&
	                    first_cost,total_fuel_cost
	  use HyGSHP_calc
	  use HyGSHP_calc_down
	  use HyGSHP_calc_conv
	  use global_print
      IMPLICIT NONE      
      INTEGER Maxdim
      PARAMETER (Maxdim = 128)      
      integer i, problem, m_or_s,s
      
!+----------------------------------------------------------------------+
!| DIRECT specific variables.                                           |
!+----------------------------------------------------------------------+
      real(DP) DIReps, DIRf
      Integer DIRmaxf, DIRmaxT
      Integer DIRalg
      Integer IError
      real(DP) :: fglobal, fglper, volper, sigmaper,volper_main,volper_starter, &
                  sigmaper_starter,sigmaper_main,x_in(300)
      real(DP) :: u(Maxdim), l(Maxdim),L_hyb_m,C_CCCT,multiplier_l,multiplier_u
      integer n
      real(DP) :: DIRx(Maxdim)

!+-----------------------------------------------------------------------+
!| Variables to pass user defined data to the function to be optimized.  |
!+-----------------------------------------------------------------------+
      INTEGER iisize, idsize, icsize
      Parameter (iisize = 300)
      Parameter (idsize = 300)
      Parameter (icsize = 30)
      INTEGER iidata(iisize)
      real(DP) ddata(idsize),PCmin
      Character*40 cdata(icsize)
      
      Integer resultfile,nom,shrink,ii,nom_main,nom_starter,shrink_main, &
              shrink_starter,cont_starter,cont_main,cont
      real(DP) per2
      real(DP) dtime,diff,tarray(2),term_fmin,restart_lim       
      
      namelist /direct_optim/ DIReps,DIRmaxf,DIRmaxT,DIRalg, &
                              logfile,fglobal,fglper, &
                              term_fmin,restart_lim,nom_main,multiplier_l, &
                              multiplier_u,volper_main,sigmaper_main, &
                              volper_starter,sigmaper_starter,nom_starter, &
                              shrink_main,shrink_starter,cont_starter,cont_main
!+----------------------------------------------------------------------+
!| Read in the problem specific data and the parameters for DIRECT.     |
!+----------------------------------------------------------------------+
      open(unit = 11, file = '..\..\namelist_files\direct_optim.txt')
      read(11,nml = direct_optim)
      close(11)
      
      n = var_vect(1)
      
      if(double == 0)then
        nom = nom_starter
      elseif(double == 1)then
        nom = nom_main
      endif
      
      if(s == 1)then
        cont = cont_starter
      else
        cont = cont_main
      endif
      
      if(m_or_s == 0)then
        volper = volper_starter
        sigmaper = sigmaper_starter
        shrink = shrink_starter
      else
        volper = volper_main
        sigmaper = sigmaper_main
        shrink = shrink_main
      endif
      
      if(guess == 1)then
        do ii = 2,n+1
            if((var_vect(ii) == 1).or.(var_vect(ii) == 2))then
                l((ii-1)) = x_in(var_vect(ii))-x_in(var_vect(ii))*0.5 
                u((ii-1)) = x_in(var_vect(ii))+x_in(var_vect(ii))*0.5 
                
!                l((ii-1)) = x_in(var_vect(ii))*0.5 
!                u((ii-1)) = x_in(var_vect(ii))*2 
            elseif(var_vect(ii) == 8)then
                l((ii-1)) = max(x_in(var_vect(ii))-x_in(var_vect(ii))*0.5,0.)
                u((ii-1)) = min(x_in(var_vect(ii))+x_in(var_vect(ii))*0.5,1.)
                
!                l((ii-1)) = max(x_in(var_vect(ii))*0.5,0.)
!                u((ii-1)) = min(x_in(var_vect(ii))*2,1.)
            else
                l((ii-1)) = x_in(var_vect(ii))-10
                u((ii-1)) = x_in(var_vect(ii))+10
            endif       
        enddo
        goto 507
      endif      
      
      do ii = 2,n
        if((var_vect(ii) == 1).or.(var_vect(ii) == 2))then
            goto 505
        elseif(var_vect(ii) > 2)then
            shrink = 5000
            goto 505
        endif
      enddo
      
 505  if(strat > 8) shrink = 5000
      
      !If not using a 1yr optimization to get the starting point
      if(double == 0)then
          if((strat <= 8).and.(var_vect(2) == 1))then
              l(1) = max(multiplier_l*L_hyb_m,0.)
              u(1) = multiplier_u*L_hyb_m
          endif
          do ii = 2,n
              if((strat <= 8).and.(var_vect(ii) == 2))then    
                  if(strat /= 0) l(ii-1) = max(multiplier_l*C_CCCT,0.)              
                  if(strat /= 0) u(ii-1) = multiplier_u*C_CCCT
                  goto 500
              endif
          enddo
 500      nom = nom_starter
      elseif(double == 1)then
          !Will need to tweak this for different CS's in the future
          !The idea is to add to the temperature to get the bounds for it
          !May make this an advanced input
          
          do ii = 2,n
            if(var_vect(ii) > 2)then
                l(ii-1) = max(DIRx(var_vect(ii)) - 3,l(ii-1))
                u(ii-1) = min(DIRx(var_vect(ii)) + 3,u(ii-1))
            endif
          enddo
                   
          if((strat <= 8).and.(var_vect(2) == 1))then
              l(1) = 0.75*L_hyb_m
              u(1) = 1.5*L_hyb_m
          endif
          
          do ii = 2,n
              if((strat <= 8).and.(var_vect(ii) == 2))then
                  if(strat /= 0) l(ii-1) = 0.5*C_CCCT              
                  if(strat /= 0) u(ii-1) = 1.5*C_CCCT
                  goto 501
              endif
          enddo
 501      nom = nom_main
          do ii = 1,n
            if(l(ii) == u(ii))then
                l(ii) = l(ii)-0.01*l(ii)
                u(ii) = u(ii)+0.01*u(ii)
            endif
          enddo
      endif
!+----------------------------------------------------------------------+
!|  Define and open the logfile and the resultfile, where we store the  |
!|  results of the run. We save the problem number, the number of       |
!|  function evaluations needed, the time used by DIRECT, the percent   |
!|  error and a flag signaling if we used DIRECT or DIRECT-l in the     |
!|  resultfile.                                                         |
!+----------------------------------------------------------------------+
507      resultfile = 29
      if(select_file == 1)then
        open(logfile, file=file15)
      elseif(select_file == 2)then
        open(logfile, file=file15a) 
      endif
      open(resultfile, file=file16)
      CALL mainheader

!+----------------------------------------------------------------------+
!| Initialize and start the timing.                                     |
!+----------------------------------------------------------------------+
        tarray(1) = 0.D0
        tarray(2) = 0.D0
        diff=dtime(tarray)
        
!+----------------------------------------------------------------------+
!| If the budget is lower zero, multiply the absolute value of the      |
!| budget by the dimension of the problem.                              |
!+----------------------------------------------------------------------+
        if (DIRmaxf .lt. 0) then
           DIRmaxf = -DIRmaxf*n
        endif
        Dirf = 0
!+----------------------------------------------------------------------+
!| Call the optimization method.                                        |
!+----------------------------------------------------------------------+
        if(conv == 1)then
            CALL DIRect(HyGSHP_sub_conv, DIRx, n, DIReps, DIRmaxf, DIRmaxT, &
                        DIRf, l, u, DIRalg, Ierror,  &
                        fglobal, fglper, volper, sigmaper, &
                        iidata, iisize, ddata, idsize, cdata, icsize,term_fmin, &
                        restart_lim,cont,nom, shrink,PCmin )
        elseif(downstream == 0)then
            CALL DIRect(HyGSHP_sub, DIRx, n, DIReps, DIRmaxf, DIRmaxT, &
                        DIRf, l, u, DIRalg, Ierror,  &
                        fglobal, fglper, volper, sigmaper, &
                        iidata, iisize, ddata, idsize, cdata, icsize,term_fmin, &
                        restart_lim,cont,nom, shrink,PCmin )
        else
            CALL DIRect(HyGSHP_sub_down, DIRx, n, DIReps, DIRmaxf, DIRmaxT, &
                        DIRf, l, u, DIRalg, Ierror,  &
                        fglobal, fglper, volper, sigmaper, &
                        iidata, iisize, ddata, idsize, cdata, icsize,term_fmin, &
                        restart_lim,cont,nom, shrink,PCmin )
        endif
!+----------------------------------------------------------------------+
!| Give out the results of the optimization.                            |
!+----------------------------------------------------------------------+
        if(PCmin < 1.)then
            write(logfile,*) 'There is no PC'
        else
            write(logfile,*) 'There is PC'    
        endif

        DIRx(n+1) = first_cost
        DIRx(n+2) = total_fuel_cost

        Write(*,100)
        Write(*,110) IError
        Write(*,120) (DIRx(i),i=1,n)
        Write(*,130) DIRf
        Write(*,140) DIRmaxf
        Write(logfile,100) 
        Write(logfile,110) IError
        Write(logfile,120) (DIRx(i),i=1,n)
        Write(logfile,130) DIRf
        Write(logfile,140) DIRmaxf
        diff=dtime(tarray)
        write(*,150) diff
        write(logfile,150) diff
!+----------------------------------------------------------------------+
!| Calculate the percent error.                                         |
!+----------------------------------------------------------------------+
        per2 = DIRf - fglobal
        if (fglobal .eq. 0.D0) then
           per2 = per2*100.D0
        else
           per2 = per2/ abs(fglobal)*100.D0
        end if
!+----------------------------------------------------------------------+
!| Save the results in the extra resultfile for use with Matlab etc.    |
!+----------------------------------------------------------------------+
      write(resultfile,200) problem, DIRmaxf, diff, per2, DIRalg
!+----------------------------------------------------------------------+
!| Close the logfile and rsultfile.                                     |
!+----------------------------------------------------------------------+
      close(logfile)
      close(resultfile)
100   FORMAT('-------------- Final result ------------------')
110   FORMAT('DIRECT termination flag : ',I3)
120   FORMAT('DIRECT minimal point    : ',20f17.2)
130   FORMAT('DIRECT minimal value    : ',f12.2)
140   FORMAT('DIRECT number of f-eval : ',I5)
150   FORMAT('Time needed             : ',f20.6,' seconds.')
200   FORMAT(I3, ' ', I6, ' ', f12.4, ' ', e10.5, ' ',I6)
      end subroutine main_direct

!+----------------------------------------------------------------------+
!| Give out a header for the main program.                              |
!+----------------------------------------------------------------------+
      SUBROUTINE mainheader
      use opt_global
      IMPLICIT None

      write(*,100)
      write(*,110)
      write(*,120)
      write(*,130)
      write(*,140)
      write(*,150)
      write(*,160)
      write(*,170)
      write(*,200)
      write(logfile,100)
      write(logfile,110)
      write(logfile,120)
      write(logfile,130)
      write(logfile,140)
      write(logfile,150)
      write(logfile,160)
      write(logfile,170)
      write(logfile,200)
100   FORMAT('+----------------------------------------+')
110   FORMAT('|    The optimization uses a modified    |')
120   FORMAT('|     version of the DIRECT algorithm    |')
130   FORMAT('|              developed by              |')
140   FORMAT('|    Owen Esslinger, Joerg Gablonsky,    |')
150   FORMAT('|              Alton Patrick             |')
160   FORMAT('|               04/15/2001               |')
170   FORMAT('|          DIRECT Version 2.0.4          |')
200   FORMAT('+----------------------------------------+')

      end subroutine mainheader
