/****************************************************************
Copyright (C) 1997, 2001 Lucent Technologies
All Rights Reserved

Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.

LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/

#include "nlp.h"
#undef f_OPNUM
#include "opcode.hd"
#include "nlc.h"
#define asl ((ASL_fg*)cur_ASL)

#define Egulp 400

typedef struct maxinfo {
	int ncond;
	int ndv;
	int nvt;
	int needT1;
	/* !! plus goo for function evaluations */
	} maxinfo;

static maxinfo bmax, cmax, omax;

 static int cwant = 2, derkind = 2, djoff = -1, owant = 2;
 static int needx0check, nvseen, output_time, *vseen;
 static list **c_ifset, **com_ifset, **o_ifset;
 int branches;
 static int Fortran1, condlevel, fwalk, iflevel, maxa, ncond, ndv, ndvmax,
	npd, nplterm, nv1, nv2, nvt, nvtmax;
 static int *cvmap;
 static fint nJ;
 static expr *Plterms;
 static char *intsk;
 efunc *r_ops[N_OPS];
 extern efunc *r_op[];
#include "r_ops1.hd"

 static real *bounds;
 extern char *progname;
 char *Half = "0.5", *One = "1.", *Negone = "-1.", *Zero = "0.";
 char *opEQ = "==", *opGE = ">=", *opGT = ">", *opLE = "<=",
	*opLT = "<", *opNE = "!=";
 char *opAND = "&&", *opOR = "||";
 char *T = "t1", *T1 = "t2";	/* temporary variable names */
 char *offlfmt1 = "&%s", *offlfmt2 = "&%s, &%s";
 static char *assign_fmt = "\t%s = %s;\n";
 static char *binop_fmt = "\t%s = %s %s %s;\n";
 static char *goto_fmt = "\tgoto L%d;\n";
 static char *ifgo_fmt = "\tif (%s %s %s) goto L%d;\n";
 static char *label_fmt = " L%d:\n";
 static char *ifstart_fmt = "\tif (%s %s %s) {\n";
 static char *elseif_fmt = "\t} else if (%s %s %s) {\n";
 static char *else_fmt = "\t} else {\n";
 static char *endif_fmt = "\t}\n";
 static char *case_fmt = "   case %d:\n";
 char *cond_fmt = "cond[%d]";
 static char *break_fmt = "\tbreak;\n";
 static char *endswitch_fmt = "\t}\n";
 static char *zerdiv_fmt = "\tzerdiv_(&%s);";
 static char *call_fmt = "\t%s(&%s);\n";
 static char *dv_fmt = "dv[%d]";
 static char *Fortstar = "";
 char *pd_fmt = "pd[%d]";
 static char *tv_fmt = "v[%d]";
 static char *x_fmt = "x[%d]";
 static char seen[N_OPS], *declare[N_OPS];
 static char *grad_fmt = "g[%d]", *jac_fmt = "J[%d]";
 static char *g0_fmt = "\tg[%d] = %s;\n", *j0_fmt = "\tJ[%d] = %s;\n";
 static char *eos = ";\n";
 static char *star = "";
 static char *Void = "void";
 static char *xcheck = "\tfint wantfg = *needfg;\n\
	if (xcheck(x) && wantfg == 2)\n\t\twantfg = 3;\n";
 static char *xcheck0 = "\txcheck(x);\n";
 static char *xcheckdcl = "(real *x)";
 static char *xkind = "\tif (!(xkind & %d)) {\n\t\txkind |= %d;\n";
 static real negone = -1.0, one = 1.0;
 extern real edagread_one;	/* &edagread_one = special derp->c value */

 static expr_nx *nums;
 static ograd **stog;

 static Adjoint *dwalk_one;
 static real *rdwalk_one;
 static Adjoint A0;
 static int *acount, *dvtfree, *rownos;
 static int densejac, krflag, ndvfree, ndvtmax, needT1, nvtfree,
	    pdts, stor_grad;

 typedef struct
pdl {
	struct pdl *next;
	int i, ts;
	} pdl;

 static pdl *pdlbusy, *pdlfree;

 static int
new_pd(VOID)
{
	pdl *p;
	if (!iflevel)
		return npd++;
	if (p = pdlfree)
		pdlfree = p->next;
	else {
		p = (pdl *)mem(sizeof(pdl));
		p->i = npd++;
		}
	p->next = pdlbusy;
	pdlbusy = p;
	p->ts = pdts;
	return p->i;
	}

 static int
#ifdef KR_headers
pdlsave(opdb)
	pdl **opdb;
#else
pdlsave(pdl **opdb)
#endif
{
	condlevel++;
	if (iflevel++)
		*opdb = pdlbusy;
	else {
		pdlfree = 0;
		*opdb = 0;
		}
	pdlbusy = 0;
	return ++pdts;
	}

 static void
pdlreset(VOID)
{
	pdl *p, *pnext;
	for(p = pdlbusy; p; p = pnext) {
		pnext = p->next;
		p->next = pdlfree;
		pdlfree = p;
		}
	pdlbusy = 0;
	}

 static void
#ifdef KR_headers
pdlrestore(opdb, mts)
	pdl *opdb;
	int mts;
#else
pdlrestore(pdl *opdb, int mts)
#endif
{
	pdl *p, *pnext;
	condlevel--;
	if (--iflevel) {
		pnext = pdlbusy;
		while(p = pnext) {
			pnext = p->next;
			p->next = opdb;
			opdb = p;
			}
		pnext = pdlfree;
		while((p = pnext) && p->ts >= mts) {
			pnext = p->next;
			p->next = opdb;
			opdb = p;
			}
		pdlfree = pnext;
		pdlbusy = opdb;
		}
	}

#define NDVTGULP 1000

 static void
dvtfree_inc(VOID)
{
	int i, j;
	i = ndvtmax;
	j = ndvtmax += NDVTGULP;
	dvtfree = (int *)Realloc(dvtfree, j*sizeof(int));
	while(i > nvtfree)
		dvtfree[--j] = dvtfree[--i];
	nvtfree += NDVTGULP;
	}

 static int
new_dv(VOID)
{ return ndvfree ? dvtfree[--ndvfree] : ndv++; }

 static void
#ifdef KR_headers
dv_free(k)
	int k;
#else
dv_free(int k)
#endif
{
	if (ndvfree >= nvtfree)
		dvtfree_inc();
	dvtfree[ndvfree++] = k;
	}

 static int
#ifdef KR_headers
new_vt(deriv)
	int deriv;
#else
new_vt(int deriv)
#endif
{
	if (deriv & 2)
		return -(new_pd() + 1);
	return nvtfree < ndvtmax ? dvtfree[nvtfree++] : nvt++;
	}

 static void
#ifdef KR_headers
vt_free(k)
	int k;
#else
vt_free(int k)
#endif
{
	if (ndvfree >= nvtfree)
		dvtfree_inc();
	dvtfree[--nvtfree] = k;
	}

 static list *list_freelist;

 static list *
#ifdef KR_headers
new_list(nxt)
	list *nxt;
#else
new_list(list *nxt)
#endif
{
	list *rv;
	if (rv = list_freelist)
		list_freelist = rv->next;
	else
		rv = (list *)mem(sizeof(list));
	rv->next = nxt;
	return rv;
	}

 static v_i *vi_freelist;

 static v_i *
#ifdef KR_headers
new_vi(v, i)
	real v;
	int i;
#else
new_vi(real v, int i)
#endif
{
	v_i *rv;
	if (rv = vi_freelist)
		vi_freelist = rv->u.next;
	else
		rv = (v_i *)mem(sizeof(v_i));
	rv->u.v = v;
	rv->i = i;
	return rv;
	}

 static void
#ifdef KR_headers
free_vi(vi)
	v_i *vi;
#else
free_vi(v_i *vi)
#endif
{
	vi->u.next = vi_freelist;
	vi_freelist = vi;
	}

 static char *
#ifdef KR_headers
fpval(r)
	real r;
#else
fpval(real r)
#endif
{
	static char buf[32];
	if (r >= Infinity)
		return "1.7e308";
	else if (r <= negInfinity)
		return "-1.7e308";
	g_fmt(buf, r);
	return buf;
	}

 void
#ifdef KR_headers
assign(a, b)
	char *a;
	char *b;
#else
assign(char *a, char *b)
#endif
{
	if (a)
		printf(assign_fmt, a, b);
	}

 void
#ifdef KR_headers
binop(a, b, op, c)
	char *a;
	char *b;
	char *op;
	char *c;
#else
binop(char *a, char *b, char *op, char *c)
#endif
{
	if (a) {
		if (a == b && !Fortran)
			printf("\t%s %s= %s;\n", a, op, c);
		else
			printf(binop_fmt, a, b, op, c);
		}
	}

 void
#ifdef KR_headers
Goto(n)
	int n;
#else
Goto(int n)
#endif
{ printf(goto_fmt, n); }

 void
#ifdef KR_headers
ifgo(a, op, b, n)
	char *a;
	char *op;
	char *b;
	int n;
#else
ifgo(char *a, char *op, char *b, int n)
#endif
{ printf(ifgo_fmt, a, op, b, n); }

 void
#ifdef KR_headers
label(n)
	int n;
#else
label(int n)
#endif
{ printf(label_fmt, n); }

 void
#ifdef KR_headers
ifstart(a, op, b)
	char *a;
	char *op;
	char *b;
#else
ifstart(char *a, char *op, char *b)
#endif
{ printf(ifstart_fmt, a, op, b); }

 void
elsestart(VOID)
{ printf(else_fmt); }

 void
#ifdef KR_headers
elseif(a, b, c)
	char *a;
	char *b;
	char *c;
#else
elseif(char *a, char *b, char *c)
#endif
{ printf(elseif_fmt, a, b, c); }

 void
endif(VOID)
{ printf(endif_fmt); }

 int
#ifdef KR_headers
Switch(v, n)
	char *v;
	int n;
#else
Switch(char *v, int n)
#endif
{
	int i, j, rv;
	if (Fortran) {
		rv = ++branches;
		branches += n;
		printf("\tgoto(");
		for(i = j = 0, --n; i < n; i++) {
			printf("%d,", rv+i);
			if (++j == 9) {
				printf("\n     x\t\t");
				j = 0;
				}
			}
		printf("%d),%s\n", rv+n, v);
		return rv;
		}
	printf("\tswitch(%s) {\n", v);
	return 0;
	}

 void
#ifdef KR_headers
Case(k)
	int k;
#else
Case(int k)
#endif
{ printf(case_fmt, k); }

 void
#ifdef KR_headers
Break(k)
	int k;
#else
Break(int k)
#endif
{ printf(break_fmt, k); }

 void
#ifdef KR_headers
endswitch(lbl)
	int lbl;
#else
endswitch(int lbl)
#endif
{ printf(endswitch_fmt, lbl); }

 void
#ifdef KR_headers
domain(s, x)
	char *s;
	char *x;
#else
domain(char *s, char *x)
#endif
{
	if (Fortran)
		printf("\tcall domain(\"%s\", %s)\n", s, x);
	else
		printf("\tdomain_(\"%s\", %s, %dL);\n", s, x, strlen(s));
	}

 void
