From 6a6bb48404ede142d1a1d6676dbe1b3e013e3723 Mon Sep 17 00:00:00 2001 From: stevenj Date: Sat, 1 Sep 2007 16:51:35 -0400 Subject: [PATCH] added first cut at Praxis, renamed constants darcs-hash:20070901205135-c8de0-23fc5194f42a4bb5c52b8ecb88ab2f404af2907f.gz --- Makefile.am | 4 +- api/Makefile.am | 2 +- api/nlopt.c | 99 +++- api/nlopt.h | 43 +- cdirect/cdirect.c | 19 +- cdirect/cdirect.h | 28 + configure.ac | 1 + praxis/Makefile.am | 6 + praxis/README | 15 + praxis/praxis.c | 1287 ++++++++++++++++++++++++++++++++++++++++++++ praxis/praxis.h | 18 + test/testopt.cpp | 8 +- 12 files changed, 1474 insertions(+), 56 deletions(-) create mode 100644 praxis/Makefile.am create mode 100644 praxis/README create mode 100644 praxis/praxis.c create mode 100644 praxis/praxis.h diff --git a/Makefile.am b/Makefile.am index 7e8c741..af4ce1c 100644 --- a/Makefile.am +++ b/Makefile.am @@ -3,13 +3,13 @@ lib_LTLIBRARIES = libnlopt.la ACLOCAL_AMFLAGS=-I ./m4 -SUBDIRS= util lbfgs subplex direct cdirect stogo api octave . test +SUBDIRS= util lbfgs subplex direct cdirect stogo praxis api octave . test EXTRA_DIST=COPYRIGHT autogen.sh nlopt.pc.in m4 libnlopt_la_SOURCES = libnlopt_la_LIBADD = lbfgs/liblbfgs.la subplex/libsubplex.la \ direct/libdirect.la cdirect/libcdirect.la stogo/libstogo.la \ -api/libapi.la util/libutil.la +praxis/libpraxis.la api/libapi.la util/libutil.la libnlopt_la_LDFLAGS = -no-undefined -version-info @SHARED_VERSION_INFO@ diff --git a/api/Makefile.am b/api/Makefile.am index 9d3dbf2..08bffdc 100644 --- a/api/Makefile.am +++ b/api/Makefile.am @@ -1,4 +1,4 @@ -AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/lbfgs -I$(top_srcdir)/subplex -I$(top_srcdir)/util +AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/lbfgs -I$(top_srcdir)/subplex -I$(top_srcdir)/praxis -I$(top_srcdir)/util include_HEADERS = nlopt.h noinst_LTLIBRARIES = libapi.la diff --git a/api/nlopt.c b/api/nlopt.c index 3e8f44c..b5c5031 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -41,15 +41,19 @@ void nlopt_version(int *major, int *minor, int *bugfix) /*************************************************************************/ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = { - "DIRECT (global)", - "DIRECT-L (global)", - "Randomized DIRECT-L (global)", - "Original DIRECT version (global)", - "Original DIRECT-L version (global)", - "Subplex (local)", - "StoGO (global)", - "StoGO with randomized search (global)", - "Low-storage BFGS (LBFGS) (local)" + "DIRECT (global, no-derivative)", + "DIRECT-L (global, no-derivative)", + "Randomized DIRECT-L (global, no-derivative)", + "Unscaled DIRECT (global, no-derivative)", + "Unscaled DIRECT-L (global, no-derivative)", + "Unscaled Randomized DIRECT-L (global, no-derivative)", + "Original DIRECT version (global, no-derivative)", + "Original DIRECT-L version (global, no-derivative)", + "Subplex (local, no-derivative)", + "StoGO (global, derivative-based)", + "StoGO with randomized search (global, derivative-based)", + "Principal-axis, praxis (local, no-derivative)", + "Low-storage BFGS (LBFGS) (local, derivative-based)" }; const char *nlopt_algorithm_name(nlopt_algorithm a) @@ -79,6 +83,7 @@ typedef struct { } nlopt_data; #include "subplex.h" +#include "praxis.h" static double f_subplex(int n, const double *x, void *data_) { @@ -115,6 +120,38 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) /*************************************************************************/ +/* for "hybrid" algorithms that combine global search with some + local search algorithm, most of the time we anticipate that people + will stick with the default local search algorithm, so we + don't add this as a parameter to nlopt_minimize. Instead, the user + can call nlopt_{set/get}_hybrid_local_algorithm to get/set the defaults. */ + +/* default local-search algorithms */ +static nlopt_algorithm local_search_alg_deriv = NLOPT_LD_LBFGS; +static nlopt_algorithm local_search_alg_nonderiv = NLOPT_LN_SUBPLEX; + +static int local_search_maxeval = -1; /* no maximum by default */ + +void nlopt_get_local_search_algorithm(nlopt_algorithm *deriv, + nlopt_algorithm *nonderiv, + int *maxeval) +{ + *deriv = local_search_alg_deriv; + *nonderiv = local_search_alg_nonderiv; + *maxeval = local_search_maxeval; +} + +void nlopt_set_local_search_algorithm(nlopt_algorithm deriv, + nlopt_algorithm nonderiv, + int maxeval) +{ + local_search_alg_deriv = deriv; + local_search_alg_nonderiv = nonderiv; + local_search_maxeval = maxeval; +} + +/*************************************************************************/ + /* same as nlopt_minimize, but xtol_abs is required to be non-NULL */ static nlopt_result nlopt_minimize_( nlopt_algorithm algorithm, @@ -161,22 +198,31 @@ static nlopt_result nlopt_minimize_( stop.start = nlopt_seconds(); switch (algorithm) { - case NLOPT_GLOBAL_DIRECT: - case NLOPT_GLOBAL_DIRECT_L: - case NLOPT_GLOBAL_DIRECT_L_RANDOMIZED: + case NLOPT_GN_DIRECT: + case NLOPT_GN_DIRECT_L: + case NLOPT_GN_DIRECT_L_RAND: return cdirect(n, f, f_data, lb, ub, x, fmin, &stop, 0.0, - (algorithm != NLOPT_GLOBAL_DIRECT) - + 3 * (algorithm == NLOPT_GLOBAL_DIRECT_L_RANDOMIZED ? 2 : (algorithm != NLOPT_GLOBAL_DIRECT)) - + 9 * (algorithm == NLOPT_GLOBAL_DIRECT_L_RANDOMIZED ? 1 : (algorithm != NLOPT_GLOBAL_DIRECT))); + (algorithm != NLOPT_GN_DIRECT) + + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT)) + + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT))); + + case NLOPT_GN_DIRECT_NOSCAL: + case NLOPT_GN_DIRECT_L_NOSCAL: + case NLOPT_GN_DIRECT_L_RAND_NOSCAL: + return cdirect_unscaled(n, f, f_data, lb, ub, x, fmin, + &stop, 0.0, + (algorithm != NLOPT_GN_DIRECT) + + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT)) + + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT))); - case NLOPT_GLOBAL_ORIG_DIRECT: - case NLOPT_GLOBAL_ORIG_DIRECT_L: + case NLOPT_GN_ORIG_DIRECT: + case NLOPT_GN_ORIG_DIRECT_L: switch (direct_optimize(f_direct, &d, n, lb, ub, x, fmin, maxeval, -1, 0.0, 0.0, pow(xtol_rel, (double) n), -1.0, stop.fmin_max, 0.0, NULL, - algorithm == NLOPT_GLOBAL_ORIG_DIRECT + algorithm == NLOPT_GN_ORIG_DIRECT ? DIRECT_ORIGINAL : DIRECT_GABLONSKY)) { case DIRECT_INVALID_BOUNDS: @@ -200,15 +246,15 @@ static nlopt_result nlopt_minimize_( break; } - case NLOPT_GLOBAL_STOGO: - case NLOPT_GLOBAL_STOGO_RANDOMIZED: + case NLOPT_GD_STOGO: + case NLOPT_GD_STOGO_RAND: if (!stogo_minimize(n, f, f_data, x, fmin, lb, ub, &stop, - algorithm == NLOPT_GLOBAL_STOGO + algorithm == NLOPT_GD_STOGO ? 0 : 2*n)) return NLOPT_FAILURE; break; - case NLOPT_LOCAL_SUBPLEX: { + case NLOPT_LN_SUBPLEX: { int iret; double *scale = (double *) malloc(sizeof(double) * n); if (!scale) return NLOPT_OUT_OF_MEMORY; @@ -238,7 +284,14 @@ static nlopt_result nlopt_minimize_( break; } - case NLOPT_LOCAL_LBFGS: { + case NLOPT_LN_PRAXIS: { + double t0 = xtol_rel, macheps = 1e-14; + double h0 = 0.1; + *fmin = praxis_(&t0, &macheps, &h0, &n, x, f_subplex, &d); + break; + } + + case NLOPT_LD_LBFGS: { int iret, *nbd = (int *) malloc(sizeof(int) * n); if (!nbd) return NLOPT_OUT_OF_MEMORY; for (i = 0; i < n; ++i) { diff --git a/api/nlopt.h b/api/nlopt.h index 8f9e30d..5038062 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -11,22 +11,36 @@ typedef double (*nlopt_func)(int n, const double *x, void *func_data); typedef enum { - /* non-gradient algorithms */ + /* Naming conventions: - NLOPT_GLOBAL_DIRECT = 0, - NLOPT_GLOBAL_DIRECT_L, - NLOPT_GLOBAL_DIRECT_L_RANDOMIZED, - NLOPT_GLOBAL_ORIG_DIRECT, - NLOPT_GLOBAL_ORIG_DIRECT_L, + NLOPT_{G/L}{D/N}_* + = global/local derivative/no-derivative optimization, + respectively + + *_RAND algorithms involve some randomization. - NLOPT_LOCAL_SUBPLEX, + *_NOSCAL algorithms are *not* scaled to a unit hypercube + (i.e. they are sensitive to the units of x) + */ - /* gradient-based algorithms */ + NLOPT_GN_DIRECT = 0, + NLOPT_GN_DIRECT_L, + NLOPT_GN_DIRECT_L_RAND, + NLOPT_GN_DIRECT_NOSCAL, + NLOPT_GN_DIRECT_L_NOSCAL, + NLOPT_GN_DIRECT_L_RAND_NOSCAL, - NLOPT_GLOBAL_STOGO, - NLOPT_GLOBAL_STOGO_RANDOMIZED, + NLOPT_GN_ORIG_DIRECT, + NLOPT_GN_ORIG_DIRECT_L, - NLOPT_LOCAL_LBFGS, + NLOPT_LN_SUBPLEX, + + NLOPT_GD_STOGO, + NLOPT_GD_STOGO_RAND, + + NLOPT_LD_LBFGS, + + NLOPT_LN_PRAXIS, NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ } nlopt_algorithm; @@ -61,6 +75,13 @@ extern void nlopt_srand_time(void); extern void nlopt_version(int *major, int *minor, int *bugfix); +extern void nlopt_get_local_search_algorithm(nlopt_algorithm *deriv, + nlopt_algorithm *nonderiv, + int *maxeval); +extern void nlopt_set_local_search_algorithm(nlopt_algorithm deriv, + nlopt_algorithm nonderiv, + int maxeval); + #ifdef __cplusplus } /* extern "C" */ #endif /* __cplusplus */ diff --git a/cdirect/cdirect.c b/cdirect/cdirect.c index fae9316..e1d2636 100644 --- a/cdirect/cdirect.c +++ b/cdirect/cdirect.c @@ -45,6 +45,7 @@ typedef struct { 1: Gablonsky DIRECT-L (pick one pt, if equal pts) 2: ~ 1, but pick points randomly if equal pts ... 2 seems to suck compared to just picking oldest pt */ + const double *lb, *ub; nlopt_stopping *stop; /* stopping criteria */ nlopt_func f; void *f_data; @@ -440,7 +441,7 @@ static nlopt_result divide_good_rects(params *p) /***************************************************************************/ /* lexographic sort order (d,f,age) of hyper-rects, for red-black tree */ -static int hyperrect_compare(double *a, double *b) +int cdirect_hyperrect_compare(double *a, double *b) { if (a[0] < b[0]) return -1; if (a[0] > b[0]) return +1; @@ -482,7 +483,7 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data, p.hull = 0; p.age = 0; - rb_tree_init(&p.rtree, hyperrect_compare); + rb_tree_init(&p.rtree, cdirect_hyperrect_compare); p.work = (double *) malloc(sizeof(double) * (2*n)); if (!p.work) goto done; @@ -532,15 +533,9 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data, coordinates to a unit hypercube ... we do this simply by wrapping cdirect() around cdirect_unscaled(). */ -typedef struct { - nlopt_func f; - void *f_data; - double *x; - const double *lb, *ub; -} uf_data; -static double uf(int n, const double *xu, double *grad, void *d_) +double cdirect_uf(int n, const double *xu, double *grad, void *d_) { - uf_data *d = (uf_data *) d_; + cdirect_uf_data *d = (cdirect_uf_data *) d_; double f; int i; for (i = 0; i < n; ++i) @@ -559,7 +554,7 @@ nlopt_result cdirect(int n, nlopt_func f, void *f_data, nlopt_stopping *stop, double magic_eps, int which_alg) { - uf_data d; + cdirect_uf_data d; nlopt_result ret; const double *xtol_abs_save; int i; @@ -576,7 +571,7 @@ nlopt_result cdirect(int n, nlopt_func f, void *f_data, } xtol_abs_save = stop->xtol_abs; stop->xtol_abs = d.x + 3*n; - ret = cdirect_unscaled(n, uf, &d, d.x+n, d.x+2*n, x, fmin, stop, + ret = cdirect_unscaled(n, cdirect_uf, &d, d.x+n, d.x+2*n, x, fmin, stop, magic_eps, which_alg); stop->xtol_abs = xtol_abs_save; for (i = 0; i < n; ++i) diff --git a/cdirect/cdirect.h b/cdirect/cdirect.h index f101cda..41148a6 100644 --- a/cdirect/cdirect.h +++ b/cdirect/cdirect.h @@ -23,6 +23,34 @@ extern nlopt_result cdirect(int n, nlopt_func f, void *f_data, nlopt_stopping *stop, double magic_eps, int which_alg); +extern nlopt_result cdirect_hybrid(int n, nlopt_func f, void *f_data, + const double *lb, const double *ub, + double *x, + double *fmin, + nlopt_stopping *stop, + nlopt_algorithm local_alg, + int local_maxeval, + int randomized_div); + +extern nlopt_result cdirect_hybrid_unscaled(int n, nlopt_func f, void *f_data, + const double *lb, const double *ub, + double *x, + double *fmin, + nlopt_stopping *stop, + nlopt_algorithm local_alg, + int local_maxeval, + int randomized_div); + +/* internal routines and data structures: */ +extern int cdirect_hyperrect_compare(double *a, double *b); +typedef struct { + nlopt_func f; + void *f_data; + double *x; + const double *lb, *ub; +} cdirect_uf_data; +extern double cdirect_uf(int n, const double *xu, double *grad, void *d_); + #ifdef __cplusplus } /* extern "C" */ #endif /* __cplusplus */ diff --git a/configure.ac b/configure.ac index 7503eda..534abab 100644 --- a/configure.ac +++ b/configure.ac @@ -139,6 +139,7 @@ AC_CONFIG_FILES([ stogo/Makefile lbfgs/Makefile subplex/Makefile + praxis/Makefile test/Makefile ]) diff --git a/praxis/Makefile.am b/praxis/Makefile.am new file mode 100644 index 0000000..9034ee1 --- /dev/null +++ b/praxis/Makefile.am @@ -0,0 +1,6 @@ +AM_CPPFLAGS = -I$(top_srcdir)/util + +noinst_LTLIBRARIES = libpraxis.la +libpraxis_la_SOURCES = praxis.c praxis.h + +EXTRA_DIST = README diff --git a/praxis/README b/praxis/README new file mode 100644 index 0000000..c51dd42 --- /dev/null +++ b/praxis/README @@ -0,0 +1,15 @@ +praxis gradient-free local optimization via the "principal-axis method", +downloaded from Netlib: + + http://netlib.org/opt/praxis + +The original Fortran code was written by Richard Brent and made +available by the Stanford Linear Accelerator Center, dated 3/1/73. +Since this code contains no copyright statements and is dated prior to +1977, under US copyright law it is in the public domain (not copyrighted). + +Converted to C via f2c and cleaned up by Steven G. Johnson +(stevenj@alum.mit.edu). C version is licensed under the MIT license +(which is in the same spirit as public domain, but disclaims warranty +and is clearer legally). + diff --git a/praxis/praxis.c b/praxis/praxis.c new file mode 100644 index 0000000..5b9a52f --- /dev/null +++ b/praxis/praxis.c @@ -0,0 +1,1287 @@ +/* See README */ + +#include +#include "nlopt-util.h" +#include "praxis.h" + +/* Common Block Declarations */ + +struct global_s { + double fx, ldt, dmin__; + int nf, nl; +}; + +struct q_s { + double v[400] /* was [20][20] */, q0[20], q1[20], qa, qb, qc, qd0, + qd1, qf1; +}; + +/* Table of constant values */ + +static int pow_ii(int x, int n) /* compute x^n, n >= 0 */ +{ + int p = 1; + while (n > 0) { + if (n & 1) p *= x; + else { n >>= 1; x *= x; } + } + return p; +} + +static void minfit_(int *m, int *n, double *machep, + double *tol, double *ab, double *q); +static void min_(int n, int j, int nits, double *d2, double *x1, double *f1, int fk, praxis_func f, void *f_data, double *x, double *t, double *machep, double *h__, struct global_s *global_1, struct q_s *q_1); +static double flin_(int n, int j, double *l, praxis_func f, void *f_data, double *x, int *nf, struct q_s *q_1); +static void sort_(int *m, int *n, double *d__, double *v); +static void quad_(int *n, praxis_func f, void *f_data, double *x, double *t, + double *machep, double *h__, struct global_s *global_1, struct q_s *q_1); + +double praxis_(double *t0, double *machep, double *h0, + int *n, double *x, praxis_func f, void *f_data) +{ + /* System generated locals */ + int i__1, i__2, i__3; + double ret_val, d__1; + + /* Global */ + struct global_s global_1; + struct q_s q_1; + + /* Local variables */ + double d__[20], h__; + int i__, j, k; + double s, t, y[20], z__[20], f1; + int k2; + double m2, m4, t2, df, dn; + int kl, ii; + double sf; + int kt; + double sl; + int im1, km1; + double dni, lds; + int ktm; + double scbd; + int idim; + int illc; + int klmk; + double ldfac, large, small, value; + double vlarge; + double vsmall; + +/* LAST MODIFIED 3/1/73 */ + +/* PRAXIS RETURNS THE MINIMUM OF THE FUNCTION F(N,X) OF N VARIABLES */ +/* USING THE PRINCIPAL AXIS METHOD. THE GRADIENT OF THE FUNCTION IS */ +/* NOT REQUIRED. */ + +/* FOR A DESCRIPTION OF THE ALGORITHM, SEE CHAPTER SEVEN OF */ +/* "ALGORITHMS FOR FINDING ZEROS AND EXTREMA OF FUNCTIONS WITHOUT */ +/* CALCULATING DERIVATIVES" BY RICHARD P BRENT. */ + +/* THE PARAMETERS ARE: */ +/* T0 IS A TOLERANCE. PRAXIS ATTEMPTS TO RETURN PRAXIS=F(X) */ +/* SUCH THAT IF X0 IS THE TRUE LOCAL MINIMUM NEAR X, THEN */ +/* NORM(X-X0) < T0 + SQUAREROOT(MACHEP)*NORM(X). */ +/* MACHEP IS THE MACHINE PRECISION, THE SMALLEST NUMBER SUCH THAT */ +/* 1 + MACHEP > 1. MACHEP SHOULD BE 16.**-13 (ABOUT */ +/* 2.22D-16) FOR REAL*8 ARITHMETIC ON THE IBM 360. */ +/* H0 IS THE MAXIMUM STEP SIZE. H0 SHOULD BE SET TO ABOUT THE */ +/* MAXIMUM DISTANCE FROM THE INITIAL GUESS TO THE MINIMUM. */ +/* (IF H0 IS SET TOO LARGE OR TOO SMALL, THE INITIAL RATE OF */ +/* CONVERGENCE MAY BE SLOW.) */ +/* N (AT LEAST TWO) IS THE NUMBER OF VARIABLES UPON WHICH */ +/* THE FUNCTION DEPENDS. */ +/* X IS AN ARRAY CONTAINING ON ENTRY A GUESS OF THE POINT OF */ +/* MINIMUM, ON RETURN THE ESTIMATED POINT OF MINIMUM. */ +/* F(N,X) IS THE FUNCTION TO BE MINIMIZED. F SHOULD BE A REAL*8 */ +/* FUNCTION DECLARED EXTERNAL IN THE CALLING PROGRAM. */ +/* THE APPROXIMATING QUADRATIC FORM IS */ +/* Q(X') = F(N,X) + (1/2) * (X'-X)-TRANSPOSE * A * (X'-X) */ +/* WHERE X IS THE BEST ESTIMATE OF THE MINIMUM AND A IS */ +/* INVERSE(V-TRANSPOSE) * D * INVERSE(V) */ +/* (V(*,*) IS THE MATRIX OF SEARCH DIRECTIONS; D(*) IS THE ARRAY */ +/* OF SECOND DIFFERENCES). IF F HAS CONTINUOUS SECOND DERIVATIVES */ +/* NEAR X0, A WILL TEND TO THE HESSIAN OF F AT X0 AS X APPROACHES X0. */ + +/* IT IS ASSUMED THAT ON FLOATING-POINT UNDERFLOW THE RESULT IS SET */ +/* TO ZERO. */ +/* THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS AFTER */ +/* THE INITIALIZATION OF MACHINE DEPENDENT NUMBERS. */ + + +/* .....IF N>20 OR IF N<20 AND YOU NEED MORE SPACE, CHANGE '20' TO THE */ +/* LARGEST VALUE OF N IN THE NEXT CARD, IN THE CARD 'IDIM=20', AND */ +/* IN THE DIMENSION STATEMENTS IN SUBROUTINES MINFIT,MIN,FLIN,QUAD. */ + + +/* .....INITIALIZATION..... */ +/* MACHINE DEPENDENT NUMBERS: */ + + /* Parameter adjustments */ + --x; + + /* Function Body */ + small = *machep * *machep; + vsmall = small * small; + large = 1. / small; + vlarge = 1. / vsmall; + m2 = sqrt(*machep); + m4 = sqrt(m2); + +/* HEURISTIC NUMBERS: */ +/* IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF */ +/* POSSIBLE), THEN SET SCBD=10. OTHERWISE SET SCBD=1. */ +/* IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED, SET ILLC=TRUE. */ +/* OTHERWISE SET ILLC=FALSE. */ +/* KTM IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT BEFORE THE */ +/* ALGORITHM TERMINATES. KTM=4 IS VERY CAUTIOUS; USUALLY KTM=1 */ +/* IS SATISFACTORY. */ + + scbd = 1.; + illc = 0 /* false */; + ktm = 1; + + ldfac = .01; + if (illc) { + ldfac = .1; + } + kt = 0; + global_1.nl = 0; + global_1.nf = 1; + global_1.fx = f(*n, &x[1], f_data); + q_1.qf1 = global_1.fx; + t = small + fabs(*t0); + t2 = t; + global_1.dmin__ = small; + h__ = *h0; + if (h__ < t * 100) { + h__ = t * 100; + } + global_1.ldt = h__; +/* .....THE FIRST SET OF SEARCH DIRECTIONS V IS THE IDENTITY MATRIX..... */ + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + i__2 = *n; + for (j = 1; j <= i__2; ++j) { +/* L10: */ + q_1.v[i__ + j * 20 - 21] = 0.; + } +/* L20: */ + q_1.v[i__ + i__ * 20 - 21] = 1.; + } + d__[0] = 0.; + q_1.qd0 = 0.; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + q_1.q0[i__ - 1] = x[i__]; +/* L30: */ + q_1.q1[i__ - 1] = x[i__]; + } + +/* .....THE MAIN LOOP STARTS HERE..... */ +L40: + sf = d__[0]; + d__[0] = 0.; + s = 0.; + +/* .....MINIMIZE ALONG THE FIRST DIRECTION V(*,1). */ +/* FX MUST BE PASSED TO MIN BY VALUE. */ + value = global_1.fx; + min_(*n, 1, 2, d__, &s, &value, 0, f,f_data, &x[1], &t, + machep, &h__, &global_1, &q_1); + if (s > 0.) { + goto L50; + } + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { +/* L45: */ + q_1.v[i__ - 1] = -q_1.v[i__ - 1]; + } +L50: + if (sf > d__[0] * .9 && sf * .9 < d__[0]) { + goto L70; + } + i__1 = *n; + for (i__ = 2; i__ <= i__1; ++i__) { +/* L60: */ + d__[i__ - 1] = 0.; + } + +/* .....THE INNER LOOP STARTS HERE..... */ +L70: + i__1 = *n; + for (k = 2; k <= i__1; ++k) { + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L75: */ + y[i__ - 1] = x[i__]; + } + sf = global_1.fx; + if (kt > 0) { + illc = 1 /* true */; + } +L80: + kl = k; + df = 0.; + +/* .....A RANDOM STEP FOLLOWS (TO AVOID RESOLUTION VALLEYS). */ +/* PRAXIS ASSUMES THAT RANDOM RETURNS A RANDOM NUMBER UNIFORMLY */ +/* DISTRIBUTED IN (0,1). */ + + if (! illc) { + goto L95; + } + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { + s = (global_1.ldt * .1 + t2 * pow_ii(10, kt)) + * nlopt_urand(-.5,.5); + /* was: (random_(n) - .5); */ + z__[i__ - 1] = s; + i__3 = *n; + for (j = 1; j <= i__3; ++j) { +/* L85: */ + x[j] += s * q_1.v[j + i__ * 20 - 21]; + } +/* L90: */ + } + global_1.fx = (*f)(*n, &x[1], f_data); + ++global_1.nf; + +/* .....MINIMIZE ALONG THE "NON-CONJUGATE" DIRECTIONS V(*,K),...,V(*,N) */ + +L95: + i__2 = *n; + for (k2 = k; k2 <= i__2; ++k2) { + sl = global_1.fx; + s = 0.; + value = global_1.fx; + min_(*n, k2, 2, &d__[k2 - 1], &s, &value, 0, f,f_data, & + x[1], &t, machep, &h__, &global_1, &q_1); + if (illc) { + goto L97; + } + s = sl - global_1.fx; + goto L99; +L97: +/* Computing 2nd power */ + d__1 = s + z__[k2 - 1]; + s = d__[k2 - 1] * (d__1 * d__1); +L99: + if (df > s) { + goto L105; + } + df = s; + kl = k2; +L105: + ; + } + if (illc || df >= (d__1 = *machep * 100 * global_1.fx, fabs(d__1))) { + goto L110; + } + +/* .....IF THERE WAS NOT MUCH IMPROVEMENT ON THE FIRST TRY, SET */ +/* ILLC=TRUE AND START THE INNER LOOP AGAIN..... */ + + illc = 1 /* true */; + goto L80; +L110: + +/* .....MINIMIZE ALONG THE "CONJUGATE" DIRECTIONS V(*,1),...,V(*,K-1) */ + + km1 = k - 1; + i__2 = km1; + for (k2 = 1; k2 <= i__2; ++k2) { + s = 0.; + value = global_1.fx; + min_(*n, k2, 2, &d__[k2 - 1], &s, &value, 0, f,f_data, & + x[1], &t, machep, &h__, &global_1, &q_1); +/* L120: */ + } + f1 = global_1.fx; + global_1.fx = sf; + lds = 0.; + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { + sl = x[i__]; + x[i__] = y[i__ - 1]; + sl -= y[i__ - 1]; + y[i__ - 1] = sl; +/* L130: */ + lds += sl * sl; + } + lds = sqrt(lds); + if (lds <= small) { + goto L160; + } + +/* .....DISCARD DIRECTION V(*,KL). */ +/* IF NO RANDOM STEP WAS TAKEN, V(*,KL) IS THE "NON-CONJUGATE" */ +/* DIRECTION ALONG WHICH THE GREATEST IMPROVEMENT WAS MADE..... */ + + klmk = kl - k; + if (klmk < 1) { + goto L141; + } + i__2 = klmk; + for (ii = 1; ii <= i__2; ++ii) { + i__ = kl - ii; + i__3 = *n; + for (j = 1; j <= i__3; ++j) { +/* L135: */ + q_1.v[j + (i__ + 1) * 20 - 21] = q_1.v[j + i__ * 20 - 21]; + } +/* L140: */ + d__[i__] = d__[i__ - 1]; + } +L141: + d__[k - 1] = 0.; + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L145: */ + q_1.v[i__ + k * 20 - 21] = y[i__ - 1] / lds; + } + +/* .....MINIMIZE ALONG THE NEW "CONJUGATE" DIRECTION V(*,K), WHICH IS */ +/* THE NORMALIZED VECTOR: (NEW X) - (0LD X)..... */ + + value = f1; + min_(*n, k, 4, &d__[k - 1], &lds, &value, 1, f,f_data, &x[1], + &t, machep, &h__, &global_1, &q_1); + if (lds > 0.) { + goto L160; + } + lds = -lds; + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L150: */ + q_1.v[i__ + k * 20 - 21] = -q_1.v[i__ + k * 20 - 21]; + } +L160: + global_1.ldt = ldfac * global_1.ldt; + if (global_1.ldt < lds) { + global_1.ldt = lds; + } + t2 = 0.; + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L165: */ +/* Computing 2nd power */ + d__1 = x[i__]; + t2 += d__1 * d__1; + } + t2 = m2 * sqrt(t2) + t; + +/* .....SEE WHETHER THE LENGTH OF THE STEP TAKEN SINCE STARTING THE */ +/* INNER LOOP EXCEEDS HALF THE TOLERANCE..... */ + + if (global_1.ldt > t2 * .5f) { + kt = -1; + } + ++kt; + if (kt > ktm) { + goto L400; + } +/* L170: */ + } +/* .....THE INNER LOOP ENDS HERE. */ + +/* TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE IN A CURVED VALLEY. */ + +/* L171: */ + quad_(n, f,f_data, &x[1], &t, machep, &h__, &global_1, &q_1); + dn = 0.; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + d__[i__ - 1] = 1. / sqrt(d__[i__ - 1]); + if (dn < d__[i__ - 1]) { + dn = d__[i__ - 1]; + } +/* L175: */ + } + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + s = d__[j - 1] / dn; + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { +/* L180: */ + q_1.v[i__ + j * 20 - 21] = s * q_1.v[i__ + j * 20 - 21]; + } + } + +/* .....SCALE THE AXES TO TRY TO REDUCE THE CONDITION NUMBER..... */ + + if (scbd <= 1.) { + goto L200; + } + s = vlarge; + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { + sl = 0.; + i__1 = *n; + for (j = 1; j <= i__1; ++j) { +/* L182: */ + sl += q_1.v[i__ + j * 20 - 21] * q_1.v[i__ + j * 20 - 21]; + } + z__[i__ - 1] = sqrt(sl); + if (z__[i__ - 1] < m4) { + z__[i__ - 1] = m4; + } + if (s > z__[i__ - 1]) { + s = z__[i__ - 1]; + } +/* L185: */ + } + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { + sl = s / z__[i__ - 1]; + z__[i__ - 1] = 1. / sl; + if (z__[i__ - 1] <= scbd) { + goto L189; + } + sl = 1. / scbd; + z__[i__ - 1] = scbd; +L189: + i__1 = *n; + for (j = 1; j <= i__1; ++j) { +/* L190: */ + q_1.v[i__ + j * 20 - 21] = sl * q_1.v[i__ + j * 20 - 21]; + } +/* L195: */ + } + +/* .....CALCULATE A NEW SET OF ORTHOGONAL DIRECTIONS BEFORE REPEATING */ +/* THE MAIN LOOP. */ +/* FIRST TRANSPOSE V FOR MINFIT: */ + +L200: + i__2 = *n; + for (i__ = 2; i__ <= i__2; ++i__) { + im1 = i__ - 1; + i__1 = im1; + for (j = 1; j <= i__1; ++j) { + s = q_1.v[i__ + j * 20 - 21]; + q_1.v[i__ + j * 20 - 21] = q_1.v[j + i__ * 20 - 21]; +/* L210: */ + q_1.v[j + i__ * 20 - 21] = s; + } +/* L220: */ + } + +/* .....CALL MINFIT TO FIND THE SINGULAR VALUE DECOMPOSITION OF V. */ +/* THIS GIVES THE PRINCIPAL VALUES AND PRINCIPAL DIRECTIONS OF THE */ +/* APPROXIMATING QUADRATIC FORM WITHOUT SQUARING THE CONDITION */ +/* NUMBER..... */ + + minfit_(&idim, n, machep, &vsmall, q_1.v, d__); + +/* .....UNSCALE THE AXES..... */ + + if (scbd <= 1.) { + goto L250; + } + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { + s = z__[i__ - 1]; + i__1 = *n; + for (j = 1; j <= i__1; ++j) { +/* L225: */ + q_1.v[i__ + j * 20 - 21] = s * q_1.v[i__ + j * 20 - 21]; + } +/* L230: */ + } + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { + s = 0.; + i__1 = *n; + for (j = 1; j <= i__1; ++j) { +/* L235: */ +/* Computing 2nd power */ + d__1 = q_1.v[j + i__ * 20 - 21]; + s += d__1 * d__1; + } + s = sqrt(s); + d__[i__ - 1] = s * d__[i__ - 1]; + s = 1 / s; + i__1 = *n; + for (j = 1; j <= i__1; ++j) { +/* L240: */ + q_1.v[j + i__ * 20 - 21] = s * q_1.v[j + i__ * 20 - 21]; + } +/* L245: */ + } + +L250: + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { + dni = dn * d__[i__ - 1]; + if (dni > large) { + goto L265; + } + if (dni < small) { + goto L260; + } + d__[i__ - 1] = 1 / (dni * dni); + goto L270; +L260: + d__[i__ - 1] = vlarge; + goto L270; +L265: + d__[i__ - 1] = vsmall; +L270: + ; + } + +/* .....SORT THE EIGENVALUES AND EIGENVECTORS..... */ + + sort_(&idim, n, d__, q_1.v); + global_1.dmin__ = d__[*n - 1]; + if (global_1.dmin__ < small) { + global_1.dmin__ = small; + } + illc = 0 /* false */; + if (m2 * d__[0] > global_1.dmin__) { + illc = 1 /* true */; + } + +/* .....THE MAIN LOOP ENDS HERE..... */ + + goto L40; + +/* .....RETURN..... */ + +L400: + ret_val = global_1.fx; + return ret_val; +} /* praxis_ */ + +static void minfit_(int *m, int *n, double *machep, + double *tol, double *ab, double *q) +{ + /* System generated locals */ + int ab_dim1, ab_offset, i__1, i__2, i__3; + double d__1, d__2; + + /* Local variables */ + double c__, e[20], f, g, h__; + int i__, j, k, l; + double s, x, y, z__; + int l2, ii, kk, kt, ll2, lp1; + double eps, temp; + +/* ...AN IMPROVED VERSION OF MINFIT (SEE GOLUB AND REINSCH, 1969) */ +/* RESTRICTED TO M=N,P=0. */ +/* THE SINGULAR VALUES OF THE ARRAY AB ARE RETURNED IN Q AND AB IS */ +/* OVERWRITTEN WITH THE ORTHOGONAL MATRIX V SUCH THAT U.DIAG(Q) = AB.V, */ +/* WHERE U IS ANOTHER ORTHOGONAL MATRIX. */ +/* ...HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM... */ + /* Parameter adjustments */ + --q; + ab_dim1 = *m; + ab_offset = 1 + ab_dim1; + ab -= ab_offset; + + /* Function Body */ + if (*n == 1) { + goto L200; + } + eps = *machep; + g = 0.; + x = 0.; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + e[i__ - 1] = g; + s = 0.; + l = i__ + 1; + i__2 = *n; + for (j = i__; j <= i__2; ++j) { +/* L1: */ +/* Computing 2nd power */ + d__1 = ab[j + i__ * ab_dim1]; + s += d__1 * d__1; + } + g = 0.; + if (s < *tol) { + goto L4; + } + f = ab[i__ + i__ * ab_dim1]; + g = sqrt(s); + if (f >= 0.) { + g = -g; + } + h__ = f * g - s; + ab[i__ + i__ * ab_dim1] = f - g; + if (l > *n) { + goto L4; + } + i__2 = *n; + for (j = l; j <= i__2; ++j) { + f = 0.; + i__3 = *n; + for (k = i__; k <= i__3; ++k) { +/* L2: */ + f += ab[k + i__ * ab_dim1] * ab[k + j * ab_dim1]; + } + f /= h__; + i__3 = *n; + for (k = i__; k <= i__3; ++k) { +/* L3: */ + ab[k + j * ab_dim1] += f * ab[k + i__ * ab_dim1]; + } + } +L4: + q[i__] = g; + s = 0.; + if (i__ == *n) { + goto L6; + } + i__3 = *n; + for (j = l; j <= i__3; ++j) { +/* L5: */ + s += ab[i__ + j * ab_dim1] * ab[i__ + j * ab_dim1]; + } +L6: + g = 0.; + if (s < *tol) { + goto L10; + } + if (i__ == *n) { + goto L16; + } + f = ab[i__ + (i__ + 1) * ab_dim1]; +L16: + g = sqrt(s); + if (f >= 0.) { + g = -g; + } + h__ = f * g - s; + if (i__ == *n) { + goto L10; + } + ab[i__ + (i__ + 1) * ab_dim1] = f - g; + i__3 = *n; + for (j = l; j <= i__3; ++j) { +/* L7: */ + e[j - 1] = ab[i__ + j * ab_dim1] / h__; + } + i__3 = *n; + for (j = l; j <= i__3; ++j) { + s = 0.; + i__2 = *n; + for (k = l; k <= i__2; ++k) { +/* L8: */ + s += ab[j + k * ab_dim1] * ab[i__ + k * ab_dim1]; + } + i__2 = *n; + for (k = l; k <= i__2; ++k) { +/* L9: */ + ab[j + k * ab_dim1] += s * e[k - 1]; + } + } +L10: + y = (d__1 = q[i__], fabs(d__1)) + (d__2 = e[i__ - 1], fabs(d__2)); +/* L11: */ + if (y > x) { + x = y; + } + } +/* ...ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS... */ + ab[*n + *n * ab_dim1] = 1.; + g = e[*n - 1]; + l = *n; + i__1 = *n; + for (ii = 2; ii <= i__1; ++ii) { + i__ = *n - ii + 1; + if (g == 0.) { + goto L23; + } + h__ = ab[i__ + (i__ + 1) * ab_dim1] * g; + i__2 = *n; + for (j = l; j <= i__2; ++j) { +/* L20: */ + ab[j + i__ * ab_dim1] = ab[i__ + j * ab_dim1] / h__; + } + i__2 = *n; + for (j = l; j <= i__2; ++j) { + s = 0.; + i__3 = *n; + for (k = l; k <= i__3; ++k) { +/* L21: */ + s += ab[i__ + k * ab_dim1] * ab[k + j * ab_dim1]; + } + i__3 = *n; + for (k = l; k <= i__3; ++k) { +/* L22: */ + ab[k + j * ab_dim1] += s * ab[k + i__ * ab_dim1]; + } + } +L23: + i__3 = *n; + for (j = l; j <= i__3; ++j) { + ab[i__ + j * ab_dim1] = 0.; +/* L24: */ + ab[j + i__ * ab_dim1] = 0.; + } + ab[i__ + i__ * ab_dim1] = 1.; + g = e[i__ - 1]; +/* L25: */ + l = i__; + } +/* ...DIAGONALIZATION OF THE BIDIAGONAL FORM... */ +/* L100: */ + eps *= x; + i__1 = *n; + for (kk = 1; kk <= i__1; ++kk) { + k = *n - kk + 1; + kt = 0; +L101: + ++kt; + if (kt <= 30) { + goto L102; + } + e[k - 1] = 0.; + /* fprintf(stderr, "QR failed\n"); */ +L102: + i__3 = k; + for (ll2 = 1; ll2 <= i__3; ++ll2) { + l2 = k - ll2 + 1; + l = l2; + if ((d__1 = e[l - 1], fabs(d__1)) <= eps) { + goto L120; + } + if (l == 1) { + goto L103; + } + if ((d__1 = q[l - 1], fabs(d__1)) <= eps) { + goto L110; + } +L103: + ; + } +/* ...CANCELLATION OF E(L) IF L>1... */ +L110: + c__ = 0.; + s = 1.; + i__3 = k; + for (i__ = l; i__ <= i__3; ++i__) { + f = s * e[i__ - 1]; + e[i__ - 1] = c__ * e[i__ - 1]; + if (fabs(f) <= eps) { + goto L120; + } + g = q[i__]; +/* ...Q(I) = H = DSQRT(G*G + F*F)... */ + if (fabs(f) < fabs(g)) { + goto L113; + } + if (f != 0.) { + goto L112; + } else { + goto L111; + } +L111: + h__ = 0.; + goto L114; +L112: +/* Computing 2nd power */ + d__1 = g / f; + h__ = fabs(f) * sqrt(d__1 * d__1 + 1); + goto L114; +L113: +/* Computing 2nd power */ + d__1 = f / g; + h__ = fabs(g) * sqrt(d__1 * d__1 + 1); +L114: + q[i__] = h__; + if (h__ != 0.) { + goto L115; + } + g = 1.; + h__ = 1.; +L115: + c__ = g / h__; +/* L116: */ + s = -f / h__; + } +/* ...TEST FOR CONVERGENCE... */ +L120: + z__ = q[k]; + if (l == k) { + goto L140; + } +/* ...SHIFT FROM BOTTOM 2*2 MINOR... */ + x = q[l]; + y = q[k - 1]; + g = e[k - 2]; + h__ = e[k - 1]; + f = ((y - z__) * (y + z__) + (g - h__) * (g + h__)) / (h__ * 2 * y); + g = sqrt(f * f + 1.); + temp = f - g; + if (f >= 0.) { + temp = f + g; + } + f = ((x - z__) * (x + z__) + h__ * (y / temp - h__)) / x; +/* ...NEXT QR TRANSFORMATION... */ + c__ = 1.; + s = 1.; + lp1 = l + 1; + if (lp1 > k) { + goto L133; + } + i__3 = k; + for (i__ = lp1; i__ <= i__3; ++i__) { + g = e[i__ - 1]; + y = q[i__]; + h__ = s * g; + g *= c__; + if (fabs(f) < fabs(h__)) { + goto L123; + } + if (f != 0.) { + goto L122; + } else { + goto L121; + } +L121: + z__ = 0.; + goto L124; +L122: +/* Computing 2nd power */ + d__1 = h__ / f; + z__ = fabs(f) * sqrt(d__1 * d__1 + 1); + goto L124; +L123: +/* Computing 2nd power */ + d__1 = f / h__; + z__ = fabs(h__) * sqrt(d__1 * d__1 + 1); +L124: + e[i__ - 2] = z__; + if (z__ != 0.) { + goto L125; + } + f = 1.; + z__ = 1.; +L125: + c__ = f / z__; + s = h__ / z__; + f = x * c__ + g * s; + g = -x * s + g * c__; + h__ = y * s; + y *= c__; + i__2 = *n; + for (j = 1; j <= i__2; ++j) { + x = ab[j + (i__ - 1) * ab_dim1]; + z__ = ab[j + i__ * ab_dim1]; + ab[j + (i__ - 1) * ab_dim1] = x * c__ + z__ * s; +/* L126: */ + ab[j + i__ * ab_dim1] = -x * s + z__ * c__; + } + if (fabs(f) < fabs(h__)) { + goto L129; + } + if (f != 0.) { + goto L128; + } else { + goto L127; + } +L127: + z__ = 0.; + goto L130; +L128: +/* Computing 2nd power */ + d__1 = h__ / f; + z__ = fabs(f) * sqrt(d__1 * d__1 + 1); + goto L130; +L129: +/* Computing 2nd power */ + d__1 = f / h__; + z__ = fabs(h__) * sqrt(d__1 * d__1 + 1); +L130: + q[i__ - 1] = z__; + if (z__ != 0.) { + goto L131; + } + f = 1.; + z__ = 1.; +L131: + c__ = f / z__; + s = h__ / z__; + f = c__ * g + s * y; +/* L132: */ + x = -s * g + c__ * y; + } +L133: + e[l - 1] = 0.; + e[k - 1] = f; + q[k] = x; + goto L101; +/* ...CONVERGENCE: Q(K) IS MADE NON-NEGATIVE... */ +L140: + if (z__ >= 0.) { + goto L150; + } + q[k] = -z__; + i__3 = *n; + for (j = 1; j <= i__3; ++j) { +/* L141: */ + ab[j + k * ab_dim1] = -ab[j + k * ab_dim1]; + } +L150: + ; + } + return; +L200: + q[1] = ab[ab_dim1 + 1]; + ab[ab_dim1 + 1] = 1.; +} /* minfit_ */ + +static void min_(int n, int j, int nits, double * + d2, double *x1, double *f1, int fk, praxis_func f, void *f_data, double * + x, double *t, double *machep, double *h__, struct global_s *global_1, struct q_s *q_1) +{ + /* System generated locals */ + int i__1; + double d__1, d__2; + + /* Local variables */ + int i__, k; + double s, d1, f0, f2, m2, m4, t2, x2, fm; + int dz; + double xm, sf1, sx1; + double temp, small; + +/* ...THE SUBROUTINE MIN MINIMIZES F FROM X IN THE DIRECTION V(*,J) UNLESS */ +/* J IS LESS THAN 1, WHEN A QUADRATIC SEARCH IS MADE IN THE PLANE */ +/* DEFINED BY Q0,Q1,X. */ +/* D2 IS EITHER ZERO OR AN APPROXIMATION TO HALF F". */ +/* ON ENTRY, X1 IS AN ESTIMATE OF THE DISTANCE FROM X TO THE MINIMUM */ +/* ALONG V(*,J) (OR, IF J=0, A CURVE). ON RETURN, X1 IS THE DISTANCE */ +/* FOUND. */ +/* IF FK=.TRUE., THEN F1 IS FLIN(X1). OTHERWISE X1 AND F1 ARE IGNORED */ +/* ON ENTRY UNLESS FINAL FX IS GREATER THAN F1. */ +/* NITS CONTROLS THE NUMBER OF TIMES AN ATTEMPT WILL BE MADE TO HALVE */ +/* THE INTERVAL. */ + /* Parameter adjustments */ + --x; + + /* Function Body */ +/* Computing 2nd power */ + d__1 = *machep; + small = d__1 * d__1; + m2 = sqrt(*machep); + m4 = sqrt(m2); + sf1 = *f1; + sx1 = *x1; + k = 0; + xm = 0.; + fm = global_1->fx; + f0 = global_1->fx; + dz = *d2 < *machep; +/* ...FIND THE STEP SIZE... */ + s = 0.; + i__1 = n; + for (i__ = 1; i__ <= i__1; ++i__) { +/* L1: */ +/* Computing 2nd power */ + d__1 = x[i__]; + s += d__1 * d__1; + } + s = sqrt(s); + temp = *d2; + if (dz) { + temp = global_1->dmin__; + } + t2 = m4 * sqrt(fabs(global_1->fx) / temp + s * global_1->ldt) + m2 * + global_1->ldt; + s = m4 * s + *t; + if (dz && t2 > s) { + t2 = s; + } + t2 = t2 > small ? t2 : small; +/* Computing MIN */ + d__1 = t2, d__2 = *h__ * .01; + t2 = d__1 < d__2 ? d__1 : d__2; + if (! (fk) || *f1 > fm) { + goto L2; + } + xm = *x1; + fm = *f1; +L2: + if (fk && fabs(*x1) >= t2) { + goto L3; + } + temp = 1.; + if (*x1 < 0.) { + temp = -1.; + } + *x1 = temp * t2; + *f1 = flin_(n, j, x1, f,f_data, &x[1], &global_1->nf, q_1); +L3: + if (*f1 > fm) { + goto L4; + } + xm = *x1; + fm = *f1; +L4: + if (! dz) { + goto L6; + } +/* ...EVALUATE FLIN AT ANOTHER POINT AND ESTIMATE THE SECOND DERIVATIVE... */ + x2 = -(*x1); + if (f0 >= *f1) { + x2 = *x1 * 2.; + } + f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1); + if (f2 > fm) { + goto L5; + } + xm = x2; + fm = f2; +L5: + *d2 = (x2 * (*f1 - f0) - *x1 * (f2 - f0)) / (*x1 * x2 * (*x1 - x2)); +/* ...ESTIMATE THE FIRST DERIVATIVE AT 0... */ +L6: + d1 = (*f1 - f0) / *x1 - *x1 * *d2; + dz = 1 /* true */; +/* ...PREDICT THE MINIMUM... */ + if (*d2 > small) { + goto L7; + } + x2 = *h__; + if (d1 >= 0.) { + x2 = -x2; + } + goto L8; +L7: + x2 = d1 * -.5 / *d2; +L8: + if (fabs(x2) <= *h__) { + goto L11; + } + if (x2 <= 0.) { + goto L9; + } else { + goto L10; + } +L9: + x2 = -(*h__); + goto L11; +L10: + x2 = *h__; +/* ...EVALUATE F AT THE PREDICTED MINIMUM... */ +L11: + f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1); + if (k >= nits || f2 <= f0) { + goto L12; + } +/* ...NO SUCCESS, SO TRY AGAIN... */ + ++k; + if (f0 < *f1 && *x1 * x2 > 0.) { + goto L4; + } + x2 *= .5; + goto L11; +/* ...INCREMENT THE ONE-DIMENSIONAL SEARCH COUNTER... */ +L12: + ++global_1->nl; + if (f2 <= fm) { + goto L13; + } + x2 = xm; + goto L14; +L13: + fm = f2; +/* ...GET A NEW ESTIMATE OF THE SECOND DERIVATIVE... */ +L14: + if ((d__1 = x2 * (x2 - *x1), fabs(d__1)) <= small) { + goto L15; + } + *d2 = (x2 * (*f1 - f0) - *x1 * (fm - f0)) / (*x1 * x2 * (*x1 - x2)); + goto L16; +L15: + if (k > 0) { + *d2 = 0.; + } +L16: + if (*d2 <= small) { + *d2 = small; + } + *x1 = x2; + global_1->fx = fm; + if (sf1 >= global_1->fx) { + goto L17; + } + global_1->fx = sf1; + *x1 = sx1; +/* ...UPDATE X FOR LINEAR BUT NOT PARABOLIC SEARCH... */ +L17: + if (j == 0) { + return; + } + i__1 = n; + for (i__ = 1; i__ <= i__1; ++i__) { +/* L18: */ + x[i__] += *x1 * q_1->v[i__ + j * 20 - 21]; + } +} /* min_ */ + +static double flin_(int n, int j, double *l, praxis_func f, void *f_data, double *x, + int *nf, struct q_s *q_1) +{ + /* System generated locals */ + int i__1; + double ret_val; + + /* Local variables */ + int i__; + double t[20]; + +/* ...FLIN IS THE FUNCTION OF ONE REAL VARIABLE L THAT IS MINIMIZED */ +/* BY THE SUBROUTINE MIN... */ + /* Parameter adjustments */ + --x; + + /* Function Body */ + if (j == 0) { + goto L2; + } +/* ...THE SEARCH IS LINEAR... */ + i__1 = n; + for (i__ = 1; i__ <= i__1; ++i__) { +/* L1: */ + t[i__ - 1] = x[i__] + *l * q_1->v[i__ + j * 20 - 21]; + } + goto L4; +/* ...THE SEARCH IS ALONG A PARABOLIC SPACE CURVE... */ +L2: + q_1->qa = *l * (*l - q_1->qd1) / (q_1->qd0 * (q_1->qd0 + q_1->qd1)); + q_1->qb = (*l + q_1->qd0) * (q_1->qd1 - *l) / (q_1->qd0 * q_1->qd1); + q_1->qc = *l * (*l + q_1->qd0) / (q_1->qd1 * (q_1->qd0 + q_1->qd1)); + i__1 = n; + for (i__ = 1; i__ <= i__1; ++i__) { +/* L3: */ + t[i__ - 1] = q_1->qa * q_1->q0[i__ - 1] + q_1->qb * x[i__] + q_1->qc * + q_1->q1[i__ - 1]; + } +/* ...THE FUNCTION EVALUATION COUNTER NF IS INCREMENTED... */ +L4: + ++(*nf); + ret_val = f(n, t, f_data); + return ret_val; +} /* flin_ */ + +static void sort_(int *m, int *n, double *d__, double *v) +{ + /* System generated locals */ + int v_dim1, v_offset, i__1, i__2; + + /* Local variables */ + int i__, j, k; + double s; + int ip1, nm1; + +/* ...SORTS THE ELEMENTS OF D(N) INTO DESCENDING ORDER AND MOVES THE */ +/* CORRESPONDING COLUMNS OF V(N,N). */ +/* M IS THE ROW DIMENSION OF V AS DECLARED IN THE CALLING PROGRAM. */ + /* Parameter adjustments */ + v_dim1 = *m; + v_offset = 1 + v_dim1; + v -= v_offset; + --d__; + + /* Function Body */ + if (*n == 1) { + return; + } + nm1 = *n - 1; + i__1 = nm1; + for (i__ = 1; i__ <= i__1; ++i__) { + k = i__; + s = d__[i__]; + ip1 = i__ + 1; + i__2 = *n; + for (j = ip1; j <= i__2; ++j) { + if (d__[j] <= s) { + goto L1; + } + k = j; + s = d__[j]; +L1: + ; + } + if (k <= i__) { + goto L3; + } + d__[k] = d__[i__]; + d__[i__] = s; + i__2 = *n; + for (j = 1; j <= i__2; ++j) { + s = v[j + i__ * v_dim1]; + v[j + i__ * v_dim1] = v[j + k * v_dim1]; +/* L2: */ + v[j + k * v_dim1] = s; + } +L3: + ; + } +} /* sort_ */ + +static void quad_(int *n, praxis_func f, void *f_data, double *x, double *t, + double *machep, double *h__, struct global_s *global_1, struct q_s *q_1) +{ + /* System generated locals */ + int i__1; + double d__1; + + /* Local variables */ + int i__; + double l, s; + double value; + +/* ...QUAD LOOKS FOR THE MINIMUM OF F ALONG A CURVE DEFINED BY Q0,Q1,X... */ + /* Parameter adjustments */ + --x; + + /* Function Body */ + s = global_1->fx; + global_1->fx = q_1->qf1; + q_1->qf1 = s; + q_1->qd1 = 0.; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + s = x[i__]; + l = q_1->q1[i__ - 1]; + x[i__] = l; + q_1->q1[i__ - 1] = s; +/* L1: */ +/* Computing 2nd power */ + d__1 = s - l; + q_1->qd1 += d__1 * d__1; + } + q_1->qd1 = sqrt(q_1->qd1); + l = q_1->qd1; + s = 0.; + if (q_1->qd0 <= 0. || q_1->qd1 <= 0. || global_1->nl < *n * 3 * *n) { + goto L2; + } + value = q_1->qf1; + min_(*n, 0, 2, &s, &l, &value, 1, f,f_data, &x[1], t, machep, + h__, global_1, q_1); + q_1->qa = l * (l - q_1->qd1) / (q_1->qd0 * (q_1->qd0 + q_1->qd1)); + q_1->qb = (l + q_1->qd0) * (q_1->qd1 - l) / (q_1->qd0 * q_1->qd1); + q_1->qc = l * (l + q_1->qd0) / (q_1->qd1 * (q_1->qd0 + q_1->qd1)); + goto L3; +L2: + global_1->fx = q_1->qf1; + q_1->qa = 0.; + q_1->qb = q_1->qa; + q_1->qc = 1.; +L3: + q_1->qd0 = q_1->qd1; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + s = q_1->q0[i__ - 1]; + q_1->q0[i__ - 1] = x[i__]; +/* L4: */ + x[i__] = q_1->qa * s + q_1->qb * x[i__] + q_1->qc * q_1->q1[i__ - 1]; + } +} /* quad_ */ diff --git a/praxis/praxis.h b/praxis/praxis.h new file mode 100644 index 0000000..b6d0616 --- /dev/null +++ b/praxis/praxis.h @@ -0,0 +1,18 @@ +#ifndef PRAXIS_H +#define PRAXIS_H + +#ifdef __cplusplus +extern "C" +{ +#endif /* __cplusplus */ + +typedef double (*praxis_func)(int n, const double *x, void *f_data); + +double praxis_(double *t0, double *machep, double *h0, + int *n, double *x, praxis_func f, void *f_data); + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ + +#endif /* PRAXIS_H */ diff --git a/test/testopt.cpp b/test/testopt.cpp index 9bbe9a2..b93cf29 100644 --- a/test/testopt.cpp +++ b/test/testopt.cpp @@ -15,7 +15,7 @@ #include "nlopt-util.h" #include "testfuncs.h" -static nlopt_algorithm algorithm = NLOPT_GLOBAL_DIRECT; +static nlopt_algorithm algorithm = NLOPT_GN_DIRECT_L; static double ftol_rel = 0, ftol_abs = 0, xtol_rel = 0, xtol_abs = 0, fmin_max_delta = -HUGE_VAL; static int maxeval = 1000, iterations = 1, center_start = 0; static double maxtime = 0.0; @@ -51,12 +51,6 @@ static int test_function(int ifunc) return 0; } func = testfuncs[ifunc]; - if (!func.has_gradient && algorithm >= NLOPT_GLOBAL_STOGO) { - fprintf(stderr, - "testopt: A function with gradients is required for %s\n", - nlopt_algorithm_name(algorithm)); - return 0; - } x = (double *) malloc(sizeof(double) * func.n * 3); if (!x) { fprintf(stderr, "testopt: Out of memory!\n"); return 0; } -- 2.30.2