// routines.f -- translated by f2c (version 20031025).

#include <ctime>

#ifdef __cplusplus
extern "C" {
#endif
#include "routines.h"
//#include "f2c.h"

// Table of constant values

static integer c__1 = 1;
static integer c__9 = 9;
static integer c__3 = 3;
static integer c__5 = 5;

// ================    L-BFGS-B (version 2.1)   ========================== 
int setulb_(const integer n, const integer m, doublereal *x, 
	const doublereal *l, const doublereal *u, const integer *nbd, 
	doublereal *f, doublereal *g, const doublereal factr, 
	const doublereal pgtol, doublereal *wa, integer *iwa, char *task, 
	const integer iprint, char *csave, logical *lsave, integer *isave, 
	doublereal *dsave, ftnlen task_len, ftnlen csave_len)
{
    // Builtin functions 
    integer s_cmp(char *, char *, ftnlen, ftnlen);

    // Local variables 
    static integer l1, l2, l3, ld, lr, lt, lz, lwa, lsg, lyg, lwn, lss, lws, 
	    lwt, lsy, lwy, lyy, lsnd, lsgo, lygo;
	    
    // Parameter adjustments 
    --iwa;
    --g;
    --nbd;
    --u;
    --l;
    --x;
    --wa;
    --lsave;
    --isave;
    --dsave;

    // Function Body 
    if (s_cmp(task, "START", (ftnlen)60, (ftnlen)5) == 0) {
		isave[1] = m * n;
	// Computing 2nd power 
		isave[2] = m * m;
	// Computing 2nd power 
		isave[3] = m * m << 2;
		isave[4] = 1;
		isave[5] = isave[4] + isave[1];
		isave[6] = isave[5] + isave[1];
		isave[7] = isave[6] + isave[2];
		isave[8] = isave[7] + isave[2];
		isave[9] = isave[8] + isave[2];
		isave[10] = isave[9] + isave[2];
		isave[11] = isave[10] + isave[3];
		isave[12] = isave[11] + isave[3];
		isave[13] = isave[12] + n;
		isave[14] = isave[13] + n;
		isave[15] = isave[14] + n;
		isave[16] = isave[15] + n;
		isave[17] = isave[16] + (m << 3);
		isave[18] = isave[17] + m;
		isave[19] = isave[18] + m;
		isave[20] = isave[19] + m;
    }
    l1 = isave[1];
    l2 = isave[2];
    l3 = isave[3];
    lws = isave[4];
    lwy = isave[5];
    lsy = isave[6];
    lss = isave[7];
    lyy = isave[8];
    lwt = isave[9];
    lwn = isave[10];
    lsnd = isave[11];
    lz = isave[12];
    lr = isave[13];
    ld = isave[14];
    lt = isave[15];
    lwa = isave[16];
    lsg = isave[17];
    lsgo = isave[18];
    lyg = isave[19];
    lygo = isave[20];
    mainlb_(n, m, &x[1], &l[1], &u[1], &nbd[1], f, &g[1], factr, pgtol, &wa[
	    lws], &wa[lwy], &wa[lsy], &wa[lss], &wa[lyy], &wa[lwt], &wa[lwn], 
	    &wa[lsnd], &wa[lz], &wa[lr], &wa[ld], &wa[lt], &wa[lwa], &wa[lsg],
	     &wa[lsgo], &wa[lyg], &wa[lygo], &iwa[1], &iwa[n + 1], &iwa[(n 
	    << 1) + 1], task, iprint, csave, &lsave[1], &isave[22], &dsave[1],
	     (ftnlen)60, (ftnlen)60);
    return 0;
} // setulb_ 

// ======================= The end of setulb ============================= 
int mainlb_(const integer n, const integer m, doublereal *x, 
	const doublereal *l, const doublereal *u, const integer *nbd, doublereal *f, doublereal 
	*g, const doublereal factr, const doublereal pgtol, doublereal *ws, doublereal *
	wy, doublereal *sy, doublereal *ss, doublereal *yy, doublereal *wt, 
	doublereal *wn, doublereal *snd, doublereal *z__, doublereal *r__, 
	doublereal *d__, doublereal *t, doublereal *wa, doublereal *sg, 
	doublereal *sgo, doublereal *yg, doublereal *ygo, integer *index, 
	integer *iwhere, integer *indx2, char *task, const integer iprint, char *
	csave, logical *lsave, integer *isave, doublereal *dsave, ftnlen 
	task_len, ftnlen csave_len)
{
    // Format strings 
    static char fmt_1002[] = "(/,\002At iterate\002,i5,4x,\002f= \002,1p,d12\
.5,4x,\002|proj g|= \002,1p,d12.5)";
    static char fmt_1003[] = "(2(1x,i4),5x,\002-\002,5x,\002-\002,3x,\002\
-\002,5x,\002-\002,5x,\002-\002,8x,\002-\002,3x,1p,2(1x,d10.3))";
    static char fmt_1001[] = "(//,\002ITERATION \002,i5)";
    static char fmt_1005[] = "(/,\002 Singular triangular system detected\
;\002,/,\002   refresh the lbfgs memory and restart the iteration.\002)";
    static char fmt_1006[] = "(/,\002 Nonpositive definiteness in Cholesky f\
actorization in formk;\002,/,\002   refresh the lbfgs memory and restart the\
 iteration.\002)";
    static char fmt_1008[] = "(/,\002 Bad direction in the line search;\002,\
/,\002   refresh the lbfgs memory and restart the iteration.\002)";
    static char fmt_1004[] = "(\002  ys=\002,1p,e10.3,\002  -gs=\002,1p,e10.\
3,\002 BFGS update SKIPPED\002)";
    static char fmt_1007[] = "(/,\002 Nonpositive definiteness in Cholesky f\
actorization in formt;\002,/,\002   refresh the lbfgs memory and restart the\
 iteration.\002)";

    // System generated locals 
    integer ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, 
	    ss_dim1, ss_offset, yy_dim1, yy_offset, wt_dim1, wt_offset, 
	    wn_dim1, wn_offset, snd_dim1, snd_offset, i__1;
    doublereal d__1, d__2;
    olist o__1;

    // Builtin functions 
    integer s_cmp(char *, char *, ftnlen, ftnlen);
    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
    integer f_open(olist *), s_wsfe(cilist *), do_fio(integer *, char *, 
	    ftnlen), e_wsfe();

    // Local variables 
    static integer i__, k;
    static doublereal gd, dr, rr, dtd;
    static integer col;
    static doublereal tol;
    static logical wrk;
    static doublereal stp, cpu1, cpu2;
    static integer head;
    static doublereal fold;
    static integer nact;
    static doublereal ddum;
    static integer info;
    static doublereal time;
    static integer nfgv, ifun, iter, nint;
    static char word[3];
    static doublereal time1, time2;
    static integer iback;
    static doublereal gdold;
    static integer nfree;
    static logical boxed;
    static integer itail;
    static doublereal theta;
    static doublereal dnorm;
    static integer nskip, iword;
    static doublereal xstep, stpmx;
    static integer ileave;
    static doublereal cachyt;
    static integer itfile;
    static doublereal epsmch;
    static logical updatd;
    static doublereal sbtime;
    static logical prjctd;
    static integer iupdat;
    static logical cnstnd;
    static doublereal sbgnrm;
    static integer nenter;
    static doublereal lnscht;
    static integer nintol;

    // Fortran I/O blocks 
    static cilist io___62 = { 0, 6, 0, fmt_1002, 0 };
    static cilist io___63 = { 0, 0, 0, fmt_1003, 0 };
    static cilist io___64 = { 0, 6, 0, fmt_1001, 0 };
    static cilist io___66 = { 0, 6, 0, fmt_1005, 0 };
    static cilist io___68 = { 0, 6, 0, fmt_1006, 0 };
    static cilist io___69 = { 0, 6, 0, fmt_1005, 0 };
    static cilist io___71 = { 0, 6, 0, fmt_1008, 0 };
    static cilist io___75 = { 0, 6, 0, fmt_1004, 0 };
    static cilist io___76 = { 0, 6, 0, fmt_1007, 0 };


    // Parameter adjustments 
    --indx2;
    --iwhere;
    --index;
    --t;
    --d__;
    --r__;
    --z__;
    --g;
    --nbd;
    --u;
    --l;
    --x;
    --ygo;
    --yg;
    --sgo;
    --sg;
    --wa;
    snd_dim1 = 2 * m;
    snd_offset = 1 + snd_dim1;
    snd -= snd_offset;
    wn_dim1 = 2 * m;
    wn_offset = 1 + wn_dim1;
    wn -= wn_offset;
    wt_dim1 = m;
    wt_offset = 1 + wt_dim1;
    wt -= wt_offset;
    yy_dim1 = m;
    yy_offset = 1 + yy_dim1;
    yy -= yy_offset;
    ss_dim1 = m;
    ss_offset = 1 + ss_dim1;
    ss -= ss_offset;
    sy_dim1 = m;
    sy_offset = 1 + sy_dim1;
    sy -= sy_offset;
    wy_dim1 = n;
    wy_offset = 1 + wy_dim1;
    wy -= wy_offset;
    ws_dim1 = n;
    ws_offset = 1 + ws_dim1;
    ws -= ws_offset;
    --lsave;
    --isave;
    --dsave;

    /* Function Body */
    if (s_cmp(task, "START", (ftnlen)60, (ftnlen)5) == 0) {
		timer_(&time1);
	/*        Generate the current machine precision. */
		epsmch = dpmeps_();
	/*        Initialize counters and scalars when task='START'. */
	/*           for the limited memory BFGS matrices: */
		col = 0;
		head = 1;
		theta = 1.;
		iupdat = 0;
		updatd = FALSE_;
	/*           for operation counts: */
		iter = 0;
		nfgv = 0;
		nint = 0;
		nintol = 0;
		nskip = 0;
		nfree = n;
	/*           for stopping tolerance: */
		tol = factr * epsmch;
	/*           for measuring running time: */
		cachyt = 0.;
		sbtime = 0.;
		lnscht = 0.;
	/*           'word' records the status of subspace solutions. */
		s_copy(word, "---", (ftnlen)3, (ftnlen)3);
	/*           'info' records the termination information. */
		info = 0;
		if (iprint >= 1) {
	/*                                open a summary file 'iterate.dat' */
		    o__1.oerr = 0;
		    o__1.ounit = 8;
		    o__1.ofnmlen = 11;
		    o__1.ofnm = "iterate.dat";
		    o__1.orl = 0;
		    o__1.osta = "unknown";
		    o__1.oacc = 0;
		    o__1.ofm = 0;
		    o__1.oblnk = 0;
		    f_open(&o__1);
		    itfile = 8;
		}
	/*        Check the input arguments for errors. */
		errclb_(n, m, factr, &l[1], &u[1], &nbd[1], task, &info, &k, (ftnlen)
			60);
		if (s_cmp(task, "ERROR", (ftnlen)5, (ftnlen)5) == 0) {
		    prn3lb_(n, &x[1], *f, task, iprint, info, &itfile, iter, nfgv, 
			    nintol, nskip, nact, sbgnrm, 0., nint, word, 
			    iback, stp, xstep, k, cachyt, sbtime, lnscht, (
			    ftnlen)60, (ftnlen)3);
		    return 0;
		}
		prn1lb_(n, m, &l[1], &u[1], &x[1], iprint, &itfile, epsmch);
	/*        Initialize iwhere & project x onto the feasible set. */
		active_(n, &l[1], &u[1], &nbd[1], &x[1], &iwhere[1], iprint, &prjctd, 
			&cnstnd, &boxed);
	/*        The end of the initialization. */
    } 
    else {
/*          restore local variables. */
	prjctd = lsave[1];
	cnstnd = lsave[2];
	boxed = lsave[3];
	updatd = lsave[4];
	nintol = isave[1];
	itfile = isave[3];
	iback = isave[4];
	nskip = isave[5];
	head = isave[6];
	col = isave[7];
	itail = isave[8];
	iter = isave[9];
	iupdat = isave[10];
	nint = isave[12];
	nfgv = isave[13];
	info = isave[14];
	ifun = isave[15];
	iword = isave[16];
	nfree = isave[17];
	nact = isave[18];
	ileave = isave[19];
	nenter = isave[20];
	theta = dsave[1];
	fold = dsave[2];
	tol = dsave[3];
	dnorm = dsave[4];
	epsmch = dsave[5];
	cpu1 = dsave[6];
	cachyt = dsave[7];
	sbtime = dsave[8];
	lnscht = dsave[9];
	time1 = dsave[10];
	gd = dsave[11];
	stpmx = dsave[12];
	sbgnrm = dsave[13];
	stp = dsave[14];
	gdold = dsave[15];
	dtd = dsave[16];
/*        After returning from the driver go to the point where execution */
/*        is to resume. */
	if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) {
	    goto L666;
	}
	if (s_cmp(task, "NEW_X", (ftnlen)5, (ftnlen)5) == 0) {
	    goto L777;
	}
	if (s_cmp(task, "FG_ST", (ftnlen)5, (ftnlen)5) == 0) {
	    goto L111;
	}
	if (s_cmp(task, "STOP", (ftnlen)4, (ftnlen)4) == 0) {
	    if (s_cmp(task + 6, "CPU", (ftnlen)3, (ftnlen)3) == 0) {
	/*                                          restore the previous iterate. */
			dcopy_(n, &t[1], 1, &x[1], 1);
			dcopy_(n, &r__[1], 1, &g[1], 1);
			*f = fold;
	    }
	    goto L999;
	}
    }
/*     Compute f0 and g0. */
    s_copy(task, "FG_START", (ftnlen)60, (ftnlen)8);
/*          return to the driver to calculate f and g; reenter at 111. */
    goto L1000;
L111:
    nfgv = 1;
/*     Compute the infinity norm of the (-) projected gradient. */
    projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm);
    if (iprint >= 1) {
		s_wsfe(&io___62);
		do_fio(&c__1, (char *)&iter, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(doublereal));
		do_fio(&c__1, (char *)&sbgnrm, (ftnlen)sizeof(doublereal));
		e_wsfe();
		io___63.ciunit = itfile;
		s_wsfe(&io___63);
		do_fio(&c__1, (char *)&iter, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&nfgv, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&sbgnrm, (ftnlen)sizeof(doublereal));
		do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(doublereal));
		e_wsfe();
    }
    if (sbgnrm <= pgtol) {
/*                                terminate the algorithm. */
	s_copy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL", (
		ftnlen)60, (ftnlen)48);
	goto L999;
    }
