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];
1754 } else if (f < fopt) {
1759 /* Set the nonzero initial elements of BMAT and the quadratic model in */
1760 /* the cases when NF is at most 2*N+1. */
1762 if (nfm <= *n << 1) {
1763 if (nfm >= 1 && nfm <= *n) {
1764 gq[nfm] = (f - fbeg) / *rhobeg;
1765 if (*npt < nf + *n) {
1766 bmat[nfm * bmat_dim1 + 1] = -one / *rhobeg;
1767 bmat[nf + nfm * bmat_dim1] = one / *rhobeg;
1768 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
1770 } else if (nfm > *n) {
1771 bmat[nf - *n + nfmm * bmat_dim1] = half / *rhobeg;
1772 bmat[nf + nfmm * bmat_dim1] = -half / *rhobeg;
1773 zmat[nfmm * zmat_dim1 + 1] = -reciq - reciq;
1774 zmat[nf - *n + nfmm * zmat_dim1] = reciq;
1775 zmat[nf + nfmm * zmat_dim1] = reciq;
1776 ih = nfmm * (nfmm + 1) / 2;
1777 temp = (fbeg - f) / *rhobeg;
1778 hq[ih] = (gq[nfmm] - temp) / *rhobeg;
1779 gq[nfmm] = half * (gq[nfmm] + temp);
1782 /* Set the off-diagonal second derivatives of the Lagrange functions and */
1783 /* the initial quadratic model. */
1786 ih = ipt * (ipt - 1) / 2 + jpt;
1793 zmat[nfmm * zmat_dim1 + 1] = recip;
1794 zmat[nf + nfmm * zmat_dim1] = recip;
1795 zmat[ipt + 1 + nfmm * zmat_dim1] = -recip;
1796 zmat[jpt + 1 + nfmm * zmat_dim1] = -recip;
1797 hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / (xipt * xjpt);
1803 /* Begin the iterative procedure, because the initial model is complete. */
1813 for (i__ = 1; i__ <= i__1; ++i__) {
1814 xopt[i__] = xpt[kopt + i__ * xpt_dim1];
1816 /* Computing 2nd power */
1818 xoptsq += d__1 * d__1;
1823 /* Generate the next trust region step and test its length. Set KNEW */
1824 /* to -1 if the purpose of the next F will be to improve the model. */
1828 rc2 = trsapp_(n, npt, &xopt[1], &xpt[xpt_offset], &gq[1], &hq[1], &pq[1], &
1829 delta, &d__[1], &w[1], &w[np], &w[np + *n], &w[np + (*n << 1)], &
1830 crvmin, &xbase[1], lb, ub);
1831 if (rc2 < 0) { rc = rc2; goto L530; }
1834 for (i__ = 1; i__ <= i__1; ++i__) {
1836 /* Computing 2nd power */
1841 d__1 = delta, d__2 = sqrt(dsq);
1842 dnorm = min(d__1,d__2);
1843 if (dnorm < half * rho) {
1845 delta = tenth * delta;
1847 if (delta <= rho * 1.5) {
1850 if (nf <= nfsav + 2) {
1853 temp = crvmin * .125 * rho * rho;
1855 d__1 = max(diffa,diffb);
1856 if (temp <= max(d__1,diffc)) {
1862 /* Shift XBASE if XOPT may be too far from XBASE. First make the changes */
1863 /* to BMAT that do not depend on ZMAT. */
1866 if (dsq <= xoptsq * .001) {
1867 tempq = xoptsq * .25;
1869 for (k = 1; k <= i__1; ++k) {
1872 for (i__ = 1; i__ <= i__2; ++i__) {
1874 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
1877 sum -= half * xoptsq;
1880 for (i__ = 1; i__ <= i__2; ++i__) {
1881 gq[i__] += temp * xpt[k + i__ * xpt_dim1];
1882 xpt[k + i__ * xpt_dim1] -= half * xopt[i__];
1883 vlag[i__] = bmat[k + i__ * bmat_dim1];
1884 w[i__] = sum * xpt[k + i__ * xpt_dim1] + tempq * xopt[i__];
1887 for (j = 1; j <= i__3; ++j) {
1889 bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] +
1890 vlag[i__] * w[j] + w[i__] * vlag[j];
1895 /* Then the revisions of BMAT that depend on ZMAT are calculated. */
1898 for (k = 1; k <= i__3; ++k) {
1901 for (i__ = 1; i__ <= i__2; ++i__) {
1902 sumz += zmat[i__ + k * zmat_dim1];
1904 w[i__] = w[*npt + i__] * zmat[i__ + k * zmat_dim1];
1907 for (j = 1; j <= i__2; ++j) {
1908 sum = tempq * sumz * xopt[j];
1910 for (i__ = 1; i__ <= i__1; ++i__) {
1912 sum += w[i__] * xpt[i__ + j * xpt_dim1];
1919 for (i__ = 1; i__ <= i__1; ++i__) {
1921 bmat[i__ + j * bmat_dim1] += sum * zmat[i__ + k *
1926 for (i__ = 1; i__ <= i__1; ++i__) {
1933 for (j = 1; j <= i__2; ++j) {
1935 bmat[ip + j * bmat_dim1] += temp * vlag[j];
1940 /* The following instructions complete the shift of XBASE, including */
1941 /* the changes to the parameters of the quadratic model. */
1945 for (j = 1; j <= i__2; ++j) {
1948 for (k = 1; k <= i__1; ++k) {
1949 w[j] += pq[k] * xpt[k + j * xpt_dim1];
1951 xpt[k + j * xpt_dim1] -= half * xopt[j];
1954 for (i__ = 1; i__ <= i__1; ++i__) {
1957 gq[j] += hq[ih] * xopt[i__];
1959 gq[i__] += hq[ih] * xopt[j];
1960 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
1962 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
1967 for (j = 1; j <= i__1; ++j) {
1968 xbase[j] += xopt[j];
1975 /* Pick the model step if KNEW is positive. A different choice of D */
1976 /* may be made later, if the choice of D by BIGLAG causes substantial */
1977 /* cancellation in DENOM. */
1981 biglag_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &zmat[
1982 zmat_offset], &idz, ndim, &knew, &dstep, &d__[1], &alpha, &
1983 vlag[1], &vlag[*npt + 1], &w[1], &w[np], &w[np + *n],
1985 if (rc2 < 0) { rc = rc2; goto L530; }
1988 /* Calculate VLAG and BETA for the current choice of D. The first NPT */
1989 /* components of W_check will be held in W. */
1992 for (k = 1; k <= i__1; ++k) {
1997 for (j = 1; j <= i__2; ++j) {
1998 suma += xpt[k + j * xpt_dim1] * d__[j];
1999 sumb += xpt[k + j * xpt_dim1] * xopt[j];
2001 sum += bmat[k + j * bmat_dim1] * d__[j];
2003 w[k] = suma * (half * suma + sumb);
2009 for (k = 1; k <= i__1; ++k) {
2012 for (i__ = 1; i__ <= i__2; ++i__) {
2014 sum += zmat[i__ + k * zmat_dim1] * w[i__];
2023 for (i__ = 1; i__ <= i__2; ++i__) {
2025 vlag[i__] += sum * zmat[i__ + k * zmat_dim1];
2031 for (j = 1; j <= i__2; ++j) {
2034 for (i__ = 1; i__ <= i__1; ++i__) {
2036 sum += w[i__] * bmat[i__ + j * bmat_dim1];
2038 bsum += sum * d__[j];
2041 for (k = 1; k <= i__1; ++k) {
2043 sum += bmat[jp + k * bmat_dim1] * d__[k];
2046 bsum += sum * d__[j];
2048 dx += d__[j] * xopt[j];
2050 beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
2053 /* If KNEW is positive and if the cancellation in DENOM is unacceptable, */
2054 /* then BIGDEN calculates an alternative model step, XNEW being used for */
2055 /* working space. */
2058 /* Computing 2nd power */
2060 temp = one + alpha * beta / (d__1 * d__1);
2061 if (fabs(temp) <= .8) {
2062 bigden_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &
2063 zmat[zmat_offset], &idz, ndim, &kopt, &knew, &d__[1], &w[
2064 1], &vlag[1], &beta, &xnew[1], &w[*ndim + 1], &w[*ndim *
2065 6 + 1], &xbase[1], lb, ub);
2069 /* Calculate the next value of the objective function. */
2073 for (i__ = 1; i__ <= i__2; ++i__) {
2074 xnew[i__] = xopt[i__] + d__[i__];
2076 x[i__] = xbase[i__] + xnew[i__];
2080 if (nlopt_stop_evals(stop)) rc = NLOPT_MAXEVAL_REACHED;
2081 else if (nlopt_stop_time(stop)) rc = NLOPT_MAXTIME_REACHED;
2082 if (rc != NLOPT_SUCCESS) goto L530;
2085 f = calfun(*n, &x[1], calfun_data);
2086 if (f < stop->minf_max) {
2087 rc = NLOPT_MINF_MAX_REACHED;
2098 /* Use the quadratic model to predict the change in F due to the step D, */
2099 /* and set DIFF to the error of this prediction. */
2104 for (j = 1; j <= i__2; ++j) {
2105 vquad += d__[j] * gq[j];
2107 for (i__ = 1; i__ <= i__1; ++i__) {
2109 temp = d__[i__] * xnew[j] + d__[j] * xopt[i__];
2114 vquad += temp * hq[ih];
2118 for (k = 1; k <= i__1; ++k) {
2120 vquad += pq[k] * w[k];
2122 diff = f - fopt - vquad;
2130 /* Update FOPT and XOPT if the new F is the least value of the objective */
2131 /* function so far. The branch when KNEW is positive occurs if D is not */
2132 /* a trust region step. */
2139 for (i__ = 1; i__ <= i__1; ++i__) {
2140 xopt[i__] = xnew[i__];
2142 /* Computing 2nd power */
2144 xoptsq += d__1 * d__1;
2146 if (nlopt_stop_ftol(stop, fopt, fsave)) {
2147 rc = NLOPT_FTOL_REACHED;
2157 /* Pick the next value of DELTA after a trust region step. */
2159 if (vquad >= zero) {
2162 ratio = (f - fsave) / vquad;
2163 if (ratio <= tenth) {
2164 delta = half * dnorm;
2165 } else if (ratio <= .7) {
2167 d__1 = half * delta;
2168 delta = max(d__1,dnorm);
2171 d__1 = half * delta, d__2 = dnorm + dnorm;
2172 delta = max(d__1,d__2);
2174 if (delta <= rho * 1.5) {
2178 /* Set KNEW to the index of the next interpolation point to be deleted. */
2181 d__2 = tenth * delta;
2182 /* Computing 2nd power */
2183 d__1 = max(d__2,rho);
2184 rhosq = d__1 * d__1;
2192 for (k = 1; k <= i__1; ++k) {
2195 for (j = 1; j <= i__2; ++j) {
2201 /* Computing 2nd power */
2202 d__1 = zmat[k + j * zmat_dim1];
2203 hdiag += temp * (d__1 * d__1);
2205 /* Computing 2nd power */
2207 temp = (d__1 = beta * hdiag + d__2 * d__2, fabs(d__1));
2210 for (j = 1; j <= i__2; ++j) {
2212 /* Computing 2nd power */
2213 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2214 distsq += d__1 * d__1;
2216 if (distsq > rhosq) {
2217 /* Computing 3rd power */
2218 d__1 = distsq / rhosq;
2219 temp *= d__1 * (d__1 * d__1);
2221 if (temp > detrat && k != ktemp) {
2231 /* Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point */
2232 /* can be moved. Begin the updating of the quadratic model, starting */
2233 /* with the explicit second derivative term. */
2236 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], &idz, ndim, &vlag[
2237 1], &beta, &knew, &w[1]);
2241 for (i__ = 1; i__ <= i__1; ++i__) {
2242 temp = pq[knew] * xpt[knew + i__ * xpt_dim1];
2244 for (j = 1; j <= i__2; ++j) {
2247 hq[ih] += temp * xpt[knew + j * xpt_dim1];
2252 /* Update the other second derivative parameters, and then the gradient */
2253 /* vector of the model. Also include the new interpolation point. */
2256 for (j = 1; j <= i__2; ++j) {
2257 temp = diff * zmat[knew + j * zmat_dim1];
2262 for (k = 1; k <= i__1; ++k) {
2264 pq[k] += temp * zmat[k + j * zmat_dim1];
2269 for (i__ = 1; i__ <= i__1; ++i__) {
2270 gq[i__] += diff * bmat[knew + i__ * bmat_dim1];
2271 /* Computing 2nd power */
2273 gqsq += d__1 * d__1;
2275 xpt[knew + i__ * xpt_dim1] = xnew[i__];
2278 /* If a trust region step makes a small change to the objective function, */
2279 /* then calculate the gradient of the least Frobenius norm interpolant at */
2280 /* XBASE, and store it in W, using VLAG for a vector of right hand sides. */
2282 if (ksave == 0 && delta == rho) {
2283 if (fabs(ratio) > .01) {
2287 for (k = 1; k <= i__1; ++k) {
2289 vlag[k] = fval[k] - fval[kopt];
2293 for (i__ = 1; i__ <= i__1; ++i__) {
2296 for (k = 1; k <= i__2; ++k) {
2298 sum += bmat[k + i__ * bmat_dim1] * vlag[k];
2305 /* Test whether to replace the new quadratic model by the least Frobenius */
2306 /* norm interpolant, making the replacement if the test is satisfied. */
2309 if (gqsq < gisq * 100.) {
2314 for (i__ = 1; i__ <= i__1; ++i__) {
2319 for (ih = 1; ih <= i__1; ++ih) {
2324 for (j = 1; j <= i__1; ++j) {
2327 for (k = 1; k <= i__2; ++k) {
2329 w[j] += vlag[k] * zmat[k + j * zmat_dim1];
2337 for (k = 1; k <= i__1; ++k) {
2340 for (j = 1; j <= i__2; ++j) {
2342 pq[k] += zmat[k + j * zmat_dim1] * w[j];
2353 /* If a trust region step has provided a sufficient decrease in F, then */
2354 /* branch for another trust region calculation. The case KSAVE>0 occurs */
2355 /* when the new function value was calculated by a model step. */
2357 if (f <= fsave + tenth * vquad) {
2364 /* Alternatively, find out if the interpolation points are close enough */
2365 /* to the best point so far. */
2369 distsq = delta * 4. * delta;
2371 for (k = 1; k <= i__2; ++k) {
2374 for (j = 1; j <= i__1; ++j) {
2376 /* Computing 2nd power */
2377 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
2387 /* If KNEW is positive, then set DSTEP, and branch back for the next */
2388 /* iteration, which will generate a "model step". */
2393 d__2 = tenth * sqrt(distsq), d__3 = half * delta;
2394 d__1 = min(d__2,d__3);
2395 dstep = max(d__1,rho);
2396 dsq = dstep * dstep;
2402 if (max(delta,dnorm) > rho) {
2406 /* The calculations with the current value of RHO are complete. Pick the */
2407 /* next values of RHO and DELTA. */
2412 ratio = rho / rhoend;
2415 } else if (ratio <= 250.) {
2416 rho = sqrt(ratio) * rhoend;
2420 delta = max(delta,rho);
2424 /* Return from the calculation, after another Newton-Raphson step, if */
2425 /* it is too short to have been tried before. */
2430 rc = NLOPT_XTOL_REACHED;
2434 for (i__ = 1; i__ <= i__2; ++i__) {
2436 x[i__] = xbase[i__] + xopt[i__];
2444 /*************************************************************************/
2447 nlopt_result newuoa(int n, int npt, double *x,
2448 const double *lb, const double *ub,
2449 double rhobeg, nlopt_stopping *stop, double *minf,
2450 newuoa_func calfun, void *calfun_data)
2452 /* Local variables */
2453 int id, np, iw, igq, ihq, ixb, ifv, ipq, ivl, ixn, ixo, ixp, ndim,
2458 /* This subroutine seeks the least value of a function of many variables, */
2459 /* by a trust region method that forms quadratic models by interpolation. */
2460 /* There can be some freedom in the interpolation conditions, which is */
2461 /* taken up by minimizing the Frobenius norm of the change to the second */
2462 /* derivative of the quadratic model, beginning with a zero matrix. The */
2463 /* arguments of the subroutine are as follows. */
2465 /* N must be set to the number of variables and must be at least two. */
2466 /* NPT is the number of interpolation conditions. Its value must be in the */
2467 /* interval [N+2,(N+1)(N+2)/2]. */
2468 /* Initial values of the variables must be set in X(1),X(2),...,X(N). They */
2469 /* will be changed to the values that give the least calculated F. */
2470 /* RHOBEG and RHOEND must be set to the initial and final values of a trust */
2471 /* region radius, so both must be positive with RHOEND<=RHOBEG. Typically */
2472 /* RHOBEG should be about one tenth of the greatest expected change to a */
2473 /* variable, and RHOEND should indicate the accuracy that is required in */
2474 /* the final values of the variables. */
2475 /* The value of IPRINT should be set to 0, 1, 2 or 3, which controls the */
2476 /* amount of printing. Specifically, there is no output if IPRINT=0 and */
2477 /* there is output only at the return if IPRINT=1. Otherwise, each new */
2478 /* value of RHO is printed, with the best vector of variables so far and */
2479 /* the corresponding value of the objective function. Further, each new */
2480 /* value of F with its variables are output if IPRINT=3. */
2481 /* MAXFUN must be set to an upper bound on the number of calls of CALFUN. */
2482 /* The array W will be used for working space. Its length must be at least */
2483 /* (NPT+13)*(NPT+N)+3*N*(N+3)/2. */
2485 /* SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set F to */
2486 /* the value of the objective function for the variables X(1),X(2),...,X(N). */
2488 /* Partition the working space array, so that different parts of it can be */
2489 /* treated separately by the subroutine that performs the main calculation. */
2491 /* Parameter adjustments */
2497 if (n < 2 || npt < n + 2 || npt > (n + 2) * np / 2) {
2498 return NLOPT_INVALID_ARGS;
2505 ifv = ixp + n * npt;
2508 ipq = ihq + n * np / 2;
2510 izmat = ibmat + ndim * n;
2511 id = izmat + npt * nptm;
2515 w = (double *) malloc(sizeof(double) * ((npt+13)*(npt+n) + 3*(n*(n+3))/2));
2516 if (!w) return NLOPT_OUT_OF_MEMORY;
2519 /* The above settings provide a partition of W for subroutine NEWUOB. */
2520 /* The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements of */
2521 /* W plus the space that is needed by the last array of NEWUOB. */
2523 ret = newuob_(&n, &npt, &x[1], &rhobeg,
2524 lb, ub, stop, minf, calfun, calfun_data,
2525 &w[ixb], &w[ixo], &w[ixn], &w[ixp], &w[ifv],
2526 &w[igq], &w[ihq], &w[ipq], &w[ibmat], &w[izmat],
2527 &ndim, &w[id], &w[ivl], &w[iw]);
2534 /*************************************************************************/
2535 /*************************************************************************/