#ifdef KR_headers
zerdiv(s)
	char *s;
#else
zerdiv(char *s)
#endif
{ printf(zerdiv_fmt, s); }

 void
#ifdef KR_headers
call(what, arg)
	char *what;
	char *arg;
#else
call(char *what, char *arg)
#endif
{ printf(call_fmt, what, arg); }

 char *
#ifdef KR_headers
call1(what, a)
	char *what;
	char *a;
#else
call1(char *what, char *a)
#endif
{
	static char buf[48];
	sprintf(buf, "%s(%s)", what, a);
	return buf;
	}

 char *
#ifdef KR_headers
call2(what, a, b)
	char *what;
	char *a;
	char *b;
#else
call2(char *what, char *a, char *b)
#endif
{
	static char buf[64];
	sprintf(buf, "%s(%s, %s)", what, a, b);
	return buf;
	}

 void
#ifdef KR_headers
introuble(who, x)
	char *who;
	char *x;
#else
introuble(char *who, char *x)
#endif
{
	if (!Fortran)
		printf("\tif (errno) in_trouble(\"%s\",%s);\n",
			who, x);
	}

 void
#ifdef KR_headers
introuble2(who, x, y)
	char *who;
	char *x;
	char *y;
#else
introuble2(char *who, char *x, char *y)
#endif
{
	if (!Fortran)
		printf("\tif (errno) in_trouble2(\"%s\",%s,%s);\n",
			who, x, y);
	}

 char *
#ifdef KR_headers
num(x)
	int x;
#else
num(int x)
#endif
{
	static char buf[16];
	sprintf(buf, "%d", x);
	return buf;
	}

 static char op_type[] = {
#include "op_type.hd"
	};

 static void
#ifdef KR_headers
dvset(dp, deriv)
	register real *dp;
	int deriv;
#else
dvset(register real *dp, int deriv)
#endif
{
	register dLR *d;
	real t;

	t = *dp;
	d = Make_dLR(dp);
	if (!deriv) {
		d->kind = dLR_UNUSED;
		return;
		}
	if (t == -1.) {
		d->kind = dLR_negone;
		d->o.vp = &negone;
		}
	else if (t == 1.) {
		d->kind = dLR_one;
		d->o.vp = &one;
		}
	else {
		d->kind = dLR_PD;
		d->o.i = new_pd();
		}
	}

 static void
#ifdef KR_headers
Lset(L, nlin)
	linpart *L;
	int nlin;
#else
Lset(linpart *L, int nlin)
#endif
{
	linpart *Le;
	double t;
	register dLR *d;

	for(Le = L + nlin; L < Le; L++) {
		t = L->fac;
		d = Make_dLR(&L->fac);
		if (t == 1.)
			d->kind = dLR_one;
		else if (t == -1.)
			d->kind = dLR_negone;
		else {
			d->kind = dLR_VP;
			*(d->o.vp = (real *)mem(sizeof(real))) = t;
			}
		}
	}

 static void
#ifdef KR_headers
mpd(e, dp, deriv)
	expr *e;
	real *dp;
	int deriv;
#else
mpd(expr *e, real *dp, int deriv)
#endif
{
	int i, j;
	register dLR *d = Make_dLR(dp);

	if (!deriv) {
		d->kind = dLR_UNUSED;
		if (e->op != f_OPNUM && e->op != f_OPVARVAL && (i = e->a) > 0)
			vt_free(i);
		}
	else if (e->op == f_OPNUM) {
		d->kind = dLR_VP;
		d->o.vp = &((expr_nx *)e)->v;
		i = 0;
		}
	else {
		i = e->a;
		if (e->op == f_OPVARVAL)
			if (i >= nv1) {
				i = (expr_v *)e - var_e;
				if (!(j = cvmap[i -= nv1]))
					j = cvmap[i] = -(++npd);
				d->o.i = -(1 + j);
				d->kind = dLR_PD;
				}
			else {
				d->kind = dLR_VARVAL;
				d->o.i = i;
				}
		else if (i < 0) {
			d->kind = dLR_PD;
			d->o.i = -(1 + i);
			}
		else if (deriv & 2) {
			d->kind = dLR_PD;
			e->a = -((d->o.i = new_pd()) + 1);
			}
		else /* DEBUG */ {
			fprintf(Stderr, "Bug in mpd!\n");
			exit(13);
			}
		}
	}

 static int ewalk(expr*, int);

 static void
#ifdef KR_headers
lwalk(ep, neg)
	expr **ep;
	int neg;
#else
lwalk(expr **ep, int neg)
#endif
{
	int i, j;
	expr *e = *ep;
 top:
	switch(Intcast e->op) {
		case OPNOT:
			neg = 1 - neg;
			e = *ep = e->L.e;
			goto top;
		case OPOR:
			e->op = neg ? f_OPAND : f_OPOR;
			goto andor;
		case OPAND:
			e->op = neg ? f_OPOR : f_OPAND;
 andor:
			lwalk(&e->L.e, neg);
			ep = &e->R.e;
			e = *ep;
			goto top;
		case LT:
			e->op = neg ? f_GE : f_LT;
			goto compare;
		case LE:
			e->op = neg ? f_GT : f_LE;
			goto compare;
		case EQ:
			e->op = neg ? f_NE : f_EQ;
			goto compare;
		case GE:
			e->op = neg ? f_LT : f_GE;
			goto compare;
		case GT:
			e->op = neg ? f_LE : f_GT;
			goto compare;
		case NE:
			e->op = neg ? f_EQ : f_NE;
 compare:
			i = ewalk(e->L.e,0);
			j = ewalk(e->R.e,0);
			if (i > 0)
				vt_free(i);
			if (j > 0)
				vt_free(j);
		}
	}

 static char *
#ifdef KR_headers
f_OPVARVAL1(e, buf)
	expr *e;
	char *buf;
#else
f_OPVARVAL1(expr *e, char *buf)
#endif
{
	int k;
	Adjoint *A;
	char *fmt;

	if ((k = e->a) < nv1) {
		if (k < 0) {
			fmt = pd_fmt;
			k = Fortran1 - k;
			}
		else {
			A = Adjp(&adjoints[k]);
			if (!A->seen) {
				A->seen = 1;
				vseen[nvseen++] = e->a;
				}
			fmt = x_fmt;
			k += Fortran;
			}
		}
	else {
		k = cvmap[(expr_v *)e - var_e - nv1];
		if (k < 0) {
			fmt = pd_fmt;
			k = Fortran1 - k;
			}
		else {
			fmt = tv_fmt;
			k += Fortran1;
			}
		}
	sprintf(buf, fmt, k);
	return buf;
	}

 char *
#ifdef KR_headers
f_OPNUM1(e, rv)
	register expr *e;
	char *rv;
#else
f_OPNUM1(register expr *e, char *rv)
#endif
{
	g_fmt(rv,((expr_nx *)e)->v);
	return rv;
	}

 static int
#ifdef KR_headers
ewalk(e, deriv)
	expr *e;
	int deriv;
#else
ewalk(expr *e, int deriv)
#endif
{
	int deriv1, i, j, k, k1, lderiv, mts, rderiv;
	expr *e1, **ep, **epe;
	expr_if *eif;
	expr_va *eva;
	expr_f *ef;
	argpair *ap, *ape;
	de *d;
	derp *dp;
	dLR *LR, **LRp;
	efunc *op;
	Adjoint *A;
	double t;
	pdl *opdb;
	static int achk[12] = { 0, 1, 2, 3, 4, 5, 6, 7, 0, 0, 10, 11 };

	k = Intcast e->op;
	e->op = r_op[k];
	deriv1 = deriv & 1;
	if (achk[k1 = op_type[k]] && e->a == nv1)
		deriv = deriv1 = 0;
	switch(k1) {
		case 11: /* OP2POW */
			e->dL = 0;
			/* no break */
		case 1: /* unary */
			j = ewalk(e->L.e, deriv1);
			i = new_vt(deriv);
			if (deriv)
				dvset(k == OPCPOW ? &e->dR : &e->dL, deriv);
			if (j > 0)
				vt_free(j);
			seen[k] = 1;
 reti:
			return e->a = i;

		case 2: /* binary */
			k1 = 1;
			switch(k) {
				case OPPLUS:
				case OPMINUS:
				case OPREM:
					e->dL = 1.;
					break;
				case OPMULT:
					k1 = 3;
					/* no break */
				default:
					e->dL = 0.;
				}
			if (lderiv = rderiv = deriv1) {
				lderiv = achk[op_type[Intcast e->L.e->op]]
					&& e->L.e->a != nv1 ? k1 : 0;
				rderiv = achk[op_type[Intcast e->R.e->op]]
					&& e->R.e->a != nv1 ? k1 : 0;
				}
			i = ewalk(e->L.e,lderiv);
			j = ewalk(e->R.e,rderiv);
			if (deriv) {
				switch(k) {
					case OPREM:
					case OPintDIV:
					case OPprecision:
					case OPround:
					case OPtrunc:
						rderiv = 0;
						break;
					case OPMULT:
						mpd(e->R.e, &e->dL, lderiv);
						mpd(e->L.e, &e->dR, rderiv);
						goto retnew;
					case OP_atan2:
						needT1 = 1;
					}
				dvset(&e->dL, lderiv);
				dvset(&e->dR, rderiv);
				k = new_vt(deriv1);
				if (i > 0)
					vt_free(i);
				if (j > 0)
					vt_free(j);
				return e->a = k;
				}
			k = i;
			i = new_vt(0);
			if (j > 0)
				vt_free(j);
			if (k > 0)
				vt_free(k);
			goto reti;

		case 3: /* vararg (min, max) */
			eva = (expr_va *)e;
			d = eva->L.d;
			condlevel++;
			if (!(i = ewalk(d->e,deriv)))
				i = new_vt(deriv);
			while(e1 = (++d)->e) {
				pdlreset();
				if ((j = ewalk(e1,deriv1)) > 0)
					vt_free(j);
				}
			condlevel--;
			if (eva->a != nv1) {
				eva->next = (expr_va *)ncond++;
				/* arrange to find this expr_va */
				/* when walking derps */
				dp = eva->R.D;
				LRp = Make_dLRp(dp->c.rp);
				*LRp = LR = (dLR *)mem(sizeof(dLR));
				LR->kind = dLR_VARARG;
				LR->o.eva = eva;
				}
			goto reti;

		case 4: /* piece-wise linear */
			if (deriv) {
				e->dL = 0;
				dvset(&e->dL, 1);
				}
			LR = Make_dLR(&e->dR);
			LR->kind = nplterm++;
			LR->o.ep = Plterms;
			Plterms = e;
 retnew:
			return e->a = new_vt(deriv);

		case 5: /* if */
			eif = (expr_if *)e;
			eif->next = (expr_if*) -1;
			lwalk(&eif->e, 0);
			mts = pdlsave(&opdb);
			if ((i = ewalk(eif->T,deriv1)) > 0)
				vt_free(i);
			pdlreset();
			if ((i = ewalk(eif->F,deriv1)) > 0)
				vt_free(i);
			pdlrestore(opdb, mts);
			if (eif->a != nv1) {
				eif->next = (expr_if *)ncond++;
				/* arrange to find this expr_if */
				/* when walking derps */
				dp = eif->D;
				LRp = Make_dLRp(dp->c.rp);
				*LRp = LR = (dLR *)mem(sizeof(dLR));
				LR->kind = dLR_IF;
				LR->o.eif = eif;
				}
			goto retnew;

		case 6: /* sumlist */
			ep = e->L.ep;
			epe = e->R.ep;
			i = ewalk(*ep++, deriv);
			if (i < 0)
				deriv = deriv1;
			j = ewalk(*ep++, deriv);
			if (i > 0) {
				if (j > 0)
					vt_free(j);
				}
			else if (!(i = j))
				i = new_vt(deriv);
			do {
				if ((j = ewalk(*ep++, deriv1)) > 0)
					vt_free(j);
				}
				while(ep < epe);
			goto reti;

		case 7: /* function call */
			ef = (expr_f *)e;
			i = k = 0;
			for(ap = ef->ap, ape = ef->ape; ap < ape; ap++)
				if (j = ewalk(ap->e,deriv))
					if (j < 0) {
						i = j;
						deriv = deriv1;
						}
					else
						k++;
			if (k)
				for(ap = ef->ap; ap < ape; ap++) {
					op = ap->e->op;
					if (op != f_OPNUM
					 && op != f_OPHOL
					 && op != f_OPVARVAL
					 && (j = ap->e->a) > 0)
						vt_free(j);
					}
			if (!i)
				goto retnew;
			goto reti;

		case 8: /* Hollerith */
			break;

		case 9: /* number */
			t = ((expr_n *)e)->v;
			((expr_nx *)e)->v = t;
			((expr_nx *)e)->next = nums;
			nums = (expr_nx *)e;
			break;

		case 10: /* variable value */
			i = (expr_v *)e - var_e;
			A = Adjp(&adjoints[e->a]);
			A->ifset = condlevel ? 1 : 0;
			if (!A->seen) {
				A->seen = 1;
				vseen[nvseen++] = i;
				}
			break;
		/*DEBUG*/default:
		/*DEBUG*/ fprintf(Stderr, "bad opnumber %d in ewalk\n", k);
		/*DEBUG*/ exit(1);
		}
	return 0;
	}

 static void
