praxis/praxis.c

   1 /* See README */
   2
   3 #include <stdlib.h>
   4 #include <math.h>
   5 #include <string.h>
   6 #include "nlopt-util.h"
   7 #include "praxis.h"
   8
   9 /* Common Block Declarations */
  10
  11 struct global_s {
  12     double fx, ldt, dmin__;
  13     int nf, nl;
  14 };
  15
  16 struct q_s {
  17      double *v; /* size n x n */
  18      double *q0, *q1, *t_flin; /* size n */
  19      double qa, qb, qc, qd0, qd1, qf1;
  20
  21      double fbest, *xbest; /* size n */
  22      nlopt_stopping *stop;
  23 };
  24
  25 /* Table of constant values */
  26
  27 static int pow_ii(int x, int n) /* compute x^n, n >= 0 */
  28 {
  29      int p = 1;
  30      while (n > 0) {
  31           if (n & 1) { n--; p *= x; }
  32           else { n >>= 1; x *= x; }
  33      }
  34      return p;
  35 }
  36
  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);
  44
  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)
  48 {
  49     /* System generated locals */
  50     int i__1, i__2, i__3;
  51     nlopt_result ret = NLOPT_SUCCESS;
  52     double d__1;
  53
  54     /* Global */
  55     struct global_s global_1;
  56     struct q_s q_1;
  57
  58     /* Local variables */
  59     double *d__, *y, *z__, *e_minfit, *prev_xbest; /* size n */
  60     double prev_fbest;
  61     double h__;
  62     int i__, j, k;
  63     double s, t_old, f1;
  64     int k2;
  65     double m2, m4, t2_old, df, dn;
  66     int kl, ii;
  67     double sf;
  68     int kt;
  69     double sl;
  70     int im1, km1;
  71     double dni, lds;
  72     int ktm;
  73     double scbd;
  74     int illc;
  75     int klmk;
  76     double ldfac, large, small, value;
  77     double vlarge;
  78     double vsmall;
  79     double *work;
  80
  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
  89 */
  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 */
  92 /*     NOT REQUIRED. */
  93
  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. */
  97
  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. */
 122
 123 /*     IT IS ASSUMED THAT ON FLOATING-POINT UNDERFLOW THE RESULT IS SET */
 124 /*     TO ZERO. */
 125 /*     THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS AFTER */
 126 /*     THE INITIALIZATION OF MACHINE DEPENDENT NUMBERS. */
 127
 128
 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 */
 133
 134 /* .....INITIALIZATION..... */
 135 /*     MACHINE DEPENDENT NUMBERS: */
 136
 137     /* Parameter adjustments */
 138     --x;
 139
 140     /* Function Body */
 141     small = machep * machep;
 142     vsmall = small * small;
 143     large = 1. / small;
 144     vlarge = 1. / vsmall;
 145     m2 = sqrt(machep);
 146     m4 = sqrt(m2);
 147
 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;
 151     q_1.v = work;
 152     q_1.q0 = q_1.v + n*n;
 153     q_1.q1 = q_1.q0 + n;
 154     q_1.t_flin = q_1.q1 + n;
 155     q_1.xbest = q_1.t_flin + n;
 156     d__ = q_1.xbest + n;
 157     y = d__ + n;
 158     z__ = y + n;
 159     e_minfit = y + n;
 160     prev_xbest = e_minfit + n;
 161
 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. */
 170
 171     scbd = 1.;
 172     illc = 0 /* false */;
 173     ktm = 1;
 174
 175     ldfac = .01;
 176     if (illc) {
 177         ldfac = .1;
 178     }
 179     kt = 0;
 180     global_1.nl = 0;
 181     global_1.nf = 1;
 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));
 185     stop->nevals++;
 186     q_1.stop = stop;
 187     q_1.qf1 = global_1.fx;
 188     if (t0 > 0)
 189          t_old = small + t0;
 190     else {
 191          t_old = 0;
 192          for (i__ = 0; i__ < n; ++i__)
 193               if (stop->xtol_abs[i__] > t_old)
 194                    t_old = stop->xtol_abs[i__];
 195          t_old += small;
 196     }
 197     t2_old = t_old;
 198     global_1.dmin__ = small;
 199     h__ = h0;
 200     if (h__ < t_old * 100) {
 201         h__ = t_old * 100;
 202     }
 203     global_1.ldt = h__;
 204 /* .....THE FIRST SET OF SEARCH DIRECTIONS V IS THE IDENTITY MATRIX..... */
 205     i__1 = n;
 206     for (i__ = 1; i__ <= i__1; ++i__) {
 207         i__2 = n;
 208         for (j = 1; j <= i__2; ++j) {
 209 /* L10: */
 210             q_1.v[i__ + j * n - (n+1)] = 0.;
 211         }
 212 /* L20: */
 213         q_1.v[i__ + i__ * n - (n+1)] = 1.;
 214     }
 215     d__[0] = 0.;
 216     q_1.qd0 = 0.;
 217     i__1 = n;
 218     for (i__ = 1; i__ <= i__1; ++i__) {
 219         q_1.q0[i__ - 1] = x[i__];
 220 /* L30: */
 221         q_1.q1[i__ - 1] = x[i__];
 222     }
 223
 224 /* .....THE MAIN LOOP STARTS HERE..... */
 225 L40:
 226     sf = d__[0];
 227     d__[0] = 0.;
 228     s = 0.;
 229
 230 /* .....MINIMIZE ALONG THE FIRST DIRECTION V(*,1). */
 231 /*     FX MUST BE PASSED TO MIN BY VALUE. */
 232     value = global_1.fx;
 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;
 236     if (s > 0.) {
 237         goto L50;
 238     }
 239     i__1 = n;
 240     for (i__ = 1; i__ <= i__1; ++i__) {
 241 /* L45: */
 242         q_1.v[i__ - 1] = -q_1.v[i__ - 1];
 243     }
 244 L50:
 245     if (sf > d__[0] * .9 && sf * .9 < d__[0]) {
 246         goto L70;
 247     }
 248     i__1 = n;
 249     for (i__ = 2; i__ <= i__1; ++i__) {
 250 /* L60: */
 251         d__[i__ - 1] = 0.;
 252     }
 253
 254 /* .....THE INNER LOOP STARTS HERE..... */
 255 L70:
 256     i__1 = n;
 257     for (k = 2; k <= i__1; ++k) {
 258         i__2 = n;
 259         for (i__ = 1; i__ <= i__2; ++i__) {
 260 /* L75: */
 261             y[i__ - 1] = x[i__];
 262         }
 263         sf = global_1.fx;
 264         if (kt > 0) {
 265             illc = 1 /* true */;
 266         }
 267 L80:
 268         kl = k;
 269         df = 0.;
 270
 271 /* .....A RANDOM STEP FOLLOWS (TO AVOID RESOLUTION VALLEYS). */
 272 /*     PRAXIS ASSUMES THAT RANDOM RETURNS A RANDOM NUMBER UNIFORMLY */
 273 /*     DISTRIBUTED IN (0,1). */
 274
 275         if (! illc) {
 276             goto L95;
 277         }
 278         i__2 = n;
 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); */
 283             z__[i__ - 1] = s;
 284             i__3 = n;
 285             for (j = 1; j <= i__3; ++j) {
 286 /* L85: */
 287                 x[j] += s * q_1.v[j + i__ * n - (n+1)];
 288             }
 289 /* L90: */
 290         }
 291         global_1.fx = (*f)(n, &x[1], f_data);
 292         ++global_1.nf;
 293
 294 /* .....MINIMIZE ALONG THE "NON-CONJUGATE" DIRECTIONS V(*,K),...,V(*,N) */
 295
 296 L95:
 297         i__2 = n;
 298         for (k2 = k; k2 <= i__2; ++k2) {
 299             sl = global_1.fx;
 300             s = 0.;
 301             value = global_1.fx;
 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;
 305             if (illc) {
 306                 goto L97;
 307             }
 308             s = sl - global_1.fx;
 309             goto L99;
 310 L97:
 311 /* Computing 2nd power */
 312             d__1 = s + z__[k2 - 1];
 313             s = d__[k2 - 1] * (d__1 * d__1);
 314 L99:
 315             if (df > s) {
 316                 goto L105;
 317             }
 318             df = s;
 319             kl = k2;
 320 L105:
 321             ;
 322         }
 323         if (illc || df >= (d__1 = machep * 100 * global_1.fx, fabs(d__1))) {
 324             goto L110;
 325         }
 326
 327 /* .....IF THERE WAS NOT MUCH IMPROVEMENT ON THE FIRST TRY, SET */
 328 /*     ILLC=TRUE AND START THE INNER LOOP AGAIN..... */
 329
 330         illc = 1 /* true */;
 331         goto L80;
 332 L110:
 333
 334 /* .....MINIMIZE ALONG THE "CONJUGATE" DIRECTIONS V(*,1),...,V(*,K-1) */
 335
 336         km1 = k - 1;
 337         i__2 = km1;
 338         for (k2 = 1; k2 <= i__2; ++k2) {
 339             s = 0.;
 340             value = global_1.fx;
 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;
 344 /* L120: */
 345         }
 346         f1 = global_1.fx;
 347         global_1.fx = sf;
 348         lds = 0.;
 349         i__2 = n;
 350         for (i__ = 1; i__ <= i__2; ++i__) {
 351             sl = x[i__];
 352             x[i__] = y[i__ - 1];
 353             sl -= y[i__ - 1];
 354             y[i__ - 1] = sl;
 355 /* L130: */
 356             lds += sl * sl;
 357         }
 358         lds = sqrt(lds);
 359         if (lds <= small) {
 360             goto L160;
 361         }
 362
 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..... */
 366
 367         klmk = kl - k;
 368         if (klmk < 1) {
 369             goto L141;
 370         }
 371         i__2 = klmk;
 372         for (ii = 1; ii <= i__2; ++ii) {
 373             i__ = kl - ii;
 374             i__3 = n;
 375             for (j = 1; j <= i__3; ++j) {
 376 /* L135: */
 377                 q_1.v[j + (i__ + 1) * n - (n+1)] = q_1.v[j + i__ * n - (n+1)];
 378             }
 379 /* L140: */
 380             d__[i__] = d__[i__ - 1];
 381         }
 382 L141:
 383         d__[k - 1] = 0.;
 384         i__2 = n;
 385         for (i__ = 1; i__ <= i__2; ++i__) {
 386 /* L145: */
 387             q_1.v[i__ + k * n - (n+1)] = y[i__ - 1] / lds;
 388         }
 389
 390 /* .....MINIMIZE ALONG THE NEW "CONJUGATE" DIRECTION V(*,K), WHICH IS */
 391 /*     THE NORMALIZED VECTOR:  (NEW X) - (0LD X)..... */
 392
 393         value = f1;
 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;
 397         if (lds > 0.) {
 398             goto L160;
 399         }
 400         lds = -lds;
 401         i__2 = n;
 402         for (i__ = 1; i__ <= i__2; ++i__) {
 403 /* L150: */
 404             q_1.v[i__ + k * n - (n+1)] = -q_1.v[i__ + k * n - (n+1)];
 405         }
 406 L160:
 407         global_1.ldt = ldfac * global_1.ldt;
 408         if (global_1.ldt < lds) {
 409             global_1.ldt = lds;
 410         }
 411         t2_old = 0.;
 412         i__2 = n;
 413         for (i__ = 1; i__ <= i__2; ++i__) {
 414 /* L165: */
 415 /* Computing 2nd power */
 416             d__1 = x[i__];
 417             t2_old += d__1 * d__1;
 418         }
 419         t2_old = m2 * sqrt(t2_old) + t_old;
 420
 421 /* .....SEE WHETHER THE LENGTH OF THE STEP TAKEN SINCE STARTING THE */
 422 /*     INNER LOOP EXCEEDS HALF THE TOLERANCE..... */
 423
 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)) {
 427             kt = -1;
 428         }
 429         ++kt;
 430         if (kt > ktm) {
 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;
 435              goto done;
 436         }
 437         prev_fbest = q_1.fbest;
 438         memcpy(prev_xbest, q_1.xbest, n * sizeof(double));
 439 /* L170: */
 440     }
 441 /* .....THE INNER LOOP ENDS HERE. */
 442
 443 /*     TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE IN A CURVED VALLEY. */
 444
 445 /* L171: */
 446     quad_(n, f,f_data, &x[1], &t_old, machep, &h__, &global_1, &q_1);
 447     dn = 0.;
 448     i__1 = n;
 449     for (i__ = 1; i__ <= i__1; ++i__) {
 450         d__[i__ - 1] = 1. / sqrt(d__[i__ - 1]);
 451         if (dn < d__[i__ - 1]) {
 452             dn = d__[i__ - 1];
 453         }
 454 /* L175: */
 455     }
 456     i__1 = n;
 457     for (j = 1; j <= i__1; ++j) {
 458         s = d__[j - 1] / dn;
 459         i__2 = n;
 460         for (i__ = 1; i__ <= i__2; ++i__) {
 461 /* L180: */
 462             q_1.v[i__ + j * n - (n+1)] = s * q_1.v[i__ + j * n - (n+1)];
 463         }
 464     }
 465
 466 /* .....SCALE THE AXES TO TRY TO REDUCE THE CONDITION NUMBER..... */
 467
 468     if (scbd <= 1.) {
 469         goto L200;
 470     }
 471     s = vlarge;
 472     i__2 = n;
 473     for (i__ = 1; i__ <= i__2; ++i__) {
 474         sl = 0.;
 475         i__1 = n;
 476         for (j = 1; j <= i__1; ++j) {
 477 /* L182: */
 478             sl += q_1.v[i__ + j * n - (n+1)] * q_1.v[i__ + j * n - (n+1)];
 479         }
 480         z__[i__ - 1] = sqrt(sl);
 481         if (z__[i__ - 1] < m4) {
 482             z__[i__ - 1] = m4;
 483         }
 484         if (s > z__[i__ - 1]) {
 485             s = z__[i__ - 1];
 486         }
 487 /* L185: */
 488     }
 489     i__2 = n;
 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) {
 494             goto L189;
 495         }
 496         sl = 1. / scbd;
 497         z__[i__ - 1] = scbd;
 498 L189:
 499         i__1 = n;
 500         for (j = 1; j <= i__1; ++j) {
 501 /* L190: */
 502             q_1.v[i__ + j * n - (n+1)] = sl * q_1.v[i__ + j * n - (n+1)];
 503         }
 504 /* L195: */
 505     }
 506
 507 /* .....CALCULATE A NEW SET OF ORTHOGONAL DIRECTIONS BEFORE REPEATING */
 508 /*     THE MAIN LOOP. */
 509 /*     FIRST TRANSPOSE V FOR MINFIT: */
 510
 511 L200:
 512     i__2 = n;
 513     for (i__ = 2; i__ <= i__2; ++i__) {
 514         im1 = i__ - 1;
 515         i__1 = im1;
 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)];
 519 /* L210: */
 520             q_1.v[j + i__ * n - (n+1)] = s;
 521         }
 522 /* L220: */
 523     }
 524
 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 */
 528 /*     NUMBER..... */
 529
 530     minfit_(n, n, machep, &vsmall, q_1.v, d__, e_minfit);
 531
 532 /* .....UNSCALE THE AXES..... */
 533
 534     if (scbd <= 1.) {
 535         goto L250;
 536     }
 537     i__2 = n;
 538     for (i__ = 1; i__ <= i__2; ++i__) {
 539         s = z__[i__ - 1];
 540         i__1 = n;
 541         for (j = 1; j <= i__1; ++j) {
 542 /* L225: */
 543             q_1.v[i__ + j * n - (n+1)] = s * q_1.v[i__ + j * n - (n+1)];
 544         }
 545 /* L230: */
 546     }
 547     i__2 = n;
 548     for (i__ = 1; i__ <= i__2; ++i__) {
 549         s = 0.;
 550         i__1 = n;
 551         for (j = 1; j <= i__1; ++j) {
 552 /* L235: */
 553 /* Computing 2nd power */
 554             d__1 = q_1.v[j + i__ * n - (n+1)];
 555             s += d__1 * d__1;
 556         }
 557         s = sqrt(s);
 558         d__[i__ - 1] = s * d__[i__ - 1];
 559         s = 1 / s;
 560         i__1 = n;
 561         for (j = 1; j <= i__1; ++j) {
 562 /* L240: */
 563             q_1.v[j + i__ * n - (n+1)] = s * q_1.v[j + i__ * n - (n+1)];
 564         }
 565 /* L245: */
 566     }
 567
 568 L250:
 569     i__2 = n;
 570     for (i__ = 1; i__ <= i__2; ++i__) {
 571         dni = dn * d__[i__ - 1];
 572         if (dni > large) {
 573             goto L265;
 574         }
 575         if (dni < small) {
 576             goto L260;
 577         }
 578         d__[i__ - 1] = 1 / (dni * dni);
 579         goto L270;
 580 L260:
 581         d__[i__ - 1] = vlarge;
 582         goto L270;
 583 L265:
 584         d__[i__ - 1] = vsmall;
 585 L270:
 586         ;
 587     }
 588
 589 /* .....SORT THE EIGENVALUES AND EIGENVECTORS..... */
 590
 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;
 595     }
 596     illc = 0 /* false */;
 597     if (m2 * d__[0] > global_1.dmin__) {
 598         illc = 1 /* true */;
 599     }
 600
 601 /* .....THE MAIN LOOP ENDS HERE..... */
 602
 603     goto L40;
 604
 605 /* .....RETURN..... */
 606
 607 done:
 608     if (ret != NLOPT_OUT_OF_MEMORY) {
 609          *minf = q_1.fbest;
 610          memcpy(&x[1], q_1.xbest, n * sizeof(double));
 611     }
 612     free(work);
 613     return ret;
 614 } /* praxis_ */
 615
 616 static void minfit_(int m, int n, double machep,
 617         double *tol, double *ab, double *q, double *ework)
 618 {
 619     /* System generated locals */
 620     int ab_dim1, ab_offset, i__1, i__2, i__3;
 621     double d__1, d__2;
 622
 623     /* Local variables */
 624     double *e; /* size n */
 625     double c__, f = 0.0, g, h__;
 626     int i__, j, k, l;
 627     double s, x, y, z__;
 628     int l2, ii, kk, kt, ll2, lp1;
 629     double eps, temp;
 630
 631     e = ework;
 632
 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 */
 640     --q;
 641     ab_dim1 = m;
 642     ab_offset = 1 + ab_dim1;
 643     ab -= ab_offset;
 644
 645     /* Function Body */
 646     if (n == 1) {
 647         goto L200;
 648     }
 649     eps = machep;
 650     g = 0.;
 651     x = 0.;
 652     i__1 = n;
 653     for (i__ = 1; i__ <= i__1; ++i__) {
 654         e[i__ - 1] = g;
 655         s = 0.;
 656         l = i__ + 1;
 657         i__2 = n;
 658         for (j = i__; j <= i__2; ++j) {
 659 /* L1: */
 660 /* Computing 2nd power */
 661             d__1 = ab[j + i__ * ab_dim1];
 662             s += d__1 * d__1;
 663         }
 664         g = 0.;
 665         if (s < *tol) {
 666             goto L4;
 667         }
 668         f = ab[i__ + i__ * ab_dim1];
 669         g = sqrt(s);
 670         if (f >= 0.) {
 671             g = -g;
 672         }
 673         h__ = f * g - s;
 674         ab[i__ + i__ * ab_dim1] = f - g;
 675         if (l > n) {
 676             goto L4;
 677         }
 678         i__2 = n;
 679         for (j = l; j <= i__2; ++j) {
 680             f = 0.;
 681             i__3 = n;
 682             for (k = i__; k <= i__3; ++k) {
 683 /* L2: */
 684                 f += ab[k + i__ * ab_dim1] * ab[k + j * ab_dim1];
 685             }
 686             f /= h__;
 687             i__3 = n;
 688             for (k = i__; k <= i__3; ++k) {
 689 /* L3: */
 690                 ab[k + j * ab_dim1] += f * ab[k + i__ * ab_dim1];
 691             }
 692         }
 693 L4:
 694         q[i__] = g;
 695         s = 0.;
 696         if (i__ == n) {
 697             goto L6;
 698         }
 699         i__3 = n;
 700         for (j = l; j <= i__3; ++j) {
 701 /* L5: */
 702             s += ab[i__ + j * ab_dim1] * ab[i__ + j * ab_dim1];
 703         }
 704 L6:
 705         g = 0.;
 706         if (s < *tol) {
 707             goto L10;
 708         }
 709         if (i__ == n) {
 710             goto L16;
 711         }
 712         f = ab[i__ + (i__ + 1) * ab_dim1];
 713 L16:
 714         g = sqrt(s);
 715         if (f >= 0.) {
 716             g = -g;
 717         }
 718         h__ = f * g - s;
 719         if (i__ == n) {
 720             goto L10;
 721         }
 722         ab[i__ + (i__ + 1) * ab_dim1] = f - g;
 723         i__3 = n;
 724         for (j = l; j <= i__3; ++j) {
 725 /* L7: */
 726             e[j - 1] = ab[i__ + j * ab_dim1] / h__;
 727         }
 728         i__3 = n;
 729         for (j = l; j <= i__3; ++j) {
 730             s = 0.;
 731             i__2 = n;
 732             for (k = l; k <= i__2; ++k) {
 733 /* L8: */
 734                 s += ab[j + k * ab_dim1] * ab[i__ + k * ab_dim1];
 735             }
 736             i__2 = n;
 737             for (k = l; k <= i__2; ++k) {
 738 /* L9: */
 739                 ab[j + k * ab_dim1] += s * e[k - 1];
 740             }
 741         }
 742 L10:
 743         y = (d__1 = q[i__], fabs(d__1)) + (d__2 = e[i__ - 1], fabs(d__2));
 744 /* L11: */
 745         if (y > x) {
 746             x = y;
 747         }
 748     }
 749 /* ...ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS... */
 750     ab[n + n * ab_dim1] = 1.;
 751     g = e[n - 1];
 752     l = n;
 753     i__1 = n;
 754     for (ii = 2; ii <= i__1; ++ii) {
 755         i__ = n - ii + 1;
 756         if (g == 0.) {
 757             goto L23;
 758         }
 759         h__ = ab[i__ + (i__ + 1) * ab_dim1] * g;
 760         i__2 = n;
 761         for (j = l; j <= i__2; ++j) {
 762 /* L20: */
 763             ab[j + i__ * ab_dim1] = ab[i__ + j * ab_dim1] / h__;
 764         }
 765         i__2 = n;
 766         for (j = l; j <= i__2; ++j) {
 767             s = 0.;
 768             i__3 = n;
 769             for (k = l; k <= i__3; ++k) {
 770 /* L21: */
 771                 s += ab[i__ + k * ab_dim1] * ab[k + j * ab_dim1];
 772             }
 773             i__3 = n;
 774             for (k = l; k <= i__3; ++k) {
 775 /* L22: */
 776                 ab[k + j * ab_dim1] += s * ab[k + i__ * ab_dim1];
 777             }
 778         }
 779 L23:
 780         i__3 = n;
 781         for (j = l; j <= i__3; ++j) {
 782             ab[i__ + j * ab_dim1] = 0.;
 783 /* L24: */
 784             ab[j + i__ * ab_dim1] = 0.;
 785         }
 786         ab[i__ + i__ * ab_dim1] = 1.;
 787         g = e[i__ - 1];
 788 /* L25: */
 789         l = i__;
 790     }
 791 /* ...DIAGONALIZATION OF THE BIDIAGONAL FORM... */
 792 /* L100: */
 793     eps *= x;
 794     i__1 = n;
 795     for (kk = 1; kk <= i__1; ++kk) {
 796         k = n - kk + 1;
 797         kt = 0;
 798 L101:
 799         ++kt;
 800         if (kt <= 30) {
 801             goto L102;
 802         }
 803         e[k - 1] = 0.;
 804         /* fprintf(stderr, "QR failed\n"); */
 805 L102:
 806         i__3 = k;
 807         for (ll2 = 1; ll2 <= i__3; ++ll2) {
 808             l2 = k - ll2 + 1;
 809             l = l2;
 810             if ((d__1 = e[l - 1], fabs(d__1)) <= eps) {
 811                 goto L120;
 812             }
 813             if (l == 1) {
 814                 goto L103;
 815             }
 816             if ((d__1 = q[l - 1], fabs(d__1)) <= eps) {
 817                 goto L110;
 818             }
 819 L103:
 820             ;
 821         }
 822 /* ...CANCELLATION OF E(L) IF L>1... */
 823 L110:
 824         c__ = 0.;
 825         s = 1.;
 826         i__3 = k;
 827         for (i__ = l; i__ <= i__3; ++i__) {
 828             f = s * e[i__ - 1];
 829             e[i__ - 1] = c__ * e[i__ - 1];
 830             if (fabs(f) <= eps) {
 831                 goto L120;
 832             }
 833             g = q[i__];
 834 /* ...Q(I) = H = DSQRT(G*G + F*F)... */
 835             if (fabs(f) < fabs(g)) {
 836                 goto L113;
 837             }
 838             if (f != 0.) {
 839                 goto L112;
 840             } else {
 841                 goto L111;
 842             }
 843 L111:
 844             h__ = 0.;
 845             goto L114;
 846 L112:
 847 /* Computing 2nd power */
 848             d__1 = g / f;
 849             h__ = fabs(f) * sqrt(d__1 * d__1 + 1);
 850             goto L114;
 851 L113:
 852 /* Computing 2nd power */
 853             d__1 = f / g;
 854             h__ = fabs(g) * sqrt(d__1 * d__1 + 1);
 855 L114:
 856             q[i__] = h__;
 857             if (h__ != 0.) {
 858                 goto L115;
 859             }
 860             g = 1.;
 861             h__ = 1.;
 862 L115:
 863             c__ = g / h__;
 864 /* L116: */
 865             s = -f / h__;
 866         }
 867 /* ...TEST FOR CONVERGENCE... */
 868 L120:
 869         z__ = q[k];
 870         if (l == k) {
 871             goto L140;
 872         }
 873 /* ...SHIFT FROM BOTTOM 2*2 MINOR... */
 874         x = q[l];
 875         y = q[k - 1];
 876         g = e[k - 2];
 877         h__ = e[k - 1];
 878         f = ((y - z__) * (y + z__) + (g - h__) * (g + h__)) / (h__ * 2 * y);
 879         g = sqrt(f * f + 1.);
 880         temp = f - g;
 881         if (f >= 0.) {
 882             temp = f + g;
 883         }
 884         f = ((x - z__) * (x + z__) + h__ * (y / temp - h__)) / x;
 885 /* ...NEXT QR TRANSFORMATION... */
 886         c__ = 1.;
 887         s = 1.;
 888         lp1 = l + 1;
 889         if (lp1 > k) {
 890             goto L133;
 891         }
 892         i__3 = k;
 893         for (i__ = lp1; i__ <= i__3; ++i__) {
 894             g = e[i__ - 1];
 895             y = q[i__];
 896             h__ = s * g;
 897             g *= c__;
 898             if (fabs(f) < fabs(h__)) {
 899                 goto L123;
 900             }
 901             if (f != 0.) {
 902                 goto L122;
 903             } else {
 904                 goto L121;
 905             }
 906 L121:
 907             z__ = 0.;
 908             goto L124;
 909 L122:
 910 /* Computing 2nd power */
 911             d__1 = h__ / f;
 912             z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
 913             goto L124;
 914 L123:
 915 /* Computing 2nd power */
 916             d__1 = f / h__;
 917             z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
 918 L124:
 919             e[i__ - 2] = z__;
 920             if (z__ != 0.) {
 921                 goto L125;
 922             }
 923             f = 1.;
 924             z__ = 1.;
 925 L125:
 926             c__ = f / z__;
 927             s = h__ / z__;
 928             f = x * c__ + g * s;
 929             g = -x * s + g * c__;
 930             h__ = y * s;
 931             y *= c__;
 932             i__2 = n;
 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;
 937 /* L126: */
 938                 ab[j + i__ * ab_dim1] = -x * s + z__ * c__;
 939             }
 940             if (fabs(f) < fabs(h__)) {
 941                 goto L129;
 942             }
 943             if (f != 0.) {
 944                 goto L128;
 945             } else {
 946                 goto L127;
 947             }
 948 L127:
 949             z__ = 0.;
 950             goto L130;
 951 L128:
 952 /* Computing 2nd power */
 953             d__1 = h__ / f;
 954             z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
 955             goto L130;
 956 L129:
 957 /* Computing 2nd power */
 958             d__1 = f / h__;
 959             z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
 960 L130:
 961             q[i__ - 1] = z__;
 962             if (z__ != 0.) {
 963                 goto L131;
 964             }
 965             f = 1.;
 966             z__ = 1.;
 967 L131:
 968             c__ = f / z__;
 969             s = h__ / z__;
 970             f = c__ * g + s * y;
 971 /* L132: */
 972             x = -s * g + c__ * y;
 973         }
 974 L133:
 975         e[l - 1] = 0.;
 976         e[k - 1] = f;
 977         q[k] = x;
 978         goto L101;
 979 /* ...CONVERGENCE:  Q(K) IS MADE NON-NEGATIVE... */
 980 L140:
 981         if (z__ >= 0.) {
 982             goto L150;
 983         }
 984         q[k] = -z__;
 985         i__3 = n;
 986         for (j = 1; j <= i__3; ++j) {
 987 /* L141: */
 988             ab[j + k * ab_dim1] = -ab[j + k * ab_dim1];
 989         }
 990 L150:
 991         ;
 992     }
 993     return;
 994 L200:
 995     q[1] = ab[ab_dim1 + 1];
 996     ab[ab_dim1 + 1] = 1.;
 997 } /* minfit_ */
 998
 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)
