dirsubrout.f

优化问题中的直接搜索法

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+-----------------------------------------------------------------------+
      SUBROUTINE DIRreplaceInf(free,freeold,f,c,thirds,length,anchor,
     +       point,c1,c2,maxfunc,maxdeep,maxdim,n,logfile, fmax)
      Implicit None
      Integer maxfunc, maxdeep, maxdim, n, free, freeold, logfile
      Double Precision  f(maxfunc,2)
      Integer anchor(-1:maxdeep)
      Integer point(maxfunc)
      Double Precision  c1(n),c2(n)
      Double Precision  c(maxfunc,MaxDim)
      Double Precision  thirds(0:maxdeep)
      Integer length(maxfunc,MaxDim)
      
      Integer LMaxDim
      PARAMETER (LMaxDim = 32)
      
      Double Precision sidelength
      Double Precision a(LmaxDim),b(LmaxDim),x(LmaxDim)
      Integer i,j,k,l, help, Isinbox
      Integer DIRgetmaxdeep
C+-----------------------------------------------------------------------+
C| JG 01/22/01 Added variable to keep track of the maximum value found.  |
C+-----------------------------------------------------------------------+
      Double Precision fmax
      
      DO 10, i = 1,free-1
         if (f(i,2) .gt. 0) then
C+-----------------------------------------------------------------------+
C| Get the maximum side length of the hyper rectangle and then set the   |
C| new side length to this lengths times the growth factor.              |
C+-----------------------------------------------------------------------+
            help = DIRgetmaxdeep(i,length,maxfunc,n)
            sidelength = thirds(help)*2.D0
C+-----------------------------------------------------------------------+
C| Set the Center and the upper and lower bounds of the rectangles.      |  
C+-----------------------------------------------------------------------+
            do 20, j = 1,n
               sidelength = thirds(length(i,j))               
               a(j) = c(i,j) - sidelength
               b(j) = c(i,j) + sidelength
20          continue
C+-----------------------------------------------------------------------+
C| The function value is reset to 'Inf', since it may have been changed  |
C| in an earlier iteration and now the feasible point which was close    |
C| is not close anymore (since the hyper rectangle surrounding the       |
C| current point may have shrunk).                                       |
C+-----------------------------------------------------------------------+            
            f(i,1) = 1.0E+6
            f(i,2) = 2.D0
C+-----------------------------------------------------------------------+
C| Check if any feasible point is near this infeasible point.            |
C+-----------------------------------------------------------------------+
            DO 30, k = 1,free-1
C+-----------------------------------------------------------------------+
C| If the point k is feasible, check if it is near.                      |
C+-----------------------------------------------------------------------+
               if (f(k,2) .eq. 0) then
C+-----------------------------------------------------------------------+
C| Copy the coordinates of the point k into x.                           |
C+-----------------------------------------------------------------------+
                  DO 40, l = 1,n
                     x(l) = c(k,l)
40                continue
C+-----------------------------------------------------------------------+
C| Check if the point k is near the infeasible point, if so, replace the |
C| value at 
C+-----------------------------------------------------------------------+
                  if (Isinbox(x,a,b,n,Lmaxdim) .eq. 1) then
                     f(i,1) = min(f(i,1), f(k,1))
                     f(i,2) = 1.D0
                  end if
               end if               
30          continue
            if (f(i,2) .eq. 1.0D0) then
               f(i,1) = f(i,1) + 1.0E-6*abs(f(i,1))
               do 200,l=1,n
                  x(l)=c(i,l)*c1(l)+c(i,l)*c2(l)
200            continue
               CALL DIRResortlist(i,anchor,f,point,length,n,maxfunc,
     +              maxdim,maxdeep,logfile)
            else
C+-----------------------------------------------------------------------+
C| JG 01/22/01                                                           |
C| Replaced fixed value for infeasible points with maximum value found,  |
C| increased by 1.                                                       |
C+-----------------------------------------------------------------------+
              if (.NOT. (fmax .eq. f(i,1))) then
                f(i,1) = max(fmax + 1.0D0,f(i,1))
              end if
            end if
         end if
