chiark / gitweb /
initial checkin
authorstevenj <stevenj@alum.mit.edu>
Wed, 22 Aug 2007 22:10:32 +0000 (18:10 -0400)
committerstevenj <stevenj@alum.mit.edu>
Wed, 22 Aug 2007 22:10:32 +0000 (18:10 -0400)
darcs-hash:20070822221032-c8de0-6fe4c3d6af0aa9a73f141cdbda9fabc08fe7805a.gz

29 files changed:
direct/AUTHORS [new file with mode: 0644]
direct/COPYRIGHT [new file with mode: 0644]
direct/DIRect.c [new file with mode: 0644]
direct/DIRparallel.c [new file with mode: 0644]
direct/DIRserial.c [new file with mode: 0644]
direct/DIRsubrout.c [new file with mode: 0644]
direct/direct-internal.h [new file with mode: 0644]
direct/direct.h [new file with mode: 0644]
direct/direct_wrap.c [new file with mode: 0644]
direct/tstc.c [new file with mode: 0644]
direct/userguide.ps.gz [new file with mode: 0644]
stogo/README [new file with mode: 0644]
stogo/global.cc [new file with mode: 0644]
stogo/global.h [new file with mode: 0644]
stogo/linalg.cc [new file with mode: 0644]
stogo/linalg.h [new file with mode: 0644]
stogo/local.cc [new file with mode: 0644]
stogo/local.h [new file with mode: 0644]
stogo/prog.cc [new file with mode: 0644]
stogo/rosen.h [new file with mode: 0644]
stogo/stogo.cc [new file with mode: 0644]
stogo/stogo.h [new file with mode: 0644]
stogo/stogo_config.h [new file with mode: 0644]
stogo/testfun.h [new file with mode: 0644]
stogo/testros.cc [new file with mode: 0644]
stogo/tools.cc [new file with mode: 0644]
stogo/tools.h [new file with mode: 0644]
stogo/tst.cc [new file with mode: 0644]
stogo/tstc.c [new file with mode: 0644]

