direct/DIRsubrout.c

   1 /* DIRsubrout.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 /* Table of constant values */
  10
  11 static integer c__1 = 1;
  12 static integer c__32 = 32;
  13 static integer c__0 = 0;
  14
  15 /* +-----------------------------------------------------------------------+ */
  16 /* | INTEGER Function DIRGetlevel                                          | */
  17 /* | Returns the level of the hyperrectangle. Depending on the value of the| */
  18 /* | global variable JONES. IF JONES equals 0, the level is given by       | */
  19 /* |               kN + p, where the rectangle has p sides with a length of| */
  20 /* |             1/3^(k+1), and N-p sides with a length of 1/3^k.          | */
  21 /* | If JONES equals 1, the level is the power of 1/3 of the length of the | */
  22 /* | longest side hyperrectangle.                                          | */
  23 /* |                                                                       | */
  24 /* | On Return :                                                           | */
  25 /* |    the maximal length                                                 | */
  26 /* |                                                                       | */
  27 /* | pos     -- the position of the midpoint in the array length           | */
  28 /* | length  -- the array with the dimensions                              | */
  29 /* | maxfunc -- the leading dimension of length                            | */
  30 /* | n     -- the dimension of the problem                                  | */
  31 /* |                                                                       | */
  32 /* +-----------------------------------------------------------------------+ */
  33 integer direct_dirgetlevel_(integer *pos, integer *length, integer *maxfunc, integer
  34         *n, integer jones)
  35 {
  36     /* System generated locals */
  37     integer length_dim1, length_offset, ret_val, i__1;
  38
  39     /* Local variables */
  40     integer i__, k, p, help;
  41
  42 /* JG 09/15/00 Added variable JONES (see above) */
  43     /* Parameter adjustments */
  44     length_dim1 = *n;
  45     length_offset = 1 + length_dim1;
  46     length -= length_offset;
  47
  48     /* Function Body */
  49     if (jones == 0) {
  50         help = length[*pos * length_dim1 + 1];
  51         k = help;
  52         p = 1;
  53         i__1 = *n;
  54         for (i__ = 2; i__ <= i__1; ++i__) {
  55             if (length[i__ + *pos * length_dim1] < k) {
  56                 k = length[i__ + *pos * length_dim1];
  57             }
  58             if (length[i__ + *pos * length_dim1] == help) {
  59                 ++p;
  60             }
  61 /* L100: */
  62         }
  63         if (k == help) {
  64             ret_val = k * *n + *n - p;
  65         } else {
  66             ret_val = k * *n + p;
  67         }
  68     } else {
  69         help = length[*pos * length_dim1 + 1];
  70         i__1 = *n;
  71         for (i__ = 2; i__ <= i__1; ++i__) {
  72             if (length[i__ + *pos * length_dim1] < help) {
  73                 help = length[i__ + *pos * length_dim1];
  74             }
  75 /* L10: */
  76         }
  77         ret_val = help;
  78     }
  79     return ret_val;
  80 } /* dirgetlevel_ */
  81
  82 /* +-----------------------------------------------------------------------+ */
  83 /* | Program       : Direct.f (subfile DIRsubrout.f)                       | */
  84 /* | Last modified : 07-16-2001                                            | */
  85 /* | Written by    : Joerg Gablonsky                                       | */
  86 /* | Subroutines used by the algorithm DIRECT.                             | */
  87 /* +-----------------------------------------------------------------------+ */
  88 /* +-----------------------------------------------------------------------+ */
  89 /* |    SUBROUTINE DIRChoose                                               | */
  90 /* |    Decide, which is the next sampling point.                          | */
  91 /* |    Changed 09/25/00 JG                                                | */
  92 /* |         Added maxdiv to call and changed S to size maxdiv.            | */
  93 /* |    Changed 01/22/01 JG                                                | */
  94 /* |         Added Ifeasiblef to call to keep track if a feasible point has| */
  95 /* |         been found.                                                   | */
  96 /* |    Changed 07/16/01 JG                                                | */
  97 /* |         Changed if statement to prevent run-time errors.              |
  98                  | */
  99 /* +-----------------------------------------------------------------------+ */
 100 /* Subroutine */ void direct_dirchoose_(integer *anchor, integer *s, integer *actdeep,
 101          doublereal *f, doublereal *minf, doublereal epsrel, doublereal epsabs, doublereal *thirds,
 102          integer *maxpos, integer *length, integer *maxfunc, const integer *maxdeep,
 103          const integer *maxdiv, integer *n, FILE *logfile,
 104         integer *cheat, doublereal *kmax, integer *ifeasiblef, integer jones)
 105 {
 106     /* System generated locals */
 107     integer s_dim1, s_offset, length_dim1, length_offset, i__1;
 108
 109     /* Local variables */
 110     integer i__, j, k;
 111     doublereal helplower;
 112     integer i___, j___;
 113     doublereal helpgreater;
 114     integer novaluedeep = 0;
 115     doublereal help2;
 116     integer novalue;
 117
 118     /* Parameter adjustments */
 119     f -= 3;
 120     ++anchor;
 121     s_dim1 = *maxdiv;
 122     s_offset = 1 + s_dim1;
 123     s -= s_offset;
 124     length_dim1 = *n;
 125     length_offset = 1 + length_dim1;
 126     length -= length_offset;
 127
 128     /* Function Body */
 129     helplower = HUGE_VAL;
 130     helpgreater = 0.;
 131     k = 1;
 132     if (*ifeasiblef >= 1) {
 133         i__1 = *actdeep;
 134         for (j = 0; j <= i__1; ++j) {
 135             if (anchor[j] > 0) {
 136                 s[k + s_dim1] = anchor[j];
 137                 s[k + (s_dim1 << 1)] = direct_dirgetlevel_(&s[k + s_dim1], &length[
 138                         length_offset], maxfunc, n, jones);
 139                 goto L12;
 140             }
 141 /* L1001: */
 142         }
 143 L12:
 144         ++k;
 145         *maxpos = 1;
 146         return;
 147     } else {
 148         i__1 = *actdeep;
 149         for (j = 0; j <= i__1; ++j) {
 150             if (anchor[j] > 0) {
 151                 s[k + s_dim1] = anchor[j];
 152                 s[k + (s_dim1 << 1)] = direct_dirgetlevel_(&s[k + s_dim1], &length[
 153                         length_offset], maxfunc, n, jones);
 154                 ++k;
 155             }
 156 /* L10: */
 157         }
 158     }
 159     novalue = 0;
 160     if (anchor[-1] > 0) {
 161         novalue = anchor[-1];
 162         novaluedeep = direct_dirgetlevel_(&novalue, &length[length_offset], maxfunc,
 163                 n, jones);
 164     }
 165     *maxpos = k - 1;
 166     i__1 = *maxdeep;
 167     for (j = k - 1; j <= i__1; ++j) {
 168         s[k + s_dim1] = 0;
 169 /* L11: */
 170     }
 171     for (j = *maxpos; j >= 1; --j) {
 172         helplower = HUGE_VAL;
 173         helpgreater = 0.;
 174         j___ = s[j + s_dim1];
 175         i__1 = j - 1;
 176         for (i__ = 1; i__ <= i__1; ++i__) {
 177             i___ = s[i__ + s_dim1];
 178 /* +-----------------------------------------------------------------------+ */
 179 /* | JG 07/16/01 Changed IF statement into two to prevent run-time errors  | */
 180 /* |             which could occur if the compiler checks the second       | */
 181 /* |             expression in an .AND. statement although the first       | */
 182 /* |             statement is already not true.                            | */
 183 /* +-----------------------------------------------------------------------+ */
 184             if (i___ > 0 && ! (i__ == j)) {
 185                 if (f[(i___ << 1) + 2] <= 1.) {
 186                     help2 = thirds[s[i__ + (s_dim1 << 1)]] - thirds[s[j + (
 187                             s_dim1 << 1)]];
 188                     help2 = (f[(i___ << 1) + 1] - f[(j___ << 1) + 1]) / help2;
 189                     if (help2 <= 0.) {
 190                          if (logfile)
 191                               fprintf(logfile, "thirds > 0, help2 <= 0\n");
 192                         goto L60;
 193                     }
 194                     if (help2 < helplower) {
 195                          if (logfile)
 196                               fprintf(logfile, "helplower = %g\n", help2);
 197                         helplower = help2;
 198                     }
 199                 }
 200             }
 201 /* L30: */
 202         }
 203         i__1 = *maxpos;
 204         for (i__ = j + 1; i__ <= i__1; ++i__) {
 205             i___ = s[i__ + s_dim1];
 206 /* +-----------------------------------------------------------------------+ */
 207 /* | JG 07/16/01 Changed IF statement into two to prevent run-time errors  | */
 208 /* |             which could occur if the compiler checks the second       | */
 209 /* |             expression in an .AND. statement although the first       | */
 210 /* |             statement is already not true.                            | */
 211 /* +-----------------------------------------------------------------------+ */
 212             if (i___ > 0 && ! (i__ == j)) {
 213                 if (f[(i___ << 1) + 2] <= 1.) {
 214                     help2 = thirds[s[i__ + (s_dim1 << 1)]] - thirds[s[j + (
 215                             s_dim1 << 1)]];
 216                     help2 = (f[(i___ << 1) + 1] - f[(j___ << 1) + 1]) / help2;
 217                     if (help2 <= 0.) {
 218                         if (logfile)
 219                              fprintf(logfile, "thirds < 0, help2 <= 0\n");
 220                         goto L60;
 221                     }
 222                     if (help2 > helpgreater) {
 223                         if (logfile)
 224                               fprintf(logfile, "helpgreater = %g\n", help2);
 225                         helpgreater = help2;
 226                     }
 227                 }
 228             }
 229 /* L31: */
 230         }
 231         if (helpgreater <= helplower) {
 232             if (*cheat == 1 && helplower > *kmax) {
 233                 helplower = *kmax;
 234             }
 235             if (f[(j___ << 1) + 1] - helplower * thirds[s[j + (s_dim1 << 1)]] >
 236                      MIN(*minf - epsrel * fabs(*minf),
 237                          *minf - epsabs)) {
 238                 if (logfile)
 239                      fprintf(logfile, "> minf - epslminfl\n");
 240                 goto L60;
 241             }
 242         } else {
 243             if (logfile)
 244                  fprintf(logfile, "helpgreater > helplower: %g  %g  %g\n",
 245                          helpgreater, helplower, helpgreater - helplower);
 246             goto L60;
 247         }
 248         goto L40;
 249 L60:
 250         s[j + s_dim1] = 0;
 251 L40:
 252         ;
 253     }
 254     if (novalue > 0) {
 255         ++(*maxpos);
 256         s[*maxpos + s_dim1] = novalue;
 257         s[*maxpos + (s_dim1 << 1)] = novaluedeep;
 258     }
 259 } /* dirchoose_ */
 260
 261 /* +-----------------------------------------------------------------------+ */
 262 /* |    SUBROUTINE DIRDoubleInsert                                         | */
 263 /* |      Routine to make sure that if there are several potential optimal | */
 264 /* |      hyperrectangles of the same level (i.e. hyperrectangles that have| */
 265 /* |      the same level and the same function value at the center), all of| */
 266 /* |      them are divided. This is the way as originally described in     | */
 267 /* |      Jones et.al.                                                     | */
 268 /* | JG 07/16/01 Added errorflag to calling sequence. We check if more     | */
 269 /* |             we reach the capacity of the array S. If this happens, we | */
 270 /* |             return to the main program with an error.                 | */
 271 /* +-----------------------------------------------------------------------+ */
 272 /* Subroutine */ void direct_dirdoubleinsert_(integer *anchor, integer *s, integer *
 273         maxpos, integer *point, doublereal *f, const integer *maxdeep, integer *
 274         maxfunc, const integer *maxdiv, integer *ierror)
 275 {
 276     /* System generated locals */
 277     integer s_dim1, s_offset, i__1;
 278
 279     /* Local variables */
 280     integer i__, oldmaxpos, pos, help, iflag, actdeep;
 281
 282 /* +-----------------------------------------------------------------------+ */
 283 /* | JG 07/16/01 Added flag to prevent run time-errors on some systems.    | */
 284 /* +-----------------------------------------------------------------------+ */
 285     /* Parameter adjustments */
 286     ++anchor;
 287     f -= 3;
 288     --point;
 289     s_dim1 = *maxdiv;
 290     s_offset = 1 + s_dim1;
 291     s -= s_offset;
 292
 293     /* Function Body */
 294     oldmaxpos = *maxpos;
 295     i__1 = oldmaxpos;
 296     for (i__ = 1; i__ <= i__1; ++i__) {
 297         if (s[i__ + s_dim1] > 0) {
 298             actdeep = s[i__ + (s_dim1 << 1)];
 299             help = anchor[actdeep];
 300             pos = point[help];
 301             iflag = 0;
 302 /* +-----------------------------------------------------------------------+ */
 303 /* | JG 07/16/01 Added flag to prevent run time-errors on some systems. On | */
 304 /* |             some systems the second conditions in an AND statement is | */
 305 /* |             evaluated even if the first one is already not true.      | */
 306 /* +-----------------------------------------------------------------------+ */
 307             while(pos > 0 && iflag == 0) {
 308                 if (f[(pos << 1) + 1] - f[(help << 1) + 1] <= 1e-13) {
 309                     if (*maxpos < *maxdiv) {
 310                         ++(*maxpos);
 311                         s[*maxpos + s_dim1] = pos;
 312                         s[*maxpos + (s_dim1 << 1)] = actdeep;
 313                         pos = point[pos];
 314                     } else {
 315 /* +-----------------------------------------------------------------------+ */
 316 /* | JG 07/16/01 Maximum number of elements possible in S has been reached!| */
 317 /* +-----------------------------------------------------------------------+ */
 318                         *ierror = -6;
 319                         return;
 320                     }
 321                 } else {
 322                     iflag = 1;
 323                 }
 324             }
 325         }
 326 /* L10: */
 327     }
 328 } /* dirdoubleinsert_ */
 329
 330 /* +-----------------------------------------------------------------------+ */
 331 /* | INTEGER Function GetmaxDeep                                           | */
 332 /* | function to get the maximal length (1/length) of the n-dimensional    | */
 333 /* | rectangle with midpoint pos.                                          | */
 334 /* |                                                                       | */
 335 /* | On Return :                                                           | */
 336 /* |    the maximal length                                                 | */
 337 /* |                                                                       | */
 338 /* | pos     -- the position of the midpoint in the array length           | */
 339 /* | length  -- the array with the dimensions                              | */
 340 /* | maxfunc -- the leading dimension of length                            | */
 341 /* | n     -- the dimension of the problem                                  | */
 342 /* |                                                                       | */
 343 /* +-----------------------------------------------------------------------+ */
 344 integer direct_dirgetmaxdeep_(integer *pos, integer *length, integer *maxfunc,
 345         integer *n)
 346 {
 347     /* System generated locals */
 348     integer length_dim1, length_offset, i__1, i__2, i__3;
 349
 350     /* Local variables */
 351     integer i__, help;
 352
 353     /* Parameter adjustments */
 354     length_dim1 = *n;
 355     length_offset = 1 + length_dim1;
 356     length -= length_offset;
 357
 358     /* Function Body */
 359     help = length[*pos * length_dim1 + 1];
 360     i__1 = *n;
 361     for (i__ = 2; i__ <= i__1; ++i__) {
 362 /* Computing MIN */
 363         i__2 = help, i__3 = length[i__ + *pos * length_dim1];
 364         help = MIN(i__2,i__3);
 365 /* L10: */
 366     }
 367     return help;
 368 } /* dirgetmaxdeep_ */
 369
 370 static integer isinbox_(doublereal *x, doublereal *a, doublereal *b, integer *n,
 371         integer *lmaxdim)
 372 {
 373     /* System generated locals */
 374     integer ret_val, i__1;
 375
 376     /* Local variables */
 377     integer outofbox, i__;
 378
 379     /* Parameter adjustments */
 380     --b;
 381     --a;
 382     --x;
 383
 384     /* Function Body */
 385     outofbox = 1;
 386     i__1 = *n;
 387     for (i__ = 1; i__ <= i__1; ++i__) {
 388         if (a[i__] > x[i__] || b[i__] < x[i__]) {
 389             outofbox = 0;
 390             goto L1010;
 391         }
 392 /* L1000: */
 393     }
 394 L1010:
 395     ret_val = outofbox;
 396     return ret_val;
 397 } /* isinbox_ */
 398
 399 /* +-----------------------------------------------------------------------+ */
 400 /* | JG Added 09/25/00                                                     | */
 401 /* |                                                                       | */
 402 /* |                       SUBROUTINE DIRResortlist                        | */
 403 /* |                                                                       | */
 404 /* | Resort the list so that the infeasible point is in the list with the  | */
 405 /* | replaced value.                                                       | */
 406 /* +-----------------------------------------------------------------------+ */
 407 /* Subroutine */ static void dirresortlist_(integer *replace, integer *anchor,
 408         doublereal *f, integer *point, integer *length, integer *n, integer *
 409         maxfunc, integer *maxdim, const integer *maxdeep, FILE *logfile,
 410                                             integer jones)
 411 {
 412     /* System generated locals */
 413     integer length_dim1, length_offset, i__1;
 414
 415     /* Local variables */
 416     integer i__, l, pos;
 417     integer start;
 418
 419 /* +-----------------------------------------------------------------------+ */
 420 /* | Get the length of the hyper rectangle with infeasible mid point and   | */
 421 /* | Index of the corresponding list.                                      | */
 422 /* +-----------------------------------------------------------------------+ */
 423 /* JG 09/25/00 Replaced with DIRgetlevel */
 424 /*      l = DIRgetmaxDeep(replace,length,maxfunc,n) */
 425     /* Parameter adjustments */
 426     --point;
 427     f -= 3;
 428     length_dim1 = *n;
 429     length_offset = 1 + length_dim1;
 430     length -= length_offset;
 431     ++anchor;
 432
 433     /* Function Body */
 434     l = direct_dirgetlevel_(replace, &length[length_offset], maxfunc, n, jones);
 435     start = anchor[l];
 436 /* +-----------------------------------------------------------------------+ */
 437 /* | If the hyper rectangle with infeasibel midpoint is already the start  | */
 438 /* | of the list, give out message, nothing to do.                         | */
 439 /* +-----------------------------------------------------------------------+ */
 440     if (*replace == start) {
 441 /*         write(logfile,*) 'No resorting of list necessarry, since new ', */
 442 /*     + 'point is already anchor of list .',l */
 443     } else {
 444 /* +-----------------------------------------------------------------------+ */
 445 /* | Take the hyper rectangle with infeasible midpoint out of the list.    | */
 446 /* +-----------------------------------------------------------------------+ */
 447         pos = start;
 448         i__1 = *maxfunc;
 449         for (i__ = 1; i__ <= i__1; ++i__) {
 450             if (point[pos] == *replace) {
 451                 point[pos] = point[*replace];
 452                 goto L20;
 453             } else {
 454                 pos = point[pos];
 455             }
 456             if (pos == 0) {
 457                 if (logfile)
 458                      fprintf(logfile, "Error in DIRREsortlist: "
 459                              "We went through the whole list\n"
 460                              "and could not find the point to replace!!\n");
 461                 goto L20;
 462             }
 463 /* L10: */
 464         }
 465 /* +-----------------------------------------------------------------------+ */
 466 /* | If the anchor of the list has a higher value than the value of a      | */
 467 /* | nearby point, put the infeasible point at the beginning of the list.  | */
 468 /* +-----------------------------------------------------------------------+ */
 469 L20:
 470         if (f[(start << 1) + 1] > f[(*replace << 1) + 1]) {
 471             anchor[l] = *replace;
 472             point[*replace] = start;
 473 /*            write(logfile,*) 'Point is replacing current anchor for ' */
 474 /*     +             , 'this list ',l,replace,start */
 475         } else {
 476 /* +-----------------------------------------------------------------------+ */
 477 /* | Insert the point into the list according to its (replaced) function   | */
 478 /* | value.                                                                | */
 479 /* +-----------------------------------------------------------------------+ */
 480             pos = start;
 481             i__1 = *maxfunc;
 482             for (i__ = 1; i__ <= i__1; ++i__) {
 483 /* +-----------------------------------------------------------------------+ */
 484 /* | The point has to be added at the end of the list.                     | */
 485 /* +-----------------------------------------------------------------------+ */
 486                 if (point[pos] == 0) {
 487                     point[*replace] = point[pos];
 488                     point[pos] = *replace;
 489 /*                  write(logfile,*) 'Point is added at the end of the ' */
 490 /*     +             , 'list ',l, replace */
 491                     goto L40;
 492                 } else {
 493                     if (f[(point[pos] << 1) + 1] > f[(*replace << 1) + 1]) {
 494                         point[*replace] = point[pos];
 495                         point[pos] = *replace;
 496 /*                     write(logfile,*) 'There are points with a higher ' */
 497 /*     +               ,'f-value in the list ',l,replace, pos */
 498                         goto L40;
 499                     }
 500                     pos = point[pos];
 501                 }
 502 /* L30: */
 503             }
 504 L40:
 505             pos = pos;
 506         }
 507     }
 508 } /* dirresortlist_ */
 509
 510 /* +-----------------------------------------------------------------------+ */
 511 /* | JG Added 09/25/00                                                     | */
 512 /* |                       SUBROUTINE DIRreplaceInf                        | */
 513 /* |                                                                       | */
 514 /* | Find out if there are infeasible points which are near feasible ones. | */
 515 /* | If this is the case, replace the function value at the center of the  | */
 516 /* | hyper rectangle by the lowest function value of a nearby function.    | */
 517 /* +-----------------------------------------------------------------------+ */
 518 /* Subroutine */ void direct_dirreplaceinf_(integer *free, integer *freeold,
 519         doublereal *f, doublereal *c__, doublereal *thirds, integer *length,
 520         integer *anchor, integer *point, doublereal *c1, doublereal *c2,
 521         integer *maxfunc, const integer *maxdeep, integer *maxdim, integer *n,
 522         FILE *logfile, doublereal *fmax, integer jones)
 523 {
 524     /* System generated locals */
 525     integer c_dim1, c_offset, length_dim1, length_offset, i__1, i__2, i__3;
 526     doublereal d__1, d__2;
 527
 528     /* Local variables */
 529     doublereal a[32], b[32];
 530     integer i__, j, k, l;
 531     doublereal x[32], sidelength;
 532     integer help;
 533
 534 /* +-----------------------------------------------------------------------+ */
 535 /* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
 536 /* +-----------------------------------------------------------------------+ */
 537     /* Parameter adjustments */
 538     --point;
 539     f -= 3;
 540     ++anchor;
 541     length_dim1 = *maxdim;
 542     length_offset = 1 + length_dim1;
 543     length -= length_offset;
 544     c_dim1 = *maxdim;
 545     c_offset = 1 + c_dim1;
 546     c__ -= c_offset;
 547     --c2;
 548     --c1;
 549
 550     /* Function Body */
 551     i__1 = *free - 1;
 552     for (i__ = 1; i__ <= i__1; ++i__) {
 553         if (f[(i__ << 1) + 2] > 0.) {
 554 /* +-----------------------------------------------------------------------+ */
 555 /* | Get the maximum side length of the hyper rectangle and then set the   | */
 556 /* | new side length to this lengths times the growth factor.              | */
 557 /* +-----------------------------------------------------------------------+ */
 558             help = direct_dirgetmaxdeep_(&i__, &length[length_offset], maxfunc, n);
 559             sidelength = thirds[help] * 2.;
 560 /* +-----------------------------------------------------------------------+ */
 561 /* | Set the Center and the upper and lower bounds of the rectangles.      | */
 562 /* +-----------------------------------------------------------------------+ */
 563             i__2 = *n;
 564             for (j = 1; j <= i__2; ++j) {
 565                 sidelength = thirds[length[i__ + j * length_dim1]];
 566                 a[j - 1] = c__[j + i__ * c_dim1] - sidelength;
 567                 b[j - 1] = c__[j + i__ * c_dim1] + sidelength;
 568 /* L20: */
 569             }
 570 /* +-----------------------------------------------------------------------+ */
 571 /* | The function value is reset to 'Inf', since it may have been changed  | */
 572 /* | in an earlier iteration and now the feasible point which was close    | */
 573 /* | is not close anymore (since the hyper rectangle surrounding the       | */
 574 /* | current point may have shrunk).                                       | */
 575 /* +-----------------------------------------------------------------------+ */
 576             f[(i__ << 1) + 1] = HUGE_VAL;
 577             f[(i__ << 1) + 2] = 2.;
 578 /* +-----------------------------------------------------------------------+ */
 579 /* | Check if any feasible point is near this infeasible point.            | */
 580 /* +-----------------------------------------------------------------------+ */
 581             i__2 = *free - 1;
 582             for (k = 1; k <= i__2; ++k) {
 583 /* +-----------------------------------------------------------------------+ */
 584 /* | If the point k is feasible, check if it is near.                      | */
 585 /* +-----------------------------------------------------------------------+ */
 586                 if (f[(k << 1) + 2] == 0.) {
 587 /* +-----------------------------------------------------------------------+ */
 588 /* | Copy the coordinates of the point k into x.                           | */
 589 /* +-----------------------------------------------------------------------+ */
 590                     i__3 = *n;
 591                     for (l = 1; l <= i__3; ++l) {
 592                         x[l - 1] = c__[l + k * c_dim1];
 593 /* L40: */
 594                     }
 595 /* +-----------------------------------------------------------------------+ */
 596 /* | Check if the point k is near the infeasible point, if so, replace the | */
 597 /* | value at */
 598 /* +-----------------------------------------------------------------------+ */
 599                     if (isinbox_(x, a, b, n, &c__32) == 1) {
 600 /* Computing MIN */
 601                          d__1 = f[(i__ << 1) + 1], d__2 = f[(k << 1) + 1];
 602                          f[(i__ << 1) + 1] = MIN(d__1,d__2);
 603                          f[(i__ << 1) + 2] = 1.;
 604                     }
 605                 }
 606 /* L30: */
 607             }
 608             if (f[(i__ << 1) + 2] == 1.) {
 609                 f[(i__ << 1) + 1] += (d__1 = f[(i__ << 1) + 1], fabs(d__1)) *
 610                         1e-6f;
 611                 i__2 = *n;
 612                 for (l = 1; l <= i__2; ++l) {
 613                     x[l - 1] = c__[l + i__ * c_dim1] * c1[l] + c__[l + i__ *
 614                             c_dim1] * c2[l];
 615 /* L200: */
 616                 }
 617                 dirresortlist_(&i__, &anchor[-1], &f[3], &point[1],
 618                                &length[length_offset], n,
 619                                maxfunc, maxdim, maxdeep, logfile, jones);
 620             } else {
 621 /* +-----------------------------------------------------------------------+ */
 622 /* | JG 01/22/01                                                           | */
 623 /* | Replaced fixed value for infeasible points with maximum value found,  | */
 624 /* | increased by 1.                                                       | */
 625 /* +-----------------------------------------------------------------------+ */
 626                 if (! (*fmax == f[(i__ << 1) + 1])) {
 627 /* Computing MAX */
 628                     d__1 = *fmax + 1., d__2 = f[(i__ << 1) + 1];
 629                     f[(i__ << 1) + 1] = MAX(d__1,d__2);
 630                 }
 631             }
 632         }
 633 /* L10: */
 634     }
 635 /* L1000: */
 636 } /* dirreplaceinf_ */
 637
 638 /* +-----------------------------------------------------------------------+ */
 639 /* |    SUBROUTINE DIRInsert                                               | */
 640 /* +-----------------------------------------------------------------------+ */
 641 /* Subroutine */ static void dirinsert_(integer *start, integer *ins, integer *point,
 642         doublereal *f, integer *maxfunc)
 643 {
 644     /* System generated locals */
 645     integer i__1;
 646
 647     /* Local variables */
 648     integer i__, help;
 649
 650 /* JG 09/17/00 Rewrote this routine. */
 651 /*      DO 10,i = 1,maxfunc */
 652 /*        IF (f(ins,1) .LT. f(point(start),1)) THEN */
 653 /*          help = point(start) */
 654 /*          point(start) = ins */
 655 /*          point(ins) = help */
 656 /*          GOTO 20 */
 657 /*        END IF */
 658 /*        IF (point(start) .EQ. 0) THEN */
 659 /*           point(start) = ins */
 660 /*           point(ins) = 0 */
 661 /*           GOTO 20 */
 662 /*        END IF */
 663 /*        start = point(start) */
 664 /* 10    CONTINUE */
 665 /* 20    END */
 666     /* Parameter adjustments */
 667     f -= 3;
 668     --point;
 669
 670     /* Function Body */
 671     i__1 = *maxfunc;
 672     for (i__ = 1; i__ <= i__1; ++i__) {
 673         if (point[*start] == 0) {
 674             point[*start] = *ins;
 675             point[*ins] = 0;
 676             return;
 677         } else if (f[(*ins << 1) + 1] < f[(point[*start] << 1) + 1]) {
 678              help = point[*start];
 679              point[*start] = *ins;
 680              point[*ins] = help;
 681              return;
 682         }
 683         *start = point[*start];
 684 /* L10: */
 685     }
 686 } /* dirinsert_ */
 687
 688 /* +-----------------------------------------------------------------------+ */
 689 /* |    SUBROUTINE DIRInsertList                                           | */
 690 /* |    Changed 02-24-2000                                                 | */
 691 /* |      Got rid of the distinction between feasible and infeasible points| */
 692 /* |      I could do this since infeasible points get set to a high        | */
 693 /* |      function value, which may be replaced by a function value of a   | */
 694 /* |      nearby function at the end of the main loop.                     | */
 695 /* +-----------------------------------------------------------------------+ */
 696 /* Subroutine */ void direct_dirinsertlist_(integer *new__, integer *anchor, integer *
 697         point, doublereal *f, integer *maxi, integer *length, integer *
 698         maxfunc, const integer *maxdeep, integer *n, integer *samp,
 699                                             integer jones)
 700 {
 701     /* System generated locals */
 702     integer length_dim1, length_offset, i__1;
 703
 704     /* Local variables */
 705     integer j;
 706     integer pos;
 707     integer pos1, pos2, deep;
 708
 709 /* JG 09/24/00 Changed this to Getlevel */
 710     /* Parameter adjustments */
 711     f -= 3;
 712     --point;
 713     ++anchor;
 714     length_dim1 = *n;
 715     length_offset = 1 + length_dim1;
 716     length -= length_offset;
 717
 718     /* Function Body */
 719     i__1 = *maxi;
 720     for (j = 1; j <= i__1; ++j) {
 721         pos1 = *new__;
 722         pos2 = point[pos1];
 723         *new__ = point[pos2];
 724 /* JG 09/24/00 Changed this to Getlevel */
 725 /*        deep = DIRGetMaxdeep(pos1,length,maxfunc,n) */
 726         deep = direct_dirgetlevel_(&pos1, &length[length_offset], maxfunc, n, jones);
 727         if (anchor[deep] == 0) {
 728             if (f[(pos2 << 1) + 1] < f[(pos1 << 1) + 1]) {
 729                 anchor[deep] = pos2;
 730                 point[pos2] = pos1;
 731                 point[pos1] = 0;
 732             } else {
 733                 anchor[deep] = pos1;
 734                 point[pos2] = 0;
 735             }
 736         } else {
 737             pos = anchor[deep];
 738             if (f[(pos2 << 1) + 1] < f[(pos1 << 1) + 1]) {
 739                 if (f[(pos2 << 1) + 1] < f[(pos << 1) + 1]) {
 740                     anchor[deep] = pos2;
 741 /* JG 08/30/00 Fixed bug. Sorting was not correct when */
 742 /*      f(1,pos2) < f(1,pos1) < f(1,pos) */
 743                     if (f[(pos1 << 1) + 1] < f[(pos << 1) + 1]) {
 744                         point[pos2] = pos1;
 745                         point[pos1] = pos;
 746                     } else {
 747                         point[pos2] = pos;
 748                         dirinsert_(&pos, &pos1, &point[1], &f[3], maxfunc);
 749                     }
 750                 } else {
 751                     dirinsert_(&pos, &pos2, &point[1], &f[3], maxfunc);
 752                     dirinsert_(&pos, &pos1, &point[1], &f[3], maxfunc);
 753                 }
 754             } else {
 755                 if (f[(pos1 << 1) + 1] < f[(pos << 1) + 1]) {
 756 /* JG 08/30/00 Fixed bug. Sorting was not correct when */
 757 /*      f(pos1,1) < f(pos2,1) < f(pos,1) */
 758                     anchor[deep] = pos1;
 759                     if (f[(pos << 1) + 1] < f[(pos2 << 1) + 1]) {
 760                         point[pos1] = pos;
 761                         dirinsert_(&pos, &pos2, &point[1], &f[3], maxfunc);
 762                     } else {
 763                         point[pos1] = pos2;
 764                         point[pos2] = pos;
 765                     }
 766                 } else {
 767                     dirinsert_(&pos, &pos1, &point[1], &f[3], maxfunc);
 768                     dirinsert_(&pos, &pos2, &point[1], &f[3], maxfunc);
 769                 }
 770             }
 771         }
 772 /* L10: */
 773     }
 774 /* JG 09/24/00 Changed this to Getlevel */
 775 /*      deep = DIRGetMaxdeep(samp,length,maxfunc,n) */
 776     deep = direct_dirgetlevel_(samp, &length[length_offset], maxfunc, n, jones);
 777     pos = anchor[deep];
 778     if (f[(*samp << 1) + 1] < f[(pos << 1) + 1]) {
 779         anchor[deep] = *samp;
 780         point[*samp] = pos;
 781     } else {
 782         dirinsert_(&pos, samp, &point[1], &f[3], maxfunc);
 783     }
 784 } /* dirinsertlist_ */
 785
 786 /* +-----------------------------------------------------------------------+ */
 787 /* |    SUBROUTINE DIRInsertList2  (Old way to do it.)                     | */
 788 /* +-----------------------------------------------------------------------+ */
 789 /* Subroutine */ static void dirinsertlist_2__(integer *start, integer *j, integer *k,
 790          integer *list2, doublereal *w, integer *maxi, integer *n)
 791 {
 792     /* System generated locals */
 793     integer list2_dim1, list2_offset, i__1;
 794
 795     /* Local variables */
 796     integer i__, pos;
 797
 798     /* Parameter adjustments */
 799     --w;
 800     list2_dim1 = *n;
 801     list2_offset = 1 + list2_dim1;
 802     list2 -= list2_offset;
 803
 804     /* Function Body */
 805     pos = *start;
 806     if (*start == 0) {
 807         list2[*j + list2_dim1] = 0;
 808         *start = *j;
 809         goto L50;
 810     }
 811     if (w[*start] > w[*j]) {
 812         list2[*j + list2_dim1] = *start;
 813         *start = *j;
 814     } else {
 815         i__1 = *maxi;
 816         for (i__ = 1; i__ <= i__1; ++i__) {
 817             if (list2[pos + list2_dim1] == 0) {
 818                 list2[*j + list2_dim1] = 0;
 819                 list2[pos + list2_dim1] = *j;
 820                 goto L50;
 821             } else {
 822                 if (w[*j] < w[list2[pos + list2_dim1]]) {
 823                     list2[*j + list2_dim1] = list2[pos + list2_dim1];
 824                     list2[pos + list2_dim1] = *j;
 825                     goto L50;
 826                 }
 827             }
 828             pos = list2[pos + list2_dim1];
 829 /* L10: */
 830         }
 831     }
 832 L50:
 833     list2[*j + (list2_dim1 << 1)] = *k;
 834 } /* dirinsertlist_2__ */
 835
 836 /* +-----------------------------------------------------------------------+ */
 837 /* |    SUBROUTINE DIRSearchmin                                            | */
 838 /* |    Search for the minimum in the list.                                ! */
 839 /* +-----------------------------------------------------------------------+ */
 840 /* Subroutine */ static void dirsearchmin_(integer *start, integer *list2, integer *
 841         pos, integer *k, integer *n)
 842 {
 843     /* System generated locals */
 844     integer list2_dim1, list2_offset;
 845
 846     /* Parameter adjustments */
 847     list2_dim1 = *n;
 848     list2_offset = 1 + list2_dim1;
 849     list2 -= list2_offset;
 850
 851     /* Function Body */
 852     *k = *start;
 853     *pos = list2[*start + (list2_dim1 << 1)];
 854     *start = list2[*start + list2_dim1];
 855 } /* dirsearchmin_ */
 856
 857 /* +-----------------------------------------------------------------------+ */
 858 /* |    SUBROUTINE DIRSamplepoints                                         | */
 859 /* |    Subroutine to sample the new points.                               | */
 860 /* +-----------------------------------------------------------------------+ */
 861 /* Subroutine */ void direct_dirsamplepoints_(doublereal *c__, integer *arrayi,
 862         doublereal *delta, integer *sample, integer *start, integer *length,
 863         FILE *logfile, doublereal *f, integer *free,
 864         integer *maxi, integer *point, doublereal *x, doublereal *l,
 865          doublereal *minf, integer *minpos, doublereal *u, integer *n,
 866         integer *maxfunc, const integer *maxdeep, integer *oops)
 867 {
 868     /* System generated locals */
 869     integer length_dim1, length_offset, c_dim1, c_offset, i__1, i__2;
 870
 871     /* Local variables */
 872     integer j, k, pos;
 873
 874
 875     /* Parameter adjustments */
 876     --u;
 877     --l;
 878     --x;
 879     --arrayi;
 880     --point;
 881     f -= 3;
 882     length_dim1 = *n;
 883     length_offset = 1 + length_dim1;
 884     length -= length_offset;
 885     c_dim1 = *n;
 886     c_offset = 1 + c_dim1;
 887     c__ -= c_offset;
 888
 889     /* Function Body */
 890     *oops = 0;
 891     pos = *free;
 892     *start = *free;
 893     i__1 = *maxi + *maxi;
 894     for (k = 1; k <= i__1; ++k) {
 895         i__2 = *n;
 896         for (j = 1; j <= i__2; ++j) {
 897             length[j + *free * length_dim1] = length[j + *sample *
 898                     length_dim1];
 899             c__[j + *free * c_dim1] = c__[j + *sample * c_dim1];
 900 /* L20: */
 901         }
 902         pos = *free;
 903         *free = point[*free];
 904         if (*free == 0) {
 905              if (logfile)
 906                   fprintf(logfile, "Error, no more free positions! "
 907                           "Increase maxfunc!\n");
 908             *oops = 1;
 909             return;
 910         }
 911 /* L10: */
 912     }
 913     point[pos] = 0;
 914     pos = *start;
 915     i__1 = *maxi;
 916     for (j = 1; j <= i__1; ++j) {
 917         c__[arrayi[j] + pos * c_dim1] = c__[arrayi[j] + *sample * c_dim1] + *
 918                 delta;
 919         pos = point[pos];
 920         c__[arrayi[j] + pos * c_dim1] = c__[arrayi[j] + *sample * c_dim1] - *
 921                 delta;
 922         pos = point[pos];
 923 /* L30: */
 924     }
 925     ASRT(pos <= 0);
 926 } /* dirsamplepoints_ */
 927
 928 /* +-----------------------------------------------------------------------+ */
 929 /* |    SUBROUTINE DIRDivide                                               | */
 930 /* |    Subroutine to divide the hyper rectangles according to the rules.  | */
 931 /* |    Changed 02-24-2000                                                 | */
 932 /* |      Replaced if statement by min (line 367)                          | */
 933 /* +-----------------------------------------------------------------------+ */
 934 /* Subroutine */ void direct_dirdivide_(integer *new__, integer *currentlength,
 935         integer *length, integer *point, integer *arrayi, integer *sample,
 936         integer *list2, doublereal *w, integer *maxi, doublereal *f, integer *
 937         maxfunc, const integer *maxdeep, integer *n)
 938 {
 939     /* System generated locals */
 940     integer length_dim1, length_offset, list2_dim1, list2_offset, i__1, i__2;
 941     doublereal d__1, d__2;
 942
 943     /* Local variables */
 944     integer i__, j, k, pos, pos2;
 945     integer start;
 946
 947
 948     /* Parameter adjustments */
 949     f -= 3;
 950     --point;
 951     --w;
 952     list2_dim1 = *n;
 953     list2_offset = 1 + list2_dim1;
 954     list2 -= list2_offset;
 955     --arrayi;
 956     length_dim1 = *n;
 957     length_offset = 1 + length_dim1;
 958     length -= length_offset;
 959
 960     /* Function Body */
 961     start = 0;
 962     pos = *new__;
 963     i__1 = *maxi;
 964     for (i__ = 1; i__ <= i__1; ++i__) {
 965         j = arrayi[i__];
 966         w[j] = f[(pos << 1) + 1];
 967         k = pos;
 968         pos = point[pos];
 969 /* Computing MIN */
 970         d__1 = f[(pos << 1) + 1], d__2 = w[j];
 971         w[j] = MIN(d__1,d__2);
 972         pos = point[pos];
 973         dirinsertlist_2__(&start, &j, &k, &list2[list2_offset], &w[1], maxi,
 974                 n);
 975 /* L10: */
 976     }
 977     ASRT(pos <= 0);
 978     i__1 = *maxi;
 979     for (j = 1; j <= i__1; ++j) {
 980         dirsearchmin_(&start, &list2[list2_offset], &pos, &k, n);
 981         pos2 = start;
 982         length[k + *sample * length_dim1] = *currentlength + 1;
 983         i__2 = *maxi - j + 1;
 984         for (i__ = 1; i__ <= i__2; ++i__) {
 985             length[k + pos * length_dim1] = *currentlength + 1;
 986             pos = point[pos];
 987             length[k + pos * length_dim1] = *currentlength + 1;
 988 /* JG 07/10/01 pos2 = 0 at the end of the 30-loop. Since we end */
 989 /*             the loop now, we do not need to reassign pos and pos2. */
 990             if (pos2 > 0) {
 991                 pos = list2[pos2 + (list2_dim1 << 1)];
 992                 pos2 = list2[pos2 + list2_dim1];
 993             }
 994 /* L30: */
 995         }
 996 /* L20: */
 997     }
 998 } /* dirdivide_ */
 999
