From: stevenj Date: Mon, 3 Sep 2007 19:45:32 +0000 (-0400) Subject: added initial f2c conversions of plip and pnet X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=2ed5ad784725b0da8e26ad11826bd37f0eb217ca;p=nlopt.git added initial f2c conversions of plip and pnet darcs-hash:20070903194532-c8de0-725aba9fd7b6e3d4cba7c7f0f89f051a7fa2e60e.gz --- diff --git a/luksan/plip.c b/luksan/plip.c new file mode 100644 index 0000000..519e2d5 --- /dev/null +++ b/luksan/plip.c @@ -0,0 +1,473 @@ +#include +#include +#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_ */ + diff --git a/luksan/pnet.c b/luksan/pnet.c new file mode 100644 index 0000000..e37b57b --- /dev/null +++ b/luksan/pnet.c @@ -0,0 +1,611 @@ +#include +#include +#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_ */ +