direct.f

优化问题中的直接搜索法

      INTEGER anchor(-1:maxdeep),S(maxdiv,2)
      INTEGER point(maxfunc), free
      DOUBLE PRECISION  c(maxfunc,MaxDim)
      DOUBLE PRECISION  thirds(0:maxdeep),levels(0:maxdeep)
      INTEGER length(maxfunc,MaxDim),t,j,actdeep
      INTEGER Minpos,file,maxpos,help,numfunc,file2
      INTEGER ArrayI(MaxDim),maxi,oops,cheat,writed
      INTEGER List2(MaxDim,2),i,actmaxdeep,oldpos
      INTEGER tstart,start,Newtosample
      DOUBLE PRECISION  w(MaxDim),kmax, delta
      INTEGER pos1
C+-----------------------------------------------------------------------+
C| JG 09/25/00 Version counter.                                          |
C+-----------------------------------------------------------------------+
      INTEGER version
      INTEGER oldmaxf,increase, freeold
C+-----------------------------------------------------------------------+
C| JG 09/24/00 Add another actdeep to keep track of the current depths   |
C|             for divisions.                                            |
C+-----------------------------------------------------------------------+
      INTEGER actdeep_div
      DOUBLE PRECISION  oldl(MaxDim), oldu(MaxDim)
C+-----------------------------------------------------------------------+
C|JG 01/13/01 Added epsfix for epsilon update. If eps < 0, we use Jones  |
C|            update formula. epsfix stores the absolute value of epsilon|
C|            then. Also added flag iepschange to store if epsilon stays |
C|            fixed or is changed.                                       |
C+-----------------------------------------------------------------------+
      DOUBLE PRECISION epsfix
      INTEGER iepschange

      INTEGER DIRGetMaxdeep, DIRgetlevel
C+-----------------------------------------------------------------------+
C| JG 01/22/01 fmax is used to keep track of the maximum value found.    |
C+-----------------------------------------------------------------------+
      DOUBLE PRECISION fmax
C+-----------------------------------------------------------------------+
C| JG 01/22/01 Ifeasiblef is used to keep track if a feasible point has  |
C|             been found so far. Ifeasiblef = 0 means a feasible point  |
C|             has been found, Ifeasiblef = 1 if not.                    |
C| JG 03/09/01 IInfeasible is used to keep track if an infeasible point  |
C|             has been found. IInfeasible > 0 means a infeasible point  |
C|             has been found, IInfeasible = 0 if not.                   |
C+-----------------------------------------------------------------------+
      INTEGER Ifeasiblef, IInfesiblef

C+-----------------------------------------------------------------------+
C+-----------------------------------------------------------------------+
C|                            Start of code.                             |
C+-----------------------------------------------------------------------+
C+-----------------------------------------------------------------------+
      writed = 0
      dwrit  = 0
      JONES  = algmethod
C+-----------------------------------------------------------------------+
C| Save the upper and lower bounds.                                      |
C+-----------------------------------------------------------------------+
      DO 150,i=1,n
        oldu(i) = u(i)
        oldl(i) = l(i)
150   CONTINUE

C+-----------------------------------------------------------------------+
C| Set version.                                                          |
C+-----------------------------------------------------------------------+
      version = 204
C+-----------------------------------------------------------------------+
C| Set parameters.                                                       |
C|    If cheat > 0, we do not allow \tilde{K} to be larger than kmax, and|
C|    set \tilde{K} to set value if necessary. Not used anymore.         |
C+-----------------------------------------------------------------------+
      cheat = 0
      kmax  = 1.D10
      mdeep = maxdeep
C+-----------------------------------------------------------------------+
C| Write the header of the logfile.                                      |
C+-----------------------------------------------------------------------+
      CALL DIRheader(logfile, version, x, n, eps, maxf, maxT, l, u,
     +               algmethod, maxfunc, maxdeep, fglobal, fglper,
     +               Ierror, epsfix, iepschange, volper, sigmaper,
     +               iidata, iisize, ddata, idsize, cdata,
     +               icsize)
C+-----------------------------------------------------------------------+
C| If an error has occured while writing the header (we do some checking |
C| of variables there), return to the main program.                      |
C+-----------------------------------------------------------------------+
      IF (Ierror .lt. 0) then
         RETURN
      END IF
C+-----------------------------------------------------------------------+
C| If the known global minimum is equal 0, we cannot divide by it.       |
C| Therefore we set it to 1. If not, we set the divisionfactor to the    |
C| absolute value of the global minimum.                                 |
C+-----------------------------------------------------------------------+
      IF (fglobal .EQ. 0.D0) then
         divfactor = 1.D0
      ELSE
         divfactor = abs(fglobal)
      END IF

C+-----------------------------------------------------------------------+
C| Start of application-specific initialisation.                         |
C+-----------------------------------------------------------------------+
      CALL DIRInitSpecific(x,n,
     +   iidata, iisize, ddata, idsize, cdata, icsize)
C+-----------------------------------------------------------------------+
C| End of application-specific initialisation.                           |
C+-----------------------------------------------------------------------+

C+-----------------------------------------------------------------------+
C| Save the budget given by the user. The variable maxf will be changed  |
C| if in the beginning no feasible points are found.                     |
C+-----------------------------------------------------------------------+
      oldmaxf = maxf
      increase = 0
C+-----------------------------------------------------------------------+
C| Initialiase the lists.                                                |
C+-----------------------------------------------------------------------+
      CALL DIRInitList(anchor,free,point,f,maxfunc,maxdeep)
