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