📄 lbfgs.c
字号:
/* opt/lbfgs.f -- translated by f2c (version 20050501).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"
#undef abs
#undef min
#undef max
#include <math.h>
#include <stdio.h>
#define abs(x) ((x) >= 0 ? (x) : -(x))
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))
/* Provide initialization function for global data argument. */
#define lb3_1 (*v3p_netlib_lbfgs_global_arg)
void
v3p_netlib_lbfgs_init(v3p_netlib_lbfgs_global_t* v3p_netlib_lbfgs_global_arg)
{
lb3_1.mp = 6;
lb3_1.lp = 6;
lb3_1.gtol = 0.9;
lb3_1.stpmin = 1e-20;
lb3_1.stpmax = 1e20;
lb3_1.stpinit = 1; /* line search default step length, added by awf */
}
/* Map persistent data. */
#define info (*v3p_netlib_lbfgs_global_arg).private_info
#define infoc (*v3p_netlib_lbfgs_global_arg).private_infoc
#define nfev (*v3p_netlib_lbfgs_global_arg).private_nfev
#define maxfev (*v3p_netlib_lbfgs_global_arg).private_maxfev
#define stp (*v3p_netlib_lbfgs_global_arg).private_stp
#define stp1 (*v3p_netlib_lbfgs_global_arg).private_stp1
#define beta (*v3p_netlib_lbfgs_global_arg).private_beta
#define ftol (*v3p_netlib_lbfgs_global_arg).private_ftol
#define gnorm (*v3p_netlib_lbfgs_global_arg).private_gnorm
#define xnorm (*v3p_netlib_lbfgs_global_arg).private_xnorm
#define inmc (*v3p_netlib_lbfgs_global_arg).private_inmc
#define iscn (*v3p_netlib_lbfgs_global_arg).private_iscn
#define iycn (*v3p_netlib_lbfgs_global_arg).private_iycn
#define iter (*v3p_netlib_lbfgs_global_arg).private_iter
#define nfun (*v3p_netlib_lbfgs_global_arg).private_nfun
#define ispt (*v3p_netlib_lbfgs_global_arg).private_ispt
#define iypt (*v3p_netlib_lbfgs_global_arg).private_iypt
#define bound (*v3p_netlib_lbfgs_global_arg).private_bound
#define point (*v3p_netlib_lbfgs_global_arg).private_point
#define finish (*v3p_netlib_lbfgs_global_arg).private_finish
#define dg (*v3p_netlib_lbfgs_global_arg).private_dg
#define fm (*v3p_netlib_lbfgs_global_arg).private_fm
#define fx (*v3p_netlib_lbfgs_global_arg).private_fx
#define fy (*v3p_netlib_lbfgs_global_arg).private_fy
#define dgm (*v3p_netlib_lbfgs_global_arg).private_dgm
#define dgx (*v3p_netlib_lbfgs_global_arg).private_dgx
#define dgy (*v3p_netlib_lbfgs_global_arg).private_dgy
#define fxm (*v3p_netlib_lbfgs_global_arg).private_fxm
#define fym (*v3p_netlib_lbfgs_global_arg).private_fym
#define stx (*v3p_netlib_lbfgs_global_arg).private_stx
#define sty (*v3p_netlib_lbfgs_global_arg).private_sty
#define dgxm (*v3p_netlib_lbfgs_global_arg).private_dgxm
#define dgym (*v3p_netlib_lbfgs_global_arg).private_dgym
#define finit (*v3p_netlib_lbfgs_global_arg).private_finit
#define width (*v3p_netlib_lbfgs_global_arg).private_width
#define stmin (*v3p_netlib_lbfgs_global_arg).private_stmin
#define stmax (*v3p_netlib_lbfgs_global_arg).private_stmax
#define stage1 (*v3p_netlib_lbfgs_global_arg).private_stage1
#define width1 (*v3p_netlib_lbfgs_global_arg).private_width1
#define ftest1 (*v3p_netlib_lbfgs_global_arg).private_ftest1
#define brackt (*v3p_netlib_lbfgs_global_arg).private_brackt
#define dginit (*v3p_netlib_lbfgs_global_arg).private_dginit
#define dgtest (*v3p_netlib_lbfgs_global_arg).private_dgtest
/* Table of constant values */
static integer c__1 = 1;
/* ---------------------------------------------------------------------- */
/* This file contains the LBFGS algorithm and supporting routines */
/* **************** */
/* LBFGS SUBROUTINE */
/* **************** */
/*< SUBROUTINE LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG) >*/
/* Subroutine */ int lbfgs_(integer *n, integer *m, doublereal *x, doublereal
*f, doublereal *g, logical *diagco, doublereal *diag, integer *iprint,
doublereal *eps, doublereal *xtol, doublereal *w, integer *iflag,
v3p_netlib_lbfgs_global_t* v3p_netlib_lbfgs_global_arg)
{
/* Initialized data */
static doublereal one = 1.; /* constant */
static doublereal zero = 0.; /* constant */
/* System generated locals */
integer i__1;
doublereal d__1;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
integer i__, cp;
doublereal sq, yr, ys=0, yy;
extern /* Subroutine */ int lb1_(
integer *iprint, integer *n, integer *m, doublereal *x, doublereal *f,
doublereal *g, v3p_netlib_lbfgs_global_t* v3p_netlib_lbfgs_global_arg);
extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
integer *);
extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
integer *, doublereal *, integer *);
extern /* Subroutine */ int mcsrch_(integer *, doublereal *, doublereal *,
doublereal *, doublereal *,
doublereal *, doublereal *,
v3p_netlib_lbfgs_global_t*);
integer npt=0;
/*< INTEGER N,M,IPRINT(2),IFLAG >*/
/*< DOUBLE PRECISION X(N),G(N),DIAG(N),W(N*(2*M+1)+2*M) >*/
/*< DOUBLE PRECISION F,EPS,XTOL >*/
/*< LOGICAL DIAGCO >*/
/* LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION */
/* JORGE NOCEDAL */
/* *** July 1990 *** */
/* This subroutine solves the unconstrained minimization problem */
/* min F(x), x= (x1,x2,...,xN), */
/* using the limited memory BFGS method. The routine is especially */
/* effective on problems involving a large number of variables. In */
/* a typical iteration of this method an approximation Hk to the */
/* inverse of the Hessian is obtained by applying M BFGS updates to */
/* a diagonal matrix Hk0, using information from the previous M steps. */
/* The user specifies the number M, which determines the amount of */
/* storage required by the routine. The user may also provide the */
/* diagonal matrices Hk0 if not satisfied with the default choice. */
/* The algorithm is described in "On the limited memory BFGS method */
/* for large scale optimization", by D. Liu and J. Nocedal, */
/* Mathematical Programming B 45 (1989) 503-528. */
/* The user is required to calculate the function value F and its */
/* gradient G. In order to allow the user complete control over */
/* these computations, reverse communication is used. The routine */
/* must be called repeatedly under the control of the parameter */
/* IFLAG. */
/* The steplength is determined at each iteration by means of the */
/* line search routine MCVSRCH, which is a slight modification of */
/* the routine CSRCH written by More' and Thuente. */
/* The calling statement is */
/* CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG) */
/* where */
/* N is an INTEGER variable that must be set by the user to the */
/* number of variables. It is not altered by the routine. */
/* Restriction: N>0. */
/* M is an INTEGER variable that must be set by the user to */
/* the number of corrections used in the BFGS update. It */
/* is not altered by the routine. Values of M less than 3 are */
/* not recommended; large values of M will result in excessive */
/* computing time. 3<= M <=7 is recommended. Restriction: M>0. */
/* X is a DOUBLE PRECISION array of length N. On initial entry */
/* it must be set by the user to the values of the initial */
/* estimate of the solution vector. On exit with IFLAG=0, it */
/* contains the values of the variables at the best point */
/* found (usually a solution). */
/* F is a DOUBLE PRECISION variable. Before initial entry and on */
/* a re-entry with IFLAG=1, it must be set by the user to */
/* contain the value of the function F at the point X. */
/* G is a DOUBLE PRECISION array of length N. Before initial */
/* entry and on a re-entry with IFLAG=1, it must be set by */
/* the user to contain the components of the gradient G at */
/* the point X. */
/* DIAGCO is a LOGICAL variable that must be set to .TRUE. if the */
/* user wishes to provide the diagonal matrix Hk0 at each */
/* iteration. Otherwise it should be set to .FALSE., in which */
/* case LBFGS will use a default value described below. If */
/* DIAGCO is set to .TRUE. the routine will return at each */
/* iteration of the algorithm with IFLAG=2, and the diagonal */
/* matrix Hk0 must be provided in the array DIAG. */
/* DIAG is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE., */
/* then on initial entry or on re-entry with IFLAG=2, DIAG */
/* it must be set by the user to contain the values of the */
/* diagonal matrix Hk0. Restriction: all elements of DIAG */
/* must be positive. */
/* IPRINT is an INTEGER array of length two which must be set by the */
/* user. */
/* IPRINT(1) specifies the frequency of the output: */
/* IPRINT(1) < 0 : no output is generated, */
/* IPRINT(1) = 0 : output only at first and last iteration, */
/* IPRINT(1) > 0 : output every IPRINT(1) iterations. */
/* IPRINT(2) specifies the type of output generated: */
/* IPRINT(2) = 0 : iteration count, number of function */
/* evaluations, function value, norm of the */
/* gradient, and steplength, */
/* IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of */
/* variables and gradient vector at the */
/* initial point, */
/* IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of */
/* variables, */
/* IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector. */
/* EPS is a positive DOUBLE PRECISION variable that must be set by */
/* the user, and determines the accuracy with which the solution */
/* is to be found. The subroutine terminates when */
/* ||G|| < EPS max(1,||X||), */
/* where ||.|| denotes the Euclidean norm. */
/* XTOL is a positive DOUBLE PRECISION variable that must be set by */
/* the user to an estimate of the machine precision (e.g. */
/* 10**(-16) on a SUN station 3/60). The line search routine will */
/* terminate if the relative width of the interval of uncertainty */
/* is less than XTOL. */
/* W is a DOUBLE PRECISION array of length N(2M+1)+2M used as */
/* workspace for LBFGS. This array must not be altered by the */
/* user. */
/* IFLAG is an INTEGER variable that must be set to 0 on initial entry */
/* to the subroutine. A return with IFLAG<0 indicates an error, */
/* and IFLAG=0 indicates that the routine has terminated without */
/* detecting errors. On a return with IFLAG=1, the user must */
/* evaluate the function F and gradient G. On a return with */
/* IFLAG=2, the user must provide the diagonal matrix Hk0. */
/* The following negative values of IFLAG, detecting an error, */
/* are possible: */
/* IFLAG=-1 The line search routine MCSRCH failed. The */
/* parameter INFO provides more detailed information */
/* (see also the documentation of MCSRCH): */
/* INFO = 0 IMPROPER INPUT PARAMETERS. */
/* INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF */
/* UNCERTAINTY IS AT MOST XTOL. */
/* INFO = 3 MORE THAN 20 FUNCTION EVALUATIONS WERE */
/* REQUIRED AT THE PRESENT ITERATION. */
/* INFO = 4 THE STEP IS TOO SMALL. */
/* INFO = 5 THE STEP IS TOO LARGE. */
/* INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS. */
/* THERE MAY NOT BE A STEP WHICH SATISFIES */
/* THE SUFFICIENT DECREASE AND CURVATURE */
/* CONDITIONS. TOLERANCES MAY BE TOO SMALL. */
/* IFLAG=-2 The i-th diagonal element of the diagonal inverse */
/* Hessian approximation, given in DIAG, is not */
/* positive. */
/* IFLAG=-3 Improper input parameters for LBFGS (N or M are */
/* not positive). */
/* ON THE DRIVER: */
/* The program that calls LBFGS must contain the declaration: */
/* COMMON: */
/* The subroutine contains one common area, which the user may wish to */
/* reference: */
/*< COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX >*/
/* MP is an INTEGER variable with default value 6. It is used as the */
/* unit number for the printing of the monitoring information */
/* controlled by IPRINT. */
/* LP is an INTEGER variable with default value 6. It is used as the */
/* unit number for the printing of error messages. This printing */
/* may be suppressed by setting LP to a non-positive value. */
/* GTOL is a DOUBLE PRECISION variable with default value 0.9, which */
/* controls the accuracy of the line search routine MCSRCH. If the */
/* function and gradient evaluations are inexpensive with respect */
/* to the cost of the iteration (which is sometimes the case when */
/* solving very large problems) it may be advantageous to set GTOL */
/* to a small value. A typical small value is 0.1. Restriction: */
/* GTOL should be greater than 1.D-04. */
/* STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which */
/* specify lower and uper bounds for the step in the line search. */
/* Their default values are 1.D-20 and 1.D+20, respectively. These */
/* values need not be modified unless the exponents are too large */
/* for the machine being used, or unless the problem is extremely */
/* badly scaled (in which case the exponents should be increased). */
/* MACHINE DEPENDENCIES */
/* The only variables that are machine-dependent are XTOL, */
/* STPMIN and STPMAX. */
/* GENERAL INFORMATION */
/* Other routines called directly: DAXPY, DDOT, LB1, MCSRCH */
/* Input/Output : No input; diagnostic messages on unit MP and */
/* error messages on unit LP. */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*< >*/
/*< >*/
/*< LOGICAL FINISH >*/
/*< SAVE >*/
/*< DATA ONE,ZERO/1.0D+0,0.0D+0/ >*/
/* Parameter adjustments */
--diag;
--g;
--x;
--w;
--iprint;
/* Function Body */
/* INITIALIZE */
/* ---------- */
/*< IF(IFLAG.EQ.0) GO TO 10 >*/
if (*iflag == 0) {
goto L10;
}
/*< GO TO (172,100) IFLAG >*/
switch (*iflag) {
case 1: goto L172;
case 2: goto L100;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -