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