/* ----------------- the beginning of the loop -------------------------- */
L222:
    if (iprint >= 99) {
		s_wsfe(&io___64);
		i__1 = iter + 1;
		do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer));
		e_wsfe();
    }
    iword = -1;

    if (! cnstnd && col > 0) {
/*                                            skip the search for GCP. */
		dcopy_(n, &x[1], 1, &z__[1], 1);
		wrk = updatd;
		nint = 0;
		goto L333;
    }
/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */

/*     Compute the Generalized Cauchy Point (GCP). */

/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
    timer_(&cpu1);
    cauchy_(n, &x[1], &l[1], &u[1], &nbd[1], &g[1], &indx2[1], &iwhere[1], &t[
	    1], &d__[1], &z__[1], m, &wy[wy_offset], &ws[ws_offset], &sy[
	    sy_offset], &wt[wt_offset], theta, col, head, &wa[1], &wa[(m 
	    << 1) + 1], &wa[(m << 2) + 1], &wa[m * 6 + 1], &nint, &sg[1], &
	    yg[1], iprint, sbgnrm, &info, &epsmch);
    if (info != 0) {
/*         singular triangular system detected; refresh the lbfgs memory. */
	if (iprint >= 1) {
	    s_wsfe(&io___66);
	    e_wsfe();
	}
	info = 0;
	col = 0;
	head = 1;
	theta = 1.;
	iupdat = 0;
	updatd = FALSE_;
	timer_(&cpu2);
	cachyt = cachyt + cpu2 - cpu1;
	goto L222;
    }
    timer_(&cpu2);
    cachyt = cachyt + cpu2 - cpu1;
    nintol += nint;
/*     Count the entering and leaving variables for iter > 0; */
/*     find the index set of free and active variables at the GCP. */
    freev_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iwhere[1], &
	    wrk, updatd, cnstnd, iprint, iter);
    nact = n - nfree;
L333:
/*     If there are no free variables or B=theta*I, then */
/*                                        skip the subspace minimization. */
    if (nfree == 0 || col == 0) goto L555;
/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */

/*     Subspace minimization. */

/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
    timer_(&cpu1);
/*     Form  the LEL^T factorization of the indefinite */
/*       matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ] */
/*                     [L_a -R_z           theta*S'AA'S ] */
/*       where     E = [-I  0] */
/*                     [ 0  I] */
    if (wrk) {
	formk_(n, nfree, &index[1], nenter, ileave, &indx2[1], iupdat, 
		updatd, &wn[wn_offset], &snd[snd_offset], m, &ws[ws_offset], &
		wy[wy_offset], &sy[sy_offset], theta, col, head, &info);
    }
    if (info != 0) {
/*          nonpositive definiteness in Cholesky factorization; */
/*          refresh the lbfgs memory and restart the iteration. */
	if (iprint >= 1) {
	    s_wsfe(&io___68);
	    e_wsfe();
	}
	info = 0;
	col = 0;
	head = 1;
	theta = 1.;
	iupdat = 0;
	updatd = FALSE_;
	timer_(&cpu2);
	sbtime = sbtime + cpu2 - cpu1;
	goto L222;
    }
/*        compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) */
/*                                                   from 'cauchy'). */
    cmprlb_(n, m, &x[1], &g[1], &ws[ws_offset], &wy[wy_offset], &sy[sy_offset]
	    , &wt[wt_offset], &z__[1], &r__[1], &wa[1], &index[1], theta, 
	    col, head, nfree, cnstnd, &info);
    if (info != 0) {
	goto L444;
    }
/*       call the direct method. */
    subsm_(n, m, nfree, &index[1], &l[1], &u[1], &nbd[1], &z__[1], &r__[1], &
	    ws[ws_offset], &wy[wy_offset], theta, col, head, &iword, &wa[1]
	    , &wn[wn_offset], iprint, &info);
L444:
    if (info != 0) {
/*          singular triangular system detected; */
/*          refresh the lbfgs memory and restart the iteration. */
	if (iprint >= 1) {
	    s_wsfe(&io___69);
	    e_wsfe();
	}
	info = 0;
	col = 0;
	head = 1;
	theta = 1.;
	iupdat = 0;
	updatd = FALSE_;
	timer_(&cpu2);
	sbtime = sbtime + cpu2 - cpu1;
	goto L222;
    }
    timer_(&cpu2);
    sbtime = sbtime + cpu2 - cpu1;
L555:
/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */

/*     Line search and optimality tests. */

/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
/*     Generate the search direction d:=z-x. */
    for (i__ = 1; i__ <= n; ++i__) {
		d__[i__] = z__[i__] - x[i__];
/* L40: */
    }
    timer_(&cpu1);
L666:
    lnsrlb_(n, &l[1], &u[1], &nbd[1], &x[1], f, &fold, &gd, &gdold, &g[1], &
	    d__[1], &r__[1], &t[1], &z__[1], &stp, &dnorm, &dtd, &xstep, &
	    stpmx, iter, &ifun, &iback, &nfgv, &info, task, boxed, cnstnd, 
	    csave, &isave[22], &dsave[17], (ftnlen)60, (ftnlen)60);
    if (info != 0 || iback >= 20) {
	/*          restore the previous iterate. */
		dcopy_(n, &t[1], 1, &x[1], 1);
		dcopy_(n, &r__[1], 1, &g[1], 1);
		*f = fold;
		if (col == 0) {
	/*             abnormal termination. */
		    if (info == 0) {
			info = -9;
	/*                restore the actual number of f and g evaluations etc. */
			--nfgv;
			--ifun;
			--iback;
		    }
		    s_copy(task, "ABNORMAL_TERMINATION_IN_LNSRCH", (ftnlen)60, (
			    ftnlen)30);
		    ++iter;
		    goto L999;
		} 
		else {
	/*             refresh the lbfgs memory and restart the iteration. */
		    if (iprint >= 1) {
				s_wsfe(&io___71);
				e_wsfe();
		    }
		    if (info == 0) --nfgv;
		    info = 0;
		    col = 0;
		    head = 1;
		    theta = 1.;
		    iupdat = 0;
		    updatd = FALSE_;
		    s_copy(task, "RESTART_FROM_LNSRCH", (ftnlen)60, (ftnlen)19);
		    timer_(&cpu2);
		    lnscht = lnscht + cpu2 - cpu1;
		    goto L222;
		}
    } 
    else if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) {
/*          return to the driver for calculating f and g; reenter at 666. */
		goto L1000;
    } 
    else {
	/*          calculate and print out the quantities related to the new X. */
		timer_(&cpu2);
		lnscht = lnscht + cpu2 - cpu1;
		++iter;
	/*        Compute the infinity norm of the projected (-)gradient. */
		projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm);
	/*        Print iteration information. */
		prn2lb_(n, &x[1], *f, &g[1], iprint, &itfile, iter, nfgv, nact, 
			sbgnrm, nint, word, iword, iback, stp, xstep, (ftnlen)3);
		goto L1000;
    }
L777:
/*     Test for termination. */
    if (sbgnrm <= pgtol) {
	/*                                terminate the algorithm. */
		s_copy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL", (
			ftnlen)60, (ftnlen)48);
		goto L999;
    }
/* Computing MAX */
    d__1 = abs(fold), d__2 = abs(*f), d__1 = max(d__1,d__2);
    ddum = max(d__1,1.);
    if (fold - *f <= tol * ddum) {
	/*                                        terminate the algorithm. */
		s_copy(task, "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH", (
			ftnlen)60, (ftnlen)47);
		if (iback >= 10) {
		    info = -5;
		}
	/*           i.e., to issue a warning if iback>10 in the line search. */
		goto L999;
    }
/*     Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. */
    for (i__ = 1; i__ <= n; ++i__) {
		r__[i__] = g[i__] - r__[i__];
/* L42: */
    }
    rr = ddot_(n, &r__[1], 1, &r__[1], 1);
    if (stp == 1.) {
		dr = gd - gdold;
		ddum = -gdold;
    } else {
		dr = (gd - gdold) * stp;
		dscal_(n, stp, &d__[1], 1);
		ddum = -gdold * stp;
    }
    if (dr <= epsmch * ddum) {
/*                            skip the L-BFGS update. */
	++nskip;
	updatd = FALSE_;
	if (iprint >= 1) {
	    s_wsfe(&io___75);
	    do_fio(&c__1, (char *)&dr, (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&ddum, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	}
	goto L888;
    }
/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */

/*     Update the L-BFGS matrix. */

/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
    updatd = TRUE_;
    ++iupdat;
/*     Update matrices WS and WY and form the middle matrix in B. */
    matupd_(n, m, &ws[ws_offset], &wy[wy_offset], &sy[sy_offset], &ss[
	    ss_offset], &d__[1], &r__[1], &itail, iupdat, &col, &head, &
	    theta, rr, dr, stp, dtd);
/*     Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; */
/*        Store T in the upper triangular of the array wt; */
/*        Cholesky factorize T to J*J' with */
/*           J' stored in the upper triangular of wt. */
    formt_(m, &wt[wt_offset], &sy[sy_offset], &ss[ss_offset], col, theta, &
	    info);
    if (info != 0) {
/*          nonpositive definiteness in Cholesky factorization; */
/*          refresh the lbfgs memory and restart the iteration. */
	if (iprint >= 1) {
	    s_wsfe(&io___76);
	    e_wsfe();
	}
	info = 0;
	col = 0;
	head = 1;
	theta = 1.;
	iupdat = 0;
	updatd = FALSE_;
	goto L222;
    }
/*     Now the inverse of the middle matrix in B is */
/*       [  D^(1/2)      O ] [ -D^(1/2)  D^(-1/2)*L' ] */
/*       [ -L*D^(-1/2)   J ] [  0        J'          ] */
L888:
/* -------------------- the end of the loop ----------------------------- */
    goto L222;
L999:
    timer_(&time2);
    time = time2 - time1;
    prn3lb_(n, &x[1], *f, task, iprint, info, &itfile, iter, nfgv, nintol, 
	    nskip, nact, sbgnrm, time, nint, word, iback, stp, xstep, 
	    k, cachyt, sbtime, lnscht, (ftnlen)60, (ftnlen)3);
L1000:
/*     Save local variables. */
    lsave[1] = prjctd;
    lsave[2] = cnstnd;
    lsave[3] = boxed;
    lsave[4] = updatd;
    isave[1] = nintol;
    isave[3] = itfile;
    isave[4] = iback;
    isave[5] = nskip;
    isave[6] = head;
    isave[7] = col;
    isave[8] = itail;
    isave[9] = iter;
    isave[10] = iupdat;
    isave[12] = nint;
    isave[13] = nfgv;
    isave[14] = info;
    isave[15] = ifun;
    isave[16] = iword;
    isave[17] = nfree;
    isave[18] = nact;
    isave[19] = ileave;
    isave[20] = nenter;
    dsave[1] = theta;
    dsave[2] = fold;
    dsave[3] = tol;
    dsave[4] = dnorm;
    dsave[5] = epsmch;
    dsave[6] = cpu1;
    dsave[7] = cachyt;
    dsave[8] = sbtime;
    dsave[9] = lnscht;
    dsave[10] = time1;
    dsave[11] = gd;
    dsave[12] = stpmx;
    dsave[13] = sbgnrm;
    dsave[14] = stp;
    dsave[15] = gdold;
    dsave[16] = dtd;
    return 0;
} /* mainlb_ */

/* ======================= The end of mainlb ============================= */
int active_(const integer n, const doublereal *l, const doublereal *u, 
	const integer *nbd, doublereal *x, integer *iwhere, const integer iprint, 
	logical *prjctd, logical *cnstnd, logical *boxed)
{
    /* Format strings */
    static char fmt_1001[] = "(/,\002At X0 \002,i9,\002 variables are exactl\
y at the bounds\002)";

    /* Builtin functions */
    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_wsle(), s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), 
	    e_wsfe();

    /* Local variables */
    static integer i__, nbdd;

    /* Fortran I/O blocks */
    static cilist io___81 = { 0, 6, 0, 0, 0 };
    static cilist io___82 = { 0, 6, 0, 0, 0 };
    static cilist io___83 = { 0, 6, 0, fmt_1001, 0 };


    /* Parameter adjustments */
    --iwhere;
    --x;
    --nbd;
    --u;
    --l;

    /* Function Body */
    nbdd = 0;
    *prjctd = FALSE_;
    *cnstnd = FALSE_;
    *boxed = TRUE_;
/*     Project the initial x to the easible set if necessary. */
    for (i__ = 1; i__ <= n; ++i__) {
		if (nbd[i__] > 0) {
		    if (nbd[i__] <= 2 && x[i__] <= l[i__]) {
			if (x[i__] < l[i__]) {
			    *prjctd = TRUE_;
			    x[i__] = l[i__];
			}
			++nbdd;
		    } else if (nbd[i__] >= 2 && x[i__] >= u[i__]) {
			if (x[i__] > u[i__]) {
			    *prjctd = TRUE_;
			    x[i__] = u[i__];
			}
			++nbdd;
		    }
		}
/* L10: */
    }
/*     Initialize iwhere and assign values to cnstnd and boxed. */
    for (i__ = 1; i__ <= n; ++i__) {
		if (nbd[i__] != 2) {
		    *boxed = FALSE_;
		}
		if (nbd[i__] == 0) {
	/*                                this variable is always free */
		    iwhere[i__] = -1;
	/*           otherwise set x(i)=mid(x(i), u(i), l(i)). */
		} 
		else {
		    *cnstnd = TRUE_;
		    if (nbd[i__] == 2 && u[i__] - l[i__] <= 0.) {
	/*                   this variable is always fixed */
			iwhere[i__] = 3;
		    } else {
			iwhere[i__] = 0;
		    }
		}
	/* L20: */
    }
    if (iprint >= 0) {
		if (*prjctd) {
		    s_wsle(&io___81);
		    do_lio(&c__9, &c__1, "The initial X is infeasible.  Restart with\
 its projection.", (ftnlen)58);
		    e_wsle();
		}
		if (! (*cnstnd)) {
		    s_wsle(&io___82);
		    do_lio(&c__9, &c__1, "This problem is unconstrained.", (ftnlen)30)
			    ;
		    e_wsle();
		}
    }
    if (iprint > 0) {
		s_wsfe(&io___83);
		do_fio(&c__1, (char *)&nbdd, (ftnlen)sizeof(integer));
		e_wsfe();
    }
    return 0;
} /* active_ */

