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_ */