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 (nlopt_stop_forced(stop))
1504          *iflag = -20;
1505     else if (*fx < stop->minf_max)
1506          *iflag = 2;
1507     else if (nlopt_stop_evals(stop))
1508          *iflag = -1;
1509     else if (nlopt_stop_time(stop))
1510          *iflag = -10;
1511     else if (dist_(ns, &s[ih * s_dim1 + 1], &s[il * s_dim1 + 1]) <= tol
1512              || small)
1513          *iflag = 0;
1514     else
1515          goto L20;
1516
1517 /*     end main loop, return best point */
1518
1519     i__1 = *ns;
1520     for (i__ = 1; i__ <= i__1; ++i__) {
1521         x[ips[i__]] = s[i__ + il * s_dim1];
1522 /* L60: */
1523     }
1524     return 0;
1525 } /* simplx_ */
1526
1527 /* subopt.f -- translated by f2c (version 19991025).
1528    You must link the resulting object file with the libraries:
1529         -lf2c -lm   (in that order)
1530 */
1531
1532 static int subopt_(integer *n)
1533 {
1534
1535
1536 /*                                         Coded by Tom Rowan */
1537 /*                            Department of Computer Sciences */
1538 /*                              University of Texas at Austin */
1539
1540 /* subopt sets options for subplx. */
1541
1542 /* input */
1543
1544 /*   n      - problem dimension */
1545
1546 /* common */
1547
1548
1549
1550
1551 /* subroutines and functions */
1552
1553 /*   fortran */
1554
1555 /* ----------------------------------------------------------- */
1556
1557 /* *********************************************************** */
1558 /* simplex method strategy parameters */
1559 /* *********************************************************** */
1560
1561 /* alpha  - reflection coefficient */
1562 /*          alpha .gt. 0 */
1563
1564     usubc_1.alpha = 1.;
1565
1566 /* beta   - contraction coefficient */
1567 /*          0 .lt. beta .lt. 1 */
1568
1569     usubc_1.beta = .5;
1570
1571 /* gamma  - expansion coefficient */
1572 /*          gamma .gt. 1 */
1573
1574     usubc_1.gamma = 2.;
1575
1576 /* delta  - shrinkage (massive contraction) coefficient */
1577 /*          0 .lt. delta .lt. 1 */
1578
1579     usubc_1.delta = .5;
1580
1581 /* *********************************************************** */
1582 /* subplex method strategy parameters */
1583 /* *********************************************************** */
1584
1585 /* psi    - simplex reduction coefficient */
1586 /*          0 .lt. psi .lt. 1 */
1587
1588     usubc_1.psi = .25;
1589
1590 /* omega  - step reduction coefficient */
1591 /*          0 .lt. omega .lt. 1 */
1592
1593     usubc_1.omega = .1;
1594
1595 /* nsmin and nsmax specify a range of subspace dimensions. */
1596 /* In addition to satisfying  1 .le. nsmin .le. nsmax .le. n, */
1597 /* nsmin and nsmax must be chosen so that n can be expressed */
1598 /* as a sum of positive integers where each of these integers */
1599 /* ns(i) satisfies   nsmin .le. ns(i) .ge. nsmax. */
1600 /* Specifically, */
1601 /*     nsmin*ceil(n/nsmax) .le. n   must be true. */
1602
1603 /* nsmin  - subspace dimension minimum */
1604
1605     usubc_1.nsmin = min(2,*n);
1606
1607 /* nsmax  - subspace dimension maximum */
1608
1609     usubc_1.nsmax = min(5,*n);
1610
1611 /* *********************************************************** */
1612 /* subplex method special cases */
1613 /* *********************************************************** */
1614 /* nelder-mead simplex method with periodic restarts */
1615 /*   nsmin = nsmax = n */
1616 /* *********************************************************** */
1617 /* nelder-mead simplex method */
1618 /*   nsmin = nsmax = n, psi = small positive */
1619 /* *********************************************************** */
1620
1621 /* irepl, ifxsw, and bonus deal with measurement replication. */
1622 /* Objective functions subject to large amounts of noise can */
1623 /* cause an optimization method to halt at a false optimum. */
1624 /* An expensive solution to this problem is to evaluate f */
1625 /* several times at each point and return the average (or max */
1626 /* or min) of these trials as the function value.  subplx */
1627 /* performs measurement replication only at the current best */
1628 /* point. The longer a point is retained as best, the more */
1629 /* accurate its function value becomes. */
1630
1631 /* The common variable nfxe contains the number of function */
1632 /* evaluations at the current best point. fxstat contains the */
1633 /* mean, max, min, and standard deviation of these trials. */
1634
1635 /* irepl  - measurement replication switch */
1636 /* irepl  = 0, 1, or 2 */
1637 /*        = 0 : no measurement replication */
1638 /*        = 1 : subplx performs measurement replication */
1639 /*        = 2 : user performs measurement replication */
1640 /*              (This is useful when optimizing on the mean, */
1641 /*              max, or min of trials is insufficient. Common */
1642 /*              variable initx is true for first function */
1643 /*              evaluation. newx is true for first trial at */
1644 /*              this point. The user uses subroutine fstats */
1645 /*              within his objective function to maintain */
1646 /*              fxstat. By monitoring newx, the user can tell */
1647 /*              whether to return the function evaluation */
1648 /*              (newx = .true.) or to use the new function */
1649 /*              evaluation to refine the function evaluation */
1650 /*              of the current best point (newx = .false.). */
1651 /*              The common variable ftest gives the function */
1652 /*              value that a new point must beat to be */
1653 /*              considered the new best point.) */
1654
1655     usubc_1.irepl = 0;
1656
1657 /* ifxsw  - measurement replication optimization switch */
1658 /* ifxsw  = 1, 2, or 3 */
1659 /*        = 1 : retain mean of trials as best function value */
1660 /*        = 2 : retain max */
1661 /*        = 3 : retain min */
1662
1663     usubc_1.ifxsw = 1;
1664
1665 /* Since the current best point will also be the most */
1666 /* accurately evaluated point whenever irepl .gt. 0, a bonus */
1667 /* should be added to the function value of the best point */
1668 /* so that the best point is not replaced by a new point */
1669 /* that only appears better because of noise. */
1670 /* subplx uses bonus to determine how many multiples of */
1671 /* fxstat(4) should be added as a bonus to the function */
1672 /* evaluation. (The bonus is adjusted automatically by */
1673 /* subplx when ifxsw or minf is changed.) */
1674
1675 /* bonus  - measurement replication bonus coefficient */
1676 /*          bonus .ge. 0 (normally, bonus = 0 or 1) */
1677 /*        = 0 : bonus not used */
1678 /*        = 1 : bonus used */
1679
1680     usubc_1.bonus = 1.;
1681
1682 /* nfstop = 0 : f(x) is not tested against fstop */
1683 /*        = 1 : if f(x) has reached fstop, subplx returns */
1684 /*              iflag = 2 */
1685 /*        = 2 : (only valid when irepl .gt. 0) */
1686 /*              if f(x) has reached fstop and */
1687 /*              nfxe .gt. nfstop, subplx returns iflag = 2 */
1688
1689     usubc_1.nfstop = 0;
1690
1691 /* fstop  - f target value */
1692 /*          Its usage is determined by the value of nfstop. */
1693
1694 /* minf   - logical switch */
1695 /*        = .true.  : subplx performs minimization */
1696 /*        = .false. : subplx performs maximization */
1697
1698     usubc_1.minf = TRUE_;
1699     return 0;
1700 } /* subopt_ */
1701
1702 /* setstp.f -- translated by f2c (version 19991025).
1703    You must link the resulting object file with the libraries:
1704         -lf2c -lm   (in that order)
1705 */
1706
1707 static double d_sign(doublereal *x, doublereal *y)
1708 {
1709      return copysign(*x, *y);
1710 }
1711
1712 static int setstp_(integer *nsubs, integer *n, doublereal *deltax,
1713                    doublereal *step)
1714 {
1715     /* System generated locals */
1716     integer i__1;
1717     doublereal d__1, d__2, d__3;
1718
1719     /* Builtin functions */
1720 /*    double d_sign(doublereal *, doublereal *); */
1721
1722     /* Local variables */
1723     static integer i__;
1724     static doublereal stpfac;
1725
1726
1727
1728 /*                                         Coded by Tom Rowan */
1729 /*                            Department of Computer Sciences */
1730 /*                              University of Texas at Austin */
1731
1732 /* setstp sets the stepsizes for the corresponding components */
1733 /* of the solution vector. */
1734
1735 /* input */
1736
1737 /*   nsubs  - number of subspaces */
1738
1739 /*   n      - number of components (problem dimension) */
1740
1741 /*   deltax - vector of change in solution vector */
1742
1743 /*   step   - stepsizes for corresponding components of */
1744 /*            solution vector */
1745
1746 /* output */
1747
1748 /*   step   - new stepsizes */
1749
1750 /* common */
1751
1752
1753
1754 /* local variables */
1755
1756
1757
1758 /* subroutines and functions */
1759
1760 /*   blas */
1761 /*   fortran */
1762
1763 /* ----------------------------------------------------------- */
1764
1765 /*     set new step */
1766
1767     /* Parameter adjustments */
1768     --step;
1769     --deltax;
1770
1771     /* Function Body */
1772     if (*nsubs > 1) {
1773 /* Computing MIN */
1774 /* Computing MAX */
1775         d__3 = dasum_(n, &deltax[1], &c__1) / dasum_(n, &step[1], &c__1);
1776         d__1 = max(d__3,usubc_1.omega), d__2 = 1. / usubc_1.omega;
1777         stpfac = min(d__1,d__2);
1778     } else {
1779         stpfac = usubc_1.psi;
1780     }
1781     dscal_(n, &stpfac, &step[1], &c__1);
1782
1783 /*     reorient simplex */
1784
1785     i__1 = *n;
1786     for (i__ = 1; i__ <= i__1; ++i__) {
1787         if (deltax[i__] != 0.) {
1788             step[i__] = d_sign(&step[i__], &deltax[i__]);
1789         } else {
1790             step[i__] = -step[i__];
1791         }
1792 /* L10: */
1793     }
1794     return 0;
1795 } /* setstp_ */
1796
1797 /* subplx.f -- translated by f2c (version 19991025).
1798    You must link the resulting object file with the libraries:
1799         -lf2c -lm   (in that order)
1800 */
1801
1802 static int subplx_(D_fp f, void *fdata, integer *n,
1803                    nlopt_stopping *stop, integer *mode,
1804                    const doublereal *scale, doublereal *x, doublereal *fx,
1805                    integer *nfe, doublereal *work, integer *iwork,
1806                    integer *iflag)
1807 {
1808     /* Initialized data */
1809
1810     static doublereal bnsfac[6] /* was [3][2] */ = { -1.,-2.,0.,1.,0.,2. };
1811
1812     /* System generated locals */
1813     integer i__1;
1814     doublereal d__1, d__2, d__3, d__4, d__5, d__6;
1815
1816     /* Local variables */
1817     static integer i__;
1818     static logical cmode;
1819     static integer istep;
1820     static doublereal xpscl;
1821     static integer nsubs, ipptr;
1822     static integer isptr;
1823     static integer ns, insfnl, ifsptr;
1824     static integer insptr;
1825     static integer istptr;
1826     static doublereal scl, dum;
1827     static integer ins;
1828     static doublereal sfx, sfx_old, *x_old;
1829
1830
1831
1832 /*                                         Coded by Tom Rowan */
1833 /*                            Department of Computer Sciences */
1834 /*                              University of Texas at Austin */
1835
1836 /* subplx uses the subplex method to solve unconstrained */
1837 /* optimization problems.  The method is well suited for */
1838 /* optimizing objective functions that are noisy or are */
1839 /* discontinuous at the solution. */
1840
1841 /* subplx sets default optimization options by calling the */
1842 /* subroutine subopt.  The user can override these defaults */
1843 /* by calling subopt prior to calling subplx, changing the */
1844 /* appropriate common variables, and setting the value of */
1845 /* mode as indicated below. */
1846
1847 /* By default, subplx performs minimization. */
1848
1849 /* input */
1850
1851 /*   f      - user supplied function f(n,x) to be optimized, */
1852 /*            declared external in calling routine */
1853
1854 /*   n      - problem dimension */
1855
1856 /*   tol    - relative error tolerance for x (tol .ge. 0.) */
1857
1858 /*   maxnfe - maximum number of function evaluations */
1859
1860 /*   mode   - integer mode switch with binary expansion */
1861 /*            (bit 1) (bit 0) : */
1862 /*            bit 0 = 0 : first call to subplx */
1863 /*                  = 1 : continuation of previous call */
1864 /*            bit 1 = 0 : use default options */
1865 /*                  = 1 : user set options */
1866
1867 /*   scale  - scale and initial stepsizes for corresponding */
1868 /*            components of x */
1869 /*            (If scale(1) .lt. 0., */
1870 /*            abs(scale(1)) is used for all components of x, */
1871 /*            and scale(2),...,scale(n) are not referenced.) */
1872
1873 /*   x      - starting guess for optimum */
1874
1875 /*   work   - double precision work array of dimension .ge. */
1876 /*            2*n + nsmax*(nsmax+4) + 1 */
1877 /*            (nsmax is set in subroutine subopt. */
1878 /*            default: nsmax = min(5,n)) */
1879
1880 /*   iwork  - integer work array of dimension .ge. */
1881 /*            n + int(n/nsmin) */
1882 /*            (nsmin is set in subroutine subopt. */
1883 /*            default: nsmin = min(2,n)) */
1884
1885 /* output */
1886
1887 /*   x      - computed optimum */
1888
1889 /*   fx     - value of f at x */
1890
1891 /*   nfe    - number of function evaluations */
1892
1893 /*   iflag  - error flag */
1894 /*            = -2 : invalid input */
1895 /*            = -1 : maxnfe exceeded */
1896 /*            =  0 : tol satisfied */
1897 /*            =  1 : limit of machine precision */
1898 /*            =  2 : fstop reached (fstop usage is determined */
1899 /*                   by values of options minf, nfstop, and */
1900 /*                   irepl. default: f(x) not tested against */
1901 /*                   fstop) */
1902 /*            iflag should not be reset between calls to */
1903 /*            subplx. */
1904
1905 /* common */
1906
1907
1908
1909
1910
1911 /* local variables */
1912
1913
1914
1915 /* subroutines and functions */
1916
1917 /*   blas */
1918 /*   fortran */
1919
1920 /* data */
1921
1922     /* Parameter adjustments */
1923     --x;
1924     --scale;
1925     --work;
1926     --iwork;
1927     x_old = work;
1928     work += *n;
1929
1930     /* Function Body */
1931 /* ----------------------------------------------------------- */
1932
1933     if (*mode % 2 == 0) {
1934
1935 /*       first call, check input */
1936
1937         if (*n < 1) {
1938             goto L120;
1939         }
1940         if (scale[1] >= 0.) {
1941             i__1 = *n;
1942             for (i__ = 1; i__ <= i__1; ++i__) {
1943                 xpscl = x[i__] + scale[i__];
1944                 if (xpscl == x[i__]) {
1945                     goto L120;
1946                 }
1947 /* L10: */
1948             }
1949         } else {
1950             scl = abs(scale[1]);
1951             i__1 = *n;
1952             for (i__ = 1; i__ <= i__1; ++i__) {
1953                 xpscl = x[i__] + scl;
1954                 if (xpscl == x[i__]) {
1955                     goto L120;
1956                 }
1957 /* L20: */
1958             }
1959         }
1960         if (*mode / 2 % 2 == 0) {
1961             subopt_(n);
1962         } else {
1963             if (usubc_1.alpha <= 0.) {
1964                 goto L120;
1965             }
1966             if (usubc_1.beta <= 0. || usubc_1.beta >= 1.) {
1967                 goto L120;
1968             }
1969             if (usubc_1.gamma <= 1.) {
1970                 goto L120;
1971             }
1972             if (usubc_1.delta <= 0. || usubc_1.delta >= 1.) {
1973                 goto L120;
1974             }
1975             if (usubc_1.psi <= 0. || usubc_1.psi >= 1.) {
1976                 goto L120;
1977             }
1978             if (usubc_1.omega <= 0. || usubc_1.omega >= 1.) {
1979                 goto L120;
1980             }
1981             if (usubc_1.nsmin < 1 || usubc_1.nsmax < usubc_1.nsmin || *n <
1982                     usubc_1.nsmax) {
1983                 goto L120;
1984             }
1985             if (*n < ((*n - 1) / usubc_1.nsmax + 1) * usubc_1.nsmin) {
1986                 goto L120;
1987             }
1988             if (usubc_1.irepl < 0 || usubc_1.irepl > 2) {
1989                 goto L120;
1990             }
1991             if (usubc_1.ifxsw < 1 || usubc_1.ifxsw > 3) {
1992                 goto L120;
1993             }
1994             if (usubc_1.bonus < 0.) {
1995                 goto L120;
1996             }
1997             if (usubc_1.nfstop < 0) {
1998                 goto L120;
1999             }
2000         }
2001
2002 /*       initialization */
2003
2004         istptr = *n + 1;
2005         isptr = istptr + *n;
2006         ifsptr = isptr + usubc_1.nsmax * (usubc_1.nsmax + 3);
2007         insptr = *n + 1;
2008         if (scale[1] > 0.) {
2009             dcopy_(n, &scale[1], &c__1, &work[1], &c__1);
2010             dcopy_(n, &scale[1], &c__1, &work[istptr], &c__1);
2011         } else {
2012             dcopy_(n, &scl, &c__0, &work[1], &c__1);
2013             dcopy_(n, &scl, &c__0, &work[istptr], &c__1);
2014         }
2015         i__1 = *n;
2016         for (i__ = 1; i__ <= i__1; ++i__) {
2017             iwork[i__] = i__;
2018 /* L30: */
2019         }
2020         *nfe = 0;
2021         usubc_1.nfxe = 1;
2022         if (usubc_1.irepl == 0) {
2023             isubc_1.fbonus = 0.;
2024         } else if (usubc_1.minf) {
2025             isubc_1.fbonus = bnsfac[usubc_1.ifxsw - 1] * usubc_1.bonus;
2026         } else {
2027             isubc_1.fbonus = bnsfac[usubc_1.ifxsw + 2] * usubc_1.bonus;
2028         }
2029         if (usubc_1.nfstop == 0) {
2030             isubc_1.sfstop = 0.;
2031         } else if (usubc_1.minf) {
2032             isubc_1.sfstop = usubc_1.fstop;
2033         } else {
2034             isubc_1.sfstop = -usubc_1.fstop;
2035         }
2036         usubc_1.ftest = 0.;
2037         cmode = FALSE_;
2038         isubc_1.new__ = TRUE_;
2039         usubc_1.initx = TRUE_;
2040         evalf_((D_fp)f, fdata, &c__0, &iwork[1], &dum, n, &x[1], &sfx, nfe);
2041         stop->nevals++;
2042         usubc_1.initx = FALSE_;
2043     } else {
2044
2045 /*       continuation of previous call */
2046
2047         if (*iflag == 2) {
2048             if (usubc_1.minf) {
2049                 isubc_1.sfstop = usubc_1.fstop;
2050             } else {
2051                 isubc_1.sfstop = -usubc_1.fstop;
2052             }
2053             cmode = TRUE_;
2054             goto L70;
2055         } else if (*iflag == -1) {
2056             cmode = TRUE_;
2057             goto L70;
2058         } else if (*iflag == 0) {
2059             cmode = FALSE_;
2060             goto L90;
2061         } else {
2062             return 0;
2063         }
2064     }
2065
2066 /*     subplex loop */
2067
2068 L40:
2069     i__1 = *n;
2070     for (i__ = 1; i__ <= i__1; ++i__) {
2071         work[i__] = (d__1 = work[i__], abs(d__1));
2072 /* L50: */
2073     }
2074     sortd_(n, &work[1], &iwork[1]);
2075     partx_(n, &iwork[1], &work[1], &nsubs, &iwork[insptr]);
2076     dcopy_(n, &x[1], &c__1, &work[1], &c__1);
2077     ins = insptr;
2078     insfnl = insptr + nsubs - 1;
2079     ipptr = 1;
2080
2081 /*       simplex loop */
2082     sfx_old = sfx;
2083     memcpy(&x_old[1], &x[1], sizeof(doublereal) * *n);
2084 L60:
2085     ns = iwork[ins];
2086 L70:
2087     simplx_((D_fp)f, fdata, n, &work[istptr], &ns, &iwork[ipptr], stop, &cmode, &x[1], &sfx, nfe, &work[isptr], &work[ifsptr], iflag);
2088     cmode = FALSE_;
2089     if (*iflag != 0) {
2090         goto L110;
2091     }
2092     if (ins < insfnl) {
2093         ++ins;
2094         ipptr += ns;
2095         goto L60;
2096     }
2097
2098 /*       end simplex loop */
2099
2100     i__1 = *n;
2101     for (i__ = 1; i__ <= i__1; ++i__) {
2102         work[i__] = x[i__] - work[i__];
2103 /* L80: */
2104     }
2105
2106 /*       check termination */
2107
2108     if (nlopt_stop_ftol(stop, sfx, sfx_old)) {
2109          *iflag = 20;
2110          goto L110;
2111     }
2112     if (nlopt_stop_x(stop, &x[1], &x_old[1])) {
2113          *iflag = 0;
2114          goto L110;
2115     }
2116
2117 L90:
2118     istep = istptr;
2119     i__1 = *n;
2120     for (i__ = 1; i__ <= i__1; ++i__) {
2121 /* Computing MAX */
2122         d__4 = (d__2 = work[i__], abs(d__2)), d__5 = (d__1 = work[istep], abs(
2123                 d__1)) * usubc_1.psi;
2124 /* Computing MAX */
2125         d__6 = (d__3 = x[i__], abs(d__3));
2126         if (max(d__4,d__5) / max(d__6,1.) > stop->xtol_rel) {
2127             setstp_(&nsubs, n, &work[1], &work[istptr]);
2128             goto L40;
2129         }
2130         ++istep;
2131 /* L100: */
2132     }
2133
2134 /*     end subplex loop */
2135
2136     *iflag = 0;
2137 L110:
2138     if (usubc_1.minf) {
2139         *fx = sfx;
2140     } else {
2141         *fx = -sfx;
2142     }
2143     return 0;
2144
2145 /*     invalid input */
2146
2147 L120:
2148     *iflag = -2;
2149     return 0;
2150 } /* subplx_ */
2151
2152 /****************************************************************************/
2153 /****************************************************************************/
2154
2155 /* front-end for subplex routines */
2156
2157 /* Wrapper around f2c'ed subplx_ routine, for multidimensinal
2158    unconstrained optimization:
2159
2160    Parameters:
2161
2162    f: function f(n,x,fdata) to be optimized
2163    n: problem dimension
2164    minf: (output) value of f at minimum
2165    x[n]: (input) starting guess position, (output) computed minimum
2166    fdata: data pointer passed to f
2167
2168    old args:
2169    tol: relative error tolerance for x
2170    maxnfe: maximum number of function evaluations
2171    minf_max, use_minf_max: if use_minf_max, stop when f <= minf_max
2172
2173    new args: nlopt_stopping *stop (stopping criteria)
2174
2175    scale[n]: (input) scale & initial stepsizes for components of x
2176              (if *scale < 0, |*scale| is used for all components)
2177
2178    Return value:
2179             = -2 : invalid input
2180             = -10 : maxtime exceeded
2181             = -1 : maxevals exceeded
2182             =  0 : tol satisfied
2183             =  1 : limit of machine precision
2184             =  2 : fstop reached (fstop usage is determined by values of
2185                    options minf, nfstop, and irepl. default: f(x) not
2186                    tested against fstop)
2187             = 20 : ftol reached
2188             = -200 : out of memory
2189 */
2190 int nlopt_subplex(subplex_func f, double *minf, double *x, int n, void *fdata,
2191             nlopt_stopping *stop,
2192             const double *scale)
2193 {
2194      int mode = 0, *iwork, nsmax, nsmin, errflag, nfe;
2195      double *work;
2196
2197      nsmax = min(5,n);
2198      nsmin = min(2,n);
2199      work = (double*) malloc(sizeof(double) * (3*n + nsmax*(nsmax+4) + 1));
2200      if (!work)
2201           return -200;
2202      iwork = (int*) malloc(sizeof(int) * (n + n/nsmin + 1));
2203      if (!iwork) {
2204           free(work);
2205           return -200;
2206      }
2207
2208      subplx_(f,fdata, &n,
2209              stop, &mode,
2210              scale, x,
2211              minf, &nfe,
2212              work, iwork, &errflag);
2213
2214      free(iwork);
2215      free(work);
2216
2217      return errflag;
2218 }