1 /* Copyright (c) 2004 M. J. D. Powell (mjdp@cam.ac.uk)
2 * Copyright (c) 2007-2008 Massachusetts Institute of Technology
4 * Permission is hereby granted, free of charge, to any person obtaining
5 * a copy of this software and associated documentation files (the
6 * "Software"), to deal in the Software without restriction, including
7 * without limitation the rights to use, copy, modify, merge, publish,
8 * distribute, sublicense, and/or sell copies of the Software, and to
9 * permit persons to whom the Software is furnished to do so, subject to
10 * the following conditions:
12 * The above copyright notice and this permission notice shall be
13 * included in all copies or substantial portions of the Software.
15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
16 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
18 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
19 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
20 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
21 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
24 /* NEWUOA derivative-free optimization algorithm by M. J. D. Powell.
25 Original Fortran code by Powell (2004). Converted via f2c, cleaned up,
26 and incorporated into NLopt by S. G. Johnson (2008). See README. */
35 #define min(a,b) ((a) <= (b) ? (a) : (b))
36 #define max(a,b) ((a) >= (b) ? (a) : (b))
38 /*************************************************************************/
43 double *xpt, *pq, *hq, *gq, *xopt;
48 static double quad_model(int n, const double *x, double *grad, void *data)
50 quad_model_data *d = (quad_model_data *) data;
51 const double *xpt = d->xpt, *pq = d->pq, *hq = d->hq, *gq = d->gq, *xopt = d->xopt;
57 /* first, set hd to be Hessian matrix times x */
58 memset(hd, 0, sizeof(double) * n);
59 /* implicit Hessian terms (a sum of outer products of xpt vectors) */
60 for (k = 0; k < npt; ++k) {
62 for (j = 0; j < n; ++j)
63 temp += xpt[k + j * npt] * (xopt[j] + x[j]);
65 for (i = 0; i < n; ++i)
66 hd[i] += temp * xpt[k + i * npt];
68 /* explicit Hessian terms (stored as compressed lower triangle hq) */
70 for (j = 0; j < n; ++j) {
71 for (i = 0; i < j; ++i) {
72 hd[j] += hq[k] * (xopt[i] + x[i]);
73 hd[i] += hq[k] * (xopt[j] + x[j]);
76 hd[j] += hq[k++] * (xopt[j] + x[j]);
79 for (i = 0; i < n; ++i) {
80 val += (gq[i] + 0.5 * hd[i]) * (xopt[i] + x[i]);
81 if (grad) grad[i] = gq[i] + hd[i];
87 /* constraint function to enforce |x|^2 <= rho*rho */
88 static double rho_constraint(int n, const double *x, double *grad, void *data)
90 double rho = *((double *) data);
91 double val = - rho*rho;
93 for (i = 0; i < n; ++i)
95 if (grad) for (i = 0; i < n; ++i) grad[i] = 2*x[i];
99 static nlopt_result trsapp_(int *n, int *npt, double *xopt,
100 double *xpt, double *gq, double *hq, double *pq,
101 double *delta, double *step, double *d__, double *g,
102 double *hd, double *hs, double *crvmin,
103 const double *xbase, const double *lb, const double *ub)
105 /* System generated locals */
106 int xpt_dim1, xpt_offset, i__1, i__2;
109 /* Local variables */
111 double dd, cf, dg, gg;
115 double ss, dhd, dhs, cth, sgk, shs, sth, qadd, half, qbeg, qred, qmin,
116 temp, qsav, qnew, zero, ggbeg, alpha, angle, reduc;
118 double ggsav, delsq, tempa, tempb;
120 double bstep, ratio, twopi;
126 /* N is the number of variables of a quadratic objective function, Q say. */
127 /* The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meanings, */
128 /* in order to define the current quadratic model Q. */
129 /* DELTA is the trust region radius, and has to be positive. */
130 /* STEP will be set to the calculated trial step. */
131 /* The arrays D, G, HD and HS will be used for working space. */
132 /* CRVMIN will be set to the least curvature of H along the conjugate */
133 /* directions that occur, except that it is set to zero if STEP goes */
134 /* all the way to the trust region boundary. */
136 /* The calculation of STEP begins with the truncated conjugate gradient */
137 /* method. If the boundary of the trust region is reached, then further */
138 /* changes to STEP may be made, each one being in the 2D space spanned */
139 /* by the current STEP and the corresponding gradient of Q. Thus STEP */
140 /* should provide a substantial reduction to Q within the trust region. */
142 /* Initialization, which includes setting HD to H times XOPT. */
144 /* Parameter adjustments */
146 xpt_offset = 1 + xpt_dim1;
159 double *slb, *sub, *xtol, minf, crv;
163 qmd.xpt = &xpt[1 + 1 * xpt_dim1];
173 for (j = 0; j < *n; ++j) {
174 slb[j] = -(sub[j] = *delta);
175 if (slb[j] < lb[j] - xbase[j] - xopt[j+1])
176 slb[j] = lb[j] - xbase[j] - xopt[j+1];
177 if (sub[j] > ub[j] - xbase[j] - xopt[j+1])
178 sub[j] = ub[j] - xbase[j] - xopt[j+1];
179 if (slb[j] > 0) slb[j] = 0;
180 if (sub[j] < 0) sub[j] = 0;
181 xtol[j] = 1e-7 * *delta; /* absolute x tolerance */
183 memset(&step[1], 0, sizeof(double) * *n);
184 ret = nlopt_minimize_constrained(NLOPT_LD_MMA, *n, quad_model, &qmd,
185 1, rho_constraint, delta, sizeof(double),
186 slb, sub, &step[1], &minf, -HUGE_VAL,
187 0., 0., 0., xtol, 1000, 0.);
188 if (rho_constraint(*n, &step[1], 0, delta) > -1e-6*(*delta)*(*delta))
191 for (j = 1; j <= *n; ++j) d__[j] = step[j] - xopt[j];
192 quad_model(*n, &d__[1], &g[1], &qmd);
194 for (j = 1; j <= *n; ++j) {
195 crv += step[j] * (g[j] - gq[j]);
196 gg += step[j] * step[j];
198 if (gg <= 1e-16 * crv)
210 twopi = atan(1.) * 8.;
211 delsq = *delta * *delta;
216 for (i__ = 1; i__ <= i__1; ++i__) {
218 d__[i__] = xopt[i__];
222 /* Prepare for the first line search. */
228 for (i__ = 1; i__ <= i__1; ++i__) {
231 g[i__] = gq[i__] + hd[i__];
234 /* Computing 2nd power */
247 /* Calculate the step to the trust region boundary and the product HD. */
252 bstep = temp / (ds + sqrt(ds * ds + dd * temp));
257 for (j = 1; j <= i__1; ++j) {
259 dhd += d__[j] * hd[j];
262 /* Update CRVMIN and set the step-length ALPHA. */
270 *crvmin = min(*crvmin,temp);
272 d__1 = alpha, d__2 = gg / dhd;
273 alpha = min(d__1,d__2);
275 qadd = alpha * (gg - half * alpha * dhd);
278 /* Update STEP and HS. */
283 for (i__ = 1; i__ <= i__1; ++i__) {
284 step[i__] += alpha * d__[i__];
285 hs[i__] += alpha * hd[i__];
287 /* Computing 2nd power */
288 d__1 = g[i__] + hs[i__];
292 /* Begin another conjugate direction iteration if required. */
295 if (qadd <= qred * .01) {
298 if (gg <= ggbeg * 1e-4) {
301 if (iterc == itermax) {
309 for (i__ = 1; i__ <= i__1; ++i__) {
310 d__[i__] = temp * d__[i__] - g[i__] - hs[i__];
311 /* Computing 2nd power */
314 ds += d__[i__] * step[i__];
316 /* Computing 2nd power */
330 /* Test whether an alternative iteration is required. */
333 if (gg <= ggbeg * 1e-4) {
339 for (i__ = 1; i__ <= i__1; ++i__) {
340 sg += step[i__] * g[i__];
342 shs += step[i__] * hs[i__];
345 angtest = sgk / sqrt(gg * delsq);
346 if (angtest <= -.99) {
350 /* Begin the alternative iteration by calculating D and HD and some */
351 /* scalar products. */
354 temp = sqrt(delsq * gg - sgk * sgk);
355 tempa = delsq / temp;
358 for (i__ = 1; i__ <= i__1; ++i__) {
360 d__[i__] = tempa * (g[i__] + hs[i__]) - tempb * step[i__];
368 for (i__ = 1; i__ <= i__1; ++i__) {
369 dg += d__[i__] * g[i__];
370 dhd += hd[i__] * d__[i__];
372 dhs += hd[i__] * step[i__];
375 /* Seek the value of the angle that minimizes Q. */
377 cf = half * (shs - dhd);
383 temp = twopi / (double) (iu + 1);
385 for (i__ = 1; i__ <= i__1; ++i__) {
386 angle = (double) i__ * temp;
389 qnew = (sg + cf * cth) * cth + (dg + dhs * cth) * sth;
394 } else if (i__ == isave + 1) {
400 if ((double) isave == zero) {
407 if (tempa != tempb) {
410 angle = half * (tempa - tempb) / (tempa + tempb);
412 angle = temp * ((double) isave + angle);
414 /* Calculate the new STEP and HS. Then test for convergence. */
418 reduc = qbeg - (sg + cf * cth) * cth - (dg + dhs * cth) * sth;
421 for (i__ = 1; i__ <= i__1; ++i__) {
422 step[i__] = cth * step[i__] + sth * d__[i__];
423 hs[i__] = cth * hs[i__] + sth * hd[i__];
425 /* Computing 2nd power */
426 d__1 = g[i__] + hs[i__];
430 ratio = reduc / qred;
431 if (iterc < itermax && ratio > .01) {
435 return NLOPT_SUCCESS;
437 /* The following instructions act as a subroutine for setting the vector */
438 /* HD to the vector D multiplied by the second derivative matrix of Q. */
439 /* They are called from three different places, which are distinguished */
440 /* by the value of ITERC. */
444 for (i__ = 1; i__ <= i__1; ++i__) {
449 for (k = 1; k <= i__1; ++k) {
452 for (j = 1; j <= i__2; ++j) {
454 temp += xpt[k + j * xpt_dim1] * d__[j];
458 for (i__ = 1; i__ <= i__2; ++i__) {
460 hd[i__] += temp * xpt[k + i__ * xpt_dim1];
465 for (j = 1; j <= i__2; ++j) {
467 for (i__ = 1; i__ <= i__1; ++i__) {
470 hd[j] += hq[ih] * d__[i__];
473 hd[i__] += hq[ih] * d__[j];
479 if (iterc <= itersw) {
486 /*************************************************************************/
489 static void bigden_(int *n, int *npt, double *xopt,
490 double *xpt, double *bmat, double *zmat, int *idz,
491 int *ndim, int *kopt, int *knew, double *d__,
492 double *w, double *vlag, double *beta, double *s,
493 double *wvec, double *prod,
494 const double *xbase, const double *lb, const double *ub)
496 /* System generated locals */
497 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
498 zmat_offset, wvec_dim1, wvec_offset, prod_dim1, prod_offset, i__1,
502 /* Local variables */
508 double ss, den[9], one, par[9], tau, sum, two, diff, half, temp;
512 double zero, alpha, angle, denex[9];
514 double tempa, tempb, tempc;
516 double ssden, dtest, quart, xoptd, twopi, xopts, denold, denmax,
517 densav, dstemp, sumold, sstemp, xoptsq;
520 /* N is the number of variables. */
521 /* NPT is the number of interpolation equations. */
522 /* XOPT is the best interpolation point so far. */
523 /* XPT contains the coordinates of the current interpolation points. */
524 /* BMAT provides the last N columns of H. */
525 /* ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H. */
526 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
527 /* KOPT is the index of the optimal interpolation point. */
528 /* KNEW is the index of the interpolation point that is going to be moved. */
529 /* D will be set to the step from XOPT to the new point, and on entry it */
530 /* should be the D that was calculated by the last call of BIGLAG. The */
531 /* length of the initial D provides a trust region bound on the final D. */
532 /* W will be set to Wcheck for the final choice of D. */
533 /* VLAG will be set to Theta*Wcheck+e_b for the final choice of D. */
534 /* BETA will be set to the value that will occur in the updating formula */
535 /* when the KNEW-th interpolation point is moved to its new position. */
536 /* S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be used */
537 /* for working space. */
539 /* D is calculated in a way that should provide a denominator with a large */
540 /* modulus in the updating formula when the KNEW-th interpolation point is */
541 /* shifted to the new position XOPT+D. */
543 /* Set some constants. */
545 /* Parameter adjustments */
547 zmat_offset = 1 + zmat_dim1;
550 xpt_offset = 1 + xpt_dim1;
554 prod_offset = 1 + prod_dim1;
557 wvec_offset = 1 + wvec_dim1;
560 bmat_offset = 1 + bmat_dim1;
573 twopi = atan(one) * 8.;
574 nptm = *npt - *n - 1;
576 /* Store the first NPT elements of the KNEW-th column of H in W(N+1) */
580 for (k = 1; k <= i__1; ++k) {
585 for (j = 1; j <= i__1; ++j) {
586 temp = zmat[*knew + j * zmat_dim1];
591 for (k = 1; k <= i__2; ++k) {
593 w[*n + k] += temp * zmat[k + j * zmat_dim1];
596 alpha = w[*n + *knew];
598 /* The initial search direction D is taken from the last call of BIGLAG, */
599 /* and the initial S is set below, usually to the direction from X_OPT */
600 /* to X_KNEW, but a different direction to an interpolation point may */
601 /* be chosen, in order to prevent S from being nearly parallel to D. */
608 for (i__ = 1; i__ <= i__2; ++i__) {
609 /* Computing 2nd power */
612 s[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__];
613 ds += d__[i__] * s[i__];
614 /* Computing 2nd power */
618 /* Computing 2nd power */
620 xoptsq += d__1 * d__1;
622 if (ds * ds > dd * .99 * ss) {
624 dtest = ds * ds / ss;
626 for (k = 1; k <= i__2; ++k) {
631 for (i__ = 1; i__ <= i__1; ++i__) {
632 diff = xpt[k + i__ * xpt_dim1] - xopt[i__];
633 dstemp += d__[i__] * diff;
635 sstemp += diff * diff;
637 if (dstemp * dstemp / sstemp < dtest) {
639 dtest = dstemp * dstemp / sstemp;
647 for (i__ = 1; i__ <= i__2; ++i__) {
649 s[i__] = xpt[ksav + i__ * xpt_dim1] - xopt[i__];
652 ssden = dd * ss - ds * ds;
656 /* Begin the iteration by overwriting S with a vector that has the */
657 /* required length and direction. */
661 temp = one / sqrt(ssden);
665 for (i__ = 1; i__ <= i__2; ++i__) {
666 s[i__] = temp * (dd * s[i__] - ds * d__[i__]);
667 xoptd += xopt[i__] * d__[i__];
669 xopts += xopt[i__] * s[i__];
672 /* Set the coefficients of the first two terms of BETA. */
674 tempa = half * xoptd * xoptd;
675 tempb = half * xopts * xopts;
676 den[0] = dd * (xoptsq + half * dd) + tempa + tempb;
677 den[1] = two * xoptd * dd;
678 den[2] = two * xopts * dd;
679 den[3] = tempa - tempb;
680 den[4] = xoptd * xopts;
681 for (i__ = 6; i__ <= 9; ++i__) {
686 /* Put the coefficients of Wcheck in WVEC. */
689 for (k = 1; k <= i__2; ++k) {
694 for (i__ = 1; i__ <= i__1; ++i__) {
695 tempa += xpt[k + i__ * xpt_dim1] * d__[i__];
696 tempb += xpt[k + i__ * xpt_dim1] * s[i__];
698 tempc += xpt[k + i__ * xpt_dim1] * xopt[i__];
700 wvec[k + wvec_dim1] = quart * (tempa * tempa + tempb * tempb);
701 wvec[k + (wvec_dim1 << 1)] = tempa * tempc;
702 wvec[k + wvec_dim1 * 3] = tempb * tempc;
703 wvec[k + (wvec_dim1 << 2)] = quart * (tempa * tempa - tempb * tempb);
705 wvec[k + wvec_dim1 * 5] = half * tempa * tempb;
708 for (i__ = 1; i__ <= i__2; ++i__) {
710 wvec[ip + wvec_dim1] = zero;
711 wvec[ip + (wvec_dim1 << 1)] = d__[i__];
712 wvec[ip + wvec_dim1 * 3] = s[i__];
713 wvec[ip + (wvec_dim1 << 2)] = zero;
715 wvec[ip + wvec_dim1 * 5] = zero;
718 /* Put the coefficents of THETA*Wcheck in PROD. */
720 for (jc = 1; jc <= 5; ++jc) {
722 if (jc == 2 || jc == 3) {
726 for (k = 1; k <= i__2; ++k) {
728 prod[k + jc * prod_dim1] = zero;
731 for (j = 1; j <= i__2; ++j) {
734 for (k = 1; k <= i__1; ++k) {
736 sum += zmat[k + j * zmat_dim1] * wvec[k + jc * wvec_dim1];
742 for (k = 1; k <= i__1; ++k) {
744 prod[k + jc * prod_dim1] += sum * zmat[k + j * zmat_dim1];
749 for (k = 1; k <= i__1; ++k) {
752 for (j = 1; j <= i__2; ++j) {
754 sum += bmat[k + j * bmat_dim1] * wvec[*npt + j + jc *
758 prod[k + jc * prod_dim1] += sum;
762 for (j = 1; j <= i__1; ++j) {
765 for (i__ = 1; i__ <= i__2; ++i__) {
767 sum += bmat[i__ + j * bmat_dim1] * wvec[i__ + jc * wvec_dim1];
770 prod[*npt + j + jc * prod_dim1] = sum;
774 /* Include in DEN the part of BETA that depends on THETA. */
777 for (k = 1; k <= i__1; ++k) {
779 for (i__ = 1; i__ <= 5; ++i__) {
780 par[i__ - 1] = half * prod[k + i__ * prod_dim1] * wvec[k + i__ *
785 den[0] = den[0] - par[0] - sum;
786 tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 1)] + prod[k + (
787 prod_dim1 << 1)] * wvec[k + wvec_dim1];
788 tempb = prod[k + (prod_dim1 << 1)] * wvec[k + (wvec_dim1 << 2)] +
789 prod[k + (prod_dim1 << 2)] * wvec[k + (wvec_dim1 << 1)];
790 tempc = prod[k + prod_dim1 * 3] * wvec[k + wvec_dim1 * 5] + prod[k +
791 prod_dim1 * 5] * wvec[k + wvec_dim1 * 3];
792 den[1] = den[1] - tempa - half * (tempb + tempc);
793 den[5] -= half * (tempb - tempc);
794 tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 3] + prod[k +
795 prod_dim1 * 3] * wvec[k + wvec_dim1];
796 tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 5] + prod[k
797 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 1)];
798 tempc = prod[k + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 2)] + prod[k
799 + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 3];
800 den[2] = den[2] - tempa - half * (tempb - tempc);
801 den[6] -= half * (tempb + tempc);
802 tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 2)] + prod[k + (
803 prod_dim1 << 2)] * wvec[k + wvec_dim1];
804 den[3] = den[3] - tempa - par[1] + par[2];
805 tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 5] + prod[k +
806 prod_dim1 * 5] * wvec[k + wvec_dim1];
807 tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 3] + prod[k
808 + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 1)];
809 den[4] = den[4] - tempa - half * tempb;
810 den[7] = den[7] - par[3] + par[4];
811 tempa = prod[k + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 5] + prod[k
812 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 2)];
814 den[8] -= half * tempa;
817 /* Extend DEN so that it holds all the coefficients of DENOM. */
820 for (i__ = 1; i__ <= 5; ++i__) {
821 /* Computing 2nd power */
822 d__1 = prod[*knew + i__ * prod_dim1];
823 par[i__ - 1] = half * (d__1 * d__1);
827 denex[0] = alpha * den[0] + par[0] + sum;
828 tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 1)];
829 tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + (prod_dim1 << 2)];
830 tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + prod_dim1 * 5];
831 denex[1] = alpha * den[1] + tempa + tempb + tempc;
832 denex[5] = alpha * den[5] + tempb - tempc;
833 tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 3];
834 tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + prod_dim1 * 5];
835 tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + (prod_dim1 << 2)];
836 denex[2] = alpha * den[2] + tempa + tempb - tempc;
837 denex[6] = alpha * den[6] + tempb + tempc;
838 tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 2)];
839 denex[3] = alpha * den[3] + tempa + par[1] - par[2];
840 tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 5];
841 denex[4] = alpha * den[4] + tempa + prod[*knew + (prod_dim1 << 1)] * prod[
842 *knew + prod_dim1 * 3];
843 denex[7] = alpha * den[7] + par[3] - par[4];
844 denex[8] = alpha * den[8] + prod[*knew + (prod_dim1 << 2)] * prod[*knew +
847 /* Seek the value of the angle that maximizes the modulus of DENOM. */
849 sum = denex[0] + denex[1] + denex[3] + denex[5] + denex[7];
854 temp = twopi / (double) (iu + 1);
857 for (i__ = 1; i__ <= i__1; ++i__) {
858 angle = (double) i__ * temp;
861 for (j = 4; j <= 8; j += 2) {
862 par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
864 par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
868 for (j = 1; j <= 9; ++j) {
870 sum += denex[j - 1] * par[j - 1];
872 if (fabs(sum) > fabs(denmax)) {
876 } else if (i__ == isave + 1) {
888 if (tempa != tempb) {
891 step = half * (tempa - tempb) / (tempa + tempb);
893 angle = temp * ((double) isave + step);
895 /* Calculate the new parameters of the denominator, the new VLAG vector */
896 /* and the new D. Then test for convergence. */
900 for (j = 4; j <= 8; j += 2) {
901 par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
903 par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
907 for (j = 1; j <= 9; ++j) {
908 *beta += den[j - 1] * par[j - 1];
910 denmax += denex[j - 1] * par[j - 1];
913 for (k = 1; k <= i__1; ++k) {
915 for (j = 1; j <= 5; ++j) {
917 vlag[k] += prod[k + j * prod_dim1] * par[j - 1];
925 for (i__ = 1; i__ <= i__1; ++i__) {
926 d__[i__] = par[1] * d__[i__] + par[2] * s[i__];
927 w[i__] = xopt[i__] + d__[i__];
928 /* Computing 2nd power */
931 tempa += d__[i__] * w[i__];
933 tempb += w[i__] * w[i__];
939 densav = max(densav,denold);
941 if (fabs(denmax) <= fabs(densav) * 1.1) {
946 /* Set S to half the gradient of the denominator with respect to D. */
947 /* Then branch for the next iteration. */
950 for (i__ = 1; i__ <= i__1; ++i__) {
951 temp = tempa * xopt[i__] + tempb * d__[i__] - vlag[*npt + i__];
953 s[i__] = tau * bmat[*knew + i__ * bmat_dim1] + alpha * temp;
956 for (k = 1; k <= i__1; ++k) {
959 for (j = 1; j <= i__2; ++j) {
961 sum += xpt[k + j * xpt_dim1] * w[j];
963 temp = (tau * w[*n + k] - alpha * vlag[k]) * sum;
965 for (i__ = 1; i__ <= i__2; ++i__) {
967 s[i__] += temp * xpt[k + i__ * xpt_dim1];
973 for (i__ = 1; i__ <= i__2; ++i__) {
974 /* Computing 2nd power */
978 ds += d__[i__] * s[i__];
980 ssden = dd * ss - ds * ds;
981 if (ssden >= dd * 1e-8 * ss) {
985 /* Set the vector W before the RETURN from the subroutine. */
988 /* SGJ, 2008: crude hack: truncate d to [lb,ub] bounds if given */
990 for (k = 1; k <= *n; ++k) {
991 if (d__[k] > ub[k-1] - xbase[k-1] - xopt[k])
992 d__[k] = ub[k-1] - xbase[k-1] - xopt[k];
993 else if (d__[k] < lb[k-1] - xbase[k-1] - xopt[k])
994 d__[k] = lb[k-1] - xbase[k-1] - xopt[k];
999 for (k = 1; k <= i__2; ++k) {
1001 for (j = 1; j <= 5; ++j) {
1003 w[k] += wvec[k + j * wvec_dim1] * par[j - 1];
1010 /*************************************************************************/
1014 int npt, ndim, iter;
1015 double *hcol, *xpt, *bmat, *xopt;
1019 /* the Lagrange function, whose absolute value biglag maximizes */
1020 static double lag(int n, const double *dx, double *grad, void *data)
1022 lag_data *d = (lag_data *) data;
1023 int i, j, npt = d->npt, ndim = d->ndim;
1024 const double *hcol = d->hcol, *xpt = d->xpt, *bmat = d->bmat, *xopt = d->xopt;
1026 for (j = 0; j < n; ++j) {
1027 val += bmat[j * ndim] * (xopt[j] + dx[j]);
1028 if (grad) grad[j] = bmat[j * ndim];
1030 for (i = 0; i < npt; ++i) {
1032 for (j = 0; j < n; ++j)
1033 dot += xpt[i + j * npt] * (xopt[j] + dx[j]);
1034 val += 0.5 * hcol[i] * (dot * dot);
1036 if (grad) for (j = 0; j < n; ++j)
1037 grad[j] += dot * xpt[i + j * npt];
1041 if (grad) for (j = 0; j < n; ++j) grad[j] = -grad[j];
1047 static nlopt_result biglag_(int *n, int *npt, double *xopt,
1048 double *xpt, double *bmat, double *zmat, int *idz,
1049 int *ndim, int *knew, double *delta, double *d__,
1050 double *alpha, double *hcol, double *gc, double *gd,
1051 double *s, double *w,
1052 const double *xbase, const double *lb, const double *ub)
1054 /* System generated locals */
1055 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1056 zmat_offset, i__1, i__2;
1059 /* Local variables */
1063 double sp, ss, cf1, cf2, cf3, cf4, cf5, dhd, cth, one, tau, sth, sum,
1066 double zero, angle, scale, denom;
1068 double delsq, tempa, tempb, twopi, taubeg, tauold, taumax;
1071 /* N is the number of variables. */
1072 /* NPT is the number of interpolation equations. */
1073 /* XOPT is the best interpolation point so far. */
1074 /* XPT contains the coordinates of the current interpolation points. */
1075 /* BMAT provides the last N columns of H. */
1076 /* ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H. */
1077 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
1078 /* KNEW is the index of the interpolation point that is going to be moved. */
1079 /* DELTA is the current trust region bound. */
1080 /* D will be set to the step from XOPT to the new point. */
1081 /* ALPHA will be set to the KNEW-th diagonal element of the H matrix. */
1082 /* HCOL, GC, GD, S and W will be used for working space. */
1084 /* The step D is calculated in a way that attempts to maximize the modulus */
1085 /* of LFUNC(XOPT+D), subject to the bound ||D|| .LE. DELTA, where LFUNC is */
1086 /* the KNEW-th Lagrange function. */
1088 /* Set some constants. */
1090 /* Parameter adjustments */
1092 zmat_offset = 1 + zmat_dim1;
1093 zmat -= zmat_offset;
1095 xpt_offset = 1 + xpt_dim1;
1099 bmat_offset = 1 + bmat_dim1;
1100 bmat -= bmat_offset;
1112 twopi = atan(one) * 8.;
1113 delsq = *delta * *delta;
1114 nptm = *npt - *n - 1;
1116 /* Set the first NPT components of HCOL to the leading elements of the */
1117 /* KNEW-th column of H. */
1121 for (k = 1; k <= i__1; ++k) {
1126 for (j = 1; j <= i__1; ++j) {
1127 temp = zmat[*knew + j * zmat_dim1];
1132 for (k = 1; k <= i__2; ++k) {
1134 hcol[k] += temp * zmat[k + j * zmat_dim1];
1137 *alpha = hcol[*knew];
1139 /* Set the unscaled initial direction D. Form the gradient of LFUNC at */
1140 /* XOPT, and multiply D by the second derivative matrix of LFUNC. */
1144 for (i__ = 1; i__ <= i__2; ++i__) {
1145 d__[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__];
1146 gc[i__] = bmat[*knew + i__ * bmat_dim1];
1149 /* Computing 2nd power */
1154 for (k = 1; k <= i__2; ++k) {
1158 for (j = 1; j <= i__1; ++j) {
1159 temp += xpt[k + j * xpt_dim1] * xopt[j];
1161 sum += xpt[k + j * xpt_dim1] * d__[j];
1163 temp = hcol[k] * temp;
1164 sum = hcol[k] * sum;
1166 for (i__ = 1; i__ <= i__1; ++i__) {
1167 gc[i__] += temp * xpt[k + i__ * xpt_dim1];
1169 gd[i__] += sum * xpt[k + i__ * xpt_dim1];
1173 /* Scale D and GD, with a sign change if required. Set S to another */
1174 /* vector in the initial two dimensional subspace. */
1180 for (i__ = 1; i__ <= i__1; ++i__) {
1181 /* Computing 2nd power */
1184 sp += d__[i__] * gc[i__];
1186 dhd += d__[i__] * gd[i__];
1188 scale = *delta / sqrt(dd);
1189 if (sp * dhd < zero) {
1193 if (sp * sp > dd * .99 * gg) {
1196 tau = scale * (fabs(sp) + half * scale * fabs(dhd));
1197 if (gg * delsq < tau * .01 * tau) {
1201 for (i__ = 1; i__ <= i__1; ++i__) {
1202 d__[i__] = scale * d__[i__];
1203 gd[i__] = scale * gd[i__];
1205 s[i__] = gc[i__] + temp * gd[i__];
1209 double minf, *dlb, *dub, *xtol;
1215 ld.xpt = &xpt[1 + 1 * xpt_dim1];
1216 ld.bmat = &bmat[*knew + 1 * bmat_dim1];
1219 dlb = &gc[1]; dub = &gd[1]; xtol = &s[1];
1220 /* make sure rounding errors don't push initial |d| > delta */
1221 for (j = 1; j <= *n; ++j) d__[j] *= 0.99999;
1222 for (j = 0; j < *n; ++j) {
1223 dlb[j] = -(dub[j] = *delta);
1224 if (dlb[j] < lb[j] - xbase[j] - xopt[j+1])
1225 dlb[j] = lb[j] - xbase[j] - xopt[j+1];
1226 if (dub[j] > ub[j] - xbase[j] - xopt[j+1])
1227 dub[j] = ub[j] - xbase[j] - xopt[j+1];
1228 if (dlb[j] > 0) dlb[j] = 0;
1229 if (dub[j] < 0) dub[j] = 0;
1230 if (d__[j+1] < dlb[j]) d__[j+1] = dlb[j];
1231 else if (d__[j+1] > dub[j]) d__[j+1] = dub[j];
1232 xtol[j] = 1e-5 * *delta;
1234 ld.flipsign = lag(*n, &d__[1], 0, &ld) > 0; /* maximize if > 0 */
1235 return nlopt_minimize_constrained(NLOPT_LD_MMA, *n, lag, &ld,
1236 1, rho_constraint, delta, sizeof(double),
1237 dlb, dub, &d__[1], &minf, -HUGE_VAL,
1238 0., 0., 0., xtol, 1000, 0.);
1241 /* Begin the iteration by overwriting S with a vector that has the */
1242 /* required length and direction, except that termination occurs if */
1243 /* the given D and S are nearly parallel. */
1251 for (i__ = 1; i__ <= i__1; ++i__) {
1252 /* Computing 2nd power */
1255 sp += d__[i__] * s[i__];
1257 /* Computing 2nd power */
1261 temp = dd * ss - sp * sp;
1262 if (temp <= dd * 1e-8 * ss) {
1267 for (i__ = 1; i__ <= i__1; ++i__) {
1268 s[i__] = (dd * s[i__] - sp * d__[i__]) / denom;
1273 /* Calculate the coefficients of the objective function on the circle, */
1274 /* beginning with the multiplication of S by the second derivative matrix. */
1277 for (k = 1; k <= i__1; ++k) {
1280 for (j = 1; j <= i__2; ++j) {
1282 sum += xpt[k + j * xpt_dim1] * s[j];
1284 sum = hcol[k] * sum;
1286 for (i__ = 1; i__ <= i__2; ++i__) {
1288 w[i__] += sum * xpt[k + i__ * xpt_dim1];
1297 for (i__ = 1; i__ <= i__2; ++i__) {
1298 cf1 += s[i__] * w[i__];
1299 cf2 += d__[i__] * gc[i__];
1300 cf3 += s[i__] * gc[i__];
1301 cf4 += d__[i__] * gd[i__];
1303 cf5 += s[i__] * gd[i__];
1306 cf4 = half * cf4 - cf1;
1308 /* Seek the value of the angle that maximizes the modulus of TAU. */
1310 taubeg = cf1 + cf2 + cf4;
1315 temp = twopi / (double) (iu + 1);
1317 for (i__ = 1; i__ <= i__2; ++i__) {
1318 angle = (double) i__ * temp;
1321 tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
1322 if (fabs(tau) > fabs(taumax)) {
1326 } else if (i__ == isave + 1) {
1339 if (tempa != tempb) {
1342 step = half * (tempa - tempb) / (tempa + tempb);
1344 angle = temp * ((double) isave + step);
1346 /* Calculate the new D and GD. Then test for convergence. */
1350 tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
1352 for (i__ = 1; i__ <= i__2; ++i__) {
1353 d__[i__] = cth * d__[i__] + sth * s[i__];
1354 gd[i__] = cth * gd[i__] + sth * w[i__];
1356 s[i__] = gc[i__] + gd[i__];
1358 if (fabs(tau) <= fabs(taubeg) * 1.1) {
1365 return NLOPT_SUCCESS;
1370 /*************************************************************************/
1373 static void update_(int *n, int *npt, double *bmat,
1374 double *zmat, int *idz, int *ndim, double *vlag,
1375 double *beta, int *knew, double *w)
1377 /* System generated locals */
1378 int bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, i__1, i__2;
1381 /* Local variables */
1382 int i__, j, ja, jb, jl, jp;
1383 double one, tau, temp;
1387 double scala, scalb_, alpha, denom, tempa, tempb, tausq;
1390 /* The arrays BMAT and ZMAT with IDZ are updated, in order to shift the */
1391 /* interpolation point that has index KNEW. On entry, VLAG contains the */
1392 /* components of the vector Theta*Wcheck+e_b of the updating formula */
1393 /* (6.11), and BETA holds the value of the parameter that has this name. */
1394 /* The vector W is used for working space. */
1396 /* Set some constants. */
1398 /* Parameter adjustments */
1400 zmat_offset = 1 + zmat_dim1;
1401 zmat -= zmat_offset;
1403 bmat_offset = 1 + bmat_dim1;
1404 bmat -= bmat_offset;
1411 nptm = *npt - *n - 1;
1413 /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */
1417 for (j = 2; j <= i__1; ++j) {
1420 } else if (zmat[*knew + j * zmat_dim1] != zero) {
1421 /* Computing 2nd power */
1422 d__1 = zmat[*knew + jl * zmat_dim1];
1423 /* Computing 2nd power */
1424 d__2 = zmat[*knew + j * zmat_dim1];
1425 temp = sqrt(d__1 * d__1 + d__2 * d__2);
1426 tempa = zmat[*knew + jl * zmat_dim1] / temp;
1427 tempb = zmat[*knew + j * zmat_dim1] / temp;
1429 for (i__ = 1; i__ <= i__2; ++i__) {
1430 temp = tempa * zmat[i__ + jl * zmat_dim1] + tempb * zmat[i__
1432 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1]
1433 - tempb * zmat[i__ + jl * zmat_dim1];
1435 zmat[i__ + jl * zmat_dim1] = temp;
1437 zmat[*knew + j * zmat_dim1] = zero;
1442 /* Put the first NPT components of the KNEW-th column of HLAG into W, */
1443 /* and calculate the parameters of the updating formula. */
1445 tempa = zmat[*knew + zmat_dim1];
1450 tempb = zmat[*knew + jl * zmat_dim1];
1453 for (i__ = 1; i__ <= i__1; ++i__) {
1454 w[i__] = tempa * zmat[i__ + zmat_dim1];
1456 w[i__] += tempb * zmat[i__ + jl * zmat_dim1];
1463 denom = alpha * *beta + tausq;
1466 /* Complete the updating of ZMAT when there is only one nonzero element */
1467 /* in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to one, */
1468 /* then the first column of ZMAT will be exchanged with another one later. */
1472 temp = sqrt((fabs(denom)));
1473 tempb = tempa / temp;
1476 for (i__ = 1; i__ <= i__1; ++i__) {
1478 zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb *
1481 if (*idz == 1 && temp < zero) {
1484 if (*idz >= 2 && temp >= zero) {
1489 /* Complete the updating of ZMAT in the alternative case. */
1492 if (*beta >= zero) {
1496 temp = zmat[*knew + jb * zmat_dim1] / denom;
1497 tempa = temp * *beta;
1499 temp = zmat[*knew + ja * zmat_dim1];
1500 scala = one / sqrt(fabs(*beta) * temp * temp + tausq);
1501 scalb_ = scala * sqrt((fabs(denom)));
1503 for (i__ = 1; i__ <= i__1; ++i__) {
1504 zmat[i__ + ja * zmat_dim1] = scala * (tau * zmat[i__ + ja *
1505 zmat_dim1] - temp * vlag[i__]);
1507 zmat[i__ + jb * zmat_dim1] = scalb_ * (zmat[i__ + jb * zmat_dim1]
1508 - tempa * w[i__] - tempb * vlag[i__]);
1510 if (denom <= zero) {
1514 if (*beta >= zero) {
1520 /* IDZ is reduced in the following case, and usually the first column */
1521 /* of ZMAT is exchanged with a later one. */
1526 for (i__ = 1; i__ <= i__1; ++i__) {
1527 temp = zmat[i__ + zmat_dim1];
1528 zmat[i__ + zmat_dim1] = zmat[i__ + *idz * zmat_dim1];
1530 zmat[i__ + *idz * zmat_dim1] = temp;
1534 /* Finally, update the matrix BMAT. */
1537 for (j = 1; j <= i__1; ++j) {
1539 w[jp] = bmat[*knew + j * bmat_dim1];
1540 tempa = (alpha * vlag[jp] - tau * w[jp]) / denom;
1541 tempb = (-(*beta) * w[jp] - tau * vlag[jp]) / denom;
1543 for (i__ = 1; i__ <= i__2; ++i__) {
1544 bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa *
1545 vlag[i__] + tempb * w[i__];
1547 bmat[jp + (i__ - *npt) * bmat_dim1] = bmat[i__ + j *
1557 /*************************************************************************/
1560 static nlopt_result newuob_(int *n, int *npt, double *x,
1562 const double *lb, const double *ub,
1563 nlopt_stopping *stop, double *minf,
1564 newuoa_func calfun, void *calfun_data,
1565 double *xbase, double *xopt, double *xnew,
1566 double *xpt, double *fval, double *gq, double *hq,
1567 double *pq, double *bmat, double *zmat, int *ndim,
1568 double *d__, double *vlag, double *w)
1570 /* System generated locals */
1571 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1572 zmat_offset, i__1, i__2, i__3;
1573 double d__1, d__2, d__3;
1575 /* Local variables */
1577 int i__, j, k, ih, nf, nh, ip, jp;
1584 double sum, fbeg, diff, half, beta;
1588 double temp, suma, sumb, fopt = HUGE_VAL, bsum, gqsq;
1590 double zero, xipt, xjpt, sumz, diffa, diffb, diffc, hdiag, alpha,
1591 delta, recip, reciq, fsave;
1592 int ksave, nfsav, itemp;
1593 double dnorm, ratio, dstep, tenth, vquad;
1598 double detrat, crvmin;
1602 nlopt_result rc = NLOPT_SUCCESS, rc2;
1604 /* SGJ, 2008: compute rhoend from NLopt stop info */
1605 rhoend = stop->xtol_rel * (*rhobeg);
1606 for (j = 0; j < *n; ++j)
1607 if (rhoend < stop->xtol_abs[j])
1608 rhoend = stop->xtol_abs[j];
1610 /* The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are identical */
1611 /* to the corresponding arguments in SUBROUTINE NEWUOA. */
1612 /* XBASE will hold a shift of origin that should reduce the contributions */
1613 /* from rounding errors to values of the model and Lagrange functions. */
1614 /* XOPT will be set to the displacement from XBASE of the vector of */
1615 /* variables that provides the least calculated F so far. */
1616 /* XNEW will be set to the displacement from XBASE of the vector of */
1617 /* variables for the current calculation of F. */
1618 /* XPT will contain the interpolation point coordinates relative to XBASE. */
1619 /* FVAL will hold the values of F at the interpolation points. */
1620 /* GQ will hold the gradient of the quadratic model at XBASE. */
1621 /* HQ will hold the explicit second derivatives of the quadratic model. */
1622 /* PQ will contain the parameters of the implicit second derivatives of */
1623 /* the quadratic model. */
1624 /* BMAT will hold the last N columns of H. */
1625 /* ZMAT will hold the factorization of the leading NPT by NPT submatrix of */
1626 /* H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, where */
1627 /* the elements of DZ are plus or minus one, as specified by IDZ. */
1628 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
1629 /* D is reserved for trial steps from XOPT. */
1630 /* VLAG will contain the values of the Lagrange functions at a new point X. */
1631 /* They are part of a product that requires VLAG to be of length NDIM. */
1632 /* The array W will be used for working space. Its length must be at least */
1633 /* 10*NDIM = 10*(NPT+N). */
1635 /* Set some constants. */
1637 /* Parameter adjustments */
1639 zmat_offset = 1 + zmat_dim1;
1640 zmat -= zmat_offset;
1642 xpt_offset = 1 + xpt_dim1;
1653 bmat_offset = 1 + bmat_dim1;
1654 bmat -= bmat_offset;
1668 /* Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */
1671 for (j = 1; j <= i__1; ++j) {
1674 for (k = 1; k <= i__2; ++k) {
1676 xpt[k + j * xpt_dim1] = zero;
1679 for (i__ = 1; i__ <= i__2; ++i__) {
1681 bmat[i__ + j * bmat_dim1] = zero;
1685 for (ih = 1; ih <= i__2; ++ih) {
1690 for (k = 1; k <= i__2; ++k) {
1693 for (j = 1; j <= i__1; ++j) {
1695 zmat[k + j * zmat_dim1] = zero;
1699 /* Begin the initialization procedure. NF becomes one more than the number */
1700 /* of function values so far. The coordinates of the displacement of the */
1701 /* next initial interpolation point from XBASE are set in XPT(NF,.). */
1703 rhosq = *rhobeg * *rhobeg;
1704 recip = one / rhosq;
1705 reciq = sqrt(half) / rhosq;
1711 if (nfm <= *n << 1) {
1712 if (nfm >= 1 && nfm <= *n) {
1713 xpt[nf + nfm * xpt_dim1] = *rhobeg;
1714 } else if (nfm > *n) {
1715 xpt[nf + nfmm * xpt_dim1] = -(*rhobeg);
1718 itemp = (nfmm - 1) / *n;
1719 jpt = nfm - itemp * *n - *n;
1727 if (fval[ipt + np] < fval[ipt + 1]) {
1731 if (fval[jpt + np] < fval[jpt + 1]) {
1734 xpt[nf + ipt * xpt_dim1] = xipt;
1735 xpt[nf + jpt * xpt_dim1] = xjpt;
1738 /* Calculate the next value of F, label 70 being reached immediately */
1739 /* after this calculation. The least function value so far and its index */
1743 for (j = 1; j <= i__1; ++j) {
1745 x[j] = xpt[nf + j * xpt_dim1] + xbase[j];
1747 if (lb && ub) { /* SGJ, 2008: make sure we are within bounds */
1748 for (j = 1; j <= i__1; ++j) {
1749 if (x[j] < lb[j-1]) x[j] = lb[j-1];
1750 else if (x[j] > ub[j-1]) x[j] = ub[j-1];
1760 } else if (f < fopt) {
1765 /* Set the nonzero initial elements of BMAT and the quadratic model in */
1766 /* the cases when NF is at most 2*N+1. */
1768 if (nfm <= *n << 1) {
1769 if (nfm >= 1 && nfm <= *n) {
1770 gq[nfm] = (f - fbeg) / *rhobeg;
1771 if (*npt < nf + *n) {
1772 bmat[nfm * bmat_dim1 + 1] = -one / *rhobeg;
1773 bmat[nf + nfm * bmat_dim1] = one / *rhobeg;
1774 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
1776 } else if (nfm > *n) {
1777 bmat[nf - *n + nfmm * bmat_dim1] = half / *rhobeg;
1778 bmat[nf + nfmm * bmat_dim1] = -half / *rhobeg;
1779 zmat[nfmm * zmat_dim1 + 1] = -reciq - reciq;
1780 zmat[nf - *n + nfmm * zmat_dim1] = reciq;
1781 zmat[nf + nfmm * zmat_dim1] = reciq;
1782 ih = nfmm * (nfmm + 1) / 2;
1783 temp = (fbeg - f) / *rhobeg;
1784 hq[ih] = (gq[nfmm] - temp) / *rhobeg;
1785 gq[nfmm] = half * (gq[nfmm] + temp);
1788 /* Set the off-diagonal second derivatives of the Lagrange functions and */
1789 /* the initial quadratic model. */
1792 ih = ipt * (ipt - 1) / 2 + jpt;
1799 zmat[nfmm * zmat_dim1 + 1] = recip;
1800 zmat[nf + nfmm * zmat_dim1] = recip;
1801 zmat[ipt + 1 + nfmm * zmat_dim1] = -recip;
1802 zmat[jpt + 1 + nfmm * zmat_dim1] = -recip;
1803 hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / (xipt * xjpt);
1809 /* Begin the iterative procedure, because the initial model is complete. */
1819 for (i__ = 1; i__ <= i__1; ++i__) {
1820 xopt[i__] = xpt[kopt + i__ * xpt_dim1];
1822 /* Computing 2nd power */
1824 xoptsq += d__1 * d__1;
1829 /* Generate the next trust region step and test its length. Set KNEW */
1830 /* to -1 if the purpose of the next F will be to improve the model. */
1834 rc2 = trsapp_(n, npt, &xopt[1], &xpt[xpt_offset], &gq[1], &hq[1], &pq[1], &
1835 delta, &d__[1], &w[1], &w[np], &w[np + *n], &w[np + (*n << 1)], &
1836 crvmin, &xbase[1], lb, ub);
1837 if (rc2 < 0) { rc = rc2; goto L530; }
1840 for (i__ = 1; i__ <= i__1; ++i__) {
1842 /* Computing 2nd power */
1847 d__1 = delta, d__2 = sqrt(dsq);
1848 dnorm = min(d__1,d__2);
1849 if (dnorm < half * rho) {
1851 delta = tenth * delta;
1853 if (delta <= rho * 1.5) {
1856 if (nf <= nfsav + 2) {
1859 temp = crvmin * .125 * rho * rho;
1861 d__1 = max(diffa,diffb);
1862 if (temp <= max(d__1,diffc)) {
1868 /* Shift XBASE if XOPT may be too far from XBASE. First make the changes */
1869 /* to BMAT that do not depend on ZMAT. */
1872 if (dsq <= xoptsq * .001) {
1873 tempq = xoptsq * .25;
1875 for (k = 1; k <= i__1; ++k) {
1878 for (i__ = 1; i__ <= i__2; ++i__) {
1880 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
1883 sum -= half * xoptsq;
1886 for (i__ = 1; i__ <= i__2; ++i__) {
1887 gq[i__] += temp * xpt[k + i__ * xpt_dim1];
1888 xpt[k + i__ * xpt_dim1] -= half * xopt[i__];
1889 vlag[i__] = bmat[k + i__ * bmat_dim1];
1890 w[i__] = sum * xpt[k + i__ * xpt_dim1] + tempq * xopt[i__];
1893 for (j = 1; j <= i__3; ++j) {
1895 bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] +
1896 vlag[i__] * w[j] + w[i__] * vlag[j];
1901 /* Then the revisions of BMAT that depend on ZMAT are calculated. */
1904 for (k = 1; k <= i__3; ++k) {
1907 for (i__ = 1; i__ <= i__2; ++i__) {
1908 sumz += zmat[i__ + k * zmat_dim1];
1910 w[i__] = w[*npt + i__] * zmat[i__ + k * zmat_dim1];
1913 for (j = 1; j <= i__2; ++j) {
1914 sum = tempq * sumz * xopt[j];
1916 for (i__ = 1; i__ <= i__1; ++i__) {
1918 sum += w[i__] * xpt[i__ + j * xpt_dim1];
1925 for (i__ = 1; i__ <= i__1; ++i__) {
1927 bmat[i__ + j * bmat_dim1] += sum * zmat[i__ + k *
1932 for (i__ = 1; i__ <= i__1; ++i__) {
1939 for (j = 1; j <= i__2; ++j) {
1941 bmat[ip + j * bmat_dim1] += temp * vlag[j];
1946 /* The following instructions complete the shift of XBASE, including */
1947 /* the changes to the parameters of the quadratic model. */
1951 for (j = 1; j <= i__2; ++j) {
1954 for (k = 1; k <= i__1; ++k) {
1955 w[j] += pq[k] * xpt[k + j * xpt_dim1];
1957 xpt[k + j * xpt_dim1] -= half * xopt[j];
1960 for (i__ = 1; i__ <= i__1; ++i__) {
1963 gq[j] += hq[ih] * xopt[i__];
1965 gq[i__] += hq[ih] * xopt[j];
1966 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
1968 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
1973 for (j = 1; j <= i__1; ++j) {
1974 xbase[j] += xopt[j];
1981 /* Pick the model step if KNEW is positive. A different choice of D */
1982 /* may be made later, if the choice of D by BIGLAG causes substantial */
1983 /* cancellation in DENOM. */
1987 biglag_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &zmat[
1988 zmat_offset], &idz, ndim, &knew, &dstep, &d__[1], &alpha, &
1989 vlag[1], &vlag[*npt + 1], &w[1], &w[np], &w[np + *n],
1991 if (rc2 < 0) { rc = rc2; goto L530; }
1994 /* Calculate VLAG and BETA for the current choice of D. The first NPT */
1995 /* components of W_check will be held in W. */
1998 for (k = 1; k <= i__1; ++k) {
2003 for (j = 1; j <= i__2; ++j) {
2004 suma += xpt[k + j * xpt_dim1] * d__[j];
2005 sumb += xpt[k + j * xpt_dim1] * xopt[j];
2007 sum += bmat[k + j * bmat_dim1] * d__[j];
2009 w[k] = suma * (half * suma + sumb);
2015 for (k = 1; k <= i__1; ++k) {
2018 for (i__ = 1; i__ <= i__2; ++i__) {
2020 sum += zmat[i__ + k * zmat_dim1] * w[i__];
2029 for (i__ = 1; i__ <= i__2; ++i__) {
2031 vlag[i__] += sum * zmat[i__ + k * zmat_dim1];
2037 for (j = 1; j <= i__2; ++j) {
2040 for (i__ = 1; i__ <= i__1; ++i__) {
2042 sum += w[i__] * bmat[i__ + j * bmat_dim1];
2044 bsum += sum * d__[j];
2047 for (k = 1; k <= i__1; ++k) {
2049 sum += bmat[jp + k * bmat_dim1] * d__[k];
2052 bsum += sum * d__[j];
2054 dx += d__[j] * xopt[j];
2056 beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2059 /* If KNEW is positive and if the cancellation in DENOM is unacceptable, */
2060 /* then BIGDEN calculates an alternative model step, XNEW being used for */
2061 /* working space. */
2064 /* Computing 2nd power */
2066 temp = one + alpha * beta / (d__1 * d__1);
2067 if (fabs(temp) <= .8) {
2068 bigden_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &
2069 zmat[zmat_offset], &idz, ndim, &kopt, &knew, &d__[1], &w[
2070 1], &vlag[1], &beta, &xnew[1], &w[*ndim + 1], &w[*ndim *
2071 6 + 1], &xbase[1], lb, ub);
2075 /* Calculate the next value of the objective function. */
2079 for (i__ = 1; i__ <= i__2; ++i__) {
2080 xnew[i__] = xopt[i__] + d__[i__];
2082 x[i__] = xbase[i__] + xnew[i__];
2086 if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
2087 else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
2088 if (rc != NLOPT_SUCCESS) goto L530;
2091 f = calfun(*n, &x[1], calfun_data);
2092 if (f < stop->minf_max) {
2093 rc = NLOPT_MINF_MAX_REACHED;
2104 /* Use the quadratic model to predict the change in F due to the step D, */
2105 /* and set DIFF to the error of this prediction. */
2110 for (j = 1; j <= i__2; ++j) {
2111 vquad += d__[j] * gq[j];
2113 for (i__ = 1; i__ <= i__1; ++i__) {
2115 temp = d__[i__] * xnew[j] + d__[j] * xopt[i__];
2120 vquad += temp * hq[ih];
2124 for (k = 1; k <= i__1; ++k) {
2126 vquad += pq[k] * w[k];
2128 diff = f - fopt - vquad;
2136 /* Update FOPT and XOPT if the new F is the least value of the objective */
2137 /* function so far. The branch when KNEW is positive occurs if D is not */
2138 /* a trust region step. */
2145 for (i__ = 1; i__ <= i__1; ++i__) {
2146 xopt[i__] = xnew[i__];
2148 /* Computing 2nd power */
2150 xoptsq += d__1 * d__1;
2152 if (nlopt_stop_ftol(stop, fopt, fsave)) {
2153 rc = NLOPT_FTOL_REACHED;
2163 /* Pick the next value of DELTA after a trust region step. */
2165 if (vquad >= zero) {
2168 ratio = (f - fsave) / vquad;
2169 if (ratio <= tenth) {
2170 delta = half * dnorm;
2171 } else if (ratio <= .7) {
2173 d__1 = half * delta;
2174 delta = max(d__1,dnorm);
2177 d__1 = half * delta, d__2 = dnorm + dnorm;
2178 delta = max(d__1,d__2);
2180 if (delta <= rho * 1.5) {
2184 /* Set KNEW to the index of the next interpolation point to be deleted. */
2187 d__2 = tenth * delta;
2188 /* Computing 2nd power */
2189 d__1 = max(d__2,rho);
2190 rhosq = d__1 * d__1;
2198 for (k = 1; k <= i__1; ++k) {
2201 for (j = 1; j <= i__2; ++j) {
2207 /* Computing 2nd power */
2208 d__1 = zmat[k + j * zmat_dim1];
2209 hdiag += temp * (d__1 * d__1);
2211 /* Computing 2nd power */
2213 temp = (d__1 = beta * hdiag + d__2 * d__2, fabs(d__1));
2216 for (j = 1; j <= i__2; ++j) {
2218 /* Computing 2nd power */
2219 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2220 distsq += d__1 * d__1;
2222 if (distsq > rhosq) {
2223 /* Computing 3rd power */
2224 d__1 = distsq / rhosq;
2225 temp *= d__1 * (d__1 * d__1);
2227 if (temp > detrat && k != ktemp) {
2237 /* Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point */
2238 /* can be moved. Begin the updating of the quadratic model, starting */
2239 /* with the explicit second derivative term. */
2242 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], &idz, ndim, &vlag[
2243 1], &beta, &knew, &w[1]);
2247 for (i__ = 1; i__ <= i__1; ++i__) {
2248 temp = pq[knew] * xpt[knew + i__ * xpt_dim1];
2250 for (j = 1; j <= i__2; ++j) {
2253 hq[ih] += temp * xpt[knew + j * xpt_dim1];
2258 /* Update the other second derivative parameters, and then the gradient */
2259 /* vector of the model. Also include the new interpolation point. */
2262 for (j = 1; j <= i__2; ++j) {
2263 temp = diff * zmat[knew + j * zmat_dim1];
2268 for (k = 1; k <= i__1; ++k) {
2270 pq[k] += temp * zmat[k + j * zmat_dim1];
2275 for (i__ = 1; i__ <= i__1; ++i__) {
2276 gq[i__] += diff * bmat[knew + i__ * bmat_dim1];
2277 /* Computing 2nd power */
2279 gqsq += d__1 * d__1;
2281 xpt[knew + i__ * xpt_dim1] = xnew[i__];
2284 /* If a trust region step makes a small change to the objective function, */
2285 /* then calculate the gradient of the least Frobenius norm interpolant at */
2286 /* XBASE, and store it in W, using VLAG for a vector of right hand sides. */
2288 if (ksave == 0 && delta == rho) {
2289 if (fabs(ratio) > .01) {
2293 for (k = 1; k <= i__1; ++k) {
2295 vlag[k] = fval[k] - fval[kopt];
2299 for (i__ = 1; i__ <= i__1; ++i__) {
2302 for (k = 1; k <= i__2; ++k) {
2304 sum += bmat[k + i__ * bmat_dim1] * vlag[k];
2311 /* Test whether to replace the new quadratic model by the least Frobenius */
2312 /* norm interpolant, making the replacement if the test is satisfied. */
2315 if (gqsq < gisq * 100.) {
2320 for (i__ = 1; i__ <= i__1; ++i__) {
2325 for (ih = 1; ih <= i__1; ++ih) {
2330 for (j = 1; j <= i__1; ++j) {
2333 for (k = 1; k <= i__2; ++k) {
2335 w[j] += vlag[k] * zmat[k + j * zmat_dim1];
2343 for (k = 1; k <= i__1; ++k) {
2346 for (j = 1; j <= i__2; ++j) {
2348 pq[k] += zmat[k + j * zmat_dim1] * w[j];
2359 /* If a trust region step has provided a sufficient decrease in F, then */
2360 /* branch for another trust region calculation. The case KSAVE>0 occurs */
2361 /* when the new function value was calculated by a model step. */
2363 if (f <= fsave + tenth * vquad) {
2370 /* Alternatively, find out if the interpolation points are close enough */
2371 /* to the best point so far. */
2375 distsq = delta * 4. * delta;
2377 for (k = 1; k <= i__2; ++k) {
2380 for (j = 1; j <= i__1; ++j) {
2382 /* Computing 2nd power */
2383 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2393 /* If KNEW is positive, then set DSTEP, and branch back for the next */
2394 /* iteration, which will generate a "model step". */
2399 d__2 = tenth * sqrt(distsq), d__3 = half * delta;
2400 d__1 = min(d__2,d__3);
2401 dstep = max(d__1,rho);
2402 dsq = dstep * dstep;
2408 if (max(delta,dnorm) > rho) {
2412 /* The calculations with the current value of RHO are complete. Pick the */
2413 /* next values of RHO and DELTA. */
2418 ratio = rho / rhoend;
2421 } else if (ratio <= 250.) {
2422 rho = sqrt(ratio) * rhoend;
2426 delta = max(delta,rho);
2430 /* Return from the calculation, after another Newton-Raphson step, if */
2431 /* it is too short to have been tried before. */
2436 rc = NLOPT_XTOL_REACHED;
2440 for (i__ = 1; i__ <= i__2; ++i__) {
2442 x[i__] = xbase[i__] + xopt[i__];
2450 /*************************************************************************/
2453 nlopt_result newuoa(int n, int npt, double *x,
2454 const double *lb, const double *ub,
2455 double rhobeg, nlopt_stopping *stop, double *minf,
2456 newuoa_func calfun, void *calfun_data)
2458 /* Local variables */
2459 int id, np, iw, igq, ihq, ixb, ifv, ipq, ivl, ixn, ixo, ixp, ndim,
2464 /* This subroutine seeks the least value of a function of many variables, */
2465 /* by a trust region method that forms quadratic models by interpolation. */
2466 /* There can be some freedom in the interpolation conditions, which is */
2467 /* taken up by minimizing the Frobenius norm of the change to the second */
2468 /* derivative of the quadratic model, beginning with a zero matrix. The */
2469 /* arguments of the subroutine are as follows. */
2471 /* N must be set to the number of variables and must be at least two. */
2472 /* NPT is the number of interpolation conditions. Its value must be in the */
2473 /* interval [N+2,(N+1)(N+2)/2]. */
2474 /* Initial values of the variables must be set in X(1),X(2),...,X(N). They */
2475 /* will be changed to the values that give the least calculated F. */
2476 /* RHOBEG and RHOEND must be set to the initial and final values of a trust */
2477 /* region radius, so both must be positive with RHOEND<=RHOBEG. Typically */
2478 /* RHOBEG should be about one tenth of the greatest expected change to a */
2479 /* variable, and RHOEND should indicate the accuracy that is required in */
2480 /* the final values of the variables. */
2481 /* The value of IPRINT should be set to 0, 1, 2 or 3, which controls the */
2482 /* amount of printing. Specifically, there is no output if IPRINT=0 and */
2483 /* there is output only at the return if IPRINT=1. Otherwise, each new */
2484 /* value of RHO is printed, with the best vector of variables so far and */
2485 /* the corresponding value of the objective function. Further, each new */
2486 /* value of F with its variables are output if IPRINT=3. */
2487 /* MAXFUN must be set to an upper bound on the number of calls of CALFUN. */
2488 /* The array W will be used for working space. Its length must be at least */
2489 /* (NPT+13)*(NPT+N)+3*N*(N+3)/2. */
2491 /* SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set F to */
2492 /* the value of the objective function for the variables X(1),X(2),...,X(N). */
2494 /* Partition the working space array, so that different parts of it can be */
2495 /* treated separately by the subroutine that performs the main calculation. */
2497 /* Parameter adjustments */
2503 if (n < 2 || npt < n + 2 || npt > (n + 2) * np / 2) {
2504 return NLOPT_INVALID_ARGS;
2511 ifv = ixp + n * npt;
2514 ipq = ihq + n * np / 2;
2516 izmat = ibmat + ndim * n;
2517 id = izmat + npt * nptm;
2521 w = (double *) malloc(sizeof(double) * ((npt+13)*(npt+n) + 3*(n*(n+3))/2));
2522 if (!w) return NLOPT_OUT_OF_MEMORY;
2525 /* The above settings provide a partition of W for subroutine NEWUOB. */
2526 /* The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements of */
2527 /* W plus the space that is needed by the last array of NEWUOB. */
2529 ret = newuob_(&n, &npt, &x[1], &rhobeg,
2530 lb, ub, stop, minf, calfun, calfun_data,
2531 &w[ixb], &w[ixo], &w[ixn], &w[ixp], &w[ifv],
2532 &w[igq], &w[ihq], &w[ipq], &w[ibmat], &w[izmat],
2533 &ndim, &w[id], &w[ivl], &w[iw]);
2540 /*************************************************************************/
2541 /*************************************************************************/