#ifdef KR_headers
max_save(m)
	maxinfo *m;
#else
max_save(maxinfo *m)
#endif
{
	if (nvtmax < nvt)
		nvtmax = nvt;
	m->nvt = nvtmax - 1;
	m->ndv = ndvmax;
	m->ncond = ncond;
	m->needT1 = needT1;
	}

 static void
#ifdef KR_headers
max_restore(m)
	maxinfo *m;
#else
max_restore(maxinfo *m)
#endif
{
	nvtmax = m->nvt + 1;
	ndvmax = m->ndv;
	ncond = m->ncond;
	needT1 = m->needT1;
	}

 static void
#ifdef KR_headers
vreset(ifset)
	list **ifset;
#else
vreset(list **ifset)
#endif
{
	Adjoint *A;
	list *lastif = 0;
	int k;

	while(nvseen) {
		k = vseen[--nvseen];
		A = Adjp(&adjoints[var_e[k].a]);
		if (A->ifset)
			(lastif = new_list(lastif))->item.i = k;
		*A = A0;
		}
	*ifset = lastif;
	}

 static void
#ifdef KR_headers
ewalkvt(e, ifset)
	expr *e;
	list **ifset;
#else
ewalkvt(expr *e, list **ifset)
#endif
{
	nvtfree = ndvtmax;
	ewalk(e,want_derivs);
	vreset(ifset);
	}

 static void
#ifdef KR_headers
zset(z)
	register int *z;
#else
zset(register int *z)
#endif
{
	register int *ze;
	if (z)
		for(ze = z + 2**z; z < ze; z += 2)
			var_e[z[1]].a = z[2];
	}

 static void
#ifdef KR_headers
ndvtreset(z)
	int *z;
#else
ndvtreset(int *z)
#endif
{
	nvt = 1;
	nvtfree = ndvtmax;
	ndv = Fortran;
	ndvfree = 0;
	zset(z);
	}

 static void
#ifdef KR_headers
comwalk(i, n)
	int i;
	int n;
#else
comwalk(int i, int n)
#endif
{
	cexp *c, *ce;
	list **ifset;

	if (i >= n)
		return;
	zset(zaC[i]);
	ifset = com_ifset + i;
	c = cexps + i;
	for(ce = cexps + n; c < ce; c++) {
		Lset(c->L, c->nlin);
		ndvtreset(zaC[i]);
		ewalkvt(c->e, ifset++);
		if (nvtmax < nvt)
			nvtmax = nvt;
		}
	}

 static void dwalk0(derp *, derp*);
 static void dwalk1(derp *, derp*);

 static void
#ifdef KR_headers
dwalk0(d, d0)
	derp *d;
	derp *d0;
#else
dwalk0(derp *d, derp *d0)
#endif
{
	int i;

	if (d == d0)
		return;
	do {
		if ((i = d->b.rp - adjoints - nv1) > 0)
			acount[i]++;
		}
		while((d = d->next) != d0);
	}

 static void
#ifdef KR_headers
cond_magic(LR)
	dLR *LR;
#else
cond_magic(dLR *LR)
#endif
{
	expr_va *eva;
	expr_if *eif;
	derp *D, *D0, Dsave;
	de *d;
	int iftype;

	if (LR->kind == dLR_IF) {
		eif = LR->o.eif;
		D = eif->D;
		iftype = 1;
		}
	else {
		eva = LR->o.eva;
		D = eva->R.D;
		iftype = 0;
		}
	Dsave = *D;
	D->c.rp = &edagread_one;
	if (iftype) {
		D->a.rp = eif->Tv.rp;
		D->next = eif->dT;
		dwalk1(D, eif->d0);
		D->a.rp = eif->Fv.rp;
		D->next = eif->dF;
		dwalk1(D, eif->d0);
		}
	else {
		D0 = eva->d0;
		for(d = eva->L.d; d->e; d++) {
			D->a.rp = d->dv.rp;
			D->next = d->d;
			dwalk1(D, D0);
			}
		*D = Dsave;
		}
	*D = Dsave;
	}

#ifdef DEBUG
 /*debug*/static int awrite, derps_in, derps_out, stats, vwrite;
#endif

 static void
#ifdef KR_headers
dwalk1(d, d0)
	derp *d;
	derp *d0;
#else
dwalk1(derp *d,  derp *d0)
#endif
{
	Adjoint *a, *b;
	dLR *c;
	int i, j, j0, neg;

	if (d == d0)
		return;
	dwalk0(d, d0);
	do {
		j = -1;
		if ((j0 = d->b.rp - adjoints - nv1) > 0)
			j = --acount[j0];
		b = Adjp(d->b.rp);
		if (!b->storage)	/* defined var used only in if */
			continue;
		i = d->a.rp - adjoints;
		if (maxa < i)
			maxa = i;
		c = dLRp(*d->c.rp);
		if (c->kind >= dLR_VARARG) {
			if (!output_time)
				cond_magic(c);
			goto jfree;
			}
		if (i < nv1 && !fwalk)
			continue;
		a = Adjp(d->a.rp);
		switch(a->storage) {
			case STOR_VP:
				a->storage = STOR_DV;
				a->o.i = new_dv();
			default:
				goto jfree;
			case 0:
				break;
			}
		if (b->storage == STOR_IMPLICIT) {
			a->neg = b->neg;
			switch(c->kind) {
			  case dLR_negone:
				a->neg = 1 - b->neg;
				/* no break */
			  case dLR_one:
				a->storage = STOR_IMPLICIT;
				break;
			  case dLR_VP:
				a->storage = STOR_VP;
				a->o.vp = c->o.vp;
				break;
			  case dLR_PD:
				a->storage = STOR_PD;
				a->o.i = c->o.i;
				break;
			  case dLR_VARVAL:
				a->storage = STOR_VARVAL;
				a->o.i = c->o.i;
				break;
			  default:
				/*DEBUG*/ fprintf(Stderr,
				"\nBad c->kind = %d with STOR_IMPLICIT\n",
						c->kind);
				/*DEBUG*/ exit(10);
			  }
 jfree:
			if (!j && b->storage == STOR_DV)
				dv_free(b->o.i);
			continue;
			}
		neg = b->neg;
		switch(c->kind) {
		  case dLR_negone:
			neg = 1 - neg;
			/* no break */
		  case dLR_one:
			if (j || fwalk && b->storage == STOR_PD)
				break;
			a->storage = b->storage;
			a->o = b->o;
			a->neg = neg;
			continue;
		  }
		a->storage = STOR_DV;
		a->o.i = !j && b->storage == STOR_DV ? b->o.i : new_dv();
		}
		while((d = d->next) != d0);
	if (b && b->storage == STOR_DV)
		dv_free(b->o.i);
	}

 static void
#ifdef KR_headers
areset(cg)
	cgrad *cg;
#else
areset(cgrad *cg)
#endif
{
	Adjoint *A;
	int k;
	while(nvseen) {
		A = Adjp(&adjoints[var_e[vseen[--nvseen]].a]);
		*A = A0;
		}

	for(k = nv1; k <= maxa; k++) {
		A = Adjp(&adjoints[k]);
		*A = A0;
		}
	maxa = 0;
	for(; cg; cg = cg->next) {
		A = Adjp(&adjoints[cg->varno]);
		*A = A0;
		}
	}

 static void
#ifdef KR_headers
dwalk(d, cg, z, L)
	register derp *d;
	cgrad *cg;
	int *z;
	list *L;
#else
dwalk(register derp *d, cgrad *cg, int *z, list *L)
#endif
{
	register real *rp;
	register derp *d1;
	int i;
	Adjoint *A;
	static cgrad *cg0;

	areset(cg0);
	cg0 = cg;
	ndv = Fortran;
	ndvfree = 0;
	zset(z);
	for(; L; L = L->next) {
		i = L->item.i;
		A = Adjp(&adjoints[var_e[i].a]);
		A->seen = 1;
		vseen[nvseen++] = i;
		A->storage = STOR_DEFV;
		A->o.i = ndv++;
		}
	if (d1 = d) {
		dwalk_one->storage = STOR_IMPLICIT;
		rp = d->b.rp;
		do {
			if (d1->b.rp == rp)
				d1->b.rp = rdwalk_one;
			}
			while(d1 = d1->next);
		dwalk1(d, 0);
		}
	if (ndvmax < ndv)
		ndvmax = ndv;
	}

 static void
#ifdef KR_headers
funnelwalk(f)
	funnel *f;
#else
funnelwalk(funnel *f)
#endif
{
	cplist *cl;
	dLR *d;
	int k;

	if (!want_derivs)
		return;
	for(; f; f = f->next) {
		k = f->ce - cexps;
		dwalk(f->fcde.d, 0, zaC[k], com_ifset[k]);
		for(cl = f->cl; cl; cl = cl->next) {
			d = Make_dLR(cl->cfa);
			d->kind = dLR_PD;
			d->o.i = npd++;
			}
		}
	}

 static void
