chiark / gitweb /
added original .for files; this will make it easier to diff the changes if upstream...
[nlopt.git] / 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