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 /**************************************************************************/
62 /* SGJ, 2008: NLopt-style interface function: */
72 const double *lb, *ub;
77 static int func_wrap(int n, int m, double *x, double *f, double *con,
81 double *xtmp = s->xtmp;
82 const double *lb = s->lb, *ub = s->ub;
84 /* in nlopt, we guarante that the function is never evaluated outside
85 the lb and ub bounds, so we need force this with xtmp ... note
86 that this leads to discontinuity in the first derivative, which
87 slows convergence if we don't enable the ENFORCE_BOUNDS feature
89 for (j = 0; j < n; ++j) {
90 if (x[j] < lb[j]) xtmp[j] = lb[j];
91 else if (x[j] > ub[j]) xtmp[j] = ub[j];
95 *f = s->f(n, xtmp, NULL, s->f_data);
96 if (nlopt_stop_forced(s->stop)) return 1;
98 for (j = 0; j < s->m_orig; ++j) {
99 nlopt_eval_constraint(con + i, NULL, s->fc+j, n, xtmp);
100 if (nlopt_stop_forced(s->stop)) return 1;
101 for (k = 0; k < s->fc[j].m; ++k)
102 con[i + k] = -con[i + k];
105 for (j = 0; j < s->p; ++j) {
106 nlopt_eval_constraint(con + i, NULL, s->h+j, n, xtmp);
107 if (nlopt_stop_forced(s->stop)) return 1;
108 for (k = 0; k < s->h[j].m; ++k)
109 con[(i + s->h[j].m) + k] = -con[i + k];
112 for (j = 0; j < n; ++j) {
113 if (!nlopt_isinf(lb[j]))
114 con[i++] = x[j] - lb[j];
115 if (!nlopt_isinf(ub[j]))
116 con[i++] = ub[j] - x[j];
125 COBYLA_MSG_NONE = 0, /* No messages */
126 COBYLA_MSG_EXIT = 1, /* Exit reasons */
127 COBYLA_MSG_ITER = 2, /* Rho and Sigma changes */
128 COBYLA_MSG_INFO = 3 /* Informational messages */
132 * A function as required by cobyla
133 * state is a void pointer provided to the function at each call
135 * n : the number of variables
136 * m : the number of constraints
137 * x : on input, then vector of variables (should not be modified)
138 * f : on output, the value of the function
139 * con : on output, the value of the constraints (vector of size m)
140 * state : on input, the value of the state variable as provided to cobyla
142 * COBYLA will try to make all the values of the constraints positive.
143 * So if you want to input a constraint j such as x[i] <= MAX, set:
144 * con[j] = MAX - x[i]
145 * The function must returns 0 if no error occurs or 1 to immediately end the
149 typedef int cobyla_function(int n, int m, double *x, double *f, double *con,
150 func_wrap_state *state);
153 * cobyla : minimize a function subject to constraints
155 * n : number of variables (>=0)
156 * m : number of constraints (>=0)
157 * x : on input, initial estimate ; on output, the solution
158 * minf : on output, minimum objective function
159 * rhobeg : a reasonable initial change to the variables
160 * stop : the NLopt stopping criteria
161 * lb, ub : lower and upper bounds on x
162 * message : see the cobyla_message enum
163 * calcfc : the function to minimize (see cobyla_function)
164 * state : used by function (see cobyla_function)
166 * The cobyla function returns the usual nlopt_result codes.
169 extern nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, nlopt_stopping *stop, const double *lb, const double *ub,
170 int message, cobyla_function *calcfc, func_wrap_state *state);
172 nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data,
173 int m, nlopt_constraint *fc,
174 int p, nlopt_constraint *h,
175 const double *lb, const double *ub, /* bounds */
176 double *x, /* in: initial guess, out: minimizer */
178 nlopt_stopping *stop,
185 s.f = f; s.f_data = f_data;
190 s.lb = lb; s.ub = ub;
192 s.xtmp = (double *) malloc(sizeof(double) * n);
193 if (!s.xtmp) return NLOPT_OUT_OF_MEMORY;
195 /* each equality constraint gives two inequality constraints */
196 m = nlopt_count_constraints(m, fc) + 2 * nlopt_count_constraints(p, h);
198 /* add constraints for lower/upper bounds (if any) */
199 for (j = 0; j < n; ++j) {
200 if (!nlopt_isinf(lb[j]))
202 if (!nlopt_isinf(ub[j]))
206 s.con_tol = (double *) malloc(sizeof(double) * m);
207 if (m && !s.con_tol) {
209 return NLOPT_OUT_OF_MEMORY;
211 for (j = 0; j < m; ++j) s.con_tol[j] = 0;
212 for (j = i = 0; i < s.m_orig; ++i) {
213 int j0 = j, jnext = j + fc[i].m;
214 for (; j < jnext; ++j) s.con_tol[j] = fc[i].tol[j - j0];
216 for (i = 0; i < s.p; ++i) {
217 int j0 = j, jnext = j + h[i].m;
218 for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - j0];
219 j0 = j; jnext = j + h[i].m;
220 for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - j0];
223 ret = cobyla(n, m, x, minf, rhobegin, stop, lb, ub, COBYLA_MSG_NONE,
226 /* make sure e.g. rounding errors didn't push us slightly out of bounds */
227 for (j = 0; j < n; ++j) {
228 if (x[j] < lb[j]) x[j] = lb[j];
229 if (x[j] > ub[j]) x[j] = ub[j];
237 /**************************************************************************/
239 /* SGJ, 2010: I find that it seems to increase robustness of the algorithm
240 if, on "simplex" steps (which are intended to improve the linear
241 independence of the simplex, not improve the objective), we multiply
242 the steps by pseudorandom numbers in [0,1). Intuitively, pseudorandom
243 steps are likely to avoid any accidental dependency in the simplex
244 vertices. However, since I still want COBYLA to be a deterministic
245 algorithm, and I don't care all that much about the quality of the
246 randomness here, I implement a simple linear congruential generator
247 rather than use nlopt_urand, and set the initial seed deterministically. */
249 #if defined(HAVE_STDINT_H)
253 #ifndef HAVE_UINT32_T
254 # if SIZEOF_UNSIGNED_LONG == 4
255 typedef unsigned long uint32_t;
256 # elif SIZEOF_UNSIGNED_INT == 4
257 typedef unsigned int uint32_t;
259 # error No 32-bit unsigned integer type
263 /* a simple linear congruential generator */
265 static uint32_t lcg_rand(uint32_t *seed)
267 return (*seed = *seed * 1103515245 + 12345);
270 static double lcg_urand(uint32_t *seed, double a, double b)
272 return a + lcg_rand(seed) * (b - a) / ((uint32_t) -1);
275 /**************************************************************************/
277 static nlopt_result cobylb(int *n, int *m, int *mpp, double *x, double *minf, double *rhobeg,
278 nlopt_stopping *stop, const double *lb, const double *ub, int *iprint, double *con, double *sim,
279 double *simi, double *datmat, double *a, double *vsig, double *veta,
280 double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
281 func_wrap_state *state);
282 static nlopt_result trstlp(int *n, int *m, double *a, double *b, double *rho,
283 double *dx, int *ifull, int *iact, double *z__, double *zdota, double *vmultc,
284 double *sdirn, double *dxnew, double *vmultd);
286 /* ------------------------------------------------------------------------ */
288 nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, nlopt_stopping *stop, const double *lb, const double *ub, int iprint,
289 cobyla_function *calcfc, func_wrap_state *state)
291 int icon, isim, isigb, idatm, iveta, isimi, ivsig, iwork, ia, idx, mpp;
297 * This subroutine minimizes an objective function F(X) subject to M
298 * inequality constraints on X, where X is a vector of variables that has
299 * N components. The algorithm employs linear approximations to the
300 * objective and constraint functions, the approximations being formed by
301 * linear interpolation at N+1 points in the space of the variables.
302 * We regard these interpolation points as vertices of a simplex. The
303 * parameter RHO controls the size of the simplex and it is reduced
304 * automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
305 * to achieve a good vector of variables for the current size, and then
306 * RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and
307 * RHOEND should be set to reasonable initial changes to and the required
308 * accuracy in the variables respectively, but this accuracy should be
309 * viewed as a subject for experimentation because it is not guaranteed.
310 * The subroutine has an advantage over many of its competitors, however,
311 * which is that it treats each constraint individually when calculating
312 * a change to the variables, instead of lumping the constraints together
313 * into a single penalty function. The name of the subroutine is derived
314 * from the phrase Constrained Optimization BY Linear Approximations.
316 * The user must set the values of N, M, RHOBEG and RHOEND, and must
317 * provide an initial vector of variables in X. Further, the value of
318 * IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
319 * printing during the calculation. Specifically, there is no output if
320 * IPRINT=0 and there is output only at the end of the calculation if
321 * IPRINT=1. Otherwise each new value of RHO and SIGMA is printed.
322 * Further, the vector of variables and some function information are
323 * given either when RHO is reduced or when each new value of F(X) is
324 * computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
325 * is a penalty parameter, it being assumed that a change to X is an
326 * improvement if it reduces the merit function
327 * F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)),
328 * where C1,C2,...,CM denote the constraint functions that should become
329 * nonnegative eventually, at least to the precision of RHOEND. In the
330 * printed output the displayed term that is multiplied by SIGMA is
331 * called MAXCV, which stands for 'MAXimum Constraint Violation'. The
332 * argument MAXFUN is an int variable that must be set by the user to a
333 * limit on the number of calls of CALCFC, the purpose of this routine being
334 * given below. The value of MAXFUN will be altered to the number of calls
335 * of CALCFC that are made. The arguments W and IACT provide real and
336 * int arrays that are used as working space. Their lengths must be at
337 * least N*(3*N+2*M+11)+4*M+6 and M+1 respectively.
339 * In order to define the objective and constraint functions, we require
340 * a subroutine that has the name and arguments
341 * SUBROUTINE CALCFC (N,M,X,F,CON)
342 * DIMENSION X(*),CON(*) .
343 * The values of N and M are fixed and have been defined already, while
344 * X is now the current vector of variables. The subroutine should return
345 * the objective and constraint functions at X in F and CON(1),CON(2),
346 * ...,CON(M). Note that we are trying to adjust X so that F(X) is as
347 * small as possible subject to the constraint functions being nonnegative.
349 * Partition the working space array W to provide the storage that is needed
350 * for the main calculation.
357 if (iprint>=1) fprintf(stderr, "cobyla: N==0.\n");
358 return NLOPT_SUCCESS;
363 if (iprint>=1) fprintf(stderr, "cobyla: N<0 or M<0.\n");
364 return NLOPT_INVALID_ARGS;
367 /* workspace allocation */
368 w = (double*) malloc((n*(3*n+2*m+11)+4*m+6)*sizeof(*w));
371 if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
372 return NLOPT_OUT_OF_MEMORY;
374 iact = (int*)malloc((m+1)*sizeof(*iact));
377 if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
379 return NLOPT_OUT_OF_MEMORY;
382 /* Parameter adjustments */
392 isimi = isim + n * n + n;
393 idatm = isimi + n * n;
394 ia = idatm + n * mpp + mpp;
395 ivsig = ia + m * n + n;
400 rc = cobylb(&n, &m, &mpp, &x[1], minf, &rhobeg, stop, &lb[1], &ub[1], &iprint,
401 &w[icon], &w[isim], &w[isimi], &w[idatm], &w[ia], &w[ivsig], &w[iveta],
402 &w[isigb], &w[idx], &w[iwork], &iact[1], calcfc, state);
404 /* Parameter adjustments (reverse) */
414 /* ------------------------------------------------------------------------- */
415 static nlopt_result cobylb(int *n, int *m, int *mpp,
416 double *x, double *minf, double *rhobeg,
417 nlopt_stopping *stop, const double *lb, const double *ub,
418 int *iprint, double *con, double *sim, double *simi,
419 double *datmat, double *a, double *vsig, double *veta,
420 double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
421 func_wrap_state *state)
423 /* System generated locals */
424 int sim_dim1, sim_offset, simi_dim1, simi_offset, datmat_dim1,
425 datmat_offset, a_dim1, a_offset, i__1, i__2, i__3;
428 /* Local variables */
429 double alpha, delta, denom, tempa, barmu;
430 double beta, cmin = 0.0, cmax = 0.0;
431 double cvmaxm, dxsign, prerem = 0.0;
432 double edgmax, pareta, prerec = 0.0, phimin, parsig = 0.0;
434 double phi, rho, sum = 0.0;
435 double ratio, vmold, parmu, error, vmnew;
436 double resmax, cvmaxp;
437 double resnew, trured;
438 double temp, wsig, f;
447 int mp, np, iz, ibrnch;
448 int nbest, ifull, iptem, jdrop;
449 nlopt_result rc = NLOPT_SUCCESS;
451 uint32_t seed = (uint32_t) (*n + *m); /* arbitrary deterministic LCG seed */
454 /* SGJ, 2008: compute rhoend from NLopt stop info */
455 rhoend = stop->xtol_rel * (*rhobeg);
456 for (j = 0; j < *n; ++j)
457 if (rhoend < stop->xtol_abs[j])
458 rhoend = stop->xtol_abs[j];
460 /* SGJ, 2008: added code to keep track of minimum feasible function val */
463 /* Set the initial values of some parameters. The last column of SIM holds */
464 /* the optimal vertex of the current simplex, and the preceding N columns */
465 /* hold the displacements from the optimal vertex to the other vertices. */
466 /* Further, SIMI holds the inverse of the matrix that is contained in the */
467 /* first N columns of SIM. */
469 /* Parameter adjustments */
471 a_offset = 1 + a_dim1 * 1;
474 simi_offset = 1 + simi_dim1 * 1;
477 sim_offset = 1 + sim_dim1 * 1;
480 datmat_offset = 1 + datmat_dim1 * 1;
481 datmat -= datmat_offset;
505 "cobyla: the initial value of RHO is %12.6E and PARMU is set to zero.\n",
510 for (i__ = 1; i__ <= i__1; ++i__) {
512 sim[i__ + np * sim_dim1] = x[i__];
514 for (j = 1; j <= i__2; ++j) {
515 sim[i__ + j * sim_dim1] = 0.;
516 simi[i__ + j * simi_dim1] = 0.;
520 /* SGJ: make sure step rhocur stays inside [lb,ub] */
521 if (x[i__] + rhocur > ub[i__]) {
522 if (x[i__] - rhocur >= lb[i__])
524 else if (ub[i__] - x[i__] > x[i__] - lb[i__])
525 rhocur = 0.5 * (ub[i__] - x[i__]);
527 rhocur = 0.5 * (x[i__] - lb[i__]);
530 sim[i__ + i__ * sim_dim1] = rhocur;
531 simi[i__ + i__ * simi_dim1] = 1.0 / rhocur;
536 /* Make the next call of the user-supplied subroutine CALCFC. These */
537 /* instructions are also used for calling CALCFC during the iterations of */
540 /* SGJ comment: why the hell couldn't he just use a damn subroutine?
541 #*&!%*@ Fortran-66 spaghetti code */
544 if (nlopt_stop_forced(stop)) rc = NLOPT_FORCED_STOP;
545 else if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
546 else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
547 if (rc != NLOPT_SUCCESS) goto L600;
550 if (calcfc(*n, *m, &x[1], &f, &con[1], state))
553 fprintf(stderr, "cobyla: user requested end of minimization.\n");
555 rc = NLOPT_FORCED_STOP;
560 feasible = 1; /* SGJ, 2010 */
563 for (k = 1; k <= i__1; ++k) {
564 d__1 = resmax, d__2 = -con[k];
565 resmax = max(d__1,d__2);
566 if (d__2 > state->con_tol[k-1])
567 feasible = 0; /* SGJ, 2010 */
571 /* SGJ, 2008: check for minf_max reached by feasible point */
572 if (f < stop->minf_max && feasible) {
573 rc = NLOPT_MINF_MAX_REACHED;
574 goto L620; /* not L600 because we want to use current x, f, resmax */
577 if (stop->nevals == *iprint - 1 || *iprint == 3) {
578 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
579 stop->nevals, f, resmax);
581 fprintf(stderr, "cobyla: X =");
582 for (i__ = 1; i__ <= i__1; ++i__) {
583 if (i__>1) fprintf(stderr, " ");
584 fprintf(stderr, "%13.6E", x[i__]);
588 for (i__ = iptemp; i__ <= i__1; ++i__) {
589 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
590 fprintf(stderr, "%15.6E", x[i__]);
593 fprintf(stderr, "\n");
601 /* Set the recently calculated function values in a column of DATMAT. This */
602 /* array has a column for each vertex of the current simplex, the entries of */
603 /* each column being the values of the constraint functions (if any) */
604 /* followed by the objective function and the greatest constraint violation */
608 for (k = 1; k <= i__1; ++k) {
609 datmat[k + jdrop * datmat_dim1] = con[k];
611 if (stop->nevals > np) {
615 /* Exchange the new vertex of the initial simplex with the optimal vertex if */
616 /* necessary. Then, if the initial simplex is not complete, pick its next */
617 /* vertex and calculate the function values there. */
620 if (datmat[mp + np * datmat_dim1] <= f) {
621 x[jdrop] = sim[jdrop + np * sim_dim1];
622 } else { /* improvement in function val */
623 double rhocur = x[jdrop] - sim[jdrop + np * sim_dim1];
624 /* SGJ: use rhocur instead of rho. In original code, rhocur == rho
625 always, but I want to change this to ensure that simplex points
626 fall within [lb,ub]. */
627 sim[jdrop + np * sim_dim1] = x[jdrop];
629 for (k = 1; k <= i__1; ++k) {
630 datmat[k + jdrop * datmat_dim1] = datmat[k + np * datmat_dim1]
632 datmat[k + np * datmat_dim1] = con[k];
635 for (k = 1; k <= i__1; ++k) {
636 sim[jdrop + k * sim_dim1] = -rhocur;
639 for (i__ = k; i__ <= i__2; ++i__) {
640 temp -= simi[i__ + k * simi_dim1];
642 simi[jdrop + k * simi_dim1] = temp;
646 if (stop->nevals <= *n) { /* evaluating initial simplex */
647 jdrop = stop->nevals;
648 /* SGJ: was += rho, but using sim[jdrop,jdrop] enforces consistency
649 if we change the stepsize above to stay in [lb,ub]. */
650 x[jdrop] += sim[jdrop + jdrop * sim_dim1];
656 /* Identify the optimal vertex of the current simplex. */
659 phimin = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
663 for (j = 1; j <= i__1; ++j) {
664 temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
669 } else if (temp == phimin && parmu == 0.) {
670 if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + nbest *
677 /* Switch the best vertex into pole position if it is not there already, */
678 /* and also update SIM, SIMI and DATMAT. */
682 for (i__ = 1; i__ <= i__1; ++i__) {
683 temp = datmat[i__ + np * datmat_dim1];
684 datmat[i__ + np * datmat_dim1] = datmat[i__ + nbest * datmat_dim1]
686 datmat[i__ + nbest * datmat_dim1] = temp;
689 for (i__ = 1; i__ <= i__1; ++i__) {
690 temp = sim[i__ + nbest * sim_dim1];
691 sim[i__ + nbest * sim_dim1] = 0.;
692 sim[i__ + np * sim_dim1] += temp;
695 for (k = 1; k <= i__2; ++k) {
696 sim[i__ + k * sim_dim1] -= temp;
697 tempa -= simi[k + i__ * simi_dim1];
699 simi[nbest + i__ * simi_dim1] = tempa;
703 /* Make an error return if SIGI is a poor approximation to the inverse of */
704 /* the leading N by N submatrix of SIG. */
708 for (i__ = 1; i__ <= i__1; ++i__) {
710 for (j = 1; j <= i__2; ++j) {
716 for (k = 1; k <= i__3; ++k) if (sim[k + j * sim_dim1] != 0) {
717 temp += simi[i__ + k * simi_dim1] * sim[k + j * sim_dim1];
719 d__1 = error, d__2 = abs(temp);
720 error = max(d__1,d__2);
725 fprintf(stderr, "cobyla: rounding errors are becoming damaging.\n");
731 /* Calculate the coefficients of the linear approximations to the objective */
732 /* and constraint functions, placing minus the objective function gradient */
733 /* after the constraint gradients in the array A. The vector W is used for */
737 for (k = 1; k <= i__2; ++k) {
738 con[k] = -datmat[k + np * datmat_dim1];
740 for (j = 1; j <= i__1; ++j) {
741 w[j] = datmat[k + j * datmat_dim1] + con[k];
744 for (i__ = 1; i__ <= i__1; ++i__) {
747 for (j = 1; j <= i__3; ++j) {
748 temp += w[j] * simi[j + i__ * simi_dim1];
753 a[i__ + k * a_dim1] = temp;
757 /* Calculate the values of sigma and eta, and set IFLAG=0 if the current */
758 /* simplex is not acceptable. */
761 parsig = alpha * rho;
764 for (j = 1; j <= i__1; ++j) {
768 for (i__ = 1; i__ <= i__2; ++i__) {
769 d__1 = simi[j + i__ * simi_dim1];
771 d__1 = sim[i__ + j * sim_dim1];
774 vsig[j] = 1. / sqrt(wsig);
775 veta[j] = sqrt(weta);
776 if (vsig[j] < parsig || veta[j] > pareta) {
781 /* If a new vertex is needed to improve acceptability, then decide which */
782 /* vertex to drop from the simplex. */
784 if (ibrnch == 1 || iflag == 1) {
790 for (j = 1; j <= i__1; ++j) {
791 if (veta[j] > temp) {
798 for (j = 1; j <= i__1; ++j) {
799 if (vsig[j] < temp) {
806 /* Calculate the step to the new vertex and its sign. */
808 temp = gamma_ * rho * vsig[jdrop];
810 for (i__ = 1; i__ <= i__1; ++i__) {
811 dx[i__] = temp * simi[jdrop + i__ * simi_dim1];
816 for (k = 1; k <= i__1; ++k) {
819 for (i__ = 1; i__ <= i__2; ++i__) {
820 sum += a[i__ + k * a_dim1] * dx[i__];
823 temp = datmat[k + np * datmat_dim1];
824 d__1 = cvmaxp, d__2 = -sum - temp;
825 cvmaxp = max(d__1,d__2);
826 d__1 = cvmaxm, d__2 = sum - temp;
827 cvmaxm = max(d__1,d__2);
831 if (parmu * (cvmaxp - cvmaxm) > sum + sum) {
835 /* Update the elements of SIM and SIMI, and set the next X. */
839 for (i__ = 1; i__ <= i__1; ++i__) {
840 /* SGJ, 2010: pseudo-randomize simplex steps (see LCG comments above) */
841 dx[i__] = dxsign * dx[i__] * lcg_urand(&seed, 0.01, 1);
842 /* SGJ: make sure dx step says in [lb,ub] */
845 double xi = sim[i__ + np * sim_dim1];
847 if (xi + dx[i__] > ub[i__])
849 if (xi + dx[i__] < lb[i__]) {
850 if (xi - dx[i__] <= ub[i__])
852 else { /* try again with halved step */
859 sim[i__ + jdrop * sim_dim1] = dx[i__];
860 temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
863 for (i__ = 1; i__ <= i__1; ++i__) {
864 simi[jdrop + i__ * simi_dim1] /= temp;
867 for (j = 1; j <= i__1; ++j) {
871 for (i__ = 1; i__ <= i__2; ++i__) {
872 temp += simi[j + i__ * simi_dim1] * dx[i__];
875 for (i__ = 1; i__ <= i__2; ++i__) {
876 simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
880 x[j] = sim[j + np * sim_dim1] + dx[j];
884 /* Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO. */
888 izdota = iz + *n * *n;
891 idxnew = isdirn + *n;
893 rc = trstlp(n, m, &a[a_offset], &con[1], &rho, &dx[1], &ifull, &iact[1], &w[
894 iz], &w[izdota], &w[ivmc], &w[isdirn], &w[idxnew], &w[ivmd]);
895 if (rc != NLOPT_SUCCESS) goto L600;
897 /* SGJ: since the bound constraints are linear, we should never get
898 a dx that lies outside the [lb,ub] constraints here, but we'll
899 enforce this anyway just to be paranoid */
901 for (i__ = 1; i__ <= i__1; ++i__) {
902 double xi = sim[i__ + np * sim_dim1];
903 if (xi + dx[i__] > ub[i__]) dx[i__] = ub[i__] - xi;
904 if (xi + dx[i__] < lb[i__]) dx[i__] = xi - lb[i__];
910 for (i__ = 1; i__ <= i__1; ++i__) {
914 if (temp < rho * .25 * rho) {
920 /* Predict the change to F and the new maximum constraint violation if the */
921 /* variables are altered from x(0) to x(0)+DX. */
926 for (k = 1; k <= i__1; ++k) {
929 for (i__ = 1; i__ <= i__2; ++i__) {
930 sum -= a[i__ + k * a_dim1] * dx[i__];
933 resnew = max(resnew,sum);
937 /* Increase PARMU if necessary and branch back if this change alters the */
938 /* optimal vertex. Otherwise PREREM and PREREC will be set to the predicted */
939 /* reductions in the merit function and the maximum constraint violation */
943 prerec = datmat[*mpp + np * datmat_dim1] - resnew;
945 barmu = sum / prerec;
947 if (parmu < barmu * 1.5) {
950 fprintf(stderr, "cobyla: increase in PARMU to %12.6E\n", parmu);
952 phi = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
955 for (j = 1; j <= i__1; ++j) {
956 temp = datmat[mp + j * datmat_dim1] + parmu * datmat[*mpp + j *
961 if (temp == phi && parmu == 0.f) {
962 if (datmat[*mpp + j * datmat_dim1] < datmat[*mpp + np *
969 prerem = parmu * prerec - sum;
971 /* Calculate the constraint and objective functions at x(*). Then find the */
972 /* actual reduction in the merit function. */
975 for (i__ = 1; i__ <= i__1; ++i__) {
976 x[i__] = sim[i__ + np * sim_dim1] + dx[i__];
981 vmold = datmat[mp + np * datmat_dim1] + parmu * datmat[*mpp + np *
983 vmnew = f + parmu * resmax;
984 trured = vmold - vmnew;
985 if (parmu == 0. && f == datmat[mp + np * datmat_dim1]) {
987 trured = datmat[*mpp + np * datmat_dim1] - resmax;
990 /* Begin the operations that decide whether x(*) should replace one of the */
991 /* vertices of the current simplex, the change being mandatory if TRURED is */
992 /* positive. Firstly, JDROP is set to the index of the vertex that is to be */
1001 for (j = 1; j <= i__1; ++j) {
1004 for (i__ = 1; i__ <= i__2; ++i__) {
1005 temp += simi[j + i__ * simi_dim1] * dx[i__];
1012 sigbar[j] = temp * vsig[j];
1015 /* Calculate the value of ell. */
1017 edgmax = delta * rho;
1020 for (j = 1; j <= i__1; ++j) {
1021 if (sigbar[j] >= parsig || sigbar[j] >= vsig[j]) {
1026 for (i__ = 1; i__ <= i__2; ++i__) {
1027 d__1 = dx[i__] - sim[i__ + j * sim_dim1];
1028 temp += d__1 * d__1;
1032 if (temp > edgmax) {
1045 /* Revise the simplex by updating the elements of SIM, SIMI and DATMAT. */
1049 for (i__ = 1; i__ <= i__1; ++i__) {
1050 sim[i__ + jdrop * sim_dim1] = dx[i__];
1051 temp += simi[jdrop + i__ * simi_dim1] * dx[i__];
1054 for (i__ = 1; i__ <= i__1; ++i__) {
1055 simi[jdrop + i__ * simi_dim1] /= temp;
1058 for (j = 1; j <= i__1; ++j) {
1062 for (i__ = 1; i__ <= i__2; ++i__) {
1063 temp += simi[j + i__ * simi_dim1] * dx[i__];
1066 for (i__ = 1; i__ <= i__2; ++i__) {
1067 simi[j + i__ * simi_dim1] -= temp * simi[jdrop + i__ *
1073 for (k = 1; k <= i__1; ++k) {
1074 datmat[k + jdrop * datmat_dim1] = con[k];
1077 /* Branch back for further iterations with the current RHO. */
1079 if (trured > 0. && trured >= prerem * .1) {
1080 /* SGJ, 2010: following a suggestion in the SAS manual (which
1081 mentions a similar modification to COBYLA, although they didn't
1082 publish their source code), increase rho if predicted reduction
1083 is sufficiently close to the actual reduction. Otherwise,
1084 COBLYA seems to easily get stuck making very small steps.
1085 Also require iflag != 0 (i.e., acceptable simplex), again
1086 following SAS suggestion (otherwise I get convergence failure
1088 if (trured >= prerem * 0.9 && trured <= prerem * 1.1 && iflag) {
1099 /* SGJ, 2008: convergence tests for function vals; note that current
1100 best val is stored in datmat[mp + np * datmat_dim1], or in f if
1101 ifull == 1, and previous best is in *minf. This seems like a
1102 sensible place to put the convergence test, as it is the same
1103 place that Powell checks the x tolerance (rho > rhoend). */
1105 double fbest = ifull == 1 ? f : datmat[mp + np * datmat_dim1];
1106 if (fbest < *minf && nlopt_stop_ftol(stop, fbest, *minf)) {
1107 rc = NLOPT_FTOL_REACHED;
1113 /* Otherwise reduce RHO if it is not at its least value and reset PARMU. */
1117 if (rho <= rhoend * 1.5) {
1123 for (k = 1; k <= i__1; ++k) {
1124 cmin = datmat[k + np * datmat_dim1];
1127 for (i__ = 1; i__ <= i__2; ++i__) {
1128 d__1 = cmin, d__2 = datmat[k + i__ * datmat_dim1];
1129 cmin = min(d__1,d__2);
1130 d__1 = cmax, d__2 = datmat[k + i__ * datmat_dim1];
1131 cmax = max(d__1,d__2);
1133 if (k <= *m && cmin < cmax * .5) {
1134 temp = max(cmax,0.) - cmin;
1138 denom = min(denom,temp);
1144 } else if (cmax - cmin < parmu * denom) {
1145 parmu = (cmax - cmin) / denom;
1149 fprintf(stderr, "cobyla: reduction in RHO to %12.6E and PARMU =%13.6E\n",
1153 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1154 stop->nevals, datmat[mp + np * datmat_dim1], datmat[*mpp + np * datmat_dim1]);
1156 fprintf(stderr, "cobyla: X =");
1158 for (i__ = 1; i__ <= i__1; ++i__) {
1159 if (i__>1) fprintf(stderr, " ");
1160 fprintf(stderr, "%13.6E", sim[i__ + np * sim_dim1]);
1164 for (i__ = iptemp; i__ <= i__1; ++i__) {
1165 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
1166 fprintf(stderr, "%15.6E", x[i__]);
1169 fprintf(stderr, "\n");
1174 rc = NLOPT_XTOL_REACHED;
1176 /* Return the best calculated values of the variables. */
1179 fprintf(stderr, "cobyla: normal return.\n");
1186 for (i__ = 1; i__ <= i__1; ++i__) {
1187 x[i__] = sim[i__ + np * sim_dim1];
1189 f = datmat[mp + np * datmat_dim1];
1190 resmax = datmat[*mpp + np * datmat_dim1];
1194 fprintf(stderr, "cobyla: NFVALS = %4d, F =%13.6E, MAXCV =%13.6E\n",
1195 stop->nevals, f, resmax);
1197 fprintf(stderr, "cobyla: X =");
1198 for (i__ = 1; i__ <= i__1; ++i__) {
1199 if (i__>1) fprintf(stderr, " ");
1200 fprintf(stderr, "%13.6E", x[i__]);
1204 for (i__ = iptemp; i__ <= i__1; ++i__) {
1205 if (!((i__-1) % 4)) fprintf(stderr, "\ncobyla: ");
1206 fprintf(stderr, "%15.6E", x[i__]);
1209 fprintf(stderr, "\n");
1214 /* ------------------------------------------------------------------------- */
1215 static nlopt_result trstlp(int *n, int *m, double *a,
1216 double *b, double *rho, double *dx, int *ifull,
1217 int *iact, double *z__, double *zdota, double *vmultc,
1218 double *sdirn, double *dxnew, double *vmultd)
1220 /* System generated locals */
1221 int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
1224 /* Local variables */
1225 double alpha, tempa;
1227 double optnew, stpful, sum, tot, acca, accb;
1228 double ratio, vsave, zdotv, zdotw, dd;
1230 double sp, ss, resold = 0.0, zdvabs, zdwabs, sumabs, resmax, optold;
1234 int iout, i__, j, k;
1238 int nact, icon = 0, mcon;
1242 /* This subroutine calculates an N-component vector DX by applying the */
1243 /* following two stages. In the first stage, DX is set to the shortest */
1244 /* vector that minimizes the greatest violation of the constraints */
1245 /* A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K=2,3,...,M, */
1246 /* subject to the Euclidean length of DX being at most RHO. If its length is */
1247 /* strictly less than RHO, then we use the resultant freedom in DX to */
1248 /* minimize the objective function */
1249 /* -A(1,M+1)*DX(1)-A(2,M+1)*DX(2)-...-A(N,M+1)*DX(N) */
1250 /* subject to no increase in any greatest constraint violation. This */
1251 /* notation allows the gradient of the objective function to be regarded as */
1252 /* the gradient of a constraint. Therefore the two stages are distinguished */
1253 /* by MCON .EQ. M and MCON .GT. M respectively. It is possible that a */
1254 /* degeneracy may prevent DX from attaining the target length RHO. Then the */
1255 /* value IFULL=0 would be set, but usually IFULL=1 on return. */
1257 /* In general NACT is the number of constraints in the active set and */
1258 /* IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT */
1259 /* contains a permutation of the remaining constraint indices. Further, Z is */
1260 /* an orthogonal matrix whose first NACT columns can be regarded as the */
1261 /* result of Gram-Schmidt applied to the active constraint gradients. For */
1262 /* J=1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th */
1263 /* column of Z with the gradient of the J-th active constraint. DX is the */
1264 /* current vector of variables and here the residuals of the active */
1265 /* constraints should be zero. Further, the active constraints have */
1266 /* nonnegative Lagrange multipliers that are held at the beginning of */
1267 /* VMULTC. The remainder of this vector holds the residuals of the inactive */
1268 /* constraints at DX, the ordering of the components of VMULTC being in */
1269 /* agreement with the permutation of the indices of the constraints that is */
1270 /* in IACT. All these residuals are nonnegative, which is achieved by the */
1271 /* shift RESMAX that makes the least residual zero. */
1273 /* Initialize Z and some other variables. The value of RESMAX will be */
1274 /* appropriate to DX=0, while ICON will be the index of a most violated */
1275 /* constraint if RESMAX is positive. Usually during the first stage the */
1276 /* vector SDIRN gives a search direction that reduces all the active */
1277 /* constraint violations by one simultaneously. */
1279 /* Parameter adjustments */
1281 z_offset = 1 + z_dim1 * 1;
1284 a_offset = 1 + a_dim1 * 1;
1301 for (i__ = 1; i__ <= i__1; ++i__) {
1303 for (j = 1; j <= i__2; ++j) {
1304 z__[i__ + j * z_dim1] = 0.;
1306 z__[i__ + i__ * z_dim1] = 1.;
1311 for (k = 1; k <= i__1; ++k) {
1312 if (b[k] > resmax) {
1318 for (k = 1; k <= i__1; ++k) {
1320 vmultc[k] = resmax - b[k];
1327 for (i__ = 1; i__ <= i__1; ++i__) {
1331 /* End the current stage of the calculation if 3 consecutive iterations */
1332 /* have either failed to reduce the best calculated value of the objective */
1333 /* function or to increase the number of active constraints since the best */
1334 /* value was calculated. This strategy prevents cycling, but there is a */
1335 /* remote possibility that it will cause premature termination. */
1346 for (i__ = 1; i__ <= i__1; ++i__) {
1347 optnew -= dx[i__] * a[i__ + mcon * a_dim1];
1350 if (icount == 0 || optnew < optold) {
1354 } else if (nact > nactx) {
1364 /* If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to */
1365 /* the active set. Apply Givens rotations so that the last N-NACT-1 columns */
1366 /* of Z are orthogonal to the gradient of the new constraint, a scalar */
1367 /* product being set to zero if its nonzero value could be due to computer */
1368 /* rounding errors. The array DXNEW is used for working space. */
1375 for (i__ = 1; i__ <= i__1; ++i__) {
1376 dxnew[i__] = a[i__ + kk * a_dim1];
1385 for (i__ = 1; i__ <= i__1; ++i__) {
1386 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1390 acca = spabs + abs(sp) * .1;
1391 accb = spabs + abs(sp) * .2;
1392 if (spabs >= acca || acca >= accb) {
1399 temp = sqrt(sp * sp + tot * tot);
1404 for (i__ = 1; i__ <= i__1; ++i__) {
1405 temp = alpha * z__[i__ + k * z_dim1] + beta * z__[i__ + kp *
1407 z__[i__ + kp * z_dim1] = alpha * z__[i__ + kp * z_dim1] -
1408 beta * z__[i__ + k * z_dim1];
1409 z__[i__ + k * z_dim1] = temp;
1416 /* Add the new constraint if this can be done without a deletion from the */
1422 vmultc[icon] = vmultc[nact];
1427 /* The next instruction is reached if a deletion has to be made from the */
1428 /* active set in order to make room for the new active constraint, because */
1429 /* the new constraint gradient is a linear combination of the gradients of */
1430 /* the old active constraints. Set the elements of VMULTD to the multipliers */
1431 /* of the linear combination. Further, set IOUT to the index of the */
1432 /* constraint to be deleted, but branch if no suitable index can be found. */
1440 for (i__ = 1; i__ <= i__1; ++i__) {
1441 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1443 zdvabs += abs(temp);
1445 acca = zdvabs + abs(zdotv) * .1;
1446 accb = zdvabs + abs(zdotv) * .2;
1447 if (zdvabs < acca && acca < accb) {
1448 temp = zdotv / zdota[k];
1449 if (temp > 0. && iact[k] <= *m) {
1450 tempa = vmultc[k] / temp;
1451 if (ratio < 0. || tempa < ratio) {
1459 for (i__ = 1; i__ <= i__1; ++i__) {
1460 dxnew[i__] -= temp * a[i__ + kw * a_dim1];
1475 /* Revise the Lagrange multipliers and reorder the active constraints so */
1476 /* that the one to be replaced is at the end of the list. Also calculate the */
1477 /* new value of ZDOTA(NACT) and branch if it is not acceptable. */
1480 for (k = 1; k <= i__1; ++k) {
1481 d__1 = 0., d__2 = vmultc[k] - ratio * vmultd[k];
1482 vmultc[k] = max(d__1,d__2);
1486 vsave = vmultc[icon];
1493 for (i__ = 1; i__ <= i__1; ++i__) {
1494 sp += z__[i__ + k * z_dim1] * a[i__ + kw * a_dim1];
1497 temp = sqrt(sp * sp + d__1 * d__1);
1498 alpha = zdota[kp] / temp;
1500 zdota[kp] = alpha * zdota[k];
1503 for (i__ = 1; i__ <= i__1; ++i__) {
1504 temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1506 z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1507 z__[i__ + kp * z_dim1];
1508 z__[i__ + k * z_dim1] = temp;
1511 vmultc[k] = vmultc[kp];
1521 for (i__ = 1; i__ <= i__1; ++i__) {
1522 temp += z__[i__ + nact * z_dim1] * a[i__ + kk * a_dim1];
1529 vmultc[nact] = ratio;
1531 /* Update IACT and ensure that the objective function continues to be */
1532 /* treated as the last active constraint when MCON>M. */
1535 iact[icon] = iact[nact];
1537 if (mcon > *m && kk != mcon) {
1541 for (i__ = 1; i__ <= i__1; ++i__) {
1542 sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1545 temp = sqrt(sp * sp + d__1 * d__1);
1546 alpha = zdota[nact] / temp;
1548 zdota[nact] = alpha * zdota[k];
1551 for (i__ = 1; i__ <= i__1; ++i__) {
1552 temp = alpha * z__[i__ + nact * z_dim1] + beta * z__[i__ + k *
1554 z__[i__ + nact * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1555 z__[i__ + nact * z_dim1];
1556 z__[i__ + k * z_dim1] = temp;
1558 iact[nact] = iact[k];
1561 vmultc[k] = vmultc[nact];
1562 vmultc[nact] = temp;
1565 /* If stage one is in progress, then set SDIRN to the direction of the next */
1566 /* change to the current vector of variables. */
1574 for (i__ = 1; i__ <= i__1; ++i__) {
1575 temp += sdirn[i__] * a[i__ + kk * a_dim1];
1578 temp /= zdota[nact];
1580 for (i__ = 1; i__ <= i__1; ++i__) {
1581 sdirn[i__] -= temp * z__[i__ + nact * z_dim1];
1585 /* Delete the constraint that has the index IACT(ICON) from the active set. */
1590 vsave = vmultc[icon];
1597 for (i__ = 1; i__ <= i__1; ++i__) {
1598 sp += z__[i__ + k * z_dim1] * a[i__ + kk * a_dim1];
1601 temp = sqrt(sp * sp + d__1 * d__1);
1602 alpha = zdota[kp] / temp;
1604 zdota[kp] = alpha * zdota[k];
1607 for (i__ = 1; i__ <= i__1; ++i__) {
1608 temp = alpha * z__[i__ + kp * z_dim1] + beta * z__[i__ + k *
1610 z__[i__ + kp * z_dim1] = alpha * z__[i__ + k * z_dim1] - beta *
1611 z__[i__ + kp * z_dim1];
1612 z__[i__ + k * z_dim1] = temp;
1615 vmultc[k] = vmultc[kp];
1625 /* If stage one is in progress, then set SDIRN to the direction of the next */
1626 /* change to the current vector of variables. */
1633 for (i__ = 1; i__ <= i__1; ++i__) {
1634 temp += sdirn[i__] * z__[i__ + (nact + 1) * z_dim1];
1637 for (i__ = 1; i__ <= i__1; ++i__) {
1638 sdirn[i__] -= temp * z__[i__ + (nact + 1) * z_dim1];
1642 /* Pick the next search direction of stage two. */
1645 temp = 1. / zdota[nact];
1647 for (i__ = 1; i__ <= i__1; ++i__) {
1648 sdirn[i__] = temp * z__[i__ + nact * z_dim1];
1651 /* Calculate the step to the boundary of the trust region or take the step */
1652 /* that reduces RESMAX to zero. The two statements below that include the */
1653 /* factor 1.0E-6 prevent some harmless underflows that occurred in a test */
1654 /* calculation. Further, we skip the step if it could be zero within a */
1655 /* reasonable tolerance for computer rounding errors. */
1662 for (i__ = 1; i__ <= i__1; ++i__) {
1663 if ((d__1 = dx[i__], abs(d__1)) >= *rho * 1e-6f) {
1667 sd += dx[i__] * sdirn[i__];
1674 temp = sqrt(ss * dd);
1675 if (abs(sd) >= temp * 1e-6f) {
1676 temp = sqrt(ss * dd + sd * sd);
1678 stpful = dd / (temp + sd);
1681 acca = step + resmax * .1;
1682 accb = step + resmax * .2;
1683 if (step >= acca || acca >= accb) {
1686 step = min(step,resmax);
1689 /* SGJ, 2010: check for error here */
1690 if (nlopt_isinf(step)) return NLOPT_ROUNDOFF_LIMITED;
1692 /* Set DXNEW to the new variables if STEP is the steplength, and reduce */
1693 /* RESMAX to the corresponding maximum residual if stage one is being done. */
1694 /* Because DXNEW will be changed during the calculation of some Lagrange */
1695 /* multipliers, it will be restored to the following value later. */
1698 for (i__ = 1; i__ <= i__1; ++i__) {
1699 dxnew[i__] = dx[i__] + step * sdirn[i__];
1705 for (k = 1; k <= i__1; ++k) {
1709 for (i__ = 1; i__ <= i__2; ++i__) {
1710 temp -= a[i__ + kk * a_dim1] * dxnew[i__];
1712 resmax = max(resmax,temp);
1716 /* Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A */
1717 /* device is included to force VMULTD(K)=0.0 if deviations from this value */
1718 /* can be attributed to computer rounding errors. First calculate the new */
1719 /* Lagrange multipliers. */
1726 for (i__ = 1; i__ <= i__1; ++i__) {
1727 temp = z__[i__ + k * z_dim1] * dxnew[i__];
1729 zdwabs += abs(temp);
1731 acca = zdwabs + abs(zdotw) * .1;
1732 accb = zdwabs + abs(zdotw) * .2;
1733 if (zdwabs >= acca || acca >= accb) {
1736 vmultd[k] = zdotw / zdota[k];
1740 for (i__ = 1; i__ <= i__1; ++i__) {
1741 dxnew[i__] -= vmultd[k] * a[i__ + kk * a_dim1];
1747 d__1 = 0., d__2 = vmultd[nact];
1748 vmultd[nact] = max(d__1,d__2);
1751 /* Complete VMULTC by finding the new constraint residuals. */
1754 for (i__ = 1; i__ <= i__1; ++i__) {
1755 dxnew[i__] = dx[i__] + step * sdirn[i__];
1760 for (k = kl; k <= i__1; ++k) {
1762 sum = resmax - b[kk];
1763 sumabs = resmax + (d__1 = b[kk], abs(d__1));
1765 for (i__ = 1; i__ <= i__2; ++i__) {
1766 temp = a[i__ + kk * a_dim1] * dxnew[i__];
1768 sumabs += abs(temp);
1770 acca = sumabs + abs(sum) * .1f;
1771 accb = sumabs + abs(sum) * .2f;
1772 if (sumabs >= acca || acca >= accb) {
1779 /* Calculate the fraction of the step from DX to DXNEW that will be taken. */
1784 for (k = 1; k <= i__1; ++k) {
1785 if (vmultd[k] < 0.) {
1786 temp = vmultc[k] / (vmultc[k] - vmultd[k]);
1794 /* Update DX, VMULTC and RESMAX. */
1798 for (i__ = 1; i__ <= i__1; ++i__) {
1799 dx[i__] = temp * dx[i__] + ratio * dxnew[i__];
1802 for (k = 1; k <= i__1; ++k) {
1803 d__1 = 0., d__2 = temp * vmultc[k] + ratio * vmultd[k];
1804 vmultc[k] = max(d__1,d__2);
1807 resmax = resold + ratio * (resmax - resold);
1810 /* If the full step is not acceptable then begin another iteration. */
1811 /* Otherwise switch to stage two or end the calculation. */
1816 if (step == stpful) {
1826 /* We employ any freedom that may be available to reduce the objective */
1827 /* function before returning a DX whose length is less than RHO. */
1835 return NLOPT_SUCCESS;