#ifdef KR_headers
cde_walk(d, n, m, ifset, c1st, z)
	cde *d;
	int n;
	maxinfo *m;
	list **ifset;
	int *c1st;
	int **z;
#else
cde_walk(cde *d, int n, maxinfo *m, list **ifset, int *c1st, int **z)
#endif
{
	cde *De;
	int i, i1, j, k;
	cexp1 *c1;

	j = *c1st++;
	for(De = d + n; d < De; d++, ifset++) {
		i = j;
		k = *c1st++;
		ndvtreset(*z++);
		while(j < k) {
			c1 = cexps1 + j++;
			Lset(c1->L, c1->nlin);
			ewalk(c1->e, want_derivs);
			}
		ewalkvt(d->e, ifset);
		while(i < k)
			if (!cvmap[i1 = i++ + ncom0])
				cvmap[i1] = nvt++;
		if (nvtmax < nvt)
			nvtmax = nvt;
		if (want_derivs)
			dwalk(d->d, 0, 0, *ifset);
		}
	max_save(m);
	}

 static void
#ifdef KR_headers
cde_dwalk(d, n, m, ifset, z, cgp)
	cde *d;
	int n;
	maxinfo *m;
	list **ifset;
	int **z;
	cgrad **cgp;
#else
cde_dwalk(cde *d, int n, maxinfo *m, list **ifset, int **z, cgrad **cgp)
#endif
{
	cde *De;

	if (!want_derivs)
		return;
	for(De = d + n; d < De; d++, ifset++) {
		dwalk(d->d, *cgp++, *z++, *ifset);
		if (m->ndv < ndv)
			m->ndv = ndv;
		}
	}

 static void derivs(derp *, derp *);
 static list **dstored;

 static void
#ifdef KR_headers
dstore(a, d)
	Adjoint *a;
	derp *d;
#else
dstore(Adjoint *a, derp *d)
#endif
{
	a->stored = 1;
	if (dstored)
		(*dstored = new_list(*dstored))->item.D = d;
	}

 static void
derprestore(VOID)
{
	Adjoint *a;
	list *L, *Lnext;

	if (L = *dstored) {
		*dstored = 0;
		do {
			Lnext = L->next;
			a = Adjp(L->item.D->a.rp);
			a->stored = 0;
			L->next = list_freelist;
			list_freelist = L;
			}
			while(L = Lnext);
		}
	}

 static void
#ifdef KR_headers
cond_dmagic(LR)
	dLR *LR;
#else
cond_dmagic(dLR *LR)
#endif
{
	expr_va *eva;
	expr_if *eif;
	derp *D, *D0;
	de *d;
	char cbuf[16];
	int i, j, n;
	list **ds0, *dstore1;

	ds0 = dstored;
	dstored = &dstore1;
	dstore1 = 0;
	if (LR->kind == dLR_IF) {
		eif = LR->o.eif;
		D = eif->D;
		D->c.rp = &edagread_one;
		sprintf(cbuf, cond_fmt, Intcast eif->next);
		if (eif->dT != eif->d0) {
			printf(ifstart_fmt, cbuf, opNE, "0");
			D->a.rp = eif->Tv.rp;
			D->next = eif->dT;
			derivs(D, eif->d0);
			derprestore();
			if (eif->dF != eif->d0) {
				printf(else_fmt);
 just_else:
				D->a.rp = eif->Fv.rp;
				D->next = eif->dF;
				derivs(D, eif->d0);
				derprestore();
				}
			}
		else {
			printf(ifstart_fmt, cbuf, opEQ, "0");
			goto just_else;
			}
#ifdef DEBUG
		--derps_in;
#endif
		printf(endif_fmt);
		D->next = eif->d0;
		}
	else {
		eva = LR->o.eva;
		D = eva->R.D;
		D->c.rp = &edagread_one;
		sprintf(cbuf, cond_fmt, Intcast eva->next);
		for(n = 0, d = eva->L.d; d->e; d++)
			n++;
		i = Switch(cbuf, n);
		j = i + n;
		D0 = eva->d0;
#ifdef DEBUG
		++derps_in;
#endif
		for(d = eva->L.d; ;) {
			Case(i++);
			D->a.rp = d->dv.rp;
			D->next = d->d;
#ifdef DEBUG
			--derps_in;
#endif
			derivs(D, D0);
			derprestore();
			if (!(++d)->e)
				break;
			Break(j);
			}
		endswitch(j);
		}
	dstored = ds0;
	}

 static void
#ifdef KR_headers
derivs(d, d0)
	derp *d;
	derp *d0;
#else
derivs(derp *d, derp *d0)
#endif
{
	Adjoint *a, *b;
	dLR *c;
	char *aop, *fmt, *ginc, *mult;
	int havenum, i, neg, num2chk, stored;
	derp *dnext;
	char bufb[32], bufc[32], bufg[32], bufi[32];
	double t, t1;
	v_i *vi;

	for(; d != d0; d = dnext) {
		dnext = d->next;
		c = dLRp(*d->c.rp);
		if (c->kind >= dLR_VARARG) {
			cond_magic(c);
			cond_dmagic(c);
			continue;
			}
#ifdef DEBUG
		++derps_in;
#endif
		b = Adjp(d->b.rp);
		if (!b->storage)	/* defined var used only in if */
			continue;
		neg = b->neg;
		a = Adjp(d->a.rp);
		if (c->kind == dLR_one
		 && a->storage == b->storage && a->o.i == b->o.i) {
			switch(a->storage) {
		  		case STOR_GRAD:
		 		case STOR_JAC:
				case STOR_DEFV:
					break;
				default:
					dstore(a,d);
					a->neg = neg;
					continue;
				}
			}
		num2chk = 0;
		stored = a->stored;
		switch(a->storage) {
		  case STOR_VI:
			i = a->o.vi->i;
			fmt = dv_fmt;
			break;
		  case STOR_PD:
			if (b->storage == STOR_IMPLICIT
			 && (stor_grad != STOR_DEFV
				|| a->storage == STOR_PD
					&& c->kind == dLR_PD
					&& a->o.i == c->o.i))
				continue;
			stored = a->stored = 1;
			if (a->neg)
				neg = 1 - neg;
			fmt = pd_fmt;
			i = a->o.i + Fortran;
			break;
		  case STOR_DV:
		  case STOR_DEFV:
			num2chk = 1;
			i = a->o.i;
			fmt = dv_fmt;
			break;
		  case STOR_GRAD:
			i = a->o.og->varno + Fortran;
			fmt = grad_fmt;
			break;
		  case STOR_JAC:
			i = (djoff >= 0
				? djoff + a->o.cg->varno*n_con
				: a->o.cg->goff) + Fortran;
			fmt = jac_fmt;
			break;
		  default:
			continue;
			}
		mult = "*";
		havenum = 0;
		switch(b->storage) {
			case STOR_PD:
				sprintf(bufb, pd_fmt, b->o.i + Fortran);
				break;
			case STOR_DV:
			case STOR_DEFV:
				sprintf(bufb, dv_fmt, b->o.i);
				break;
			case STOR_VARVAL:
				sprintf(bufb, x_fmt, b->o.i + Fortran);
				break;
			case STOR_VI:
				t = b->o.vi->u.v;
				goto have_t;
			case STOR_VP:
				t = *b->o.vp;
 have_t:
				havenum = 1;
				if (t < 0) {
					t = -t;
					neg = 1 - neg;
					}
				g_fmt(bufb, t);
				break;
			case STOR_IMPLICIT:
				havenum = 1;
				t = 1.;
				mult = "";
				bufb[0] = 0;
				break;
			default:
				/*DEBUG*/ fprintf(Stderr,
				"\nBad b->storage = %d\n", b->storage);
				/*DEBUG*/ exit(12);
				}
		switch(c->kind) {
			case dLR_VP:
				if (havenum) {
					t *= *c->o.vp;
					if (!stored && num2chk) {
						havenum = 2;
						break;
						}
					mult = "";
					bufb[0] = 0;
					}
				else
					t = *c->o.vp;
 num_fmt:
				if (t < 0) {
					t = -t;
					neg = 1 - neg;
					}
				g_fmt(bufc, t);
				break;
			case dLR_PD:
				sprintf(bufc, pd_fmt, c->o.i + Fortran);
				break;
			case dLR_VARVAL:
				sprintf(bufc, x_fmt, c->o.i + Fortran);
				break;
			case dLR_negone:
				neg = 1 - neg;
			case dLR_one:
				if (havenum) {
					if (!stored) {
						havenum = 2;
						if (*mult)
							break;
						}
					goto num_fmt;
					}
				if (*mult) {
					mult = "";
					bufc[0] = 0;
					}
				else
					strcpy(bufc, One);
				break;
			default:
				/*DEBUG*/fprintf(Stderr,
				"\nBad c->kind = %d\n", c->kind);
				/*DEBUG*/ exit(11);
			}
		ginc = 0;
		sprintf(bufg, fmt, i);
		if (stored) {
			aop = neg ? "-" : "+";
			if (Fortran)
				printf("\t%s = %s %s ", bufg, bufg, aop);
			else
				printf("\t%s %s= ", bufg, aop);
			}
		else {
			if (neg && havenum)
				t = -t;
			switch(a->storage) {
				case STOR_DV:
				case STOR_DEFV:
					if (havenum == 2) {
						a->storage = STOR_VI;
						a->o.vi = new_vi(t, a->o.i);
						a->neg = 0;
						continue;
						}
					break;
				case STOR_VI:
					if (havenum == 2) {
						a->o.vi->u.v += t;
						continue;
						}
					a->storage = STOR_DV;
					a->o.i = (vi = a->o.vi)->i;
					if (vi->u.v)
						g_fmt(ginc = bufi, vi->u.v);
					sprintf(bufg, dv_fmt, vi->i);
					free_vi(vi);
					break;
				case STOR_GRAD:
					t1 = a->o.og->coef;
					goto t1_test;
				case STOR_JAC:
					t1 = a->o.cg->coef;
 t1_test:
					if (t1)
						g_fmt(ginc = bufi, t1);
				}
			dstore(a,d);
			if (!*mult && !neg
			 && (!bufb[0] && !strcmp(bufc,bufg))
			  || !bufc[0] && !strcmp(bufb,bufg))
				continue;
			if (!strcmp(bufg,bufb) && !Fortran
			 && *mult == '*' && !neg && !ginc) {
				printf("\t%s *= %s;\n", bufg, bufc);
#ifdef DEBUG
				++derps_out;
				if (awrite)
					goto a_write;
#endif
				continue;
				}
			printf("\t%s %s", bufg, neg ? "= -" : "= ");
			}
		printf("%s%s%s", bufb, mult, bufc);
		if (ginc)
			printf(" + %s", ginc);
		printf(eos);
#ifdef DEBUG
		++derps_out;
		if (awrite) {
 a_write:
		if (!*bufb)
			strcpy(bufb, One);
		if (!*bufc)
			strcpy(bufc, One);
		  if (Fortran) {
			printf("\twrite(*,*) '%d: %s', %s, ' * ', %s%s",
				derps_in, neg ? "-(" : "", bufb, bufc,
				neg ? ",')'" : "");
			if (ginc)
				printf(", ' +  %s'", ginc);
			printf(", ': a[%d] = ', %s\n",
				d->a.rp - adjoints, bufg);
		  } else {
			printf("fprintf(Stderr, \"%d: %s", derps_in,
				neg ? "-(%g * %g)" : "%g * %g");
			if (ginc)
				printf(" + %s", ginc);
			printf(": a[%d] = %%g\\n\", %s, %s, %s);\n",
				d->a.rp - adjoints, bufb, bufc, bufg);
		  }
			}
#endif
		}
	}

 static void
