1 /* Copyright (c) 2004 M. J. D. Powell (mjdp@cam.ac.uk)
2 * Copyright (c) 2007-2014 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 MIN2(a,b) ((a) <= (b) ? (a) : (b))
36 #define MAX2(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 = 0.0, cf, dg, gg = 0.0;
115 double ss, dhd, dhs, cth, sgk, shs = 0.0, sth, qadd, half, qbeg, qred = 0.0, qmin,
116 temp, qsav, qnew, zero, ggbeg = 0.0, alpha, angle, reduc;
118 double ggsav, delsq, tempa = 0.0, tempb = 0.0;
120 double bstep = 0.0, 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 = MIN2(*crvmin,temp);
272 d__1 = alpha, d__2 = gg / dhd;
273 alpha = MIN2(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 nlopt_result 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 (sstemp == 0) return NLOPT_ROUNDOFF_LIMITED;
638 if (dstemp * dstemp / sstemp < dtest) {
640 dtest = dstemp * dstemp / sstemp;
648 for (i__ = 1; i__ <= i__2; ++i__) {
650 s[i__] = xpt[ksav + i__ * xpt_dim1] - xopt[i__];
653 ssden = dd * ss - ds * ds;
657 /* Begin the iteration by overwriting S with a vector that has the */
658 /* required length and direction. */
662 if (ssden < 0) return NLOPT_ROUNDOFF_LIMITED;
663 temp = one / sqrt(ssden);
667 for (i__ = 1; i__ <= i__2; ++i__) {
668 s[i__] = temp * (dd * s[i__] - ds * d__[i__]);
669 xoptd += xopt[i__] * d__[i__];
671 if (nlopt_isinf(s[i__])) return NLOPT_ROUNDOFF_LIMITED;
672 xopts += xopt[i__] * s[i__];
675 /* Set the coefficients of the first two terms of BETA. */
677 tempa = half * xoptd * xoptd;
678 tempb = half * xopts * xopts;
679 den[0] = dd * (xoptsq + half * dd) + tempa + tempb;
680 den[1] = two * xoptd * dd;
681 den[2] = two * xopts * dd;
682 den[3] = tempa - tempb;
683 den[4] = xoptd * xopts;
684 for (i__ = 6; i__ <= 9; ++i__) {
689 /* Put the coefficients of Wcheck in WVEC. */
692 for (k = 1; k <= i__2; ++k) {
697 for (i__ = 1; i__ <= i__1; ++i__) {
698 tempa += xpt[k + i__ * xpt_dim1] * d__[i__];
699 tempb += xpt[k + i__ * xpt_dim1] * s[i__];
701 tempc += xpt[k + i__ * xpt_dim1] * xopt[i__];
703 wvec[k + wvec_dim1] = quart * (tempa * tempa + tempb * tempb);
704 wvec[k + (wvec_dim1 << 1)] = tempa * tempc;
705 wvec[k + wvec_dim1 * 3] = tempb * tempc;
706 wvec[k + (wvec_dim1 << 2)] = quart * (tempa * tempa - tempb * tempb);
708 wvec[k + wvec_dim1 * 5] = half * tempa * tempb;
711 for (i__ = 1; i__ <= i__2; ++i__) {
713 wvec[ip + wvec_dim1] = zero;
714 wvec[ip + (wvec_dim1 << 1)] = d__[i__];
715 wvec[ip + wvec_dim1 * 3] = s[i__];
716 wvec[ip + (wvec_dim1 << 2)] = zero;
718 wvec[ip + wvec_dim1 * 5] = zero;
721 /* Put the coefficents of THETA*Wcheck in PROD. */
723 for (jc = 1; jc <= 5; ++jc) {
725 if (jc == 2 || jc == 3) {
729 for (k = 1; k <= i__2; ++k) {
731 prod[k + jc * prod_dim1] = zero;
734 for (j = 1; j <= i__2; ++j) {
737 for (k = 1; k <= i__1; ++k) {
739 sum += zmat[k + j * zmat_dim1] * wvec[k + jc * wvec_dim1];
745 for (k = 1; k <= i__1; ++k) {
747 prod[k + jc * prod_dim1] += sum * zmat[k + j * zmat_dim1];
752 for (k = 1; k <= i__1; ++k) {
755 for (j = 1; j <= i__2; ++j) {
757 sum += bmat[k + j * bmat_dim1] * wvec[*npt + j + jc *
761 prod[k + jc * prod_dim1] += sum;
765 for (j = 1; j <= i__1; ++j) {
768 for (i__ = 1; i__ <= i__2; ++i__) {
770 sum += bmat[i__ + j * bmat_dim1] * wvec[i__ + jc * wvec_dim1];
773 prod[*npt + j + jc * prod_dim1] = sum;
777 /* Include in DEN the part of BETA that depends on THETA. */
780 for (k = 1; k <= i__1; ++k) {
782 for (i__ = 1; i__ <= 5; ++i__) {
783 par[i__ - 1] = half * prod[k + i__ * prod_dim1] * wvec[k + i__ *
788 den[0] = den[0] - par[0] - sum;
789 tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 1)] + prod[k + (
790 prod_dim1 << 1)] * wvec[k + wvec_dim1];
791 tempb = prod[k + (prod_dim1 << 1)] * wvec[k + (wvec_dim1 << 2)] +
792 prod[k + (prod_dim1 << 2)] * wvec[k + (wvec_dim1 << 1)];
793 tempc = prod[k + prod_dim1 * 3] * wvec[k + wvec_dim1 * 5] + prod[k +
794 prod_dim1 * 5] * wvec[k + wvec_dim1 * 3];
795 den[1] = den[1] - tempa - half * (tempb + tempc);
796 den[5] -= half * (tempb - tempc);
797 tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 3] + prod[k +
798 prod_dim1 * 3] * wvec[k + wvec_dim1];
799 tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 5] + prod[k
800 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 1)];
801 tempc = prod[k + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 2)] + prod[k
802 + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 3];
803 den[2] = den[2] - tempa - half * (tempb - tempc);
804 den[6] -= half * (tempb + tempc);
805 tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 2)] + prod[k + (
806 prod_dim1 << 2)] * wvec[k + wvec_dim1];
807 den[3] = den[3] - tempa - par[1] + par[2];
808 tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 5] + prod[k +
809 prod_dim1 * 5] * wvec[k + wvec_dim1];
810 tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 3] + prod[k
811 + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 1)];
812 den[4] = den[4] - tempa - half * tempb;
813 den[7] = den[7] - par[3] + par[4];
814 tempa = prod[k + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 5] + prod[k
815 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 2)];
817 den[8] -= half * tempa;
820 /* Extend DEN so that it holds all the coefficients of DENOM. */
823 for (i__ = 1; i__ <= 5; ++i__) {
824 /* Computing 2nd power */
825 d__1 = prod[*knew + i__ * prod_dim1];
826 par[i__ - 1] = half * (d__1 * d__1);
830 denex[0] = alpha * den[0] + par[0] + sum;
831 tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 1)];
832 tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + (prod_dim1 << 2)];
833 tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + prod_dim1 * 5];
834 denex[1] = alpha * den[1] + tempa + tempb + tempc;
835 denex[5] = alpha * den[5] + tempb - tempc;
836 tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 3];
837 tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + prod_dim1 * 5];
838 tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + (prod_dim1 << 2)];
839 denex[2] = alpha * den[2] + tempa + tempb - tempc;
840 denex[6] = alpha * den[6] + tempb + tempc;
841 tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 2)];
842 denex[3] = alpha * den[3] + tempa + par[1] - par[2];
843 tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 5];
844 denex[4] = alpha * den[4] + tempa + prod[*knew + (prod_dim1 << 1)] * prod[
845 *knew + prod_dim1 * 3];
846 denex[7] = alpha * den[7] + par[3] - par[4];
847 denex[8] = alpha * den[8] + prod[*knew + (prod_dim1 << 2)] * prod[*knew +
850 /* Seek the value of the angle that maximizes the modulus of DENOM. */
852 sum = denex[0] + denex[1] + denex[3] + denex[5] + denex[7];
857 temp = twopi / (double) (iu + 1);
860 for (i__ = 1; i__ <= i__1; ++i__) {
861 angle = (double) i__ * temp;
864 for (j = 4; j <= 8; j += 2) {
865 par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
867 par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
871 for (j = 1; j <= 9; ++j) {
873 sum += denex[j - 1] * par[j - 1];
875 if (fabs(sum) > fabs(denmax)) {
879 } else if (i__ == isave + 1) {
891 if (tempa != tempb) {
894 step = half * (tempa - tempb) / (tempa + tempb);
896 angle = temp * ((double) isave + step);
898 /* Calculate the new parameters of the denominator, the new VLAG vector */
899 /* and the new D. Then test for convergence. */
903 for (j = 4; j <= 8; j += 2) {
904 par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
906 par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
910 for (j = 1; j <= 9; ++j) {
911 *beta += den[j - 1] * par[j - 1];
913 denmax += denex[j - 1] * par[j - 1];
916 for (k = 1; k <= i__1; ++k) {
918 for (j = 1; j <= 5; ++j) {
920 vlag[k] += prod[k + j * prod_dim1] * par[j - 1];
928 for (i__ = 1; i__ <= i__1; ++i__) {
929 d__[i__] = par[1] * d__[i__] + par[2] * s[i__];
930 w[i__] = xopt[i__] + d__[i__];
931 /* Computing 2nd power */
934 tempa += d__[i__] * w[i__];
936 tempb += w[i__] * w[i__];
942 densav = MAX2(densav,denold);
944 if (fabs(denmax) <= fabs(densav) * 1.1) {
949 /* Set S to half the gradient of the denominator with respect to D. */
950 /* Then branch for the next iteration. */
953 for (i__ = 1; i__ <= i__1; ++i__) {
954 temp = tempa * xopt[i__] + tempb * d__[i__] - vlag[*npt + i__];
956 s[i__] = tau * bmat[*knew + i__ * bmat_dim1] + alpha * temp;
959 for (k = 1; k <= i__1; ++k) {
962 for (j = 1; j <= i__2; ++j) {
964 sum += xpt[k + j * xpt_dim1] * w[j];
966 if (nlopt_isinf(tau * w[*n + k]) ||
967 nlopt_isinf(alpha * vlag[k])) return NLOPT_ROUNDOFF_LIMITED;
968 temp = (tau * w[*n + k] - alpha * vlag[k]) * sum;
970 for (i__ = 1; i__ <= i__2; ++i__) {
972 s[i__] += temp * xpt[k + i__ * xpt_dim1];
978 for (i__ = 1; i__ <= i__2; ++i__) {
979 /* Computing 2nd power */
983 ds += d__[i__] * s[i__];
985 ssden = dd * ss - ds * ds;
986 if (ssden >= dd * 1e-8 * ss) {
990 /* Set the vector W before the RETURN from the subroutine. */
993 /* SGJ, 2008: crude hack: truncate d to [lb,ub] bounds if given */
995 for (k = 1; k <= *n; ++k) {
996 if (d__[k] > ub[k-1] - xbase[k-1] - xopt[k])
997 d__[k] = ub[k-1] - xbase[k-1] - xopt[k];
998 else if (d__[k] < lb[k-1] - xbase[k-1] - xopt[k])
999 d__[k] = lb[k-1] - xbase[k-1] - xopt[k];
1004 for (k = 1; k <= i__2; ++k) {
1006 for (j = 1; j <= 5; ++j) {
1008 w[k] += wvec[k + j * wvec_dim1] * par[j - 1];
1012 return NLOPT_SUCCESS;
1015 /*************************************************************************/
1019 int npt, ndim, iter;
1020 double *hcol, *xpt, *bmat, *xopt;
1024 /* the Lagrange function, whose absolute value biglag maximizes */
1025 static double lag(int n, const double *dx, double *grad, void *data)
1027 lag_data *d = (lag_data *) data;
1028 int i, j, npt = d->npt, ndim = d->ndim;
1029 const double *hcol = d->hcol, *xpt = d->xpt, *bmat = d->bmat, *xopt = d->xopt;
1031 for (j = 0; j < n; ++j) {
1032 val += bmat[j * ndim] * (xopt[j] + dx[j]);
1033 if (grad) grad[j] = bmat[j * ndim];
1035 for (i = 0; i < npt; ++i) {
1037 for (j = 0; j < n; ++j)
1038 dot += xpt[i + j * npt] * (xopt[j] + dx[j]);
1039 val += 0.5 * hcol[i] * (dot * dot);
1041 if (grad) for (j = 0; j < n; ++j)
1042 grad[j] += dot * xpt[i + j * npt];
1046 if (grad) for (j = 0; j < n; ++j) grad[j] = -grad[j];
1052 static nlopt_result biglag_(int *n, int *npt, double *xopt,
1053 double *xpt, double *bmat, double *zmat, int *idz,
1054 int *ndim, int *knew, double *delta, double *d__,
1055 double *alpha, double *hcol, double *gc, double *gd,
1056 double *s, double *w,
1057 const double *xbase, const double *lb, const double *ub)
1059 /* System generated locals */
1060 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1061 zmat_offset, i__1, i__2;
1064 /* Local variables */
1068 double sp, ss, cf1, cf2, cf3, cf4, cf5, dhd, cth, one, tau, sth, sum,
1071 double zero, angle, scale, denom;
1073 double delsq, tempa, tempb = 0.0, twopi, taubeg, tauold, taumax;
1076 /* N is the number of variables. */
1077 /* NPT is the number of interpolation equations. */
1078 /* XOPT is the best interpolation point so far. */
1079 /* XPT contains the coordinates of the current interpolation points. */
1080 /* BMAT provides the last N columns of H. */
1081 /* ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H. */
1082 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
1083 /* KNEW is the index of the interpolation point that is going to be moved. */
1084 /* DELTA is the current trust region bound. */
1085 /* D will be set to the step from XOPT to the new point. */
1086 /* ALPHA will be set to the KNEW-th diagonal element of the H matrix. */
1087 /* HCOL, GC, GD, S and W will be used for working space. */
1089 /* The step D is calculated in a way that attempts to maximize the modulus */
1090 /* of LFUNC(XOPT+D), subject to the bound ||D|| .LE. DELTA, where LFUNC is */
1091 /* the KNEW-th Lagrange function. */
1093 /* Set some constants. */
1095 /* Parameter adjustments */
1097 zmat_offset = 1 + zmat_dim1;
1098 zmat -= zmat_offset;
1100 xpt_offset = 1 + xpt_dim1;
1104 bmat_offset = 1 + bmat_dim1;
1105 bmat -= bmat_offset;
1117 twopi = atan(one) * 8.;
1118 delsq = *delta * *delta;
1119 nptm = *npt - *n - 1;
1121 /* Set the first NPT components of HCOL to the leading elements of the */
1122 /* KNEW-th column of H. */
1126 for (k = 1; k <= i__1; ++k) {
1131 for (j = 1; j <= i__1; ++j) {
1132 temp = zmat[*knew + j * zmat_dim1];
1137 for (k = 1; k <= i__2; ++k) {
1139 hcol[k] += temp * zmat[k + j * zmat_dim1];
1140 if (nlopt_isinf(hcol[k])) return NLOPT_ROUNDOFF_LIMITED;
1143 *alpha = hcol[*knew];
1145 /* Set the unscaled initial direction D. Form the gradient of LFUNC at */
1146 /* XOPT, and multiply D by the second derivative matrix of LFUNC. */
1150 for (i__ = 1; i__ <= i__2; ++i__) {
1151 d__[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__];
1152 gc[i__] = bmat[*knew + i__ * bmat_dim1];
1155 /* Computing 2nd power */
1160 for (k = 1; k <= i__2; ++k) {
1164 for (j = 1; j <= i__1; ++j) {
1165 temp += xpt[k + j * xpt_dim1] * xopt[j];
1167 sum += xpt[k + j * xpt_dim1] * d__[j];
1169 temp = hcol[k] * temp;
1170 sum = hcol[k] * sum;
1172 for (i__ = 1; i__ <= i__1; ++i__) {
1173 gc[i__] += temp * xpt[k + i__ * xpt_dim1];
1175 gd[i__] += sum * xpt[k + i__ * xpt_dim1];
1179 /* Scale D and GD, with a sign change if required. Set S to another */
1180 /* vector in the initial two dimensional subspace. */
1186 for (i__ = 1; i__ <= i__1; ++i__) {
1187 /* Computing 2nd power */
1190 sp += d__[i__] * gc[i__];
1192 dhd += d__[i__] * gd[i__];
1194 scale = *delta / sqrt(dd);
1195 if (sp * dhd < zero) {
1199 if (sp * sp > dd * .99 * gg) {
1202 tau = scale * (fabs(sp) + half * scale * fabs(dhd));
1203 if (gg * delsq < tau * .01 * tau) {
1207 for (i__ = 1; i__ <= i__1; ++i__) {
1208 d__[i__] = scale * d__[i__];
1209 gd[i__] = scale * gd[i__];
1211 s[i__] = gc[i__] + temp * gd[i__];
1215 double minf, *dlb, *dub, *xtol;
1221 ld.xpt = &xpt[1 + 1 * xpt_dim1];
1222 ld.bmat = &bmat[*knew + 1 * bmat_dim1];
1225 dlb = &gc[1]; dub = &gd[1]; xtol = &s[1];
1226 /* make sure rounding errors don't push initial |d| > delta */
1227 for (j = 1; j <= *n; ++j) d__[j] *= 0.99999;
1228 for (j = 0; j < *n; ++j) {
1229 dlb[j] = -(dub[j] = *delta);
1230 if (dlb[j] < lb[j] - xbase[j] - xopt[j+1])
1231 dlb[j] = lb[j] - xbase[j] - xopt[j+1];
1232 if (dub[j] > ub[j] - xbase[j] - xopt[j+1])
1233 dub[j] = ub[j] - xbase[j] - xopt[j+1];
1234 if (dlb[j] > 0) dlb[j] = 0;
1235 if (dub[j] < 0) dub[j] = 0;
1236 if (d__[j+1] < dlb[j]) d__[j+1] = dlb[j];
1237 else if (d__[j+1] > dub[j]) d__[j+1] = dub[j];
1238 xtol[j] = 1e-5 * *delta;
1240 ld.flipsign = lag(*n, &d__[1], 0, &ld) > 0; /* maximize if > 0 */
1241 return nlopt_minimize_constrained(NLOPT_LD_MMA, *n, lag, &ld,
1242 1, rho_constraint, delta, sizeof(double),
1243 dlb, dub, &d__[1], &minf, -HUGE_VAL,
1244 0., 0., 0., xtol, 1000, 0.);
1247 /* Begin the iteration by overwriting S with a vector that has the */
1248 /* required length and direction, except that termination occurs if */
1249 /* the given D and S are nearly parallel. */
1257 for (i__ = 1; i__ <= i__1; ++i__) {
1258 /* Computing 2nd power */
1261 sp += d__[i__] * s[i__];
1263 /* Computing 2nd power */
1267 temp = dd * ss - sp * sp;
1268 if (temp <= dd * 1e-8 * ss) {
1273 for (i__ = 1; i__ <= i__1; ++i__) {
1274 s[i__] = (dd * s[i__] - sp * d__[i__]) / denom;
1279 /* Calculate the coefficients of the objective function on the circle, */
1280 /* beginning with the multiplication of S by the second derivative matrix. */
1283 for (k = 1; k <= i__1; ++k) {
1286 for (j = 1; j <= i__2; ++j) {
1288 sum += xpt[k + j * xpt_dim1] * s[j];
1290 sum = hcol[k] * sum;
1292 for (i__ = 1; i__ <= i__2; ++i__) {
1294 w[i__] += sum * xpt[k + i__ * xpt_dim1];
1303 for (i__ = 1; i__ <= i__2; ++i__) {
1304 cf1 += s[i__] * w[i__];
1305 cf2 += d__[i__] * gc[i__];
1306 cf3 += s[i__] * gc[i__];
1307 cf4 += d__[i__] * gd[i__];
1309 cf5 += s[i__] * gd[i__];
1312 cf4 = half * cf4 - cf1;
1314 /* Seek the value of the angle that maximizes the modulus of TAU. */
1316 taubeg = cf1 + cf2 + cf4;
1321 temp = twopi / (double) (iu + 1);
1323 for (i__ = 1; i__ <= i__2; ++i__) {
1324 angle = (double) i__ * temp;
1327 tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
1328 if (fabs(tau) > fabs(taumax)) {
1332 } else if (i__ == isave + 1) {
1345 if (tempa != tempb) {
1348 step = half * (tempa - tempb) / (tempa + tempb);
1350 angle = temp * ((double) isave + step);
1352 /* Calculate the new D and GD. Then test for convergence. */
1356 tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
1358 for (i__ = 1; i__ <= i__2; ++i__) {
1359 d__[i__] = cth * d__[i__] + sth * s[i__];
1360 gd[i__] = cth * gd[i__] + sth * w[i__];
1362 s[i__] = gc[i__] + gd[i__];
1364 if (fabs(tau) <= fabs(taubeg) * 1.1) {
1371 return NLOPT_SUCCESS;
1376 /*************************************************************************/
1379 static void update_(int *n, int *npt, double *bmat,
1380 double *zmat, int *idz, int *ndim, double *vlag,
1381 double *beta, int *knew, double *w)
1383 /* System generated locals */
1384 int bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, i__1, i__2;
1387 /* Local variables */
1388 int i__, j, ja, jb, jl, jp;
1389 double one, tau, temp;
1393 double scala, scalb_, alpha, denom, tempa, tempb = 0.0, tausq;
1396 /* The arrays BMAT and ZMAT with IDZ are updated, in order to shift the */
1397 /* interpolation point that has index KNEW. On entry, VLAG contains the */
1398 /* components of the vector Theta*Wcheck+e_b of the updating formula */
1399 /* (6.11), and BETA holds the value of the parameter that has this name. */
1400 /* The vector W is used for working space. */
1402 /* Set some constants. */
1404 /* Parameter adjustments */
1406 zmat_offset = 1 + zmat_dim1;
1407 zmat -= zmat_offset;
1409 bmat_offset = 1 + bmat_dim1;
1410 bmat -= bmat_offset;
1417 nptm = *npt - *n - 1;
1419 /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */
1423 for (j = 2; j <= i__1; ++j) {
1426 } else if (zmat[*knew + j * zmat_dim1] != zero) {
1427 /* Computing 2nd power */
1428 d__1 = zmat[*knew + jl * zmat_dim1];
1429 /* Computing 2nd power */
1430 d__2 = zmat[*knew + j * zmat_dim1];
1431 temp = sqrt(d__1 * d__1 + d__2 * d__2);
1432 tempa = zmat[*knew + jl * zmat_dim1] / temp;
1433 tempb = zmat[*knew + j * zmat_dim1] / temp;
1435 for (i__ = 1; i__ <= i__2; ++i__) {
1436 temp = tempa * zmat[i__ + jl * zmat_dim1] + tempb * zmat[i__
1438 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1]
1439 - tempb * zmat[i__ + jl * zmat_dim1];
1441 zmat[i__ + jl * zmat_dim1] = temp;
1443 zmat[*knew + j * zmat_dim1] = zero;
1448 /* Put the first NPT components of the KNEW-th column of HLAG into W, */
1449 /* and calculate the parameters of the updating formula. */
1451 tempa = zmat[*knew + zmat_dim1];
1456 tempb = zmat[*knew + jl * zmat_dim1];
1459 for (i__ = 1; i__ <= i__1; ++i__) {
1460 w[i__] = tempa * zmat[i__ + zmat_dim1];
1462 w[i__] += tempb * zmat[i__ + jl * zmat_dim1];
1469 denom = alpha * *beta + tausq;
1472 /* Complete the updating of ZMAT when there is only one nonzero element */
1473 /* in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to one, */
1474 /* then the first column of ZMAT will be exchanged with another one later. */
1478 temp = sqrt((fabs(denom)));
1479 tempb = tempa / temp;
1482 for (i__ = 1; i__ <= i__1; ++i__) {
1484 zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb *
1487 if (*idz == 1 && temp < zero) {
1490 if (*idz >= 2 && temp >= zero) {
1495 /* Complete the updating of ZMAT in the alternative case. */
1498 if (*beta >= zero) {
1502 temp = zmat[*knew + jb * zmat_dim1] / denom;
1503 tempa = temp * *beta;
1505 temp = zmat[*knew + ja * zmat_dim1];
1506 scala = one / sqrt(fabs(*beta) * temp * temp + tausq);
1507 scalb_ = scala * sqrt((fabs(denom)));
1509 for (i__ = 1; i__ <= i__1; ++i__) {
1510 zmat[i__ + ja * zmat_dim1] = scala * (tau * zmat[i__ + ja *
1511 zmat_dim1] - temp * vlag[i__]);
1513 zmat[i__ + jb * zmat_dim1] = scalb_ * (zmat[i__ + jb * zmat_dim1]
1514 - tempa * w[i__] - tempb * vlag[i__]);
1516 if (denom <= zero) {
1520 if (*beta >= zero) {
1526 /* IDZ is reduced in the following case, and usually the first column */
1527 /* of ZMAT is exchanged with a later one. */
1532 for (i__ = 1; i__ <= i__1; ++i__) {
1533 temp = zmat[i__ + zmat_dim1];
1534 zmat[i__ + zmat_dim1] = zmat[i__ + *idz * zmat_dim1];
1536 zmat[i__ + *idz * zmat_dim1] = temp;
1540 /* Finally, update the matrix BMAT. */
1543 for (j = 1; j <= i__1; ++j) {
1545 w[jp] = bmat[*knew + j * bmat_dim1];
1546 tempa = (alpha * vlag[jp] - tau * w[jp]) / denom;
1547 tempb = (-(*beta) * w[jp] - tau * vlag[jp]) / denom;
1549 for (i__ = 1; i__ <= i__2; ++i__) {
1550 bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa *
1551 vlag[i__] + tempb * w[i__];
1553 bmat[jp + (i__ - *npt) * bmat_dim1] = bmat[i__ + j *
1563 /*************************************************************************/
1566 static nlopt_result newuob_(int *n, int *npt, double *x,
1568 const double *lb, const double *ub,
1569 nlopt_stopping *stop, double *minf,
1570 newuoa_func calfun, void *calfun_data,
1571 double *xbase, double *xopt, double *xnew,
1572 double *xpt, double *fval, double *gq, double *hq,
1573 double *pq, double *bmat, double *zmat, int *ndim,
1574 double *d__, double *vlag, double *w)
1576 /* System generated locals */
1577 int xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1,
1578 zmat_offset, i__1, i__2, i__3;
1579 double d__1, d__2, d__3;
1581 /* Local variables */
1583 int i__, j, k, ih, nf, nh, ip, jp;
1588 double dsq, rho = 0.0;
1589 int ipt = 0, jpt = 0;
1590 double sum, fbeg = 0.0, diff, half, beta;
1594 double temp, suma, sumb, fopt = HUGE_VAL, bsum, gqsq;
1596 double zero, xipt = 0.0, xjpt = 0.0, sumz, diffa = 0.0, diffb = 0.0, diffc = 0.0, hdiag, alpha = 0.0,
1597 delta, recip, reciq, fsave;
1598 int ksave, nfsav = 0, itemp;
1599 double dnorm = 0.0, ratio = 0.0, dstep, tenth, vquad;
1604 double detrat, crvmin = 0.0;
1606 double xoptsq = 0.0;
1608 nlopt_result rc = NLOPT_SUCCESS, rc2;
1610 /* SGJ, 2008: compute rhoend from NLopt stop info */
1611 rhoend = stop->xtol_rel * (*rhobeg);
1612 for (j = 0; j < *n; ++j)
1613 if (rhoend < stop->xtol_abs[j])
1614 rhoend = stop->xtol_abs[j];
1616 /* The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are identical */
1617 /* to the corresponding arguments in SUBROUTINE NEWUOA. */
1618 /* XBASE will hold a shift of origin that should reduce the contributions */
1619 /* from rounding errors to values of the model and Lagrange functions. */
1620 /* XOPT will be set to the displacement from XBASE of the vector of */
1621 /* variables that provides the least calculated F so far. */
1622 /* XNEW will be set to the displacement from XBASE of the vector of */
1623 /* variables for the current calculation of F. */
1624 /* XPT will contain the interpolation point coordinates relative to XBASE. */
1625 /* FVAL will hold the values of F at the interpolation points. */
1626 /* GQ will hold the gradient of the quadratic model at XBASE. */
1627 /* HQ will hold the explicit second derivatives of the quadratic model. */
1628 /* PQ will contain the parameters of the implicit second derivatives of */
1629 /* the quadratic model. */
1630 /* BMAT will hold the last N columns of H. */
1631 /* ZMAT will hold the factorization of the leading NPT by NPT submatrix of */
1632 /* H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, where */
1633 /* the elements of DZ are plus or minus one, as specified by IDZ. */
1634 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
1635 /* D is reserved for trial steps from XOPT. */
1636 /* VLAG will contain the values of the Lagrange functions at a new point X. */
1637 /* They are part of a product that requires VLAG to be of length NDIM. */
1638 /* The array W will be used for working space. Its length must be at least */
1639 /* 10*NDIM = 10*(NPT+N). */
1641 /* Set some constants. */
1643 /* Parameter adjustments */
1645 zmat_offset = 1 + zmat_dim1;
1646 zmat -= zmat_offset;
1648 xpt_offset = 1 + xpt_dim1;
1659 bmat_offset = 1 + bmat_dim1;
1660 bmat -= bmat_offset;
1674 /* Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */
1677 for (j = 1; j <= i__1; ++j) {
1680 for (k = 1; k <= i__2; ++k) {
1682 xpt[k + j * xpt_dim1] = zero;
1685 for (i__ = 1; i__ <= i__2; ++i__) {
1687 bmat[i__ + j * bmat_dim1] = zero;
1691 for (ih = 1; ih <= i__2; ++ih) {
1696 for (k = 1; k <= i__2; ++k) {
1699 for (j = 1; j <= i__1; ++j) {
1701 zmat[k + j * zmat_dim1] = zero;
1705 /* Begin the initialization procedure. NF becomes one more than the number */
1706 /* of function values so far. The coordinates of the displacement of the */
1707 /* next initial interpolation point from XBASE are set in XPT(NF,.). */
1709 rhosq = *rhobeg * *rhobeg;
1710 recip = one / rhosq;
1711 reciq = sqrt(half) / rhosq;
1717 if (nfm <= *n << 1) {
1718 if (nfm >= 1 && nfm <= *n) {
1719 xpt[nf + nfm * xpt_dim1] = *rhobeg;
1720 } else if (nfm > *n) {
1721 xpt[nf + nfmm * xpt_dim1] = -(*rhobeg);
1724 itemp = (nfmm - 1) / *n;
1725 jpt = nfm - itemp * *n - *n;
1733 if (fval[ipt + np] < fval[ipt + 1]) {
1737 if (fval[jpt + np] < fval[jpt + 1]) {
1740 xpt[nf + ipt * xpt_dim1] = xipt;
1741 xpt[nf + jpt * xpt_dim1] = xjpt;
1744 /* Calculate the next value of F, label 70 being reached immediately */
1745 /* after this calculation. The least function value so far and its index */
1749 for (j = 1; j <= i__1; ++j) {
1751 x[j] = xpt[nf + j * xpt_dim1] + xbase[j];
1753 if (lb && ub) { /* SGJ, 2008: make sure we are within bounds */
1754 for (j = 1; j <= i__1; ++j) {
1755 if (x[j] < lb[j-1]) x[j] = lb[j-1];
1756 else if (x[j] > ub[j-1]) x[j] = ub[j-1];
1766 } else if (f < fopt) {
1771 /* Set the nonzero initial elements of BMAT and the quadratic model in */
1772 /* the cases when NF is at most 2*N+1. */
1774 if (nfm <= *n << 1) {
1775 if (nfm >= 1 && nfm <= *n) {
1776 gq[nfm] = (f - fbeg) / *rhobeg;
1777 if (*npt < nf + *n) {
1778 bmat[nfm * bmat_dim1 + 1] = -one / *rhobeg;
1779 bmat[nf + nfm * bmat_dim1] = one / *rhobeg;
1780 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
1782 } else if (nfm > *n) {
1783 bmat[nf - *n + nfmm * bmat_dim1] = half / *rhobeg;
1784 bmat[nf + nfmm * bmat_dim1] = -half / *rhobeg;
1785 zmat[nfmm * zmat_dim1 + 1] = -reciq - reciq;
1786 zmat[nf - *n + nfmm * zmat_dim1] = reciq;
1787 zmat[nf + nfmm * zmat_dim1] = reciq;
1788 ih = nfmm * (nfmm + 1) / 2;
1789 temp = (fbeg - f) / *rhobeg;
1790 hq[ih] = (gq[nfmm] - temp) / *rhobeg;
1791 gq[nfmm] = half * (gq[nfmm] + temp);
1794 /* Set the off-diagonal second derivatives of the Lagrange functions and */
1795 /* the initial quadratic model. */
1798 ih = ipt * (ipt - 1) / 2 + jpt;
1805 zmat[nfmm * zmat_dim1 + 1] = recip;
1806 zmat[nf + nfmm * zmat_dim1] = recip;
1807 zmat[ipt + 1 + nfmm * zmat_dim1] = -recip;
1808 zmat[jpt + 1 + nfmm * zmat_dim1] = -recip;
1809 hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / (xipt * xjpt);
1815 /* Begin the iterative procedure, because the initial model is complete. */
1825 for (i__ = 1; i__ <= i__1; ++i__) {
1826 xopt[i__] = xpt[kopt + i__ * xpt_dim1];
1828 /* Computing 2nd power */
1830 xoptsq += d__1 * d__1;
1835 /* Generate the next trust region step and test its length. Set KNEW */
1836 /* to -1 if the purpose of the next F will be to improve the model. */
1840 rc2 = trsapp_(n, npt, &xopt[1], &xpt[xpt_offset], &gq[1], &hq[1], &pq[1], &
1841 delta, &d__[1], &w[1], &w[np], &w[np + *n], &w[np + (*n << 1)], &
1842 crvmin, &xbase[1], lb, ub);
1843 if (rc2 < 0) { rc = rc2; goto L530; }
1846 for (i__ = 1; i__ <= i__1; ++i__) {
1848 /* Computing 2nd power */
1853 d__1 = delta, d__2 = sqrt(dsq);
1854 dnorm = MIN2(d__1,d__2);
1855 if (dnorm < half * rho) {
1857 delta = tenth * delta;
1859 if (delta <= rho * 1.5) {
1862 if (nf <= nfsav + 2) {
1865 temp = crvmin * .125 * rho * rho;
1867 d__1 = MAX2(diffa,diffb);
1868 if (temp <= MAX2(d__1,diffc)) {
1874 /* Shift XBASE if XOPT may be too far from XBASE. First make the changes */
1875 /* to BMAT that do not depend on ZMAT. */
1878 if (dsq <= xoptsq * .001) {
1879 tempq = xoptsq * .25;
1881 for (k = 1; k <= i__1; ++k) {
1884 for (i__ = 1; i__ <= i__2; ++i__) {
1886 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
1889 sum -= half * xoptsq;
1892 for (i__ = 1; i__ <= i__2; ++i__) {
1893 gq[i__] += temp * xpt[k + i__ * xpt_dim1];
1894 xpt[k + i__ * xpt_dim1] -= half * xopt[i__];
1895 vlag[i__] = bmat[k + i__ * bmat_dim1];
1896 w[i__] = sum * xpt[k + i__ * xpt_dim1] + tempq * xopt[i__];
1899 for (j = 1; j <= i__3; ++j) {
1901 bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] +
1902 vlag[i__] * w[j] + w[i__] * vlag[j];
1907 /* Then the revisions of BMAT that depend on ZMAT are calculated. */
1910 for (k = 1; k <= i__3; ++k) {
1913 for (i__ = 1; i__ <= i__2; ++i__) {
1914 sumz += zmat[i__ + k * zmat_dim1];
1916 w[i__] = w[*npt + i__] * zmat[i__ + k * zmat_dim1];
1919 for (j = 1; j <= i__2; ++j) {
1920 sum = tempq * sumz * xopt[j];
1922 for (i__ = 1; i__ <= i__1; ++i__) {
1924 sum += w[i__] * xpt[i__ + j * xpt_dim1];
1931 for (i__ = 1; i__ <= i__1; ++i__) {
1933 bmat[i__ + j * bmat_dim1] += sum * zmat[i__ + k *
1938 for (i__ = 1; i__ <= i__1; ++i__) {
1945 for (j = 1; j <= i__2; ++j) {
1947 bmat[ip + j * bmat_dim1] += temp * vlag[j];
1952 /* The following instructions complete the shift of XBASE, including */
1953 /* the changes to the parameters of the quadratic model. */
1957 for (j = 1; j <= i__2; ++j) {
1960 for (k = 1; k <= i__1; ++k) {
1961 w[j] += pq[k] * xpt[k + j * xpt_dim1];
1963 xpt[k + j * xpt_dim1] -= half * xopt[j];
1966 for (i__ = 1; i__ <= i__1; ++i__) {
1969 gq[j] += hq[ih] * xopt[i__];
1971 gq[i__] += hq[ih] * xopt[j];
1972 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
1974 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
1979 for (j = 1; j <= i__1; ++j) {
1980 xbase[j] += xopt[j];
1987 /* Pick the model step if KNEW is positive. A different choice of D */
1988 /* may be made later, if the choice of D by BIGLAG causes substantial */
1989 /* cancellation in DENOM. */
1993 biglag_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &zmat[
1994 zmat_offset], &idz, ndim, &knew, &dstep, &d__[1], &alpha, &
1995 vlag[1], &vlag[*npt + 1], &w[1], &w[np], &w[np + *n],
1997 if (rc2 < 0) { rc = rc2; goto L530; }
2000 /* Calculate VLAG and BETA for the current choice of D. The first NPT */
2001 /* components of W_check will be held in W. */
2004 for (k = 1; k <= i__1; ++k) {
2009 for (j = 1; j <= i__2; ++j) {
2010 suma += xpt[k + j * xpt_dim1] * d__[j];
2011 sumb += xpt[k + j * xpt_dim1] * xopt[j];
2013 sum += bmat[k + j * bmat_dim1] * d__[j];
2015 w[k] = suma * (half * suma + sumb);
2021 for (k = 1; k <= i__1; ++k) {
2024 for (i__ = 1; i__ <= i__2; ++i__) {
2026 sum += zmat[i__ + k * zmat_dim1] * w[i__];
2035 for (i__ = 1; i__ <= i__2; ++i__) {
2037 vlag[i__] += sum * zmat[i__ + k * zmat_dim1];
2043 for (j = 1; j <= i__2; ++j) {
2046 for (i__ = 1; i__ <= i__1; ++i__) {
2048 sum += w[i__] * bmat[i__ + j * bmat_dim1];
2050 bsum += sum * d__[j];
2053 for (k = 1; k <= i__1; ++k) {
2055 sum += bmat[jp + k * bmat_dim1] * d__[k];
2058 bsum += sum * d__[j];
2060 dx += d__[j] * xopt[j];
2062 beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2065 /* If KNEW is positive and if the cancellation in DENOM is unacceptable, */
2066 /* then BIGDEN calculates an alternative model step, XNEW being used for */
2067 /* working space. */
2070 /* Computing 2nd power */
2072 if (d__1 == 0) { rc = NLOPT_ROUNDOFF_LIMITED; goto L530; }
2073 temp = one + alpha * beta / (d__1 * d__1);
2074 if (fabs(temp) <= .8) {
2076 bigden_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &
2077 zmat[zmat_offset], &idz, ndim, &kopt, &knew, &d__[1], &w[
2078 1], &vlag[1], &beta, &xnew[1], &w[*ndim + 1], &w[*ndim *
2079 6 + 1], &xbase[1], lb, ub);
2080 if (rc2 < 0) { rc = rc2; goto L530; }
2084 /* Calculate the next value of the objective function. */
2088 for (i__ = 1; i__ <= i__2; ++i__) {
2089 xnew[i__] = xopt[i__] + d__[i__];
2091 x[i__] = xbase[i__] + xnew[i__];
2093 if (lb && ub) { /* SGJ, 2008: make sure we are within bounds,
2094 since roundoff errors can push us slightly outside */
2095 for (j = 1; j <= i__1; ++j) {
2096 if (x[j] < lb[j-1]) x[j] = lb[j-1];
2097 else if (x[j] > ub[j-1]) x[j] = ub[j-1];
2102 if (nlopt_stop_forced(stop)) rc = NLOPT_FORCED_STOP;
2103 else if (stop->nevals > 0) {
2104 if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
2105 else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
2107 if (rc != NLOPT_SUCCESS) goto L530;
2110 f = calfun(*n, &x[1], calfun_data);
2111 if (f < stop->minf_max) {
2112 rc = NLOPT_MINF_MAX_REACHED;
2123 /* Use the quadratic model to predict the change in F due to the step D, */
2124 /* and set DIFF to the error of this prediction. */
2129 for (j = 1; j <= i__2; ++j) {
2130 vquad += d__[j] * gq[j];
2132 for (i__ = 1; i__ <= i__1; ++i__) {
2134 temp = d__[i__] * xnew[j] + d__[j] * xopt[i__];
2139 vquad += temp * hq[ih];
2143 for (k = 1; k <= i__1; ++k) {
2145 vquad += pq[k] * w[k];
2147 diff = f - fopt - vquad;
2155 /* Update FOPT and XOPT if the new F is the least value of the objective */
2156 /* function so far. The branch when KNEW is positive occurs if D is not */
2157 /* a trust region step. */
2164 for (i__ = 1; i__ <= i__1; ++i__) {
2165 xopt[i__] = xnew[i__];
2167 /* Computing 2nd power */
2169 xoptsq += d__1 * d__1;
2171 if (nlopt_stop_ftol(stop, fopt, fsave)) {
2172 rc = NLOPT_FTOL_REACHED;
2182 /* Pick the next value of DELTA after a trust region step. */
2184 if (vquad >= zero) {
2187 ratio = (f - fsave) / vquad;
2188 if (ratio <= tenth) {
2189 delta = half * dnorm;
2190 } else if (ratio <= .7) {
2192 d__1 = half * delta;
2193 delta = MAX2(d__1,dnorm);
2196 d__1 = half * delta, d__2 = dnorm + dnorm;
2197 delta = MAX2(d__1,d__2);
2199 if (delta <= rho * 1.5) {
2203 /* Set KNEW to the index of the next interpolation point to be deleted. */
2206 d__2 = tenth * delta;
2207 /* Computing 2nd power */
2208 d__1 = MAX2(d__2,rho);
2209 rhosq = d__1 * d__1;
2217 for (k = 1; k <= i__1; ++k) {
2220 for (j = 1; j <= i__2; ++j) {
2226 /* Computing 2nd power */
2227 d__1 = zmat[k + j * zmat_dim1];
2228 hdiag += temp * (d__1 * d__1);
2230 /* Computing 2nd power */
2232 temp = (d__1 = beta * hdiag + d__2 * d__2, fabs(d__1));
2235 for (j = 1; j <= i__2; ++j) {
2237 /* Computing 2nd power */
2238 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2239 distsq += d__1 * d__1;
2241 if (distsq > rhosq) {
2242 /* Computing 3rd power */
2243 d__1 = distsq / rhosq;
2244 temp *= d__1 * (d__1 * d__1);
2246 if (temp > detrat && k != ktemp) {
2256 /* Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point */
2257 /* can be moved. Begin the updating of the quadratic model, starting */
2258 /* with the explicit second derivative term. */
2261 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], &idz, ndim, &vlag[
2262 1], &beta, &knew, &w[1]);
2266 for (i__ = 1; i__ <= i__1; ++i__) {
2267 temp = pq[knew] * xpt[knew + i__ * xpt_dim1];
2269 for (j = 1; j <= i__2; ++j) {
2272 hq[ih] += temp * xpt[knew + j * xpt_dim1];
2277 /* Update the other second derivative parameters, and then the gradient */
2278 /* vector of the model. Also include the new interpolation point. */
2281 for (j = 1; j <= i__2; ++j) {
2282 temp = diff * zmat[knew + j * zmat_dim1];
2287 for (k = 1; k <= i__1; ++k) {
2289 pq[k] += temp * zmat[k + j * zmat_dim1];
2294 for (i__ = 1; i__ <= i__1; ++i__) {
2295 gq[i__] += diff * bmat[knew + i__ * bmat_dim1];
2296 /* Computing 2nd power */
2298 gqsq += d__1 * d__1;
2300 xpt[knew + i__ * xpt_dim1] = xnew[i__];
2303 /* If a trust region step makes a small change to the objective function, */
2304 /* then calculate the gradient of the least Frobenius norm interpolant at */
2305 /* XBASE, and store it in W, using VLAG for a vector of right hand sides. */
2307 if (ksave == 0 && delta == rho) {
2308 if (fabs(ratio) > .01) {
2312 for (k = 1; k <= i__1; ++k) {
2314 vlag[k] = fval[k] - fval[kopt];
2318 for (i__ = 1; i__ <= i__1; ++i__) {
2321 for (k = 1; k <= i__2; ++k) {
2323 sum += bmat[k + i__ * bmat_dim1] * vlag[k];
2330 /* Test whether to replace the new quadratic model by the least Frobenius */
2331 /* norm interpolant, making the replacement if the test is satisfied. */
2334 if (gqsq < gisq * 100.) {
2339 for (i__ = 1; i__ <= i__1; ++i__) {
2344 for (ih = 1; ih <= i__1; ++ih) {
2349 for (j = 1; j <= i__1; ++j) {
2352 for (k = 1; k <= i__2; ++k) {
2354 w[j] += vlag[k] * zmat[k + j * zmat_dim1];
2362 for (k = 1; k <= i__1; ++k) {
2365 for (j = 1; j <= i__2; ++j) {
2367 pq[k] += zmat[k + j * zmat_dim1] * w[j];
2378 /* If a trust region step has provided a sufficient decrease in F, then */
2379 /* branch for another trust region calculation. The case KSAVE>0 occurs */
2380 /* when the new function value was calculated by a model step. */
2382 if (f <= fsave + tenth * vquad) {
2389 /* Alternatively, find out if the interpolation points are close enough */
2390 /* to the best point so far. */
2394 distsq = delta * 4. * delta;
2396 for (k = 1; k <= i__2; ++k) {
2399 for (j = 1; j <= i__1; ++j) {
2401 /* Computing 2nd power */
2402 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2412 /* If KNEW is positive, then set DSTEP, and branch back for the next */
2413 /* iteration, which will generate a "model step". */
2418 d__2 = tenth * sqrt(distsq), d__3 = half * delta;
2419 d__1 = MIN2(d__2,d__3);
2420 dstep = MAX2(d__1,rho);
2421 dsq = dstep * dstep;
2427 if (MAX2(delta,dnorm) > rho) {
2431 /* The calculations with the current value of RHO are complete. Pick the */
2432 /* next values of RHO and DELTA. */
2437 ratio = rho / rhoend;
2440 } else if (ratio <= 250.) {
2441 rho = sqrt(ratio) * rhoend;
2445 delta = MAX2(delta,rho);
2449 /* Return from the calculation, after another Newton-Raphson step, if */
2450 /* it is too short to have been tried before. */
2455 rc = NLOPT_XTOL_REACHED;
2459 for (i__ = 1; i__ <= i__2; ++i__) {
2461 x[i__] = xbase[i__] + xopt[i__];
2469 /*************************************************************************/
2472 nlopt_result newuoa(int n, int npt, double *x,
2473 const double *lb, const double *ub,
2474 double rhobeg, nlopt_stopping *stop, double *minf,
2475 newuoa_func calfun, void *calfun_data)
2477 /* Local variables */
2478 int id, np, iw, igq, ihq, ixb, ifv, ipq, ivl, ixn, ixo, ixp, ndim,
2483 /* This subroutine seeks the least value of a function of many variables, */
2484 /* by a trust region method that forms quadratic models by interpolation. */
2485 /* There can be some freedom in the interpolation conditions, which is */
2486 /* taken up by minimizing the Frobenius norm of the change to the second */
2487 /* derivative of the quadratic model, beginning with a zero matrix. The */
2488 /* arguments of the subroutine are as follows. */
2490 /* N must be set to the number of variables and must be at least two. */
2491 /* NPT is the number of interpolation conditions. Its value must be in the */
2492 /* interval [N+2,(N+1)(N+2)/2]. */
2493 /* Initial values of the variables must be set in X(1),X(2),...,X(N). They */
2494 /* will be changed to the values that give the least calculated F. */
2495 /* RHOBEG and RHOEND must be set to the initial and final values of a trust */
2496 /* region radius, so both must be positive with RHOEND<=RHOBEG. Typically */
2497 /* RHOBEG should be about one tenth of the greatest expected change to a */
2498 /* variable, and RHOEND should indicate the accuracy that is required in */
2499 /* the final values of the variables. */
2500 /* The value of IPRINT should be set to 0, 1, 2 or 3, which controls the */
2501 /* amount of printing. Specifically, there is no output if IPRINT=0 and */
2502 /* there is output only at the return if IPRINT=1. Otherwise, each new */
2503 /* value of RHO is printed, with the best vector of variables so far and */
2504 /* the corresponding value of the objective function. Further, each new */
2505 /* value of F with its variables are output if IPRINT=3. */
2506 /* MAXFUN must be set to an upper bound on the number of calls of CALFUN. */
2507 /* The array W will be used for working space. Its length must be at least */
2508 /* (NPT+13)*(NPT+N)+3*N*(N+3)/2. */
2510 /* SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set F to */
2511 /* the value of the objective function for the variables X(1),X(2),...,X(N). */
2513 /* Partition the working space array, so that different parts of it can be */
2514 /* treated separately by the subroutine that performs the main calculation. */
2516 /* Parameter adjustments */
2523 nlopt_stop_msg(stop, "dimension %d must be >= 2", n);
2524 return NLOPT_INVALID_ARGS;
2526 if (npt < n + 2 || npt > (n + 2) * np / 2) {
2527 nlopt_stop_msg(stop, "invalid # of interpolation conditions %d", npt);
2528 return NLOPT_INVALID_ARGS;
2535 ifv = ixp + n * npt;
2538 ipq = ihq + n * np / 2;
2540 izmat = ibmat + ndim * n;
2541 id = izmat + npt * nptm;
2545 w = (double *) malloc(sizeof(double) * ((npt+13)*(npt+n) + 3*(n*(n+3))/2));
2546 if (!w) return NLOPT_OUT_OF_MEMORY;
2549 /* The above settings provide a partition of W for subroutine NEWUOB. */
2550 /* The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements of */
2551 /* W plus the space that is needed by the last array of NEWUOB. */
2553 ret = newuob_(&n, &npt, &x[1], &rhobeg,
2554 lb, ub, stop, minf, calfun, calfun_data,
2555 &w[ixb], &w[ixo], &w[ixn], &w[ixp], &w[ifv],
2556 &w[igq], &w[ihq], &w[ipq], &w[ibmat], &w[izmat],
2557 &ndim, &w[id], &w[ivl], &w[iw]);
2564 /*************************************************************************/
2565 /*************************************************************************/