Attachment 14459 Details for Bug 73524 – After opening this file, eclipse hangs

After opening this file, eclipse hangs

routines.c (text/plain), 106.10 KB, created by Boris von Loesch

on 2004-09-09 09:46:48 EDT

(hide)

Description:

Filename:

MIME Type:

Creator: Boris von Loesch

Created: 2004-09-09 09:46:48 EDT

Size: 106.10 KB

patch

obsolete

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

Actions: View

Attachments on bug 73524: 14459