#ifdef KR_headers
co_derivs(d, L, gj0_fmt)
	derp *d;
	list *L;
	char *gj0_fmt;
#else
co_derivs(derp *d, list *L, char *gj0_fmt)
#endif
{
	char buf[32];
	Adjoint *a;
	int j, k;
	ograd *og;

	for(j = Fortran; L; L = L->next, j++) {
		a = Adjp(&adjoints[var_e[k = L->item.i].a]);
		if (k < nv1)
			if (stor_grad == STOR_DEFV) {
				sprintf(buf, dv_fmt, j + Fortran);
				assign(buf, Zero);
				}
			else {
				og = stog[k];
				g_fmt(buf, stor_grad == STOR_GRAD ? og->coef
						: ((cgrad *)og)->coef);
				printf(gj0_fmt, k + Fortran, buf);
				}
		else {
			a->storage = STOR_DEFV;
			a->o.i = j;
			sprintf(buf, dv_fmt, j);
			assign(buf, Zero);
			}
		a->stored = 1;
		}
	derivs(d,0);
	}

 static char *
#ifdef KR_headers
vprod(t, k)
	real t;
	int k;
#else
vprod(real t, int k)
#endif
{
	static char buf[64];
	int i;

	k += Fortran;
	if (t == 1.)
		sprintf(buf, x_fmt, k);
	else if (t == -1.) {
		buf[0] = '-';
		sprintf(buf+1, x_fmt, k);
		}
	else {
		i = g_fmt(buf, t);
		buf[i++] = '*';
		sprintf(buf+i, x_fmt, k);
		}
	return buf;
	}

 static void
#ifdef KR_headers
plcommon(deriv)
	int deriv;
#else
plcommon(int deriv)
#endif
{
	register plterm *p;
	register expr *e;
	register dLR *LR;

	for(e = Plterms; e; e = LR->o.ep) {
		p = e->L.p;
		LR = dLRp(e->dR);
		printf("\tcommon /pltcom/ bs%d\n\
	double precision bs%d(%d)\n", LR->kind, LR->kind, 2*p->n - 1);
		}
	printf(
	 "\tcommon /xkindc/ xkind\n\tinteger xkind\n\tsave /xkindc/\n");

	if (npd)
		printf(
	"\tcommon /pdcomn/ pd\n\tdouble precision pd(%d)\n\tsave /pdcomn/\n",
			npd);
	if (!deriv && needx0check)
		printf("\tlogical xchk\n");
	}

 static char *
#ifdef KR_headers
rv_output(rv, eval, og)
	char *rv;
	char *eval;
	ograd *og;
#else
rv_output(char *rv, char *eval, ograd *og)
#endif
{
	char *s;

	for(; og; og = og->next)
		if (og->coef)
			break;
	if (og) {
		s = vprod(og->coef, og->varno);
		if (strcmp(eval, Zero))
			binop(rv, eval, "+", s);
		else
			assign(rv, s);
		while(og = og->next)
			if (og->coef)
				binop(rv, rv, "+",
					vprod(og->coef, og->varno));
		return rv;
		}
	if (Fortran)
		assign(rv, eval);
	return eval;
	}

 static void
#ifdef KR_headers
obj_derivs(i)
	int i;
#else
obj_derivs(int i)
#endif
{
	ograd *og, *og0;
	Adjoint *a;
	cde *c;
	char buf[32];

	if (!want_derivs)
		return;
	c = obj_de + i;
	og0 = Ograd[i];
	dwalk(c->d, (cgrad*)og0, zao[i], o_ifset[i]);	/* added 20010312 */
	for(og = og0; og; og = og->next) {
		a = Adjp(&adjoints[og->varno]);
		a->storage = STOR_GRAD;
		a->o.og = og;
		a->stored = 0;
		stog[og->varno] = og;
		}
	stor_grad = STOR_GRAD;
	co_derivs(c->d, o_ifset[i], g0_fmt);
	for(og = og0; og; og = og->next) {
		a = Adjp(&adjoints[og->varno]);
		if (!a->stored && (!intsk || !intsk[og->varno])) {
			g_fmt(buf, og->coef);
			printf(g0_fmt, og->varno + Fortran, buf);
			}
		}
	}

 static char *
#ifdef KR_headers
con_linadd(i, s)
	int i;
	char *s;
#else
con_linadd(int i, char *s)
#endif
{
	cgrad *cg;
	char *s1;

	for(cg = Cgrad[i]; cg; cg = cg->next)
		if (cg->coef) {
			s1 = vprod(cg->coef, cg->varno);
			if (strcmp(s,Zero))
				binop(T, s, "+", s1);
			else
				assign(T, s1);
			s = T;
			while(cg = cg->next)
				if (cg->coef)
					binop(s, s, "+",
						vprod(cg->coef, cg->varno));
			break;
			}
	return s;
	}

 static void
#ifdef KR_headers
con_derivs(i)
	int i;
#else
con_derivs(int i)
#endif
{
	cgrad *cg, *cg0;
	Adjoint *a;
	char buf[32];
	int i1, j;
	cde *c = con_de + i;

	if (densejac)
		djoff = i;
	printf("\n%s   /*** derivatives for constraint %d ***/\n\n",
		Fortstar, i+1);
	dwalk(c->d, cg0 = Cgrad[i], zac[i], c_ifset[i]);
	for(cg = cg0; cg; cg = cg->next) {
		a = Adjp(&adjoints[cg->varno]);
		a->storage = STOR_JAC;
		a->stored = 0;
		a->o.cg = cg;
		stog[cg->varno] = (ograd *)cg;
		}
	stor_grad = STOR_JAC;
	co_derivs(con_de[i].d, c_ifset[i], j0_fmt);
	for(cg = cg0; cg; cg = cg->next) {
		a = Adjp(&adjoints[cg->varno]);
		if (!a->stored && (!intsk || !intsk[cg->varno])) {
			g_fmt(buf, cg->coef);
			j = densejac ? i + cg->varno*n_con : cg->goff;
			printf(j0_fmt, j + Fortran, buf);
			}
		}
	if (densejac) {
		i1 = i + 1;
		for(cg = cg0; cg; cg = cg->next)
			rownos[cg->varno] = i1;
		for(j = 0; j < nv1; j++)
			if (rownos[j] != i1)
				printf(j0_fmt, i + j*n_con + Fortran, Zero);
		}
	}

 static int
nzcgrad(VOID)
{
	cgrad *cg;
	int i;
	for(i = 0; i < n_con; i++)
		for(cg = Cgrad[i]; cg; cg = cg->next)
			if (cg->coef)
				return 1;
	return 0;
	}

 static char *
#ifdef KR_headers
cv_name(L, buf)
	linpart *L;
	char *buf;
#else
cv_name(linpart *L, char *buf)
#endif
{
	expr_v ev;
	expr *ep;

	ep = (expr *)((char *)L->v.rp - ((char *)&ev.v - (char *)&ev));
	return f_OPVARVAL1(ep, buf);
	}

 static void
#ifdef KR_headers
com_out(e, L, nlin, k0)
	expr *e;
	linpart *L;
	int nlin;
	int k0;
#else
com_out(expr *e, linpart *L, int nlin, int k0)
#endif
{
	char buf[32], bufg[32], res[32], vn[32];
	char *s;
	double t;
	linpart *Le;
	int asg, j, k, op;
	efuncb *eb;
	dLR *d;

	printf("\n%s\t/*** defined variable %d ***/\n\n", star, k0+1);
	if ((j = cvmap[k0]) < 0) {
		j = Fortran1 - j;
		s = pd_fmt;
		}
	else {
		eb = (efuncb *)e->op;
		if (eb != (efuncb *)OPNUM && eb != (efuncb *)OPVARVAL)
			e->a = j;
		j += Fortran1;
		s = tv_fmt;
		}
	sprintf(res, s, j);
	s = callb(e,buf);
	if (!L) {
		if (strcmp(res,s))
			assign(res, s);
#ifdef DEBUG
		goto done;
#else
		return;
#endif
		}
	asg = !strcmp(s, Zero);
	for(Le = L + nlin; L < Le; L++, asg = 0) {
		d = dLRp(L->fac);
		op = '+';
		switch(k = d->kind) {
			case dLR_negone:
				op = '-';
			case dLR_one:
				break;
			case dLR_VP:
				t = *d->o.vp;
				if (t < 0. && !asg) {
					t = -t;
					op = '-';
					}
				g_fmt(bufg, t);
				break;
			default:/*DEBUG*/ fprintf(Stderr,
				"Bad d->kind = %d in com_walk\n", d->kind);
				/*DEBUG*/ exit(14);
			}
		if (asg)
			printf(op == '-' ? "\t%s = -" : "\t%s = ", res);
		else if (Fortran || res != s)
			printf("\t%s = %s %c ", res, s, op);
		else
			printf("\t%s %c= ", res, op);
		if (k == dLR_VP)
			printf("%s*", bufg);
		printf("%s%s", cv_name(L,vn), eos);
		s = res;
		}
#ifdef DEBUG
 done:
	if (vwrite)
	printf("fprintf(Stderr, \"defvar[%d] = %%g\\n\", %s);\n", k0, res);
#endif
	}

 static char *
#ifdef KR_headers
putout(e, i, j, what, k, z)
	expr *e;
	int i;
	int j;
	char *what;
	int k;
	int *z;
#else
putout(expr *e, int i, int j, char *what, int k, int *z)
#endif
{
	static char buf[32];
	cexp1 *c1;

	++k;
	ndvtreset(z);
	if (i < j) {
		if (what)
		printf("\n\n%s\t/*** defined variables for %s %d ***/\n",
			star, what, k);
		for(c1 = cexps1 + i; i < j; i++, c1++)
			com_out(c1->e, c1->L, c1->nlin, i + ncom0);
		}
	printf(what ? "\n%s  /***  %s %d  ***/\n\n"
		    : "\n%s  /***  objective ***/\n\n", star, what, k);
	return callb(e,buf);
	}

 static void
#ifdef KR_headers
x0check(i, j, k, deriv)
	int i;
	int j;
	int k;
	int deriv;
#else
x0check(int i, int j, int k, int deriv)
#endif
{
	cexp *c;

	if (!needx0check)
		return;
	printf(deriv ? xcheck : xcheck0);
	if (i < j) {
		printf(xkind, k, k, k);
		for(c = cexps + i; i < j; i++, c++)
			com_out(c->e, c->L, c->nlin, i);
		printf("\t%s", endif_fmt);
		}
	}

 static void