/* ======================= The end of active ============================= */
int bmv_(const integer m, const doublereal *sy, const doublereal *wt, 
	const integer col, const doublereal *v, doublereal *p, integer *info)
{
    /* System generated locals */
    integer sy_dim1, sy_offset, wt_dim1, wt_offset, i__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, k, i2;
    static doublereal sum;

    /* Parameter adjustments */
    wt_dim1 = m;
    wt_offset = 1 + wt_dim1;
    wt -= wt_offset;
    sy_dim1 = m;
    sy_offset = 1 + sy_dim1;
    sy -= sy_offset;
    --p;
    --v;

    /* Function Body */
    if (col == 0) return 0;

/*     PART I: solve [  D^(1/2)      O ] [ p1 ] = [ v1 ] */
/*                   [ -L*D^(-1/2)   J ] [ p2 ]   [ v2 ]. */
/* 	solve Jp2=v2+LD^(-1)v1. */
    p[col + 1] = v[col + 1];
    for (i__ = 2; i__ <= col; ++i__) {
		i2 = col + i__;
		sum = 0.;
		i__2 = i__ - 1;
		for (k = 1; k <= i__2; ++k) {
		    sum += sy[i__ + k * sy_dim1] * v[k] / sy[k + k * sy_dim1];
		}
		p[i2] = v[i2] + sum;
    }
/*     Solve the triangular system */
    dtrsl_(&wt[wt_offset], m, col, &p[col + 1], 11, info);
    if (*info != 0) return 0;

/*     	solve D^(1/2)p1=v1. */
    for (i__ = 1; i__ <= col; ++i__) {
		p[i__] = v[i__] / sqrt(sy[i__ + i__ * sy_dim1]);
    }
/*     PART II: solve [ -D^(1/2)   D^(-1/2)*L'  ] [ p1 ] = [ p1 ] */
/*                    [  0         J'           ] [ p2 ]   [ p2 ]. */
/*       solve J^Tp2=p2. */
    dtrsl_(&wt[wt_offset], m, col, &p[col + 1], 1, info);
    if (*info != 0) return 0;
    
/*       compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2) */
/*                 =-D^(-1/2)p1+D^(-1)L'p2. */
    for (i__ = 1; i__ <= col; ++i__) {
		p[i__] = -p[i__] / sqrt(sy[i__ + i__ * sy_dim1]);
    }
    for (i__ = 1; i__ <= col; ++i__) {
		sum = 0.;
		for (k = i__ + 1; k <= col; ++k) {
		    sum += sy[k + i__ * sy_dim1] * p[col + k] / sy[i__ + i__ * sy_dim1];
		}
		p[i__] += sum;
    }
    return 0;
} /* bmv_ */

// ======================== The end of bmv =============================== 
int cauchy_(const integer n, const doublereal *x, const doublereal *l, 
	const doublereal *u, const integer *nbd, doublereal *g, integer *iorder, integer *
	iwhere, doublereal *t, doublereal *d__, doublereal *xcp, const integer m, 
	const doublereal *wy, const doublereal *ws, const doublereal *sy, const doublereal *wt, 
	const doublereal theta, const integer col, const integer head, doublereal *p, 
	doublereal *c__, doublereal *wbp, doublereal *v, integer *nint, 
	const doublereal *sg, const doublereal *yg, const integer iprint, const 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;
    doublereal d__1;

    /* Builtin functions */
    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_wsle(), s_wsfe(cilist *), e_wsfe(), do_fio(integer *, char *, 
	    ftnlen);

    /* Local variables */
    static integer i__, j;
    static doublereal f1, f2, dt, tj, tl, tu, tj0;
    static integer ibp;
    static doublereal dtm;
    static doublereal wmc, wmp, wmw;
    static integer col2;
    static doublereal dibp;
    static integer iter;
    static doublereal zibp, tsum, dibp2;
    static logical bnded;
    static doublereal neggi;
    static integer nfree;
    static doublereal bkmin;
    static integer nleft;
    static doublereal f2_org__;
    static integer nbreak, ibkmin;
    static integer pointr;
    static logical xlower, xupper;

    /* 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 };

    
    // 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;
    wt -= wt_offset;
    sy_dim1 = m;
    sy_offset = 1 + sy_dim1;
    sy -= sy_offset;
    ws_dim1 = n;
    ws_offset = 1 + ws_dim1;
    ws -= ws_offset;
    wy_dim1 = n;
    wy_offset = 1 + wy_dim1;
    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], 1, &xcp[1], 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. */
    for (i__ = 1; i__ <= col2; ++i__) p[i__] = 0.;

/*     In the following loop we determine for each variable its bound 
        status and its breakpoint, and update p accordingly. 
        Smallest breakpoint is identified. */
    for (i__ = 1; i__ <= n; ++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). 
		    for (j = 1; j <= col; ++j) {
				p[j] += wy[i__ + pointr * wy_dim1] * neggi;
				p[col + j] += ws[i__ + pointr * ws_dim1] * neggi;
				pointr = pointr % m + 1;
		    }
		    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], 1);
    }
/*     Initialize GCP xcp = x. */
    dcopy_(n, &x[1], 1, &xcp[1], 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);
		    for (i__ = 1; i__ <= n; ++i__) {
				do_fio(&c__1, (char *)&xcp[i__], (ftnlen)sizeof(doublereal));
		    }
		    e_wsfe();
		}
		return 0;
    }
/*     Initialize c = W'(xcp - x) = 0. */
    for (j = 1; j <= col2; ++j) c__[j] = 0.;

/*     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], 1, &p[1], 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). */
		}
		hpsolb_(nleft, &t[1], &iorder[1], iter - 2);
		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], 1, &c__[1], 1);
	/*           choose wbp, */
	/*           the row of W corresponding to the breakpoint encountered. */
		pointr = head;
		for (j = 1; j <= col; ++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], 1, &v[1], 1);
		wmp = ddot_(col2, &p[1], 1, &v[1], 1);
		wmw = ddot_(col2, &wbp[1], 1, &v[1], 1);
	/*           update p = p - dibp*wbp. */
		daxpy_(col2, -dibp, &wbp[1], 1, &p[1], 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], 1, &xcp[1], 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], 1, &c__[1], 1);
    }
    if (iprint > 100) {
		s_wsfe(&io___131);
		for (i__ = 1; i__ <= n; ++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_ */

/* ====================== The end of cauchy ============================== */
int cmprlb_(const integer n, const integer m, const doublereal *x, 
	const doublereal *g, const doublereal *ws, const doublereal *wy, const doublereal *sy, 
	const doublereal *wt, doublereal *z__, doublereal *r__, doublereal *wa, 
	const integer *index, const doublereal theta, const integer col, const integer head, 
	const integer nfree, const logical cnstnd, integer *info)
{
    /* System generated locals */
    integer ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, 
	    wt_dim1, wt_offset;

    /* Local variables */
    static integer i__, j, k;
    static doublereal a1, a2;
    static integer pointr;

    // Parameter adjustments 
    --index;
    --r__;
    --z__;
    --g;
    --x;
    --wa;
    wt_dim1 = m;
    wt_offset = 1 + wt_dim1;
    wt -= wt_offset;
    sy_dim1 = m;
    sy_offset = 1 + sy_dim1;
    sy -= sy_offset;
    wy_dim1 = n;
    wy_offset = 1 + wy_dim1;
    wy -= wy_offset;
    ws_dim1 = n;
    ws_offset = 1 + ws_dim1;
    ws -= ws_offset;

    /* Function Body */
    if (! (cnstnd) && col > 0) {
		for (i__ = 1; i__ <= n; ++i__) r__[i__] = -g[i__];
    } 
    else {
		for (i__ = 1; i__ <= nfree; ++i__) {
		    k = index[i__];
		    r__[i__] = -theta * (z__[k] - x[k]) - g[k];
		}
		bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wa[(m << 1) + 1], &wa[
			1], info);
		if (*info != 0) {
		    *info = -8;
		    return 0;
		}
		pointr = head;
		for (j = 1; j <= col; ++j) {
		    a1 = wa[j];
		    a2 = theta * wa[col + j];
		    for (i__ = 1; i__ <= nfree; ++i__) {
				k = index[i__];
				r__[i__] = r__[i__] + wy[k + pointr * wy_dim1] * a1 + ws[k + 
					pointr * ws_dim1] * a2;
		    }
		    pointr = pointr % m + 1;
		}
    }
    return 0;
} /* cmprlb_ */

/* ======================= The end of cmprlb ============================= */
int errclb_(const integer n, const integer m, const doublereal factr, 
	const doublereal *l, const doublereal *u, const integer *nbd, 
	char *task, integer *info, integer *k, ftnlen task_len)
{
    /* Builtin functions */
    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);

    /* Local variables */
    static integer i__;

    /* Parameter adjustments */
    --nbd;
    --u;
    --l;

    /* Function Body */
    if (n <= 0) s_copy(task, "ERROR: N .LE. 0", (ftnlen)60, (ftnlen)15);
    if (m <= 0) s_copy(task, "ERROR: M .LE. 0", (ftnlen)60, (ftnlen)15);
    if (factr < 0.) s_copy(task, "ERROR: FACTR .LT. 0", (ftnlen)60, (ftnlen)19);

/*     Check the validity of the arrays nbd(i), u(i), and l(i). */
    for (i__ = 1; i__ <= n; ++i__) {
		if (nbd[i__] < 0 || nbd[i__] > 3) {
	/*                                                   return */
		    s_copy(task, "ERROR: INVALID NBD", (ftnlen)60, (ftnlen)18);
		    *info = -6;
		    *k = i__;
		}
		if (nbd[i__] == 2) {
		    if (l[i__] > u[i__]) {
	/*                                    return */
			s_copy(task, "ERROR: NO FEASIBLE SOLUTION", (ftnlen)60, (
				ftnlen)27);
			*info = -7;
			*k = i__;
		    }
		}
    }
    return 0;
} /* errclb_ */

/* ======================= The end of errclb ============================= */
int formk_(const integer n, const integer nsub, integer *ind, 
	const integer nenter, integer const ileave, integer *indx2, 
	const integer iupdat, const logical updatd, doublereal *wn, 
	doublereal *wn1, const integer m, const doublereal *ws, 
	const doublereal *wy, const doublereal *sy, const doublereal theta, 
	const integer col, const integer head, integer *info)
{
    /* System generated locals */
    integer wn_dim1, wn_offset, wn1_dim1, wn1_offset, ws_dim1, ws_offset, 
	    wy_dim1, wy_offset, sy_dim1, sy_offset, i__1;

    // Local variables 
    static integer i__, k, k1, m2, is, js, iy, jy, is1, js1, col2, dend, pend;
    static integer upcl;
    static doublereal temp1, temp2, temp3, temp4;
    static integer ipntr, jpntr, dbegin, pbegin;

    // Parameter adjustments 
    --indx2;
    --ind;
    sy_dim1 = m;
    sy_offset = 1 + sy_dim1;
    sy -= sy_offset;
    wy_dim1 = n;
    wy_offset = 1 + wy_dim1;
    wy -= wy_offset;
    ws_dim1 = n;
    ws_offset = 1 + ws_dim1;
    ws -= ws_offset;
    wn1_dim1 = 2 * m;
    wn1_offset = 1 + wn1_dim1;
    wn1 -= wn1_offset;
    wn_dim1 = 2 * m;
    wn_offset = 1 + wn_dim1;
    wn -= wn_offset;

    /* Function Body */
    if (updatd) {
		if (iupdat > m) {
	/*                                 shift old part of WN1. */
		    for (jy = 1; jy <= m - 1; ++jy) {
				js = m + jy;
				dcopy_(m - jy, &wn1[jy + 1 + (jy + 1) * wn1_dim1], 1, &wn1[
					jy + jy * wn1_dim1], 1);
				dcopy_(m - jy, &wn1[js + 1 + (js + 1) * wn1_dim1], 1, &wn1[
					js + js * wn1_dim1], 1);
				dcopy_(m - 1, &wn1[m + 2 + (jy + 1) * wn1_dim1], 1, &wn1[
					m + 1 + jy * wn1_dim1], 1);
		    }
		}
	/*          put new rows in blocks (1,1), (2,1) and (2,2). */
		pbegin = 1;
		pend = nsub;
		dbegin = nsub + 1;
		dend = n;
		iy = col;
		is = m + col;
		ipntr = head + col - 1;
		if (ipntr > m) ipntr -= m;
		jpntr = head;
		for (jy = 1; jy <= col; ++jy) {
		    js = m + jy;
		    temp1 = 0.;
		    temp2 = 0.;
		    temp3 = 0.;
	//             compute element jy of row 'col' of Y'ZZ'Y 
		    for (k = pbegin; k <= pend; ++k) {
				k1 = ind[k];
				temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
		    }
	//             compute elements jy of row 'col' of L_a and S'AA'S 
		    for (k = dbegin; k <= dend; ++k) {
				k1 = ind[k];
				temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
				temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
		    }
		    wn1[iy + jy * wn1_dim1] = temp1;
		    wn1[is + js * wn1_dim1] = temp2;
		    wn1[is + jy * wn1_dim1] = temp3;
		    jpntr = jpntr % m + 1;
		}
	/*          put new column in block (2,1). */
		jy = col;
		jpntr = head + col - 1;
		if (jpntr > m) jpntr -= m;
		ipntr = head;
		for (i__ = 1; i__ <= col; ++i__) {
		    is = m + i__;
		    temp3 = 0.;
	/*             compute element i of column 'col' of R_z */
		    for (k = pbegin; k <= pend; ++k) {
				k1 = ind[k];
				temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
		    }
		    ipntr = ipntr % m + 1;
		    wn1[is + jy * wn1_dim1] = temp3;
		}
		upcl = col - 1;
    } 
    else upcl = col;
/*       modify the old parts in blocks (1,1) and (2,2) due to changes 
       in the set of free variables. */
    ipntr = head;
    for (iy = 1; iy <= upcl; ++iy) {
		is = m + iy;
		jpntr = head;
		for (jy = 1; jy <= iy; ++jy) {
		    js = m + jy;
		    temp1 = 0.;
		    temp2 = 0.;
		    temp3 = 0.;
		    temp4 = 0.;
		    for (k = 1; k <= nenter; ++k) {
				k1 = indx2[k];
				temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
				temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
		    }
		    for (k = ileave; k <= n; ++k) {
				k1 = indx2[k];
				temp3 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
				temp4 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
		    }
		    wn1[iy + jy * wn1_dim1] = wn1[iy + jy * wn1_dim1] + temp1 - temp3;
		    wn1[is + js * wn1_dim1] = wn1[is + js * wn1_dim1] - temp2 + temp4;
		    jpntr = jpntr % m + 1;
		}
		ipntr = ipntr % m + 1;
    }
/*       modify the old parts in block (2,1). */
    ipntr = head;
    i__1 = m + upcl;
    for (is = m + 1; is <= i__1; ++is) {
		jpntr = head;
		for (jy = 1; jy <= upcl; ++jy) {
		    temp1 = 0.;
		    temp3 = 0.;
		    for (k = 1; k <= nenter; ++k) {
				k1 = indx2[k];
				temp1 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
		    }
		    for (k = ileave; k <= n; ++k) {
				k1 = indx2[k];
				temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
		    }
		    if (is <= jy + m) {
				wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] + temp1 - 
					temp3;
		    } 
		    else {
				wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] - temp1 + 
					temp3;
		    }
		    jpntr = jpntr % m + 1;
		}
		ipntr = ipntr % m + 1;
    }