1002 {
1003     /* System generated locals */
1004     int i__1;
1005     double d__1, d__2;
1006
1007     /* Local variables */
1008     int i__, k;
1009     double s, d1, f0, f2, m2, m4, t2, x2, fm;
1010     int dz;
1011     double xm, sf1, sx1;
1012     double temp, small;
1013     nlopt_result ret = NLOPT_SUCCESS;
1014
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 */
1021 /*   FOUND. */
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 */
1025 /*   THE INTERVAL. */
1026     /* Parameter adjustments */
1027     --x;
1028
1029     /* Function Body */
1030 /* Computing 2nd power */
1031     d__1 = machep;
1032     small = d__1 * d__1;
1033     m2 = sqrt(machep);
1034     m4 = sqrt(m2);
1035     sf1 = *f1;
1036     sx1 = *x1;
1037     k = 0;
1038     xm = 0.;
1039     fm = global_1->fx;
1040     f0 = global_1->fx;
1041     dz = *d2 < machep;
1042 /* ...FIND THE STEP SIZE... */
1043     s = 0.;
1044     i__1 = n;
1045     for (i__ = 1; i__ <= i__1; ++i__) {
1046 /* L1: */
1047 /* Computing 2nd power */
1048         d__1 = x[i__];
1049         s += d__1 * d__1;
1050     }
1051     s = sqrt(s);
1052     temp = *d2;
1053     if (dz) {
1054         temp = global_1->dmin__;
1055     }
1056     t2 = m4 * sqrt(fabs(global_1->fx) / temp + s * global_1->ldt) + m2 *
1057             global_1->ldt;
1058     s = m4 * s + *t_old;
1059     if (dz && t2 > s) {
1060         t2 = s;
1061     }
1062     t2 = t2 > small ? t2 : small;
1063 /* Computing MIN */
1064     d__1 = t2, d__2 = *h__ * .01;
1065     t2 = d__1 < d__2 ? d__1 : d__2;
1066     if (! (fk) || *f1 > fm) {
1067         goto L2;
1068     }
1069     xm = *x1;
1070     fm = *f1;
1071 L2:
1072     if (fk && fabs(*x1) >= t2) {
1073         goto L3;
1074     }
1075     temp = 1.;
1076     if (*x1 < 0.) {
1077         temp = -1.;
1078     }
1079     *x1 = temp * 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;
1082 L3:
1083     if (*f1 > fm) {
1084         goto L4;
1085     }
1086     xm = *x1;
1087     fm = *f1;
1088 L4:
1089     if (! dz) {
1090         goto L6;
1091     }
1092 /* ...EVALUATE FLIN AT ANOTHER POINT AND ESTIMATE THE SECOND DERIVATIVE... */
1093     x2 = -(*x1);
1094     if (f0 >= *f1) {
1095         x2 = *x1 * 2.;
1096     }
1097     f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1, &ret);
1098     if (ret != NLOPT_SUCCESS) return ret;
1099     if (f2 > fm) {
1100         goto L5;
1101     }
1102     xm = x2;
1103     fm = f2;
1104 L5:
1105     *d2 = (x2 * (*f1 - f0) - *x1 * (f2 - f0)) / (*x1 * x2 * (*x1 - x2));
1106 /* ...ESTIMATE THE FIRST DERIVATIVE AT 0... */
1107 L6:
1108     d1 = (*f1 - f0) / *x1 - *x1 * *d2;
1109     dz = 1 /* true */;
1110 /* ...PREDICT THE MINIMUM... */
1111     if (*d2 > small) {
1112         goto L7;
1113     }
1114     x2 = *h__;
1115     if (d1 >= 0.) {
1116         x2 = -x2;
1117     }
1118     goto L8;
1119 L7:
1120     x2 = d1 * -.5 / *d2;
1121 L8:
1122     if (fabs(x2) <= *h__) {
1123         goto L11;
1124     }
1125     if (x2 <= 0.) {
1126         goto L9;
1127     } else {
1128         goto L10;
1129     }
1130 L9:
1131     x2 = -(*h__);
1132     goto L11;
1133 L10:
1134     x2 = *h__;
1135 /* ...EVALUATE F AT THE PREDICTED MINIMUM... */
1136 L11:
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) {
1140         goto L12;
1141     }
1142 /* ...NO SUCCESS, SO TRY AGAIN... */
1143     ++k;
1144     if (f0 < *f1 && *x1 * x2 > 0.) {
1145         goto L4;
1146     }
1147     x2 *= .5;
1148     goto L11;
1149 /* ...INCREMENT THE ONE-DIMENSIONAL SEARCH COUNTER... */
1150 L12:
1151     ++global_1->nl;
1152     if (f2 <= fm) {
1153         goto L13;
1154     }
1155     x2 = xm;
1156     goto L14;
1157 L13:
1158     fm = f2;
1159 /* ...GET A NEW ESTIMATE OF THE SECOND DERIVATIVE... */
1160 L14:
1161     if ((d__1 = x2 * (x2 - *x1), fabs(d__1)) <= small) {
1162         goto L15;
1163     }
1164     *d2 = (x2 * (*f1 - f0) - *x1 * (fm - f0)) / (*x1 * x2 * (*x1 - x2));
1165     goto L16;
1166 L15:
1167     if (k > 0) {
1168         *d2 = 0.;
1169     }
1170 L16:
1171     if (*d2 <= small) {
1172         *d2 = small;
1173     }
1174     *x1 = x2;
1175     global_1->fx = fm;
1176     if (sf1 >= global_1->fx) {
1177         goto L17;
1178     }
1179     global_1->fx = sf1;
1180     *x1 = sx1;
1181 /* ...UPDATE X FOR LINEAR BUT NOT PARABOLIC SEARCH... */
1182 L17:
1183     if (j == 0) {
1184         return NLOPT_SUCCESS;
1185     }
1186     i__1 = n;
1187     for (i__ = 1; i__ <= i__1; ++i__) {
1188 /* L18: */
1189         x[i__] += *x1 * q_1->v[i__ + j * n - (n+1)];
1190     }
1191     return NLOPT_SUCCESS;
1192 } /* min_ */
1193
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)
1196 {
1197     /* System generated locals */
1198     int i__1;
1199     double ret_val;
1200
1201     /* Local variables */
1202     nlopt_stopping *stop = q_1->stop;
1203     int i__;
1204     double *t; /* size n */
1205
1206     t = q_1->t_flin;
1207
1208 /* ...FLIN IS THE FUNCTION OF ONE REAL VARIABLE L THAT IS MINIMIZED */
1209 /*   BY THE SUBROUTINE MIN... */
1210     /* Parameter adjustments */
1211     --x;
1212
1213     /* Function Body */
1214     if (j == 0) {
1215         goto L2;
1216     }
1217 /* ...THE SEARCH IS LINEAR... */
1218     i__1 = n;
1219     for (i__ = 1; i__ <= i__1; ++i__) {
1220 /* L1: */
1221         t[i__ - 1] = x[i__] + *l * q_1->v[i__ + j * n - (n+1)];
1222     }
1223     goto L4;
1224 /* ...THE SEARCH IS ALONG A PARABOLIC SPACE CURVE... */
1225 L2:
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));
1229     i__1 = n;
1230     for (i__ = 1; i__ <= i__1; ++i__) {
1231 /* L3: */
1232         t[i__ - 1] = q_1->qa * q_1->q0[i__ - 1] + q_1->qb * x[i__] + q_1->qc *
1233                 q_1->q1[i__ - 1];
1234     }
1235 /* ...THE FUNCTION EVALUATION COUNTER NF IS INCREMENTED... */
1236 L4:
1237     ++(*nf);
1238     ret_val = f(n, t, f_data);
1239     stop->nevals++;
1240     if (ret_val < q_1->fbest) {
1241          q_1->fbest = ret_val;
1242          memcpy(q_1->xbest, t, n * sizeof(double));
1243     }
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;
1248     return ret_val;
1249 } /* flin_ */
1250
1251 static void sort_(int m, int n, double *d__, double *v)
1252 {
1253     /* System generated locals */
1254     int v_dim1, v_offset, i__1, i__2;
1255
1256     /* Local variables */
1257     int i__, j, k;
1258     double s;
1259     int ip1, nm1;
1260
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 */
1265     v_dim1 = m;
1266     v_offset = 1 + v_dim1;
1267     v -= v_offset;
1268     --d__;
1269
1270     /* Function Body */
1271     if (n == 1) {
1272         return;
1273     }
1274     nm1 = n - 1;
1275     i__1 = nm1;
1276     for (i__ = 1; i__ <= i__1; ++i__) {
1277         k = i__;
1278         s = d__[i__];
1279         ip1 = i__ + 1;
1280         i__2 = n;
1281         for (j = ip1; j <= i__2; ++j) {
1282             if (d__[j] <= s) {
1283                 goto L1;
1284             }
1285             k = j;
1286             s = d__[j];
1287 L1:
1288             ;
1289         }
1290         if (k <= i__) {
1291             goto L3;
1292         }
1293         d__[k] = d__[i__];
1294         d__[i__] = s;
1295         i__2 = n;
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];
1299 /* L2: */
1300             v[j + k * v_dim1] = s;
1301         }
1302 L3:
1303         ;
1304     }
1305 } /* sort_ */
1306
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)
1309 {
1310     /* System generated locals */
1311     int i__1;
1312     double d__1;
1313
1314     /* Local variables */
1315     int i__;
1316     double l, s;
1317     double value;
1318
1319 /* ...QUAD LOOKS FOR THE MINIMUM OF F ALONG A CURVE DEFINED BY Q0,Q1,X... */
1320     /* Parameter adjustments */
1321     --x;
1322
1323     /* Function Body */
1324     s = global_1->fx;
1325     global_1->fx = q_1->qf1;
1326     q_1->qf1 = s;
1327     q_1->qd1 = 0.;
1328     i__1 = n;
1329     for (i__ = 1; i__ <= i__1; ++i__) {
1330         s = x[i__];
1331         l = q_1->q1[i__ - 1];
1332         x[i__] = l;
1333         q_1->q1[i__ - 1] = s;
1334 /* L1: */
1335 /* Computing 2nd power */
1336         d__1 = s - l;
1337         q_1->qd1 += d__1 * d__1;
1338     }
1339     q_1->qd1 = sqrt(q_1->qd1);
1340     l = q_1->qd1;
1341     s = 0.;
1342     if (q_1->qd0 <= 0. || q_1->qd1 <= 0. || global_1->nl < n * 3 * n) {
1343         goto L2;
1344     }
1345     value = q_1->qf1;
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));
1351     goto L3;
1352 L2:
1353     global_1->fx = q_1->qf1;
1354     q_1->qa = 0.;
1355     q_1->qb = q_1->qa;
1356     q_1->qc = 1.;
1357 L3:
1358     q_1->qd0 = q_1->qd1;
1359     i__1 = n;
1360     for (i__ = 1; i__ <= i__1; ++i__) {
1361         s = q_1->q0[i__ - 1];
1362         q_1->q0[i__ - 1] = x[i__];
1363 /* L4: */
1364         x[i__] = q_1->qa * s + q_1->qb * x[i__] + q_1->qc * q_1->q1[i__ - 1];
1365     }
1366 } /* quad_ */