#ifdef KR_headers
funnel_set(f, gj0_fmt)
	funnel *f; char *gj0_fmt;
#else
funnel_set(funnel *f, char *gj0_fmt)
#endif
{
	cplist *cl;
	dLR *d;
	Adjoint *a;
	int k;
	list *Lp;
	char *b, *b1, buf[32], bufb[32];

	stor_grad = STOR_DEFV;
	fwalk = 1;
	do {
		printf("\n%s\t/*** funnel ***/\n\n", Fortstar);
		k = f->ce - cexps;
		dwalk(f->fcde.d, 0, zaC[k], Lp = com_ifset[k]);
		co_derivs(f->fcde.d, Lp, gj0_fmt);
		cl = f->cl;
		do {
			d = dLRp(*cl->cfa); /*!! was Make_dLR(cl->cfa) */
			sprintf(buf, pd_fmt, d->o.i + Fortran);
			a = Adjp(cl->ca.rp);
			b = bufb;
			switch(a->storage) {
			  case STOR_IMPLICIT:
				b = a->neg ? "-1" : "1";
				break;
			  case STOR_VP:
				g_fmt(b, *a->o.vp);
				break;
			  case STOR_PD:
				b1 = b;
				if (a->neg)
					*b1++ = '-';
				sprintf(b1, pd_fmt, a->o.i + Fortran);
				break;
			  case STOR_DV:
			  case STOR_DEFV:
				sprintf(b, dv_fmt, a->o.i);
				break;
			  default:/*DEBUG*/ fprintf(Stderr,
				"Bad a->storage in funnel_set\n");
				/*DEBUG*/exit(15);
				}
			*a = A0;
			assign(buf, b);
			}
			while(cl = cl->next);
		}
		while(f = f->next);
	fwalk = 0;
	}

 static void
end(VOID)
{
	printf(Fortran ? "\tend\n" : "}\n");
	branches = 0;
	}

 void
#ifdef KR_headers
fobj_output(deriv)
	int deriv;
#else
fobj_output(int deriv)
#endif
{
	int *c1, i, lbl;
	char buf[24], *eval;
	ograd *og;
	static char *rv[2] = { "feval0", "feval" };

	sprintf(buf, n_obj ? "%d" : "*", o_vars);
	printf(deriv ?
	"\n\tdouble precision function feval(nobj, needfg, x, g)\n\
	integer nobj, needfg\n\
	double precision x(%s), g(%s)\n"
	: "\n\tdouble precision function feval0(nobj, x)\n\
	integer nobj\n\
	double precision x(%s)\n", buf, buf);

	if (n_obj > 1)
		printf("\tinteger nobj1\n");
	else if (!n_obj) {
		printf("\tend\n");
		return;
		}

	want_derivs = deriv;
	for(og = 0, i = 0; i < n_obj; i++) {
		for(og = Ograd[i]; og; og = og->next)
			if (og->coef)
				goto break2;
		}
 break2:
	cde_dwalk(obj_de, n_obj, &omax, o_ifset, zao, (cgrad**)Ograd);
	if (needx0check)
		printf("%s\tlogical xcheck\n\texternal xcheck\n",
			deriv ? "\tinteger wantfg\n" : "");
	if (omax.nvt) {
		printf("\n\tdouble precision v(%d)", omax.nvt);
		if (omax.ndv > 1 && deriv)
			printf(", dv(%d)", omax.ndv-1);
		printf(og ? ", rv\n" : "\n");
		if (omax.ncond > 1 && want_derivs)
			printf("\tinteger cond(%d)\n", omax.ncond-1);
		}
	else if (og)
		printf("\n\tdouble precision rv\n");

	if (omax.needT1)
		printf("\tdouble precision t1, t2\n");

	plcommon(deriv);

	x0check(combc, ncom0, 2, deriv);

	if ((f_b || f_o) && deriv) {
		printf("\tif (wantfg .ge. 2) then\n");
		if (f_b)
			printf("\t\tcall funelb(x)\n");
		if (f_o)
			funnel_set(f_o, g0_fmt);
		printf("\t\tendif\n");
		}

	if (n_obj > 1) {
		printf("\tnobj1 = nobj + 1\n");
		lbl = Switch("nobj1", n_obj);
		}
	c1 = o_cexp1st;
	for(i = 0;; c1++) {
		if (n_obj > 1)
			Case(lbl+i);
		if (deriv)
			printf("\tif (mod(wantfg,2) .eq. 1) then\n");
		eval = putout(obj_de[i].e, c1[0], c1[1],
				n_obj > 1 ? "objective" : NULL, i, zao[i]);
		rv_output(rv[deriv], eval, Ograd[i]);
		if (deriv) {
			printf("\n\tendif\n\n\tif (wantfg .gt. 1) then\n");
			obj_derivs(i);
			printf("\n\tendif\n");
			}
		if (++i >= n_obj)
			break;
		printf("\treturn\n");
		}
	end();
	}

 void
#ifdef KR_headers
obj_output(deriv)
	int deriv;
#else
obj_output(int deriv)
#endif
{
	int *c1, i;
	ograd *og;
	static char rv[] = "rv";
	static char *header[4] = {
		"0_(fint *nobj, real *x)",
		"_(fint *nobj, fint *needfg, real *x, real *g)",
		"0_(nobj, x)\n\tfint *nobj;\n real *x;",
		"_(nobj, needfg, x, g)\n fint *nobj, *needfg;\n real *x, *g;"
		};
	char *eval, *s;

	printf(" real\nfeval%s\n", header[deriv + krflag]);

	if (!n_obj) {
		/*{*/ printf("{ return 0.; }\n");
		return;
		}
	want_derivs = deriv;
	printf("{");
	for(og = 0, i = 0; i < n_obj; i++) {
		for(og = Ograd[i]; og; og = og->next)
			if (og->coef)
				goto break2;
		}
 break2:
	cde_dwalk(obj_de, n_obj, &omax, o_ifset, zao, (cgrad**)Ograd);
	if (omax.nvt) {
		printf("\n\treal v[%d]", omax.nvt);
		if (omax.ndv)
			printf(", dv[%d]", omax.ndv);
		printf("%s;\n", og ? ", rv" : "");
		}
	else if (og)
		printf("\n\treal rv;\n");

	if (omax.ncond)
		printf("\tstatic int cond[%d];\n", omax.ncond);

	if (omax.needT1)
		printf("\treal t1, t2;\n");

	x0check(combc, ncom0, 2, deriv);

	if ((f_b || f_o) && deriv) {
		printf("\tif (wantfg & 2) {\n");
		if (f_b)
			printf("\t\tfunnelb(x);\n");
		if (f_o)
			funnel_set(f_o, g0_fmt);
		printf("\t\t}\n");
		}

	printf(n_obj > 1 ? "\n\tswitch(*nobj) {\n" : "\n");
	c1 = o_cexp1st;
	for(i = 0; i < n_obj; i++, c1++) {
		if (n_obj > 1)
			printf("\n  case %d:\n", i);
		if (deriv)
			printf("\tif (wantfg & 1) {\n");
		eval = putout(obj_de[i].e, c1[0], c1[1],
			n_obj > 1 ? "objective" : NULL, i, zao[i]);
		s = rv_output(rv, eval, Ograd[i]);
		if (deriv) {
			printf("\t;}\n\n\tif (wantfg & 2) {\n");
			obj_derivs(i);
			printf("\t}\n");
			}
		printf("\n\treturn %s;\n", s);
		}
	if (n_obj > 1)
		printf("\n\t}\n");
	end();
	}

 void
#ifdef KR_headers
fcon_output(deriv)
	int deriv;
#else
fcon_output(int deriv)
#endif
{
	int *c1, i;
	char *s;

	printf(deriv ? "\n\tsubroutine ceval(needfg, x, c, J)\n\
	integer needfg\n\tdouble precision x(%d), c(%d), J(%ld)\n"
	: "\n\tsubroutine ceval0(x, c)\n\
	double precision x(%d), c(%d)\n",
		c_vars ? c_vars : 1, n_con ? n_con : 1, nJ ? nJ : 1);

	if (!n_con) {
		printf("\tend\n");
		return;
		}

	want_derivs = deriv;
	cde_dwalk(con_de, n_con, &cmax, c_ifset, zac, Cgrad);
	if (needx0check)
		printf("%s\tlogical xcheck\n\texternal xcheck\n",
			deriv ? "\tinteger wantfg\n" : "");

	if (cmax.nvt) {
		printf("\n\tdouble precision v(%d)", cmax.nvt);
		printf(cmax.ndv > 1  && deriv ? ", dv(%d)\n" : "\n",
			cmax.ndv-1);
		}
	if (cmax.ncond > 1)
		printf("\tinteger cond(%d)\n", cmax.ncond-1);

	if (cmax.needT1)
		printf("\tdouble precision t1, t2\n");
	else if (nzcgrad())
		printf("\tdouble precision t1\n");

	plcommon(deriv);

	x0check(comb, combc, 1, deriv);

	if (deriv)
		printf("\n\tif (mod(wantfg,2) .gt. 0) then\n");
	c1 = c_cexp1st;
	for(i = 0; i < n_con; i++, c1++) {
		s = putout(con_de[i].e, c1[0], c1[1], "constraint", i, zac[i]);
		printf("\tc(%d) = %s\n", i+1, con_linadd(i,s));
		}
	if (deriv) {
		printf("\tendif\n\n\tif (wantfg .gt. 1) then\n");
		if (f_b)
			printf("\t\tcall funelb(x)\n");
		if (f_c)
			funnel_set(f_c, j0_fmt);
		for(i = 0; i < n_con; i++)
			con_derivs(i);
		printf("\tendif\n");
		djoff = -1;
		}
	end();
	}

 void
#ifdef KR_headers
con_output(deriv)
	int deriv;
#else
con_output(int deriv)
#endif
{
	int *c1, i;
	char *s;
	static char *header[4] = {
		"0_(real *x, real *c)",
		"_(fint *needfg, real *x, real *c, real *J)",
		"0_(x, c)\n real *x, *c;",
		"_(needfg, x, c, J)\n fint *needfg;\n real *x, *c, *J;"
		};
	printf("\n void\nceval%s\n{", header[deriv + krflag]);

	if (!n_con) {
		printf("}\n");
		return;
		} /*}*/

	want_derivs = deriv;
	cde_dwalk(con_de, n_con, &cmax, c_ifset, zac, Cgrad);
	if (cmax.nvt) {
		printf("\n\treal v[%d]", cmax.nvt);
		printf(cmax.ndv && deriv ? ", dv[%d];\n" : ";\n", cmax.ndv);
		}
	if (cmax.ncond)
		printf("\tstatic int cond[%d];\n", cmax.ncond);

	if (cmax.needT1)
		printf("\treal t1, t2;\n");
	else if (nzcgrad())
		printf("\treal t1;\n");

	x0check(comb, combc, 1, deriv);
	if (deriv)
		printf("\n\tif (wantfg & 1) {\n");
	c1 = c_cexp1st;
	for(i = 0; i < n_con; i++, c1++) {
		s = putout(con_de[i].e, c1[0], c1[1], "constraint", i, zac[i]);
		printf("\tc[%d] = %s;\n", i, con_linadd(i,s));
		}
	if (deriv) {
		printf("\t;}\n\tif (wantfg & 2) {\n");
		if (f_b)
			printf("\tfunnelb(x);\n");
		if (f_c)
			funnel_set(f_c, j0_fmt);
		for(i = 0; i < n_con; i++)
			con_derivs(i);
		printf("\t}\n");
		djoff = -1;
		}
	end();
	}

 void
