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 min(a,b) ((a) <= (b) ? (a) : (b))
58 #define max(a,b) ((a) >= (b) ? (a) : (b))
59 #define abs(x) ((x) >= 0 ? (x) : -(x))
61 #define U(n) ((unsigned) (n))
63 /**************************************************************************/
64 /* SGJ, 2008: NLopt-style interface function: */
75 double *con_tol, *scale;
79 static int func_wrap(int ni, int mi, double *x, double *f, double *con,
84 double *xtmp = s->xtmp;
85 const double *lb = s->lb, *ub = s->ub;
87 (void) mi; /* unused */
89 /* in nlopt, we guarante that the function is never evaluated outside
90 the lb and ub bounds, so we need force this with xtmp ... note
91 that this leads to discontinuity in the first derivative, which
92 slows convergence if we don't enable the ENFORCE_BOUNDS feature
94 for (j = 0; j < n; ++j) {
95 if (x[j] < lb[j]) xtmp[j] = lb[j];
96 else if (x[j] > ub[j]) xtmp[j] = ub[j];
99 nlopt_unscale(n, s->scale, xtmp, xtmp);
101 *f = s->f(n, xtmp, NULL, s->f_data);
102 if (nlopt_stop_forced(s->stop)) return 1;
104 for (j = 0; j < s->m_orig; ++j) {
105 nlopt_eval_constraint(con + i, NULL, s->fc+j, n, xtmp);
106 if (nlopt_stop_forced(s->stop)) return 1;
107 for (k = 0; k < s->fc[j].m; ++k)
108 con[i + k] = -con[i + k];
111 for (j = 0; j < s->p; ++j) {
112 nlopt_eval_constraint(con + i, NULL, s->h+j, n, xtmp);
113 if (nlopt_stop_forced(s->stop)) return 1;
114 for (k = 0; k < s->h[j].m; ++k)
115 con[(i + s->h[j].m) + k] = -con[i + k];
118 for (j = 0; j < n; ++j) {
119 if (!nlopt_isinf(lb[j]))
120 con[i++] = x[j] - lb[j];
121 if (!nlopt_isinf(ub[j]))
122 con[i++] = ub[j] - x[j];
131 COBYLA_MSG_NONE = 0, /* No messages */
132 COBYLA_MSG_EXIT = 1, /* Exit reasons */
133 COBYLA_MSG_ITER = 2, /* Rho and Sigma changes */
134 COBYLA_MSG_INFO = 3 /* Informational messages */
138 * A function as required by cobyla
139 * state is a void pointer provided to the function at each call
141 * n : the number of variables
142 * m : the number of constraints
143 * x : on input, then vector of variables (should not be modified)
144 * f : on output, the value of the function
145 * con : on output, the value of the constraints (vector of size m)
146 * state : on input, the value of the state variable as provided to cobyla
148 * COBYLA will try to make all the values of the constraints positive.
149 * So if you want to input a constraint j such as x[i] <= MAX, set:
150 * con[j] = MAX - x[i]
151 * The function must returns 0 if no error occurs or 1 to immediately end the
155 typedef int cobyla_function(int n, int m, double *x, double *f, double *con,
156 func_wrap_state *state);
159 * cobyla : minimize a function subject to constraints
161 * n : number of variables (>=0)
162 * m : number of constraints (>=0)
163 * x : on input, initial estimate ; on output, the solution
164 * minf : on output, minimum objective function
165 * rhobeg : a reasonable initial change to the variables
166 * stop : the NLopt stopping criteria
167 * lb, ub : lower and upper bounds on x
168 * message : see the cobyla_message enum
169 * calcfc : the function to minimize (see cobyla_function)
170 * state : used by function (see cobyla_function)
172 * The cobyla function returns the usual nlopt_result codes.
175 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,
176 int message, cobyla_function *calcfc, func_wrap_state *state);
178 nlopt_result cobyla_minimize(unsigned n, nlopt_func f, void *f_data,
179 unsigned m, nlopt_constraint *fc,
180 unsigned p, nlopt_constraint *h,
181 const double *lb, const double *ub, /* bounds */
182 double *x, /* in: initial guess, out: minimizer */
184 nlopt_stopping *stop,
190 double rhobeg, rhoend;
192 s.f = f; s.f_data = f_data;
198 s.lb = s.ub = s.xtmp = s.con_tol = s.scale = NULL;
200 s.scale = nlopt_compute_rescaling(n, dx);
201 if (!s.scale) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
203 s.lb = nlopt_new_rescaled(n, s.scale, lb);
204 if (!s.lb) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
205 s.ub = nlopt_new_rescaled(n, s.scale, ub);
206 if (!s.ub) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
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 = 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] / s.scale[j])
216 rhoend = stop->xtol_abs[j] / 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 (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
567 else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
568 if (rc != NLOPT_SUCCESS) goto L600;
571 if (calcfc(*n, *m, &x[1], &f, &con[1], state))
574 fprintf(stderr, "cobyla: user requested end of minimization.\n");
576 rc = NLOPT_FORCED_STOP;
581 feasible = 1; /* SGJ, 2010 */
584 for (k = 1; k <= i__1; ++k) {
585 d__1 = resmax, d__2 = -con[k];
586 resmax = max(d__1,d__2);
587 if (d__2 > state->con_tol[k-1])
588 feasible = 0; /* SGJ, 2010 */
592 /* SGJ, 2008: check for minf_max reached by feasible point */
593 if (f < stop->minf_max && feasible) {
594 rc = NLOPT_MINF_MAX_REACHED;
595 goto L620; /* not L600 because we want to use current x, f, resmax */
598 if (stop->nevals == *iprint - 1 || *iprint == 3) {
599 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
600 stop->nevals, f, resmax);
602 fprintf(stderr, "cobyla: X =");
603 for (i__ = 1; i__ <= i__1; ++i__) {
604 if (i__>1) fprintf(stderr, " ");
605 fprintf(stderr, "%13.6E", x[i__]);
609 for (i__ = iptemp; i__ <= i__1; ++i__) {
610 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
611 fprintf(stderr, "%15.6E", x[i__]);
614 fprintf(stderr, "\n");
622 /* Set the recently calculated function values in a column of DATMAT. This */
623 /* array has a column for each vertex of the current simplex, the entries of */
624 /* each column being the values of the constraint functions (if any) */
625 /* followed by the objective function and the greatest constraint violation */
629 for (k = 1; k <= i__1; ++k) {
630 datmat[k + jdrop * datmat_dim1] = con[k];
632 if (stop->nevals > np) {
636 /* Exchange the new vertex of the initial simplex with the optimal vertex if */
637 /* necessary. Then, if the initial simplex is not complete, pick its next */
638 /* vertex and calculate the function values there. */
641 if (datmat[mp + np * datmat_dim1] <= f) {
642 x[jdrop] = sim[jdrop + np * sim_dim1];
643 } else { /* improvement in function val */
644 double rhocur = x[jdrop] - sim[jdrop + np * sim_dim1];
645 /* SGJ: use rhocur instead of rho. In original code, rhocur == rho
646 always, but I want to change this to ensure that simplex points
647 fall within [lb,ub]. */
648 sim[jdrop + np * sim_dim1] = x[jdrop];
650 for (k = 1; k <= i__1; ++k) {
651 datmat[k + jdrop * datmat_dim1] = datmat[k + np * datmat_dim1]
653 datmat[k + np * datmat_dim1] = con[k];
656 for (k = 1; k <= i__1; ++k) {
657 sim[jdrop + k * sim_dim1] = -rhocur;
660 for (i__ = k; i__ <= i__2; ++i__) {
661 temp -= simi[i__ + k * simi_dim1];
663 simi[jdrop + k * simi_dim1] = temp;
667 if (stop->nevals <= *n) { /* evaluating initial simplex */
668 jdrop = stop->nevals;
669 /* SGJ: was += rho, but using sim[jdrop,jdrop] enforces consistency
670 if we change the stepsize above to stay in [lb,ub]. */
671 x[jdrop] += sim[jdrop + jdrop * sim_dim1];
677 /* Identify the optimal vertex of the current simplex. */
680 phimin = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
684 for (j = 1; j <= i__1; ++j) {
685 temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
690 } else if (temp == phimin && parmu == 0.) {
691 if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + nbest *
698 /* Switch the best vertex into pole position if it is not there already, */
699 /* and also update SIM, SIMI and DATMAT. */
703 for (i__ = 1; i__ <= i__1; ++i__) {
704 temp = datmat[i__ + np * datmat_dim1];
705 datmat[i__ + np * datmat_dim1] = datmat[i__ + nbest * datmat_dim1]
707 datmat[i__ + nbest * datmat_dim1] = temp;
710 for (i__ = 1; i__ <= i__1; ++i__) {
711 temp = sim[i__ + nbest * sim_dim1];
712 sim[i__ + nbest * sim_dim1] = 0.;
713 sim[i__ + np * sim_dim1] += temp;
716 for (k = 1; k <= i__2; ++k) {
717 sim[i__ + k * sim_dim1] -= temp;
718 tempa -= simi[k + i__ * simi_dim1];
720 simi[nbest + i__ * simi_dim1] = tempa;
724 /* Make an error return if SIGI is a poor approximation to the inverse of */
725 /* the leading N by N submatrix of SIG. */
729 for (i__ = 1; i__ <= i__1; ++i__) {
731 for (j = 1; j <= i__2; ++j) {
737 for (k = 1; k <= i__3; ++k) if (sim[k + j * sim_dim1] != 0) {
738 temp += simi[i__ + k * simi_dim1] * sim[k + j * sim_dim1];
740 d__1 = error, d__2 = abs(temp);
741 error = max(d__1,d__2);
746 fprintf(stderr, "cobyla: rounding errors are becoming damaging.\n");
752 /* Calculate the coefficients of the linear approximations to the objective */
753 /* and constraint functions, placing minus the objective function gradient */
754 /* after the constraint gradients in the array A. The vector W is used for */
758 for (k = 1; k <= i__2; ++k) {
759 con[k] = -datmat[k + np * datmat_dim1];
761 for (j = 1; j <= i__1; ++j) {
762 w[j] = datmat[k + j * datmat_dim1] + con[k];
765 for (i__ = 1; i__ <= i__1; ++i__) {
768 for (j = 1; j <= i__3; ++j) {
769 temp += w[j] * simi[j + i__ * simi_dim1];
774 a[i__ + k * a_dim1] = temp;
778 /* Calculate the values of sigma and eta, and set IFLAG=0 if the current */
779 /* simplex is not acceptable. */
782 parsig = alpha * rho;
785 for (j = 1; j <= i__1; ++j) {
789 for (i__ = 1; i__ <= i__2; ++i__) {
790 d__1 = simi[j + i__ * simi_dim1];
792 d__1 = sim[i__ + j * sim_dim1];
795 vsig[j] = 1. / sqrt(wsig);
796 veta[j] = sqrt(weta);
797 if (vsig[j] < parsig || veta[j] > pareta) {
802 /* If a new vertex is needed to improve acceptability, then decide which */
803 /* vertex to drop from the simplex. */
805 if (ibrnch == 1 || iflag == 1) {
811 for (j = 1; j <= i__1; ++j) {
812 if (veta[j] > temp) {
819 for (j = 1; j <= i__1; ++j) {
820 if (vsig[j] < temp) {
827 /* Calculate the step to the new vertex and its sign. */
829 temp = gamma_ * rho * vsig[jdrop];
831 for (i__ = 1; i__ <= i__1; ++i__) {
832 dx[i__] = temp * simi[jdrop + i__ * simi_dim1];
837 for (k = 1; k <= i__1; ++k) {
840 for (i__ = 1; i__ <= i__2; ++i__) {
841 sum += a[i__ + k * a_dim1] * dx[i__];
844 temp = datmat[k + np * datmat_dim1];
845 d__1 = cvmaxp, d__2 = -sum - temp;
846 cvmaxp = max(d__1,d__2);
847 d__1 = cvmaxm, d__2 = sum - temp;
848 cvmaxm = max(d__1,d__2);
852 if (parmu * (cvmaxp - cvmaxm) > sum + sum) {
856 /* Update the elements of SIM and SIMI, and set the next X. */
860 for (i__ = 1; i__ <= i__1; ++i__) {
861 /* SGJ, 2010: pseudo-randomize simplex steps (see LCG comments above) */
862 dx[i__] = dxsign * dx[i__] * lcg_urand(&seed, 0.01, 1);
863 /* SGJ: make sure dx step says in [lb,ub] */
866 double xi = sim[i__ + np * sim_dim1];
868 if (xi + dx[i__] > ub[i__])
870 if (xi + dx[i__] < lb[i__]) {
871 if (xi - dx[i__] <= ub[i__])
873 else { /* try again with halved step */
880 sim[i__ + jdrop * sim_dim1] = dx[i__];
881 temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
884 for (i__ = 1; i__ <= i__1; ++i__) {
885 simi[jdrop + i__ * simi_dim1] /= temp;
888 for (j = 1; j <= i__1; ++j) {
892 for (i__ = 1; i__ <= i__2; ++i__) {
893 temp += simi[j + i__ * simi_dim1] * dx[i__];
896 for (i__ = 1; i__ <= i__2; ++i__) {
897 simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
901 x[j] = sim[j + np * sim_dim1] + dx[j];
905 /* Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO. */
909 izdota = iz + *n * *n;
912 idxnew = isdirn + *n;
914 rc = trstlp(n, m, &a[a_offset], &con[1], &rho, &dx[1], &ifull, &iact[1], &w[
915 iz], &w[izdota], &w[ivmc], &w[isdirn], &w[idxnew], &w[ivmd]);
916 if (rc != NLOPT_SUCCESS) goto L600;
918 /* SGJ: since the bound constraints are linear, we should never get
919 a dx that lies outside the [lb,ub] constraints here, but we'll
920 enforce this anyway just to be paranoid */
922 for (i__ = 1; i__ <= i__1; ++i__) {
923 double xi = sim[i__ + np * sim_dim1];
924 if (xi + dx[i__] > ub[i__]) dx[i__] = ub[i__] - xi;
925 if (xi + dx[i__] < lb[i__]) dx[i__] = xi - lb[i__];
931 for (i__ = 1; i__ <= i__1; ++i__) {
935 if (temp < rho * .25 * rho) {
941 /* Predict the change to F and the new maximum constraint violation if the */
942 /* variables are altered from x(0) to x(0)+DX. */
947 for (k = 1; k <= i__1; ++k) {
950 for (i__ = 1; i__ <= i__2; ++i__) {
951 sum -= a[i__ + k * a_dim1] * dx[i__];
954 resnew = max(resnew,sum);
958 /* Increase PARMU if necessary and branch back if this change alters the */
959 /* optimal vertex. Otherwise PREREM and PREREC will be set to the predicted */
960 /* reductions in the merit function and the maximum constraint violation */
964 prerec = datmat[*mpp + np * datmat_dim1] - resnew;
966 barmu = sum / prerec;
968 if (parmu < barmu * 1.5) {
971 fprintf(stderr, "cobyla: increase in PARMU to %12.6E\n", parmu);
973 phi = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
976 for (j = 1; j <= i__1; ++j) {
977 temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
982 if (temp == phi && parmu == 0.f) {
983 if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + np *
990 prerem = parmu * prerec - sum;
992 /* Calculate the constraint and objective functions at x(*). Then find the */
993 /* actual reduction in the merit function. */
996 for (i__ = 1; i__ <= i__1; ++i__) {
997 x[i__] = sim[i__ + np * sim_dim1] + dx[i__];
1002 vmold = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
1004 vmnew = f + parmu * resmax;
1005 trured = vmold - vmnew;
1006 if (parmu == 0. && f == datmat[mp + np * datmat_dim1]) {
1008 trured = datmat[*mpp + np * datmat_dim1] - resmax;
1011 /* Begin the operations that decide whether x(*) should replace one of the */
1012 /* vertices of the current simplex, the change being mandatory if TRURED is */
1013 /* positive. Firstly, JDROP is set to the index of the vertex that is to be */
1017 if (trured <= 0.f) {
1022 for (j = 1; j <= i__1; ++j) {
1025 for (i__ = 1; i__ <= i__2; ++i__) {
1026 temp += simi[j + i__ * simi_dim1] * dx[i__];
1033 sigbar[j] = temp * vsig[j];
1036 /* Calculate the value of ell. */
1038 edgmax = delta * rho;
1041 for (j = 1; j <= i__1; ++j) {
1042 if (sigbar[j] >= parsig || sigbar[j] >= vsig[j]) {
1047 for (i__ = 1; i__ <= i__2; ++i__) {
1048 d__1 = dx[i__] - sim[i__ + j * sim_dim1];
1049 temp += d__1 * d__1;
1053 if (temp > edgmax) {
1066 /* Revise the simplex by updating the elements of SIM, SIMI and DATMAT. */
1070 for (i__ = 1; i__ <= i__1; ++i__) {
1071 sim[i__ + jdrop * sim_dim1] = dx[i__];
1072 temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
1075 for (i__ = 1; i__ <= i__1; ++i__) {
1076 simi[jdrop + i__ * simi_dim1] /= temp;
1079 for (j = 1; j <= i__1; ++j) {
1083 for (i__ = 1; i__ <= i__2; ++i__) {
1084 temp += simi[j + i__ * simi_dim1] * dx[i__];
1087 for (i__ = 1; i__ <= i__2; ++i__) {
1088 simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
1094 for (k = 1; k <= i__1; ++k) {
1095 datmat[k + jdrop * datmat_dim1] = con[k];
1098 /* Branch back for further iterations with the current RHO. */
1100 if (trured > 0. && trured >= prerem * .1) {
1101 /* SGJ, 2010: following a suggestion in the SAS manual (which
1102 mentions a similar modification to COBYLA, although they didn't
1103 publish their source code), increase rho if predicted reduction
1104 is sufficiently close to the actual reduction. Otherwise,
1105 COBLYA seems to easily get stuck making very small steps.
1106 Also require iflag != 0 (i.e., acceptable simplex), again
1107 following SAS suggestion (otherwise I get convergence failure
1109 if (trured >= prerem * 0.9 && trured <= prerem * 1.1 && iflag) {
1120 /* SGJ, 2008: convergence tests for function vals; note that current
1121 best val is stored in datmat[mp + np * datmat_dim1], or in f if
1122 ifull == 1, and previous best is in *minf. This seems like a
1123 sensible place to put the convergence test, as it is the same
1124 place that Powell checks the x tolerance (rho > rhoend). */
1126 double fbest = ifull == 1 ? f : datmat[mp + np * datmat_dim1];
1127 if (fbest < *minf && nlopt_stop_ftol(stop, fbest, *minf)) {
1128 rc = NLOPT_FTOL_REACHED;
1134 /* Otherwise reduce RHO if it is not at its least value and reset PARMU. */
1138 if (rho <= rhoend * 1.5) {
1144 for (k = 1; k <= i__1; ++k) {
1145 cmin = datmat[k + np * datmat_dim1];
1148 for (i__ = 1; i__ <= i__2; ++i__) {
1149 d__1 = cmin, d__2 = datmat[k + i__ * datmat_dim1];
1150 cmin = min(d__1,d__2);
1151 d__1 = cmax, d__2 = datmat[k + i__ * datmat_dim1];
1152 cmax = max(d__1,d__2);
1154 if (k <= *m && cmin < cmax * .5) {
1155 temp = max(cmax,0.) - cmin;
1159 denom = min(denom,temp);
1165 } else if (cmax - cmin < parmu * denom) {
1166 parmu = (cmax - cmin) / denom;
1170 fprintf(stderr, "cobyla: reduction in RHO to %12.6E and PARMU =%13.6E\n",
1174 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1175 stop->nevals, datmat[mp + np * datmat_dim1], datmat[*mpp + np * datmat_dim1]);
1177 fprintf(stderr, "cobyla: X =");
1179 for (i__ = 1; i__ <= i__1; ++i__) {
1180 if (i__>1) fprintf(stderr, " ");
1181 fprintf(stderr, "%13.6E", sim[i__ + np * sim_dim1]);
1185 for (i__ = iptemp; i__ <= i__1; ++i__) {
1186 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
1187 fprintf(stderr, "%15.6E", x[i__]);
1190 fprintf(stderr, "\n");
1195 rc = NLOPT_XTOL_REACHED;
1197 /* Return the best calculated values of the variables. */
1200 fprintf(stderr, "cobyla: normal return.\n");
1207 for (i__ = 1; i__ <= i__1; ++i__) {
1208 x[i__] = sim[i__ + np * sim_dim1];
1210 f = datmat[mp + np * datmat_dim1];
1211 resmax = datmat[*mpp + np * datmat_dim1];
1215 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1216 stop->nevals, f, resmax);
1218 fprintf(stderr, "cobyla: X =");
1219 for (i__ = 1; i__ <= i__1; ++i__) {
1220 if (i__>1) fprintf(stderr, " ");
1221 fprintf(stderr, "%13.6E", x[i__]);
1225 for (i__ = iptemp; i__ <= i__1; ++i__) {
1226 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
1227 fprintf(stderr, "%15.6E", x[i__]);
1230 fprintf(stderr, "\n");
1235 /* ------------------------------------------------------------------------- */
1236 static nlopt_result trstlp(int *n, int *m, double *a,
1237 double *b, double *rho, double *dx, int *ifull,
1238 int *iact, double *z__, double *zdota, double *vmultc,
1239 double *sdirn, double *dxnew, double *vmultd)
1241 /* System generated locals */
1242 int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
1245 /* Local variables */
1246 double alpha, tempa;
1248 double optnew, stpful, sum, tot, acca, accb;
1249 double ratio, vsave, zdotv, zdotw, dd;
1251 double sp, ss, resold = 0.0, zdvabs, zdwabs, sumabs, resmax, optold;
1255 int iout, i__, j, k;
1259 int nact, icon = 0, mcon;
1263 /* This subroutine calculates an N-component vector DX by applying the */
1264 /* following two stages. In the first stage, DX is set to the shortest */
1265 /* vector that minimizes the greatest violation of the constraints */
1266 /* A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K=2,3,...,M, */
1267 /* subject to the Euclidean length of DX being at most RHO. If its length is */
1268 /* strictly less than RHO, then we use the resultant freedom in DX to */
1269 /* minimize the objective function */
1270 /* -A(1,M+1)*DX(1)-A(2,M+1)*DX(2)-...-A(N,M+1)*DX(N) */
1271 /* subject to no increase in any greatest constraint violation. This */
1272 /* notation allows the gradient of the objective function to be regarded as */
1273 /* the gradient of a constraint. Therefore the two stages are distinguished */
1274 /* by MCON .EQ. M and MCON .GT. M respectively. It is possible that a */
1275 /* degeneracy may prevent DX from attaining the target length RHO. Then the */
1276 /* value IFULL=0 would be set, but usually IFULL=1 on return. */
1278 /* In general NACT is the number of constraints in the active set and */
1279 /* IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT */
1280 /* contains a permutation of the remaining constraint indices. Further, Z is */
1281 /* an orthogonal matrix whose first NACT columns can be regarded as the */
1282 /* result of Gram-Schmidt applied to the active constraint gradients. For */
1283 /* J=1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th */
1284 /* column of Z with the gradient of the J-th active constraint. DX is the */
1285 /* current vector of variables and here the residuals of the active */
1286 /* constraints should be zero. Further, the active constraints have */
1287 /* nonnegative Lagrange multipliers that are held at the beginning of */
1288 /* VMULTC. The remainder of this vector holds the residuals of the inactive */
1289 /* constraints at DX, the ordering of the components of VMULTC being in */
1290 /* agreement with the permutation of the indices of the constraints that is */
1291 /* in IACT. All these residuals are nonnegative, which is achieved by the */
1292 /* shift RESMAX that makes the least residual zero. */
1294 /* Initialize Z and some other variables. The value of RESMAX will be */
1295 /* appropriate to DX=0, while ICON will be the index of a most violated */
1296 /* constraint if RESMAX is positive. Usually during the first stage the */
1297 /* vector SDIRN gives a search direction that reduces all the active */
1298 /* constraint violations by one simultaneously. */
1300 /* Parameter adjustments */
1302 z_offset = 1 + z_dim1 * 1;
1305 a_offset = 1 + a_dim1 * 1;
1322 for (i__ = 1; i__ <= i__1; ++i__) {
1324 for (j = 1; j <= i__2; ++j) {
1325 z__[i__ + j * z_dim1] = 0.;
1327 z__[i__ + i__ * z_dim1] = 1.;
1332 for (k = 1; k <= i__1; ++k) {
1333 if (b[k] > resmax) {
1339 for (k = 1; k <= i__1; ++k) {
1341 vmultc[k] = resmax - b[k];
1348 for (i__ = 1; i__ <= i__1; ++i__) {
1352 /* End the current stage of the calculation if 3 consecutive iterations */
1353 /* have either failed to reduce the best calculated value of the objective */
1354 /* function or to increase the number of active constraints since the best */
1355 /* value was calculated. This strategy prevents cycling, but there is a */
1356 /* remote possibility that it will cause premature termination. */
1367 for (i__ = 1; i__ <= i__1; ++i__) {
1368 optnew -= dx[i__] * a[i__ + mcon * a_dim1];
1371 if (icount == 0 || optnew < optold) {
1375 } else if (nact > nactx) {
1385 /* If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to */
1386 /* the active set. Apply Givens rotations so that the last N-NACT-1 columns */
1387 /* of Z are orthogonal to the gradient of the new constraint, a scalar */
1388 /* product being set to zero if its nonzero value could be due to computer */
1389 /* rounding errors. The array DXNEW is used for working space. */
1396 for (i__ = 1; i__ <= i__1; ++i__) {
1397 dxnew[i__] = a[i__ + kk * a_dim1];
1406 for (i__ = 1; i__ <= i__1; ++i__) {
1407 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1411 acca = spabs + abs(sp) * .1;
1412 accb = spabs + abs(sp) * .2;
1413 if (spabs >= acca || acca >= accb) {
1420 temp = sqrt(sp * sp + tot * tot);
1425 for (i__ = 1; i__ <= i__1; ++i__) {
1426 temp = alpha * z__[i__ + k * z_dim1] + beta * z__[i__ + kp *
1428 z__[i__ + kp * z_dim1] = alpha * z__[i__ + kp * z_dim1] -
1429 beta * z__[i__ + k * z_dim1];
1430 z__[i__ + k * z_dim1] = temp;
1437 /* Add the new constraint if this can be done without a deletion from the */
1443 vmultc[icon] = vmultc[nact];
1448 /* The next instruction is reached if a deletion has to be made from the */
1449 /* active set in order to make room for the new active constraint, because */
1450 /* the new constraint gradient is a linear combination of the gradients of */
1451 /* the old active constraints. Set the elements of VMULTD to the multipliers */
1452 /* of the linear combination. Further, set IOUT to the index of the */
1453 /* constraint to be deleted, but branch if no suitable index can be found. */
1461 for (i__ = 1; i__ <= i__1; ++i__) {
1462 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1464 zdvabs += abs(temp);
1466 acca = zdvabs + abs(zdotv) * .1;
1467 accb = zdvabs + abs(zdotv) * .2;
1468 if (zdvabs < acca && acca < accb) {
1469 temp = zdotv / zdota[k];
1470 if (temp > 0. && iact[k] <= *m) {
1471 tempa = vmultc[k] / temp;
1472 if (ratio < 0. || tempa < ratio) {
1480 for (i__ = 1; i__ <= i__1; ++i__) {
1481 dxnew[i__] -= temp * a[i__ + kw * a_dim1];
1496 /* Revise the Lagrange multipliers and reorder the active constraints so */
1497 /* that the one to be replaced is at the end of the list. Also calculate the */
1498 /* new value of ZDOTA(NACT) and branch if it is not acceptable. */
1501 for (k = 1; k <= i__1; ++k) {
1502 d__1 = 0., d__2 = vmultc[k] - ratio * vmultd[k];
1503 vmultc[k] = max(d__1,d__2);
1507 vsave = vmultc[icon];
1514 for (i__ = 1; i__ <= i__1; ++i__) {
1515 sp += z__[i__ + k * z_dim1] * a[i__ + kw * a_dim1];
1518 temp = sqrt(sp * sp + d__1 * d__1);
1519 alpha = zdota[kp] / temp;
1521 zdota[kp] = alpha * zdota[k];
1524 for (i__ = 1; i__ <= i__1; ++i__) {
1525 temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1527 z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1528 z__[i__ + kp * z_dim1];
1529 z__[i__ + k * z_dim1] = temp;
1532 vmultc[k] = vmultc[kp];
1542 for (i__ = 1; i__ <= i__1; ++i__) {
1543 temp += z__[i__ + nact * z_dim1] * a[i__ + kk * a_dim1];
1550 vmultc[nact] = ratio;
1552 /* Update IACT and ensure that the objective function continues to be */
1553 /* treated as the last active constraint when MCON>M. */
1556 iact[icon] = iact[nact];
1558 if (mcon > *m && kk != mcon) {
1562 for (i__ = 1; i__ <= i__1; ++i__) {
1563 sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1566 temp = sqrt(sp * sp + d__1 * d__1);
1567 alpha = zdota[nact] / temp;
1569 zdota[nact] = alpha * zdota[k];
1572 for (i__ = 1; i__ <= i__1; ++i__) {
1573 temp = alpha * z__[i__ + nact * z_dim1] + beta * z__[i__ + k *
1575 z__[i__ + nact * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1576 z__[i__ + nact * z_dim1];
1577 z__[i__ + k * z_dim1] = temp;
1579 iact[nact] = iact[k];
1582 vmultc[k] = vmultc[nact];
1583 vmultc[nact] = temp;
1586 /* If stage one is in progress, then set SDIRN to the direction of the next */
1587 /* change to the current vector of variables. */
1595 for (i__ = 1; i__ <= i__1; ++i__) {
1596 temp += sdirn[i__] * a[i__ + kk * a_dim1];
1599 temp /= zdota[nact];
1601 for (i__ = 1; i__ <= i__1; ++i__) {
1602 sdirn[i__] -= temp * z__[i__ + nact * z_dim1];
1606 /* Delete the constraint that has the index IACT(ICON) from the active set. */
1611 vsave = vmultc[icon];
1618 for (i__ = 1; i__ <= i__1; ++i__) {
1619 sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1622 temp = sqrt(sp * sp + d__1 * d__1);
1623 alpha = zdota[kp] / temp;
1625 zdota[kp] = alpha * zdota[k];
1628 for (i__ = 1; i__ <= i__1; ++i__) {
1629 temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1631 z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1632 z__[i__ + kp * z_dim1];
1633 z__[i__ + k * z_dim1] = temp;
1636 vmultc[k] = vmultc[kp];
1646 /* If stage one is in progress, then set SDIRN to the direction of the next */
1647 /* change to the current vector of variables. */
1654 for (i__ = 1; i__ <= i__1; ++i__) {
1655 temp += sdirn[i__] * z__[i__ + (nact + 1) * z_dim1];
1658 for (i__ = 1; i__ <= i__1; ++i__) {
1659 sdirn[i__] -= temp * z__[i__ + (nact + 1) * z_dim1];
1663 /* Pick the next search direction of stage two. */
1666 temp = 1. / zdota[nact];
1668 for (i__ = 1; i__ <= i__1; ++i__) {
1669 sdirn[i__] = temp * z__[i__ + nact * z_dim1];
1672 /* Calculate the step to the boundary of the trust region or take the step */
1673 /* that reduces RESMAX to zero. The two statements below that include the */
1674 /* factor 1.0E-6 prevent some harmless underflows that occurred in a test */
1675 /* calculation. Further, we skip the step if it could be zero within a */
1676 /* reasonable tolerance for computer rounding errors. */
1683 for (i__ = 1; i__ <= i__1; ++i__) {
1684 if ((d__1 = dx[i__], abs(d__1)) >= *rho * 1e-6f) {
1688 sd += dx[i__] * sdirn[i__];
1695 temp = sqrt(ss * dd);
1696 if (abs(sd) >= temp * 1e-6f) {
1697 temp = sqrt(ss * dd + sd * sd);
1699 stpful = dd / (temp + sd);
1702 acca = step + resmax * .1;
1703 accb = step + resmax * .2;
1704 if (step >= acca || acca >= accb) {
1707 step = min(step,resmax);
1710 /* SGJ, 2010: check for error here */
1711 if (nlopt_isinf(step)) return NLOPT_ROUNDOFF_LIMITED;
1713 /* Set DXNEW to the new variables if STEP is the steplength, and reduce */
1714 /* RESMAX to the corresponding maximum residual if stage one is being done. */
1715 /* Because DXNEW will be changed during the calculation of some Lagrange */
1716 /* multipliers, it will be restored to the following value later. */
1719 for (i__ = 1; i__ <= i__1; ++i__) {
1720 dxnew[i__] = dx[i__] + step * sdirn[i__];
1726 for (k = 1; k <= i__1; ++k) {
1730 for (i__ = 1; i__ <= i__2; ++i__) {
1731 temp -= a[i__ + kk * a_dim1] * dxnew[i__];
1733 resmax = max(resmax,temp);
1737 /* Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A */
1738 /* device is included to force VMULTD(K)=0.0 if deviations from this value */
1739 /* can be attributed to computer rounding errors. First calculate the new */
1740 /* Lagrange multipliers. */
1747 for (i__ = 1; i__ <= i__1; ++i__) {
1748 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1750 zdwabs += abs(temp);
1752 acca = zdwabs + abs(zdotw) * .1;
1753 accb = zdwabs + abs(zdotw) * .2;
1754 if (zdwabs >= acca || acca >= accb) {
1757 vmultd[k] = zdotw / zdota[k];
1761 for (i__ = 1; i__ <= i__1; ++i__) {
1762 dxnew[i__] -= vmultd[k] * a[i__ + kk * a_dim1];
1768 d__1 = 0., d__2 = vmultd[nact];
1769 vmultd[nact] = max(d__1,d__2);
1772 /* Complete VMULTC by finding the new constraint residuals. */
1775 for (i__ = 1; i__ <= i__1; ++i__) {
1776 dxnew[i__] = dx[i__] + step * sdirn[i__];
1781 for (k = kl; k <= i__1; ++k) {
1783 sum = resmax - b[kk];
1784 sumabs = resmax + (d__1 = b[kk], abs(d__1));
1786 for (i__ = 1; i__ <= i__2; ++i__) {
1787 temp = a[i__ + kk * a_dim1] * dxnew[i__];
1789 sumabs += abs(temp);
1791 acca = sumabs + abs(sum) * .1f;
1792 accb = sumabs + abs(sum) * .2f;
1793 if (sumabs >= acca || acca >= accb) {
1800 /* Calculate the fraction of the step from DX to DXNEW that will be taken. */
1805 for (k = 1; k <= i__1; ++k) {
1806 if (vmultd[k] < 0.) {
1807 temp = vmultc[k] / (vmultc[k] - vmultd[k]);
1815 /* Update DX, VMULTC and RESMAX. */
1819 for (i__ = 1; i__ <= i__1; ++i__) {
1820 dx[i__] = temp * dx[i__] + ratio * dxnew[i__];
1823 for (k = 1; k <= i__1; ++k) {
1824 d__1 = 0., d__2 = temp * vmultc[k] + ratio * vmultd[k];
1825 vmultc[k] = max(d__1,d__2);
1828 resmax = resold + ratio * (resmax - resold);
1831 /* If the full step is not acceptable then begin another iteration. */
1832 /* Otherwise switch to stage two or end the calculation. */
1837 if (step == stpful) {
1847 /* We employ any freedom that may be available to reduce the objective */
1848 /* function before returning a DX whose length is less than RHO. */
1856 return NLOPT_SUCCESS;