dirsubrout.f

优化问题中的直接搜索法

C+-----------------------------------------------------------------------+
C| Program       : Direct.f (subfile DIRsubrout.f)                       |
C| Last modified : 07-16-2001                                            |
C| Written by    : Joerg Gablonsky                                       |
C| Subroutines used by the algorithm DIRECT.                             |
C+-----------------------------------------------------------------------+

C+-----------------------------------------------------------------------+
C|    SUBROUTINE DIRChoose                                               |
C|    Decide, which is the next sampling point.                          |
C|    Changed 09/25/00 JG                                                |
C|         Added maxdiv to call and changed S to size maxdiv.            |
C|    Changed 01/22/01 JG                                                |
C|         Added Ifeasiblef to call to keep track if a feasible point has|
C|         been found.                                                   |
C|    Changed 07/16/01 JG                                                |
C|         Changed if statement to prevent run-time errors.              |                                                   |
C+-----------------------------------------------------------------------+
      SUBROUTINE DIRChoose(anchor,S,actdeep,f,fmin,eps,thirds,maxpos,
     +      length,maxfunc,maxdeep,maxdiv,n,logfile,dwrit,cheat,kmax,
     +      Ifeasiblef)

      IMPLICIT None
      Double Precision MaxLower
      Parameter (MaxLower = 1.D20)
      Integer maxfunc,maxdeep,n, maxdiv
      Integer Anchor(-1:maxdeep),S(maxdiv,2),length(maxfunc,n)
      Integer actdeep,dwrit
      Integer DIRGetlevel
      Double Precision f(maxfunc,2),fmin,eps,thirds(0:maxfunc)
      Integer maxpos,i,j,k,i_,j_,logfile,cheat
      Double Precision help2,helplower,helpgreater,kmax
      Integer novalue,novaluedeep
      Integer Ifeasiblef

        helplower = MaxLower
        helpgreater = 0.D0
        k = 1
        if (Ifeasiblef .ge. 1) then
           DO 1001,j=0,actdeep
              IF (anchor(j) .GT. 0) THEN
                 S(k,1) = anchor(j)
                 S(k,2) = DIRGetlevel(S(k,1),length,maxfunc,n)
                 goto 12
              END IF
1001       CONTINUE
12         k = k + 1
           maxpos = 1
           return
        else
           DO 10,j=0,actdeep
              IF (anchor(j) .GT. 0) THEN
                 S(k,1) = anchor(j)
                 S(k,2) = DIRGetlevel(S(k,1),length,maxfunc,n)
                 k = k + 1
              END IF
10         CONTINUE
        END IF
        novalue = 0
        if (anchor(-1) .gt. 0) then
          novalue = anchor(-1)
          novaluedeep = DIRGetlevel(novalue,length,maxfunc,n)
        end if
        maxpos = k - 1
       DO 11,j=k-1,maxdeep
         S(k,1) = 0
11     CONTINUE
       DO 40,j=maxpos,1,-1
         helplower = Maxlower
         helpgreater = 0.D0
         j_ = S(j,1)
         DO 30,i=1,j-1
           i_ = S(i,1)
C+-----------------------------------------------------------------------+
C| JG 07/16/01 Changed IF statement into two to prevent run-time errors  |
C|             which could occur if the compiler checks the second       |
C|             expression in an .AND. statement although the first       |
C|             statement is already not true.                            |
C+-----------------------------------------------------------------------+
           IF ((i_ .GT. 0) .AND. .NOT. (i .EQ. j)) THEN
             IF (f(i_,2) .le. 1.D0) THEN
               help2 = thirds(S(i,2)) - thirds(S(j,2))
               help2 = (f(i_,1) - f(j_,1))/help2
               IF (help2 .LE. 0.D0) THEN
                 IF (dwrit .EQ. 2) THEN
                   Write(logfile,*) "thirds > 0,help2 <= 0"
                 END IF
                 GOTO 60
               END IF
               IF (help2 .LT. helplower) THEN
                 IF (dwrit .EQ. 2) THEN
                   Write(logfile,*) "helplower = ",help2
                 END IF
                 helplower = help2
               END IF
         	   END IF
           END IF
