subplex/subplex.c

   1 /*
   2 Downloaded from http://www.netlib.org/opt/subplex.tgz
   3
   4 README file for SUBPLEX
   5
   6 NAME
   7      subplex - subspace-searching simplex method for unconstrained
   8      optimization
   9
  10 DESCRIPTION
  11      Subplex is a subspace-searching simplex method for the
  12      unconstrained optimization of general multivariate functions.
  13      Like the Nelder-Mead simplex method it generalizes, the subplex
  14      method is well suited for optimizing noisy objective functions.
  15      The number of function evaluations required for convergence
  16      typically increases only linearly with the problem size, so for
  17      most applications the subplex method is much more efficient than
  18      the simplex method.
  19
  20 INSTALLATION
  21      To build subplex on UNIX systems, edit the Makefile as necessary
  22      and type:
  23
  24              make
  25
  26      This will create a linkable library named subplex.a and a
  27      demonstration executable named demo.
  28
  29 EXAMPLE
  30      To run subplex on a simple objective function type:
  31
  32              demo < demo.in
  33
  34      To run subplex on other problems, edit a copy of the sample driver
  35      demo.f as necessary.
  36
  37 AUTHOR
  38      Tom Rowan
  39      Oak Ridge National Laboratory
  40      Mathematical Sciences Section
  41      P.O. Box 2008, Bldg. 6012
  42      Oak Ridge, TN 37831-6367
  43
  44      Phone:   (423) 574-3131
  45      Fax  :   (423) 574-0680
  46      Email:   na.rowan@na-net.ornl.gov
  47
  48 REFERENCE
  49      T. Rowan, "Functional Stability Analysis of Numerical Algorithms",
  50      Ph.D. thesis, Department of Computer Sciences, University of Texas
  51      at Austin, 1990.
  52
  53 COMMENTS
  54      Please send comments, suggestions, or bug reports to
  55      na.rowan@na-net.ornl.gov.
  56  */
  57
  58 #include <stdio.h>
  59 #include <stdlib.h>
  60 #include <math.h>
  61 #include <string.h>
  62
  63 #include "subplex.h"
  64
  65 typedef double doublereal;
  66 typedef int logical;
  67 typedef int integer;
  68
  69 #define TRUE_ 1
  70 #define FALSE_ 0
  71
  72 typedef subplex_func D_fp;
  73
  74 #define max(a,b) ((a) > (b) ? (a) : (b))
  75 #define min(a,b) ((a) < (b) ? (a) : (b))
  76 #define abs(x) fabs(x)
  77
  78 /****************************************************************************/
  79 /****************************************************************************/
  80
  81 /* dasum.f -- translated by f2c (version 19991025).
  82    You must link the resulting object file with the libraries:
  83         -lf2c -lm   (in that order)
  84 */
  85
  86 static doublereal dasum_(integer *n, doublereal *dx, integer *incx)
  87 {
  88     /* System generated locals */
  89     integer i__1;
  90     doublereal ret_val, d__1, d__2, d__3, d__4, d__5, d__6;
  91
  92     /* Local variables */
  93     integer i__, m;
  94     doublereal dtemp;
  95     integer ix, mp1;
  96
  97
  98 /*     takes the sum of the absolute values. */
  99 /*     uses unrolled loops for increment equal to one. */
 100 /*     jack dongarra, linpack, 3/11/78. */
 101 /*     modified to correct problem with negative increment, 8/21/90. */
 102
 103
 104     /* Parameter adjustments */
 105     --dx;
 106
 107     /* Function Body */
 108     ret_val = 0.;
 109     dtemp = 0.;
 110     if (*n <= 0) {
 111         return ret_val;
 112     }
 113     if (*incx == 1) {
 114         goto L20;
 115     }
 116
 117 /*        code for increment not equal to 1 */
 118
 119     ix = 1;
 120     if (*incx < 0) {
 121         ix = (-(*n) + 1) * *incx + 1;
 122     }
 123     i__1 = *n;
 124     for (i__ = 1; i__ <= i__1; ++i__) {
 125         dtemp += (d__1 = dx[ix], abs(d__1));
 126         ix += *incx;
 127 /* L10: */
 128     }
 129     ret_val = dtemp;
 130     return ret_val;
 131
 132 /*        code for increment equal to 1 */
 133
 134
 135 /*        clean-up loop */
 136
 137 L20:
 138     m = *n % 6;
 139     if (m == 0) {
 140         goto L40;
 141     }
 142     i__1 = m;
 143     for (i__ = 1; i__ <= i__1; ++i__) {
 144         dtemp += (d__1 = dx[i__], abs(d__1));
 145 /* L30: */
 146     }
 147     if (*n < 6) {
 148         goto L60;
 149     }
 150 L40:
 151     mp1 = m + 1;
 152     i__1 = *n;
 153     for (i__ = mp1; i__ <= i__1; i__ += 6) {
 154         dtemp = dtemp + (d__1 = dx[i__], abs(d__1)) + (d__2 = dx[i__ + 1],
 155                 abs(d__2)) + (d__3 = dx[i__ + 2], abs(d__3)) + (d__4 = dx[i__
 156                 + 3], abs(d__4)) + (d__5 = dx[i__ + 4], abs(d__5)) + (d__6 =
 157                 dx[i__ + 5], abs(d__6));
 158 /* L50: */
 159     }
 160 L60:
 161     ret_val = dtemp;
 162     return ret_val;
 163 } /* dasum_ */
 164
 165 /* daxpy.f -- translated by f2c (version 19991025).
 166    You must link the resulting object file with the libraries:
 167         -lf2c -lm   (in that order)
 168 */
 169
 170 static int daxpy_(integer *n, doublereal *da, doublereal *dx,
 171         integer *incx, doublereal *dy, integer *incy)
 172 {
 173     /* System generated locals */
 174     integer i__1;
 175
 176     /* Local variables */
 177     integer i__, m, ix, iy, mp1;
 178
 179
 180 /*     constant times a vector plus a vector. */
 181 /*     uses unrolled loops for increments equal to one. */
 182 /*     jack dongarra, linpack, 3/11/78. */
 183
 184
 185     /* Parameter adjustments */
 186     --dy;
 187     --dx;
 188
 189     /* Function Body */
 190     if (*n <= 0) {
 191         return 0;
 192     }
 193     if (*da == 0.) {
 194         return 0;
 195     }
 196     if (*incx == 1 && *incy == 1) {
 197         goto L20;
 198     }
 199
 200 /*        code for unequal increments or equal increments */
 201 /*          not equal to 1 */
 202
 203     ix = 1;
 204     iy = 1;
 205     if (*incx < 0) {
 206         ix = (-(*n) + 1) * *incx + 1;
 207     }
 208     if (*incy < 0) {
 209         iy = (-(*n) + 1) * *incy + 1;
 210     }
 211     i__1 = *n;
 212     for (i__ = 1; i__ <= i__1; ++i__) {
 213         dy[iy] += *da * dx[ix];
 214         ix += *incx;
 215         iy += *incy;
 216 /* L10: */
 217     }
 218     return 0;
 219
 220 /*        code for both increments equal to 1 */
 221
 222
 223 /*        clean-up loop */
 224
 225 L20:
 226     m = *n % 4;
 227     if (m == 0) {
 228         goto L40;
 229     }
 230     i__1 = m;
 231     for (i__ = 1; i__ <= i__1; ++i__) {
 232         dy[i__] += *da * dx[i__];
 233 /* L30: */
 234     }
 235     if (*n < 4) {
 236         return 0;
 237     }
 238 L40:
 239     mp1 = m + 1;
 240     i__1 = *n;
 241     for (i__ = mp1; i__ <= i__1; i__ += 4) {
 242         dy[i__] += *da * dx[i__];
 243         dy[i__ + 1] += *da * dx[i__ + 1];
 244         dy[i__ + 2] += *da * dx[i__ + 2];
 245         dy[i__ + 3] += *da * dx[i__ + 3];
 246 /* L50: */
 247     }
 248     return 0;
 249 } /* daxpy_ */
 250
 251 /* dcopy.f -- translated by f2c (version 19991025).
 252    You must link the resulting object file with the libraries:
 253         -lf2c -lm   (in that order)
 254 */
 255
 256 static int dcopy_(integer *n, const doublereal *dx, integer *incx,
 257         doublereal *dy, integer *incy)
 258 {
 259     /* System generated locals */
 260     integer i__1;
 261
 262     /* Local variables */
 263     integer i__, m, ix, iy, mp1;
 264
 265
 266 /*     copies a vector, x, to a vector, y. */
 267 /*     uses unrolled loops for increments equal to one. */
 268 /*     jack dongarra, linpack, 3/11/78. */
 269
 270
 271     /* Parameter adjustments */
 272     --dy;
 273     --dx;
 274
 275     /* Function Body */
 276     if (*n <= 0) {
 277         return 0;
 278     }
 279     if (*incx == 1 && *incy == 1) {
 280         goto L20;
 281     }
 282
 283 /*        code for unequal increments or equal increments */
 284 /*          not equal to 1 */
 285
 286     ix = 1;
 287     iy = 1;
 288     if (*incx < 0) {
 289         ix = (-(*n) + 1) * *incx + 1;
 290     }
 291     if (*incy < 0) {
 292         iy = (-(*n) + 1) * *incy + 1;
 293     }
 294     i__1 = *n;
 295     for (i__ = 1; i__ <= i__1; ++i__) {
 296         dy[iy] = dx[ix];
 297         ix += *incx;
 298         iy += *incy;
 299 /* L10: */
 300     }
 301     return 0;
 302
 303 /*        code for both increments equal to 1 */
 304
 305
 306 /*        clean-up loop */
 307
 308 L20:
 309     m = *n % 7;
 310     if (m == 0) {
 311         goto L40;
 312     }
 313     i__1 = m;
 314     for (i__ = 1; i__ <= i__1; ++i__) {
 315         dy[i__] = dx[i__];
 316 /* L30: */
 317     }
 318     if (*n < 7) {
 319         return 0;
 320     }
 321 L40:
 322     mp1 = m + 1;
 323     i__1 = *n;
 324     for (i__ = mp1; i__ <= i__1; i__ += 7) {
 325         dy[i__] = dx[i__];
 326         dy[i__ + 1] = dx[i__ + 1];
 327         dy[i__ + 2] = dx[i__ + 2];
 328         dy[i__ + 3] = dx[i__ + 3];
 329         dy[i__ + 4] = dx[i__ + 4];
 330         dy[i__ + 5] = dx[i__ + 5];
 331         dy[i__ + 6] = dx[i__ + 6];
 332 /* L50: */
 333     }
 334     return 0;
 335 } /* dcopy_ */
 336
 337 /* dscal.f -- translated by f2c (version 19991025).
 338    You must link the resulting object file with the libraries:
 339         -lf2c -lm   (in that order)
 340 */
 341
 342 static int dscal_(integer *n, doublereal *da, doublereal *dx,
 343         integer *incx)
 344 {
 345     /* System generated locals */
 346     integer i__1;
 347
 348     /* Local variables */
 349     integer i__, m, ix, mp1;
 350
 351
 352 /*     scales a vector by a constant. */
 353 /*     uses unrolled loops for increment equal to one. */
 354 /*     jack dongarra, linpack, 3/11/78. */
 355 /*     modified to correct problem with negative increment, 8/21/90. */
 356
 357
 358     /* Parameter adjustments */
 359     --dx;
 360
 361     /* Function Body */
 362     if (*n <= 0) {
 363         return 0;
 364     }
 365     if (*incx == 1) {
 366         goto L20;
 367     }
 368
 369 /*        code for increment not equal to 1 */
 370
 371     ix = 1;
 372     if (*incx < 0) {
 373         ix = (-(*n) + 1) * *incx + 1;
 374     }
 375     i__1 = *n;
 376     for (i__ = 1; i__ <= i__1; ++i__) {
 377         dx[ix] = *da * dx[ix];
 378         ix += *incx;
 379 /* L10: */
 380     }
 381     return 0;
 382
 383 /*        code for increment equal to 1 */
 384
 385
 386 /*        clean-up loop */
 387
 388 L20:
 389     m = *n % 5;
 390     if (m == 0) {
 391         goto L40;
 392     }
 393     i__1 = m;
 394     for (i__ = 1; i__ <= i__1; ++i__) {
 395         dx[i__] = *da * dx[i__];
 396 /* L30: */
 397     }
 398     if (*n < 5) {
 399         return 0;
 400     }
 401 L40:
 402     mp1 = m + 1;
 403     i__1 = *n;
 404     for (i__ = mp1; i__ <= i__1; i__ += 5) {
 405         dx[i__] = *da * dx[i__];
 406         dx[i__ + 1] = *da * dx[i__ + 1];
 407         dx[i__ + 2] = *da * dx[i__ + 2];
 408         dx[i__ + 3] = *da * dx[i__ + 3];
 409         dx[i__ + 4] = *da * dx[i__ + 4];
 410 /* L50: */
 411     }
 412     return 0;
 413 } /* dscal_ */
 414
 415 /* dist.f -- translated by f2c (version 19991025).
 416    You must link the resulting object file with the libraries:
 417         -lf2c -lm   (in that order)
 418 */
 419
 420 static doublereal dist_(integer *n, doublereal *x, doublereal *y)
 421 {
 422     /* System generated locals */
 423     integer i__1;
 424     doublereal ret_val, d__1;
 425
 426     /* Local variables */
 427     integer i__;
 428     doublereal scale, absxmy, sum;
 429
 430
 431
 432 /*                                         Coded by Tom Rowan */
 433 /*                            Department of Computer Sciences */
 434 /*                              University of Texas at Austin */
 435
 436 /* dist calculates the distance between the points x,y. */
 437
 438 /* input */
 439
 440 /*   n      - number of components */
 441
 442 /*   x      - point in n-space */
 443
 444 /*   y      - point in n-space */
 445
 446 /* local variables */
 447
 448
 449 /* subroutines and functions */
 450
 451 /*   fortran */
 452
 453 /* ----------------------------------------------------------- */
 454
 455     /* Parameter adjustments */
 456     --y;
 457     --x;
 458
 459     /* Function Body */
 460     absxmy = (d__1 = x[1] - y[1], abs(d__1));
 461     if (absxmy <= 1.) {
 462         sum = absxmy * absxmy;
 463         scale = 1.;
 464     } else {
 465         sum = 1.;
 466         scale = absxmy;
 467     }
 468     i__1 = *n;
 469     for (i__ = 2; i__ <= i__1; ++i__) {
 470         absxmy = (d__1 = x[i__] - y[i__], abs(d__1));
 471         if (absxmy <= scale) {
 472 /* Computing 2nd power */
 473             d__1 = absxmy / scale;
 474             sum += d__1 * d__1;
 475         } else {
 476 /* Computing 2nd power */
 477             d__1 = scale / absxmy;
 478             sum = sum * (d__1 * d__1) + 1.;
 479             scale = absxmy;
 480         }
 481 /* L10: */
 482     }
 483     ret_val = scale * sqrt(sum);
 484     return ret_val;
 485 } /* dist_ */
 486
 487 /* calcc.f -- translated by f2c (version 19991025).
 488    You must link the resulting object file with the libraries:
 489         -lf2c -lm   (in that order)
 490 */
 491
 492 /* Table of constant values */
 493
 494 static doublereal c_b3 = 0.;
 495 static integer c__0 = 0;
 496 static integer c__1 = 1;
 497 static doublereal c_b7 = 1.;
 498
 499 static int calcc_(integer *ns, doublereal *s, integer *ih, integer *
 500         inew, logical *updatc, doublereal *c__)
 501 {
 502     /* System generated locals */
 503     integer s_dim1, s_offset, i__1;
 504     doublereal d__1;
 505
 506     /* Local variables */
 507     integer i__, j;
 508
 509 /*                                         Coded by Tom Rowan */
 510 /*                            Department of Computer Sciences */
 511 /*                              University of Texas at Austin */
 512
 513 /* calcc calculates the centroid of the simplex without the */
 514 /* vertex with highest function value. */
 515
 516 /* input */
 517
 518 /*   ns     - subspace dimension */
 519
 520 /*   s      - double precision work space of dimension .ge. */
 521 /*            ns*(ns+3) used to store simplex */
 522
 523 /*   ih     - index to vertex with highest function value */
 524
 525 /*   inew   - index to new point */
 526
 527 /*   updatc - logical switch */
 528 /*            = .true.  : update centroid */
 529 /*            = .false. : calculate centroid from scratch */
 530
 531 /*   c      - centroid of the simplex without vertex with */
 532 /*            highest function value */
 533
 534 /* output */
 535
 536 /*   c      - new centroid */
 537
 538 /* local variables */
 539
 540
 541 /* subroutines and functions */
 542
 543 /*   blas */
 544
 545 /* ----------------------------------------------------------- */
 546
 547     /* Parameter adjustments */
 548     --c__;
 549     s_dim1 = *ns;
 550     s_offset = 1 + s_dim1 * 1;
 551     s -= s_offset;
 552
 553     /* Function Body */
 554     if (*updatc) {
 555         if (*ih == *inew) {
 556             return 0;
 557         }
 558         i__1 = *ns;
 559         for (i__ = 1; i__ <= i__1; ++i__) {
 560             c__[i__] += (s[i__ + *inew * s_dim1] - s[i__ + *ih * s_dim1]) / *
 561                     ns;
 562 /* L10: */
 563         }
 564     } else {
 565         dcopy_(ns, &c_b3, &c__0, &c__[1], &c__1);
 566         i__1 = *ns + 1;
 567         for (j = 1; j <= i__1; ++j) {
 568             if (j != *ih) {
 569                 daxpy_(ns, &c_b7, &s[j * s_dim1 + 1], &c__1, &c__[1], &c__1);
 570             }
 571 /* L20: */
 572         }
 573         d__1 = 1. / *ns;
 574         dscal_(ns, &d__1, &c__[1], &c__1);
 575     }
 576     return 0;
 577 } /* calcc_ */
 578
 579 /* order.f -- translated by f2c (version 19991025).
 580    You must link the resulting object file with the libraries:
 581         -lf2c -lm   (in that order)
 582 */
 583
 584 static int order_(integer *npts, doublereal *fs, integer *il,
 585         integer *is, integer *ih)
 586 {
 587     /* System generated locals */
 588     integer i__1;
 589
 590     /* Local variables */
 591     integer i__, j, il0;
 592
 593
 594
 595 /*                                         Coded by Tom Rowan */
 596 /*                            Department of Computer Sciences */
 597 /*                              University of Texas at Austin */
 598
 599 /* order determines the indices of the vertices with the */
 600 /* lowest, second highest, and highest function values. */
 601
 602 /* input */
 603
 604 /*   npts   - number of points in simplex */
 605
 606 /*   fs     - double precision vector of function values of */
 607 /*            simplex */
 608
 609 /*   il     - index to vertex with lowest function value */
 610
 611 /* output */
 612
 613 /*   il     - new index to vertex with lowest function value */
 614
 615 /*   is     - new index to vertex with second highest */
 616 /*            function value */
 617
 618 /*   ih     - new index to vertex with highest function value */
 619
 620 /* local variables */
 621
 622
 623 /* subroutines and functions */
 624
 625 /*   fortran */
 626
 627 /* ----------------------------------------------------------- */
 628
 629     /* Parameter adjustments */
 630     --fs;
 631
 632     /* Function Body */
 633     il0 = *il;
 634     j = il0 % *npts + 1;
 635     if (fs[j] >= fs[*il]) {
 636         *ih = j;
 637         *is = il0;
 638     } else {
 639         *ih = il0;
 640         *is = j;
 641         *il = j;
 642     }
 643     i__1 = il0 + *npts - 2;
 644     for (i__ = il0 + 1; i__ <= i__1; ++i__) {
 645         j = i__ % *npts + 1;
 646         if (fs[j] >= fs[*ih]) {
 647             *is = *ih;
 648             *ih = j;
 649         } else if (fs[j] > fs[*is]) {
 650             *is = j;
 651         } else if (fs[j] < fs[*il]) {
 652             *il = j;
 653         }
 654 /* L10: */
 655     }
 656     return 0;
 657 } /* order_ */
 658
 659 /* partx.f -- translated by f2c (version 19991025).
 660    You must link the resulting object file with the libraries:
 661         -lf2c -lm   (in that order)
 662 */
 663
 664 /* Common Block Declarations */
 665
 666 static struct {
 667     doublereal alpha, beta, gamma, delta, psi, omega;
 668     integer nsmin, nsmax, irepl, ifxsw;
 669     doublereal bonus, fstop;
 670     integer nfstop, nfxe;
 671     doublereal fxstat[4], ftest;
 672     logical minf, initx, newx;
 673 } usubc_;
 674
 675 #define usubc_1 usubc_
 676
 677 static int partx_(integer *n, integer *ip, doublereal *absdx,
 678         integer *nsubs, integer *nsvals)
 679 {
 680     /* System generated locals */
 681     integer i__1;
 682
 683     /* Local variables */
 684     static integer i__, nleft, nused;
 685     static doublereal as1max, gapmax, asleft, as1, as2;
 686     static integer ns1, ns2;
 687     static doublereal gap;
 688
 689
 690
 691 /*                                         Coded by Tom Rowan */
 692 /*                            Department of Computer Sciences */
 693 /*                              University of Texas at Austin */
 694
 695 /* partx partitions the vector x by grouping components of */
 696 /* similar magnitude of change. */
 697
 698 /* input */
 699
 700 /*   n      - number of components (problem dimension) */
 701
 702 /*   ip     - permutation vector */
 703
 704 /*   absdx  - vector of magnitude of change in x */
 705
 706 /*   nsvals - integer array dimensioned .ge. int(n/nsmin) */
 707
 708 /* output */
 709
 710 /*   nsubs  - number of subspaces */
 711
 712 /*   nsvals - integer array of subspace dimensions */
 713
 714 /* common */
 715
 716
 717
 718 /* local variables */
 719
 720
 721
 722 /* subroutines and functions */
 723
 724 /*   fortran */
 725
 726 /* ----------------------------------------------------------- */
 727
 728     /* Parameter adjustments */
 729     --absdx;
 730     --ip;
 731     --nsvals;
 732
 733     /* Function Body */
 734     *nsubs = 0;
 735     nused = 0;
 736     nleft = *n;
 737     asleft = absdx[1];
 738     i__1 = *n;
 739     for (i__ = 2; i__ <= i__1; ++i__) {
 740         asleft += absdx[i__];
 741 /* L10: */
 742     }
 743 L20:
 744     if (nused < *n) {
 745         ++(*nsubs);
 746         as1 = 0.;
 747         i__1 = usubc_1.nsmin - 1;
 748         for (i__ = 1; i__ <= i__1; ++i__) {
 749             as1 += absdx[ip[nused + i__]];
 750 /* L30: */
 751         }
 752         gapmax = -1.;
 753         i__1 = min(usubc_1.nsmax,nleft);
 754         for (ns1 = usubc_1.nsmin; ns1 <= i__1; ++ns1) {
 755             as1 += absdx[ip[nused + ns1]];
 756             ns2 = nleft - ns1;
 757             if (ns2 > 0) {
 758                 if (ns2 >= ((ns2 - 1) / usubc_1.nsmax + 1) * usubc_1.nsmin) {
 759                     as2 = asleft - as1;
 760                     gap = as1 / ns1 - as2 / ns2;
 761                     if (gap > gapmax) {
 762                         gapmax = gap;
 763                         nsvals[*nsubs] = ns1;
 764                         as1max = as1;
 765                     }
 766                 }
 767             } else {
 768                 if (as1 / ns1 > gapmax) {
 769                     nsvals[*nsubs] = ns1;
 770                     return 0;
 771                 }
 772             }
 773 /* L40: */
 774         }
 775         nused += nsvals[*nsubs];
 776         nleft = *n - nused;
 777         asleft -= as1max;
 778         goto L20;
 779     }
 780     return 0;
 781 } /* partx_ */
 782
 783 /* sortd.f -- translated by f2c (version 19991025).
 784    You must link the resulting object file with the libraries:
 785         -lf2c -lm   (in that order)
 786 */
 787
 788 static int sortd_(integer *n, doublereal *xkey, integer *ix)
 789 {
 790     /* System generated locals */
 791     integer i__1;
 792
 793     /* Local variables */
 794     integer ixip1, i__, ilast, iswap, ifirst, ixi;
 795
 796
 797
 798 /*                                         Coded by Tom Rowan */
 799 /*                            Department of Computer Sciences */
 800 /*                              University of Texas at Austin */
 801
 802 /* sortd uses the shakersort method to sort an array of keys */
 803 /* in decreasing order. The sort is performed implicitly by */
 804 /* modifying a vector of indices. */
 805
 806 /* For nearly sorted arrays, sortd requires O(n) comparisons. */
 807 /* for completely unsorted arrays, sortd requires O(n**2) */
 808 /* comparisons and will be inefficient unless n is small. */
 809
 810 /* input */
 811
 812 /*   n      - number of components */
 813
 814 /*   xkey   - double precision vector of keys */
 815
 816 /*   ix     - integer vector of indices */
 817
 818 /* output */
 819
 820 /*   ix     - indices satisfy xkey(ix(i)) .ge. xkey(ix(i+1)) */
 821 /*            for i = 1,...,n-1 */
 822
 823 /* local variables */
 824
 825
 826 /* ----------------------------------------------------------- */
 827
 828     /* Parameter adjustments */
 829     --ix;
 830     --xkey;
 831
 832     /* Function Body */
 833     ifirst = 1;
 834     iswap = 1;
 835     ilast = *n - 1;
 836 L10:
 837     if (ifirst <= ilast) {
 838         i__1 = ilast;
 839         for (i__ = ifirst; i__ <= i__1; ++i__) {
 840             ixi = ix[i__];
 841             ixip1 = ix[i__ + 1];
 842             if (xkey[ixi] < xkey[ixip1]) {
 843                 ix[i__] = ixip1;
 844                 ix[i__ + 1] = ixi;
 845                 iswap = i__;
 846             }
 847 /* L20: */
 848         }
 849         ilast = iswap - 1;
 850         i__1 = ifirst;
 851         for (i__ = ilast; i__ >= i__1; --i__) {
 852             ixi = ix[i__];
 853             ixip1 = ix[i__ + 1];
 854             if (xkey[ixi] < xkey[ixip1]) {
 855                 ix[i__] = ixip1;
 856                 ix[i__ + 1] = ixi;
 857                 iswap = i__;
 858             }
 859 /* L30: */
 860         }
 861         ifirst = iswap + 1;
 862         goto L10;
 863     }
 864     return 0;
 865 } /* sortd_ */
 866
 867 /* newpt.f -- translated by f2c (version 19991025).
 868    You must link the resulting object file with the libraries:
 869         -lf2c -lm   (in that order)
 870 */
 871
 872 static int newpt_(integer *ns, doublereal *coef, doublereal *xbase,
 873         doublereal *xold, logical *new__, doublereal *xnew, logical *small)
 874 {
 875     /* System generated locals */
 876     integer i__1;
 877
 878     /* Local variables */
 879     integer i__;
 880     logical eqold;
 881     doublereal xoldi;
 882     logical eqbase;
 883
 884
 885
 886 /*                                         Coded by Tom Rowan */
 887 /*                            Department of Computer Sciences */
 888 /*                              University of Texas at Austin */
 889
 890 /* newpt performs reflections, expansions, contractions, and */
 891 /* shrinkages (massive contractions) by computing: */
 892
 893 /* xbase + coef * (xbase - xold) */
 894
 895 /* The result is stored in xnew if new .eq. .true., */
 896 /* in xold otherwise. */
 897
 898 /* use :  coef .gt. 0 to reflect */
 899 /*        coef .lt. 0 to expand, contract, or shrink */
 900
 901 /* input */
 902
 903 /*   ns     - number of components (subspace dimension) */
 904
 905 /*   coef   - one of four simplex method coefficients */
 906
 907 /*   xbase  - double precision ns-vector representing base */
 908 /*            point */
 909
 910 /*   xold   - double precision ns-vector representing old */
 911 /*            point */
 912
 913 /*   new    - logical switch */
 914 /*            = .true.  : store result in xnew */
 915 /*            = .false. : store result in xold, xnew is not */
 916 /*                        referenced */
 917
 918 /* output */
 919
 920 /*   xold   - unchanged if new .eq. .true., contains new */
 921 /*            point otherwise */
 922
 923 /*   xnew   - double precision ns-vector representing new */
 924 /*            point if  new .eq. .true., not referenced */
 925 /*            otherwise */
 926
 927 /*   small  - logical flag */
 928 /*            = .true.  : coincident points */
 929 /*            = .false. : otherwise */
 930
 931 /* local variables */
 932
 933
 934 /* subroutines and functions */
 935
 936 /*   fortran */
 937
 938 /* ----------------------------------------------------------- */
 939
 940     /* Parameter adjustments */
 941     --xold;
 942     --xbase;
 943     --xnew;
 944
 945     /* Function Body */
 946     eqbase = TRUE_;
 947     eqold = TRUE_;
 948     if (*new__) {
 949         i__1 = *ns;
 950         for (i__ = 1; i__ <= i__1; ++i__) {
 951             xnew[i__] = xbase[i__] + *coef * (xbase[i__] - xold[i__]);
 952             eqbase = eqbase && xnew[i__] == xbase[i__];
 953             eqold = eqold && xnew[i__] == xold[i__];
 954 /* L10: */
 955         }
 956     } else {
 957         i__1 = *ns;
 958         for (i__ = 1; i__ <= i__1; ++i__) {
 959             xoldi = xold[i__];
 960             xold[i__] = xbase[i__] + *coef * (xbase[i__] - xold[i__]);
 961             eqbase = eqbase && xold[i__] == xbase[i__];
 962             eqold = eqold && xold[i__] == xoldi;
 963 /* L20: */
 964         }
 965     }
 966     *small = eqbase || eqold;
 967     return 0;
 968 } /* newpt_ */
 969
 970 /* start.f -- translated by f2c (version 19991025).
 971    You must link the resulting object file with the libraries:
 972         -lf2c -lm   (in that order)
 973 */
 974
 975 static int start_(integer *n, doublereal *x, doublereal *step,
 976         integer *ns, integer *ips, doublereal *s, logical *small)
 977 {
 978     /* System generated locals */
 979     integer s_dim1, s_offset, i__1;
 980
 981     /* Local variables */
 982     integer i__, j;
 983
 984
 985 /*                                         Coded by Tom Rowan */
 986 /*                            Department of Computer Sciences */
 987 /*                              University of Texas at Austin */
 988
 989 /* start creates the initial simplex for simplx minimization. */
 990
 991 /* input */
 992
 993 /*   n      - problem dimension */
 994
 995 /*   x      - current best point */
 996
 997 /*   step   - stepsizes for corresponding components of x */
 998
 999 /*   ns     - subspace dimension */
1000
1001 /*   ips    - permutation vector */
1002
1003
1004 /* output */
1005
1006 /*   s      - first ns+1 columns contain initial simplex */
1007
1008 /*   small  - logical flag */
1009 /*            = .true.  : coincident points */
1010 /*            = .false. : otherwise */
1011
1012 /* local variables */
1013
1014
1015 /* subroutines and functions */
1016
1017 /*   blas */
1018 /*   fortran */
1019
1020 /* ----------------------------------------------------------- */
1021
1022     /* Parameter adjustments */
1023     --ips;
1024     --step;
1025     --x;
1026     s_dim1 = *ns;
1027     s_offset = 1 + s_dim1 * 1;
1028     s -= s_offset;
1029
1030     /* Function Body */
1031     i__1 = *ns;
1032     for (i__ = 1; i__ <= i__1; ++i__) {
1033         s[i__ + s_dim1] = x[ips[i__]];
1034 /* L10: */
1035     }
1036     i__1 = *ns + 1;
1037     for (j = 2; j <= i__1; ++j) {
1038         dcopy_(ns, &s[s_dim1 + 1], &c__1, &s[j * s_dim1 + 1], &c__1);
1039         s[j - 1 + j * s_dim1] = s[j - 1 + s_dim1] + step[ips[j - 1]];
1040 /* L20: */
1041     }
1042
1043 /* check for coincident points */
1044
1045     i__1 = *ns + 1;
1046     for (j = 2; j <= i__1; ++j) {
1047         if (s[j - 1 + j * s_dim1] == s[j - 1 + s_dim1]) {
1048             goto L40;
1049         }
1050 /* L30: */
1051     }
1052     *small = FALSE_;
1053     return 0;
1054
1055 /* coincident points */
1056
1057 L40:
1058     *small = TRUE_;
1059     return 0;
1060 } /* start_ */
1061
1062 /* fstats.f -- translated by f2c (version 19991025).
1063    You must link the resulting object file with the libraries:
1064         -lf2c -lm   (in that order)
1065 */
1066
1067 static int fstats_(doublereal *fx, integer *ifxwt, logical *reset)
1068 {
1069     /* System generated locals */
1070     doublereal d__1, d__2, d__3;
1071
1072     /* Local variables */
1073     static doublereal fscale;
1074     static integer nsv;
1075     static doublereal f1sv;
1076
1077
1078
1079 /*                                         Coded by Tom Rowan */
1080 /*                            Department of Computer Sciences */
1081 /*                              University of Texas at Austin */
1082
1083 /* fstats modifies the common /usubc/ variables nfxe,fxstat. */
1084
1085 /* input */
1086
1087 /*   fx     - most recent evaluation of f at best x */
1088
1089 /*   ifxwt  - integer weight for fx */
1090
1091 /*   reset  - logical switch */
1092 /*            = .true.  : initialize nfxe,fxstat */
1093 /*            = .false. : update nfxe,fxstat */
1094
1095 /* common */
1096
1097
1098
1099 /* local variables */
1100
1101
1102
1103 /* subroutines and functions */
1104
1105 /*   fortran */
1106
1107 /* ----------------------------------------------------------- */
1108
1109     if (*reset) {
1110         usubc_1.nfxe = *ifxwt;
1111         usubc_1.fxstat[0] = *fx;
1112         usubc_1.fxstat[1] = *fx;
1113         usubc_1.fxstat[2] = *fx;
1114         usubc_1.fxstat[3] = 0.;
1115     } else {
1116         nsv = usubc_1.nfxe;
1117         f1sv = usubc_1.fxstat[0];
1118         usubc_1.nfxe += *ifxwt;
1119         usubc_1.fxstat[0] += *ifxwt * (*fx - usubc_1.fxstat[0]) /
1120                 usubc_1.nfxe;
1121         usubc_1.fxstat[1] = max(usubc_1.fxstat[1],*fx);
1122         usubc_1.fxstat[2] = min(usubc_1.fxstat[2],*fx);
1123 /* Computing MAX */
1124         d__1 = abs(usubc_1.fxstat[1]), d__2 = abs(usubc_1.fxstat[2]), d__1 =
1125                 max(d__1,d__2);
1126         fscale = max(d__1,1.);
1127 /* Computing 2nd power */
1128         d__1 = usubc_1.fxstat[3] / fscale;
1129 /* Computing 2nd power */
1130         d__2 = (usubc_1.fxstat[0] - f1sv) / fscale;
1131 /* Computing 2nd power */
1132         d__3 = (*fx - usubc_1.fxstat[0]) / fscale;
1133         usubc_1.fxstat[3] = fscale * sqrt(((nsv - 1) * (d__1 * d__1) + nsv * (
1134                 d__2 * d__2) + *ifxwt * (d__3 * d__3)) / (usubc_1.nfxe - 1));
1135     }
1136     return 0;
1137 } /* fstats_ */
1138
1139 /* evalf.f -- translated by f2c (version 19991025).
1140    You must link the resulting object file with the libraries:
1141         -lf2c -lm   (in that order)
1142 */
1143
1144 /* Common Block Declarations */
1145
1146 static struct {
1147     doublereal fbonus, sfstop, sfbest;
1148     logical new__;
1149 } isubc_;
1150
1151 #define isubc_1 isubc_
1152
1153 static logical c_true = TRUE_;
1154 static logical c_false = FALSE_;
1155
1156 static int evalf_(D_fp f,void*fdata, integer *ns, integer *ips, doublereal *xs,
1157          integer *n, doublereal *x, doublereal *sfx, integer *nfe)
1158 {
1159     /* System generated locals */
1160     integer i__1;
1161
1162     /* Local variables */
1163     static integer i__;
1164     static doublereal fx;
1165     static logical newbst;
1166
1167 /*                                         Coded by Tom Rowan */
1168 /*                            Department of Computer Sciences */
1169 /*                              University of Texas at Austin */
1170
1171 /* evalf evaluates the function f at a point defined by x */
1172 /* with ns of its components replaced by those in xs. */
1173
1174 /* input */
1175
1176 /*   f      - user supplied function f(n,x) to be optimized */
1177
1178 /*   ns     - subspace dimension */
1179
1180 /*   ips    - permutation vector */
1181
1182 /*   xs     - double precision ns-vector to be mapped to x */
1183
1184 /*   n      - problem dimension */
1185
1186 /*   x      - double precision n-vector */
1187
1188 /*   nfe    - number of function evaluations */
1189
1190 /* output */
1191
1192 /*   sfx    - signed value of f evaluated at x */
1193
1194 /*   nfe    - incremented number of function evaluations */
1195
1196 /* common */
1197
1198
1199
1200
1201
1202 /* local variables */
1203
1204
1205
1206 /* subroutines and functions */
1207
1208
1209 /* ----------------------------------------------------------- */
1210
1211     /* Parameter adjustments */
1212     --ips;
1213     --xs;
1214     --x;
1215
1216     /* Function Body */
1217     i__1 = *ns;
1218     for (i__ = 1; i__ <= i__1; ++i__) {
1219         x[ips[i__]] = xs[i__];
1220 /* L10: */
1221     }
1222     usubc_1.newx = isubc_1.new__ || usubc_1.irepl != 2;
1223     fx = (*f)(*n, &x[1], fdata);
1224     if (usubc_1.irepl == 0) {
1225         if (usubc_1.minf) {
1226             *sfx = fx;
1227         } else {
1228             *sfx = -fx;
1229         }
1230     } else if (isubc_1.new__) {
1231         if (usubc_1.minf) {
1232             *sfx = fx;
1233             newbst = fx < usubc_1.ftest;
1234         } else {
1235             *sfx = -fx;
1236             newbst = fx > usubc_1.ftest;
1237         }
1238         if (usubc_1.initx || newbst) {
1239             if (usubc_1.irepl == 1) {
1240                 fstats_(&fx, &c__1, &c_true);
1241             }
1242             usubc_1.ftest = fx;
1243             isubc_1.sfbest = *sfx;
1244         }
1245     } else {
1246         if (usubc_1.irepl == 1) {
1247             fstats_(&fx, &c__1, &c_false);
1248             fx = usubc_1.fxstat[usubc_1.ifxsw - 1];
1249         }
1250         usubc_1.ftest = fx + isubc_1.fbonus * usubc_1.fxstat[3];
1251         if (usubc_1.minf) {
1252             *sfx = usubc_1.ftest;
1253             isubc_1.sfbest = fx;
1254         } else {
1255             *sfx = -usubc_1.ftest;
1256             isubc_1.sfbest = -fx;
1257         }
1258     }
1259     ++(*nfe);
1260     return 0;
1261 } /* evalf_ */
1262
1263 /* simplx.f -- translated by f2c (version 19991025).
1264    You must link the resulting object file with the libraries:
1265         -lf2c -lm   (in that order)
1266 */
1267
1268 static int simplx_(D_fp f, void *fdata, integer *n, doublereal *step, integer *
1269         ns, integer *ips, nlopt_stopping *stop, logical *cmode, doublereal *x,
1270         doublereal *fx, integer *nfe, doublereal *s, doublereal *fs, integer *
1271         iflag)
1272 {
1273     /* System generated locals */
1274     integer s_dim1, s_offset, i__1;
1275     doublereal d__1, d__2;
1276
1277     /* Local variables */
1278     static integer inew;
1279     static integer npts;
1280     static integer i__, j;
1281     static integer icent;
1282     static logical small;
1283     static integer itemp;
1284     static doublereal fc, fe;
1285     static integer ih, il;
1286     static doublereal fr;
1287     static integer is;
1288     static logical updatc;
1289     static doublereal dum, tol;
1290
1291
1292
1293 /*                                         Coded by Tom Rowan */
1294 /*                            Department of Computer Sciences */
1295 /*                              University of Texas at Austin */
1296
1297 /* simplx uses the Nelder-Mead simplex method to minimize the */
1298 /* function f on a subspace. */
1299
1300 /* input */
1301
1302 /*   f      - function to be minimized, declared external in */
1303 /*            calling routine */
1304
1305 /*   n      - problem dimension */
1306
1307 /*   step   - stepsizes for corresponding components of x */
1308
1309 /*   ns     - subspace dimension */
1310
1311 /*   ips    - permutation vector */
1312
1313 /*   maxnfe - maximum number of function evaluations */
1314
1315 /*   cmode  - logical switch */
1316 /*            = .true.  : continuation of previous call */
1317 /*            = .false. : first call */
1318
1319 /*   x      - starting guess for minimum */
1320
1321 /*   fx     - value of f at x */
1322
1323 /*   nfe    - number of function evaluations */
1324
1325 /*   s      - double precision work array of dimension .ge. */
1326 /*            ns*(ns+3) used to store simplex */
1327
1328 /*   fs     - double precision work array of dimension .ge. */
1329 /*            ns+1 used to store function values of simplex */
1330 /*            vertices */
1331
1332 /* output */
1333
1334 /*   x      - computed minimum */
1335
1336 /*   fx     - value of f at x */
1337
1338 /*   nfe    - incremented number of function evaluations */
1339
1340 /*   iflag  - error flag */
1341 /*            = -1 : maxnfe exceeded */
1342 /*            =  0 : simplex reduced by factor of psi */
1343 /*            =  1 : limit of machine precision */
1344 /*            =  2 : reached fstop */
1345
1346 /* common */
1347
1348
1349
1350
1351
1352 /* local variables */
1353
1354
1355
1356 /* subroutines and functions */
1357
1358 /*   blas */
1359 /*   fortran */
1360
1361 /* ----------------------------------------------------------- */
1362
1363     /* Parameter adjustments */
1364     --x;
1365     --step;
1366     --fs;
1367     s_dim1 = *ns;
1368     s_offset = 1 + s_dim1 * 1;
1369     s -= s_offset;
1370     --ips;
1371
1372     /* Function Body */
1373     if (*cmode) {
1374         goto L50;
1375     }
1376     npts = *ns + 1;
1377     icent = *ns + 2;
1378     itemp = *ns + 3;
1379     updatc = FALSE_;
1380     start_(n, &x[1], &step[1], ns, &ips[1], &s[s_offset], &small);
1381     if (small) {
1382         *iflag = 1;
1383         return 0;
1384     }
1385     if (usubc_1.irepl > 0) {
1386         isubc_1.new__ = FALSE_;
1387         evalf_((D_fp)f,fdata, ns, &ips[1], &s[s_dim1 + 1], n, &x[1], &fs[1], nfe);
1388         stop->nevals++;
1389     } else {
1390         fs[1] = *fx;
1391     }
1392     isubc_1.new__ = TRUE_;
1393     i__1 = npts;
1394     for (j = 2; j <= i__1; ++j) {
1395         evalf_((D_fp)f, fdata,ns, &ips[1], &s[j * s_dim1 + 1], n, &x[1], &fs[j],
1396                 nfe);
1397         stop->nevals++;
1398 /* L10: */
1399     }
1400     il = 1;
1401     order_(&npts, &fs[1], &il, &is, &ih);
1402     tol = usubc_1.psi * dist_(ns, &s[ih * s_dim1 + 1], &s[il * s_dim1 + 1]);
1403
1404 /*     main loop */
1405
1406 L20:
1407     calcc_(ns, &s[s_offset], &ih, &inew, &updatc, &s[icent * s_dim1 + 1]);
1408     updatc = TRUE_;
1409     inew = ih;
1410
1411 /*       reflect */
1412
1413     newpt_(ns, &usubc_1.alpha, &s[icent * s_dim1 + 1], &s[ih * s_dim1 + 1], &
1414             c_true, &s[itemp * s_dim1 + 1], &small);
1415     if (small) {
1416         goto L40;
1417     }
1418     evalf_((D_fp)f,fdata, ns, &ips[1], &s[itemp * s_dim1 + 1], n, &x[1], &fr, nfe);
1419     stop->nevals++;
1420     if (fr < fs[il]) {
1421
1422 /*         expand */
1423
1424         d__1 = -usubc_1.gamma;
1425         newpt_(ns, &d__1, &s[icent * s_dim1 + 1], &s[itemp * s_dim1 + 1], &
1426                 c_true, &s[ih * s_dim1 + 1], &small);
1427         if (small) {
1428             goto L40;
1429         }
1430         evalf_((D_fp)f,fdata, ns, &ips[1], &s[ih * s_dim1 + 1], n, &x[1], &fe, nfe);
1431         stop->nevals++;
1432         if (fe < fr) {
1433             fs[ih] = fe;
1434         } else {
1435             dcopy_(ns, &s[itemp * s_dim1 + 1], &c__1, &s[ih * s_dim1 + 1], &
1436                     c__1);
1437             fs[ih] = fr;
1438         }
1439     } else if (fr < fs[is]) {
1440
1441 /*         accept reflected point */
1442
1443         dcopy_(ns, &s[itemp * s_dim1 + 1], &c__1, &s[ih * s_dim1 + 1], &c__1);
1444         fs[ih] = fr;
1445     } else {
1446
1447 /*         contract */
1448
1449         if (fr > fs[ih]) {
1450             d__1 = -usubc_1.beta;
1451             newpt_(ns, &d__1, &s[icent * s_dim1 + 1], &s[ih * s_dim1 + 1], &
1452                     c_true, &s[itemp * s_dim1 + 1], &small);
1453         } else {
1454             d__1 = -usubc_1.beta;
1455             newpt_(ns, &d__1, &s[icent * s_dim1 + 1], &s[itemp * s_dim1 + 1],
1456                     &c_false, &dum, &small);
1457         }
1458         if (small) {
1459             goto L40;
1460         }
1461         evalf_((D_fp)f,fdata, ns, &ips[1], &s[itemp * s_dim1 + 1], n, &x[1], &fc,
1462                 nfe);
1463         stop->nevals++;
1464 /* Computing MIN */
1465         d__1 = fr, d__2 = fs[ih];
1466         if (fc < min(d__1,d__2)) {
1467             dcopy_(ns, &s[itemp * s_dim1 + 1], &c__1, &s[ih * s_dim1 + 1], &
1468                     c__1);
1469             fs[ih] = fc;
1470         } else {
1471
1472 /*           shrink simplex */
1473
1474             i__1 = npts;
1475             for (j = 1; j <= i__1; ++j) {
1476                 if (j != il) {
1477                     d__1 = -usubc_1.delta;
1478                     newpt_(ns, &d__1, &s[il * s_dim1 + 1], &s[j * s_dim1 + 1],
1479                              &c_false, &dum, &small);
1480                     if (small) {
1481                         goto L40;
1482                     }
1483                     evalf_((D_fp)f,fdata, ns, &ips[1], &s[j * s_dim1 + 1], n, &x[1],
1484                              &fs[j], nfe);
1485                     stop->nevals++;
1486                 }
1487 /* L30: */
1488             }
1489         }
1490         updatc = FALSE_;
1491     }
1492     order_(&npts, &fs[1], &il, &is, &ih);
1493
1494 /*       check termination */
1495
1496 L40:
1497     if (usubc_1.irepl == 0) {
1498         *fx = fs[il];
1499     } else {
1500         *fx = isubc_1.sfbest;
1501     }
1502 L50:
1503     if (*fx < stop->minf_max)
1504          *iflag = 2;
1505     else if (nlopt_stop_evals(stop))
1506          *iflag = -1;
1507     else if (nlopt_stop_time(stop))
1508          *iflag = -10;
1509     else if (dist_(ns, &s[ih * s_dim1 + 1], &s[il * s_dim1 + 1]) <= tol
1510              || small)
1511          *iflag = 0;
1512     else
1513          goto L20;
1514
1515 /*     end main loop, return best point */
1516
1517     i__1 = *ns;
1518     for (i__ = 1; i__ <= i__1; ++i__) {
1519         x[ips[i__]] = s[i__ + il * s_dim1];
1520 /* L60: */
1521     }
1522     return 0;
1523 } /* simplx_ */
1524
1525 /* subopt.f -- translated by f2c (version 19991025).
1526    You must link the resulting object file with the libraries:
1527         -lf2c -lm   (in that order)
1528 */
1529
1530 static int subopt_(integer *n)
1531 {
1532
1533
1534 /*                                         Coded by Tom Rowan */
1535 /*                            Department of Computer Sciences */
1536 /*                              University of Texas at Austin */
1537
1538 /* subopt sets options for subplx. */
1539
1540 /* input */
1541
1542 /*   n      - problem dimension */
1543
1544 /* common */
1545
1546
1547
1548
1549 /* subroutines and functions */
1550
1551 /*   fortran */
1552
1553 /* ----------------------------------------------------------- */
1554
1555 /* *********************************************************** */
1556 /* simplex method strategy parameters */
1557 /* *********************************************************** */
1558
1559 /* alpha  - reflection coefficient */
1560 /*          alpha .gt. 0 */
1561
1562     usubc_1.alpha = 1.;
1563
1564 /* beta   - contraction coefficient */
1565 /*          0 .lt. beta .lt. 1 */
1566
1567     usubc_1.beta = .5;
1568
1569 /* gamma  - expansion coefficient */
1570 /*          gamma .gt. 1 */
1571
1572     usubc_1.gamma = 2.;
1573
1574 /* delta  - shrinkage (massive contraction) coefficient */
1575 /*          0 .lt. delta .lt. 1 */
1576
1577     usubc_1.delta = .5;
1578
1579 /* *********************************************************** */
1580 /* subplex method strategy parameters */
1581 /* *********************************************************** */
1582
1583 /* psi    - simplex reduction coefficient */
1584 /*          0 .lt. psi .lt. 1 */
1585
1586     usubc_1.psi = .25;
1587
1588 /* omega  - step reduction coefficient */
1589 /*          0 .lt. omega .lt. 1 */
1590
1591     usubc_1.omega = .1;
1592
1593 /* nsmin and nsmax specify a range of subspace dimensions. */
1594 /* In addition to satisfying  1 .le. nsmin .le. nsmax .le. n, */
1595 /* nsmin and nsmax must be chosen so that n can be expressed */
1596 /* as a sum of positive integers where each of these integers */
1597 /* ns(i) satisfies   nsmin .le. ns(i) .ge. nsmax. */
1598 /* Specifically, */
1599 /*     nsmin*ceil(n/nsmax) .le. n   must be true. */
1600
1601 /* nsmin  - subspace dimension minimum */
1602
1603     usubc_1.nsmin = min(2,*n);
1604
1605 /* nsmax  - subspace dimension maximum */
1606
1607     usubc_1.nsmax = min(5,*n);
1608
1609 /* *********************************************************** */
1610 /* subplex method special cases */
1611 /* *********************************************************** */
1612 /* nelder-mead simplex method with periodic restarts */
1613 /*   nsmin = nsmax = n */
1614 /* *********************************************************** */
1615 /* nelder-mead simplex method */
1616 /*   nsmin = nsmax = n, psi = small positive */
1617 /* *********************************************************** */
1618
1619 /* irepl, ifxsw, and bonus deal with measurement replication. */
1620 /* Objective functions subject to large amounts of noise can */
1621 /* cause an optimization method to halt at a false optimum. */
1622 /* An expensive solution to this problem is to evaluate f */
1623 /* several times at each point and return the average (or max */
1624 /* or min) of these trials as the function value.  subplx */
1625 /* performs measurement replication only at the current best */
1626 /* point. The longer a point is retained as best, the more */
1627 /* accurate its function value becomes. */
1628
1629 /* The common variable nfxe contains the number of function */
1630 /* evaluations at the current best point. fxstat contains the */
1631 /* mean, max, min, and standard deviation of these trials. */
1632
1633 /* irepl  - measurement replication switch */
1634 /* irepl  = 0, 1, or 2 */
1635 /*        = 0 : no measurement replication */
1636 /*        = 1 : subplx performs measurement replication */
1637 /*        = 2 : user performs measurement replication */
1638 /*              (This is useful when optimizing on the mean, */
1639 /*              max, or min of trials is insufficient. Common */
1640 /*              variable initx is true for first function */
1641 /*              evaluation. newx is true for first trial at */
1642 /*              this point. The user uses subroutine fstats */
1643 /*              within his objective function to maintain */
1644 /*              fxstat. By monitoring newx, the user can tell */
1645 /*              whether to return the function evaluation */
1646 /*              (newx = .true.) or to use the new function */
1647 /*              evaluation to refine the function evaluation */
1648 /*              of the current best point (newx = .false.). */
1649 /*              The common variable ftest gives the function */
1650 /*              value that a new point must beat to be */
1651 /*              considered the new best point.) */
1652
1653     usubc_1.irepl = 0;
1654
1655 /* ifxsw  - measurement replication optimization switch */
1656 /* ifxsw  = 1, 2, or 3 */
1657 /*        = 1 : retain mean of trials as best function value */
1658 /*        = 2 : retain max */
1659 /*        = 3 : retain min */
1660
1661     usubc_1.ifxsw = 1;
1662
1663 /* Since the current best point will also be the most */
1664 /* accurately evaluated point whenever irepl .gt. 0, a bonus */
1665 /* should be added to the function value of the best point */
1666 /* so that the best point is not replaced by a new point */
1667 /* that only appears better because of noise. */
1668 /* subplx uses bonus to determine how many multiples of */
1669 /* fxstat(4) should be added as a bonus to the function */
1670 /* evaluation. (The bonus is adjusted automatically by */
1671 /* subplx when ifxsw or minf is changed.) */
1672
1673 /* bonus  - measurement replication bonus coefficient */
1674 /*          bonus .ge. 0 (normally, bonus = 0 or 1) */
1675 /*        = 0 : bonus not used */
1676 /*        = 1 : bonus used */
1677
1678     usubc_1.bonus = 1.;
1679
1680 /* nfstop = 0 : f(x) is not tested against fstop */
1681 /*        = 1 : if f(x) has reached fstop, subplx returns */
1682 /*              iflag = 2 */
1683 /*        = 2 : (only valid when irepl .gt. 0) */
1684 /*              if f(x) has reached fstop and */
1685 /*              nfxe .gt. nfstop, subplx returns iflag = 2 */
1686
1687     usubc_1.nfstop = 0;
1688
1689 /* fstop  - f target value */
1690 /*          Its usage is determined by the value of nfstop. */
1691
1692 /* minf   - logical switch */
1693 /*        = .true.  : subplx performs minimization */
1694 /*        = .false. : subplx performs maximization */
1695
1696     usubc_1.minf = TRUE_;
1697     return 0;
1698 } /* subopt_ */
1699
1700 /* setstp.f -- translated by f2c (version 19991025).
1701    You must link the resulting object file with the libraries:
1702         -lf2c -lm   (in that order)
1703 */
1704
1705 static double d_sign(doublereal *x, doublereal *y)
1706 {
1707      return copysign(*x, *y);
1708 }
1709
1710 static int setstp_(integer *nsubs, integer *n, doublereal *deltax,
1711                    doublereal *step)
1712 {
1713     /* System generated locals */
1714     integer i__1;
1715     doublereal d__1, d__2, d__3;
1716
1717     /* Builtin functions */
1718 /*    double d_sign(doublereal *, doublereal *); */
1719
1720     /* Local variables */
1721     static integer i__;
1722     static doublereal stpfac;
1723
1724
1725
1726 /*                                         Coded by Tom Rowan */
1727 /*                            Department of Computer Sciences */
1728 /*                              University of Texas at Austin */
1729
1730 /* setstp sets the stepsizes for the corresponding components */
1731 /* of the solution vector. */
1732
1733 /* input */
1734
1735 /*   nsubs  - number of subspaces */
1736
1737 /*   n      - number of components (problem dimension) */
1738
1739 /*   deltax - vector of change in solution vector */
1740
1741 /*   step   - stepsizes for corresponding components of */
1742 /*            solution vector */
1743
1744 /* output */
1745
1746 /*   step   - new stepsizes */
1747
1748 /* common */
1749
1750
1751
1752 /* local variables */
1753
1754
1755
1756 /* subroutines and functions */
1757
1758 /*   blas */
1759 /*   fortran */
1760
1761 /* ----------------------------------------------------------- */
1762
1763 /*     set new step */
1764
1765     /* Parameter adjustments */
1766     --step;
1767     --deltax;
1768
1769     /* Function Body */
1770     if (*nsubs > 1) {
1771 /* Computing MIN */
1772 /* Computing MAX */
1773         d__3 = dasum_(n, &deltax[1], &c__1) / dasum_(n, &step[1], &c__1);
1774         d__1 = max(d__3,usubc_1.omega), d__2 = 1. / usubc_1.omega;
1775         stpfac = min(d__1,d__2);
1776     } else {
1777         stpfac = usubc_1.psi;
1778     }
1779     dscal_(n, &stpfac, &step[1], &c__1);
1780
1781 /*     reorient simplex */
1782
1783     i__1 = *n;
1784     for (i__ = 1; i__ <= i__1; ++i__) {
1785         if (deltax[i__] != 0.) {
1786             step[i__] = d_sign(&step[i__], &deltax[i__]);
1787         } else {
1788             step[i__] = -step[i__];
1789         }
1790 /* L10: */
1791     }
1792     return 0;
1793 } /* setstp_ */
1794
1795 /* subplx.f -- translated by f2c (version 19991025).
1796    You must link the resulting object file with the libraries:
1797         -lf2c -lm   (in that order)
1798 */
1799
1800 static int subplx_(D_fp f, void *fdata, integer *n,
1801                    nlopt_stopping *stop, integer *mode,
1802                    const doublereal *scale, doublereal *x, doublereal *fx,
1803                    integer *nfe, doublereal *work, integer *iwork,
1804                    integer *iflag)
1805 {
1806     /* Initialized data */
1807
1808     static doublereal bnsfac[6] /* was [3][2] */ = { -1.,-2.,0.,1.,0.,2. };
1809
1810     /* System generated locals */
1811     integer i__1;
1812     doublereal d__1, d__2, d__3, d__4, d__5, d__6;
1813
1814     /* Local variables */
1815     static integer i__;
1816     static logical cmode;
1817     static integer istep;
1818     static doublereal xpscl;
1819     static integer nsubs, ipptr;
1820     static integer isptr;
1821     static integer ns, insfnl, ifsptr;
1822     static integer insptr;
1823     static integer istptr;
1824     static doublereal scl, dum;
1825     static integer ins;
1826     static doublereal sfx, sfx_old, *x_old;
1827
1828
1829
1830 /*                                         Coded by Tom Rowan */
1831 /*                            Department of Computer Sciences */
1832 /*                              University of Texas at Austin */
1833
1834 /* subplx uses the subplex method to solve unconstrained */
1835 /* optimization problems.  The method is well suited for */
1836 /* optimizing objective functions that are noisy or are */
1837 /* discontinuous at the solution. */
1838
1839 /* subplx sets default optimization options by calling the */
1840 /* subroutine subopt.  The user can override these defaults */
1841 /* by calling subopt prior to calling subplx, changing the */
1842 /* appropriate common variables, and setting the value of */
1843 /* mode as indicated below. */
1844
1845 /* By default, subplx performs minimization. */
1846
1847 /* input */
1848
1849 /*   f      - user supplied function f(n,x) to be optimized, */
1850 /*            declared external in calling routine */
1851
1852 /*   n      - problem dimension */
1853
1854 /*   tol    - relative error tolerance for x (tol .ge. 0.) */
1855
1856 /*   maxnfe - maximum number of function evaluations */
1857
1858 /*   mode   - integer mode switch with binary expansion */
1859 /*            (bit 1) (bit 0) : */
1860 /*            bit 0 = 0 : first call to subplx */
1861 /*                  = 1 : continuation of previous call */
1862 /*            bit 1 = 0 : use default options */
1863 /*                  = 1 : user set options */
1864
1865 /*   scale  - scale and initial stepsizes for corresponding */
1866 /*            components of x */
1867 /*            (If scale(1) .lt. 0., */
1868 /*            abs(scale(1)) is used for all components of x, */
1869 /*            and scale(2),...,scale(n) are not referenced.) */
1870
1871 /*   x      - starting guess for optimum */
1872
1873 /*   work   - double precision work array of dimension .ge. */
1874 /*            2*n + nsmax*(nsmax+4) + 1 */
1875 /*            (nsmax is set in subroutine subopt. */
1876 /*            default: nsmax = min(5,n)) */
1877
1878 /*   iwork  - integer work array of dimension .ge. */
1879 /*            n + int(n/nsmin) */
1880 /*            (nsmin is set in subroutine subopt. */
1881 /*            default: nsmin = min(2,n)) */
1882
1883 /* output */
1884
1885 /*   x      - computed optimum */
1886
1887 /*   fx     - value of f at x */
1888
1889 /*   nfe    - number of function evaluations */
1890
1891 /*   iflag  - error flag */
1892 /*            = -2 : invalid input */
1893 /*            = -1 : maxnfe exceeded */
1894 /*            =  0 : tol satisfied */
1895 /*            =  1 : limit of machine precision */
1896 /*            =  2 : fstop reached (fstop usage is determined */
1897 /*                   by values of options minf, nfstop, and */
1898 /*                   irepl. default: f(x) not tested against */
1899 /*                   fstop) */
1900 /*            iflag should not be reset between calls to */
1901 /*            subplx. */
1902
1903 /* common */
1904
1905
1906
1907
1908
1909 /* local variables */
1910
1911
1912
1913 /* subroutines and functions */
1914
1915 /*   blas */
1916 /*   fortran */
1917
1918 /* data */
1919
1920     /* Parameter adjustments */
1921     --x;
1922     --scale;
1923     --work;
1924     --iwork;
1925     x_old = work;
1926     work += *n;
1927
1928     /* Function Body */
1929 /* ----------------------------------------------------------- */
1930
1931     if (*mode % 2 == 0) {
1932
1933 /*       first call, check input */
1934
1935         if (*n < 1) {
1936             goto L120;
1937         }
1938         if (scale[1] >= 0.) {
1939             i__1 = *n;
1940             for (i__ = 1; i__ <= i__1; ++i__) {
1941                 xpscl = x[i__] + scale[i__];
1942                 if (xpscl == x[i__]) {
1943                     goto L120;
1944                 }
1945 /* L10: */
1946             }
1947         } else {
1948             scl = abs(scale[1]);
1949             i__1 = *n;
1950             for (i__ = 1; i__ <= i__1; ++i__) {
1951                 xpscl = x[i__] + scl;
1952                 if (xpscl == x[i__]) {
1953                     goto L120;
1954                 }
1955 /* L20: */
1956             }
1957         }
1958         if (*mode / 2 % 2 == 0) {
1959             subopt_(n);
1960         } else {
1961             if (usubc_1.alpha <= 0.) {
1962                 goto L120;
1963             }
1964             if (usubc_1.beta <= 0. || usubc_1.beta >= 1.) {
1965                 goto L120;
1966             }
1967             if (usubc_1.gamma <= 1.) {
1968                 goto L120;
1969             }
1970             if (usubc_1.delta <= 0. || usubc_1.delta >= 1.) {
1971                 goto L120;
1972             }
1973             if (usubc_1.psi <= 0. || usubc_1.psi >= 1.) {
1974                 goto L120;
1975             }
1976             if (usubc_1.omega <= 0. || usubc_1.omega >= 1.) {
1977                 goto L120;
1978             }
1979             if (usubc_1.nsmin < 1 || usubc_1.nsmax < usubc_1.nsmin || *n <
1980                     usubc_1.nsmax) {
1981                 goto L120;
1982             }
1983             if (*n < ((*n - 1) / usubc_1.nsmax + 1) * usubc_1.nsmin) {
1984                 goto L120;
1985             }
1986             if (usubc_1.irepl < 0 || usubc_1.irepl > 2) {
1987                 goto L120;
1988             }
1989             if (usubc_1.ifxsw < 1 || usubc_1.ifxsw > 3) {
1990                 goto L120;
1991             }
1992             if (usubc_1.bonus < 0.) {
1993                 goto L120;
1994             }
1995             if (usubc_1.nfstop < 0) {
1996                 goto L120;
1997             }
1998         }
1999
2000 /*       initialization */
2001
2002         istptr = *n + 1;
2003         isptr = istptr + *n;
2004         ifsptr = isptr + usubc_1.nsmax * (usubc_1.nsmax + 3);
2005         insptr = *n + 1;
2006         if (scale[1] > 0.) {
2007             dcopy_(n, &scale[1], &c__1, &work[1], &c__1);
2008             dcopy_(n, &scale[1], &c__1, &work[istptr], &c__1);
2009         } else {
2010             dcopy_(n, &scl, &c__0, &work[1], &c__1);
2011             dcopy_(n, &scl, &c__0, &work[istptr], &c__1);
2012         }
2013         i__1 = *n;
2014         for (i__ = 1; i__ <= i__1; ++i__) {
2015             iwork[i__] = i__;
2016 /* L30: */
2017         }
2018         *nfe = 0;
2019         usubc_1.nfxe = 1;
2020         if (usubc_1.irepl == 0) {
2021             isubc_1.fbonus = 0.;
2022         } else if (usubc_1.minf) {
2023             isubc_1.fbonus = bnsfac[usubc_1.ifxsw - 1] * usubc_1.bonus;
2024         } else {
2025             isubc_1.fbonus = bnsfac[usubc_1.ifxsw + 2] * usubc_1.bonus;
2026         }
2027         if (usubc_1.nfstop == 0) {
2028             isubc_1.sfstop = 0.;
2029         } else if (usubc_1.minf) {
2030             isubc_1.sfstop = usubc_1.fstop;
2031         } else {
2032             isubc_1.sfstop = -usubc_1.fstop;
2033         }
2034         usubc_1.ftest = 0.;
2035         cmode = FALSE_;
2036         isubc_1.new__ = TRUE_;
2037         usubc_1.initx = TRUE_;
2038         evalf_((D_fp)f, fdata, &c__0, &iwork[1], &dum, n, &x[1], &sfx, nfe);
2039         stop->nevals++;
2040         usubc_1.initx = FALSE_;
2041     } else {
2042
2043 /*       continuation of previous call */
2044
2045         if (*iflag == 2) {
2046             if (usubc_1.minf) {
2047                 isubc_1.sfstop = usubc_1.fstop;
2048             } else {
2049                 isubc_1.sfstop = -usubc_1.fstop;
2050             }
2051             cmode = TRUE_;
2052             goto L70;
2053         } else if (*iflag == -1) {
2054             cmode = TRUE_;
2055             goto L70;
2056         } else if (*iflag == 0) {
2057             cmode = FALSE_;
2058             goto L90;
2059         } else {
2060             return 0;
2061         }
2062     }
2063
2064 /*     subplex loop */
2065
2066 L40:
2067     i__1 = *n;
2068     for (i__ = 1; i__ <= i__1; ++i__) {
2069         work[i__] = (d__1 = work[i__], abs(d__1));
2070 /* L50: */
2071     }
2072     sortd_(n, &work[1], &iwork[1]);
2073     partx_(n, &iwork[1], &work[1], &nsubs, &iwork[insptr]);
2074     dcopy_(n, &x[1], &c__1, &work[1], &c__1);
2075     ins = insptr;
2076     insfnl = insptr + nsubs - 1;
2077     ipptr = 1;
2078
2079 /*       simplex loop */
2080     sfx_old = sfx;
2081     memcpy(&x_old[1], &x[1], sizeof(doublereal) * *n);
2082 L60:
2083     ns = iwork[ins];
2084 L70:
2085     simplx_((D_fp)f, fdata, n, &work[istptr], &ns, &iwork[ipptr], stop, &cmode, &x[1], &sfx, nfe, &work[isptr], &work[ifsptr], iflag);
2086     cmode = FALSE_;
2087     if (*iflag != 0) {
2088         goto L110;
2089     }
2090     if (ins < insfnl) {
2091         ++ins;
2092         ipptr += ns;
2093         goto L60;
2094     }
2095
2096 /*       end simplex loop */
2097
2098     i__1 = *n;
2099     for (i__ = 1; i__ <= i__1; ++i__) {
2100         work[i__] = x[i__] - work[i__];
2101 /* L80: */
2102     }
2103
2104 /*       check termination */
2105
2106     if (nlopt_stop_ftol(stop, sfx, sfx_old)) {
2107          *iflag = 20;
2108          goto L110;
2109     }
2110     if (nlopt_stop_x(stop, &x[1], &x_old[1])) {
2111          *iflag = 0;
2112          goto L110;
2113     }
2114
2115 L90:
2116     istep = istptr;
2117     i__1 = *n;
2118     for (i__ = 1; i__ <= i__1; ++i__) {
2119 /* Computing MAX */
2120         d__4 = (d__2 = work[i__], abs(d__2)), d__5 = (d__1 = work[istep], abs(
2121                 d__1)) * usubc_1.psi;
2122 /* Computing MAX */
2123         d__6 = (d__3 = x[i__], abs(d__3));
2124         if (max(d__4,d__5) / max(d__6,1.) > stop->xtol_rel) {
2125             setstp_(&nsubs, n, &work[1], &work[istptr]);
2126             goto L40;
2127         }
2128         ++istep;
2129 /* L100: */
2130     }
2131
2132 /*     end subplex loop */
2133
2134     *iflag = 0;
2135 L110:
2136     if (usubc_1.minf) {
2137         *fx = sfx;
2138     } else {
2139         *fx = -sfx;
2140     }
2141     return 0;
2142
2143 /*     invalid input */
2144
2145 L120:
2146     *iflag = -2;
2147     return 0;
2148 } /* subplx_ */
2149
2150 /****************************************************************************/
2151 /****************************************************************************/
2152
2153 /* front-end for subplex routines */
2154
2155 /* Wrapper around f2c'ed subplx_ routine, for multidimensinal
2156    unconstrained optimization:
2157
2158    Parameters:
2159
2160    f: function f(n,x,fdata) to be optimized
2161    n: problem dimension
2162    minf: (output) value of f at minimum
2163    x[n]: (input) starting guess position, (output) computed minimum
2164    fdata: data pointer passed to f
2165
2166    old args:
2167    tol: relative error tolerance for x
2168    maxnfe: maximum number of function evaluations
2169    minf_max, use_minf_max: if use_minf_max, stop when f <= minf_max
2170
2171    new args: nlopt_stopping *stop (stopping criteria)
2172
2173    scale[n]: (input) scale & initial stepsizes for components of x
2174              (if *scale < 0, |*scale| is used for all components)
2175
2176    Return value:
2177             = -2 : invalid input
2178             = -10 : maxtime exceeded
2179             = -1 : maxevals exceeded
2180             =  0 : tol satisfied
2181             =  1 : limit of machine precision
2182             =  2 : fstop reached (fstop usage is determined by values of
2183                    options minf, nfstop, and irepl. default: f(x) not
2184                    tested against fstop)
2185             = 20 : ftol reached
2186             = -200 : out of memory
2187 */
2188 int nlopt_subplex(subplex_func f, double *minf, double *x, int n, void *fdata,
2189             nlopt_stopping *stop,
2190             const double *scale)
2191 {
2192      int mode = 0, *iwork, nsmax, nsmin, errflag, nfe;
2193      double *work;
2194
2195      nsmax = min(5,n);
2196      nsmin = min(2,n);
2197      work = (double*) malloc(sizeof(double) * (3*n + nsmax*(nsmax+4) + 1));
2198      if (!work)
2199           return -200;
2200      iwork = (int*) malloc(sizeof(int) * (n + n/nsmin + 1));
2201      if (!iwork) {
2202           free(work);
2203           return -200;
2204      }
2205
2206      subplx_(f,fdata, &n,
2207              stop, &mode,
2208              scale, x,
2209              minf, &nfe,
2210              work, iwork, &errflag);
2211
2212      free(iwork);
2213      free(work);
2214
2215      return errflag;
2216 }