--- /dev/null
+#include <stdlib.h>
+#include <math.h>
+#include "luksan.h"
+
+#define TRUE_ 1
+#define FALSE_ 0
+
+/* Common Block Declarations */
+
+struct {
+ int nres, ndec, nin, nit, nfv, nfg, nfh;
+} stat_;
+
+#define stat_1 stat_
+
+/* *********************************************************************** */
+/* SUBROUTINE PLIP ALL SYSTEMS 01/09/22 */
+/* PURPOSE : */
+/* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT */
+/* USE THE SHIFTED LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE */
+/* PRODUCT FORM UPDATES. */
+
+/* PARAMETERS : */
+/* II NF NUMBER OF VARIABLES. */
+/* II NB CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. */
+/* NB>0-SIMPLE BOUNDS ACCEPTED. */
+/* RI X(NF) VECTOR OF VARIABLES. */
+/* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE */
+/* X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I). */
+/* IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND */
+/* XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. */
+/* RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. */
+/* RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. */
+/* RA GF(NF) GRADIENT OF THE OBJECTIVE FUNCTION. */
+/* RO S(NF) DIRECTION VECTOR. */
+/* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. */
+/* RI GO(NF) GRADIENTS DIFFERENCE. */
+/* RA SO(NF) AUXILIARY VECTOR. */
+/* RA XM(NF*MF) AUXILIARY VECTOR. */
+/* RA XR(MF) AUXILIARY VECTOR. */
+/* RA GR(MF) AUXILIARY VECTOR. */
+/* RI XMAX MAXIMUM STEPSIZE. */
+/* RI TOLX TOLERANCE FOR CHANGE OF VARIABLES. */
+/* RI TOLF TOLERANCE FOR CHANGE OF FUNCTION VALUES. */
+/* RI TOLB TOLERANCE FOR THE FUNCTION VALUE. */
+/* RI TOLG TOLERANCE FOR THE GRADIENT NORM. */
+/* RI FMIN ESTIMATION OF THE MINIMUM FUNCTION VALUE. */
+/* RO GMAX MAXIMUM PARTIAL DERIVATIVE. */
+/* RO F VALUE OF THE OBJECTIVE FUNCTION. */
+/* II MIT MAXIMUM NUMBER OF ITERATIONS. */
+/* II MFV MAXIMUM NUMBER OF FUNCTION EVALUATIONS. */
+/* II IEST ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. */
+/* IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN. */
+/* II MET METHOD USED. MET=1-RANK-ONE METHOD. MET=2-RANK-TWO */
+/* METHOD. */
+/* II MF NUMBER OF LIMITED MEMORY STEPS. */
+/* II IPRNT PRINT SPECIFICATION. IPRNT=0-NO PRINT. */
+/* ABS(IPRNT)=1-PRINT OF FINAL RESULTS. */
+/* ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS. */
+/* IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL */
+/* RESULTS. */
+/* IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. */
+/* ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN */
+/* MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
+/* ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN */
+/* MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
+/* ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. */
+/* ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. */
+/* ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, */
+/* BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. */
+/* ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. */
+/* ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. */
+
+/* VARIABLES IN COMMON /STAT/ (STATISTICS) : */
+/* IO NRES NUMBER OF RESTARTS. */
+/* IO NDEC NUMBER OF MATRIX DECOMPOSITION. */
+/* IO NIN NUMBER OF INNER ITERATIONS. */
+/* IO NIT NUMBER OF ITERATIONS. */
+/* IO NFV NUMBER OF FUNCTION EVALUATIONS. */
+/* IO NFG NUMBER OF GRADIENT EVALUATIONS. */
+/* IO NFH NUMBER OF HESSIAN EVALUATIONS. */
+
+/* SUBPROGRAMS USED : */
+/* S PCBS04 ELIMINATION OF BOX CONSTRAINT VIOLATIONS. */
+/* S PS1L01 STEPSIZE SELECTION USING LINE SEARCH. */
+/* S PULSP3 SHIFTED VARIABLE METRIC UPDATE. */
+/* S PULVP3 SHIFTED LIMITED-MEMORY VARIABLE METRIC UPDATE. */
+/* S PYADC0 ADDITION OF A BOX CONSTRAINT. */
+/* S PYFUT1 TEST ON TERMINATION. */
+/* S PYRMC0 DELETION OF A BOX CONSTRAINT. */
+/* S PYTRCD COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC */
+/* UPDATE. */
+/* S PYTRCG COMPUTATION OF THE PROJECTED GRADIENT. */
+/* S PYTRCS COMPUTATION OF THE PROJECTED DIRECTION VECTOR. */
+/* S MXDRMM MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR */
+/* MATRIX A BY A VECTOR X. */
+/* S MXDCMD MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR */
+/* MATRIX A BY A VECTOR X AND ADDITION OF THE SCALED VECTOR */
+/* ALF*Y. */
+/* S MXUCOP COPYING OF A VECTOR. */
+/* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. */
+/* RF MXUDOT DOT PRODUCT OF TWO VECTORS. */
+/* S MXUNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. */
+/* S MXUZER VECTOR ELEMENTS CORRESPONDING TO ACTIVE BOUNDS ARE SET */
+/* TO ZERO. */
+/* S MXVCOP COPYING OF A VECTOR. */
+
+/* EXTERNAL SUBROUTINES : */
+/* SE OBJ COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. */
+/* CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER */
+/* OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS */
+/* THE VALUE OF THE OBJECTIVE FUNCTION. */
+/* SE DOBJ COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. */
+/* CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER */
+/* OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF) */
+/* IS THE GRADIENT OF THE OBJECTIVE FUNCTION. */
+
+/* METHOD : */
+/* HYBRID METHOD WITH SPARSE MARWIL UPDATES FOR SPARSE LEAST SQUARES */
+/* PROBLEMS. */
+
+static void plip_(int *nf, int *nb, double *x, int *
+ ix, double *xl, double *xu, double *gf, double *s,
+ double *xo, double *go, double *so, double *xm,
+ double *xr, double *gr, double *xmax, double *tolx,
+ double *tolf, double *tolb, double *tolg, double *
+ fmin, double *gmax, double *f, int *mit, int *mfv,
+ int *iest, int *met, int *mf, int *iprnt, int *
+ iterm)
+{
+ /* System generated locals */
+ int i__1;
+ double d__1, d__2;
+
+ /* Builtin functions */
+
+ /* Local variables */
+ int i__, n;
+ double p, r__;
+ int kd, ld;
+ double fo, fp;
+ int nn;
+ double po, pp, ro, rp;
+ int kbf, mec, mfg;
+ extern static void obj_(int *, double *, double *);
+ double par;
+ int mes, kit;
+ double alf1, alf2, eta0, eta9, par1, par2;
+ int mes1, mes2, mes3, met3;
+ double eps8, eps9;
+ extern static void dobj_(int *, double *, double *);
+ int meta, mred, nred, iold;
+ double fmax, dmax__;
+ extern static void luksan_ps1l01__(double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, int *, int *, int *,
+ int *, int *, int *, int *, int *, int *,
+ int *, int *, int *, int *);
+ int inew;
+ double told;
+ int ites;
+ double rmin, rmax, umax, tolp, tols;
+ int isys;
+ extern static void luksan_pcbs04__(int *, double *,
+ int *, double *, double *, double *, int *);
+ int ires1, ires2;
+ extern static void luksan_pyadc0__(int *, int *,
+ double *, int *, double *, double *, int *);
+ int iterd, iterh, mtesf, ntesf;
+ double gnorm;
+ extern static void luksan_pyrmc0__(int *, int *, int
+ *, double *, double *, double *, double *,
+ double *, int *, int *);
+ int iters, irest, inits, kters, maxst;
+ double snorm;
+ int mtesx, ntesx;
+ extern static void luksan_pulsp3__(int *, int *, int
+ *, double *, double *, double *, double *,
+ double *, double *, double *, int *, int *),
+ luksan_pyfut1__(int *, double *, double *, double
+ *, double *, double *, double *, double *,
+ double *, double *, int *, int *, int *,
+ int *, int *, int *, int *, int *, int *,
+ int *, int *, int *, int *, int *, int *,
+ int *, int *, int *), luksan_pulvp3__(int *,
+ int *, double *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, int *, int *, int *,
+ int *), luksan_mxdcmd__(int *, int *, double *,
+ double *, double *, double *, double *),
+ luksan_mxuneg__(int *, double *, double *, int *,
+ int *), luksan_mxdrmm__(int *, int *, double *,
+ double *, double *), luksan_pytrcd__(int *,
+ double *, int *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, int *, int *, int *,
+ int *), luksan_pytrcg__(int *, int *, int *,
+ double *, double *, double *, int *, int *),
+ luksan_mxudir__(int *, double *, double *, double
+ *, double *, int *, int *), luksan_mxucop__(int *,
+ double *, double *, int *, int *),
+ luksan_mxvcop__(int *, double *, double *),
+ luksan_pytrcs__(int *, double *, int *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, int *), luksan_mxuzer__(int *, double *,
+ int *, int *);
+ extern double mxudot_(int *, double *, double *, int *
+ , int *);
+
+
+/* INITIATION */
+
+ /* Parameter adjustments */
+ --gr;
+ --xr;
+ --xm;
+ --so;
+ --go;
+ --xo;
+ --s;
+ --gf;
+ --xu;
+ --xl;
+ --ix;
+ --x;
+
+ /* Function Body */
+ kbf = 0;
+ if (*nb > 0) {
+ kbf = 2;
+ }
+ stat_1.nres = 0;
+ stat_1.ndec = 0;
+ stat_1.nin = 0;
+ stat_1.nit = 0;
+ stat_1.nfv = 0;
+ stat_1.nfg = 0;
+ stat_1.nfh = 0;
+ isys = 0;
+ ites = 1;
+ mtesx = 2;
+ mtesf = 2;
+ inits = 2;
+ *iterm = 0;
+ iterd = 0;
+ iters = 2;
+ iterh = 0;
+ kters = 3;
+ irest = 0;
+ ires1 = 999;
+ ires2 = 0;
+ mred = 10;
+ meta = 1;
+ met3 = 4;
+ mec = 4;
+ mes = 4;
+ mes1 = 2;
+ mes2 = 2;
+ mes3 = 2;
+ eta0 = 1e-15;
+ eta9 = 1e120;
+ eps8 = 1.;
+ eps9 = 1e-8;
+ alf1 = 1e-10;
+ alf2 = 1e10;
+ rmax = eta9;
+ dmax__ = eta9;
+ fmax = 1e20;
+ if (*iest <= 0) {
+ *fmin = -1e60;
+ }
+ if (*iest > 0) {
+ *iest = 1;
+ }
+ if (*xmax <= 0.) {
+ *xmax = 1e16;
+ }
+ if (*tolx <= 0.) {
+ *tolx = 1e-16;
+ }
+ if (*tolf <= 0.) {
+ *tolf = 1e-14;
+ }
+ if (*tolg <= 0.) {
+ *tolg = 1e-6;
+ }
+ if (*tolb <= 0.) {
+ *tolb = *fmin + 1e-16;
+ }
+ told = 1e-4;
+ tols = 1e-4;
+ tolp = .9;
+ if (*met <= 0) {
+ *met = 2;
+ }
+ if (*mit <= 0) {
+ *mit = 9000;
+ }
+ if (*mfv <= 0) {
+ *mfv = 9000;
+ }
+ mfg = *mfv;
+ kd = 1;
+ ld = -1;
+ kit = -(ires1 * *nf + ires2);
+ fo = *fmin;
+
+/* INITIAL OPERATIONS WITH SIMPLE BOUNDS */
+
+ if (kbf > 0) {
+ i__1 = *nf;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
+ xu[i__] = xl[i__];
+ ix[i__] = 5;
+ } else if (ix[i__] == 5 || ix[i__] == 6) {
+ xl[i__] = x[i__];
+ xu[i__] = x[i__];
+ ix[i__] = 5;
+ }
+/* L2: */
+ }
+ luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
+ luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
+ }
+ if (*iterm != 0) {
+ goto L11190;
+ }
+ obj_(nf, &x[1], f);
+ ++stat_1.nfv;
+ dobj_(nf, &x[1], &gf[1]);
+ ++stat_1.nfg;
+L11120:
+ luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
+ luksan_pyfut1__(nf, f, &fo, &umax, gmax, &dmax__, tolx, tolf, tolb, tolg,
+ &kd, &stat_1.nit, &kit, mit, &stat_1.nfv, mfv, &stat_1.nfg, &mfg,
+ &ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
+ iters, iterm);
+ if (*iterm != 0) {
+ goto L11190;
+ }
+ if (kbf > 0 && rmax > 0.) {
+ luksan_pyrmc0__(nf, &n, &ix[1], &gf[1], &eps8, &umax, gmax, &rmax, &
+ iold, &irest);
+ }
+L11130:
+ if (irest > 0) {
+ nn = 0;
+ par = 1.;
+ ld = min(ld,1);
+ if (kit < stat_1.nit) {
+ ++stat_1.nres;
+ kit = stat_1.nit;
+ } else {
+ *iterm = -10;
+ if (iters < 0) {
+ *iterm = iters - 5;
+ }
+ }
+ }
+ if (*iterm != 0) {
+ goto L11190;
+ }
+
+/* DIRECTION DETERMINATION */
+
+ gnorm = sqrt(mxudot_(nf, &gf[1], &gf[1], &ix[1], &kbf));
+
+/* NEWTON LIKE STEP */
+
+ luksan_mxuneg__(nf, &gf[1], &s[1], &ix[1], &kbf);
+ luksan_mxdrmm__(nf, &nn, &xm[1], &s[1], &gr[1]);
+ luksan_mxdcmd__(nf, &nn, &xm[1], &gr[1], &par, &s[1], &s[1]);
+ luksan_mxuzer__(nf, &s[1], &ix[1], &kbf);
+ iterd = 1;
+ snorm = sqrt(mxudot_(nf, &s[1], &s[1], &ix[1], &kbf));
+
+/* TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH */
+
+ if (kd > 0) {
+ p = mxudot_(nf, &gf[1], &s[1], &ix[1], &kbf);
+ }
+ if (iterd < 0) {
+ *iterm = iterd;
+ } else {
+
+/* TEST ON DESCENT DIRECTION */
+
+ if (snorm <= 0.) {
+ irest = max(irest,1);
+ } else if (p + told * gnorm * snorm <= 0.) {
+ irest = 0;
+ } else {
+
+/* UNIFORM DESCENT CRITERION */
+
+ irest = max(irest,1);
+ }
+ if (irest == 0) {
+
+/* PREPARATION OF LINE SEARCH */
+
+ nred = 0;
+ rmin = alf1 * gnorm / snorm;
+/* Computing MIN */
+ d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
+ rmax = min(d__1,d__2);
+ }
+ }
+ if (*iterm != 0) {
+ goto L11190;
+ }
+ if (irest != 0) {
+ goto L11130;
+ }
+ luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
+ &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
+ if (rmax == 0.) {
+ goto L11175;
+ }
+L11170:
+ luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, fmin, &fmax, &rmin,
+ &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1.nit, &kit, &
+ nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
+ if (isys == 0) {
+ goto L11174;
+ }
+ luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
+ luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
+ obj_(nf, &x[1], f);
+ ++stat_1.nfv;
+ dobj_(nf, &x[1], &gf[1]);
+ ++stat_1.nfg;
+ p = mxudot_(nf, &gf[1], &s[1], &ix[1], &kbf);
+ goto L11170;
+L11174:
+ if (iters <= 0) {
+ r__ = 0.;
+ *f = fo;
+ p = po;
+ luksan_mxvcop__(nf, &xo[1], &x[1]);
+ luksan_mxvcop__(nf, &go[1], &gf[1]);
+ irest = max(irest,1);
+ ld = kd;
+ goto L11130;
+ }
+ luksan_mxuneg__(nf, &go[1], &s[1], &ix[1], &kbf);
+ luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
+ p, &po, &dmax__, &kbf, &kd, &ld, &iters);
+ luksan_mxucop__(nf, &gf[1], &so[1], &ix[1], &kbf);
+ if (nn < *mf) {
+ luksan_pulsp3__(nf, &nn, mf, &xm[1], &gr[1], &xo[1], &go[1], &r__, &
+ po, &par, &iterh, &met3);
+ } else {
+ luksan_pulvp3__(nf, &nn, &xm[1], &xr[1], &gr[1], &s[1], &so[1], &xo[1]
+ , &go[1], &r__, &po, &par, &iterh, &mec, &met3, met);
+ }
+L11175:
+ if (iterh != 0) {
+ irest = max(irest,1);
+ }
+ if (kbf > 0) {
+ luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
+ }
+ goto L11120;
+L11190:
+ return;
+} /* plip_ */
+
--- /dev/null
+#include <stdlib.h>
+#include <math.h>
+#include "luksan.h"
+
+#define TRUE_ 1
+#define FALSE_ 0
+
+/* Common Block Declarations */
+
+struct {
+ int nres, ndec, nin, nit, nfv, nfg, nfh;
+} stat_;
+
+#define stat_1 stat_
+
+/* Table of constant values */
+
+static double c_b7 = 0.;
+
+/* *********************************************************************** */
+/* SUBROUTINE PNET ALL SYSTEMS 01/09/22 */
+/* PURPOSE : */
+/* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT */
+/* USE THE LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
+/* RECURRENCES. */
+
+/* PARAMETERS : */
+/* II NF NUMBER OF VARIABLES. */
+/* II NB CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED. */
+/* NB>0-SIMPLE BOUNDS ACCEPTED. */
+/* RI X(NF) VECTOR OF VARIABLES. */
+/* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE */
+/* X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I). */
+/* IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND */
+/* XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. */
+/* RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. */
+/* RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. */
+/* RO GF(NF) GRADIENT OF THE OBJECTIVE FUNCTION. */
+/* RA GN(NF) OLD GRADIENT OF THE OBJECTIVE FUNCTION. */
+/* RO S(NF) DIRECTION VECTOR. */
+/* RA XO(NF) ARRAY CONTAINING INCREMENTS OF VARIABLES. */
+/* RA GO(NF) ARRAY CONTAINING INCREMENTS OF GRADIENTS. */
+/* RA XS(NF) AUXILIARY VECTOR. */
+/* RA GS(NF) AUXILIARY VECTOR. */
+/* RA XM(NF*MF) ARRAY CONTAINING INCREMENTS OF VARIABLES. */
+/* RA GM(NF*MF) ARRAY CONTAINING INCREMENTS OF GRADIENTS. */
+/* RA U1(MF) AUXILIARY VECTOR. */
+/* RA U2(MF) AUXILIARY VECTOR. */
+/* RI XMAX MAXIMUM STEPSIZE. */
+/* RI TOLX TOLERANCE FOR CHANGE OF VARIABLES. */
+/* RI TOLF TOLERANCE FOR CHANGE OF FUNCTION VALUES. */
+/* RI TOLB TOLERANCE FOR THE FUNCTION VALUE. */
+/* RI TOLG TOLERANCE FOR THE GRADIENT NORM. */
+/* RI FMIN ESTIMATION OF THE MINIMUM FUNCTION VALUE. */
+/* RO GMAX MAXIMUM PARTIAL DERIVATIVE. */
+/* RO F VALUE OF THE OBJECTIVE FUNCTION. */
+/* II MIT MAXIMUM NUMBER OF ITERATIONS. */
+/* II MFV MAXIMUM NUMBER OF FUNCTION EVALUATIONS. */
+/* II MFG MAXIMUM NUMBER OF GRADIENT EVALUATIONS. */
+/* II IEST ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. */
+/* IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN. */
+/* II MOS1 CHOICE OF RESTARTS AFTER A CONSTRAINT CHANGE. */
+/* MOS1=1-RESTARTS ARE SUPPRESSED. MOS1=2-RESTARTS WITH */
+/* STEEPEST DESCENT DIRECTIONS ARE USED. */
+/* II MOS1 CHOICE OF DIRECTION VECTORS AFTER RESTARTS. MOS1=1-THE */
+/* NEWTON DIRECTIONS ARE USED. MOS1=2-THE STEEPEST DESCENT */
+/* DIRECTIONS ARE USED. */
+/* II MOS2 CHOICE OF PRECONDITIONING STRATEGY. MOS2=1-PRECONDITIONING */
+/* IS NOT USED. MOS2=2-PRECONDITIONING BY THE LIMITED MEMORY */
+/* BFGS METHOD IS USED. */
+/* II MF THE NUMBER OF LIMITED-MEMORY VARIABLE METRIC UPDATES */
+/* IN EACH ITERATION (THEY USE 2*MF STORED VECTORS). */
+/* II IPRNT PRINT SPECIFICATION. IPRNT=0-NO PRINT. */
+/* ABS(IPRNT)=1-PRINT OF FINAL RESULTS. */
+/* ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS. */
+/* IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL */
+/* RESULTS. */
+/* IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. */
+/* ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN */
+/* MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
+/* ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN */
+/* MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS. */
+/* ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB. */
+/* ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG. */
+/* ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, */
+/* BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. */
+/* ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. */
+/* ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. */
+
+/* VARIABLES IN COMMON /STAT/ (STATISTICS) : */
+/* IO NRES NUMBER OF RESTARTS. */
+/* IO NDEC NUMBER OF MATRIX DECOMPOSITION. */
+/* IO NIN NUMBER OF INNER ITERATIONS. */
+/* IO NIT NUMBER OF ITERATIONS. */
+/* IO NFV NUMBER OF FUNCTION EVALUATIONS. */
+/* IO NFG NUMBER OF GRADIENT EVALUATIONS. */
+/* IO NFH NUMBER OF HESSIAN EVALUATIONS. */
+
+/* SUBPROGRAMS USED : */
+/* S PCBS04 ELIMINATION OF BOX CONSTRAINT VIOLATIONS. */
+/* S PS1L01 STEPSIZE SELECTION USING LINE SEARCH. */
+/* S PYADC0 ADDITION OF A BOX CONSTRAINT. */
+/* S PYFUT1 TEST ON TERMINATION. */
+/* S PYRMC0 DELETION OF A BOX CONSTRAINT. */
+/* S PYTRCD COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC */
+/* UPDATE. */
+/* S PYTRCG COMPUTATION OF THE PROJECTED GRADIENT. */
+/* S PYTRCS COMPUTATION OF THE PROJECTED DIRECTION VECTOR. */
+/* S MXDRCB BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
+/* OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
+/* S MXDRCF FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION */
+/* OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE. */
+/* S MXDRSU SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B. */
+/* SHIFT OF ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN */
+/* THE LIMITED MEMORY BFGS METHOD. */
+/* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. */
+/* RF MXUDOT DOT PRODUCT OF TWO VECTORS. */
+/* S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. */
+/* S MXVCOP COPYING OF A VECTOR. */
+/* S MXVSCL SCALING OF A VECTOR. */
+/* S MXVSET INITIATINON OF A VECTOR. */
+/* S MXVDIF DIFFERENCE OF TWO VECTORS. */
+
+/* EXTERNAL SUBROUTINES : */
+/* SE OBJ COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION. */
+/* CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER */
+/* OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS */
+/* THE VALUE OF THE OBJECTIVE FUNCTION. */
+/* SE DOBJ COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION. */
+/* CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER */
+/* OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF) */
+/* IS THE GRADIENT OF THE OBJECTIVE FUNCTION. */
+
+/* METHOD : */
+/* LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG */
+/* RECURRENCES. */
+
+static void pnet_(int *nf, int *nb, double *x, int *
+ ix, double *xl, double *xu, double *gf, double *gn,
+ double *s, double *xo, double *go, double *xs,
+ double *gs, double *xm, double *gm, double *u1,
+ double *u2, double *xmax, double *tolx, double *tolf,
+ double *tolb, double *tolg, double *fmin, double *
+ gmax, double *f, int *mit, int *mfv, int *mfg,
+ int *iest, int *mos1, int *mos2, int *mf, int *
+ iprnt, int *iterm)
+{
+ /* System generated locals */
+ int i__1;
+ double d__1, d__2;
+
+ /* Builtin functions */
+
+ /* Local variables */
+ double a, b;
+ int i__, n;
+ double p, r__;
+ int kd, ld;
+ double fo, fp, po, pp, ro, rp;
+ int mx, kbf;
+ double alf;
+ extern static void obj_(int *, double *, double *);
+ double par;
+ int mes, kit;
+ double rho, eps;
+ int mmx;
+ double alf1, alf2, eta0, eta9, par1, par2;
+ int mes1, mes2, mes3;
+ double rho1, rho2, eps8, eps9;
+ extern static void dobj_(int *, double *, double *);
+ int mred, iold, nred;
+ double fmax, dmax__;
+ extern static void luksan_ps1l01__(double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, int *, int *, int *,
+ int *, int *, int *, int *, int *, int *,
+ int *, int *, int *, int *);
+ int inew;
+ double told;
+ int ites;
+ double rmin, rmax, umax, tolp, tols;
+ int isys;
+ extern static void luksan_pcbs04__(int *, double *,
+ int *, double *, double *, double *, int *);
+ int ires1, ires2;
+ extern static void luksan_pyadc0__(int *, int *,
+ double *, int *, double *, double *, int *);
+ int iterd, mtesf, ntesf;
+ double gnorm;
+ extern static void luksan_pyrmc0__(int *, int *, int
+ *, double *, double *, double *, double *,
+ double *, int *, int *);
+ int iters, irest, inits, kters, maxst;
+ double snorm;
+ int mtesx, ntesx;
+ extern static void luksan_pyfut1__(int *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, double *, int *,
+ int *, int *, int *, int *, int *, int *,
+ int *, int *, int *, int *, int *, int *,
+ int *, int *, int *, int *, int *),
+ luksan_mxdrcb__(int *, int *, double *, double *,
+ double *, double *, double *, int *, int *),
+ luksan_mxdrcf__(int *, int *, double *, double *,
+ double *, double *, double *, int *, int *),
+ luksan_mxvdif__(int *, double *, double *, double
+ *), luksan_mxvneg__(int *, double *, double *),
+ luksan_pytrcd__(int *, double *, int *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, double *, int *,
+ int *, int *, int *), luksan_pytrcg__(int *,
+ int *, int *, double *, double *, double *,
+ int *, int *), luksan_mxudir__(int *, double *,
+ double *, double *, double *, int *, int *),
+ luksan_mxvcop__(int *, double *, double *),
+ luksan_mxvscl__(int *, double *, double *, double
+ *), luksan_mxdrsu__(int *, int *, double *,
+ double *, double *), luksan_pytrcs__(int *,
+ double *, int *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, double *,
+ double *, double *, double *, int *),
+ luksan_mxvset__(int *, double *, double *);
+ extern double mxudot_(int *, double *, double *, int *
+ , int *);
+
+
+/* INITIATION */
+
+ /* Parameter adjustments */
+ --u2;
+ --u1;
+ --gm;
+ --xm;
+ --gs;
+ --xs;
+ --go;
+ --xo;
+ --s;
+ --gn;
+ --gf;
+ --xu;
+ --xl;
+ --ix;
+ --x;
+
+ /* Function Body */
+ kbf = 0;
+ if (*nb > 0) {
+ kbf = 2;
+ }
+ stat_1.nres = 0;
+ stat_1.ndec = 0;
+ stat_1.nin = 0;
+ stat_1.nit = 0;
+ stat_1.nfv = 0;
+ stat_1.nfg = 0;
+ stat_1.nfh = 0;
+ isys = 0;
+ ites = 1;
+ mtesx = 2;
+ mtesf = 2;
+ inits = 2;
+ *iterm = 0;
+ iterd = 0;
+ iters = 2;
+ kters = 3;
+ irest = 0;
+ ires1 = 999;
+ ires2 = 0;
+ mred = 10;
+ mes = 4;
+ mes1 = 2;
+ mes2 = 2;
+ mes3 = 2;
+ eps = .8;
+ eta0 = 1e-15;
+ eta9 = 1e120;
+ eps8 = 1.;
+ eps9 = 1e-8;
+ alf1 = 1e-10;
+ alf2 = 1e10;
+ rmax = eta9;
+ dmax__ = eta9;
+ fmax = 1e20;
+ if (*iest <= 0) {
+ *fmin = -1e60;
+ }
+ if (*iest > 0) {
+ *iest = 1;
+ }
+ if (*xmax <= 0.) {
+ *xmax = 1e16;
+ }
+ if (*tolx <= 0.) {
+ *tolx = 1e-16;
+ }
+ if (*tolf <= 0.) {
+ *tolf = 1e-14;
+ }
+ if (*tolg <= 0.) {
+ *tolg = 1e-6;
+ }
+ if (*tolb <= 0.) {
+ *tolb = *fmin + 1e-16;
+ }
+ told = 1e-4;
+ tols = 1e-4;
+ tolp = .9;
+ if (*mit <= 0) {
+ *mit = 5000;
+ }
+ if (*mfv <= 0) {
+ *mfv = 5000;
+ }
+ if (*mfg <= 0) {
+ *mfg = 30000;
+ }
+ if (*mos1 <= 0) {
+ *mos1 = 1;
+ }
+ if (*mos2 <= 0) {
+ *mos2 = 1;
+ }
+ kd = 1;
+ ld = -1;
+ kit = -(ires1 * *nf + ires2);
+ fo = *fmin;
+
+/* INITIAL OPERATIONS WITH SIMPLE BOUNDS */
+
+ if (kbf > 0) {
+ i__1 = *nf;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if ((ix[i__] == 3 || ix[i__] == 4) && xu[i__] <= xl[i__]) {
+ xu[i__] = xl[i__];
+ ix[i__] = 5;
+ } else if (ix[i__] == 5 || ix[i__] == 6) {
+ xl[i__] = x[i__];
+ xu[i__] = x[i__];
+ ix[i__] = 5;
+ }
+/* L2: */
+ }
+ luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
+ luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
+ }
+ obj_(nf, &x[1], f);
+ ++stat_1.nfv;
+ dobj_(nf, &x[1], &gf[1]);
+ ++stat_1.nfg;
+ ld = kd;
+L11020:
+ luksan_pytrcg__(nf, nf, &ix[1], &gf[1], &umax, gmax, &kbf, &iold);
+ luksan_mxvcop__(nf, &gf[1], &gn[1]);
+ luksan_pyfut1__(nf, f, &fo, &umax, gmax, &dmax__, tolx, tolf, tolb, tolg,
+ &kd, &stat_1.nit, &kit, mit, &stat_1.nfv, mfv, &stat_1.nfg, mfg, &
+ ntesx, &mtesx, &ntesf, &mtesf, &ites, &ires1, &ires2, &irest, &
+ iters, iterm);
+ if (*iterm != 0) {
+ goto L11080;
+ }
+ if (kbf > 0) {
+ luksan_pyrmc0__(nf, &n, &ix[1], &gn[1], &eps8, &umax, gmax, &rmax, &
+ iold, &irest);
+ if (umax > eps8 * *gmax) {
+ irest = max(irest,1);
+ }
+ }
+ luksan_mxvcop__(nf, &x[1], &xo[1]);
+L11040:
+
+/* DIRECTION DETERMINATION */
+
+ if (irest != 0) {
+ if (kit < stat_1.nit) {
+ mx = 0;
+ ++stat_1.nres;
+ kit = stat_1.nit;
+ } else {
+ *iterm = -10;
+ if (iters < 0) {
+ *iterm = iters - 5;
+ }
+ goto L11080;
+ }
+ if (*mos1 > 1) {
+ luksan_mxvneg__(nf, &gn[1], &s[1]);
+ gnorm = sqrt(mxudot_(nf, &gn[1], &gn[1], &ix[1], &kbf));
+ snorm = gnorm;
+ goto L12560;
+ }
+ }
+ rho1 = mxudot_(nf, &gn[1], &gn[1], &ix[1], &kbf);
+ gnorm = sqrt(rho1);
+/* Computing MIN */
+ d__1 = eps, d__2 = sqrt(gnorm);
+ par = min(d__1,d__2);
+ if (par > .01) {
+/* Computing MIN */
+ d__1 = par, d__2 = 1. / (double) stat_1.nit;
+ par = min(d__1,d__2);
+ }
+ par *= par;
+
+/* CG INITIATION */
+
+ rho = rho1;
+ snorm = 0.;
+ luksan_mxvset__(nf, &c_b7, &s[1]);
+ luksan_mxvneg__(nf, &gn[1], &gs[1]);
+ luksan_mxvcop__(nf, &gs[1], &xs[1]);
+ if (*mos2 > 1) {
+ if (mx == 0) {
+ b = 0.;
+ } else {
+ b = mxudot_(nf, &xm[1], &gm[1], &ix[1], &kbf);
+ }
+ if (b > 0.) {
+ u1[1] = 1. / b;
+ luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
+ ix[1], &kbf);
+ a = mxudot_(nf, &gm[1], &gm[1], &ix[1], &kbf);
+ if (a > 0.) {
+ d__1 = b / a;
+ luksan_mxvscl__(nf, &d__1, &xs[1], &xs[1]);
+ }
+ luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &xs[1], &
+ ix[1], &kbf);
+ }
+ }
+ rho = mxudot_(nf, &gs[1], &xs[1], &ix[1], &kbf);
+/* SIG=RHO */
+ mmx = *nf + 3;
+ nred = 0;
+L12520:
+ ++nred;
+ if (nred > mmx) {
+ goto L12550;
+ }
+ fo = *f;
+ pp = sqrt(eta0 / mxudot_(nf, &xs[1], &xs[1], &ix[1], &kbf));
+ ld = 0;
+ luksan_mxudir__(nf, &pp, &xs[1], &xo[1], &x[1], &ix[1], &kbf);
+ dobj_(nf, &x[1], &gf[1]);
+ ++stat_1.nfg;
+ ld = kd;
+ luksan_mxvdif__(nf, &gf[1], &gn[1], &go[1]);
+ *f = fo;
+ d__1 = 1. / pp;
+ luksan_mxvscl__(nf, &d__1, &go[1], &go[1]);
+ alf = mxudot_(nf, &xs[1], &go[1], &ix[1], &kbf);
+ if (alf <= 1. / eta9) {
+/* IF (ALF.LE.1.0D-8*SIG) THEN */
+
+/* CG FAILS (THE MATRIX IS NOT POSITIVE DEFINITE) */
+
+ if (nred == 1) {
+ luksan_mxvneg__(nf, &gn[1], &s[1]);
+ snorm = gnorm;
+ }
+ iterd = 0;
+ goto L12560;
+ } else {
+ iterd = 2;
+ }
+
+/* CG STEP */
+
+ alf = rho / alf;
+ luksan_mxudir__(nf, &alf, &xs[1], &s[1], &s[1], &ix[1], &kbf);
+ d__1 = -alf;
+ luksan_mxudir__(nf, &d__1, &go[1], &gs[1], &gs[1], &ix[1], &kbf);
+ rho2 = mxudot_(nf, &gs[1], &gs[1], &ix[1], &kbf);
+ snorm = sqrt(mxudot_(nf, &s[1], &s[1], &ix[1], &kbf));
+ if (rho2 <= par * rho1) {
+ goto L12560;
+ }
+ if (nred >= mmx) {
+ goto L12550;
+ }
+ if (*mos2 > 1) {
+ if (b > 0.) {
+ luksan_mxvcop__(nf, &gs[1], &go[1]);
+ luksan_mxdrcb__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
+ ix[1], &kbf);
+ if (a > 0.) {
+ d__1 = b / a;
+ luksan_mxvscl__(nf, &d__1, &go[1], &go[1]);
+ }
+ luksan_mxdrcf__(nf, &mx, &xm[1], &gm[1], &u1[1], &u2[1], &go[1], &
+ ix[1], &kbf);
+ rho2 = mxudot_(nf, &gs[1], &go[1], &ix[1], &kbf);
+ alf = rho2 / rho;
+ luksan_mxudir__(nf, &alf, &xs[1], &go[1], &xs[1], &ix[1], &kbf);
+ } else {
+ alf = rho2 / rho;
+ luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
+ }
+ } else {
+ alf = rho2 / rho;
+ luksan_mxudir__(nf, &alf, &xs[1], &gs[1], &xs[1], &ix[1], &kbf);
+ }
+ rho = rho2;
+/* SIG=RHO2+ALF*ALF*SIG */
+ goto L12520;
+L12550:
+
+/* AN INEXACT SOLUTION IS OBTAINED */
+
+L12560:
+
+/* ------------------------------ */
+/* END OF DIRECTION DETERMINATION */
+/* ------------------------------ */
+
+ luksan_mxvcop__(nf, &xo[1], &x[1]);
+ luksan_mxvcop__(nf, &gn[1], &gf[1]);
+ if (kd > 0) {
+ p = mxudot_(nf, &gn[1], &s[1], &ix[1], &kbf);
+ }
+ if (iterd < 0) {
+ *iterm = iterd;
+ } else {
+
+/* TEST ON DESCENT DIRECTION */
+
+ if (snorm <= 0.) {
+ irest = max(irest,1);
+ } else if (p + told * gnorm * snorm <= 0.) {
+ irest = 0;
+ } else {
+
+/* UNIFORM DESCENT CRITERION */
+
+ irest = max(irest,1);
+ }
+ if (irest == 0) {
+
+/* PREPARATION OF LINE SEARCH */
+
+ nred = 0;
+ rmin = alf1 * gnorm / snorm;
+/* Computing MIN */
+ d__1 = alf2 * gnorm / snorm, d__2 = *xmax / snorm;
+ rmax = min(d__1,d__2);
+ }
+ }
+ ld = kd;
+ if (*iterm != 0) {
+ goto L11080;
+ }
+ if (irest != 0) {
+ goto L11040;
+ }
+ luksan_pytrcs__(nf, &x[1], &ix[1], &xo[1], &xl[1], &xu[1], &gf[1], &go[1],
+ &s[1], &ro, &fp, &fo, f, &po, &p, &rmax, &eta9, &kbf);
+ if (rmax == 0.) {
+ goto L11075;
+ }
+L11060:
+ luksan_ps1l01__(&r__, &rp, f, &fo, &fp, &p, &po, &pp, fmin, &fmax, &rmin,
+ &rmax, &tols, &tolp, &par1, &par2, &kd, &ld, &stat_1.nit, &kit, &
+ nred, &mred, &maxst, iest, &inits, &iters, &kters, &mes, &isys);
+ if (isys == 0) {
+ goto L11064;
+ }
+ luksan_mxudir__(nf, &r__, &s[1], &xo[1], &x[1], &ix[1], &kbf);
+ luksan_pcbs04__(nf, &x[1], &ix[1], &xl[1], &xu[1], &eps9, &kbf);
+ obj_(nf, &x[1], f);
+ ++stat_1.nfv;
+ dobj_(nf, &x[1], &gf[1]);
+ ++stat_1.nfg;
+ ld = kd;
+ p = mxudot_(nf, &gf[1], &s[1], &ix[1], &kbf);
+ goto L11060;
+L11064:
+ if (iters <= 0) {
+ r__ = 0.;
+ *f = fo;
+ p = po;
+ luksan_mxvcop__(nf, &xo[1], &x[1]);
+ luksan_mxvcop__(nf, &go[1], &gf[1]);
+ irest = max(irest,1);
+ ld = kd;
+ goto L11040;
+ }
+ luksan_pytrcd__(nf, &x[1], &ix[1], &xo[1], &gf[1], &go[1], &r__, f, &fo, &
+ p, &po, &dmax__, &kbf, &kd, &ld, &iters);
+ if (*mos2 > 1) {
+/* Computing MIN */
+ i__1 = mx + 1;
+ mx = min(i__1,*mf);
+ luksan_mxdrsu__(nf, &mx, &xm[1], &gm[1], &u1[1]);
+ luksan_mxvcop__(nf, &xo[1], &xm[1]);
+ luksan_mxvcop__(nf, &go[1], &gm[1]);
+ }
+L11075:
+ if (kbf > 0) {
+ luksan_pyadc0__(nf, &n, &x[1], &ix[1], &xl[1], &xu[1], &inew);
+ if (inew > 0) {
+ irest = max(irest,1);
+ }
+ }
+ goto L11020;
+L11080:
+ return;
+} /* pnet_ */
+