30       CONTINUE	     
         DO 31,i=j+1,maxpos
           i_ = S(i,1)
C+-----------------------------------------------------------------------+
C| JG 07/16/01 Changed IF statement into two to prevent run-time errors  |
C|             which could occur if the compiler checks the second       |
C|             expression in an .AND. statement although the first       |
C|             statement is already not true.                            |
C+-----------------------------------------------------------------------+
           IF ((i_ .GT. 0) .AND. .NOT. (i .EQ. j)) THEN
             IF (f(i_,2) .le. 1.D0) THEN 
               help2 = thirds(S(i,2)) - thirds(S(j,2))
               help2 = (f(i_,1) - f(j_,1))/help2
               IF (help2 .LE. 0.D0) THEN
                 IF (dwrit .EQ. 2) THEN
                   Write(logfile,*) "thirds < 0,help2 <= 0"
                 END IF
                 GOTO 60
               END IF
               IF (help2 .GT. helpgreater) THEN
                 IF (dwrit .EQ. 2) THEN
                   Write(logfile,*) "helpgreater = ",help2
                 END IF
                 helpgreater = help2
               END IF
             END IF
           END IF
31       CONTINUE
         IF ((helplower .GT. Maxlower) .AND. 
     +       (helpgreater .gt. 0)) THEN
           helplower = helpgreater
           helpgreater = helpgreater - 1.D0
         END IF
         IF (helpgreater .LE. helplower) THEN
           IF ((cheat .EQ. 1) .AND. (helplower .GT. Kmax)) THEN
             helplower = Kmax
           END IF
           IF ((f(j_,1) - helplower * thirds(S(j,2))) .GT. 
     +        (fmin - eps*abs(fmin))) THEN
              IF (dwrit .EQ. 2) THEN
                Write(logfile,*) "> fmin - eps|fmin|"
              END IF
              GOTO 60
            END IF
         ELSE
           IF (dwrit .EQ. 2) THEN
           Write(logfile,*) "helpgreater > helplower",helpgreater,
     +            helplower,helpgreater - helplower
           END IF
           GOTO 60
         END IF
         GOTO 40
60       S(j,1) = 0
40     CONTINUE
       if (novalue .gt. 0) then
          maxpos = maxpos + 1
          S(maxpos,1) = novalue
          S(maxpos,2) = novaluedeep
       end if
      END

C+-----------------------------------------------------------------------+
C| INTEGER Function GetmaxDeep                                           |
C| function to get the maximal length (1/length) of the n-dimensional    |
C| rectangle with midpoint pos.                                          |
C|                                                                       |
C| On Return :                                                           |
C|    the maximal length                                                 |
C|                                                                       |
C| pos     -- the position of the midpoint in the array length           |
C| length  -- the array with the dimensions                              |
C| maxfunc -- the leading dimension of length                            |
C| n	   -- the dimension of the problem                                  |
C|                                                                       |
C+-----------------------------------------------------------------------+

      INTEGER Function DIRGetmaxDeep(pos,length,maxfunc,n)
      IMPLICIT None
      INTEGER pos,maxfunc,n,length(maxfunc,n),help,i

      help = length(pos,1)
      DO 10,i = 2,n
        help = min(help, length(pos,i))
10    CONTINUE            
      DIRGetMaxDeep = help
      END

C+-----------------------------------------------------------------------+
C| INTEGER Function DIRGetlevel                                          |
C| Returns the level of the hyperrectangle. Depending on the value of the|
C| global variable JONES. IF JONES equals 0, the level is given by       |
C|               kN + p, where the rectangle has p sides with a length of|
C|             1/3^(k+1), and N-p sides with a length of 1/3^k.          |
C| If JONES equals 1, the level is the power of 1/3 of the length of the |
C| longest side hyperrectangle.                                          |
C|                                                                       |
C| On Return :                                                           |
C|    the maximal length                                                 |
C|                                                                       |
C| pos     -- the position of the midpoint in the array length           |
C| length  -- the array with the dimensions                              |
C| maxfunc -- the leading dimension of length                            |
C| n	   -- the dimension of the problem                                  |
C|                                                                       |
C+-----------------------------------------------------------------------+

      INTEGER Function DIRGetlevel(pos,length,maxfunc,n)
      IMPLICIT None
      INTEGER pos,maxfunc,n,length(maxfunc,n),help,i,p,k