10    continue
1000  format(20f18.8) 
      END 
      
C+-----------------------------------------------------------------------+
C| JG Added 09/25/00                                                     |
C|                                                                       |
C|                       SUBROUTINE DIRResortlist                        |
C|                                                                       |
C| Resort the list so that the infeasible point is in the list with the  |
C| replaced value.                                                       |
C+-----------------------------------------------------------------------+
      SUBROUTINE DIRResortlist(replace,anchor,f,point,length,n,maxfunc,
     +           maxdim,maxdeep,logfile)
      Implicit None
      Integer maxfunc, maxdim, maxdeep, n, logfile
      Integer replace
      Double Precision  f(maxfunc,2)
      Integer anchor(-1:maxdeep)
      Integer point(maxfunc)
      Integer length(maxfunc,MaxDim)
      Integer start, l, i, pos
      Integer DIRgetlevel
      
C+-----------------------------------------------------------------------+
C| Get the length of the hyper rectangle with infeasible mid point and   |
C| Index of the corresponding list.                                      |
C+-----------------------------------------------------------------------+
C JG 09/25/00 Replaced with DIRgetlevel
C      l = DIRgetmaxDeep(replace,length,maxfunc,n)
      l = DIRgetlevel(replace,length,maxfunc,n)
      start = anchor(l)
C+-----------------------------------------------------------------------+
C| If the hyper rectangle with infeasibel midpoint is already the start  |
C| of the list, give out message, nothing to do.                         |
C+-----------------------------------------------------------------------+
      if (replace .eq. start) then
C         write(logfile,*) 'No resorting of list necessarry, since new ',
C     + 'point is already anchor of list .',l
      else
C+-----------------------------------------------------------------------+
C| Take the hyper rectangle with infeasible midpoint out of the list.    |
C+-----------------------------------------------------------------------+
         pos = start
         do 10, i = 1,maxfunc
            if (point(pos) .eq. replace) then
               point(pos) = point(replace)
               goto 20
            else
               pos = point(pos)
            end if
            if (pos .eq. 0) then
               write(logfile,*) 'Error in DIRREsortlist: We went ',
     + 'through the whole list and could not find the point to 
     +  replace!!'
               goto 20
            end if
10       continue
C+-----------------------------------------------------------------------+
C| If the anchor of the list has a higher value than the value of a      |
C| nearby point, put the infeasible point at the beginning of the list.  |
C+-----------------------------------------------------------------------+
20       if (f(start,1) .gt. f(replace,1)) then
            anchor(l) = replace
            point(replace) = start
C            write(logfile,*) 'Point is replacing current anchor for '
C     +             , 'this list ',l,replace,start
         else
C+-----------------------------------------------------------------------+
C| Insert the point into the list according to its (replaced) function   |
C| value.                                                                |
C+-----------------------------------------------------------------------+
            pos = start
            do 30, i = 1,maxfunc
C+-----------------------------------------------------------------------+
C| The point has to be added at the end of the list.                     |
C+-----------------------------------------------------------------------+
               if (point(pos) .eq. 0) then
                  point(replace) = point(pos)
                  point(pos) = replace
C                  write(logfile,*) 'Point is added at the end of the '
C     +             , 'list ',l, replace
                  goto 40
               else
                  if (f(point(pos),1) .gt. f(replace,1)) then
                     point(replace) = point(pos)
                     point(pos) = replace
C                     write(logfile,*) 'There are points with a higher '
C     +               ,'f-value in the list ',l,replace, pos
                     goto 40
                  end if
                  pos = point(pos)
               end if
30          continue
40          pos = pos
         end if
      end if   
      end
C+-----------------------------------------------------------------------+
C|    SUBROUTINE DIRInsertList                                           |
C|    Changed 02-24-2000                                                 |
C|      Got rid of the distinction between feasible and infeasible points|
C|      I could do this since infeasible points get set to a high        |
C|      function value, which may be replaced by a function value of a   |
C|      nearby function at the end of the main loop.                     |
C+-----------------------------------------------------------------------+
      SUBROUTINE DIRInsertList(new,anchor,point,f,maxI,
     +                      length,maxfunc,maxdeep,n,samp)
      IMPLICIT None
      INTEGER maxfunc,maxdeep,j,maxI,n,samp
      INTEGER pos1,pos2,pos,new,deep,anchor(-1:maxdeep)
      INTEGER point(maxfunc),length(maxfunc,n)
