luksan/pnet.c

   1 #include <limits.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <string.h>
   5 #include "luksan.h"
   6
   7 #define MAX2(a,b) ((a) > (b) ? (a) : (b))
   8 #define MIN2(a,b) ((a) < (b) ? (a) : (b))
   9
  10 /* Table of constant values */
  11
  12 static double c_b7 = 0.;
  13
  14 /* *********************************************************************** */
  15 /* SUBROUTINE PNET               ALL SYSTEMS                   01/09/22 */
  16 /* PURPOSE : */
  17 /* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT */
  18 /* USE THE LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
  19 /* RECURRENCES. */
  20
  21 /* PARAMETERS : */
  22 /*  II  NF  NUMBER OF VARIABLES. */
  23 /*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. */
  24 /*         NB>0-SIMPLE BOUNDS ACCEPTED. */
  25 /*  RI  X(NF)  VECTOR OF VARIABLES. */
  26 /*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE */
  27 /*         X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I). */
  28 /*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND */
  29 /*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. */
  30 /*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. */
  31 /*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. */
  32 /*  RO  GF(NF)  GRADIENT OF THE OBJECTIVE FUNCTION. */
  33 /*  RA  GN(NF)  OLD GRADIENT OF THE OBJECTIVE FUNCTION. */
  34 /*  RO  S(NF)  DIRECTION VECTOR. */
  35 /*  RA  XO(NF)  ARRAY CONTAINING INCREMENTS OF VARIABLES. */
  36 /*  RA  GO(NF)  ARRAY CONTAINING INCREMENTS OF GRADIENTS. */
  37 /*  RA  XS(NF)  AUXILIARY VECTOR. */
  38 /*  RA  GS(NF)  AUXILIARY VECTOR. */
  39 /*  RA  XM(NF*MF)  ARRAY CONTAINING INCREMENTS OF VARIABLES. */
  40 /*  RA  GM(NF*MF)  ARRAY CONTAINING INCREMENTS OF GRADIENTS. */
  41 /*  RA  U1(MF)  AUXILIARY VECTOR. */
  42 /*  RA  U2(MF)  AUXILIARY VECTOR. */
  43 /*  RI  XMAX  MAXIMUM STEPSIZE. */
  44 /*  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES. */
  45 /*  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES. */
  46 /*  RI  TOLB  TOLERANCE FOR THE FUNCTION VALUE. */
  47 /*  RI  TOLG  TOLERANCE FOR THE GRADIENT NORM. */
  48 /*  RI  MINF_EST  ESTIMATION OF THE MINIMUM FUNCTION VALUE. */
  49 /*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE. */
  50 /*  RO  F  VALUE OF THE OBJECTIVE FUNCTION. */
  51 /*  II  MIT  MAXIMUM NUMBER OF ITERATIONS. */
  52 /*  II  MFV  MAXIMUM NUMBER OF FUNCTION EVALUATIONS. */
  53 /*  II  MFG  MAXIMUM NUMBER OF GRADIENT EVALUATIONS. */
  54 /*  II  IEST  ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. */
  55 /*         IEST=1-MINIMUM IS ESTIMATED BY THE VALUE MINF_EST. */
  56 /*  II  MOS1  CHOICE OF RESTARTS AFTER A CONSTRAINT CHANGE. */
  57 /*         MOS1=1-RESTARTS ARE SUPPRESSED. MOS1=2-RESTARTS WITH */
  58 /*         STEEPEST DESCENT DIRECTIONS ARE USED. */
  59 /*  II  MOS1  CHOICE OF DIRECTION VECTORS AFTER RESTARTS. MOS1=1-THE */
  60 /*         NEWTON DIRECTIONS ARE USED. MOS1=2-THE STEEPEST DESCENT */
  61 /*         DIRECTIONS ARE USED. */
  62 /*  II  MOS2  CHOICE OF PRECONDITIONING STRATEGY. MOS2=1-PRECONDITIONING */
  63 /*         IS NOT USED. MOS2=2-PRECONDITIONING BY THE LIMITED MEMORY */
  64 /*         BFGS METHOD IS USED. */
  65 /*  II  MF  THE NUMBER OF LIMITED-MEMORY VARIABLE METRIC UPDATES */
  66 /*         IN EACH ITERATION (THEY USE 2*MF STORED VECTORS). */
  67 /*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. */
  68 /*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN */
  69 /*                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
  70 /*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN */
  71 /*                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
  72 /*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. */
  73 /*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. */
  74 /*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, */
  75 /*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. */
  76 /*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. */
  77 /*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. */
  78
  79 /* VARIABLES IN COMMON /STAT/ (STATISTICS) : */
  80 /*  IO  NRES  NUMBER OF RESTARTS. */
  81 /*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION. */
  82 /*  IO  NIN  NUMBER OF INNER ITERATIONS. */
  83 /*  IO  NIT  NUMBER OF ITERATIONS. */
  84 /*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS. */
  85 /*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS. */
  86 /*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS. */
  87
  88 /* SUBPROGRAMS USED : */
  89 /*  S   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS. */
  90 /*  S   PS1L01  STEPSIZE SELECTION USING LINE SEARCH. */
  91 /*  S   PYADC0  ADDITION OF A BOX CONSTRAINT. */
  92 /*  S   PYFUT1  TEST ON TERMINATION. */
  93 /*  S   PYRMC0  DELETION OF A BOX CONSTRAINT. */
  94 /*  S   PYTRCD  COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC */
  95 /*         UPDATE. */
  96 /*  S   PYTRCG  COMPUTATION OF THE PROJECTED GRADIENT. */
  97 /*  S   PYTRCS  COMPUTATION OF THE PROJECTED DIRECTION VECTOR. */
  98 /*  S   MXDRCB BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
  99 /*         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
 100 /*  S   MXDRCF FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
 101 /*         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
 102 /*  S   MXDRSU SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B. */
 103 /*         SHIFT OF ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN */
 104 /*         THE LIMITED MEMORY BFGS METHOD. */
 105 /*  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR. */
 106 /*  RF  MXUDOT  DOT PRODUCT OF TWO VECTORS. */
 107 /*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN. */
 108 /*  S   MXVCOP  COPYING OF A VECTOR. */
 109 /*  S   MXVSCL  SCALING OF A VECTOR. */
 110 /*  S   MXVSET  INITIATINON OF A VECTOR. */
 111 /*  S   MXVDIF  DIFFERENCE OF TWO VECTORS. */
 112
 113 /* EXTERNAL SUBROUTINES : */
 114 /*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. */
 115 /*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER */
 116 /*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS */
 117 /*         THE VALUE OF THE OBJECTIVE FUNCTION. */
 118 /*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. */
 119 /*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER */
 120 /*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF) */
 121 /*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION. */
 122 /* -- OBJ and DOBJ are replaced by a single function, objgrad, in NLopt */
 123
 124 /* METHOD : */
 125 /* LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
 126 /* RECURRENCES. */
 127
 128 static void pnet_(int *nf, int *nb, double *x, int *
 129                   ix, double *xl, double *xu, double *gf, double *gn,
 130                   double *s, double *xo, double *go, double *xs,
 131                   double *gs, double *xm, double *gm, double *u1,
 132                   double *u2, double *xmax, double *tolx, double *tolf,
 133                   double *tolb, double *tolg, nlopt_stopping *stop,
 134                   double *minf_est, double *
 135                   gmax, double *f, int *mit, int *mfv, int *mfg,
 136                   int *iest, int *mos1, int *mos2, int *mf,
 137                   int *iterm, stat_common *stat_1,
 138                   nlopt_func objgrad, void *objgrad_data)
 139 {
 140     /* System generated locals */
 141     int i__1;
 142     double d__1, d__2;
 143
 144     /* Builtin functions */
 145
 146     /* Local variables */
 147     double a, b;
 148     int i__, n;
 149     double p, r__;
 150     int kd, ld;
 151     double fo, fp, po, pp, ro, rp;
 152     int mx, kbf;
 153     double alf;
 154     double par;
 155     int mes, kit;
 156     double rho, eps;
 157     int mmx;
 158     double alf1, alf2, eta0, eta9, par1, par2;
 159     int mes1, mes2, mes3;
 160     double rho1, rho2, eps8, eps9;
 161     int mred, iold, nred;
 162     double maxf, dmax__;
 163     int xstop = 0;
 164     int inew;
 165     double told;
 166     int ites;
 167     double rmin, rmax, umax, tolp, tols;
 168     int isys;
 169     int ires1, ires2;
 170     int iterd, mtesf, ntesf;
 171     double gnorm;
 172     int iters, irest, inits, kters, maxst;
 173     double snorm;
 174     int mtesx, ntesx;
 175
 176 /*     INITIATION */
 177
 178     /* Parameter adjustments */
 179     --u2;
 180     --u1;
 181     --gm;
 182     --xm;
 183     --gs;
 184     --xs;
 185     --go;
 186     --xo;
 187     --s;
 188     --gn;
 189     --gf;
 190     --xu;
 191     --xl;
 192     --ix;
 193     --x;
 194
 195     /* Function Body */
 196     kbf = 0;
 197     if (*nb > 0) {
 198         kbf = 2;
 199     }
 200     stat_1->nres = 0;
 201     stat_1->ndec = 0;
 202     stat_1->nin = 0;
 203     stat_1->nit = 0;
 204     stat_1->nfg = 0;
 205     stat_1->nfh = 0;
 206     isys = 0;
 207     ites = 1;
 208     mtesx = 2;
 209     mtesf = 2;
 210     inits = 2;
 211     *iterm = 0;
 212     iterd = 0;
 213     iters = 2;
 214     kters = 3;
 215     irest = 0;
 216     ires1 = 999;
 217     ires2 = 0;
 218     mred = 10;
 219     mes = 4;
 220     mes1 = 2;
 221     mes2 = 2;
 222     mes3 = 2;
 223     eps = .8;
 224     eta0 = 1e-15;
 225     eta9 = 1e120;
 226     eps8 = 1.;
 227     eps9 = 1e-8;
 228     alf1 = 1e-10;
 229     alf2 = 1e10;
 230     rmax = eta9;
 231     dmax__ = eta9;
 232     maxf = 1e20;
 233     if (*iest <= 0) {
 234          *minf_est = -HUGE_VAL; /* changed from -1e60 by SGJ */
 235     }
 236     if (*iest > 0) {
 237         *iest = 1;
 238     }
 239     if (*xmax <= 0.) {
 240         *xmax = 1e16;
 241     }
 242     if (*tolx <= 0.) {
 243         *tolx = 1e-16;
 244     }
 245     if (*tolf <= 0.) {
 246         *tolf = 1e-14;
 247     }
 248     if (*tolg <= 0.) {
 249          *tolg = 1e-8; /* SGJ: was 1e-6, but this sometimes stops too soon */
 250     }
 251 #if 0
 252     /* removed by SGJ: this check prevented us from using minf_max <= 0,
 253        which doesn't make sense.  Instead, if you don't want to have a
 254        lower limit, you should set minf_max = -HUGE_VAL */
 255     if (*tolb <= 0.) {
 256         *tolb = *minf_est + 1e-16;
 257     }
 258 #endif
 259     told = 1e-4;
 260     tols = 1e-4;
 261     tolp = .9;
 262     /* changed by SGJ: default is no limit (INT_MAX) on # iterations/fevals */
 263     if (*mit <= 0) {
 264         *mit = INT_MAX;
 265     }
 266     if (*mfv <= 0) {
 267         *mfv = INT_MAX;
 268     }
 269     if (*mfg <= 0) {
 270         *mfg = INT_MAX;
 271     }
 272     if (*mos1 <= 0) {
 273         *mos1 = 1;
 274     }
 275     if (*mos2 <= 0) {
 276         *mos2 = 1;
 277     }
 278     kd = 1;
 279     ld = -1;
 280     kit = -(ires1 * *nf + ires2);
 281     fo = *minf_est;
 282
 283 /*     INITIAL OPERATIONS WITH SIMPLE BOUNDS */
 284
 285     if (kbf > 0) {
 286         i__1 = *nf;
 287         for (i__ = 1; i__ <= i__1; ++i__) {
 288             if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
 289                 xu[i__] = xl[i__];
 290                 ix[i__] = 5;
 291             } else if (ix[i__] == 5 || ix[i__] == 6) {
 292                 xl[i__] = x[i__];
 293                 xu[i__] = x[i__];
 294                 ix[i__] = 5;
 295             }
 296 /* L2: */
 297         }
 298         luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
 299         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
 300     }
 301     *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
 302     ++stop->nevals;
 303     ++stat_1->nfg;
 304     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11080; }
 305     ld = kd;
 306 L11020:
 307     luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
 308     luksan_mxvcop__(nf, &gf[1], &gn[1]);
 309     luksan_pyfut1__(nf, f, &fo, &umax, gmax, xstop, stop, tolg,
 310             &kd, &stat_1->nit, &kit, mit, &stat_1->nfg, mfg, &
 311             ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
 312             iters, iterm);
 313     if (*iterm != 0) {
 314         goto L11080;
 315     }
 316     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11080; }
 317     if (kbf > 0) {
 318         luksan_pyrmc0__(nf, &n, &ix[1], &gn[1], &eps8, &umax, gmax, &rmax, &
 319                 iold, &irest);
 320         if (umax > eps8 * *gmax) {
 321             irest = MAX2(irest,1);
 322         }
 323     }
 324     luksan_mxvcop__(nf, &x[1], &xo[1]);
 325 L11040:
 326
 327 /*     DIRECTION DETERMINATION */
 328
 329     if (irest != 0) {
 330         if (kit < stat_1->nit) {
 331             mx = 0;
 332             ++stat_1->nres;
 333             kit = stat_1->nit;
 334         } else {
 335             *iterm = -10;
 336             if (iters < 0) {
 337                 *iterm = iters - 5;
 338             }
 339             goto L11080;
 340         }
 341         if (*mos1 > 1) {
 342             luksan_mxvneg__(nf, &gn[1], &s[1]);
 343             gnorm = sqrt(luksan_mxudot__(nf, &gn[1], &gn[1], &ix[1], &kbf));
 344             snorm = gnorm;
 345             goto L12560;
 346         }
 347     }
 348     rho1 = luksan_mxudot__(nf, &gn[1], &gn[1], &ix[1], &kbf);
 349     gnorm = sqrt(rho1);
 350 /* Computing MIN */
 351     d__1 = eps, d__2 = sqrt(gnorm);
 352     par = MIN2(d__1,d__2);
 353     if (par > .01) {
 354 /* Computing MIN */
 355         d__1 = par, d__2 = 1. / (double) stat_1->nit;
 356         par = MIN2(d__1,d__2);
 357     }
 358     par *= par;
 359
 360 /*     CG INITIATION */
 361
 362     rho = rho1;
 363     snorm = 0.;
 364     luksan_mxvset__(nf, &c_b7, &s[1]);
 365     luksan_mxvneg__(nf, &gn[1], &gs[1]);
 366     luksan_mxvcop__(nf, &gs[1], &xs[1]);
 367     if (*mos2 > 1) {
 368         if (mx == 0) {
 369             b = 0.;
 370         } else {
 371             b = luksan_mxudot__(nf, &xm[1], &gm[1], &ix[1], &kbf);
 372         }
 373         if (b > 0.) {
 374             u1[1] = 1. / b;
 375             luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
 376                     ix[1], &kbf);
 377             a = luksan_mxudot__(nf, &gm[1], &gm[1], &ix[1], &kbf);
 378             if (a > 0.) {
 379                 d__1 = b / a;
 380                 luksan_mxvscl__(nf, &d__1, &xs[1], &xs[1]);
 381             }
 382             luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
 383                     ix[1], &kbf);
 384         }
 385     }
 386     rho = luksan_mxudot__(nf, &gs[1], &xs[1], &ix[1], &kbf);
 387 /*      SIG=RHO */
 388     mmx = *nf + 3;
 389     nred = 0;
 390 L12520:
 391     ++nred;
 392     if (nred > mmx) {
 393         goto L12550;
 394     }
 395     fo = *f;
 396     pp = sqrt(eta0 / luksan_mxudot__(nf, &xs[1], &xs[1], &ix[1], &kbf));
 397     ld = 0;
 398     luksan_mxudir__(nf, &pp, &xs[1], &xo[1], &x[1], &ix[1], &kbf);
 399     objgrad(*nf, &x[1], &gf[1], objgrad_data);
 400     ++stop->nevals;
 401     ++stat_1->nfg;
 402     ld = kd;
 403     luksan_mxvdif__(nf, &gf[1], &gn[1], &go[1]);
 404     *f = fo;
 405     d__1 = 1. / pp;
 406     luksan_mxvscl__(nf, &d__1, &go[1], &go[1]);
 407     alf = luksan_mxudot__(nf, &xs[1], &go[1], &ix[1], &kbf);
 408     if (alf <= 1. / eta9) {
 409 /*      IF (ALF.LE.1.0D-8*SIG) THEN */
 410
 411 /*     CG FAILS (THE MATRIX IS NOT POSITIVE DEFINITE) */
 412
 413         if (nred == 1) {
 414             luksan_mxvneg__(nf, &gn[1], &s[1]);
 415             snorm = gnorm;
 416         }
 417         iterd = 0;
 418         goto L12560;
 419     } else {
 420         iterd = 2;
 421     }
 422
 423 /*     CG STEP */
 424
 425     alf = rho / alf;
 426     luksan_mxudir__(nf, &alf, &xs[1], &s[1], &s[1], &ix[1], &kbf);
 427     d__1 = -alf;
 428     luksan_mxudir__(nf, &d__1, &go[1], &gs[1], &gs[1], &ix[1], &kbf);
 429     rho2 = luksan_mxudot__(nf, &gs[1], &gs[1], &ix[1], &kbf);
 430     snorm = sqrt(luksan_mxudot__(nf, &s[1], &s[1], &ix[1], &kbf));
 431     if (rho2 <= par * rho1) {
 432         goto L12560;
 433     }
 434     if (nred >= mmx) {
 435         goto L12550;
 436     }
 437     if (*mos2 > 1) {
 438         if (b > 0.) {
 439             luksan_mxvcop__(nf, &gs[1], &go[1]);
 440             luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
 441                     ix[1], &kbf);
 442             if (a > 0.) {
 443                 d__1 = b / a;
 444                 luksan_mxvscl__(nf, &d__1, &go[1], &go[1]);
 445             }
 446             luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
 447                     ix[1], &kbf);
 448             rho2 = luksan_mxudot__(nf, &gs[1], &go[1], &ix[1], &kbf);
 449             alf = rho2 / rho;
 450             luksan_mxudir__(nf, &alf, &xs[1], &go[1], &xs[1], &ix[1], &kbf);
 451         } else {
 452             alf = rho2 / rho;
 453             luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
 454         }
 455     } else {
 456         alf = rho2 / rho;
 457         luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
 458     }
 459     rho = rho2;
 460 /*      SIG=RHO2+ALF*ALF*SIG */
 461     goto L12520;
 462 L12550:
 463
 464 /*     AN INEXACT SOLUTION IS OBTAINED */
 465
 466 L12560:
 467
 468 /*     ------------------------------ */
 469 /*     END OF DIRECTION DETERMINATION */
 470 /*     ------------------------------ */
 471
 472     luksan_mxvcop__(nf, &xo[1], &x[1]);
 473     luksan_mxvcop__(nf, &gn[1], &gf[1]);
 474     if (kd > 0) {
 475         p = luksan_mxudot__(nf, &gn[1], &s[1], &ix[1], &kbf);
 476     }
 477     if (iterd < 0) {
 478         *iterm = iterd;
 479     } else {
 480
 481 /*     TEST ON DESCENT DIRECTION */
 482
 483         if (snorm <= 0.) {
 484             irest = MAX2(irest,1);
 485         } else if (p + told * gnorm * snorm <= 0.) {
 486             irest = 0;
 487         } else {
 488
 489 /*     UNIFORM DESCENT CRITERION */
 490
 491             irest = MAX2(irest,1);
 492         }
 493         if (irest == 0) {
 494
 495 /*     PREPARATION OF LINE SEARCH */
 496
 497             nred = 0;
 498             rmin = alf1 * gnorm / snorm;
 499 /* Computing MIN */
 500             d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
 501             rmax = MIN2(d__1,d__2);
 502         }
 503     }
 504     ld = kd;
 505     if (*iterm != 0) {
 506         goto L11080;
 507     }
 508     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11080; }
 509     if (irest != 0) {
 510         goto L11040;
 511     }
 512     luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
 513              &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
 514     if (rmax == 0.) {
 515         goto L11075;
 516     }
 517 L11060:
 518     luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, minf_est, &maxf, &rmin,
 519             &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1->nit, &kit, &
 520             nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
 521     if (isys == 0) {
 522         goto L11064;
 523     }
 524     luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
 525     luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
 526     *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
 527     ++stop->nevals;
 528     ++stat_1->nfg;
 529     ld = kd;
 530     p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
 531     goto L11060;
 532 L11064:
 533     if (iters <= 0) {
 534         r__ = 0.;
 535         *f = fo;
 536         p = po;
 537         luksan_mxvcop__(nf, &xo[1], &x[1]);
 538         luksan_mxvcop__(nf, &go[1], &gf[1]);
 539         irest = MAX2(irest,1);
 540         ld = kd;
 541         goto L11040;
 542     }
 543     luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
 544             p, &po, &dmax__, &kbf, &kd, &ld, &iters);
 545     xstop = nlopt_stop_dx(stop, &x[1], &xo[1]);
 546     if (*mos2 > 1) {
 547 /* Computing MIN */
 548         i__1 = mx + 1;
 549         mx = MIN2(i__1,*mf);
 550         luksan_mxdrsu__(nf, &mx, &xm[1], &gm[1], &u1[1]);
 551         luksan_mxvcop__(nf, &xo[1], &xm[1]);
 552         luksan_mxvcop__(nf, &go[1], &gm[1]);
 553     }
 554 L11075:
 555     if (kbf > 0) {
 556         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
 557         if (inew > 0) {
 558             irest = MAX2(irest,1);
 559         }
 560     }
 561     goto L11020;
 562 L11080:
 563     return;
 564 } /* pnet_ */
 565
 566 /* NLopt wrapper around pnet_, handling dynamic allocation etc. */
 567 nlopt_result luksan_pnet(int n, nlopt_func f, void *f_data,
 568                          const double *lb, const double *ub, /* bounds */
 569                          double *x, /* in: initial guess, out: minimizer */
 570                          double *minf,
 571                          nlopt_stopping *stop,
 572                          int mf, /* subspace dimension (0 for default) */
 573                          int mos1, int mos2) /* 1 or 2 */
 574 {
 575      int i, *ix, nb = 1;
 576      double *work;
 577      double *xl, *xu, *gf, *gn, *s, *xo, *go, *xs, *gs, *xm, *gm, *u1, *u2;
 578      double gmax, minf_est;
 579      double xmax = 0; /* no maximum */
 580      double tolg = 0; /* default gradient tolerance */
 581      int iest = 0; /* we have no estimate of min function value */
 582      int mit = 0, mfg = 0; /* default no limit on #iterations */
 583      int mfv = stop->maxeval;
 584      stat_common stat;
 585      int iterm;
 586
 587      ix = (int*) malloc(sizeof(int) * n);
 588      if (!ix) return NLOPT_OUT_OF_MEMORY;
 589
 590      /* FIXME: what should we set mf to?  The example program tlis.for
 591         sets it to zero as far as I can tell, but it seems to greatly
 592         improve convergence to make it > 0.  The computation time
 593         per iteration, and of course the memory, seem to go as O(n * mf),
 594         and we'll assume that the main limiting factor is the memory.
 595         We'll assume that at least MEMAVAIL memory, or 4*n memory, whichever
 596         is bigger, is available. */
 597      if (mf <= 0) {
 598           mf = MAX2(MEMAVAIL/n, 4);
 599           if (stop->maxeval && stop->maxeval <= mf)
 600                mf = MAX2(stop->maxeval - 5, 1); /* mf > maxeval seems not good */
 601      }
 602
 603  retry_alloc:
 604      work = (double*) malloc(sizeof(double) * (n * 9 + MAX2(n,n*mf)*2 +
 605                                                MAX2(n,mf)*2));
 606      if (!work) {
 607           if (mf > 0) {
 608                mf = 0; /* allocate minimal memory */
 609                goto retry_alloc;
 610           }
 611           free(ix);
 612           return NLOPT_OUT_OF_MEMORY;
 613      }
 614
 615      xl = work; xu = xl + n;
 616      gf = xu + n; gn = gf + n; s = gn + n;
 617      xo = s + n; go = xo + n; xs = go + n; gs = xs + n;
 618      xm = gs + n; gm = xm + MAX2(n*mf,n);
 619      u1 = gm + MAX2(n*mf,n); u2 = u1 + MAX2(n,mf);
 620
 621      for (i = 0; i < n; ++i) {
 622           int lbu = lb[i] <= -0.99 * HUGE_VAL; /* lb unbounded */
 623           int ubu = ub[i] >= 0.99 * HUGE_VAL;  /* ub unbounded */
 624           ix[i] = lbu ? (ubu ? 0 : 2) : (ubu ? 1 : (lb[i] == ub[i] ? 5 : 3));
 625           xl[i] = lb[i];
 626           xu[i] = ub[i];
 627      }
 628
 629      /* ?  xo does not seem to be initialized in the
 630         original Fortran code, but it is used upon
 631         input to pnet if mf > 0 ... perhaps ALLOCATE initializes
 632         arrays to zero by default? */
 633      memset(xo, 0, sizeof(double) * MAX2(n,n*mf));
 634
 635      pnet_(&n, &nb, x, ix, xl, xu,
 636            gf, gn, s, xo, go, xs, gs, xm, gm, u1, u2,
 637            &xmax,
 638
 639            /* fixme: pass tol_rel and tol_abs and use NLopt check */
 640            &stop->xtol_rel,
 641            &stop->ftol_rel,
 642            &stop->minf_max,
 643            &tolg,
 644            stop,
 645
 646            &minf_est, &gmax,
 647            minf,
 648            &mit, &mfv, &mfg,
 649            &iest,
 650            &mos1, &mos2,
 651            &mf,
 652            &iterm, &stat,
 653            f, f_data);
 654
 655      free(work);
 656      free(ix);
 657
 658      switch (iterm) {
 659          case 1: return NLOPT_XTOL_REACHED;
 660          case 2: return NLOPT_FTOL_REACHED;
 661          case 3: return NLOPT_MINF_MAX_REACHED;
 662          case 4: return NLOPT_SUCCESS; /* gradient tolerance reached */
 663          case 6: return NLOPT_SUCCESS;
 664          case 12: case 13: return NLOPT_MAXEVAL_REACHED;
 665          case 100: return NLOPT_MAXTIME_REACHED;
 666          case -999: return NLOPT_FORCED_STOP;
 667          default: return NLOPT_FAILURE;
 668      }
 669 }