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
pnet_(int * nf,int * nb,double * x,int * ix,double * xl,double * xu,double * gf,double * gn,double * s,double * xo,double * go,double * xs,double * gs,double * xm,double * gm,double * u1,double * u2,double * xmax,double * tolx,double * tolf,double * tolb,double * tolg,nlopt_stopping * stop,double * minf_est,double * gmax,double * f,int * mit,int * mfv,int * mfg,int * iest,int * mos1,int * mos2,int * mf,int * iterm,stat_common * stat_1,nlopt_func objgrad,void * objgrad_data)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 = 0.0, b = 0.0;
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 double rho1, rho2, eps8, eps9;
160 int mred, iold, nred;
161 double maxf, dmax__;
162 int xstop = 0;
163 int inew;
164 double told;
165 int ites;
166 double rmin, rmax, umax, tolp, tols;
167 int isys;
168 int ires1, ires2;
169 int iterd, mtesf, ntesf;
170 double gnorm;
171 int iters, irest, inits, kters, maxst;
172 double snorm;
173 int mtesx, ntesx;
174 ps1l01_state state;
175
176 (void) tolb;
177
178 /* INITIATION */
179
180 /* Parameter adjustments */
181 --u2;
182 --u1;
183 --gm;
184 --xm;
185 --gs;
186 --xs;
187 --go;
188 --xo;
189 --s;
190 --gn;
191 --gf;
192 --xu;
193 --xl;
194 --ix;
195 --x;
196
197 /* Function Body */
198 kbf = 0;
199 if (*nb > 0) {
200 kbf = 2;
201 }
202 stat_1->nres = 0;
203 stat_1->ndec = 0;
204 stat_1->nin = 0;
205 stat_1->nit = 0;
206 stat_1->nfg = 0;
207 stat_1->nfh = 0;
208 isys = 0;
209 ites = 1;
210 mtesx = 2;
211 mtesf = 2;
212 inits = 2;
213 *iterm = 0;
214 iterd = 0;
215 iters = 2;
216 kters = 3;
217 irest = 0;
218 ires1 = 999;
219 ires2 = 0;
220 mred = 10;
221 mes = 4;
222 eps = .8;
223 eta0 = 1e-15;
224 eta9 = 1e120;
225 eps8 = 1.;
226 eps9 = 1e-8;
227 alf1 = 1e-10;
228 alf2 = 1e10;
229 rmax = eta9;
230 dmax__ = eta9;
231 maxf = 1e20;
232 if (*iest <= 0) {
233 *minf_est = -HUGE_VAL; /* changed from -1e60 by SGJ */
234 }
235 if (*iest > 0) {
236 *iest = 1;
237 }
238 if (*xmax <= 0.) {
239 *xmax = 1e16;
240 }
241 if (*tolx <= 0.) {
242 *tolx = 1e-16;
243 }
244 if (*tolf <= 0.) {
245 *tolf = 1e-14;
246 }
247 if (*tolg <= 0.) {
248 *tolg = 1e-8; /* SGJ: was 1e-6, but this sometimes stops too soon */
249 }
250 #if 0
251 /* removed by SGJ: this check prevented us from using minf_max <= 0,
252 which doesn't make sense. Instead, if you don't want to have a
253 lower limit, you should set minf_max = -HUGE_VAL */
254 if (*tolb <= 0.) {
255 *tolb = *minf_est + 1e-16;
256 }
257 #endif
258 told = 1e-4;
259 tols = 1e-4;
260 tolp = .9;
261 /* changed by SGJ: default is no limit (INT_MAX) on # iterations/fevals */
262 if (*mit <= 0) {
263 *mit = INT_MAX;
264 }
265 if (*mfv <= 0) {
266 *mfv = INT_MAX;
267 }
268 if (*mfg <= 0) {
269 *mfg = INT_MAX;
270 }
271 if (*mos1 <= 0) {
272 *mos1 = 1;
273 }
274 if (*mos2 <= 0) {
275 *mos2 = 1;
276 }
277 kd = 1;
278 ld = -1;
279 kit = -(ires1 * *nf + ires2);
280 fo = *minf_est;
281
282 /* INITIAL OPERATIONS WITH SIMPLE BOUNDS */
283
284 if (kbf > 0) {
285 i__1 = *nf;
286 for (i__ = 1; i__ <= i__1; ++i__) {
287 if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
288 xu[i__] = xl[i__];
289 ix[i__] = 5;
290 } else if (ix[i__] == 5 || ix[i__] == 6) {
291 xl[i__] = x[i__];
292 xu[i__] = x[i__];
293 ix[i__] = 5;
294 }
295 /* L2: */
296 }
297 luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
298 luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
299 }
300 *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
301 ++(*stop->nevals_p);
302 ++stat_1->nfg;
303 if (nlopt_stop_time(stop)) { *iterm = 100; goto L11080; }
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 = MAX2(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 = MIN2(d__1,d__2);
352 if (par > .01) {
353 /* Computing MIN */
354 d__1 = par, d__2 = 1. / (double) stat_1->nit;
355 par = MIN2(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_p);
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 = MAX2(irest,1);
484 } else if (p + told * gnorm * snorm <= 0.) {
485 irest = 0;
486 } else {
487
488 /* UNIFORM DESCENT CRITERION */
489
490 irest = MAX2(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 = MIN2(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,
520 &isys, &state);
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_p);
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. */
luksan_pnet(int n,nlopt_func f,void * f_data,const double * lb,const double * ub,double * x,double * minf,nlopt_stopping * stop,int mf,int mos1,int mos2)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 if (mf <= 0) {
591 mf = MAX2(MEMAVAIL/n, 10);
592 if (stop->maxeval && stop->maxeval <= mf)
593 mf = MAX2(stop->maxeval, 1);
594 }
595
596 retry_alloc:
597 work = (double*) malloc(sizeof(double) * (n * 9 + MAX2(n,n*mf)*2 +
598 MAX2(n,mf)*2));
599 if (!work) {
600 if (mf > 0) {
601 mf = 0; /* allocate minimal memory */
602 goto retry_alloc;
603 }
604 free(ix);
605 return NLOPT_OUT_OF_MEMORY;
606 }
607
608 xl = work; xu = xl + n;
609 gf = xu + n; gn = gf + n; s = gn + n;
610 xo = s + n; go = xo + n; xs = go + n; gs = xs + n;
611 xm = gs + n; gm = xm + MAX2(n*mf,n);
612 u1 = gm + MAX2(n*mf,n); u2 = u1 + MAX2(n,mf);
613
614 for (i = 0; i < n; ++i) {
615 int lbu = lb[i] <= -0.99 * HUGE_VAL; /* lb unbounded */
616 int ubu = ub[i] >= 0.99 * HUGE_VAL; /* ub unbounded */
617 ix[i] = lbu ? (ubu ? 0 : 2) : (ubu ? 1 : (lb[i] == ub[i] ? 5 : 3));
618 xl[i] = lb[i];
619 xu[i] = ub[i];
620 }
621
622 /* ? xo does not seem to be initialized in the
623 original Fortran code, but it is used upon
624 input to pnet if mf > 0 ... perhaps ALLOCATE initializes
625 arrays to zero by default? */
626 memset(xo, 0, sizeof(double) * MAX2(n,n*mf));
627
628 pnet_(&n, &nb, x, ix, xl, xu,
629 gf, gn, s, xo, go, xs, gs, xm, gm, u1, u2,
630 &xmax,
631
632 /* fixme: pass tol_rel and tol_abs and use NLopt check */
633 &stop->xtol_rel,
634 &stop->ftol_rel,
635 &stop->minf_max,
636 &tolg,
637 stop,
638
639 &minf_est, &gmax,
640 minf,
641 &mit, &mfv, &mfg,
642 &iest,
643 &mos1, &mos2,
644 &mf,
645 &iterm, &stat,
646 f, f_data);
647
648 free(work);
649 free(ix);
650
651 switch (iterm) {
652 case 1: return NLOPT_XTOL_REACHED;
653 case 2: return NLOPT_FTOL_REACHED;
654 case 3: return NLOPT_MINF_MAX_REACHED;
655 case 4: return NLOPT_SUCCESS; /* gradient tolerance reached */
656 case 6: return NLOPT_SUCCESS;
657 case 12: case 13: return NLOPT_MAXEVAL_REACHED;
658 case 100: return NLOPT_MAXTIME_REACHED;
659 case -999: return NLOPT_FORCED_STOP;
660 default: return NLOPT_FAILURE;
661 }
662 }
663