chiark / gitweb /
Use trusty
[nlopt.git] / direct / DIRect.c
index f6862643491bb4da60b86980f83c93176d39e1c2..f0a7c1a95cd7f6bd311d3f26a1a80f0e21d729d3 100644 (file)
@@ -1,4 +1,4 @@
-/* DIRect.f -- translated by f2c (version 20050501).
+/* DIRect-transp.f -- translated by f2c (version 20050501).
    
    f2c output hand-cleaned by SGJ (August 2007). 
 */
@@ -44,7 +44,7 @@
 /* |   Lipschitz continues. However, DIRECT has proven to be effective on  | */
 /* |   more complex problems than these.                                   | */
 /* +-----------------------------------------------------------------------+ */
-/* Subroutine */ void direct_direct_(fp fcn, doublereal *x, integer *n, doublereal *eps, doublereal epsabs, integer *maxf, integer *maxt, doublereal *fmin, doublereal *l, 
+/* Subroutine */ void direct_direct_(fp fcn, doublereal *x, integer *n, doublereal *eps, doublereal epsabs, integer *maxf, integer *maxt, double starttime, double maxtime, int *force_stop, doublereal *minf, doublereal *l, 
        doublereal *u, integer *algmethod, integer *ierror, FILE *logfile, 
        doublereal *fglobal, doublereal *fglper, doublereal *volper, 
        doublereal *sigmaper, void *fcn_data)
     integer i__1, i__2;
     doublereal d__1;
 
-    /* constants (FIXME: change to variable?) */
-    const integer MAXFUNC = 90000;
-    const integer MAXDEEP = 600;
-    const integer MAXDIV = 3000;
+    /* changed by SGJ to be dynamically allocated ... would be
+       even better to use realloc, below, to grow these as needed */
+    integer MAXFUNC = *maxf <= 0 ? 101000 : (*maxf + 1000 + *maxf / 2);
+    integer MAXDEEP = *maxt <= 0 ? MAXFUNC/5: *maxt + 1000;
+    const integer MAXDIV = 5000;
 
     /* Local variables */
     integer increase;
     /* FIXME: change sizes dynamically? */
 #define MY_ALLOC(p, t, n) p = (t *) malloc(sizeof(t) * (n)); \
                           if (!(p)) { *ierror = -100; goto cleanup; }
+
+    /* Note that I've transposed c__, length, and f relative to the 
+       original Fortran code.  e.g. length was length(maxfunc,n) 
+       in Fortran [ or actually length(maxfunc, maxdims), but by
+       using malloc I can just allocate n ], corresponding to
+       length[n][maxfunc] in C, but I've changed the code to access
+       it as length[maxfunc][n].  That is, the maxfunc direction
+       is the discontiguous one.  This makes it easier to resize
+       dynamically (by adding contiguous rows) using realloc, without
+       having to move data around manually. */
     MY_ALLOC(c__, doublereal, MAXFUNC * (*n));
+    MY_ALLOC(length, integer, MAXFUNC * (*n));
     MY_ALLOC(f, doublereal, MAXFUNC * 2);
+    MY_ALLOC(point, integer, MAXFUNC);
+    if (*maxf <= 0) *maxf = MAXFUNC - 1000;
+
     MY_ALLOC(s, integer, MAXDIV * 2);
+
+    MY_ALLOC(anchor, integer, MAXDEEP + 2);
+    MY_ALLOC(levels, doublereal, MAXDEEP + 1);
+    MY_ALLOC(thirds, doublereal, MAXDEEP + 1);    
+    if (*maxt <= 0) *maxt = MAXDEEP;
+
     MY_ALLOC(w, doublereal, (*n));
     MY_ALLOC(oldl, doublereal, (*n));
     MY_ALLOC(oldu, doublereal, (*n));
     MY_ALLOC(list2, integer, (*n) * 2);
-    MY_ALLOC(point, integer, MAXFUNC);
-    MY_ALLOC(anchor, integer, MAXDEEP + 2);
-    MY_ALLOC(length, integer, MAXFUNC * (*n));
     MY_ALLOC(arrayi, integer, (*n));
-    MY_ALLOC(levels, doublereal, MAXDEEP + 1);
-    MY_ALLOC(thirds, doublereal, MAXDEEP + 1);    
 
 /* +-----------------------------------------------------------------------+ */
 /* |    SUBROUTINE Direct                                                  | */
 /* |     eps -- Exceeding value. If eps > 0, we use the same epsilon for   | */
 /* |            all iterations. If eps < 0, we use the update formula from | */
 /* |            Jones:                                                     | */
-/* |               eps = max(1.D-4*abs(fmin),epsfix),                      | */
+/* |               eps = max(1.D-4*abs(minf),epsfix),                      | */
 /* |            where epsfix = abs(eps), the absolute value of eps which is| */
 /* |            passed to the function.                                    | */
 /* |    maxf -- The maximum number of function evaluations.                | */
 /* |  fglper -- Terminate the optimization when the percent error          | */
 /* |                100(f_min - fglobal)/max(1,abs(fglobal)) < fglper.     | */
 /* |  volper -- Terminate the optimization when the volume of the          | */
-/* |            hyperrectangle S with f(c(S)) = fmin is less then volper   | */
+/* |            hyperrectangle S with f(c(S)) = minf is less then volper   | */
 /* |            percent of the volume of the original hyperrectangle.      | */
 /* |sigmaper -- Terminate the optimization when the measure of the         | */
-/* |            hyperrectangle S with f(c(S)) = fmin is less then sigmaper.| */
+/* |            hyperrectangle S with f(c(S)) = minf is less then sigmaper.| */
 /* |                                                                       | */
 /* | User data that is passed through without being changed:               | */
 /* |  fcn_data - opaque pointer to any user data                           | */
 /* |            X should be a good approximation to the global minimum     | */
 /* |            for the function within the hyper-box.                     | */
 /* |                                                                       | */
-/* |    fmin -- The value of the function at x.                            | */
+/* |    minf -- The value of the function at x.                            | */
 /* |  Ierror -- Error flag. If Ierror is lower 0, an error has occured. The| */
 /* |            values of Ierror mean                                      | */
 /* |            Fatal errors :                                             | */
 /* |              2   Number of iterations is equal to maxT.               | */
 /* |              3   The best function value found is within fglper of    | */
 /* |                  the (known) global optimum, that is                  | */
-/* |                   100(fmin - fglobal/max(1,|fglobal|))  < fglper.     | */
+/* |                   100(minf - fglobal/max(1,|fglobal|))  < fglper.     | */
 /* |                  Note that this termination signal only occurs when   | */
 /* |                  the global optimal value is known, that is, a test   | */
 /* |                  function is optimized.                               | */
-/* |              4   The volume of the hyperrectangle with fmin at its    | */
+/* |              4   The volume of the hyperrectangle with minf at its    | */
 /* |                  center is less than volper percent of the volume of  | */
 /* |                  the original hyperrectangle.                         | */
-/* |              5   The measure of the hyperrectangle with fmin at its   | */
+/* |              5   The measure of the hyperrectangle with minf at its   | */
 /* |                  center is less than sigmaper.                        | */
 /* |                                                                       | */
 /* | SUBROUTINEs used :                                                    | */
 /* +-----------------------------------------------------------------------+ */
     direct_dirinit_(f, fcn, c__, length, &actdeep, point, anchor, &ifree,
            logfile, arrayi, &maxi, list2, w, &x[1], &l[1], &u[1], 
-           fmin, &minpos, thirds, levels, &MAXFUNC, &MAXDEEP, n, n, &
-           fmax, &ifeasiblef, &iinfesiblef, ierror, fcn_data, jones);
+           minf, &minpos, thirds, levels, &MAXFUNC, &MAXDEEP, n, n, &
+           fmax, &ifeasiblef, &iinfesiblef, ierror, fcn_data, jones,
+                   starttime, maxtime, force_stop);
 /* +-----------------------------------------------------------------------+ */
 /* | Added error checking.                                                 | */
 /* +-----------------------------------------------------------------------+ */
                 fprintf(logfile, "WARNING: Error occured in routine DIRsamplef..\n");
            goto cleanup;
        }