C+-----------------------------------------------------------------------+
C| Call the routine to initialise the mapping of x from the n-dimensional|
C| unit cube to the hypercube given by u and l. If an error occured,     |
C| give out a error message and return to the main program with the error|
C| flag set.                                                             |
C| JG 07/16/01 Changed call to remove unused data.                       |
C+-----------------------------------------------------------------------+
      CALL DIRpreprc(u,l,n,l,u,oops)
      IF (oops .GT. 0) THEN
        Write(*,10005)
        Write(logfile,10005)
        IError = -3
        Return
      END IF
      tstart = 2
C+-----------------------------------------------------------------------+
C| Initialise the algorithm DIRECT.                                      |
C+-----------------------------------------------------------------------+
C+-----------------------------------------------------------------------+
C| Added variable to keep track of the maximum value found.              |
C+-----------------------------------------------------------------------+
      CALL DIRInit(f,fcn,c,length,actdeep,point,anchor,free,
     +   dwrit,logfile,ArrayI,maxI,List2,w,x,l,u,fmin,minpos,
     +   thirds,levels,maxfunc,maxdeep,n,MaxDim,fmax,Ifeasiblef,
     +   IInfesiblef, Ierror,
     +   iidata, iisize, ddata, idsize, cdata, icsize)
C+-----------------------------------------------------------------------+
C| Added error checking.                                                 |
C+-----------------------------------------------------------------------+
      IF (Ierror .lt. 0) then
         IF (Ierror .eq. -4) THEN
            Write(*,10006)
            Write(logfile,10006)
            return
         END IF
         IF (Ierror .eq. -5) THEN
            Write(*,10007)
            Write(logfile,10007)
            return
         END IF
      END IF
      
      numfunc = 1 + maxI + maxI
      actmaxdeep = 1
      oldpos = 0
      tstart = 2
C+-----------------------------------------------------------------------+
C| If no feasible point has been found, give out the iteration, the      |
C| number of function evaluations and a warning. Otherwise, give out     |
C| the iteration, the number of function evaluations done and fmin.      |
C+-----------------------------------------------------------------------+
      IF (Ifeasiblef .gt. 0) then
        write(*,10012) tstart-1,numfunc
        write(logfile,10012) t,numfunc
      ELSE
        Write(*,10002) numfunc, fmin, fmax
        Write(logfile,10003) tstart-1,numfunc,fmin
      END IF
C+-----------------------------------------------------------------------+
C+-----------------------------------------------------------------------+
C| Main loop!                                                            |
C+-----------------------------------------------------------------------+
C+-----------------------------------------------------------------------+
      DO 10,t=tstart,MaxT
C+-----------------------------------------------------------------------+
C| Choose the sample points. The indices of the sample points are stored |
C| in the list S.                                                        |
C+-----------------------------------------------------------------------+
        actdeep = actmaxdeep
        CALL DIRChoose(anchor,S,maxdeep,f,fmin,eps,levels,maxpos,
     +     length,maxfunc,maxdeep,maxdiv,n,logfile,dwrit,cheat,kmax,
     +     Ifeasiblef)
C+-----------------------------------------------------------------------+
C| Add other hyperrectangles to S, which have the same level and the same|
C| function value at the center as the ones found above (that are stored |
C| in S). This is only done if we use the original DIRECT algorithm.     |
C| JG 07/16/01 Added Errorflag.                                          |
C+-----------------------------------------------------------------------+
        IF (algmethod .EQ. 0) THEN
          CALL DIRDoubleInsert(anchor, S, maxpos, point, f,
     +           maxdeep, maxfunc, maxdiv, Ierror)
          IF (Ierror .eq. -6) THEN
              Write(*,10020)
              Write(*,10021)
              Write(*,10022)
              Write(*,10023)
              Write(logfile,10020)
              Write(logfile,10021)
              Write(logfile,10022)
              Write(logfile,10023)
              return
          END IF
        ENDIF
        oldpos = minpos
C+-----------------------------------------------------------------------+
C| Initialise the number of sample points in this outer loop.            |
C+-----------------------------------------------------------------------+
        Newtosample = 0
        DO 20,j=1,maxpos
           actdeep = S(j,2)
C+-----------------------------------------------------------------------+
C| If the actual index is a point to sample, do it.                      |
C+-----------------------------------------------------------------------+
           IF (S(j,1) .GT. 0) THEN
C+-----------------------------------------------------------------------+
C| JG 09/24/00 Calculate the value delta used for sampling points.       |
C+-----------------------------------------------------------------------+
              actdeep_div = DIRGetmaxdeep(S(j,1),length,maxfunc,n)
              delta = thirds(actdeep_div+1)
              actdeep = S(j,2)
C+-----------------------------------------------------------------------+
C| If the current dept of division is only one under the maximal allowed |
C| dept, stop the computation.                                           |
C+-----------------------------------------------------------------------+
              IF (actdeep+1 .GE. mdeep) THEN
                 Write(*,10004)
                 write(logfile,10004)
                 Ierror = -6
                 GOTO 100
              END IF
              actmaxdeep = max(actdeep,actmaxdeep)
              help = S(j,1)
              IF (.NOT. (anchor(actdeep) .EQ. help)) THEN
                pos1 = anchor(actdeep)
                DO WHILE (.NOT. (point(pos1) .EQ. help))
                  pos1 = point(pos1)
                END DO
                point(pos1) = point(help)
              ELSE
                anchor(actdeep) = point(help)
              END IF
              IF (actdeep .lt. 0) then