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; }
207 s.xtmp = (double *) malloc(sizeof(double) * n);
208 if (!s.xtmp) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
210 /* SGJ, 2008: compute rhoend from NLopt stop info */
211 rhobeg = dx[0] / s.scale[0];
212 rhoend = stop->xtol_rel * (rhobeg);
213 for (j = 0; j < n; ++j)
214 if (rhoend < stop->xtol_abs[j] / s.scale[j])
215 rhoend = stop->xtol_abs[j] / s.scale[j];
217 /* each equality constraint gives two inequality constraints */
218 m = nlopt_count_constraints(m, fc) + 2 * nlopt_count_constraints(p, h);
220 /* add constraints for lower/upper bounds (if any) */
221 for (j = 0; j < n; ++j) {
222 if (!nlopt_isinf(lb[j]))
224 if (!nlopt_isinf(ub[j]))
228 s.con_tol = (double *) malloc(sizeof(double) * m);
229 if (m && !s.con_tol) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
231 for (j = 0; j < m; ++j) s.con_tol[j] = 0;
232 for (j = i = 0; i < s.m_orig; ++i) {
233 unsigned ji = j, jnext = j + fc[i].m;
234 for (; j < jnext; ++j) s.con_tol[j] = fc[i].tol[j - ji];
236 for (i = 0; i < s.p; ++i) {
237 unsigned ji = j, jnext = j + h[i].m;
238 for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - ji];
239 ji = j; jnext = j + h[i].m;
240 for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - ji];
243 nlopt_rescale(n, s.scale, x, x);
244 ret = cobyla((int) n, (int) m, x, minf, rhobeg, rhoend,
245 stop, s.lb, s.ub, COBYLA_MSG_NONE,
247 nlopt_unscale(n, s.scale, x, x);
249 /* make sure e.g. rounding errors didn't push us slightly out of bounds */
250 for (j = 0; j < n; ++j) {
251 if (x[j] < lb[j]) x[j] = lb[j];
252 if (x[j] > ub[j]) x[j] = ub[j];
264 /**************************************************************************/
266 /* SGJ, 2010: I find that it seems to increase robustness of the algorithm
267 if, on "simplex" steps (which are intended to improve the linear
268 independence of the simplex, not improve the objective), we multiply
269 the steps by pseudorandom numbers in [0,1). Intuitively, pseudorandom
270 steps are likely to avoid any accidental dependency in the simplex
271 vertices. However, since I still want COBYLA to be a deterministic
272 algorithm, and I don't care all that much about the quality of the
273 randomness here, I implement a simple linear congruential generator
274 rather than use nlopt_urand, and set the initial seed deterministically. */
276 #if defined(HAVE_STDINT_H)
280 #ifndef HAVE_UINT32_T
281 # if SIZEOF_UNSIGNED_LONG == 4
282 typedef unsigned long uint32_t;
283 # elif SIZEOF_UNSIGNED_INT == 4
284 typedef unsigned int uint32_t;
286 # error No 32-bit unsigned integer type
290 /* a simple linear congruential generator */
292 static uint32_t lcg_rand(uint32_t *seed)
294 return (*seed = *seed * 1103515245 + 12345);
297 static double lcg_urand(uint32_t *seed, double a, double b)
299 return a + lcg_rand(seed) * (b - a) / ((uint32_t) -1);
302 /**************************************************************************/
304 static nlopt_result cobylb(int *n, int *m, int *mpp, double *x, double *minf, double *rhobeg, double rhoend,
305 nlopt_stopping *stop, const double *lb, const double *ub, int *iprint, double *con, double *sim,
306 double *simi, double *datmat, double *a, double *vsig, double *veta,
307 double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
308 func_wrap_state *state);
309 static nlopt_result trstlp(int *n, int *m, double *a, double *b, double *rho,
310 double *dx, int *ifull, int *iact, double *z__, double *zdota, double *vmultc,
311 double *sdirn, double *dxnew, double *vmultd);
313 /* ------------------------------------------------------------------------ */
315 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,
316 cobyla_function *calcfc, func_wrap_state *state)
318 int icon, isim, isigb, idatm, iveta, isimi, ivsig, iwork, ia, idx, mpp;
324 * This subroutine minimizes an objective function F(X) subject to M
325 * inequality constraints on X, where X is a vector of variables that has
326 * N components. The algorithm employs linear approximations to the
327 * objective and constraint functions, the approximations being formed by
328 * linear interpolation at N+1 points in the space of the variables.
329 * We regard these interpolation points as vertices of a simplex. The
330 * parameter RHO controls the size of the simplex and it is reduced
331 * automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
332 * to achieve a good vector of variables for the current size, and then
333 * RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and
334 * RHOEND should be set to reasonable initial changes to and the required
335 * accuracy in the variables respectively, but this accuracy should be
336 * viewed as a subject for experimentation because it is not guaranteed.
337 * The subroutine has an advantage over many of its competitors, however,
338 * which is that it treats each constraint individually when calculating
339 * a change to the variables, instead of lumping the constraints together
340 * into a single penalty function. The name of the subroutine is derived
341 * from the phrase Constrained Optimization BY Linear Approximations.
343 * The user must set the values of N, M, RHOBEG and RHOEND, and must
344 * provide an initial vector of variables in X. Further, the value of
345 * IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
346 * printing during the calculation. Specifically, there is no output if
347 * IPRINT=0 and there is output only at the end of the calculation if
348 * IPRINT=1. Otherwise each new value of RHO and SIGMA is printed.
349 * Further, the vector of variables and some function information are
350 * given either when RHO is reduced or when each new value of F(X) is
351 * computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
352 * is a penalty parameter, it being assumed that a change to X is an
353 * improvement if it reduces the merit function
354 * F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)),
355 * where C1,C2,...,CM denote the constraint functions that should become
356 * nonnegative eventually, at least to the precision of RHOEND. In the
357 * printed output the displayed term that is multiplied by SIGMA is
358 * called MAXCV, which stands for 'MAXimum Constraint Violation'. The
359 * argument MAXFUN is an int variable that must be set by the user to a
360 * limit on the number of calls of CALCFC, the purpose of this routine being
361 * given below. The value of MAXFUN will be altered to the number of calls
362 * of CALCFC that are made. The arguments W and IACT provide real and
363 * int arrays that are used as working space. Their lengths must be at
364 * least N*(3*N+2*M+11)+4*M+6 and M+1 respectively.
366 * In order to define the objective and constraint functions, we require
367 * a subroutine that has the name and arguments
368 * SUBROUTINE CALCFC (N,M,X,F,CON)
369 * DIMENSION X(*),CON(*) .
370 * The values of N and M are fixed and have been defined already, while
371 * X is now the current vector of variables. The subroutine should return
372 * the objective and constraint functions at X in F and CON(1),CON(2),
373 * ...,CON(M). Note that we are trying to adjust X so that F(X) is as
374 * small as possible subject to the constraint functions being nonnegative.
376 * Partition the working space array W to provide the storage that is needed
377 * for the main calculation.
384 if (iprint>=1) fprintf(stderr, "cobyla: N==0.\n");
385 return NLOPT_SUCCESS;
390 if (iprint>=1) fprintf(stderr, "cobyla: N<0 or M<0.\n");
391 return NLOPT_INVALID_ARGS;
394 /* workspace allocation */
395 w = (double*) malloc(U(n*(3*n+2*m+11)+4*m+6)*sizeof(*w));
398 if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
399 return NLOPT_OUT_OF_MEMORY;
401 iact = (int*)malloc(U(m+1)*sizeof(*iact));
404 if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
406 return NLOPT_OUT_OF_MEMORY;
409 /* Parameter adjustments */
419 isimi = isim + n * n + n;
420 idatm = isimi + n * n;
421 ia = idatm + n * mpp + mpp;
422 ivsig = ia + m * n + n;
427 rc = cobylb(&n, &m, &mpp, &x[1], minf, &rhobeg, rhoend, stop, &lb[1], &ub[1], &iprint,
428 &w[icon], &w[isim], &w[isimi], &w[idatm], &w[ia], &w[ivsig], &w[iveta],
429 &w[isigb], &w[idx], &w[iwork], &iact[1], calcfc, state);
431 /* Parameter adjustments (reverse) */
441 /* ------------------------------------------------------------------------- */
442 static nlopt_result cobylb(int *n, int *m, int *mpp,
443 double *x, double *minf, double *rhobeg, double rhoend,
444 nlopt_stopping *stop, const double *lb, const double *ub,
445 int *iprint, double *con, double *sim, double *simi,
446 double *datmat, double *a, double *vsig, double *veta,
447 double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
448 func_wrap_state *state)
450 /* System generated locals */
451 int sim_dim1, sim_offset, simi_dim1, simi_offset, datmat_dim1,
452 datmat_offset, a_dim1, a_offset, i__1, i__2, i__3;
455 /* Local variables */
456 double alpha, delta, denom, tempa, barmu;
457 double beta, cmin = 0.0, cmax = 0.0;
458 double cvmaxm, dxsign, prerem = 0.0;
459 double edgmax, pareta, prerec = 0.0, phimin, parsig = 0.0;
461 double phi, rho, sum = 0.0;
462 double ratio, vmold, parmu, error, vmnew;
463 double resmax, cvmaxp;
464 double resnew, trured;
465 double temp, wsig, f;
474 int mp, np, iz, ibrnch;
475 int nbest, ifull, iptem, jdrop;
476 nlopt_result rc = NLOPT_SUCCESS;
477 uint32_t seed = (uint32_t) (*n + *m); /* arbitrary deterministic LCG seed */
480 /* SGJ, 2008: added code to keep track of minimum feasible function val */
483 /* Set the initial values of some parameters. The last column of SIM holds */
484 /* the optimal vertex of the current simplex, and the preceding N columns */
485 /* hold the displacements from the optimal vertex to the other vertices. */
486 /* Further, SIMI holds the inverse of the matrix that is contained in the */
487 /* first N columns of SIM. */
489 /* Parameter adjustments */
491 a_offset = 1 + a_dim1 * 1;
494 simi_offset = 1 + simi_dim1 * 1;
497 sim_offset = 1 + sim_dim1 * 1;
500 datmat_offset = 1 + datmat_dim1 * 1;
501 datmat -= datmat_offset;
525 "cobyla: the initial value of RHO is %12.6E and PARMU is set to zero.\n",
530 for (i__ = 1; i__ <= i__1; ++i__) {
532 sim[i__ + np * sim_dim1] = x[i__];
534 for (j = 1; j <= i__2; ++j) {
535 sim[i__ + j * sim_dim1] = 0.;
536 simi[i__ + j * simi_dim1] = 0.;
540 /* SGJ: make sure step rhocur stays inside [lb,ub] */
541 if (x[i__] + rhocur > ub[i__]) {
542 if (x[i__] - rhocur >= lb[i__])
544 else if (ub[i__] - x[i__] > x[i__] - lb[i__])
545 rhocur = 0.5 * (ub[i__] - x[i__]);
547 rhocur = 0.5 * (x[i__] - lb[i__]);
550 sim[i__ + i__ * sim_dim1] = rhocur;
551 simi[i__ + i__ * simi_dim1] = 1.0 / rhocur;
556 /* Make the next call of the user-supplied subroutine CALCFC. These */
557 /* instructions are also used for calling CALCFC during the iterations of */
560 /* SGJ comment: why the hell couldn't he just use a damn subroutine?
561 #*&!%*@ Fortran-66 spaghetti code */
564 if (nlopt_stop_forced(stop)) rc = NLOPT_FORCED_STOP;
565 else if (stop->nevals > 0) {
566 if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
567 else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
569 if (rc != NLOPT_SUCCESS) goto L600;
572 if (calcfc(*n, *m, &x[1], &f, &con[1], state))
575 fprintf(stderr, "cobyla: user requested end of minimization.\n");
577 rc = NLOPT_FORCED_STOP;
582 feasible = 1; /* SGJ, 2010 */
585 for (k = 1; k <= i__1; ++k) {
586 d__1 = resmax, d__2 = -con[k];
587 resmax = MAX2(d__1,d__2);
588 if (d__2 > state->con_tol[k-1])
589 feasible = 0; /* SGJ, 2010 */
593 /* SGJ, 2008: check for minf_max reached by feasible point */
594 if (f < stop->minf_max && feasible) {
595 rc = NLOPT_MINF_MAX_REACHED;
596 goto L620; /* not L600 because we want to use current x, f, resmax */
599 if (stop->nevals == *iprint - 1 || *iprint == 3) {
600 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
601 stop->nevals, f, resmax);
603 fprintf(stderr, "cobyla: X =");
604 for (i__ = 1; i__ <= i__1; ++i__) {
605 if (i__>1) fprintf(stderr, " ");
606 fprintf(stderr, "%13.6E", x[i__]);
610 for (i__ = iptemp; i__ <= i__1; ++i__) {
611 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
612 fprintf(stderr, "%15.6E", x[i__]);
615 fprintf(stderr, "\n");
623 /* Set the recently calculated function values in a column of DATMAT. This */
624 /* array has a column for each vertex of the current simplex, the entries of */
625 /* each column being the values of the constraint functions (if any) */
626 /* followed by the objective function and the greatest constraint violation */
630 for (k = 1; k <= i__1; ++k) {
631 datmat[k + jdrop * datmat_dim1] = con[k];
633 if (stop->nevals > np) {
637 /* Exchange the new vertex of the initial simplex with the optimal vertex if */
638 /* necessary. Then, if the initial simplex is not complete, pick its next */
639 /* vertex and calculate the function values there. */
642 if (datmat[mp + np * datmat_dim1] <= f) {
643 x[jdrop] = sim[jdrop + np * sim_dim1];
644 } else { /* improvement in function val */
645 double rhocur = x[jdrop] - sim[jdrop + np * sim_dim1];
646 /* SGJ: use rhocur instead of rho. In original code, rhocur == rho
647 always, but I want to change this to ensure that simplex points
648 fall within [lb,ub]. */
649 sim[jdrop + np * sim_dim1] = x[jdrop];
651 for (k = 1; k <= i__1; ++k) {
652 datmat[k + jdrop * datmat_dim1] = datmat[k + np * datmat_dim1]
654 datmat[k + np * datmat_dim1] = con[k];
657 for (k = 1; k <= i__1; ++k) {
658 sim[jdrop + k * sim_dim1] = -rhocur;
661 for (i__ = k; i__ <= i__2; ++i__) {
662 temp -= simi[i__ + k * simi_dim1];
664 simi[jdrop + k * simi_dim1] = temp;
668 if (stop->nevals <= *n) { /* evaluating initial simplex */
669 jdrop = stop->nevals;
670 /* SGJ: was += rho, but using sim[jdrop,jdrop] enforces consistency
671 if we change the stepsize above to stay in [lb,ub]. */
672 x[jdrop] += sim[jdrop + jdrop * sim_dim1];
678 /* Identify the optimal vertex of the current simplex. */
681 phimin = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
685 for (j = 1; j <= i__1; ++j) {
686 temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
691 } else if (temp == phimin && parmu == 0.) {
692 if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + nbest *
699 /* Switch the best vertex into pole position if it is not there already, */
700 /* and also update SIM, SIMI and DATMAT. */
704 for (i__ = 1; i__ <= i__1; ++i__) {
705 temp = datmat[i__ + np * datmat_dim1];
706 datmat[i__ + np * datmat_dim1] = datmat[i__ + nbest * datmat_dim1]
708 datmat[i__ + nbest * datmat_dim1] = temp;
711 for (i__ = 1; i__ <= i__1; ++i__) {
712 temp = sim[i__ + nbest * sim_dim1];
713 sim[i__ + nbest * sim_dim1] = 0.;
714 sim[i__ + np * sim_dim1] += temp;
717 for (k = 1; k <= i__2; ++k) {
718 sim[i__ + k * sim_dim1] -= temp;
719 tempa -= simi[k + i__ * simi_dim1];
721 simi[nbest + i__ * simi_dim1] = tempa;
725 /* Make an error return if SIGI is a poor approximation to the inverse of */
726 /* the leading N by N submatrix of SIG. */
730 for (i__ = 1; i__ <= i__1; ++i__) {
732 for (j = 1; j <= i__2; ++j) {
738 for (k = 1; k <= i__3; ++k) if (sim[k + j * sim_dim1] != 0) {
739 temp += simi[i__ + k * simi_dim1] * sim[k + j * sim_dim1];
741 d__1 = error, d__2 = fabs(temp);
742 error = MAX2(d__1,d__2);
747 fprintf(stderr, "cobyla: rounding errors are becoming damaging.\n");
753 /* Calculate the coefficients of the linear approximations to the objective */
754 /* and constraint functions, placing minus the objective function gradient */
755 /* after the constraint gradients in the array A. The vector W is used for */
759 for (k = 1; k <= i__2; ++k) {
760 con[k] = -datmat[k + np * datmat_dim1];
762 for (j = 1; j <= i__1; ++j) {
763 w[j] = datmat[k + j * datmat_dim1] + con[k];
766 for (i__ = 1; i__ <= i__1; ++i__) {
769 for (j = 1; j <= i__3; ++j) {
770 temp += w[j] * simi[j + i__ * simi_dim1];
775 a[i__ + k * a_dim1] = temp;
779 /* Calculate the values of sigma and eta, and set IFLAG=0 if the current */
780 /* simplex is not acceptable. */
783 parsig = alpha * rho;
786 for (j = 1; j <= i__1; ++j) {
790 for (i__ = 1; i__ <= i__2; ++i__) {
791 d__1 = simi[j + i__ * simi_dim1];
793 d__1 = sim[i__ + j * sim_dim1];
796 vsig[j] = 1. / sqrt(wsig);
797 veta[j] = sqrt(weta);
798 if (vsig[j] < parsig || veta[j] > pareta) {
803 /* If a new vertex is needed to improve acceptability, then decide which */
804 /* vertex to drop from the simplex. */
806 if (ibrnch == 1 || iflag == 1) {
812 for (j = 1; j <= i__1; ++j) {
813 if (veta[j] > temp) {
820 for (j = 1; j <= i__1; ++j) {
821 if (vsig[j] < temp) {
828 /* Calculate the step to the new vertex and its sign. */
830 temp = gamma_ * rho * vsig[jdrop];
832 for (i__ = 1; i__ <= i__1; ++i__) {
833 dx[i__] = temp * simi[jdrop + i__ * simi_dim1];
838 for (k = 1; k <= i__1; ++k) {
841 for (i__ = 1; i__ <= i__2; ++i__) {
842 sum += a[i__ + k * a_dim1] * dx[i__];
845 temp = datmat[k + np * datmat_dim1];
846 d__1 = cvmaxp, d__2 = -sum - temp;
847 cvmaxp = MAX2(d__1,d__2);
848 d__1 = cvmaxm, d__2 = sum - temp;
849 cvmaxm = MAX2(d__1,d__2);
853 if (parmu * (cvmaxp - cvmaxm) > sum + sum) {
857 /* Update the elements of SIM and SIMI, and set the next X. */
861 for (i__ = 1; i__ <= i__1; ++i__) {
862 /* SGJ, 2010: pseudo-randomize simplex steps (see LCG comments above) */
863 dx[i__] = dxsign * dx[i__] * lcg_urand(&seed, 0.01, 1);
864 /* SGJ: make sure dx step says in [lb,ub] */
867 double xi = sim[i__ + np * sim_dim1];
869 if (xi + dx[i__] > ub[i__])
871 if (xi + dx[i__] < lb[i__]) {
872 if (xi - dx[i__] <= ub[i__])
874 else { /* try again with halved step */
881 sim[i__ + jdrop * sim_dim1] = dx[i__];
882 temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
885 for (i__ = 1; i__ <= i__1; ++i__) {
886 simi[jdrop + i__ * simi_dim1] /= temp;
889 for (j = 1; j <= i__1; ++j) {
893 for (i__ = 1; i__ <= i__2; ++i__) {
894 temp += simi[j + i__ * simi_dim1] * dx[i__];
897 for (i__ = 1; i__ <= i__2; ++i__) {
898 simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
902 x[j] = sim[j + np * sim_dim1] + dx[j];
906 /* Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO. */
910 izdota = iz + *n * *n;
913 idxnew = isdirn + *n;
915 rc = trstlp(n, m, &a[a_offset], &con[1], &rho, &dx[1], &ifull, &iact[1], &w[
916 iz], &w[izdota], &w[ivmc], &w[isdirn], &w[idxnew], &w[ivmd]);
917 if (rc != NLOPT_SUCCESS) goto L600;
919 /* SGJ: since the bound constraints are linear, we should never get
920 a dx that lies outside the [lb,ub] constraints here, but we'll
921 enforce this anyway just to be paranoid */
923 for (i__ = 1; i__ <= i__1; ++i__) {
924 double xi = sim[i__ + np * sim_dim1];
925 if (xi + dx[i__] > ub[i__]) dx[i__] = ub[i__] - xi;
926 if (xi + dx[i__] < lb[i__]) dx[i__] = xi - lb[i__];
932 for (i__ = 1; i__ <= i__1; ++i__) {
936 if (temp < rho * .25 * rho) {
942 /* Predict the change to F and the new maximum constraint violation if the */
943 /* variables are altered from x(0) to x(0)+DX. */
948 for (k = 1; k <= i__1; ++k) {
951 for (i__ = 1; i__ <= i__2; ++i__) {
952 sum -= a[i__ + k * a_dim1] * dx[i__];
955 resnew = MAX2(resnew,sum);
959 /* Increase PARMU if necessary and branch back if this change alters the */
960 /* optimal vertex. Otherwise PREREM and PREREC will be set to the predicted */
961 /* reductions in the merit function and the maximum constraint violation */
965 prerec = datmat[*mpp + np * datmat_dim1] - resnew;
967 barmu = sum / prerec;
969 if (parmu < barmu * 1.5) {
972 fprintf(stderr, "cobyla: increase in PARMU to %12.6E\n", parmu);
974 phi = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
977 for (j = 1; j <= i__1; ++j) {
978 temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
983 if (temp == phi && parmu == 0.f) {
984 if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + np *
991 prerem = parmu * prerec - sum;
993 /* Calculate the constraint and objective functions at x(*). Then find the */
994 /* actual reduction in the merit function. */
997 for (i__ = 1; i__ <= i__1; ++i__) {
998 x[i__] = sim[i__ + np * sim_dim1] + dx[i__];
1003 vmold = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
1005 vmnew = f + parmu * resmax;
1006 trured = vmold - vmnew;
1007 if (parmu == 0. && f == datmat[mp + np * datmat_dim1]) {
1009 trured = datmat[*mpp + np * datmat_dim1] - resmax;
1012 /* Begin the operations that decide whether x(*) should replace one of the */
1013 /* vertices of the current simplex, the change being mandatory if TRURED is */
1014 /* positive. Firstly, JDROP is set to the index of the vertex that is to be */
1018 if (trured <= 0.f) {
1023 for (j = 1; j <= i__1; ++j) {
1026 for (i__ = 1; i__ <= i__2; ++i__) {
1027 temp += simi[j + i__ * simi_dim1] * dx[i__];
1034 sigbar[j] = temp * vsig[j];
1037 /* Calculate the value of ell. */
1039 edgmax = delta * rho;
1042 for (j = 1; j <= i__1; ++j) {
1043 if (sigbar[j] >= parsig || sigbar[j] >= vsig[j]) {
1048 for (i__ = 1; i__ <= i__2; ++i__) {
1049 d__1 = dx[i__] - sim[i__ + j * sim_dim1];
1050 temp += d__1 * d__1;
1054 if (temp > edgmax) {
1067 /* Revise the simplex by updating the elements of SIM, SIMI and DATMAT. */
1071 for (i__ = 1; i__ <= i__1; ++i__) {
1072 sim[i__ + jdrop * sim_dim1] = dx[i__];
1073 temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
1076 for (i__ = 1; i__ <= i__1; ++i__) {
1077 simi[jdrop + i__ * simi_dim1] /= temp;
1080 for (j = 1; j <= i__1; ++j) {
1084 for (i__ = 1; i__ <= i__2; ++i__) {
1085 temp += simi[j + i__ * simi_dim1] * dx[i__];
1088 for (i__ = 1; i__ <= i__2; ++i__) {
1089 simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
1095 for (k = 1; k <= i__1; ++k) {
1096 datmat[k + jdrop * datmat_dim1] = con[k];
1099 /* Branch back for further iterations with the current RHO. */
1101 if (trured > 0. && trured >= prerem * .1) {
1102 /* SGJ, 2010: following a suggestion in the SAS manual (which
1103 mentions a similar modification to COBYLA, although they didn't
1104 publish their source code), increase rho if predicted reduction
1105 is sufficiently close to the actual reduction. Otherwise,
1106 COBLYA seems to easily get stuck making very small steps.
1107 Also require iflag != 0 (i.e., acceptable simplex), again
1108 following SAS suggestion (otherwise I get convergence failure
1110 if (trured >= prerem * 0.9 && trured <= prerem * 1.1 && iflag) {
1121 /* SGJ, 2008: convergence tests for function vals; note that current
1122 best val is stored in datmat[mp + np * datmat_dim1], or in f if
1123 ifull == 1, and previous best is in *minf. This seems like a
1124 sensible place to put the convergence test, as it is the same
1125 place that Powell checks the x tolerance (rho > rhoend). */
1127 double fbest = ifull == 1 ? f : datmat[mp + np * datmat_dim1];
1128 if (fbest < *minf && nlopt_stop_ftol(stop, fbest, *minf)) {
1129 rc = NLOPT_FTOL_REACHED;
1135 /* Otherwise reduce RHO if it is not at its least value and reset PARMU. */
1139 if (rho <= rhoend * 1.5) {
1145 for (k = 1; k <= i__1; ++k) {
1146 cmin = datmat[k + np * datmat_dim1];
1149 for (i__ = 1; i__ <= i__2; ++i__) {
1150 d__1 = cmin, d__2 = datmat[k + i__ * datmat_dim1];
1151 cmin = MIN2(d__1,d__2);
1152 d__1 = cmax, d__2 = datmat[k + i__ * datmat_dim1];
1153 cmax = MAX2(d__1,d__2);
1155 if (k <= *m && cmin < cmax * .5) {
1156 temp = MAX2(cmax,0.) - cmin;
1160 denom = MIN2(denom,temp);
1166 } else if (cmax - cmin < parmu * denom) {
1167 parmu = (cmax - cmin) / denom;
1171 fprintf(stderr, "cobyla: reduction in RHO to %12.6E and PARMU =%13.6E\n",
1175 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1176 stop->nevals, datmat[mp + np * datmat_dim1], datmat[*mpp + np * datmat_dim1]);
1178 fprintf(stderr, "cobyla: X =");
1180 for (i__ = 1; i__ <= i__1; ++i__) {
1181 if (i__>1) fprintf(stderr, " ");
1182 fprintf(stderr, "%13.6E", sim[i__ + np * sim_dim1]);
1186 for (i__ = iptemp; i__ <= i__1; ++i__) {
1187 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
1188 fprintf(stderr, "%15.6E", x[i__]);
1191 fprintf(stderr, "\n");
1196 rc = NLOPT_XTOL_REACHED;
1198 /* Return the best calculated values of the variables. */
1201 fprintf(stderr, "cobyla: normal return.\n");
1208 for (i__ = 1; i__ <= i__1; ++i__) {
1209 x[i__] = sim[i__ + np * sim_dim1];
1211 f = datmat[mp + np * datmat_dim1];
1212 resmax = datmat[*mpp + np * datmat_dim1];
1216 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1217 stop->nevals, f, resmax);
1219 fprintf(stderr, "cobyla: X =");
1220 for (i__ = 1; i__ <= i__1; ++i__) {
1221 if (i__>1) fprintf(stderr, " ");
1222 fprintf(stderr, "%13.6E", x[i__]);
1226 for (i__ = iptemp; i__ <= i__1; ++i__) {
1227 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
1228 fprintf(stderr, "%15.6E", x[i__]);
1231 fprintf(stderr, "\n");
1236 /* ------------------------------------------------------------------------- */
1237 static nlopt_result trstlp(int *n, int *m, double *a,
1238 double *b, double *rho, double *dx, int *ifull,
1239 int *iact, double *z__, double *zdota, double *vmultc,
1240 double *sdirn, double *dxnew, double *vmultd)
1242 /* System generated locals */
1243 int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
1246 /* Local variables */
1247 double alpha, tempa;
1249 double optnew, stpful, sum, tot, acca, accb;
1250 double ratio, vsave, zdotv, zdotw, dd;
1252 double sp, ss, resold = 0.0, zdvabs, zdwabs, sumabs, resmax, optold;
1256 int iout, i__, j, k;
1260 int nact, icon = 0, mcon;
1264 /* This subroutine calculates an N-component vector DX by applying the */
1265 /* following two stages. In the first stage, DX is set to the shortest */
1266 /* vector that minimizes the greatest violation of the constraints */
1267 /* A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K=2,3,...,M, */
1268 /* subject to the Euclidean length of DX being at most RHO. If its length is */
1269 /* strictly less than RHO, then we use the resultant freedom in DX to */
1270 /* minimize the objective function */
1271 /* -A(1,M+1)*DX(1)-A(2,M+1)*DX(2)-...-A(N,M+1)*DX(N) */
1272 /* subject to no increase in any greatest constraint violation. This */
1273 /* notation allows the gradient of the objective function to be regarded as */
1274 /* the gradient of a constraint. Therefore the two stages are distinguished */
1275 /* by MCON .EQ. M and MCON .GT. M respectively. It is possible that a */
1276 /* degeneracy may prevent DX from attaining the target length RHO. Then the */
1277 /* value IFULL=0 would be set, but usually IFULL=1 on return. */
1279 /* In general NACT is the number of constraints in the active set and */
1280 /* IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT */
1281 /* contains a permutation of the remaining constraint indices. Further, Z is */
1282 /* an orthogonal matrix whose first NACT columns can be regarded as the */
1283 /* result of Gram-Schmidt applied to the active constraint gradients. For */
1284 /* J=1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th */
1285 /* column of Z with the gradient of the J-th active constraint. DX is the */
1286 /* current vector of variables and here the residuals of the active */
1287 /* constraints should be zero. Further, the active constraints have */
1288 /* nonnegative Lagrange multipliers that are held at the beginning of */
1289 /* VMULTC. The remainder of this vector holds the residuals of the inactive */
1290 /* constraints at DX, the ordering of the components of VMULTC being in */
1291 /* agreement with the permutation of the indices of the constraints that is */
1292 /* in IACT. All these residuals are nonnegative, which is achieved by the */
1293 /* shift RESMAX that makes the least residual zero. */
1295 /* Initialize Z and some other variables. The value of RESMAX will be */
1296 /* appropriate to DX=0, while ICON will be the index of a most violated */
1297 /* constraint if RESMAX is positive. Usually during the first stage the */
1298 /* vector SDIRN gives a search direction that reduces all the active */
1299 /* constraint violations by one simultaneously. */
1301 /* Parameter adjustments */
1303 z_offset = 1 + z_dim1 * 1;
1306 a_offset = 1 + a_dim1 * 1;
1323 for (i__ = 1; i__ <= i__1; ++i__) {
1325 for (j = 1; j <= i__2; ++j) {
1326 z__[i__ + j * z_dim1] = 0.;
1328 z__[i__ + i__ * z_dim1] = 1.;
1333 for (k = 1; k <= i__1; ++k) {
1334 if (b[k] > resmax) {
1340 for (k = 1; k <= i__1; ++k) {
1342 vmultc[k] = resmax - b[k];
1349 for (i__ = 1; i__ <= i__1; ++i__) {
1353 /* End the current stage of the calculation if 3 consecutive iterations */
1354 /* have either failed to reduce the best calculated value of the objective */
1355 /* function or to increase the number of active constraints since the best */
1356 /* value was calculated. This strategy prevents cycling, but there is a */
1357 /* remote possibility that it will cause premature termination. */
1368 for (i__ = 1; i__ <= i__1; ++i__) {
1369 optnew -= dx[i__] * a[i__ + mcon * a_dim1];
1372 if (icount == 0 || optnew < optold) {
1376 } else if (nact > nactx) {
1386 /* If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to */
1387 /* the active set. Apply Givens rotations so that the last N-NACT-1 columns */
1388 /* of Z are orthogonal to the gradient of the new constraint, a scalar */
1389 /* product being set to zero if its nonzero value could be due to computer */
1390 /* rounding errors. The array DXNEW is used for working space. */
1397 for (i__ = 1; i__ <= i__1; ++i__) {
1398 dxnew[i__] = a[i__ + kk * a_dim1];
1407 for (i__ = 1; i__ <= i__1; ++i__) {
1408 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1410 spabs += fabs(temp);
1412 acca = spabs + fabs(sp) * .1;
1413 accb = spabs + fabs(sp) * .2;
1414 if (spabs >= acca || acca >= accb) {
1421 temp = sqrt(sp * sp + tot * tot);
1426 for (i__ = 1; i__ <= i__1; ++i__) {
1427 temp = alpha * z__[i__ + k * z_dim1] + beta * z__[i__ + kp *
1429 z__[i__ + kp * z_dim1] = alpha * z__[i__ + kp * z_dim1] -
1430 beta * z__[i__ + k * z_dim1];
1431 z__[i__ + k * z_dim1] = temp;
1438 /* Add the new constraint if this can be done without a deletion from the */
1444 vmultc[icon] = vmultc[nact];
1449 /* The next instruction is reached if a deletion has to be made from the */
1450 /* active set in order to make room for the new active constraint, because */
1451 /* the new constraint gradient is a linear combination of the gradients of */
1452 /* the old active constraints. Set the elements of VMULTD to the multipliers */
1453 /* of the linear combination. Further, set IOUT to the index of the */
1454 /* constraint to be deleted, but branch if no suitable index can be found. */
1462 for (i__ = 1; i__ <= i__1; ++i__) {
1463 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1465 zdvabs += fabs(temp);
1467 acca = zdvabs + fabs(zdotv) * .1;
1468 accb = zdvabs + fabs(zdotv) * .2;
1469 if (zdvabs < acca && acca < accb) {
1470 temp = zdotv / zdota[k];
1471 if (temp > 0. && iact[k] <= *m) {
1472 tempa = vmultc[k] / temp;
1473 if (ratio < 0. || tempa < ratio) {
1481 for (i__ = 1; i__ <= i__1; ++i__) {
1482 dxnew[i__] -= temp * a[i__ + kw * a_dim1];
1497 /* Revise the Lagrange multipliers and reorder the active constraints so */
1498 /* that the one to be replaced is at the end of the list. Also calculate the */
1499 /* new value of ZDOTA(NACT) and branch if it is not acceptable. */
1502 for (k = 1; k <= i__1; ++k) {
1503 d__1 = 0., d__2 = vmultc[k] - ratio * vmultd[k];
1504 vmultc[k] = MAX2(d__1,d__2);
1508 vsave = vmultc[icon];
1515 for (i__ = 1; i__ <= i__1; ++i__) {
1516 sp += z__[i__ + k * z_dim1] * a[i__ + kw * a_dim1];
1519 temp = sqrt(sp * sp + d__1 * d__1);
1520 alpha = zdota[kp] / temp;
1522 zdota[kp] = alpha * zdota[k];
1525 for (i__ = 1; i__ <= i__1; ++i__) {
1526 temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1528 z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1529 z__[i__ + kp * z_dim1];
1530 z__[i__ + k * z_dim1] = temp;
1533 vmultc[k] = vmultc[kp];
1543 for (i__ = 1; i__ <= i__1; ++i__) {
1544 temp += z__[i__ + nact * z_dim1] * a[i__ + kk * a_dim1];
1551 vmultc[nact] = ratio;
1553 /* Update IACT and ensure that the objective function continues to be */
1554 /* treated as the last active constraint when MCON>M. */
1557 iact[icon] = iact[nact];
1559 if (mcon > *m && kk != mcon) {
1563 for (i__ = 1; i__ <= i__1; ++i__) {
1564 sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1567 temp = sqrt(sp * sp + d__1 * d__1);
1568 alpha = zdota[nact] / temp;
1570 zdota[nact] = alpha * zdota[k];
1573 for (i__ = 1; i__ <= i__1; ++i__) {
1574 temp = alpha * z__[i__ + nact * z_dim1] + beta * z__[i__ + k *
1576 z__[i__ + nact * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1577 z__[i__ + nact * z_dim1];
1578 z__[i__ + k * z_dim1] = temp;
1580 iact[nact] = iact[k];
1583 vmultc[k] = vmultc[nact];
1584 vmultc[nact] = temp;
1587 /* If stage one is in progress, then set SDIRN to the direction of the next */
1588 /* change to the current vector of variables. */
1596 for (i__ = 1; i__ <= i__1; ++i__) {
1597 temp += sdirn[i__] * a[i__ + kk * a_dim1];
1600 temp /= zdota[nact];
1602 for (i__ = 1; i__ <= i__1; ++i__) {
1603 sdirn[i__] -= temp * z__[i__ + nact * z_dim1];
1607 /* Delete the constraint that has the index IACT(ICON) from the active set. */
1612 vsave = vmultc[icon];
1619 for (i__ = 1; i__ <= i__1; ++i__) {
1620 sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1623 temp = sqrt(sp * sp + d__1 * d__1);
1624 alpha = zdota[kp] / temp;
1626 zdota[kp] = alpha * zdota[k];
1629 for (i__ = 1; i__ <= i__1; ++i__) {
1630 temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1632 z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1633 z__[i__ + kp * z_dim1];
1634 z__[i__ + k * z_dim1] = temp;
1637 vmultc[k] = vmultc[kp];
1647 /* If stage one is in progress, then set SDIRN to the direction of the next */
1648 /* change to the current vector of variables. */
1655 for (i__ = 1; i__ <= i__1; ++i__) {
1656 temp += sdirn[i__] * z__[i__ + (nact + 1) * z_dim1];
1659 for (i__ = 1; i__ <= i__1; ++i__) {
1660 sdirn[i__] -= temp * z__[i__ + (nact + 1) * z_dim1];
1664 /* Pick the next search direction of stage two. */
1667 temp = 1. / zdota[nact];
1669 for (i__ = 1; i__ <= i__1; ++i__) {
1670 sdirn[i__] = temp * z__[i__ + nact * z_dim1];
1673 /* Calculate the step to the boundary of the trust region or take the step */
1674 /* that reduces RESMAX to zero. The two statements below that include the */
1675 /* factor 1.0E-6 prevent some harmless underflows that occurred in a test */
1676 /* calculation. Further, we skip the step if it could be zero within a */
1677 /* reasonable tolerance for computer rounding errors. */
1684 for (i__ = 1; i__ <= i__1; ++i__) {
1685 if ((d__1 = dx[i__], fabs(d__1)) >= *rho * 1e-6f) {
1689 sd += dx[i__] * sdirn[i__];
1696 temp = sqrt(ss * dd);
1697 if (fabs(sd) >= temp * 1e-6f) {
1698 temp = sqrt(ss * dd + sd * sd);
1700 stpful = dd / (temp + sd);
1703 acca = step + resmax * .1;
1704 accb = step + resmax * .2;
1705 if (step >= acca || acca >= accb) {
1708 step = MIN2(step,resmax);
1711 /* SGJ, 2010: check for error here */
1712 if (nlopt_isinf(step)) return NLOPT_ROUNDOFF_LIMITED;
1714 /* Set DXNEW to the new variables if STEP is the steplength, and reduce */
1715 /* RESMAX to the corresponding maximum residual if stage one is being done. */
1716 /* Because DXNEW will be changed during the calculation of some Lagrange */
1717 /* multipliers, it will be restored to the following value later. */
1720 for (i__ = 1; i__ <= i__1; ++i__) {
1721 dxnew[i__] = dx[i__] + step * sdirn[i__];
1727 for (k = 1; k <= i__1; ++k) {
1731 for (i__ = 1; i__ <= i__2; ++i__) {
1732 temp -= a[i__ + kk * a_dim1] * dxnew[i__];
1734 resmax = MAX2(resmax,temp);
1738 /* Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A */
1739 /* device is included to force VMULTD(K)=0.0 if deviations from this value */
1740 /* can be attributed to computer rounding errors. First calculate the new */
1741 /* Lagrange multipliers. */
1748 for (i__ = 1; i__ <= i__1; ++i__) {
1749 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1751 zdwabs += fabs(temp);
1753 acca = zdwabs + fabs(zdotw) * .1;
1754 accb = zdwabs + fabs(zdotw) * .2;
1755 if (zdwabs >= acca || acca >= accb) {
1758 vmultd[k] = zdotw / zdota[k];
1762 for (i__ = 1; i__ <= i__1; ++i__) {
1763 dxnew[i__] -= vmultd[k] * a[i__ + kk * a_dim1];
1769 d__1 = 0., d__2 = vmultd[nact];
1770 vmultd[nact] = MAX2(d__1,d__2);
1773 /* Complete VMULTC by finding the new constraint residuals. */
1776 for (i__ = 1; i__ <= i__1; ++i__) {
1777 dxnew[i__] = dx[i__] + step * sdirn[i__];
1782 for (k = kl; k <= i__1; ++k) {
1784 sum = resmax - b[kk];
1785 sumabs = resmax + (d__1 = b[kk], fabs(d__1));
1787 for (i__ = 1; i__ <= i__2; ++i__) {
1788 temp = a[i__ + kk * a_dim1] * dxnew[i__];
1790 sumabs += fabs(temp);
1792 acca = sumabs + fabs(sum) * .1f;
1793 accb = sumabs + fabs(sum) * .2f;
1794 if (sumabs >= acca || acca >= accb) {
1801 /* Calculate the fraction of the step from DX to DXNEW that will be taken. */
1806 for (k = 1; k <= i__1; ++k) {
1807 if (vmultd[k] < 0.) {
1808 temp = vmultc[k] / (vmultc[k] - vmultd[k]);
1816 /* Update DX, VMULTC and RESMAX. */
1820 for (i__ = 1; i__ <= i__1; ++i__) {
1821 dx[i__] = temp * dx[i__] + ratio * dxnew[i__];
1824 for (k = 1; k <= i__1; ++k) {
1825 d__1 = 0., d__2 = temp * vmultc[k] + ratio * vmultd[k];
1826 vmultc[k] = MAX2(d__1,d__2);
1829 resmax = resold + ratio * (resmax - resold);
1832 /* If the full step is not acceptable then begin another iteration. */
1833 /* Otherwise switch to stage two or end the calculation. */
1838 if (step == stpful) {
1848 /* We employ any freedom that may be available to reduce the objective */
1849 /* function before returning a DX whose length is less than RHO. */
1857 return NLOPT_SUCCESS;