/*     Form the upper triangle of WN = [D+Y' ZZ'Y/theta   -L_a'+R_z' ] 
                                       [-L_a +R_z        S'AA'S*theta] */
    m2 = m << 1;
    for (iy = 1; iy <= col; ++iy) {
		is = col + iy;
		is1 = m + iy;
		for (jy = 1; jy <= iy; ++jy) {
		    js = col + jy;
		    js1 = m + jy;
		    wn[jy + iy * wn_dim1] = wn1[iy + jy * wn1_dim1] / theta;
		    wn[js + is * wn_dim1] = wn1[is1 + js1 * wn1_dim1] * theta;
		}
		for (jy = 1; jy <= iy - 1; ++jy) {
		    wn[jy + is * wn_dim1] = -wn1[is1 + jy * wn1_dim1];
		}
		for (jy = iy; jy <= col; ++jy) {
		    wn[jy + is * wn_dim1] = wn1[is1 + jy * wn1_dim1];
		}
		wn[iy + iy * wn_dim1] += sy[iy + iy * sy_dim1];
    }
/*     Form the upper triangle of WN= [  LL'            L^-1(-L_a'+R_z')] 
                                      [(-L_a +R_z)L'^-1   S'AA'S*theta  ] 
        first Cholesky factor (1,1) block of wn to get LL' 
                          with L' stored in the upper triangle of wn. */
    dpofa_(&wn[wn_offset], m2, col, info);
    if (*info != 0) {
		*info = -1;
		return 0;
    }
/*        then form L^-1(-L_a'+R_z') in the (1,2) block. */
    col2 = col << 1;
    for (js = col + 1; js <= col2; ++js) {
		dtrsl_(&wn[wn_offset], m2, col, &wn[js * wn_dim1 + 1], 11, info);
    }
/*     Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the 
        upper triangle of (2,2) block of wn. */
    for (is = col + 1; is <= col2; ++is) {
		for (js = is; js <= col2; ++js) {
		    wn[is + js * wn_dim1] += ddot_(col, &wn[is * wn_dim1 + 1], 1, 
			    &wn[js * wn_dim1 + 1], 1);
		}
    }
/*     Cholesky factorization of (2,2) block of wn. */
    dpofa_(&wn[col + 1 + (col + 1) * wn_dim1], m2, col, info);
    if (*info != 0) {
		*info = -2;
		return 0;
    }
    return 0;
} /* formk_ */

/* ======================= The end of formk ============================== */
int formt_(const integer m, doublereal *wt, const doublereal *sy, 
	const doublereal *ss, const integer col, const doublereal theta, integer *info)
{
    /* System generated locals */
    integer wt_dim1, wt_offset, sy_dim1, sy_offset, ss_dim1, ss_offset;

    /* Local variables */
    static integer i__, j, k, k1;
    static doublereal ddum;


	// Parameter adjustments 
    ss_dim1 = m;
    ss_offset = 1 + ss_dim1;
    ss -= ss_offset;
    sy_dim1 = m;
    sy_offset = 1 + sy_dim1;
    sy -= sy_offset;
    wt_dim1 = m;
    wt_offset = 1 + wt_dim1;
    wt -= wt_offset;

    // Function Body 
    for (j = 1; j <= col; ++j) {
		wt[j * wt_dim1 + 1] = theta * ss[j * ss_dim1 + 1];
    }
    for (i__ = 2; i__ <= col; ++i__) {
		for (j = i__; j <= col; ++j) {
		    k1 = min(i__,j) - 1;
		    ddum = 0.;
		    for (k = 1; k <= k1; ++k) {
				ddum += sy[i__ + k * sy_dim1] * sy[j + k * sy_dim1] / sy[k + 
					k * sy_dim1];
		    }
		    wt[i__ + j * wt_dim1] = ddum + theta * ss[i__ + j * ss_dim1];
		}
    }
/*     Cholesky factorize T to J*J' with 
        J' stored in the upper triangle of wt. */
    dpofa_(&wt[wt_offset], m, col, info);
    if (*info != 0) *info = -3;
    return 0;
} /* formt_ */

/* ======================= The end of formt ============================== */
int freev_(const integer n, integer *nfree, integer *index, 
	integer *nenter, integer *ileave, integer *indx2, integer *iwhere, 
	logical *wrk, const logical updatd, const logical cnstnd, const integer iprint, 
	const integer iter)
{
    /* System generated locals */
    integer i__1;

    /* Builtin functions */
    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_wsle();

    /* Local variables */
    static integer i__, k, iact;

    /* Fortran I/O blocks */
    static cilist io___169 = { 0, 6, 0, 0, 0 };
    static cilist io___170 = { 0, 6, 0, 0, 0 };
    static cilist io___171 = { 0, 6, 0, 0, 0 };
    static cilist io___173 = { 0, 6, 0, 0, 0 };


    /* Parameter adjustments */
    --iwhere;
    --indx2;
    --index;

    /* Function Body */
    *nenter = 0;
    *ileave = n + 1;
    if (iter > 0 && cnstnd) {
/*                           count the entering and leaving variables. */
	for (i__ = 1; i__ <= *nfree; ++i__) {
	    k = index[i__];
	    if (iwhere[k] > 0) {
			--(*ileave);
			indx2[*ileave] = k;
			if (iprint >= 100) {
			    s_wsle(&io___169);
			    do_lio(&c__9, &c__1, "Variable ", (ftnlen)9);
			    do_lio(&c__3, &c__1, (char *)&k, (ftnlen)sizeof(integer));
			    do_lio(&c__9, &c__1, " leaves the set of free variables", 
				    (ftnlen)33);
			    e_wsle();
			}
	    }
	}
	for (i__ = *nfree + 1; i__ <= n; ++i__) {
	    k = index[i__];
	    if (iwhere[k] <= 0) {
			++(*nenter);
			indx2[*nenter] = k;
			if (iprint >= 100) {
			    s_wsle(&io___170);
			    do_lio(&c__9, &c__1, "Variable ", (ftnlen)9);
			    do_lio(&c__3, &c__1, (char *)&k, (ftnlen)sizeof(integer));
			    do_lio(&c__9, &c__1, " enters the set of free variables", 
				    (ftnlen)33);
			    e_wsle();
			}
	    }
	}
	if (iprint >= 99) {
	    s_wsle(&io___171);
	    i__1 = n + 1 - *ileave;
	    do_lio(&c__3, &c__1, (char *)&i__1, (ftnlen)sizeof(integer));
	    do_lio(&c__9, &c__1, " variables leave; ", (ftnlen)18);
	    do_lio(&c__3, &c__1, (char *)&(*nenter), (ftnlen)sizeof(integer));
	    do_lio(&c__9, &c__1, " variables enter", (ftnlen)16);
	    e_wsle();
	}
    }
    *wrk = *ileave < n + 1 || *nenter > 0 || updatd;
/*     Find the index set of free and active variables at the GCP. */
    *nfree = 0;
    iact = n + 1;
    for (i__ = 1; i__ <= n; ++i__) {
		if (iwhere[i__] <= 0) {
		    ++(*nfree);
		    index[*nfree] = i__;
		} else {
		    --iact;
		    index[iact] = i__;
		}
    }
    if (iprint >= 99) {
		s_wsle(&io___173);
		do_lio(&c__3, &c__1, (char *)&(*nfree), (ftnlen)sizeof(integer));
		do_lio(&c__9, &c__1, " variables are free at GCP ", (ftnlen)27);
		i__1 = iter + 1;
		do_lio(&c__3, &c__1, (char *)&i__1, (ftnlen)sizeof(integer));
		e_wsle();
    }
    return 0;
} /* freev_ */

/* ======================= The end of freev ============================== */
int hpsolb_(const integer n, doublereal *t, integer *iorder, 
	const integer iheap)
{
    /* Local variables */
    static integer i__, j, k;
    static doublereal out, ddum;
    static integer indxin, indxou;

    /* Parameter adjustments */
    --iorder;
    --t;

    /* Function Body */
    if (iheap == 0) {
	/*        Rearrange the elements t(1) to t(n) to form a heap. */
		for (k = 2; k <= n; ++k) {
		    ddum = t[k];
		    indxin = iorder[k];
	/*           Add ddum to the heap. */
		    i__ = k;
	L10:
		    if (i__ > 1) {
				j = i__ / 2;
				if (ddum < t[j]) {
				    t[i__] = t[j];
				    iorder[i__] = iorder[j];
				    i__ = j;
				    goto L10;
				}
		    }
		    t[i__] = ddum;
		    iorder[i__] = indxin;
		}
    }
/*     Assign to 'out' the value of t(1), the least member of the heap, 
        and rearrange the remaining members to form a heap as
        elements 1 to n-1 of t. */
    if (n > 1) {
		i__ = 1;
		out = t[1];
		indxou = iorder[1];
		ddum = t[n];
		indxin = iorder[n];
	/*        Restore the heap */
	L30:
		j = i__ + i__;
		if (j <= n - 1) {
		    if (t[j + 1] < t[j]) ++j;
		    if (t[j] < ddum) {
				t[i__] = t[j];
				iorder[i__] = iorder[j];
				i__ = j;
				goto L30;
		    }
		}
		t[i__] = ddum;
		iorder[i__] = indxin;
	//     Put the least member in t(n). 
		t[n] = out;
		iorder[n] = indxou;
    }
    return 0;
} /* hpsolb_ */

/* ====================== The end of hpsolb ============================== */
int lnsrlb_(const integer n, const doublereal *l, const doublereal *u, 
	const integer *nbd, doublereal *x, doublereal *f, doublereal *fold, 
	doublereal *gd, doublereal *gdold, doublereal *g, doublereal *d__, 
	doublereal *r__, doublereal *t, doublereal *z__, doublereal *stp, 
	doublereal *dnorm, doublereal *dtd, doublereal *xstep, doublereal *
	stpmx, const integer iter, integer *ifun, integer *iback, integer *nfgv, 
	integer *info, char *task, const logical boxed, const logical cnstnd, char *
	csave, integer *isave, doublereal *dsave, ftnlen task_len, ftnlen 
	csave_len)
{
    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    integer s_cmp(char *, char *, ftnlen, ftnlen);
    double sqrt(doublereal);
    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);

    /* Local variables */
    static integer i__;
    static doublereal a1, a2;

    /* Parameter adjustments */
    --z__;
    --t;
    --r__;
    --d__;
    --g;
    --x;
    --nbd;
    --u;
    --l;
    --isave;
    --dsave;

    /* Function Body */
    if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) goto L556;

    *dtd = ddot_(n, &d__[1], 1, &d__[1], 1);
    *dnorm = sqrt(*dtd);
/*     Determine the maximum step length. */
    *stpmx = 1e10;
    if (cnstnd) {
		if (iter == 0) *stpmx = 1.;
		else {
		    for (i__ = 1; i__ <= n; ++i__) {
				a1 = d__[i__];
				if (nbd[i__] != 0) {
				    if (a1 < 0. && nbd[i__] <= 2) {
						a2 = l[i__] - x[i__];
						if (a2 >= 0.) *stpmx = 0.;
						else if (a1 * *stpmx < a2) *stpmx = a2 / a1;
				    }
				    else if (a1 > 0. && nbd[i__] >= 2) {
						a2 = u[i__] - x[i__];
						if (a2 <= 0.) *stpmx = 0.;
						else if (a1 * *stpmx > a2) *stpmx = a2 / a1;
				    }
				}
		    }
		}
    }
    if (iter == 0 && ! (boxed)) {
	/* Computing MIN */
		d__1 = 1. / *dnorm;
		*stp = min(d__1,*stpmx);
    } 
    else *stp = 1.;
    dcopy_(n, &x[1], 1, &t[1], 1);
    dcopy_(n, &g[1], 1, &r__[1], 1);
    *fold = *f;
    *ifun = 0;
    *iback = 0;
    s_copy(csave, "START", (ftnlen)60, (ftnlen)5);
L556:
    *gd = ddot_(n, &g[1], 1, &d__[1], 1);
    if (*ifun == 0) {
		*gdold = *gd;
		if (*gd >= 0.) {
	/*                               the directional derivative >=0. 
	                               Line search is impossible. */
		    *info = -4;
		    return 0;
		}
    }
    dcsrch_(f, gd, stp, .001, .9, .1, 0., *stpmx, csave, &
	    isave[1], &dsave[1], (ftnlen)60);
    *xstep = *stp * *dnorm;
    if (s_cmp(csave, "CONV", (ftnlen)4, (ftnlen)4) != 0 && s_cmp(csave, "WARN"
	    , (ftnlen)4, (ftnlen)4) != 0) {
		s_copy(task, "FG_LNSRCH", (ftnlen)60, (ftnlen)9);
		++(*ifun);
		++(*nfgv);
		*iback = *ifun - 1;
		if (*stp == 1.) {
		    dcopy_(n, &z__[1], 1, &x[1], 1);
		} 
		else {
		    for (i__ = 1; i__ <= n; ++i__) {
				x[i__] = *stp * d__[i__] + t[i__];
		    }
		}
    } 
    else {
		s_copy(task, "NEW_X", (ftnlen)60, (ftnlen)5);
    }
    return 0;
} /* lnsrlb_ */

