luksan/pssubs.c

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