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, 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__ = 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 /* 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 = MAX2(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__ = MIN2(d__1,d__2);
 332     } else if (init1 == 3) {
 333 /* Computing MIN */
 334         d__1 = 1., d__2 = rtemp * 2. / *po;
 335         *r__ = MIN2(d__1,d__2);
 336     } else if (init1 == 4) {
 337         *r__ = rtemp * 2. / *po;
 338     }
 339     *r__ = MAX2(*r__,*rmin);
 340     *r__ = MIN2(*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__ = MIN2(*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 = MAX2(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((MAX2(d__1,d__2))) / (sqrt((MAX2(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((MAX2(d__1,d__2))) / (sqrt((MAX2(d__3,d__4))) + 1.) *
 555                     pom;
 556         }
 557 /* Computing MAX */
 558         d__1 = *sig, d__2 = pom * .2;
 559         *sig = MAX2(d__1,d__2);
 560 /* Computing MIN */
 561         d__1 = *sig, d__2 = pom * .8;
 562         *sig = MIN2(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((MAX2(d__1,d__2))) / (sqrt((MAX2(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((MAX2(d__1,d__2))) / (sqrt((MAX2(d__3,d__4))) + 1.) *
 698                     pom;
 699         }
 700 /* Computing MAX */
 701         d__1 = *sig, d__2 = pom * .2;
 702         *sig = MAX2(d__1,d__2);
 703 /* Computing MIN */
 704         d__1 = *sig, d__2 = pom * .8;
 705         *sig = MIN2(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 *fo, double *umax,
 870                      double *gmax, int xstop, /* double *dmax__,  */
 871                      const nlopt_stopping *stop,
 872                      double *tolg, int *kd, int *nit, int *kit, int *mit,
 873                      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     if (*iterm < 0) {
 884         return;
 885     }
 886     if (*ites <= 0) {
 887         goto L1;
 888     }
 889     if (*iters == 0) {
 890         goto L1;
 891     }
 892     if (*nit <= 0) {
 893 /* Computing MIN */
 894         d__1 = sqrt((fabs(*f))), d__2 = fabs(*f) / 10.;
 895         *fo = *f + MIN2(d__1,d__2);
 896     }
 897     if (nlopt_stop_forced(stop)) {
 898         *iterm = -999;
 899         return;
 900     }
 901     if (*f <= stop->minf_max /* *tolb */) {
 902         *iterm = 3;
 903         return;
 904     }
 905     if (*kd > 0) {
 906         if (*gmax <= *tolg && *umax <= *tolg) {
 907             *iterm = 4;
 908             return;
 909         }
 910     }
 911     if (*nit <= 0) {
 912         *ntesx = 0;
 913         *ntesf = 0;
 914     }
 915     if (xstop) /* (*dmax__ <= *tolx) */ {
 916         *iterm = 1;
 917         ++(*ntesx);
 918         if (*ntesx >= *mtesx) {
 919             return;
 920         }
 921     } else {
 922         *ntesx = 0;
 923     }
 924     if (nlopt_stop_ftol(stop, *f, *fo)) {
 925         *iterm = 2;
 926         ++(*ntesf);
 927         if (*ntesf >= *mtesf) {
 928             return;
 929         }
 930     } else {
 931         *ntesf = 0;
 932     }
 933 L1:
 934     if (*nit >= *mit) {
 935         *iterm = 11;
 936         return;
 937     }
 938     if (nlopt_stop_evals(stop)) /* (*nfv >= *mfv) */ {
 939         *iterm = 12;
 940         return;
 941     }
 942     if (*nfg >= *mfg) {
 943         *iterm = 13;
 944         return;
 945     }
 946     *iterm = 0;
 947     if (*n > 0 && *nit - *kit >= *ires1 * *n + *ires2) {
 948         *irest = MAX2(*irest,1);
 949     }
 950     ++(*nit);
 951     return;
 952 } /* luksan_pyfut1__ */
 953
 954 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 955 /* SUBROUTINE PYRMC0                ALL SYSTEMS                98/12/01
 956 * PURPOSE :
 957 * OLD SIMPLE BOUND IS REMOVED FROM THE ACTIVE SET. TRANSFORMED
 958 * GRADIENT OF THE OBJECTIVE FUNCTION IS UPDATED.
 959 *
 960 * PARAMETERS :
 961 *  II  NF  DECLARED NUMBER OF VARIABLES.
 962 *  II  N  REDUCED NUMBER OF VARIABLES.
 963 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
 964 *  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
 965 *  RI  EPS8  TOLERANCE FOR CONSTRAINT TO BE REMOVED.
 966 *  RI  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
 967 *  RI  GMAX  NORM OF THE TRANSFORMED GRADIENT.
 968 *  RO  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
 969 *  II  IOLD  NUMBER OF REMOVED CONSTRAINTS.
 970 *  IU  IREST  RESTART INDICATOR.
 971 */
 972 void luksan_pyrmc0__(int *nf, int *n, int *ix,
 973         double *g, double *eps8, double *umax, double *gmax,
 974         double *rmax, int *iold, int *irest)
 975 {
 976     /* System generated locals */
 977     int i__1, i__2, i__3;
 978
 979     /* Local variables */
 980     int i__, ixi;
 981
 982     /* Parameter adjustments */
 983     --g;
 984     --ix;
 985
 986     /* Function Body */
 987     if (*n == 0 || *rmax > 0.) {
 988         if (*umax > *eps8 * *gmax) {
 989             *iold = 0;
 990             i__1 = *nf;
 991             for (i__ = 1; i__ <= i__1; ++i__) {
 992                 ixi = ix[i__];
 993                 if (ixi >= 0) {
 994                 } else if (ixi <= -5) {
 995                 } else if ((ixi == -1 || ixi == -3) && -g[i__] <= 0.) {
 996                 } else if ((ixi == -2 || ixi == -4) && g[i__] <= 0.) {
 997                 } else {
 998                     ++(*iold);
 999 /* Computing MIN */
1000                     i__3 = (i__2 = ix[i__], iabs(i__2));
1001                     ix[i__] = MIN2(i__3,3);
1002                     if (*rmax == 0.) {
1003                         goto L2;
1004                     }
1005                 }
1006 /* L1: */
1007             }
1008 L2:
1009             if (*iold > 1) {
1010                 *irest = MAX2(*irest,1);
1011             }
1012         }
1013     }
1014     return;
1015 } /* luksan_pyrmc0__ */
1016
1017 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1018 /* SUBROUTINE PYTRCD             ALL SYSTEMS                   98/12/01
1019 * PURPOSE :
1020 * VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
1021 * AND SCALED AND REDUCED. TEST VALUE DMAX IS DETERMINED.
1022 *
1023 * PARAMETERS :
1024 *  II  NF DECLARED NUMBER OF VARIABLES.
1025 *  RI  X(NF)  VECTOR OF VARIABLES.
1026 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
1027 *  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
1028 *  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
1029 *  RU  GO(NF)  GRADIENTS DIFFERENCE.
1030 *  RO  R  VALUE OF THE STEPSIZE PARAMETER.
1031 *  RO  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
1032 *  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
1033 *  RO  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
1034 *  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
1035 *  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
1036 *  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
1037 *         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
1038 *  IO  KD  DEGREE OF REQUIRED DERIVATIVES.
1039 *  IO  LD  DEGREE OF COMPUTED DERIVATIVES.
1040 *  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
1041 *         ITERS=0 FOR ZERO STEP.
1042 *
1043 * SUBPROGRAMS USED :
1044 *  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
1045 *  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
1046 *         SUBSTRACTED ONE.
1047 */
1048 void luksan_pytrcd__(int *nf, double *x, int *ix,
1049         double *xo, double *g, double *go, double *r__,
1050         double *f, double *fo, double *p, double *po,
1051         double *dmax__, int *kbf, int *kd, int *ld, int *
1052         iters)
1053 {
1054     /* System generated locals */
1055     int i__1;
1056     double d__1, d__2, d__3, d__4, d__5;
1057
1058     /* Local variables */
1059     int i__;
1060
1061     /* Parameter adjustments */
1062     --go;
1063     --g;
1064     --xo;
1065     --ix;
1066     --x;
1067
1068     /* Function Body */
1069     if (*iters > 0) {
1070         luksan_mxvdif__(nf, &x[1], &xo[1], &xo[1]);
1071         luksan_mxvdif__(nf, &g[1], &go[1], &go[1]);
1072         *po = *r__ * *po;
1073         *p = *r__ * *p;
1074     } else {
1075         *f = *fo;
1076         *p = *po;
1077         luksan_mxvsav__(nf, &x[1], &xo[1]);
1078         luksan_mxvsav__(nf, &g[1], &go[1]);
1079         *ld = *kd;
1080     }
1081     *dmax__ = 0.;
1082     i__1 = *nf;
1083     for (i__ = 1; i__ <= i__1; ++i__) {
1084         if (*kbf > 0) {
1085             if (ix[i__] < 0) {
1086                 xo[i__] = 0.;
1087                 go[i__] = 0.;
1088                 goto L1;
1089             }
1090         }
1091 /* Computing MAX */
1092 /* Computing MAX */
1093         d__5 = (d__2 = x[i__], fabs(d__2));
1094         d__3 = *dmax__, d__4 = (d__1 = xo[i__], fabs(d__1)) / MAX2(d__5,1.);
1095         *dmax__ = MAX2(d__3,d__4);
1096 L1:
1097         ;
1098     }
1099     return;
1100 } /* luksan_pytrcd__ */
1101
1102 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1103 /* SUBROUTINE PYTRCG                ALL SYSTEMS                99/12/01
1104 * PURPOSE :
1105 *  GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. TEST VALUES
1106 *  GMAX AND UMAX ARE COMPUTED.
1107 *
1108 * PARAMETERS :
1109 *  II  NF DECLARED NUMBER OF VARIABLES.
1110 *  II  N  ACTUAL NUMBER OF VARIABLES.
1111 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
1112 *  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
1113 *  RI  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
1114 *  RI  GMAX  NORM OF THE TRANSFORMED GRADIENT.
1115 *  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
1116 *         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
1117 *  II  IOLD  INDEX OF THE REMOVED CONSTRAINT.
1118 *
1119 * SUBPROGRAMS USED :
1120 *  RF  MXVMAX  L-INFINITY NORM OF A VECTOR.
1121 */
1122 void luksan_pytrcg__(int *nf, int *n, int *ix,
1123         double *g, double *umax, double *gmax, int *kbf,
1124         int *iold)
1125 {
1126     /* System generated locals */
1127     int i__1;
1128     double d__1, d__2;
1129
1130     /* Local variables */
1131     int i__;
1132     double temp;
1133
1134     /* Parameter adjustments */
1135     --g;
1136     --ix;
1137
1138     /* Function Body */
1139     if (*kbf > 0) {
1140         *gmax = 0.;
1141         *umax = 0.;
1142         *iold = 0;
1143         i__1 = *nf;
1144         for (i__ = 1; i__ <= i__1; ++i__) {
1145             temp = g[i__];
1146             if (ix[i__] >= 0) {
1147 /* Computing MAX */
1148                 d__1 = *gmax, d__2 = fabs(temp);
1149                 *gmax = MAX2(d__1,d__2);
1150             } else if (ix[i__] <= -5) {
1151             } else if ((ix[i__] == -1 || ix[i__] == -3) && *umax + temp >= 0.)
1152                      {
1153             } else if ((ix[i__] == -2 || ix[i__] == -4) && *umax - temp >= 0.)
1154                      {
1155             } else {
1156                 *iold = i__;
1157                 *umax = fabs(temp);
1158             }
1159 /* L1: */
1160         }
1161     } else {
1162         *umax = 0.;
1163         *gmax = luksan_mxvmax__(nf, &g[1]);
1164     }
1165     *n = *nf;
1166     return;
1167 } /* luksan_pytrcg__ */
1168
1169 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1170 /* SUBROUTINE PYTRCS             ALL SYSTEMS                   98/12/01
1171 * PURPOSE :
1172 * SCALED AND REDUCED DIRECTION VECTOR IS BACK TRANSFORMED. VECTORS
1173 * X,G AND VALUES F,P ARE SAVED.
1174 *
1175 * PARAMETERS :
1176 *  II  NF DECLARED NUMBER OF VARIABLES.
1177 *  RI  X(NF)  VECTOR OF VARIABLES.
1178 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
1179 *  RO  XO(NF)  SAVED VECTOR OF VARIABLES.
1180 *  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
1181 *  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
1182 *  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
1183 *  RO  GO(NF)  SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
1184 *  RO  S(NF)  DIRECTION VECTOR.
1185 *  RO  RO  SAVED VALUE OF THE STEPSIZE PARAMETER.
1186 *  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
1187 *  RU  FO  SAVED VALUE OF THE OBJECTIVE FUNCTION.
1188 *  RI  F  VALUE OF THE OBJECTIVE FUNCTION.
1189 *  RO  PO  SAVED VALUE OF THE DIRECTIONAL DERIVATIVE.
1190 *  RI  P  VALUE OF THE DIRECTIONAL DERIVATIVE.
1191 *  RO  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
1192 *  RI  ETA9  MAXIMUM FOR REAL NUMBERS.
1193 *  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
1194 *         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
1195 *
1196 * SUBPROGRAMS USED :
1197 *  S   MXVCOP  COPYING OF A VECTOR.
1198 */
1199 void luksan_pytrcs__(int *nf, double *x, int *ix,
1200         double *xo, double *xl, double *xu, double *g,
1201         double *go, double *s, double *ro, double *fp,
1202         double *fo, double *f, double *po, double *p,
1203         double *rmax, double *eta9, int *kbf)
1204 {
1205     /* System generated locals */
1206     int i__1;
1207     double d__1, d__2;
1208
1209     /* Local variables */
1210     int i__;
1211
1212     /* Parameter adjustments */
1213     --s;
1214     --go;
1215     --g;
1216     --xu;
1217     --xl;
1218     --xo;
1219     --ix;
1220     --x;
1221
1222     /* Function Body */
1223     *fp = *fo;
1224     *ro = 0.;
1225     *fo = *f;
1226     *po = *p;
1227     luksan_mxvcop__(nf, &x[1], &xo[1]);
1228     luksan_mxvcop__(nf, &g[1], &go[1]);
1229     if (*kbf > 0) {
1230         i__1 = *nf;
1231         for (i__ = 1; i__ <= i__1; ++i__) {
1232             if (ix[i__] < 0) {
1233                 s[i__] = 0.;
1234             } else {
1235                 if (ix[i__] == 1 || ix[i__] >= 3) {
1236                     if (s[i__] < -1. / *eta9) {
1237 /* Computing MIN */
1238                         d__1 = *rmax, d__2 = (xl[i__] - x[i__]) / s[i__];
1239                         *rmax = MIN2(d__1,d__2);
1240                     }
1241                 }
1242                 if (ix[i__] == 2 || ix[i__] >= 3) {
1243                     if (s[i__] > 1. / *eta9) {
1244 /* Computing MIN */
1245                         d__1 = *rmax, d__2 = (xu[i__] - x[i__]) / s[i__];
1246                         *rmax = MIN2(d__1,d__2);
1247                     }
1248                 }
1249             }
1250 /* L1: */
1251         }
1252     }
1253     return;
1254 } /* luksan_pytrcs__ */
1255