luksan/plip.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3 #include "luksan.h"
   4
   5 #define TRUE_ 1
   6 #define FALSE_ 0
   7
   8 /* Common Block Declarations */
   9
  10 struct {
  11     int nres, ndec, nin, nit, nfv, nfg, nfh;
  12 } stat_;
  13
  14 #define stat_1 stat_
  15
  16 /* *********************************************************************** */
  17 /* SUBROUTINE PLIP               ALL SYSTEMS                   01/09/22 */
  18 /* PURPOSE : */
  19 /* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT */
  20 /* USE THE SHIFTED LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE */
  21 /* PRODUCT FORM UPDATES. */
  22
  23 /* PARAMETERS : */
  24 /*  II  NF  NUMBER OF VARIABLES. */
  25 /*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. */
  26 /*         NB>0-SIMPLE BOUNDS ACCEPTED. */
  27 /*  RI  X(NF)  VECTOR OF VARIABLES. */
  28 /*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE */
  29 /*         X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I). */
  30 /*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND */
  31 /*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. */
  32 /*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. */
  33 /*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. */
  34 /*  RA  GF(NF)  GRADIENT OF THE OBJECTIVE FUNCTION. */
  35 /*  RO  S(NF)  DIRECTION VECTOR. */
  36 /*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE. */
  37 /*  RI  GO(NF)  GRADIENTS DIFFERENCE. */
  38 /*  RA  SO(NF)  AUXILIARY VECTOR. */
  39 /*  RA  XM(NF*MF)  AUXILIARY VECTOR. */
  40 /*  RA  XR(MF)  AUXILIARY VECTOR. */
  41 /*  RA  GR(MF)  AUXILIARY VECTOR. */
  42 /*  RI  XMAX  MAXIMUM STEPSIZE. */
  43 /*  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES. */
  44 /*  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES. */
  45 /*  RI  TOLB  TOLERANCE FOR THE FUNCTION VALUE. */
  46 /*  RI  TOLG  TOLERANCE FOR THE GRADIENT NORM. */
  47 /*  RI  FMIN  ESTIMATION OF THE MINIMUM FUNCTION VALUE. */
  48 /*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE. */
  49 /*  RO  F  VALUE OF THE OBJECTIVE FUNCTION. */
  50 /*  II  MIT  MAXIMUM NUMBER OF ITERATIONS. */
  51 /*  II  MFV  MAXIMUM NUMBER OF FUNCTION EVALUATIONS. */
  52 /*  II  IEST  ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. */
  53 /*         IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN. */
  54 /*  II  MET  METHOD USED. MET=1-RANK-ONE METHOD. MET=2-RANK-TWO */
  55 /*         METHOD. */
  56 /*  II  MF  NUMBER OF LIMITED MEMORY STEPS. */
  57 /*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT. */
  58 /*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS. */
  59 /*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS. */
  60 /*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL */
  61 /*         RESULTS. */
  62 /*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. */
  63 /*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN */
  64 /*                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
  65 /*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN */
  66 /*                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
  67 /*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. */
  68 /*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. */
  69 /*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, */
  70 /*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. */
  71 /*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. */
  72 /*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. */
  73
  74 /* VARIABLES IN COMMON /STAT/ (STATISTICS) : */
  75 /*  IO  NRES  NUMBER OF RESTARTS. */
  76 /*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION. */
  77 /*  IO  NIN  NUMBER OF INNER ITERATIONS. */
  78 /*  IO  NIT  NUMBER OF ITERATIONS. */
  79 /*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS. */
  80 /*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS. */
  81 /*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS. */
  82
  83 /* SUBPROGRAMS USED : */
  84 /*  S   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS. */
  85 /*  S   PS1L01  STEPSIZE SELECTION USING LINE SEARCH. */
  86 /*  S   PULSP3  SHIFTED VARIABLE METRIC UPDATE. */
  87 /*  S   PULVP3  SHIFTED LIMITED-MEMORY VARIABLE METRIC UPDATE. */
  88 /*  S   PYADC0  ADDITION OF A BOX CONSTRAINT. */
  89 /*  S   PYFUT1  TEST ON TERMINATION. */
  90 /*  S   PYRMC0  DELETION OF A BOX CONSTRAINT. */
  91 /*  S   PYTRCD  COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC */
  92 /*         UPDATE. */
  93 /*  S   PYTRCG  COMPUTATION OF THE PROJECTED GRADIENT. */
  94 /*  S   PYTRCS  COMPUTATION OF THE PROJECTED DIRECTION VECTOR. */
  95 /*  S   MXDRMM  MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR */
  96 /*         MATRIX A BY A VECTOR X. */
  97 /*  S   MXDCMD  MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR */
  98 /*         MATRIX A BY A VECTOR X AND ADDITION OF THE SCALED VECTOR */
  99 /*         ALF*Y. */
 100 /*  S   MXUCOP  COPYING OF A VECTOR. */
 101 /*  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR. */
 102 /*  RF  MXUDOT  DOT PRODUCT OF TWO VECTORS. */
 103 /*  S   MXUNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN. */
 104 /*  S   MXUZER  VECTOR ELEMENTS CORRESPONDING TO ACTIVE BOUNDS ARE SET */
 105 /*         TO ZERO. */
 106 /*  S   MXVCOP  COPYING OF A VECTOR. */
 107
 108 /* EXTERNAL SUBROUTINES : */
 109 /*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. */
 110 /*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER */
 111 /*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS */
 112 /*         THE VALUE OF THE OBJECTIVE FUNCTION. */
 113 /*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. */
 114 /*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER */
 115 /*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF) */
 116 /*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION. */
 117
 118 /* METHOD : */
 119 /* HYBRID METHOD WITH SPARSE MARWIL UPDATES FOR SPARSE LEAST SQUARES */
 120 /* PROBLEMS. */
 121
 122 static void plip_(int *nf, int *nb, double *x, int *
 123         ix, double *xl, double *xu, double *gf, double *s,
 124         double *xo, double *go, double *so, double *xm,
 125         double *xr, double *gr, double *xmax, double *tolx,
 126         double *tolf, double *tolb, double *tolg, double *
 127         fmin, double *gmax, double *f, int *mit, int *mfv,
 128         int *iest, int *met, int *mf, int *iprnt, int *
 129         iterm)
 130 {
 131     /* System generated locals */
 132     int i__1;
 133     double d__1, d__2;
 134
 135     /* Builtin functions */
 136
 137     /* Local variables */
 138     int i__, n;
 139     double p, r__;
 140     int kd, ld;
 141     double fo, fp;
 142     int nn;
 143     double po, pp, ro, rp;
 144     int kbf, mec, mfg;
 145     extern static void obj_(int *, double *, double *);
 146     double par;
 147     int mes, kit;
 148     double alf1, alf2, eta0, eta9, par1, par2;
 149     int mes1, mes2, mes3, met3;
 150     double eps8, eps9;
 151     extern static void dobj_(int *, double *, double *);
 152     int meta, mred, nred, iold;
 153     double fmax, dmax__;
 154     extern static void luksan_ps1l01__(double *, double *,
 155             double *, double *, double *, double *,
 156             double *, double *, double *, double *,
 157             double *, double *, double *, double *,
 158             double *, double *, int *, int *, int *,
 159             int *, int *, int *, int *, int *, int *,
 160             int *, int *, int *, int *);
 161     int inew;
 162     double told;
 163     int ites;
 164     double rmin, rmax, umax, tolp, tols;
 165     int isys;
 166     extern static void luksan_pcbs04__(int *, double *,
 167             int *, double *, double *, double *, int *);
 168     int ires1, ires2;
 169     extern static void luksan_pyadc0__(int *, int *,
 170             double *, int *, double *, double *, int *);
 171     int iterd, iterh, mtesf, ntesf;
 172     double gnorm;
 173     extern static void luksan_pyrmc0__(int *, int *, int
 174             *, double *, double *, double *, double *,
 175             double *, int *, int *);
 176     int iters, irest, inits, kters, maxst;
 177     double snorm;
 178     int mtesx, ntesx;
 179     extern static void luksan_pulsp3__(int *, int *, int
 180             *, double *, double *, double *, double *,
 181             double *, double *, double *, int *, int *),
 182             luksan_pyfut1__(int *, double *, double *, double
 183             *, double *, double *, double *, double *,
 184             double *, double *, int *, int *, int *,
 185             int *, int *, int *, int *, int *, int *,
 186             int *, int *, int *, int *, int *, int *,
 187             int *, int *, int *), luksan_pulvp3__(int *,
 188             int *, double *, double *, double *, double *,
 189              double *, double *, double *, double *,
 190             double *, double *, int *, int *, int *,
 191             int *), luksan_mxdcmd__(int *, int *, double *,
 192             double *, double *, double *, double *),
 193             luksan_mxuneg__(int *, double *, double *, int *,
 194             int *), luksan_mxdrmm__(int *, int *, double *,
 195             double *, double *), luksan_pytrcd__(int *,
 196             double *, int *, double *, double *, double *,
 197              double *, double *, double *, double *,
 198             double *, double *, int *, int *, int *,
 199             int *), luksan_pytrcg__(int *, int *, int *,
 200             double *, double *, double *, int *, int *),
 201             luksan_mxudir__(int *, double *, double *, double
 202             *, double *, int *, int *), luksan_mxucop__(int *,
 203              double *, double *, int *, int *),
 204             luksan_mxvcop__(int *, double *, double *),
 205             luksan_pytrcs__(int *, double *, int *, double *,
 206             double *, double *, double *, double *,
 207             double *, double *, double *, double *,
 208             double *, double *, double *, double *,
 209             double *, int *), luksan_mxuzer__(int *, double *,
 210              int *, int *);
 211     extern double mxudot_(int *, double *, double *, int *
 212             , int *);
 213
 214
 215 /*     INITIATION */
 216
 217     /* Parameter adjustments */
 218     --gr;
 219     --xr;
 220     --xm;
 221     --so;
 222     --go;
 223     --xo;
 224     --s;
 225     --gf;
 226     --xu;
 227     --xl;
 228     --ix;
 229     --x;
 230
 231     /* Function Body */
 232     kbf = 0;
 233     if (*nb > 0) {
 234         kbf = 2;
 235     }
 236     stat_1.nres = 0;
 237     stat_1.ndec = 0;
 238     stat_1.nin = 0;
 239     stat_1.nit = 0;
 240     stat_1.nfv = 0;
 241     stat_1.nfg = 0;
 242     stat_1.nfh = 0;
 243     isys = 0;
 244     ites = 1;
 245     mtesx = 2;
 246     mtesf = 2;
 247     inits = 2;
 248     *iterm = 0;
 249     iterd = 0;
 250     iters = 2;
 251     iterh = 0;
 252     kters = 3;
 253     irest = 0;
 254     ires1 = 999;
 255     ires2 = 0;
 256     mred = 10;
 257     meta = 1;
 258     met3 = 4;
 259     mec = 4;
 260     mes = 4;
 261     mes1 = 2;
 262     mes2 = 2;
 263     mes3 = 2;
 264     eta0 = 1e-15;
 265     eta9 = 1e120;
 266     eps8 = 1.;
 267     eps9 = 1e-8;
 268     alf1 = 1e-10;
 269     alf2 = 1e10;
 270     rmax = eta9;
 271     dmax__ = eta9;
 272     fmax = 1e20;
 273     if (*iest <= 0) {
 274         *fmin = -1e60;
 275     }
 276     if (*iest > 0) {
 277         *iest = 1;
 278     }
 279     if (*xmax <= 0.) {
 280         *xmax = 1e16;
 281     }
 282     if (*tolx <= 0.) {
 283         *tolx = 1e-16;
 284     }
 285     if (*tolf <= 0.) {
 286         *tolf = 1e-14;
 287     }
 288     if (*tolg <= 0.) {
 289         *tolg = 1e-6;
 290     }
 291     if (*tolb <= 0.) {
 292         *tolb = *fmin + 1e-16;
 293     }
 294     told = 1e-4;
 295     tols = 1e-4;
 296     tolp = .9;
 297     if (*met <= 0) {
 298         *met = 2;
 299     }
 300     if (*mit <= 0) {
 301         *mit = 9000;
 302     }
 303     if (*mfv <= 0) {
 304         *mfv = 9000;
 305     }
 306     mfg = *mfv;
 307     kd = 1;
 308     ld = -1;
 309     kit = -(ires1 * *nf + ires2);
 310     fo = *fmin;
 311
 312 /*     INITIAL OPERATIONS WITH SIMPLE BOUNDS */
 313
 314     if (kbf > 0) {
 315         i__1 = *nf;
 316         for (i__ = 1; i__ <= i__1; ++i__) {
 317             if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
 318                 xu[i__] = xl[i__];
 319                 ix[i__] = 5;
 320             } else if (ix[i__] == 5 || ix[i__] == 6) {
 321                 xl[i__] = x[i__];
 322                 xu[i__] = x[i__];
 323                 ix[i__] = 5;
 324             }
 325 /* L2: */
 326         }
 327         luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
 328         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
 329     }
 330     if (*iterm != 0) {
 331         goto L11190;
 332     }
 333     obj_(nf, &x[1], f);
 334     ++stat_1.nfv;
 335     dobj_(nf, &x[1], &gf[1]);
 336     ++stat_1.nfg;
 337 L11120:
 338     luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
 339     luksan_pyfut1__(nf, f, &fo, &umax, gmax, &dmax__, tolx, tolf, tolb, tolg,
 340             &kd, &stat_1.nit, &kit, mit, &stat_1.nfv, mfv, &stat_1.nfg, &mfg,
 341             &ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
 342             iters, iterm);
 343     if (*iterm != 0) {
 344         goto L11190;
 345     }
 346     if (kbf > 0 && rmax > 0.) {
 347         luksan_pyrmc0__(nf, &n, &ix[1], &gf[1], &eps8, &umax, gmax, &rmax, &
 348                 iold, &irest);
 349     }
 350 L11130:
 351     if (irest > 0) {
 352         nn = 0;
 353         par = 1.;
 354         ld = min(ld,1);
 355         if (kit < stat_1.nit) {
 356             ++stat_1.nres;
 357             kit = stat_1.nit;
 358         } else {
 359             *iterm = -10;
 360             if (iters < 0) {
 361                 *iterm = iters - 5;
 362             }
 363         }
 364     }
 365     if (*iterm != 0) {
 366         goto L11190;
 367     }
 368
 369 /*     DIRECTION DETERMINATION */
 370
 371     gnorm = sqrt(mxudot_(nf, &gf[1], &gf[1], &ix[1], &kbf));
 372
 373 /*     NEWTON LIKE STEP */
 374
 375     luksan_mxuneg__(nf, &gf[1], &s[1], &ix[1], &kbf);
 376     luksan_mxdrmm__(nf, &nn, &xm[1], &s[1], &gr[1]);
 377     luksan_mxdcmd__(nf, &nn, &xm[1], &gr[1], &par, &s[1], &s[1]);
 378     luksan_mxuzer__(nf, &s[1], &ix[1], &kbf);
 379     iterd = 1;
 380     snorm = sqrt(mxudot_(nf, &s[1], &s[1], &ix[1], &kbf));
 381
 382 /*     TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH */
 383
 384     if (kd > 0) {
 385         p = mxudot_(nf, &gf[1], &s[1], &ix[1], &kbf);
 386     }
 387     if (iterd < 0) {
 388         *iterm = iterd;
 389     } else {
 390
 391 /*     TEST ON DESCENT DIRECTION */
 392
 393         if (snorm <= 0.) {
 394             irest = max(irest,1);
 395         } else if (p + told * gnorm * snorm <= 0.) {
 396             irest = 0;
 397         } else {
 398
 399 /*     UNIFORM DESCENT CRITERION */
 400
 401             irest = max(irest,1);
 402         }
 403         if (irest == 0) {
 404
 405 /*     PREPARATION OF LINE SEARCH */
 406
 407             nred = 0;
 408             rmin = alf1 * gnorm / snorm;
 409 /* Computing MIN */
 410             d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
 411             rmax = min(d__1,d__2);
 412         }
 413     }
 414     if (*iterm != 0) {
 415         goto L11190;
 416     }
 417     if (irest != 0) {
 418         goto L11130;
 419     }
 420     luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
 421              &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
 422     if (rmax == 0.) {
 423         goto L11175;
 424     }
 425 L11170:
 426     luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, fmin, &fmax, &rmin,
 427             &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1.nit, &kit, &
 428             nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
 429     if (isys == 0) {
 430         goto L11174;
 431     }
 432     luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
 433     luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
 434     obj_(nf, &x[1], f);
 435     ++stat_1.nfv;
 436     dobj_(nf, &x[1], &gf[1]);
 437     ++stat_1.nfg;
 438     p = mxudot_(nf, &gf[1], &s[1], &ix[1], &kbf);
 439     goto L11170;
 440 L11174:
 441     if (iters <= 0) {
 442         r__ = 0.;
 443         *f = fo;
 444         p = po;
 445         luksan_mxvcop__(nf, &xo[1], &x[1]);
 446         luksan_mxvcop__(nf, &go[1], &gf[1]);
 447         irest = max(irest,1);
 448         ld = kd;
 449         goto L11130;
 450     }
 451     luksan_mxuneg__(nf, &go[1], &s[1], &ix[1], &kbf);
 452     luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
 453             p, &po, &dmax__, &kbf, &kd, &ld, &iters);
 454     luksan_mxucop__(nf, &gf[1], &so[1], &ix[1], &kbf);
 455     if (nn < *mf) {
 456         luksan_pulsp3__(nf, &nn, mf, &xm[1], &gr[1], &xo[1], &go[1], &r__, &
 457                 po, &par, &iterh, &met3);
 458     } else {
 459         luksan_pulvp3__(nf, &nn, &xm[1], &xr[1], &gr[1], &s[1], &so[1], &xo[1]
 460                 , &go[1], &r__, &po, &par, &iterh, &mec, &met3, met);
 461     }
 462 L11175:
 463     if (iterh != 0) {
 464         irest = max(irest,1);
 465     }
 466     if (kbf > 0) {
 467         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
 468     }
 469     goto L11120;
 470 L11190:
 471     return;
 472 } /* plip_ */
 473