luksan/pnet.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3 #include "luksan.h"
   4
   5 #define TRUE_ 1
   6 #define FALSE_ 0
   7
   8 /* Common Block Declarations */
   9
  10 struct {
  11     int nres, ndec, nin, nit, nfv, nfg, nfh;
  12 } stat_;
  13
  14 #define stat_1 stat_
  15
  16 /* Table of constant values */
  17
  18 static double c_b7 = 0.;
  19
  20 /* *********************************************************************** */
  21 /* SUBROUTINE PNET               ALL SYSTEMS                   01/09/22 */
  22 /* PURPOSE : */
  23 /* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT */
  24 /* USE THE LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
  25 /* RECURRENCES. */
  26
  27 /* PARAMETERS : */
  28 /*  II  NF  NUMBER OF VARIABLES. */
  29 /*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. */
  30 /*         NB>0-SIMPLE BOUNDS ACCEPTED. */
  31 /*  RI  X(NF)  VECTOR OF VARIABLES. */
  32 /*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE */
  33 /*         X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I). */
  34 /*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND */
  35 /*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. */
  36 /*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. */
  37 /*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. */
  38 /*  RO  GF(NF)  GRADIENT OF THE OBJECTIVE FUNCTION. */
  39 /*  RA  GN(NF)  OLD GRADIENT OF THE OBJECTIVE FUNCTION. */
  40 /*  RO  S(NF)  DIRECTION VECTOR. */
  41 /*  RA  XO(NF)  ARRAY CONTAINING INCREMENTS OF VARIABLES. */
  42 /*  RA  GO(NF)  ARRAY CONTAINING INCREMENTS OF GRADIENTS. */
  43 /*  RA  XS(NF)  AUXILIARY VECTOR. */
  44 /*  RA  GS(NF)  AUXILIARY VECTOR. */
  45 /*  RA  XM(NF*MF)  ARRAY CONTAINING INCREMENTS OF VARIABLES. */
  46 /*  RA  GM(NF*MF)  ARRAY CONTAINING INCREMENTS OF GRADIENTS. */
  47 /*  RA  U1(MF)  AUXILIARY VECTOR. */
  48 /*  RA  U2(MF)  AUXILIARY VECTOR. */
  49 /*  RI  XMAX  MAXIMUM STEPSIZE. */
  50 /*  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES. */
  51 /*  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES. */
  52 /*  RI  TOLB  TOLERANCE FOR THE FUNCTION VALUE. */
  53 /*  RI  TOLG  TOLERANCE FOR THE GRADIENT NORM. */
  54 /*  RI  FMIN  ESTIMATION OF THE MINIMUM FUNCTION VALUE. */
  55 /*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE. */
  56 /*  RO  F  VALUE OF THE OBJECTIVE FUNCTION. */
  57 /*  II  MIT  MAXIMUM NUMBER OF ITERATIONS. */
  58 /*  II  MFV  MAXIMUM NUMBER OF FUNCTION EVALUATIONS. */
  59 /*  II  MFG  MAXIMUM NUMBER OF GRADIENT EVALUATIONS. */
  60 /*  II  IEST  ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. */
  61 /*         IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN. */
  62 /*  II  MOS1  CHOICE OF RESTARTS AFTER A CONSTRAINT CHANGE. */
  63 /*         MOS1=1-RESTARTS ARE SUPPRESSED. MOS1=2-RESTARTS WITH */
  64 /*         STEEPEST DESCENT DIRECTIONS ARE USED. */
  65 /*  II  MOS1  CHOICE OF DIRECTION VECTORS AFTER RESTARTS. MOS1=1-THE */
  66 /*         NEWTON DIRECTIONS ARE USED. MOS1=2-THE STEEPEST DESCENT */
  67 /*         DIRECTIONS ARE USED. */
  68 /*  II  MOS2  CHOICE OF PRECONDITIONING STRATEGY. MOS2=1-PRECONDITIONING */
  69 /*         IS NOT USED. MOS2=2-PRECONDITIONING BY THE LIMITED MEMORY */
  70 /*         BFGS METHOD IS USED. */
  71 /*  II  MF  THE NUMBER OF LIMITED-MEMORY VARIABLE METRIC UPDATES */
  72 /*         IN EACH ITERATION (THEY USE 2*MF STORED VECTORS). */
  73 /*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT. */
  74 /*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS. */
  75 /*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS. */
  76 /*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL */
  77 /*         RESULTS. */
  78 /*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. */
  79 /*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN */
  80 /*                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
  81 /*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN */
  82 /*                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
  83 /*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. */
  84 /*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. */
  85 /*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, */
  86 /*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. */
  87 /*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. */
  88 /*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. */
  89
  90 /* VARIABLES IN COMMON /STAT/ (STATISTICS) : */
  91 /*  IO  NRES  NUMBER OF RESTARTS. */
  92 /*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION. */
  93 /*  IO  NIN  NUMBER OF INNER ITERATIONS. */
  94 /*  IO  NIT  NUMBER OF ITERATIONS. */
  95 /*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS. */
  96 /*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS. */
  97 /*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS. */
  98
  99 /* SUBPROGRAMS USED : */
 100 /*  S   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS. */
 101 /*  S   PS1L01  STEPSIZE SELECTION USING LINE SEARCH. */
 102 /*  S   PYADC0  ADDITION OF A BOX CONSTRAINT. */
 103 /*  S   PYFUT1  TEST ON TERMINATION. */
 104 /*  S   PYRMC0  DELETION OF A BOX CONSTRAINT. */
 105 /*  S   PYTRCD  COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC */
 106 /*         UPDATE. */
 107 /*  S   PYTRCG  COMPUTATION OF THE PROJECTED GRADIENT. */
 108 /*  S   PYTRCS  COMPUTATION OF THE PROJECTED DIRECTION VECTOR. */
 109 /*  S   MXDRCB BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
 110 /*         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
 111 /*  S   MXDRCF FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
 112 /*         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
 113 /*  S   MXDRSU SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B. */
 114 /*         SHIFT OF ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN */
 115 /*         THE LIMITED MEMORY BFGS METHOD. */
 116 /*  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR. */
 117 /*  RF  MXUDOT  DOT PRODUCT OF TWO VECTORS. */
 118 /*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN. */
 119 /*  S   MXVCOP  COPYING OF A VECTOR. */
 120 /*  S   MXVSCL  SCALING OF A VECTOR. */
 121 /*  S   MXVSET  INITIATINON OF A VECTOR. */
 122 /*  S   MXVDIF  DIFFERENCE OF TWO VECTORS. */
 123
 124 /* EXTERNAL SUBROUTINES : */
 125 /*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. */
 126 /*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER */
 127 /*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS */
 128 /*         THE VALUE OF THE OBJECTIVE FUNCTION. */
 129 /*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. */
 130 /*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER */
 131 /*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF) */
 132 /*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION. */
 133
 134 /* METHOD : */
 135 /* LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
 136 /* RECURRENCES. */
 137
 138 static void pnet_(int *nf, int *nb, double *x, int *
 139         ix, double *xl, double *xu, double *gf, double *gn,
 140         double *s, double *xo, double *go, double *xs,
 141         double *gs, double *xm, double *gm, double *u1,
 142         double *u2, double *xmax, double *tolx, double *tolf,
 143         double *tolb, double *tolg, double *fmin, double *
 144         gmax, double *f, int *mit, int *mfv, int *mfg,
 145         int *iest, int *mos1, int *mos2, int *mf, int *
 146         iprnt, int *iterm)
 147 {
 148     /* System generated locals */
 149     int i__1;
 150     double d__1, d__2;
 151
 152     /* Builtin functions */
 153
 154     /* Local variables */
 155     double a, b;
 156     int i__, n;
 157     double p, r__;
 158     int kd, ld;
 159     double fo, fp, po, pp, ro, rp;
 160     int mx, kbf;
 161     double alf;
 162     extern static void obj_(int *, double *, double *);
 163     double par;
 164     int mes, kit;
 165     double rho, eps;
 166     int mmx;
 167     double alf1, alf2, eta0, eta9, par1, par2;
 168     int mes1, mes2, mes3;
 169     double rho1, rho2, eps8, eps9;
 170     extern static void dobj_(int *, double *, double *);
 171     int mred, iold, nred;
 172     double fmax, dmax__;
 173     extern static void luksan_ps1l01__(double *, double *,
 174             double *, double *, double *, double *,
 175             double *, double *, double *, double *,
 176             double *, double *, double *, double *,
 177             double *, double *, int *, int *, int *,
 178             int *, int *, int *, int *, int *, int *,
 179             int *, int *, int *, int *);
 180     int inew;
 181     double told;
 182     int ites;
 183     double rmin, rmax, umax, tolp, tols;
 184     int isys;
 185     extern static void luksan_pcbs04__(int *, double *,
 186             int *, double *, double *, double *, int *);
 187     int ires1, ires2;
 188     extern static void luksan_pyadc0__(int *, int *,
 189             double *, int *, double *, double *, int *);
 190     int iterd, mtesf, ntesf;
 191     double gnorm;
 192     extern static void luksan_pyrmc0__(int *, int *, int
 193             *, double *, double *, double *, double *,
 194             double *, int *, int *);
 195     int iters, irest, inits, kters, maxst;
 196     double snorm;
 197     int mtesx, ntesx;
 198     extern static void luksan_pyfut1__(int *, double *,
 199             double *, double *, double *, double *,
 200             double *, double *, double *, double *, int *,
 201              int *, int *, int *, int *, int *, int *,
 202              int *, int *, int *, int *, int *, int *,
 203              int *, int *, int *, int *, int *),
 204             luksan_mxdrcb__(int *, int *, double *, double *,
 205             double *, double *, double *, int *, int *),
 206             luksan_mxdrcf__(int *, int *, double *, double *,
 207             double *, double *, double *, int *, int *),
 208             luksan_mxvdif__(int *, double *, double *, double
 209             *), luksan_mxvneg__(int *, double *, double *),
 210             luksan_pytrcd__(int *, double *, int *, double *,
 211             double *, double *, double *, double *,
 212             double *, double *, double *, double *, int *,
 213              int *, int *, int *), luksan_pytrcg__(int *,
 214             int *, int *, double *, double *, double *,
 215             int *, int *), luksan_mxudir__(int *, double *,
 216             double *, double *, double *, int *, int *),
 217             luksan_mxvcop__(int *, double *, double *),
 218             luksan_mxvscl__(int *, double *, double *, double
 219             *), luksan_mxdrsu__(int *, int *, double *,
 220             double *, double *), luksan_pytrcs__(int *,
 221             double *, int *, double *, double *, double *,
 222              double *, double *, double *, double *,
 223             double *, double *, double *, double *,
 224             double *, double *, double *, int *),
 225             luksan_mxvset__(int *, double *, double *);
 226     extern double mxudot_(int *, double *, double *, int *
 227             , int *);
 228
 229
 230 /*     INITIATION */
 231
 232     /* Parameter adjustments */
 233     --u2;
 234     --u1;
 235     --gm;
 236     --xm;
 237     --gs;
 238     --xs;
 239     --go;
 240     --xo;
 241     --s;
 242     --gn;
 243     --gf;
 244     --xu;
 245     --xl;
 246     --ix;
 247     --x;
 248
 249     /* Function Body */
 250     kbf = 0;
 251     if (*nb > 0) {
 252         kbf = 2;
 253     }
 254     stat_1.nres = 0;
 255     stat_1.ndec = 0;
 256     stat_1.nin = 0;
 257     stat_1.nit = 0;
 258     stat_1.nfv = 0;
 259     stat_1.nfg = 0;
 260     stat_1.nfh = 0;
 261     isys = 0;
 262     ites = 1;
 263     mtesx = 2;
 264     mtesf = 2;
 265     inits = 2;
 266     *iterm = 0;
 267     iterd = 0;
 268     iters = 2;
 269     kters = 3;
 270     irest = 0;
 271     ires1 = 999;
 272     ires2 = 0;
 273     mred = 10;
 274     mes = 4;
 275     mes1 = 2;
 276     mes2 = 2;
 277     mes3 = 2;
 278     eps = .8;
 279     eta0 = 1e-15;
 280     eta9 = 1e120;
 281     eps8 = 1.;
 282     eps9 = 1e-8;
 283     alf1 = 1e-10;
 284     alf2 = 1e10;
 285     rmax = eta9;
 286     dmax__ = eta9;
 287     fmax = 1e20;
 288     if (*iest <= 0) {
 289         *fmin = -1e60;
 290     }
 291     if (*iest > 0) {
 292         *iest = 1;
 293     }
 294     if (*xmax <= 0.) {
 295         *xmax = 1e16;
 296     }
 297     if (*tolx <= 0.) {
 298         *tolx = 1e-16;
 299     }
 300     if (*tolf <= 0.) {
 301         *tolf = 1e-14;
 302     }
 303     if (*tolg <= 0.) {
 304         *tolg = 1e-6;
 305     }
 306     if (*tolb <= 0.) {
 307         *tolb = *fmin + 1e-16;
 308     }
 309     told = 1e-4;
 310     tols = 1e-4;
 311     tolp = .9;
 312     if (*mit <= 0) {
 313         *mit = 5000;
 314     }
 315     if (*mfv <= 0) {
 316         *mfv = 5000;
 317     }
 318     if (*mfg <= 0) {
 319         *mfg = 30000;
 320     }
 321     if (*mos1 <= 0) {
 322         *mos1 = 1;
 323     }
 324     if (*mos2 <= 0) {
 325         *mos2 = 1;
 326     }
 327     kd = 1;
 328     ld = -1;
 329     kit = -(ires1 * *nf + ires2);
 330     fo = *fmin;
 331
 332 /*     INITIAL OPERATIONS WITH SIMPLE BOUNDS */
 333
 334     if (kbf > 0) {
 335         i__1 = *nf;
 336         for (i__ = 1; i__ <= i__1; ++i__) {
 337             if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
 338                 xu[i__] = xl[i__];
 339                 ix[i__] = 5;
 340             } else if (ix[i__] == 5 || ix[i__] == 6) {
 341                 xl[i__] = x[i__];
 342                 xu[i__] = x[i__];
 343                 ix[i__] = 5;
 344             }
 345 /* L2: */
 346         }
 347         luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
 348         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
 349     }
 350     obj_(nf, &x[1], f);
 351     ++stat_1.nfv;
 352     dobj_(nf, &x[1], &gf[1]);
 353     ++stat_1.nfg;
 354     ld = kd;
 355 L11020:
 356     luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
 357     luksan_mxvcop__(nf, &gf[1], &gn[1]);
 358     luksan_pyfut1__(nf, f, &fo, &umax, gmax, &dmax__, tolx, tolf, tolb, tolg,
 359             &kd, &stat_1.nit, &kit, mit, &stat_1.nfv, mfv, &stat_1.nfg, mfg, &
 360             ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
 361             iters, iterm);
 362     if (*iterm != 0) {
 363         goto L11080;
 364     }
 365     if (kbf > 0) {
 366         luksan_pyrmc0__(nf, &n, &ix[1], &gn[1], &eps8, &umax, gmax, &rmax, &
 367                 iold, &irest);
 368         if (umax > eps8 * *gmax) {
 369             irest = max(irest,1);
 370         }
 371     }
 372     luksan_mxvcop__(nf, &x[1], &xo[1]);
 373 L11040:
 374
 375 /*     DIRECTION DETERMINATION */
 376
 377     if (irest != 0) {
 378         if (kit < stat_1.nit) {
 379             mx = 0;
 380             ++stat_1.nres;
 381             kit = stat_1.nit;
 382         } else {
 383             *iterm = -10;
 384             if (iters < 0) {
 385                 *iterm = iters - 5;
 386             }
 387             goto L11080;
 388         }
 389         if (*mos1 > 1) {
 390             luksan_mxvneg__(nf, &gn[1], &s[1]);
 391             gnorm = sqrt(mxudot_(nf, &gn[1], &gn[1], &ix[1], &kbf));
 392             snorm = gnorm;
 393             goto L12560;
 394         }
 395     }
 396     rho1 = mxudot_(nf, &gn[1], &gn[1], &ix[1], &kbf);
 397     gnorm = sqrt(rho1);
 398 /* Computing MIN */
 399     d__1 = eps, d__2 = sqrt(gnorm);
 400     par = min(d__1,d__2);
 401     if (par > .01) {
 402 /* Computing MIN */
 403         d__1 = par, d__2 = 1. / (double) stat_1.nit;
 404         par = min(d__1,d__2);
 405     }
 406     par *= par;
 407
 408 /*     CG INITIATION */
 409
 410     rho = rho1;
 411     snorm = 0.;
 412     luksan_mxvset__(nf, &c_b7, &s[1]);
 413     luksan_mxvneg__(nf, &gn[1], &gs[1]);
 414     luksan_mxvcop__(nf, &gs[1], &xs[1]);
 415     if (*mos2 > 1) {
 416         if (mx == 0) {
 417             b = 0.;
 418         } else {
 419             b = mxudot_(nf, &xm[1], &gm[1], &ix[1], &kbf);
 420         }
 421         if (b > 0.) {
 422             u1[1] = 1. / b;
 423             luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
 424                     ix[1], &kbf);
 425             a = mxudot_(nf, &gm[1], &gm[1], &ix[1], &kbf);
 426             if (a > 0.) {
 427                 d__1 = b / a;
 428                 luksan_mxvscl__(nf, &d__1, &xs[1], &xs[1]);
 429             }
 430             luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
 431                     ix[1], &kbf);
 432         }
 433     }
 434     rho = mxudot_(nf, &gs[1], &xs[1], &ix[1], &kbf);
 435 /*      SIG=RHO */
 436     mmx = *nf + 3;
 437     nred = 0;
 438 L12520:
 439     ++nred;
 440     if (nred > mmx) {
 441         goto L12550;
 442     }
 443     fo = *f;
 444     pp = sqrt(eta0 / mxudot_(nf, &xs[1], &xs[1], &ix[1], &kbf));
 445     ld = 0;
 446     luksan_mxudir__(nf, &pp, &xs[1], &xo[1], &x[1], &ix[1], &kbf);
 447     dobj_(nf, &x[1], &gf[1]);
 448     ++stat_1.nfg;
 449     ld = kd;
 450     luksan_mxvdif__(nf, &gf[1], &gn[1], &go[1]);
 451     *f = fo;
 452     d__1 = 1. / pp;
 453     luksan_mxvscl__(nf, &d__1, &go[1], &go[1]);
 454     alf = mxudot_(nf, &xs[1], &go[1], &ix[1], &kbf);
 455     if (alf <= 1. / eta9) {
 456 /*      IF (ALF.LE.1.0D-8*SIG) THEN */
 457
 458 /*     CG FAILS (THE MATRIX IS NOT POSITIVE DEFINITE) */
 459
 460         if (nred == 1) {
 461             luksan_mxvneg__(nf, &gn[1], &s[1]);
 462             snorm = gnorm;
 463         }
 464         iterd = 0;
 465         goto L12560;
 466     } else {
 467         iterd = 2;
 468     }
 469
 470 /*     CG STEP */
 471
 472     alf = rho / alf;
 473     luksan_mxudir__(nf, &alf, &xs[1], &s[1], &s[1], &ix[1], &kbf);
 474     d__1 = -alf;
 475     luksan_mxudir__(nf, &d__1, &go[1], &gs[1], &gs[1], &ix[1], &kbf);
 476     rho2 = mxudot_(nf, &gs[1], &gs[1], &ix[1], &kbf);
 477     snorm = sqrt(mxudot_(nf, &s[1], &s[1], &ix[1], &kbf));
 478     if (rho2 <= par * rho1) {
 479         goto L12560;
 480     }
 481     if (nred >= mmx) {
 482         goto L12550;
 483     }
 484     if (*mos2 > 1) {
 485         if (b > 0.) {
 486             luksan_mxvcop__(nf, &gs[1], &go[1]);
 487             luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
 488                     ix[1], &kbf);
 489             if (a > 0.) {
 490                 d__1 = b / a;
 491                 luksan_mxvscl__(nf, &d__1, &go[1], &go[1]);
 492             }
 493             luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
 494                     ix[1], &kbf);
 495             rho2 = mxudot_(nf, &gs[1], &go[1], &ix[1], &kbf);
 496             alf = rho2 / rho;
 497             luksan_mxudir__(nf, &alf, &xs[1], &go[1], &xs[1], &ix[1], &kbf);
 498         } else {
 499             alf = rho2 / rho;
 500             luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
 501         }
 502     } else {
 503         alf = rho2 / rho;
 504         luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
 505     }
 506     rho = rho2;
 507 /*      SIG=RHO2+ALF*ALF*SIG */
 508     goto L12520;
 509 L12550:
 510
 511 /*     AN INEXACT SOLUTION IS OBTAINED */
 512
 513 L12560:
 514
 515 /*     ------------------------------ */
 516 /*     END OF DIRECTION DETERMINATION */
 517 /*     ------------------------------ */
 518
 519     luksan_mxvcop__(nf, &xo[1], &x[1]);
 520     luksan_mxvcop__(nf, &gn[1], &gf[1]);
 521     if (kd > 0) {
 522         p = mxudot_(nf, &gn[1], &s[1], &ix[1], &kbf);
 523     }
 524     if (iterd < 0) {
 525         *iterm = iterd;
 526     } else {
 527
 528 /*     TEST ON DESCENT DIRECTION */
 529
 530         if (snorm <= 0.) {
 531             irest = max(irest,1);
 532         } else if (p + told * gnorm * snorm <= 0.) {
 533             irest = 0;
 534         } else {
 535
 536 /*     UNIFORM DESCENT CRITERION */
 537
 538             irest = max(irest,1);
 539         }
 540         if (irest == 0) {
 541
 542 /*     PREPARATION OF LINE SEARCH */
 543
 544             nred = 0;
 545             rmin = alf1 * gnorm / snorm;
 546 /* Computing MIN */
 547             d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
 548             rmax = min(d__1,d__2);
 549         }
 550     }
 551     ld = kd;
 552     if (*iterm != 0) {
 553         goto L11080;
 554     }
 555     if (irest != 0) {
 556         goto L11040;
 557     }
 558     luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
 559              &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
 560     if (rmax == 0.) {
 561         goto L11075;
 562     }
 563 L11060:
 564     luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, fmin, &fmax, &rmin,
 565             &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1.nit, &kit, &
 566             nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
 567     if (isys == 0) {
 568         goto L11064;
 569     }
 570     luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
 571     luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
 572     obj_(nf, &x[1], f);
 573     ++stat_1.nfv;
 574     dobj_(nf, &x[1], &gf[1]);
 575     ++stat_1.nfg;
 576     ld = kd;
 577     p = mxudot_(nf, &gf[1], &s[1], &ix[1], &kbf);
 578     goto L11060;
 579 L11064:
 580     if (iters <= 0) {
 581         r__ = 0.;
 582         *f = fo;
 583         p = po;
 584         luksan_mxvcop__(nf, &xo[1], &x[1]);
 585         luksan_mxvcop__(nf, &go[1], &gf[1]);
 586         irest = max(irest,1);
 587         ld = kd;
 588         goto L11040;
 589     }
 590     luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
 591             p, &po, &dmax__, &kbf, &kd, &ld, &iters);
 592     if (*mos2 > 1) {
 593 /* Computing MIN */
 594         i__1 = mx + 1;
 595         mx = min(i__1,*mf);
 596         luksan_mxdrsu__(nf, &mx, &xm[1], &gm[1], &u1[1]);
 597         luksan_mxvcop__(nf, &xo[1], &xm[1]);
 598         luksan_mxvcop__(nf, &go[1], &gm[1]);
 599     }
 600 L11075:
 601     if (kbf > 0) {
 602         luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
 603         if (inew > 0) {
 604             irest = max(irest,1);
 605         }
 606     }
 607     goto L11020;
 608 L11080:
 609     return;
 610 } /* pnet_ */
 611