C JG 09/24/00 Changed this to Getlevel
      INTEGER DIRGetlevel
      Double Precision  f(maxfunc,2)

      DO 10,j = 1,maxI
        pos1 = new
        pos2 = point(pos1)
        new = point(pos2)
C JG 09/24/00 Changed this to Getlevel
C        deep = DIRGetMaxdeep(pos1,length,maxfunc,n)
        deep = DIRGetlevel(pos1,length,maxfunc,n)
        IF (anchor(deep) .EQ. 0) THEN
           IF (f(pos2,1) .LT. f(pos1,1)) THEN
              anchor(deep) = pos2
              point(pos2) = pos1
              point(pos1) = 0
           ELSE
              anchor(deep) = pos1
              point(pos2) = 0
           END IF
        ELSE
           pos = anchor(deep)
           IF (f(pos2,1) .LT. f(pos1,1)) THEN
              IF (f(pos2,1) .LT. f(pos,1)) THEN
                 anchor(deep) = pos2
C JG 08/30/00 Fixed bug. Sorting was not correct when 
C      f(pos2,1) < f(pos1,1) < f(pos,1)
                 IF (f(pos1,1) .LT. f(pos,1)) THEN
                    point(pos2) = pos1
                    point(pos1) = pos
                 ELSE
                    point(pos2) = pos
                    CALL DIRInsert(pos,pos1,point,f,maxfunc)
                 END IF
              ELSE
                 CALL DIRInsert(pos,pos2,point,f,maxfunc)
                 CALL DIRInsert(pos,pos1,point,f,maxfunc)
              END IF
           ELSE
              IF (f(pos1,1) .LT. f(pos,1)) THEN
C JG 08/30/00 Fixed bug. Sorting was not correct when 
C      f(pos1,1) < f(pos2,1) < f(pos,1)
                 anchor(deep) = pos1
                 IF (f(pos,1) .LT. f(pos2,1)) THEN
                    point(pos1) = pos
                    CALL DIRInsert(pos,pos2,point,f,maxfunc)
                 ELSE
                    point(pos1) = pos2
                    point(pos2) = pos
                 END IF
              ELSE
                 CALL DIRInsert(pos,pos1,point,f,maxfunc)
                 CALL DIRInsert(pos,pos2,point,f,maxfunc)
              END IF              
           END IF
        END IF
10    CONTINUE
C JG 09/24/00 Changed this to Getlevel
C      deep = DIRGetMaxdeep(samp,length,maxfunc,n)
      deep = DIRGetlevel(samp,length,maxfunc,n)
      pos = anchor(deep)
      IF (f(samp,1) .LT. f(pos,1)) THEN
         anchor(deep) = samp
         point(samp) = pos
      ELSE
         CALL DIRInsert(pos,samp,point,f,maxfunc)
      END IF
      END

C+-----------------------------------------------------------------------+
C|    SUBROUTINE DIRInsertList2  (Old way to do it.)                     |
C+-----------------------------------------------------------------------+
      SUBROUTINE DIRInsertList_2(start,j,k,List2,w,maxI,n)
      IMPLICIT None
      INTEGER start,n,j,k
      INTEGER List2(n,2)
      Double Precision w(n)
      INTEGER pos,i,maxI

        pos = start
        IF (start .EQ. 0) THEN
          List2(j,1) = 0
          start = j
          GOTO 50
        END IF
        IF (w(start) .GT. w(j)) THEN
          List2(j,1) = start
          start = j
        ELSE
          DO 10,i=1,maxI
            IF (List2(pos,1) .EQ. 0) THEN
              List2(j,1) = 0
              List2(pos,1) =  j
              GOTO 50
            ELSE
              IF (w(j) .LT. w(List2(pos,1))) THEN
                List2(j,1) = List2(pos,1)