1000 /* +-----------------------------------------------------------------------+ */
1001 /* |                                                                       | */
1002 /* |                       SUBROUTINE DIRINFCN                             | */
1003 /* |                                                                       | */
1004 /* | Subroutine DIRinfcn unscales the variable x for use in the            | */
1005 /* | user-supplied function evaluation subroutine fcn. After fcn returns   | */
1006 /* | to DIRinfcn, DIRinfcn then rescales x for use by DIRECT.              | */
1007 /* |                                                                       | */
1008 /* | On entry                                                              | */
1009 /* |                                                                       | */
1010 /* |        fcn -- The argument containing the name of the user-supplied   | */
1011 /* |               subroutine that returns values for the function to be   | */
1012 /* |               minimized.                                              | */
1013 /* |                                                                       | */
1014 /* |          x -- A double-precision vector of length n. The point at     | */
1015 /* |               which the derivative is to be evaluated.                | */
1016 /* |                                                                       | */
1017 /* |        xs1 -- A double-precision vector of length n. Used for         | */
1018 /* |               scaling and unscaling the vector x by DIRinfcn.         | */
1019 /* |                                                                       | */
1020 /* |        xs2 -- A double-precision vector of length n. Used for         | */
1021 /* |               scaling and unscaling the vector x by DIRinfcn.         | */
1022 /* |                                                                       | */
1023 /* |          n -- An integer. The dimension of the problem.               | */
1024 /* |       kret -- An Integer. If kret =  1, the point is infeasible,      | */
1025 /* |                              kret = -1, bad problem set up,           | */
1026 /* |                              kret =  0, feasible.                     | */
1027 /* |                                                                       | */
1028 /* | On return                                                             | */
1029 /* |                                                                       | */
1030 /* |          f -- A double-precision scalar.                              | */
1031 /* |                                                                       | */
1032 /* | Subroutines and Functions                                             | */
1033 /* |                                                                       | */
1034 /* | The subroutine whose name is passed through the argument fcn.         | */
1035 /* |                                                                       | */
1036 /* +-----------------------------------------------------------------------+ */
1037 /* Subroutine */ void direct_dirinfcn_(fp fcn, doublereal *x, doublereal *c1,
1038         doublereal *c2, integer *n, doublereal *f, integer *flag__,
1039                                        void *fcn_data)
1040 {
1041     /* System generated locals */
1042     integer i__1;
1043
1044     /* Local variables */
1045     integer i__;
1046
1047 /* +-----------------------------------------------------------------------+ */
1048 /* | Variables to pass user defined data to the function to be optimized.  | */
1049 /* +-----------------------------------------------------------------------+ */
1050 /* +-----------------------------------------------------------------------+ */
1051 /* | Unscale the variable x.                                               | */
1052 /* +-----------------------------------------------------------------------+ */
1053     /* Parameter adjustments */
1054     --c2;
1055     --c1;
1056     --x;
1057
1058     /* Function Body */
1059     i__1 = *n;
1060     for (i__ = 1; i__ <= i__1; ++i__) {
1061         x[i__] = (x[i__] + c2[i__]) * c1[i__];
1062 /* L20: */
1063     }
1064 /* +-----------------------------------------------------------------------+ */
1065 /* | Call the function-evaluation subroutine fcn.                          | */
1066 /* +-----------------------------------------------------------------------+ */
1067     *flag__ = 0;
1068     *f = fcn(*n, &x[1], flag__, fcn_data);
1069 /* +-----------------------------------------------------------------------+ */
1070 /* | Rescale the variable x.                                               | */
1071 /* +-----------------------------------------------------------------------+ */
1072     i__1 = *n;
1073     for (i__ = 1; i__ <= i__1; ++i__) {
1074         x[i__] = x[i__] / c1[i__] - c2[i__];
1075 /* L30: */
1076     }
1077 } /* dirinfcn_ */
1078
1079 /* +-----------------------------------------------------------------------+ */
1080 /* |    SUBROUTINE DIRGet_I                                                | */
1081 /* +-----------------------------------------------------------------------+ */
1082 /* Subroutine */ void direct_dirget_i__(integer *length, integer *pos, integer *
1083         arrayi, integer *maxi, integer *n, integer *maxfunc)
1084 {
1085     /* System generated locals */
1086     integer length_dim1, length_offset, i__1;
1087
1088     /* Local variables */
1089     integer i__, j, help;
1090
1091     /* Parameter adjustments */
1092     --arrayi;
1093     length_dim1 = *n;
1094     length_offset = 1 + length_dim1;
1095     length -= length_offset;
1096
1097     /* Function Body */
1098     j = 1;
1099     help = length[*pos * length_dim1 + 1];
1100     i__1 = *n;
1101     for (i__ = 2; i__ <= i__1; ++i__) {
1102         if (length[i__ + *pos * length_dim1] < help) {
1103             help = length[i__ + *pos * length_dim1];
1104         }
1105 /* L10: */
1106     }
1107     i__1 = *n;
1108     for (i__ = 1; i__ <= i__1; ++i__) {
1109         if (length[i__ + *pos * length_dim1] == help) {
1110             arrayi[j] = i__;
1111             ++j;
1112         }
1113 /* L20: */
1114     }
1115     *maxi = j - 1;
1116 } /* dirget_i__ */
1117
1118 /* +-----------------------------------------------------------------------+ */
1119 /* |    SUBROUTINE DIRInit                                                 | */
1120 /* |    Initialise all needed variables and do the first run of the        | */
1121 /* |    algorithm.                                                         | */
1122 /* |    Changed 02/24/2000                                                 | */
1123 /* |       Changed fcn Double precision to fcn external!                   | */
1124 /* |    Changed 09/15/2000                                                 | */
1125 /* |       Added distinction between Jones way to characterize rectangles  | */
1126 /* |       and our way. Common variable JONES controls which way we use.   | */
1127 /* |          JONES = 0    Jones way (Distance from midpoint to corner)    | */
1128 /* |          JONES = 1    Our way (Length of longest side)                | */
1129 /* |    Changed 09/24/00                                                   | */
1130 /* |       Added array levels. Levels contain the values to characterize   | */
1131 /* |       the hyperrectangles.                                            | */
1132 /* |    Changed 01/22/01                                                   | */
1133 /* |       Added variable fmax to keep track of maximum value found.       | */
1134 /* |       Added variable Ifeasiblef to keep track if feasibel point has   | */
1135 /* |       been found.                                                     | */
1136 /* |    Changed 01/23/01                                                   | */
1137 /* |       Added variable Ierror to keep track of errors.                  | */
1138 /* +-----------------------------------------------------------------------+ */
1139 /* Subroutine */ void direct_dirinit_(doublereal *f, fp fcn, doublereal *c__,
1140         integer *length, integer *actdeep, integer *point, integer *anchor,
1141         integer *free, FILE *logfile, integer *arrayi,
1142         integer *maxi, integer *list2, doublereal *w, doublereal *x,
1143         doublereal *l, doublereal *u, doublereal *minf, integer *minpos,
1144         doublereal *thirds, doublereal *levels, integer *maxfunc, const integer *
1145         maxdeep, integer *n, integer *maxor, doublereal *fmax, integer *
1146         ifeasiblef, integer *iinfeasible, integer *ierror, void *fcndata,
1147         integer jones, double starttime, double maxtime, int *force_stop)
1148 {
1149     /* System generated locals */
1150     integer c_dim1, c_offset, length_dim1, length_offset, list2_dim1,
1151             list2_offset, i__1, i__2;
1152
1153     /* Local variables */
1154     integer i__, j;
1155     integer new__, help, oops;
1156     doublereal help2, delta;
1157     doublereal costmin;
1158
1159 /* +-----------------------------------------------------------------------+ */
1160 /* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
1161 /* +-----------------------------------------------------------------------+ */
1162 /* +-----------------------------------------------------------------------+ */
1163 /* | JG 01/22/01 Added variable Ifeasiblef to keep track if feasibel point | */
1164 /* |             has been found.                                           | */
1165 /* | JG 01/23/01 Added variable Ierror to keep track of errors.            | */
1166 /* | JG 03/09/01 Added IInfeasible to keep track if an infeasible point has| */
1167 /* |             been found.                                               | */
1168 /* +-----------------------------------------------------------------------+ */
1169 /* JG 09/15/00 Added variable JONES (see above) */
1170 /* +-----------------------------------------------------------------------+ */
1171 /* | Variables to pass user defined data to the function to be optimized.  | */
1172 /* +-----------------------------------------------------------------------+ */
1173     /* Parameter adjustments */
1174     --point;
1175     f -= 3;
1176     ++anchor;
1177     --u;
1178     --l;
1179     --x;
1180     --w;
1181     list2_dim1 = *maxor;
1182     list2_offset = 1 + list2_dim1;
1183     list2 -= list2_offset;
1184     --arrayi;
1185     length_dim1 = *n;
1186     length_offset = 1 + length_dim1;
1187     length -= length_offset;
1188     c_dim1 = *maxor;
1189     c_offset = 1 + c_dim1;
1190     c__ -= c_offset;
1191
1192     /* Function Body */
1193     *minf = HUGE_VAL;
1194     costmin = *minf;
1195 /* JG 09/15/00 If Jones way of characterising rectangles is used, */
1196 /*             initialise thirds to reflect this. */
1197     if (jones == 0) {
1198         i__1 = *n - 1;
1199         for (j = 0; j <= i__1; ++j) {
1200             w[j + 1] = sqrt(*n - j + j / 9.) * .5;
1201 /* L5: */
1202         }
1203         help2 = 1.;
1204         i__1 = *maxdeep / *n;
1205         for (i__ = 1; i__ <= i__1; ++i__) {
1206             i__2 = *n - 1;
1207             for (j = 0; j <= i__2; ++j) {
1208                 levels[(i__ - 1) * *n + j] = w[j + 1] / help2;
1209 /* L8: */
1210             }
1211             help2 *= 3.;
1212 /* L10: */
1213         }
1214     } else {
1215 /* JG 09/15/00 Initialiase levels to contain 1/j */
1216         help2 = 3.;
1217         i__1 = *maxdeep;
1218         for (i__ = 1; i__ <= i__1; ++i__) {
1219             levels[i__] = 1. / help2;
1220             help2 *= 3.;
1221 /* L11: */
1222         }
1223         levels[0] = 1.;
1224     }
1225     help2 = 3.;
1226     i__1 = *maxdeep;
1227     for (i__ = 1; i__ <= i__1; ++i__) {
1228         thirds[i__] = 1. / help2;
1229         help2 *= 3.;
1230 /* L21: */
1231     }
1232     thirds[0] = 1.;
1233     i__1 = *n;
1234     for (i__ = 1; i__ <= i__1; ++i__) {
1235         c__[i__ + c_dim1] = .5;
1236         x[i__] = .5;
1237         length[i__ + length_dim1] = 0;
1238 /* L20: */
1239     }
1240     direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &f[3], &help, fcndata);
1241     if (force_stop && *force_stop) {
1242          *ierror = -102;
1243          return;
1244     }
1245     f[4] = (doublereal) help;
1246     *iinfeasible = help;
1247     *fmax = f[3];
1248 /* 09/25/00 Added this */
1249 /*      if (f(1,1) .ge. 1.E+6) then */
1250     if (f[4] > 0.) {
1251         f[3] = HUGE_VAL;
1252         *fmax = f[3];
1253         *ifeasiblef = 1;
1254     } else {
1255         *ifeasiblef = 0;
1256     }
1257 /* JG 09/25/00 Remove IF */
1258     *minf = f[3];
1259     costmin = f[3];
1260     *minpos = 1;
1261     *actdeep = 2;
1262     point[1] = 0;
1263     *free = 2;
1264     delta = thirds[1];
1265     if (nlopt_stop_time_(starttime, maxtime)) {
1266          *ierror = DIRECT_MAXTIME_EXCEEDED;
1267          return;
1268     }
1269     direct_dirget_i__(&length[length_offset], &c__1, &arrayi[1], maxi, n, maxfunc);
1270     new__ = *free;
1271     direct_dirsamplepoints_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &
1272             length[length_offset], logfile, &f[3], free, maxi, &
1273             point[1], &x[1], &l[1], minf, minpos, &u[1], n,
1274             maxfunc, maxdeep, &oops);
1275 /* +-----------------------------------------------------------------------+ */
1276 /* | JG 01/23/01 Added error checking.                                     | */
1277 /* +-----------------------------------------------------------------------+ */
1278     if (oops > 0) {
1279         *ierror = -4;
1280         return;
1281     }
1282 /* +-----------------------------------------------------------------------+ */
1283 /* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
1284 /* |             Added variable to keep track if feasible point was found. | */
1285 /* +-----------------------------------------------------------------------+ */
1286     direct_dirsamplef_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &length[
1287             length_offset], logfile, &f[3], free, maxi, &point[
1288             1], fcn, &x[1], &l[1], minf, minpos, &u[1], n, maxfunc,
1289             maxdeep, &oops, fmax, ifeasiblef, iinfeasible, fcndata,
1290             force_stop);
1291     if (force_stop && *force_stop) {
1292          *ierror = -102;
1293          return;
1294     }
1295     if (nlopt_stop_time_(starttime, maxtime)) {
1296          *ierror = DIRECT_MAXTIME_EXCEEDED;
1297          return;
1298     }
1299 /* +-----------------------------------------------------------------------+ */
1300 /* | JG 01/23/01 Added error checking.                                     | */
1301 /* +-----------------------------------------------------------------------+ */
1302     if (oops > 0) {
1303         *ierror = -5;
1304         return;
1305     }
1306     direct_dirdivide_(&new__, &c__0, &length[length_offset], &point[1], &arrayi[1], &
1307             c__1, &list2[list2_offset], &w[1], maxi, &f[3], maxfunc,
1308             maxdeep, n);
1309     direct_dirinsertlist_(&new__, &anchor[-1], &point[1], &f[3], maxi, &
1310             length[length_offset], maxfunc, maxdeep, n, &c__1, jones);
1311 } /* dirinit_ */
1312
1313 /* +-----------------------------------------------------------------------+ */
1314 /* |    SUBROUTINE DIRInitList                                             | */
1315 /* |    Initialise the list.                                               | */
1316 /* +-----------------------------------------------------------------------+ */
1317 /* Subroutine */ void direct_dirinitlist_(integer *anchor, integer *free, integer *
1318         point, doublereal *f, integer *maxfunc, const integer *maxdeep)
1319 {
1320     /* System generated locals */
1321     integer i__1;
1322
1323     /* Local variables */
1324     integer i__;
1325
1326 /*   f -- values of functions. */
1327 /*   anchor -- anchors of lists with deep i */
1328 /*   point -- lists */
1329 /*   free  -- first free position */
1330     /* Parameter adjustments */
1331     f -= 3;
1332     --point;
1333     ++anchor;
1334
1335     /* Function Body */
1336     i__1 = *maxdeep;
1337     for (i__ = -1; i__ <= i__1; ++i__) {
1338         anchor[i__] = 0;
1339 /* L10: */
1340     }
1341     i__1 = *maxfunc;
1342     for (i__ = 1; i__ <= i__1; ++i__) {
1343         f[(i__ << 1) + 1] = 0.;
1344         f[(i__ << 1) + 2] = 0.;
1345         point[i__] = i__ + 1;
1346 /*       point(i) = 0 */
1347 /* L20: */
1348     }
1349     point[*maxfunc] = 0;
1350     *free = 1;
1351 } /* dirinitlist_ */
1352
1353 /* +-----------------------------------------------------------------------+ */
1354 /* |                                                                       | */
1355 /* |                       SUBROUTINE DIRPREPRC                            | */
1356 /* |                                                                       | */
1357 /* | Subroutine DIRpreprc uses an afine mapping to map the hyper-box given | */
1358 /* | by the constraints on the variable x onto the n-dimensional unit cube.| */
1359 /* | This mapping is done using the following equation:                    | */
1360 /* |                                                                       | */
1361 /* |               x(i)=x(i)/(u(i)-l(i))-l(i)/(u(i)-l(i)).                 | */
1362 /* |                                                                       | */
1363 /* | DIRpreprc checks if the bounds l and u are well-defined. That is, if  | */
1364 /* |                                                                       | */
1365 /* |               l(i) < u(i) forevery i.                                 | */
1366 /* |                                                                       | */
1367 /* | On entry                                                              | */
1368 /* |                                                                       | */
1369 /* |          u -- A double-precision vector of length n. The vector       | */
1370 /* |               containing the upper bounds for the n independent       | */
1371 /* |               variables.                                              | */
1372 /* |                                                                       | */
1373 /* |          l -- A double-precision vector of length n. The vector       | */
1374 /* |               containing the lower bounds for the n independent       | */
1375 /* |               variables.                                              | */
1376 /* |                                                                       | */
1377 /* |          n -- An integer. The dimension of the problem.               | */
1378 /* |                                                                       | */
1379 /* | On return                                                             | */
1380 /* |                                                                       | */
1381 /* |        xs1 -- A double-precision vector of length n, used for scaling | */
1382 /* |               and unscaling the vector x.                             | */
1383 /* |                                                                       | */
1384 /* |        xs2 -- A double-precision vector of length n, used for scaling | */
1385 /* |               and unscaling the vector x.                             | */
1386 /* |                                                                       | */
1387 /* |                                                                       | */
1388 /* |       oops -- An integer. If an upper bound is less than a lower      | */
1389 /* |               bound or if the initial point is not in the             | */
1390 /* |               hyper-box oops is set to 1 and iffco terminates.        | */
1391 /* |                                                                       | */
1392 /* +-----------------------------------------------------------------------+ */
1393 /* Subroutine */ void direct_dirpreprc_(doublereal *u, doublereal *l, integer *n,
1394         doublereal *xs1, doublereal *xs2, integer *oops)
1395 {
1396     /* System generated locals */
1397     integer i__1;
1398
1399     /* Local variables */
1400     integer i__;
1401     doublereal help;
1402
1403     /* Parameter adjustments */
1404     --xs2;
1405     --xs1;
1406     --l;
1407     --u;
1408
1409     /* Function Body */
1410     *oops = 0;
1411     i__1 = *n;
1412     for (i__ = 1; i__ <= i__1; ++i__) {
1413 /* +-----------------------------------------------------------------------+ */
1414 /* | Check if the hyper-box is well-defined.                               | */
1415 /* +-----------------------------------------------------------------------+ */
1416         if (u[i__] <= l[i__]) {
1417             *oops = 1;
1418             return;
1419         }
1420 /* L20: */
1421     }
1422 /* +-----------------------------------------------------------------------+ */
1423 /* | Scale the initial iterate so that it is in the unit cube.             | */
1424 /* +-----------------------------------------------------------------------+ */
1425     i__1 = *n;
1426     for (i__ = 1; i__ <= i__1; ++i__) {
1427         help = u[i__] - l[i__];
1428         xs2[i__] = l[i__] / help;
1429         xs1[i__] = help;
1430 /* L50: */
1431     }
1432 } /* dirpreprc_ */
1433
1434 /* Subroutine */ void direct_dirheader_(FILE *logfile, integer *version,
1435         doublereal *x, integer *n, doublereal *eps, integer *maxf, integer *
1436         maxt, doublereal *l, doublereal *u, integer *algmethod, integer *
1437         maxfunc, const integer *maxdeep, doublereal *fglobal, doublereal *fglper,
1438         integer *ierror, doublereal *epsfix, integer *iepschange, doublereal *
1439         volper, doublereal *sigmaper)
1440 {
1441     /* System generated locals */
1442     integer i__1;
1443
1444     /* Local variables */
1445     integer imainver, i__, numerrors, isubsubver, ihelp, isubver;
1446
1447
1448 /* +-----------------------------------------------------------------------+ */
1449 /* | Variables to pass user defined data to the function to be optimized.  | */
1450 /* +-----------------------------------------------------------------------+ */
1451     /* Parameter adjustments */
1452     --u;
1453     --l;
1454     --x;
1455
1456     /* Function Body */
1457     if (logfile)
1458          fprintf(logfile, "------------------- Log file ------------------\n");
1459
1460     numerrors = 0;
1461     *ierror = 0;
1462     imainver = *version / 100;
1463     ihelp = *version - imainver * 100;
1464     isubver = ihelp / 10;
1465     ihelp -= isubver * 10;
1466     isubsubver = ihelp;
1467 /* +-----------------------------------------------------------------------+ */
1468 /* | JG 01/13/01 Added check for epsilon. If epsilon is smaller 0, we use  | */
1469 /* |             the update formula from Jones. We then set the flag       | */
1470 /* |             iepschange to 1, and store the absolute value of eps in   | */
1471 /* |             epsfix. epsilon is then changed after each iteration.     | */
1472 /* +-----------------------------------------------------------------------+ */
1473     if (*eps < 0.) {
1474         *iepschange = 1;
1475         *epsfix = -(*eps);
1476         *eps = -(*eps);
1477     } else {
1478         *iepschange = 0;
1479         *epsfix = 1e100;
1480     }
1481
1482 /* +-----------------------------------------------------------------------+ */
1483 /* | JG 07/16/01 Removed printout of contents in cdata(1).                 | */
1484 /* +-----------------------------------------------------------------------+ */
1485 /*      write(logfile,*) cdata(1) */
1486
1487     if (logfile) {
1488          fprintf(logfile, "DIRECT Version %d.%d.%d\n"
1489                  " Problem dimension n: %d\n"
1490                  " Eps value: %e\n"
1491                  " Maximum number of f-evaluations (maxf): %d\n"
1492                  " Maximum number of iterations (MaxT): %d\n"
1493                  " Value of f_global: %e\n"
1494                  " Global percentage wanted: %e\n"
1495                  " Volume percentage wanted: %e\n"
1496                  " Measure percentage wanted: %e\n",
1497                  imainver, isubver, isubsubver, *n, *eps, *maxf, *maxt,
1498                  *fglobal, *fglper, *volper, *sigmaper);
1499          fprintf(logfile, *iepschange == 1
1500                  ? "Epsilon is changed using the Jones formula.\n"
1501                  : "Epsilon is constant.\n");
1502          fprintf(logfile, *algmethod == 0
1503                  ? "Jones original DIRECT algorithm is used.\n"
1504                  : "Our modification of the DIRECT algorithm is used.\n");
1505     }
1506
1507     i__1 = *n;
1508     for (i__ = 1; i__ <= i__1; ++i__) {
1509         if (u[i__] <= l[i__]) {
1510             *ierror = -1;
1511             if (logfile)
1512                  fprintf(logfile, "WARNING: bounds on variable x%d: "
1513                          "%g <= xi <= %g\n", i__, l[i__], u[i__]);
1514             ++numerrors;
1515         } else {
1516             if (logfile)
1517                  fprintf(logfile, "Bounds on variable x%d: "
1518                          "%g <= xi <= %g\n", i__, l[i__], u[i__]);
1519         }
1520 /* L1010: */
1521     }
1522 /* +-----------------------------------------------------------------------+ */
1523 /* | If there are to many function evaluations or to many iteration, note  | */
1524 /* | this and set the error flag accordingly. Note: If more than one error | */
1525 /* | occurred, we give out an extra message.                               | */
1526 /* +-----------------------------------------------------------------------+ */
1527     if (*maxf + 20 > *maxfunc) {
1528         if (logfile)
1529              fprintf(logfile,
1530 "WARNING: The maximum number of function evaluations (%d) is higher than\n"
1531 "         the constant maxfunc (%d).  Increase maxfunc in subroutine DIRECT\n"
1532 "         or decrease the maximum number of function evaluations.\n",
1533                      *maxf, *maxfunc);
1534         ++numerrors;
1535         *ierror = -2;
1536     }
1537     if (*ierror < 0) {
1538         if (logfile) fprintf(logfile, "----------------------------------\n");
1539         if (numerrors == 1) {
1540              if (logfile)
1541                   fprintf(logfile, "WARNING: One error in the input!\n");
1542         } else {
1543              if (logfile)
1544                   fprintf(logfile, "WARNING: %d errors in the input!\n",
1545                           numerrors);
1546         }
1547     }
1548     if (logfile) fprintf(logfile, "----------------------------------\n");
1549     if (*ierror >= 0) {
1550          if (logfile)
1551               fprintf(logfile, "Iteration # of f-eval. minf\n");
1552     }
1553 /* L10005: */
1554 } /* dirheader_ */
1555
1556 /* Subroutine */ void direct_dirsummary_(FILE *logfile, doublereal *x, doublereal *
1557         l, doublereal *u, integer *n, doublereal *minf, doublereal *fglobal,
1558         integer *numfunc, integer *ierror)
1559 {
1560     /* Local variables */
1561     integer i__;
1562
1563     /* Parameter adjustments */
1564     --u;
1565     --l;
1566     --x;
1567
1568     /* Function Body */
1569     if (logfile) {
1570          fprintf(logfile, "-----------------------Summary------------------\n"
1571                  "Final function value: %g\n"
1572                  "Number of function evaluations: %d\n", *minf, *numfunc);
1573          if (*fglobal > -1e99)
1574               fprintf(logfile, "Final function value is within %g%% of global optimum\n", 100*(*minf - *fglobal) / MAX(1.0, fabs(*fglobal)));
1575          fprintf(logfile, "Index, final solution, x(i)-l(i), u(i)-x(i)\n");
1576          for (i__ = 1; i__ <= *n; ++i__)
1577               fprintf(logfile, "%d, %g, %g, %g\n", i__, x[i__],
1578                       x[i__]-l[i__], u[i__] - x[i__]);
1579          fprintf(logfile, "-----------------------------------------------\n");
1580
1581     }
1582 } /* dirsummary_ */