8 /* Common Block Declarations */
11 int nres, ndec, nin, nit, nfv, nfg, nfh;
16 /* Table of constant values */
18 static double c_b7 = 0.;
20 /* *********************************************************************** */
21 /* SUBROUTINE PNET ALL SYSTEMS 01/09/22 */
23 /* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT */
24 /* USE THE LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
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 */
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. */
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. */
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 */
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. */
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. */
135 /* LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
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 *
148 /* System generated locals */
152 /* Builtin functions */
154 /* Local variables */
159 double fo, fp, po, pp, ro, rp;
162 extern static void obj_(int *, double *, double *);
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;
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 *);
183 double rmin, rmax, umax, tolp, tols;
185 extern static void luksan_pcbs04__(int *, double *,
186 int *, double *, double *, double *, int *);
188 extern static void luksan_pyadc0__(int *, int *,
189 double *, int *, double *, double *, int *);
190 int iterd, mtesf, ntesf;
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;
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 *
232 /* Parameter adjustments */
307 *tolb = *fmin + 1e-16;
329 kit = -(ires1 * *nf + ires2);
332 /* INITIAL OPERATIONS WITH SIMPLE BOUNDS */
336 for (i__ = 1; i__ <= i__1; ++i__) {
337 if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
340 } else if (ix[i__] == 5 || ix[i__] == 6) {
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);
352 dobj_(nf, &x[1], &gf[1]);
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, &
366 luksan_pyrmc0__(nf, &n, &ix[1], &gn[1], &eps8, &umax, gmax, &rmax, &
368 if (umax > eps8 * *gmax) {
369 irest = max(irest,1);
372 luksan_mxvcop__(nf, &x[1], &xo[1]);
375 /* DIRECTION DETERMINATION */
378 if (kit < stat_1.nit) {
390 luksan_mxvneg__(nf, &gn[1], &s[1]);
391 gnorm = sqrt(mxudot_(nf, &gn[1], &gn[1], &ix[1], &kbf));
396 rho1 = mxudot_(nf, &gn[1], &gn[1], &ix[1], &kbf);
399 d__1 = eps, d__2 = sqrt(gnorm);
400 par = min(d__1,d__2);
403 d__1 = par, d__2 = 1. / (double) stat_1.nit;
404 par = min(d__1,d__2);
412 luksan_mxvset__(nf, &c_b7, &s[1]);
413 luksan_mxvneg__(nf, &gn[1], &gs[1]);
414 luksan_mxvcop__(nf, &gs[1], &xs[1]);
419 b = mxudot_(nf, &xm[1], &gm[1], &ix[1], &kbf);
423 luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
425 a = mxudot_(nf, &gm[1], &gm[1], &ix[1], &kbf);
428 luksan_mxvscl__(nf, &d__1, &xs[1], &xs[1]);
430 luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
434 rho = mxudot_(nf, &gs[1], &xs[1], &ix[1], &kbf);
444 pp = sqrt(eta0 / mxudot_(nf, &xs[1], &xs[1], &ix[1], &kbf));
446 luksan_mxudir__(nf, &pp, &xs[1], &xo[1], &x[1], &ix[1], &kbf);
447 dobj_(nf, &x[1], &gf[1]);
450 luksan_mxvdif__(nf, &gf[1], &gn[1], &go[1]);
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 */
458 /* CG FAILS (THE MATRIX IS NOT POSITIVE DEFINITE) */
461 luksan_mxvneg__(nf, &gn[1], &s[1]);
473 luksan_mxudir__(nf, &alf, &xs[1], &s[1], &s[1], &ix[1], &kbf);
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) {
486 luksan_mxvcop__(nf, &gs[1], &go[1]);
487 luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
491 luksan_mxvscl__(nf, &d__1, &go[1], &go[1]);
493 luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
495 rho2 = mxudot_(nf, &gs[1], &go[1], &ix[1], &kbf);
497 luksan_mxudir__(nf, &alf, &xs[1], &go[1], &xs[1], &ix[1], &kbf);
500 luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
504 luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
507 /* SIG=RHO2+ALF*ALF*SIG */
511 /* AN INEXACT SOLUTION IS OBTAINED */
515 /* ------------------------------ */
516 /* END OF DIRECTION DETERMINATION */
517 /* ------------------------------ */
519 luksan_mxvcop__(nf, &xo[1], &x[1]);
520 luksan_mxvcop__(nf, &gn[1], &gf[1]);
522 p = mxudot_(nf, &gn[1], &s[1], &ix[1], &kbf);
528 /* TEST ON DESCENT DIRECTION */
531 irest = max(irest,1);
532 } else if (p + told * gnorm * snorm <= 0.) {
536 /* UNIFORM DESCENT CRITERION */
538 irest = max(irest,1);
542 /* PREPARATION OF LINE SEARCH */
545 rmin = alf1 * gnorm / snorm;
547 d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
548 rmax = min(d__1,d__2);
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);
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);
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);
574 dobj_(nf, &x[1], &gf[1]);
577 p = mxudot_(nf, &gf[1], &s[1], &ix[1], &kbf);
584 luksan_mxvcop__(nf, &xo[1], &x[1]);
585 luksan_mxvcop__(nf, &go[1], &gf[1]);
586 irest = max(irest,1);
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);
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]);
602 luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
604 irest = max(irest,1);