src/algs/luksan/pssubs.c

   1 #include <math.h>
   2 #include "luksan.h"
   3
   4 #define FALSE_ 0
   5 #define MAX2(a,b) ((a) > (b) ? (a) : (b))
   6 #define MIN2(a,b) ((a) < (b) ? (a) : (b))
   7 #define iabs(a) ((a) < 0 ? -(a) : (a))
   8
   9 /*     subroutines extracted from pssubs.for */
  10 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
  11 /* SUBROUTINE PCBS04             ALL SYSTEMS                   98/12/01
  12 * PURPOSE :
  13 * INITIATION OF THE VECTOR CONTAINING TYPES OF CONSTRAINTS.
  14 *
  15 * PARAMETERS :
  16 *  II  NF  NUMBER OF VARIABLES.
  17 *  RI  X(NF)  VECTOR OF VARIABLES.
  18 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
  19 *  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
  20 *  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
  21 *  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
  22 *  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
  23 *         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
  24 */
  25 void luksan_pcbs04__(int *nf, double *x, int *ix,
  26         double *xl, double *xu, double *eps9, int *kbf)
  27 {
  28     /* System generated locals */
  29     int i__1, i__2;
  30     double d__1, d__2;
  31
  32     /* Local variables */
  33     int i__, ixi;
  34     double temp;
  35
  36     /* Parameter adjustments */
  37     --xu;
  38     --xl;
  39     --ix;
  40     --x;
  41
  42     /* Function Body */
  43     if (*kbf > 0) {
  44         i__1 = *nf;
  45         for (i__ = 1; i__ <= i__1; ++i__) {
  46             temp = 1.;
  47             ixi = (i__2 = ix[i__], iabs(i__2));
  48 /* Computing MAX */
  49             d__2 = (d__1 = xl[i__], fabs(d__1));
  50             if ((ixi == 1 || ixi == 3 || ixi == 4) && x[i__] <= xl[i__] + *
  51                     eps9 * MAX2(d__2,temp)) {
  52                 x[i__] = xl[i__];
  53             }
  54 /* Computing MAX */
  55             d__2 = (d__1 = xu[i__], fabs(d__1));
  56             if ((ixi == 2 || ixi == 3 || ixi == 4) && x[i__] >= xu[i__] - *
  57                     eps9 * MAX2(d__2,temp)) {
  58                 x[i__] = xu[i__];
  59             }
  60 /* L1: */
  61         }
  62     }
  63     return;
  64 } /* luksan_pcbs04__ */
  65
  66 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
  67 /* SUBROUTINE PNINT1                ALL SYSTEMS                91/12/01 */
  68 /* PURPOSE : */
  69 /* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITH DIRECTIONAL */
  70 /* DERIVATIVES. */
  71
  72 /* PARAMETERS : */
  73 /*  RI  RL  LOWER VALUE OF THE STEPSIZE PARAMETER. */
  74 /*  RI  RU  UPPER VALUE OF THE STEPSIZE PARAMETER. */
  75 /*  RI  FL  VALUE OF THE OBJECTIVE FUNCTION FOR R=RL. */
  76 /*  RI  FU  VALUE OF THE OBJECTIVE FUNCTION FOR R=RU. */
  77 /*  RI  PL  DIRECTIONAL DERIVATIVE FOR R=RL. */
  78 /*  RI  PU  DIRECTIONAL DERIVATIVE FOR R=RU. */
  79 /*  RO  R  VALUE OF THE STEPSIZE PARAMETER OBTAINED. */
  80 /*  II  MODE  MODE OF LINE SEARCH. */
  81 /*  II  MTYP  METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-QUADRATIC */
  82 /*         INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE). */
  83 /*         MTYP=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL */
  84 /*         DERIVATIVES). MTYP=4-CUBIC INTERPOLATION. MTYP=5-CONIC */
  85 /*         INTERPOLATION. */
  86 /*  IO  MERR  ERROR INDICATOR. MERR=0 FOR NORMAL RETURN. */
  87
  88 /* METHOD : */
  89 /* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS. */
  90
  91 void luksan_pnint1__(double *rl, double *ru, double *fl,
  92         double *fu, double *pl, double *pu, double *r__,
  93         int *mode, int *mtyp, int *merr)
  94 {
  95     /* System generated locals */
  96     double d__1, d__2;
  97
  98     /* Local variables */
  99     double a = 0.0, b = 0.0, c__, d__, den = 0.0, dis;
 100     int ntyp;
 101
 102     *merr = 0;
 103     if (*mode <= 0) {
 104         return;
 105     }
 106     if (*pl >= 0.) {
 107         *merr = 2;
 108         return;
 109     } else if (*ru <= *rl) {
 110         *merr = 3;
 111         return;
 112     }
 113     for (ntyp = *mtyp; ntyp >= 1; --ntyp) {
 114         if (ntyp == 1) {
 115
 116 /*     BISECTION */
 117
 118             if (*mode == 1) {
 119                 *r__ = *ru * 4.;
 120                 return;
 121             } else {
 122                 *r__ = (*rl + *ru) * .5;
 123                 return;
 124             }
 125         } else if (ntyp == *mtyp) {
 126             a = (*fu - *fl) / (*pl * (*ru - *rl));
 127             b = *pu / *pl;
 128         }
 129         if (ntyp == 2) {
 130
 131 /*     QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL */
 132 /*     DERIVATIVE */
 133
 134             den = (1. - a) * 2.;
 135         } else if (ntyp == 3) {
 136
 137 /*     QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL */
 138 /*     DERIVATIVES */
 139
 140             den = 1. - b;
 141         } else if (ntyp == 4) {
 142
 143 /*     CUBIC EXTRAPOLATION OR INTERPOLATION */
 144
 145             c__ = b - a * 2. + 1.;
 146             d__ = b - a * 3. + 2.;
 147             dis = d__ * d__ - c__ * 3.;
 148             if (dis < 0.) {
 149                 goto L1;
 150             }
 151             den = d__ + sqrt(dis);
 152         } else if (ntyp == 5) {
 153
 154 /*     CONIC EXTRAPOLATION OR INTERPOLATION */
 155
 156             dis = a * a - b;
 157             if (dis < 0.) {
 158                 goto L1;
 159             }
 160             den = a + sqrt(dis);
 161             if (den <= 0.) {
 162                 goto L1;
 163             }
 164 /* Computing 3rd power */
 165             d__1 = 1. / den;
 166             den = 1. - b * (d__1 * (d__1 * d__1));
 167         }
 168         if (*mode == 1 && den > 0. && den < 1.) {
 169
 170 /*     EXTRAPOLATION ACCEPTED */
 171
 172             *r__ = *rl + (*ru - *rl) / den;
 173 /* Computing MAX */
 174             d__1 = *r__, d__2 = *ru * 1.1;
 175             *r__ = MAX2(d__1,d__2);
 176 /* Computing MIN */
 177             d__1 = *r__, d__2 = *ru * 1e3;
 178             *r__ = MIN2(d__1,d__2);
 179             return;
 180         } else if (*mode == 2 && den > 1.) {
 181
 182 /*     INTERPOLATION ACCEPTED */
 183
 184             *r__ = *rl + (*ru - *rl) / den;
 185             if (*rl == 0.) {
 186 /* Computing MAX */
 187                 d__1 = *r__, d__2 = *rl + (*ru - *rl) * .01;
 188                 *r__ = MAX2(d__1,d__2);
 189             } else {
 190 /* Computing MAX */
 191                 d__1 = *r__, d__2 = *rl + (*ru - *rl) * .1;
 192                 *r__ = MAX2(d__1,d__2);
 193             }
 194 /* Computing MIN */
 195             d__1 = *r__, d__2 = *rl + (*ru - *rl) * .9;
 196             *r__ = MIN2(d__1,d__2);
 197             return;
 198         }
 199 L1:
 200         ;
 201     }
 202     return;
 203 } /* luksan_pnint1__ */
 204
 205 /* save and restore state, replacing old non-reeentrant implementation
 206    that used static local variables */
 207 #define SS(var) state->var = var
 208 #define SAVE_STATE SS(fl); SS(fu); SS(pl); SS(rl); SS(pu); SS(ru); \
 209                    SS(mes1); SS(mes2); SS(mes3); SS(mode); SS(mtyp)
 210 #define RS(var) var = state->var
 211 #define RESTORE_STATE RS(fl); RS(fu); RS(pl); RS(rl); RS(pu); RS(ru); \
 212                       RS(mes1); RS(mes2); RS(mes3); RS(mode); RS(mtyp)
 213
 214
 215 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 216 /* SUBROUTINE PS1L01                ALL SYSTEMS                97/12/01
 217 * PURPOSE :
 218 *  STANDARD LINE SEARCH WITH DIRECTIONAL DERIVATIVES.
 219 *
 220 * PARAMETERS :
 221 *  RO  R  VALUE OF THE STEPSIZE PARAMETER.
 222 *  RO  RP  PREVIOUS VALUE OF THE STEPSIZE PARAMETER.
 223 *  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
 224 *  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
 225 *  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
 226 *  RO  P  VALUE OF THE DIRECTIONAL DERIVATIVE.
 227 *  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
 228 *  RO  PP  PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
 229 *  RI  FMIN  LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
 230 *  RI  MAXF  UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
 231 *  RI  RMIN  MINIMUM VALUE OF THE STEPSIZE PARAMETER
 232 *  RI  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER
 233 *  RI  TOLS  TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
 234 *         CHANGE OF THE FUNCTION VALUE).
 235 *  RI  TOLP  TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
 236 *         CHANGE OF THE DIRECTIONAL DERIVATIVE).
 237 *  RO  PAR1  PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
 238 *         UPDATES.
 239 *  RO  PAR2  PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
 240 *         UPDATES.
 241 *  II  KD  DEGREE OF REQUIRED DERIVATIVES.
 242 *  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES OF OBJECTIVE
 243 *  II  NIT  ACTUAL NUMBER OF ITERATIONS.
 244 *  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
 245 *  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
 246 *  II  MRED  MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
 247 *  IO  MAXST  MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
 248 *         STEPSIZE WAS NOT OR WAS REACHED.
 249 *  II  IEST  LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND
 250 *         IS NOT OR IS GIVEN.
 251 *  II  INITS  CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE
 252 *         IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL
 253 *         STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED
 254 *         INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL
 255 *         STEPSIZE.
 256 *  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
 257 *         LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
 258 *         STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
 259 *         ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
 260 *         ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
 261 *         ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
 262 *         DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
 263 *  II  KTERS  TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH.
 264 *         KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE.
 265 *         KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE.
 266 *         KTERS=6-FIRST STEPSIZE.
 267 *  II  MES  METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC
 268 *         INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
 269 *         MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
 270 *         DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC
 271 *         INTERPOLATION.
 272 *  IU  ISYS  CONTROL PARAMETER.
 273 *
 274 * SUBPROGRAM USED :
 275 *  S   PNINT1  EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL
 276 *         DERIVATIVES.
 277 *
 278 * METHOD :
 279 * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH STANDARD TERMINATION
 280 * CRITERIA.
 281 */
 282 void luksan_ps1l01__(double *r__, double *rp,
 283         double *f, double *fo, double *fp, double *p,
 284         double *po, double *pp, double *minf, double *maxf,
 285         double *rmin, double *rmax, double *tols, double *
 286         tolp, double *par1, double *par2, int *kd, int *ld,
 287         int *nit, int *kit, int *nred, int *mred, int *
 288         maxst, int *iest, int *inits, int *iters, int *kters,
 289         int *mes, int *isys, ps1l01_state *state)
 290 {
 291     /* System generated locals */
 292     double d__1, d__2;
 293
 294     /* Local variables */
 295     unsigned l1, l2, l3, m1, l5, m2, l7, m3;
 296     double fl, fu, pl, rl, pu, ru;
 297     int mes1, mes2, mes3, mode;
 298     int merr;
 299     int mtyp;
 300     int init1;
 301     double rtemp;
 302
 303     RESTORE_STATE;
 304
 305     if (*isys == 1) {
 306         goto L3;
 307     }
 308     mes1 = 2;
 309     mes2 = 2;
 310     mes3 = 2;
 311     *iters = 0;
 312     if (*po >= 0.) {
 313         *r__ = 0.;
 314         *iters = -2;
 315         goto L4;
 316     }
 317     if (*rmax <= 0.) {
 318         *iters = 0;
 319         goto L4;
 320     }
 321
 322 /*     INITIAL STEPSIZE SELECTION */
 323
 324     if (*inits > 0) {
 325         rtemp = *minf - *f;
 326     } else if (*iest == 0) {
 327         rtemp = *f - *fp;
 328     } else {
 329 /* Computing MAX */
 330         d__1 = *f - *fp, d__2 = *minf - *f;
 331         rtemp = MAX2(d__1,d__2);
 332     }
 333     init1 = iabs(*inits);
 334     *rp = 0.;
 335     *fp = *fo;
 336     *pp = *po;
 337     if (init1 == 0) {
 338     } else if (init1 == 1 || (*inits >= 1 && *iest == 0)) {
 339         *r__ = 1.;
 340     } else if (init1 == 2) {
 341 /* Computing MIN */
 342         d__1 = 1., d__2 = rtemp * 4. / *po;
 343         *r__ = MIN2(d__1,d__2);
 344     } else if (init1 == 3) {
 345 /* Computing MIN */
 346         d__1 = 1., d__2 = rtemp * 2. / *po;
 347         *r__ = MIN2(d__1,d__2);
 348     } else if (init1 == 4) {
 349         *r__ = rtemp * 2. / *po;
 350     }
 351     *r__ = MAX2(*r__,*rmin);
 352     *r__ = MIN2(*r__,*rmax);
 353     mode = 0;
 354     ru = 0.;
 355     fu = *fo;
 356     pu = *po;
 357
 358 /*     NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) */
 359
 360 L2:
 361     luksan_pnint1__(&rl, &ru, &fl, &fu, &pl, &pu, r__, &mode, &mtyp, &merr);
 362     if (merr > 0) {
 363         *iters = -merr;
 364         goto L4;
 365     } else if (mode == 1) {
 366         --(*nred);
 367         *r__ = MIN2(*r__,*rmax);
 368     } else if (mode == 2) {
 369         ++(*nred);
 370     }
 371
 372 /*     COMPUTATION OF THE NEW FUNCTION VALUE AND THE NEW DIRECTIONAL */
 373 /*     DERIVATIVE */
 374
 375     *kd = 1;
 376     *ld = -1;
 377     *isys = 1;
 378     SAVE_STATE;
 379     return;
 380 L3:
 381     if (mode == 0) {
 382         *par1 = *p / *po;
 383         *par2 = *f - *fo;
 384     }
 385     if (*iters != 0) {
 386         goto L4;
 387     }
 388     if (*f <= *minf) {
 389         *iters = 7;
 390         goto L4;
 391     } else {
 392         l1 = *r__ <= *rmin && *nit != *kit;
 393         l2 = *r__ >= *rmax;
 394         l3 = *f - *fo <= *tols * *r__ * *po;
 395         l5 = *p >= *tolp * *po || (mes2 == 2 && mode == 2);
 396         l7 = mes2 <= 2 || mode != 0;
 397         m1 = FALSE_;
 398         m2 = FALSE_;
 399         m3 = l3;
 400         if (mes3 >= 1) {
 401             m1 = fabs(*p) <= fabs(*po) * .01 && *fo - *f >= fabs(*fo) *
 402                     9.9999999999999994e-12;
 403             l3 = l3 || m1;
 404         }
 405         if (mes3 >= 2) {
 406             m2 = fabs(*p) <= fabs(*po) * .5 && (d__1 = *fo - *f, fabs(d__1)) <=
 407                     fabs(*fo) * 2.0000000000000001e-13;
 408             l3 = l3 || m2;
 409         }
 410         *maxst = 0;
 411         if (l2) {
 412             *maxst = 1;
 413         }
 414     }
 415
 416 /*     TEST ON TERMINATION */
 417
 418     if (l1 && ! l3) {
 419         *iters = 0;
 420         goto L4;
 421     } else if (l2 && l3 && ! l5) {
 422         *iters = 7;
 423         goto L4;
 424     } else if (m3 && mes1 == 3) {
 425         *iters = 5;
 426         goto L4;
 427     } else if (l3 && l5 && l7) {
 428         *iters = 4;
 429         goto L4;
 430     } else if (*kters < 0 || (*kters == 6 && l7)) {
 431         *iters = 6;
 432         goto L4;
 433     } else if (iabs(*nred) >= *mred) {
 434         *iters = -1;
 435         goto L4;
 436     } else {
 437         *rp = *r__;
 438         *fp = *f;
 439         *pp = *p;
 440         mode = MAX2(mode,1);
 441         mtyp = iabs(*mes);
 442         if (*f >= *maxf) {
 443             mtyp = 1;
 444         }
 445     }
 446     if (mode == 1) {
 447
 448 /*     INTERVAL CHANGE AFTER EXTRAPOLATION */
 449
 450         rl = ru;
 451         fl = fu;
 452         pl = pu;
 453         ru = *r__;
 454         fu = *f;
 455         pu = *p;
 456         if (! l3) {
 457             *nred = 0;
 458             mode = 2;
 459         } else if (mes1 == 1) {
 460             mtyp = 1;
 461         }
 462     } else {
 463
 464 /*     INTERVAL CHANGE AFTER INTERPOLATION */
 465
 466         if (! l3) {
 467             ru = *r__;
 468             fu = *f;
 469             pu = *p;
 470         } else {
 471             rl = *r__;
 472             fl = *f;
 473             pl = *p;
 474         }
 475     }
 476     goto L2;
 477 L4:
 478     *isys = 0;
 479     SAVE_STATE;
 480     return;
 481 } /* luksan_ps1l01__ */
 482
 483 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 484 /* SUBROUTINE PULSP3                ALL SYSTEMS                02/12/01
 485 * PURPOSE :
 486 * LIMITED STORAGE VARIABLE METRIC UPDATE.
 487 *
 488 * PARAMETERS :
 489 *  II  N  NUMBER OF VARIABLES (NUMBER OF ROWS OF XM).
 490 *  II  M  NUMBER OF COLUMNS OF XM.
 491 *  II  MF  MAXIMUM NUMBER OF COLUMNS OF XM.
 492 *  RI  XM(N*M)  RECTANGULAR MATRIX IN THE PRODUCT FORM SHIFTED BROYDEN
 493 *         METHOD (STORED COLUMNWISE): H-SIGMA*I=XM*TRANS(XM)
 494 *  RO  GR(M)  MATRIX TRANS(XM)*GO.
 495 *  RU  XO(N)  VECTORS OF VARIABLES DIFFERENCE XO AND VECTOR XO-TILDE.
 496 *  RU  GO(N)  GRADIENT DIFFERENCE GO AND VECTOR XM*TRANS(XM)*GO.
 497 *  RI  R  STEPSIZE PARAMETER.
 498 *  RI  PO  OLD DIRECTIONAL DERIVATIVE (MULTIPLIED BY R)
 499 *  RU  SIG  SCALING PARAMETER (ZETA AND SIGMA).
 500 *  IO  ITERH  TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
 501 *         ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
 502 *  II  MET3  CHOICE OF SIGMA (1-CONSTANT, 2-QUADRATIC EQUATION).
 503 *
 504 * SUBPROGRAMS USED :
 505 *  S   MXDRMM  MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR
 506 *         MATRIX BY A VECTOR.
 507 *  S   MXDCMU  UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX.
 508 *         WITH CONTROLLING OF POSITIVE DEFINITENESS.
 509 *  S   MXVDIR  VECTOR AUGMENTED BY A SCALED VECTOR.
 510 *  RF  MXVDOT  DOT PRODUCT OF VECTORS.
 511 *  S   MXVSCL  SCALING OF A VECTOR.
 512 *
 513 * METHOD :
 514 * SHIFTED BFGS METHOD IN THE PRODUCT FORM.
 515 */
 516 void luksan_pulsp3__(int *n, int *m, int *mf,
 517         double *xm, double *gr, double *xo, double *go,
 518         double *r__, double *po, double *sig, int *iterh,
 519         int *met3)
 520 {
 521     /* System generated locals */
 522     double d__1, d__2, d__3, d__4;
 523
 524     /* Builtin functions */
 525
 526     /* Local variables */
 527     double a, b, c__, aa, bb, ah, den, par, pom;
 528
 529     /* Parameter adjustments */
 530     --go;
 531     --xo;
 532     --gr;
 533     --xm;
 534
 535     /* Function Body */
 536     if (*m >= *mf) {
 537         return;
 538     }
 539     b = luksan_mxvdot__(n, &xo[1], &go[1]);
 540     if (b <= 0.) {
 541         *iterh = 2;
 542         goto L22;
 543     }
 544     luksan_mxdrmm__(n, m, &xm[1], &go[1], &gr[1]);
 545     ah = luksan_mxvdot__(n, &go[1], &go[1]);
 546     aa = luksan_mxvdot__(m, &gr[1], &gr[1]);
 547     a = aa + ah * *sig;
 548     c__ = -(*r__) * *po;
 549
 550 /*     DETERMINATION OF THE PARAMETER SIG (SHIFT) */
 551
 552     par = 1.;
 553     pom = b / ah;
 554     if (a > 0.) {
 555         den = luksan_mxvdot__(n, &xo[1], &xo[1]);
 556         if (*met3 <= 4) {
 557 /* Computing MAX */
 558             d__1 = 0., d__2 = 1. - aa / a;
 559 /* Computing MAX */
 560             d__3 = 0., d__4 = 1. - b * b / (den * ah);
 561             *sig = sqrt((MAX2(d__1,d__2))) / (sqrt((MAX2(d__3,d__4))) + 1.) *
 562                     pom;
 563         } else {
 564 /* Computing MAX */
 565             d__1 = 0., d__2 = *sig * ah / a;
 566 /* Computing MAX */
 567             d__3 = 0., d__4 = 1. - b * b / (den * ah);
 568             *sig = sqrt((MAX2(d__1,d__2))) / (sqrt((MAX2(d__3,d__4))) + 1.) *
 569                     pom;
 570         }
 571 /* Computing MAX */
 572         d__1 = *sig, d__2 = pom * .2;
 573         *sig = MAX2(d__1,d__2);
 574 /* Computing MIN */
 575         d__1 = *sig, d__2 = pom * .8;
 576         *sig = MIN2(d__1,d__2);
 577     } else {
 578         *sig = pom * .25;
 579     }
 580
 581 /*     COMPUTATION OF SHIFTED XO AND SHIFTED B */
 582
 583     bb = b - ah * *sig;
 584     d__1 = -(*sig);
 585     luksan_mxvdir__(n, &d__1, &go[1], &xo[1], &xo[1]);
 586
 587 /*     BFGS-BASED SHIFTED BFGS UPDATE */
 588
 589     pom = 1.;
 590     d__1 = -1. / bb;
 591     luksan_mxdcmu__(n, m, &xm[1], &d__1, &xo[1], &gr[1]);
 592     d__1 = sqrt(par / bb);
 593     luksan_mxvscl__(n, &d__1, &xo[1], &xm[*n * *m + 1]);
 594     ++(*m);
 595 L22:
 596     *iterh = 0;
 597     return;
 598 } /* luksan_pulsp3__ */
 599
 600 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 601 /* SUBROUTINE PULVP3                ALL SYSTEMS                03/12/01
 602 * PURPOSE :
 603 * RANK-TWO LIMITED-STORAGE VARIABLE-METRIC METHODS IN THE PRODUCT FORM.
 604 *
 605 * PARAMETERS :
 606 *  II  N  NUMBER OF VARIABLES (NUMBER OF ROWS OF XM).
 607 *  II  M  NUMBER OF COLUMNS OF XM.
 608 *  RI  XM(N*M)  RECTANGULAR MATRIX IN THE PRODUCT FORM SHIFTED BROYDEN
 609 *         METHOD (STORED COLUMNWISE): H-SIGMA*I=XM*TRANS(XM)
 610 *  RO  XR(M)  VECTOR TRANS(XM)*H**(-1)*XO.
 611 *  RO  GR(M)  MATRIX TRANS(XM)*GO.
 612 *  RA  S(N)  AUXILIARY VECTORS (H**(-1)*XO AND U).
 613 *  RA  SO(N)  AUXILIARY VECTORS ((H-SIGMA*I)*H**(-1)*XO AND V).
 614 *  RU  XO(N)  VECTORS OF VARIABLES DIFFERENCE XO AND VECTOR XO-TILDE.
 615 *  RU  GO(N)  GRADIENT DIFFERENCE GO AND VECTOR XM*TRANS(XM)*GO.
 616 *  RI  R  STEPSIZE PARAMETER.
 617 *  RI  PO  OLD DIRECTIONAL DERIVATIVE (MULTIPLIED BY R)
 618 *  RU  SIG  SCALING PARAMETER (ZETA AND SIGMA).
 619 *  IO  ITERH  TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
 620 *         ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
 621 *  II  MET2  CHOICE OF THE CORRECTION PARAMETER (1-THE UNIT VALUE,
 622 *         2-THE BALANCING VALUE, 3-THE SQUARE ROOT, 4-THE GEOMETRIC
 623 *         MEAN).
 624 *  II  MET3  CHOICE OF THE SHIFT PARAMETER (4-THE FIRST FORMULA,
 625 *         5-THE SECOND FORMULA).
 626 *  II  MET5  CHOICE OF THE METHOD (1-RANK-ONE METHOD, 2-RANK-TWO
 627 *         METHOD).
 628 *
 629 * SUBPROGRAMS USED :
 630 *  S   MXDRMM  MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR
 631 *         MATRIX BY A VECTOR.
 632 *  S   MXDCMU  UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX.
 633 *         WITH CONTROLLING OF POSITIVE DEFINITENESS. RANK-ONE FORMULA.
 634 *  S   MXDCMV  UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX.
 635 *         WITH CONTROLLING OF POSITIVE DEFINITENESS. RANK-TWO FORMULA.
 636 *  S   MXVDIR  VECTOR AUGMENTED BY A SCALED VECTOR.
 637 *  RF  MXVDOT  DOT PRODUCT OF VECTORS.
 638 *  S   MXVLIN  LINEAR COMBINATION OF TWO VECTORS.
 639 *  S   MXVSCL  SCALING OF A VECTOR.
 640 *
 641 * METHOD :
 642 * RANK-ONE LIMITED-STORAGE VARIABLE-METRIC METHOD IN THE PRODUCT FORM.
 643 */
 644 void luksan_pulvp3__(int *n, int *m, double *xm,
 645         double *xr, double *gr, double *s, double *so,
 646         double *xo, double *go, double *r__, double *po,
 647         double *sig, int *iterh, int *met2, int *met3,
 648         int *met5)
 649 {
 650     /* System generated locals */
 651     double d__1, d__2, d__3, d__4;
 652
 653     /* Builtin functions */
 654
 655     /* Local variables */
 656     double a, b, c__, aa, bb, cc, ah, den, par, pom, zet;
 657
 658     /* Parameter adjustments */
 659     --go;
 660     --xo;
 661     --so;
 662     --s;
 663     --gr;
 664     --xr;
 665     --xm;
 666
 667     /* Function Body */
 668     zet = *sig;
 669
 670 /*     COMPUTATION OF B */
 671
 672     b = luksan_mxvdot__(n, &xo[1], &go[1]);
 673     if (b <= 0.) {
 674         *iterh = 2;
 675         goto L22;
 676     }
 677
 678 /*     COMPUTATION OF GR=TRANS(XM)*GO, XR=TRANS(XM)*H**(-1)*XO */
 679 /*     AND S=H**(-1)*XO, SO=(H-SIGMA*I)*H**(-1)*XO. COMPUTATION */
 680 /*     OF AA=GR*GR, BB=GR*XR, CC=XR*XR. COMPUTATION OF A AND C. */
 681
 682     luksan_mxdrmm__(n, m, &xm[1], &go[1], &gr[1]);
 683     luksan_mxvscl__(n, r__, &s[1], &s[1]);
 684     luksan_mxdrmm__(n, m, &xm[1], &s[1], &xr[1]);
 685     d__1 = -(*sig);
 686     luksan_mxvdir__(n, &d__1, &s[1], &xo[1], &so[1]);
 687     ah = luksan_mxvdot__(n, &go[1], &go[1]);
 688     aa = luksan_mxvdot__(m, &gr[1], &gr[1]);
 689     bb = luksan_mxvdot__(m, &gr[1], &xr[1]);
 690     cc = luksan_mxvdot__(m, &xr[1], &xr[1]);
 691     a = aa + ah * *sig;
 692     c__ = -(*r__) * *po;
 693
 694 /*     DETERMINATION OF THE PARAMETER SIG (SHIFT) */
 695
 696     pom = b / ah;
 697     if (a > 0.) {
 698         den = luksan_mxvdot__(n, &xo[1], &xo[1]);
 699         if (*met3 <= 4) {
 700 /* Computing MAX */
 701             d__1 = 0., d__2 = 1. - aa / a;
 702 /* Computing MAX */
 703             d__3 = 0., d__4 = 1. - b * b / (den * ah);
 704             *sig = sqrt((MAX2(d__1,d__2))) / (sqrt((MAX2(d__3,d__4))) + 1.) *
 705                     pom;
 706         } else {
 707 /* Computing MAX */
 708             d__1 = 0., d__2 = *sig * ah / a;
 709 /* Computing MAX */
 710             d__3 = 0., d__4 = 1. - b * b / (den * ah);
 711             *sig = sqrt((MAX2(d__1,d__2))) / (sqrt((MAX2(d__3,d__4))) + 1.) *
 712                     pom;
 713         }
 714 /* Computing MAX */
 715         d__1 = *sig, d__2 = pom * .2;
 716         *sig = MAX2(d__1,d__2);
 717 /* Computing MIN */
 718         d__1 = *sig, d__2 = pom * .8;
 719         *sig = MIN2(d__1,d__2);
 720     } else {
 721         *sig = pom * .25;
 722     }
 723
 724 /*     COMPUTATION OF SHIFTED XO AND SHIFTED B */
 725
 726     b -= ah * *sig;
 727     d__1 = -(*sig);
 728     luksan_mxvdir__(n, &d__1, &go[1], &xo[1], &xo[1]);
 729
 730 /*     COMPUTATION OF THE PARAMETER RHO (CORRECTION) */
 731
 732     if (*met2 <= 1) {
 733         par = 1.;
 734     } else if (*met2 == 2) {
 735         par = *sig * ah / b;
 736     } else if (*met2 == 3) {
 737         par = sqrt(1. - aa / a);
 738     } else if (*met2 == 4) {
 739         par = sqrt(sqrt(1. - aa / a) * (*sig * ah / b));
 740     } else {
 741         par = zet / (zet + *sig);
 742     }
 743
 744 /*     COMPUTATION OF THE PARAMETER THETA (BFGS) */
 745
 746     d__1 = sqrt(par * b / cc);
 747     pom = copysign(d__1, bb);
 748
 749 /*     COMPUTATION OF Q AND P */
 750
 751     if (*met5 == 1) {
 752
 753 /*     RANK ONE UPDATE OF XM */
 754
 755         luksan_mxvdir__(m, &pom, &xr[1], &gr[1], &xr[1]);
 756         luksan_mxvlin__(n, &par, &xo[1], &pom, &so[1], &s[1]);
 757         d__1 = -1. / (par * b + pom * bb);
 758         luksan_mxdcmu__(n, m, &xm[1], &d__1, &s[1], &xr[1]);
 759     } else {
 760
 761 /*     RANK TWO UPDATE OF XM */
 762
 763         d__1 = par / pom - bb / b;
 764         luksan_mxvdir__(n, &d__1, &xo[1], &so[1], &s[1]);
 765         d__1 = -1. / b;
 766         d__2 = -1. / cc;
 767         luksan_mxdcmv__(n, m, &xm[1], &d__1, &xo[1], &gr[1], &d__2, &s[1], &xr[1]);
 768     }
 769 L22:
 770     *iterh = 0;
 771     return;
 772 } /* luksan_pulvp3__ */
 773
 774 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 775 /* SUBROUTINE PYADC0                ALL SYSTEMS                98/12/01
 776 * PURPOSE :
 777 * NEW SIMPLE BOUNDS ARE ADDED TO THE ACTIVE SET.
 778 *
 779 * PARAMETERS :
 780 *  II  NF  DECLARED NUMBER OF VARIABLES.
 781 *  II  N  REDUCED NUMBER OF VARIABLES.
 782 *  RI  X(NF)  VECTOR OF VARIABLES.
 783 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
 784 *  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
 785 *  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
 786 *  IO  INEW  NUMBER OF ACTIVE CONSTRAINTS.
 787 */
 788 void luksan_pyadc0__(int *nf, int *n, double *x,
 789         int *ix, double *xl, double *xu, int *inew)
 790 {
 791     /* System generated locals */
 792     int i__1;
 793
 794     /* Local variables */
 795     int i__, ii, ixi;
 796
 797     /* Parameter adjustments */
 798     --ix;
 799     --x;
 800     --xl;
 801     --xu;
 802
 803     /* Function Body */
 804     *n = *nf;
 805     *inew = 0;
 806     i__1 = *nf;
 807     for (i__ = 1; i__ <= i__1; ++i__) {
 808         ii = ix[i__];
 809         ixi = iabs(ii);
 810         if (ixi >= 5) {
 811             ix[i__] = -ixi;
 812         } else if ((ixi == 1 || ixi == 3 || ixi == 4) && x[i__] <= xl[i__]) {
 813             x[i__] = xl[i__];
 814             if (ixi == 4) {
 815                 ix[i__] = -3;
 816             } else {
 817                 ix[i__] = -ixi;
 818             }
 819             --(*n);
 820             if (ii > 0) {
 821                 ++(*inew);
 822             }
 823         } else if ((ixi == 2 || ixi == 3 || ixi == 4) && x[i__] >= xu[i__]) {
 824             x[i__] = xu[i__];
 825             if (ixi == 3) {
 826                 ix[i__] = -4;
 827             } else {
 828                 ix[i__] = -ixi;
 829             }
 830             --(*n);
 831             if (ii > 0) {
 832                 ++(*inew);
 833             }
 834         }
 835 /* L1: */
 836     }
 837     return;
 838 } /* luksan_pyadc0__ */
 839
 840 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 841 /* SUBROUTINE PYFUT1                ALL SYSTEMS                98/12/01
 842 * PURPOSE :
 843 * TERMINATION CRITERIA AND TEST ON RESTART.
 844 *
 845 * PARAMETERS :
 846 *  II  N  ACTUAL NUMBER OF VARIABLES.
 847 *  RI  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
 848 *  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
 849 *  RI  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
 850 *  RO  GMAX  NORM OF THE TRANSFORMED GRADIENT.
 851 *  RI  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
 852 *  RI  TOLX  LOWER BOUND FOR STEPLENGTH.
 853 *  RI  TOLF  LOWER BOUND FOR FUNCTION DECREASE.
 854 *  RI  TOLB  LOWER BOUND FOR FUNCTION VALUE.
 855 *  RI  TOLG  LOWER BOUND FOR GRADIENT.
 856 *  II  KD  DEGREE OF REQUIRED DERIVATIVES.
 857 *  IU  NIT  ACTUAL NUMBER OF ITERATIONS.
 858 *  II  KIT  NUMBER OF THE ITERATION AFTER RESTART.
 859 *  II  MIT  MAXIMUM NUMBER OF ITERATIONS.
 860 *  IU  NFV  ACTUAL NUMBER OF COMPUTED FUNCTION VALUES.
 861 *  II  MFV  MAXIMUM NUMBER OF COMPUTED FUNCTION VALUES.
 862 *  IU  NFG  ACTUAL NUMBER OF COMPUTED GRADIENT VALUES.
 863 *  II  MFG  MAXIMUM NUMBER OF COMPUTED GRADIENT VALUES.
 864 *  IU  NTESX  ACTUAL NUMBER OF TESTS ON STEPLENGTH.
 865 *  II  MTESX  MAXIMUM NUMBER OF TESTS ON STEPLENGTH.
 866 *  IU  NTESF  ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE.
 867 *  II  MTESF  MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE.
 868 *  II  IRES1  RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
 869 *         IRES1*N+IRES2 ITERATIONS.
 870 *  II  IRES2  RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
 871 *         IRES1*N+IRES2 ITERATIONS.
 872 *  IU  IREST  RESTART INDICATOR. RESTART IS PERFORMED IF IREST>0.
 873 *  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
 874 *         ITERS=0 FOR ZERO STEP.
 875 *  IO  ITERM  TERMINATION INDICATOR. ITERM=1-TERMINATION AFTER MTESX
 876 *         UNSUFFICIENT STEPLENGTHS. ITERM=2-TERMINATION AFTER MTESF
 877 *         UNSUFFICIENT FUNCTION DECREASES. ITERM=3-TERMINATION ON LOWER
 878 *         BOUND FOR FUNCTION VALUE. ITERM=4-TERMINATION ON LOWER BOUND
 879 *         FOR GRADIENT. ITERM=11-TERMINATION AFTER MAXIMUM NUMBER OF
 880 *         ITERATIONS. ITERM=12-TERMINATION AFTER MAXIMUM NUMBER OF
 881 *         COMPUTED FUNCTION VALUES.
 882 */
 883 void luksan_pyfut1__(int *n, double *f, double *fo, double *umax,
 884                      double *gmax, int xstop, /* double *dmax__,  */
 885                      const nlopt_stopping *stop,
 886                      double *tolg, int *kd, int *nit, int *kit, int *mit,
 887                      int *nfg, int *mfg, int *ntesx,
 888         int *mtesx, int *ntesf, int *mtesf, int *ites,
 889         int *ires1, int *ires2, int *irest, int *iters,
 890         int *iterm)
 891 {
 892     /* System generated locals */
 893     double d__1, d__2;
 894
 895     /* Builtin functions */
 896
 897     if (*iterm < 0) {
 898         return;
 899     }
 900     if (*ites <= 0) {
 901         goto L1;
 902     }
 903     if (*iters == 0) {
 904         goto L1;
 905     }
 906     if (*nit <= 0) {
 907 /* Computing MIN */
 908         d__1 = sqrt((fabs(*f))), d__2 = fabs(*f) / 10.;
 909         *fo = *f + MIN2(d__1,d__2);
 910     }
 911     if (nlopt_stop_forced(stop)) {
 912         *iterm = -999;
 913         return;
 914     }
 915     if (*f <= stop->minf_max /* *tolb */) {
 916         *iterm = 3;
 917         return;
 918     }
 919     if (*kd > 0) {
 920         if (*gmax <= *tolg && *umax <= *tolg) {
 921             *iterm = 4;
 922             return;
 923         }
 924     }
 925     if (*nit <= 0) {
 926         *ntesx = 0;
 927         *ntesf = 0;
 928     }
 929     if (xstop) /* (*dmax__ <= *tolx) */ {
 930         *iterm = 1;
 931         ++(*ntesx);
 932         if (*ntesx >= *mtesx) {
 933             return;
 934         }
 935     } else {
 936         *ntesx = 0;
 937     }
 938     if (nlopt_stop_ftol(stop, *f, *fo)) {
 939         *iterm = 2;
 940         ++(*ntesf);
 941         if (*ntesf >= *mtesf) {
 942             return;
 943         }
 944     } else {
 945         *ntesf = 0;
 946     }
 947 L1:
 948     if (*nit >= *mit) {
 949         *iterm = 11;
 950         return;
 951     }
 952     if (nlopt_stop_evals(stop)) /* (*nfv >= *mfv) */ {
 953         *iterm = 12;
 954         return;
 955     }
 956     if (*nfg >= *mfg) {
 957         *iterm = 13;
 958         return;
 959     }
 960     *iterm = 0;
 961     if (*n > 0 && *nit - *kit >= *ires1 * *n + *ires2) {
 962         *irest = MAX2(*irest,1);
 963     }
 964     ++(*nit);
 965     return;
 966 } /* luksan_pyfut1__ */
 967
 968 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 969 /* SUBROUTINE PYRMC0                ALL SYSTEMS                98/12/01
 970 * PURPOSE :
 971 * OLD SIMPLE BOUND IS REMOVED FROM THE ACTIVE SET. TRANSFORMED
 972 * GRADIENT OF THE OBJECTIVE FUNCTION IS UPDATED.
 973 *
 974 * PARAMETERS :
 975 *  II  NF  DECLARED NUMBER OF VARIABLES.
 976 *  II  N  REDUCED NUMBER OF VARIABLES.
 977 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
 978 *  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
 979 *  RI  EPS8  TOLERANCE FOR CONSTRAINT TO BE REMOVED.
 980 *  RI  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
 981 *  RI  GMAX  NORM OF THE TRANSFORMED GRADIENT.
 982 *  RO  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
 983 *  II  IOLD  NUMBER OF REMOVED CONSTRAINTS.
 984 *  IU  IREST  RESTART INDICATOR.
 985 */
 986 void luksan_pyrmc0__(int *nf, int *n, int *ix,
 987         double *g, double *eps8, double *umax, double *gmax,
 988         double *rmax, int *iold, int *irest)
 989 {
 990     /* System generated locals */
 991     int i__1, i__2, i__3;
 992
 993     /* Local variables */
 994     int i__, ixi;
 995
 996     /* Parameter adjustments */
 997     --g;
 998     --ix;
 999
1000     /* Function Body */
1001     if (*n == 0 || *rmax > 0.) {
1002         if (*umax > *eps8 * *gmax) {
1003             *iold = 0;
1004             i__1 = *nf;
1005             for (i__ = 1; i__ <= i__1; ++i__) {
1006                 ixi = ix[i__];
1007                 if (ixi >= 0) {
1008                 } else if (ixi <= -5) {
1009                 } else if ((ixi == -1 || ixi == -3) && -g[i__] <= 0.) {
1010                 } else if ((ixi == -2 || ixi == -4) && g[i__] <= 0.) {
1011                 } else {
1012                     ++(*iold);
1013 /* Computing MIN */
1014                     i__3 = (i__2 = ix[i__], iabs(i__2));
1015                     ix[i__] = MIN2(i__3,3);
1016                     if (*rmax == 0.) {
1017                         goto L2;
1018                     }
1019                 }
1020 /* L1: */
1021             }
1022 L2:
1023             if (*iold > 1) {
1024                 *irest = MAX2(*irest,1);
1025             }
1026         }
1027     }
1028     return;
1029 } /* luksan_pyrmc0__ */
1030
1031 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1032 /* SUBROUTINE PYTRCD             ALL SYSTEMS                   98/12/01
1033 * PURPOSE :
1034 * VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
1035 * AND SCALED AND REDUCED. TEST VALUE DMAX IS DETERMINED.
1036 *
1037 * PARAMETERS :
1038 *  II  NF DECLARED NUMBER OF VARIABLES.
1039 *  RI  X(NF)  VECTOR OF VARIABLES.
1040 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
1041 *  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
1042 *  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
1043 *  RU  GO(NF)  GRADIENTS DIFFERENCE.
1044 *  RO  R  VALUE OF THE STEPSIZE PARAMETER.
1045 *  RO  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
1046 *  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
1047 *  RO  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
1048 *  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
1049 *  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
1050 *  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
1051 *         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
1052 *  IO  KD  DEGREE OF REQUIRED DERIVATIVES.
1053 *  IO  LD  DEGREE OF COMPUTED DERIVATIVES.
1054 *  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
1055 *         ITERS=0 FOR ZERO STEP.
1056 *
1057 * SUBPROGRAMS USED :
1058 *  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
1059 *  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
1060 *         SUBSTRACTED ONE.
1061 */
1062 void luksan_pytrcd__(int *nf, double *x, int *ix,
1063         double *xo, double *g, double *go, double *r__,
1064         double *f, double *fo, double *p, double *po,
1065         double *dmax__, int *kbf, int *kd, int *ld, int *
1066         iters)
1067 {
1068     /* System generated locals */
1069     int i__1;
1070     double d__1, d__2, d__3, d__4, d__5;
1071
1072     /* Local variables */
1073     int i__;
1074
1075     /* Parameter adjustments */
1076     --go;
1077     --g;
1078     --xo;
1079     --ix;
1080     --x;
1081
1082     /* Function Body */
1083     if (*iters > 0) {
1084         luksan_mxvdif__(nf, &x[1], &xo[1], &xo[1]);
1085         luksan_mxvdif__(nf, &g[1], &go[1], &go[1]);
1086         *po = *r__ * *po;
1087         *p = *r__ * *p;
1088     } else {
1089         *f = *fo;
1090         *p = *po;
1091         luksan_mxvsav__(nf, &x[1], &xo[1]);
1092         luksan_mxvsav__(nf, &g[1], &go[1]);
1093         *ld = *kd;
1094     }
1095     *dmax__ = 0.;
1096     i__1 = *nf;
1097     for (i__ = 1; i__ <= i__1; ++i__) {
1098         if (*kbf > 0) {
1099             if (ix[i__] < 0) {
1100                 xo[i__] = 0.;
1101                 go[i__] = 0.;
1102                 goto L1;
1103             }
1104         }
1105 /* Computing MAX */
1106 /* Computing MAX */
1107         d__5 = (d__2 = x[i__], fabs(d__2));
1108         d__3 = *dmax__, d__4 = (d__1 = xo[i__], fabs(d__1)) / MAX2(d__5,1.);
1109         *dmax__ = MAX2(d__3,d__4);
1110 L1:
1111         ;
1112     }
1113     return;
1114 } /* luksan_pytrcd__ */
1115
1116 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1117 /* SUBROUTINE PYTRCG                ALL SYSTEMS                99/12/01
1118 * PURPOSE :
1119 *  GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. TEST VALUES
1120 *  GMAX AND UMAX ARE COMPUTED.
1121 *
1122 * PARAMETERS :
1123 *  II  NF DECLARED NUMBER OF VARIABLES.
1124 *  II  N  ACTUAL NUMBER OF VARIABLES.
1125 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
1126 *  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
1127 *  RI  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
1128 *  RI  GMAX  NORM OF THE TRANSFORMED GRADIENT.
1129 *  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
1130 *         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
1131 *  II  IOLD  INDEX OF THE REMOVED CONSTRAINT.
1132 *
1133 * SUBPROGRAMS USED :
1134 *  RF  MXVMAX  L-INFINITY NORM OF A VECTOR.
1135 */
1136 void luksan_pytrcg__(int *nf, int *n, int *ix,
1137         double *g, double *umax, double *gmax, int *kbf,
1138         int *iold)
1139 {
1140     /* System generated locals */
1141     int i__1;
1142     double d__1, d__2;
1143
1144     /* Local variables */
1145     int i__;
1146     double temp;
1147
1148     /* Parameter adjustments */
1149     --g;
1150     --ix;
1151
1152     /* Function Body */
1153     if (*kbf > 0) {
1154         *gmax = 0.;
1155         *umax = 0.;
1156         *iold = 0;
1157         i__1 = *nf;
1158         for (i__ = 1; i__ <= i__1; ++i__) {
1159             temp = g[i__];
1160             if (ix[i__] >= 0) {
1161 /* Computing MAX */
1162                 d__1 = *gmax, d__2 = fabs(temp);
1163                 *gmax = MAX2(d__1,d__2);
1164             } else if (ix[i__] <= -5) {
1165             } else if ((ix[i__] == -1 || ix[i__] == -3) && *umax + temp >= 0.)
1166                      {
1167             } else if ((ix[i__] == -2 || ix[i__] == -4) && *umax - temp >= 0.)
1168                      {
1169             } else {
1170                 *iold = i__;
1171                 *umax = fabs(temp);
1172             }
1173 /* L1: */
1174         }
1175     } else {
1176         *umax = 0.;
1177         *gmax = luksan_mxvmax__(nf, &g[1]);
1178     }
1179     *n = *nf;
1180     return;
1181 } /* luksan_pytrcg__ */
1182
1183 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1184 /* SUBROUTINE PYTRCS             ALL SYSTEMS                   98/12/01
1185 * PURPOSE :
1186 * SCALED AND REDUCED DIRECTION VECTOR IS BACK TRANSFORMED. VECTORS
1187 * X,G AND VALUES F,P ARE SAVED.
1188 *
1189 * PARAMETERS :
1190 *  II  NF DECLARED NUMBER OF VARIABLES.
1191 *  RI  X(NF)  VECTOR OF VARIABLES.
1192 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
1193 *  RO  XO(NF)  SAVED VECTOR OF VARIABLES.
1194 *  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
1195 *  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
1196 *  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
1197 *  RO  GO(NF)  SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
1198 *  RO  S(NF)  DIRECTION VECTOR.
1199 *  RO  RO  SAVED VALUE OF THE STEPSIZE PARAMETER.
1200 *  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
1201 *  RU  FO  SAVED VALUE OF THE OBJECTIVE FUNCTION.
1202 *  RI  F  VALUE OF THE OBJECTIVE FUNCTION.
1203 *  RO  PO  SAVED VALUE OF THE DIRECTIONAL DERIVATIVE.
1204 *  RI  P  VALUE OF THE DIRECTIONAL DERIVATIVE.
1205 *  RO  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
1206 *  RI  ETA9  MAXIMUM FOR REAL NUMBERS.
1207 *  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
1208 *         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
1209 *
1210 * SUBPROGRAMS USED :
1211 *  S   MXVCOP  COPYING OF A VECTOR.
1212 */
1213 void luksan_pytrcs__(int *nf, double *x, int *ix,
1214         double *xo, double *xl, double *xu, double *g,
1215         double *go, double *s, double *ro, double *fp,
1216         double *fo, double *f, double *po, double *p,
1217         double *rmax, double *eta9, int *kbf)
1218 {
1219     /* System generated locals */
1220     int i__1;
1221     double d__1, d__2;
1222
1223     /* Local variables */
1224     int i__;
1225
1226     /* Parameter adjustments */
1227     --s;
1228     --go;
1229     --g;
1230     --xu;
1231     --xl;
1232     --xo;
1233     --ix;
1234     --x;
1235
1236     /* Function Body */
1237     *fp = *fo;
1238     *ro = 0.;
1239     *fo = *f;
1240     *po = *p;
1241     luksan_mxvcop__(nf, &x[1], &xo[1]);
1242     luksan_mxvcop__(nf, &g[1], &go[1]);
1243     if (*kbf > 0) {
1244         i__1 = *nf;
1245         for (i__ = 1; i__ <= i__1; ++i__) {
1246             if (ix[i__] < 0) {
1247                 s[i__] = 0.;
1248             } else {
1249                 if (ix[i__] == 1 || ix[i__] >= 3) {
1250                     if (s[i__] < -1. / *eta9) {
1251 /* Computing MIN */
1252                         d__1 = *rmax, d__2 = (xl[i__] - x[i__]) / s[i__];
1253                         *rmax = MIN2(d__1,d__2);
1254                     }
1255                 }
1256                 if (ix[i__] == 2 || ix[i__] >= 3) {
1257                     if (s[i__] > 1. / *eta9) {
1258 /* Computing MIN */
1259                         d__1 = *rmax, d__2 = (xu[i__] - x[i__]) / s[i__];
1260                         *rmax = MIN2(d__1,d__2);
1261                     }
1262                 }
1263             }
1264 /* L1: */
1265         }
1266     }
1267     return;
1268 } /* luksan_pytrcs__ */
1269