/* ======================= The end of lnsrlb ============================= */
int matupd_(const integer n, const integer m, doublereal *ws, 
	doublereal *wy, doublereal *sy, doublereal *ss, doublereal *d__, 
	doublereal *r__, integer *itail, const integer iupdat, integer *col, 
	integer *head, doublereal *theta, const doublereal rr, const doublereal dr, 
	const doublereal stp, const doublereal dtd)
{
    /* System generated locals */
    integer ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, 
	    ss_dim1, ss_offset;

    /* Local variables */
    static integer j;
    static integer pointr;

    /* Parameter adjustments */
    --r__;
    --d__;
    ss_dim1 = m;
    ss_offset = 1 + ss_dim1;
    ss -= ss_offset;
    sy_dim1 = m;
    sy_offset = 1 + sy_dim1;
    sy -= sy_offset;
    wy_dim1 = n;
    wy_offset = 1 + wy_dim1;
    wy -= wy_offset;
    ws_dim1 = n;
    ws_offset = 1 + ws_dim1;
    ws -= ws_offset;

    // Function Body 
    if (iupdat <= m) {
		*col = iupdat;
		*itail = (*head + iupdat - 2) % m + 1;
    } 
    else {
		*itail = *itail % m + 1;
		*head = *head % m + 1;
    }
/*     Update matrices WS and WY. */
    dcopy_(n, &d__[1], 1, &ws[*itail * ws_dim1 + 1], 1);
    dcopy_(n, &r__[1], 1, &wy[*itail * wy_dim1 + 1], 1);
/*     Set theta=yy/ys. */
    *theta = rr / dr;
/*     Form the middle matrix in B. */
/*        update the upper triangle of SS, */
/*                                         and the lower triangle of SY: */
    if (iupdat > m) {
	/*                              move old information */
		for (j = 1; j <= *col - 1; ++j) {
		    dcopy_(j, &ss[(j + 1) * ss_dim1 + 2], 1, &ss[j * ss_dim1 + 1]
			    , 1);
		    dcopy_(*col - j, &sy[j + 1 + (j + 1) * sy_dim1], 1, &sy[j + j * 
			    sy_dim1], 1);
		}
    }
/*        add new information: the last row of SY */
/*                                             and the last column of SS: */
    pointr = *head;
    for (j = 1; j <= *col - 1; ++j) {
		sy[*col + j * sy_dim1] = ddot_(n, &d__[1], 1, &wy[pointr * 
			wy_dim1 + 1], 1);
		ss[j + *col * ss_dim1] = ddot_(n, &ws[pointr * ws_dim1 + 1], 1, &
			d__[1], 1);
		pointr = pointr % m + 1;
    }
    if (stp == 1.) ss[*col + *col * ss_dim1] = dtd;
    else ss[*col + *col * ss_dim1] = stp * stp * dtd;
    sy[*col + *col * sy_dim1] = dr;
    return 0;
} /* matupd_ */

/* ======================= The end of matupd ============================= */
int prn1lb_(const integer n, const integer m, const doublereal *l, 
	const doublereal *u, const doublereal *x, 
	const integer iprint, const integer *itfile, const doublereal epsmch)
{
    /* Format strings */
    static char fmt_7001[] = "(\002RUNNING THE L-BFGS-B CODE\002,/,/,\002   \
        * * *\002,/,/,\002Machine precision =\002,1p,d10.3)";
    static char fmt_2001[] = "(\002RUNNING THE L-BFGS-B CODE\002,/,/,\002it \
   = iteration number\002,/,\002nf    = number of function evaluations\002,/,\
\002nint  = number of segments explored during the Cauchy search\002,/,\002n\
act  = number of active bounds at the generalized Cauchy point\002,/,\002sub\
   = manner in which the subspace minimization terminated:\002,/,\002       \
 con = converged, bnd = a bound was reached\002,/,\002itls  = number of iter\
ations performed in the line search\002,/,\002stepl = step length used\002,/,\
\002tstep = norm of the displacement (total step)\002,/,\002projg = norm of \
the projected gradient\002,/,\002f     = function value\002,/,/,\002        \
   * * *\002,/,/,\002Machine precision =\002,1p,d10.3)";
    static char fmt_9001[] = "(/,3x,\002it\002,3x,\002nf\002,2x,\002nint\002\
,2x,\002nact\002,2x,\002sub\002,2x,\002itls\002,2x,\002stepl\002,4x,\002tstep\
\002,5x,\002projg\002,8x,\002f\002)";
    static char fmt_1004[] = "(/,a4,1p,6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))";

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(), 
	    s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_wsle();

    /* Local variables */
    static integer i__;

    /* Fortran I/O blocks */
    static cilist io___186 = { 0, 6, 0, fmt_7001, 0 };
    static cilist io___187 = { 0, 6, 0, 0, 0 };
    static cilist io___188 = { 0, 0, 0, fmt_2001, 0 };
    static cilist io___189 = { 0, 0, 0, 0, 0 };
    static cilist io___190 = { 0, 0, 0, fmt_9001, 0 };
    static cilist io___191 = { 0, 6, 0, fmt_1004, 0 };
    static cilist io___193 = { 0, 6, 0, fmt_1004, 0 };
    static cilist io___194 = { 0, 6, 0, fmt_1004, 0 };

    /* Parameter adjustments */
    --x;
    --u;
    --l;

    /* Function Body */
    if (iprint >= 0) {
		s_wsfe(&io___186);
		do_fio(&c__1, (char *)&epsmch, (ftnlen)sizeof(doublereal));
		e_wsfe();
		s_wsle(&io___187);
		do_lio(&c__9, &c__1, "N = ", (ftnlen)4);
		do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
		do_lio(&c__9, &c__1, "    M = ", (ftnlen)8);
		do_lio(&c__3, &c__1, (char *)&m, (ftnlen)sizeof(integer));
		e_wsle();
		if (iprint >= 1) {
		    io___188.ciunit = *itfile;
		    s_wsfe(&io___188);
		    do_fio(&c__1, (char *)&epsmch, (ftnlen)sizeof(doublereal));
		    e_wsfe();
		    io___189.ciunit = *itfile;
		    s_wsle(&io___189);
		    do_lio(&c__9, &c__1, "N = ", (ftnlen)4);
		    do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
		    do_lio(&c__9, &c__1, "    M = ", (ftnlen)8);
		    do_lio(&c__3, &c__1, (char *)&m, (ftnlen)sizeof(integer));
		    e_wsle();
		    io___190.ciunit = *itfile;
		    s_wsfe(&io___190);
		    e_wsfe();
		    if (iprint > 100) {
				s_wsfe(&io___191);
				do_fio(&c__1, "L =", (ftnlen)3);
				for (i__ = 1; i__ <= n; ++i__) {
				    do_fio(&c__1, (char *)&l[i__], (ftnlen)sizeof(doublereal));
				}
				e_wsfe();
				s_wsfe(&io___193);
				do_fio(&c__1, "X0 =", (ftnlen)4);
				for (i__ = 1; i__ <= n; ++i__) {
				    do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
				}
				e_wsfe();
				s_wsfe(&io___194);
				do_fio(&c__1, "U =", (ftnlen)3);
				for (i__ = 1; i__ <= n; ++i__) {
				    do_fio(&c__1, (char *)&u[i__], (ftnlen)sizeof(doublereal));
				}
				e_wsfe();
		    }
		}
    }
    return 0;
} /* prn1lb_ */

/* ======================= The end of prn1lb ============================= */
int prn2lb_(const integer n, const doublereal *x, const doublereal f, 
	const doublereal *g, const integer iprint, integer *itfile, const integer iter, 
	const integer nfgv, const integer nact, const doublereal sbgnrm, 
	const integer nint, char *word, const integer iword, const integer iback, 
	const doublereal stp, const doublereal xstep, ftnlen word_len)
{
    /* Format strings */
    static char fmt_2001[] = "(/,\002At iterate\002,i5,4x,\002f= \002,1p,d12\
.5,4x,\002|proj g|= \002,1p,d12.5)";
    static char fmt_1004[] = "(/,a4,1p,6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))";
    static char fmt_3001[] = "(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),1\
p,2(1x,d10.3))";

    /* Builtin functions */
    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_wsle(), s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), 
	    e_wsfe();

    /* Local variables */
    static integer i__, imod;

    /* Fortran I/O blocks */
    static cilist io___195 = { 0, 6, 0, 0, 0 };
    static cilist io___196 = { 0, 6, 0, fmt_2001, 0 };
    static cilist io___197 = { 0, 6, 0, fmt_1004, 0 };
    static cilist io___199 = { 0, 6, 0, fmt_1004, 0 };
    static cilist io___201 = { 0, 6, 0, fmt_2001, 0 };
    static cilist io___202 = { 0, 0, 0, fmt_3001, 0 };

    // Parameter adjustments 
    --g;
    --x;

    // Function Body 
    if (iword == 0) {
/*                            the subspace minimization converged. */
		s_copy(word, "con", (ftnlen)3, (ftnlen)3);
    } 
    else if (iword == 1) {
/*                          the subspace minimization stopped at a bound. */
		s_copy(word, "bnd", (ftnlen)3, (ftnlen)3);
    } 
    else if (iword == 5) {
/*                             the truncated Newton step has been used. */
		s_copy(word, "TNT", (ftnlen)3, (ftnlen)3);
    } 
    else {
		s_copy(word, "---", (ftnlen)3, (ftnlen)3);
    }
    if (iprint >= 99) {
		s_wsle(&io___195);
		do_lio(&c__9, &c__1, "LINE SEARCH", (ftnlen)11);
		do_lio(&c__3, &c__1, (char *)&iback, (ftnlen)sizeof(integer));
		do_lio(&c__9, &c__1, " times; norm of step = ", (ftnlen)23);
		do_lio(&c__5, &c__1, (char *)&xstep, (ftnlen)sizeof(doublereal));
		e_wsle();
		s_wsfe(&io___196);
		do_fio(&c__1, (char *)&iter, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&f, (ftnlen)sizeof(doublereal));
		do_fio(&c__1, (char *)&sbgnrm, (ftnlen)sizeof(doublereal));
		e_wsfe();
		if (iprint > 100) {
		    s_wsfe(&io___197);
		    do_fio(&c__1, "X =", (ftnlen)3);
		    for (i__ = 1; i__ <= n; ++i__) {
				do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
		    }
		    e_wsfe();
		    s_wsfe(&io___199);
		    do_fio(&c__1, "G =", (ftnlen)3);
		    for (i__ = 1; i__ <= n; ++i__) {
				do_fio(&c__1, (char *)&g[i__], (ftnlen)sizeof(doublereal));
		    }
		    e_wsfe();
		}
    } 
    else if (iprint > 0) {
		imod = iter % iprint;
		if (imod == 0) {
		    s_wsfe(&io___201);
		    do_fio(&c__1, (char *)&iter, (ftnlen)sizeof(integer));
		    do_fio(&c__1, (char *)&f, (ftnlen)sizeof(doublereal));
		    do_fio(&c__1, (char *)&sbgnrm, (ftnlen)sizeof(doublereal));
		    e_wsfe();
		}
    }
    if (iprint >= 1) {
		io___202.ciunit = *itfile;
		s_wsfe(&io___202);
		do_fio(&c__1, (char *)&iter, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&nfgv, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&nint, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&nact, (ftnlen)sizeof(integer));
		do_fio(&c__1, word, (ftnlen)3);
		do_fio(&c__1, (char *)&iback, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&stp, (ftnlen)sizeof(doublereal));
		do_fio(&c__1, (char *)&xstep, (ftnlen)sizeof(doublereal));
		do_fio(&c__1, (char *)&sbgnrm, (ftnlen)sizeof(doublereal));
		do_fio(&c__1, (char *)&f, (ftnlen)sizeof(doublereal));
		e_wsfe();
    }
    return 0;
} /* prn2lb_ */

