dirsubrout.f

优化问题中的直接搜索法

                List2(pos,1) = j
                GOTO 50
              END IF
            END IF
            pos = List2(pos,1)
10	   CONTINUE
        END IF
50     List2(j,2) = k 

      END

C+-----------------------------------------------------------------------+
C|    SUBROUTINE DIRSearchmin                                            |
C|    Search for the minimum in the list.                                !
C+-----------------------------------------------------------------------+
      SUBROUTINE DIRSearchmin(start,List2,pos,k,n)
      IMPLICIT None
      Integer start,pos,k,n
      INTEGER List2(n,2)

        k = start
        pos = List2(start,2)
        start = List2(start,1)
      END

C+-----------------------------------------------------------------------+
C|    SUBROUTINE DIRInit                                                 |
C|    Initialise all needed variables and do the first run of the        |
C|    algorithm.                                                         |
C|    Changed 02/24/2000                                                 |
C|       Changed fcn Double precision to fcn external!                   |
C|    Changed 09/15/2000                                                 |
C|       Added distinction between Jones way to characterize rectangles  |
C|       and our way. Common variable JONES controls which way we use.   |
C|          JONES = 0    Jones way (Distance from midpoint to corner)    |
C|          JONES = 1    Our way (Length of longest side)                |
C|    Changed 09/24/00                                                   |
C|       Added array levels. Levels contain the values to characterize   |
C|       the hyperrectangles.                                            |
C|    Changed 01/22/01                                                   |
C|       Added variable fmax to keep track of maximum value found.       |
C|       Added variable Ifeasiblef to keep track if feasibel point has   |
C|       been found.                                                     |
C|    Changed 01/23/01                                                   |
C|       Added variable Ierror to keep track of errors.                  |
C+-----------------------------------------------------------------------+
      SUBROUTINE 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,maxor,fmax,Ifeasiblef,IInfeasible, 
     + Ierror,
     + iidata, iisize, ddata, idsize, cdata, icsize)
      IMPLICIT None
      Integer maxfunc,maxdeep,n,maxor
      Double Precision  f(maxfunc,2),c(maxfunc,maxor),fmin
      Double Precision  x(n),delta, thirds(0:maxdeep)
      Double Precision  levels(0:maxdeep)
      Integer length(maxfunc,maxor),actdeep,minpos,i,oops,j
      Integer point(maxfunc),anchor(-1:maxdeep),free
      Integer ArrayI(maxor),maxI,new,dwrit,logfile,List2(maxor,2)
      Double Precision  w(maxor)
      External fcn
      Double Precision help2,l(n),u(n)
      Double Precision costmin
C+-----------------------------------------------------------------------+
C| JG 01/22/01 Added variable to keep track of the maximum value found.  |
C+-----------------------------------------------------------------------+
      Double Precision fmax
C+-----------------------------------------------------------------------+
C| JG 01/22/01 Added variable Ifeasiblef to keep track if feasibel point |
C|             has been found.                                           |
C| JG 01/23/01 Added variable Ierror to keep track of errors.            |
C| JG 03/09/01 Added IInfeasible to keep track if an infeasible point has|
C|             been found.                                               |
C+-----------------------------------------------------------------------+
      Integer Ifeasiblef, Ierror, IInfeasible
C JG 09/15/00 Added variable JONES (see above)
      Integer JONES, help
      COMMON /directcontrol/ JONES
C+-----------------------------------------------------------------------+
C| Variables to pass user defined data to the function to be optimized.  |
C+-----------------------------------------------------------------------+
      INTEGER iisize, idsize, icsize
      INTEGER iidata(iisize)
      Double Precision ddata(idsize)
      Character*40 cdata(icsize)
        
      fmin = 1.D20
      costmin = fmin
C JG 09/15/00 If Jones way of characterising rectangles is used, 
C             initialise thirds to reflect this.
      IF (JONES .eq. 0) THEN
        DO 5,j = 0,n-1
          w(j+1) = 0.5D0 * dsqrt(n - j + j/9.D0)
