luksan/plip.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 /* *********************************************************************** */
  11 /* SUBROUTINE PLIP               ALL SYSTEMS                   01/09/22 */
  12 /* PURPOSE : */
  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. */
  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 /*  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 */
  49 /*         METHOD. */
  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. */
  62
  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. */
  71
  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 */
  81 /*         UPDATE. */
  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 */
  88 /*         ALF*Y. */
  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 */
  94 /*         TO ZERO. */
  95 /*  S   MXVCOP  COPYING OF A VECTOR. */
  96
  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 */
 107
 108 /* METHOD : */
 109 /* HYBRID METHOD WITH SPARSE MARWIL UPDATES FOR SPARSE LEAST SQUARES */
 110 /* PROBLEMS. */
 111
 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)
 122 {
 123     /* System generated locals */
 124     int i__1;
 125     double d__1, d__2;
 126
 127     /* Builtin functions */
 128
 129     /* Local variables */
 130     int i__, n;
 131     double p, r__;
 132     int kd, ld;
 133     double fo, fp;
 134     int nn;
 135     double po, pp, ro, rp;
 136     int kbf, mec, mfg;
 137     double par;
 138     int mes, kit;
 139     double alf1, alf2, eta0, eta9, par1, par2;
 140     int mes1, mes2, mes3, met3;
 141     double eps8, eps9;
 142     int meta, mred, nred, iold;
 143     double maxf, dmax__;
 144     int inew;
 145     double told;
 146     int ites;
 147     double rmin, rmax, umax, tolp, tols;
 148     int isys;
 149     int ires1, ires2;
 150     int iterd, iterh, mtesf, ntesf;
 151     double gnorm;
 152     int iters, irest, inits, kters, maxst;
 153     double snorm;
 154     int mtesx, ntesx;
 155
 156 /*     INITIATION */
 157
 158     /* Parameter adjustments */
 159     --gr;
 160     --xr;
 161     --xm;
 162     --so;
 163     --go;
 164     --xo;
 165     --s;
 166     --gf;
 167     --xu;
 168     --xl;
 169     --ix;
 170     --x;
 171
 172     /* Function Body */
 173     kbf = 0;
 174     if (*nb > 0) {
 175         kbf = 2;
 176     }
 177     stat_1->nres = 0;
 178     stat_1->ndec = 0;
 179     stat_1->nin = 0;
 180     stat_1->nit = 0;
 181     stop->nevals = 0;
 182     stat_1->nfg = 0;
 183     stat_1->nfh = 0;
 184     isys = 0;
 185     ites = 1;
 186     mtesx = 2;
 187     mtesf = 2;
 188     inits = 2;
 189     *iterm = 0;
 190     iterd = 0;
 191     iters = 2;
 192     iterh = 0;
 193     kters = 3;
 194     irest = 0;
 195     ires1 = 999;
 196     ires2 = 0;
 197     mred = 10;
 198     meta = 1;
 199     met3 = 4;
 200     mec = 4;
 201     mes = 4;
 202     mes1 = 2;
 203     mes2 = 2;
 204     mes3 = 2;
 205     eta0 = 1e-15;
 206     eta9 = 1e120;
 207     eps8 = 1.;
 208     eps9 = 1e-8;
 209     alf1 = 1e-10;
 210     alf2 = 1e10;
 211     rmax = eta9;
 212     dmax__ = eta9;
 213     maxf = 1e20;
 214     if (*iest <= 0) {
 215          *minf_est = -HUGE_VAL; /* changed from -1e60 by SGJ */
 216     }
 217     if (*iest > 0) {
 218         *iest = 1;
 219     }
 220     if (*xmax <= 0.) {
 221         *xmax = 1e16;
 222     }
 223     if (*tolx <= 0.) {
 224         *tolx = 1e-16;
 225     }
 226     if (*tolf <= 0.) {
 227         *tolf = 1e-14;
 228     }
 229     if (*tolg <= 0.) {
 230          *tolg = 1e-8; /* SGJ: was 1e-6, but this sometimes stops too soon */
 231     }
 232 #if 0
 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 */
 236     if (*tolb <= 0.) {
 237         *tolb = *minf_est + 1e-16;
 238     }
 239 #endif
 240     told = 1e-4;
 241     tols = 1e-4;
 242     tolp = .9;
 243     if (*met <= 0) {
 244         *met = 2;
 245     }
 246     /* changed by SGJ: default is no limit (INT_MAX) on # iterations/fevals */
 247     if (*mit <= 0) {
 248         *mit = INT_MAX;
 249     }
 250     if (*mfv <= 0) {
 251         *mfv = INT_MAX;
 252     }
 253     mfg = *mfv;
 254     kd = 1;
 255     ld = -1;
 256     kit = -(ires1 * *nf + ires2);
 257     fo = *minf_est;
 258
 259 /*     INITIAL OPERATIONS WITH SIMPLE BOUNDS */
 260
 261     if (kbf > 0) {
 262         i__1 = *nf;
 263         for (i__ = 1; i__ <= i__1; ++i__) {
 264             if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
 265                 xu[i__] = xl[i__];
 266                 ix[i__] = 5;
 267             } else if (ix[i__] == 5 || ix[i__] == 6) {
 268                 xl[i__] = x[i__];
 269                 xu[i__] = x[i__];
 270                 ix[i__] = 5;
 271             }
 272 /* L2: */
 273         }
 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);
 276     }
 277     if (*iterm != 0) {
 278         goto L11190;
 279     }
 280     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
 281     *f = objgrad(*nf, &x[1], &gf[1], objgrad_data);
 282     ++stop->nevals;
 283     ++stat_1->nfg;
 284 L11120:
 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, &
 289             iters, iterm);
 290     if (*iterm != 0) {
 291         goto L11190;
 292     }
 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, &
 296                 iold, &irest);
 297     }
 298 L11130:
 299     if (irest > 0) {
 300         nn = 0;
 301         par = 1.;
 302         ld = min(ld,1);
 303         if (kit < stat_1->nit) {
 304             ++stat_1->nres;
 305             kit = stat_1->nit;
 306         } else {
 307             *iterm = -10;
 308             if (iters < 0) {
 309                 *iterm = iters - 5;
 310             }
 311         }
 312     }
 313     if (*iterm != 0) {
 314         goto L11190;
 315     }
 316     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
 317
 318 /*     DIRECTION DETERMINATION */
 319
 320     gnorm = sqrt(luksan_mxudot__(nf, &gf[1], &gf[1], &ix[1], &kbf));
 321
 322 /*     NEWTON LIKE STEP */
 323
 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);
 328     iterd = 1;
 329     snorm = sqrt(luksan_mxudot__(nf, &s[1], &s[1], &ix[1], &kbf));
 330
 331 /*     TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH */
 332
 333     if (kd > 0) {
 334         p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
 335     }
 336     if (iterd < 0) {
 337         *iterm = iterd;
 338     } else {
 339
 340 /*     TEST ON DESCENT DIRECTION */
 341
 342         if (snorm <= 0.) {
 343             irest = max(irest,1);
 344         } else if (p + told * gnorm * snorm <= 0.) {
 345             irest = 0;
 346         } else {
 347
 348 /*     UNIFORM DESCENT CRITERION */
 349
 350             irest = max(irest,1);
 351         }
 352         if (irest == 0) {
 353
 354 /*     PREPARATION OF LINE SEARCH */
 355
 356             nred = 0;
 357             rmin = alf1 * gnorm / snorm;
 358 /* Computing MIN */
 359             d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
 360             rmax = min(d__1,d__2);
 361         }
 362     }
 363     if (*iterm != 0) {
 364         goto L11190;
 365     }
 366     if (nlopt_stop_time(stop)) { *iterm = 100; goto L11190; }
 367     if (irest != 0) {
 368         goto L11130;
 369     }
 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);
 372     if (rmax == 0.) {
 373         goto L11175;
 374     }
 375 L11170:
 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);
 379     if (isys == 0) {
 380         goto L11174;
 381     }
 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);
 385     ++stop->nevals;
 386     ++stat_1->nfg;
 387     p = luksan_mxudot__(nf, &gf[1], &s[1], &ix[1], &kbf);
 388     goto L11170;
 389 L11174:
 390     if (iters <= 0) {
 391         r__ = 0.;
 392         *f = fo;
 393         p = po;
 394         luksan_mxvcop__(nf, &xo[1], &x[1]);
 395         luksan_mxvcop__(nf, &go[1], &gf[1]);
 396         irest = max(irest,1);
 397         ld = kd;
 398         goto L11130;
 399     }
 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);
 404     if (nn < *mf) {
 405         luksan_pulsp3__(nf, &nn, mf, &xm[1], &gr[1], &xo[1], &go[1], &r__, &
 406                 po, &par, &iterh, &met3);
 407     } else {
 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);
 410     }
 411 L11175:
 412     if (iterh != 0) {
 413         irest = max(irest,1);
 414     }
 415     if (kbf > 0) {
 416         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
 417     }
 418     goto L11120;
 419 L11190:
 420     return;
 421 } /* plip_ */
 422