/* ======================= The end of prn2lb ============================= */
int prn3lb_(const integer n, const doublereal *x, const doublereal f, char *
	task, const integer iprint, const integer info, integer *itfile, const integer iter, 
	const integer nfgv, const integer nintol, const integer nskip, const integer nact, 
	const doublereal sbgnrm, const doublereal time, const integer nint, char *word, 
	const integer iback, const doublereal stp, const doublereal xstep, const integer k, 
	const doublereal cachyt, const doublereal sbtime, const doublereal lnscht, ftnlen 
	task_len, ftnlen word_len)
{
    /* Format strings */
    static char fmt_3003[] = "(/,\002           * * *\002,/,/,\002Tit   = to\
tal number of iterations\002,/,\002Tnf   = total number of function evaluati\
ons\002,/,\002Tnint = total number of segments explored during\002,\002 Cauc\
hy searches\002,/,\002Skip  = number of BFGS updates skipped\002,/,\002Nact \
 = number of active bounds at final generalized\002,\002 Cauchy point\002,/\
,\002Projg = norm of the final projected gradient\002,/,\002F     = final fu\
nction value\002,/,/,\002           * * *\002)";
    static char fmt_3004[] = "(/,3x,\002N\002,3x,\002Tit\002,2x,\002Tnf\002,\
2x,\002Tnint\002,2x,\002Skip\002,2x,\002Nact\002,5x,\002Projg\002,8x,\002\
F\002)";
    static char fmt_3005[] = "(i5,2(1x,i4),(1x,i6),(2x,i4),(1x,i5),1p,2(2x,d\
10.3))";
    static char fmt_1004[] = "(/,a4,1p,6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))";
    static char fmt_3009[] = "(/,a60)";
    static char fmt_9011[] = "(/,\002 Matrix in 1st Cholesky factorization i\
n formk is not Pos. Def.\002)";
    static char fmt_9012[] = "(/,\002 Matrix in 2st Cholesky factorization i\
n formk is not Pos. Def.\002)";
    static char fmt_9013[] = "(/,\002 Matrix in the Cholesky factorization i\
n formt is not Pos. Def.\002)";
    static char fmt_9014[] = "(/,\002 Derivative >= 0, backtracking line sea\
rch impossible.\002,/,\002   Previous x, f and g restored.\002,/,\002 Possib\
le causes: 1 error in function or gradient evaluation;\002,/,\002           \
       2 rounding errors dominate computation.\002)";
    static char fmt_9015[] = "(/,\002 Warning:  more than 10 function and gr\
adient\002,/,\002   evaluations in the last line search.  Termination\002,/\
,\002   may possibly be caused by a bad search direction.\002)";
    static char fmt_9018[] = "(/,\002 The triangular system is singular.\002)"
	    ;
    static char fmt_9019[] = "(/,\002 Line search cannot locate an adequate \
point after 20 function\002,/,\002  and gradient evaluations.  Previous x, f\
 and g restored.\002,/,\002 Possible causes: 1 error in function or gradient\
 evaluation;\002,/,\002                  2 rounding error dominate computati\
on.\002)";
    static char fmt_3007[] = "(/,\002 Cauchy                time\002,1p,e10.\
3,\002 seconds.\002,/\002 Subspace minimization time\002,1p,e10.3,\002 secon\
ds.\002,/\002 Line search           time\002,1p,e10.3,\002 seconds.\002)";
    static char fmt_3008[] = "(/,\002 Total User time\002,1p,e10.3,\002 seco\
nds.\002,/)";
    static char fmt_3002[] = "(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),6\
x,\002-\002,10x,\002-\002)";

    /* Builtin functions */
    integer s_cmp(char *, char *, ftnlen, ftnlen), s_wsfe(cilist *), e_wsfe(),
	     do_fio(integer *, char *, ftnlen), s_wsle(cilist *), do_lio(
	    integer *, integer *, char *, ftnlen), e_wsle();

    /* Local variables */
    static integer i__;

    /* Fortran I/O blocks */
    static cilist io___203 = { 0, 6, 0, fmt_3003, 0 };
    static cilist io___204 = { 0, 6, 0, fmt_3004, 0 };
    static cilist io___205 = { 0, 6, 0, fmt_3005, 0 };
    static cilist io___206 = { 0, 6, 0, fmt_1004, 0 };
    static cilist io___208 = { 0, 6, 0, 0, 0 };
    static cilist io___209 = { 0, 6, 0, fmt_3009, 0 };
    static cilist io___210 = { 0, 6, 0, fmt_9011, 0 };
    static cilist io___211 = { 0, 6, 0, fmt_9012, 0 };
    static cilist io___212 = { 0, 6, 0, fmt_9013, 0 };
    static cilist io___213 = { 0, 6, 0, fmt_9014, 0 };
    static cilist io___214 = { 0, 6, 0, fmt_9015, 0 };
    static cilist io___215 = { 0, 6, 0, 0, 0 };
    static cilist io___216 = { 0, 6, 0, 0, 0 };
    static cilist io___217 = { 0, 6, 0, fmt_9018, 0 };
    static cilist io___218 = { 0, 6, 0, fmt_9019, 0 };
    static cilist io___219 = { 0, 6, 0, fmt_3007, 0 };
    static cilist io___220 = { 0, 6, 0, fmt_3008, 0 };
    static cilist io___221 = { 0, 0, 0, fmt_3002, 0 };
    static cilist io___222 = { 0, 0, 0, fmt_3009, 0 };
    static cilist io___223 = { 0, 0, 0, fmt_9011, 0 };
    static cilist io___224 = { 0, 0, 0, fmt_9012, 0 };
    static cilist io___225 = { 0, 0, 0, fmt_9013, 0 };
    static cilist io___226 = { 0, 0, 0, fmt_9014, 0 };
    static cilist io___227 = { 0, 0, 0, fmt_9015, 0 };
    static cilist io___228 = { 0, 0, 0, fmt_9018, 0 };
    static cilist io___229 = { 0, 0, 0, fmt_9019, 0 };
    static cilist io___230 = { 0, 0, 0, fmt_3008, 0 };

    /* Parameter adjustments */
    --x;

    /* Function Body */
    if (s_cmp(task, "ERROR", (ftnlen)5, (ftnlen)5) == 0) goto L999;
    if (iprint >= 0) {
		s_wsfe(&io___203);
		e_wsfe();
		s_wsfe(&io___204);
		e_wsfe();
		s_wsfe(&io___205);
		do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&iter, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&nfgv, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&nintol, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&nskip, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&nact, (ftnlen)sizeof(integer));
		do_fio(&c__1, (char *)&sbgnrm, (ftnlen)sizeof(doublereal));
		do_fio(&c__1, (char *)&f, (ftnlen)sizeof(doublereal));
		e_wsfe();
		if (iprint >= 100) {
		    s_wsfe(&io___206);
		    do_fio(&c__1, "X =", (ftnlen)3);
		    for (i__ = 1; i__ <= n; ++i__) {
				do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
		    }
		    e_wsfe();
		}
		if (iprint >= 1) {
		    s_wsle(&io___208);
		    do_lio(&c__9, &c__1, " F =", (ftnlen)4);
		    do_lio(&c__5, &c__1, (char *)&f, (ftnlen)sizeof(doublereal));
		    e_wsle();
		}
    }
L999:
    if (iprint >= 0) {
		s_wsfe(&io___209);
		do_fio(&c__1, task, (ftnlen)60);
		e_wsfe();
		if (info != 0) {
		    if (info == -1) {
				s_wsfe(&io___210);
				e_wsfe();
		    }
		    else if (info == -2) {
				s_wsfe(&io___211);
				e_wsfe();
		    }
		    else if (info == -3) {
				s_wsfe(&io___212);
				e_wsfe();
		    }
		    else if (info == -4) {
				s_wsfe(&io___213);
				e_wsfe();
		    }
		    else if (info == -5) {
				s_wsfe(&io___214);
				e_wsfe();
		    }
		    else if (info == -6) {
				s_wsle(&io___215);
				do_lio(&c__9, &c__1, " Input nbd(", (ftnlen)11);
				do_lio(&c__3, &c__1, (char *)&k, (ftnlen)sizeof(integer));
				do_lio(&c__9, &c__1, ") is invalid.", (ftnlen)13);
				e_wsle();
		    }
		    else if (info == -7) {
				s_wsle(&io___216);
				do_lio(&c__9, &c__1, " l(", (ftnlen)3);
				do_lio(&c__3, &c__1, (char *)&k, (ftnlen)sizeof(integer));
				do_lio(&c__9, &c__1, ") > u(", (ftnlen)6);
				do_lio(&c__3, &c__1, (char *)&k, (ftnlen)sizeof(integer));
				do_lio(&c__9, &c__1, ").  No feasible solution.", (ftnlen)25);
				e_wsle();
		    }
		    else if (info == -8) {
				s_wsfe(&io___217);
				e_wsfe();
		    }
		    else if (info == -9) {
				s_wsfe(&io___218);
				e_wsfe();
		    }
		}
		if (iprint >= 1) {
		    s_wsfe(&io___219);
		    do_fio(&c__1, (char *)&cachyt, (ftnlen)sizeof(doublereal));
		    do_fio(&c__1, (char *)&sbtime, (ftnlen)sizeof(doublereal));
		    do_fio(&c__1, (char *)&lnscht, (ftnlen)sizeof(doublereal));
		    e_wsfe();
		}
		s_wsfe(&io___220);
		do_fio(&c__1, (char *)&time, (ftnlen)sizeof(doublereal));
		e_wsfe();
		if (iprint >= 1) {
		    if (info == -4 || info == -9) {
				io___221.ciunit = *itfile;
				s_wsfe(&io___221);
				do_fio(&c__1, (char *)&iter, (ftnlen)sizeof(integer));
				do_fio(&c__1, (char *)&nfgv, (ftnlen)sizeof(integer));
				do_fio(&c__1, (char *)&nint, (ftnlen)sizeof(integer));
				do_fio(&c__1, (char *)&nact, (ftnlen)sizeof(integer));
				do_fio(&c__1, word, (ftnlen)3);
				do_fio(&c__1, (char *)&iback, (ftnlen)sizeof(integer));
				do_fio(&c__1, (char *)&stp, (ftnlen)sizeof(doublereal));
				do_fio(&c__1, (char *)&xstep, (ftnlen)sizeof(doublereal));
				e_wsfe();
		    }
		    io___222.ciunit = *itfile;
		    s_wsfe(&io___222);
		    do_fio(&c__1, task, (ftnlen)60);
		    e_wsfe();
		    if (info != 0) {
				if (info == -1) {
				    io___223.ciunit = *itfile;
				    s_wsfe(&io___223);
				    e_wsfe();
				}
				else if (info == -2) {
				    io___224.ciunit = *itfile;
				    s_wsfe(&io___224);
				    e_wsfe();
				}
				else if (info == -3) {
				    io___225.ciunit = *itfile;
				    s_wsfe(&io___225);
				    e_wsfe();
				}
				else if (info == -4) {
				    io___226.ciunit = *itfile;
				    s_wsfe(&io___226);
				    e_wsfe();
				}
				else if (info == -5) {
				    io___227.ciunit = *itfile;
				    s_wsfe(&io___227);
				    e_wsfe();
				}
				else if (info == -8) {
				    io___228.ciunit = *itfile;
				    s_wsfe(&io___228);
				    e_wsfe();
				}
				else if (info == -9) {
				    io___229.ciunit = *itfile;
				    s_wsfe(&io___229);
				    e_wsfe();
				}
		    }
		    io___230.ciunit = *itfile;
		    s_wsfe(&io___230);
		    do_fio(&c__1, (char *)&time, (ftnlen)sizeof(doublereal));
		    e_wsfe();
		}
    }
    return 0;
} /* prn3lb_ */

/* ======================= The end of prn3lb ============================= */
int projgr_(const integer n, const doublereal *l, const doublereal *u, 
	const integer *nbd, const doublereal *x, const doublereal *g, doublereal *sbgnrm)
{
    /* System generated locals */
    doublereal d__1, d__2;

    /* Local variables */
    static integer i__;
    static doublereal gi;

    // Parameter adjustments
    --g;
    --x;
    --nbd;
    --u;
    --l;

    /* Function Body */
    *sbgnrm = 0.;
    for (i__ = 1; i__ <= n; ++i__) {
		gi = g[i__];
		if (nbd[i__] != 0) {
		    if (gi < 0.) {
				if (nbd[i__] >= 2) {
		/* Computing MAX */
				    d__1 = x[i__] - u[i__];
				    gi = max(d__1,gi);
				}
		    } 
		    else {
				if (nbd[i__] <= 2) {
		/* Computing MIN */
				    d__1 = x[i__] - l[i__];
				    gi = min(d__1,gi);
				}
		    }
		}
	/* Computing MAX */
		d__1 = *sbgnrm, d__2 = abs(gi);
		*sbgnrm = max(d__1,d__2);
    }
    return 0;
} /* projgr_ */

/* ======================= The end of projgr ============================= */
int subsm_(const integer n, const integer m, const integer nsub, const integer *
	ind, const doublereal *l, const doublereal *u, const integer *nbd, doublereal *x, 
	doublereal *d__, const doublereal *ws, const doublereal *wy, const doublereal theta, 
	const integer col, const integer head, integer *iword, doublereal *wv, 
	const doublereal *wn, const integer iprint, integer *info)
{
    /* Format strings */
    static char fmt_1001[] = "(/,\002----------------SUBSM entered----------\
-------\002,/)";
    static char fmt_1002[] = "(\002ALPHA = \002,f7.5,\002 backtrack to the B\
OX\002)";
    static char fmt_1003[] = "(\002Subspace solution X =  \002,/,(4x,1p,6(1x\
,d11.4)))";
    static char fmt_1004[] = "(/,\002----------------exit SUBSM ------------\
--------\002,/)";

    /* System generated locals */
    integer ws_dim1, ws_offset, wy_dim1, wy_offset, wn_dim1, wn_offset;

    /* Builtin functions */
    integer s_wsfe(cilist *), e_wsfe(), do_fio(integer *, char *, ftnlen), 
	    s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_wsle();

    /* Local variables */
    static integer i__, j, k, m2;
    static doublereal dk;
    static integer js, jy, ibd, col2;
    static doublereal temp1, temp2, alpha;
    static integer pointr;

    /* Fortran I/O blocks */
    static cilist io___233 = { 0, 6, 0, fmt_1001, 0 };
    static cilist io___247 = { 0, 6, 0, fmt_1002, 0 };
    static cilist io___248 = { 0, 6, 0, 0, 0 };
    static cilist io___249 = { 0, 6, 0, fmt_1003, 0 };
    static cilist io___250 = { 0, 6, 0, fmt_1004, 0 };


    /* Parameter adjustments */
    --d__;
    --x;
    --nbd;
    --u;
    --l;
    wn_dim1 = 2 * m;
    wn_offset = 1 + wn_dim1;
    wn -= wn_offset;
    --wv;
    wy_dim1 = n;
    wy_offset = 1 + wy_dim1;
    wy -= wy_offset;
    ws_dim1 = n;
    ws_offset = 1 + ws_dim1;
    ws -= ws_offset;
    --ind;

    /* Function Body */
    if (nsub <= 0) return 0;
    if (iprint >= 99) {
		s_wsfe(&io___233);
		e_wsfe();
    }
/*     Compute wv = W'Zd. */
    pointr = head;
    for (i__ = 1; i__ <= col; ++i__) {
		temp1 = 0.;
		temp2 = 0.;
		for (j = 1; j <= nsub; ++j) {
		    k = ind[j];
		    temp1 += wy[k + pointr * wy_dim1] * d__[j];
		    temp2 += ws[k + pointr * ws_dim1] * d__[j];
		}
		wv[i__] = temp1;
		wv[col + i__] = theta * temp2;
		pointr = pointr % m + 1;
    }
/*     Compute wv:=K^(-1)wv. */
    m2 = m << 1;
    col2 = col << 1;
    dtrsl_(&wn[wn_offset], m2, col2, &wv[1], 11, info);
    if (*info != 0) return 0;
    for (i__ = 1; i__ <= col; ++i__) {
		wv[i__] = -wv[i__];
    }
    dtrsl_(&wn[wn_offset], m2, col2, &wv[1], 1, info);
    if (*info != 0) return 0;

/*     Compute d = (1/theta)d + (1/theta**2)Z'W wv. */
    pointr = head;
    for (jy = 1; jy <= col; ++jy) {
		js = col + jy;
		for (i__ = 1; i__ <= nsub; ++i__) {
		    k = ind[i__];
		    d__[i__] = d__[i__] + wy[k + pointr * wy_dim1] * wv[jy] / theta 
			    + ws[k + pointr * ws_dim1] * wv[js];
		}
		pointr = pointr % m + 1;
    }
    for (i__ = 1; i__ <= nsub; ++i__) d__[i__] /= theta;

/*     Backtrack to the feasible region. */
    alpha = 1.;
    temp1 = alpha;
    for (i__ = 1; i__ <= nsub; ++i__) {
		k = ind[i__];
		dk = d__[i__];
		if (nbd[k] != 0) {
		    if (dk < 0. && nbd[k] <= 2) {
				temp2 = l[k] - x[k];
				if (temp2 >= 0.) temp1 = 0.;
				else if (dk * alpha < temp2) temp1 = temp2 / dk;
		    } 
		    else if (dk > 0. && nbd[k] >= 2) {
				temp2 = u[k] - x[k];
				if (temp2 <= 0.) temp1 = 0.;
				else if (dk * alpha > temp2) temp1 = temp2 / dk;
		    }
		    if (temp1 < alpha) {
				alpha = temp1;
				ibd = i__;
		    }
		}
    }
    if (alpha < 1.) {
		dk = d__[ibd];
		k = ind[ibd];
		if (dk > 0.) {
		    x[k] = u[k];
		    d__[ibd] = 0.;
		} 
		else if (dk < 0.) {
		    x[k] = l[k];
		    d__[ibd] = 0.;
		}
    }
    for (i__ = 1; i__ <= nsub; ++i__) {
		k = ind[i__];
		x[k] += alpha * d__[i__];
    }
    if (iprint >= 99) {
		if (alpha < 1.) {
		    s_wsfe(&io___247);
		    do_fio(&c__1, (char *)&alpha, (ftnlen)sizeof(doublereal));
		    e_wsfe();
		} 
		else {
		    s_wsle(&io___248);
		    do_lio(&c__9, &c__1, "SM solution inside the box", (ftnlen)26);
		    e_wsle();
		}
		if (iprint > 100) {
		    s_wsfe(&io___249);
		    for (i__ = 1; i__ <= n; ++i__) {
				do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
		    }
		    e_wsfe();
		}
    }
    if (alpha < 1.) *iword = 1;
    else *iword = 0;
    if (iprint >= 99) {
		s_wsfe(&io___250);
		e_wsfe();
    }
    return 0;
} /* subsm_ */
/* ====================== The end of subsm =============================== */
int dcsrch_(doublereal *f, doublereal *g, doublereal *stp, 
	const doublereal ftol, const doublereal gtol, const doublereal xtol, 
	const doublereal stpmin, const doublereal stpmax, 
	char *task, integer *isave, doublereal *dsave, ftnlen task_len)
{
    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    integer s_cmp(char *, char *, ftnlen, ftnlen);
    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);

    /* Local variables */
    static doublereal fm, gm, fx, fy, gx, gy, fxm, fym, gxm, gym, stx, sty;
    static integer stage;
    static doublereal finit, ginit, width, ftest, gtest, stmin, stmax, width1;
    static logical brackt;

