Public feed
#380018
static int cauchy_(
integer *n, doublereal *x, doublereal *l, doublereal *u, integer *nbd,
doublereal *g, integer *iorder, integer *iwhere, doublereal *t,
doublereal *d__, doublereal *xcp, integer *m, doublereal *wy, doublereal *ws,
doublereal *sy, doublereal *wt, doublereal *theta, integer *col, integer *head,
doublereal *p, doublereal *c__, doublereal *wbp, doublereal *v, integer *nint,
doublereal *sg, doublereal *yg, integer *iprint, doublereal *sbgnrm,
integer *info, doublereal *epsmch)
{
/* Format strings */
/*
static char fmt_3010[] = "(/,\002---------------- CAUCHY entered--------\
-----------\002)";
static char fmt_1010[] = "(\002Cauchy X = \002,/,(4x,1p,6(1x,d11.4)))";
static char fmt_4011[] = "(/,\002Piece \002,i3,\002 --f1, f2 at start\
point \002,1p,2(1x,d11.4))";
static char fmt_5010[] = "(\002Distance to the next break point = \002,\
1p,d11.4)";
static char fmt_6010[] = "(\002Distance to the stationary point = \002,\
1p,d11.4)";
static char fmt_4010[] = "(\002Piece \002,i3,\002 --f1, f2 at start p\
oint \002,1p,2(1x,d11.4))";
static char fmt_2010[] = "(/,\002---------------- exit CAUCHY-----------\
-----------\002,/)";
*/
/* System generated locals */
integer wy_dim1, wy_offset, ws_dim1, ws_offset, sy_dim1, sy_offset,
wt_dim1, wt_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
static doublereal dibp;
static integer iter;
static doublereal zibp, tsum, dibp2;
static integer i__, j;
static logical bnded;
static doublereal neggi;
static integer nfree;
static doublereal bkmin;
static integer nleft;
static doublereal f1, f2, f2_org__, dt, tj, tl;
static integer nbreak, ibkmin;
static doublereal tu;
static integer pointr;
static doublereal tj0;
static logical xlower, xupper;
static integer ibp;
static doublereal dtm;
static doublereal wmc, wmp, wmw;
static integer col2;
/* Fortran I/O blocks */
/*
static cilist io___88 = { 0, 6, 0, 0, 0 };
static cilist io___96 = { 0, 6, 0, fmt_3010, 0 };
static cilist io___105 = { 0, 6, 0, fmt_1010, 0 };
static cilist io___110 = { 0, 6, 0, 0, 0 };
static cilist io___117 = { 0, 6, 0, fmt_4011, 0 };
static cilist io___118 = { 0, 6, 0, fmt_5010, 0 };
static cilist io___119 = { 0, 6, 0, fmt_6010, 0 };
static cilist io___122 = { 0, 6, 0, 0, 0 };
static cilist io___127 = { 0, 6, 0, 0, 0 };
static cilist io___128 = { 0, 6, 0, 0, 0 };
static cilist io___129 = { 0, 6, 0, fmt_4010, 0 };
static cilist io___130 = { 0, 6, 0, fmt_6010, 0 };
static cilist io___131 = { 0, 6, 0, fmt_1010, 0 };
static cilist io___132 = { 0, 6, 0, fmt_2010, 0 };
*/
/* ************ */
/* Subroutine cauchy */
/* For given x, l, u, g (with sbgnrm > 0), and a limited memory */
/* BFGS matrix B defined in terms of matrices WY, WS, WT, and */
/* scalars head, col, and theta, this subroutine computes the */
/* generalized Cauchy point (GCP), defined as the first local */
/* minimizer of the quadratic */
/* Q(x + s) = g's + 1/2 s'Bs */
/* along the projected gradient direction P(x-tg,l,u). */
/* The routine returns the GCP in xcp. */
/* n is an integer variable. */
/* On entry n is the dimension of the problem. */
/* On exit n is unchanged. */
/* x is a double precision array of dimension n. */
/* On entry x is the starting point for the GCP computation. */
/* On exit x is unchanged. */
/* l is a double precision array of dimension n. */
/* On entry l is the lower bound of x. */
/* On exit l is unchanged. */
/* u is a double precision array of dimension n. */
/* On entry u is the upper bound of x. */
/* On exit u is unchanged. */
/* nbd is an integer array of dimension n. */
/* On entry nbd represents the type of bounds imposed on the */
/* variables, and must be specified as follows: */
/* nbd(i)=0 if x(i) is unbounded, */
/* 1 if x(i) has only a lower bound, */
/* 2 if x(i) has both lower and upper bounds, and */
/* 3 if x(i) has only an upper bound. */
/* On exit nbd is unchanged. */
/* g is a double precision array of dimension n. */
/* On entry g is the gradient of f(x). g must be a nonzero vector. */
/* On exit g is unchanged. */
/* iorder is an integer working array of dimension n. */
/* iorder will be used to store the breakpoints in the piecewise */
/* linear path and free variables encountered. On exit, */
/* iorder(1),...,iorder(nleft) are indices of breakpoints */
/* which have not been encountered; */
/* iorder(nleft+1),...,iorder(nbreak) are indices of */
/* encountered breakpoints; and */
/* iorder(nfree),...,iorder(n) are indices of variables which */
/* have no bound constraits along the search direction. */
/* iwhere is an integer array of dimension n. */
/* On entry iwhere indicates only the permanently fixed (iwhere=3) */
/* or free (iwhere= -1) components of x. */
/* On exit iwhere records the status of the current x variables. */
/* iwhere(i)=-3 if x(i) is free and has bounds, but is not moved */
/* 0 if x(i) is free and has bounds, and is moved */
/* 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) */
/* 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) */
/* 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) */
/* -1 if x(i) is always free, i.e., it has no bounds. */
/* t is a double precision working array of dimension n. */
/* t will be used to store the break points. */
/* d is a double precision array of dimension n used to store */
/* the Cauchy direction P(x-tg)-x. */
/* xcp is a double precision array of dimension n used to return the */
/* GCP on exit. */
/* m is an integer variable. */
/* On entry m is the maximum number of variable metric corrections */
/* used to define the limited memory matrix. */
/* On exit m is unchanged. */
/* ws, wy, sy, and wt are double precision arrays. */
/* On entry they store information that defines the */
/* limited memory BFGS matrix: */
/* ws(n,m) stores S, a set of s-vectors; */
/* wy(n,m) stores Y, a set of y-vectors; */
/* sy(m,m) stores S'Y; */
/* wt(m,m) stores the */
/* Cholesky factorization of (theta*S'S+LD^(-1)L'). */
/* On exit these arrays are unchanged. */
/* theta is a double precision variable. */
/* On entry theta is the scaling factor specifying B_0 = theta I. */
/* On exit theta is unchanged. */
/* col is an integer variable. */
/* On entry col is the actual number of variable metric */
/* corrections stored so far. */
/* On exit col is unchanged. */
/* head is an integer variable. */
/* On entry head is the location of the first s-vector (or y-vector) */
/* in S (or Y). */
/* On exit col is unchanged. */
/* p is a double precision working array of dimension 2m. */
/* p will be used to store the vector p = W^(T)d. */
/* c is a double precision working array of dimension 2m. */
/* c will be used to store the vector c = W^(T)(xcp-x). */
/* wbp is a double precision working array of dimension 2m. */
/* wbp will be used to store the row of W corresponding */
/* to a breakpoint. */
/* v is a double precision working array of dimension 2m. */
/* nint is an integer variable. */
/* On exit nint records the number of quadratic segments explored */
/* in searching for the GCP. */
/* sg and yg are double precision arrays of dimension m. */
/* On entry sg and yg store S'g and Y'g correspondingly. */
/* On exit they are unchanged. */
/* iprint is an INTEGER variable that must be set by the user. */
/* It controls the frequency and type of output generated: */
/* iprint<0 no output is generated; */
/* iprint=0 print only one line at the last iteration; */
/* 0<iprint<99 print also f and |proj g| every iprint iterations; */
/* iprint=99 print details of every iteration except n-vectors; */
/* iprint=100 print also the changes of active set and final x; */
/* iprint>100 print details of every iteration including x and g; */
/* When iprint > 0, the file iterate.dat will be created to */
/* summarize the iteration. */
/* sbgnrm is a double precision variable. */
/* On entry sbgnrm is the norm of the projected gradient at x. */
/* On exit sbgnrm is unchanged. */
/* info is an integer variable. */
/* On entry info is 0. */
/* On exit info = 0 for normal return, */
/* = nonzero for abnormal return when the the system */
/* used in routine bmv is singular. */
/* Subprograms called: */
/* L-BFGS-B Library ... hpsolb, bmv. */
/* Linpack ... dscal dcopy, daxpy. */
/* References: */
/* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
/* memory algorithm for bound constrained optimization'', */
/* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
/* [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN */
/* Subroutines for Large Scale Bound Constrained Optimization'' */
/* Tech. Report, NAM-11, EECS Department, Northwestern University, */
/* 1994. */
/* (Postscript files of these papers are available via anonymous */
/* ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
/* * * * */
/* NEOS, November 1994. (Latest revision June 1996.) */
/* Optimization Technology Center. */
/* Argonne National Laboratory and Northwestern University. */
/* Written by */
/* Ciyou Zhu */
/* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
/* ************ */
/* Check the status of the variables, reset iwhere(i) if necessary; */
/* compute the Cauchy direction d and the breakpoints t; initialize */
/* the derivative f1 and the vector p = W'd (for theta = 1). */
/* Parameter adjustments */
--xcp;
--d__;
--t;
--iwhere;
--iorder;
--g;
--nbd;
--u;
--l;
--x;
--yg;
--sg;
--v;
--wbp;
--c__;
--p;
wt_dim1 = *m;
wt_offset = 1 + wt_dim1 * 1;
wt -= wt_offset;
sy_dim1 = *m;
sy_offset = 1 + sy_dim1 * 1;
sy -= sy_offset;
ws_dim1 = *n;
ws_offset = 1 + ws_dim1 * 1;
ws -= ws_offset;
wy_dim1 = *n;
wy_offset = 1 + wy_dim1 * 1;
wy -= wy_offset;
/* Function Body */
if (*sbgnrm <= 0.) {
if (*iprint >= 0) {
/*
s_wsle(&io___88);
do_lio(&c__9, &c__1, "Subgnorm = 0. GCP = X.", (ftnlen)23);
e_wsle();
*/
}
dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
return 0;
}
bnded = TRUE_;
nfree = *n + 1;
nbreak = 0;
ibkmin = 0;
bkmin = 0.;
col2 = *col << 1;
f1 = 0.;
if (*iprint >= 99) {
/*
s_wsfe(&io___96);
e_wsfe();
*/
}
/* We set p to zero and build it up as we determine d. */
i__1 = col2;
for (i__ = 1; i__ <= i__1; ++i__) {
p[i__] = 0.;
/* L20: */
}
/* In the following loop we determine for each variable its bound */
/* status and its breakpoint, and update p accordingly. */
/* Smallest breakpoint is identified. */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
neggi = -g[i__];
if (iwhere[i__] != 3 && iwhere[i__] != -1) {
/* if x(i) is not a constant and has bounds, */
/* compute the difference between x(i) and its bounds. */
if (nbd[i__] <= 2) {
tl = x[i__] - l[i__];
}
if (nbd[i__] >= 2) {
tu = u[i__] - x[i__];
}
/* If a variable is close enough to a bound */
/* we treat it as at bound. */
xlower = nbd[i__] <= 2 && tl <= 0.;
xupper = nbd[i__] >= 2 && tu <= 0.;
/* reset iwhere(i). */
iwhere[i__] = 0;
if (xlower) {
if (neggi <= 0.) {
iwhere[i__] = 1;
}
} else if (xupper) {
if (neggi >= 0.) {
iwhere[i__] = 2;
}
} else {
if (abs(neggi) <= 0.) {
iwhere[i__] = -3;
}
}
}
pointr = *head;
if (iwhere[i__] != 0 && iwhere[i__] != -1) {
d__[i__] = 0.;
} else {
d__[i__] = neggi;
f1 -= neggi * neggi;
/* calculate p := p - W'e_i* (g_i). */
i__2 = *col;
for (j = 1; j <= i__2; ++j) {
p[j] += wy[i__ + pointr * wy_dim1] * neggi;
p[*col + j] += ws[i__ + pointr * ws_dim1] * neggi;
pointr = pointr % *m + 1;
/* L40: */
}
if (nbd[i__] <= 2 && nbd[i__] != 0 && neggi < 0.) {
/* x(i) + d(i) is bounded; compute t(i). */
++nbreak;
iorder[nbreak] = i__;
t[nbreak] = tl / (-neggi);
if (nbreak == 1 || t[nbreak] < bkmin) {
bkmin = t[nbreak];
ibkmin = nbreak;
}
} else if (nbd[i__] >= 2 && neggi > 0.) {
/* x(i) + d(i) is bounded; compute t(i). */
++nbreak;
iorder[nbreak] = i__;
t[nbreak] = tu / neggi;
if (nbreak == 1 || t[nbreak] < bkmin) {
bkmin = t[nbreak];
ibkmin = nbreak;
}
} else {
/* x(i) + d(i) is not bounded. */
--nfree;
iorder[nfree] = i__;
if (abs(neggi) > 0.) {
bnded = FALSE_;
}
}
}
/* L50: */
}
/* The indices of the nonzero components of d are now stored */
/* in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n). */
/* The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin. */
if (*theta != 1.) {
/* complete the initialization of p for theta not= one. */
dscal_(col, theta, &p[*col + 1], &c__1);
}
/* Initialize GCP xcp = x. */
dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
if (nbreak == 0 && nfree == *n + 1) {
/* is a zero vector, return with the initial xcp as GCP. */
if (*iprint > 100) {
/*
s_wsfe(&io___105);
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
do_fio(&c__1, (char *)&xcp[i__], (ftnlen)sizeof(doublereal));
}
e_wsfe();
*/
}
return 0;
}
/* Initialize c = W'(xcp - x) = 0. */
i__1 = col2;
for (j = 1; j <= i__1; ++j) {
c__[j] = 0.;
/* L60: */
}
/* Initialize derivative f2. */
f2 = -(*theta) * f1;
f2_org__ = f2;
if (*col > 0) {
bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &p[1], &v[1], info);
if (*info != 0) {
return 0;
}
f2 -= ddot_(&col2, &v[1], &c__1, &p[1], &c__1);
}
dtm = -f1 / f2;
tsum = 0.;
*nint = 1;
if (*iprint >= 99) {
/*
s_wsle(&io___110);
do_lio(&c__9, &c__1, "There are ", (ftnlen)10);
do_lio(&c__3, &c__1, (char *)&nbreak, (ftnlen)sizeof(integer));
do_lio(&c__9, &c__1, " breakpoints ", (ftnlen)14);
e_wsle();
*/
}
/* If there are no breakpoints, locate the GCP and return. */
if (nbreak == 0) {
goto L888;
}
nleft = nbreak;
iter = 1;
tj = 0.;
/* ------------------- the beginning of the loop ------------------------- */
L777:
/* Find the next smallest breakpoint; */
/* compute dt = t(nleft) - t(nleft + 1). */
tj0 = tj;
if (iter == 1) {
/* Since we already have the smallest breakpoint we need not do */
/* heapsort yet. Often only one breakpoint is used and the */
/* cost of heapsort is avoided. */
tj = bkmin;
ibp = iorder[ibkmin];
} else {
if (iter == 2) {
/* Replace the already used smallest breakpoint with the */
/* breakpoint numbered nbreak > nlast, before heapsort call. */
if (ibkmin != nbreak) {
t[ibkmin] = t[nbreak];
iorder[ibkmin] = iorder[nbreak];
}
/* Update heap structure of breakpoints */
/* (if iter=2, initialize heap). */
}
i__1 = iter - 2;
hpsolb_(&nleft, &t[1], &iorder[1], &i__1);
tj = t[nleft];
ibp = iorder[nleft];
}
dt = tj - tj0;
if (dt != 0. && *iprint >= 100) {
/*
s_wsfe(&io___117);
do_fio(&c__1, (char *)&(*nint), (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&f1, (ftnlen)sizeof(doublereal));
do_fio(&c__1, (char *)&f2, (ftnlen)sizeof(doublereal));
e_wsfe();
s_wsfe(&io___118);
do_fio(&c__1, (char *)&dt, (ftnlen)sizeof(doublereal));
e_wsfe();
s_wsfe(&io___119);
do_fio(&c__1, (char *)&dtm, (ftnlen)sizeof(doublereal));
e_wsfe();
*/
}
/* If a minimizer is within this interval, locate the GCP and return. */
if (dtm < dt) {
goto L888;
}
/* Otherwise fix one variable and */
/* reset the corresponding component of d to zero. */
tsum += dt;
--nleft;
++iter;
dibp = d__[ibp];
d__[ibp] = 0.;
if (dibp > 0.) {
zibp = u[ibp] - x[ibp];
xcp[ibp] = u[ibp];
iwhere[ibp] = 2;
} else {
zibp = l[ibp] - x[ibp];
xcp[ibp] = l[ibp];
iwhere[ibp] = 1;
}
if (*iprint >= 100) {
/*
s_wsle(&io___122);
do_lio(&c__9, &c__1, "Variable ", (ftnlen)10);
do_lio(&c__3, &c__1, (char *)&ibp, (ftnlen)sizeof(integer));
do_lio(&c__9, &c__1, " is fixed.", (ftnlen)11);
e_wsle();
*/
}
if (nleft == 0 && nbreak == *n) {
/* all n variables are fixed, */
/* return with xcp as GCP. */
dtm = dt;
goto L999;
}
/* Update the derivative information. */
++(*nint);
/* Computing 2nd power */
d__1 = dibp;
dibp2 = d__1 * d__1;
/* Update f1 and f2. */
/* temporarily set f1 and f2 for col=0. */
f1 = f1 + dt * f2 + dibp2 - *theta * dibp * zibp;
f2 -= *theta * dibp2;
if (*col > 0) {
/* update c = c + dt*p. */
daxpy_(&col2, &dt, &p[1], &c__1, &c__[1], &c__1);
/* choose wbp, */
/* the row of W corresponding to the breakpoint encountered. */
pointr = *head;
i__1 = *col;
for (j = 1; j <= i__1; ++j) {
wbp[j] = wy[ibp + pointr * wy_dim1];
wbp[*col + j] = *theta * ws[ibp + pointr * ws_dim1];
pointr = pointr % *m + 1;
/* L70: */
}
/* compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'. */
bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wbp[1], &v[1], info);
if (*info != 0) {
return 0;
}
wmc = ddot_(&col2, &c__[1], &c__1, &v[1], &c__1);
wmp = ddot_(&col2, &p[1], &c__1, &v[1], &c__1);
wmw = ddot_(&col2, &wbp[1], &c__1, &v[1], &c__1);
/* update p = p - dibp*wbp. */
d__1 = -dibp;
daxpy_(&col2, &d__1, &wbp[1], &c__1, &p[1], &c__1);
/* complete updating f1 and f2 while col > 0. */
f1 += dibp * wmc;
f2 = f2 + dibp * 2. * wmp - dibp2 * wmw;
}
/* Computing MAX */
d__1 = *epsmch * f2_org__;
f2 = max(d__1,f2);
if (nleft > 0) {
dtm = -f1 / f2;
goto L777;
/* to repeat the loop for unsearched intervals. */
} else if (bnded) {
f1 = 0.;
f2 = 0.;
dtm = 0.;
} else {
dtm = -f1 / f2;
}
/* ------------------- the end of the loop ------------------------------- */
L888:
if (*iprint >= 99) {
/*
s_wsle(&io___127);
e_wsle();
s_wsle(&io___128);
do_lio(&c__9, &c__1, "GCP found in this segment", (ftnlen)25);
e_wsle();
s_wsfe(&io___129);
do_fio(&c__1, (char *)&(*nint), (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&f1, (ftnlen)sizeof(doublereal));
do_fio(&c__1, (char *)&f2, (ftnlen)sizeof(doublereal));
e_wsfe();
s_wsfe(&io___130);
do_fio(&c__1, (char *)&dtm, (ftnlen)sizeof(doublereal));
e_wsfe();
*/
}
if (dtm <= 0.) {
dtm = 0.;
}
tsum += dtm;
/* Move free variables (i.e., the ones w/o breakpoints) and */
/* the variables whose breakpoints haven't been reached. */
daxpy_(n, &tsum, &d__[1], &c__1, &xcp[1], &c__1);
L999:
/* Update c = c + dtm*p = W'(x^c - x) */
/* which will be used in computing r = Z'(B(x^c - x) + g). */
if (*col > 0) {
daxpy_(&col2, &dtm, &p[1], &c__1, &c__[1], &c__1);
}
if (*iprint > 100) {
/*
s_wsfe(&io___131);
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
do_fio(&c__1, (char *)&xcp[i__], (ftnlen)sizeof(doublereal));
}
e_wsfe();
*/
}
if (*iprint >= 99) {
/*
s_wsfe(&io___132);
e_wsfe();
*/
}
return 0;
} /* cauchy_ */