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