/*     Initialization block. */
    /* Parameter adjustments */
    --dsave;
    --isave;

    /* Function Body */
    if (s_cmp(task, "START", (ftnlen)5, (ftnlen)5) == 0) {
	/*        Check the input arguments for errors. */
		if (*stp < stpmin) {
		    s_copy(task, "ERROR: STP .LT. STPMIN", task_len, (ftnlen)22);
		}
		if (*stp > stpmax) {
		    s_copy(task, "ERROR: STP .GT. STPMAX", task_len, (ftnlen)22);
		}
		if (*g >= 0.) {
		    s_copy(task, "ERROR: INITIAL G .GE. ZERO", task_len, (ftnlen)26);
		}
		if (ftol < 0.) {
		    s_copy(task, "ERROR: FTOL .LT. ZERO", task_len, (ftnlen)21);
		}
		if (gtol < 0.) {
		    s_copy(task, "ERROR: GTOL .LT. ZERO", task_len, (ftnlen)21);
		}
		if (xtol < 0.) {
		    s_copy(task, "ERROR: XTOL .LT. ZERO", task_len, (ftnlen)21);
		}
		if (stpmin < 0.) {
		    s_copy(task, "ERROR: STPMIN .LT. ZERO", task_len, (ftnlen)23);
		}
		if (stpmax < stpmin) {
		    s_copy(task, "ERROR: STPMAX .LT. STPMIN", task_len, (ftnlen)25);
		}
	/*        Exit if there are errors on input. */
		if (s_cmp(task, "ERROR", (ftnlen)5, (ftnlen)5) == 0) {
		    return 0;
		}
	/*        Initialize local variables. */
		brackt = FALSE_;
		stage = 1;
		finit = *f;
		ginit = *g;
		gtest = ftol * ginit;
		width = stpmax - stpmin;
		width1 = width / .5;
	/*        The variables stx, fx, gx contain the values of the step, 
	          function, and derivative at the best step. 
	          The variables sty, fy, gy contain the value of the step, 
	          function, and derivative at sty.
	          The variables stp, f, g contain the values of the step, 
	          function, and derivative at stp. */
		stx = 0.;
		fx = finit;
		gx = ginit;
		sty = 0.;
		fy = finit;
		gy = ginit;
		stmin = 0.;
		stmax = *stp + *stp * 4.;
		s_copy(task, "FG", task_len, (ftnlen)2);
		goto L1000;
	} 
	else {
	/*        Restore local variables. */
		if (isave[1] == 1) brackt = TRUE_;
		else brackt = FALSE_;
		stage = isave[2];
		ginit = dsave[1];
		gtest = dsave[2];
		gx = dsave[3];
		gy = dsave[4];
		finit = dsave[5];
		fx = dsave[6];
		fy = dsave[7];
		stx = dsave[8];
		sty = dsave[9];
		stmin = dsave[10];
		stmax = dsave[11];
		width = dsave[12];
		width1 = dsave[13];
    }
/*     If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the 
       algorithm enters the second stage. */
    ftest = finit + *stp * gtest;
    if (stage == 1 && *f <= ftest && *g >= 0.) stage = 2;
//     Test for warnings. 
    if (brackt && (*stp <= stmin || *stp >= stmax)) {
		s_copy(task, "WARNING: ROUNDING ERRORS PREVENT PROGRESS", task_len, (
			ftnlen)41);
    }
    if (brackt && stmax - stmin <= xtol * stmax) {
		s_copy(task, "WARNING: XTOL TEST SATISFIED", task_len, (ftnlen)28);
    }
    if (*stp == stpmax && *f <= ftest && *g <= gtest) {
		s_copy(task, "WARNING: STP = STPMAX", task_len, (ftnlen)21);
    }
    if (*stp == stpmin && (*f > ftest || *g >= gtest)) {
		s_copy(task, "WARNING: STP = STPMIN", task_len, (ftnlen)21);
    }
/*     Test for convergence. */
    if (*f <= ftest && abs(*g) <= gtol * (-ginit)) {
		s_copy(task, "CONVERGENCE", task_len, (ftnlen)11);
    }
/*     Test for termination. */
    if (s_cmp(task, "WARN", (ftnlen)4, (ftnlen)4) == 0 || s_cmp(task, "CONV", 
	    (ftnlen)4, (ftnlen)4) == 0) {
		goto L1000;
    }
/*     A modified function is used to predict the step during the 
       first stage if a lower function value has been obtained but 
       the decrease is not sufficient. */
    if (stage == 1 && *f <= fx && *f > ftest) {
	/*        Define the modified function and derivative values. */
		fm = *f - *stp * gtest;
		fxm = fx - stx * gtest;
		fym = fy - sty * gtest;
		gm = *g - gtest;
		gxm = gx - gtest;
		gym = gy - gtest;
	/*        Call dcstep to update stx, sty, and to compute the new step. */
		dcstep_(&stx, &fxm, &gxm, &sty, &fym, &gym, stp, fm, gm, &brackt, 
			stmin, stmax);
	/*        Reset the function and derivative values for f. */
		fx = fxm + stx * gtest;
		fy = fym + sty * gtest;
		gx = gxm + gtest;
		gy = gym + gtest;
    } 
    else {
/*       Call dcstep to update stx, sty, and to compute the new step. */
		dcstep_(&stx, &fx, &gx, &sty, &fy, &gy, stp, *f, *g, &brackt, stmin, 
			stmax);
    }
/*     Decide if a bisection step is needed. */
    if (brackt) {
		if ((d__1 = sty - stx, abs(d__1)) >= width1 * .66) {
		    *stp = stx + (sty - stx) * .5;
		}
		width1 = width;
		width = (d__1 = sty - stx, abs(d__1));
    }
/*     Set the minimum and maximum steps allowed for stp. */
    if (brackt) {
		stmin = min(stx,sty);
		stmax = max(stx,sty);
    } 
    else {
		stmin = *stp + (*stp - stx) * 1.1;
		stmax = *stp + (*stp - stx) * 4.;
    }
/*     Force the step to be within the bounds stpmax and stpmin. */
    *stp = max(*stp,stpmin);
    *stp = min(*stp,stpmax);
/*     If further progress is not possible, let stp be the best */
/*     point obtained during the search. */
    if (brackt && (*stp <= stmin || *stp >= stmax) || brackt && stmax - stmin 
	    <= xtol * stmax) {
		*stp = stx;
    }
/*     Obtain another function and derivative. */
    s_copy(task, "FG", task_len, (ftnlen)2);
L1000:
/*     Save local variables. */
    if (brackt) isave[1] = 1;
    else isave[1] = 0;
    isave[2] = stage;
    dsave[1] = ginit;
    dsave[2] = gtest;
    dsave[3] = gx;
    dsave[4] = gy;
    dsave[5] = finit;
    dsave[6] = fx;
    dsave[7] = fy;
    dsave[8] = stx;
    dsave[9] = sty;
    dsave[10] = stmin;
    dsave[11] = stmax;
    dsave[12] = width;
    dsave[13] = width1;
    return 0;
} /* dcsrch_ */

/* ====================== The end of dcsrch ============================== */
int dcstep_(doublereal *stx, doublereal *fx, doublereal *dx, 
	doublereal *sty, doublereal *fy, doublereal *dy, doublereal *stp, 
	const doublereal fp, const doublereal dp, logical *brackt, const doublereal stpmin, 
	const doublereal stpmax)
{
    /* System generated locals */
    doublereal d__1, d__2, d__3;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static doublereal p, q, r__, s, sgnd, stpc, stpf, stpq, gamma, theta;

    sgnd = dp * (*dx / abs(*dx));
/*     First case: A higher function value. The minimum is bracketed. 
       If the cubic step is closer to stx than the quadratic step, the 
       cubic step is taken, otherwise the average of the cubic and 
       quadratic steps is taken. */
    if (fp > *fx) {
		theta = (*fx - fp) * 3. / (*stp - *stx) + *dx + dp;
	/* Computing MAX */
		d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(dp);
		s = max(d__1,d__2);
	/* Computing 2nd power */
		d__1 = theta / s;
		gamma = s * sqrt(d__1 * d__1 - *dx / s * (dp / s));
		if (*stp < *stx) gamma = -gamma;
	
		p = gamma - *dx + theta;
		q = gamma - *dx + gamma + dp;
		r__ = p / q;
		stpc = *stx + r__ * (*stp - *stx);
		stpq = *stx + *dx / ((*fx - fp) / (*stp - *stx) + *dx) / 2. * (*stp 
			- *stx);
		if ((d__1 = stpc - *stx, abs(d__1)) < (d__2 = stpq - *stx, abs(d__2)))
			 {
		    stpf = stpc;
		} 
		else stpf = stpc + (stpq - stpc) / 2.;
		*brackt = TRUE_;
	/*     Second case: A lower function value and derivatives of opposite 
	       sign. The minimum is bracketed. If the cubic step is farther from 
	       stp than the secant step, the cubic step is taken, otherwise the 
	       secant step is taken. */
    } 
    else if (sgnd < 0.) {
		theta = (*fx - fp) * 3. / (*stp - *stx) + *dx + dp;
	/* Computing MAX */
		d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(dp);
		s = max(d__1,d__2);
	/* Computing 2nd power */
		d__1 = theta / s;
		gamma = s * sqrt(d__1 * d__1 - *dx / s * (dp / s));
		if (*stp > *stx) gamma = -gamma;
		p = gamma - dp + theta;
		q = gamma - dp + gamma + *dx;
		r__ = p / q;
		stpc = *stp + r__ * (*stx - *stp);
		stpq = *stp + dp / (dp - *dx) * (*stx - *stp);
		if ((d__1 = stpc - *stp, abs(d__1)) > (d__2 = stpq - *stp, abs(d__2)))
			 {
		    stpf = stpc;
		} 
		else stpf = stpq;
		*brackt = TRUE_;
	/*     Third case: A lower function value, derivatives of the same sign, 
	     and the magnitude of the derivative decreases. */
    } 
    else if (abs(dp) < abs(*dx)) {
	/*        The cubic step is computed only if the cubic tends to infinity 
	          in the direction of the step or if the minimum of the cubic 
	          is beyond stp. Otherwise the cubic step is defined to be the 
	          secant step. */
		theta = (*fx - fp) * 3. / (*stp - *stx) + *dx + dp;
	/* Computing MAX */
		d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(dp);
		s = max(d__1,d__2);
	/*        The case gamma = 0 only arises if the cubic does not tend 
	          to infinity in the direction of the step. 
	   Computing MAX 
	   Computing 2nd power */
		d__3 = theta / s;
		d__1 = 0., d__2 = d__3 * d__3 - *dx / s * (dp / s);
		gamma = s * sqrt((max(d__1,d__2)));
		if (*stp > *stx) gamma = -gamma;
		p = gamma - dp + theta;
		q = gamma + (*dx - dp) + gamma;
		r__ = p / q;
		if (r__ < 0. && gamma != 0.) stpc = *stp + r__ * (*stx - *stp);
		else if (*stp > *stx) stpc = stpmax;
		else stpc = stpmin;
		stpq = *stp + dp / (dp - *dx) * (*stx - *stp);
		if (*brackt) {
	/*           A minimizer has been bracketed. If the cubic step is 
	             closer to stp than the secant step, the cubic step is 
	             taken, otherwise the secant step is taken. */
		    if ((d__1 = stpc - *stp, abs(d__1)) < (d__2 = stpq - *stp, abs(
			    d__2))) {
				stpf = stpc;
		    } 
		    else stpf = stpq;
		    if (*stp > *stx) {
		// Computing MIN 
				d__1 = *stp + (*sty - *stp) * .66;
				stpf = min(d__1,stpf);
		    } 
		    else {
		// Computing MAX 
				d__1 = *stp + (*sty - *stp) * .66;
				stpf = max(d__1,stpf);
		    }
		} 
		else {
	/*           A minimizer has not been bracketed. If the cubic step is 
	             farther from stp than the secant step, the cubic step is 
	             taken, otherwise the secant step is taken. */
		    if ((d__1 = stpc - *stp, abs(d__1)) > (d__2 = stpq - *stp, abs(
			    d__2))) {
				stpf = stpc;
		    } 
		    else stpf = stpq;
		    stpf = min(stpmax,stpf);
		    stpf = max(stpmin,stpf);
		}
	/*     Fourth case: A lower function value, derivatives of the same sign, 
	       and the magnitude of the derivative does not decrease. If the 
	       minimum is not bracketed, the step is either stpmin or stpmax, 
	       otherwise the cubic step is taken. */
    } 
    else {
		if (*brackt) {
		    theta = (fp - *fy) * 3. / (*sty - *stp) + *dy + dp;
	/* Computing MAX */
		    d__1 = abs(theta), d__2 = abs(*dy), d__1 = max(d__1,d__2), d__2 = 
			    abs(dp);
		    s = max(d__1,d__2);
	/* Computing 2nd power */
		    d__1 = theta / s;
		    gamma = s * sqrt(d__1 * d__1 - *dy / s * (dp / s));
		    if (*stp > *sty) gamma = -gamma;
		    p = gamma - dp + theta;
		    q = gamma - dp + gamma + *dy;
		    r__ = p / q;
		    stpc = *stp + r__ * (*sty - *stp);
		    stpf = stpc;
		} 
		else if (*stp > *stx) stpf = stpmax;
		else stpf = stpmin;
    }
//     Update the interval which contains a minimizer. 
    if (fp > *fx) {
		*sty = *stp;
		*fy = fp;
		*dy = dp;
    } 
    else {
		if (sgnd < 0.) {
		    *sty = *stx;
		    *fy = *fx;
		    *dy = *dx;
		}
		*stx = *stp;
		*fx = fp;
		*dx = dp;
    }
//     Compute the new step. 
    *stp = stpf;
    return 0;
} /* dcstep_ */

