7 #define max(a,b) ((a) > (b) ? (a) : (b))
8 #define min(a,b) ((a) < (b) ? (a) : (b))
10 /* *********************************************************************** */
11 /* SUBROUTINE PLIP ALL SYSTEMS 01/09/22 */
13 /* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT */
14 /* USE THE SHIFTED LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE */
15 /* PRODUCT FORM UPDATES. */
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 /* RA 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 SO(NF) AUXILIARY VECTOR. */
33 /* RA XM(NF*MF) AUXILIARY VECTOR. */
34 /* RA XR(MF) AUXILIARY VECTOR. */
35 /* RA GR(MF) AUXILIARY VECTOR. */
36 /* RI XMAX MAXIMUM STEPSIZE. */
37 /* RI TOLX TOLERANCE FOR CHANGE OF VARIABLES. */
38 /* RI TOLF TOLERANCE FOR CHANGE OF FUNCTION VALUES. */
39 /* RI TOLB TOLERANCE FOR THE FUNCTION VALUE. */
40 /* RI TOLG TOLERANCE FOR THE GRADIENT NORM. */
41 /* RI MINF_EST ESTIMATION OF THE MINIMUM FUNCTION VALUE. */
42 /* RO GMAX MAXIMUM PARTIAL DERIVATIVE. */
43 /* RO F VALUE OF THE OBJECTIVE FUNCTION. */
44 /* II MIT MAXIMUM NUMBER OF ITERATIONS. */
45 /* II MFV MAXIMUM NUMBER OF FUNCTION EVALUATIONS. */
46 /* II IEST ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. */
47 /* IEST=1-MINIMUM IS ESTIMATED BY THE VALUE MINF_EST. */
48 /* II MET METHOD USED. MET=1-RANK-ONE METHOD. MET=2-RANK-TWO */
50 /* II MF NUMBER OF LIMITED MEMORY STEPS. */
51 /* IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. */
52 /* ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN */
53 /* MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
54 /* ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN */
55 /* MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
56 /* ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. */
57 /* ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. */
58 /* ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, */
59 /* BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. */
60 /* ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. */
61 /* ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. */
63 /* VARIABLES IN COMMON /STAT/ (STATISTICS) : */
64 /* IO NRES NUMBER OF RESTARTS. */
65 /* IO NDEC NUMBER OF MATRIX DECOMPOSITION. */
66 /* IO NIN NUMBER OF INNER ITERATIONS. */
67 /* IO NIT NUMBER OF ITERATIONS. */
68 /* IO NFV NUMBER OF FUNCTION EVALUATIONS. */
69 /* IO NFG NUMBER OF GRADIENT EVALUATIONS. */
70 /* IO NFH NUMBER OF HESSIAN EVALUATIONS. */
72 /* SUBPROGRAMS USED : */
73 /* S PCBS04 ELIMINATION OF BOX CONSTRAINT VIOLATIONS. */
74 /* S PS1L01 STEPSIZE SELECTION USING LINE SEARCH. */
75 /* S PULSP3 SHIFTED VARIABLE METRIC UPDATE. */
76 /* S PULVP3 SHIFTED LIMITED-MEMORY VARIABLE METRIC UPDATE. */
77 /* S PYADC0 ADDITION OF A BOX CONSTRAINT. */
78 /* S PYFUT1 TEST ON TERMINATION. */
79 /* S PYRMC0 DELETION OF A BOX CONSTRAINT. */
80 /* S PYTRCD COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC */
82 /* S PYTRCG COMPUTATION OF THE PROJECTED GRADIENT. */
83 /* S PYTRCS COMPUTATION OF THE PROJECTED DIRECTION VECTOR. */
84 /* S MXDRMM MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR */
85 /* MATRIX A BY A VECTOR X. */
86 /* S MXDCMD MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR */
87 /* MATRIX A BY A VECTOR X AND ADDITION OF THE SCALED VECTOR */
89 /* S MXUCOP COPYING OF A VECTOR. */
90 /* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. */
91 /* RF MXUDOT DOT PRODUCT OF TWO VECTORS. */
92 /* S MXUNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. */
93 /* S MXUZER VECTOR ELEMENTS CORRESPONDING TO ACTIVE BOUNDS ARE SET */
95 /* S MXVCOP COPYING OF A VECTOR. */
97 /* EXTERNAL SUBROUTINES : */
98 /* SE OBJ COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. */
99 /* CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER */
100 /* OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS */
101 /* THE VALUE OF THE OBJECTIVE FUNCTION. */
102 /* SE DOBJ COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. */
103 /* CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER */
104 /* OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF) */
105 /* IS THE GRADIENT OF THE OBJECTIVE FUNCTION. */
106 /* -- OBJ and DOBJ are replaced by a single function, objgrad, in NLopt */
109 /* HYBRID METHOD WITH SPARSE MARWIL UPDATES FOR SPARSE LEAST SQUARES */
112 static void plip_(int *nf, int *nb, double *x, int *
113 ix, double *xl, double *xu, double *gf, double *s,
114 double *xo, double *go, double *so, double *xm,
115 double *xr, double *gr, double *xmax, double *tolx,
116 double *tolf, double *tolb, double *tolg,
117 nlopt_stopping *stop, double *
118 minf_est, double *gmax, double *f, int *mit, int *mfv,
119 int *iest, int *met, int *mf,
120 int *iterm, stat_common *stat_1,
121 nlopt_func objgrad, void *objgrad_data)
123 /* System generated locals */
127 /* Builtin functions */
129 /* Local variables */
135 double po, pp, ro, rp;
139 double alf1, alf2, eta0, eta9, par1, par2;
140 int mes1, mes2, mes3, met3;
142 int meta, mred, nred, iold;
147 double rmin, rmax, umax, tolp, tols;
150 int iterd, iterh, mtesf, ntesf;
152 int iters, irest, inits, kters, maxst;
158 /* Parameter adjustments */
215 *minf_est = -HUGE_VAL; /* changed from -1e60 by SGJ */
230 *tolg = 1e-8; /* SGJ: was 1e-6, but this sometimes stops too soon */
233 /* removed by SGJ: this check prevented us from using minf_max <= 0,
234 which doesn't make sense. Instead, if you don't want to have a
235 lower limit, you should set minf_max = -HUGE_VAL */
237 *tolb = *minf_est + 1e-16;
246 /* changed by SGJ: default is no limit (INT_MAX) on # iterations/fevals */
256 kit = -(ires1 * *nf + ires2);
259 /* INITIAL OPERATIONS WITH SIMPLE BOUNDS */
263 for (i__ = 1; i__ <= i__1; ++i__) {
264 if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
267 } else if (ix[i__] == 5 || ix[i__] == 6) {
274 luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
275 luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
280 if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
281 *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
285 luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
286 luksan_pyfut1__(nf, f, &fo, &umax, gmax, &dmax__, tolx, tolf, tolb, tolg,
287 &kd, &stat_1->nit, &kit, mit, &stop->nevals, mfv, &stat_1->nfg, &mfg,
288 &ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
293 if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
294 if (kbf > 0 && rmax > 0.) {
295 luksan_pyrmc0__(nf, &n, &ix[1], &gf[1], &eps8, &umax, gmax, &rmax, &
303 if (kit < stat_1->nit) {
316 if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
318 /* DIRECTION DETERMINATION */
320 gnorm = sqrt(luksan_mxudot__(nf, &gf[1], &gf[1], &ix[1], &kbf));
322 /* NEWTON LIKE STEP */
324 luksan_mxuneg__(nf, &gf[1], &s[1], &ix[1], &kbf);
325 luksan_mxdrmm__(nf, &nn, &xm[1], &s[1], &gr[1]);
326 luksan_mxdcmd__(nf, &nn, &xm[1], &gr[1], &par, &s[1], &s[1]);
327 luksan_mxuzer__(nf, &s[1], &ix[1], &kbf);
329 snorm = sqrt(luksan_mxudot__(nf, &s[1], &s[1], &ix[1], &kbf));
331 /* TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH */
334 p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
340 /* TEST ON DESCENT DIRECTION */
343 irest = max(irest,1);
344 } else if (p + told * gnorm * snorm <= 0.) {
348 /* UNIFORM DESCENT CRITERION */
350 irest = max(irest,1);
354 /* PREPARATION OF LINE SEARCH */
357 rmin = alf1 * gnorm / snorm;
359 d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
360 rmax = min(d__1,d__2);
366 if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
370 luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
371 &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
376 luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, minf_est, &maxf, &rmin,
377 &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1->nit, &kit, &
378 nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
382 luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
383 luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
384 *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
387 p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
394 luksan_mxvcop__(nf, &xo[1], &x[1]);
395 luksan_mxvcop__(nf, &go[1], &gf[1]);
396 irest = max(irest,1);
400 luksan_mxuneg__(nf, &go[1], &s[1], &ix[1], &kbf);
401 luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
402 p, &po, &dmax__, &kbf, &kd, &ld, &iters);
403 luksan_mxucop__(nf, &gf[1], &so[1], &ix[1], &kbf);
405 luksan_pulsp3__(nf, &nn, mf, &xm[1], &gr[1], &xo[1], &go[1], &r__, &
406 po, &par, &iterh, &met3);
408 luksan_pulvp3__(nf, &nn, &xm[1], &xr[1], &gr[1], &s[1], &so[1], &xo[1]
409 , &go[1], &r__, &po, &par, &iterh, &mec, &met3, met);
413 irest = max(irest,1);
416 luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
423 /* NLopt wrapper around plip_, handling dynamic allocation etc. */
424 nlopt_result luksan_plip(int n, nlopt_func f, void *f_data,
425 const double *lb, const double *ub, /* bounds */
426 double *x, /* in: initial guess, out: minimizer */
428 nlopt_stopping *stop,
429 int method) /* 1 or 2, see below */
432 double *work, *xl, *xu, *gf, *s, *xo, *go, *so, *xm, *xr, *gr;
433 double gmax, minf_est;
434 double xmax = 0; /* no maximum */
435 double tolg = 0; /* default gradient tolerance */
436 int iest = 0; /* we have no estimate of min function value */
437 int mit = 0; /* default no limit on #iterations */
438 int mfv = stop->maxeval;
443 ix = (int*) malloc(sizeof(int) * n);
444 if (!ix) return NLOPT_OUT_OF_MEMORY;
446 /* FIXME: what should we set mf to? The example program tlis.for
447 sets it to zero as far as I can tell, but it seems to greatly
448 improve convergence to make it > 0. The computation time
449 per iteration, and of course the memory, seem to go as O(n * mf),
450 and we'll assume that the main limiting factor is the memory.
451 We'll assume that at least MEMAVAIL memory, or 4*n memory, whichever
452 is bigger, is available. */
453 mf = max(MEMAVAIL/n, 4);
454 if (stop->maxeval && stop->maxeval <= mf)
455 mf = max(stop->maxeval - 5, 1); /* mf > maxeval seems not good */
458 work = (double*) malloc(sizeof(double) * (n * 7 + max(n,n*mf) +
462 mf = 0; /* allocate minimal memory */
466 return NLOPT_OUT_OF_MEMORY;
469 xl = work; xu = xl + n;
470 gf = xu + n; s = gf + n; xo = s + n; go = xo + n; so = go + n;
472 xr = xm + max(n*mf,n); gr = xr + max(n,mf);
474 for (i = 0; i < n; ++i) {
475 int lbu = lb[i] <= -0.99 * HUGE_VAL; /* lb unbounded */
476 int ubu = ub[i] >= 0.99 * HUGE_VAL; /* ub unbounded */
477 ix[i] = lbu ? (ubu ? 0 : 2) : (ubu ? 1 : (lb[i] == ub[i] ? 5 : 3));
482 /* ? xo does not seem to be initialized in the
483 original Fortran code, but it is used upon
484 input to plip if mf > 0 ... perhaps ALLOCATE initializes
485 arrays to zero by default? */
486 memset(xo, 0, sizeof(double) * max(n,n*mf));
488 plip_(&n, &nb, x, ix, xl, xu,
489 gf, s, xo, go, so, xm, xr, gr,
492 /* fixme: pass tol_rel and tol_abs and use NLopt check */
503 &method, /* 1 == rank-one method VAR1, 2 == rank-two method VAR2 */
512 case 1: return NLOPT_XTOL_REACHED;
513 case 2: return NLOPT_FTOL_REACHED;
514 case 3: return NLOPT_MINF_MAX_REACHED;
515 case 4: return NLOPT_SUCCESS; /* gradient tolerance reached */
516 case 6: return NLOPT_SUCCESS;
517 case 12: case 13: return NLOPT_MAXEVAL_REACHED;
518 default: return NLOPT_FAILURE;