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