5       CONTINUE
        help2 = 1.D0
        DO 10,i = 1,maxdeep/n
          DO 8, j = 0, n-1
            levels((i-1)*n+j) = w(j+1) / help2
8         CONTINUE
          help2 = help2 * 3.D0
10      CONTINUE
      ELSE
C JG 09/15/00 Initialiase levels to contain 1/j
        help2 = 3.D0
        DO 11,i = 1,maxdeep
          levels(i) = 1.D0 / help2
          help2 = help2 * 3.D0
11      CONTINUE
        levels(0) = 1.D0
      ENDIF
      help2 = 3.D0
      DO 21,i = 1,maxdeep
        thirds(i) = 1.D0 / help2
        help2 = help2 * 3.D0
21    CONTINUE
      thirds(0) = 1.D0        
      DO 20,i=1,n
        c(1,i) = 0.5D0
        x(i) = 0.5D0
        length(1,i) = 0
20    CONTINUE
      CALL DIRinfcn(fcn,x,l,u,n,f(1,1),help,
     +                iidata, iisize, ddata, idsize, cdata, icsize)
      f(1,2) = help
      IInfeasible = help
      fmax = f(1,1)
C 09/25/00 Added this
C      if (f(1,1) .ge. 1.E+6) then
      if (f(1,2) .gt. 0.D0) then
        f(1,1) = 1.D6
        fmax = f(1,1)
        Ifeasiblef = 1
      else
        Ifeasiblef = 0
      end if

C JG 09/25/00 Remove IF
      fmin = f(1,1)
      costmin = f(1,1)
      minpos = 1
      actdeep = 2
      point(1) = 0
      free = 2
      delta = thirds(1)
      CALL DIRGet_I(length,1,ArrayI,maxI,n,maxfunc)
      new = free
      CALL DIRSamplepoints(c,ArrayI,delta,1,new,length,
     +           dwrit,logfile,f,free,maxI,point,fcn,x,l,
     +           fmin,minpos,u,n,
     +           maxfunc,maxdeep,oops)
C+-----------------------------------------------------------------------+
C| JG 01/23/01 Added error checking.                                     |
C+-----------------------------------------------------------------------+
      IF (oops .GT. 0) THEN
         IError = -4
         return
      END IF
C+-----------------------------------------------------------------------+
C| JG 01/22/01 Added variable to keep track of the maximum value found.  |
C|             Added variable to keep track if feasible point was found. |
C+-----------------------------------------------------------------------+
      CALL DIRSamplef(c,ArrayI,delta,1,new,length,
     +         dwrit,logfile,f,free,maxI,point,fcn,x,l,
     +         fmin,minpos,u,n,maxfunc,maxdeep,oops,fmax,Ifeasiblef,
     +         IInfeasible,
     +         iidata, iisize, ddata, idsize, cdata, icsize)
C+-----------------------------------------------------------------------+
C| JG 01/23/01 Added error checking.                                     |
C+-----------------------------------------------------------------------+
      IF (oops .GT. 0) THEN
         IError = -5
         return
      END IF
      CALL DIRDivide(new,0,length,point,
     +  ArrayI,1,List2,w,maxI,f,maxfunc,maxdeep,n)
      CALL DIRInsertList(new,anchor,point,f,maxI,length,
     +                    maxfunc,maxdeep,n,1)
      END

C+-----------------------------------------------------------------------+
C|    SUBROUTINE DIRDivide                                               |
C|    Subroutine to divide the hyper rectangles according to the rules.  |
C|    Changed 02-24-2000                                                 |
C|      Replaced if statement by min (line 367)                          |
C+-----------------------------------------------------------------------+
      SUBROUTINE DIRDivide(new,currentlength,length,point,
     +	   ArrayI,sample,List2,w,maxI,f,maxfunc,maxdeep,n)
      IMPLICIT None
      INTEGER start,new,maxfunc,maxdeep,n,sample
      INTEGER currentlength,length(maxfunc,n)
      INTEGER point(maxfunc)
      INTEGER List2(n,2),maxI,ArrayI(n)
      Double Precision f(maxfunc,2),w(n)
      INTEGER pos,i,j,k,pos2

        start = 0
        pos = new
        DO 10,i=1,maxI
          j = ArrayI(i)
          w(j) = f(pos,1)
          k = pos
          pos = point(pos)
          w(j) = min(f(pos,1),w(j))
          pos = point(pos)
          CALL DIRInsertList_2(start,j,k,list2,w,maxI,n)
