chiark / gitweb /
fix support for maxtime in Luksan algorithms and ORIG_DIRECT; thanks to Jurgen Werner...
[nlopt.git] / 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 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
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, eta0, eta9, par1, par2;
131     int mes1, mes2, mes3;
132     double eps8, eps9;
133     int mred, iold, nred;
134     double maxf, dmax__;
135     int xstop = 0;
136     int inew;
137     double told;
138     int ites;
139     double rmin, rmax, umax, tolp, tols;
140     int isys;
141     int ires1, ires2;
142     int iterd, mtesf, ntesf;
143     double gnorm;
144     int iters, irest, inits, kters, maxst;
145     double snorm;
146     int mtesx, ntesx;
147
148 /*     INITIATION */
149
150     /* Parameter adjustments */
151     --vo;
152     --uo;
153     --go;
154     --xo;
155     --s;
156     --gf;
157     --xu;
158     --xl;
159     --ix;
160     --x;
161
162     /* Function Body */
163     kbf = 0;
164     if (*nb > 0) {
165         kbf = 2;
166     }
167     stat_1->nres = 0;
168     stat_1->ndec = 0;
169     stat_1->nin = 0;
170     stat_1->nit = 0;
171     stat_1->nfg = 0;
172     stat_1->nfh = 0;
173     isys = 0;
174     ites = 1;
175     mtesx = 2;
176     mtesf = 2;
177     inits = 2;
178     *iterm = 0;
179     iterd = 0;
180     iters = 2;
181     kters = 3;
182     irest = 0;
183     ires1 = 999;
184     ires2 = 0;
185     mred = 10;
186     mes = 4;
187     mes1 = 2;
188     mes2 = 2;
189     mes3 = 2;
190     eta0 = 1e-15;
191     eta9 = 1e120;
192     eps8 = 1.;
193     eps9 = 1e-8;
194     alf1 = 1e-10;
195     alf2 = 1e10;
196     rmax = eta9;
197     dmax__ = eta9;
198     maxf = 1e20;
199     if (*iest <= 0) {
200          *minf_est = -HUGE_VAL; /* changed from -1e60 by SGJ */
201     }
202     if (*iest > 0) {
203         *iest = 1;
204     }
205     if (*xmax <= 0.) {
206         *xmax = 1e16;
207     }
208     if (*tolx <= 0.) {
209         *tolx = 1e-16;
210     }
211     if (*tolf <= 0.) {
212         *tolf = 1e-14;
213     }
214     if (*tolg <= 0.) {
215          *tolg = 1e-8; /* SGJ: was 1e-6, but this sometimes stops too soon */
216     }
217 #if 0
218     /* removed by SGJ: this check prevented us from using minf_max <= 0,
219        which doesn't make sense.  Instead, if you don't want to have a
220        lower limit, you should set minf_max = -HUGE_VAL */
221     if (*tolb <= 0.) {
222         *tolb = *minf_est + 1e-16;
223     }
224 #endif
225     told = 1e-4;
226     tols = 1e-4;
227     tolp = .8;
228     /* changed by SGJ: default is no limit (INT_MAX) on # iterations/fevals */
229     if (*mit <= 0) {
230         *mit = INT_MAX;
231     }
232     if (*mfv <= 0) {
233         *mfv = INT_MAX;
234     }
235     mfg = *mfv;
236     kd = 1;
237     ld = -1;
238     kit = -(ires1 * *nf + ires2);
239     fo = *minf_est;
240
241 /*     INITIAL OPERATIONS WITH SIMPLE BOUNDS */
242
243     if (kbf > 0) {
244         i__1 = *nf;
245         for (i__ = 1; i__ <= i__1; ++i__) {
246             if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
247                 xu[i__] = xl[i__];
248                 ix[i__] = 5;
249             } else if (ix[i__] == 5 || ix[i__] == 6) {
250                 xl[i__] = x[i__];
251                 xu[i__] = x[i__];
252                 ix[i__] = 5;
253             }
254 /* L2: */
255         }
256         luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
257         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
258     }
259     if (*iterm != 0) {
260         goto L11190;
261     }
262     *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
263     ++stop->nevals;
264     ++stat_1->nfg;
265     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
266 L11120:
267     luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
268     luksan_pyfut1__(nf, f, &fo, &umax, gmax, xstop, stop, tolg, 
269             &kd, &stat_1->nit, &kit, mit, &stat_1->nfg, &mfg, 
270             &ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
271             iters, iterm);
272     if (*iterm != 0) {
273         goto L11190;
274     }
275     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
276     if (kbf > 0 && rmax > 0.) {
277         luksan_pyrmc0__(nf, &n, &ix[1], &gf[1], &eps8, &umax, gmax, &rmax, &
278                 iold, &irest);
279     }
280 L11130:
281
282 /*     DIRECTION DETERMINATION */
283
284     gnorm = sqrt(luksan_mxudot__(nf, &gf[1], &gf[1], &ix[1], &kbf));
285     if (irest != 0) {
286         goto L12620;
287     }
288 /* Computing MIN */
289     i__1 = stat_1->nit - kit;
290     k = MIN2(i__1,*mf);
291     if (k <= 0) {
292         irest = MAX2(irest,1);
293         goto L12620;
294     }
295
296 /*     DETERMINATION OF THE PARAMETER B */
297
298     b = luksan_mxudot__(nf, &xo[1], &go[1], &ix[1], &kbf);
299     if (b <= 0.) {
300         irest = MAX2(irest,1);
301         goto L12620;
302     }
303     uo[1] = 1. / b;
304     luksan_mxuneg__(nf, &gf[1], &s[1], &ix[1], &kbf);
305     luksan_mxdrcb__(nf, &k, &xo[1], &go[1], &uo[1], &vo[1], &s[1], &ix[1], &
306             kbf);
307     a = luksan_mxudot__(nf, &go[1], &go[1], &ix[1], &kbf);
308     if (a > 0.) {
309         d__1 = b / a;
310         luksan_mxvscl__(nf, &d__1, &s[1], &s[1]);
311     }
312     luksan_mxdrcf__(nf, &k, &xo[1], &go[1], &uo[1], &vo[1], &s[1], &ix[1], &
313             kbf);
314     snorm = sqrt(luksan_mxudot__(nf, &s[1], &s[1], &ix[1], &kbf));
315 /* Computing MIN */
316     i__1 = k + 1;
317     k = MIN2(i__1,*mf);
318     luksan_mxdrsu__(nf, &k, &xo[1], &go[1], &uo[1]);
319 L12620:
320     iterd = 0;
321     if (irest != 0) {
322
323 /*     STEEPEST DESCENT DIRECTION */
324
325         luksan_mxuneg__(nf, &gf[1], &s[1], &ix[1], &kbf);
326         snorm = gnorm;
327         if (kit < stat_1->nit) {
328             ++stat_1->nres;
329             kit = stat_1->nit;
330         } else {
331              *iterm = -10;
332             if (iters < 0) {
333                 *iterm = iters - 5;
334             }
335         }
336     }
337
338 /*     TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH */
339
340     if (kd > 0) {
341         p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
342     }
343     if (iterd < 0) {
344         *iterm = iterd;
345     } else {
346
347 /*     TEST ON DESCENT DIRECTION */
348
349         if (snorm <= 0.) {
350             irest = MAX2(irest,1);
351         } else if (p + told * gnorm * snorm <= 0.) {
352             irest = 0;
353         } else {
354
355 /*     UNIFORM DESCENT CRITERION */
356
357             irest = MAX2(irest,1);
358         }
359         if (irest == 0) {
360
361 /*     PREPARATION OF LINE SEARCH */
362
363             nred = 0;
364             rmin = alf1 * gnorm / snorm;
365 /* Computing MIN */
366             d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
367             rmax = MIN2(d__1,d__2);
368         }
369     }
370     if (*iterm != 0) {
371         goto L11190;
372     }
373     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
374     if (irest != 0) {
375         goto L11130;
376     }
377     luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
378              &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
379     if (rmax == 0.) {
380         goto L11175;
381     }
382 L11170:
383     luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, minf_est, &maxf, &rmin, 
384             &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1->nit, &kit, &
385             nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
386     if (isys == 0) {
387         goto L11174;
388     }
389     luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
390     luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
391     *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
392     ++stop->nevals;
393     ++stat_1->nfg;
394     p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
395     goto L11170;
396 L11174:
397     if (iters <= 0) {
398         r__ = 0.;
399         *f = fo;
400         p = po;
401         luksan_mxvcop__(nf, &xo[1], &x[1]);
402         luksan_mxvcop__(nf, &go[1], &gf[1]);
403         irest = MAX2(irest,1);
404         ld = kd;
405         goto L11130;
406     }
407     luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
408             p, &po, &dmax__, &kbf, &kd, &ld, &iters);
409     xstop = nlopt_stop_dx(stop, &x[1], &xo[1]);
410 L11175:
411     if (kbf > 0) {
412         luksan_mxvine__(nf, &ix[1]);
413         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
414     }
415     goto L11120;
416 L11190:
417     return;
418 } /* plis_ */
419
420 /* NLopt wrapper around plis_, handling dynamic allocation etc. */
421 nlopt_result luksan_plis(int n, nlopt_func f, void *f_data,
422                   const double *lb, const double *ub, /* bounds */
423                   double *x, /* in: initial guess, out: minimizer */
424                   double *minf,
425                   nlopt_stopping *stop)
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      int mf;
438
439      ix = (int*) malloc(sizeof(int) * n);
440      if (!ix) return NLOPT_OUT_OF_MEMORY;
441
442      /* FIXME: what should we set mf to?  The example program tlis.for
443         sets it to zero as far as I can tell, but it seems to greatly
444         improve convergence to make it > 0.  The computation time
445         per iteration, and of course the memory, seem to go as O(n * mf),
446         and we'll assume that the main limiting factor is the memory.
447         We'll assume that at least MEMAVAIL memory, or 4*n memory, whichever
448         is bigger, is available. */
449      mf = MAX2(MEMAVAIL/n, 4);
450      if (stop->maxeval && stop->maxeval <= mf)
451           mf = MAX2(stop->maxeval - 5, 1); /* mf > maxeval seems not good */
452
453  retry_alloc:
454      work = (double*) malloc(sizeof(double) * (n * 4 + MAX2(n,n*mf)*2 + 
455                                                MAX2(n,mf)*2));
456      if (!work) { 
457           if (mf > 0) {
458                mf = 0; /* allocate minimal memory */
459                goto retry_alloc;
460           }
461           free(ix);
462           return NLOPT_OUT_OF_MEMORY;
463      }
464
465      xl = work; xu = xl + n; gf = xu + n; s = gf + n; 
466      xo = s + n; go = xo + MAX2(n,n*mf);
467      uo = go + MAX2(n,n*mf); vo = uo + MAX2(n,mf);
468
469      for (i = 0; i < n; ++i) {
470           int lbu = lb[i] <= -0.99 * HUGE_VAL; /* lb unbounded */
471           int ubu = ub[i] >= 0.99 * HUGE_VAL;  /* ub unbounded */
472           ix[i] = lbu ? (ubu ? 0 : 2) : (ubu ? 1 : (lb[i] == ub[i] ? 5 : 3));
473           xl[i] = lb[i];
474           xu[i] = ub[i];
475      }
476
477      /* ?  xo does not seem to be initialized in the
478         original Fortran code, but it is used upon
479         input to plis if mf > 0 ... perhaps ALLOCATE initializes
480         arrays to zero by default? */
481      memset(xo, 0, sizeof(double) * MAX2(n,n*mf));
482
483      plis_(&n, &nb, x, ix, xl, xu, 
484            gf, s, xo, go, uo, vo,
485            &xmax,
486
487            /* fixme: pass tol_rel and tol_abs and use NLopt check */
488            &stop->xtol_rel,
489            &stop->ftol_rel,
490            &stop->minf_max,
491            &tolg,
492            stop,
493
494            &minf_est, &gmax,
495            minf,
496            &mit, &mfv,
497            &iest,
498            &mf,
499            &iterm, &stat,
500            f, f_data);
501
502      free(work);
503      free(ix);
504
505      switch (iterm) {
506          case 1: return NLOPT_XTOL_REACHED;
507          case 2: return NLOPT_FTOL_REACHED;
508          case 3: return NLOPT_MINF_MAX_REACHED;
509          case 4: return NLOPT_SUCCESS; /* gradient tolerance reached */
510          case 6: return NLOPT_SUCCESS;
511          case 12: case 13: return NLOPT_MAXEVAL_REACHED;
512          case 100: return NLOPT_MAXTIME_REACHED;
513          case -999: return NLOPT_FORCED_STOP;
514          default: return NLOPT_FAILURE;
515      }
516 }