luksan/plis.c

   1 #include <limits.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <string.h>
   5 #include "luksan.h"
   6
   7 #define max(a,b) ((a) > (b) ? (a) : (b))
   8 #define min(a,b) ((a) < (b) ? (a) : (b))
   9
  10 /* Common Block Declarations */
  11
  12 typedef struct {
  13      int nres, ndec, nin, nit;
  14      /* int nfv;   -- now stored in stop->nevals */
  15      int nfg, nfh;
  16 } stat_common;
  17
  18 /* *********************************************************************** */
  19 /* SUBROUTINE PLIS               ALL SYSTEMS                   01/09/22 */
  20 /* PURPOSE : */
  21 /* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT */
  22 /* USE THE LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
  23 /* RECURRENCES. */
  24
  25 /* PARAMETERS : */
  26 /*  II  NF  NUMBER OF VARIABLES. */
  27 /*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. */
  28 /*         NB>0-SIMPLE BOUNDS ACCEPTED. */
  29 /*  RI  X(NF)  VECTOR OF VARIABLES. */
  30 /*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE */
  31 /*         X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I). */
  32 /*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND */
  33 /*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. */
  34 /*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. */
  35 /*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. */
  36 /*  RO  GF(NF)  GRADIENT OF THE OBJECTIVE FUNCTION. */
  37 /*  RO  S(NF)  DIRECTION VECTOR. */
  38 /*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE. */
  39 /*  RI  GO(NF)  GRADIENTS DIFFERENCE. */
  40 /*  RA  UO(NF)  AUXILIARY VECTOR. */
  41 /*  RA  VO(NF)  AUXILIARY VECTOR. */
  42 /*  RI  XMAX  MAXIMUM STEPSIZE. */
  43 /*  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES. */
  44 /*  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES. */
  45 /*  RI  TOLB  TOLERANCE FOR THE FUNCTION VALUE. */
  46 /*  RI  TOLG  TOLERANCE FOR THE GRADIENT NORM. */
  47 /*  RI  MINF_EST  ESTIMATION OF THE MINIMUM FUNCTION VALUE. */
  48 /*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE. */
  49 /*  RO  F  VALUE OF THE OBJECTIVE FUNCTION. */
  50 /*  II  MIT  MAXIMUM NUMBER OF ITERATIONS. */
  51 /*  II  MFV  MAXIMUM NUMBER OF FUNCTION EVALUATIONS. */
  52 /*  II  IEST  ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. */
  53 /*         IEST=1-MINIMUM IS ESTIMATED BY THE VALUE MINF_EST. */
  54 /*  II  MF  NUMBER OF LIMITED MEMORY STEPS. */
  55 /*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. */
  56 /*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN */
  57 /*                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
  58 /*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN */
  59 /*                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
  60 /*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. */
  61 /*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. */
  62 /*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, */
  63 /*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. */
  64 /*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. */
  65 /*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. */
  66
  67 /* VARIABLES IN COMMON /STAT/ (STATISTICS) : */
  68 /*  IO  NRES  NUMBER OF RESTARTS. */
  69 /*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION. */
  70 /*  IO  NIN  NUMBER OF INNER ITERATIONS. */
  71 /*  IO  NIT  NUMBER OF ITERATIONS. */
  72 /*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS. */
  73 /*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS. */
  74 /*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS. */
  75
  76 /* SUBPROGRAMS USED : */
  77 /*  S   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS. */
  78 /*  S   PS1L01  STEPSIZE SELECTION USING LINE SEARCH. */
  79 /*  S   PYADC0  ADDITION OF A BOX CONSTRAINT. */
  80 /*  S   PYFUT1  TEST ON TERMINATION. */
  81 /*  S   PYRMC0  DELETION OF A BOX CONSTRAINT. */
  82 /*  S   PYTRCD  COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC */
  83 /*         UPDATE. */
  84 /*  S   PYTRCG  COMPUTATION OF THE PROJECTED GRADIENT. */
  85 /*  S   PYTRCS  COMPUTATION OF THE PROJECTED DIRECTION VECTOR. */
  86 /*  S   MXDRCB BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
  87 /*         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
  88 /*  S   MXDRCF FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
  89 /*         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
  90 /*  S   MXDRSU SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B. */
  91 /*         SHIFT OF ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN */
  92 /*         THE LIMITED MEMORY BFGS METHOD. */
  93 /*  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR. */
  94 /*  RF  MXUDOT  DOT PRODUCT OF TWO VECTORS. */
  95 /*  S   MXUNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN. */
  96 /*  S   MXVCOP  COPYING OF A VECTOR. */
  97 /*  S   MXVSCL  SCALING OF A VECTOR. */
  98
  99 /* EXTERNAL SUBROUTINES : */
 100 /*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. */
 101 /*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER */
 102 /*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS */
 103 /*         THE VALUE OF THE OBJECTIVE FUNCTION. */
 104 /*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. */
 105 /*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER */
 106 /*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF) */
 107 /*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION. */
 108 /* -- OBJ and DOBJ are replaced by a single function, objgrad, in NLopt */
 109
 110 /* METHOD : */
 111 /* LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
 112 /* RECURRENCES. */
 113
 114 static void plis_(int *nf, int *nb, double *x, int *
 115                   ix, double *xl, double *xu, double *gf, double *s,
 116                   double *xo, double *go, double *uo, double *vo,
 117                   double *xmax, double *tolx, double *tolf, double *
 118                   tolb, double *tolg, nlopt_stopping *stop,
 119                   double *minf_est, double *gmax,
 120                   double *f, int *mit, int *mfv, int *iest, int *mf,
 121                   int *iterm, stat_common *stat_1,
 122                   nlopt_func objgrad, void *objgrad_data)
 123 {
 124     /* System generated locals */
 125     int i__1;
 126     double d__1, d__2;
 127
 128     /* Builtin functions */
 129
 130     /* Local variables */
 131     double a, b;
 132     int i__, k, n;
 133     double p, r__;
 134     int kd, ld;
 135     double fo, fp, po, pp, ro, rp;
 136     int kbf, mfg;
 137     int mes, kit;
 138     double alf1, alf2, eta0, eta9, par1, par2;
 139     int mes1, mes2, mes3;
 140     double eps8, eps9;
 141     int mred, iold, nred;
 142     double maxf, dmax__;
 143     int inew;
 144     double told;
 145     int ites;
 146     double rmin, rmax, umax, tolp, tols;
 147     int isys;
 148     int ires1, ires2;
 149     int iterd, mtesf, ntesf;
 150     double gnorm;
 151     int iters, irest, inits, kters, maxst;
 152     double snorm;
 153     int mtesx, ntesx;
 154
 155 /*     INITIATION */
 156
 157     /* Parameter adjustments */
 158     --vo;
 159     --uo;
 160     --go;
 161     --xo;
 162     --s;
 163     --gf;
 164     --xu;
 165     --xl;
 166     --ix;
 167     --x;
 168
 169     /* Function Body */
 170     kbf = 0;
 171     if (*nb > 0) {
 172         kbf = 2;
 173     }
 174     stat_1->nres = 0;
 175     stat_1->ndec = 0;
 176     stat_1->nin = 0;
 177     stat_1->nit = 0;
 178     stat_1->nfg = 0;
 179     stat_1->nfh = 0;
 180     isys = 0;
 181     ites = 1;
 182     mtesx = 2;
 183     mtesf = 2;
 184     inits = 2;
 185     *iterm = 0;
 186     iterd = 0;
 187     iters = 2;
 188     kters = 3;
 189     irest = 0;
 190     ires1 = 999;
 191     ires2 = 0;
 192     mred = 10;
 193     mes = 4;
 194     mes1 = 2;
 195     mes2 = 2;
 196     mes3 = 2;
 197     eta0 = 1e-15;
 198     eta9 = 1e120;
 199     eps8 = 1.;
 200     eps9 = 1e-8;
 201     alf1 = 1e-10;
 202     alf2 = 1e10;
 203     rmax = eta9;
 204     dmax__ = eta9;
 205     maxf = 1e20;
 206     if (*iest <= 0) {
 207         *minf_est = -HUGE_VAL;
 208     }
 209     if (*iest > 0) {
 210         *iest = 1;
 211     }
 212     if (*xmax <= 0.) {
 213         *xmax = 1e16;
 214     }
 215     if (*tolx <= 0.) {
 216         *tolx = 1e-16;
 217     }
 218     if (*tolf <= 0.) {
 219         *tolf = 1e-14;
 220     }
 221     if (*tolg <= 0.) {
 222          *tolg = 1e-8; /* SGJ: was 1e-6, but this sometimes stops too soon */
 223     }
 224 #if 0
 225     /* removed by SGJ: this check prevented us from using minf_max <= 0,
 226        which doesn't make sense.  Instead, if you don't want to have a
 227        lower limit, you should set minf_max = -HUGE_VAL */
 228     if (*tolb <= 0.) {
 229         *tolb = *minf_est + 1e-16;
 230     }
 231 #endif
 232     told = 1e-4;
 233     tols = 1e-4;
 234     tolp = .8;
 235     /* changed by SGJ: default is no limit (INT_MAX) on # iterations/fevals */
 236     if (*mit <= 0) {
 237         *mit = INT_MAX;
 238     }
 239     if (*mfv <= 0) {
 240         *mfv = INT_MAX;
 241     }
 242     mfg = *mfv;
 243     kd = 1;
 244     ld = -1;
 245     kit = -(ires1 * *nf + ires2);
 246     fo = *minf_est;
 247
 248 /*     INITIAL OPERATIONS WITH SIMPLE BOUNDS */
 249
 250     if (kbf > 0) {
 251         i__1 = *nf;
 252         for (i__ = 1; i__ <= i__1; ++i__) {
 253             if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
 254                 xu[i__] = xl[i__];
 255                 ix[i__] = 5;
 256             } else if (ix[i__] == 5 || ix[i__] == 6) {
 257                 xl[i__] = x[i__];
 258                 xu[i__] = x[i__];
 259                 ix[i__] = 5;
 260             }
 261 /* L2: */
 262         }
 263         luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
 264         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
 265     }
 266     if (*iterm != 0) {
 267         goto L11190;
 268     }
 269     *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
 270     ++stop->nevals;
 271     ++stat_1->nfg;
 272 L11120:
 273     luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
 274     luksan_pyfut1__(nf, f, &fo, &umax, gmax, &dmax__, tolx, tolf, tolb, tolg,
 275             &kd, &stat_1->nit, &kit, mit, &stop->nevals, mfv, &stat_1->nfg, &mfg,
 276             &ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
 277             iters, iterm);
 278     if (*iterm != 0) {
 279         goto L11190;
 280     }
 281     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
 282     if (kbf > 0 && rmax > 0.) {
 283         luksan_pyrmc0__(nf, &n, &ix[1], &gf[1], &eps8, &umax, gmax, &rmax, &
 284                 iold, &irest);
 285     }
 286 L11130:
 287
 288 /*     DIRECTION DETERMINATION */
 289
 290     gnorm = sqrt(luksan_mxudot__(nf, &gf[1], &gf[1], &ix[1], &kbf));
 291     if (irest != 0) {
 292         goto L12620;
 293     }
 294 /* Computing MIN */
 295     i__1 = stat_1->nit - kit;
 296     k = min(i__1,*mf);
 297     if (k <= 0) {
 298         irest = max(irest,1);
 299         goto L12620;
 300     }
 301
 302 /*     DETERMINATION OF THE PARAMETER B */
 303
 304     b = luksan_mxudot__(nf, &xo[1], &go[1], &ix[1], &kbf);
 305     if (b <= 0.) {
 306         irest = max(irest,1);
 307         goto L12620;
 308     }
 309     uo[1] = 1. / b;
 310     luksan_mxuneg__(nf, &gf[1], &s[1], &ix[1], &kbf);
 311     luksan_mxdrcb__(nf, &k, &xo[1], &go[1], &uo[1], &vo[1], &s[1], &ix[1], &
 312             kbf);
 313     a = luksan_mxudot__(nf, &go[1], &go[1], &ix[1], &kbf);
 314     if (a > 0.) {
 315         d__1 = b / a;
 316         luksan_mxvscl__(nf, &d__1, &s[1], &s[1]);
 317     }
 318     luksan_mxdrcf__(nf, &k, &xo[1], &go[1], &uo[1], &vo[1], &s[1], &ix[1], &
 319             kbf);
 320     snorm = sqrt(luksan_mxudot__(nf, &s[1], &s[1], &ix[1], &kbf));
 321 /* Computing MIN */
 322     i__1 = k + 1;
 323     k = min(i__1,*mf);
 324     luksan_mxdrsu__(nf, &k, &xo[1], &go[1], &uo[1]);
 325 L12620:
 326     iterd = 0;
 327     if (irest != 0) {
 328
 329 /*     STEEPEST DESCENT DIRECTION */
 330
 331         luksan_mxuneg__(nf, &gf[1], &s[1], &ix[1], &kbf);
 332         snorm = gnorm;
 333         if (kit < stat_1->nit) {
 334             ++stat_1->nres;
 335             kit = stat_1->nit;
 336         } else {
 337              *iterm = -10;
 338             if (iters < 0) {
 339                 *iterm = iters - 5;
 340             }
 341         }
 342     }
 343
 344 /*     TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH */
 345
 346     if (kd > 0) {
 347         p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
 348     }
 349     if (iterd < 0) {
 350         *iterm = iterd;
 351     } else {
 352
 353 /*     TEST ON DESCENT DIRECTION */
 354
 355         if (snorm <= 0.) {
 356             irest = max(irest,1);
 357         } else if (p + told * gnorm * snorm <= 0.) {
 358             irest = 0;
 359         } else {
 360
 361 /*     UNIFORM DESCENT CRITERION */
 362
 363             irest = max(irest,1);
 364         }
 365         if (irest == 0) {
 366
 367 /*     PREPARATION OF LINE SEARCH */
 368
 369             nred = 0;
 370             rmin = alf1 * gnorm / snorm;
 371 /* Computing MIN */
 372             d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
 373             rmax = min(d__1,d__2);
 374         }
 375     }
 376     if (*iterm != 0) {
 377         goto L11190;
 378     }
 379     if (irest != 0) {
 380         goto L11130;
 381     }
 382     luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
 383              &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
 384     if (rmax == 0.) {
 385         goto L11175;
 386     }
 387 L11170:
 388     luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, minf_est, &maxf, &rmin,
 389             &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1->nit, &kit, &
 390             nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
 391     if (isys == 0) {
 392         goto L11174;
 393     }
 394     luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
 395     luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
 396     *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
 397     ++stop->nevals;
 398     ++stat_1->nfg;
 399     p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
 400     goto L11170;
 401 L11174:
 402     if (iters <= 0) {
 403         r__ = 0.;
 404         *f = fo;
 405         p = po;
 406         luksan_mxvcop__(nf, &xo[1], &x[1]);
 407         luksan_mxvcop__(nf, &go[1], &gf[1]);
 408         irest = max(irest,1);
 409         ld = kd;
 410         goto L11130;
 411     }
 412     luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
 413             p, &po, &dmax__, &kbf, &kd, &ld, &iters);
 414 L11175:
 415     if (kbf > 0) {
 416         luksan_mxvine__(nf, &ix[1]);
 417         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
 418     }
 419     goto L11120;
 420 L11190:
 421     return;
 422 } /* plis_ */
 423
 424 /* number of double variables that can be stored in scratch memory
 425    ... it's >= 2007, and this is in the context of scientific computation,
 426    so assume that at least 10M are available, and that sizeof(double)==8 */
 427 #define MEMAVAIL 1310720
 428
 429 /* NLopt wrapper around plis_, handling dynamic allocation etc. */
 430 nlopt_result luksan_plis(int n, nlopt_func f, void *f_data,
 431                   const double *lb, const double *ub, /* bounds */
 432                   double *x, /* in: initial guess, out: minimizer */
 433                   double *minf,
 434                   nlopt_stopping *stop)
 435 {
 436      int i, *ix;
 437      double *work, *xl, *xu, *xo, *gf, *s, *go, *uo, *vo;
 438      double gmax, minf_est;
 439      double xmax = 0; /* no maximum */
 440      double tolg = 0; /* default gradient tolerance */
 441      int iest = 0; /* we have no estimate of min function value */
 442      int mit = 0; /* default no limit on #iterations */
 443      int mfv = stop->maxeval;
 444      stat_common stat;
 445      int iterm;
 446      int mf;
 447
 448      ix = (int*) malloc(sizeof(int) * n);
 449      if (!ix) return NLOPT_OUT_OF_MEMORY;
 450
 451      /* FIXME: what should we set mf to?  The example program tlis.for
 452         sets it to zero as far as I can tell, but it seems to greatly
 453         improve convergence to make it > 0.  The computation time
 454         per iteration, and of course the memory, seem to go as O(n * mf),
 455         and we'll assume that the main limiting factor is the memory.
 456         We'll assume that at least MEMAVAIL memory, or 4*n memory, whichever
 457         is bigger, is available. */
 458      mf = max(MEMAVAIL/n, 4);
 459      if (stop->maxeval && stop->maxeval <= mf)
 460           mf = max(stop->maxeval - 5, 1); /* mf > maxeval seems not good */
 461
 462  retry_alloc:
 463      work = (double*) malloc(sizeof(double) * (n * 4 + max(n,n*mf)*2 +
 464                                                max(n,mf)*2));
 465      if (!work) {
 466           if (mf > 0) {
 467                mf /= 4;
 468                goto retry_alloc;
 469           }
 470           free(ix);
 471           return NLOPT_OUT_OF_MEMORY;
 472      }
 473
 474      xl = work; xu = xl + n; gf = xu + n; s = gf + n;
 475      xo = s + n; go = xo + max(n,n*mf);
 476      uo = go + max(n,n*mf); vo = uo + max(n,mf);
 477
 478      for (i = 0; i < n; ++i) {
 479           int lbu = lb[i] <= -0.99 * HUGE_VAL; /* lb unbounded */
 480           int ubu = ub[i] >= 0.99 * HUGE_VAL;  /* ub unbounded */
 481           ix[i] = lbu ? (ubu ? 0 : 2) : (ubu ? 1 : (lb[i] == ub[i] ? 5 : 3));
 482           xl[i] = lb[i];
 483           xu[i] = ub[i];
 484      }
 485
 486      /* ?  xo does not seem to be initialized in the
 487         original Fortran code, but it is used upon
 488         input to plis if mf > 0 ... perhaps ALLOCATE initializes
 489         arrays to zero by default? */
 490      memset(xo, 0, sizeof(double) * max(n,n*mf));
 491
 492      plis_(&n, &n, x, ix, xl, xu,
 493            gf, s, xo, go, uo, vo,
 494            &xmax,
 495
 496            /* fixme: pass tol_rel and tol_abs and use NLopt check */
 497            &stop->xtol_rel,
 498            &stop->ftol_rel,
 499            &stop->minf_max,
 500            &tolg,
 501            stop,
 502
 503            &minf_est, &gmax,
 504            minf,
 505            &mit, &mfv,
 506            &iest,
 507            &mf,
 508            &iterm, &stat,
 509            f, f_data);
 510
 511      free(work);
 512      free(ix);
 513
 514      switch (iterm) {
 515          case 1: return NLOPT_XTOL_REACHED;
 516          case 2: return NLOPT_FTOL_REACHED;
 517          case 3: return NLOPT_MINF_MAX_REACHED;
 518          case 4: return NLOPT_SUCCESS; /* gradient tolerance reached */
 519          case 6: return NLOPT_SUCCESS;
 520          case 12: case 13: return NLOPT_MAXEVAL_REACHED;
 521          default: return NLOPT_FAILURE;
 522      }
 523 }