output(VOID)
{
	int i, j;
	plterm *p;
	expr *e;
	real *b, *be;
	dLR *LR;
	char buf[32], *x0;
	cexp *c;
	static char *fhead[2] = { "(real *x)", "(x) real *x;" };

#ifdef DEBUG
	if (awrite | vwrite)
		printf("#include \"stdio.h\"\n");
#endif
	printf("#include \"math.h\"\n#include \"errno.h\"\n\
#ifndef fint\n\
#ifndef Long\n\
#include \"arith.h\"	/* for Long */\n\
#ifndef Long\n\
#define Long long\n\
#endif\n\
#endif\n\
#define fint Long\n\
#endif\n\
#ifndef real\n\
#define real double\n\
#endif\n");
	if (!krflag)
		printf("#ifdef __cplusplus\nextern \"C\" {\n#endif\n");
	if (krflag < 2)
		printf(" %s\n %s\n %s\n %s\n %s\n %s\n %s\n %s\n",
			"real acosh_(real *);",
			"real asinh_(real *);",
			"real acoshd_(real *, real *);",
			"real asinhd_(real *, real *);",
			"void in_trouble(char *, real);",
			"void in_trouble2(char *, real, real);",
			"void domain_(char *, real *, fint);",
			"void zerdiv_(real *);");
	printf(" fint auxcom_[1] = { %d /* nlc */ };\n", nlc);
	printf(" fint funcom_[%d] = {\n\
	%d /* nvar */,\n\
	%d /* nobj */,\n\
	%d /* ncon */,\n\
	%d /* nzc */,\n\
	%d /* densejac */",
		n_con && !densejac ? nzc + nv1 + n_obj + 6 : n_obj + 5,
		nv1, n_obj, n_con, densejac ? nv1*n_con : nzc, densejac);

	if (n_obj) {
		printf(",\n\n\t/* objtype (0 = minimize, 1 = maximize) */\n");
		for(i = 0; i < n_obj; i++)
			printf("%s\n\t%d", i ? "," : "", objtype[i]);
		}

	if (n_con && !densejac) {
		printf(",\n\n\t/* colstarts */\n");
		for(i = 0; i <= n_var; i++)
			printf("\t%d,\n", A_colstarts[i]);
		printf("\n\t/* rownos */\n\t1");
		for(i = 1; i < nzc; i++)
			printf(",\n\t%d", rownos[i]);
		}
	printf(" };\n\n");

	for(i = j = 0; i < N_OPS; i++)
		if (seen[i] && declare[i])
			printf(j++ ? ", %s_()" : " extern real %s_()",
				declare[i]);
	if (j)
		printf(";\n");
	for(e = Plterms; e; e = LR->o.ep) {
		LR = dLRp(e->dR);
		p = e->L.p;
		i = 2*p->n - 1;
		printf(" real bs%d[%d] = {\n", LR->kind, i);
		for(b = p->bs, be = b + i; b < be; b++) {
			g_fmt(buf, *b);
			printf("\t%s%s\n", buf,
				b + 1 < be ? "," : "};\n");
			}
		}

	printf(
" real boundc_[1+%d+%d] /* Infinity, variable bounds, constraint bounds */ = {\n\t\t1.7e308",
		2*nv1, 2*n_con);
	b = bounds;
	be = b + 2*(n_con + nv1);
	while(b < be)
		printf(",\n\t\t%s", fpval(*b++));
	printf("};\n\n");

	printf(" real x0comn_[%d] = {\n", nv1);
	for(i = 0; i < nv1; i++)
		printf("%s\t\t%s", i ? ",\n" : "", fpval(X0[i]));
	printf(" };\n\n");

	if (npd)
		printf(" static real pd[%d];\n", npd);

	if (f_b && derkind & 2) {
		printf("\n static void\nfunnelb%s\n{\n", fhead[krflag>>1]);
		if ((i = bmax.ndv) > 0)
			printf("\treal dv[%d];\n", i);
		funnel_set(f_b, "Botch<%d> = %s;");
		printf("\t}\n");
		}

	if (needx0check) {
		printf("static real old_x[%d];\nstatic int xkind = -1;\n\n\
 static int\nxcheck%s\n{\n\treal", nv1, xcheckdcl);
		if (comb > 0) {
			printf(" *x0 = x,");
			x0 = "x0";
			}
		else
			x0 = "x";
		printf(" *x1 = old_x, *xe = x + %d;\n", nv1);
		if ((i = bmax.nvt) > 0)
			printf("\treal v[%d];\n", i);
		printf("\terrno = 0;\n\
	if (xkind >= 0) {\n\t\twhile(*%s++ == *x1++)\n\
		\tif (%s == xe)\n\t\t\t\treturn 0;\n\t\t--%s, --x1;\n\t\t}\n\
	do *x1++ = *%s++;\n\t\twhile(%s < xe);\n\txkind = 0;\n",
			x0,x0,x0,x0,x0);
		for(i = 0, c = cexps; i < comb; c++, i++)
			com_out(c->e, c->L, c->nlin, i);
		printf("\treturn 1;\n\t}\n");
		}
	for(i = 1; i < 3; i++) {
		if (owant & i)
			obj_output(i-1);
		if (cwant & i)
			con_output(i-1);
		}
	if (!krflag)
		printf("#ifdef __cplusplus\n\t}\n#endif\n");
	}

 void
foutput(VOID)
{
	int colrow, i, j, nb, ncr;
	plterm *p;
	real *b, *be;
	expr *e;
	dLR *LR;
	char buf[32];
	cexp *c;

	ncr = n_obj;
	if (colrow = n_con && !densejac)
		ncr += nzc + n_var + 1;
	printf(
	 "\tblock data\n\tcommon /funcom/ nvar, nobj, ncon, nzc, densej%s\n",
		ncr ? ", colrow" : "");
	printf("\tinteger nvar, nobj, ncon, nzc, densej");
	if (ncr)
		printf(", colrow(%d)", ncr);

	nb = 2*(n_con + o_vars);
	printf("\n\n\tcommon /boundc/ bounds\n\
	double precision bounds(%d)\n\
	common /x0comn/ x0\n\tdouble precision x0(%d)\n",
		nb + 1, nv1);
	printf("\tcommon /auxcom/ nlc\n\tinteger nlc\n");
	plcommon(1);

	printf("\tdata nvar/%d/, nobj/%d/, ncon/%d/, nzc/%d/, densej/%d/\n",
		nv1, n_obj, n_con, colrow ? nzc : n_con*nv1, densejac);
	if (n_obj) {
		printf("\n*\t*** objtype (0 = minimize, 1 = maximize) ***\n");
		for(i = 0; i < n_obj; i++)
			printf("\tdata colrow(%d)/%d/\n", i+1, objtype[i]);
		}
	if (colrow) {
		printf("\n*\t*** colstarts ***\n");
		j = n_obj + 1;
		for(i = 0; i <= n_var; i++)
			printf("\tdata colrow(%d)/%d/\n", i+j, A_colstarts[i]+1);
		j += n_var + 1;
		printf("\n*\t*** rownos ***\n\tdata colrow(%d)/1/\n", j);
		for(i = 1; i < nzc; i++)
			printf("\tdata colrow(%d)/%d/\n", i+j, rownos[i]);
		}
	for(i = 0; i < nv1; i++)
		printf("\tdata x0(%d)/%s/\n", i+1, fpval(X0[i]));
	for(e = Plterms; e; e = LR->o.ep) {
		LR = dLRp(e->dR);
		p = e->L.p;
		i = 2*p->n - 1;
		j = 0;
		for(b = p->bs, be = b + i; b < be; b++) {
			g_fmt(buf, *b);
			printf("\tdata bs%d(%d)/%s/\n", LR->kind, ++j, buf);
			}
		}

	b = bounds;
	be = b + nb;
	i = 2;
	printf("\tdata bounds(1)/1.7e+308/\n");
	while(b < be)
		printf("\tdata bounds(%d)/%s/\n", i++, fpval(*b++));
	printf("\tdata nlc/%d/\n", nlc);
	printf("\tdata xkind/-1/\n\tend\n\n");

	if (f_b && derkind & 2) {
		printf("\tsubroutine funelb(x)\n");
		printf("\tdouble precision x(%d)\n", n_var);
		plcommon(1);
		if ((i = bmax.ndv) > 0)
			printf("\tdouble precision dv(%d)\n", i);
		funnel_set(f_b, "Botch<%d> = %s");
		printf("\tend\n");
		}

	if (needx0check) {
		printf("\tlogical function xcheck(x)\n\
	double precision x(%d)\n\
	double precision oldx(%d)\n", nv1, nv1);
		if ((i = bmax.nvt) > 0)
			printf("\tdouble precision v(%d)\n", i);
		if (comb > 0)
			plcommon(1);
		else
			printf(
		  "\tcommon /xkindc/ xkind\n\tinteger xkind\n\tsave /xkindc/\n");
		printf("\tif (xkind .lt. 0) then\n\t\ti = 1\n\t\tgoto 20\n\t\tendif\n\
	do 10 i = 1, %d\n\t\tif (x(i) .ne. oldx(i)) goto 20\n 10\t\tcontinue\n\
	xcheck = .false.\n\treturn\n 20 \tdo 30 i = i, %d\n\
 30\t\toldx(i) = x(i)\n\txkind = 0\n", nv1, nv1);
		for(i = 0, c = cexps; i < comb; c++, i++)
			com_out(c->e, c->L, c->nlin, i);
		printf("\txcheck = .true.\n\tend\n");
		}

	for(i = 1; i < 3; i++) {
		if (owant & i)
			fobj_output(i-1);
		if (cwant & i)
			fcon_output(i-1);
		}
	}

 static void
#ifdef KR_headers
usage(rc)
	int rc;
#else
usage(int rc)
#endif
{
	static char *op[] = {
		"-1 { produce feval0, ceval0 (function values only) }",
		"-2 { produce feval, ceval (functions and gradients, default) }",
		"-3 { produce feval0, ceval0, feval, ceval }",
		"-4 { produce feval, ceval0 }",
		"-5 { produce feval0, ceval }",
#ifdef DEBUG
		"-a { write adjoint operations on Stderr (DEBUG option) }",
#endif
		"-d { dense Jacobian }",
		"-f { Fortran output }",
		"-i { no derivatives for integer variables }",
		"-k { K&R C output (default = ANSI C) }",
#ifdef DEBUG
		"-s { show statistics on Stderr (DEBUG option) }",
		"-v { write def var assignments on Stderr (DEBUG option) }",
#endif
		0};
	char **o;

	fprintf(Stderr, "usage: %s [options] file[.nl]\noptions:\n", progname);
	for(o = op; *o; o++)
		fprintf(Stderr, "\t%s\n", *o);
	exit(rc);
	}

 static void
