Public feed
#380019
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_ */