+       if (*ierror == -102) goto L100;
     }
+    else if (*ierror == DIRECT_MAXTIME_EXCEEDED) goto L100;
     numfunc = maxi + 1 + maxi;
     actmaxdeep = 1;
     oldpos = 0;
 /* +-----------------------------------------------------------------------+ */
 /* | If no feasible point has been found, give out the iteration, the      | */
 /* | number of function evaluations and a warning. Otherwise, give out     | */
-/* | the iteration, the number of function evaluations done and fmin.      | */
+/* | the iteration, the number of function evaluations done and minf.      | */
 /* +-----------------------------------------------------------------------+ */
     if (ifeasiblef > 0) {
         if (logfile)
     } else {
         if (logfile)
              fprintf(logfile, "%d, %d, %g, %g\n", 
-                     tstart-1, numfunc, *fmin, fmax);
+                     tstart-1, numfunc, *minf, fmax);
     }
 /* +-----------------------------------------------------------------------+ */
 /* +-----------------------------------------------------------------------+ */
 /* | in the list S.                                                        | */
 /* +-----------------------------------------------------------------------+ */
        actdeep = actmaxdeep;
-       direct_dirchoose_(anchor, s, &MAXDEEP, f, fmin, *eps, epsabs, levels, &maxpos, length, 
+       direct_dirchoose_(anchor, s, &MAXDEEP, f, minf, *eps, epsabs, levels, &maxpos, length, 
                &MAXFUNC, &MAXDEEP, &MAXDIV, n, logfile, &cheat, &
                kmax, &ifeasiblef, jones);
 /* +-----------------------------------------------------------------------+ */
        newtosample = 0;
        i__2 = maxpos;
        for (j = 1; j <= i__2; ++j) {
-           actdeep = s[j + 2999];
+           actdeep = s[j + MAXDIV-1];
 /* +-----------------------------------------------------------------------+ */
 /* | If the actual index is a point to sample, do it.                      | */
 /* +-----------------------------------------------------------------------+ */
                actdeep_div__ = direct_dirgetmaxdeep_(&s[j - 1], length, &MAXFUNC, 
                        n);
                delta = thirds[actdeep_div__ + 1];
-               actdeep = s[j + 2999];
+               actdeep = s[j + MAXDIV-1];
 /* +-----------------------------------------------------------------------+ */
 /* | If the current dept of division is only one under the maximal allowed | */
 /* | dept, stop the computation.                                           | */
                    anchor[actdeep + 1] = point[help - 1];
                }
                if (actdeep < 0) {
-                   actdeep = (integer) f[help - 1];
+                   actdeep = (integer) f[(help << 1) - 2];
                }
 /* +-----------------------------------------------------------------------+ */
 /* | Get the Directions in which to decrease the intervall-length.         | */
 /* +-----------------------------------------------------------------------+ */
                direct_dirsamplepoints_(c__, arrayi, &delta, &help, &start, length, 
                        logfile, f, &ifree, &maxi, point, &x[
-                       1], &l[1], fmin, &minpos, &u[1], n, &MAXFUNC, &
+                       1], &l[1], minf, &minpos, &u[1], n, &MAXFUNC, &
                        MAXDEEP, &oops);
                if (oops > 0) {
                    if (logfile)
 /* +-----------------------------------------------------------------------+ */
                direct_dirsamplef_(c__, arrayi, &delta, &help, &start, length,
                            logfile, f, &ifree, &maxi, point, fcn, &x[
-                       1], &l[1], fmin, &minpos, &u[1], n, &MAXFUNC, &
+                       1], &l[1], minf, &minpos, &u[1], n, &MAXFUNC, &
                        MAXDEEP, &oops, &fmax, &ifeasiblef, &iinfesiblef, 
-                       fcn_data);
+                                  fcn_data, force_stop);
+               if (force_stop && *force_stop) {
+                    *ierror = -102;
+                    goto L100;
+               }
+               if (nlopt_stop_time_(starttime, maxtime)) {
+                    *ierror = DIRECT_MAXTIME_EXCEEDED;
+                    goto L100;
+               }
                if (oops > 0) {
                    if (logfile)
                         fprintf(logfile, "WARNING: Error occured in routine DIRsamplef.\n");
        if (oldpos < minpos) {
            if (logfile)
                 fprintf(logfile, "%d, %d, %g, %g\n",
-                        t, numfunc, *fmin, fmax);
+                        t, numfunc, *minf, fmax);
        }
 /* +-----------------------------------------------------------------------+ */
 /* | If no feasible point has been found, give out the iteration, the      | */
 /* +-----------------------------------------------------------------------+ */
 /* +-----------------------------------------------------------------------+ */
 /* | JG 01/22/01 Calculate the index for the hyperrectangle at which       | */
