/*************************************************************************/
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)
} nlopt_data;
#include "subplex.h"
+#include "praxis.h"
static double f_subplex(int n, const double *x, 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,
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:
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;
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) {
--- /dev/null
+/* See README */
+
+#include <math.h>
+#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_ */
~ianmdlvl / nlopt.git / commitdiff