diff --git a/direct/AUTHORS b/direct/AUTHORS
new file mode 100644 (file)
index 0000000..b676a64
--- /dev/null
@@ -0,0 +1,2 @@
+C conversion: Steven G. Johnson (stevenj@alum.mit.edu)
+Original Fortran code: Joerg.M.Gablonsky (jmgablon@mailandnews.com)
diff --git a/direct/COPYRIGHT b/direct/COPYRIGHT
new file mode 100644 (file)
index 0000000..4b4a9c3
--- /dev/null
@@ -0,0 +1,28 @@
+This code is based on the DIRECT 2.0.4 Fortran code by Gablonsky et al. at
+       http://www4.ncsu.edu/~ctk/SOFTWARE/DIRECTv204.tar.gz
+The C version was initially converted via f2c and then cleaned up and
+reorganized by Steven G. Johnson (stevenj@alum.mit.edu), August 2007.
+
+******** Copyright and license for the original Fortran DIRECT code ********
+Copyright (c) 1999, 2000, 2001 North Carolina State University
+
+This program is distributed under the MIT License (see
+http://www.opensource.org/licenses/mit-license.php):
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/direct/DIRect.c b/direct/DIRect.c
new file mode 100644 (file)
index 0000000..f208332
--- /dev/null
@@ -0,0 +1,743 @@
+/* DIRect.f -- translated by f2c (version 20050501).
+   
+   f2c output hand-cleaned by SGJ (August 2007). 
+*/
+
+#include <math.h>
+#include "direct-internal.h"
+
+/* Common Block Declarations */
+
+/* Table of constant values */
+
+static integer c__90000 = 90000;
+static integer c__600 = 600;
+static integer c__64 = 64;
+static integer c__3000 = 3000;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Program       : Direct.f                                              | */
+/* | Last modified : 07-16-2001                                            | */
+/* | Written by    : Joerg Gablonsky (jmgablon@unity.ncsu.edu)             | */
+/* |                 North Carolina State University                       | */
+/* |                 Dept. of Mathematics                                  | */
+/* | DIRECT is a method to solve problems of the form:                     | */
+/* |              min f: Q --> R,                                          | */
+/* | where f is the function to be minimized and Q is an n-dimensional     | */
+/* | hyperrectangle given by the the following equation:                   | */
+/* |                                                                       | */
+/* |       Q={ x : l(i) <= x(i) <= u(i), i = 1,...,n }.                    | */
+/* | Note: This version of DIRECT can also handle hidden constraints. By   | */
+/* |       this we mean that the function may not be defined over the whole| */
+/* |       hyperrectangle Q, but only over a subset, which is not given    | */
+/* |       analytically.                                                   | */
+/* |                                                                       | */
+/* | We now give a brief outline of the algorithm:                         | */
+/* |                                                                       | */
+/* |   The algorithm starts with mapping the hyperrectangle Q to the       | */
+/* |   n-dimensional unit hypercube. DIRECT then samples the function at   | */
+/* |   the center of this hypercube and at 2n more points, 2 in each       | */
+/* |   coordinate direction. Uisng these function values, DIRECT then      | */
+/* |   divides the domain into hyperrectangles, each having exactly one of | */
+/* |   the sampling points as its center. In each iteration, DIRECT chooses| */
+/* |   some of the existing hyperrectangles to be further divided.         | */
+/* |   We provide two different strategies of how to decide which          | */
+/* |   hyperrectangles DIRECT divides and several different convergence    | */
+/* |   criteria.                                                           | */
+/* |                                                                       | */
+/* |   DIRECT was designed to solve problems where the function f is       | */
+/* |   Lipschitz continues. However, DIRECT has proven to be effective on  | */
+/* |   more complex problems than these.                                   | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_direct_(fp fcn, doublereal *x, integer *n, doublereal *eps, doublereal epsabs, integer *maxf, integer *maxt, doublereal *fmin, doublereal *l, 
+       doublereal *u, integer *algmethod, integer *ierror, FILE *logfile, 
+       doublereal *fglobal, doublereal *fglper, doublereal *volper, 
+       doublereal *sigmaper, void *fcn_data)
+{
+    /* System generated locals */
+    integer i__1, i__2;
+    doublereal d__1;
+
+    /* Local variables */
+    integer increase;
+    doublereal *c__ = 0        /* was [90000][64] */, *f = 0   /* 
+           was [90000][2] */;
+    integer i__, j, *s = 0     /* was [3000][2] */, t;
+    doublereal *w = 0;
+    doublereal divfactor;
+    integer ifeasiblef, iepschange, actmaxdeep;
+    integer actdeep_div__, iinfesiblef;
+    integer pos1, newtosample;
+    integer ifree, help;
+    doublereal *oldl = 0, fmax;
+    integer maxi;
+    doublereal kmax, *oldu = 0;
+    integer oops, *list2 = 0   /* was [64][2] */, cheat;
+    doublereal delta;
+    integer mdeep, *point = 0, start;
+    integer *anchor = 0, *length = 0   /* was [90000][64] */, *arrayi = 0;
+    doublereal *levels = 0, *thirds = 0;
+    integer writed;
+    doublereal epsfix;
+    integer oldpos, minpos, maxpos, tstart, actdeep, ifreeold, oldmaxf;
+    integer numfunc, version;
+    integer jones;
+
+    /* FIXME: change sizes dynamically? */
+#define MY_ALLOC(p, t, n) p = (t *) malloc(sizeof(t) * (n)); \
+                          if (!(p)) { *ierror = -100; goto cleanup; }
+    MY_ALLOC(c__, doublereal, 5760000);
+    MY_ALLOC(f, doublereal, 180000);
+    MY_ALLOC(s, integer, 6000);
+    MY_ALLOC(w, doublereal, 64);
+    MY_ALLOC(oldl, doublereal, 64);
+    MY_ALLOC(oldu, doublereal, 64);
+    MY_ALLOC(list2, integer, 128);
+    MY_ALLOC(point, integer, 90000);
+    MY_ALLOC(anchor, integer, 602);
+    MY_ALLOC(length, integer, 5760000);
+    MY_ALLOC(arrayi, integer, 64);
+    MY_ALLOC(levels, doublereal, 601);
+    MY_ALLOC(thirds, doublereal, 601);    
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE Direct                                                  | */
+/* | On entry                                                              | */
+/* |     fcn -- The argument containing the name of the user-supplied      | */
+/* |            SUBROUTINE that returns values for the function to be      | */
+/* |            minimized.                                                 | */
+/* |       n -- The dimension of the problem.                              | */
+/* |     eps -- Exceeding value. If eps > 0, we use the same epsilon for   | */
+/* |            all iterations. If eps < 0, we use the update formula from | */
+/* |            Jones:                                                     | */
+/* |               eps = max(1.D-4*abs(fmin),epsfix),                      | */
+/* |            where epsfix = abs(eps), the absolute value of eps which is| */
+/* |            passed to the function.                                    | */
+/* |    maxf -- The maximum number of function evaluations.                | */
+/* |    maxT -- The maximum number of iterations.                          | */
+/* |            Direct stops when either the maximum number of iterations  | */
+/* |            is reached or more than maxf function-evalutions were made.| */
+/* |       l -- The lower bounds of the hyperbox.                          | */
+/* |       u -- The upper bounds of the hyperbox.                          | */
+/* |algmethod-- Choose the method, that is either use the original method  | */
+/* |            as described by Jones et.al. (0) or use our modification(1)| */
+/* | logfile -- File-Handle for the logfile. DIRECT expects this file to be| */
+/* |            opened and closed by the user outside of DIRECT. We moved  | */
+/* |            this to the outside so the user can add extra informations | */
+/* |            to this file before and after the call to DIRECT.          | */
+/* | fglobal -- Function value of the global optimum. If this value is not | */
+/* |            known (that is, we solve a real problem, not a testproblem)| */
+/* |            set this value to -1.D100 and fglper (see below) to 0.D0.  | */
+/* |  fglper -- Terminate the optimization when the percent error          | */
+/* |                100(f_min - fglobal)/max(1,abs(fglobal)) < fglper.     | */
+/* |  volper -- Terminate the optimization when the volume of the          | */
+/* |            hyperrectangle S with f(c(S)) = fmin is less then volper   | */
+/* |            percent of the volume of the original hyperrectangle.      | */
+/* |sigmaper -- Terminate the optimization when the measure of the         | */
+/* |            hyperrectangle S with f(c(S)) = fmin is less then sigmaper.| */
+/* |                                                                       | */
+/* | User data that is passed through without being changed:               | */
+/* |  fcn_data - opaque pointer to any user data                           | */
+/* |                                                                       | */
+/* | On return                                                             | */
+/* |                                                                       | */
+/* |       x -- The final point obtained in the optimization process.      | */
+/* |            X should be a good approximation to the global minimum     | */
+/* |            for the function within the hyper-box.                     | */
+/* |                                                                       | */
+/* |    fmin -- The value of the function at x.                            | */
+/* |  Ierror -- Error flag. If Ierror is lower 0, an error has occured. The| */
+/* |            values of Ierror mean                                      | */
+/* |            Fatal errors :                                             | */
+/* |             -1   u(i) <= l(i) for some i.                             | */
+/* |             -2   maxf is too large.                                   | */
+/* |             -3   Initialization in DIRpreprc failed.                  | */
+/* |             -4   Error in DIRSamplepoints, that is there was an error | */
+/* |                  in the creation of the sample points.                | */
+/* |             -5   Error in DIRSamplef, that is an error occured while  | */
+/* |                  the function was sampled.                            | */
+/* |             -6   Error in DIRDoubleInsert, that is an error occured   | */
+/* |                  DIRECT tried to add all hyperrectangles with the same| */
+/* |                  size and function value at the center. Either        | */
+/* |                  increase maxdiv or use our modification (Jones = 1). | */
+/* |            Termination values :                                       | */
+/* |              1   Number of function evaluations done is larger then   | */
+/* |                  maxf.                                                | */
+/* |              2   Number of iterations is equal to maxT.               | */
+/* |              3   The best function value found is within fglper of    | */
+/* |                  the (known) global optimum, that is                  | */
+/* |                   100(fmin - fglobal/max(1,|fglobal|))  < fglper.     | */
+/* |                  Note that this termination signal only occurs when   | */
+/* |                  the global optimal value is known, that is, a test   | */
+/* |                  function is optimized.                               | */
+/* |              4   The volume of the hyperrectangle with fmin at its    | */
+/* |                  center is less than volper percent of the volume of  | */
+/* |                  the original hyperrectangle.                         | */
+/* |              5   The measure of the hyperrectangle with fmin at its   | */
+/* |                  center is less than sigmaper.                        | */
+/* |                                                                       | */
+/* | SUBROUTINEs used :                                                    | */
+/* |                                                                       | */
+/* | DIRheader, DIRInitSpecific, DIRInitList, DIRpreprc, DIRInit, DIRChoose| */
+/* | DIRDoubleInsert, DIRGet_I, DIRSamplepoints, DIRSamplef, DIRDivide     | */
+/* | DIRInsertList, DIRreplaceInf, DIRWritehistbox, DIRsummary, Findareas  | */
+/* |                                                                       | */
+/* | Functions used :                                                      | */
+/* |                                                                       | */
+/* | DIRgetMaxdeep, DIRgetlevel                                            | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Parameters                                                            | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | The maximum of function evaluations allowed.                          | */
+/* | The maximum dept of the algorithm.                                    | */
+/* | The maximum number of divisions allowed.                              | */
+/* | The maximal dimension of the problem.                                 | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Global Variables.                                                     | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | EXTERNAL Variables.                                                   | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | User Variables.                                                       | */
+/* | These can be used to pass user defined data to the function to be     | */
+/* | optimized.                                                            | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Place to define, if needed, some application-specific variables.      | */
+/* | Note: You should try to use the arrays defined above for this.        | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | End of application - specific variables !                             | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Internal variables :                                                  | */
+/* |       f -- values of functions.                                       | */
+/* |divfactor-- Factor used for termination with known global minimum.     | */
+/* |  anchor -- anchors of lists with deepness i, -1 is anchor for list of | */
+/* |            NaN - values.                                              | */
+/* |       S -- List of potentially optimal points.                        | */
+/* |   point -- lists                                                      | */
+/* |    ifree -- first free position                                        | */
+/* |       c -- midpoints of arrays                                        | */
+/* |  thirds -- Precalculated values of 1/3^i.                             | */
+/* |  levels -- Length of intervals.                                       | */
+/* |  length -- Length of intervall (index)                                | */
+/* |       t -- actual iteration                                           | */
+/* |       j -- loop-variable                                              | */
+/* | actdeep -- the actual minimal interval-length index                   | */
+/* |  Minpos -- position of the actual minimum                             | */
+/* |    file -- The filehandle for a datafile.                             | */
+/* |  maxpos -- The number of intervalls, which are truncated.             | */
+/* |    help -- A help variable.                                           | */
+/* | numfunc -- The actual number of function evaluations.                 | */
+/* |   file2 -- The filehandle for an other datafile.                      | */
+/* |  ArrayI -- Array with the indexes of the sides with maximum length.   | */
+/* |    maxi -- Number of directions with maximal side length.             | */
+/* |    oops -- Flag which shows if anything went wrong in the             | */
+/* |            initialisation.                                            | */
+/* |   cheat -- Obsolete. If equal 1, we don't allow Ktilde > kmax.        | */
+/* |  writed -- If writed=1, store final division to plot with Matlab.     | */
+/* |   List2 -- List of indicies of intervalls, which are to be truncated. | */
+/* |       i -- Another loop-variable.                                     | */
+/* |actmaxdeep-- The actual maximum (minimum) of possible Interval length. | */
+/* |  oldpos -- The old index of the minimum. Used to print only, if there | */
+/* |            is a new minimum found.                                    | */
+/* |  tstart -- The start of the outer loop.                               | */
+/* |   start -- The postion of the starting point in the inner loop.       | */
+/* | Newtosample -- The total number of points to sample in the inner loop.| */
+/* |       w -- Array used to divide the intervalls                        | */
+/* |    kmax -- Obsolete. If cheat = 1, Ktilde was not allowed to be larger| */
+/* |            than kmax. If Ktilde > kmax, we set ktilde = kmax.         | */
+/* |   delta -- The distance to new points from center of old hyperrec.    | */
+/* |    pos1 -- Help variable used as an index.                            | */
+/* | version -- Store the version number of DIRECT.                        | */
+/* | oldmaxf -- Store the original function budget.                        | */
+/* |increase -- Flag used to keep track if function budget was increased   | */
+/* |            because no feasible point was found.                       | */
+/* | ifreeold -- Keep track which index was free before. Used with          | */
+/* |            SUBROUTINE DIRReplaceInf.                                  | */
+/* |actdeep_div-- Keep track of the current depths for divisions.          | */
+/* |    oldl -- Array used to store the original bounds of the domain.     | */
+/* |    oldu -- Array used to store the original bounds of the domain.     | */
+/* |  epsfix -- If eps < 0, we use Jones update formula. epsfix stores the | */
+/* |            absolute value of epsilon.                                 | */
+/* |iepschange-- flag iepschange to store if epsilon stays fixed or is     | */
+/* |             changed.                                                  | */
+/* |DIRgetMaxdeep-- Function to calculate the level of a hyperrectangle.   | */
+/* |DIRgetlevel-- Function to calculate the level and stage of a hyperrec. | */
+/* |    fmax -- Keep track of the maximum value of the function found.     | */
+/* |Ifeasiblef-- Keep track if a feasible point has  been found so far.    | */
+/* |             Ifeasiblef = 0 means a feasible point has been found,     | */
+/* |             Ifeasiblef = 1 no feasible point has been found.          | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 09/25/00 Version counter.                                          | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 09/24/00 Add another actdeep to keep track of the current depths   | */
+/* |             for divisions.                                            | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* |JG 01/13/01 Added epsfix for epsilon update. If eps < 0, we use Jones  | */
+/* |            update formula. epsfix stores the absolute value of epsilon| */
+/* |            then. Also added flag iepschange to store if epsilon stays | */
+/* |            fixed or is changed.                                       | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 fmax is used to keep track of the maximum value found.    | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Ifeasiblef is used to keep track if a feasible point has  | */
+/* |             been found so far. Ifeasiblef = 0 means a feasible point  | */
+/* |             has been found, Ifeasiblef = 1 if not.                    | */
+/* | JG 03/09/01 IInfeasible is used to keep track if an infeasible point  | */
+/* |             has been found. IInfeasible > 0 means a infeasible point  | */
+/* |             has been found, IInfeasible = 0 if not.                   | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* |                            Start of code.                             | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+    /* Parameter adjustments */
+    --u;
+    --l;
+    --x;
+
+    /* Function Body */
+    writed = 0;
+    jones = *algmethod;
+/* +-----------------------------------------------------------------------+ */
+/* | Save the upper and lower bounds.                                      | */
+/* +-----------------------------------------------------------------------+ */
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       oldu[i__ - 1] = u[i__];
+       oldl[i__ - 1] = l[i__];
+/* L150: */
+    }
+/* +-----------------------------------------------------------------------+ */
+/* | Set version.                                                          | */
+/* +-----------------------------------------------------------------------+ */
+    version = 204;
+/* +-----------------------------------------------------------------------+ */
+/* | Set parameters.                                                       | */
+/* |    If cheat > 0, we do not allow \tilde{K} to be larger than kmax, and| */
+/* |    set \tilde{K} to set value if necessary. Not used anymore.         | */
+/* +-----------------------------------------------------------------------+ */
+    cheat = 0;
+    kmax = 1e10;
+    mdeep = 600;
+/* +-----------------------------------------------------------------------+ */
+/* | Write the header of the logfile.                                      | */
+/* +-----------------------------------------------------------------------+ */
+    direct_dirheader_(logfile, &version, &x[1], n, eps, maxf, maxt, &l[1], &u[1], 
+           algmethod, &c__90000, &c__600, fglobal, fglper, ierror, &epsfix, &
+                     iepschange, volper, sigmaper);
+/* +-----------------------------------------------------------------------+ */
+/* | If an error has occured while writing the header (we do some checking | */
+/* | of variables there), return to the main program.                      | */
+/* +-----------------------------------------------------------------------+ */
+    if (*ierror < 0) {
+       goto cleanup;
+    }
+/* +-----------------------------------------------------------------------+ */
+/* | If the known global minimum is equal 0, we cannot divide by it.       | */
+/* | Therefore we set it to 1. If not, we set the divisionfactor to the    | */
+/* | absolute value of the global minimum.                                 | */
+/* +-----------------------------------------------------------------------+ */
+    if (*fglobal == 0.) {
+       divfactor = 1.;
+    } else {
+       divfactor = fabs(*fglobal);
+    }
+/* +-----------------------------------------------------------------------+ */
+/* | Save the budget given by the user. The variable maxf will be changed  | */
+/* | if in the beginning no feasible points are found.                     | */
+/* +-----------------------------------------------------------------------+ */
+    oldmaxf = *maxf;
+    increase = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | Initialiase the lists.                                                | */
+/* +-----------------------------------------------------------------------+ */
+    direct_dirinitlist_(anchor, &ifree, point, f, &c__90000, &c__600);
+/* +-----------------------------------------------------------------------+ */
+/* | Call the routine to initialise the mapping of x from the n-dimensional| */
+/* | unit cube to the hypercube given by u and l. If an error occured,     | */
+/* | give out a error message and return to the main program with the error| */
+/* | flag set.                                                             | */
+/* | JG 07/16/01 Changed call to remove unused data.                       | */
+/* +-----------------------------------------------------------------------+ */
+    direct_dirpreprc_(&u[1], &l[1], n, &l[1], &u[1], &oops);
+    if (oops > 0) {
+       if (logfile)
+            fprintf(logfile,"WARNING: Initialization in DIRpreprc failed.\n");
+       *ierror = -3;
+       goto cleanup;
+    }
+    tstart = 2;
+/* +-----------------------------------------------------------------------+ */
+/* | Initialise the algorithm DIRECT.                                      | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Added variable to keep track of the maximum value found.              | */
+/* +-----------------------------------------------------------------------+ */
+    direct_dirinit_(f, fcn, c__, length, &actdeep, point, anchor, &ifree,
+           logfile, arrayi, &maxi, list2, w, &x[1], &l[1], &u[1], 
+           fmin, &minpos, thirds, levels, &c__90000, &c__600, n, &c__64, &
+           fmax, &ifeasiblef, &iinfesiblef, ierror, fcn_data, jones);
+/* +-----------------------------------------------------------------------+ */
+/* | Added error checking.                                                 | */
+/* +-----------------------------------------------------------------------+ */
+    if (*ierror < 0) {
+       if (*ierror == -4) {
+           if (logfile)
+                fprintf(logfile, "WARNING: Error occured in routine DIRsamplepoints.\n");
+           goto cleanup;
+       }
+       if (*ierror == -5) {
+           if (logfile)
+                fprintf(logfile, "WARNING: Error occured in routine DIRsamplef..\n");
+           goto cleanup;
+       }
+    }
+    numfunc = maxi + 1 + maxi;
+    actmaxdeep = 1;
+    oldpos = 0;
+    tstart = 2;
+/* +-----------------------------------------------------------------------+ */
+/* | If no feasible point has been found, give out the iteration, the      | */
+/* | number of function evaluations and a warning. Otherwise, give out     | */
+/* | the iteration, the number of function evaluations done and fmin.      | */
+/* +-----------------------------------------------------------------------+ */
+    if (ifeasiblef > 0) {
+        if (logfile)
+             fprintf(logfile, "No feasible point found in %d iterations "
+                     "and %d function evaluations.\n", tstart-1, numfunc);
+    } else {
+        if (logfile)
+             fprintf(logfile, "%d, %d, %g, %g\n", 
+                     tstart-1, numfunc, *fmin, fmax);
+    }
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Main loop!                                                            | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+    i__1 = *maxt;
+    for (t = tstart; t <= i__1; ++t) {
+/* +-----------------------------------------------------------------------+ */
+/* | Choose the sample points. The indices of the sample points are stored | */
+/* | in the list S.                                                        | */
+/* +-----------------------------------------------------------------------+ */
+       actdeep = actmaxdeep;
+       direct_dirchoose_(anchor, s, &c__600, f, fmin, *eps, epsabs, levels, &maxpos, length, 
+               &c__90000, &c__600, &c__3000, n, logfile, &cheat, &
+               kmax, &ifeasiblef, jones);
+/* +-----------------------------------------------------------------------+ */
+/* | Add other hyperrectangles to S, which have the same level and the same| */
+/* | function value at the center as the ones found above (that are stored | */
+/* | in S). This is only done if we use the original DIRECT algorithm.     | */
+/* | JG 07/16/01 Added Errorflag.                                          | */
+/* +-----------------------------------------------------------------------+ */
+       if (*algmethod == 0) {
+            direct_dirdoubleinsert_(anchor, s, &maxpos, point, f, &c__600, &c__90000,
+                    &c__3000, ierror);
+           if (*ierror == -6) {
+               if (logfile)
+                    fprintf(logfile,
+"WARNING: Capacity of array S in DIRDoubleInsert reached. Increase maxdiv.\n"
+"This means that there are a lot of hyperrectangles with the same function\n"
+"value at the center. We suggest to use our modification instead (Jones = 1)\n"
+                         );
+               goto cleanup;
+           }
+       }
+       oldpos = minpos;
+/* +-----------------------------------------------------------------------+ */
+/* | Initialise the number of sample points in this outer loop.            | */
+/* +-----------------------------------------------------------------------+ */
+       newtosample = 0;
+       i__2 = maxpos;
+       for (j = 1; j <= i__2; ++j) {
+           actdeep = s[j + 2999];
+/* +-----------------------------------------------------------------------+ */
+/* | If the actual index is a point to sample, do it.                      | */
+/* +-----------------------------------------------------------------------+ */
+           if (s[j - 1] > 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | JG 09/24/00 Calculate the value delta used for sampling points.       | */
+/* +-----------------------------------------------------------------------+ */
+               actdeep_div__ = direct_dirgetmaxdeep_(&s[j - 1], length, &c__90000, 
+                       n);
+               delta = thirds[actdeep_div__ + 1];
+               actdeep = s[j + 2999];
+/* +-----------------------------------------------------------------------+ */
+/* | If the current dept of division is only one under the maximal allowed | */
+/* | dept, stop the computation.                                           | */
+/* +-----------------------------------------------------------------------+ */
+               if (actdeep + 1 >= mdeep) {
+                   if (logfile)
+                        fprintf(logfile, "WARNING: Maximum number of levels reached. Increase maxdeep.\n");
+                   *ierror = -6;
+                   goto L100;
+               }
+               actmaxdeep = MAX(actdeep,actmaxdeep);
+               help = s[j - 1];
+               if (! (anchor[actdeep + 1] == help)) {
+                   pos1 = anchor[actdeep + 1];
+                   while(! (point[pos1 - 1] == help)) {
+                       pos1 = point[pos1 - 1];
+                   }
+                   point[pos1 - 1] = point[help - 1];
+               } else {
+                   anchor[actdeep + 1] = point[help - 1];
+               }
+               if (actdeep < 0) {
+                   actdeep = (integer) f[help - 1];
+               }
+/* +-----------------------------------------------------------------------+ */
+/* | Get the Directions in which to decrease the intervall-length.         | */
+/* +-----------------------------------------------------------------------+ */
+               direct_dirget_i__(length, &help, arrayi, &maxi, n, &c__90000);
+/* +-----------------------------------------------------------------------+ */
+/* | Sample the function. To do this, we first calculate the points where  | */
+/* | we need to sample the function. After checking for errors, we then do | */
+/* | the actual evaluation of the function, again followed by checking for | */
+/* | errors.                                                               | */
+/* +-----------------------------------------------------------------------+ */
+               direct_dirsamplepoints_(c__, arrayi, &delta, &help, &start, length, 
+                       logfile, f, &ifree, &maxi, point, &x[
+                       1], &l[1], fmin, &minpos, &u[1], n, &c__90000, &
+                       c__600, &oops);
+               if (oops > 0) {
+                   if (logfile)
+                        fprintf(logfile, "WARNING: Error occured in routine DIRsamplepoints.\n");
+                   *ierror = -4;
+                   goto cleanup;
+               }
+               newtosample += maxi;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
+/* +-----------------------------------------------------------------------+ */
+               direct_dirsamplef_(c__, arrayi, &delta, &help, &start, length,
+                           logfile, f, &ifree, &maxi, point, fcn, &x[
+                       1], &l[1], fmin, &minpos, &u[1], n, &c__90000, &
+                       c__600, &oops, &fmax, &ifeasiblef, &iinfesiblef, 
+                       fcn_data);
+               if (oops > 0) {
+                   if (logfile)
+                        fprintf(logfile, "WARNING: Error occured in routine DIRsamplef.\n");
+                   *ierror = -5;
+                   goto cleanup;
+               }
+/* +-----------------------------------------------------------------------+ */
+/* | Divide the intervalls.                                                | */
+/* +-----------------------------------------------------------------------+ */
+               direct_dirdivide_(&start, &actdeep_div__, length, point, arrayi, &
+                       help, list2, w, &maxi, f, &c__90000, &c__600, n);
+/* +-----------------------------------------------------------------------+ */
+/* | Insert the new intervalls into the list (sorted).                     | */
+/* +-----------------------------------------------------------------------+ */
+               direct_dirinsertlist_(&start, anchor, point, f, &maxi, length, &
+                       c__90000, &c__600, n, &help, jones);
+/* +-----------------------------------------------------------------------+ */
+/* | Increase the number of function evaluations.                          | */
+/* +-----------------------------------------------------------------------+ */
+               numfunc = numfunc + maxi + maxi;
+           }
+/* +-----------------------------------------------------------------------+ */
+/* | End of main loop.                                                     | */
+/* +-----------------------------------------------------------------------+ */
+/* L20: */
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | If there is a new minimum, show the actual iteration, the number of   | */
+/* | function evaluations, the minimum value of f (so far) and the position| */
+/* | in the array.                                                         | */
+/* +-----------------------------------------------------------------------+ */
+       if (oldpos < minpos) {
+           if (logfile)
+                fprintf(logfile, "%d, %d, %g, %g\n",
+                        t, numfunc, *fmin, fmax);
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | If no feasible point has been found, give out the iteration, the      | */
+/* | number of function evaluations and a warning.                         | */
+/* +-----------------------------------------------------------------------+ */
+       if (ifeasiblef > 0) {
+           if (logfile)
+                fprintf(logfile, "No feasible point found in %d iterations "
+                        "and %d function evaluations\n", t, numfunc);
+       }
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* |                       Termination Checks                              | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Calculate the index for the hyperrectangle at which       | */
+/* |             fmin is assumed. We then calculate the volume of this     | */
+/* |             hyperrectangle and store it in delta. This delta can be   | */
+/* |             used to stop DIRECT once the volume is below a certain    | */
+/* |             percentage of the original volume. Since the original     | */
+/* |             is 1 (scaled), we can stop once delta is below a certain  | */
+/* |             percentage, given by volper.                              | */
+/* +-----------------------------------------------------------------------+ */
+       *ierror = jones;
+       jones = 0;
+       actdeep_div__ = direct_dirgetlevel_(&minpos, length, &c__90000, n, jones);
+       jones = *ierror;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Use precalculated values to calculate volume.             | */
+/* +-----------------------------------------------------------------------+ */
+       delta = thirds[actdeep_div__] * 100;
+       if (delta <= *volper) {
+           *ierror = 4;
+           if (logfile)
+                fprintf(logfile, "DIRECT stopped: Volume of S_min is "
+                        "%g%% < %g%% of the original volume.\n",
+                        delta, *volper);
+           goto L100;
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/23/01 Calculate the measure for the hyperrectangle at which     | */
+/* |             fmin is assumed. If this measure is smaller then sigmaper,| */
+/* |             we stop DIRECT.                                           | */
+/* +-----------------------------------------------------------------------+ */
+       actdeep_div__ = direct_dirgetlevel_(&minpos, length, &c__90000, n, jones);
+       delta = levels[actdeep_div__];
+       if (delta <= *sigmaper) {
+           *ierror = 5;
+           if (logfile)
+                fprintf(logfile, "DIRECT stopped: Measure of S_min "
+                        "= %g < %g.\n", delta, *sigmaper);
+           goto L100;
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | If the best found function value is within fglper of the (known)      | */
+/* | global minimum value, terminate. This only makes sense if this optimal| */
+/* | value is known, that is, in test problems.                            | */
+/* +-----------------------------------------------------------------------+ */
+       if ((*fmin - *fglobal) * 100 / divfactor <= *fglper) {
+           *ierror = 3;
+           if (logfile)
+                fprintf(logfile, "DIRECT stopped: fmin within fglper of global minimum.\n");
+           goto L100;
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | Find out if there are infeasible points which are near feasible ones. | */
+/* | If this is the case, replace the function value at the center of the  | */
+/* | hyper rectangle by the lowest function value of a nearby function.    | */
+/* | If no infeasible points exist (IInfesiblef = 0), skip this.           | */
+/* +-----------------------------------------------------------------------+ */
+       if (iinfesiblef > 0) {
+            direct_dirreplaceinf_(&ifree, &ifreeold, f, c__, thirds, length, anchor, 
+                   point, &u[1], &l[1], &c__90000, &c__600, &c__64, n, 
+                   logfile, &fmax, jones);
+       }
+       ifreeold = ifree;
+/* +-----------------------------------------------------------------------+ */
+/* | If iepschange = 1, we use the epsilon change formula from Jones.      | */
+/* +-----------------------------------------------------------------------+ */
+       if (iepschange == 1) {
+/* Computing MAX */
+           d__1 = fabs(*fmin) * 1e-4;
+           *eps = MAX(d__1,epsfix);
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | If no feasible point has been found yet, set the maximum number of    | */
+/* | function evaluations to the number of evaluations already done plus   | */
+/* | the budget given by the user.                                         | */
+/* | If the budget has already be increased, increase it again. If a       | */
+/* | feasible point has been found, remark that and reset flag. No further | */
+/* | increase is needed.                                                   | */
+/* +-----------------------------------------------------------------------+ */
+       if (increase == 1) {
+           *maxf = numfunc + oldmaxf;
+           if (ifeasiblef == 0) {
+               if (logfile)
+                    fprintf(logfile, "DIRECT found a feasible point.  The "
+                            "adjusted budget is now set to %d.\n", *maxf);
+               increase = 0;
+           }
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | Check if the number of function evaluations done is larger than the   | */
+/* | allocated budget. If this is the case, check if a feasible point was  | */
+/* | found. If this is a case, terminate. If no feasible point was found,  | */
+/* | increase the budget and set flag increase.                            | */
+/* +-----------------------------------------------------------------------+ */
+       if (numfunc > *maxf) {
+           if (ifeasiblef == 0) {
+               *ierror = 1;
+               if (logfile)
+                    fprintf(logfile, "DIRECT stopped: numfunc >= maxf.\n");
+               goto L100;
+           } else {
+               increase = 1;
+               if (logfile)
+                     fprintf(logfile, 
+"DIRECT could not find a feasible point after %d function evaluations.\n"
+"DIRECT continues until a feasible point is found.\n", numfunc);
+               *maxf = numfunc + oldmaxf;
+           }
+       }
+/* L10: */
+    }
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | End of main loop.                                                     | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | The algorithm stopped after maxT iterations.                          | */
+/* +-----------------------------------------------------------------------+ */
+    *ierror = 2;
+    if (logfile)
+        fprintf(logfile, "DIRECT stopped: maxT iterations.\n");
+
+L100:
+/* +-----------------------------------------------------------------------+ */
+/* | Store the position of the minimum in x.                               | */
+/* +-----------------------------------------------------------------------+ */
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       x[i__] = c__[minpos + i__ * 90000 - 90001] * l[i__] + l[i__] * u[i__];
+       u[i__] = oldu[i__ - 1];
+       l[i__] = oldl[i__ - 1];
+/* L50: */
+    }
+/* +-----------------------------------------------------------------------+ */
+/* | Store the number of function evaluations in maxf.                     | */
+/* +-----------------------------------------------------------------------+ */
+    *maxf = numfunc;
+/* +-----------------------------------------------------------------------+ */
+/* | Give out a summary of the run.                                        | */
+/* +-----------------------------------------------------------------------+ */
+    direct_dirsummary_(logfile, &x[1], &l[1], &u[1], n, fmin, fglobal, &numfunc, 
+           ierror);
+/* +-----------------------------------------------------------------------+ */
+/* | Format statements.                                                    | */
+/* +-----------------------------------------------------------------------+ */
+
+ cleanup:
+#define MY_FREE(p) if (p) free(p)
+    MY_FREE(c__);
+    MY_FREE(f);
+    MY_FREE(s);
+    MY_FREE(w);
+    MY_FREE(oldl);
+    MY_FREE(oldu);
+    MY_FREE(list2);
+    MY_FREE(point);
+    MY_FREE(anchor);
+    MY_FREE(length);
+    MY_FREE(arrayi);
+    MY_FREE(levels);
+    MY_FREE(thirds);
+} /* direct_ */
+
diff --git a/direct/DIRparallel.c b/direct/DIRparallel.c
new file mode 100644 (file)
index 0000000..aba092e
--- /dev/null
@@ -0,0 +1,381 @@
+/* DIRparallel.f -- translated by f2c (version 20050501).
+
+   f2c output hand-cleaned by SGJ (August 2007).
+*/
+
+#include "direct-internal.h"
+
+/* Table of constant values */
+
+static integer c__0 = 0;
+static integer c_n1 = -1;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Program       : Direct.f (subfile DIRseriell.f)                       | */
+/* | Last modified : 02-22-01                                              | */
+/* | Written by    : Joerg Gablonsky                                       | */
+/* | Subroutines, which differ depending on the serial or parallel version.| */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Parallel Direct. This routine replaces the normal main routine DIRect.| */
+/* | In it, we find out if this pe is the master or slave. If it is the    | */
+/* | master, it calls the serial DIRect main routine. The only routine that| */
+/* | has to change for parallel Direct is DIRSamplef, where the actual     | */
+/* | sampling of the function is done. If we are on the slave, wait for    | */
+/* | either the coordinates of a point to sample the function or the       | */
+/* | termination signal.                                                   | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ int direct_pardirect_(fp fcn, doublereal *x, integer *n, 
+       doublereal *eps, integer *maxf, integer *maxt, doublereal *fmin, 
+       doublereal *l, doublereal *u, integer *algmethod, integer *ierror, 
+       FILE *logfile, doublereal *fglobal, doublereal *fglper, doublereal 
+       *volper, doublereal *sigmaper, void *fcn_data)
+{
+    /* System generated locals */
+    integer i__1;
+
+    /* Local variables */
+    integer i__, k;
+    integer tid;
+    integer flag__;
+    doublereal fval;
+    integer tids[360], kret;
+    integer mytid;
+    doublereal fscale;
+    integer nprocs;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Parameters                                                            | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | The maximum of function evaluations allowed.                          | */
+/* | The maximum dept of the algorithm.                                    | */
+/* | The maximum number of divisions allowed.                              | */
+/* | The maximal dimension of the problem.                                 | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Global Variables.                                                     | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | External Variables.                                                   | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | User Variables.                                                       | */
+/* | These can be used to pass user defined data to the function to be     | */
+/* | optimized.                                                            | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Parallel programming variables                                        | */
+/* +-----------------------------------------------------------------------+ */
+/*       maxprocs should be >= the number of processes used for DIRECT */
+/* +-----------------------------------------------------------------------+ */
+/* | End of parallel programming variables                                 | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Internal variables                                                    | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 02/28/01 Begin of parallel additions                               | */
+/* | DETERMINE MASTER PROCESSOR. GET TIDS OF ALL PROCESSORS.               | */
+/* +-----------------------------------------------------------------------+ */
+    /* Parameter adjustments */
+    --u;
+    --l;
+    --x;
+
+    /* Function Body */
+    getmytidif_(&mytid);
+    getnprocsif_(&nprocs);
+    gettidif_(&c__0, tids);
+/* +-----------------------------------------------------------------------+ */
+/* | If I am the master get the other tids and start running DIRECT.       | */
+/* | Otherwise, branch off to do function evaluations.                     | */
+/* +-----------------------------------------------------------------------+ */
+    if (mytid == tids[0]) {
+       i__1 = nprocs - 1;
+       for (i__ = 1; i__ <= i__1; ++i__) {
+           gettidif_(&i__, &tids[i__]);
+/* L46: */
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | Call Direct main routine. This routine calls DIRSamplef for the       | */
+/* | function evaluations, which are then done in parallel.                | */
+/* +-----------------------------------------------------------------------+ */
+       direct_direct_(fcn, &x[1], n, eps, maxf, maxt, fmin, &l[1], &u[1], 
+               algmethod, ierror, logfile, fglobal, fglper, volper, sigmaper,
+               fcn_data);
+/* +-----------------------------------------------------------------------+ */
+/* | Send exit message to rest of pe's.                                    | */
+/* +-----------------------------------------------------------------------+ */
+       flag__ = 0;
+       i__1 = nprocs;
+       for (tid = 2; tid <= i__1; ++tid) {
+           mastersendif_(&tids[tid - 1], &tids[tid - 1], n, &flag__, &flag__,
+                    &x[1], &u[1], &l[1], &x[1]);
+/* L200: */
+       }
+    } else {
+/* +-----------------------------------------------------------------------+ */
+/* | This is what the slaves do!!                                          | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* |   Receive the first point from the master processor.                  | */
+/* +-----------------------------------------------------------------------+ */
+       slaverecvif_(tids, &c_n1, n, &flag__, &k, &fscale, &u[1], &l[1], &x[1]
+               );
+/* +-----------------------------------------------------------------------+ */
+/* | Repeat until master signals to stop.                                  | */
+/* +-----------------------------------------------------------------------+ */
+       while(flag__ > 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | Evaluate f(x).                                                        | */
+/* +-----------------------------------------------------------------------+ */
+            direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &fval, &kret, &
+                     fcn_data);
+/* +-----------------------------------------------------------------------+ */
+/* | Send result and wait for next point / message with signal to stop.    | */
+/* +-----------------------------------------------------------------------+ */
+           slavesendif_(tids, &mytid, &k, &mytid, &fval, &kret);
+           slaverecvif_(tids, &c_n1, n, &flag__, &k, &fscale, &u[1], &l[1], &
+                   x[1]);
+       }
+    }
+    return 0;
+} /* pardirect_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | Subroutine for sampling. This sampling is done in parallel, the master| */
+/* | prozessor is also evaluating the function sometimes.                  | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirsamplef_(doublereal *c__, integer *arrayi, doublereal 
+       *delta, integer *sample, integer *new__, integer *length, 
+       FILE *logfile, doublereal *f, integer *free, integer *maxi, 
+       integer *point, fp fcn, doublereal *x, doublereal *l, doublereal *
+       fmin, integer *minpos, doublereal *u, integer *n, integer *maxfunc, 
+       integer *maxdeep, integer *oops, doublereal *fmax, integer *
+       ifeasiblef, integer *iinfesiblef, void *fcn_data)
+{
+    /* System generated locals */
+    integer length_dim1, length_offset, c_dim1, c_offset, f_dim1, f_offset, 
+           i__1;
+    doublereal d__1;
+
+    /* Local variables */
+    integer i__, j, k, helppoint, tid, pos;
+    integer flag__, tids[360], kret, npts;
+    doublereal fhelp;
+    integer oldpos, nprocs, datarec;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 fcn must be declared external.                            | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Removed fcn.                                              | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
+/* |             Added variable to keep track if feasible point was found. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized.  | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Parallel programming variables.                                       | */
+/* +-----------------------------------------------------------------------+ */
+/* JG 09/05/00 Increase this if more processors are used. */
+/* +-----------------------------------------------------------------------+ */
+/* | Find out the id's of all processors.                                  | */
+/* +-----------------------------------------------------------------------+ */
+    /* Parameter adjustments */
+    --u;
+    --l;
+    --x;
+    --arrayi;
+    --point;
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+    c_dim1 = *maxfunc;
+    c_offset = 1 + c_dim1;
+    c__ -= c_offset;
+
+    /* Function Body */
+    getnprocsif_(&nprocs);
+    i__1 = nprocs - 1;
+    for (i__ = 0; i__ <= i__1; ++i__) {
+       gettidif_(&i__, &tids[i__]);
+/* L46: */
+    }
+/* +-----------------------------------------------------------------------+ */
+/* | Set the pointer to the first function to be evaluated,                | */
+/* | store this position also in helppoint.                                | */
+/* +-----------------------------------------------------------------------+ */
+    pos = *new__;
+    helppoint = pos;
+/* +-----------------------------------------------------------------------+ */
+/* | Iterate over all points, where the function should be                 | */
+/* | evaluated.                                                            | */
+/* +-----------------------------------------------------------------------+ */
+    flag__ = 1;
+    npts = *maxi + *maxi;
+    k = 1;
+    while(k <= npts && k < nprocs) {
+/* +-----------------------------------------------------------------------+ */
+/* | tid is the id of the prozessor the next points is send to.            | */
+/* +-----------------------------------------------------------------------+ */
+       tid = k + 1;
+/* +-----------------------------------------------------------------------+ */
+/* | Copy the position into the helparray x.                               | */
+/* +-----------------------------------------------------------------------+ */
+       i__1 = *n;
+       for (i__ = 1; i__ <= i__1; ++i__) {
+           x[i__] = c__[pos + i__ * c_dim1];
+/* L60: */
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | Send the point.                                                       | */
+/* +-----------------------------------------------------------------------+ */
+       mastersendif_(&tids[tid - 1], &tids[tid - 1], n, &flag__, &pos, &x[1],
+                &u[1], &l[1], &x[1]);
+       ++k;
+       pos = point[pos];
+/* +-----------------------------------------------------------------------+ */
+/* | Get the next point.                                                   | */
+/* +-----------------------------------------------------------------------+ */
+    }
+/* +-----------------------------------------------------------------------+ */
+/* |  Get data until it is all received.                                   | */
+/* +-----------------------------------------------------------------------+ */
+    datarec = 0;
+    while(datarec < npts) {
+       if ((doublereal) datarec / (doublereal) nprocs - datarec / nprocs < 
+               1e-5 && k <= npts) {
+           i__1 = *n;
+           for (i__ = 1; i__ <= i__1; ++i__) {
+               x[i__] = c__[pos + i__ * c_dim1];
+/* L165: */
+           }
+           direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &fhelp, &kret,
+                     fcn_data);
+           oldpos = pos;
+           f[oldpos + f_dim1] = fhelp;
+           ++datarec;
+/* +-----------------------------------------------------------------------+ */
+/* | Remember if an infeasible point has been found.                       | */
+/* +-----------------------------------------------------------------------+ */
+           *iinfesiblef = MAX(*iinfesiblef,kret);
+           if (kret == 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | if the function evaluation was O.K., set the flag in                  | */
+/* | f(pos,2).                                                             | */
+/* +-----------------------------------------------------------------------+ */
+               f[oldpos + (f_dim1 << 1)] = 0.;
+               *ifeasiblef = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
+/* +-----------------------------------------------------------------------+ */
+/* Computing MAX */
+               d__1 = f[pos + f_dim1];
+               *fmax = MAX(d__1,*fmax);
+           }
+/* +-----------------------------------------------------------------------+ */
+/* | Remember if an infeasible point has been found.                       | */
+/* +-----------------------------------------------------------------------+ */
+           *iinfesiblef = MAX(*iinfesiblef,kret);
+           if (kret == 1) {
+/* +-----------------------------------------------------------------------+ */
+/* | If the function could not be evaluated at the given point,            | */
+/* | set flag to mark this (f(pos,2) and store the maximum                 | */
+/* | box-sidelength in f(pos,1).                                           | */
+/* +-----------------------------------------------------------------------+ */
+               f[oldpos + (f_dim1 << 1)] = 2.;
+               f[oldpos + f_dim1] = *fmax;
+           }
+/* +-----------------------------------------------------------------------+ */
+/* | If the function could not be evaluated due to a failure in            | */
+/* | the setup, mark this.                                                 | */
+/* +-----------------------------------------------------------------------+ */
+           if (kret == -1) {
+               f[oldpos + (f_dim1 << 1)] = -1.;
+           }
+           ++k;
+           pos = point[pos];
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | Recover where to store the value.                                     | */
+/* +-----------------------------------------------------------------------+ */
+       masterrecvif_(&c_n1, &c_n1, &oldpos, &tid, &fhelp, &kret);
+       f[oldpos + f_dim1] = fhelp;
+       ++datarec;
+/* +-----------------------------------------------------------------------+ */
+/* | Remember if an infeasible point has been found.                       | */
+/* +-----------------------------------------------------------------------+ */
+       *iinfesiblef = MAX(*iinfesiblef,kret);
+       if (kret == 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | if the function evaluation was O.K., set the flag in                  | */
+/* | f(pos,2).                                                             | */
+/* +-----------------------------------------------------------------------+ */
+           f[oldpos + (f_dim1 << 1)] = 0.;
+           *ifeasiblef = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
+/* +-----------------------------------------------------------------------+ */
+/* Computing MAX */
+           d__1 = f[oldpos + f_dim1];
+           *fmax = MAX(d__1,*fmax);
+       }
+       if (kret == 1) {
+/* +-----------------------------------------------------------------------+ */
+/* | If the function could not be evaluated at the given point,            | */
+/* | set flag to mark this (f(pos,2) and store the maximum                 | */
+/* | box-sidelength in f(pos,1).                                           | */
+/* +-----------------------------------------------------------------------+ */
+           f[oldpos + (f_dim1 << 1)] = 2.;
+           f[oldpos + f_dim1] = *fmax;
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | If the function could not be evaluated due to a failure in            | */
+/* | the setup, mark this.                                                 | */
+/* +-----------------------------------------------------------------------+ */
+       if (kret == -1) {
+           f[oldpos + (f_dim1 << 1)] = -1.;
+       }
+/* +-----------------------------------------------------------------------+ */
+/* |         Send data until it is all sent.                               | */
+/* +-----------------------------------------------------------------------+ */
+       if (k <= npts) {
+/* +-----------------------------------------------------------------------+ */
+/* | Copy the position into the helparray x.                               | */
+/* +-----------------------------------------------------------------------+ */
+           i__1 = *n;
+           for (i__ = 1; i__ <= i__1; ++i__) {
+               x[i__] = c__[pos + i__ * c_dim1];
+/* L160: */
+           }
+           mastersendif_(&tid, &tid, n, &flag__, &pos, &x[1], &u[1], &l[1], &
+                   x[1]);
+           ++k;
+           pos = point[pos];
+       }
+    }
+    pos = helppoint;
+/* +-----------------------------------------------------------------------+ */
+/* | Iterate over all evaluated points and see, if the minimal             | */
+/* | value of the function has changed. If this has happend,               | */
+/* | store the minimal value and its position in the array.                | */
+/* | Attention: Only valied values are checked!!                           | */
+/* +-----------------------------------------------------------------------+ */
+    i__1 = *maxi + *maxi;
+    for (j = 1; j <= i__1; ++j) {
+       if (f[pos + f_dim1] < *fmin && f[pos + (f_dim1 << 1)] == 0.) {
+           *fmin = f[pos + f_dim1];
+           *minpos = pos;
+       }
+       pos = point[pos];
+/* L50: */
+    }
+} /* dirsamplef_ */
diff --git a/direct/DIRserial.c b/direct/DIRserial.c
new file mode 100644 (file)
index 0000000..bc45a51
--- /dev/null
@@ -0,0 +1,145 @@
+/* DIRserial.f -- translated by f2c (version 20050501).
+
+   f2c output hand-cleaned by SGJ (August 2007).
+*/
+
+#include "direct-internal.h"
+
+/* +-----------------------------------------------------------------------+ */
+/* | Program       : Direct.f (subfile DIRserial.f)                        | */
+/* | Last modified : 04-12-2001                                            | */
+/* | Written by    : Joerg Gablonsky                                       | */
+/* | SUBROUTINEs, which differ depENDing on the serial or parallel version.| */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE for sampling.                                              | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirsamplef_(doublereal *c__, integer *arrayi, doublereal 
+       *delta, integer *sample, integer *new__, integer *length, 
+       FILE *logfile, doublereal *f, integer *free, integer *maxi, 
+       integer *point, fp fcn, doublereal *x, doublereal *l, doublereal *
+       fmin, integer *minpos, doublereal *u, integer *n, integer *maxfunc, 
+       integer *maxdeep, integer *oops, doublereal *fmax, integer *
+       ifeasiblef, integer *iinfesiblef, void *fcn_data)
+{
+    /* System generated locals */
+    integer length_dim1, length_offset, c_dim1, c_offset, f_dim1, f_offset, 
+           i__1, i__2;
+    doublereal d__1;
+
+    /* Local variables */
+    integer i__, j, helppoint, pos, kret;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 fcn must be declared external.                            | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Removed fcn.                                              | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
+/* |             Added variable to keep track IF feasible point was found. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized.  | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Set the pointer to the first function to be evaluated,                | */
+/* | store this position also in helppoint.                                | */
+/* +-----------------------------------------------------------------------+ */
+    /* Parameter adjustments */
+    --u;
+    --l;
+    --x;
+    --arrayi;
+    --point;
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+    c_dim1 = *maxfunc;
+    c_offset = 1 + c_dim1;
+    c__ -= c_offset;
+
+    /* Function Body */
+    pos = *new__;
+    helppoint = pos;
+/* +-----------------------------------------------------------------------+ */
+/* | Iterate over all points, where the function should be                 | */
+/* | evaluated.                                                            | */
+/* +-----------------------------------------------------------------------+ */
+    i__1 = *maxi + *maxi;
+    for (j = 1; j <= i__1; ++j) {
+/* +-----------------------------------------------------------------------+ */
+/* | Copy the position into the helparrayy x.                              | */
+/* +-----------------------------------------------------------------------+ */
+       i__2 = *n;
+       for (i__ = 1; i__ <= i__2; ++i__) {
+           x[i__] = c__[pos + i__ * c_dim1];
+/* L60: */
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | Call the function.                                                    | */
+/* +-----------------------------------------------------------------------+ */
+       direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &f[pos + f_dim1], &kret, 
+                 fcn_data);
+/* +-----------------------------------------------------------------------+ */
+/* | Remember IF an infeasible point has been found.                       | */
+/* +-----------------------------------------------------------------------+ */
+       *iinfesiblef = MAX(*iinfesiblef,kret);
+       if (kret == 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | IF the function evaluation was O.K., set the flag in                  | */
+/* | f(pos,2). Also mark that a feasible point has been found.             | */
+/* +-----------------------------------------------------------------------+ */
+           f[pos + (f_dim1 << 1)] = 0.;
+           *ifeasiblef = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
+/* +-----------------------------------------------------------------------+ */
+/* Computing MAX */
+           d__1 = f[pos + f_dim1];
+           *fmax = MAX(d__1,*fmax);
+       }
+       if (kret >= 1) {
+/* +-----------------------------------------------------------------------+ */
+/* |  IF the function could not be evaluated at the given point,            | */
+/* | set flag to mark this (f(pos,2) and store the maximum                 | */
+/* | box-sidelength in f(pos,1).                                           | */
+/* +-----------------------------------------------------------------------+ */
+           f[pos + (f_dim1 << 1)] = 2.;
+           f[pos + f_dim1] = *fmax;
+       }
+/* +-----------------------------------------------------------------------+ */
+/* |  IF the function could not be evaluated due to a failure in            | */
+/* | the setup, mark this.                                                 | */
+/* +-----------------------------------------------------------------------+ */
+       if (kret == -1) {
+           f[pos + (f_dim1 << 1)] = -1.;
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | Set the position to the next point, at which the function             | */
+/* | should e evaluated.                                                   | */
+/* +-----------------------------------------------------------------------+ */
+       pos = point[pos];
+/* L40: */
+    }
+    pos = helppoint;
+/* +-----------------------------------------------------------------------+ */
+/* | Iterate over all evaluated points and see, IF the minimal             | */
+/* | value of the function has changed.  IF this has happEND,               | */
+/* | store the minimal value and its position in the array.                | */
+/* | Attention: Only valied values are checked!!                           | */
+/* +-----------------------------------------------------------------------+ */
+    i__1 = *maxi + *maxi;
+    for (j = 1; j <= i__1; ++j) {
+       if (f[pos + f_dim1] < *fmin && f[pos + (f_dim1 << 1)] == 0.) {
+           *fmin = f[pos + f_dim1];
+           *minpos = pos;
+       }
+       pos = point[pos];
+/* L50: */
+    }
+} /* dirsamplef_ */
diff --git a/direct/DIRsubrout.c b/direct/DIRsubrout.c
new file mode 100644 (file)
index 0000000..4d8a45e
--- /dev/null
@@ -0,0 +1,1715 @@
+/* DIRsubrout.f -- translated by f2c (version 20050501).
+
+   f2c output hand-cleaned by SGJ (August 2007).
+*/
+
+#include <math.h>
+#include "direct-internal.h"
+
+/* Table of constant values */
+
+static integer c__1 = 1;
+static integer c__32 = 32;
+static integer c__0 = 0;
+
+/* +-----------------------------------------------------------------------+ */
+/* | INTEGER Function DIRGetlevel                                          | */
+/* | Returns the level of the hyperrectangle. Depending on the value of the| */
+/* | global variable JONES. IF JONES equals 0, the level is given by       | */
+/* |               kN + p, where the rectangle has p sides with a length of| */
+/* |             1/3^(k+1), and N-p sides with a length of 1/3^k.          | */
+/* | If JONES equals 1, the level is the power of 1/3 of the length of the | */
+/* | longest side hyperrectangle.                                          | */
+/* |                                                                       | */
+/* | On Return :                                                           | */
+/* |    the maximal length                                                 | */
+/* |                                                                       | */
+/* | pos     -- the position of the midpoint in the array length           | */
+/* | length  -- the array with the dimensions                              | */
+/* | maxfunc -- the leading dimension of length                            | */
+/* | n    -- the dimension of the problem                                  | */
+/* |                                                                       | */
+/* +-----------------------------------------------------------------------+ */
+integer direct_dirgetlevel_(integer *pos, integer *length, integer *maxfunc, integer 
+       *n, integer jones)
+{
+    /* System generated locals */
+    integer length_dim1, length_offset, ret_val, i__1;
+
+    /* Local variables */
+    integer i__, k, p, help;
+
+/* JG 09/15/00 Added variable JONES (see above) */
+    /* Parameter adjustments */
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+
+    /* Function Body */
+    if (jones == 0) {
+       help = length[*pos + length_dim1];
+       k = help;
+       p = 1;
+       i__1 = *n;
+       for (i__ = 2; i__ <= i__1; ++i__) {
+           if (length[*pos + i__ * length_dim1] < k) {
+               k = length[*pos + i__ * length_dim1];
+           }
+           if (length[*pos + i__ * length_dim1] == help) {
+               ++p;
+           }
+/* L100: */
+       }
+       if (k == help) {
+           ret_val = k * *n + *n - p;
+       } else {
+           ret_val = k * *n + p;
+       }
+    } else {
+       help = length[*pos + length_dim1];
+       i__1 = *n;
+       for (i__ = 2; i__ <= i__1; ++i__) {
+           if (length[*pos + i__ * length_dim1] < help) {
+               help = length[*pos + i__ * length_dim1];
+           }
+/* L10: */
+       }
+       ret_val = help;
+    }
+    return ret_val;
+} /* dirgetlevel_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | Program       : Direct.f (subfile DIRsubrout.f)                       | */
+/* | Last modified : 07-16-2001                                            | */
+/* | Written by    : Joerg Gablonsky                                       | */
+/* | Subroutines used by the algorithm DIRECT.                             | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRChoose                                               | */
+/* |    Decide, which is the next sampling point.                          | */
+/* |    Changed 09/25/00 JG                                                | */
+/* |         Added maxdiv to call and changed S to size maxdiv.            | */
+/* |    Changed 01/22/01 JG                                                | */
+/* |         Added Ifeasiblef to call to keep track if a feasible point has| */
+/* |         been found.                                                   | */
+/* |    Changed 07/16/01 JG                                                | */
+/* |         Changed if statement to prevent run-time errors.              |                                  
+                 | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirchoose_(integer *anchor, integer *s, integer *actdeep,
+        doublereal *f, doublereal *fmin, doublereal epsrel, doublereal epsabs, doublereal *thirds,
+        integer *maxpos, integer *length, integer *maxfunc, integer *maxdeep,
+        integer *maxdiv, integer *n, FILE *logfile,
+       integer *cheat, doublereal *kmax, integer *ifeasiblef, integer jones)
+{
+    /* System generated locals */
+    integer s_dim1, s_offset, length_dim1, length_offset, f_dim1, f_offset, 
+           i__1;
+
+    /* Local variables */
+    integer i__, j, k;
+    doublereal helplower;
+    integer i___, j___;
+    doublereal helpgreater;
+    integer novaluedeep;
+    doublereal help2;
+    integer novalue;
+
+    /* Parameter adjustments */
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    ++anchor;
+    s_dim1 = *maxdiv;
+    s_offset = 1 + s_dim1;
+    s -= s_offset;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+
+    /* Function Body */
+    helplower = 1e20;
+    helpgreater = 0.;
+    k = 1;
+    if (*ifeasiblef >= 1) {
+       i__1 = *actdeep;
+       for (j = 0; j <= i__1; ++j) {
+           if (anchor[j] > 0) {
+               s[k + s_dim1] = anchor[j];
+               s[k + (s_dim1 << 1)] = direct_dirgetlevel_(&s[k + s_dim1], &length[
+                       length_offset], maxfunc, n, jones);
+               goto L12;
+           }
+/* L1001: */
+       }
+L12:
+       ++k;
+       *maxpos = 1;
+       return;
+    } else {
+       i__1 = *actdeep;
+       for (j = 0; j <= i__1; ++j) {
+           if (anchor[j] > 0) {
+               s[k + s_dim1] = anchor[j];
+               s[k + (s_dim1 << 1)] = direct_dirgetlevel_(&s[k + s_dim1], &length[
+                       length_offset], maxfunc, n, jones);
+               ++k;
+           }
+/* L10: */
+       }
+    }
+    novalue = 0;
+    if (anchor[-1] > 0) {
+       novalue = anchor[-1];
+       novaluedeep = direct_dirgetlevel_(&novalue, &length[length_offset], maxfunc, 
+               n, jones);
+    }
+    *maxpos = k - 1;
+    i__1 = *maxdeep;
+    for (j = k - 1; j <= i__1; ++j) {
+       s[k + s_dim1] = 0;
+/* L11: */
+    }
+    for (j = *maxpos; j >= 1; --j) {
+       helplower = 1e20;
+       helpgreater = 0.;
+       j___ = s[j + s_dim1];
+       i__1 = j - 1;
+       for (i__ = 1; i__ <= i__1; ++i__) {
+           i___ = s[i__ + s_dim1];
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Changed IF statement into two to prevent run-time errors  | */
+/* |             which could occur if the compiler checks the second       | */
+/* |             expression in an .AND. statement although the first       | */
+/* |             statement is already not true.                            | */
+/* +-----------------------------------------------------------------------+ */
+           if (i___ > 0 && ! (i__ == j)) {
+               if (f[i___ + (f_dim1 << 1)] <= 1.) {
+                   help2 = thirds[s[i__ + (s_dim1 << 1)]] - thirds[s[j + (
+                           s_dim1 << 1)]];
+                   help2 = (f[i___ + f_dim1] - f[j___ + f_dim1]) / help2;
+                   if (help2 <= 0.) {
+                        if (logfile)
+                             fprintf(logfile, "thirds > 0, help2 <= 0\n");
+                       goto L60;
+                   }
+                   if (help2 < helplower) {
+                        if (logfile)
+                             fprintf(logfile, "helplower = %g\n", help2);
+                       helplower = help2;
+                   }
+               }
+           }
+/* L30: */
+       }
+       i__1 = *maxpos;
+       for (i__ = j + 1; i__ <= i__1; ++i__) {
+           i___ = s[i__ + s_dim1];
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Changed IF statement into two to prevent run-time errors  | */
+/* |             which could occur if the compiler checks the second       | */
+/* |             expression in an .AND. statement although the first       | */
+/* |             statement is already not true.                            | */
+/* +-----------------------------------------------------------------------+ */
+           if (i___ > 0 && ! (i__ == j)) {
+               if (f[i___ + (f_dim1 << 1)] <= 1.) {
+                   help2 = thirds[s[i__ + (s_dim1 << 1)]] - thirds[s[j + (
+                           s_dim1 << 1)]];
+                   help2 = (f[i___ + f_dim1] - f[j___ + f_dim1]) / help2;
+                   if (help2 <= 0.) {
+                       if (logfile)
+                            fprintf(logfile, "thirds < 0, help2 <= 0\n");
+                       goto L60;
+                   }
+                   if (help2 > helpgreater) {
+                       if (logfile)
+                              fprintf(logfile, "helpgreater = %g\n", help2);
+                       helpgreater = help2;
+                   }
+               }
+           }
+/* L31: */
+       }
+       if (helplower > 1e20 && helpgreater > 0.) {
+           helplower = helpgreater;
+           helpgreater += -1.;
+       }
+       if (helpgreater <= helplower) {
+           if (*cheat == 1 && helplower > *kmax) {
+               helplower = *kmax;
+           }
+           if (f[j___ + f_dim1] - helplower * thirds[s[j + (s_dim1 << 1)]] > 
+                    MIN(*fmin - epsrel * fabs(*fmin), 
+                        *fmin - epsabs)) {
+               if (logfile)
+                    fprintf(logfile, "> fmin - epslfminl\n");
+               goto L60;
+           }
+       } else {
+           if (logfile)
+                fprintf(logfile, "helpgreater > helplower: %g  %g  %g\n",
+                        helpgreater, helplower, helpgreater - helplower);
+           goto L60;
+       }
+       goto L40;
+L60:
+       s[j + s_dim1] = 0;
+L40:
+       ;
+    }
+    if (novalue > 0) {
+       ++(*maxpos);
+       s[*maxpos + s_dim1] = novalue;
+       s[*maxpos + (s_dim1 << 1)] = novaluedeep;
+    }
+} /* dirchoose_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRDoubleInsert                                         | */
+/* |      Routine to make sure that if there are several potential optimal | */
+/* |      hyperrectangles of the same level (i.e. hyperrectangles that have| */
+/* |      the same level and the same function value at the center), all of| */
+/* |      them are divided. This is the way as originally described in     | */
+/* |      Jones et.al.                                                     | */
+/* | JG 07/16/01 Added errorflag to calling sequence. We check if more     | */
+/* |             we reach the capacity of the array S. If this happens, we | */
+/* |             return to the main program with an error.                 | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirdoubleinsert_(integer *anchor, integer *s, integer *
+       maxpos, integer *point, doublereal *f, integer *maxdeep, integer *
+       maxfunc, integer *maxdiv, integer *ierror)
+{
+    /* System generated locals */
+    integer s_dim1, s_offset, f_dim1, f_offset, i__1;
+
+    /* Local variables */
+    integer i__, oldmaxpos, pos, help, iflag, actdeep;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Added flag to prevent run time-errors on some systems.    | */
+/* +-----------------------------------------------------------------------+ */
+    /* Parameter adjustments */
+    ++anchor;
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    --point;
+    s_dim1 = *maxdiv;
+    s_offset = 1 + s_dim1;
+    s -= s_offset;
+
+    /* Function Body */
+    oldmaxpos = *maxpos;
+    i__1 = oldmaxpos;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       if (s[i__ + s_dim1] > 0) {
+           actdeep = s[i__ + (s_dim1 << 1)];
+           help = anchor[actdeep];
+           pos = point[help];
+           iflag = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Added flag to prevent run time-errors on some systems. On | */
+/* |             some systems the second conditions in an AND statement is | */
+/* |             evaluated even if the first one is already not true.      | */
+/* +-----------------------------------------------------------------------+ */
+           while(pos > 0 && iflag == 0) {
+               if (f[pos + f_dim1] - f[help + f_dim1] <= 1e-13) {
+                   if (*maxpos < *maxdiv) {
+                       ++(*maxpos);
+                       s[*maxpos + s_dim1] = pos;
+                       s[*maxpos + (s_dim1 << 1)] = actdeep;
+                       pos = point[pos];
+                   } else {
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Maximum number of elements possible in S has been reached!| */
+/* +-----------------------------------------------------------------------+ */
+                       *ierror = -6;
+                       return;
+                   }
+               } else {
+                   iflag = 1;
+               }
+           }
+       }
+/* L10: */
+    }
+} /* dirdoubleinsert_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | INTEGER Function GetmaxDeep                                           | */
+/* | function to get the maximal length (1/length) of the n-dimensional    | */
+/* | rectangle with midpoint pos.                                          | */
+/* |                                                                       | */
+/* | On Return :                                                           | */
+/* |    the maximal length                                                 | */
+/* |                                                                       | */
+/* | pos     -- the position of the midpoint in the array length           | */
+/* | length  -- the array with the dimensions                              | */
+/* | maxfunc -- the leading dimension of length                            | */
+/* | n    -- the dimension of the problem                                  | */
+/* |                                                                       | */
+/* +-----------------------------------------------------------------------+ */
+integer direct_dirgetmaxdeep_(integer *pos, integer *length, integer *maxfunc, 
+       integer *n)
+{
+    /* System generated locals */
+    integer length_dim1, length_offset, i__1, i__2, i__3;
+
+    /* Local variables */
+    integer i__, help;
+
+    /* Parameter adjustments */
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+
+    /* Function Body */
+    help = length[*pos + length_dim1];
+    i__1 = *n;
+    for (i__ = 2; i__ <= i__1; ++i__) {
+/* Computing MIN */
+       i__2 = help, i__3 = length[*pos + i__ * length_dim1];
+       help = MIN(i__2,i__3);
+/* L10: */
+    }
+    return help;
+} /* dirgetmaxdeep_ */
+
+static integer isinbox_(doublereal *x, doublereal *a, doublereal *b, integer *n, 
+       integer *lmaxdim)
+{
+    /* System generated locals */
+    integer ret_val, i__1;
+
+    /* Local variables */
+    integer outofbox, i__;
+
+    /* Parameter adjustments */
+    --b;
+    --a;
+    --x;
+
+    /* Function Body */
+    outofbox = 1;
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       if (a[i__] > x[i__] || b[i__] < x[i__]) {
+           outofbox = 0;
+           goto L1010;
+       }
+/* L1000: */
+    }
+L1010:
+    ret_val = outofbox;
+    return ret_val;
+} /* isinbox_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG Added 09/25/00                                                     | */
+/* |                                                                       | */
+/* |                       SUBROUTINE DIRResortlist                        | */
+/* |                                                                       | */
+/* | Resort the list so that the infeasible point is in the list with the  | */
+/* | replaced value.                                                       | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirresortlist_(integer *replace, integer *anchor, 
+       doublereal *f, integer *point, integer *length, integer *n, integer *
+       maxfunc, integer *maxdim, integer *maxdeep, FILE *logfile,
+                                           integer jones)
+{
+    /* System generated locals */
+    integer f_dim1, f_offset, length_dim1, length_offset, i__1;
+
+    /* Local variables */
+    integer i__, l, pos;
+    integer start;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Get the length of the hyper rectangle with infeasible mid point and   | */
+/* | Index of the corresponding list.                                      | */
+/* +-----------------------------------------------------------------------+ */
+/* JG 09/25/00 Replaced with DIRgetlevel */
+/*      l = DIRgetmaxDeep(replace,length,maxfunc,n) */
+    /* Parameter adjustments */
+    --point;
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+    ++anchor;
+
+    /* Function Body */
+    l = direct_dirgetlevel_(replace, &length[length_offset], maxfunc, n, jones);
+    start = anchor[l];
+/* +-----------------------------------------------------------------------+ */
+/* | If the hyper rectangle with infeasibel midpoint is already the start  | */
+/* | of the list, give out message, nothing to do.                         | */
+/* +-----------------------------------------------------------------------+ */
+    if (*replace == start) {
+/*         write(logfile,*) 'No resorting of list necessarry, since new ', */
+/*     + 'point is already anchor of list .',l */
+    } else {
+/* +-----------------------------------------------------------------------+ */
+/* | Take the hyper rectangle with infeasible midpoint out of the list.    | */
+/* +-----------------------------------------------------------------------+ */
+       pos = start;
+       i__1 = *maxfunc;
+       for (i__ = 1; i__ <= i__1; ++i__) {
+           if (point[pos] == *replace) {
+               point[pos] = point[*replace];
+               goto L20;
+           } else {
+               pos = point[pos];
+           }
+           if (pos == 0) {
+               if (logfile)
+                    fprintf(logfile, "Error in DIRREsortlist: "
+                            "We went through the whole list\n"
+                            "and could not find the point to replace!!\n");
+               goto L20;
+           }
+/* L10: */
+       }
+/* +-----------------------------------------------------------------------+ */
+/* | If the anchor of the list has a higher value than the value of a      | */
+/* | nearby point, put the infeasible point at the beginning of the list.  | */
+/* +-----------------------------------------------------------------------+ */
+L20:
+       if (f[start + f_dim1] > f[*replace + f_dim1]) {
+           anchor[l] = *replace;
+           point[*replace] = start;
+/*            write(logfile,*) 'Point is replacing current anchor for ' */
+/*     +             , 'this list ',l,replace,start */
+       } else {
+/* +-----------------------------------------------------------------------+ */
+/* | Insert the point into the list according to its (replaced) function   | */
+/* | value.                                                                | */
+/* +-----------------------------------------------------------------------+ */
+           pos = start;
+           i__1 = *maxfunc;
+           for (i__ = 1; i__ <= i__1; ++i__) {
+/* +-----------------------------------------------------------------------+ */
+/* | The point has to be added at the end of the list.                     | */
+/* +-----------------------------------------------------------------------+ */
+               if (point[pos] == 0) {
+                   point[*replace] = point[pos];
+                   point[pos] = *replace;
+/*                  write(logfile,*) 'Point is added at the end of the ' */
+/*     +             , 'list ',l, replace */
+                   goto L40;
+               } else {
+                   if (f[point[pos] + f_dim1] > f[*replace + f_dim1]) {
+                       point[*replace] = point[pos];
+                       point[pos] = *replace;
+/*                     write(logfile,*) 'There are points with a higher ' */
+/*     +               ,'f-value in the list ',l,replace, pos */
+                       goto L40;
+                   }
+                   pos = point[pos];
+               }
+/* L30: */
+           }
+L40:
+           pos = pos;
+       }
+    }
+} /* dirresortlist_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG Added 09/25/00                                                     | */
+/* |                       SUBROUTINE DIRreplaceInf                        | */
+/* |                                                                       | */
+/* | Find out if there are infeasible points which are near feasible ones. | */
+/* | If this is the case, replace the function value at the center of the  | */
+/* | hyper rectangle by the lowest function value of a nearby function.    | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirreplaceinf_(integer *free, integer *freeold, 
+       doublereal *f, doublereal *c__, doublereal *thirds, integer *length, 
+       integer *anchor, integer *point, doublereal *c1, doublereal *c2, 
+       integer *maxfunc, integer *maxdeep, integer *maxdim, integer *n, 
+       FILE *logfile, doublereal *fmax, integer jones)
+{
+    /* System generated locals */
+    integer f_dim1, f_offset, c_dim1, c_offset, length_dim1, length_offset, 
+           i__1, i__2, i__3;
+    doublereal d__1, d__2;
+
+    /* Local variables */
+    doublereal a[32], b[32];
+    integer i__, j, k, l;
+    doublereal x[32], sidelength;
+    integer help;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
+/* +-----------------------------------------------------------------------+ */
+    /* Parameter adjustments */
+    --point;
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    ++anchor;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+    c_dim1 = *maxfunc;
+    c_offset = 1 + c_dim1;
+    c__ -= c_offset;
+    --c2;
+    --c1;
+
+    /* Function Body */
+    i__1 = *free - 1;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       if (f[i__ + (f_dim1 << 1)] > 0.) {
+/* +-----------------------------------------------------------------------+ */
+/* | Get the maximum side length of the hyper rectangle and then set the   | */
+/* | new side length to this lengths times the growth factor.              | */
+/* +-----------------------------------------------------------------------+ */
+           help = direct_dirgetmaxdeep_(&i__, &length[length_offset], maxfunc, n);
+           sidelength = thirds[help] * 2.;
+/* +-----------------------------------------------------------------------+ */
+/* | Set the Center and the upper and lower bounds of the rectangles.      | */
+/* +-----------------------------------------------------------------------+ */
+           i__2 = *n;
+           for (j = 1; j <= i__2; ++j) {
+               sidelength = thirds[length[i__ + j * length_dim1]];
+               a[j - 1] = c__[i__ + j * c_dim1] - sidelength;
+               b[j - 1] = c__[i__ + j * c_dim1] + sidelength;
+/* L20: */
+           }
+/* +-----------------------------------------------------------------------+ */
+/* | The function value is reset to 'Inf', since it may have been changed  | */
+/* | in an earlier iteration and now the feasible point which was close    | */
+/* | is not close anymore (since the hyper rectangle surrounding the       | */
+/* | current point may have shrunk).                                       | */
+/* +-----------------------------------------------------------------------+ */
+           f[i__ + f_dim1] = 1e6f;
+           f[i__ + (f_dim1 << 1)] = 2.;
+/* +-----------------------------------------------------------------------+ */
+/* | Check if any feasible point is near this infeasible point.            | */
+/* +-----------------------------------------------------------------------+ */
+           i__2 = *free - 1;
+           for (k = 1; k <= i__2; ++k) {
+/* +-----------------------------------------------------------------------+ */
+/* | If the point k is feasible, check if it is near.                      | */
+/* +-----------------------------------------------------------------------+ */
+               if (f[k + (f_dim1 << 1)] == 0.) {
+/* +-----------------------------------------------------------------------+ */
+/* | Copy the coordinates of the point k into x.                           | */
+/* +-----------------------------------------------------------------------+ */
+                   i__3 = *n;
+                   for (l = 1; l <= i__3; ++l) {
+                       x[l - 1] = c__[k + l * c_dim1];
+/* L40: */
+                   }
+/* +-----------------------------------------------------------------------+ */
+/* | Check if the point k is near the infeasible point, if so, replace the | */
+/* | value at */
+/* +-----------------------------------------------------------------------+ */
+                   if (isinbox_(x, a, b, n, &c__32) == 1) {
+/* Computing MIN */
+                       d__1 = f[i__ + f_dim1], d__2 = f[k + f_dim1];
+                       f[i__ + f_dim1] = MIN(d__1,d__2);
+                       f[i__ + (f_dim1 << 1)] = 1.;
+                   }
+               }
+/* L30: */
+           }
+           if (f[i__ + (f_dim1 << 1)] == 1.) {
+               f[i__ + f_dim1] += (d__1 = f[i__ + f_dim1], fabs(d__1)) * 
+                       1e-6f;
+               i__2 = *n;
+               for (l = 1; l <= i__2; ++l) {
+                   x[l - 1] = c__[i__ + l * c_dim1] * c1[l] + c__[i__ + l * 
+                           c_dim1] * c2[l];
+/* L200: */
+               }
+               dirresortlist_(&i__, &anchor[-1], &f[f_offset], &point[1], &
+                       length[length_offset], n, maxfunc, maxdim, maxdeep, 
+                       logfile, jones);
+           } else {
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01                                                           | */
+/* | Replaced fixed value for infeasible points with maximum value found,  | */
+/* | increased by 1.                                                       | */
+/* +-----------------------------------------------------------------------+ */
+               if (! (*fmax == f[i__ + f_dim1])) {
+/* Computing MAX */
+                   d__1 = *fmax + 1., d__2 = f[i__ + f_dim1];
+                   f[i__ + f_dim1] = MAX(d__1,d__2);
+               }
+           }
+       }
+/* L10: */
+    }
+/* L1000: */
+} /* dirreplaceinf_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRInsert                                               | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirinsert_(integer *start, integer *ins, integer *point, 
+       doublereal *f, integer *maxfunc)
+{
+    /* System generated locals */
+    integer f_dim1, f_offset, i__1;
+
+    /* Local variables */
+    integer i__, help;
+
+/* JG 09/17/00 Rewrote this routine. */
+/*      DO 10,i = 1,maxfunc */
+/*        IF (f(ins,1) .LT. f(point(start),1)) THEN */
+/*          help = point(start) */
+/*          point(start) = ins */
+/*          point(ins) = help */
+/*          GOTO 20 */
+/*        END IF */
+/*        IF (point(start) .EQ. 0) THEN */
+/*           point(start) = ins */
+/*           point(ins) = 0 */
+/*           GOTO 20 */
+/*        END IF */
+/*        start = point(start) */
+/* 10    CONTINUE */
+/* 20    END */
+    /* Parameter adjustments */
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    --point;
+
+    /* Function Body */
+    i__1 = *maxfunc;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       if (point[*start] == 0) {
+           point[*start] = *ins;
+           point[*ins] = 0;
+           return;
+       } else if (f[*ins + f_dim1] < f[point[*start] + f_dim1]) {
+            help = point[*start];
+            point[*start] = *ins;
+            point[*ins] = help;
+            return;
+       }
+       *start = point[*start];
+/* L10: */
+    }
+} /* dirinsert_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRInsertList                                           | */
+/* |    Changed 02-24-2000                                                 | */
+/* |      Got rid of the distinction between feasible and infeasible points| */
+/* |      I could do this since infeasible points get set to a high        | */
+/* |      function value, which may be replaced by a function value of a   | */
+/* |      nearby function at the end of the main loop.                     | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinsertlist_(integer *new__, integer *anchor, integer *
+       point, doublereal *f, integer *maxi, integer *length, integer *
+       maxfunc, integer *maxdeep, integer *n, integer *samp,
+                                           integer jones)
+{
+    /* System generated locals */
+    integer length_dim1, length_offset, f_dim1, f_offset, i__1;
+
+    /* Local variables */
+    integer j;
+    integer pos;
+    integer pos1, pos2, deep;
+
+/* JG 09/24/00 Changed this to Getlevel */
+    /* Parameter adjustments */
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    --point;
+    ++anchor;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+
+    /* Function Body */
+    i__1 = *maxi;
+    for (j = 1; j <= i__1; ++j) {
+       pos1 = *new__;
+       pos2 = point[pos1];
+       *new__ = point[pos2];
+/* JG 09/24/00 Changed this to Getlevel */
+/*        deep = DIRGetMaxdeep(pos1,length,maxfunc,n) */
+       deep = direct_dirgetlevel_(&pos1, &length[length_offset], maxfunc, n, jones);
+       if (anchor[deep] == 0) {
+           if (f[pos2 + f_dim1] < f[pos1 + f_dim1]) {
+               anchor[deep] = pos2;
+               point[pos2] = pos1;
+               point[pos1] = 0;
+           } else {
+               anchor[deep] = pos1;
+               point[pos2] = 0;
+           }
+       } else {
+           pos = anchor[deep];
+           if (f[pos2 + f_dim1] < f[pos1 + f_dim1]) {
+               if (f[pos2 + f_dim1] < f[pos + f_dim1]) {
+                   anchor[deep] = pos2;
+/* JG 08/30/00 Fixed bug. Sorting was not correct when */
+/*      f(pos2,1) < f(pos1,1) < f(pos,1) */
+                   if (f[pos1 + f_dim1] < f[pos + f_dim1]) {
+                       point[pos2] = pos1;
+                       point[pos1] = pos;
+                   } else {
+                       point[pos2] = pos;
+                       dirinsert_(&pos, &pos1, &point[1], &f[f_offset], 
+                               maxfunc);
+                   }
+               } else {
+                   dirinsert_(&pos, &pos2, &point[1], &f[f_offset], maxfunc);
+                   dirinsert_(&pos, &pos1, &point[1], &f[f_offset], maxfunc);
+               }
+           } else {
+               if (f[pos1 + f_dim1] < f[pos + f_dim1]) {
+/* JG 08/30/00 Fixed bug. Sorting was not correct when */
+/*      f(pos1,1) < f(pos2,1) < f(pos,1) */
+                   anchor[deep] = pos1;
+                   if (f[pos + f_dim1] < f[pos2 + f_dim1]) {
+                       point[pos1] = pos;
+                       dirinsert_(&pos, &pos2, &point[1], &f[f_offset], 
+                               maxfunc);
+                   } else {
+                       point[pos1] = pos2;
+                       point[pos2] = pos;
+                   }
+               } else {
+                   dirinsert_(&pos, &pos1, &point[1], &f[f_offset], maxfunc);
+                   dirinsert_(&pos, &pos2, &point[1], &f[f_offset], maxfunc);
+               }
+           }
+       }
+/* L10: */
+    }
+/* JG 09/24/00 Changed this to Getlevel */
+/*      deep = DIRGetMaxdeep(samp,length,maxfunc,n) */
+    deep = direct_dirgetlevel_(samp, &length[length_offset], maxfunc, n, jones);
+    pos = anchor[deep];
+    if (f[*samp + f_dim1] < f[pos + f_dim1]) {
+       anchor[deep] = *samp;
+       point[*samp] = pos;
+    } else {
+       dirinsert_(&pos, samp, &point[1], &f[f_offset], maxfunc);
+    }
+} /* dirinsertlist_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRInsertList2  (Old way to do it.)                     | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirinsertlist_2__(integer *start, integer *j, integer *k,
+        integer *list2, doublereal *w, integer *maxi, integer *n)
+{
+    /* System generated locals */
+    integer list2_dim1, list2_offset, i__1;
+
+    /* Local variables */
+    integer i__, pos;
+
+    /* Parameter adjustments */
+    --w;
+    list2_dim1 = *n;
+    list2_offset = 1 + list2_dim1;
+    list2 -= list2_offset;
+
+    /* Function Body */
+    pos = *start;
+    if (*start == 0) {
+       list2[*j + list2_dim1] = 0;
+       *start = *j;
+       goto L50;
+    }
+    if (w[*start] > w[*j]) {
+       list2[*j + list2_dim1] = *start;
+       *start = *j;
+    } else {
+       i__1 = *maxi;
+       for (i__ = 1; i__ <= i__1; ++i__) {
+           if (list2[pos + list2_dim1] == 0) {
+               list2[*j + list2_dim1] = 0;
+               list2[pos + list2_dim1] = *j;
+               goto L50;
+           } else {
+               if (w[*j] < w[list2[pos + list2_dim1]]) {
+                   list2[*j + list2_dim1] = list2[pos + list2_dim1];
+                   list2[pos + list2_dim1] = *j;
+                   goto L50;
+               }
+           }
+           pos = list2[pos + list2_dim1];
+/* L10: */
+       }
+    }
+L50:
+    list2[*j + (list2_dim1 << 1)] = *k;
+} /* dirinsertlist_2__ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRSearchmin                                            | */
+/* |    Search for the minimum in the list.                                ! */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirsearchmin_(integer *start, integer *list2, integer *
+       pos, integer *k, integer *n)
+{
+    /* System generated locals */
+    integer list2_dim1, list2_offset;
+
+    /* Parameter adjustments */
+    list2_dim1 = *n;
+    list2_offset = 1 + list2_dim1;
+    list2 -= list2_offset;
+
+    /* Function Body */
+    *k = *start;
+    *pos = list2[*start + (list2_dim1 << 1)];
+    *start = list2[*start + list2_dim1];
+} /* dirsearchmin_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRSamplepoints                                         | */
+/* |    Subroutine to sample the new points.                               | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirsamplepoints_(doublereal *c__, integer *arrayi, 
+       doublereal *delta, integer *sample, integer *start, integer *length, 
+       FILE *logfile, doublereal *f, integer *free, 
+       integer *maxi, integer *point, doublereal *x, doublereal *l,
+        doublereal *fmin, integer *minpos, doublereal *u, integer *n, 
+       integer *maxfunc, integer *maxdeep, integer *oops)
+{
+    /* System generated locals */
+    integer length_dim1, length_offset, c_dim1, c_offset, f_dim1, f_offset, 
+           i__1, i__2;
+
+    /* Local variables */
+    integer j, k, pos;
+
+
+    /* Parameter adjustments */
+    --u;
+    --l;
+    --x;
+    --arrayi;
+    --point;
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+    c_dim1 = *maxfunc;
+    c_offset = 1 + c_dim1;
+    c__ -= c_offset;
+
+    /* Function Body */
+    *oops = 0;
+    pos = *free;
+    *start = *free;
+    i__1 = *maxi + *maxi;
+    for (k = 1; k <= i__1; ++k) {
+       i__2 = *n;
+       for (j = 1; j <= i__2; ++j) {
+           length[*free + j * length_dim1] = length[*sample + j * 
+                   length_dim1];
+           c__[*free + j * c_dim1] = c__[*sample + j * c_dim1];
+/* L20: */
+       }
+       pos = *free;
+       *free = point[*free];
+       if (*free == 0) {
+            if (logfile)
+                 fprintf(logfile, "Error, no more free positions! "
+                         "Increase maxfunc!\n");
+           *oops = 1;
+           return;
+       }
+/* L10: */
+    }
+    point[pos] = 0;
+    pos = *start;
+    i__1 = *maxi;
+    for (j = 1; j <= i__1; ++j) {
+       c__[pos + arrayi[j] * c_dim1] = c__[*sample + arrayi[j] * c_dim1] + *
+               delta;
+       pos = point[pos];
+       c__[pos + arrayi[j] * c_dim1] = c__[*sample + arrayi[j] * c_dim1] - *
+               delta;
+       pos = point[pos];
+/* L30: */
+    }
+    ASRT(pos <= 0);
+} /* dirsamplepoints_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRDivide                                               | */
+/* |    Subroutine to divide the hyper rectangles according to the rules.  | */
+/* |    Changed 02-24-2000                                                 | */
+/* |      Replaced if statement by min (line 367)                          | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirdivide_(integer *new__, integer *currentlength, 
+       integer *length, integer *point, integer *arrayi, integer *sample, 
+       integer *list2, doublereal *w, integer *maxi, doublereal *f, integer *
+       maxfunc, integer *maxdeep, integer *n)
+{
+    /* System generated locals */
+    integer length_dim1, length_offset, list2_dim1, list2_offset, f_dim1, 
+           f_offset, i__1, i__2;
+    doublereal d__1, d__2;
+
+    /* Local variables */
+    integer i__, j, k, pos, pos2;
+    integer start;
+
+
+    /* Parameter adjustments */
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    --point;
+    --w;
+    list2_dim1 = *n;
+    list2_offset = 1 + list2_dim1;
+    list2 -= list2_offset;
+    --arrayi;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+
+    /* Function Body */
+    start = 0;
+    pos = *new__;
+    i__1 = *maxi;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       j = arrayi[i__];
+       w[j] = f[pos + f_dim1];
+       k = pos;
+       pos = point[pos];
+/* Computing MIN */
+       d__1 = f[pos + f_dim1], d__2 = w[j];
+       w[j] = MIN(d__1,d__2);
+       pos = point[pos];
+       dirinsertlist_2__(&start, &j, &k, &list2[list2_offset], &w[1], maxi, 
+               n);
+/* L10: */
+    }
+    ASRT(pos <= 0);
+    i__1 = *maxi;
+    for (j = 1; j <= i__1; ++j) {
+       dirsearchmin_(&start, &list2[list2_offset], &pos, &k, n);
+       pos2 = start;
+       length[*sample + k * length_dim1] = *currentlength + 1;
+       i__2 = *maxi - j + 1;
+       for (i__ = 1; i__ <= i__2; ++i__) {
+           length[pos + k * length_dim1] = *currentlength + 1;
+           pos = point[pos];
+           length[pos + k * length_dim1] = *currentlength + 1;
+/* JG 07/10/01 pos2 = 0 at the end of the 30-loop. Since we end */
+/*             the loop now, we do not need to reassign pos and pos2. */
+           if (pos2 > 0) {
+               pos = list2[pos2 + (list2_dim1 << 1)];
+               pos2 = list2[pos2 + list2_dim1];
+           }
+/* L30: */
+       }
+/* L20: */
+    }
+} /* dirdivide_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |                                                                       | */
+/* |                       SUBROUTINE DIRINFCN                             | */
+/* |                                                                       | */
+/* | Subroutine DIRinfcn unscales the variable x for use in the            | */
+/* | user-supplied function evaluation subroutine fcn. After fcn returns   | */
+/* | to DIRinfcn, DIRinfcn then rescales x for use by DIRECT.              | */
+/* |                                                                       | */
+/* | On entry                                                              | */
+/* |                                                                       | */
+/* |        fcn -- The argument containing the name of the user-supplied   | */
+/* |               subroutine that returns values for the function to be   | */
+/* |               minimized.                                              | */
+/* |                                                                       | */
+/* |          x -- A double-precision vector of length n. The point at     | */
+/* |               which the derivative is to be evaluated.                | */
+/* |                                                                       | */
+/* |        xs1 -- A double-precision vector of length n. Used for         | */
+/* |               scaling and unscaling the vector x by DIRinfcn.         | */
+/* |                                                                       | */
+/* |        xs2 -- A double-precision vector of length n. Used for         | */
+/* |               scaling and unscaling the vector x by DIRinfcn.         | */
+/* |                                                                       | */
+/* |          n -- An integer. The dimension of the problem.               | */
+/* |       kret -- An Integer. If kret =  1, the point is infeasible,      | */
+/* |                              kret = -1, bad problem set up,           | */
+/* |                              kret =  0, feasible.                     | */
+/* |                                                                       | */
+/* | On return                                                             | */
+/* |                                                                       | */
+/* |          f -- A double-precision scalar.                              | */
+/* |                                                                       | */
+/* | Subroutines and Functions                                             | */
+/* |                                                                       | */
+/* | The subroutine whose name is passed through the argument fcn.         | */
+/* |                                                                       | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinfcn_(fp fcn, doublereal *x, doublereal *c1, 
+       doublereal *c2, integer *n, doublereal *f, integer *flag__, 
+                                      void *fcn_data)
+{
+    /* System generated locals */
+    integer i__1;
+
+    /* Local variables */
+    integer i__;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized.  | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Unscale the variable x.                                               | */
+/* +-----------------------------------------------------------------------+ */
+    /* Parameter adjustments */
+    --c2;
+    --c1;
+    --x;
+
+    /* Function Body */
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       x[i__] = (x[i__] + c2[i__]) * c1[i__];
+/* L20: */
+    }
+/* +-----------------------------------------------------------------------+ */
+/* | Call the function-evaluation subroutine fcn.                          | */
+/* +-----------------------------------------------------------------------+ */
+    *flag__ = 0;
+    *f = fcn(*n, &x[1], flag__, fcn_data);
+/* +-----------------------------------------------------------------------+ */
+/* | Rescale the variable x.                                               | */
+/* +-----------------------------------------------------------------------+ */
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       x[i__] = x[i__] / c1[i__] - c2[i__];
+/* L30: */
+    }
+} /* dirinfcn_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRGet_I                                                | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirget_i__(integer *length, integer *pos, integer *
+       arrayi, integer *maxi, integer *n, integer *maxfunc)
+{
+    /* System generated locals */
+    integer length_dim1, length_offset, i__1;
+
+    /* Local variables */
+    integer i__, j, help;
+
+    /* Parameter adjustments */
+    --arrayi;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+
+    /* Function Body */
+    j = 1;
+    help = length[*pos + length_dim1];
+    i__1 = *n;
+    for (i__ = 2; i__ <= i__1; ++i__) {
+       if (length[*pos + i__ * length_dim1] < help) {
+           help = length[*pos + i__ * length_dim1];
+       }
+/* L10: */
+    }
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       if (length[*pos + i__ * length_dim1] == help) {
+           arrayi[j] = i__;
+           ++j;
+       }
+/* L20: */
+    }
+    *maxi = j - 1;
+} /* dirget_i__ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRInit                                                 | */
+/* |    Initialise all needed variables and do the first run of the        | */
+/* |    algorithm.                                                         | */
+/* |    Changed 02/24/2000                                                 | */
+/* |       Changed fcn Double precision to fcn external!                   | */
+/* |    Changed 09/15/2000                                                 | */
+/* |       Added distinction between Jones way to characterize rectangles  | */
+/* |       and our way. Common variable JONES controls which way we use.   | */
+/* |          JONES = 0    Jones way (Distance from midpoint to corner)    | */
+/* |          JONES = 1    Our way (Length of longest side)                | */
+/* |    Changed 09/24/00                                                   | */
+/* |       Added array levels. Levels contain the values to characterize   | */
+/* |       the hyperrectangles.                                            | */
+/* |    Changed 01/22/01                                                   | */
+/* |       Added variable fmax to keep track of maximum value found.       | */
+/* |       Added variable Ifeasiblef to keep track if feasibel point has   | */
+/* |       been found.                                                     | */
+/* |    Changed 01/23/01                                                   | */
+/* |       Added variable Ierror to keep track of errors.                  | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinit_(doublereal *f, fp fcn, doublereal *c__, 
+       integer *length, integer *actdeep, integer *point, integer *anchor, 
+       integer *free, FILE *logfile, integer *arrayi, 
+       integer *maxi, integer *list2, doublereal *w, doublereal *x, 
+       doublereal *l, doublereal *u, doublereal *fmin, integer *minpos, 
+       doublereal *thirds, doublereal *levels, integer *maxfunc, integer *
+       maxdeep, integer *n, integer *maxor, doublereal *fmax, integer *
+       ifeasiblef, integer *iinfeasible, integer *ierror, void *fcndata,
+                                     integer jones)
+{
+    /* System generated locals */
+    integer f_dim1, f_offset, c_dim1, c_offset, length_dim1, length_offset, 
+           list2_dim1, list2_offset, i__1, i__2;
+
+    /* Local variables */
+    integer i__, j;
+    integer new__, help, oops;
+    doublereal help2, delta;
+    doublereal costmin;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable Ifeasiblef to keep track if feasibel point | */
+/* |             has been found.                                           | */
+/* | JG 01/23/01 Added variable Ierror to keep track of errors.            | */
+/* | JG 03/09/01 Added IInfeasible to keep track if an infeasible point has| */
+/* |             been found.                                               | */
+/* +-----------------------------------------------------------------------+ */
+/* JG 09/15/00 Added variable JONES (see above) */
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized.  | */
+/* +-----------------------------------------------------------------------+ */
+    /* Parameter adjustments */
+    --point;
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    ++anchor;
+    --u;
+    --l;
+    --x;
+    --w;
+    list2_dim1 = *maxor;
+    list2_offset = 1 + list2_dim1;
+    list2 -= list2_offset;
+    --arrayi;
+    length_dim1 = *maxfunc;
+    length_offset = 1 + length_dim1;
+    length -= length_offset;
+    c_dim1 = *maxfunc;
+    c_offset = 1 + c_dim1;
+    c__ -= c_offset;
+
+    /* Function Body */
+    *fmin = 1e20;
+    costmin = *fmin;
+/* JG 09/15/00 If Jones way of characterising rectangles is used, */
+/*             initialise thirds to reflect this. */
+    if (jones == 0) {
+       i__1 = *n - 1;
+       for (j = 0; j <= i__1; ++j) {
+           w[j + 1] = sqrt(*n - j + j / 9.) * .5;
+/* L5: */
+       }
+       help2 = 1.;
+       i__1 = *maxdeep / *n;
+       for (i__ = 1; i__ <= i__1; ++i__) {
+           i__2 = *n - 1;
+           for (j = 0; j <= i__2; ++j) {
+               levels[(i__ - 1) * *n + j] = w[j + 1] / help2;
+/* L8: */
+           }
+           help2 *= 3.;
+/* L10: */
+       }
+    } else {
+/* JG 09/15/00 Initialiase levels to contain 1/j */
+       help2 = 3.;
+       i__1 = *maxdeep;
+       for (i__ = 1; i__ <= i__1; ++i__) {
+           levels[i__] = 1. / help2;
+           help2 *= 3.;
+/* L11: */
+       }
+       levels[0] = 1.;
+    }
+    help2 = 3.;
+    i__1 = *maxdeep;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       thirds[i__] = 1. / help2;
+       help2 *= 3.;
+/* L21: */
+    }
+    thirds[0] = 1.;
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       c__[i__ * c_dim1 + 1] = .5;
+       x[i__] = .5;
+       length[i__ * length_dim1 + 1] = 0;
+/* L20: */
+    }
+    direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &f[f_dim1 + 1], &help, fcndata);
+    f[(f_dim1 << 1) + 1] = (doublereal) help;
+    *iinfeasible = help;
+    *fmax = f[f_dim1 + 1];
+/* 09/25/00 Added this */
+/*      if (f(1,1) .ge. 1.E+6) then */
+    if (f[(f_dim1 << 1) + 1] > 0.) {
+       f[f_dim1 + 1] = 1e6;
+       *fmax = f[f_dim1 + 1];
+       *ifeasiblef = 1;
+    } else {
+       *ifeasiblef = 0;
+    }
+/* JG 09/25/00 Remove IF */
+    *fmin = f[f_dim1 + 1];
+    costmin = f[f_dim1 + 1];
+    *minpos = 1;
+    *actdeep = 2;
+    point[1] = 0;
+    *free = 2;
+    delta = thirds[1];
+    direct_dirget_i__(&length[length_offset], &c__1, &arrayi[1], maxi, n, maxfunc);
+    new__ = *free;
+    direct_dirsamplepoints_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &
+           length[length_offset], logfile, &f[f_offset], free, maxi, &
+           point[1], &x[1], &l[1], fmin, minpos, &u[1], n, 
+           maxfunc, maxdeep, &oops);
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/23/01 Added error checking.                                     | */
+/* +-----------------------------------------------------------------------+ */
+    if (oops > 0) {
+       *ierror = -4;
+       return;
+    }
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found.  | */
+/* |             Added variable to keep track if feasible point was found. | */
+/* +-----------------------------------------------------------------------+ */
+    direct_dirsamplef_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &length[
+           length_offset], logfile, &f[f_offset], free, maxi, &point[
+           1], fcn, &x[1], &l[1], fmin, minpos, &u[1], n, maxfunc, 
+           maxdeep, &oops, fmax, ifeasiblef, iinfeasible, fcndata);
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/23/01 Added error checking.                                     | */
+/* +-----------------------------------------------------------------------+ */
+    if (oops > 0) {
+       *ierror = -5;
+       return;
+    }
+    direct_dirdivide_(&new__, &c__0, &length[length_offset], &point[1], &arrayi[1], &
+           c__1, &list2[list2_offset], &w[1], maxi, &f[f_offset], maxfunc, 
+           maxdeep, n);
+    direct_dirinsertlist_(&new__, &anchor[-1], &point[1], &f[f_offset], maxi, &
+           length[length_offset], maxfunc, maxdeep, n, &c__1, jones);
+} /* dirinit_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRInitList                                             | */
+/* |    Initialise the list.                                               | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinitlist_(integer *anchor, integer *free, integer *
+       point, doublereal *f, integer *maxfunc, integer *maxdeep)
+{
+    /* System generated locals */
+    integer f_dim1, f_offset, i__1;
+
+    /* Local variables */
+    integer i__;
+
+/*   f -- values of functions. */
+/*   anchor -- anchors of lists with deep i */
+/*   point -- lists */
+/*   free  -- first free position */
+    /* Parameter adjustments */
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    --point;
+    ++anchor;
+
+    /* Function Body */
+    i__1 = *maxdeep;
+    for (i__ = -1; i__ <= i__1; ++i__) {
+       anchor[i__] = 0;
+/* L10: */
+    }
+    i__1 = *maxfunc;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       f[i__ + f_dim1] = 0.;
+       f[i__ + (f_dim1 << 1)] = 0.;
+       point[i__] = i__ + 1;
+/*       point(i) = 0 */
+/* L20: */
+    }
+    point[*maxfunc] = 0;
+    *free = 1;
+} /* dirinitlist_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRSort3                                                | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirsort3_(integer *pos1, integer *pos2, integer *pos3, 
+       doublereal *f, integer *maxfunc)
+{
+    /* System generated locals */
+    integer f_dim1, f_offset;
+
+    /* Local variables */
+    integer help;
+
+    /* Parameter adjustments */
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+
+    /* Function Body */
+    if (f[*pos1 + f_dim1] < f[*pos2 + f_dim1]) {
+       if (f[*pos1 + f_dim1] < f[*pos3 + f_dim1]) {
+           if (f[*pos3 + f_dim1] < f[*pos2 + f_dim1]) {
+               help = *pos2;
+               *pos2 = *pos3;
+               *pos3 = help;
+           }
+       } else {
+           help = *pos1;
+           *pos1 = *pos3;
+           *pos3 = *pos2;
+           *pos2 = help;
+       }
+    } else {
+       if (f[*pos2 + f_dim1] < f[*pos3 + f_dim1]) {
+           if (f[*pos3 + f_dim1] < f[*pos1 + f_dim1]) {
+               help = *pos1;
+               *pos1 = *pos2;
+               *pos2 = *pos3;
+               *pos3 = help;
+           } else {
+               help = *pos1;
+               *pos1 = *pos2;
+               *pos2 = help;
+           }
+       } else {
+           help = *pos1;
+           *pos1 = *pos3;
+           *pos3 = help;
+       }
+    }
+} /* dirsort3_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |    SUBROUTINE DIRInsert3                                              | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinsert3_(integer *pos1, integer *pos2, integer *pos3, 
+       integer *deep, integer *anchor, integer *point, integer *free, 
+       doublereal *f, doublereal *fmin, integer *minpos, integer *maxfunc, 
+       integer *maxdeep)
+{
+    /* System generated locals */
+    integer f_dim1, f_offset;
+
+    /* Local variables */
+    integer pos;
+
+    /* Parameter adjustments */
+    f_dim1 = *maxfunc;
+    f_offset = 1 + f_dim1;
+    f -= f_offset;
+    --point;
+    ++anchor;
+
+    /* Function Body */
+    dirsort3_(pos1, pos2, pos3, &f[f_offset], maxfunc);
+    if (anchor[*deep] == 0) {
+       anchor[*deep] = *pos1;
+       point[*pos1] = *pos2;
+       point[*pos2] = *pos3;
+       point[*pos3] = 0;
+    } else {
+       pos = anchor[*deep];
+       if (f[*pos1 + f_dim1] < f[pos + f_dim1]) {
+           anchor[*deep] = *pos1;
+           point[*pos1] = pos;
+       } else {
+           dirinsert_(&pos, pos1, &point[1], &f[f_offset], maxfunc);
+       }
+       dirinsert_(&pos, pos2, &point[1], &f[f_offset], maxfunc);
+       dirinsert_(&pos, pos3, &point[1], &f[f_offset], maxfunc);
+    }
+    if (f[*pos1 + f_dim1] < *fmin && f[*pos1 + (f_dim1 << 1)] == 0.) {
+       *fmin = f[*pos1 + f_dim1];
+       *minpos = *pos1;
+    }
+} /* dirinsert3_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* |                                                                       | */
+/* |                       SUBROUTINE DIRPREPRC                            | */
+/* |                                                                       | */
+/* | Subroutine DIRpreprc uses an afine mapping to map the hyper-box given | */
+/* | by the constraints on the variable x onto the n-dimensional unit cube.| */
+/* | This mapping is done using the following equation:                    | */
+/* |                                                                       | */
+/* |               x(i)=x(i)/(u(i)-l(i))-l(i)/(u(i)-l(i)).                 | */
+/* |                                                                       | */
+/* | DIRpreprc checks if the bounds l and u are well-defined. That is, if  | */
+/* |                                                                       | */
+/* |               l(i) < u(i) forevery i.                                 | */
+/* |                                                                       | */
+/* | On entry                                                              | */
+/* |                                                                       | */
+/* |          u -- A double-precision vector of length n. The vector       | */
+/* |               containing the upper bounds for the n independent       | */
+/* |               variables.                                              | */
+/* |                                                                       | */
+/* |          l -- A double-precision vector of length n. The vector       | */
+/* |               containing the lower bounds for the n independent       | */
+/* |               variables.                                              | */
+/* |                                                                       | */
+/* |          n -- An integer. The dimension of the problem.               | */
+/* |                                                                       | */
+/* | On return                                                             | */
+/* |                                                                       | */
+/* |        xs1 -- A double-precision vector of length n, used for scaling | */
+/* |               and unscaling the vector x.                             | */
+/* |                                                                       | */
+/* |        xs2 -- A double-precision vector of length n, used for scaling | */
+/* |               and unscaling the vector x.                             | */
+/* |                                                                       | */
+/* |                                                                       | */
+/* |       oops -- An integer. If an upper bound is less than a lower      | */
+/* |               bound or if the initial point is not in the             | */
+/* |               hyper-box oops is set to 1 and iffco terminates.        | */
+/* |                                                                       | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirpreprc_(doublereal *u, doublereal *l, integer *n, 
+       doublereal *xs1, doublereal *xs2, integer *oops)
+{
+    /* System generated locals */
+    integer i__1;
+
+    /* Local variables */
+    integer i__;
+    doublereal help;
+
+    /* Parameter adjustments */
+    --xs2;
+    --xs1;
+    --l;
+    --u;
+
+    /* Function Body */
+    *oops = 0;
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+/* +-----------------------------------------------------------------------+ */
+/* | Check if the hyper-box is well-defined.                               | */
+/* +-----------------------------------------------------------------------+ */
+       if (u[i__] <= l[i__]) {
+           *oops = 1;
+           return;
+       }
+/* L20: */
+    }
+/* +-----------------------------------------------------------------------+ */
+/* | Scale the initial iterate so that it is in the unit cube.             | */
+/* +-----------------------------------------------------------------------+ */
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       help = u[i__] - l[i__];
+       xs2[i__] = l[i__] / help;
+       xs1[i__] = help;
+/* L50: */
+    }
+} /* dirpreprc_ */
+
+/* Subroutine */ void direct_dirheader_(FILE *logfile, integer *version, 
+       doublereal *x, integer *n, doublereal *eps, integer *maxf, integer *
+       maxt, doublereal *l, doublereal *u, integer *algmethod, integer *
+       maxfunc, integer *maxdeep, doublereal *fglobal, doublereal *fglper, 
+       integer *ierror, doublereal *epsfix, integer *iepschange, doublereal *
+       volper, doublereal *sigmaper)
+{
+    /* System generated locals */
+    integer i__1;
+
+    /* Local variables */
+    integer imainver, i__, numerrors, isubsubver, ihelp, isubver;
+
+
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized.  | */
+/* +-----------------------------------------------------------------------+ */
+    /* Parameter adjustments */
+    --u;
+    --l;
+    --x;
+
+    /* Function Body */
+    if (logfile)
+        fprintf(logfile, "------------------- Log file ------------------\n");
+
+    numerrors = 0;
+    *ierror = 0;
+    imainver = *version / 100;
+    ihelp = *version - imainver * 100;
+    isubver = ihelp / 10;
+    ihelp -= isubver * 10;
+    isubsubver = ihelp;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/13/01 Added check for epsilon. If epsilon is smaller 0, we use  | */
+/* |             the update formula from Jones. We then set the flag       | */
+/* |             iepschange to 1, and store the absolute value of eps in   | */
+/* |             epsfix. epsilon is then changed after each iteration.     | */
+/* +-----------------------------------------------------------------------+ */
+    if (*eps < 0.) {
+       *iepschange = 1;
+       *epsfix = -(*eps);
+       *eps = -(*eps);
+    } else {
+       *iepschange = 0;
+       *epsfix = 1e100;
+    }
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Removed printout of contents in cdata(1).                 | */
+/* +-----------------------------------------------------------------------+ */
+/*      write(logfile,*) cdata(1) */
+
+    if (logfile) {
+        fprintf(logfile, "DIRECT Version %d.%d.%d\n"
+                " Problem dimension n: %d\n"
+                " Eps value: %e\n"
+                " Maximum number of f-evaluations (maxf): %d\n"
+                " Maximum number of iterations (MaxT): %d\n"
+                " Value of f_global: %e\n"
+                " Global percentage wanted: %e\n"
+                " Volume percentage wanted: %e\n"
+                " Measure percentage wanted: %e\n",
+                imainver, isubver, isubsubver, *n, *eps, *maxf, *maxt,
+                *fglobal, *fglper, *volper, *sigmaper);
+        fprintf(logfile, *iepschange == 1 
+                ? "Epsilon is changed using the Jones formula.\n" 
+                : "Epsilon is constant.\n");
+        fprintf(logfile, *algmethod == 0
+                ? "Jones original DIRECT algorithm is used.\n"
+                : "Our modification of the DIRECT algorithm is used.\n");
+    }
+
+    i__1 = *n;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+       if (u[i__] <= l[i__]) {
+           *ierror = -1;
+           if (logfile)
+                fprintf(logfile, "WARNING: bounds on variable x%d: "
+                        "%g <= xi <= %g\n", i__, l[i__], u[i__]);
+           ++numerrors;
+       } else {
+           if (logfile)
+                fprintf(logfile, "Bounds on variable x%d: "
+                        "%g <= xi <= %g\n", i__, l[i__], u[i__]);
+       }
+/* L1010: */
+    }
+/* +-----------------------------------------------------------------------+ */
+/* | If there are to many function evaluations or to many iteration, note  | */
+/* | this and set the error flag accordingly. Note: If more than one error | */
+/* | occurred, we give out an extra message.                               | */
+/* +-----------------------------------------------------------------------+ */
+    if (*maxf + 20 > *maxfunc) {
+       if (logfile)
+            fprintf(logfile,
+"WARNING: The maximum number of function evaluations (%d) is higher than\n"
+"         the constant maxfunc (%d).  Increase maxfunc in subroutine DIRECT\n"
+"         or decrease the maximum number of function evaluations.\n",
+                    *maxf, *maxfunc);
+       ++numerrors;
+       *ierror = -2;
+    }
+    if (*ierror < 0) {
+       if (logfile) fprintf(logfile, "----------------------------------\n");
+       if (numerrors == 1) {
+            if (logfile) 
+                 fprintf(logfile, "WARNING: One error in the input!\n");
+       } else {
+            if (logfile) 
+                 fprintf(logfile, "WARNING: %d errors in the input!\n",
+                         numerrors);
+       }
+    }
+    if (logfile) fprintf(logfile, "----------------------------------\n");
+    if (*ierror >= 0) {
+        if (logfile)
+             fprintf(logfile, "Iteration # of f-eval. fmin\n");
+    }
+/* L10005: */
+} /* dirheader_ */
+
+/* Subroutine */ void direct_dirsummary_(FILE *logfile, doublereal *x, doublereal *
+       l, doublereal *u, integer *n, doublereal *fmin, doublereal *fglobal, 
+       integer *numfunc, integer *ierror)
+{
+    /* Local variables */
+    integer i__;
+
+    /* Parameter adjustments */
+    --u;
+    --l;
+    --x;
+
+    /* Function Body */
+    if (logfile) {
+        fprintf(logfile, "-----------------------Summary------------------\n"
+                "Final function value: %g\n"
+                "Number of function evaluations: %d\n", *fmin, *numfunc);
+        if (*fglobal > -1e99)
+             fprintf(logfile, "Final function value is within %g%% of global optimum\n", 100*(*fmin - *fglobal) / MAX(1.0, fabs(*fglobal)));
+        fprintf(logfile, "Index, final solution, x(i)-l(i), u(i)-x(i)\n");
+        for (i__ = 1; i__ <= *n; ++i__)
+             fprintf(logfile, "%d, %g, %g, %g\n", i__, x[i__], 
+                     x[i__]-l[i__], u[i__] - x[i__]);
+        fprintf(logfile, "-----------------------------------------------\n");
+             
+    }
+} /* dirsummary_ */
+
+/* Subroutine */ void direct_dirmaxf_to_high1__(integer *maxf, integer *maxfunc, 
+                                       FILE *logfile)
+{
+     if (logfile) {
+         fprintf(logfile, 
+"The maximum number of function evaluations (%d)\n"
+"is higher than the constant maxfunc (%d).  Increase maxfunc\n"
+"in the SUBROUTINE DIRECT or decrease the maximum number\n"
+                 "of function evaluations.\n", *maxf, *maxfunc);
+     }
+} /* dirmaxf_to_high1__ */
+
+/* Subroutine */ void direct_dirmaxt_to_high1__(integer *maxt, integer *maxdeep, 
+                                       FILE *logfile)
+{
+     if (logfile) {
+         fprintf(logfile,
+"The maximum number of iterations (%d) is higher than the constant\n"
+"maxdeep (%d).  Increase maxdeep or decrease the number of iterations.\n",
+                 *maxt, *maxdeep);
+     }
+} /* dirmaxt_to_high1__ */
+
diff --git a/direct/direct-internal.h b/direct/direct-internal.h
new file mode 100644 (file)
index 0000000..84fd63e
--- /dev/null
@@ -0,0 +1,130 @@
+#ifndef DIRECT_INTERNAL_H
+#define DIRECT_INTERNAL_H
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "direct.h"
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif /* __cplusplus */
+
+typedef int integer;
+typedef double doublereal;
+typedef direct_objective_func fp;
+
+#define ASRT(c) if (!(c)) { fprintf(stderr, "DIRECT assertion failure at " __FILE__ ":%d -- " #c "\n", __LINE__); exit(EXIT_FAILURE); }
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+/* DIRect.c */
+extern void direct_dirsamplef_(
+     doublereal *c__, integer *arrayi, doublereal
+     *delta, integer *sample, integer *new__, integer *length,
+     FILE *logfile, doublereal *f, integer *free, integer *maxi,
+     integer *point, fp fcn, doublereal *x, doublereal *l, doublereal *
+     fmin, integer *minpos, doublereal *u, integer *n, integer *maxfunc,
+     integer *maxdeep, integer *oops, doublereal *fmax, integer *
+     ifeasiblef, integer *iinfesiblef, void *fcndata);
+
+/* DIRsubrout.c */
+
+extern void direct_dirheader_(
+     FILE *logfile, integer *version,
+     doublereal *x, integer *n, doublereal *eps, integer *maxf, integer *
+     maxt, doublereal *l, doublereal *u, integer *algmethod, integer *
+     maxfunc, integer *maxdeep, doublereal *fglobal, doublereal *fglper,
+     integer *ierror, doublereal *epsfix, integer *iepschange, doublereal *
+     volper, doublereal *sigmaper);
+extern void direct_dirinit_(
+     doublereal *f, fp fcn, doublereal *c__,
+     integer *length, integer *actdeep, integer *point, integer *anchor,
+     integer *free, FILE *logfile, integer *arrayi,
+     integer *maxi, integer *list2, doublereal *w, doublereal *x,
+     doublereal *l, doublereal *u, doublereal *fmin, integer *minpos,
+     doublereal *thirds, doublereal *levels, integer *maxfunc, integer *
+     maxdeep, integer *n, integer *maxor, doublereal *fmax, integer *
+     ifeasiblef, integer *iinfeasible, integer *ierror, void *fcndata,
+     integer jones);
+extern void direct_dirinitlist_(
+     integer *anchor, integer *free, integer *
+     point, doublereal *f, integer *maxfunc, integer *maxdeep);
+extern void direct_dirpreprc_(doublereal *u, doublereal *l, integer *n, 
+                             doublereal *xs1, doublereal *xs2, integer *oops);
+extern void direct_dirchoose_(
+     integer *anchor, integer *s, integer *actdeep,
+     doublereal *f, doublereal *fmin, doublereal epsrel, doublereal epsabs, doublereal *thirds,
+     integer *maxpos, integer *length, integer *maxfunc, integer *maxdeep,
+     integer *maxdiv, integer *n, FILE *logfile,
+     integer *cheat, doublereal *kmax, integer *ifeasiblef, integer jones);
+extern void direct_dirdoubleinsert_(
+     integer *anchor, integer *s, integer *maxpos, integer *point, 
+     doublereal *f, integer *maxdeep, integer *maxfunc, 
+     integer *maxdiv, integer *ierror);
+extern integer direct_dirgetmaxdeep_(integer *pos, integer *length, integer *maxfunc,
+                             integer *n);
+extern void direct_dirget_i__(
+     integer *length, integer *pos, integer *arrayi, integer *maxi, 
+     integer *n, integer *maxfunc);
+extern void direct_dirsamplepoints_(
+     doublereal *c__, integer *arrayi, 
+     doublereal *delta, integer *sample, integer *start, integer *length, 
+     FILE *logfile, doublereal *f, integer *free, 
+     integer *maxi, integer *point, doublereal *x, doublereal *l,
+     doublereal *fmin, integer *minpos, doublereal *u, integer *n, 
+     integer *maxfunc, integer *maxdeep, integer *oops);
+extern void direct_dirdivide_(
+     integer *new__, integer *currentlength, 
+     integer *length, integer *point, integer *arrayi, integer *sample, 
+     integer *list2, doublereal *w, integer *maxi, doublereal *f, 
+     integer *maxfunc, integer *maxdeep, integer *n);
+extern void direct_dirinsertlist_(
+     integer *new__, integer *anchor, integer *point, doublereal *f, 
+     integer *maxi, integer *length, integer *maxfunc, 
+     integer *maxdeep, integer *n, integer *samp, integer jones);
+extern void direct_dirreplaceinf_(
+     integer *free, integer *freeold, 
+     doublereal *f, doublereal *c__, doublereal *thirds, integer *length, 
+     integer *anchor, integer *point, doublereal *c1, doublereal *c2, 
+     integer *maxfunc, integer *maxdeep, integer *maxdim, integer *n, 
+     FILE *logfile, doublereal *fmax, integer jones);
+extern void direct_dirsummary_(
+     FILE *logfile, doublereal *x, doublereal *l, doublereal *u, 
+     integer *n, doublereal *fmin, doublereal *fglobal, 
+     integer *numfunc, integer *ierror);
+extern integer direct_dirgetlevel_(
+     integer *pos, integer *length, 
+     integer *maxfunc, integer *n, integer jones);
+extern void direct_dirinfcn_(
+     fp fcn, doublereal *x, doublereal *c1, 
+     doublereal *c2, integer *n, doublereal *f, integer *flag__, 
+     void *fcn_data);
+
+/* DIRserial.c / DIRparallel.c */
+extern void direct_dirsamplef_(
+     doublereal *c__, integer *arrayi, doublereal 
+     *delta, integer *sample, integer *new__, integer *length, 
+     FILE *logfile, doublereal *f, integer *free, integer *maxi, 
+     integer *point, fp fcn, doublereal *x, doublereal *l, doublereal *
+     fmin, integer *minpos, doublereal *u, integer *n, integer *maxfunc, 
+     integer *maxdeep, integer *oops, doublereal *fmax, integer *
+     ifeasiblef, integer *iinfesiblef, void *fcn_data);
+
+/* DIRect.c */
+
+extern void direct_direct_(
+     fp fcn, doublereal *x, integer *n, doublereal *eps, doublereal epsabs,
+     integer *maxf, integer *maxt, doublereal *fmin, doublereal *l, 
+     doublereal *u, integer *algmethod, integer *ierror, FILE *logfile, 
+     doublereal *fglobal, doublereal *fglper, doublereal *volper, 
+     doublereal *sigmaper, void *fcn_data);
+
+#ifdef __cplusplus
+}  /* extern "C" */
+#endif /* __cplusplus */
+
+#endif /* DIRECT_INTERNAL_H */
diff --git a/direct/direct.h b/direct/direct.h
new file mode 100644 (file)
index 0000000..41ed8c4
--- /dev/null
@@ -0,0 +1,59 @@
+#ifndef DIRECT_H
+#define DIRECT_H
+
+#include <math.h>
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif /* __cplusplus */
+
+typedef double (*direct_objective_func)(int n, const double *x,
+                                       int *undefined_flag, 
+                                       void *data);
+
+typedef enum {
+     DIRECT_ORIGINAL, DIRECT_GABLONSKY
+} direct_algorithm;
+
+typedef enum {
+     DIRECT_INVALID_BOUNDS = -1,
+     DIRECT_MAXFEVAL_TOOBIG = -2,
+     DIRECT_INIT_FAILED = -3,
+     DIRECT_SAMPLEPOINTS_FAILED = -4,
+     DIRECT_SAMPLE_FAILED = -5,
+     DIRECT_MAXFEVAL_EXCEEDED = 1,
+     DIRECT_MAXITER_EXCEEDED = 2,
+     DIRECT_GLOBAL_FOUND = 3,
+     DIRECT_VOLTOL = 4,
+     DIRECT_SIGMATOL = 5,
+
+     DIRECT_OUT_OF_MEMORY = -100,
+     DIRECT_INVALID_ARGS = -101
+} direct_return_code;
+
+#define DIRECT_UNKNOWN_FGLOBAL (-HUGE_VAL)
+#define DIRECT_UNKNOWN_FGLOBAL_RELTOL (0.0)
+
+extern direct_return_code direct_optimize(
+     direct_objective_func f, void *f_data,
+     int dimension,
+     const double *lower_bounds, const double *upper_bounds,
+
+     double *x, double *fmin, 
+
+     int max_feval, int max_iter,
+     double reltol, double abstol,
+     double volume_reltol, double sigma_reltol,
+
+     double fglobal,
+     double fglobal_reltol,
+
+     FILE *logfile,
+     direct_algorithm algorithm);
+
+#ifdef __cplusplus
+}  /* extern "C" */
+#endif /* __cplusplus */
+
+#endif /* DIRECT_H */
diff --git a/direct/direct_wrap.c b/direct/direct_wrap.c
new file mode 100644 (file)
index 0000000..406a44d
--- /dev/null
@@ -0,0 +1,97 @@
+/* C-style API for DIRECT functions.  SGJ (August 2007). */
+
+#include "direct-internal.h"
+
+/* Perform global minimization using (Gablonsky implementation of) DIRECT
+   algorithm.   Arguments:
+
+   f, f_data: the objective function and any user data
+       -- the objective function f(n, x, undefined_flag, data) takes 4 args:
+              int n: the dimension, same as dimension arg. to direct_optimize
+              const double *x: array x[n] of point to evaluate
+             int *undefined_flag: set to 1 on return if x violates constraints
+                                  or don't touch otherwise
+              void *data: same as f_data passed to direct_optimize
+          return value = value of f(x)
+
+   dimension: the number of minimization variable dimensions
+   lower_bounds, upper_bounds: arrays of length dimension of variable bounds
+
+   x: an array of length dimension, set to optimum variables upon return
+   fmin: on return, set to minimum f value
+
+   max_feval, max_iter: maximum number of function evaluations & DIRECT iters
+   reltol, abstol: relative and absolute tolerances (0 if none)
+   volume_reltol: relative tolerance on hypercube volume (0 if none)
+   sigma_reltol: relative tolerance on hypercube "measure" (??) (0 if none)
+
+   fglobal: the global minimum of f, if known ahead of time
+       -- this is mainly for benchmarking, in most cases it
+          is not known and you should pass DIRECT_UNKNOWN_FGLOBAL
+   fglobal_reltol: relative tolerance on how close we should find fglobal
+       -- ignored if fglobal is DIRECT_UNKNOWN_FGLOBAL
+
+   logfile: an output file to write diagnostic info to (NULL for no I/O)
+
+   algorithm: whether to use the original DIRECT algorithm (DIRECT_ORIGINAL)
+              or Gablonsky's "improved" version (DIRECT_GABLONSKY)
+*/
+direct_return_code direct_optimize(
+     direct_objective_func f, void *f_data,
+     int dimension,
+     const double *lower_bounds, const double *upper_bounds,
+
+     double *x, double *fmin, 
+
+     int max_feval, int max_iter,
+     double reltol, double abstol,
+     double volume_reltol, double sigma_reltol,
+
+     double fglobal,
+     double fglobal_reltol,
+
+     FILE *logfile,
+     direct_algorithm algorithm)
+{
+     integer algmethod = algorithm == DIRECT_GABLONSKY;
+     integer ierror;
+     doublereal *l, *u;
+     int i;
+
+     /* convert to percentages: */
+     volume_reltol *= 100;
+     sigma_reltol *= 100;
+     fglobal_reltol *= 100;
+
+     /* make sure these are ignored if <= 0 */
+     if (volume_reltol <= 0) volume_reltol = -1;
+     if (sigma_reltol <= 0) sigma_reltol = -1;
+
+     if (fglobal == DIRECT_UNKNOWN_FGLOBAL)
+         fglobal_reltol = DIRECT_UNKNOWN_FGLOBAL_RELTOL;
+
+     if (dimension < 1) return DIRECT_INVALID_ARGS;
+
+     l = (doublereal *) malloc(sizeof(doublereal) * dimension * 2);
+     if (!l) return DIRECT_OUT_OF_MEMORY;
+     u = l + dimension;
+     for (i = 0; i < dimension; ++i) {
+         l[i] = lower_bounds[i];
+         u[i] = upper_bounds[i];
+     }
+     
+     direct_direct_(f, x, &dimension, &reltol, abstol,
+                   &max_feval, &max_iter,
+                   fmin,
+                   l, u,
+                   &algmethod,
+                   &ierror,
+                   logfile,
+                   &fglobal, &fglobal_reltol,
+                   &volume_reltol, &sigma_reltol,
+                   f_data);
+
+     free(l);
+
+     return (direct_return_code) ierror;
+}
diff --git a/direct/tstc.c b/direct/tstc.c
new file mode 100644 (file)
index 0000000..23a161a
--- /dev/null
@@ -0,0 +1,41 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "direct.h"
+
+/* has two global minima at (0.09,-0.71) and (-0.09,0.71), plus
+   4 additional local minima */
+static int cnt=0;
+double tst_obj(int n, const double *xy, int *undefined_flag, void *unused)
+{
+  double x, y, f;
+  x = xy[0];
+  y = xy[1];
+  f = ((x*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3) + x*y + (y*y)*(-4+4*(y*y)));
+  printf("feval:, %d, %g, %g, %g\n", ++cnt, x,y, f);
+  return f;
+}
+
+int main(int argc, char **argv)
+{
+  int n = 2;
+  double x[2], l[2], u[2];
+  long int maxits = 0;
+  int info;
+  double fmin;
+
+  maxits = argc < 2 ? 100 : atoi(argv[1]);
+
+  l[0] = -3; l[1] = -3;
+  u[0] = 3; u[1] = 3;
+
+  info = direct_optimize(tst_obj, NULL, n, l, u, x, &fmin,
+                        maxits, 500,
+                        0, 0, 0, 0, DIRECT_UNKNOWN_FGLOBAL, 0,
+                        stdout, DIRECT_GABLONSKY);
+
+  printf("min f = %g at (%g,%g) after %d evals, return value %d\n",
+        fmin, x[0], x[1], cnt, info);
+
+  return EXIT_SUCCESS;
+}
diff --git a/direct/userguide.ps.gz b/direct/userguide.ps.gz
new file mode 100644 (file)
index 0000000..c3a1655
Binary files /dev/null and b/direct/userguide.ps.gz differ
diff --git a/stogo/README b/stogo/README
new file mode 100644 (file)
index 0000000..71c6c62
--- /dev/null
@@ -0,0 +1,21 @@
+This code is modified from the original StoGO Global Optimization library
+by Madsen et al., downloaded from:
+
+       http://www2.imm.dtu.dk/~km/GlobOpt/opt.html
+
+It was modified to allow C-callable wrappers and some other niceties by
+Steven G. Johnson (stevenj@alum.mit.edu) in 2007.
+
+StoGO uses a gradient-based direct-search branch-and-bound algorithm,
+described in:
+
+S. Gudmundsson, "Parallel Global Optimization," M.Sc. Thesis, IMM,
+       Technical University of Denmark, 1998.
+
+K. Madsen, S. Zertchaninov, and A. Zilinskas, "Global Optimization
+       using Branch-and-Bound," Submitted to the Journal of Global
+       Optimization, 1998.  [ does not seem to have been published? ]
+
+S. Zertchaninov and K. Madsen, "A C++ Programme for Global Optimization,"
+       IMM-REP-1998-04, Department of Mathematical Modelling,
+       Technical University of Denmark, DK-2800 Lyngby, Denmark, 1998.
diff --git a/stogo/global.cc b/stogo/global.cc
new file mode 100644 (file)
index 0000000..9496358
--- /dev/null
@@ -0,0 +1,315 @@
+/*
+   Multi Dimensional Global Search.
+
+   Author: Steinn Gudmundsson
+   Email: steinng@hotmail.com
+
+   This program is supplied without any warranty whatsoever.
+
+   NB The RNGs seed should be initialized using some timer
+*/
+
+#include <iostream>
+#include <time.h>
+
+#include <iterator>
+#include <algorithm>
+#include <stack>
+
+#include "stogo_config.h"
+#include "global.h"
+#include "local.h"
+
+// Timer stuff
+time_t   StartTime;
+
+double MacEpsilon ;
+int FC=0, GC=0 ;
+
+Global::Global(RTBox D, Pobj o, Pgrad g, GlobalParams P): Domain(D) {
+
+  dim=Domain.GetDim();
+  Objective=o;
+  Gradient=g;
+
+  // Initialize parameters
+  maxtime=P.maxtime;
+  maxeval = P.maxeval;
+  numeval = 0;
+  eps_cl=P.eps_cl; mu=P.mu; rshift=P.rshift;
+  det_pnts=P.det_pnts; rnd_pnts=P.rnd_pnts;
+  fbound=DBL_MAX;
+}
+
+#if 0 // not necessary; default copy is sufficient 
+Global& Global::operator=(const Global &G) {
+  // Copy the problem info and parameter settings
+  Domain=G.Domain; Objective=G.Objective;  Gradient=G.Gradient;
+  maxtime=G.maxtime;
+  maxeval = G.maxeval;
+  numeval = G.numeval;
+  eps_cl=G.eps_cl; mu=G.mu; rshift=G.rshift;
+  det_pnts=G.det_pnts; rnd_pnts=G.rnd_pnts;
+  return *this;
+}
+#endif
+
+void Global::FillRegular(RTBox SampleBox, RTBox box) {
+  // Generation of regular sampling points
+  double w;
+  int i, flag, dir;
+  Trial tmpTrial(dim);
+  RVector m(dim), x(dim);
+
+  if (det_pnts>0) {
+    // Add midpoint
+    box.Midpoint(m) ;
+    tmpTrial.xvals=m ; tmpTrial.objval=DBL_MAX ;
+    SampleBox.AddTrial(tmpTrial) ;
+    // Add the rest
+    i=1 ; flag=1 ; dir=0 ;
+    x=m ; w=box.Width(dir) ;
+    while (i<det_pnts) {
+      x(dir)=m(dir)+flag*rshift*w ;
+      tmpTrial.xvals=x ; 
+      SampleBox.AddTrial(tmpTrial) ;
+      flag=-flag;
+      if (flag==1 && dir<dim) {
+       x(dir)=m(dir) ;
+       dir++ ;
+       w=box.Width(dir) ;
+      }
+      i++ ;
+    }
+  }
+}
+
+void Global::FillRandom(RTBox SampleBox, RTBox box) {
+  // Generation of stochastic sampling points
+  Trial tmpTrial(dim);
+
+  tmpTrial.objval=DBL_MAX;
+  for (int i=1 ; i<=rnd_pnts ; i++) {
+    for (int dir=0 ; dir<dim ; dir++)
+      tmpTrial.xvals(dir) =
+       box.lb(dir)+(box.ub(dir)-box.lb(dir))*(double(rand())/RAND_MAX) ;
+    SampleBox.AddTrial(tmpTrial) ;
+  }
+}
+
+double Global::NewtonTest(RTBox box, int axis, RCRVector x_av, int *noutside) {
+  // Perform the Newton test
+
+  int info,nout=0;
+  Trial tmpTrial(dim);
+  TBox SampleBox(dim) ;
+  double maxgrad=0 ;
+
+  // Create sampling points
+  FillRegular(SampleBox, box);
+  FillRandom(SampleBox, box);
+
+  // Perform the actual sampling
+  while ( !SampleBox.EmptyBox() ) {
+    SampleBox.RemoveTrial(tmpTrial) ;
+    info = local(tmpTrial, box, Domain, eps_cl, &maxgrad, *this,
+                axis, x_av) ;
+    // What should we do when info=LS_Unstable?
+    if (info == LS_Out)
+      nout++;
+    if (info == LS_New ) {
+      box.AddTrial(tmpTrial) ;
+
+      if (tmpTrial.objval<=fbound+mu && tmpTrial.objval<=box.fmin+mu) {
+       cout << "Found a candidate, x=" << tmpTrial.xvals;
+       cout << " F=" <<tmpTrial.objval << " FC=" << FC << endl;
+       SolSet.push_back(tmpTrial);
+      }
+#ifdef GS_DEBUG
+      cout << "Found a stationary point, X= " << tmpTrial.xvals;
+      cout <<" objval=" << tmpTrial.objval << endl;
+#endif
+    }
+
+    if (!InTime())
+      break;
+  }
+  *noutside=nout;
+  return maxgrad;
+}
+
+void Global::ReduceOrSubdivide(RTBox box, int axis, RCRVector x_av) {
+  TBox B1(dim), B2(dim);
+  Trial tmpTrial(dim);
+  double maxgrad;
+  int ns,nout;
+
+  // Monotonicity test has not been implemented yet
+  maxgrad=NewtonTest(box, axis, x_av, &nout);
+  ns=box.NStationary() ;
+  if (ns==0) {
+    // All iterates outside
+    // NB result=Intersection(B,boundary(Domain))
+    Garbage.push(box) ;
+  }
+  else
+    if (ns==1 && nout==0) {
+      // All iterates converge to same point
+      Garbage.push(box) ;
+    }
+    else
+      if ( (ns>1) && (box.LowerBound(maxgrad)>fbound) ) {
+       // Several stationary points found and lower bound > fbound
+       Garbage.push(box) ;
+      }
+      else {
+       // Subdivision
+       B1.ClearBox() ; B2.ClearBox() ;
+       box.split(B1,B2) ;
+       CandSet.push(B1) ; CandSet.push(B2) ;
+      }
+
+  // Update fbound
+  if (box.fmin < fbound) {
+    fbound=box.fmin ;
+#ifdef GS_DEBUG
+    cout <<"*** Improving fbound, fbound=" << fbound << endl;
+#endif
+  }
+}
+
+void Global::Search(int axis, RCRVector x_av){
+  Trial tmpTrial(dim) ;
+  TBox box(dim), B1(dim), B2(dim);
+  RVector m(dim), x(dim);
+  int inner_iter, outer_iter;
+
+  MacEpsilon=eps(); // Get machine precision
+  if (det_pnts>2*dim+1) {
+    det_pnts=2*dim+1;
+    cout << "Warning: Reducing det_pnts to " << det_pnts << endl;
+  }
+
+  // Initialize timer
+  time(&StartTime);
+
+  // Clear priority_queues
+  while (!Garbage.empty())
+    Garbage.pop();
+  while (!CandSet.empty())
+    CandSet.pop();
+
+  box=Domain;
+  CandSet.push(box);
+  int done=0 ; outer_iter=0 ;
+
+  while (!done) {
+    outer_iter++ ;
+
+    // Inner loop
+    inner_iter=0 ;
+    while (!CandSet.empty()) {
+      inner_iter++ ;
+      // Get best box from Candidate set
+      box=CandSet.top() ; CandSet.pop() ;
+
+#ifdef GS_DEBUG
+      cout << "Iteration..." << inner_iter << " #CS=" << CandSet.size()+1 ;
+      cout << " Processing " << box.NStationary() << " trials in the box " <<box;
+#endif
+      ReduceOrSubdivide(box, axis, x_av);
+
+      if (!InTime()) {
+        done=TRUE;
+        cout << "The program has run out of time or function evaluations\n";
+        break;
+      }
+
+    } // inner while-loop
+    cout << endl << "*** Inner loop completed ***" << endl ;
+    
+    // Reduce SolSet if necessary
+    SolSet.erase(remove_if(SolSet.begin(), SolSet.end(),
+                          TrialGT(fbound+mu)),SolSet.end());
+    if (InTime()) {
+      cout << "Current set of minimizers (" << SolSet.size() << ")" << endl ;
+      DispMinimizers() ;
+
+      while (!Garbage.empty()) {
+        box=Garbage.top() ;
+        Garbage.pop() ;
+        // Split box
+        B1.ClearBox() ; B2.ClearBox() ;
+        box.split(B1,B2) ;
+        // Add boxes to Candidate set
+        CandSet.push(B1) ; CandSet.push(B2) ;
+      }
+    }
+  } // Outer while-loop
+
+  cout << "Number of outer iterations : " << outer_iter << endl;
+  cout << "Number of unexplored boxes : " << CandSet.size() << endl;
+  cout << "Number of boxes in garbage : " << Garbage.size() << endl;
+  cout << "Number of elements in SolSet : " << SolSet.size() << endl;
+  cout << "Number of function evaluations : " << FC << endl;
+  cout << "Number of gradient evaluations : " << GC << endl;
+
+  if (axis != -1) {
+    // Return minimizer when doing the AV method
+    tmpTrial=SolSet.back();
+    x_av(axis)=tmpTrial.xvals(0);
+  }
+}
+
+/************* Various utility functions ****************/
+long int Global::GetTime()
+{
+ time_t ctime; time(&ctime);
+ return (long int)difftime(ctime,StartTime);
+}
+
+bool Global::InTime()
+{
+ return (!maxtime || GetTime()<maxtime) && (!maxeval || numeval<maxeval);
+}
+
+double Global::GetMinValue() {
+  return fbound;
+}
+
+void Global::SetMinValue(double new_fb) {
+  fbound=new_fb;
+}
+
+void Global::SetDomain(RTBox box) {
+  Domain=box;
+}
+
+void Global::GetDomain(RTBox box) {
+  box=Domain;
+}
+
+void Global::DispMinimizers() {
+  copy(SolSet.begin(), SolSet.end(), ostream_iterator<Trial>(cout));
+}
+
+double Global::OneMinimizer(RCRVector x) {
+  if (NoMinimizers()) return 0.0;
+  for (int i=0;i<x.GetLength();i++) x(i) = SolSet.front().xvals(i);
+  return SolSet.front().objval;
+}
+
+bool Global::NoMinimizers() {
+  return SolSet.empty();
+}
+
+void Global::ClearSolSet() {
+  SolSet.erase(SolSet.begin(), SolSet.end()) ;
+}
+
+void Global::AddPoint(RCRVector x, double f) {
+  Trial T(dim);
+  T.xvals=x; T.objval=f;
+  Domain.AddTrial(T);
+  SolSet.push_back(T);
+}
diff --git a/stogo/global.h b/stogo/global.h
new file mode 100644 (file)
index 0000000..a6b568e
--- /dev/null
@@ -0,0 +1,83 @@
+#ifndef GLOBAL_H
+#define GLOBAL_H
+
+#include <queue>
+//#include "function.h"
+#include "tools.h"
+using namespace std;
+
+typedef void dom(RTBox) ;
+typedef dom* Pdom ;
+
+typedef double obj(RCRVector) ;
+typedef obj* Pobj ;
+
+typedef void grad(RCRVector,RVector&) ;
+typedef grad* Pgrad ;
+
+typedef enum { OBJECTIVE_ONLY, GRADIENT_ONLY, OBJECTIVE_AND_GRADIENT } whichO;
+
+typedef double objgrad(RCRVector,RCRVector,whichO) ;
+typedef objgrad* Pobjgrad ;
+
+class GlobalParams {
+public:
+  long int maxtime, maxeval;
+  double eps_cl, mu, rshift;
+  int det_pnts, rnd_pnts;
+};
+
+class Global: public GlobalParams {
+public:
+  // Problem specification
+  int dim ;
+  Pobj  Objective ;
+  Pgrad Gradient ;
+  long int numeval;
+
+  virtual double ObjectiveGradient(RCRVector xy, RVector&grad, whichO which){
+       ++numeval;
+       switch (which) {
+          case OBJECTIVE_AND_GRADIENT:
+               Gradient(xy, grad);
+          case OBJECTIVE_ONLY:
+               return Objective(xy);
+          case GRADIENT_ONLY:
+               Gradient(xy, grad);
+               return 0.0;
+       }
+  }
+                                  
+  Global(RTBox, Pobj, Pgrad, GlobalParams);
+//  Global& operator=(const Global &);
+
+  void Search(int, RCRVector);
+  void DispMinimizers();
+  double OneMinimizer(RCRVector);
+  bool NoMinimizers();
+  void SetDomain(RTBox);
+  void GetDomain(RTBox);
+  double GetMinValue();
+  void SetMinValue(double);
+  void ClearSolSet();
+  void AddPoint(RCRVector, double);
+
+  long int GetTime();
+  bool InTime();
+
+private:
+  list<Trial> SolSet;
+  list<Trial>::const_iterator titr;
+  priority_queue<TBox> CandSet;
+  priority_queue<TBox> Garbage;
+
+  double fbound;
+  TBox Domain;
+
+  void FillRegular(RTBox, RTBox);
+  void FillRandom(RTBox, RTBox);
+  double NewtonTest(RTBox, int, RCRVector, int*);
+  void ReduceOrSubdivide(RTBox, int, RCRVector);
+};
+#endif
+
diff --git a/stogo/linalg.cc b/stogo/linalg.cc
new file mode 100644 (file)
index 0000000..04916cf
--- /dev/null
@@ -0,0 +1,216 @@
+/*
+   Temporary implementation of vector and matrix classes
+   No attempt is made to check if the function arguments are valid
+*/
+
+#include <iostream>
+#include <math.h>         // for sqrt()
+
+#include "linalg.h"
+
+double eps() {
+  /* Returns the machine precision : (min { x >= 0 : 1 + x > 1 }) 
+     NB This routine should be replaced by LAPACK_LAMCH */
+  double Current, Last, OnePlusCurrent ;
+  Current = 1.0 ;
+  do
+    {
+      Last = Current ;
+      Current /= 2.0 ;
+      OnePlusCurrent = 1.0 + Current ;
+    } while (OnePlusCurrent > 1.0) ;
+  return Last ;
+}
+
+
+RVector::RVector() {
+ // Constructor
+ len=0;
+ elements=0; (*this)=0.;
+}
+
+RVector::RVector(int n) {
+ // Constructor
+ len=n;
+ elements=new double[len]; (*this)=0.;
+}
+
+RVector::RVector(RCRVector vect)
+{
+ // Constructor + Copy
+ len=vect.len;
+ elements=new double[len]; (*this)=vect;
+}
+
+RCRVector RVector::operator=(RCRVector vect)
+{
+ // Copy
+ for (int i=0;i<len;i++) elements[i]=vect.elements[i];
+ return *this;
+}
+
+RCRVector RVector::operator=(double num) {
+ // Assignment
+ for (int i=0;i<len;i++) elements[i]=num;
+ return *this;
+}
+
+double RVector::nrm2() {
+  double sum=0 ;
+  for (int i = 0; i < len; i++)
+    sum+=elements[i]*elements[i] ;
+  return sqrt(sum) ;
+}
+
+void scal(double alpha, RCRVector x) {
+  int n=x.len ;
+  for (int i = 0; i < n; i++)
+    x.elements[i] = alpha*x.elements[i] ;
+}
+
+double norm2(RCRVector x) {
+  // Euclidian norm
+  double sum=0 ;
+  double* pa=x.elements ;
+  int n=x.len ;
+  for (int i = 0; i < n; i++) {
+    sum+=(*pa)*(*pa) ;
+    pa++ ;
+  }
+  return sqrt(sum) ;
+}
+
+double normInf(RCRVector x) {
+  // Infinity norm
+  int n=x.len ;
+  double tmp=DBL_MIN ;
+  for (int i=0 ; i<n ; i++)
+    tmp=max(tmp,fabs(x.elements[i])) ;
+  return tmp ;
+}
+
+double dot(RCRVector x, RCRVector y) {
+  // dot <- x'y
+  int n=x.len ;
+  double sum=0 ;
+  for (int i=0 ; i<n ; i++)
+    sum+= x.elements[i]*y.elements[i] ;
+  return sum ;
+}
+
+void copy(RCRVector x, RCRVector y) {
+  // y <- x
+  double *px=x.elements ;
+  double *py=y.elements ;
+  int n=x.len ;
+  for (int i=0 ; i<n ; i++)
+    (*py++)=(*px++) ;
+}
+
+void axpy(double alpha, RCRVector x, RCRVector y) {
+  // y <- alpha*x + y
+  double *px=x.elements ;
+  double *py=y.elements ;
+  int n=x.len ;
+  for (int i=0 ; i<n ; i++) {
+    *py=alpha*(*px) + *py ;
+    px++ ; py++ ;
+  }
+}
+
+void gemv(char trans, double alpha, RCRMatrix A,RCRVector x,
+         double beta, RCRVector y) {
+  // Matrix-vector multiplication
+  int i, j, dim=A.Dim;
+  double sum ;
+
+  if (trans=='N') {
+    // y=alpha*A*x + beta*y
+    for (i=0;i<dim;i++) {
+      sum=0.;
+      for (j=0;j<dim;j++)
+        sum+=A.Vals[j+i*dim]*x.elements[j]*alpha;
+      y.elements[i]=y.elements[i]*beta + sum ;
+    }
+  }
+  else {
+    // y=alpha*transpose(A)*x +  beta*y
+    for (i=0;i<dim;i++) {
+      sum=0.;
+      for (j=0;j<dim;j++) {
+        sum+=A.Vals[i+j*dim]*x.elements[j]*alpha ;
+      }
+      y.elements[i]=y.elements[i]*beta + sum ;
+    }
+  }
+}
+
+ostream & operator << (ostream & os, const RVector & v) {
+  os << '[';
+  for (int i = 0; i < v.len; i++) {
+    if (i>0) os << "," ;
+    os << v.elements[i] ;
+  }
+  return os << ']';
+}
+
+/*************************** Matrix Class ***********************/
+
+RMatrix::RMatrix() {
+ // Constructor
+ Dim=0 ; Vals=0 ; (*this)=0 ;
+}
+
+RMatrix::RMatrix(int dim) {
+ // Constructor
+ Dim=dim;
+ Vals=new double[long(Dim)*long(Dim)]; (*this)=0.;
+}
+
+RMatrix::RMatrix(RCRMatrix matr) {
+ // Constructor + Copy 
+ Dim=matr.Dim;
+ Vals=new double[long(Dim)*long(Dim)]; (*this)=matr;
+}
+
+RCRMatrix RMatrix::operator=(RCRMatrix mat)
+{ // Assignment, A=B
+ long int Len=long(Dim)*long(Dim);
+ for (int i=0;i<Len;i++) Vals[i]=mat.Vals[i] ;
+ return *this;
+}
+
+RCRMatrix RMatrix::operator=(double num) {
+ long int Len=long(Dim)*long(Dim);
+ for (long int i=0;i<Len;i++) Vals[i]=num;
+ return *this;
+}
+
+double& RMatrix::operator()(int vidx,int hidx) {
+ return Vals[vidx*Dim+hidx];
+}
+
+void ger(double alpha, RCRVector x,RCRVector y, RCRMatrix A) {
+  // Rank one update : A=alpha*xy'+A
+  int dim=x.len;
+  double* pa=A.Vals ;
+
+  for (int i=0;i<dim;i++)
+    for (int j=0;j<dim;j++) {
+      *pa=alpha*x.elements[i]*y.elements[j] + *pa;
+      pa++ ;
+    }
+}
+
+ostream & operator << (ostream & os, const RMatrix & A) {
+  int n=A.Dim ;
+  double* pa=A.Vals ;
+  os << endl ;
+  for (int i = 0; i < n; i++) {
+    for (int j=0 ; j< n ; j++) {
+      os << (*(pa++)) << " " ;
+    }
+    os << endl ;
+  }
+  return os ;
+}
diff --git a/stogo/linalg.h b/stogo/linalg.h
new file mode 100644 (file)
index 0000000..8e889a7
--- /dev/null
@@ -0,0 +1,86 @@
+/*
+   Temporary implementation of vector and matrix classes
+   This is more or less borrowed from Serguei's program
+*/
+
+#ifndef LINALG_H
+#define LINALG_H
+
+#include <iostream>
+using namespace std;
+#include <math.h>         // for sqrt()
+#include <float.h>
+
+typedef const class RVector CRVector;
+typedef CRVector& RCRVector;
+typedef const class RMatrix CRMatrix ;
+typedef CRMatrix& RCRMatrix;
+
+double eps() ;
+
+#define max(A,B)    ((A) > (B) ? (A):(B))
+#define min(A,B)    ((A) < (B) ? (A):(B))
+
+/********************* Class RVector *********************/
+
+class RVector{
+protected:
+  double*  elements;  // array of values
+  int      len;       // size of array
+
+ public:
+  RVector() ;
+  RVector(int);       // Constructor
+  RVector(RCRVector); // copy constructor
+  ~RVector() { delete[] elements; elements=0 ; len=0; }
+
+  RCRVector operator=(double) ;
+  RCRVector operator=(RCRVector);
+
+  double & operator () (int i) const {return elements[i] ; }
+  double nrm2() ; // Euclidian norm
+
+  double *raw_data() { return elements; }
+  const double *raw_data_const() const { return elements; }
+
+  friend ostream & operator << (ostream &, const RVector &);
+
+  friend double norm2(RCRVector) ;
+  friend double normInf(RCRVector) ;
+  friend double dot(RCRVector, RCRVector) ;
+  friend void scal(double, RCRVector) ;
+  friend void copy(RCRVector, RCRVector) ;
+  friend void axpy(double, RCRVector, RCRVector) ;
+  friend void gemv(char,double, RCRMatrix, RCRVector, double, RCRVector);
+  friend void ger(double alpha, RCRVector, RCRVector, RCRMatrix);
+
+  int GetLength() const { return len; }; // get vector size
+};
+
+/******************* Class RMatrix *************************/
+
+class RMatrix
+{
+ protected:
+  double*  Vals; // array of values
+  int       Dim; // dimension
+
+ public:
+   RMatrix() ;
+   RMatrix(int); // dimension
+  ~RMatrix() { delete[] Vals;  Vals=0 ; Dim=0; }
+  RMatrix(RCRMatrix); // copy constructor
+  RCRMatrix operator=(double num) ;
+  RCRMatrix operator=(RCRMatrix) ; // (needed for template stuff)
+
+  double& operator()(int vidx,int hidx) ;
+  friend ostream & operator << (ostream &, const RMatrix &);
+
+  friend void gemv(char,double, RCRMatrix, RCRVector, double, RCRVector);
+  friend void ger(double alpha,RCRVector,RCRVector,RCRMatrix);
+
+  int       GetDim() { return Dim; }; // get dimension
+};
+
+#endif
diff --git a/stogo/local.cc b/stogo/local.cc
new file mode 100644 (file)
index 0000000..615a439
--- /dev/null
@@ -0,0 +1,322 @@
+
+/*
+   Local search - A trust region algorithm with BFGS update.
+*/
+
+#include <iostream>
+#include <stdlib.h>
+
+#include "stogo_config.h"
+#include "global.h"
+#include "local.h"
+#include "tools.h"
+
+int local(Trial &T, TBox &box, TBox &domain, double eps_cl, double *mgr,
+          Global &glob, int axis, RCRVector x_av) {
+
+  int k_max, info, outside ;
+  int k, i, good_enough, iTmp ;
+  int n=box.GetDim();
+
+  double maxgrad, delta, f, f_new, tmp ;
+  double alpha, gamma, beta, d2, s2, nom, den, ro ;
+  double nrm_sd, nrm_hn, snrm_hn, nrm_dl ;
+  RVector x(n), g(n), h_sd(n), h_dl(n), h_n(n), x_new(n), g_new(n) ;
+  RVector s(n),y(n),z(n),w(n) ; // Temporary vectors
+  RMatrix B(n), H(n) ;          // Hessian and it's inverse
+
+  k_max = max_iter*n ;
+  x=T.xvals ;
+
+#ifdef LS_DEBUG
+  cout << "Local Search, x=" << x << endl;
+#endif
+
+  if (box.OutsideBox(x, domain) != 0) {
+    cout << "Starting point is not inside the boundary. Exiting...\n" ;
+    exit(1) ;
+    return LS_Out ;
+  }
+
+  // Check if we are close to a stationary point located previously
+  if (box.CloseToMin(x, &tmp, eps_cl)) {
+#ifdef LS_DEBUG
+     cout << "Close to a previously located stationary point, exiting" << endl;
+#endif
+     T.objval=tmp;
+     return LS_Old ;
+   } 
+
+  // Initially B and H are equal to the identity matrix
+  B=0 ; H=0 ;
+  for (i=0 ; i<n ; i++) {
+    B(i,i)=1 ;
+    H(i,i)=1 ;
+  }
+
+  RVector g_av(x_av.GetLength());
+  if (axis==-1) {
+    f=glob.ObjectiveGradient(x,g,OBJECTIVE_AND_GRADIENT);
+  }
+  else {
+    x_av(axis)=x(0);
+    f=glob.ObjectiveGradient(x_av,g_av,OBJECTIVE_AND_GRADIENT);
+    g(0)=g_av(axis);
+  }
+  FC++;GC++;
+
+  if (axis == -1) {
+    // Skipping AV
+#ifdef INI3
+    // Elaborate scheme to initalize delta
+    delta=delta_coef*norm2(g) ;
+    copy(g,z) ;
+    axpy(1.0,x,z) ;
+    if (!box.InsideBox(z)) {
+      if (box.Intersection(x,g,z)==TRUE) {
+       axpy(-1.0,x,z) ;
+       delta=min(delta,delta_coef*norm2(z)) ;
+      }
+      else {
+       // Algorithm broke down, use INI1
+        delta = (1.0/7)*box.ShortestSide(&iTmp) ;
+      }
+    }
+#endif
+#ifdef INI2
+    // Use INI2 scheme
+    delta = box.ClosestSide(x)*delta_coef ;
+    if (delta<MacEpsilon)
+      // Patch to avoid trust region with radius close to zero
+      delta = (1.0/7)*box.ShortestSide(&iTmp) ;
+#endif
+#ifdef INI1
+    delta = delta_coef*box.ShortestSide(&iTmp) ;
+#endif
+  }
+  else {
+    // Use a simple scheme for the 1D minimization (INI1)
+    delta = (1.0/7.0)*box.ShortestSide(&iTmp) ;
+  }
+
+  k=0 ; good_enough = 0 ; info=LS_New ; outside=0 ;
+  maxgrad=*mgr ;
+  while (good_enough == 0) {
+    k++ ;
+    if (k>k_max) {
+#ifdef LS_DEBUG
+      cout << "Maximum number of iterations reached\n" ;
+#endif
+      info=LS_MaxIter ;
+      break ;
+    }
+
+    // Update maximal gradient value
+    maxgrad=max(maxgrad,normInf(g)) ;
+
+    // Steepest descent, h_sd = -g
+    copy(g,h_sd) ;
+    scal(-1.0,h_sd) ;
+    nrm_sd=norm2(h_sd) ;
+
+    if (nrm_sd < epsilon) {
+      // Stop criterion (gradient) fullfilled
+#ifdef LS_DEBUG
+      cout << "Gradient small enough" << endl ;
+#endif
+      good_enough = 1 ;
+      break ;
+    }
+
+    // Compute Newton step, h_n = -H*g
+    gemv('N',-1.0, H, g, 0.0, h_n) ;
+    nrm_hn = norm2(h_n) ;
+
+    if (nrm_hn < delta) {
+      // Pure Newton step
+      copy(h_n, h_dl) ;
+#ifdef LS_DEBUG
+      cout << "[Newton step]      " ;
+#endif
+    }
+    else {
+      gemv('N',1.0,B,g,0.0,z) ;
+      tmp=dot(g,z) ;
+      if (tmp==0) {
+       info = LS_Unstable ;
+       break ;
+      }
+      alpha=(nrm_sd*nrm_sd)/tmp ; // Normalization (N38,eq. 3.30)
+      scal(alpha,h_sd) ;
+      nrm_sd=fabs(alpha)*nrm_sd ;
+
+      if (nrm_sd >= delta) {
+       gamma = delta/nrm_sd ; // Normalization (N38, eq. 3.33)
+       copy(h_sd,h_dl) ;
+       scal(gamma,h_dl) ;
+#ifdef LS_DEBUG
+       cout << "[Steepest descent]  " ;
+#endif
+      }
+      else {
+       // Combination of Newton and SD steps
+       d2 = delta*delta ; 
+       copy(h_sd,s) ; 
+       s2=nrm_sd*nrm_sd ;
+       nom = d2 - s2 ;
+       snrm_hn=nrm_hn*nrm_hn ;
+       tmp = dot(h_n,s) ;
+        den = tmp-s2 + sqrt((tmp-d2)*(tmp-d2)+(snrm_hn-d2)*(d2-s2)) ;
+       if (den==0) {
+         info = LS_Unstable ;
+         break ;
+       }
+       // Normalization (N38, eq. 3.31)
+       beta = nom/den ; 
+       copy(h_n,h_dl) ;
+       scal(beta,h_dl) ;
+       axpy((1-beta),h_sd,h_dl) ;
+#ifdef LS_DEBUG
+       cout << "[Mixed step]        " ;
+#endif
+      }
+    }
+    nrm_dl=norm2(h_dl) ;
+
+    //x_new = x+h_dl ;
+    copy(x,x_new) ;
+    axpy(1.0,h_dl,x_new) ;
+
+    // Check if x_new is inside the box
+    iTmp=box.OutsideBox(x_new, domain) ;
+    if (iTmp == 1) {
+#ifdef LS_DEBUG
+      cout << "x_new is outside the box " << endl ;
+#endif
+      outside++ ;
+      if (outside>max_outside_steps) {
+       // Previous point was also outside, exit
+       break ;
+      }
+    }
+    else if (iTmp == 2) {
+#ifdef LS_DEBUG
+      cout << " x_new is outside the domain" << endl ;
+#endif
+      info=LS_Out ;
+      break ;
+    }
+    else {
+      outside=0 ;  
+    }
+
+    // Compute the gain
+    if (axis==-1)
+      f_new=glob.ObjectiveGradient(x_new,x_new,OBJECTIVE_ONLY);
+    else {
+      x_av(axis)=x_new(0);
+      f_new=glob.ObjectiveGradient(x_av,x_new,OBJECTIVE_ONLY);
+    }
+    FC++;
+    gemv('N',0.5,B,h_dl,0.0,z);
+    ro = (f_new-f) / (dot(g,h_dl) + dot(h_dl,z)); // Quadratic model
+    if (ro > 0.75) {
+      delta = delta*2;
+    }
+    if (ro < 0.25) {
+      delta = delta/3;
+    }
+    if (ro > 0) {
+      // Update the Hessian and it's inverse using the BFGS formula
+      if (axis==-1)
+       glob.ObjectiveGradient(x_new,g_new,GRADIENT_ONLY);
+      else {
+       x_av(axis)=x_new(0);
+       glob.ObjectiveGradient(x_av,g_av,GRADIENT_ONLY);
+       g_new(0)=g_av(axis);
+      }
+      GC++;
+
+      // y=g_new-g
+      copy(g_new,y);
+      axpy(-1.0,g,y);
+
+      // Check curvature condition
+      alpha=dot(y,h_dl);
+      if (alpha <= sqrt(MacEpsilon)*nrm_dl*norm2(y)) {
+#ifdef LS_DEBUG
+       cout << "Curvature condition violated " ;
+#endif
+      }
+      else {
+       // Update Hessian
+       gemv('N',1.0,B,h_dl,0.0,z) ; // z=Bh_dl
+       beta=-1/dot(h_dl,z) ;
+       ger(1/alpha,y,y,B) ;
+       ger(beta,z,z,B) ;
+
+        // Update Hessian inverse
+        gemv('N',1.0,H,y,0.0,z) ; // z=H*y
+        gemv('T',1.0,H,y,0.0,w) ; // w=y'*H
+       beta=dot(y,z) ;
+       beta=(1+beta/alpha)/alpha ;
+
+       // It should be possible to do this updating more efficiently, by
+       // exploiting the fact that (h_dl*y'*H) = transpose(H*y*h_dl')
+       ger(beta,h_dl,h_dl,H) ;
+       ger(-1/alpha,z,h_dl,H) ;
+       ger(-1/alpha,h_dl,w,H) ;
+      }
+
+      if (nrm_dl < norm2(x)*epsilon) {
+       // Stop criterion (iteration progress) fullfilled
+#ifdef LS_DEBUG
+       cout << "Progress is marginal" ;
+#endif
+       good_enough = 1 ;
+      }
+
+      // Check if we are close to a stationary point located previously
+      if (box.CloseToMin(x_new, &f_new, eps_cl)) {
+       // Note that x_new and f_new may be overwritten on exit from CloseToMin
+#ifdef LS_DEBUG
+       cout << "Close to a previously located stationary point, exiting" << endl;
+#endif
+       info = LS_Old ;
+       good_enough = 1 ;
+      }
+
+      // Update x, g and f
+      copy(x_new,x) ; copy(g_new,g) ; f=f_new ;
+
+#ifdef LS_DEBUG
+      cout << " x=" << x << endl ;
+#endif
+
+    }
+    else {
+#ifdef LS_DEBUG
+      cout << "Step is no good, ro=" << ro << " delta=" << delta << endl ;
+#endif
+    }
+    
+  } // wend
+
+  // Make sure the routine returns correctly...
+  // Check if last iterate is outside the boundary
+  if (box.OutsideBox(x, domain) != 0) {
+    info=LS_Out; f=DBL_MAX;
+  }
+
+  if (info == LS_Unstable) {
+    cout << "Local search became unstable. No big deal but exiting anyway\n" ;
+    exit(1);
+  }
+
+  T.xvals=x ; T.objval=f ;
+  *mgr=maxgrad ;
+  if (outside>0)
+    return LS_Out ;
+  else
+    return info ;
+}
diff --git a/stogo/local.h b/stogo/local.h
new file mode 100644 (file)
index 0000000..21136dd
--- /dev/null
@@ -0,0 +1,25 @@
+/*
+    Definitions of various variables used in the local search
+*/
+
+#ifndef LOCAL_H
+#define LOCAL_H
+
+#include "tools.h"
+#include "global.h"
+
+extern int FC, GC ;
+
+// Results of local search
+enum {LS_Unstable, LS_MaxIter, LS_Old, LS_New,LS_Out} ;
+
+const double delta_coef = 1.0/2.0; // Initialization of trust region
+const double epsilon = 1.0E-4 ;    // Stopping criterion, var 1e-4
+const int max_outside_steps=1 ;    // Maximum number of steps outside the box
+const int max_iter=50 ;            // Max iterations = max_iter*dim. of problem
+
+extern double MacEpsilon ;   //  min {x >= 0 : 1 + x > 1}
+
+int local(Trial &, TBox &, TBox &, double, double*, Global&, int, RCRVector) ;
+
+#endif
diff --git a/stogo/prog.cc b/stogo/prog.cc
new file mode 100644 (file)
index 0000000..1cb941a
--- /dev/null
@@ -0,0 +1,301 @@
+/*
+       A simple program to test the global optimizer.
+*/
+
+#include "global.h"
+#include "tools.h"
+#include "testfun.h"
+
+#define STRLEN_MAX 80
+
+int main() {
+  bool AVfail, AVflag;
+  int testfnc, dim, axis, i;
+  double AVbest;
+  Pdom  Dom;
+  Pgrad Grad;
+  Pobj  Obj;
+
+  cout << "Global Optimizer v0.1" << endl;
+  cout << "1.  Rosenbrock (min=0)" << endl;
+  cout << "2.  McCormick (min=-1.91)" << endl;
+  cout << "3.  Box & Betts (min=0)" << endl;
+  cout << "4.  Paviani (min=-45.7)" << endl;
+  cout << "5.  Extended Beale (min=0)" << endl;
+  cout << "6.  Three-hump Camel (min=0)" << endl;
+  cout << "7.  Jacobsen & Pedersen (min=3)" << endl;
+  cout << "8.  Poly1 (min=0)" << endl;
+  cout << "9.  Six-hump Camel (min=-1.03) " << endl;
+  cout << "10. One-dimensional (min=-23.58)" << endl;
+  cout << "11. Kowalik (min=0.0003075)" << endl;
+  cout << "12. Hansen (min=-176.54)" << endl;
+  cout << "13. Shubert (min=-24.06)" << endl;
+  cout << "14. Griewank (min=0)" << endl;
+  cout << "15. Levy4 (min=-21.502)" << endl;
+  cout << "16. Levy5 (min=-11.504)" << endl;
+  cout << "17. Levy6 (min=-11.504)" << endl;
+  cout << "18. Levy7 (min=-11.504)" << endl;
+  cout << "19. Rastrigin (min=0)" << endl;
+  cout << "20. Trid" << endl;
+  cout << "21. Perm(4,50)" << endl;
+  cout << "22. Perm(4,0.5)" << endl;
+  cout << "23. Powersum(8,18,44,114)" << endl;
+  cout << "24. Schwefel (min=-12569.5)" << endl;
+  cout << "25. Shekel (-10.0)" << endl;
+  cout << endl << "Select test function :";
+  cin >> testfnc;
+
+  switch (testfnc) {
+  case 10:
+    dim=1;
+    Dom=Domain_OneDim;
+    Obj=Objective_OneDim;
+    Grad=Gradient_OneDim;
+    break;
+  case 1:
+    dim=2;
+    Dom=Domain_Rosenbrock;
+    Obj=Objective_Rosenbrock;
+    Grad=Gradient_Rosenbrock;
+    break;
+  case 8:
+    dim=2;
+    Dom=Domain_Poly1;
+    Obj=Objective_Poly1;
+    Grad=Gradient_Poly1;
+    break;
+  case 7:
+    dim=2;
+    Dom=Domain_Jacobsen;
+    Obj=Objective_Jacobsen;
+    Grad=Gradient_Jacobsen;
+    break;
+  case 6:
+    dim=2;
+    Dom=Domain_Camel3;
+    Obj=Objective_Camel3;
+    Grad=Gradient_Camel3;
+    break;
+  case 5:
+    dim=2;
+    Dom=Domain_Beale3;
+    Obj=Objective_Beale3;
+    Grad=Gradient_Beale3;
+    break;
+  case 4:
+    dim=10;
+    Dom=Domain_Paviani;
+    Obj=Objective_Paviani;
+    Grad=Gradient_Paviani;
+    break;
+  case 12:
+    dim=2;
+    Dom=Domain_Hansen;
+    Obj=Objective_Hansen;
+    Grad=Gradient_Hansen;
+    break;
+  case 13:
+    dim=2;
+    Dom=Domain_Shubert;
+    Obj=Objective_Shubert;
+    Grad=Gradient_Shubert;
+    break;
+  case 11:
+    dim=4;
+    Dom=Domain_Kowalik;
+    Obj=Objective_Kowalik;
+    Grad=Gradient_Kowalik;
+    break;
+  case 9:
+    dim=2;
+    Dom=Domain_Camel6;
+    Obj=Objective_Camel6;
+    Grad=Gradient_Camel6;
+    break;
+  case 2:
+    dim=2;
+    Dom=Domain_McCormick;
+    Obj=Objective_McCormick;
+    Grad=Gradient_McCormick;
+    break;
+  case 3:
+    dim=3;
+    Dom=Domain_BoxBetts; 
+    Obj=Objective_BoxBetts;
+    Grad=Gradient_BoxBetts; 
+    break;
+  case 14:
+    dim=10;
+    Dom=Domain_Griewank;
+    Obj=Objective_Griewank;
+    Grad=Gradient_Griewank;
+    break;
+  case 15:
+    dim=4;
+    Dom=Domain_Levy;
+    Obj=Objective_Levy;
+    Grad=Gradient_Levy;  
+    break;
+  case 16:
+    dim=5;
+    Dom=Domain_Levy; 
+    Obj=Objective_Levy; 
+    Grad=Gradient_Levy; 
+    break;
+ case 17:
+    dim=6;
+    Dom=Domain_Levy;
+    Obj=Objective_Levy;
+    Grad=Gradient_Levy;
+    break;
+ case 18:
+    dim=7;
+    Dom=Domain_Levy;
+    Obj=Objective_Levy;
+    Grad=Gradient_Levy;
+    break;
+  case 19:
+    cout << "Enter problem dimension ";
+    int rast_dim;   
+    cin >> rast_dim;
+    dim=rast_dim;
+    Dom=Domain_Rastrigin;   
+    Obj=Objective_Rastrigin;
+    Grad=Gradient_Rastrigin;
+    break;
+  case 20:
+    cout << "Enter problem dimension (two or larger) ";
+    int trid_dim;
+    cin >> trid_dim;
+    dim=trid_dim; 
+    Dom=Domain_Trid;
+    Obj=Objective_Trid;
+    Grad=Gradient_Trid;
+    break;
+  case 21:
+    dim=4;
+    Dom=Domain_Perm;
+    Obj=Objective_Perm_4_50;
+    Grad=Gradient_Perm_4_50;
+    break;
+ case 22:
+    dim=4;
+    Dom=Domain_Perm;
+    Obj=Objective_Perm_4_05;
+    Grad=Gradient_Perm_4_05;
+    break;
+ case 23:  
+    dim=4;
+    Dom=Domain_Powersum;
+    Obj=Objective_Powersum;
+    Grad=Gradient_Powersum;
+    break;
+ case 24:
+    cout << "Enter problem dimension ";
+    int schwef_dim;
+    cin >> schwef_dim;
+    dim=schwef_dim;
+    Dom=Domain_Schwefel;
+    Obj=Objective_Schwefel;
+    Grad=Gradient_Schwefel;
+    break;
+  case 25:
+    dim=4;
+    Dom=Domain_Shekel;
+    Obj=Objective_Shekel;
+    Grad=Gradient_Shekel;
+    break;
+  default:
+    cout << "Error : Function not defined" << endl;
+    exit(1);
+  }
+
+  cout << "Dimension=" <<dim << endl;
+  TBox D(dim);
+  Dom(D);
+  cout << "Domain=";
+  for (i=0; i<dim; i++)
+    cout << "[" << D.lb(i) << "," << D.ub(i) << "]";
+  cout << endl << endl;
+
+  GlobalParams params;
+  cout << "Enter time limit (seconds) ";
+  cin >> params.maxtime;
+  if (params.maxtime<1) {   
+    cout << "Warning: time limit set to 1 second\n";
+  }
+  params.maxeval = 0;
+  cout << "Use factory settings (y/n) ";
+  char str[STRLEN_MAX]; cin >> str;
+  if (str[0]=='y') {
+    params.det_pnts=2*dim+1; params.rnd_pnts=0;
+    params.eps_cl=0.1; params.rshift=0.3;
+    params.mu=1.0E-4; AVflag=FALSE;
+  }
+  else {
+    cout << "Number of deterministic points ";
+    cin >> params.det_pnts;
+    cout << "Numer of stochastic points ";
+    cin >> params.rnd_pnts;
+    cout << "Radius of attraction ";
+    cin >> params.eps_cl;
+    cout << "Parameter rshift ";
+    cin >> params.rshift;
+    cout << "Parameter mu ";
+    cin >> params.mu;
+    cout << "Use the AV initialization (y/n) ";
+    cin >> str;
+    if (str[0]=='y') AVflag=TRUE; else AVflag=FALSE;
+  }
+
+  Global Problem(D,Obj, Grad, params);
+  RVector x_av(dim);
+  if (AVflag==TRUE) {
+    cout << "Enter time limit for each coordinate direction (seconds) ";
+    cin >> params.maxtime;
+    if (params.maxtime<1) {
+      cout << "Warning: time limit set to 1 second\n";
+    }
+    params.maxeval = 0;
+    params.det_pnts=3;
+    TBox I(1);
+    Global AV(I, Obj, Grad, params);
+
+    x_av=0.0; AVfail=FALSE;
+    for (axis=0; axis<Problem.dim; axis++) {
+      cout << "### axis=" << axis << " ###" << endl;
+      I.lb=(D.lb)(axis); I.ub=(D.ub)(axis);
+      AV.SetDomain(I);
+      AV.ClearSolSet();
+      AV.Search(axis, x_av);
+
+      if (AV.NoMinimizers()) {
+       cout << "AV failed with axis=" << axis << endl;
+       AVfail=TRUE; break;
+      }
+    }
+
+    if (AVfail==FALSE) {
+      AVbest=AV.GetMinValue();
+      cout << "### AV Located x=" << x_av << " fbound=" << AVbest << endl;
+      RVector AVx(Problem.dim);
+      AVx=x_av;
+
+      // Use result from AV for new fbound
+      Problem.SetMinValue(AVbest);
+
+      // Add the best point found to the initial box (domain)
+      Problem.AddPoint(x_av, AVbest);      
+    }
+  }
+
+  // Perform the main search
+  cout << "### Starting main search ###" << endl;
+  Problem.Search(-1, x_av);
+
+  cout << "Optimization terminated. Current set of minimizers is" << endl;
+  if (Problem.NoMinimizers() && AVflag==FALSE)
+    cout << "### No improvement found ###" << endl;
+  else
+    Problem.DispMinimizers();
+}
diff --git a/stogo/rosen.h b/stogo/rosen.h
new file mode 100644 (file)
index 0000000..f8080a3
--- /dev/null
@@ -0,0 +1,18 @@
+#include "linalg.h" 
+#include "tools.h"
+#include "stogo_config.h"
+
+void Domain_Rosenbrock(RTBox box) {
+  box.lb=-10.0 ; box.ub=10.0;
+}
+
+double Objective_Rosenbrock(RCRVector x) {
+   double a=x(1)-x(0)*x(0);
+   double b=1-x(0);
+   return 100*a*a + b*b;
+}
+
+void Gradient_Rosenbrock(RCRVector x, RCRVector grad) {
+  grad(0)=200*(x(1)-x(0)*x(0))*(-2*x(0))-2*(1-x(0));
+  grad(1)=200*(x(1)-x(0)*x(0));   
+}
diff --git a/stogo/stogo.cc b/stogo/stogo.cc
new file mode 100644 (file)
index 0000000..1a5b32f
--- /dev/null
@@ -0,0 +1,60 @@
+// A C-callable front-end to the StoGO global-optimization library.
+//  -- Steven G. Johnson
+
+#include "stogo.h"
+#include "global.h"
+
+class MyGlobal : public Global {
+protected:
+  objective_func my_func;
+  void *my_data;
+
+public:
+
+  MyGlobal(RTBox D, GlobalParams P, objective_func func, void *data) : Global(D, 0, 0, P), my_func(func), my_data(data) {}
+  
+  virtual double ObjectiveGradient(RCRVector xy, RVector &grad, whichO which){
+    ++numeval;
+    switch (which) {
+    case GRADIENT_ONLY:
+    case OBJECTIVE_AND_GRADIENT:
+      return my_func(xy.GetLength(), xy.raw_data_const(), grad.raw_data(), my_data);
+    case OBJECTIVE_ONLY:
+      return my_func(xy.GetLength(), xy.raw_data_const(), NULL, my_data);
+    }
+  }
+};
+
+int stogo_minimize(int n,
+                  objective_func fgrad, void *data,
+                  double *x, double *fmin,
+                  const double *l, const double *u,
+                  long int maxeval, long int maxtime)
+{
+  GlobalParams params;
+
+  // FIXME: WTF do these parameters mean?
+  params.det_pnts=2*n+1; params.rnd_pnts=0;
+  params.eps_cl=0.1; params.rshift=0.3;
+  params.mu=1.0E-4;
+
+  params.maxtime = maxtime;
+  params.maxeval = maxeval;
+
+  TBox D(n);
+  for (int i = 0; i < n; ++i) {
+    D.lb(i) = l[i];
+    D.ub(i) = u[i];
+  }
+
+  MyGlobal Problem(D, params, fgrad, data);
+  RVector dummyvec(n);
+  Problem.Search(-1, dummyvec);
+
+  if (Problem.NoMinimizers())
+    return 0;
+  
+  *fmin = Problem.OneMinimizer(dummyvec);
+  for (int i = 0; i < n; ++i) x[i] = dummyvec(i);
+  return 1;
+}
diff --git a/stogo/stogo.h b/stogo/stogo.h
new file mode 100644 (file)
index 0000000..2e39580
--- /dev/null
@@ -0,0 +1,60 @@
+/* A C-callable front-end to the StoGO global-optimization library.
+   -- Steven G. Johnson   */
+
+#ifndef STOGO_H
+#define STOGO_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef double (*objective_func)(int n, const double *x, double *grad,
+                                void *data);
+
+/* search for the global minimum of the function fgrad(n, x, grad, data)
+   inside a simple n-dimensional hyperrectangle.
+
+   Input:
+
+       n: dimension of search space (number of decision variables)
+
+       fgrad: the objective function of the form fgrad(n, x, grad, data),
+              returning the objective function at x, where
+                 n: dimension of search space
+                x: pointer to array of length n, point to evaluate
+                 grad: if non-NULL, an array of length n which
+                       should on return be the gradient d(fgrad)/dx
+                       [ if NULL, grad should be ignored ]
+                 data: arbitrary data pointer, whose value is the
+                      data argument of stogo_minimize
+
+       data: arbitrary pointer to any auxiliary data needed by fgrad
+
+       l, u: arrays of length n giving the lower and upper bounds of the
+             search space
+
+       maxeval: if nonzero, a maximum number of fgrad evaluations
+
+       maxtime: if nonzero, a maximum time (in seconds)
+
+   Output:
+
+      fmin: the minimum value of the objective function found
+
+      x: pointer to array of length n, giving the location of the minimum
+
+   Return value: 0 if no minimum found, 1 otherwise.
+
+ */
+
+int stogo_minimize(int n,
+                   objective_func fgrad, void *data,
+                   double *x, double *fmin,
+                   const double *l, const double *u,
+                   long int maxeval, long int maxtime);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/stogo/stogo_config.h b/stogo/stogo_config.h
new file mode 100644 (file)
index 0000000..1e94f37
--- /dev/null
@@ -0,0 +1,19 @@
+#ifndef STOGO_CONFIG_H
+#define STOGO_CONFIG_H
+
+// Compiler flags to enable/disable various options in the code
+
+// To obtain extra information on the global search process:
+//#define GS_DEBUG
+
+// To obtain info on the local search process:
+//#define LS_DEBUG     
+
+// Initialization strategy in the local search (INI1,INI2 or INI3)
+#define INI2
+
+// Use the more pessimistic strategy for lower bound estimation
+//#define LB2
+
+
+#endif
diff --git a/stogo/testfun.h b/stogo/testfun.h
new file mode 100644 (file)
index 0000000..5aadafe
--- /dev/null
@@ -0,0 +1,539 @@
+#ifndef TESTFUN_H
+#define TESTFUN_H
+
+#include "linalg.h"
+#include "tools.h"
+#include "stogo_config.h"
+
+const double pi=fabs(acos(-1.));
+
+/* The Matrix a and vector c are needed in the Shekel function */
+static double a[10][4]={ { 4,4,4,4 } ,
+                        { 1,1,1,1 } ,
+                        { 8,8,8,8 } ,
+                        { 6,6,6,6 } ,
+                        { 3,7,3,7 } ,
+                        { 2,9,2,9 } ,
+                        { 5,5,3,3 } ,
+                        { 8,1,8,1 } ,
+                        { 6,2,6,2 } ,
+                        {7,3.6,7,3.6} };
+static double c[10]= { .1 , .2 , .2 , .4 , .4 , .6 , .3, .7 , .5 , .5 };
+
+void Domain_Shekel(RTBox box) {
+  box.lb=0.0 ; box.ub=10.0 ;
+}
+
+double Objective_Shekel(RCRVector x) {
+  int n=x.GetLength() ;
+  double R=0.0, S;
+  for(int i=0;i<10;i++) {
+    S=0;
+    for(int j=0;j<n;j++) S+=pow(x(j)-a[i][j],2);
+    R-=1/(S+c[i]);
+  }
+  return R;
+}
+
+void Gradient_Shekel(RCRVector x, RVector &grad) {
+  int n=x.GetLength() ;
+  double R;
+  for(int k=0;k<n;k++) {
+    R=0.0;
+    for(int i=0;i<10;i++) {
+      R+=(2.0*x(k)-2.0*a[i][k])/(pow(pow(x(0)-a[i][0],2.0)+pow(x(1)-a[i][1],2.0)
+                                +pow(x(2)-a[i][2],2.0)+pow(x(3)-a[i][3],2.0)+c[i],2.0));
+    }
+    grad(k)=R;
+  }
+}
+
+
+/******************** Unimodal functions ******************/
+void Domain_Rosenbrock(RTBox box) {
+  box.lb=-10.0 ; box.ub=10.0 ;
+}
+
+double Objective_Rosenbrock(RCRVector x) {
+   double a=x(1)-x(0)*x(0) ;
+   double b=1-x(0) ;
+   return 100*a*a + b*b ;
+}
+
+void Gradient_Rosenbrock(RCRVector x, RVector &grad) {
+  grad(0)=200*(x(1)-x(0)*x(0))*(-2*x(0))-2*(1-x(0)) ;
+  grad(1)=200*(x(1)-x(0)*x(0)) ;
+}
+
+
+void Domain_McCormick(RTBox box) {
+  box.lb(0)=-1.5 ; box.ub(0)=4.0 ;
+  box.lb(1)=-3.0 ; box.ub(1)=4.0 ;
+}
+
+double Objective_McCormick(RCRVector x) {
+  return sin(x(0)+x(1)) + pow(x(0)-x(1),2.0) - 1.5*x(0) + 2.5*x(1) + 1.0 ;
+}
+
+void Gradient_McCormick(RCRVector x, RVector &grad) {
+  grad(0)=cos(x(0)+x(1)) + 2*(x(0)-x(1)) - 1.5 ;
+  grad(1)=cos(x(0)+x(1)) - 2*(x(0)-x(1)) + 2.5 ;
+}
+
+
+void Domain_BoxBetts(RTBox box) {
+  box.lb(0)=0.9 ; box.ub(0)=1.2 ;
+  box.lb(1)=9.0 ; box.ub(1)=11.2 ;
+  box.lb(2)=0.9 ; box.ub(2)=1.2 ;
+}
+
+double Objective_BoxBetts(RCRVector x) {
+  double x0=x(0),x1=x(1),x2=x(2) ;
+  double sum=0.0 ;
+  for (int i=1 ; i<=10 ; i++)
+    sum+=pow(exp(-0.1*i*x0)-exp(-0.1*i*x1)-(exp(-0.1*i)-exp(-1.0*i))*x2,2.0);
+  return sum ;
+}
+
+void Gradient_BoxBetts(RCRVector x, RVector &grad) {
+  double x0=x(0),x1=x(1),x2=x(2) ;
+  double g0=0.0, g1=0.0, g2=0.0 ;
+  for (int i=1 ; i<=10 ; i++) {
+    g0 += -0.2*(exp(-0.1*i*x0)-exp(-0.1*i*x1)
+         -(exp(-0.1*i)-exp(-1.0*i))*x2)*i*exp(-0.1*i*x0);
+    g1 += 0.2*(exp(-0.1*i*x0)-exp(-0.1*i*x1)-(exp(-0.1*i)
+         -exp(-1.0*i))*x2)*i*exp(-0.1*i*x1);
+    g2 += 2.0*(exp(-0.1*i*x0)-exp(-0.1*i*x1)
+         -(exp(-0.1*i)-exp(-1.0*i))*x2)*(-exp(-0.1*i)+exp(-1.0*i));
+  }
+  grad(0)=g0 ; grad(1)=g1 ; grad(2)=g2 ;
+}
+
+
+void Domain_Paviani(RTBox box) {
+ box.lb=2.001 ; box.ub=9.999 ;
+}
+
+double Objective_Paviani(RCRVector x) { 
+  double a,b,sum=0.0, mul=1.0 ;
+  int n=x.GetLength() ;
+  for (int i=0 ; i<n ; i++) {
+    a=log(x(i)-2.0) ; b=log(10.0-x(i)) ;
+    sum+= a*a + b*b ;
+    mul*= x(i) ;
+  }
+  return sum - pow(mul,0.2) ;
+}
+
+void Gradient_Paviani(RCRVector x, RVector &grad) {
+  double sum, mul=1.0 ;
+  int n=10 ;
+  for (int i=0 ; i<n ; i++) {
+    mul*= x(i) ;
+  }
+
+  for (int j=0 ; j<n ; j++) {
+    sum=2*log(x(j)-2.0)/(x(j)-2.0) - 2*log(10.0-x(j))/(10.0-x(j)) ;
+    grad(j) = sum - 0.2*(mul/x(j))/pow(mul,0.8) ;
+  }
+}
+
+
+void Domain_Beale3(RTBox box) { // NB Make sure that there is only one minima
+ box.lb=-0.5 ; box.ub=4.0 ;
+}
+
+double Objective_Beale3(RCRVector x) {
+  double x1=x(0), x2=x(1) ;
+  double b1=1.500 + (x2 - 1.)*x1 ;
+  double b2=2.250 + (x2*x2 - 1.)*x1 ;
+  double b3=2.625 + (x2*x2*x2 - 1.)*x1 ;
+  return b1*b1 + b2*b2 + b3*b3 ;
+}
+
+void Gradient_Beale3(RCRVector x, RVector &grad) {
+  double x1=x(0), x2=x(1) ;
+  grad(0) = 2*(1.500 + (x2 - 1.)*x1)*(x2 - 1.)
+          + 2*(2.250 + (x2*x2 - 1.)*x1)*(x2*x2 - 1.)
+          + 2*(2.625 + (x2*x2*x2 - 1.)*x1)*(x2*x2*x2 - 1.) ;
+
+  grad(1) = 2*(1.500 + (x2 - 1.)*x1)*x1
+          + 4*(2.250 + (x2*x2  - 1.)*x1)*x2*x1
+          + 6*(2.625 + (x2*x2*x2 - 1.)*x1)*x2*x2*x1 ;
+}
+
+/************** Difficult unimodal problems ****************/
+void Domain_Trid(RTBox box) {
+  int n=(box.lb).GetLength() ;
+  double tmp=pow(n,2.0);
+  box.lb=-tmp ; box.ub=tmp ;  // [-n^2,n^2]
+}
+
+double Objective_Trid(RCRVector x) {
+  int n=x.GetLength();
+  double sum1=0.0, sum2=0.0;
+  for (int i=1 ; i<=n ; i++)
+    sum1+=pow(x(i-1)-1,2.0);
+  for (int i=2 ; i<=n ; i++)
+    sum2+=x(i-1)*x(i-2);
+  return sum1 - sum2;
+}
+
+void Gradient_Trid(RCRVector x, RVector &grad) {
+  int n=x.GetLength();
+  grad(0)=2*(x(0)-1)-x(1) ;
+  for (int i=1 ; i<=n-2 ; i++)
+    grad(i)=2*(x(i)-1) - (x(i-1) + x(i+1)) ;
+  grad(n-1)=2*(x(n-1)-1) - x(n-2) ;
+}
+
+/******************** Multimodal functions ****************/
+
+void Domain_OneDim(RTBox box) {
+ box.lb=-25.0 ; box.ub=25.0 ;
+}
+
+double Objective_OneDim(RCRVector x) {
+  return x(0)*sin(x(0)) ;
+}
+
+void Gradient_OneDim(RCRVector x, RVector &grad) {
+  grad(0) = x(0)*cos(x(0)) + sin(x(0)) ;
+}
+
+void Domain_Poly1(RTBox box) {
+// box.lb=-5.0 ; box.ub=5.0 ;
+ box.lb=-4.0;box.ub=4.01;
+}
+
+double Objective_Poly1(RCRVector x) {
+  double a=(x(0)*x(0) + x(1)*x(1) - 11) ;
+  double b=(x(0)+x(1)*x(1) - 7) ;
+  return a*a + b*b ;
+}
+
+void Gradient_Poly1(RCRVector x, RVector &grad) {
+  double a=(x(0)*x(0) + x(1)*x(1) - 11) ;
+  double b=(x(0)+x(1)*x(1) - 7) ;
+  grad(0) = 4*x(0)* a + 2*b ;
+  grad(1) = 4*x(1)* a + 4*x(1)*b ;
+}
+
+void Domain_Jacobsen(RTBox box) {
+ box.lb=-100.0 ; box.ub=100.0 ;
+}
+
+double Objective_Jacobsen(RCRVector x) {
+  double a=(x(0)*x(0) + x(1) - 11) ;
+  double b=(x(0)+x(1)*x(1) - 7) ;
+  return a*a + b*b + 3;
+}
+
+void Gradient_Jacobsen(RCRVector x, RVector &grad) {
+  double a=(x(0)*x(0) + x(1) - 11) ;
+  double b=(x(0)+x(1)*x(1) - 7) ;
+  grad(0) = 4*x(0)* a + 2*b ;
+  grad(1) = 2*x(1)* a + 4*x(1)*b ;
+}
+
+void Domain_Camel3(RTBox box) {
+ box.lb=-1.0 ; box.ub=2.0 ;
+}
+
+double Objective_Camel3(RCRVector x) {
+  double a=x(0)*x(0) ;
+  return 12*a - 6.3*a*a + a*a*a + 6*x(1)*(x(1)-x(0)) ;
+}
+
+void Gradient_Camel3(RCRVector x, RVector &grad) {
+  double a=x(0)*x(0) ;
+  grad(0) = 24*x(0) - 25.2*a*x(0) + 6*a*a*x(0) - 6*x(1) ;
+  grad(1) = 12*x(1) - 6*x(0) ;
+}
+
+
+void Domain_Hansen(RTBox box) {
+  box.lb=-10.0 ; box.ub=10.0 ;
+}
+
+double Objective_Hansen(RCRVector x) {
+  return (cos(1.0)+2.0*cos(x(0)+2.0)+3.0*cos(2.0*x(0)+3.0)+4.0*cos(3.0*x(0)
+        +4.0)+5.0*cos(4.0*x(0)+5.0))*(cos(2.0*x(1)+1.0)+2.0*cos(3.0*x(1)+2.0)   
+        +3.0*cos(4.0*x(1)+3.0)+4.0*cos(5.0*x(1)+4.0)+5.0*cos(6.0*x(1)+5.0));
+}
+
+void Gradient_Hansen(RCRVector x, RVector &grad) {
+   grad(0) = (-2.0*sin(x(0)+2.0)-6.0*sin(2.0*x(0)+3.0)-12.0*sin(3.0*x(0)+4.0)
+           -20.0*sin(4.0*x(0)+5.0))*(cos(2.0*x(1)+1.0)+2.0*cos(3.0*x(1)+2.0)
+           +3.0*cos(4.0*x(1)+3.0)+4.0*cos(5.0*x(1)+4.0)+5.0*cos(6.0*x(1)+5.0));
+
+   grad(1) = (cos(1.0)+2.0*cos(x(0)+2.0)+3.0*cos(2.0*x(0)+3.0)+4.0*cos(3.0*x(0)
+             +4.0)+5.0*cos(4.0*x(0)+5.0))*(-2.0*sin(2.0*x(1)+1.0)
+             -6.0*sin(3.0*x(1)+2.0)-12.0*sin(4.0*x(1)+3.0)
+             -20.0*sin(5.0*x(1)+4.0)-30.0*sin(6.0*x(1)+5.0));
+}
+
+
+void Domain_Shubert(RTBox box) {
+  box.lb=-10.0 ; box.ub=10.0 ;
+}
+
+double Objective_Shubert(RCRVector x) {
+   return -sin(2.0*x(0)+1.0)-2.0*sin(3.0*x(0)+2.0)-3.0*sin(4.0*x(0)+3.0)
+          -4.0*sin(5.0*x(0)+4.0)-5.0*sin(6.0*x(0)+5.0)-sin(2.0*x(1)+1.0)
+          -2.0*sin(3.0*x(1)+2.0)-3.0*sin(4.0*x(1)+3.0)-4.0*sin(5.0*x(1)+4.0)
+          -5.0*sin(6.0*x(1)+5.0);
+}
+
+void Gradient_Shubert(RCRVector x, RVector &grad) {
+   grad(0) = -2.0*cos(2.0*x(0)+1.0)-6.0*cos(3.0*x(0)+2.0)-12.0*cos(4.0*x(0)+3.0)
+             -20.0*cos(5.0*x(0)+4.0)-30.0*cos(6.0*x(0)+5.0);
+   grad(1) = -2.0*cos(2.0*x(1)+1.0)-6.0*cos(3.0*x(1)+2.0)-12.0*cos(4.0*x(1)+3.0)
+             -20.0*cos(5.0*x(1)+4.0)-30.0*cos(6.0*x(1)+5.0);
+}
+
+
+void Domain_Griewank(RTBox box) {
+  box.lb=-500 ; box.ub=700 ;
+}
+
+double Objective_Griewank(RCRVector x) {
+  double sum=0 ;
+  double prod=1 ;
+  for (int i=0 ; i<10 ; i++) {
+    sum+=x(i)*x(i) ;
+    prod*=cos(x(i)/sqrt(double(i+1))) ;
+  }
+  return sum/4000.0-prod + 1 ;
+}
+
+void Gradient_Griewank(RCRVector x, RVector &grad) {
+  double prod=1 ;
+  for (int i=0 ; i<10 ; i++) {
+    prod*=cos(x(i)/sqrt(double(i+1))) ;
+  }
+  for (int i=0 ; i<10 ; i++) {
+    grad(i)=x(i)/2000.0 + 1/sqrt(double(i+1))
+           *prod/cos(x(i)/sqrt(double(i+1)))*sin(x(i)/sqrt(double(i+1))) ;
+  }
+}
+
+
+void Domain_Levy(RTBox box) {
+  int n=(box.lb).GetLength() ;
+  switch (n) {
+  case 4 :
+    box.lb=-10.0 ; box.ub=10.0 ;
+    break ;
+  default:
+    box.lb=-5.0 ; box.ub=5.0 ;
+  }
+}
+
+double Objective_Levy(RCRVector x) {
+  int n=x.GetLength();
+  double sum=0.0;
+
+  for (int i=0 ; i<=n-2 ; i++)
+    sum+=pow(x(i)-1,2.0)*(1+pow(sin(3*pi*x(i+1)),2.0));
+  return pow(sin(3*pi*x(0)),2.0) + sum + (x(n-1)-1)*(1+pow(sin(2*pi*x(n-1)),2.0));
+}
+
+
+void Gradient_Levy(RCRVector x, RVector &grad) {
+  int n=x.GetLength();
+
+  grad(0)=6*sin(3*pi*x(0))*cos(3*pi*x(0))*pi + 2*(x(0)-1)*(1+pow(sin(3*pi*x(1)),2.0));
+
+  for (int i=1 ; i<=n-2 ; i++)
+    grad(i)=6*pow(x(i-1)-1,2.0)*sin(3*pi*x(i))*cos(3*pi*x(i))*pi
+      + 2*(x(i)-1)*(1+pow(sin(3*pi*x(i+1)),2.0)) ;
+
+  grad(n-1)=6*pow(x(n-2)-1,2.0)*sin(3*pi*x(n-1))*cos(3*pi*x(n-1))*pi
+    + 1 + pow(sin(2*pi*x(n-1)),2.0) + 4*(x(n-1)-1)*sin(2*pi*x(n-1))*cos(2*pi*x(n-1))*pi;
+}
+
+
+void Domain_Kowalik(RTBox box) {
+  box.lb=0 ; box.ub=5.0 ;
+}
+
+double Objective_Kowalik(RCRVector x) {
+  // First element in a and b is not used
+  double a[]={999,0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246};
+  double b[]={999,1./0.25,1./0.5,1.,1./2,1./4,1./6,1./8,1./10,1./12,1./14,1./16};
+  double x1=x(0),x2=x(1),x3=x(2),x4=x(3);
+
+  double s=0;
+  for (int i=1 ; i<=11 ; i++)
+    s+=pow(a[i]- (x1*(b[i]*b[i] + b[i]*x2))/(b[i]*b[i] + b[i]*x3 + x4),2.0) ;
+  return s;
+}
+
+void Gradient_Kowalik(RCRVector x, RVector &grad) {
+  double a[]={999,0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246};
+  double b[]={999,1./0.25,1./0.5,1.,1./2,1./4,1./6,1./8,1./10,1./12,1./14,1./16};
+  double x1=x(0),x2=x(1),x3=x(2),x4=x(3);
+  double g1=0,g2=0,g3=0,g4=0,tmp;
+  for (int i=1 ; i<=11 ; i++) {
+    tmp=a[i]- (x1*(b[i]*b[i] + b[i]*x2))/(b[i]*b[i] + b[i]*x3 + x4) ;
+    g1=g1-2*tmp*(b[i]*b[i]+b[i]*x2)/(b[i]*b[i]+b[i]*x3+x4) ;
+    g2=g2-2*tmp*(x1*b[i])/(b[i]*b[i]+b[i]*x3+x4) ;
+    g3=g3+2*tmp*x1*(b[i]*b[i]+b[i]*x2)*b[i]/pow(b[i]*b[i]+b[i]*x3+x4,2.0) ;
+    g4=g4+2*tmp*x1*(b[i]*b[i]+b[i]*x2)/pow(b[i]*b[i]+b[i]*x3+x4,2.0) ;
+  }
+  grad(0)=g1; grad(1)=g2; grad(2)=g3; grad(3)=g4;
+}
+
+
+void Domain_Camel6(RTBox box) {
+ box.lb=-5.0 ; box.ub=10.0 ;
+}
+
+double Objective_Camel6(RCRVector x) {
+  double x1=x(0),x2=x(1) ;
+  return 4.0*x1*x1-0.21E1*pow(x1,4.0)+pow(x1,6.0)/3+x1*x2-4.0*x2*x2 
+    + 4.0*pow(x2,4.0);
+}
+
+void Gradient_Camel6(RCRVector x, RVector &grad) {
+  double x1=x(0),x2=x(1) ;
+  grad(0) = 8.0*x1-0.84E1*x1*x1*x1+2.0*pow(x1,5.0)+x2;
+  grad(1) = x1-8.0*x2+16.0*x2*x2*x2;
+}
+
+/** Multimodal function with a huge number of local optima **/
+
+void Domain_Rastrigin(RTBox box) {
+  box.lb=-4.12 ; box.ub=6.12;
+}
+
+double Objective_Rastrigin(RCRVector x) {
+  double sum=0.0;
+  int n=x.GetLength();
+
+  for (int i=1 ; i<=n ; i++)
+     sum+=pow(x(i-1),2.0)-10*cos(2*pi*x(i-1))+10 ;
+  return sum ;
+}
+
+void Gradient_Rastrigin(RCRVector x, RVector &grad) {
+  int n=x.GetLength();
+
+  for (int i=1 ; i<=n ; i++)
+    grad(i-1)= 2*x(i-1) + 20*sin(2*pi*x(i-1))*pi ;
+}
+
+void Domain_Schwefel(RTBox box) {
+  box.lb=-500.0;box.ub=500.0;
+}
+
+double Objective_Schwefel(RCRVector x) {
+  double sum=0.0;
+  int n=x.GetLength();
+
+  for (int i=0 ; i<n ; i++)
+    sum+=x(i)*sin(sqrt(fabs(x(i))));
+  return -sum;
+}
+
+void Gradient_Schwefel(RCRVector x, RVector &grad) {
+  int n=x.GetLength();
+
+  for (int i=0; i<n ; i++) {
+    if (x(i)>=0)
+      grad(i)=-( sin(sqrt(x(i)))+sqrt(x(i))/2.0*cos(sqrt(x(i))) );
+    else
+      grad(i)=-( sin(sqrt(-x(i)))+sqrt(-x(i))/2.0*cos(sqrt(-x(i))) );
+  }
+}
+
+/******* Difficult multimodal problem, PERM(n) *********/
+
+void Domain_Perm(RTBox box) {
+  int n=(box.lb).GetLength();
+  box.lb=double(-n) ; box.ub=double(n) ;
+}
+
+double ObjPerm(RCRVector x, double beta) {
+  int n=x.GetLength();
+  double s1=0.0;
+  for (int k=1; k<=n ; k++) {
+    double s2=0.0;
+    for (int i=1 ; i<=n; i++)
+      //  s2+=(i^k+beta)*((x(i-1)/i)^k-1);
+      s2+=(pow(1.0*i,1.0*k)+beta)*(pow(x(i-1)/i,1.0*k)-1) ;
+    s1+=s2*s2;
+  }
+  return s1;
+}
+
+void GradPerm(RCRVector x, RVector &grad, double beta) {
+  int n=x.GetLength();
+  for (int j=1 ; j<=n ; j++) {
+    double s1=0.0;
+    for (int k=1 ; k<=n ; k++) {
+      double s2=0.0;
+      for (int i=1 ; i<=n; i++)
+       //      s2+=(i^k+beta)*((x(i-1)/i)^k-1);
+       s2+=(pow(1.0*i,1.0*k)+beta)*(pow(x(i-1)/i,1.0*k)-1);
+      //    s1+=2*s2*(j^k+beta)/(x(j-1)*j^k)*k*x(j-1)^k;
+      s1+=2*s2*(pow(1.0*j,1.0*k)+beta)/pow(1.0*j,1.0*k)*k*pow(x(j-1),k-1.0);
+    }
+    grad(j-1)=s1;
+  }
+}
+
+double Objective_Perm_4_50(RCRVector x) {
+  return ObjPerm(x,50);
+}
+
+void Gradient_Perm_4_50(RCRVector x, RVector &grad) {
+  GradPerm(x,grad,50);
+}
+
+double Objective_Perm_4_05(RCRVector x) {
+  return ObjPerm(x,0.5);
+}
+
+void Gradient_Perm_4_05(RCRVector x, RVector &grad) {
+  GradPerm(x,grad,0.5);
+}
+
+/******************* Powersum **************/
+void Domain_Powersum(RTBox box) {
+  box.lb=0.0 ; box.ub=4.0 ;
+}
+
+double Objective_Powersum(RCRVector x) {
+  int n=x.GetLength();
+  double b[]={8,18,44,114};
+
+  double s1=0.0;
+  for (int k=1 ; k<=n ; k++) {
+    double s2=0.0;
+    for (int i=1 ; i<=n ; i++)
+      s2+=pow(x(i-1),1.0*k);
+    s1+=pow(s2-b[k-1],2.0);
+  }
+  return s1;
+}
+
+void Gradient_Powersum(RCRVector x, RVector &grad) {
+  int n=x.GetLength();
+  double b[]={8,18,44,114};
+
+  for (int j=0 ; j<n ; j++) {
+    double s1=0.0;
+    for (int k=1 ; k<=n ; k++) {
+      double s2=0.0;
+      for (int i=1 ; i<=n ; i++)
+       s2+=pow(x(i-1),1.0*k);
+      s1+=2*(s2-b[k-1])*k*pow(x(j),k-1.0);
+    }
+    grad(j)=s1;
+  }
+}
+
+#endif
diff --git a/stogo/testros.cc b/stogo/testros.cc
new file mode 100644 (file)
index 0000000..7c0dea7
--- /dev/null
@@ -0,0 +1,25 @@
+#include "global.h"
+#include "tools.h"
+#include "rosen.h"
+
+int main() {
+  int dim=2;
+
+  Pdom Dom=Domain_Rosenbrock;
+  Pobj Obj=Objective_Rosenbrock;
+  Pgrad Grad=Gradient_Rosenbrock;
+
+  GlobalParams params;
+  params.det_pnts=2*dim+1; params.rnd_pnts=0;
+  params.eps_cl=0.1; params.rshift=0.3;
+  params.mu=1.0E-4; params.maxtime=5;
+
+  TBox D(dim);
+  Dom(D);
+  Global Problem(D,Obj, Grad, params);
+  RVector dummyvec(dim);
+  Problem.Search(-1, dummyvec);
+
+  cout << "Minimizers found\n";
+  Problem.DispMinimizers();
+}
diff --git a/stogo/tools.cc b/stogo/tools.cc
new file mode 100644 (file)
index 0000000..3def645
--- /dev/null
@@ -0,0 +1,423 @@
+
+#include <float.h>
+#include <iostream>
+
+#include "stogo_config.h"
+#include "tools.h"
+
+Trial::Trial():xvals(0) {
+  objval=DBL_MAX;
+}
+
+Trial::Trial(int n):xvals(n) {
+  objval=DBL_MAX;
+}
+
+Trial::Trial(RCTrial tr):xvals(tr.xvals) {
+  objval=tr.objval ;
+}
+
+RCTrial Trial::operator=(RCTrial tr) {
+  xvals=tr.xvals ;
+  objval=tr.objval ;
+  return *this ;
+}
+
+ostream & operator << (ostream & os, const Trial & T) {
+  os << T.xvals << "  " << "(" << T.objval << ")" << endl ;
+  return os;
+}
+
+/********************* Vanilla Box **********************/
+VBox::VBox():lb(0),ub(0) {} // Constructor
+
+VBox::VBox(int n):lb(n),ub(n) {} // Constructor
+
+VBox::VBox(RCVBox box):lb(box.lb), ub(box.ub) {} // Constructor
+
+RCVBox VBox::operator=(RCVBox box) {
+  // Copy: Box1=Box2
+  lb=box.lb; ub=box.ub;
+  return *this;
+}
+
+int VBox::GetDim() {
+  return lb.GetLength();
+}
+
+double VBox::Width(int i) {
+  // Returns the width of the i-th interval, i between [0,...,dim-1]
+  return ub(i)-lb(i);
+}
+
+void VBox::Midpoint(RCRVector x) {
+  // Returns the midpoint of Box in x
+  int n=GetDim();
+  for (int i=0 ; i<n ; i++)
+    x(i)=fabs(ub(i)-lb(i))/2 + lb(i);
+}
+
+ostream & operator << (ostream & os, const VBox & B) {
+  int n=(B.lb).GetLength() ;
+  for (int i=0 ; i<n ;i++)
+    os << '[' << B.lb(i) << "," << B.ub(i) << "]";
+  return os;
+}
+
+/************************ Trial Box ***********************/
+TBox::TBox() : VBox() {
+  // Constructor
+  fmin=DBL_MAX;
+}
+
+TBox::TBox(int n) : VBox(n) {
+  // Constructor
+  fmin=DBL_MAX;
+}
+
+TBox::TBox(RCTBox box) : VBox(box) {
+  // Constructor + Copy
+  fmin=box.fmin;
+  TList=box.TList ;
+}
+
+RCTBox TBox::operator=(RCTBox box) {
+  // Copy
+  // NB We should 'somehow' use VBox to copy lb and ub
+  // Note that assignment operators are _not_ inherited
+  lb=box.lb ; ub=box.ub ;
+  fmin=box.fmin ;
+  TList=box.TList ;
+  return *this ;
+}
+
+double TBox::GetMin() {
+  return fmin;
+}
+
+bool TBox::EmptyBox() {
+  // Returns TRUE if the list of Trials is empty
+  return TList.empty() ;
+}
+
+void TBox::AddTrial(RCTrial T) {
+  // Add a Trial to the (back of) TList
+  TList.push_back(T) ;
+  if (T.objval<fmin)
+    fmin=T.objval ;
+}
+
+void TBox::RemoveTrial(Trial &T) {
+  // Remove a trial from the (back of) box
+  T=TList.back() ;
+  TList.pop_back() ;
+}
+
+void TBox::GetLastTrial(Trial &T) {
+// Return a trial from the (back of) box
+  T=TList.back() ;
+}
+
+list<Trial>::const_iterator TBox::FirstTrial() {
+  return TList.begin();
+}
+
+list<Trial>::const_iterator TBox::LastTrial() {
+  return TList.end();
+}
+
+void TBox::GetTrial(list<Trial>::const_iterator itr, Trial &T) {
+  T.xvals=(*itr).xvals;
+  T.objval=(*itr).objval;
+}
+
+void TBox::ClearBox() {
+  TList.erase(TList.begin(), TList.end());
+  fmin=DBL_MAX;
+}
+
+bool TBox::CloseToMin(RVector &vec, double *objval, double eps_cl) {
+  // Returns TRUE if 'vec' is close to some of the trials in the box,
+  // in this case, 'vec' and 'objval' are overwritten by the Trial data
+  // otherwise 'vec' and 'objval' are not affected.
+  //
+  // Here "close" means norm2(T - some trial in box)<=eps_cl
+  //
+  // NB It might be better to accept Trial as argument instead of vector
+
+  int n=GetDim();
+  RVector x(n), y(n);
+  list<Trial>::const_iterator itr;
+  for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
+    y=vec ; // NB Should be possible to avoid this copying
+    x=(*itr).xvals;
+    axpy(-1,x,y);
+    if (norm2(y)<=eps_cl) {
+      vec=x;
+      *objval=(*itr).objval;
+      return TRUE;
+    }
+  }
+  return FALSE;
+}
+
+unsigned int TBox::NStationary() {
+  // Returns the number of trials in a box
+  return TList.size() ;
+}
+
+void TBox::split(RTBox B1, RTBox B2) {
+  list<Trial>::const_iterator itr;
+  double w,m,tmp;
+  double fm1=DBL_MAX, fm2=DBL_MAX;
+  int i, k, ns;
+  int n=GetDim();
+
+  B1.lb=lb; B1.ub=ub;
+  B2.lb=lb; B2.ub=ub;
+  w=LongestSide(&i);
+  ns=TList.size();
+  switch (ns) {
+  case 0: case 1:
+    // Bisection
+    w=ub(i)-lb(i); // Length of interval
+    m=lb(i)+w/2;   // Midpoint
+    B1.ub(i)=m; B2.lb(i)=m;
+    break;
+  default:
+    // More elaborate split when there are more than one trials in B
+    // See Serguies and Kaj's tech. report, page 11
+    // Compute the center point of all stationary points
+    RVector center(n), x(n), dispers(n);
+    center=0; dispers=0;
+    for ( itr = TList.begin(); itr != TList.end(); ++itr )
+      axpy(1.0, (*itr).xvals, center);
+    scal((double)(1.0/ns),center);
+
+    // Compute the relative deviations
+    for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
+      for (k = 0; k<n; k++) {
+       x=(*itr).xvals;
+       dispers(k)=dispers(k)+pow(center(k)-x(k),2.0);
+      }
+    }
+    scal((double)(1.0/ns),dispers);
+
+    // i=arg max(disp)
+    tmp=dispers(0);i=0;
+    for (k=1; k<n; k++) {
+      if (dispers(k)>tmp) {
+       tmp=dispers(k);i=k;
+      }
+    }
+    B1.ub(i)=center(i) ; B2.lb(i)=center(i);
+    break;
+  }
+
+  // Split the set of trials accordingly
+  for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
+    if ( B1.InsideBox((*itr).xvals) ) {
+      fm1=min(fm1,(*itr).objval);
+      B1.AddTrial(*itr);
+    }
+    else {
+      B2.AddTrial(*itr);
+      fm2=min(fm2,(*itr).objval);
+    }
+  }
+  // Set fmin of B1 and B2
+  B1.fmin=fm1 ; B2.fmin=fm2;
+}
+
+void TBox::dispTrials() {
+  // Display information about box
+#ifdef KCC
+  copy(TList.begin(), TList.end(), ostream_iterator<Trial,char>(cout));
+#else
+  copy(TList.begin(), TList.end(), ostream_iterator<Trial>(cout));
+#endif
+}
+
+ostream & operator << (ostream & os, const TBox & B) {
+  int n=(B.lb).GetLength() ;
+  for (int i=0 ; i<n ;i++)
+    os << '[' << B.lb(i) << "," << B.ub(i) << "]";
+  os << "   fmin= " << B.fmin << endl;
+  return os;
+}
+
+bool TBox::InsideBox(RCRVector x) {
+  // Returns TRUE if the point X lies inside BOX, FALSE otherwise
+  int n=GetDim();
+  for (int i=0 ; i<n ; i++)
+    if (x(i)<lb(i) || x(i)>ub(i)) return FALSE;
+  return TRUE;
+}
+
+int TBox::OutsideBox(RCRVector x, RCTBox domain) {
+  // The function returns
+  //    0 if X is inside both 'box' and 'domain'
+  //    1 if X is inside 'domain' but outside 'box'
+  //    2 if X is outside both 'domain' and 'box
+
+  int n=GetDim();
+  int ins_box=1, ins_dom=1, outs=999;
+  for (int i=0 ; i<n ; i++) {
+    if (x(i)<lb(i) || x(i)>ub(i))
+      ins_box=0;
+    if (x(i)<domain.lb(i) || x(i)>domain.ub(i)) {
+      ins_dom=0; break;
+    }
+  }
+  if (ins_box==1 && ins_dom==1)
+    outs=0;
+  if (ins_box==0 && ins_dom==1)
+    outs=1;
+  if (ins_box==0 && ins_dom==0)
+    outs=2;
+  if (outs==999) {
+    // Something has gone wrong!
+    cout << "Error in OutsideBox, exiting\n";
+    exit(1);
+  }
+  return outs;
+}
+
+double TBox::ShortestSide(int *idx) {
+  // Returns the shortest side of the box and it's index
+  int n=GetDim(),j=0;
+  double tmp=ub(0)-lb(0);
+  for (int i=1 ; i<n ;i++)
+    if ( (ub(i)-lb(i))<tmp ) {
+      tmp=ub(i)-lb(i);
+      j=i;
+    }
+  *idx=j;
+  return tmp;
+}
+
+double TBox::LongestSide(int *idx) {
+  // Returns the longest side of the box and it's index
+  int n=GetDim(),j=0;
+  double tmp=ub(0)-lb(0);
+  for (int i=1 ; i<n ;i++)
+    if ( (ub(i)-lb(i))>tmp ) {
+      tmp=ub(i)-lb(i);
+      j=i;
+    }
+  *idx=j;
+  return tmp;
+}
+
+double TBox::ClosestSide(RCRVector x) {
+  // Returns the shortest distance from point X to the box B.
+
+  //   Warning: The output of this functon is nonsense if the
+  //   point X lies outside B. Should we try to detect this case?
+  
+  double dist, tmp ;
+  int n=GetDim();
+  dist=DBL_MAX;
+  for (int i=0 ; i<n ; i++) {
+    tmp = min( x(i)-lb(i), ub(i)-x(i) );
+    dist = min(dist, tmp);
+  }
+  return dist;
+}
+
+double TBox::FarthestSide(RCRVector x) {
+  // Returns the longest distance from point X to the box B.
+  //   Same comment apply here as in ClosestSide(X)
+    
+  double dist, tmp;
+  int n=GetDim();
+  dist=DBL_MIN;
+  for (int i=0 ; i<n ; i++) {
+    tmp = max( x(i)-lb(i), ub(i)-x(i) );
+    dist = max(dist, tmp);
+  }
+  return dist;
+}
+
+bool TBox::Intersection(RCRVector x, RCRVector h, RCRVector z) {
+// Intersection of a line with a box
+// The line is described by a vector 'h' and a point 'x' and
+// the point of intersection is returned in 'z'
+// Note that h+x must lie outside of box
+//
+// Known problem:
+//   Due to round of errors the algorithm can fail to find an intersection
+//   The caller is notified and should act accordingly
+//
+//  The routine returns FALSE if no intersection was found, TRUE otherwise
+
+  int n=GetDim();
+  RVector tmpV(n);
+  bool done;
+  int i, j, k, isect;
+  double alpha, gamma;
+
+  i=0; done=FALSE;
+  while (i<n && done==FALSE) {
+    if (h(i)==0) {
+      z(i)=x(i);
+      break;
+    }
+    for (k=1; k<=2; k++) {
+      if (k==1)
+       alpha=lb(i);
+      else
+       alpha=ub(i);
+      gamma=(alpha-x(i))/h(i);
+      z(i)=alpha;
+      isect=1;
+      for (j=0; j<n; j++) {
+       if (j != i) {
+         z(j)=x(j)+gamma*h(j);
+         if (z(j)<lb(j) || z(j)>ub(j)) {
+           isect=0;
+           break;
+         }
+       }
+      }
+      copy(z,tmpV); axpy(-1.0,x,tmpV);  // tmpV=z-x
+      if (isect==1 && dot(tmpV,h)>0) {
+       done=TRUE; break;
+      }
+    }
+    i++;
+  }
+  return done;
+}
+
+double TBox::LowerBound(double maxgrad) {
+// Lower bound estimation
+  double lb=fmin ;
+  double f1,f2,est ;
+  list<Trial>::const_iterator itr1,itr2 ;
+
+  int n=GetDim();
+  RVector x1(n), x2(n) ;
+#ifndef LB2
+  for ( itr1 = TList.begin(); itr1 != TList.end(); ++itr1 ) {
+    itr2=itr1 ;
+    while (++itr2 != TList.end()) {
+      x1=(*itr1).xvals ; f1=(*itr1).objval ;
+      x2=(*itr2).xvals ; f2=(*itr2).objval ;
+      axpy(-1.0,x2,x1) ; // x1=x1-x2
+      est=0.5*(f1+f2-maxgrad*norm2(x1)) ;
+      lb=min(lb,est) ;
+      // cout << "est=" << est << " x1=" << x1 << " x2=" << x2 << endl ;
+    }
+  }
+#endif
+#ifdef LB2
+  for ( itr1 = TList.begin(); itr1 != TList.end(); ++itr1 ) {
+    // Use also max distance to border
+    x1=(*itr1).xvals; f1=(*itr1).objval;
+    est=f1-FarthestSide(x1)*maxgrad;
+    lb=min(lb,est);
+  }
+#endif
+  return lb ;
+}
diff --git a/stogo/tools.h b/stogo/tools.h
new file mode 100644 (file)
index 0000000..81d8dc3
--- /dev/null
@@ -0,0 +1,119 @@
+
+//   Various datastructures and functions used by the global optimizer
+
+#ifndef TOOLS_H
+#define TOOLS_H
+
+#include <float.h>
+#include <iostream>
+
+#include <algorithm>
+#include <iterator>
+#include <list>
+
+#include "linalg.h"
+
+#ifndef FALSE
+const int FALSE=(1==0); // boolean FALSE
+#endif
+#ifndef FALSE
+const int TRUE=(1==1); // boolean TRUE
+#endif
+
+typedef const class Trial CTrial;
+typedef CTrial& RCTrial;
+typedef CTrial* PCTrial;
+
+class Trial{
+public:
+  RVector xvals;
+  double objval;
+
+  Trial();
+  Trial(int);
+  Trial(RCTrial); // Copy constructor
+  RCTrial operator=(RCTrial) ; // assignment operator
+  friend ostream & operator << (ostream &, RCTrial) ;
+};
+
+class TrialGT : public unary_function<Trial, bool>
+// Predicate for Trial (needed for remove_if)
+{
+public:
+  explicit TrialGT(double val) : _val(val) {}
+  bool operator()(Trial& foo) { 
+    return foo.objval > _val;
+  }
+private:
+  double _val;
+};
+
+typedef class VBox& RVBox;
+typedef const class VBox CVBox;
+typedef CVBox* PCVBox;
+typedef CVBox& RCVBox;
+
+class VBox{
+public:
+  RVector lb, ub;
+
+  VBox();        // Construct a box
+  VBox(int);
+  VBox(RCVBox);  // Copy constructor
+  RCVBox operator=(RCVBox);      // assignment operator
+
+  int GetDim();                  // Returns the dimension of the box
+  double Width(int) ;            // Returns the width of the i-th interval
+  void Midpoint(RCRVector);      // Returns the midpoint
+
+  friend ostream & operator << (ostream &, const VBox &);
+};
+
+typedef class TBox& RTBox;
+typedef const class TBox CTBox;
+typedef CTBox* PCTBox;
+typedef CTBox& RCTBox;
+
+class TBox: public VBox {
+public:
+  double fmin;   // Smallest function value found so far
+  list<Trial> TList; // List of trials
+
+  TBox();        // Construct a box
+  TBox(int);
+  TBox(RCTBox);  // Copy constructor
+
+  RCTBox operator=(RCTBox);      // assignment operator
+
+  double GetMin();               // Returns 'fmin'
+  bool EmptyBox();               // Returns TRUE if Box contains no trials
+  void AddTrial(RCTrial);        // Add a trial to the (back of) box
+  void RemoveTrial(Trial &);     // Remove a trial from the (back of) box
+  void GetLastTrial(Trial &);    // Return a trial from the back of the box
+
+  list<Trial>::const_iterator FirstTrial();
+  list<Trial>::const_iterator LastTrial();
+
+  void GetTrial(list<Trial>::const_iterator, Trial&);
+  void ClearBox();            
+  bool CloseToMin(RVector&, double*, double);
+
+  unsigned int NStationary();      // Returns the number of stationary points
+
+  void split(RTBox, RTBox);     // Split a box
+  void dispTrials();
+
+  bool InsideBox(RCRVector);
+  int  OutsideBox(RCRVector, RCTBox);
+  double ShortestSide(int*);    // Returns the shortest side
+  double LongestSide(int*);     // Returns the longest side
+  double ClosestSide(RCRVector x);
+  double FarthestSide(RCRVector x);
+  bool Intersection(RCRVector, RCRVector, RCRVector);
+  double LowerBound(double);
+
+  bool operator<(const TBox & x) const {return (fmin>x.fmin);}
+  friend ostream & operator << (ostream &, const TBox &);
+};
+
+#endif
diff --git a/stogo/tst.cc b/stogo/tst.cc
new file mode 100644 (file)
index 0000000..b2721ac
--- /dev/null
@@ -0,0 +1,62 @@
+#include <unistd.h>
+
+#include "global.h"
+#include "tools.h"
+
+#include "linalg.h" 
+#include "tools.h"
+#include "stogo_config.h"
+
+void Domain_Mine(RTBox box) {
+  box.lb=-3.0 ; box.ub=3.0;
+}
+
+static int cnt = 0;
+double Objective_Mine(RCRVector xy) {
+  double x = xy(0);
+  double y = xy(1);
+  double f = ((x*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3) + x*y + (y*y)*(-4+4*(y*y)));
+  printf("feval:, %d, %g, %g, %g\n", ++cnt, x,y, f);
+  return f;
+}
+
+void Gradient_Mine(RCRVector xy, RVector &grad) {
+  double x = xy(0);
+  double y = xy(1);
+  grad(0) = /* df/dx */
+    ((2*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3)
+     + (x*x)*(-4.2*x+4*(x*(x*x))/3)
+     + y);
+  grad(1) = /* df/dy */
+    (x + (2*y)*(-4+4*(y*y)) + (y*y)*(8*(y)));
+}
+
+
+int main(int argc, char **argv) {
+  int dim=2;
+
+  Pdom Dom=Domain_Mine;
+  Pobj Obj=Objective_Mine;
+  Pgrad Grad=Gradient_Mine;
+
+  GlobalParams params;
+  params.det_pnts=2*dim+1; params.rnd_pnts=0;
+  params.eps_cl=0.1; params.rshift=0.3;
+  params.mu=1.0E-4; params.maxtime=0;
+
+  params.maxeval = argc < 2 ? 100 : atoi(argv[1]);
+
+  TBox D(dim);
+  Dom(D);
+  Global Problem(D,Obj, Grad, params);
+  RVector dummyvec(dim);
+  Problem.Search(-1, dummyvec);
+
+  cout << "Minimizers found\n";
+  Problem.DispMinimizers();
+
+  cout << "dummyvec: " << dummyvec << "\n";
+  
+  double val = Problem.OneMinimizer(dummyvec);
+  cout << "one minimizer: " << val << ": " << dummyvec << "\n";
+}
diff --git a/stogo/tstc.c b/stogo/tstc.c
new file mode 100644 (file)
index 0000000..c321b3e
--- /dev/null
@@ -0,0 +1,46 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "stogo.h"
+
+/* has two global minima at (0.09,-0.71) and (-0.09,0.71), plus
+   4 additional local minima */
+static int cnt=0;
+double tst_obj(int n, const double *xy, double *g, void *unused)
+{
+  double x, y, f;
+  x = xy[0];
+  y = xy[1];
+  f = ((x*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3) + x*y + (y*y)*(-4+4*(y*y)));
+  printf("feval:, %d, %g, %g, %g\n", ++cnt, x,y, f);
+  if (g) {
+       g[0] = /* df/dx */
+           ((2*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3)
+            + (x*x)*(-4.2*x+4*(x*(x*x))/3)
+            + y);
+       g[1] = /* df/dy */
+           (x + (2*y)*(-4+4*(y*y)) + (y*y)*(8*(y)));
+  }
+  return f;
+}
+
+int main(int argc, char **argv)
+{
+  int n = 2;
+  double x[2], l[2], u[2];
+  long int maxits = 0, maxtim = 0;
+  int info;
+  double fmin;
+
+  maxits = argc < 2 ? 100 : atoi(argv[1]);
+
+  l[0] = -3; l[1] = -3;
+  u[0] = 3; u[1] = 3;
+
+  info = stogo_minimize(n, tst_obj, NULL, x, &fmin, l, u, maxits, maxtim);
+
+  printf("min f = %g at (%g,%g) after %d evals, return value %d\n",
+        fmin, x[0], x[1], cnt, info);
+
+  return EXIT_SUCCESS;
+}