chiark / gitweb /
Use trusty
[nlopt.git] / direct / DIRect.c
index f2083328b428298133f98bafc180c3d9c771bdbb..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). 
 */
 
 /* Table of constant values */
 
-static integer c__90000 = 90000;
-static integer c__600 = 600;
-static integer c__64 = 64;
-static integer c__3000 = 3000;
-
 /* +-----------------------------------------------------------------------+ */
 /* | Program       : Direct.f                                              | */
 /* | Last modified : 07-16-2001                                            | */
@@ -49,7 +44,7 @@ static integer c__3000 = 3000;
 /* |   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)
@@ -58,6 +53,12 @@ static integer c__3000 = 3000;
     integer i__1, i__2;
     doublereal d__1;
 
+    /* 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;
     doublereal *c__ = 0        /* was [90000][64] */, *f = 0   /* 
@@ -86,19 +87,34 @@ static integer c__3000 = 3000;
     /* FIXME: change sizes dynamically? */
 #define MY_ALLOC(p, t, n) p = (t *) malloc(sizeof(t) * (n)); \
                           if (!(p)) { *ierror = -100; goto cleanup; }
-    MY_ALLOC(c__, doublereal, 5760000);
-    MY_ALLOC(f, doublereal, 180000);
-    MY_ALLOC(s, integer, 6000);
-    MY_ALLOC(w, doublereal, 64);
-    MY_ALLOC(oldl, doublereal, 64);
-    MY_ALLOC(oldu, doublereal, 64);
-    MY_ALLOC(list2, integer, 128);
-    MY_ALLOC(point, integer, 90000);
-    MY_ALLOC(anchor, integer, 602);
-    MY_ALLOC(length, integer, 5760000);
-    MY_ALLOC(arrayi, integer, 64);
-    MY_ALLOC(levels, doublereal, 601);
-    MY_ALLOC(thirds, doublereal, 601);    
+
+    /* 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(arrayi, integer, (*n));
 
 /* +-----------------------------------------------------------------------+ */
 /* |    SUBROUTINE Direct                                                  | */
@@ -110,7 +126,7 @@ static integer c__3000 = 3000;
 /* |     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.                | */
@@ -131,10 +147,10 @@ static integer c__3000 = 3000;
 /* |  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                           | */
@@ -145,7 +161,7 @@ static integer c__3000 = 3000;
 /* |            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 :                                             | */
@@ -166,14 +182,14 @@ static integer c__3000 = 3000;
 /* |              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 :                                                    | */
@@ -331,12 +347,12 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
     cheat = 0;
     kmax = 1e10;
-    mdeep = 600;
+    mdeep = MAXDEEP;
 /* +-----------------------------------------------------------------------+ */
 /* | Write the header of the logfile.                                      | */
 /* +-----------------------------------------------------------------------+ */
     direct_dirheader_(logfile, &version, &x[1], n, eps, maxf, maxt, &l[1], &u[1], 
-           algmethod, &c__90000, &c__600, fglobal, fglper, ierror, &epsfix, &
+           algmethod, &MAXFUNC, &MAXDEEP, fglobal, fglper, ierror, &epsfix, &
                      iepschange, volper, sigmaper);
 /* +-----------------------------------------------------------------------+ */
 /* | If an error has occured while writing the header (we do some checking | */
@@ -364,7 +380,7 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
 /* | Initialiase the lists.                                                | */
 /* +-----------------------------------------------------------------------+ */
-    direct_dirinitlist_(anchor, &ifree, point, f, &c__90000, &c__600);
+    direct_dirinitlist_(anchor, &ifree, point, f, &MAXFUNC, &MAXDEEP);
 /* +-----------------------------------------------------------------------+ */
 /* | Call the routine to initialise the mapping of x from the n-dimensional| */
 /* | unit cube to the hypercube given by u and l. If an error occured,     | */
@@ -388,8 +404,9 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
     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, &c__90000, &c__600, n, &c__64, &
-           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.                                                 | */
 /* +-----------------------------------------------------------------------+ */
@@ -404,7 +421,9 @@ static integer c__3000 = 3000;
                 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;
@@ -412,7 +431,7 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
 /* | 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)
@@ -421,7 +440,7 @@ static integer c__3000 = 3000;
     } else {
         if (logfile)
              fprintf(logfile, "%d, %d, %g, %g\n", 
-                     tstart-1, numfunc, *fmin, fmax);
+                     tstart-1, numfunc, *minf, fmax);
     }
 /* +-----------------------------------------------------------------------+ */
 /* +-----------------------------------------------------------------------+ */
@@ -435,8 +454,8 @@ static integer c__3000 = 3000;
 /* | in the list S.                                                        | */
 /* +-----------------------------------------------------------------------+ */
        actdeep = actmaxdeep;
