direct/DIRect.c

   1 /* DIRect.f -- translated by f2c (version 20050501).
   2
   3    f2c output hand-cleaned by SGJ (August 2007).
   4 */
   5
   6 #include <math.h>
   7 #include "direct-internal.h"
   8
   9 /* Common Block Declarations */
  10
  11 /* Table of constant values */
  12
  13 static integer c__90000 = 90000;
  14 static integer c__600 = 600;
  15 static integer c__64 = 64;
  16 static integer c__3000 = 3000;
  17
  18 /* +-----------------------------------------------------------------------+ */
  19 /* | Program       : Direct.f                                              | */
  20 /* | Last modified : 07-16-2001                                            | */
  21 /* | Written by    : Joerg Gablonsky (jmgablon@unity.ncsu.edu)             | */
  22 /* |                 North Carolina State University                       | */
  23 /* |                 Dept. of Mathematics                                  | */
  24 /* | DIRECT is a method to solve problems of the form:                     | */
  25 /* |              min f: Q --> R,                                          | */
  26 /* | where f is the function to be minimized and Q is an n-dimensional     | */
  27 /* | hyperrectangle given by the the following equation:                   | */
  28 /* |                                                                       | */
  29 /* |       Q={ x : l(i) <= x(i) <= u(i), i = 1,...,n }.                    | */
  30 /* | Note: This version of DIRECT can also handle hidden constraints. By   | */
  31 /* |       this we mean that the function may not be defined over the whole| */
  32 /* |       hyperrectangle Q, but only over a subset, which is not given    | */
  33 /* |       analytically.                                                   | */
  34 /* |                                                                       | */
  35 /* | We now give a brief outline of the algorithm:                         | */
  36 /* |                                                                       | */
  37 /* |   The algorithm starts with mapping the hyperrectangle Q to the       | */
  38 /* |   n-dimensional unit hypercube. DIRECT then samples the function at   | */
  39 /* |   the center of this hypercube and at 2n more points, 2 in each       | */
  40 /* |   coordinate direction. Uisng these function values, DIRECT then      | */
  41 /* |   divides the domain into hyperrectangles, each having exactly one of | */
  42 /* |   the sampling points as its center. In each iteration, DIRECT chooses| */
  43 /* |   some of the existing hyperrectangles to be further divided.         | */
  44 /* |   We provide two different strategies of how to decide which          | */
  45 /* |   hyperrectangles DIRECT divides and several different convergence    | */
  46 /* |   criteria.                                                           | */
  47 /* |                                                                       | */
  48 /* |   DIRECT was designed to solve problems where the function f is       | */
  49 /* |   Lipschitz continues. However, DIRECT has proven to be effective on  | */
  50 /* |   more complex problems than these.                                   | */
  51 /* +-----------------------------------------------------------------------+ */
  52 /* Subroutine */ void direct_direct_(fp fcn, doublereal *x, integer *n, doublereal *eps, doublereal epsabs, integer *maxf, integer *maxt, doublereal *fmin, doublereal *l,
  53         doublereal *u, integer *algmethod, integer *ierror, FILE *logfile,
  54         doublereal *fglobal, doublereal *fglper, doublereal *volper,
  55         doublereal *sigmaper, void *fcn_data)
  56 {
  57     /* System generated locals */
  58     integer i__1, i__2;
  59     doublereal d__1;
  60
  61     /* Local variables */
  62     integer increase;
  63     doublereal *c__ = 0 /* was [90000][64] */, *f = 0   /*
  64             was [90000][2] */;
  65     integer i__, j, *s = 0      /* was [3000][2] */, t;
  66     doublereal *w = 0;
  67     doublereal divfactor;
  68     integer ifeasiblef, iepschange, actmaxdeep;
  69     integer actdeep_div__, iinfesiblef;
  70     integer pos1, newtosample;
  71     integer ifree, help;
  72     doublereal *oldl = 0, fmax;
  73     integer maxi;
  74     doublereal kmax, *oldu = 0;
  75     integer oops, *list2 = 0    /* was [64][2] */, cheat;
  76     doublereal delta;
  77     integer mdeep, *point = 0, start;
  78     integer *anchor = 0, *length = 0    /* was [90000][64] */, *arrayi = 0;
  79     doublereal *levels = 0, *thirds = 0;
  80     integer writed;
  81     doublereal epsfix;
  82     integer oldpos, minpos, maxpos, tstart, actdeep, ifreeold, oldmaxf;
  83     integer numfunc, version;
  84     integer jones;
  85
  86     /* FIXME: change sizes dynamically? */
  87 #define MY_ALLOC(p, t, n) p = (t *) malloc(sizeof(t) * (n)); \
  88                           if (!(p)) { *ierror = -100; goto cleanup; }
  89     MY_ALLOC(c__, doublereal, 5760000);
  90     MY_ALLOC(f, doublereal, 180000);
  91     MY_ALLOC(s, integer, 6000);
  92     MY_ALLOC(w, doublereal, 64);
  93     MY_ALLOC(oldl, doublereal, 64);
  94     MY_ALLOC(oldu, doublereal, 64);
  95     MY_ALLOC(list2, integer, 128);
  96     MY_ALLOC(point, integer, 90000);
  97     MY_ALLOC(anchor, integer, 602);
  98     MY_ALLOC(length, integer, 5760000);
  99     MY_ALLOC(arrayi, integer, 64);
 100     MY_ALLOC(levels, doublereal, 601);
 101     MY_ALLOC(thirds, doublereal, 601);
 102
 103 /* +-----------------------------------------------------------------------+ */
 104 /* |    SUBROUTINE Direct                                                  | */
 105 /* | On entry                                                              | */
 106 /* |     fcn -- The argument containing the name of the user-supplied      | */
 107 /* |            SUBROUTINE that returns values for the function to be      | */
 108 /* |            minimized.                                                 | */
 109 /* |       n -- The dimension of the problem.                              | */
 110 /* |     eps -- Exceeding value. If eps > 0, we use the same epsilon for   | */
 111 /* |            all iterations. If eps < 0, we use the update formula from | */
 112 /* |            Jones:                                                     | */
 113 /* |               eps = max(1.D-4*abs(fmin),epsfix),                      | */
 114 /* |            where epsfix = abs(eps), the absolute value of eps which is| */
 115 /* |            passed to the function.                                    | */
 116 /* |    maxf -- The maximum number of function evaluations.                | */
 117 /* |    maxT -- The maximum number of iterations.                          | */
 118 /* |            Direct stops when either the maximum number of iterations  | */
 119 /* |            is reached or more than maxf function-evalutions were made.| */
 120 /* |       l -- The lower bounds of the hyperbox.                          | */
 121 /* |       u -- The upper bounds of the hyperbox.                          | */
 122 /* |algmethod-- Choose the method, that is either use the original method  | */
 123 /* |            as described by Jones et.al. (0) or use our modification(1)| */
 124 /* | logfile -- File-Handle for the logfile. DIRECT expects this file to be| */
 125 /* |            opened and closed by the user outside of DIRECT. We moved  | */
 126 /* |            this to the outside so the user can add extra informations | */
 127 /* |            to this file before and after the call to DIRECT.          | */
 128 /* | fglobal -- Function value of the global optimum. If this value is not | */
 129 /* |            known (that is, we solve a real problem, not a testproblem)| */
 130 /* |            set this value to -1.D100 and fglper (see below) to 0.D0.  | */
 131 /* |  fglper -- Terminate the optimization when the percent error          | */
 132 /* |                100(f_min - fglobal)/max(1,abs(fglobal)) < fglper.     | */
 133 /* |  volper -- Terminate the optimization when the volume of the          | */
 134 /* |            hyperrectangle S with f(c(S)) = fmin is less then volper   | */
 135 /* |            percent of the volume of the original hyperrectangle.      | */
 136 /* |sigmaper -- Terminate the optimization when the measure of the         | */
 137 /* |            hyperrectangle S with f(c(S)) = fmin is less then sigmaper.| */
 138 /* |                                                                       | */
 139 /* | User data that is passed through without being changed:               | */
 140 /* |  fcn_data - opaque pointer to any user data                           | */
 141 /* |                                                                       | */
 142 /* | On return                                                             | */
 143 /* |                                                                       | */
 144 /* |       x -- The final point obtained in the optimization process.      | */
 145 /* |            X should be a good approximation to the global minimum     | */
 146 /* |            for the function within the hyper-box.                     | */
 147 /* |                                                                       | */
 148 /* |    fmin -- The value of the function at x.                            | */
 149 /* |  Ierror -- Error flag. If Ierror is lower 0, an error has occured. The| */
 150 /* |            values of Ierror mean                                      | */
 151 /* |            Fatal errors :                                             | */
 152 /* |             -1   u(i) <= l(i) for some i.                             | */
 153 /* |             -2   maxf is too large.                                   | */
 154 /* |             -3   Initialization in DIRpreprc failed.                  | */
 155 /* |             -4   Error in DIRSamplepoints, that is there was an error | */
 156 /* |                  in the creation of the sample points.                | */
 157 /* |             -5   Error in DIRSamplef, that is an error occured while  | */
 158 /* |                  the function was sampled.                            | */
 159 /* |             -6   Error in DIRDoubleInsert, that is an error occured   | */
 160 /* |                  DIRECT tried to add all hyperrectangles with the same| */
 161 /* |                  size and function value at the center. Either        | */
 162 /* |                  increase maxdiv or use our modification (Jones = 1). | */
 163 /* |            Termination values :                                       | */
 164 /* |              1   Number of function evaluations done is larger then   | */
 165 /* |                  maxf.                                                | */
 166 /* |              2   Number of iterations is equal to maxT.               | */
 167 /* |              3   The best function value found is within fglper of    | */
 168 /* |                  the (known) global optimum, that is                  | */
 169 /* |                   100(fmin - fglobal/max(1,|fglobal|))  < fglper.     | */
 170 /* |                  Note that this termination signal only occurs when   | */
 171 /* |                  the global optimal value is known, that is, a test   | */
 172 /* |                  function is optimized.                               | */
 173 /* |              4   The volume of the hyperrectangle with fmin at its    | */
 174 /* |                  center is less than volper percent of the volume of  | */
 175 /* |                  the original hyperrectangle.                         | */
 176 /* |              5   The measure of the hyperrectangle with fmin at its   | */
 177 /* |                  center is less than sigmaper.                        | */
 178 /* |                                                                       | */
 179 /* | SUBROUTINEs used :                                                    | */
 180 /* |                                                                       | */
 181 /* | DIRheader, DIRInitSpecific, DIRInitList, DIRpreprc, DIRInit, DIRChoose| */
 182 /* | DIRDoubleInsert, DIRGet_I, DIRSamplepoints, DIRSamplef, DIRDivide     | */
 183 /* | DIRInsertList, DIRreplaceInf, DIRWritehistbox, DIRsummary, Findareas  | */
 184 /* |                                                                       | */
 185 /* | Functions used :                                                      | */
 186 /* |                                                                       | */
 187 /* | DIRgetMaxdeep, DIRgetlevel                                            | */
 188 /* +-----------------------------------------------------------------------+ */
 189 /* +-----------------------------------------------------------------------+ */
 190 /* | Parameters                                                            | */
 191 /* +-----------------------------------------------------------------------+ */
 192 /* +-----------------------------------------------------------------------+ */
 193 /* | The maximum of function evaluations allowed.                          | */
 194 /* | The maximum dept of the algorithm.                                    | */
 195 /* | The maximum number of divisions allowed.                              | */
 196 /* | The maximal dimension of the problem.                                 | */
 197 /* +-----------------------------------------------------------------------+ */
 198 /* +-----------------------------------------------------------------------+ */
 199 /* | Global Variables.                                                     | */
 200 /* +-----------------------------------------------------------------------+ */
 201 /* +-----------------------------------------------------------------------+ */
 202 /* | EXTERNAL Variables.                                                   | */
 203 /* +-----------------------------------------------------------------------+ */
 204 /* +-----------------------------------------------------------------------+ */
 205 /* | User Variables.                                                       | */
 206 /* | These can be used to pass user defined data to the function to be     | */
 207 /* | optimized.                                                            | */
 208 /* +-----------------------------------------------------------------------+ */
 209 /* +-----------------------------------------------------------------------+ */
 210 /* | Place to define, if needed, some application-specific variables.      | */
 211 /* | Note: You should try to use the arrays defined above for this.        | */
 212 /* +-----------------------------------------------------------------------+ */
 213 /* +-----------------------------------------------------------------------+ */
 214 /* | End of application - specific variables !                             | */
 215 /* +-----------------------------------------------------------------------+ */
 216 /* +-----------------------------------------------------------------------+ */
 217 /* | Internal variables :                                                  | */
 218 /* |       f -- values of functions.                                       | */
 219 /* |divfactor-- Factor used for termination with known global minimum.     | */
 220 /* |  anchor -- anchors of lists with deepness i, -1 is anchor for list of | */
 221 /* |            NaN - values.                                              | */
 222 /* |       S -- List of potentially optimal points.                        | */
 223 /* |   point -- lists                                                      | */
 224 /* |    ifree -- first free position                                        | */
 225 /* |       c -- midpoints of arrays                                        | */
 226 /* |  thirds -- Precalculated values of 1/3^i.                             | */
 227 /* |  levels -- Length of intervals.                                       | */
 228 /* |  length -- Length of intervall (index)                                | */
 229 /* |       t -- actual iteration                                           | */
 230 /* |       j -- loop-variable                                              | */
 231 /* | actdeep -- the actual minimal interval-length index                   | */
 232 /* |  Minpos -- position of the actual minimum                             | */
 233 /* |    file -- The filehandle for a datafile.                             | */
 234 /* |  maxpos -- The number of intervalls, which are truncated.             | */
 235 /* |    help -- A help variable.                                           | */
 236 /* | numfunc -- The actual number of function evaluations.                 | */
 237 /* |   file2 -- The filehandle for an other datafile.                      | */
 238 /* |  ArrayI -- Array with the indexes of the sides with maximum length.   | */
 239 /* |    maxi -- Number of directions with maximal side length.             | */
 240 /* |    oops -- Flag which shows if anything went wrong in the             | */
 241 /* |            initialisation.                                            | */
 242 /* |   cheat -- Obsolete. If equal 1, we don't allow Ktilde > kmax.        | */
 243 /* |  writed -- If writed=1, store final division to plot with Matlab.     | */
 244 /* |   List2 -- List of indicies of intervalls, which are to be truncated. | */
 245 /* |       i -- Another loop-variable.                                     | */
 246 /* |actmaxdeep-- The actual maximum (minimum) of possible Interval length. | */
 247 /* |  oldpos -- The old index of the minimum. Used to print only, if there | */
 248 /* |            is a new minimum found.                                    | */
 249 /* |  tstart -- The start of the outer loop.                               | */
 250 /* |   start -- The postion of the starting point in the inner loop.       | */
 251 /* | Newtosample -- The total number of points to sample in the inner loop.| */
 252 /* |       w -- Array used to divide the intervalls                        | */
 253 /* |    kmax -- Obsolete. If cheat = 1, Ktilde was not allowed to be larger| */
 254 /* |            than kmax. If Ktilde > kmax, we set ktilde = kmax.         | */
 255 /* |   delta -- The distance to new points from center of old hyperrec.    | */
 256 /* |    pos1 -- Help variable used as an index.                            | */
 257 /* | version -- Store the version number of DIRECT.                        | */
 258 /* | oldmaxf -- Store the original function budget.                        | */
 259 /* |increase -- Flag used to keep track if function budget was increased   | */
 260 /* |            because no feasible point was found.                       | */
 261 /* | ifreeold -- Keep track which index was free before. Used with          | */
 262 /* |            SUBROUTINE DIRReplaceInf.                                  | */
 263 /* |actdeep_div-- Keep track of the current depths for divisions.          | */
 264 /* |    oldl -- Array used to store the original bounds of the domain.     | */
 265 /* |    oldu -- Array used to store the original bounds of the domain.     | */
 266 /* |  epsfix -- If eps < 0, we use Jones update formula. epsfix stores the | */
 267 /* |            absolute value of epsilon.                                 | */
 268 /* |iepschange-- flag iepschange to store if epsilon stays fixed or is     | */
 269 /* |             changed.                                                  | */
 270 /* |DIRgetMaxdeep-- Function to calculate the level of a hyperrectangle.   | */
 271 /* |DIRgetlevel-- Function to calculate the level and stage of a hyperrec. | */
 272 /* |    fmax -- Keep track of the maximum value of the function found.     | */
 273 /* |Ifeasiblef-- Keep track if a feasible point has  been found so far.    | */
 274 /* |             Ifeasiblef = 0 means a feasible point has been found,     | */
 275 /* |             Ifeasiblef = 1 no feasible point has been found.          | */
 276 /* +-----------------------------------------------------------------------+ */
 277 /* +-----------------------------------------------------------------------+ */
 278 /* | JG 09/25/00 Version counter.                                          | */
 279 /* +-----------------------------------------------------------------------+ */
 280 /* +-----------------------------------------------------------------------+ */
 281 /* | JG 09/24/00 Add another actdeep to keep track of the current depths   | */
 282 /* |             for divisions.                                            | */
 283 /* +-----------------------------------------------------------------------+ */
 284 /* +-----------------------------------------------------------------------+ */
 285 /* |JG 01/13/01 Added epsfix for epsilon update. If eps < 0, we use Jones  | */
 286 /* |            update formula. epsfix stores the absolute value of epsilon| */
 287 /* |            then. Also added flag iepschange to store if epsilon stays | */
 288 /* |            fixed or is changed.                                       | */
 289 /* +-----------------------------------------------------------------------+ */
 290 /* +-----------------------------------------------------------------------+ */
 291 /* | JG 01/22/01 fmax is used to keep track of the maximum value found.    | */
 292 /* +-----------------------------------------------------------------------+ */
 293 /* +-----------------------------------------------------------------------+ */
 294 /* | JG 01/22/01 Ifeasiblef is used to keep track if a feasible point has  | */
 295 /* |             been found so far. Ifeasiblef = 0 means a feasible point  | */
 296 /* |             has been found, Ifeasiblef = 1 if not.                    | */
 297 /* | JG 03/09/01 IInfeasible is used to keep track if an infeasible point  | */
 298 /* |             has been found. IInfeasible > 0 means a infeasible point  | */
 299 /* |             has been found, IInfeasible = 0 if not.                   | */
 300 /* +-----------------------------------------------------------------------+ */
 301 /* +-----------------------------------------------------------------------+ */
 302 /* +-----------------------------------------------------------------------+ */
 303 /* |                            Start of code.                             | */
 304 /* +-----------------------------------------------------------------------+ */
 305 /* +-----------------------------------------------------------------------+ */
 306     /* Parameter adjustments */
 307     --u;
 308     --l;
 309     --x;
 310
 311     /* Function Body */
 312     writed = 0;
 313     jones = *algmethod;
 314 /* +-----------------------------------------------------------------------+ */
 315 /* | Save the upper and lower bounds.                                      | */
 316 /* +-----------------------------------------------------------------------+ */
 317     i__1 = *n;
 318     for (i__ = 1; i__ <= i__1; ++i__) {
 319         oldu[i__ - 1] = u[i__];
 320         oldl[i__ - 1] = l[i__];
 321 /* L150: */
 322     }
 323 /* +-----------------------------------------------------------------------+ */
 324 /* | Set version.                                                          | */
 325 /* +-----------------------------------------------------------------------+ */
 326     version = 204;
 327 /* +-----------------------------------------------------------------------+ */
 328 /* | Set parameters.                                                       | */
 329 /* |    If cheat > 0, we do not allow \tilde{K} to be larger than kmax, and| */
 330 /* |    set \tilde{K} to set value if necessary. Not used anymore.         | */
 331 /* +-----------------------------------------------------------------------+ */
 332     cheat = 0;
 333     kmax = 1e10;
 334     mdeep = 600;
 335 /* +-----------------------------------------------------------------------+ */
 336 /* | Write the header of the logfile.                                      | */
 337 /* +-----------------------------------------------------------------------+ */
 338     direct_dirheader_(logfile, &version, &x[1], n, eps, maxf, maxt, &l[1], &u[1],
 339             algmethod, &c__90000, &c__600, fglobal, fglper, ierror, &epsfix, &
 340                       iepschange, volper, sigmaper);
 341 /* +-----------------------------------------------------------------------+ */
 342 /* | If an error has occured while writing the header (we do some checking | */
 343 /* | of variables there), return to the main program.                      | */
 344 /* +-----------------------------------------------------------------------+ */
 345     if (*ierror < 0) {
 346         goto cleanup;
 347     }
 348 /* +-----------------------------------------------------------------------+ */
 349 /* | If the known global minimum is equal 0, we cannot divide by it.       | */
 350 /* | Therefore we set it to 1. If not, we set the divisionfactor to the    | */
 351 /* | absolute value of the global minimum.                                 | */
 352 /* +-----------------------------------------------------------------------+ */
 353     if (*fglobal == 0.) {
 354         divfactor = 1.;
 355     } else {
 356         divfactor = fabs(*fglobal);
 357     }
 358 /* +-----------------------------------------------------------------------+ */
 359 /* | Save the budget given by the user. The variable maxf will be changed  | */
 360 /* | if in the beginning no feasible points are found.                     | */
 361 /* +-----------------------------------------------------------------------+ */
 362     oldmaxf = *maxf;
 363     increase = 0;
 364 /* +-----------------------------------------------------------------------+ */
 365 /* | Initialiase the lists.                                                | */
 366 /* +-----------------------------------------------------------------------+ */
 367     direct_dirinitlist_(anchor, &ifree, point, f, &c__90000, &c__600);
 368 /* +-----------------------------------------------------------------------+ */
 369 /* | Call the routine to initialise the mapping of x from the n-dimensional| */
 370 /* | unit cube to the hypercube given by u and l. If an error occured,     | */
 371 /* | give out a error message and return to the main program with the error| */
 372 /* | flag set.                                                             | */
 373 /* | JG 07/16/01 Changed call to remove unused data.                       | */
 374 /* +-----------------------------------------------------------------------+ */
 375     direct_dirpreprc_(&u[1], &l[1], n, &l[1], &u[1], &oops);
 376     if (oops > 0) {
 377         if (logfile)
 378              fprintf(logfile,"WARNING: Initialization in DIRpreprc failed.\n");
 379         *ierror = -3;
 380         goto cleanup;
 381     }
 382     tstart = 2;
 383 /* +-----------------------------------------------------------------------+ */
 384 /* | Initialise the algorithm DIRECT.                                      | */
 385 /* +-----------------------------------------------------------------------+ */
 386 /* +-----------------------------------------------------------------------+ */
 387 /* | Added variable to keep track of the maximum value found.              | */
 388 /* +-----------------------------------------------------------------------+ */
 389     direct_dirinit_(f, fcn, c__, length, &actdeep, point, anchor, &ifree,
 390             logfile, arrayi, &maxi, list2, w, &x[1], &l[1], &u[1],
 391             fmin, &minpos, thirds, levels, &c__90000, &c__600, n, &c__64, &
 392             fmax, &ifeasiblef, &iinfesiblef, ierror, fcn_data, jones);
 393 /* +-----------------------------------------------------------------------+ */
 394 /* | Added error checking.                                                 | */
 395 /* +-----------------------------------------------------------------------+ */
 396     if (*ierror < 0) {
 397         if (*ierror == -4) {
 398             if (logfile)
 399                  fprintf(logfile, "WARNING: Error occured in routine DIRsamplepoints.\n");
 400             goto cleanup;
 401         }
 402         if (*ierror == -5) {
 403             if (logfile)
 404                  fprintf(logfile, "WARNING: Error occured in routine DIRsamplef..\n");
 405             goto cleanup;
 406         }
 407     }
 408     numfunc = maxi + 1 + maxi;
 409     actmaxdeep = 1;
 410     oldpos = 0;
 411     tstart = 2;
 412 /* +-----------------------------------------------------------------------+ */
 413 /* | If no feasible point has been found, give out the iteration, the      | */
 414 /* | number of function evaluations and a warning. Otherwise, give out     | */
 415 /* | the iteration, the number of function evaluations done and fmin.      | */
 416 /* +-----------------------------------------------------------------------+ */
 417     if (ifeasiblef > 0) {
 418          if (logfile)
 419               fprintf(logfile, "No feasible point found in %d iterations "
 420                       "and %d function evaluations.\n", tstart-1, numfunc);
 421     } else {
 422          if (logfile)
 423               fprintf(logfile, "%d, %d, %g, %g\n",
 424                       tstart-1, numfunc, *fmin, fmax);
 425     }
 426 /* +-----------------------------------------------------------------------+ */
 427 /* +-----------------------------------------------------------------------+ */
 428 /* | Main loop!                                                            | */
 429 /* +-----------------------------------------------------------------------+ */
 430 /* +-----------------------------------------------------------------------+ */
 431     i__1 = *maxt;
 432     for (t = tstart; t <= i__1; ++t) {
 433 /* +-----------------------------------------------------------------------+ */
 434 /* | Choose the sample points. The indices of the sample points are stored | */
 435 /* | in the list S.                                                        | */
 436 /* +-----------------------------------------------------------------------+ */
 437         actdeep = actmaxdeep;
 438         direct_dirchoose_(anchor, s, &c__600, f, fmin, *eps, epsabs, levels, &maxpos, length,
 439                 &c__90000, &c__600, &c__3000, n, logfile, &cheat, &
 440                 kmax, &ifeasiblef, jones);
 441 /* +-----------------------------------------------------------------------+ */
 442 /* | Add other hyperrectangles to S, which have the same level and the same| */
 443 /* | function value at the center as the ones found above (that are stored | */
 444 /* | in S). This is only done if we use the original DIRECT algorithm.     | */
 445 /* | JG 07/16/01 Added Errorflag.                                          | */
 446 /* +-----------------------------------------------------------------------+ */
 447         if (*algmethod == 0) {
 448              direct_dirdoubleinsert_(anchor, s, &maxpos, point, f, &c__600, &c__90000,
 449                      &c__3000, ierror);
 450             if (*ierror == -6) {
 451                 if (logfile)
 452                      fprintf(logfile,
 453 "WARNING: Capacity of array S in DIRDoubleInsert reached. Increase maxdiv.\n"
 454 "This means that there are a lot of hyperrectangles with the same function\n"
 455 "value at the center. We suggest to use our modification instead (Jones = 1)\n"
 456                           );
 457                 goto cleanup;
 458             }
 459         }
 460         oldpos = minpos;
 461 /* +-----------------------------------------------------------------------+ */
 462 /* | Initialise the number of sample points in this outer loop.            | */
 463 /* +-----------------------------------------------------------------------+ */
 464         newtosample = 0;
 465         i__2 = maxpos;
 466         for (j = 1; j <= i__2; ++j) {
 467             actdeep = s[j + 2999];
 468 /* +-----------------------------------------------------------------------+ */
 469 /* | If the actual index is a point to sample, do it.                      | */
 470 /* +-----------------------------------------------------------------------+ */
 471             if (s[j - 1] > 0) {
 472 /* +-----------------------------------------------------------------------+ */
 473 /* | JG 09/24/00 Calculate the value delta used for sampling points.       | */
 474 /* +-----------------------------------------------------------------------+ */
 475                 actdeep_div__ = direct_dirgetmaxdeep_(&s[j - 1], length, &c__90000,
 476                         n);
 477                 delta = thirds[actdeep_div__ + 1];
 478                 actdeep = s[j + 2999];
 479 /* +-----------------------------------------------------------------------+ */
 480 /* | If the current dept of division is only one under the maximal allowed | */
 481 /* | dept, stop the computation.                                           | */
 482 /* +-----------------------------------------------------------------------+ */
 483                 if (actdeep + 1 >= mdeep) {
 484                     if (logfile)
 485                          fprintf(logfile, "WARNING: Maximum number of levels reached. Increase maxdeep.\n");
 486                     *ierror = -6;
 487                     goto L100;
 488                 }
 489                 actmaxdeep = MAX(actdeep,actmaxdeep);
 490                 help = s[j - 1];
 491                 if (! (anchor[actdeep + 1] == help)) {
 492                     pos1 = anchor[actdeep + 1];
 493                     while(! (point[pos1 - 1] == help)) {
 494                         pos1 = point[pos1 - 1];
 495                     }
 496                     point[pos1 - 1] = point[help - 1];
 497                 } else {
 498                     anchor[actdeep + 1] = point[help - 1];
 499                 }
 500                 if (actdeep < 0) {
 501                     actdeep = (integer) f[help - 1];
 502                 }
 503 /* +-----------------------------------------------------------------------+ */
 504 /* | Get the Directions in which to decrease the intervall-length.         | */
 505 /* +-----------------------------------------------------------------------+ */
 506                 direct_dirget_i__(length, &help, arrayi, &maxi, n, &c__90000);
 507 /* +-----------------------------------------------------------------------+ */
 508 /* | Sample the function. To do this, we first calculate the points where  | */
 509 /* | we need to sample the function. After checking for errors, we then do | */
 510 /* | the actual evaluation of the function, again followed by checking for | */
 511 /* | errors.                                                               | */
 512 /* +-----------------------------------------------------------------------+ */
 513                 direct_dirsamplepoints_(c__, arrayi, &delta, &help, &start, length,
 514                         logfile, f, &ifree, &maxi, point, &x[
 515                         1], &l[1], fmin, &minpos, &u[1], n, &c__90000, &
 516                         c__600, &oops);
 517                 if (oops > 0) {
 518                     if (logfile)
 519                          fprintf(logfile, "WARNING: Error occured in routine DIRsamplepoints.\n");
 520                     *ierror = -4;
 521                     goto cleanup;
 522                 }
 523                 newtosample += maxi;
 524 /* +-----------------------------------------------------------------------+ */
 525 /* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
 526 /* +-----------------------------------------------------------------------+ */
 527                 direct_dirsamplef_(c__, arrayi, &delta, &help, &start, length,
 528                             logfile, f, &ifree, &maxi, point, fcn, &x[
 529                         1], &l[1], fmin, &minpos, &u[1], n, &c__90000, &
 530                         c__600, &oops, &fmax, &ifeasiblef, &iinfesiblef,
 531                         fcn_data);
 532                 if (oops > 0) {
 533                     if (logfile)
 534                          fprintf(logfile, "WARNING: Error occured in routine DIRsamplef.\n");
 535                     *ierror = -5;
 536                     goto cleanup;
 537                 }
 538 /* +-----------------------------------------------------------------------+ */
 539 /* | Divide the intervalls.                                                | */
 540 /* +-----------------------------------------------------------------------+ */
 541                 direct_dirdivide_(&start, &actdeep_div__, length, point, arrayi, &
 542                         help, list2, w, &maxi, f, &c__90000, &c__600, n);
 543 /* +-----------------------------------------------------------------------+ */
 544 /* | Insert the new intervalls into the list (sorted).                     | */
 545 /* +-----------------------------------------------------------------------+ */
 546                 direct_dirinsertlist_(&start, anchor, point, f, &maxi, length, &
 547                         c__90000, &c__600, n, &help, jones);
 548 /* +-----------------------------------------------------------------------+ */
 549 /* | Increase the number of function evaluations.                          | */
 550 /* +-----------------------------------------------------------------------+ */
 551                 numfunc = numfunc + maxi + maxi;
 552             }
 553 /* +-----------------------------------------------------------------------+ */
 554 /* | End of main loop.                                                     | */
 555 /* +-----------------------------------------------------------------------+ */
 556 /* L20: */
 557         }
 558 /* +-----------------------------------------------------------------------+ */
 559 /* | If there is a new minimum, show the actual iteration, the number of   | */
 560 /* | function evaluations, the minimum value of f (so far) and the position| */
 561 /* | in the array.                                                         | */
 562 /* +-----------------------------------------------------------------------+ */
 563         if (oldpos < minpos) {
 564             if (logfile)
 565                  fprintf(logfile, "%d, %d, %g, %g\n",
 566                          t, numfunc, *fmin, fmax);
 567         }
 568 /* +-----------------------------------------------------------------------+ */
 569 /* | If no feasible point has been found, give out the iteration, the      | */
 570 /* | number of function evaluations and a warning.                         | */
 571 /* +-----------------------------------------------------------------------+ */
 572         if (ifeasiblef > 0) {
 573             if (logfile)
 574                  fprintf(logfile, "No feasible point found in %d iterations "
 575                          "and %d function evaluations\n", t, numfunc);
 576         }
 577 /* +-----------------------------------------------------------------------+ */
 578 /* +-----------------------------------------------------------------------+ */
 579 /* |                       Termination Checks                              | */
 580 /* +-----------------------------------------------------------------------+ */
 581 /* +-----------------------------------------------------------------------+ */
 582 /* | JG 01/22/01 Calculate the index for the hyperrectangle at which       | */
 583 /* |             fmin is assumed. We then calculate the volume of this     | */
 584 /* |             hyperrectangle and store it in delta. This delta can be   | */
 585 /* |             used to stop DIRECT once the volume is below a certain    | */
 586 /* |             percentage of the original volume. Since the original     | */
 587 /* |             is 1 (scaled), we can stop once delta is below a certain  | */
 588 /* |             percentage, given by volper.                              | */
 589 /* +-----------------------------------------------------------------------+ */
 590         *ierror = jones;
 591         jones = 0;
 592         actdeep_div__ = direct_dirgetlevel_(&minpos, length, &c__90000, n, jones);
 593         jones = *ierror;
 594 /* +-----------------------------------------------------------------------+ */
 595 /* | JG 07/16/01 Use precalculated values to calculate volume.             | */
 596 /* +-----------------------------------------------------------------------+ */
 597         delta = thirds[actdeep_div__] * 100;
 598         if (delta <= *volper) {
 599             *ierror = 4;
 600             if (logfile)
 601                  fprintf(logfile, "DIRECT stopped: Volume of S_min is "
 602                          "%g%% < %g%% of the original volume.\n",
 603                          delta, *volper);
 604             goto L100;
 605         }
 606 /* +-----------------------------------------------------------------------+ */
 607 /* | JG 01/23/01 Calculate the measure for the hyperrectangle at which     | */
 608 /* |             fmin is assumed. If this measure is smaller then sigmaper,| */
 609 /* |             we stop DIRECT.                                           | */
 610 /* +-----------------------------------------------------------------------+ */
 611         actdeep_div__ = direct_dirgetlevel_(&minpos, length, &c__90000, n, jones);
 612         delta = levels[actdeep_div__];
 613         if (delta <= *sigmaper) {
 614             *ierror = 5;
 615             if (logfile)
 616                  fprintf(logfile, "DIRECT stopped: Measure of S_min "
 617                          "= %g < %g.\n", delta, *sigmaper);
 618             goto L100;
 619         }
 620 /* +-----------------------------------------------------------------------+ */
 621 /* | If the best found function value is within fglper of the (known)      | */
 622 /* | global minimum value, terminate. This only makes sense if this optimal| */
 623 /* | value is known, that is, in test problems.                            | */
 624 /* +-----------------------------------------------------------------------+ */
 625         if ((*fmin - *fglobal) * 100 / divfactor <= *fglper) {
 626             *ierror = 3;
 627             if (logfile)
 628                  fprintf(logfile, "DIRECT stopped: fmin within fglper of global minimum.\n");
 629             goto L100;
 630         }
 631 /* +-----------------------------------------------------------------------+ */
 632 /* | Find out if there are infeasible points which are near feasible ones. | */
 633 /* | If this is the case, replace the function value at the center of the  | */
 634 /* | hyper rectangle by the lowest function value of a nearby function.    | */
 635 /* | If no infeasible points exist (IInfesiblef = 0), skip this.           | */
 636 /* +-----------------------------------------------------------------------+ */
 637         if (iinfesiblef > 0) {
 638              direct_dirreplaceinf_(&ifree, &ifreeold, f, c__, thirds, length, anchor,
 639                     point, &u[1], &l[1], &c__90000, &c__600, &c__64, n,
 640                     logfile, &fmax, jones);
 641         }
 642         ifreeold = ifree;
 643 /* +-----------------------------------------------------------------------+ */
 644 /* | If iepschange = 1, we use the epsilon change formula from Jones.      | */
 645 /* +-----------------------------------------------------------------------+ */
 646         if (iepschange == 1) {
 647 /* Computing MAX */
 648             d__1 = fabs(*fmin) * 1e-4;
 649             *eps = MAX(d__1,epsfix);
 650         }
 651 /* +-----------------------------------------------------------------------+ */
 652 /* | If no feasible point has been found yet, set the maximum number of    | */
 653 /* | function evaluations to the number of evaluations already done plus   | */
 654 /* | the budget given by the user.                                         | */
 655 /* | If the budget has already be increased, increase it again. If a       | */
 656 /* | feasible point has been found, remark that and reset flag. No further | */
 657 /* | increase is needed.                                                   | */
 658 /* +-----------------------------------------------------------------------+ */
 659         if (increase == 1) {
 660             *maxf = numfunc + oldmaxf;
 661             if (ifeasiblef == 0) {
 662                 if (logfile)
 663                      fprintf(logfile, "DIRECT found a feasible point.  The "
 664                              "adjusted budget is now set to %d.\n", *maxf);
 665                 increase = 0;
 666             }
 667         }
 668 /* +-----------------------------------------------------------------------+ */
 669 /* | Check if the number of function evaluations done is larger than the   | */
 670 /* | allocated budget. If this is the case, check if a feasible point was  | */
 671 /* | found. If this is a case, terminate. If no feasible point was found,  | */
 672 /* | increase the budget and set flag increase.                            | */
 673 /* +-----------------------------------------------------------------------+ */
 674         if (numfunc > *maxf) {
 675             if (ifeasiblef == 0) {
 676                 *ierror = 1;
 677                 if (logfile)
 678                      fprintf(logfile, "DIRECT stopped: numfunc >= maxf.\n");
 679                 goto L100;
 680             } else {
 681                 increase = 1;
 682                 if (logfile)
 683                      fprintf(logfile,
 684 "DIRECT could not find a feasible point after %d function evaluations.\n"
 685 "DIRECT continues until a feasible point is found.\n", numfunc);
 686                 *maxf = numfunc + oldmaxf;
 687             }
 688         }
 689 /* L10: */
 690     }
 691 /* +-----------------------------------------------------------------------+ */
 692 /* +-----------------------------------------------------------------------+ */
 693 /* | End of main loop.                                                     | */
 694 /* +-----------------------------------------------------------------------+ */
 695 /* +-----------------------------------------------------------------------+ */
 696 /* +-----------------------------------------------------------------------+ */
 697 /* | The algorithm stopped after maxT iterations.                          | */
 698 /* +-----------------------------------------------------------------------+ */
 699     *ierror = 2;
 700     if (logfile)
 701          fprintf(logfile, "DIRECT stopped: maxT iterations.\n");
 702
 703 L100:
 704 /* +-----------------------------------------------------------------------+ */
 705 /* | Store the position of the minimum in x.                               | */
 706 /* +-----------------------------------------------------------------------+ */
 707     i__1 = *n;
 708     for (i__ = 1; i__ <= i__1; ++i__) {
 709         x[i__] = c__[minpos + i__ * 90000 - 90001] * l[i__] + l[i__] * u[i__];
 710         u[i__] = oldu[i__ - 1];
 711         l[i__] = oldl[i__ - 1];
 712 /* L50: */
 713     }
 714 /* +-----------------------------------------------------------------------+ */
 715 /* | Store the number of function evaluations in maxf.                     | */
 716 /* +-----------------------------------------------------------------------+ */
 717     *maxf = numfunc;
 718 /* +-----------------------------------------------------------------------+ */
 719 /* | Give out a summary of the run.                                        | */
 720 /* +-----------------------------------------------------------------------+ */
 721     direct_dirsummary_(logfile, &x[1], &l[1], &u[1], n, fmin, fglobal, &numfunc,
 722             ierror);
 723 /* +-----------------------------------------------------------------------+ */
 724 /* | Format statements.                                                    | */
 725 /* +-----------------------------------------------------------------------+ */
 726
 727  cleanup:
 728 #define MY_FREE(p) if (p) free(p)
 729     MY_FREE(c__);
 730     MY_FREE(f);
 731     MY_FREE(s);
 732     MY_FREE(w);
 733     MY_FREE(oldl);
 734     MY_FREE(oldu);
 735     MY_FREE(list2);
 736     MY_FREE(point);
 737     MY_FREE(anchor);
 738     MY_FREE(length);
 739     MY_FREE(arrayi);
 740     MY_FREE(levels);
 741     MY_FREE(thirds);
 742 } /* direct_ */
 743