6 #include "nlopt-util.h"
9 /* Common Block Declarations */
12 double fx, ldt, dmin__;
17 double *v; /* size n x n */
18 double *q0, *q1, *t_flin; /* size n */
19 double qa, qb, qc, qd0, qd1, qf1;
21 double fbest, *xbest; /* size n */
25 /* Table of constant values */
27 static int pow_ii(int x, int n) /* compute x^n, n >= 0 */
31 if (n & 1) { n--; p *= x; }
32 else { n >>= 1; x *= x; }
37 static void minfit_(int m, int n, double machep,
38 double *tol, double *ab, double *q, double *ework);
39 static nlopt_result min_(int n, int j, int nits, double *d2, double *x1, double *f1, int fk, praxis_func f, void *f_data, double *x, double *t_old, double machep, double *h__, struct global_s *global_1, struct q_s *q_1);
40 static double flin_(int n, int j, double *l, praxis_func f, void *f_data, double *x, int *nf, struct q_s *q_1, nlopt_result *ret);
41 static void sort_(int m, int n, double *d__, double *v);
42 static void quad_(int n, praxis_func f, void *f_data, double *x,
43 double *t_old, double machep, double *h__, struct global_s *global_1, struct q_s *q_1);
45 nlopt_result praxis_(double t0, double machep, double h0,
46 int n, double *x, praxis_func f, void *f_data,
47 nlopt_stopping *stop, double *minf)
49 /* System generated locals */
51 nlopt_result ret = NLOPT_SUCCESS;
55 struct global_s global_1;
59 double *d__, *y, *z__, *e_minfit, *prev_xbest; /* size n */
65 double m2, m4, t2_old, df, dn;
76 double ldfac, large, small, value;
81 /* LAST MODIFIED 3/1/73 */
82 /* Modified August 2007 by S. G. Johnson:
83 after conversion to C via f2c and some manual cleanup,
84 eliminating write statements, I additionally:
85 - modified the routine to use NLopt stop criteria
86 - allocate arrays dynamically to allow n > 20
87 - return the NLopt return code, where the min.
88 function value is now given by the parameter minf
90 /* PRAXIS RETURNS THE MINIMUM OF THE FUNCTION F(N,X) OF N VARIABLES */
91 /* USING THE PRINCIPAL AXIS METHOD. THE GRADIENT OF THE FUNCTION IS */
94 /* FOR A DESCRIPTION OF THE ALGORITHM, SEE CHAPTER SEVEN OF */
95 /* "ALGORITHMS FOR FINDING ZEROS AND EXTREMA OF FUNCTIONS WITHOUT */
96 /* CALCULATING DERIVATIVES" BY RICHARD P BRENT. */
98 /* THE PARAMETERS ARE: */
99 /* T0 IS A TOLERANCE. PRAXIS ATTEMPTS TO RETURN PRAXIS=F(X) */
100 /* SUCH THAT IF X0 IS THE TRUE LOCAL MINIMUM NEAR X, THEN */
101 /* NORM(X-X0) < T0 + SQUAREROOT(MACHEP)*NORM(X). */
102 /* MACHEP IS THE MACHINE PRECISION, THE SMALLEST NUMBER SUCH THAT */
103 /* 1 + MACHEP > 1. MACHEP SHOULD BE 16.**-13 (ABOUT */
104 /* 2.22D-16) FOR REAL*8 ARITHMETIC ON THE IBM 360. */
105 /* H0 IS THE MAXIMUM STEP SIZE. H0 SHOULD BE SET TO ABOUT THE */
106 /* MAXIMUM DISTANCE FROM THE INITIAL GUESS TO THE MINIMUM. */
107 /* (IF H0 IS SET TOO LARGE OR TOO SMALL, THE INITIAL RATE OF */
108 /* CONVERGENCE MAY BE SLOW.) */
109 /* N (AT LEAST TWO) IS THE NUMBER OF VARIABLES UPON WHICH */
110 /* THE FUNCTION DEPENDS. */
111 /* X IS AN ARRAY CONTAINING ON ENTRY A GUESS OF THE POINT OF */
112 /* MINIMUM, ON RETURN THE ESTIMATED POINT OF MINIMUM. */
113 /* F(N,X) IS THE FUNCTION TO BE MINIMIZED. F SHOULD BE A REAL*8 */
114 /* FUNCTION DECLARED EXTERNAL IN THE CALLING PROGRAM. */
115 /* THE APPROXIMATING QUADRATIC FORM IS */
116 /* Q(X') = F(N,X) + (1/2) * (X'-X)-TRANSPOSE * A * (X'-X) */
117 /* WHERE X IS THE BEST ESTIMATE OF THE MINIMUM AND A IS */
118 /* INVERSE(V-TRANSPOSE) * D * INVERSE(V) */
119 /* (V(*,*) IS THE MATRIX OF SEARCH DIRECTIONS; D(*) IS THE ARRAY */
120 /* OF SECOND DIFFERENCES). IF F HAS CONTINUOUS SECOND DERIVATIVES */
121 /* NEAR X0, A WILL TEND TO THE HESSIAN OF F AT X0 AS X APPROACHES X0. */
123 /* IT IS ASSUMED THAT ON FLOATING-POINT UNDERFLOW THE RESULT IS SET */
125 /* THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS AFTER */
126 /* THE INITIALIZATION OF MACHINE DEPENDENT NUMBERS. */
129 /* .....IF N>20 OR IF N<20 AND YOU NEED MORE SPACE, CHANGE '20' TO THE */
130 /* LARGEST VALUE OF N IN THE NEXT CARD, IN THE CARD 'IDIM=20', AND */
131 /* IN THE DIMENSION STATEMENTS IN SUBROUTINES MINFIT,MIN,FLIN,QUAD. */
132 /* ...changed by S. G. Johnson, 2007, to use malloc */
134 /* .....INITIALIZATION..... */
135 /* MACHINE DEPENDENT NUMBERS: */
137 /* Parameter adjustments */
141 small = machep * machep;
142 vsmall = small * small;
144 vlarge = 1. / vsmall;
148 /* new: dynamic allocation of temporary arrays */
149 work = (double *) malloc(sizeof(double) * (n*n + n*9));
150 if (!work) return NLOPT_OUT_OF_MEMORY;
152 q_1.q0 = q_1.v + n*n;
154 q_1.t_flin = q_1.q1 + n;
155 q_1.xbest = q_1.t_flin + n;
160 prev_xbest = e_minfit + n;
162 /* HEURISTIC NUMBERS: */
163 /* IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF */
164 /* POSSIBLE), THEN SET SCBD=10. OTHERWISE SET SCBD=1. */
165 /* IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED, SET ILLC=TRUE. */
166 /* OTHERWISE SET ILLC=FALSE. */
167 /* KTM IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT BEFORE THE */
168 /* ALGORITHM TERMINATES. KTM=4 IS VERY CAUTIOUS; USUALLY KTM=1 */
169 /* IS SATISFACTORY. */
172 illc = 0 /* false */;
182 prev_fbest = q_1.fbest = global_1.fx = f(n, &x[1], f_data);
183 memcpy(q_1.xbest, &x[1], n*sizeof(double));
184 memcpy(prev_xbest, &x[1], n*sizeof(double));
187 q_1.qf1 = global_1.fx;
192 for (i__ = 0; i__ < n; ++i__)
193 if (stop->xtol_abs[i__] > t_old)
194 t_old = stop->xtol_abs[i__];
198 global_1.dmin__ = small;
200 if (h__ < t_old * 100) {
204 /* .....THE FIRST SET OF SEARCH DIRECTIONS V IS THE IDENTITY MATRIX..... */
206 for (i__ = 1; i__ <= i__1; ++i__) {
208 for (j = 1; j <= i__2; ++j) {
210 q_1.v[i__ + j * n - (n+1)] = 0.;
213 q_1.v[i__ + i__ * n - (n+1)] = 1.;
218 for (i__ = 1; i__ <= i__1; ++i__) {
219 q_1.q0[i__ - 1] = x[i__];
221 q_1.q1[i__ - 1] = x[i__];
224 /* .....THE MAIN LOOP STARTS HERE..... */
230 /* .....MINIMIZE ALONG THE FIRST DIRECTION V(*,1). */
231 /* FX MUST BE PASSED TO MIN BY VALUE. */
233 ret = min_(n, 1, 2, d__, &s, &value, 0, f,f_data, &x[1],
234 &t_old, machep, &h__, &global_1, &q_1);
235 if (ret != NLOPT_SUCCESS) goto done;
240 for (i__ = 1; i__ <= i__1; ++i__) {
242 q_1.v[i__ - 1] = -q_1.v[i__ - 1];
245 if (sf > d__[0] * .9 && sf * .9 < d__[0]) {
249 for (i__ = 2; i__ <= i__1; ++i__) {
254 /* .....THE INNER LOOP STARTS HERE..... */
257 for (k = 2; k <= i__1; ++k) {
259 for (i__ = 1; i__ <= i__2; ++i__) {
271 /* .....A RANDOM STEP FOLLOWS (TO AVOID RESOLUTION VALLEYS). */
272 /* PRAXIS ASSUMES THAT RANDOM RETURNS A RANDOM NUMBER UNIFORMLY */
273 /* DISTRIBUTED IN (0,1). */
279 for (i__ = 1; i__ <= i__2; ++i__) {
280 s = (global_1.ldt * .1 + t2_old * pow_ii(10, kt))
281 * nlopt_urand(-.5,.5);
282 /* was: (random_(n) - .5); */
285 for (j = 1; j <= i__3; ++j) {
287 x[j] += s * q_1.v[j + i__ * n - (n+1)];
291 global_1.fx = (*f)(n, &x[1], f_data);
294 /* .....MINIMIZE ALONG THE "NON-CONJUGATE" DIRECTIONS V(*,K),...,V(*,N) */
298 for (k2 = k; k2 <= i__2; ++k2) {
302 ret = min_(n, k2, 2, &d__[k2 - 1], &s, &value, 0, f,f_data, &
303 x[1], &t_old, machep, &h__, &global_1, &q_1);
304 if (ret != NLOPT_SUCCESS) goto done;
308 s = sl - global_1.fx;
311 /* Computing 2nd power */
312 d__1 = s + z__[k2 - 1];
313 s = d__[k2 - 1] * (d__1 * d__1);
323 if (illc || df >= (d__1 = machep * 100 * global_1.fx, fabs(d__1))) {
327 /* .....IF THERE WAS NOT MUCH IMPROVEMENT ON THE FIRST TRY, SET */
328 /* ILLC=TRUE AND START THE INNER LOOP AGAIN..... */
334 /* .....MINIMIZE ALONG THE "CONJUGATE" DIRECTIONS V(*,1),...,V(*,K-1) */
338 for (k2 = 1; k2 <= i__2; ++k2) {
341 ret = min_(n, k2, 2, &d__[k2 - 1], &s, &value, 0, f,f_data, &
342 x[1], &t_old, machep, &h__, &global_1, &q_1);
343 if (ret != NLOPT_SUCCESS) goto done;
350 for (i__ = 1; i__ <= i__2; ++i__) {
363 /* .....DISCARD DIRECTION V(*,KL). */
364 /* IF NO RANDOM STEP WAS TAKEN, V(*,KL) IS THE "NON-CONJUGATE" */
365 /* DIRECTION ALONG WHICH THE GREATEST IMPROVEMENT WAS MADE..... */
372 for (ii = 1; ii <= i__2; ++ii) {
375 for (j = 1; j <= i__3; ++j) {
377 q_1.v[j + (i__ + 1) * n - (n+1)] = q_1.v[j + i__ * n - (n+1)];
380 d__[i__] = d__[i__ - 1];
385 for (i__ = 1; i__ <= i__2; ++i__) {
387 q_1.v[i__ + k * n - (n+1)] = y[i__ - 1] / lds;
390 /* .....MINIMIZE ALONG THE NEW "CONJUGATE" DIRECTION V(*,K), WHICH IS */
391 /* THE NORMALIZED VECTOR: (NEW X) - (0LD X)..... */
394 ret = min_(n, k, 4, &d__[k - 1], &lds, &value, 1, f,f_data, &x[1],
395 &t_old, machep, &h__, &global_1, &q_1);
396 if (ret != NLOPT_SUCCESS) goto done;
402 for (i__ = 1; i__ <= i__2; ++i__) {
404 q_1.v[i__ + k * n - (n+1)] = -q_1.v[i__ + k * n - (n+1)];
407 global_1.ldt = ldfac * global_1.ldt;
408 if (global_1.ldt < lds) {
413 for (i__ = 1; i__ <= i__2; ++i__) {
415 /* Computing 2nd power */
417 t2_old += d__1 * d__1;
419 t2_old = m2 * sqrt(t2_old) + t_old;
421 /* .....SEE WHETHER THE LENGTH OF THE STEP TAKEN SINCE STARTING THE */
422 /* INNER LOOP EXCEEDS HALF THE TOLERANCE..... */
424 if (global_1.ldt > t2_old * .5f
425 && !nlopt_stop_f(stop, q_1.fbest, prev_fbest)
426 && !nlopt_stop_x(stop, q_1.xbest, prev_xbest)) {
431 if (nlopt_stop_f(stop, q_1.fbest, prev_fbest))
432 ret = NLOPT_FTOL_REACHED;
433 else if (nlopt_stop_x(stop, q_1.xbest, prev_xbest))
434 ret = NLOPT_XTOL_REACHED;
437 prev_fbest = q_1.fbest;
438 memcpy(prev_xbest, q_1.xbest, n * sizeof(double));
441 /* .....THE INNER LOOP ENDS HERE. */
443 /* TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE IN A CURVED VALLEY. */
446 quad_(n, f,f_data, &x[1], &t_old, machep, &h__, &global_1, &q_1);
449 for (i__ = 1; i__ <= i__1; ++i__) {
450 d__[i__ - 1] = 1. / sqrt(d__[i__ - 1]);
451 if (dn < d__[i__ - 1]) {
457 for (j = 1; j <= i__1; ++j) {
460 for (i__ = 1; i__ <= i__2; ++i__) {
462 q_1.v[i__ + j * n - (n+1)] = s * q_1.v[i__ + j * n - (n+1)];
466 /* .....SCALE THE AXES TO TRY TO REDUCE THE CONDITION NUMBER..... */
473 for (i__ = 1; i__ <= i__2; ++i__) {
476 for (j = 1; j <= i__1; ++j) {
478 sl += q_1.v[i__ + j * n - (n+1)] * q_1.v[i__ + j * n - (n+1)];
480 z__[i__ - 1] = sqrt(sl);
481 if (z__[i__ - 1] < m4) {
484 if (s > z__[i__ - 1]) {
490 for (i__ = 1; i__ <= i__2; ++i__) {
491 sl = s / z__[i__ - 1];
492 z__[i__ - 1] = 1. / sl;
493 if (z__[i__ - 1] <= scbd) {
500 for (j = 1; j <= i__1; ++j) {
502 q_1.v[i__ + j * n - (n+1)] = sl * q_1.v[i__ + j * n - (n+1)];
507 /* .....CALCULATE A NEW SET OF ORTHOGONAL DIRECTIONS BEFORE REPEATING */
509 /* FIRST TRANSPOSE V FOR MINFIT: */
513 for (i__ = 2; i__ <= i__2; ++i__) {
516 for (j = 1; j <= i__1; ++j) {
517 s = q_1.v[i__ + j * n - (n+1)];
518 q_1.v[i__ + j * n - (n+1)] = q_1.v[j + i__ * n - (n+1)];
520 q_1.v[j + i__ * n - (n+1)] = s;
525 /* .....CALL MINFIT TO FIND THE SINGULAR VALUE DECOMPOSITION OF V. */
526 /* THIS GIVES THE PRINCIPAL VALUES AND PRINCIPAL DIRECTIONS OF THE */
527 /* APPROXIMATING QUADRATIC FORM WITHOUT SQUARING THE CONDITION */
530 minfit_(n, n, machep, &vsmall, q_1.v, d__, e_minfit);
532 /* .....UNSCALE THE AXES..... */
538 for (i__ = 1; i__ <= i__2; ++i__) {
541 for (j = 1; j <= i__1; ++j) {
543 q_1.v[i__ + j * n - (n+1)] = s * q_1.v[i__ + j * n - (n+1)];
548 for (i__ = 1; i__ <= i__2; ++i__) {
551 for (j = 1; j <= i__1; ++j) {
553 /* Computing 2nd power */
554 d__1 = q_1.v[j + i__ * n - (n+1)];
558 d__[i__ - 1] = s * d__[i__ - 1];
561 for (j = 1; j <= i__1; ++j) {
563 q_1.v[j + i__ * n - (n+1)] = s * q_1.v[j + i__ * n - (n+1)];
570 for (i__ = 1; i__ <= i__2; ++i__) {
571 dni = dn * d__[i__ - 1];
578 d__[i__ - 1] = 1 / (dni * dni);
581 d__[i__ - 1] = vlarge;
584 d__[i__ - 1] = vsmall;
589 /* .....SORT THE EIGENVALUES AND EIGENVECTORS..... */
591 sort_(n, n, d__, q_1.v);
592 global_1.dmin__ = d__[n - 1];
593 if (global_1.dmin__ < small) {
594 global_1.dmin__ = small;
596 illc = 0 /* false */;
597 if (m2 * d__[0] > global_1.dmin__) {
601 /* .....THE MAIN LOOP ENDS HERE..... */
605 /* .....RETURN..... */
608 if (ret != NLOPT_OUT_OF_MEMORY) {
610 memcpy(&x[1], q_1.xbest, n * sizeof(double));
616 static void minfit_(int m, int n, double machep,
617 double *tol, double *ab, double *q, double *ework)
619 /* System generated locals */
620 int ab_dim1, ab_offset, i__1, i__2, i__3;
623 /* Local variables */
624 double *e; /* size n */
625 double c__, f = 0.0, g, h__;
628 int l2, ii, kk, kt, ll2, lp1;
633 /* ...AN IMPROVED VERSION OF MINFIT (SEE GOLUB AND REINSCH, 1969) */
634 /* RESTRICTED TO M=N,P=0. */
635 /* THE SINGULAR VALUES OF THE ARRAY AB ARE RETURNED IN Q AND AB IS */
636 /* OVERWRITTEN WITH THE ORTHOGONAL MATRIX V SUCH THAT U.DIAG(Q) = AB.V, */
637 /* WHERE U IS ANOTHER ORTHOGONAL MATRIX. */
638 /* ...HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM... */
639 /* Parameter adjustments */
642 ab_offset = 1 + ab_dim1;
653 for (i__ = 1; i__ <= i__1; ++i__) {
658 for (j = i__; j <= i__2; ++j) {
660 /* Computing 2nd power */
661 d__1 = ab[j + i__ * ab_dim1];
668 f = ab[i__ + i__ * ab_dim1];
674 ab[i__ + i__ * ab_dim1] = f - g;
679 for (j = l; j <= i__2; ++j) {
682 for (k = i__; k <= i__3; ++k) {
684 f += ab[k + i__ * ab_dim1] * ab[k + j * ab_dim1];
688 for (k = i__; k <= i__3; ++k) {
690 ab[k + j * ab_dim1] += f * ab[k + i__ * ab_dim1];
700 for (j = l; j <= i__3; ++j) {
702 s += ab[i__ + j * ab_dim1] * ab[i__ + j * ab_dim1];
712 f = ab[i__ + (i__ + 1) * ab_dim1];
722 ab[i__ + (i__ + 1) * ab_dim1] = f - g;
724 for (j = l; j <= i__3; ++j) {
726 e[j - 1] = ab[i__ + j * ab_dim1] / h__;
729 for (j = l; j <= i__3; ++j) {
732 for (k = l; k <= i__2; ++k) {
734 s += ab[j + k * ab_dim1] * ab[i__ + k * ab_dim1];
737 for (k = l; k <= i__2; ++k) {
739 ab[j + k * ab_dim1] += s * e[k - 1];
743 y = (d__1 = q[i__], fabs(d__1)) + (d__2 = e[i__ - 1], fabs(d__2));
749 /* ...ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS... */
750 ab[n + n * ab_dim1] = 1.;
754 for (ii = 2; ii <= i__1; ++ii) {
759 h__ = ab[i__ + (i__ + 1) * ab_dim1] * g;
761 for (j = l; j <= i__2; ++j) {
763 ab[j + i__ * ab_dim1] = ab[i__ + j * ab_dim1] / h__;
766 for (j = l; j <= i__2; ++j) {
769 for (k = l; k <= i__3; ++k) {
771 s += ab[i__ + k * ab_dim1] * ab[k + j * ab_dim1];
774 for (k = l; k <= i__3; ++k) {
776 ab[k + j * ab_dim1] += s * ab[k + i__ * ab_dim1];
781 for (j = l; j <= i__3; ++j) {
782 ab[i__ + j * ab_dim1] = 0.;
784 ab[j + i__ * ab_dim1] = 0.;
786 ab[i__ + i__ * ab_dim1] = 1.;
791 /* ...DIAGONALIZATION OF THE BIDIAGONAL FORM... */
795 for (kk = 1; kk <= i__1; ++kk) {
804 /* fprintf(stderr, "QR failed\n"); */
807 for (ll2 = 1; ll2 <= i__3; ++ll2) {
810 if ((d__1 = e[l - 1], fabs(d__1)) <= eps) {
816 if ((d__1 = q[l - 1], fabs(d__1)) <= eps) {
822 /* ...CANCELLATION OF E(L) IF L>1... */
827 for (i__ = l; i__ <= i__3; ++i__) {
829 e[i__ - 1] = c__ * e[i__ - 1];
830 if (fabs(f) <= eps) {
834 /* ...Q(I) = H = DSQRT(G*G + F*F)... */
835 if (fabs(f) < fabs(g)) {
847 /* Computing 2nd power */
849 h__ = fabs(f) * sqrt(d__1 * d__1 + 1);
852 /* Computing 2nd power */
854 h__ = fabs(g) * sqrt(d__1 * d__1 + 1);
867 /* ...TEST FOR CONVERGENCE... */
873 /* ...SHIFT FROM BOTTOM 2*2 MINOR... */
878 f = ((y - z__) * (y + z__) + (g - h__) * (g + h__)) / (h__ * 2 * y);
879 g = sqrt(f * f + 1.);
884 f = ((x - z__) * (x + z__) + h__ * (y / temp - h__)) / x;
885 /* ...NEXT QR TRANSFORMATION... */
893 for (i__ = lp1; i__ <= i__3; ++i__) {
898 if (fabs(f) < fabs(h__)) {
910 /* Computing 2nd power */
912 z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
915 /* Computing 2nd power */
917 z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
929 g = -x * s + g * c__;
933 for (j = 1; j <= i__2; ++j) {
934 x = ab[j + (i__ - 1) * ab_dim1];
935 z__ = ab[j + i__ * ab_dim1];
936 ab[j + (i__ - 1) * ab_dim1] = x * c__ + z__ * s;
938 ab[j + i__ * ab_dim1] = -x * s + z__ * c__;
940 if (fabs(f) < fabs(h__)) {
952 /* Computing 2nd power */
954 z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
957 /* Computing 2nd power */
959 z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
972 x = -s * g + c__ * y;
979 /* ...CONVERGENCE: Q(K) IS MADE NON-NEGATIVE... */
986 for (j = 1; j <= i__3; ++j) {
988 ab[j + k * ab_dim1] = -ab[j + k * ab_dim1];
995 q[1] = ab[ab_dim1 + 1];
996 ab[ab_dim1 + 1] = 1.;
999 static nlopt_result min_(int n, int j, int nits, double *
1000 d2, double *x1, double *f1, int fk, praxis_func f, void *f_data, double *
1001 x, double *t_old, double machep, double *h__, struct global_s *global_1, struct q_s *q_1)
1003 /* System generated locals */
1007 /* Local variables */
1009 double s, d1, f0, f2, m2, m4, t2, x2, fm;
1011 double xm, sf1, sx1;
1013 nlopt_result ret = NLOPT_SUCCESS;
1015 /* ...THE SUBROUTINE MIN MINIMIZES F FROM X IN THE DIRECTION V(*,J) UNLESS */
1016 /* J IS LESS THAN 1, WHEN A QUADRATIC SEARCH IS MADE IN THE PLANE */
1017 /* DEFINED BY Q0,Q1,X. */
1018 /* D2 IS EITHER ZERO OR AN APPROXIMATION TO HALF F". */
1019 /* ON ENTRY, X1 IS AN ESTIMATE OF THE DISTANCE FROM X TO THE MINIMUM */
1020 /* ALONG V(*,J) (OR, IF J=0, A CURVE). ON RETURN, X1 IS THE DISTANCE */
1022 /* IF FK=.TRUE., THEN F1 IS FLIN(X1). OTHERWISE X1 AND F1 ARE IGNORED */
1023 /* ON ENTRY UNLESS FINAL FX IS GREATER THAN F1. */
1024 /* NITS CONTROLS THE NUMBER OF TIMES AN ATTEMPT WILL BE MADE TO HALVE */
1026 /* Parameter adjustments */
1030 /* Computing 2nd power */
1032 small = d__1 * d__1;
1042 /* ...FIND THE STEP SIZE... */
1045 for (i__ = 1; i__ <= i__1; ++i__) {
1047 /* Computing 2nd power */
1054 temp = global_1->dmin__;
1056 t2 = m4 * sqrt(fabs(global_1->fx) / temp + s * global_1->ldt) + m2 *
1058 s = m4 * s + *t_old;
1062 t2 = t2 > small ? t2 : small;
1064 d__1 = t2, d__2 = *h__ * .01;
1065 t2 = d__1 < d__2 ? d__1 : d__2;
1066 if (! (fk) || *f1 > fm) {
1072 if (fk && fabs(*x1) >= t2) {
1080 *f1 = flin_(n, j, x1, f,f_data, &x[1], &global_1->nf, q_1, &ret);
1081 if (ret != NLOPT_SUCCESS) return ret;
1092 /* ...EVALUATE FLIN AT ANOTHER POINT AND ESTIMATE THE SECOND DERIVATIVE... */
1097 f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1, &ret);
1098 if (ret != NLOPT_SUCCESS) return ret;
1105 *d2 = (x2 * (*f1 - f0) - *x1 * (f2 - f0)) / (*x1 * x2 * (*x1 - x2));
1106 /* ...ESTIMATE THE FIRST DERIVATIVE AT 0... */
1108 d1 = (*f1 - f0) / *x1 - *x1 * *d2;
1110 /* ...PREDICT THE MINIMUM... */
1120 x2 = d1 * -.5 / *d2;
1122 if (fabs(x2) <= *h__) {
1135 /* ...EVALUATE F AT THE PREDICTED MINIMUM... */
1137 f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1, &ret);
1138 if (ret != NLOPT_SUCCESS) return ret;
1139 if (k >= nits || f2 <= f0) {
1142 /* ...NO SUCCESS, SO TRY AGAIN... */
1144 if (f0 < *f1 && *x1 * x2 > 0.) {
1149 /* ...INCREMENT THE ONE-DIMENSIONAL SEARCH COUNTER... */
1159 /* ...GET A NEW ESTIMATE OF THE SECOND DERIVATIVE... */
1161 if ((d__1 = x2 * (x2 - *x1), fabs(d__1)) <= small) {
1164 *d2 = (x2 * (*f1 - f0) - *x1 * (fm - f0)) / (*x1 * x2 * (*x1 - x2));
1176 if (sf1 >= global_1->fx) {
1181 /* ...UPDATE X FOR LINEAR BUT NOT PARABOLIC SEARCH... */
1184 return NLOPT_SUCCESS;
1187 for (i__ = 1; i__ <= i__1; ++i__) {
1189 x[i__] += *x1 * q_1->v[i__ + j * n - (n+1)];
1191 return NLOPT_SUCCESS;
1194 static double flin_(int n, int j, double *l, praxis_func f, void *f_data, double *x,
1195 int *nf, struct q_s *q_1, nlopt_result *ret)
1197 /* System generated locals */
1201 /* Local variables */
1202 nlopt_stopping *stop = q_1->stop;
1204 double *t; /* size n */
1208 /* ...FLIN IS THE FUNCTION OF ONE REAL VARIABLE L THAT IS MINIMIZED */
1209 /* BY THE SUBROUTINE MIN... */
1210 /* Parameter adjustments */
1217 /* ...THE SEARCH IS LINEAR... */
1219 for (i__ = 1; i__ <= i__1; ++i__) {
1221 t[i__ - 1] = x[i__] + *l * q_1->v[i__ + j * n - (n+1)];
1224 /* ...THE SEARCH IS ALONG A PARABOLIC SPACE CURVE... */
1226 q_1->qa = *l * (*l - q_1->qd1) / (q_1->qd0 * (q_1->qd0 + q_1->qd1));
1227 q_1->qb = (*l + q_1->qd0) * (q_1->qd1 - *l) / (q_1->qd0 * q_1->qd1);
1228 q_1->qc = *l * (*l + q_1->qd0) / (q_1->qd1 * (q_1->qd0 + q_1->qd1));
1230 for (i__ = 1; i__ <= i__1; ++i__) {
1232 t[i__ - 1] = q_1->qa * q_1->q0[i__ - 1] + q_1->qb * x[i__] + q_1->qc *
1235 /* ...THE FUNCTION EVALUATION COUNTER NF IS INCREMENTED... */
1238 ret_val = f(n, t, f_data);
1240 if (ret_val < q_1->fbest) {
1241 q_1->fbest = ret_val;
1242 memcpy(q_1->xbest, t, n * sizeof(double));
1244 if (nlopt_stop_forced(stop)) *ret = NLOPT_FORCED_STOP;
1245 else if (nlopt_stop_evals(stop)) *ret = NLOPT_MAXEVAL_REACHED;
1246 else if (nlopt_stop_time(stop)) *ret = NLOPT_MAXTIME_REACHED;
1247 else if (ret_val <= stop->minf_max) *ret = NLOPT_MINF_MAX_REACHED;
1251 static void sort_(int m, int n, double *d__, double *v)
1253 /* System generated locals */
1254 int v_dim1, v_offset, i__1, i__2;
1256 /* Local variables */
1261 /* ...SORTS THE ELEMENTS OF D(N) INTO DESCENDING ORDER AND MOVES THE */
1262 /* CORRESPONDING COLUMNS OF V(N,N). */
1263 /* M IS THE ROW DIMENSION OF V AS DECLARED IN THE CALLING PROGRAM. */
1264 /* Parameter adjustments */
1266 v_offset = 1 + v_dim1;
1276 for (i__ = 1; i__ <= i__1; ++i__) {
1281 for (j = ip1; j <= i__2; ++j) {
1296 for (j = 1; j <= i__2; ++j) {
1297 s = v[j + i__ * v_dim1];
1298 v[j + i__ * v_dim1] = v[j + k * v_dim1];
1300 v[j + k * v_dim1] = s;
1307 static void quad_(int n, praxis_func f, void *f_data, double *x, double *t_old,
1308 double machep, double *h__, struct global_s *global_1, struct q_s *q_1)
1310 /* System generated locals */
1314 /* Local variables */
1319 /* ...QUAD LOOKS FOR THE MINIMUM OF F ALONG A CURVE DEFINED BY Q0,Q1,X... */
1320 /* Parameter adjustments */
1325 global_1->fx = q_1->qf1;
1329 for (i__ = 1; i__ <= i__1; ++i__) {
1331 l = q_1->q1[i__ - 1];
1333 q_1->q1[i__ - 1] = s;
1335 /* Computing 2nd power */
1337 q_1->qd1 += d__1 * d__1;
1339 q_1->qd1 = sqrt(q_1->qd1);
1342 if (q_1->qd0 <= 0. || q_1->qd1 <= 0. || global_1->nl < n * 3 * n) {
1346 min_(n, 0, 2, &s, &l, &value, 1, f,f_data, &x[1], t_old, machep,
1347 h__, global_1, q_1);
1348 q_1->qa = l * (l - q_1->qd1) / (q_1->qd0 * (q_1->qd0 + q_1->qd1));
1349 q_1->qb = (l + q_1->qd0) * (q_1->qd1 - l) / (q_1->qd0 * q_1->qd1);
1350 q_1->qc = l * (l + q_1->qd0) / (q_1->qd1 * (q_1->qd0 + q_1->qd1));
1353 global_1->fx = q_1->qf1;
1358 q_1->qd0 = q_1->qd1;
1360 for (i__ = 1; i__ <= i__1; ++i__) {
1361 s = q_1->q0[i__ - 1];
1362 q_1->q0[i__ - 1] = x[i__];
1364 x[i__] = q_1->qa * s + q_1->qb * x[i__] + q_1->qc * q_1->q1[i__ - 1];