direct.f

优化问题中的直接搜索法

                actdeep = f(help,1)
              END IF
C+-----------------------------------------------------------------------+
C| Get the Directions in which to decrease the intervall-length.         |
C+-----------------------------------------------------------------------+
              CALL DIRGet_I(length,help,ArrayI,maxI,n,maxfunc)
C+-----------------------------------------------------------------------+
C| Sample the function. To do this, we first calculate the points where  |
C| we need to sample the function. After checking for errors, we then do |
C| the actual evaluation of the function, again followed by checking for |
C| errors.                                                               |
C+-----------------------------------------------------------------------+
              CALL DIRSamplepoints(c,ArrayI,delta,help,
     +             start,length,dwrit,logfile,f,free,maxI,point,
     +             fcn,x,l,fmin,minpos,u,n,maxfunc,maxdeep,oops)
              IF (oops .GT. 0) THEN
                Write(*,10006)
                Write(logfile,10006)
                IError = -4
                return
              END IF
              Newtosample = newtosample + maxI
C+-----------------------------------------------------------------------+
C| JG 01/22/01 Added variable to keep track of the maximum value found.  |
C+-----------------------------------------------------------------------+
              CALL DIRSamplef(c,ArrayI,delta,help,start,
     +            length,dwrit,logfile,f,free,maxI,point,fcn,x,l,
     +            fmin,minpos,u,n,maxfunc,maxdeep,oops,fmax,
     +            Ifeasiblef,IInfesiblef,
     +            iidata, iisize, ddata, idsize, cdata, icsize)
              IF (oops .GT. 0) THEN
                Write(*,10007)
                Write(logfile,10007)
                IError = -5
                return
              END IF
C+-----------------------------------------------------------------------+
C| Divide the intervalls.                                                |
C+-----------------------------------------------------------------------+
              CALL DIRDivide(start,actdeep_div,length,point,
     +             ArrayI,help,List2,w,maxI,f,maxfunc,maxdeep,n)
C+-----------------------------------------------------------------------+
C| Insert the new intervalls into the list (sorted).                     |
C+-----------------------------------------------------------------------+
              CALL DIRInsertList(start,anchor,point,f,maxI,length,
     +                    maxfunc,maxdeep,n,help)
C+-----------------------------------------------------------------------+
C| Increase the number of function evaluations.                          |
C+-----------------------------------------------------------------------+
              numfunc = numfunc + maxI + maxI
           END IF
C+-----------------------------------------------------------------------+
C| End of main loop.                                                     |
C+-----------------------------------------------------------------------+
20      CONTINUE
C+-----------------------------------------------------------------------+
C| If there is a new minimum, show the actual iteration, the number of   |
C| function evaluations, the minimum value of f (so far) and the position|
C| in the array.                                                         |
C+-----------------------------------------------------------------------+
        IF (oldpos .LT. minpos) THEN
          Write(*,10002) numfunc,fmin, fmax
          Write(logfile,10003) t,numfunc,fmin
        END IF
C+-----------------------------------------------------------------------+
C| If no feasible point has been found, give out the iteration, the      |
C| number of function evaluations and a warning.                         |
C+-----------------------------------------------------------------------+
        IF (Ifeasiblef .gt. 0) then
          write(*,10012) t,numfunc
          write(logfile,10012) t,numfunc
        END IF
C+-----------------------------------------------------------------------+
C+-----------------------------------------------------------------------+
C|                       Termination Checks                              |
C+-----------------------------------------------------------------------+

C+-----------------------------------------------------------------------+
C| JG 01/22/01 Calculate the index for the hyperrectangle at which       |
C|             fmin is assumed. We then calculate the volume of this     |
C|             hyperrectangle and store it in delta. This delta can be   |
C|             used to stop DIRECT once the volume is below a certain    |
C|             percentage of the original volume. Since the original     |
C|             is 1 (scaled), we can stop once delta is below a certain  |
C|             percentage, given by volper.                              |
C+-----------------------------------------------------------------------+
        Ierror = Jones
        Jones = 0
        actdeep_div = DIRGetlevel(minpos,length,maxfunc,n)
        Jones = Ierror
C+-----------------------------------------------------------------------+
C| JG 07/16/01 Use precalculated values to calculate volume.             |
C+-----------------------------------------------------------------------+
        delta = thirds(actdeep_div)*100
        IF (delta .LE. volper) THEN
           Ierror = 4
           Write(*,10011) delta, volper
           Write(logfile,10011) delta, volper
           GOTO 100
        END IF
C+-----------------------------------------------------------------------+
C| JG 01/23/01 Calculate the measure for the hyperrectangle at which     |
C|             fmin is assumed. If this measure is smaller then sigmaper,|
C|             we stop DIRECT.                                           |
C+-----------------------------------------------------------------------+
        actdeep_div = DIRGetlevel(minpos,length,maxfunc,n)
        delta = levels(actdeep_div)
        IF (delta .LE. sigmaper) THEN
           Ierror = 5
           Write(*,10013) delta, sigmaper
           Write(logfile,10013) delta, sigmaper
           GOTO 100
        END IF
C+-----------------------------------------------------------------------+
C| If the best found function value is within fglper of the (known)      |
C| global minimum value, terminate. This only makes sense if this optimal|
C| value is known, that is, in test problems.                            |
C+-----------------------------------------------------------------------+
        IF ((100*(fmin - fglobal)/divfactor) .LE. fglper)
     +   THEN
           Ierror = 3
           Write(*,10010)
           Write(logfile,10010)
           GOTO 100
        END IF
