1 /* cobyla : contrained optimization by linear approximation */
4 * Copyright (c) 1992, Michael J. D. Powell (M.J.D.Powell@damtp.cam.ac.uk)
5 * Copyright (c) 2004, Jean-Sebastien Roy (js@jeannot.org)
7 * Permission is hereby granted, free of charge, to any person obtaining a
8 * copy of this software and associated documentation files (the
9 * "Software"), to deal in the Software without restriction, including
10 * without limitation the rights to use, copy, modify, merge, publish,
11 * distribute, sublicense, and/or sell copies of the Software, and to
12 * permit persons to whom the Software is furnished to do so, subject to
13 * the following conditions:
15 * The above copyright notice and this permission notice shall be included
16 * in all copies or substantial portions of the Software.
18 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
19 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
20 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
21 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
22 * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
23 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
24 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
28 * This software is a C version of COBYLA2, a contrained optimization by linear
29 * approximation package developed by Michael J. D. Powell in Fortran.
31 * The original source code can be found at :
32 * http://plato.la.asu.edu/topics/problems/nlores.html
35 static char const rcsid[] =
36 "@(#) $Jeannot: cobyla.c,v 1.11 2004/04/18 09:51:36 js Exp $";
44 /* SGJ, 2008: modified COBYLA code to take explicit account of bound
45 constraints. Since bound constraints are linear, these should
46 already be handled exactly when COBYLA optimizes it linear model.
47 However, they were not handled properly when COBYLA created its
48 initial simplex, or when COBYLA updated unacceptable simplices.
49 Since NLopt guarantees that the objective will not be evaluated
50 outside the bound constraints, this required us to handle such
51 points by putting a slope discontinuity into the objective &
52 constraints (below), which slows convergence considerably for
53 smooth functions. Instead, handling them explicitly prevents this
55 #define ENFORCE_BOUNDS 1
57 #define MIN2(a,b) ((a) <= (b) ? (a) : (b))
58 #define MAX2(a,b) ((a) >= (b) ? (a) : (b))
60 #define U(n) ((unsigned) (n))
62 /**************************************************************************/
63 /* SGJ, 2008: NLopt-style interface function: */
74 double *con_tol, *scale;
78 static int func_wrap(int ni, int mi, double *x, double *f, double *con,
83 double *xtmp = s->xtmp;
84 const double *lb = s->lb, *ub = s->ub;
86 (void) mi; /* unused */
88 /* in nlopt, we guarante that the function is never evaluated outside
89 the lb and ub bounds, so we need force this with xtmp ... note
90 that this leads to discontinuity in the first derivative, which
91 slows convergence if we don't enable the ENFORCE_BOUNDS feature
93 for (j = 0; j < n; ++j) {
94 if (x[j] < lb[j]) xtmp[j] = lb[j];
95 else if (x[j] > ub[j]) xtmp[j] = ub[j];
98 nlopt_unscale(n, s->scale, xtmp, xtmp);
100 *f = s->f(n, xtmp, NULL, s->f_data);
101 if (nlopt_stop_forced(s->stop)) return 1;
103 for (j = 0; j < s->m_orig; ++j) {
104 nlopt_eval_constraint(con + i, NULL, s->fc+j, n, xtmp);
105 if (nlopt_stop_forced(s->stop)) return 1;
106 for (k = 0; k < s->fc[j].m; ++k)
107 con[i + k] = -con[i + k];
110 for (j = 0; j < s->p; ++j) {
111 nlopt_eval_constraint(con + i, NULL, s->h+j, n, xtmp);
112 if (nlopt_stop_forced(s->stop)) return 1;
113 for (k = 0; k < s->h[j].m; ++k)
114 con[(i + s->h[j].m) + k] = -con[i + k];
117 for (j = 0; j < n; ++j) {
118 if (!nlopt_isinf(lb[j]))
119 con[i++] = x[j] - lb[j];
120 if (!nlopt_isinf(ub[j]))
121 con[i++] = ub[j] - x[j];
130 COBYLA_MSG_NONE = 0, /* No messages */
131 COBYLA_MSG_EXIT = 1, /* Exit reasons */
132 COBYLA_MSG_ITER = 2, /* Rho and Sigma changes */
133 COBYLA_MSG_INFO = 3 /* Informational messages */
137 * A function as required by cobyla
138 * state is a void pointer provided to the function at each call
140 * n : the number of variables
141 * m : the number of constraints
142 * x : on input, then vector of variables (should not be modified)
143 * f : on output, the value of the function
144 * con : on output, the value of the constraints (vector of size m)
145 * state : on input, the value of the state variable as provided to cobyla
147 * COBYLA will try to make all the values of the constraints positive.
148 * So if you want to input a constraint j such as x[i] <= MAX, set:
149 * con[j] = MAX - x[i]
150 * The function must returns 0 if no error occurs or 1 to immediately end the
154 typedef int cobyla_function(int n, int m, double *x, double *f, double *con,
155 func_wrap_state *state);
158 * cobyla : minimize a function subject to constraints
160 * n : number of variables (>=0)
161 * m : number of constraints (>=0)
162 * x : on input, initial estimate ; on output, the solution
163 * minf : on output, minimum objective function
164 * rhobeg : a reasonable initial change to the variables
165 * stop : the NLopt stopping criteria
166 * lb, ub : lower and upper bounds on x
167 * message : see the cobyla_message enum
168 * calcfc : the function to minimize (see cobyla_function)
169 * state : used by function (see cobyla_function)
171 * The cobyla function returns the usual nlopt_result codes.
174 extern nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, double rhoend, nlopt_stopping *stop, const double *lb, const double *ub,
175 int message, cobyla_function *calcfc, func_wrap_state *state);
177 nlopt_result cobyla_minimize(unsigned n, nlopt_func f, void *f_data,
178 unsigned m, nlopt_constraint *fc,
179 unsigned p, nlopt_constraint *h,
180 const double *lb, const double *ub, /* bounds */
181 double *x, /* in: initial guess, out: minimizer */
183 nlopt_stopping *stop,
189 double rhobeg, rhoend;
191 s.f = f; s.f_data = f_data;
197 s.lb = s.ub = s.xtmp = s.con_tol = s.scale = NULL;
199 s.scale = nlopt_compute_rescaling(n, dx);
200 if (!s.scale) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
202 s.lb = nlopt_new_rescaled(n, s.scale, lb);
203 if (!s.lb) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
204 s.ub = nlopt_new_rescaled(n, s.scale, ub);
205 if (!s.ub) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
206 nlopt_reorder_bounds(n, s.lb, s.ub);
208 s.xtmp = (double *) malloc(sizeof(double) * n);
209 if (!s.xtmp) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
211 /* SGJ, 2008: compute rhoend from NLopt stop info */
212 rhobeg = fabs(dx[0] / s.scale[0]);
213 rhoend = stop->xtol_rel * (rhobeg);
214 for (j = 0; j < n; ++j)
215 if (rhoend < stop->xtol_abs[j] / fabs(s.scale[j]))
216 rhoend = stop->xtol_abs[j] / fabs(s.scale[j]);
218 /* each equality constraint gives two inequality constraints */
219 m = nlopt_count_constraints(m, fc) + 2 * nlopt_count_constraints(p, h);
221 /* add constraints for lower/upper bounds (if any) */
222 for (j = 0; j < n; ++j) {
223 if (!nlopt_isinf(lb[j]))
225 if (!nlopt_isinf(ub[j]))
229 s.con_tol = (double *) malloc(sizeof(double) * m);
230 if (m && !s.con_tol) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
232 for (j = 0; j < m; ++j) s.con_tol[j] = 0;
233 for (j = i = 0; i < s.m_orig; ++i) {
234 unsigned ji = j, jnext = j + fc[i].m;
235 for (; j < jnext; ++j) s.con_tol[j] = fc[i].tol[j - ji];
237 for (i = 0; i < s.p; ++i) {
238 unsigned ji = j, jnext = j + h[i].m;
239 for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - ji];
240 ji = j; jnext = j + h[i].m;
241 for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - ji];
244 nlopt_rescale(n, s.scale, x, x);
245 ret = cobyla((int) n, (int) m, x, minf, rhobeg, rhoend,
246 stop, s.lb, s.ub, COBYLA_MSG_NONE,
248 nlopt_unscale(n, s.scale, x, x);
250 /* make sure e.g. rounding errors didn't push us slightly out of bounds */
251 for (j = 0; j < n; ++j) {
252 if (x[j] < lb[j]) x[j] = lb[j];
253 if (x[j] > ub[j]) x[j] = ub[j];
265 /**************************************************************************/
267 /* SGJ, 2010: I find that it seems to increase robustness of the algorithm
268 if, on "simplex" steps (which are intended to improve the linear
269 independence of the simplex, not improve the objective), we multiply
270 the steps by pseudorandom numbers in [0,1). Intuitively, pseudorandom
271 steps are likely to avoid any accidental dependency in the simplex
272 vertices. However, since I still want COBYLA to be a deterministic
273 algorithm, and I don't care all that much about the quality of the
274 randomness here, I implement a simple linear congruential generator
275 rather than use nlopt_urand, and set the initial seed deterministically. */
277 #if defined(HAVE_STDINT_H)
281 #ifndef HAVE_UINT32_T
282 # if SIZEOF_UNSIGNED_LONG == 4
283 typedef unsigned long uint32_t;
284 # elif SIZEOF_UNSIGNED_INT == 4
285 typedef unsigned int uint32_t;
287 # error No 32-bit unsigned integer type
291 /* a simple linear congruential generator */
293 static uint32_t lcg_rand(uint32_t *seed)
295 return (*seed = *seed * 1103515245 + 12345);
298 static double lcg_urand(uint32_t *seed, double a, double b)
300 return a + lcg_rand(seed) * (b - a) / ((uint32_t) -1);
303 /**************************************************************************/
305 static nlopt_result cobylb(int *n, int *m, int *mpp, double *x, double *minf, double *rhobeg, double rhoend,
306 nlopt_stopping *stop, const double *lb, const double *ub, int *iprint, double *con, double *sim,
307 double *simi, double *datmat, double *a, double *vsig, double *veta,
308 double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
309 func_wrap_state *state);
310 static nlopt_result trstlp(int *n, int *m, double *a, double *b, double *rho,
311 double *dx, int *ifull, int *iact, double *z__, double *zdota, double *vmultc,
312 double *sdirn, double *dxnew, double *vmultd);
314 /* ------------------------------------------------------------------------ */
316 nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, double rhoend, nlopt_stopping *stop, const double *lb, const double *ub, int iprint,
317 cobyla_function *calcfc, func_wrap_state *state)
319 int icon, isim, isigb, idatm, iveta, isimi, ivsig, iwork, ia, idx, mpp;
325 * This subroutine minimizes an objective function F(X) subject to M
326 * inequality constraints on X, where X is a vector of variables that has
327 * N components. The algorithm employs linear approximations to the
328 * objective and constraint functions, the approximations being formed by
329 * linear interpolation at N+1 points in the space of the variables.
330 * We regard these interpolation points as vertices of a simplex. The
331 * parameter RHO controls the size of the simplex and it is reduced
332 * automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
333 * to achieve a good vector of variables for the current size, and then
334 * RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and
335 * RHOEND should be set to reasonable initial changes to and the required
336 * accuracy in the variables respectively, but this accuracy should be
337 * viewed as a subject for experimentation because it is not guaranteed.
338 * The subroutine has an advantage over many of its competitors, however,
339 * which is that it treats each constraint individually when calculating
340 * a change to the variables, instead of lumping the constraints together
341 * into a single penalty function. The name of the subroutine is derived
342 * from the phrase Constrained Optimization BY Linear Approximations.
344 * The user must set the values of N, M, RHOBEG and RHOEND, and must
345 * provide an initial vector of variables in X. Further, the value of
346 * IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
347 * printing during the calculation. Specifically, there is no output if
348 * IPRINT=0 and there is output only at the end of the calculation if
349 * IPRINT=1. Otherwise each new value of RHO and SIGMA is printed.
350 * Further, the vector of variables and some function information are
351 * given either when RHO is reduced or when each new value of F(X) is
352 * computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
353 * is a penalty parameter, it being assumed that a change to X is an
354 * improvement if it reduces the merit function
355 * F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)),
356 * where C1,C2,...,CM denote the constraint functions that should become
357 * nonnegative eventually, at least to the precision of RHOEND. In the
358 * printed output the displayed term that is multiplied by SIGMA is
359 * called MAXCV, which stands for 'MAXimum Constraint Violation'. The
360 * argument MAXFUN is an int variable that must be set by the user to a
361 * limit on the number of calls of CALCFC, the purpose of this routine being
362 * given below. The value of MAXFUN will be altered to the number of calls
363 * of CALCFC that are made. The arguments W and IACT provide real and
364 * int arrays that are used as working space. Their lengths must be at
365 * least N*(3*N+2*M+11)+4*M+6 and M+1 respectively.
367 * In order to define the objective and constraint functions, we require
368 * a subroutine that has the name and arguments
369 * SUBROUTINE CALCFC (N,M,X,F,CON)
370 * DIMENSION X(*),CON(*) .
371 * The values of N and M are fixed and have been defined already, while
372 * X is now the current vector of variables. The subroutine should return
373 * the objective and constraint functions at X in F and CON(1),CON(2),
374 * ...,CON(M). Note that we are trying to adjust X so that F(X) is as
375 * small as possible subject to the constraint functions being nonnegative.
377 * Partition the working space array W to provide the storage that is needed
378 * for the main calculation.
385 if (iprint>=1) fprintf(stderr, "cobyla: N==0.\n");
386 return NLOPT_SUCCESS;
391 if (iprint>=1) fprintf(stderr, "cobyla: N<0 or M<0.\n");
392 return NLOPT_INVALID_ARGS;
395 /* workspace allocation */
396 w = (double*) malloc(U(n*(3*n+2*m+11)+4*m+6)*sizeof(*w));
399 if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
400 return NLOPT_OUT_OF_MEMORY;
402 iact = (int*)malloc(U(m+1)*sizeof(*iact));
405 if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
407 return NLOPT_OUT_OF_MEMORY;
410 /* Parameter adjustments */
420 isimi = isim + n * n + n;
421 idatm = isimi + n * n;
422 ia = idatm + n * mpp + mpp;
423 ivsig = ia + m * n + n;
428 rc = cobylb(&n, &m, &mpp, &x[1], minf, &rhobeg, rhoend, stop, &lb[1], &ub[1], &iprint,
429 &w[icon], &w[isim], &w[isimi], &w[idatm], &w[ia], &w[ivsig], &w[iveta],
430 &w[isigb], &w[idx], &w[iwork], &iact[1], calcfc, state);
432 /* Parameter adjustments (reverse) */
442 /* ------------------------------------------------------------------------- */
443 static nlopt_result cobylb(int *n, int *m, int *mpp,
444 double *x, double *minf, double *rhobeg, double rhoend,
445 nlopt_stopping *stop, const double *lb, const double *ub,
446 int *iprint, double *con, double *sim, double *simi,
447 double *datmat, double *a, double *vsig, double *veta,
448 double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
449 func_wrap_state *state)
451 /* System generated locals */
452 int sim_dim1, sim_offset, simi_dim1, simi_offset, datmat_dim1,
453 datmat_offset, a_dim1, a_offset, i__1, i__2, i__3;
456 /* Local variables */
457 double alpha, delta, denom, tempa, barmu;
458 double beta, cmin = 0.0, cmax = 0.0;
459 double cvmaxm, dxsign, prerem = 0.0;
460 double edgmax, pareta, prerec = 0.0, phimin, parsig = 0.0;
462 double phi, rho, sum = 0.0;
463 double ratio, vmold, parmu, error, vmnew;
464 double resmax, cvmaxp;
465 double resnew, trured;
466 double temp, wsig, f;
475 int mp, np, iz, ibrnch;
476 int nbest, ifull, iptem, jdrop;
477 nlopt_result rc = NLOPT_SUCCESS;
478 uint32_t seed = (uint32_t) (*n + *m); /* arbitrary deterministic LCG seed */
481 /* SGJ, 2008: added code to keep track of minimum feasible function val */
484 /* Set the initial values of some parameters. The last column of SIM holds */
485 /* the optimal vertex of the current simplex, and the preceding N columns */
486 /* hold the displacements from the optimal vertex to the other vertices. */
487 /* Further, SIMI holds the inverse of the matrix that is contained in the */
488 /* first N columns of SIM. */
490 /* Parameter adjustments */
492 a_offset = 1 + a_dim1 * 1;
495 simi_offset = 1 + simi_dim1 * 1;
498 sim_offset = 1 + sim_dim1 * 1;
501 datmat_offset = 1 + datmat_dim1 * 1;
502 datmat -= datmat_offset;
526 "cobyla: the initial value of RHO is %12.6E and PARMU is set to zero.\n",
531 for (i__ = 1; i__ <= i__1; ++i__) {
533 sim[i__ + np * sim_dim1] = x[i__];
535 for (j = 1; j <= i__2; ++j) {
536 sim[i__ + j * sim_dim1] = 0.;
537 simi[i__ + j * simi_dim1] = 0.;
541 /* SGJ: make sure step rhocur stays inside [lb,ub] */
542 if (x[i__] + rhocur > ub[i__]) {
543 if (x[i__] - rhocur >= lb[i__])
545 else if (ub[i__] - x[i__] > x[i__] - lb[i__])
546 rhocur = 0.5 * (ub[i__] - x[i__]);
548 rhocur = 0.5 * (x[i__] - lb[i__]);
551 sim[i__ + i__ * sim_dim1] = rhocur;
552 simi[i__ + i__ * simi_dim1] = 1.0 / rhocur;
557 /* Make the next call of the user-supplied subroutine CALCFC. These */
558 /* instructions are also used for calling CALCFC during the iterations of */
561 /* SGJ comment: why the hell couldn't he just use a damn subroutine?
562 #*&!%*@ Fortran-66 spaghetti code */
565 if (nlopt_stop_forced(stop)) rc = NLOPT_FORCED_STOP;
566 else if (stop->nevals > 0) {
567 if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
568 else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
570 if (rc != NLOPT_SUCCESS) goto L600;
573 if (calcfc(*n, *m, &x[1], &f, &con[1], state))
576 fprintf(stderr, "cobyla: user requested end of minimization.\n");
578 rc = NLOPT_FORCED_STOP;
583 feasible = 1; /* SGJ, 2010 */
586 for (k = 1; k <= i__1; ++k) {
587 d__1 = resmax, d__2 = -con[k];
588 resmax = MAX2(d__1,d__2);
589 if (d__2 > state->con_tol[k-1])
590 feasible = 0; /* SGJ, 2010 */
594 /* SGJ, 2008: check for minf_max reached by feasible point */
595 if (f < stop->minf_max && feasible) {
596 rc = NLOPT_MINF_MAX_REACHED;
597 goto L620; /* not L600 because we want to use current x, f, resmax */
600 if (stop->nevals == *iprint - 1 || *iprint == 3) {
601 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
602 stop->nevals, f, resmax);
604 fprintf(stderr, "cobyla: X =");
605 for (i__ = 1; i__ <= i__1; ++i__) {
606 if (i__>1) fprintf(stderr, " ");
607 fprintf(stderr, "%13.6E", x[i__]);
611 for (i__ = iptemp; i__ <= i__1; ++i__) {
612 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
613 fprintf(stderr, "%15.6E", x[i__]);
616 fprintf(stderr, "\n");
624 /* Set the recently calculated function values in a column of DATMAT. This */
625 /* array has a column for each vertex of the current simplex, the entries of */
626 /* each column being the values of the constraint functions (if any) */
627 /* followed by the objective function and the greatest constraint violation */
631 for (k = 1; k <= i__1; ++k) {
632 datmat[k + jdrop * datmat_dim1] = con[k];
634 if (stop->nevals > np) {
638 /* Exchange the new vertex of the initial simplex with the optimal vertex if */
639 /* necessary. Then, if the initial simplex is not complete, pick its next */
640 /* vertex and calculate the function values there. */
643 if (datmat[mp + np * datmat_dim1] <= f) {
644 x[jdrop] = sim[jdrop + np * sim_dim1];
645 } else { /* improvement in function val */
646 double rhocur = x[jdrop] - sim[jdrop + np * sim_dim1];
647 /* SGJ: use rhocur instead of rho. In original code, rhocur == rho
648 always, but I want to change this to ensure that simplex points
649 fall within [lb,ub]. */
650 sim[jdrop + np * sim_dim1] = x[jdrop];
652 for (k = 1; k <= i__1; ++k) {
653 datmat[k + jdrop * datmat_dim1] = datmat[k + np * datmat_dim1]
655 datmat[k + np * datmat_dim1] = con[k];
658 for (k = 1; k <= i__1; ++k) {
659 sim[jdrop + k * sim_dim1] = -rhocur;
662 for (i__ = k; i__ <= i__2; ++i__) {
663 temp -= simi[i__ + k * simi_dim1];
665 simi[jdrop + k * simi_dim1] = temp;
669 if (stop->nevals <= *n) { /* evaluating initial simplex */
670 jdrop = stop->nevals;
671 /* SGJ: was += rho, but using sim[jdrop,jdrop] enforces consistency
672 if we change the stepsize above to stay in [lb,ub]. */
673 x[jdrop] += sim[jdrop + jdrop * sim_dim1];
679 /* Identify the optimal vertex of the current simplex. */
682 phimin = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
686 for (j = 1; j <= i__1; ++j) {
687 temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
692 } else if (temp == phimin && parmu == 0.) {
693 if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + nbest *
700 /* Switch the best vertex into pole position if it is not there already, */
701 /* and also update SIM, SIMI and DATMAT. */
705 for (i__ = 1; i__ <= i__1; ++i__) {
706 temp = datmat[i__ + np * datmat_dim1];
707 datmat[i__ + np * datmat_dim1] = datmat[i__ + nbest * datmat_dim1]
709 datmat[i__ + nbest * datmat_dim1] = temp;
712 for (i__ = 1; i__ <= i__1; ++i__) {
713 temp = sim[i__ + nbest * sim_dim1];
714 sim[i__ + nbest * sim_dim1] = 0.;
715 sim[i__ + np * sim_dim1] += temp;
718 for (k = 1; k <= i__2; ++k) {
719 sim[i__ + k * sim_dim1] -= temp;
720 tempa -= simi[k + i__ * simi_dim1];
722 simi[nbest + i__ * simi_dim1] = tempa;
726 /* Make an error return if SIGI is a poor approximation to the inverse of */
727 /* the leading N by N submatrix of SIG. */
731 for (i__ = 1; i__ <= i__1; ++i__) {
733 for (j = 1; j <= i__2; ++j) {
739 for (k = 1; k <= i__3; ++k) if (sim[k + j * sim_dim1] != 0) {
740 temp += simi[i__ + k * simi_dim1] * sim[k + j * sim_dim1];
742 d__1 = error, d__2 = fabs(temp);
743 error = MAX2(d__1,d__2);
748 fprintf(stderr, "cobyla: rounding errors are becoming damaging.\n");
750 rc = NLOPT_ROUNDOFF_LIMITED;
754 /* Calculate the coefficients of the linear approximations to the objective */
755 /* and constraint functions, placing minus the objective function gradient */
756 /* after the constraint gradients in the array A. The vector W is used for */
760 for (k = 1; k <= i__2; ++k) {
761 con[k] = -datmat[k + np * datmat_dim1];
763 for (j = 1; j <= i__1; ++j) {
764 w[j] = datmat[k + j * datmat_dim1] + con[k];
767 for (i__ = 1; i__ <= i__1; ++i__) {
770 for (j = 1; j <= i__3; ++j) {
771 temp += w[j] * simi[j + i__ * simi_dim1];
776 a[i__ + k * a_dim1] = temp;
780 /* Calculate the values of sigma and eta, and set IFLAG=0 if the current */
781 /* simplex is not acceptable. */
784 parsig = alpha * rho;
787 for (j = 1; j <= i__1; ++j) {
791 for (i__ = 1; i__ <= i__2; ++i__) {
792 d__1 = simi[j + i__ * simi_dim1];
794 d__1 = sim[i__ + j * sim_dim1];
797 vsig[j] = 1. / sqrt(wsig);
798 veta[j] = sqrt(weta);
799 if (vsig[j] < parsig || veta[j] > pareta) {
804 /* If a new vertex is needed to improve acceptability, then decide which */
805 /* vertex to drop from the simplex. */
807 if (ibrnch == 1 || iflag == 1) {
813 for (j = 1; j <= i__1; ++j) {
814 if (veta[j] > temp) {
821 for (j = 1; j <= i__1; ++j) {
822 if (vsig[j] < temp) {
829 /* Calculate the step to the new vertex and its sign. */
831 temp = gamma_ * rho * vsig[jdrop];
833 for (i__ = 1; i__ <= i__1; ++i__) {
834 dx[i__] = temp * simi[jdrop + i__ * simi_dim1];
839 for (k = 1; k <= i__1; ++k) {
842 for (i__ = 1; i__ <= i__2; ++i__) {
843 sum += a[i__ + k * a_dim1] * dx[i__];
846 temp = datmat[k + np * datmat_dim1];
847 d__1 = cvmaxp, d__2 = -sum - temp;
848 cvmaxp = MAX2(d__1,d__2);
849 d__1 = cvmaxm, d__2 = sum - temp;
850 cvmaxm = MAX2(d__1,d__2);
854 if (parmu * (cvmaxp - cvmaxm) > sum + sum) {
858 /* Update the elements of SIM and SIMI, and set the next X. */
862 for (i__ = 1; i__ <= i__1; ++i__) {
863 /* SGJ, 2010: pseudo-randomize simplex steps (see LCG comments above) */
864 dx[i__] = dxsign * dx[i__] * lcg_urand(&seed, 0.01, 1);
865 /* SGJ: make sure dx step says in [lb,ub] */
868 double xi = sim[i__ + np * sim_dim1];
870 if (xi + dx[i__] > ub[i__])
872 if (xi + dx[i__] < lb[i__]) {
873 if (xi - dx[i__] <= ub[i__])
875 else { /* try again with halved step */
882 sim[i__ + jdrop * sim_dim1] = dx[i__];
883 temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
886 for (i__ = 1; i__ <= i__1; ++i__) {
887 simi[jdrop + i__ * simi_dim1] /= temp;
890 for (j = 1; j <= i__1; ++j) {
894 for (i__ = 1; i__ <= i__2; ++i__) {
895 temp += simi[j + i__ * simi_dim1] * dx[i__];
898 for (i__ = 1; i__ <= i__2; ++i__) {
899 simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
903 x[j] = sim[j + np * sim_dim1] + dx[j];
907 /* Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO. */
911 izdota = iz + *n * *n;
914 idxnew = isdirn + *n;
916 rc = trstlp(n, m, &a[a_offset], &con[1], &rho, &dx[1], &ifull, &iact[1], &w[
917 iz], &w[izdota], &w[ivmc], &w[isdirn], &w[idxnew], &w[ivmd]);
918 if (rc != NLOPT_SUCCESS) goto L600;
920 /* SGJ: since the bound constraints are linear, we should never get
921 a dx that lies outside the [lb,ub] constraints here, but we'll
922 enforce this anyway just to be paranoid */
924 for (i__ = 1; i__ <= i__1; ++i__) {
925 double xi = sim[i__ + np * sim_dim1];
926 if (xi + dx[i__] > ub[i__]) dx[i__] = ub[i__] - xi;
927 if (xi + dx[i__] < lb[i__]) dx[i__] = xi - lb[i__];
933 for (i__ = 1; i__ <= i__1; ++i__) {
937 if (temp < rho * .25 * rho) {
943 /* Predict the change to F and the new maximum constraint violation if the */
944 /* variables are altered from x(0) to x(0)+DX. */
949 for (k = 1; k <= i__1; ++k) {
952 for (i__ = 1; i__ <= i__2; ++i__) {
953 sum -= a[i__ + k * a_dim1] * dx[i__];
956 resnew = MAX2(resnew,sum);
960 /* Increase PARMU if necessary and branch back if this change alters the */
961 /* optimal vertex. Otherwise PREREM and PREREC will be set to the predicted */
962 /* reductions in the merit function and the maximum constraint violation */
966 prerec = datmat[*mpp + np * datmat_dim1] - resnew;
968 barmu = sum / prerec;
970 if (parmu < barmu * 1.5) {
973 fprintf(stderr, "cobyla: increase in PARMU to %12.6E\n", parmu);
975 phi = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
978 for (j = 1; j <= i__1; ++j) {
979 temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
984 if (temp == phi && parmu == 0.f) {
985 if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + np *
992 prerem = parmu * prerec - sum;
994 /* Calculate the constraint and objective functions at x(*). Then find the */
995 /* actual reduction in the merit function. */
998 for (i__ = 1; i__ <= i__1; ++i__) {
999 x[i__] = sim[i__ + np * sim_dim1] + dx[i__];
1004 vmold = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
1006 vmnew = f + parmu * resmax;
1007 trured = vmold - vmnew;
1008 if (parmu == 0. && f == datmat[mp + np * datmat_dim1]) {
1010 trured = datmat[*mpp + np * datmat_dim1] - resmax;
1013 /* Begin the operations that decide whether x(*) should replace one of the */
1014 /* vertices of the current simplex, the change being mandatory if TRURED is */
1015 /* positive. Firstly, JDROP is set to the index of the vertex that is to be */
1019 if (trured <= 0.f) {
1024 for (j = 1; j <= i__1; ++j) {
1027 for (i__ = 1; i__ <= i__2; ++i__) {
1028 temp += simi[j + i__ * simi_dim1] * dx[i__];
1035 sigbar[j] = temp * vsig[j];
1038 /* Calculate the value of ell. */
1040 edgmax = delta * rho;
1043 for (j = 1; j <= i__1; ++j) {
1044 if (sigbar[j] >= parsig || sigbar[j] >= vsig[j]) {
1049 for (i__ = 1; i__ <= i__2; ++i__) {
1050 d__1 = dx[i__] - sim[i__ + j * sim_dim1];
1051 temp += d__1 * d__1;
1055 if (temp > edgmax) {
1068 /* Revise the simplex by updating the elements of SIM, SIMI and DATMAT. */
1072 for (i__ = 1; i__ <= i__1; ++i__) {
1073 sim[i__ + jdrop * sim_dim1] = dx[i__];
1074 temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
1077 for (i__ = 1; i__ <= i__1; ++i__) {
1078 simi[jdrop + i__ * simi_dim1] /= temp;
1081 for (j = 1; j <= i__1; ++j) {
1085 for (i__ = 1; i__ <= i__2; ++i__) {
1086 temp += simi[j + i__ * simi_dim1] * dx[i__];
1089 for (i__ = 1; i__ <= i__2; ++i__) {
1090 simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
1096 for (k = 1; k <= i__1; ++k) {
1097 datmat[k + jdrop * datmat_dim1] = con[k];
1100 /* Branch back for further iterations with the current RHO. */
1102 if (trured > 0. && trured >= prerem * .1) {
1103 /* SGJ, 2010: following a suggestion in the SAS manual (which
1104 mentions a similar modification to COBYLA, although they didn't
1105 publish their source code), increase rho if predicted reduction
1106 is sufficiently close to the actual reduction. Otherwise,
1107 COBLYA seems to easily get stuck making very small steps.
1108 Also require iflag != 0 (i.e., acceptable simplex), again
1109 following SAS suggestion (otherwise I get convergence failure
1111 if (trured >= prerem * 0.9 && trured <= prerem * 1.1 && iflag) {
1122 /* SGJ, 2008: convergence tests for function vals; note that current
1123 best val is stored in datmat[mp + np * datmat_dim1], or in f if
1124 ifull == 1, and previous best is in *minf. This seems like a
1125 sensible place to put the convergence test, as it is the same
1126 place that Powell checks the x tolerance (rho > rhoend). */
1128 double fbest = ifull == 1 ? f : datmat[mp + np * datmat_dim1];
1129 if (fbest < *minf && nlopt_stop_ftol(stop, fbest, *minf)) {
1130 rc = NLOPT_FTOL_REACHED;
1136 /* Otherwise reduce RHO if it is not at its least value and reset PARMU. */
1140 if (rho <= rhoend * 1.5) {
1146 for (k = 1; k <= i__1; ++k) {
1147 cmin = datmat[k + np * datmat_dim1];
1150 for (i__ = 1; i__ <= i__2; ++i__) {
1151 d__1 = cmin, d__2 = datmat[k + i__ * datmat_dim1];
1152 cmin = MIN2(d__1,d__2);
1153 d__1 = cmax, d__2 = datmat[k + i__ * datmat_dim1];
1154 cmax = MAX2(d__1,d__2);
1156 if (k <= *m && cmin < cmax * .5) {
1157 temp = MAX2(cmax,0.) - cmin;
1161 denom = MIN2(denom,temp);
1167 } else if (cmax - cmin < parmu * denom) {
1168 parmu = (cmax - cmin) / denom;
1172 fprintf(stderr, "cobyla: reduction in RHO to %12.6E and PARMU =%13.6E\n",
1176 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1177 stop->nevals, datmat[mp + np * datmat_dim1], datmat[*mpp + np * datmat_dim1]);
1179 fprintf(stderr, "cobyla: X =");
1181 for (i__ = 1; i__ <= i__1; ++i__) {
1182 if (i__>1) fprintf(stderr, " ");
1183 fprintf(stderr, "%13.6E", sim[i__ + np * sim_dim1]);
1187 for (i__ = iptemp; i__ <= i__1; ++i__) {
1188 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
1189 fprintf(stderr, "%15.6E", x[i__]);
1192 fprintf(stderr, "\n");
1196 else /* rho <= rhoend */
1197 rc = rhoend > 0 ? NLOPT_XTOL_REACHED : NLOPT_ROUNDOFF_LIMITED;
1199 /* Return the best calculated values of the variables. */
1202 fprintf(stderr, "cobyla: normal return.\n");
1209 for (i__ = 1; i__ <= i__1; ++i__) {
1210 x[i__] = sim[i__ + np * sim_dim1];
1212 f = datmat[mp + np * datmat_dim1];
1213 resmax = datmat[*mpp + np * datmat_dim1];
1217 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1218 stop->nevals, f, resmax);
1220 fprintf(stderr, "cobyla: X =");
1221 for (i__ = 1; i__ <= i__1; ++i__) {
1222 if (i__>1) fprintf(stderr, " ");
1223 fprintf(stderr, "%13.6E", x[i__]);
1227 for (i__ = iptemp; i__ <= i__1; ++i__) {
1228 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
1229 fprintf(stderr, "%15.6E", x[i__]);
1232 fprintf(stderr, "\n");
1237 /* ------------------------------------------------------------------------- */
1238 static nlopt_result trstlp(int *n, int *m, double *a,
1239 double *b, double *rho, double *dx, int *ifull,
1240 int *iact, double *z__, double *zdota, double *vmultc,
1241 double *sdirn, double *dxnew, double *vmultd)
1243 /* System generated locals */
1244 int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
1247 /* Local variables */
1248 double alpha, tempa;
1250 double optnew, stpful, sum, tot, acca, accb;
1251 double ratio, vsave, zdotv, zdotw, dd;
1253 double sp, ss, resold = 0.0, zdvabs, zdwabs, sumabs, resmax, optold;
1257 int iout, i__, j, k;
1261 int nact, icon = 0, mcon;
1265 /* This subroutine calculates an N-component vector DX by applying the */
1266 /* following two stages. In the first stage, DX is set to the shortest */
1267 /* vector that minimizes the greatest violation of the constraints */
1268 /* A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K=2,3,...,M, */
1269 /* subject to the Euclidean length of DX being at most RHO. If its length is */
1270 /* strictly less than RHO, then we use the resultant freedom in DX to */
1271 /* minimize the objective function */
1272 /* -A(1,M+1)*DX(1)-A(2,M+1)*DX(2)-...-A(N,M+1)*DX(N) */
1273 /* subject to no increase in any greatest constraint violation. This */
1274 /* notation allows the gradient of the objective function to be regarded as */
1275 /* the gradient of a constraint. Therefore the two stages are distinguished */
1276 /* by MCON .EQ. M and MCON .GT. M respectively. It is possible that a */
1277 /* degeneracy may prevent DX from attaining the target length RHO. Then the */
1278 /* value IFULL=0 would be set, but usually IFULL=1 on return. */
1280 /* In general NACT is the number of constraints in the active set and */
1281 /* IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT */
1282 /* contains a permutation of the remaining constraint indices. Further, Z is */
1283 /* an orthogonal matrix whose first NACT columns can be regarded as the */
1284 /* result of Gram-Schmidt applied to the active constraint gradients. For */
1285 /* J=1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th */
1286 /* column of Z with the gradient of the J-th active constraint. DX is the */
1287 /* current vector of variables and here the residuals of the active */
1288 /* constraints should be zero. Further, the active constraints have */
1289 /* nonnegative Lagrange multipliers that are held at the beginning of */
1290 /* VMULTC. The remainder of this vector holds the residuals of the inactive */
1291 /* constraints at DX, the ordering of the components of VMULTC being in */
1292 /* agreement with the permutation of the indices of the constraints that is */
1293 /* in IACT. All these residuals are nonnegative, which is achieved by the */
1294 /* shift RESMAX that makes the least residual zero. */
1296 /* Initialize Z and some other variables. The value of RESMAX will be */
1297 /* appropriate to DX=0, while ICON will be the index of a most violated */
1298 /* constraint if RESMAX is positive. Usually during the first stage the */
1299 /* vector SDIRN gives a search direction that reduces all the active */
1300 /* constraint violations by one simultaneously. */
1302 /* Parameter adjustments */
1304 z_offset = 1 + z_dim1 * 1;
1307 a_offset = 1 + a_dim1 * 1;
1324 for (i__ = 1; i__ <= i__1; ++i__) {
1326 for (j = 1; j <= i__2; ++j) {
1327 z__[i__ + j * z_dim1] = 0.;
1329 z__[i__ + i__ * z_dim1] = 1.;
1334 for (k = 1; k <= i__1; ++k) {
1335 if (b[k] > resmax) {
1341 for (k = 1; k <= i__1; ++k) {
1343 vmultc[k] = resmax - b[k];
1350 for (i__ = 1; i__ <= i__1; ++i__) {
1354 /* End the current stage of the calculation if 3 consecutive iterations */
1355 /* have either failed to reduce the best calculated value of the objective */
1356 /* function or to increase the number of active constraints since the best */
1357 /* value was calculated. This strategy prevents cycling, but there is a */
1358 /* remote possibility that it will cause premature termination. */
1369 for (i__ = 1; i__ <= i__1; ++i__) {
1370 optnew -= dx[i__] * a[i__ + mcon * a_dim1];
1373 if (icount == 0 || optnew < optold) {
1377 } else if (nact > nactx) {
1387 /* If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to */
1388 /* the active set. Apply Givens rotations so that the last N-NACT-1 columns */
1389 /* of Z are orthogonal to the gradient of the new constraint, a scalar */
1390 /* product being set to zero if its nonzero value could be due to computer */
1391 /* rounding errors. The array DXNEW is used for working space. */
1398 for (i__ = 1; i__ <= i__1; ++i__) {
1399 dxnew[i__] = a[i__ + kk * a_dim1];
1408 for (i__ = 1; i__ <= i__1; ++i__) {
1409 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1411 spabs += fabs(temp);
1413 acca = spabs + fabs(sp) * .1;
1414 accb = spabs + fabs(sp) * .2;
1415 if (spabs >= acca || acca >= accb) {
1422 temp = sqrt(sp * sp + tot * tot);
1427 for (i__ = 1; i__ <= i__1; ++i__) {
1428 temp = alpha * z__[i__ + k * z_dim1] + beta * z__[i__ + kp *
1430 z__[i__ + kp * z_dim1] = alpha * z__[i__ + kp * z_dim1] -
1431 beta * z__[i__ + k * z_dim1];
1432 z__[i__ + k * z_dim1] = temp;
1439 /* Add the new constraint if this can be done without a deletion from the */
1445 vmultc[icon] = vmultc[nact];
1450 /* The next instruction is reached if a deletion has to be made from the */
1451 /* active set in order to make room for the new active constraint, because */
1452 /* the new constraint gradient is a linear combination of the gradients of */
1453 /* the old active constraints. Set the elements of VMULTD to the multipliers */
1454 /* of the linear combination. Further, set IOUT to the index of the */
1455 /* constraint to be deleted, but branch if no suitable index can be found. */
1463 for (i__ = 1; i__ <= i__1; ++i__) {
1464 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1466 zdvabs += fabs(temp);
1468 acca = zdvabs + fabs(zdotv) * .1;
1469 accb = zdvabs + fabs(zdotv) * .2;
1470 if (zdvabs < acca && acca < accb) {
1471 temp = zdotv / zdota[k];
1472 if (temp > 0. && iact[k] <= *m) {
1473 tempa = vmultc[k] / temp;
1474 if (ratio < 0. || tempa < ratio) {
1482 for (i__ = 1; i__ <= i__1; ++i__) {
1483 dxnew[i__] -= temp * a[i__ + kw * a_dim1];
1498 /* Revise the Lagrange multipliers and reorder the active constraints so */
1499 /* that the one to be replaced is at the end of the list. Also calculate the */
1500 /* new value of ZDOTA(NACT) and branch if it is not acceptable. */
1503 for (k = 1; k <= i__1; ++k) {
1504 d__1 = 0., d__2 = vmultc[k] - ratio * vmultd[k];
1505 vmultc[k] = MAX2(d__1,d__2);
1509 vsave = vmultc[icon];
1516 for (i__ = 1; i__ <= i__1; ++i__) {
1517 sp += z__[i__ + k * z_dim1] * a[i__ + kw * a_dim1];
1520 temp = sqrt(sp * sp + d__1 * d__1);
1521 alpha = zdota[kp] / temp;
1523 zdota[kp] = alpha * zdota[k];
1526 for (i__ = 1; i__ <= i__1; ++i__) {
1527 temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1529 z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1530 z__[i__ + kp * z_dim1];
1531 z__[i__ + k * z_dim1] = temp;
1534 vmultc[k] = vmultc[kp];
1544 for (i__ = 1; i__ <= i__1; ++i__) {
1545 temp += z__[i__ + nact * z_dim1] * a[i__ + kk * a_dim1];
1552 vmultc[nact] = ratio;
1554 /* Update IACT and ensure that the objective function continues to be */
1555 /* treated as the last active constraint when MCON>M. */
1558 iact[icon] = iact[nact];
1560 if (mcon > *m && kk != mcon) {
1564 for (i__ = 1; i__ <= i__1; ++i__) {
1565 sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1568 temp = sqrt(sp * sp + d__1 * d__1);
1569 alpha = zdota[nact] / temp;
1571 zdota[nact] = alpha * zdota[k];
1574 for (i__ = 1; i__ <= i__1; ++i__) {
1575 temp = alpha * z__[i__ + nact * z_dim1] + beta * z__[i__ + k *
1577 z__[i__ + nact * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1578 z__[i__ + nact * z_dim1];
1579 z__[i__ + k * z_dim1] = temp;
1581 iact[nact] = iact[k];
1584 vmultc[k] = vmultc[nact];
1585 vmultc[nact] = temp;
1588 /* If stage one is in progress, then set SDIRN to the direction of the next */
1589 /* change to the current vector of variables. */
1597 for (i__ = 1; i__ <= i__1; ++i__) {
1598 temp += sdirn[i__] * a[i__ + kk * a_dim1];
1601 temp /= zdota[nact];
1603 for (i__ = 1; i__ <= i__1; ++i__) {
1604 sdirn[i__] -= temp * z__[i__ + nact * z_dim1];
1608 /* Delete the constraint that has the index IACT(ICON) from the active set. */
1613 vsave = vmultc[icon];
1620 for (i__ = 1; i__ <= i__1; ++i__) {
1621 sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1624 temp = sqrt(sp * sp + d__1 * d__1);
1625 alpha = zdota[kp] / temp;
1627 zdota[kp] = alpha * zdota[k];
1630 for (i__ = 1; i__ <= i__1; ++i__) {
1631 temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1633 z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1634 z__[i__ + kp * z_dim1];
1635 z__[i__ + k * z_dim1] = temp;
1638 vmultc[k] = vmultc[kp];
1648 /* If stage one is in progress, then set SDIRN to the direction of the next */
1649 /* change to the current vector of variables. */
1656 for (i__ = 1; i__ <= i__1; ++i__) {
1657 temp += sdirn[i__] * z__[i__ + (nact + 1) * z_dim1];
1660 for (i__ = 1; i__ <= i__1; ++i__) {
1661 sdirn[i__] -= temp * z__[i__ + (nact + 1) * z_dim1];
1665 /* Pick the next search direction of stage two. */
1668 temp = 1. / zdota[nact];
1670 for (i__ = 1; i__ <= i__1; ++i__) {
1671 sdirn[i__] = temp * z__[i__ + nact * z_dim1];
1674 /* Calculate the step to the boundary of the trust region or take the step */
1675 /* that reduces RESMAX to zero. The two statements below that include the */
1676 /* factor 1.0E-6 prevent some harmless underflows that occurred in a test */
1677 /* calculation. Further, we skip the step if it could be zero within a */
1678 /* reasonable tolerance for computer rounding errors. */
1685 for (i__ = 1; i__ <= i__1; ++i__) {
1686 if ((d__1 = dx[i__], fabs(d__1)) >= *rho * 1e-6f) {
1690 sd += dx[i__] * sdirn[i__];
1697 temp = sqrt(ss * dd);
1698 if (fabs(sd) >= temp * 1e-6f) {
1699 temp = sqrt(ss * dd + sd * sd);
1701 stpful = dd / (temp + sd);
1704 acca = step + resmax * .1;
1705 accb = step + resmax * .2;
1706 if (step >= acca || acca >= accb) {
1709 step = MIN2(step,resmax);
1712 /* SGJ, 2010: check for error here */
1713 if (nlopt_isinf(step)) return NLOPT_ROUNDOFF_LIMITED;
1715 /* Set DXNEW to the new variables if STEP is the steplength, and reduce */
1716 /* RESMAX to the corresponding maximum residual if stage one is being done. */
1717 /* Because DXNEW will be changed during the calculation of some Lagrange */
1718 /* multipliers, it will be restored to the following value later. */
1721 for (i__ = 1; i__ <= i__1; ++i__) {
1722 dxnew[i__] = dx[i__] + step * sdirn[i__];
1728 for (k = 1; k <= i__1; ++k) {
1732 for (i__ = 1; i__ <= i__2; ++i__) {
1733 temp -= a[i__ + kk * a_dim1] * dxnew[i__];
1735 resmax = MAX2(resmax,temp);
1739 /* Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A */
1740 /* device is included to force VMULTD(K)=0.0 if deviations from this value */
1741 /* can be attributed to computer rounding errors. First calculate the new */
1742 /* Lagrange multipliers. */
1749 for (i__ = 1; i__ <= i__1; ++i__) {
1750 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1752 zdwabs += fabs(temp);
1754 acca = zdwabs + fabs(zdotw) * .1;
1755 accb = zdwabs + fabs(zdotw) * .2;
1756 if (zdwabs >= acca || acca >= accb) {
1759 vmultd[k] = zdotw / zdota[k];
1763 for (i__ = 1; i__ <= i__1; ++i__) {
1764 dxnew[i__] -= vmultd[k] * a[i__ + kk * a_dim1];
1770 d__1 = 0., d__2 = vmultd[nact];
1771 vmultd[nact] = MAX2(d__1,d__2);
1774 /* Complete VMULTC by finding the new constraint residuals. */
1777 for (i__ = 1; i__ <= i__1; ++i__) {
1778 dxnew[i__] = dx[i__] + step * sdirn[i__];
1783 for (k = kl; k <= i__1; ++k) {
1785 sum = resmax - b[kk];
1786 sumabs = resmax + (d__1 = b[kk], fabs(d__1));
1788 for (i__ = 1; i__ <= i__2; ++i__) {
1789 temp = a[i__ + kk * a_dim1] * dxnew[i__];
1791 sumabs += fabs(temp);
1793 acca = sumabs + fabs(sum) * .1f;
1794 accb = sumabs + fabs(sum) * .2f;
1795 if (sumabs >= acca || acca >= accb) {
1802 /* Calculate the fraction of the step from DX to DXNEW that will be taken. */
1807 for (k = 1; k <= i__1; ++k) {
1808 if (vmultd[k] < 0.) {
1809 temp = vmultc[k] / (vmultc[k] - vmultd[k]);
1817 /* Update DX, VMULTC and RESMAX. */
1821 for (i__ = 1; i__ <= i__1; ++i__) {
1822 dx[i__] = temp * dx[i__] + ratio * dxnew[i__];
1825 for (k = 1; k <= i__1; ++k) {
1826 d__1 = 0., d__2 = temp * vmultc[k] + ratio * vmultd[k];
1827 vmultc[k] = MAX2(d__1,d__2);
1830 resmax = resold + ratio * (resmax - resold);
1833 /* If the full step is not acceptable then begin another iteration. */
1834 /* Otherwise switch to stage two or end the calculation. */
1839 if (step == stpful) {
1849 /* We employ any freedom that may be available to reduce the objective */
1850 /* function before returning a DX whose length is less than RHO. */
1858 return NLOPT_SUCCESS;