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