10     CONTINUE
       IF (pos .GT. 0) THEN
           Write(*,*) "Error Divide"
           STOP
       END IF
       DO 20,j=1,maxI
         CALL DIRSearchmin(start,List2,pos,k,n)
         pos2 = start
         length(sample,k) = currentlength + 1
         DO 30,i=1,maxI-j+1
           length(pos,k) = currentlength + 1
           pos = point(pos)
           length(pos,k) = currentlength + 1
C JG 07/10/01 pos2 = 0 at the end of the 30-loop. Since we end
C             the loop now, we do not need to reassign pos and pos2.
           if (pos2 .gt. 0) then
             pos = List2(pos2,2)
             pos2 = List2(pos2,1)
           end if
30       CONTINUE
20     CONTINUE
      END

C+-----------------------------------------------------------------------+
C|    SUBROUTINE DIRSamplepoints                                         |
C|    Subroutine to sample the new points.                               |
C+-----------------------------------------------------------------------+
      SUBROUTINE DIRSamplepoints(c,ArrayI,delta,sample,start,length,
     +           dwrit,logfile,f,free,maxI,point,fcn,x,l,fmin,minpos,
     +           u,n,
     +           maxfunc,maxdeep,oops) 
      IMPLICIT None
      INTEGER n,maxfunc,maxdeep,oops
      INTEGER maxI,ArrayI(n),sample
      INTEGER length(maxfunc,n),free,point(maxfunc)
      Double Precision c(maxfunc,n),delta,x(n),l(n)
      Double Precision u(n),f(maxfunc,2)
      Double Precision fmin
      INTEGER pos,j,k,dwrit,logfile,minpos
      Integer start
      External fcn

      oops = 0
      pos = free
      start = free
      DO 10,k=1,maxI+maxI
        DO 20,j=1,n
          length(free,j) = length(sample,j)
          c(free,j) = c(sample,j)
20      CONTINUE
        pos = free
        free = point(free)
        IF (free .EQ. 0) THEN
           Write(*,1000) 
           Write(*,1001) 
           IF (dwrit .EQ. 2) THEN
             Write(logfile,1000)
             Write(logfile,1001)
           END IF
           oops = 1
           RETURN
1000  FORMAT("Error, no more free positions !")
1001  FORMAT("Increase maxfunc !")
        END IF
10    CONTINUE
      point(pos) = 0
      pos = start
      DO 30,j=1,maxI
         c(pos,ArrayI(j)) = c(sample,ArrayI(j)) + delta
         pos = point(pos)
         c(pos,ArrayI(j)) = c(sample,ArrayI(j)) - delta
         pos = point(pos)
30    CONTINUE
      IF (pos .GT. 0) THEN
          Write(*,2000)
          IF (dwrit .EQ. 2) THEN
             Write(logfile,2000)
           END IF
          STOP
2000      FORMAT("Error ! ") 
      END IF
      END


C+-----------------------------------------------------------------------+
C|    SUBROUTINE DIRGet_I                                                |
C+-----------------------------------------------------------------------+
      SUBROUTINE DIRGet_I(length,pos,ArrayI,maxi,n,maxfunc)
      IMPLICIT None
      Integer maxfunc,n,maxi,pos
      Integer length(maxfunc,n),ArrayI(n),i,help,j

      j = 1
      help = length(pos,1)
      DO 10,i = 2,n
        IF (length(pos,i) .LT. help) THEN
           help = length(pos,i)