#ifdef KR_headers
cant(s1, s2)
	char *s1;
	char *s2;
#else
cant(char *s1, char *s2)
#endif
{
	fprintf(Stderr, "Can't open %s", s1);
	if (s2)
		fprintf(Stderr, " or %s.nl", s2);
	fputc('\n', Stderr);
	exit(1);
	}

 static void
get_rownos(VOID)
{
	int i = n_con, i1, j, j1;
	cgrad *cg;
	memset((char *)rownos, 0, nzc*sizeof(int));
	while(i1 = i)
		for(cg = Cgrad[--i]; cg; cg = cg->next)
			rownos[cg->goff] = i1;
	if (!Fortran)
		for(i = 0; i <= n_var; i++)
			A_colstarts[i]++;
	if (intsk) {
		i1 = j = 0;
		--intsk;
		for(i = 1; i <= n_var; i++) {
			j1 = A_colstarts[i] - 1;
			if (!intsk[i])
				for(; j < j1; j++)
					rownos[i1++] = rownos[j];
			A_colstarts[i] = i1 + 1;
			j = j1;
			}
		++intsk;
		nzc = i1;
		if (nlvbi && (nlvci < nlvc || nlvoi < nlvo)
		 || nlvci && nlvoi < nlvo) {
			for(i = 0; i < n_var; i++)
				A_colstarts[i] = A_colstarts[i+1];
			i = n_con;
			while(--i >= 0)
			  for(cg = Cgrad[i]; cg; cg = cg->next)
			    if (!intsk[cg->varno])
				cg->goff = --A_colstarts[cg->varno] - 1;
			}
		}
	}

 static void
#ifdef KR_headers
nlvzap(i, j)
	int i;
	int j;
#else
nlvzap(int i, int j)
#endif
{
	memset(intsk + i - j, 1, j);
	}

 int
#ifdef KR_headers
main(argc, argv)
	int argc;
	char **argv;
#else
main(int argc, char **argv)
#endif
{
	FILE *nl;
	char *s, *s0;
	int i, ncom, nv;
	fint L;
	expr_nx *enx;

	ASL_alloc(ASL_read_fg);
	progname = *argv;
	g_fmt_decpt = 1;
 nextarg:
	while ((s = *++argv) && *s == '-')
		for(s0 = s, --argc;;)
		    switch(*++s) {
			default:
				fprintf(Stderr, "invalid option %s\n",s0);
				usage(1);
			case 0:
				goto nextarg;
			case '1':
				want_derivs = 0;
				cwant = owant = derkind = 1;
				break;
			case '2':
			case '3':
				want_derivs = 1;
				cwant = owant = derkind = *s - '0';
				break;
			case '4':
				want_derivs = cwant = 1;
				derkind = owant = 2;
				break;
			case '5':
				want_derivs = owant = 1;
				derkind = cwant = 2;
				break;
#ifdef DEBUG
			case 'a':		/*debug*/
				awrite = 1;	/*debug*/
				break;
			case 's':		/*debug*/
				stats = 1;	/*debug*/
				break;
			case 'v':		/*debug*/
				vwrite = 1;	/*debug*/
				break;
#endif
			case 'd':
				densejac = 1;
				break;
			case 'f':
				Fortran = 1;
				g_fmt_decpt = 2;
				g_fmt_E = 'd';
				Half = ".5d0";
				One = "1d0";
				Zero = "0.d+00";
				Negone = "-1d0";
				opEQ = ".eq.";
				opGE = ".ge.";
				opGT = ".gt.";
				opLE = ".le.";
				opLT = ".lt.";
				opNE = ".ne.";
				opAND = ".and.";
				opOR = ".or.";
				offlfmt1 = "%s";
				offlfmt2 = "%s, %s";
				assign_fmt = "\t%s = %s\n";
				binop_fmt = "\t%s = %s %s %s\n";
				goto_fmt = "\tgo to %d\n";
				ifgo_fmt = "\tif (%s %s %s) go to %d\n";
				label_fmt = "%5d continue\n";
				ifstart_fmt = "\tif (%s %s %s) then\n";
				elseif_fmt = "\telse if (%s %s %s) then\n";
				else_fmt = "\telse\n";
				endif_fmt = "\tendif\n";
				call_fmt = "\t%s(%s)\n";
				case_fmt = label_fmt;
				cond_fmt = "cond(%d)";
				break_fmt = goto_fmt;
				endswitch_fmt = label_fmt;
				zerdiv_fmt = "\tcall zerdiv(%s)\n";
				dv_fmt = "dv(%d)";
				Fortstar = "*";
				pd_fmt = "pd(%d)";
				tv_fmt = "v(%d)";
				x_fmt = "x(%d)";
				grad_fmt = "g(%d)";
				jac_fmt = "J(%d)";
				g0_fmt = "\tg(%d) = %s\n";
				j0_fmt = "\tJ(%d) = %s\n";
				star = "*";
				eos = "\n";
				xcheck = "\twantfg = needfg\n\
	if (xcheck(x) .and. wantfg .eq. 2) wantfg = 3\n";
				xcheck0 = "\txchk = xcheck(x)\n";
				xkind =
		"\tif (mod(xkind,2*%d) .lt. %d) then\n\t\txkind = xkind + %d\n";
				break;
			case 'i':
				skip_int_derivs = 1;
				break;
			case 'k':
				krflag = 2;
				Void = "";
				xcheckdcl = "(x) real *x;";
				break;
			case '?':
				usage(0);
			}
	if (argc > 2 || !s)
		usage(1);
	return_nofile = 1;
	nl = jacdim0(s, L = strlen(s));
	if (!nl) {
		if ((L -= 3) > 0 && !strcmp(s + L, ".nl")) {
			s[L] = 0;
			nl = jacdim0(s, L);
			if (!nl)
				cant(s, s);
			}
		else
			cant(s, 0);
		}
	for(i = 0; i < N_OPS; i++)
		r_ops[i] = (efunc *)i;

	nv1 = c_vars > o_vars ? c_vars : o_vars;
	ncom = (i = comb + comc + como) + comc1 + como1;
	nv2 = nv1 + ncom;

	c_cexp1st = (int *)Malloc((n_con + n_obj + 2)*sizeof(int));
	o_cexp1st = c_cexp1st + n_con + 1;
	zac = (int **)Malloc((n_con + n_obj + i)*sizeof(int*));
	zao = zac + n_con;
	zaC = zao + n_obj;

	if (n_con)
		if (densejac) {
			rownos = (int *)Malloc(nv1*sizeof(int));
			memset((char *)rownos, 0, nv1*sizeof(int));
			}
		else {
			rownos = (int *)Malloc((nzc + nv1 + 1)*sizeof(int));
			A_colstarts = rownos + nzc;
			}

	LUv = bounds = (real *)Malloc((3*nv1+2*n_con)*sizeof(real));
	LUrhs = LUv + 2*nv1;
	X0 = LUrhs + 2*n_con;

	size_expr_n = sizeof(expr_nx);
	fg_read(nl,0);

	needx0check = comb > 0 || derkind & 2;
	c_ifset = (list **)Malloc((n_con + n_obj + ncom0)*sizeof(list *));
	o_ifset = c_ifset + n_con;
	com_ifset = o_ifset + n_obj;
	op_type[OP1POW] = 2;
	op_type[OP2POW] = 11;

	stog = (ograd **)Malloc(nv1*sizeof(ograd *));

	declare[OP_asinh] = "asinh";
	declare[OP_asinh+1] = "asinhd";
	declare[OP_acosh] = "acosh";
	declare[OP_acosh+1] = "acoshd";
	declare[OPPLTERM] = "plterm";

	dvtfree = (int *)Malloc(NDVTGULP*sizeof(int));
	ndvtmax = nvtfree = NDVTGULP;

	Fortran1 = Fortran - 1;

	if (skip_int_derivs) {
		intsk = (char *)Malloc(nv1);
		memset(intsk, 0, nv1);
		if (nlvbi)
			nlvzap(nlvb, nlvbi);
		if (nlvci)
			nlvzap(nlvb+nlvc, nlvci);
		if (nlvoi)
			nlvzap(nlvb+nlvc+nlvo, nlvoi);
		}
	if (n_con && !densejac)
		get_rownos();

	vseen = (int *)Malloc((nv2 + ncom)*sizeof(int));
	cvmap = vseen + nv2;
	npd = 0;
	for(i = 0; i < ncom0; i++)
		cvmap[i] = -(++npd);
	if (ncom1)
		memset((char *)&cvmap[ncom0], 0, ncom1*sizeof(int));

	ndv = ncond = Fortran;

	dvset(&edagread_one, 1);
#ifdef X64_bit_pointers
	{Adjoint *Ap = (Adjoint *)Malloc(amax*sizeof(Adjoint));
	memset((char *)Ap, 0, amax*sizeof(Adjoint));
	for(i = 0; i < amax; i++)
		*(Adjoint **)&adjoints[i] = Ap++;
	}
#else
	memset((char *)adjoints, 0, amax*sizeof(real));
#endif
	rdwalk_one = &adjoints[nv1];
	dwalk_one = Adjp(rdwalk_one);
	if ((i = amax - nv1) > 0) {
		acount = (int *)Malloc(i*sizeof(int));
		memset((char *)acount, 0, i*sizeof(int));
		}

	comwalk(0,comb);
	if (f_b)
		funnelwalk(f_b);
	max_save(&bmax);
	comwalk(comb, combc);
	if (f_c)
		funnelwalk(f_c);
	cde_walk(con_de, n_con, &cmax, c_ifset, c_cexp1st, zac);
	max_restore(&bmax);
	comwalk(combc, ncom0);
	if (f_o)
		funnelwalk(f_o);
	cde_walk(obj_de, n_obj, &omax, o_ifset, o_cexp1st, zao);

	nv = nv1 + ncom;
	for(i = 0; i < nv; i++)
		var_e[i].op = (efunc *)f_OPVARVAL1;

	for(enx = nums; enx; enx = enx->next)
		enx->op = f_OPNUM1;

	if (n_obj || n_con) {
		output_time = 1;
		if (Fortran)
			foutput();
		else
			output();
		}
#ifdef DEBUG
	if (stats) {
		*stub_end = 0;
		fprintf(Stderr, "%s: derps_in %d derps_out %d nderps %d\n",
			filename, derps_in, derps_out, nderps);
		}
#endif
	return 0;
	}

 char *
#ifdef KR_headers
e_val(e, buf)
	expr *e;
	char *buf;
#else
e_val(expr *e, char *buf)
#endif
{
	int i;
	if ((i = e->a) >= 0)
		sprintf(buf, tv_fmt, Fortran1 + i);
	else
		sprintf(buf, pd_fmt, Fortran1 - i);
	return buf;
	}
