提交 2ed5ad78 编写于 作者: S stevenj

added initial f2c conversions of plip and pnet

darcs-hash:20070903194532-c8de0-725aba9fd7b6e3d4cba7c7f0f89f051a7fa2e60e.gz
上级 cab3a453
#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_ */
#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_ */
Markdown is supported
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册