-       direct_dirchoose_(anchor, s, &c__600, f, fmin, *eps, epsabs, levels, &maxpos, length, 
-               &c__90000, &c__600, &c__3000, n, logfile, &cheat, &
+       direct_dirchoose_(anchor, s, &MAXDEEP, f, minf, *eps, epsabs, levels, &maxpos, length, 
+               &MAXFUNC, &MAXDEEP, &MAXDIV, n, logfile, &cheat, &
                kmax, &ifeasiblef, jones);
 /* +-----------------------------------------------------------------------+ */
 /* | Add other hyperrectangles to S, which have the same level and the same| */
@@ -445,8 +464,8 @@ static integer c__3000 = 3000;
 /* | JG 07/16/01 Added Errorflag.                                          | */
 /* +-----------------------------------------------------------------------+ */
        if (*algmethod == 0) {
-            direct_dirdoubleinsert_(anchor, s, &maxpos, point, f, &c__600, &c__90000,
-                    &c__3000, ierror);
+            direct_dirdoubleinsert_(anchor, s, &maxpos, point, f, &MAXDEEP, &MAXFUNC,
+                    &MAXDIV, ierror);
            if (*ierror == -6) {
                if (logfile)
                     fprintf(logfile,
@@ -464,7 +483,7 @@ static integer c__3000 = 3000;
        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.                      | */
 /* +-----------------------------------------------------------------------+ */
@@ -472,10 +491,10 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
 /* | JG 09/24/00 Calculate the value delta used for sampling points.       | */
 /* +-----------------------------------------------------------------------+ */
-               actdeep_div__ = direct_dirgetmaxdeep_(&s[j - 1], length, &c__90000
+               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.                                           | */
@@ -498,12 +517,12 @@ static integer c__3000 = 3000;
                    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_dirget_i__(length, &help, arrayi, &maxi, n, &c__90000);
+               direct_dirget_i__(length, &help, arrayi, &maxi, n, &MAXFUNC);
 /* +-----------------------------------------------------------------------+ */
 /* | Sample the function. To do this, we first calculate the points where  | */
 /* | we need to sample the function. After checking for errors, we then do | */
@@ -512,8 +531,8 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
                direct_dirsamplepoints_(c__, arrayi, &delta, &help, &start, length, 
                        logfile, f, &ifree, &maxi, point, &x[
-                       1], &l[1], fmin, &minpos, &u[1], n, &c__90000, &
-                       c__600, &oops);
+                       1], &l[1], minf, &minpos, &u[1], n, &MAXFUNC, &
+                       MAXDEEP, &oops);
                if (oops > 0) {
                    if (logfile)
                         fprintf(logfile, "WARNING: Error occured in routine DIRsamplepoints.\n");
@@ -526,9 +545,17 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
                direct_dirsamplef_(c__, arrayi, &delta, &help, &start, length,
                            logfile, f, &ifree, &maxi, point, fcn, &x[
-                       1], &l[1], fmin, &minpos, &u[1], n, &c__90000, &
-                       c__600, &oops, &fmax, &ifeasiblef, &iinfesiblef, 
-                       fcn_data);
+                       1], &l[1], minf, &minpos, &u[1], n, &MAXFUNC, &
+                       MAXDEEP, &oops, &fmax, &ifeasiblef, &iinfesiblef, 
+                                  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");
@@ -539,12 +566,12 @@ static integer c__3000 = 3000;
 /* | Divide the intervalls.                                                | */
 /* +-----------------------------------------------------------------------+ */
                direct_dirdivide_(&start, &actdeep_div__, length, point, arrayi, &
-                       help, list2, w, &maxi, f, &c__90000, &c__600, n);
+                       help, list2, w, &maxi, f, &MAXFUNC, &MAXDEEP, n);
 /* +-----------------------------------------------------------------------+ */
 /* | Insert the new intervalls into the list (sorted).                     | */
 /* +-----------------------------------------------------------------------+ */
                direct_dirinsertlist_(&start, anchor, point, f, &maxi, length, &
-                       c__90000, &c__600, n, &help, jones);
+                       MAXFUNC, &MAXDEEP, n, &help, jones);
 /* +-----------------------------------------------------------------------+ */
 /* | Increase the number of function evaluations.                          | */
 /* +-----------------------------------------------------------------------+ */
@@ -563,7 +590,7 @@ static integer c__3000 = 3000;
        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      | */
@@ -580,7 +607,7 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
 /* +-----------------------------------------------------------------------+ */
 /* | 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     | */
@@ -589,7 +616,7 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
        *ierror = jones;
        jones = 0;
-       actdeep_div__ = direct_dirgetlevel_(&minpos, length, &c__90000, n, jones);
+       actdeep_div__ = direct_dirgetlevel_(&minpos, length, &MAXFUNC, n, jones);
        jones = *ierror;
 /* +-----------------------------------------------------------------------+ */
 /* | JG 07/16/01 Use precalculated values to calculate volume.             | */
@@ -605,10 +632,10 @@ static integer c__3000 = 3000;
        }
 /* +-----------------------------------------------------------------------+ */
 /* | 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, &c__90000, n, jones);
+       actdeep_div__ = direct_dirgetlevel_(&minpos, length, &MAXFUNC, n, jones);
        delta = levels[actdeep_div__];
        if (delta <= *sigmaper) {
            *ierror = 5;
@@ -622,10 +649,10 @@ static integer c__3000 = 3000;
 /* | 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;
        }
 /* +-----------------------------------------------------------------------+ */
@@ -636,7 +663,7 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
        if (iinfesiblef > 0) {
             direct_dirreplaceinf_(&ifree, &ifreeold, f, c__, thirds, length, anchor, 
-                   point, &u[1], &l[1], &c__90000, &c__600, &c__64, n, 
+                   point, &u[1], &l[1], &MAXFUNC, &MAXDEEP, n, n, 
                    logfile, &fmax, jones);
        }
        ifreeold = ifree;
@@ -645,7 +672,7 @@ static integer c__3000 = 3000;
 /* +-----------------------------------------------------------------------+ */
        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.                                                    | */