-/* |             fmin is assumed. We then calculate the volume of this     | */
+/* |             minf is assumed. We then calculate the volume of this     | */
 /* |             hyperrectangle and store it in delta. This delta can be   | */
 /* |             used to stop DIRECT once the volume is below a certain    | */
 /* |             percentage of the original volume. Since the original     | */
        }
 /* +-----------------------------------------------------------------------+ */
 /* | JG 01/23/01 Calculate the measure for the hyperrectangle at which     | */
-/* |             fmin is assumed. If this measure is smaller then sigmaper,| */
+/* |             minf is assumed. If this measure is smaller then sigmaper,| */
 /* |             we stop DIRECT.                                           | */
 /* +-----------------------------------------------------------------------+ */
        actdeep_div__ = direct_dirgetlevel_(&minpos, length, &MAXFUNC, n, jones);
 /* | global minimum value, terminate. This only makes sense if this optimal| */
 /* | value is known, that is, in test problems.                            | */
 /* +-----------------------------------------------------------------------+ */
-       if ((*fmin - *fglobal) * 100 / divfactor <= *fglper) {
+       if ((*minf - *fglobal) * 100 / divfactor <= *fglper) {
            *ierror = 3;
            if (logfile)
-                fprintf(logfile, "DIRECT stopped: fmin within fglper of global minimum.\n");
+                fprintf(logfile, "DIRECT stopped: minf within fglper of global minimum.\n");
            goto L100;
        }
 /* +-----------------------------------------------------------------------+ */
 /* +-----------------------------------------------------------------------+ */
        if (iepschange == 1) {
 /* Computing MAX */
-           d__1 = fabs(*fmin) * 1e-4;
+           d__1 = fabs(*minf) * 1e-4;
            *eps = MAX(d__1,epsfix);
        }
 /* +-----------------------------------------------------------------------+ */
@@ -706,7 +733,7 @@ L100:
 /* +-----------------------------------------------------------------------+ */
     i__1 = *n;
     for (i__ = 1; i__ <= i__1; ++i__) {
-       x[i__] = c__[minpos + i__ * 90000 - 90001] * l[i__] + l[i__] * u[i__];
+       x[i__] = c__[i__ + minpos * i__1 - i__1-1] * l[i__] + l[i__] * u[i__];
        u[i__] = oldu[i__ - 1];
        l[i__] = oldl[i__ - 1];
 /* L50: */
@@ -718,7 +745,7 @@ L100:
 /* +-----------------------------------------------------------------------+ */
 /* | Give out a summary of the run.                                        | */
 /* +-----------------------------------------------------------------------+ */
-    direct_dirsummary_(logfile, &x[1], &l[1], &u[1], n, fmin, fglobal, &numfunc, 
+    direct_dirsummary_(logfile, &x[1], &l[1], &u[1], n, minf, fglobal, &numfunc, 
            ierror);
 /* +-----------------------------------------------------------------------+ */
 /* | Format statements.                                                    | */