C JG 09/15/00 Added variable JONES (see above)
      Integer JONES
      COMMON /directcontrol/ JONES

      IF (JONES .eq. 0) THEN
        help = length(pos,1)
        k = help
        p = 1
        DO 100,i = 2,n
          IF (length(pos,i) .LT. k) THEN
            k = length(pos,i)
          END IF
          IF (length(pos,i) .EQ. help) THEN
             p = p + 1
          END IF
100     CONTINUE
        IF (k .EQ. help) THEN            
           DIRGetLevel = k*n + n-p
        ELSE
           DIRGetLevel = k*n + p
        END IF
      ELSE
        help = length(pos,1)
        DO 10,i = 2,n
          IF (length(pos,i) .LT. help) THEN
            help = length(pos,i)
          END IF
10      CONTINUE            
        DIRGetLevel = help
      END IF
      END

C+-----------------------------------------------------------------------+
C|    SUBROUTINE DIRDoubleInsert                                         |
C|      Routine to make sure that if there are several potential optimal |
C|      hyperrectangles of the same level (i.e. hyperrectangles that have|
C|      the same level and the same function value at the center), all of|
C|      them are divided. This is the way as originally described in     |
C|      Jones et.al.                                                     |
C| JG 07/16/01 Added errorflag to calling sequence. We check if more     |
C|             we reach the capacity of the array S. If this happens, we |
C|             return to the main program with an error.                 |
C+-----------------------------------------------------------------------+
      SUBROUTINE DIRDoubleInsert(anchor, S, maxpos, point, f, 
     +           maxdeep, maxfunc, maxdiv, IError)
      IMPLICIT None
      Integer maxpos, maxdeep, maxfunc, maxdiv
      Integer anchor(-1:maxdeep),S(maxdiv,2)
      Integer point(maxfunc)
      Double Precision  f(maxfunc,2)
      
C+-----------------------------------------------------------------------+
C| JG 07/16/01 Added flag to prevent run time-errors on some systems.    |
C+-----------------------------------------------------------------------+
      Integer iflag, IError
      Integer oldmaxpos, i, pos, help, actdeep
      
      oldmaxpos = maxpos
      DO 10, i = 1,oldmaxpos
        IF (S(i,1) .GT. 0) THEN
          actdeep = S(i,2)
          help = anchor(actdeep)
          pos = point(help)
          iflag = 0
C+-----------------------------------------------------------------------+
C| JG 07/16/01 Added flag to prevent run time-errors on some systems. On |
C|             some systems the second conditions in an AND statement is |
C|             evaluated even if the first one is already not true.      |
C+-----------------------------------------------------------------------+
          DO WHILE ((pos .GT. 0) .AND. (iflag .eq. 0))
             if (f(pos,1) - f(help,1) .LE. 1.D-13) then
               if (maxpos .lt. maxdiv) then
                 maxpos = maxpos + 1
                 S(maxpos,1) = pos
                 S(maxpos,2) = actdeep
                 pos = point(pos)
               else
C+-----------------------------------------------------------------------+
C| JG 07/16/01 Maximum number of elements possible in S has been reached!|
C+-----------------------------------------------------------------------+
                 IError = -6
                 return
               end if
             else
               iflag = 1
             end if
          END DO
        END IF 
10    CONTINUE
      END

C+-----------------------------------------------------------------------+
C| JG Added 09/25/00                                                     |
C|                       SUBROUTINE DIRreplaceInf                        |
C|                                                                       |
C| Find out if there are infeasible points which are near feasible ones. |