C+-----------------------------------------------------------------------+
C| Find out if there are infeasible points which are near feasible ones. |
C| If this is the case, replace the function value at the center of the  |
C| hyper rectangle by the lowest function value of a nearby function.    |
C| If no infeasible points exist (IInfesiblef = 0), skip this.           |
C+-----------------------------------------------------------------------+
        IF (IInfesiblef .gt. 0) THEN
           CALL DIRreplaceInf(free,freeold,f,c,thirds,length,anchor,
     +       point,u,l,maxfunc,maxdeep,maxdim,n,logfile, fmax)
        ENDIF
        freeold = free
C+-----------------------------------------------------------------------+
C| If iepschange = 1, we use the epsilon change formula from Jones.      |
C+-----------------------------------------------------------------------+
        IF (iepschange .eq. 1) then
           eps = max(1.D-4*abs(fmin),epsfix)
        END IF
C+-----------------------------------------------------------------------+
C| If no feasible point has been found yet, set the maximum number of    |
C| function evaluations to the number of evaluations already done plus   |
C| the budget given by the user.                                         |
C| If the budget has already be increased, increase it again. If a       |
C| feasible point has been found, remark that and reset flag. No further |
C| increase is needed.                                                   |
C+-----------------------------------------------------------------------+
        IF (increase .eq. 1) then        
           maxf = numfunc + oldmaxf
           IF (Ifeasiblef .eq. 0) then
             write(logfile,10031) maxf
             increase = 0
           END IF
        END IF
C+-----------------------------------------------------------------------+
C| Check if the number of function evaluations done is larger than the   |
C| allocated budget. If this is the case, check if a feasible point was  |
C| found. If this is a case, terminate. If no feasible point was found,  |
C| increase the budget and set flag increase.                            |
C+-----------------------------------------------------------------------+
        IF (numfunc .GT. maxf) THEN
           IF (Ifeasiblef .eq. 0) then
              Ierror = 1
              Write(*,10008)
              Write(logfile,10008)
              GOTO 100
           ELSE
              increase = 1
              write(logfile,10030) numfunc
              maxf = numfunc+ oldmaxf
           END IF
        END IF
10    CONTINUE
C+-----------------------------------------------------------------------+
C+-----------------------------------------------------------------------+
C| End of main loop.                                                     |
C+-----------------------------------------------------------------------+
C+-----------------------------------------------------------------------+

C+-----------------------------------------------------------------------+
C| The algorithm stopped after maxT iterations.                          |
C+-----------------------------------------------------------------------+
      Ierror = 2
      Write(*,10009)
      Write(logfile,10009)

100   CONTINUE
C+-----------------------------------------------------------------------+
C| Store the position of the minimum in x.                               |
C+-----------------------------------------------------------------------+
      DO 50,i=1,n
         x(i) = c(Minpos,i)*l(i)+l(i)*u(i)
         u(i) = oldu(i)
         l(i) = oldl(i)
50    CONTINUE
C+-----------------------------------------------------------------------+
C| Store the number of function evaluations in maxf.                     |
C+-----------------------------------------------------------------------+
      maxf = numfunc

C+-----------------------------------------------------------------------+
C| If needed, save the final division in a file for use with Matlab.     |
C+-----------------------------------------------------------------------+
      writed = 0
      IF (writed .EQ. 1) THEN
        file  = 12
        file2 =  0
        CALL DIRWritehistbox(point,f,thirds,c,anchor,actdeep,file,
     +               l,u,file2,maxfunc,maxdeep,n,MaxDim,length)
      END IF

C+-----------------------------------------------------------------------+
C| Give out a summary of the run.                                        |
C+-----------------------------------------------------------------------+
      CALL DIRsummary(logfile,x,l,u,n,fmin,fglobal, numfunc, Ierror)

C+-----------------------------------------------------------------------+
C| Format statements.                                                    |
C+-----------------------------------------------------------------------+
10002  FORMAT(i5," & ",f18.10," & ",f18.10," \\\\ ")
10003  FORMAT(i5,'       ', i5,'     ',f18.10)
10004  FORMAT('WARNING : Maximum number of levels reached. Increase
     +        maxdeep.')
10005  FORMAT('WARNING : Initialisation in DIRpreprc failed.')
10006  FORMAT('WARNING : Error occured in routine DIRsamplepoints.')
10007  FORMAT('WARNING : Error occured in routine DIRsamplef.')
10008  FORMAT('DIRECT stopped: numfunc >= maxf.')
10009  FORMAT('DIRECT stopped: maxT iterations.')
10010  FORMAT('DIRECT stopped: fmin within fglper of global minimum.')
10011  FORMAT('DIRECT stopped: Volume of S_min is ',d8.2,
     +       '% < ',d8.2,'% of the original volume.')
10012  FORMAT('No feasible point found in ',I4,' iterations ',
     +       'and ',I5,' function evaluations.')
10013  FORMAT('DIRECT stopped: Measure of S_min = ',d8.2,' < '
     +       ,d8.2,'.')
10020  FORMAT('WARNING : Capacity of array S in DIRDoubleInsert'
     +        ' reached. Increase maxdiv.')
10021  FORMAT('This means that there are a lot of hyperrectangles')
10022  FORMAT('with the same function value at the center. We')
10023  FORMAT('suggest to use our modification instead (Jones = 1)')

10030 FORMAT('DIRECT could not find a feasible ',
     +        'point after ', I5, ' function evaluations. ',
     +        'DIRECT continues until a feasible point is found.')
10031 FORMAT('DIRECT found a feasible point. ',
     +         'The adjusted budget is now set to ',I5,'.')
      END