chiark / gitweb /
c563c96d165e8151f0de8570df9e4649d7e2563b
[nlopt.git] / 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 max(a,b) ((a) > (b) ? (a) : (b))
8 #define min(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     ld = kd;
305 L11020:
306     luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
307     luksan_mxvcop__(nf, &gf[1], &gn[1]);
308     luksan_pyfut1__(nf, f, &fo, &umax, gmax, xstop, stop, tolg, 
309             &kd, &stat_1->nit, &kit, mit, &stat_1->nfg, mfg, &
310             ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
311             iters, iterm);
312     if (*iterm != 0) {
313         goto L11080;
314     }
315     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11080; }
316     if (kbf > 0) {
317         luksan_pyrmc0__(nf, &n, &ix[1], &gn[1], &eps8, &umax, gmax, &rmax, &
318                 iold, &irest);
319         if (umax > eps8 * *gmax) {
320             irest = max(irest,1);
321         }
322     }
323     luksan_mxvcop__(nf, &x[1], &xo[1]);
324 L11040:
325
326 /*     DIRECTION DETERMINATION */
327
328     if (irest != 0) {
329         if (kit < stat_1->nit) {
330             mx = 0;
331             ++stat_1->nres;
332             kit = stat_1->nit;
333         } else {
334             *iterm = -10;
335             if (iters < 0) {
336                 *iterm = iters - 5;
337             }
338             goto L11080;
339         }
340         if (*mos1 > 1) {
341             luksan_mxvneg__(nf, &gn[1], &s[1]);
342             gnorm = sqrt(luksan_mxudot__(nf, &gn[1], &gn[1], &ix[1], &kbf));
343             snorm = gnorm;
344             goto L12560;
345         }
346     }
347     rho1 = luksan_mxudot__(nf, &gn[1], &gn[1], &ix[1], &kbf);
348     gnorm = sqrt(rho1);
349 /* Computing MIN */
350     d__1 = eps, d__2 = sqrt(gnorm);
351     par = min(d__1,d__2);
352     if (par > .01) {
353 /* Computing MIN */
354         d__1 = par, d__2 = 1. / (double) stat_1->nit;
355         par = min(d__1,d__2);
356     }
357     par *= par;
358
359 /*     CG INITIATION */
360
361     rho = rho1;
362     snorm = 0.;
363     luksan_mxvset__(nf, &c_b7, &s[1]);
364     luksan_mxvneg__(nf, &gn[1], &gs[1]);
365     luksan_mxvcop__(nf, &gs[1], &xs[1]);
366     if (*mos2 > 1) {
367         if (mx == 0) {
368             b = 0.;
369         } else {
370             b = luksan_mxudot__(nf, &xm[1], &gm[1], &ix[1], &kbf);
371         }
372         if (b > 0.) {
373             u1[1] = 1. / b;
374             luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
375                     ix[1], &kbf);
376             a = luksan_mxudot__(nf, &gm[1], &gm[1], &ix[1], &kbf);
377             if (a > 0.) {
378                 d__1 = b / a;
379                 luksan_mxvscl__(nf, &d__1, &xs[1], &xs[1]);
380             }
381             luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
382                     ix[1], &kbf);
383         }
384     }
385     rho = luksan_mxudot__(nf, &gs[1], &xs[1], &ix[1], &kbf);
386 /*      SIG=RHO */
387     mmx = *nf + 3;
388     nred = 0;
389 L12520:
390     ++nred;
391     if (nred > mmx) {
392         goto L12550;
393     }
394     fo = *f;
395     pp = sqrt(eta0 / luksan_mxudot__(nf, &xs[1], &xs[1], &ix[1], &kbf));
396     ld = 0;
397     luksan_mxudir__(nf, &pp, &xs[1], &xo[1], &x[1], &ix[1], &kbf);
398     objgrad(*nf, &x[1], &gf[1], objgrad_data);
399     ++stop->nevals;
400     ++stat_1->nfg;
401     ld = kd;
402     luksan_mxvdif__(nf, &gf[1], &gn[1], &go[1]);
403     *f = fo;
404     d__1 = 1. / pp;
405     luksan_mxvscl__(nf, &d__1, &go[1], &go[1]);
406     alf = luksan_mxudot__(nf, &xs[1], &go[1], &ix[1], &kbf);
407     if (alf <= 1. / eta9) {
408 /*      IF (ALF.LE.1.0D-8*SIG) THEN */
409
410 /*     CG FAILS (THE MATRIX IS NOT POSITIVE DEFINITE) */
411
412         if (nred == 1) {
413             luksan_mxvneg__(nf, &gn[1], &s[1]);
414             snorm = gnorm;
415         }
416         iterd = 0;
417         goto L12560;
418     } else {
419         iterd = 2;
420     }
421
422 /*     CG STEP */
423
424     alf = rho / alf;
425     luksan_mxudir__(nf, &alf, &xs[1], &s[1], &s[1], &ix[1], &kbf);
426     d__1 = -alf;
427     luksan_mxudir__(nf, &d__1, &go[1], &gs[1], &gs[1], &ix[1], &kbf);
428     rho2 = luksan_mxudot__(nf, &gs[1], &gs[1], &ix[1], &kbf);
429     snorm = sqrt(luksan_mxudot__(nf, &s[1], &s[1], &ix[1], &kbf));
430     if (rho2 <= par * rho1) {
431         goto L12560;
432     }
433     if (nred >= mmx) {
434         goto L12550;
435     }
436     if (*mos2 > 1) {
437         if (b > 0.) {
438             luksan_mxvcop__(nf, &gs[1], &go[1]);
439             luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
440                     ix[1], &kbf);
441             if (a > 0.) {
442                 d__1 = b / a;
443                 luksan_mxvscl__(nf, &d__1, &go[1], &go[1]);
444             }
445             luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
446                     ix[1], &kbf);
447             rho2 = luksan_mxudot__(nf, &gs[1], &go[1], &ix[1], &kbf);
448             alf = rho2 / rho;
449             luksan_mxudir__(nf, &alf, &xs[1], &go[1], &xs[1], &ix[1], &kbf);
450         } else {
451             alf = rho2 / rho;
452             luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
453         }
454     } else {
455         alf = rho2 / rho;
456         luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
457     }
458     rho = rho2;
459 /*      SIG=RHO2+ALF*ALF*SIG */
460     goto L12520;
461 L12550:
462
463 /*     AN INEXACT SOLUTION IS OBTAINED */
464
465 L12560:
466
467 /*     ------------------------------ */
468 /*     END OF DIRECTION DETERMINATION */
469 /*     ------------------------------ */
470
471     luksan_mxvcop__(nf, &xo[1], &x[1]);
472     luksan_mxvcop__(nf, &gn[1], &gf[1]);
473     if (kd > 0) {
474         p = luksan_mxudot__(nf, &gn[1], &s[1], &ix[1], &kbf);
475     }
476     if (iterd < 0) {
477         *iterm = iterd;
478     } else {
479
480 /*     TEST ON DESCENT DIRECTION */
481
482         if (snorm <= 0.) {
483             irest = max(irest,1);
484         } else if (p + told * gnorm * snorm <= 0.) {
485             irest = 0;
486         } else {
487
488 /*     UNIFORM DESCENT CRITERION */
489
490             irest = max(irest,1);
491         }
492         if (irest == 0) {
493
494 /*     PREPARATION OF LINE SEARCH */
495
496             nred = 0;
497             rmin = alf1 * gnorm / snorm;
498 /* Computing MIN */
499             d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
500             rmax = min(d__1,d__2);
501         }
502     }
503     ld = kd;
504     if (*iterm != 0) {
505         goto L11080;
506     }
507     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11080; }
508     if (irest != 0) {
509         goto L11040;
510     }
511     luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
512              &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
513     if (rmax == 0.) {
514         goto L11075;
515     }
516 L11060:
517     luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, minf_est, &maxf, &rmin, 
518             &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1->nit, &kit, &
519             nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
520     if (isys == 0) {
521         goto L11064;
522     }
523     luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
524     luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
525     *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
526     ++stop->nevals;
527     ++stat_1->nfg;
528     ld = kd;
529     p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
530     goto L11060;
531 L11064:
532     if (iters <= 0) {
533         r__ = 0.;
534         *f = fo;
535         p = po;
536         luksan_mxvcop__(nf, &xo[1], &x[1]);
537         luksan_mxvcop__(nf, &go[1], &gf[1]);
538         irest = max(irest,1);
539         ld = kd;
540         goto L11040;
541     }
542     luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
543             p, &po, &dmax__, &kbf, &kd, &ld, &iters);
544     xstop = nlopt_stop_dx(stop, &x[1], &xo[1]);
545     if (*mos2 > 1) {
546 /* Computing MIN */
547         i__1 = mx + 1;
548         mx = min(i__1,*mf);
549         luksan_mxdrsu__(nf, &mx, &xm[1], &gm[1], &u1[1]);
550         luksan_mxvcop__(nf, &xo[1], &xm[1]);
551         luksan_mxvcop__(nf, &go[1], &gm[1]);
552     }
553 L11075:
554     if (kbf > 0) {
555         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
556         if (inew > 0) {
557             irest = max(irest,1);
558         }
559     }
560     goto L11020;
561 L11080:
562     return;
563 } /* pnet_ */
564
565 /* NLopt wrapper around pnet_, handling dynamic allocation etc. */
566 nlopt_result luksan_pnet(int n, nlopt_func f, void *f_data,
567                          const double *lb, const double *ub, /* bounds */
568                          double *x, /* in: initial guess, out: minimizer */
569                          double *minf, 
570                          nlopt_stopping *stop,
571                          int mos1, int mos2) /* 1 or 2 */
572 {
573      int i, *ix, nb = 1;
574      double *work;
575      double *xl, *xu, *gf, *gn, *s, *xo, *go, *xs, *gs, *xm, *gm, *u1, *u2;
576      double gmax, minf_est;
577      double xmax = 0; /* no maximum */
578      double tolg = 0; /* default gradient tolerance */
579      int iest = 0; /* we have no estimate of min function value */
580      int mit = 0, mfg = 0; /* default no limit on #iterations */
581      int mfv = stop->maxeval;
582      stat_common stat;
583      int iterm;
584      int mf;
585
586      ix = (int*) malloc(sizeof(int) * n);
587      if (!ix) return NLOPT_OUT_OF_MEMORY;
588
589      /* FIXME: what should we set mf to?  The example program tlis.for
590         sets it to zero as far as I can tell, but it seems to greatly
591         improve convergence to make it > 0.  The computation time
592         per iteration, and of course the memory, seem to go as O(n * mf),
593         and we'll assume that the main limiting factor is the memory.
594         We'll assume that at least MEMAVAIL memory, or 4*n memory, whichever
595         is bigger, is available. */
596      mf = max(MEMAVAIL/n, 4);
597      if (stop->maxeval && stop->maxeval <= mf)
598           mf = max(stop->maxeval - 5, 1); /* mf > maxeval seems not good */
599
600  retry_alloc:
601      work = (double*) malloc(sizeof(double) * (n * 9 + max(n,n*mf)*2 + 
602                                                max(n,mf)*2));
603      if (!work) { 
604           if (mf > 0) {
605                mf = 0; /* allocate minimal memory */
606                goto retry_alloc;
607           }
608           free(ix);
609           return NLOPT_OUT_OF_MEMORY;
610      }
611
612      xl = work; xu = xl + n;
613      gf = xu + n; gn = gf + n; s = gn + n; 
614      xo = s + n; go = xo + n; xs = go + n; gs = xs + n;
615      xm = gs + n; gm = xm + max(n*mf,n);
616      u1 = gm + max(n*mf,n); u2 = u1 + max(n,mf);
617
618      for (i = 0; i < n; ++i) {
619           int lbu = lb[i] <= -0.99 * HUGE_VAL; /* lb unbounded */
620           int ubu = ub[i] >= 0.99 * HUGE_VAL;  /* ub unbounded */
621           ix[i] = lbu ? (ubu ? 0 : 2) : (ubu ? 1 : (lb[i] == ub[i] ? 5 : 3));
622           xl[i] = lb[i];
623           xu[i] = ub[i];
624      }
625
626      /* ?  xo does not seem to be initialized in the
627         original Fortran code, but it is used upon
628         input to pnet if mf > 0 ... perhaps ALLOCATE initializes
629         arrays to zero by default? */
630      memset(xo, 0, sizeof(double) * max(n,n*mf));
631
632      pnet_(&n, &nb, x, ix, xl, xu, 
633            gf, gn, s, xo, go, xs, gs, xm, gm, u1, u2,
634            &xmax,
635
636            /* fixme: pass tol_rel and tol_abs and use NLopt check */
637            &stop->xtol_rel,
638            &stop->ftol_rel,
639            &stop->minf_max,
640            &tolg,
641            stop,
642
643            &minf_est, &gmax,
644            minf,
645            &mit, &mfv, &mfg,
646            &iest,
647            &mos1, &mos2,
648            &mf,
649            &iterm, &stat,
650            f, f_data);
651
652      free(work);
653      free(ix);
654
655      switch (iterm) {
656          case 1: return NLOPT_XTOL_REACHED;
657          case 2: return NLOPT_FTOL_REACHED;
658          case 3: return NLOPT_MINF_MAX_REACHED;
659          case 4: return NLOPT_SUCCESS; /* gradient tolerance reached */
660          case 6: return NLOPT_SUCCESS;
661          case 12: case 13: return NLOPT_MAXEVAL_REACHED;
662          default: return NLOPT_FAILURE;
663      }
664 }