/* ====================== The end of dcstep ============================== */
int timer_(doublereal *ttime)
{
/*    static real temp;
    static real tarray[2];*/

        
/*    temp = etime_(tarray);
    *ttime = (doublereal) tarray[0];*/
    *ttime = (clock() * 1000)/(CLOCKS_PER_SEC);
    return 0;
} /* timer_ */

/* ====================== The end of timer =============================== */
doublereal dnrm2_(const integer n, const doublereal *x, const integer incx)
{
    /* System generated locals */
    integer i__1, i__2;
    doublereal ret_val, d__1, d__2, d__3;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__;
    static doublereal scale;

    // Parameter adjustments 
    --x;

    /* Function Body */
    ret_val = 0.;
    scale = 0.;
    i__1 = n;
    i__2 = incx;
    for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
	/* Computing MAX */
		d__2 = scale, d__3 = (d__1 = x[i__], abs(d__1));
		scale = max(d__2,d__3);
    }
    if (scale == 0.) return ret_val;
    i__2 = n;
    i__1 = incx;
    for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
	/* Computing 2nd power */
		d__1 = x[i__] / scale;
		ret_val += d__1 * d__1;
    }
    ret_val = scale * sqrt(ret_val);
    return ret_val;
} /* dnrm2_ */

/* ====================== The end of dnrm2 =============================== */
doublereal dpmeps_()
{
    /* Initialized data */

    const static doublereal zero = 0.;
    const static doublereal one = 1.;
    const static doublereal two = 2.;

    /* System generated locals */
    integer i__1;
    doublereal ret_val;

    /* Local variables */
    static doublereal a, b;
    static integer i__, it;
    static doublereal beta;
    static integer irnd;
    static doublereal temp, temp1, betah;
    static integer ibeta, negep;
    static doublereal tempa;
    static integer itemp;
    static doublereal betain;

    a = one;
    b = one;
L10:
    a += a;
    temp = a + one;
    temp1 = temp - a;
    if (temp1 - one == zero) goto L10;
L20:
    b += b;
    temp = a + b;
    itemp = (integer) (temp - a);
    if (itemp == 0) goto L20;
    ibeta = itemp;
    beta = (doublereal) ibeta;
/*     determine it, irnd. */
    it = 0;
    b = one;
L30:
    ++it;
    b *= beta;
    temp = b + one;
    temp1 = temp - b;
    if (temp1 - one == zero) goto L30;
    irnd = 0;
    betah = beta / two;
    temp = a + betah;
    if (temp - a != zero) irnd = 1;
    tempa = a + beta;
    temp = tempa + betah;
    if (irnd == 0 && temp - tempa != zero) irnd = 2;
/*     determine dpmeps. */
    negep = it + 3;
    betain = one / beta;
    a = one;
    i__1 = negep;
    for (i__ = 1; i__ <= i__1; ++i__) a *= betain;
L50:
    temp = one + a;
    if (temp - one != zero) goto L60;
    a *= beta;
    goto L50;
L60:
    ret_val = a;
    if (ibeta == 2 || irnd == 0) goto L70;
    a = a * (one + a) / two;
    temp = one + a;
    if (temp - one != zero) ret_val = a;
L70:
    return ret_val;
} /* dpmeps_ */

/* ====================== The end of dpmeps ============================== */
int daxpy_(const integer n, const doublereal da, const doublereal *dx, 
	const integer incx, doublereal *dy, const integer incy)
{

    /* Local variables */
    static integer i__, m, ix, iy, mp1;

    /* Parameter adjustments */
    --dy;
    --dx;

    /* Function Body */
    if (n <= 0) return 0;
    if (da == 0.) return 0;
    if (incx == 1 && incy == 1) goto L20;

/*        code for unequal increments or equal increments 
          not equal to 1 */

    ix = 1;
    iy = 1;
    if (incx < 0) ix = (-n + 1) * incx + 1;
    if (incy < 0) iy = (-n + 1) * incy + 1;
    for (i__ = 1; i__ <= n; ++i__) {
		dy[iy] += da * dx[ix];
		ix += incx;
		iy += incy;
    }
    return 0;

/*        code for both increments equal to 1 */

/*        clean-up loop */

L20:
    m = n % 4;
    if (m == 0) goto L40;
    for (i__ = 1; i__ <= m; ++i__) {
		dy[i__] += da * dx[i__];
    }
    if (n < 4) return 0;

L40:
    mp1 = m + 1;
    for (i__ = mp1; i__ <= n; i__ += 4) {
		dy[i__] += da * dx[i__];
		dy[i__ + 1] += da * dx[i__ + 1];
		dy[i__ + 2] += da * dx[i__ + 2];
		dy[i__ + 3] += da * dx[i__ + 3];
    }
    return 0;
} /* daxpy_ */

/* ====================== The end of daxpy =============================== */
int dcopy_(const integer n, const doublereal *dx, const integer incx, 
	doublereal *dy, const integer incy)
{
    /* Local variables */
    static integer i__, m, ix, iy, mp1;

    /* Parameter adjustments */
    --dy;
    --dx;

    /* Function Body */
    if (n <= 0) return 0;
    if (incx == 1 && incy == 1) goto L20;

/*        code for unequal increments or equal increments */
/*          not equal to 1 */

    ix = 1;
    iy = 1;
    if (incx < 0) ix = (-n + 1) * incx + 1;
    if (incy < 0) iy = (-n + 1) * incy + 1;
    for (i__ = 1; i__ <= n; ++i__) {
		dy[iy] = dx[ix];
		ix += incx;
		iy += incy;
    }
    return 0;

/*        code for both increments equal to 1 */

/*        clean-up loop */

L20:
    m = n % 7;
    if (m == 0) goto L40;
    for (i__ = 1; i__ <= m; ++i__) {
		dy[i__] = dx[i__];
    }
    if (n < 7) return 0;
L40:
    mp1 = m + 1;
    for (i__ = mp1; i__ <= n; i__ += 7) {
		dy[i__] = dx[i__];
		dy[i__ + 1] = dx[i__ + 1];
		dy[i__ + 2] = dx[i__ + 2];
		dy[i__ + 3] = dx[i__ + 3];
		dy[i__ + 4] = dx[i__ + 4];
		dy[i__ + 5] = dx[i__ + 5];
		dy[i__ + 6] = dx[i__ + 6];
    }
    return 0;
} /* dcopy_ */

/* ====================== The end of dcopy =============================== */
doublereal ddot_(const integer n, const doublereal *dx, const integer incx, 
	const doublereal *dy, const integer incy)
{
    /* System generated locals */
    doublereal ret_val;

    /* Local variables */
    static integer i__, m, ix, iy, mp1;
    static doublereal dtemp;

    /* Parameter adjustments */
    --dy;
    --dx;

    /* Function Body */
    ret_val = 0.;
    dtemp = 0.;
    if (n <= 0) return ret_val;
    if (incx == 1 && incy == 1) goto L20;

/*        code for unequal increments or equal increments */
/*          not equal to 1 */

    ix = 1;
    iy = 1;
    if (incx < 0) ix = (-n + 1) * incx + 1;
    if (incy < 0) iy = (-n + 1) * incy + 1;
    for (i__ = 1; i__ <= n; ++i__) {
		dtemp += dx[ix] * dy[iy];
		ix += incx;
		iy += incy;
    }
    ret_val = dtemp;
    return ret_val;

/*        code for both increments equal to 1 */
/*        clean-up loop */

L20:
    m = n % 5;
    if (m == 0) goto L40;
    for (i__ = 1; i__ <= m; ++i__) {
		dtemp += dx[i__] * dy[i__];
    }
    if (n < 5) goto L60;
L40:
    mp1 = m + 1;
    for (i__ = mp1; i__ <= n; i__ += 5) {
		dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
			i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ + 
			4] * dy[i__ + 4];
    }
L60:
    ret_val = dtemp;
    return ret_val;
} /* ddot_ */

/* ====================== The end of ddot ================================ */
int dpofa_(doublereal *a, const integer lda, const integer n, 
	integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer j, k;
    static doublereal s, t;
    static integer jm1;
    
    /* Parameter adjustments */
    a_dim1 = lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;

    /* Function Body */
    for (j = 1; j <= n; ++j) {
		*info = j;
		s = 0.;
		jm1 = j - 1;
		if (jm1 < 1) goto L20;
		for (k = 1; k <= jm1; ++k) {
		    t = a[k + j * a_dim1] - ddot_(k - 1, &a[k * a_dim1 + 1], 1, &
			    a[j * a_dim1 + 1], 1);
		    t /= a[k + k * a_dim1];
		    a[k + j * a_dim1] = t;
		    s += t * t;
		}
	L20:
		s = a[j + j * a_dim1] - s;
	/*     ......exit */
		if (s <= 0.) goto L40;
		a[j + j * a_dim1] = sqrt(s);
    }
    *info = 0;
L40:
    return 0;
} /* dpofa_ */

/* ====================== The end of dpofa =============================== */
int dscal_(const integer n, const doublereal da, doublereal *dx, 
	const integer incx)
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer i__, m, mp1, nincx;

    /* Parameter adjustments */
    --dx;

    /* Function Body */
    if (n <= 0 || incx <= 0) return 0;
    if (incx == 1) goto L20;

/*        code for increment not equal to 1 */

    nincx = n * incx;
    i__1 = nincx;
    i__2 = incx;
    for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
		dx[i__] = da * dx[i__];
    }
    return 0;

/*        code for increment equal to 1 */
/*        clean-up loop */

L20:
    m = n % 5;
    if (m == 0) goto L40;
    for (i__ = 1; i__ <= m; ++i__) {
		dx[i__] = da * dx[i__];
    }
    if (n < 5) return 0;
L40:
    mp1 = m + 1;
    for (i__ = mp1; i__ <= n; i__ += 5) {
		dx[i__] = da * dx[i__];
		dx[i__ + 1] = da * dx[i__ + 1];
		dx[i__ + 2] = da * dx[i__ + 2];
		dx[i__ + 3] = da * dx[i__ + 3];
		dx[i__ + 4] = da * dx[i__ + 4];
    }
    return 0;
} /* dscal_ */

/* ====================== The end of dscal =============================== */
int dtrsl_(const doublereal *t, const integer ldt, const integer n, 
	doublereal *b, const integer job, integer *info)
{
    /* System generated locals */
    integer t_dim1, t_offset;

    /* Local variables */
    static integer j, jj, case__;
    static doublereal temp;


    /* Parameter adjustments */
    t_dim1 = ldt;
    t_offset = 1 + t_dim1;
    t -= t_offset;
    --b;

    /* Function Body */
    for (*info = 1; *info <= n; ++(*info)) {
/*     ......exit */
		if (t[*info + *info * t_dim1] == 0.) goto L150;
    }
    *info = 0;

/*        determine the task and go to it. */

    case__ = 1;
    if (job % 10 != 0) case__ = 2;
    if (job % 100 / 10 != 0) case__ += 2;
    switch (case__) {
		case 1:  goto L20;
		case 2:  goto L50;
		case 3:  goto L80;
		case 4:  goto L110;
    }

/*        solve t*x=b for t lower triangular */

L20:
    b[1] /= t[t_dim1 + 1];
    if (n < 2) goto L40;
    for (j = 2; j <= n; ++j) {
		temp = -b[j - 1];
		daxpy_(n - j + 1, temp, &t[j + (j - 1) * t_dim1], 1, &b[j], 1);
		b[j] /= t[j + j * t_dim1];
    }
L40:
    goto L140;

/*        solve t*x=b for t upper triangular. */

L50:
    b[n] /= t[n + n * t_dim1];
    if (n < 2) goto L70;
    for (jj = 2; jj <= n; ++jj) {
		j = n - jj + 1;
		temp = -b[j + 1];
		daxpy_(j, temp, &t[(j + 1) * t_dim1 + 1], 1, &b[1], 1);
		b[j] /= t[j + j * t_dim1];
    }
L70:
    goto L140;

/*        solve trans(t)*x=b for t lower triangular. */

L80:
    b[n] /= t[n + n * t_dim1];
    if (n < 2) goto L100;
    for (jj = 2; jj <= n; ++jj) {
		j = n - jj + 1;
		b[j] -= ddot_(jj - 1, &t[j + 1 + j * t_dim1], 1, &b[j + 1], 1);
		b[j] /= t[j + j * t_dim1];
    }
L100:
    goto L140;

/*        solve trans(t)*x=b for t upper triangular. */

L110:
    b[1] /= t[t_dim1 + 1];
    if (n < 2) goto L130;
    for (j = 2; j <= n; ++j) {
		b[j] -= ddot_(j - 1, &t[j * t_dim1 + 1], 1, &b[1], 1);
		b[j] /= t[j + j * t_dim1];
    }
L130:
L140:
L150:
    return 0;
} /* dtrsl_ */

#ifdef __cplusplus
	}
#endif