/* See README */
+#include <stdlib.h>
#include <math.h>
+#include <string.h>
#include "nlopt-util.h"
#include "praxis.h"
};
struct q_s {
- double v[400] /* was [20][20] */, q0[20], q1[20], qa, qb, qc, qd0,
- qd1, qf1;
+ double *v; /* size n x n */
+ double *q0, *q1, *t_flin; /* size n */
+ double qa, qb, qc, qd0, qd1, qf1;
+
+ double fbest, *xbest; /* size n */
+ nlopt_stopping *stop;
};
/* Table of constant values */
{
int p = 1;
while (n > 0) {
- if (n & 1) p *= x;
+ if (n & 1) { n--; 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)
+static void minfit_(int m, int n, double machep,
+ double *tol, double *ab, double *q, double *ework);
+static nlopt_result 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_old, 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, nlopt_result *ret);
+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_old, double machep, double *h__, struct global_s *global_1, struct q_s *q_1);
+
+nlopt_result praxis_(double t0, double machep, double h0,
+ int n, double *x, praxis_func f, void *f_data,
+ nlopt_stopping *stop, double *minf)
{
/* System generated locals */
int i__1, i__2, i__3;
- double ret_val, d__1;
+ nlopt_result ret = NLOPT_SUCCESS;
+ double d__1;
/* Global */
struct global_s global_1;
struct q_s q_1;
/* Local variables */
- double d__[20], h__;
+ double *d__, *y, *z__, *e_minfit, *prev_xbest; /* size n */
+ double prev_fbest;
+ double h__;
int i__, j, k;
- double s, t, y[20], z__[20], f1;
+ double s, t_old, f1;
int k2;
- double m2, m4, t2, df, dn;
+ double m2, m4, t2_old, df, dn;
int kl, ii;
double sf;
int kt;
double dni, lds;
int ktm;
double scbd;
- int idim;
int illc;
int klmk;
double ldfac, large, small, value;
double vlarge;
double vsmall;
+ double *work;
/* LAST MODIFIED 3/1/73 */
-
+/* Modified August 2007 by S. G. Johnson:
+ after conversion to C via f2c and some manual cleanup,
+ eliminating write statements, I additionally:
+ - modified the routine to use NLopt stop criteria
+ - allocate arrays dynamically to allow n > 20
+ - return the NLopt return code, where the min.
+ function value is now given by the parameter minf
+*/
/* 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. */
/* .....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. */
-
+ /* ...changed by S. G. Johnson, 2007, to use malloc */
/* .....INITIALIZATION..... */
/* MACHINE DEPENDENT NUMBERS: */
--x;
/* Function Body */
- small = *machep * *machep;
+ small = machep * machep;
vsmall = small * small;
large = 1. / small;
vlarge = 1. / vsmall;
- m2 = sqrt(*machep);
+ m2 = sqrt(machep);
m4 = sqrt(m2);
+ /* new: dynamic allocation of temporary arrays */
+ work = (double *) malloc(sizeof(double) * (n*n + n*9));
+ if (!work) return NLOPT_OUT_OF_MEMORY;
+ q_1.v = work;
+ q_1.q0 = q_1.v + n*n;
+ q_1.q1 = q_1.q0 + n;
+ q_1.t_flin = q_1.q1 + n;
+ q_1.xbest = q_1.t_flin + n;
+ d__ = q_1.xbest + n;
+ y = d__ + n;
+ z__ = y + n;
+ e_minfit = y + n;
+ prev_xbest = e_minfit + n;
+
/* HEURISTIC NUMBERS: */
/* IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF */
/* POSSIBLE), THEN SET SCBD=10. OTHERWISE SET SCBD=1. */
kt = 0;
global_1.nl = 0;
global_1.nf = 1;
- global_1.fx = f(*n, &x[1], f_data);
+ prev_fbest = q_1.fbest = global_1.fx = f(n, &x[1], f_data);
+ memcpy(q_1.xbest, &x[1], n*sizeof(double));
+ memcpy(prev_xbest, &x[1], n*sizeof(double));
+ stop->nevals++;
+ q_1.stop = stop;
q_1.qf1 = global_1.fx;
- t = small + fabs(*t0);
- t2 = t;
+ if (t0 > 0)
+ t_old = small + t0;
+ else {
+ t_old = 0;
+ for (i__ = 0; i__ < n; ++i__)
+ if (stop->xtol_abs[i__] > t_old)
+ t_old = stop->xtol_abs[i__];
+ t_old += small;
+ }
+ t2_old = t_old;
global_1.dmin__ = small;
- h__ = *h0;
- if (h__ < t * 100) {
- h__ = t * 100;
+ h__ = h0;
+ if (h__ < t_old * 100) {
+ h__ = t_old * 100;
}
global_1.ldt = h__;
/* .....THE FIRST SET OF SEARCH DIRECTIONS V IS THE IDENTITY MATRIX..... */
- i__1 = *n;
+ i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
- i__2 = *n;
+ i__2 = n;
for (j = 1; j <= i__2; ++j) {
/* L10: */
- q_1.v[i__ + j * 20 - 21] = 0.;
+ q_1.v[i__ + j * n - (n+1)] = 0.;
}
/* L20: */
- q_1.v[i__ + i__ * 20 - 21] = 1.;
+ q_1.v[i__ + i__ * n - (n+1)] = 1.;
}
d__[0] = 0.;
q_1.qd0 = 0.;
- i__1 = *n;
+ i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
q_1.q0[i__ - 1] = x[i__];
/* L30: */
/* .....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);
+ ret = min_(n, 1, 2, d__, &s, &value, 0, f,f_data, &x[1],
+ &t_old, machep, &h__, &global_1, &q_1);
+ if (ret != NLOPT_SUCCESS) goto done;
if (s > 0.) {
goto L50;
}
- i__1 = *n;
+ i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* L45: */
q_1.v[i__ - 1] = -q_1.v[i__ - 1];
if (sf > d__[0] * .9 && sf * .9 < d__[0]) {
goto L70;
}
- i__1 = *n;
+ i__1 = n;
for (i__ = 2; i__ <= i__1; ++i__) {
/* L60: */
d__[i__ - 1] = 0.;
/* .....THE INNER LOOP STARTS HERE..... */
L70:
- i__1 = *n;
+ i__1 = n;
for (k = 2; k <= i__1; ++k) {
- i__2 = *n;
+ i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
/* L75: */
y[i__ - 1] = x[i__];
if (! illc) {
goto L95;
}
- i__2 = *n;
+ i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
- s = (global_1.ldt * .1 + t2 * pow_ii(10, kt))
+ s = (global_1.ldt * .1 + t2_old * pow_ii(10, kt))
* nlopt_urand(-.5,.5);
/* was: (random_(n) - .5); */
z__[i__ - 1] = s;
- i__3 = *n;
+ i__3 = n;
for (j = 1; j <= i__3; ++j) {
/* L85: */
- x[j] += s * q_1.v[j + i__ * 20 - 21];
+ x[j] += s * q_1.v[j + i__ * n - (n+1)];
}
/* L90: */
}
- global_1.fx = (*f)(*n, &x[1], f_data);
+ 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;
+ 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);
+ ret = min_(n, k2, 2, &d__[k2 - 1], &s, &value, 0, f,f_data, &
+ x[1], &t_old, machep, &h__, &global_1, &q_1);
+ if (ret != NLOPT_SUCCESS) goto done;
if (illc) {
goto L97;
}
L105:
;
}
- if (illc || df >= (d__1 = *machep * 100 * global_1.fx, fabs(d__1))) {
+ if (illc || df >= (d__1 = machep * 100 * global_1.fx, fabs(d__1))) {
goto L110;
}
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);
+ ret = min_(n, k2, 2, &d__[k2 - 1], &s, &value, 0, f,f_data, &
+ x[1], &t_old, machep, &h__, &global_1, &q_1);
+ if (ret != NLOPT_SUCCESS) goto done;
/* L120: */
}
f1 = global_1.fx;
global_1.fx = sf;
lds = 0.;
- i__2 = *n;
+ i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
sl = x[i__];
x[i__] = y[i__ - 1];
i__2 = klmk;
for (ii = 1; ii <= i__2; ++ii) {
i__ = kl - ii;
- i__3 = *n;
+ 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];
+ q_1.v[j + (i__ + 1) * n - (n+1)] = q_1.v[j + i__ * n - (n+1)];
}
/* L140: */
d__[i__] = d__[i__ - 1];
}
L141:
d__[k - 1] = 0.;
- i__2 = *n;
+ i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
/* L145: */
- q_1.v[i__ + k * 20 - 21] = y[i__ - 1] / lds;
+ q_1.v[i__ + k * n - (n+1)] = 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);
+ ret = min_(n, k, 4, &d__[k - 1], &lds, &value, 1, f,f_data, &x[1],
+ &t_old, machep, &h__, &global_1, &q_1);
+ if (ret != NLOPT_SUCCESS) goto done;
if (lds > 0.) {
goto L160;
}
lds = -lds;
- i__2 = *n;
+ 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];
+ q_1.v[i__ + k * n - (n+1)] = -q_1.v[i__ + k * n - (n+1)];
}
L160:
global_1.ldt = ldfac * global_1.ldt;
if (global_1.ldt < lds) {
global_1.ldt = lds;
}
- t2 = 0.;
- i__2 = *n;
+ t2_old = 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_old += d__1 * d__1;
}
- t2 = m2 * sqrt(t2) + t;
+ t2_old = m2 * sqrt(t2_old) + t_old;
/* .....SEE WHETHER THE LENGTH OF THE STEP TAKEN SINCE STARTING THE */
/* INNER LOOP EXCEEDS HALF THE TOLERANCE..... */
- if (global_1.ldt > t2 * .5f) {
+ if (global_1.ldt > t2_old * .5f
+ && !nlopt_stop_f(stop, q_1.fbest, prev_fbest)
+ && !nlopt_stop_x(stop, q_1.xbest, prev_xbest)) {
kt = -1;
}
++kt;
if (kt > ktm) {
- goto L400;
+ if (nlopt_stop_f(stop, q_1.fbest, prev_fbest))
+ ret = NLOPT_FTOL_REACHED;
+ else if (nlopt_stop_x(stop, q_1.xbest, prev_xbest))
+ ret = NLOPT_XTOL_REACHED;
+ goto done;
}
+ prev_fbest = q_1.fbest;
+ memcpy(prev_xbest, q_1.xbest, n * sizeof(double));
/* 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);
+ quad_(n, f,f_data, &x[1], &t_old, machep, &h__, &global_1, &q_1);
dn = 0.;
- i__1 = *n;
+ i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
d__[i__ - 1] = 1. / sqrt(d__[i__ - 1]);
if (dn < d__[i__ - 1]) {
}
/* L175: */
}
- i__1 = *n;
+ i__1 = n;
for (j = 1; j <= i__1; ++j) {
s = d__[j - 1] / dn;
- i__2 = *n;
+ 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];
+ q_1.v[i__ + j * n - (n+1)] = s * q_1.v[i__ + j * n - (n+1)];
}
}
goto L200;
}
s = vlarge;
- i__2 = *n;
+ i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
sl = 0.;
- i__1 = *n;
+ 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];
+ sl += q_1.v[i__ + j * n - (n+1)] * q_1.v[i__ + j * n - (n+1)];
}
z__[i__ - 1] = sqrt(sl);
if (z__[i__ - 1] < m4) {
}
/* L185: */
}
- i__2 = *n;
+ i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
sl = s / z__[i__ - 1];
z__[i__ - 1] = 1. / sl;
sl = 1. / scbd;
z__[i__ - 1] = scbd;
L189:
- i__1 = *n;
+ 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];
+ q_1.v[i__ + j * n - (n+1)] = sl * q_1.v[i__ + j * n - (n+1)];
}
/* L195: */
}
/* FIRST TRANSPOSE V FOR MINFIT: */
L200:
- i__2 = *n;
+ 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];
+ s = q_1.v[i__ + j * n - (n+1)];
+ q_1.v[i__ + j * n - (n+1)] = q_1.v[j + i__ * n - (n+1)];
/* L210: */
- q_1.v[j + i__ * 20 - 21] = s;
+ q_1.v[j + i__ * n - (n+1)] = s;
}
/* L220: */
}
/* APPROXIMATING QUADRATIC FORM WITHOUT SQUARING THE CONDITION */
/* NUMBER..... */
- minfit_(&idim, n, machep, &vsmall, q_1.v, d__);
+ minfit_(n, n, machep, &vsmall, q_1.v, d__, e_minfit);
/* .....UNSCALE THE AXES..... */
if (scbd <= 1.) {
goto L250;
}
- i__2 = *n;
+ i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
s = z__[i__ - 1];
- i__1 = *n;
+ 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];
+ q_1.v[i__ + j * n - (n+1)] = s * q_1.v[i__ + j * n - (n+1)];
}
/* L230: */
}
- i__2 = *n;
+ i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
s = 0.;
- i__1 = *n;
+ i__1 = n;
for (j = 1; j <= i__1; ++j) {
/* L235: */
/* Computing 2nd power */
- d__1 = q_1.v[j + i__ * 20 - 21];
+ d__1 = q_1.v[j + i__ * n - (n+1)];
s += d__1 * d__1;
}
s = sqrt(s);
d__[i__ - 1] = s * d__[i__ - 1];
s = 1 / s;
- i__1 = *n;
+ 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];
+ q_1.v[j + i__ * n - (n+1)] = s * q_1.v[j + i__ * n - (n+1)];
}
/* L245: */
}
L250:
- i__2 = *n;
+ i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
dni = dn * d__[i__ - 1];
if (dni > large) {
/* .....SORT THE EIGENVALUES AND EIGENVECTORS..... */
- sort_(&idim, n, d__, q_1.v);
- global_1.dmin__ = d__[*n - 1];
+ sort_(n, n, d__, q_1.v);
+ global_1.dmin__ = d__[n - 1];
if (global_1.dmin__ < small) {
global_1.dmin__ = small;
}
/* .....RETURN..... */
-L400:
- ret_val = global_1.fx;
- return ret_val;
+done:
+ if (ret != NLOPT_OUT_OF_MEMORY) {
+ *minf = q_1.fbest;
+ memcpy(&x[1], q_1.xbest, n * sizeof(double));
+ }
+ free(work);
+ return ret;
} /* praxis_ */
-static void minfit_(int *m, int *n, double *machep,
- double *tol, double *ab, double *q)
+static void minfit_(int m, int n, double machep,
+ double *tol, double *ab, double *q, double *ework)
{
/* 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__;
+ double *e; /* size n */
+ double c__, f = 0.0, g, h__;
int i__, j, k, l;
double s, x, y, z__;
int l2, ii, kk, kt, ll2, lp1;
double eps, temp;
+
+ e = ework;
/* ...AN IMPROVED VERSION OF MINFIT (SEE GOLUB AND REINSCH, 1969) */
/* RESTRICTED TO M=N,P=0. */
/* ...HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM... */
/* Parameter adjustments */
--q;
- ab_dim1 = *m;
+ ab_dim1 = m;
ab_offset = 1 + ab_dim1;
ab -= ab_offset;
/* Function Body */
- if (*n == 1) {
+ if (n == 1) {
goto L200;
}
- eps = *machep;
+ eps = machep;
g = 0.;
x = 0.;
- i__1 = *n;
+ i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
e[i__ - 1] = g;
s = 0.;
l = i__ + 1;
- i__2 = *n;
+ i__2 = n;
for (j = i__; j <= i__2; ++j) {
/* L1: */
/* Computing 2nd power */
}
h__ = f * g - s;
ab[i__ + i__ * ab_dim1] = f - g;
- if (l > *n) {
+ if (l > n) {
goto L4;
}
- i__2 = *n;
+ i__2 = n;
for (j = l; j <= i__2; ++j) {
f = 0.;
- i__3 = *n;
+ 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;
+ 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) {
+ if (i__ == n) {
goto L6;
}
- i__3 = *n;
+ i__3 = n;
for (j = l; j <= i__3; ++j) {
/* L5: */
s += ab[i__ + j * ab_dim1] * ab[i__ + j * ab_dim1];
if (s < *tol) {
goto L10;
}
- if (i__ == *n) {
+ if (i__ == n) {
goto L16;
}
f = ab[i__ + (i__ + 1) * ab_dim1];
g = -g;
}
h__ = f * g - s;
- if (i__ == *n) {
+ if (i__ == n) {
goto L10;
}
ab[i__ + (i__ + 1) * ab_dim1] = f - g;
- i__3 = *n;
+ i__3 = n;
for (j = l; j <= i__3; ++j) {
/* L7: */
e[j - 1] = ab[i__ + j * ab_dim1] / h__;
}
- i__3 = *n;
+ i__3 = n;
for (j = l; j <= i__3; ++j) {
s = 0.;
- i__2 = *n;
+ 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;
+ i__2 = n;
for (k = l; k <= i__2; ++k) {
/* L9: */
ab[j + k * ab_dim1] += s * e[k - 1];
}
}
/* ...ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS... */
- ab[*n + *n * ab_dim1] = 1.;
- g = e[*n - 1];
- l = *n;
- i__1 = *n;
+ 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;
+ i__ = n - ii + 1;
if (g == 0.) {
goto L23;
}
h__ = ab[i__ + (i__ + 1) * ab_dim1] * g;
- i__2 = *n;
+ 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;
+ i__2 = n;
for (j = l; j <= i__2; ++j) {
s = 0.;
- i__3 = *n;
+ 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;
+ 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;
+ i__3 = n;
for (j = l; j <= i__3; ++j) {
ab[i__ + j * ab_dim1] = 0.;
/* L24: */
/* ...DIAGONALIZATION OF THE BIDIAGONAL FORM... */
/* L100: */
eps *= x;
- i__1 = *n;
+ i__1 = n;
for (kk = 1; kk <= i__1; ++kk) {
- k = *n - kk + 1;
+ k = n - kk + 1;
kt = 0;
L101:
++kt;
g = -x * s + g * c__;
h__ = y * s;
y *= c__;
- i__2 = *n;
+ i__2 = n;
for (j = 1; j <= i__2; ++j) {
x = ab[j + (i__ - 1) * ab_dim1];
z__ = ab[j + i__ * ab_dim1];
goto L150;
}
q[k] = -z__;
- i__3 = *n;
+ i__3 = n;
for (j = 1; j <= i__3; ++j) {
/* L141: */
ab[j + k * ab_dim1] = -ab[j + k * ab_dim1];
ab[ab_dim1 + 1] = 1.;
} /* minfit_ */
-static void min_(int n, int j, int nits, double *
+static nlopt_result 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)
+ x, double *t_old, double machep, double *h__, struct global_s *global_1, struct q_s *q_1)
{
/* System generated locals */
int i__1;
int dz;
double xm, sf1, sx1;
double temp, small;
+ nlopt_result ret = NLOPT_SUCCESS;
/* ...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 */
/* Function Body */
/* Computing 2nd power */
- d__1 = *machep;
+ d__1 = machep;
small = d__1 * d__1;
- m2 = sqrt(*machep);
+ m2 = sqrt(machep);
m4 = sqrt(m2);
sf1 = *f1;
sx1 = *x1;
xm = 0.;
fm = global_1->fx;
f0 = global_1->fx;
- dz = *d2 < *machep;
+ dz = *d2 < machep;
/* ...FIND THE STEP SIZE... */
s = 0.;
i__1 = n;
}
t2 = m4 * sqrt(fabs(global_1->fx) / temp + s * global_1->ldt) + m2 *
global_1->ldt;
- s = m4 * s + *t;
+ s = m4 * s + *t_old;
if (dz && t2 > s) {
t2 = s;
}
temp = -1.;
}
*x1 = temp * t2;
- *f1 = flin_(n, j, x1, f,f_data, &x[1], &global_1->nf, q_1);
+ *f1 = flin_(n, j, x1, f,f_data, &x[1], &global_1->nf, q_1, &ret);
+ if (ret != NLOPT_SUCCESS) return ret;
L3:
if (*f1 > fm) {
goto L4;
if (f0 >= *f1) {
x2 = *x1 * 2.;
}
- f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1);
+ f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1, &ret);
+ if (ret != NLOPT_SUCCESS) return ret;
if (f2 > fm) {
goto L5;
}
x2 = *h__;
/* ...EVALUATE F AT THE PREDICTED MINIMUM... */
L11:
- f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1);
+ f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1, &ret);
+ if (ret != NLOPT_SUCCESS) return ret;
if (k >= nits || f2 <= f0) {
goto L12;
}
/* ...UPDATE X FOR LINEAR BUT NOT PARABOLIC SEARCH... */
L17:
if (j == 0) {
- return;
+ return NLOPT_SUCCESS;
}
i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* L18: */
- x[i__] += *x1 * q_1->v[i__ + j * 20 - 21];
+ x[i__] += *x1 * q_1->v[i__ + j * n - (n+1)];
}
+ return NLOPT_SUCCESS;
} /* 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)
+ int *nf, struct q_s *q_1, nlopt_result *ret)
{
/* System generated locals */
int i__1;
double ret_val;
/* Local variables */
+ nlopt_stopping *stop = q_1->stop;
int i__;
- double t[20];
+ double *t; /* size n */
+
+ t = q_1->t_flin;
/* ...FLIN IS THE FUNCTION OF ONE REAL VARIABLE L THAT IS MINIMIZED */
/* BY THE SUBROUTINE MIN... */
i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* L1: */
- t[i__ - 1] = x[i__] + *l * q_1->v[i__ + j * 20 - 21];
+ t[i__ - 1] = x[i__] + *l * q_1->v[i__ + j * n - (n+1)];
}
goto L4;
/* ...THE SEARCH IS ALONG A PARABOLIC SPACE CURVE... */
L4:
++(*nf);
ret_val = f(n, t, f_data);
+ stop->nevals++;
+ if (ret_val < q_1->fbest) {
+ q_1->fbest = ret_val;
+ memcpy(q_1->xbest, t, n * sizeof(double));
+ }
+ if (nlopt_stop_forced(stop)) *ret = NLOPT_FORCED_STOP;
+ else if (nlopt_stop_evals(stop)) *ret = NLOPT_MAXEVAL_REACHED;
+ else if (nlopt_stop_time(stop)) *ret = NLOPT_MAXTIME_REACHED;
+ else if (ret_val <= stop->minf_max) *ret = NLOPT_MINF_MAX_REACHED;
return ret_val;
} /* flin_ */
-static void sort_(int *m, int *n, double *d__, double *v)
+static void sort_(int m, int n, double *d__, double *v)
{
/* System generated locals */
int v_dim1, v_offset, i__1, i__2;
/* 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_dim1 = m;
v_offset = 1 + v_dim1;
v -= v_offset;
--d__;
/* Function Body */
- if (*n == 1) {
+ if (n == 1) {
return;
}
- nm1 = *n - 1;
+ nm1 = n - 1;
i__1 = nm1;
for (i__ = 1; i__ <= i__1; ++i__) {
k = i__;
s = d__[i__];
ip1 = i__ + 1;
- i__2 = *n;
+ i__2 = n;
for (j = ip1; j <= i__2; ++j) {
if (d__[j] <= s) {
goto L1;
}
d__[k] = d__[i__];
d__[i__] = s;
- i__2 = *n;
+ 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];
}
} /* 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)
+static void quad_(int n, praxis_func f, void *f_data, double *x, double *t_old,
+ double machep, double *h__, struct global_s *global_1, struct q_s *q_1)
{
/* System generated locals */
int i__1;
global_1->fx = q_1->qf1;
q_1->qf1 = s;
q_1->qd1 = 0.;
- i__1 = *n;
+ i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
s = x[i__];
l = q_1->q1[i__ - 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) {
+ 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,
+ min_(n, 0, 2, &s, &l, &value, 1, f,f_data, &x[1], t_old, 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 = 1.;
L3:
q_1->qd0 = q_1->qd1;
- i__1 = *n;
+ i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
s = q_1->q0[i__ - 1];
q_1->q0[i__ - 1] = x[i__];
~ianmdlvl / nlopt.git / blobdiff