From 23f47c1b08efe756bc93f19a034f14fa5a0593ab Mon Sep 17 00:00:00 2001 From: Christian Valente Date: Tue, 14 Nov 2023 19:12:42 +0200 Subject: [PATCH] 20231111 solvers.tgz: minor tweaks that should not affect results. solvers2.tgz: minor tweaks as in solvers.tgz, plus a bug fix to pfghread.c that affected Hessians of imported functions having several arguments involving variables. --- CMakeLists.txt | 1 + src/solvers/asldate.c | 2 +- src/solvers/changes | 14 +- src/solvers/dtoa.c | 8 +- src/solvers/fgh_read.c | 2 +- src/solvers/pfg_read.c | 14 +- src/solvers/sphes.c | 1 - src/solvers/xsum0.out | 8 +- src/solvers2/asldate.c | 2 +- src/solvers2/dtoa.c | 8 +- src/solvers2/fgh_read.c | 1726 ++++++++++++++ src/solvers2/pfg_read.c | 5007 +++++++++++++++++++++++++++++++++++++++ src/solvers2/pfghread.c | 2 +- src/solvers2/sphes.c | 2 - 14 files changed, 6772 insertions(+), 25 deletions(-) create mode 100644 src/solvers2/fgh_read.c create mode 100644 src/solvers2/pfg_read.c diff --git a/CMakeLists.txt b/CMakeLists.txt index 1d14fd6..d4235aa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -247,6 +247,7 @@ add_prefix(ASL2_SOURCES ${ASL2_SOURCE_DIR}/ degree.c derprop.c dtoa1.c duthes.c dynlink.c eval1.c eval2.c ewalloc1.c ewalloc2.c f_read.c fg_read.c fg_write.c + fgh_read.c pfg_read.c fpecatch.c fpinit.c fullhes.c func_add.c funcadd1.c g_fmt.c genrowno.c getenv.c getstub.c htcl.c indic_cons.c jac0dim.c diff --git a/src/solvers/asldate.c b/src/solvers/asldate.c index 64065c5..5f6da21 100644 --- a/src/solvers/asldate.c +++ b/src/solvers/asldate.c @@ -1 +1 @@ -long ASLdate_ASL = 20230330; +long ASLdate_ASL = 20231111; diff --git a/src/solvers/changes b/src/solvers/changes index 69efbaa..f107bc6 100644 --- a/src/solvers/changes +++ b/src/solvers/changes @@ -3334,8 +3334,18 @@ piecewise-linear terms. pshvprod.c in solvers2.tgz: fix a bug in computing Hessians involving piecewise-linear terms. -20230329 +20230330 jac0dim.c: change return_nofile (= asl->i.return_nofile_) so -return_nofile = 0 or 1 work as before and (return_nofile & 2) causes +(return_nofile & 1) works as before and (return_nofile & 2) causes a quiet return of 0 when the given file does not have a valid .nl header (the first 10 lines). The default return_nofile is still 0. + +20230818 + solvers2.tgz: fix an error-message glitch in fg_write.c; omit +source file jac2dim.c. + +20231111 + solvers.tgz: minor tweaks that should not affect results. + solvers2.tgz: minor tweaks as in solvers.tgz, plus a bug fix to +pfghread.c that affected Hessians of imported functions having +several arguments involving variables. diff --git a/src/solvers/dtoa.c b/src/solvers/dtoa.c index 076d39c..20089ea 100644 --- a/src/solvers/dtoa.c +++ b/src/solvers/dtoa.c @@ -1869,7 +1869,7 @@ mult(Bigint *a, Bigint *b MTd) #else #ifdef Pack_32 for(; xb < xbe; xb++, xc0++) { - if (y = *xb & 0xffff) { + if ((y = *xb & 0xffff)) { x = xa; xc = xc0; carry = 0; @@ -1883,7 +1883,7 @@ mult(Bigint *a, Bigint *b MTd) while(x < xae); *xc = carry; } - if (y = *xb >> 16) { + if ((y = *xb >> 16)) { x = xa; xc = xc0; carry = 0; @@ -5235,9 +5235,11 @@ dtoa_r(double dd, int mode, int ndigits, int *decpt, int *sign, char **rve, char #ifndef SET_INEXACT #ifdef Check_FLT_ROUNDS try_quick = Rounding == 1; +#else + try_quick = 1; #endif #endif /*SET_INEXACT*/ -#endif +#endif /*USE_BF96*/ if (mode > 5) { mode -= 4; diff --git a/src/solvers/fgh_read.c b/src/solvers/fgh_read.c index 7ba9a75..d7d82af 100644 --- a/src/solvers/fgh_read.c +++ b/src/solvers/fgh_read.c @@ -98,7 +98,7 @@ sorry_CLP(EdRead *R, const char *what) #define memadj(x) (((x) + (sizeof(long)-1)) & ~(sizeof(long)-1)) #endif - extern char* f_OPHOL(expr* e A_ASL); + extern char* f_OPHOL; extern efunc f_OPPLTERM, f_OPVARVAL, f_OPFUNCALL; extern sfunc f_OPIFSYM; diff --git a/src/solvers/pfg_read.c b/src/solvers/pfg_read.c index 356157e..4624c86 100644 --- a/src/solvers/pfg_read.c +++ b/src/solvers/pfg_read.c @@ -1950,13 +1950,13 @@ linterms(Static *S, cexp *c, real scale) static void dvwalk(Static *S, int i) { - int n; - ograd *og, *og0; + cexp *c; + dv_info *dvi; + expr *e; + int n, split; linarg *tl, **tp; + ograd *og, *og0; ps_func f; - dv_info *dvi; - cexp *c; - int split; ASLTYPE *asl = S->asl; Termno++; @@ -1981,8 +1981,10 @@ dvwalk(Static *S, int i) tlistgen(S, &f); del_Elemtemp(S, last_psb_elem); } + else if ((e = c->e)) + og = awalk(S, e); else - og = awalk(S, c->e); + og = 0; og0 = og; if (c->nlin) og = af_sum(S, og, linterms(S,c,1.)); diff --git a/src/solvers/sphes.c b/src/solvers/sphes.c index a50065f..3cd3747 100644 --- a/src/solvers/sphes.c +++ b/src/solvers/sphes.c @@ -1219,7 +1219,6 @@ sphes_ASL(ASL *a, SputInfo **pspi, real *H, int nobj, real *ow, real *y) *rtodo++ = 0; /* reset */ while((uhw = uhwi)) { uhwi = uhwi->next; - si = s; ogp = uhw->ogp; r = uhw->r; ogpe = ogp + r->n; diff --git a/src/solvers/xsum0.out b/src/solvers/xsum0.out index 789d76c..9d3d915 100644 --- a/src/solvers/xsum0.out +++ b/src/solvers/xsum0.out @@ -9,7 +9,7 @@ arithchk.c 63b0185 6532 asl.h 1fe3527a 42925 asl_pfg.h 6483336 1268 asl_pfgh.h 49b5a53 1288 -asldate.c 158d058e 29 +asldate.c f734fb13 29 atof.c 1fabc7d3 1747 auxinfo.c 1a3c8690 1488 avltree.c 10aeaa58 11346 @@ -32,7 +32,7 @@ conval.c f1e4dc7d 4569 degree.c 3eafa26 4337 derprop.c e75c750f 1361 details.c0 fa7ec7b3 182 -dtoa.c fd6feb77 158408 +dtoa.c 576045b 158447 dtoa1.c 1928cdca 3352 duthes.c e13b176c 3874 dvalue.hd 9606f1 1125 @@ -103,7 +103,7 @@ op_type.hd f803bc0b 1287 op_typeb.hd 1c1e0c07 1287 opcode.hd 12bc8de3 1348 opnos.hd 1a7c2f8d 78 -pfg_read.c 10bf2455 98287 +pfg_read.c 1dca51d4 98321 pfghread.c b815512 1107 printf.c 147de217 24293 pshvprod.c 1dcadea7 32089 @@ -125,7 +125,7 @@ rops2.c 25c9f1 20459 sigcatch.c fe35ba8b 2374 sjac0dim.c ff686adb 1546 sos_add.c fb5da0e8 23060 -sphes.c e4e0bfbe 27806 +sphes.c f70a4208 27795 sprintf.c e6ce1d45 3050 sscanf.c f5a6f021 3327 stderr.c fc99b4f4 2096 diff --git a/src/solvers2/asldate.c b/src/solvers2/asldate.c index 5d67b5d..5f6da21 100644 --- a/src/solvers2/asldate.c +++ b/src/solvers2/asldate.c @@ -1 +1 @@ -long ASLdate_ASL = 20230818; +long ASLdate_ASL = 20231111; diff --git a/src/solvers2/dtoa.c b/src/solvers2/dtoa.c index 076d39c..20089ea 100644 --- a/src/solvers2/dtoa.c +++ b/src/solvers2/dtoa.c @@ -1869,7 +1869,7 @@ mult(Bigint *a, Bigint *b MTd) #else #ifdef Pack_32 for(; xb < xbe; xb++, xc0++) { - if (y = *xb & 0xffff) { + if ((y = *xb & 0xffff)) { x = xa; xc = xc0; carry = 0; @@ -1883,7 +1883,7 @@ mult(Bigint *a, Bigint *b MTd) while(x < xae); *xc = carry; } - if (y = *xb >> 16) { + if ((y = *xb >> 16)) { x = xa; xc = xc0; carry = 0; @@ -5235,9 +5235,11 @@ dtoa_r(double dd, int mode, int ndigits, int *decpt, int *sign, char **rve, char #ifndef SET_INEXACT #ifdef Check_FLT_ROUNDS try_quick = Rounding == 1; +#else + try_quick = 1; #endif #endif /*SET_INEXACT*/ -#endif +#endif /*USE_BF96*/ if (mode > 5) { mode -= 4; diff --git a/src/solvers2/fgh_read.c b/src/solvers2/fgh_read.c new file mode 100644 index 0000000..d7d82af --- /dev/null +++ b/src/solvers2/fgh_read.c @@ -0,0 +1,1726 @@ +/******************************************************************* +Copyright (C) 2016 AMPL Optimization, Inc.; written by David M. Gay. + +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. + +The author and AMPL Optimization, Inc. disclaim all warranties with +regard to this software, including all implied warranties of +merchantability and fitness. In no event shall the author 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 "jac2dim.h" +#include "opnos.hd" + +#ifdef __cplusplus +extern "C" { +#endif + +extern real con2ival_nomap_ASL(ASL*, int, real*, fint*); +extern void con2grd_nomap_ASL(ASL*, int, real *, real*, fint*); + +#define Egulp 400 + + static int com11, n_com1, ncom_togo, nvar0, nvinc, nvref; + static int *vrefnext, *vrefx; + static expr_if *if2list_end; + static expr_va *varg2list_end; + +#define NFHASH 23 + + static derp *last_d; + static expr *last_e; + +#define efunc efunc2 + + static relo *relolist, *relo2list; + static expr_if *iflist, *if2list; + static expr_va *varglist, *varg2list; + static real one = 1.; + static int nderp; +#undef nzc + static int amax1, imap_len, k_seen, lasta, lasta0, lasta00, lastj, + max_var, nocopy, nv0, nv01, nv0b, nv0c, nv1, nzc, nzclim; + static int co_first = 1; + static int *imap, *zc, *zci; + + static expr *(*holread) ANSI((EdRead*)); + extern real f_OPCPOW ANSI((expr* A_ASL)); + + static ASL_fgh *asl; + + static void +ed_reset(ASL *a) +{ + + a->i.memLast = a->i.memNext = 0; + k_seen = 0; + nvref = 0; + vrefnext = vrefx = 0; + relolist = relo2list = 0; + last_d = 0; + iflist = if2list = if2list_end = 0; + varglist = varg2list = varg2list_end = 0; + imap = 0; + amax1 = imap_len = lastj = nocopy = 0; + com11 = lasta = lasta0 = lasta00 = n_com1 = nderp = 0; + last_e = 0; + co_first = 1; + } + + static void +fscream(EdRead *R, const char *name, int nargs, const char *kind) +{ + scream(R, ASL_readerr_argerr, + "line %ld: attempt to call %s with %d %sargs\n", + R->Line, name, nargs, kind); + } + + static void +sorry_CLP(EdRead *R, const char *what) +{ + fprintf(Stderr, + "Sorry, %s cannot handle %s.\n", + progname ? progname : "", what); + exit_ASL(R,ASL_readerr_CLP); + } + +#ifdef Double_Align +#define memadj(x) x +#else +#define memadj(x) (((x) + (sizeof(long)-1)) & ~(sizeof(long)-1)) +#endif + + extern char* f_OPHOL; + extern efunc f_OPPLTERM, f_OPVARVAL, f_OPFUNCALL; + extern sfunc f_OPIFSYM; + + static void +new_derp(int a, int b, real *c) +{ + derp *d; + if (a == nv1) + return; + nderp++; + d = (derp *)mem(sizeof(derp)); + d->next = last_d; + last_d = d; + d->a.i = a; + d->b.i = b; + d->c.rp = c; + } + + static derp * +new_relo(expr *e, derp *Dnext, int *ap) +{ + relo *r; + derp *d; + + if (last_d != Dnext) { + *ap = e->a; + for(d = last_d; d->next != Dnext; d = d->next); + d->next = 0; + } + else { + last_d = 0; + new_derp(e->a, *ap = lasta++, &one); + } + if (!last_d) + return 0; + r = (relo *)mem(sizeof(relo)); + r->next = relolist; + r->next2 = relo2list; + relo2list = relolist = r; + r->D = r->Dcond = last_d; + r->Dnext = Dnext; + return r->D; + } + + static relo * +new_relo1(derp *Dnext) +{ + relo *r; + + r = (relo *)mem(sizeof(relo)); + r->next = relolist; + relolist = r; + r->D = 0; + r->Dnext = Dnext; + return r; + } + + static expr * +new_expr(int opcode, expr *L, expr *R, int deriv) +{ + expr *rv; + extern efunc f_OP1POW, f_OP2POW, f_OPPOW; + efunc *o; + int L1, R1, i; + + o = r_ops[opcode]; + if (o == f_OPPOW) { + if (R->op == f_OPNUM) + if (((expr_n *)R)->v == 2.) { + o = f_OP2POW; + R = 0; + } + else + o = f_OP1POW; + else if (L->op == f_OPNUM) + o = f_OPCPOW; + } + rv = (expr *)mem(sizeof(expr)); + rv->op = o; + rv->L.e = L; + rv->R.e = R; + if (deriv) { + L1 = L && L->op != f_OPNUM; + R1 = R && R->op != f_OPNUM; + if (L1 | R1) { + rv->a = lasta++; + if (L1) + new_derp(L->a, rv->a, &rv->dL); + if (R1) + new_derp(R->a, rv->a, &rv->dR); + rv->bak = last_e; + last_e = rv; + if (R) + rv->dL2 = rv->dLR = rv->dR2 = 0; + else if (o == f_OP2POW) + rv->dL2 = 2; + else + rv->dL2 = 0; + if (L1) + if (R1) + switch(opcode) { + case PLUS: i = Hv_plusLR; break; + case MINUS: i = Hv_minusLR; break; + case MULT: i = Hv_timesLR; break; + case DIV: i = Hv_divLR; break; + default: i = Hv_binaryLR; + } + else switch(opcode) { + case PLUS: + case MINUS: i = Hv_plusL; break; + case MULT: i = Hv_timesL; break; + case UMINUS: i = Hv_negate; break; + default: i = Hv_unary; + } + else + switch(opcode) { + case PLUS: i = Hv_plusR; break; + case MINUS: i = Hv_minusR; break; + case MULT: i = Hv_timesR; break; + default: i = Hv_binaryR; + } + rv->dO.i = i; + } + } + return rv; + } + + static char op_type[] = { +#include "op_type.hd" + }; + + static expr * +eread(EdRead *R, int deriv) +{ + arglist *al; + argpair *ap, *da, *sap; + char **sa; + char *dig; + de *d; + derp *dsave; + efunc *op; + expr *L, *arg, **args, **args1, **argse, *rv; + expr_f *rvf; + expr_if *rvif; + expr_n *rvn; + expr_va *rva; + func_info *fi; + int *at, *nn, *nn0; + int a0, a1, i, i1, j, j1, k, kd, kd2, ks, numargs, symargs; + int (*Xscanf)(EdRead*, const char*, ...); + long L1; + plterm *p; + real *b, **fh, *hes, r, *ra; + unsigned int Ls; + static real dvalue[] = { +#include "dvalue.hd" + }; + + Xscanf = xscanf; + switch(edag_peek(R)) { + case 'f': + if (Xscanf(R, "%d %d", &i, &j) != 2) + badline(R); + fi = funcs[i]; + if (fi->nargs >= 0) { + if (fi->nargs != j) { + bad_nargs: + fscream(R, fi->name, j, ""); + } + } + else if (-(1+fi->nargs) > j) + goto bad_nargs; + rvf = (expr_f *)mem(sizeof(expr_f) + + (j-1)*sizeof(expr *)); + rvf->op = f_OPFUNCALL; + rvf->fi = fi; + rvf->dO.i = Hv_func; + args = rvf->args; + argse = args + j; + k = ks = symargs = numargs = 0; + while(args < argse) { + arg = *args++ = eread(R, deriv); + if ((op = arg->op) == (efunc*)f_OPHOL) + symargs++; + else if (op == (efunc*)f_OPIFSYM) + ks++; + else { + numargs++; + if (op != f_OPNUM) + k++; + } + } + symargs += ks; + if (symargs && !(fi->ftype & 1)) + fscream(R, fi->name, symargs, "symbolic "); + if (deriv) { + kd2 = (kd = k) ? numargs*(numargs+1) >> 1 : 0; + rvf->a = lasta++; + } + else { + kd = kd2 = 0; + rvf->a = nv1; + } + Ls = sizeof(arglist) + + (k + kd + ks)*sizeof(argpair) + + kd*kd*sizeof(real *) + + (numargs+kd+kd2)*sizeof(real) + + symargs*sizeof(char *) + + j*sizeof(int); + if (kd) + Ls += numargs*sizeof(real); + ra = (real *)mem(Ls); + dig = 0; + if (kd < numargs && kd) + dig = (char*)mem(numargs); + b = kd ? ra + numargs : ra; + hes = b + numargs; + al = rvf->al = (arglist *)(hes + kd2); + al->n = numargs + symargs; + al->nr = numargs; + al->ra = ra; + if (kd) { + al->derivs = b; + al->hes = hes; + memset(b, 0, (numargs + kd2)*sizeof(real)); + } + else + al->derivs = al->hes = 0; + al->dig = dig; + al->funcinfo = fi->funcinfo; + al->AE = asl->i.ae; + al->sa = (Const char**)(sa = (char **)(al + 1)); + ap = rvf->ap = (argpair *)(sa + symargs); + rvf->da = da = ap + k; + sap = rvf->sap = da + kd; + rvf->fh = fh = (real **)(sap + ks); + at = al->at = (int *)(fh + kd*kd); + symargs = numargs = i = 0; + nn = nn0 = (int *)b; + for(args = rvf->args; args < argse; at++) { + arg = *args++; + if ((op = arg->op) == (efunc*)f_OPHOL) { + *at = --symargs; + *sa++ = ((expr_h *)arg)->sym; + } + else if (op == (efunc*)f_OPIFSYM) { + *at = --symargs; + sap->e = arg; + (sap++)->u.s = sa++; + } + else { + *at = numargs++; + j = 1; + if (op == f_OPNUM) + *ra = ((expr_n *)arg)->v; + else { + ap->e = arg; + (ap++)->u.v = ra; + if (kd) { + j = 0; + new_derp(arg->a, + rvf->a, b); + *b = 0; + da->e = arg; + (da++)->u.v = b; + *nn++ = i; + } + } + if (dig) + *dig++ = j; + b++; + ra++; + i++; + } + } + rvf->ape = ap; + rvf->sape = sap; + rvf->dae = da; + if (kd) { + rvf->bak = last_e; + last_e = (expr *)rvf; + rvf->dO.i = Hv_func; + for(i1 = 0; i1 < kd; i1++) { + i = nn0[i1]; + for(j1 = 0; j1 < kd; j1++) { + j = nn0[j1]; + *fh++ = &hes[i >= j ? (i*(i+1)>>1)+j + : (j*(j+1)>>1)+i]; + } + } + } + return (expr *)rvf; + + case 'h': + return holread(R); + + case 's': + if (Xscanf(R, "%hd", &L1) != 1) + badline(R); + r = L1; + goto have_r; + + case 'l': + if (Xscanf(R, "%ld", &L1) != 1) + badline(R); + r = L1; + goto have_r; + + case 'n': + if (Xscanf(R, "%lf", &r) != 1) + badline(R); + have_r: + rvn = (expr_n *)mem(sizeof(expr_n)); + rvn->op = f_OPNUM_ASL; + rvn->v = r; + return (expr *)rvn; + + case 'o': + break; + + case 'v': + if (Xscanf(R,"%d",&k) != 1 || k < 0) + badline(R); + if (k >= nvar0) + k += nvinc; + if (k > max_var) + badline(R); + if (k < nv01 && deriv && !zc[k]++) + zci[nzc++] = k; + return (expr *)(var_e + k); + + default: + badline(R); + } + + if (Xscanf(R, asl->i.opfmt, &k) != 1 || k < 0 || k >= sizeof(op_type)) + badline(R); + switch(op_type[k]) { + + case 1: /* unary */ + rv = new_expr(k, eread(R, deriv), 0, deriv); + rv->dL = dvalue[k]; /* for UMINUS, FLOOR, CEIL */ + return rv; + + case 2: /* binary */ + if (dvalue[k] == 11) + deriv = 0; + L = eread(R, deriv); + rv = new_expr(k, L, eread(R, deriv), deriv); + rv->dL = 1.; + rv->dR = dvalue[k]; /* for PLUS, MINUS, REM */ + return rv; + + case 3: /* vararg (min, max) */ + i = -1; + Xscanf(R, "%d", &i); + if (i <= 0) + badline(R); + rva = (expr_va *)mem(sizeof(expr_va)); + rva->op = r_ops[k]; + rva->L.d = d = + (de *)mem(i*sizeof(de) + sizeof(expr *)); + rva->next = varglist; + varglist = varg2list = rva; + if (!last_d) { + new_derp(lasta, lasta, &one); + lasta++; + } + rva->d0 = dsave = last_d; + rva->bak = last_e; + a0 = a1 = lasta; + for(j = 0; i > 0; i--, d++) { + last_d = dsave; + last_e = 0; + d->e = L = eread(R, deriv); + d->ee = last_e; + if (L->op == f_OPNUM || L->a == nv1) { + d->d = dsave; + d->dv.i = nv1; + } + else { + if (deriv) + d->d = new_relo(L, dsave, + &d->dv.i); + if (a1 < lasta) + a1 = lasta; + lasta = a0; + } + } + d->e = 0; /* sentinnel expr * */ + rva->a = lasta = a1; + last_d = dsave; + if (deriv) { + new_derp(0, lasta++, &one); + /* f_MINLIST or f_MAXLIST will replace the 0 */ + rva->R.D = last_d; + nocopy = 1; + last_e = (expr *)rva; + rva->dO.i = Hv_vararg; + } + else { + rva->R.D = 0; + last_e = rva->bak; + } + return (expr *)rva; + + case 4: /* piece-wise linear */ + i = -1; + Xscanf(R, "%d", &i); + if (i <= 1) + badline(R); + j = 2*i - 1; + p = (plterm *)mem(sizeof(plterm) + (j-1)*sizeof(real)); + p->n = i; + b = p->bs; + do { + switch(edag_peek(R)) { + case 's': + if (Xscanf(R,"%hd",&L1) != 1) + badline(R); + r = L1; + break; + case 'l': + if (Xscanf(R,"%ld",&L1) != 1) + badline(R); + r = L1; + break; + case 'n': + if (Xscanf(R,"%lf",&r) == 1) + break; + default: + badline(R); + } + *b++ = r; + } + while(--j > 0); + if (b[-2] <= 0.) + p->z = 2*i - 2; + else { + b = p->bs + 1; + while(*b <= 0.) + b += 2; + p->z = (b - p->bs) - 1; + } + if (edag_peek(R) != 'v' + || Xscanf(R, "%d", &k) != 1 + || k < 0 || k >= max_var) + badline(R); + if (k >= nvar0) + k += nvinc; + rv = (expr *)mem(sizeof(expr)); + rv->op = f_OPPLTERM; + rv->L.p = p; + rv->R.e = (expr *)(var_e + k); + if (deriv) { + new_derp(k, rv->a = lasta++, &rv->dL); + rv->bak = last_e; + last_e = rv; + rv->dO.i = Hv_plterm; + } + return rv; + + case 5: /* if */ + rvif = (expr_if *)mem(sizeof(expr_if)); + rvif->op = r_ops[k]; + rvif->next = iflist; + iflist = if2list = rvif; + if (!last_d) { + new_derp(lasta, lasta, &one); + lasta++; + } + rvif->d0 = dsave = last_d; + rvif->bak = last_e; + rvif->e = eread(R, 0); + last_e = 0; + a0 = lasta; + rvif->T = L = eread(R, deriv); + rvif->Te = last_e; + j = 0; + if (L->op == f_OPNUM) { + rvif->dT = dsave; + rvif->Tv.i = nv1; + } + else if ((j = deriv)) + rvif->dT = new_relo(L, dsave, &rvif->Tv.i); + a1 = lasta; + lasta = a0; + last_d = dsave; + last_e = 0; + rvif->F = L = eread(R, deriv); + rvif->Fe = last_e; + if (L->op == f_OPNUM) { + rvif->dF = dsave; + rvif->Fv.i = nv1; + } + else if (j) + rvif->dF = new_relo(L, dsave, &rvif->Fv.i); + if (lasta < a1) + lasta = a1; + last_d = dsave; + if (j) { + new_derp(0, rvif->a = lasta++, &one); + rvif->D = last_d; + nocopy = 1; + last_e = (expr *)rvif; + rvif->dO.i = Hv_if; + } + else { + rvif->a = nv1; + rvif->D = 0; + last_e = rvif->bak; + } + return (expr *)rvif; + + case 11: /* OPCOUNT */ + deriv = 0; + /* no break */ + case 6: /* sumlist */ + i = 0; + Xscanf(R, "%d", &i); + if (i <= 2 && (op_type[k] == 6 || i < 1)) + badline(R); + rv = (expr *)mem(sizeof(expr) - sizeof(real) + + (deriv ? i+i+1 : i)*sizeof(expr *)); + rv->op = r_ops[k]; + rv->a = deriv ? lasta++ : nv1; + rv->L.ep = args = (expr **)&rv->dR; + if (deriv) { + j = 0; + do { + *args++ = L = eread(R, deriv); + if (L->op != f_OPNUM) { + new_derp(L->a, rv->a, &one); + j++; + } + } + while(--i > 0); + if (j) { + rv->bak = last_e; + last_e = rv; + rv->dO.i = Hv_sumlist; + argse = args; + args1 = rv->L.ep; + while(args1 < args) { + arg = *args1++; + if (arg->op != f_OPNUM) + *argse++ = arg; + } + *argse = 0; + } + } + else do + *args++ = eread(R, deriv); + while(--i > 0); + rv->R.ep = args; + return rv; + } + badline(R); + return 0; + } + + static list * +new_list(list *nxt) +{ + list *rv = (list *)mem(sizeof(list)); + rv->next = nxt; + return rv; + } + + static list * +crefs(void) +{ + int i; + list *rv = 0; + + while(nzc > 0) { + if ((i = zci[--nzc]) >= nv0) { + rv = new_list(rv); + rv->item.i = i; + } + zc[i] = 0; + } + return rv; + } + + static funnel * +funnelfix(funnel *f) +{ + cexp *ce; + funnel *fnext, *fprev; + + for(fprev = 0; f; f = fnext) { + fnext = f->next; + f->next = fprev; + fprev = f; + ce = f->ce; + ce->z.i = ce->d->b.i; + } + return fprev; + } + + static derp * +derpadjust(derp *d0, int a, derp *dnext) +{ + derp *d, *d1; + int *r, *re; + relo *rl; + expr_if *il, *ile; + expr_va *vl, *vle; + de *de1; + + if (!(d = d0)) + return dnext; + r = imap + lasta0; + re = imap + lasta; + while(r < re) + *r++ = a++; + if (amax < a) + amax = a; + r = imap; + for(;; d = d1) { + d->a.i = r[d->a.i]; + d->b.i = r[d->b.i]; + if (!(d1 = d->next)) + break; + } + d->next = dnext; + if ((rl = relo2list)) { + relo2list = 0; + do { + d = rl->Dcond; + do { + d->a.i = r[d->a.i]; + d->b.i = r[d->b.i]; + } + while((d = d->next)); + } + while((rl = rl->next2)); + } + if (if2list != if2list_end) { + ile = if2list_end; + if2list_end = il = if2list; + do { + il->Tv.i = r[il->Tv.i]; + il->Fv.i = r[il->Fv.i]; + } + while((il = il->next) != ile); + } + if (varg2list != varg2list_end) { + vle = varg2list_end; + varg2list_end = vl = varg2list; + do { + for(de1 = vl->L.d; de1->e; de1++) + de1->dv.i = r[de1->dv.i]; + } + while((vl = vl->next) != vle); + } + return d0; + } + + static derp * +derpcopy(cexp *ce, derp *dnext) +{ + derp *d, *dprev; + int *map; + derp d00; + + if (!(d = ce->d)) + return dnext; + map = imap; + for(dprev = &d00; d; d = d->next) { + new_derp(map[d->a.i], map[d->b.i], d->c.rp); + dprev = dprev->next = last_d; + } + dprev->next = dnext; + return d00.next; + } + + static void +imap_alloc(void) +{ + int i, *r, *re; + + if (imap) { + imap_len += lasta; + imap = (int *)Realloc(imap, imap_len*Sizeof(int)); + return; + } + imap_len = amax1 > lasta ? amax1 : lasta; + imap_len += 100; + r = imap = (int *)Malloc(imap_len*Sizeof(int)); + for(i = 0, re = r + nv1+1; r < re;) + *r++ = i++; + } + + static int +compar(const void *a, const void *b, void *v) +{ + Not_Used(v); + return *(int*)a - *(int*)b; + } + + static void +comsubs(int alen, cde *d) +{ + list *L; + int a, i, j, k; + int *r, *re; + cexp *ce; + derp *D, *dnext; + relo *R; + + D = last_d; + a = lasta00; + dnext = 0; + R = 0; + for(i = j = 0; i < nzc; i++) + if ((k = zci[i]) >= nv0) + zci[j++] = k; + else + zc[k] = 0; + if ((nzc = j)) { + for(i = 0; i < nzc; i++) + for(L = cexps[zci[i]-nv0].cref; L; L = L->next) + if (!zc[L->item.i]++) + zci[nzc++] = L->item.i; + if (nzc > 1) { + if (nzc < nzclim) + qsortv(zci, nzc, sizeof(int), compar, NULL); + else for(i = nv0, j = 0; i < max_var; i++) + if (zc[i]) + zci[j++] = i; + } + if (nzc > 0) { + R = new_relo1(dnext); + i = 0; + do { + j = zci[i]; + zc[j] = 0; + ce = &cexps[j - nv0]; + if (ce->funneled) + imap[var_e[j].a] = a++; + else { + r = imap + ce->z.i; + re = r + ce->zlen; + while(r < re) + *r++ = a++; + } + dnext = R->D = derpcopy(ce, R->D); + } + while(++i < nzc); + nzc = 0; + } + } + if (D || R) { + if (!R) + R = new_relo1(dnext); + D = R->D = derpadjust(D, a, R->D); + if (d->e->op != f_OPVARVAL) + d->e->a = imap[d->e->a]; + } + d->d = D; + a += alen; + d->zaplen = (a > lasta00 ? a - nv1 : 0)*sizeof(real); + if (amax < a) + amax = a; + } + + static void +co_read(EdRead *R, cde *d, int wd) +{ + int alen; + + d->com11 = com11; + d->n_com1 = n_com1; + com11 += n_com1; + n_com1 = 0; + + if (amax1 < lasta) + amax1 = lasta; + if (co_first) { + co_first = 0; + if (imap_len < lasta) + imap_alloc(); + f_b = funnelfix(f_b); + f_c = funnelfix(f_c); + f_o = funnelfix(f_o); + } + if (!lastj) { + lasta = lasta0; + last_d = 0; + } + lastj = 0; + last_e = 0; + d->e = eread(R, wd); + d->ee = last_e; + alen = lasta - lasta0; + if (imap_len < lasta) + imap_alloc(); + comsubs(alen, d); + } + + static linpart * +linpt_read(EdRead *R, int nlin) +{ + int (*Xscanf)(EdRead*, const char*, ...); + linpart *L, *rv; + + if (nlin <= 0) + return 0; + Xscanf = xscanf; + L = rv = (linpart *)mem(nlin*sizeof(linpart)); + do { + if (Xscanf(R, "%d %lf", &L->v.i, &L->fac) != 2) + badline(R); + L++; + } + while(--nlin > 0); + return rv; + } + + static int +funnelkind(cexp *ce, int *ip) +{ + int i, j, k, nzc0, rv; + int *vr, *vre; + + ce->vref = 0; + if (!(nzc0 = nzc)) + return 0; + for(i = k = rv = 0; i < nzc; i++) + if ((j = zci[i]) < nv0) { + if (k >= maxfwd) + goto done; + vrefx[k++] = j; + } + else { + if (!(vr = cexps[j-nv0].vref)) + goto done; + vre = vr + *vr; + while(++vr <= vre) + if (!zc[*vr]++) + zci[nzc++] = *vr; + } + if (k >= nvref) { + nvref = (maxfwd + 1)*(ncom_togo < vrefGulp + ? ncom_togo : vrefGulp); + vrefnext = (int *)M1alloc(nvref*Sizeof(int)); + } + vr = ce->vref = vrefnext; + *vr = k; + vrefnext += i = k + 1; + nvref -= i; + for(i = 0; i < k; i++) + *++vr = vrefx[i]; + if (nderp > 3*k && !nocopy) { + *ip = k; + return 2; + } + else { + done: + if (nocopy || nderp > 3*nzc0) + rv = 1; + } + while(nzc > nzc0) + zc[zci[--nzc]] = 0; + return rv; + } + + static void +cexp_read(EdRead *R, int k, int nlin) +{ + int a, fk, i, j, la0, na; + funnel *f, **fp; + linpart *L, *Le; + expr *e; + cplist *cl, *cl0; + cexp *ce; + + ce = cexps + k - nv0; + L = ce->L = linpt_read(R, ce->nlin = nlin); + nocopy = 0; + last_d = 0; + last_e = 0; + ce->z.i = la0 = lasta; + nderps += nderp; + nderp = 0; + e = ce->e = eread(R, 1); + if (la0 == lasta) { + a = lasta++; + if (e->op != f_OPNUM) + new_derp(e->a, a, &edagread_one); + } + else + a = e->a; + var_e[k].a = a; + ce->zlen = lasta - la0; + for(Le = L + nlin; L < Le; L++) { + new_derp(i = L->v.i, a, &L->fac); + if (!zc[i]++) + zci[nzc++] = i; + } + i = 0; /* only needed to shut up an erroneous warning */ + if ((fk = funnelkind(ce, &i))) { + /* arrange to funnel */ + fp = k < nv0b ? &f_b : k < nv0c ? &f_c : &f_o; + ce->funneled = f = (funnel *)mem(sizeof(funnel)); + f->next = *fp; + *fp = f; + f->ce = ce; + if (imap_len < lasta) + imap_alloc(); + if (fk == 1) { + f->fulld = last_d; + a = lasta00; + for(i = nzc; --i >= 0; ) + if ((j = zci[i]) >= nv0) + imap[var_e[j].a] = a++; + if ((na = ce->zlen) || a > lasta00) + na += a - nv1; + f->fcde.zaplen = na*sizeof(real); + i = nzc; + derpadjust(last_d, a, 0); + } + else { + f->fulld = 0; + f->fcde.e = e; + comsubs(ce->zlen, &f->fcde); + memcpy(zci, vrefx, i*sizeof(int)); + } + last_d = 0; + cl0 = 0; + while(i > 0) + if ((a = var_e[zci[--i]].a) != nv1) { + new_derp(a, lasta0, 0); + cl = (cplist *)mem(sizeof(cplist)); + cl->next = cl0; + cl0 = cl; + cl->ca.i = imap[last_d->a.i]; + cl->cfa = last_d->c.rp = + (real *)mem(sizeof(real)); + } + f->cl = cl0; + var_e[k].a = lasta0++; + lasta = lasta0; + } + lasta0 = lasta; + ce->d = last_d; + ce->ee = last_e; + ce->cref = crefs(); + --ncom_togo; + } + + static void +cexp1_read(EdRead *R, int j, int k, int nlin) +{ + linpart *L, *Le; + cexp1 *ce = cexps1 + (k - nv01); + expr *e; + expr_v *v; + int la0; + + n_com1++; + L = ce->L = linpt_read(R, ce->nlin = nlin); + + if (!lastj) { + last_d = 0; + if (amax1 < lasta) + amax1 = lasta; + lasta = lasta0; + lastj = j; + } + last_e = 0; + la0 = lasta; + e = ce->e = eread(R, 1); + ce->ee = last_e; + v = var_e + k; + if (la0 == lasta) { + j = lasta++; + if (e->op != f_OPNUM) + new_derp(e->a, j, &edagread_one); + } + else + j = e->a; + v->a = j; + v->dO.r = 0; + for(Le = L + nlin; L < Le; L++) + new_derp(L->v.i, j, &L->fac); + } + + static expr * +hvadjust(expr *e) +{ + expr *e0; + + for(e0 = 0; e; e = e->bak) { + e->fwd = e0; + e0 = e; + e->a = e->dO.i; + } + return e0; + } + + static void +co_adjust(cde *d, int n) +{ + cde *de1; + + for(de1 = d + n; d < de1; d++) + d->ef = hvadjust(d->ee); + } + + static void +ifadjust(expr_if *e) +{ + for(; e; e = e->next) { + e->Tv.rp = &adjoints[e->Tv.i]; + e->Fv.rp = &adjoints[e->Fv.i]; + e->Tf = hvadjust(e->Te); + e->Ff = hvadjust(e->Fe); + } + } + + static void +vargadjust(expr_va *e) +{ + de *d; + + for(; e; e = e->next) { + for(d = e->L.d; d->e; d++) { + d->dv.rp = &adjoints[d->dv.i]; + d->ef = hvadjust(d->ee); + } + } + } + + static void +funneladj1(funnel *f) +{ + real *a = adjoints; + derp *d; + cplist *cl; + + for(a = adjoints; f; f = f->next) { + if ((d = f->fulld)) { + f->fcde.d = d; + do { + d->a.rp = &a[d->a.i]; + d->b.rp = &a[d->b.i]; + } + while((d = d->next)); + } + for(cl = f->cl; cl; cl = cl->next) + cl->ca.rp = &a[cl->ca.i]; + } + } + + static void +funneladjust(VOID) +{ + cexp *c, *ce; + linpart *L, *Le; + c = cexps; + for(ce = c + ncom0; c < ce; c++) { + if ((L = c->L)) + for(Le = L + c->nlin; L < Le; L++) + L->v.vp = (Char*)&var_e[L->v.i]; + c->ef = hvadjust(c->ee); + } + + funneladj1(f_b); + funneladj1(f_c); + funneladj1(f_o); + } + + static void +com1adjust(VOID) +{ + cexp1 *c, *ce; + linpart *L, *Le; + + for(c = cexps1, ce = c + ncom1; c < ce; c++) { + for(L = c->L, Le = L + c->nlin; L < Le; L++) + L->v.vp = (Char*)&var_e[L->v.i]; + c->ef = hvadjust(c->ee); + } + } + + static void +adjust_compl_rhs(VOID) +{ + cde *C; + expr *e; + int *Cvar, i, j, nc, stride; + real *L, *U, t, t1; + + L = LUrhs; + if ((U = Urhsx)) + stride = 1; + else { + U = L + 1; + stride = 2; + } + C = con_de; + Cvar = cvar; + nc = n_con; + for(i = nlc; i < nc; i++) + if (Cvar[i] && (e = C[i].e) && e->op == f_OPNUM + && (t = ((expr_n*)e)->v) != 0.) { + t1 = t; + if (L[j = stride*i] > negInfinity) { + L[j] -= t; + t1 = 0.; + } + if (U[j] < Infinity) { + U[j] -= t; + t1 = 0.; + } + ((expr_n*)e)->v = t1; + } + } + + static void +adjust(ASL *a0, int flags) +{ + derp *d, **dp; + real *a = adjoints; + relo *r; + + for(r = relolist; r; r = r->next) { + dp = &r->D; + while((d = *dp)) { + d->a.rp = &a[d->a.i]; + d->b.rp = &a[d->b.i]; + dp = &d->next; + } + *dp = r->Dnext; + } + ifadjust(iflist); + vargadjust(varglist); + if (ncom0) + funneladjust(); + com1adjust(); + co_adjust(con_de, n_con); + co_adjust(obj_de, n_obj); + if (k_seen) { + if (!A_vals) + goff_comp_ASL(a0); + else if (Fortran) + colstart_inc_ASL(a0); + } + if (n_cc > nlcc && nlc < n_con + && !(flags & ASL_no_linear_cc_rhs_adjust)) + adjust_compl_rhs(); + } + + static void +br_read(EdRead *R, int nc, real *L, real *U, int *Cvar, int nv) +{ + ASL *asl; + int i, inc, j, k; + int (*Xscanf)(EdRead*, const char*, ...); + + asl = R->asl; + + if (U) + inc = 1; + else { + U = L + 1; + inc = 2; + } + Xscanf = xscanf; + Xscanf(R, ""); /* purge line */ + for(i = 0; i < nc; i++, L += inc, U += inc) { + switch(edag_peek(R) - '0') { + case 0: + if (Xscanf(R,"%lf %lf",L,U)!= 2) + badline(R); + break; + case 1: + if (Xscanf(R, "%lf", U) != 1) + badline(R); + *L = negInfinity; + break; + case 2: + if (Xscanf(R, "%lf", L) != 1) + badline(R); + *U = Infinity; + break; + case 3: + *L = negInfinity; + *U = Infinity; + Xscanf(R, ""); /* purge line */ + break; + case 4: + if (Xscanf(R, "%lf", L) != 1) + badline(R); + *U = *L; + break; + case 5: + if (Cvar) { + if (Xscanf(R, "%d %d", &k, &j) == 2 + && j > 0 && j <= nv) { + Cvar[i] = j; + *L = k & 2 ? negInfinity : 0.; + *U = k & 1 ? Infinity : 0.; + break; + } + } + default: + badline(R); + } + } + } + + static expr * +aholread(EdRead *R) +{ + int i, k; + expr_h *rvh; + char *s1; + FILE *nl = R->nl; + + k = getc(nl); + if (k < '1' || k > '9') + badline(R); + i = k - '0'; + while((k = getc(nl)) != ':') { + if (k < '0' || k > '9') + badline(R); + i = 10*i + k - '0'; + } + rvh = (expr_h *)mem(memadj(sizeof(expr_h) + i)); + for(s1 = rvh->sym;;) { + if ((k = getc(nl)) < 0) { + fprintf(Stderr, + "Premature end of file in aholread, line %ld of %s\n", + R->Line, R->asl->i.filename_); + exit_ASL(R,1); + } + if (k == '\n') { + R->Line++; + if (!i) + break; + } + if (--i < 0) + badline(R); + *s1++ = k; + } + *s1 = 0; + rvh->op = (efunc*)f_OPHOL; + rvh->a = nv1; + return (expr *)rvh; + } + + static expr * +bholread(EdRead *R) +{ + int i; + expr_h *rvh; + char *s; + + if (xscanf(R, "%d", &i) != 1) + badline(R); + rvh = (expr_h *)mem(memadj(sizeof(expr_h) + i)); + s = rvh->sym; + if (fread(s, i, 1, R->nl) != 1) + badline(R); + s[i] = 0; + rvh->op = (efunc*)f_OPHOL; + rvh->a = nv1; + for(;;) switch(*s++) { + case 0: goto break2; /* so we return at end of fcn */ + case '\n': R->Line++; + } + break2: + return (expr *)rvh; + } + + static void +hv2init_ignore(ASL *asl, int hid_limit, int nobj, real *ow, real *y) {} + + int +fgh_read_ASL(ASL *a, FILE *nl, int flags) +{ + EdRead ER, *R; + Jmp_buf JB; + cgrad *cg, **cgp; + char fname[128]; + expr_v *e; + func_info *fi; + int i, j, k, *ka, maxfwd1, nc, nc0, nc1, nco, ncom, nlcon, nlin; + int no, nv, nvc, nvo, nvr, nxv, readall; + int (*Xscanf)(EdRead*, const char*, ...); + ograd *og, **ogp; + real t; + size_t *kaz, nz; + unsigned x; + + ASL_CHECK(a, ASL_read_fgh, "fgh_read"); + flagsave_ASL(a, flags); /* includes allocation of LUv, LUrhs, A_vals or Cgrad, etc. */ + asl = (ASL_fgh*)a; +#define asl ((ASL_fgh*)a) + + ed_reset(a); + R = EdReadInit_ASL(&ER, a, nl, 0); + if (flags & ASL_return_read_err) { + a->i.err_jmp_ = &JB; + i = setjmp(JB.jb); + if (i) { + a->i.err_jmp_ = 0; + return i; + } + } + nlcon = a->i.n_lcon_; + if (nlcon && !(flags & ASL_allow_CLP)) { + if (a->i.err_jmp_) + return ASL_readerr_CLP; + sorry_CLP(R, "logical constraints"); + } + + readall = flags & ASL_keep_all_suffixes; + if (nfunc) + func_add(a); + if (binary_nl) + holread = bholread; + else + holread = aholread; + + ncom = comb + comc + como + comc1 + como1; + nc0 = n_con; + nc = nc1 = nc0 + a->i.nsufext[ASL_Sufkind_con]; + no = n_obj; + nvc = c_vars; + nvo = o_vars; + nco = nc + no + nlcon; + if (no < 0 || nco <= 0) + scream(R, ASL_readerr_corrupt, + "ed2read: nc = %d, no = %d, nlcon = %d\n", + nc0, no, nlcon); + nxv = a->i.nsufext[ASL_Sufkind_var]; + nvr = n_var; /* nv for reading */ + nv0 = nv1 = nvr + nxv; + max_var = nv = nv1 + ncom; + combc = comb + comc; + ncom0 = ncom_togo = combc + como; + nzclim = ncom0 >> 3; + ncom1 = comc1 + como1; + nv0b = nv1 + comb; + nv0c = nv0b + comc; + nv01 = nv1 + ncom0; + amax = lasta = lasta0 = lasta00 = nv1 + 1; + ka = 0; + kaz = 0; + if ((maxfwd1 = maxfwd + 1) > 1) + nvref = maxfwd1*((ncom0 < vrefGulp ? ncom0 : vrefGulp) + 1); + x = nco*sizeof(cde) + no*sizeof(ograd *) + + nv*(sizeof(expr_v) + 2*sizeof(int)) + + ncom0*sizeof(cexp) + + ncom1*sizeof(cexp1) + + nfunc*sizeof(func_info *) + + nvref*sizeof(int) + + no; + nvar0 = a->i.n_var0; + if (!(nvinc = a->i.nvinc = a->i.n_var_ - nvar0 + nxv)) + nvar0 += ncom0 + ncom1; + if (pi0) { + memset(pi0, 0, nc*sizeof(real)); + if (havepi0) + memset(havepi0, 0, nc); + } + if (X0) + memset(X0, 0, nvr*sizeof(real)); + if (havex0) + memset(havex0, 0, nvr); + e = var_e = (expr_v *)M1zapalloc(x); + var_ex = e + nv0; + var_ex1 = var_ex + ncom0; + con_de = (cde *)(e + nv); + lcon_de = con_de + nc; + for(k = 0; k < nv; e++) { + e->op = f_OPVARVAL; + e->a = k++; + } + obj_de = lcon_de + nlcon; + Ograd = (ograd **)(obj_de + no); + cexps = (cexp *)(Ograd + no); + cexpsc = cexps + comb; + cexpso = cexpsc + comc; + cexps1 = (cexp1 *)(cexpse = cexps + ncom0); + funcs = (func_info **)(cexps1 + ncom1); + zc = (int *)(funcs + nfunc); + zci = zc + nv; + vrefx = zci + nv; + objtype = (char *)(vrefx + nvref); + if (nvref) { + vrefnext = vrefx + maxfwd1; + nvref -= maxfwd1; + } + if (n_cc && !cvar) + cvar = (int*)M1alloc(nc*sizeof(int)); + if (cvar) + memset(cvar, 0, nc*sizeof(int)); + last_d = 0; + nz = 0; + Xscanf = xscanf; + for(;;) { + ER.can_end = 1; + i = edag_peek(R); + if (i == EOF) { + free(imap); + adjoints = (real *)M1zapalloc(amax*Sizeof(real)); + adjoints_nv1 = &adjoints[nv1]; + adjust(a, flags); + nzjac = nz; + if (!Lastx) + Lastx = (real *)M1alloc(nv0*sizeof(real)); + fclose(nl); + nderps += nderp; + a->p.Objval = a->p.Objval_nomap = obj2val_ASL; + a->p.Objgrd = a->p.Objgrd_nomap = obj2grd_ASL; + a->p.Conval = con2val_ASL; + a->p.Jacval = jac2val_ASL; + a->p.Hvcomp = a->p.Hvcomp_nomap = hv2comp_ASL; + a->p.Hvcompd = hv2compd_ASL; + a->p.Hvcomps = hv2comps_ASL; + a->p.Hvinit = a->p.Hvinit_nomap = hv2init_ignore; + a->p.Conival = con2ival_ASL; + a->p.Conival_nomap = con2ival_nomap_ASL; + a->p.Congrd = con2grd_ASL; + a->p.Congrd_nomap = con2grd_nomap_ASL; + a->p.Lconval = lcon2val_ASL; + a->p.Xknown = x2known_ASL; + return prob_adj_ASL(a); + } + ER.can_end = 0; + k = -1; + switch(i) { + case 'C': + Xscanf(R, "%d", &k); + if (k < 0 || k >= nc0) + badline(R); + co_read(R, con_de + k, 1); + break; + case 'F': + if (Xscanf(R, "%d %d %d %127s", + &i, &j, &k, fname) != 4 + || i < 0 || i >= nfunc) + badline(R); + if ((fi = func_lookup(a, fname,0))) { + if (fi->nargs != k && fi->nargs >= 0 + && (k >= 0 || fi->nargs < -(k+1))) + scream(R, ASL_readerr_argerr, + "function %s: disagreement of nargs: %d and %d\n", + fname,fi->nargs, k); + } + else { + fi = (func_info *)mem(sizeof(func_info)); + fi->ftype = j; + fi->nargs = k; + fi->funcp = 0; + fi->name = (Const char *)strcpy((char*) + mem(memadj(strlen(fname)+1)), + fname); + } + if (!fi->funcp && !(fi->funcp = dynlink(fname))) + scream(R, ASL_readerr_unavail, + "function %s not available\n", + fname); + funcs[i] = fi; + break; + case 'G': + if (Xscanf(R, "%d %d", &j, &k) != 2 + || j < 0 || j >= no || k <= 0 || k > nvo) + badline(R); + ogp = Ograd + j; + while(k--) { + *ogp = og = (ograd *)mem(sizeof(ograd)); + ogp = &og->next; + if (Xscanf(R, "%d %lf", &i, &og->coef) != 2) + badline(R); + og->varno = i; + } + *ogp = 0; + break; + case 'J': + if (Xscanf(R, "%d %d", &j, &k) != 2 + || j < 0 || j >= nc0 || k <= 0 || k > nvc) + badline(R); + nz += k; + if (A_vals) { + j += Fortran; + if (ka) { + while(k--) { + if (Xscanf(R, "%d %lf", + &i, &t) != 2) + badline(R); + i = ka[i]++; + A_vals[i] = t; + A_rownos[i] = j; + } + } + else { + while(k--) { + if (Xscanf(R, "%d %lf", + &i, &t) != 2) + badline(R); + i = kaz[i]++; + A_vals[i] = t; + A_rownos[i] = j; + } + } + break; + } + cgp = Cgrad + j; + j = 0; + while(k--) { + *cgp = cg = (cgrad *)mem(sizeof(cgrad)); + cgp = &cg->next; + if (k_seen) { + if (Xscanf(R, "%d %lf", &i, &cg->coef) != 2) + badline(R); + } + else + if (Xscanf(R, "%d %d %lf", &i, &j, + &cg->coef) != 3) + badline(R); + cg->varno = i; + cg->goff = j; + } + *cgp = 0; + break; + case 'L': + Xscanf(R, "%d", &k); + if (k < 0 || k >= nlcon) + badline(R); + co_read(R, lcon_de + k, 0); + break; + case 'O': + if (Xscanf(R, "%d %d", &k, &j) != 2 + || k < 0 || k >= no) + badline(R); + objtype[k] = j; + co_read(R, obj_de + k, 1); + break; + case 'V': + if (Xscanf(R, "%d %d %d", &k, &nlin, &j) != 3) + badline(R); + if (k >= nvar0) + k += nvinc; + if (k < nv0 || k >= nv) + badline(R); + if (j) + cexp1_read(R, j, k, nlin); + else + cexp_read(R, k, nlin); + break; + case 'S': + Suf_read_ASL(R, readall); + break; + case 'r': + br_read(R, asl->i.n_con0, LUrhs, Urhsx, cvar, nv0); + break; + case 'b': + br_read(R, asl->i.n_var0, LUv, Uvx, 0, 0); + break; + case 'K': + case 'k': + k_seen++; + if (ka_read_ASL(a, R, i, &ka, &kaz)) + badline(R); + break; + case 'x': + if (!Xscanf(R,"%d",&k) + || k < 0 || k > nvr) + badline(R); + if (!X0 && want_xpi0 & 1) { + x = nv1*sizeof(real); + if (want_xpi0 & 4) + x += nv1; + X0 = (real *)M1zapalloc(x); + if (want_xpi0 & 4) + havex0 = (char*)(X0 + nv1); + } + while(k--) { + if (Xscanf(R, "%d %lf", &j, &t) != 2) + badline(R); + if (X0) { + X0[j] = t; + if (havex0) + havex0[j] = 1; + } + } + break; + case 'd': + if (!Xscanf(R,"%d",&k) + || k < 0 || k > nc0) + badline(R); + if (!pi0 && want_xpi0 & 2) { + x = nc*sizeof(real); + if (want_xpi0 & 4) + x += nc; + pi0 = (real *)M1zapalloc(x); + if (want_xpi0 & 4) + havepi0 = (char*)(pi0 + nc); + } + while(k--) { + if (Xscanf(R, "%d %lf", &j, &t) != 2 + || j < 0 || j >= nc0) + badline(R); + if (pi0) { + pi0[j] = t; + if (havepi0) + havepi0[j] = 1; + } + } + break; + default: + badline(R); + } + } + } +#ifdef __cplusplus + } +#endif diff --git a/src/solvers2/pfg_read.c b/src/solvers2/pfg_read.c new file mode 100644 index 0000000..4624c86 --- /dev/null +++ b/src/solvers2/pfg_read.c @@ -0,0 +1,5007 @@ +/******************************************************************* +Copyright (C) 2017, 2020 AMPL Optimization, Inc.; written by David M. Gay. + +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. + +The author and AMPL Optimization, Inc. disclaim all warranties with +regard to this software, including all implied warranties of +merchantability and fitness. In no event shall the author 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. +*******************************************************************/ + +#ifdef DEBUG +#include "assert.h" +#else +#define assert(x) /*nothing*/ +#endif + +#ifdef PSHVREAD +#include "jacpdim.h" +#define efunc efunc2 +#include "opnos.hd" +#define PSHV(x) x +#define asltype ASL_read_pfgh +#define who "pfgh_read" +#define ASLTYPE ASL_pfgh +#else +#define PSHV(x) /*nothing*/ +#include "asl_pfg.h" +#define asltype ASL_read_pfg +#define who "pfg_read" +#define ASLTYPE ASL_pfg +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +extern real conpival_nomap_ASL(ASL*, int, real*, fint*); +extern void conpgrd_nomap_ASL(ASL*, int, real *, real*, fint*); + +#define GAP_MAX 10 + +#ifdef PSHVREAD + static void hes_setup ANSI((ASL*, int, int, int, int, int)); +/* #define Static to avoid confusion in pi. */ +#define Static Static_pshv +#else +#define Static Static_psed +#endif +typedef struct Static Static; + static ograd *cotermwalk ANSI((Static*, expr **ep, ps_func *f, int wantg, int omitdv)); + static efunc *OPNUM, *OPVARVAL; + +#define NFHASH 23 + +#undef f_OPNUM +#include "r_opn0.hd" + +#undef nzc + +#if defined(IEEE_MC68k) || defined(IBM) + enum { W0 = 0, W1 = 1 }; +#else + enum { W0 = 1, W1 = 0 }; +#endif + +/*******/ + + typedef struct +expr_vx { + expr_v v; + linarg *la; + int a0; + int a1; + } expr_vx; + + typedef struct +Elemtemp { + unsigned int esize; + int nmax; + int k; + void **mp; + } Elemtemp; + + typedef struct +PSfind { + ps_func *f; + Elemtemp *b, *g; + } PSfind; + + typedef struct EU EU; + struct EU { efunc_n *op; union {real v; EU *p; } u; }; + + struct +Static { + ASLTYPE *asl; + ASL *a; + Elemtemp *_last_psb_elem; + derp *_last_d; + expr *(*_holread) ANSI((EdRead*)); + expr *_expr_free, **_slscratch; /* used in awalk(sumlist) */ + expr *_last_e; + expr_if *_iflist, *_if2list, *_if2list_end; + EU *_expr_n_free; + expr_v **_larvlist; + expr_v **_varp; + expr_va *_varglist, *_varg2list, *_varg2list_end; + int *_imap, *_vrefnext, *_vrefx; + int *_zc, *zc1, *_zci, *zci1, *_zl; + size_t _lthashmask, _nrange, _rangehashmask; + int klthash, krangehash; + int _allJ; + int _amax1; + int _cexp_k; + int _cexp_n; + int _conno; + int _groupno; + int _imap_len; + int _k_Elemtemp; + int _k_seen; + int _kimap; + int _kzc; + int _lasta; + int _lasta0; + int _lasta00; + int _max_var; + int _ncom; + int _ncom_togo; + int _nderp; + int _noa; + int _nocopy; + int _nndv; + int _nsce; + int _nv0b; + int _nv0c; + int _nv0r; + int _nv0x; + int _nv1; + int _nvref; + int _nzc; + int _nzclim; + int _slmax; + int _slscratchlev; + int _termno; /* current term number */ + int _wantCgroups; + int _wantOgroups; + int _zc_lim; + int nvar0, nvinc; + int size_exprn; + la_ref *_laref_free; + linarg **_lthash; /* hash table */ + linarg *_ltfree; /* free list */ + linarg *_tlist; /* linargs in this term */ + ograd *_freeog; /* free list */ + range **_rangehash; + real *_rnz; + real _lt_scale; + relo *_relolist, *_relo2list; + EdRead *R; + }; + +#define allJ S->_allJ +#define amax1 S->_amax1 +#define cexp_k S->_cexp_k +#define cexp_n S->_cexp_n +#define Conno S->_conno +#define expr_free S->_expr_free +#define expr_n_free S->_expr_n_free +#define freeog S->_freeog +#define Groupno S->_groupno +#define holread S->_holread +#define if2list S->_if2list +#define if2list_end S->_if2list_end +#define iflist S->_iflist +#define imap S->_imap +#define imap_len S->_imap_len +#define k_Elemtemp S->_k_Elemtemp +#define k_seen S->_k_seen +#define kimap S->_kimap +#define kzc S->_kzc +#define laref_free S->_laref_free +#define larvlist S->_larvlist +#define last_d S->_last_d +#define last_e S->_last_e +#define last_psb_elem S->_last_psb_elem +#define lasta S->_lasta +#define lasta0 S->_lasta0 +#define lasta00 S->_lasta00 +#define lt_scale S->_lt_scale +#define ltfree S->_ltfree +#define lthash S->_lthash +#define lthashmask S->_lthashmask +#define max_var S->_max_var +#define max_var1 asl->P.max_var1_ +#define Ncom S->_ncom +#define ncom_togo S->_ncom_togo +#define nderp S->_nderp +#define noa S->_noa +#define nocopy S->_nocopy +#define nndv S->_nndv +#define nrange S->_nrange +#define nsce S->_nsce +#define nv0b S->_nv0b +#define nv0c S->_nv0c +#define nv0r S->_nv0r +#define nv0x S->_nv0x +#define nv1 S->_nv1 +#define nvref S->_nvref +#define nzc S->_nzc +#define nzclim S->_nzclim +#define rangehash S->_rangehash +#define rangehashmask S->_rangehashmask +#define relo2list S->_relo2list +#define relolist S->_relolist +#define rnz S->_rnz +#define slmax S->_slmax +#define slscratch S->_slscratch +#define slscratchlev S->_slscratchlev +#define Termno S->_termno +#define tlist S->_tlist +#define varg2list S->_varg2list +#define varg2list_end S->_varg2list_end +#define varglist S->_varglist +#define varp S->_varp +#define vrefnext S->_vrefnext +#define vrefx S->_vrefx +#define wantCgroups S->_wantCgroups +#define wantOgroups S->_wantOgroups +#define zc S->_zc +#define zc_lim S->_zc_lim +#define zci S->_zci +#define zl S->_zl + +static list *crefs ANSI((Static*)); + + static Static * +S_init(Static *S, ASLTYPE *asl) +{ + memset(S, 0, sizeof(Static)); + S->asl = asl; + S->a = (ASL*)asl; + return S; + } + + static void * +new_mblkzap(ASLTYPE *asl, int k) +{ + void *rv = new_mblk_ASL((ASL*)asl,k); + memset(rv, 0, sizeof(void*)<i.memLast = asl->i.memNext = 0; + memset(asl->mblk_free, 0, MBLK_KMAX*sizeof(char*)); + memset(&asl->P.merge + 1, 0, sizeof(ps_info) - sizeof(long)); +#ifdef PSHVREAD + asl->P.hop_free = 0; + asl->P.hes_setup_called = 0; +#endif + } + +#ifdef Double_Align +#define memadj(x) x +#else +#define memadj(x) (((x) + (sizeof(long)-1)) & ~(sizeof(long)-1)) +#endif +#ifdef DEBUG + static expr *opzork; + static derp *dzork; + static ograd *ogzork; + static expr *ezork; + static EU *enzork; + static int dzork1, dzork2, exprzork, exprzork1, izork = -1; +#endif + + static Elemtemp * +new_Elemtemp(Static *S, unsigned int esize, void **mp) +{ + Elemtemp *e; + int k; + ASL *asl = S->a; + + e = (Elemtemp *)new_mblk(k_Elemtemp); + e->esize = esize; + e->mp = mp; + e->k = k = htcl(8*esize); + *mp = new_mblk(k); + e->nmax = (sizeof(void*) << k) / esize; + return e; + } + + static void +del_Elemtemp(Static *S, Elemtemp *e) +{ + ASL *asl = S->a; + del_mblk(e->k, *e->mp); + del_mblk(k_Elemtemp, e); + } + + static void +upgrade_Elemtemp(Static *S, Elemtemp *e) +{ + void *m, *m0; + int k; + ASL *asl = S->a; + + k = e->k++; + memcpy(m = new_mblk(e->k), m0 = *e->mp, e->esize * e->nmax); + del_mblk(k++, m0); + *e->mp = m; + e->nmax = (sizeof(void*) << k) / e->esize; + } + + static void +free_laref(Static *S, la_ref **L) +{ + la_ref *L1, *L2; + + if ((L1 = *L)) { + while((L2 = L1->next)) + L1 = L2; + L1->next = laref_free; + laref_free = *L; + *L = 0; + } + } + + static la_ref * +new_laref(Static *S, la_ref *nxt) +{ + la_ref *rv; + + if ((rv = laref_free)) + laref_free = rv->next; + else + rv = (la_ref *)mem_ASL(S->a, sizeof(la_ref)); + rv->next = nxt; + return rv; + } + + static list * +new_list(ASL *asl, list *nxt) +{ + list *rv = (list *)mem(sizeof(list)); + rv->next = nxt; + return rv; + } + + static void +fscream(EdRead *R, const char *name, int nargs, const char *kind) +{ + scream(R, ASL_readerr_argerr, + "line %ld: attempt to call %s with %d %sargs\n", + R->Line, name, nargs, kind); + } + + static void +new_derp(Static *S, int a, int b, real *c) +{ + derp *d; + nderp++; + d = (derp *)mem_ASL(S->a, sizeof(derp)); + d->next = last_d; + last_d = d; + d->a.i = a; + d->b.i = b; + d->c.rp = c; +#ifdef DEBUG + if (d == dzork) + printf(""); + if (++dzork1 == dzork2) + printf(""); +#endif + } + + static derp * +new_relo(Static *S, expr *e, derp *Dnext, int *ap) +{ + relo *r; + derp *d; + + r = (relo *)mem_ASL(S->a, sizeof(relo)); + r->next = relolist; + r->next2 = relo2list; + relo2list = relolist = r; + if (last_d != Dnext) { + *ap = e->a; + for(d = last_d; d->next != Dnext; d = d->next); + d->next = 0; + } + else { + last_d = 0; + new_derp(S, e->a, *ap = lasta++, &edagread_one); + } + if (!last_d) + return 0; + r->D = r->Dcond = last_d; + r->Dnext = Dnext; + return r->D; + } + + static relo * +new_relo1(Static *S, derp *Dnext) +{ + relo *r; + + r = (relo *)mem_ASL(S->a, sizeof(relo)); + r->next = relolist; + relolist = r; + r->D = 0; + r->Dnext = Dnext; + return r; + } + + static void +efree(Static *S, expr *e) +{ + EU *en; + expr **args, **argse, *e1; + top: + switch(optypeb[Intcast e->op]) { + + case 2: /* binary */ + efree(S, e->R.e); + /* no break */ + + case 1: /* unary */ + e1 = e->L.e; + e->L.e = expr_free; + expr_free = e; + e = e1; + goto top; + + case 6: /* sumlist */ + args = e->L.ep; + argse = e->R.ep; + while(args < argse) + efree(S, *args++); + e->L.e = expr_free; + expr_free = e; + break; + + case 9: + en = (EU*)e; + en->u.p = expr_n_free; + expr_n_free = en; + } + } + + static expr * +new_expr(Static *S, int o, expr *L, expr *R) +{ + expr *rv; + +#ifdef DEBUG + if (++exprzork1 == exprzork) + printf("exprzork1 = %d\n", exprzork1);; +#endif + if ((rv = expr_free)) + expr_free = rv->L.e; + else + rv = (expr *)mem_ASL(S->a, sizeof(expr)); + PSHV(rv->dL2 = 0); + if (o == f_OPPOW) { + if (Intcast R->op == f_OPNUM) { + if (((expr_n *)R)->v == 2.) { + o = f_OP2POW; + R = 0; + PSHV(rv->dL2 = 2); + } + else { + o = f_OP1POW; + rv->dR = ((expr_n *)R)->v; + /* rv->R.e may be used for a back pointer */ + } + } + else if (Intcast L->op == f_OPNUM) + o = f_OPCPOW; + } + rv->op = (efunc *)(size_t)o; + rv->L.e = L; + rv->R.e = R; + rv->a = 0; +#ifdef DEBUG + if (rv == ezork) + printf(""); +#endif + return rv; + } + + static void +dexpr(Static *S, expr *e, expr *L, expr *R) +{ + int L1, R1; +#ifdef PSHVREAD + int i, opcode; +#endif + + e->a = noa; + L1 = L && L->op != OPNUM && L->a != noa; + R1 = R && R->op != OPNUM && R->a != noa; + if (L1 | R1) { + if (L1) + new_derp(S, L->a, lasta, &e->dL); + if (R1) + new_derp(S, R->a, lasta, &e->dR); + e->a = lasta++; +#ifdef PSHVREAD + e->bak = last_e; + last_e = e; + opcode = Intcast e->op; + if (R) + e->dLR = e->dR2 = 0; + if (L1) + if (R1) + switch(opcode) { + case PLUS: i = Hv_plusLR; break; + case MINUS: i = Hv_minusLR; break; + case MULT: i = Hv_timesLR; break; + case DIV: i = Hv_divLR; break; + default: i = Hv_binaryLR; + } + else switch(opcode) { + case PLUS: + case MINUS: i = Hv_plusL; break; + case DIV: + case MULT: i = Hv_timesL; break; + case UMINUS: i = Hv_negate; break; + default: i = Hv_unary; + } + else + switch(opcode) { + case PLUS: i = Hv_plusR; break; + case MINUS: i = Hv_minusR; break; + case MULT: i = Hv_timesR; break; + default: i = Hv_binaryR; + } + e->dO.i = i; +#endif + } + } + + static real dvalue[] = { +#include "dvalue.hd" + }; + + static expr_n * +new_expr_n(Static *S, real t) +{ + EU *en; + + if ((en = expr_n_free)) + expr_n_free = en->u.p; + else + en = (EU *)mem_ASL(S->a, S->size_exprn); +#ifdef DEBUG + if (en == enzork) + printf(""); +#endif + en->u.v = t; + en->op = (efunc_n *)f_OPNUM; + return (expr_n*)en; + } + + static expr * +eread(EdRead *R) +{ + ASLTYPE *asl; + Static *S; + arglist *al; + argpair *ap, *sap; + char **sa; + de *d; + expr *L, *arg, **args, **argse, *rv; + expr_f *rvf; + expr_if *rvif; + expr_va *rva; + fint L1; + func_info *fi; + int *at, i, j, k, ks, numargs, op, symargs; + int (*Xscanf)(EdRead*, const char*, ...); + plterm *p; + real *b, *ra; + real r; + + S = (Static *)R->S; + asl = S->asl; + Xscanf = xscanf; + switch(edag_peek(R)) { + case 'f': + if (Xscanf(R, "%d %d", &i, &j) != 2) + badline(R); + fi = funcs[i]; + if (fi->nargs >= 0) { + if (fi->nargs != j) { + bad_nargs: + fscream(R, fi->name, j, ""); + } + } + else if (-(1+fi->nargs) > j) + goto bad_nargs; + rvf = (expr_f *)mem(sizeof(expr_f) + + (j-1)*sizeof(expr *)); + rvf->op = (efunc *)f_OPFUNCALL; + rvf->fi = fi; + PSHV(rvf->dO.i = Hv_func); + args = rvf->args; + argse = args + j; + k = ks = symargs = numargs = 0; + while(args < argse) { + arg = *args++ = eread(R); + if ((op = Intcast arg->op) == f_OPHOL) + symargs++; + else if (op == f_OPIFSYM) + ks++; + else { + numargs++; + if (op != f_OPNUM) + k++; + } + } + symargs += ks; + if (symargs && !(fi->ftype & 1)) + fscream(R, fi->name, symargs, "symbolic "); + ra = (real *)mem(sizeof(arglist) + + (k+ks)*sizeof(argpair) + + numargs*sizeof(real) + + symargs*sizeof(char *) + + j*sizeof(int)); + al = rvf->al = (arglist *)(ra + numargs); + al->n = numargs + symargs; + al->nr = numargs; + al->ra = ra; + al->derivs = al->hes = 0; + al->funcinfo = fi->funcinfo; + al->AE = asl->i.ae; + al->sa = (Const char**)(sa = (char **)(al + 1)); + ap = rvf->ap = (argpair *)(sa + symargs); + sap = rvf->sap = ap + k; + at = al->at = (int *)(sap + ks); + symargs = numargs = 0; + for(args = rvf->args; args < argse; at++) { + arg = *args++; + if ((op = Intcast arg->op) == f_OPHOL) { + *at = --symargs; + *sa++ = ((expr_h *)arg)->sym; + } + else if (op == f_OPIFSYM) { + *at = --symargs; + sap->e = arg; + (sap++)->u.s = sa++; + } + else { + *at = numargs++; + if (op == f_OPNUM) + *ra = ((expr_n *)arg)->v; + else { + ap->e = arg; + (ap++)->u.v = ra; + } + ra++; + } + } + rvf->ape = ap; + rvf->sape = sap; + return (expr *)rvf; + + case 'h': + return holread(R); + case 's': + if (Xscanf(R, "%hd", &L1) != 1) + badline(R); + r = L1; + goto have_r; + + case 'l': + if (Xscanf(R, "%ld", &L1) != 1) + badline(R); + r = L1; + goto have_r; + + case 'n': + if (Xscanf(R, "%lf", &r) != 1) + badline(R); + have_r: + return (expr *)new_expr_n(S, r); + + case 'o': + break; + + case 'v': + if (Xscanf(R,"%d",&k) != 1 || k < 0) + badline(R); + if (k >= S->nvar0) + k += S->nvinc; + if (k > max_var) + badline(R); + return (expr *)(var_e + k); + + default: + badline(R); + } + + if (Xscanf(R, asl->i.opfmt, &k) != 1 || k < 0 || k >= N_OPS) + badline(R); + switch(optype[k]) { + + case 1: /* unary */ + rv = new_expr(S, k, eread(R), 0); + return rv; + + case 2: /* binary */ + L = eread(R); + rv = new_expr(S, k, L, eread(R)); + return rv; + + case 3: /* vararg (min, max) */ + i = -1; + Xscanf(R, "%d", &i); + if (i <= 0) + badline(R); + rva = (expr_va *)mem(sizeof(expr_va)); + PSHV(rva->val = 0;) + rva->op = (efunc *)(size_t)k; + rva->L.d = d = (de *)mem(i*sizeof(de) + sizeof(expr *)); + rva->next = varglist; + varglist = rva; + for(j = 0; i > 0; i--, d++) + d->e = eread(R); + d->e = 0; /* sentinnel expr * */ + return (expr *)rva; + + case 4: /* piece-wise linear */ + i = -1; + Xscanf(R, "%d", &i); + if (i <= 1) + badline(R); + plterms++; + j = 2*i - 1; + p = (plterm *)mem(sizeof(plterm) + (j-1)*sizeof(real)); + p->n = i; + b = p->bs; + do { + switch(edag_peek(R)) { + case 's': + if (Xscanf(R,"%hd",&L1) != 1) + badline(R); + r = L1; + break; + case 'l': + if (Xscanf(R,"%ld",&L1) != 1) + badline(R); + r = L1; + break; + case 'n': + if (Xscanf(R,"%lf",&r) == 1) + break; + default: + badline(R); + } + *b++ = r; + } + while(--j > 0); + if (b[-2] <= 0.) + p->z = 2*i - 2; + else { + b = p->bs + 1; + while(*b <= 0.) + b += 2; + p->z = (b - p->bs) - 1; + } + rv = (expr *)mem(sizeof(expr)); + rv->op = (efunc *)(size_t)k; + rv->L.p = p; + rv->R.e = eread(R); + return rv; + + case 5: /* if */ + rvif = (expr_if *)mem(sizeof(expr_if)); + PSHV(rvif->val = 0;) + PSHV(rvif->Te = rvif->Fe = 0;) /* for bogus .nl files */ + rvif->op = (efunc *)(size_t)k; + rvif->next = iflist; + iflist = rvif; + rvif->e = eread(R); + rvif->T = L = eread(R); + rvif->F = L = eread(R); + return (expr *)rvif; + + case 11: /* OPCOUNT */ + case 6: /* sumlist */ + i = 0; + Xscanf(R, "%d", &i); + if (i <= 2 && (optype[k] == 6 || i < 1)) + badline(R); + if (slmax < i) + slmax = i; + args = (expr**)new_mblk(htcl(i*sizeof(expr*))); + rv = new_expr(S, k, (expr*)args, (expr*)(args+i)); + do + *args++ = eread(R); + while(--i > 0); + return rv; + } + badline(R); + return 0; + } + + static void +ogfree(Static *S, ograd *og) +{ + ograd *og1; + + if ((og1 = og)) { +#ifdef DEBUG + if (ogzork) { + do if (og == ogzork) + printf(""); + while((og = og->next)); + og = og1; + } +#endif + while(og1->next) + og1 = og1->next; + og1->next = freeog; + freeog = og; + } + } + + static real +ogfree1(Static *S, ograd **ogp) +{ + ograd *og = *ogp; + *ogp = og->next; + og->next = freeog; + freeog = og; + return og->coef; + } + + static ograd * +new_ograd(Static *S, ograd *next, int varno, real coef) +{ + ograd *rv; + + if ((rv = freeog)) + freeog = rv->next; + else + rv = (ograd *)mem_ASL(S->a, sizeof(ograd)); + rv->next = next; + rv->varno = varno; + rv->coef = coef; +#ifdef DEBUG + if (rv == ogzork) + printf(""); +#endif + return rv; + } + + static linarg * +lahash(Static *S, linarg *la) +{ + ASLTYPE *asl; + Ulong x; + int k, nnz; + linarg *la1, *la2, **lap, **lap1, **q, **q0, **qe; + ograd *og, *og1; + size_t mask; + union {real r; Ulong L[2];} u; + + asl = S->asl; + x = nnz = la->nnz; + for(og = la->nz; og; og = og->next) { + u.r = og->coef; + x = (x << 1 | (x & 0x80000000) >> 31) ^ + (og->varno*101 + u.L[W0] + 257*u.L[W1]); + } + for(lap = <hash[x & lthashmask]; (la1 = *lap); lap = &la1->hnext) + if (la1->nnz == nnz) { + og = la->nz; + og1 = la1->nz; + for(;;) { + if (!og) { + if (!og1) + return la1; + break; + } + if (!og1) + break; + if (og->varno != og1->varno + || og->coef != og1->coef) + break; + og = og->next; + og1 = og1->next; + } + } + *lap = la; + la->hnext = 0; + if (++S->asl->P.nlttot > lthashmask) { + mask = lthashmask; + k = S->klthash; + q = q0 = lthash; + qe = q + mask; + lthashmask = mask = (mask << 1) | 1; + lap = lthash = (linarg**)new_mblkzap(asl, S->klthash = k + 1); + while(q <= qe) { + for(la2 = *q++; la2; la2 = la1) { + la1 = la2->hnext; + x = la2->nnz; + for(og = la2->nz; og; og = og->next) { + u.r = og->coef; + x = (x << 1 | (x & 0x80000000) >> 31) ^ + (og->varno*101 + u.L[W0] + 257*u.L[W1]); + } + lap1 = lap + (x & mask); + la2->hnext = *lap1; + *lap1 = la2; + } + } + del_mblk(k, q0); + } + return la; + } + + static range * +new_range(Static *S, range *r, range **rp) +{ + ASLTYPE *asl; + int k, len, uilen; + range *r1, *r2, *r3, **rp1, **rq, **rq0, **rqe; + size_t mask; + + asl = S->asl; + uilen = r->nv*sizeof(int); + len = sizeof(range) + uilen; + r1 = (range*)mem(len); + r1->nintv = 0; + r1->n = r->n; + r1->nv = r->nv; + r1->chksum = r->chksum; + r1->refs = 0; + r1->lasttermno = -1; + r1->hnext = r1->hunext = 0; + if (uilen) + memcpy(r1->ui = (int*)(r1+1), r->ui, uilen); + r1->lap = (linarg**)new_mblk(htcl(len = r->n*sizeof(linarg*))); + memcpy(r1->lap, r->lap, len); + r2 = r1->rlist.next = asl->P.rlist.next; + r1->rlist.prev = r2->rlist.prev; + r2->rlist.prev = asl->P.rlist.next = r1; + *rp = r1; + if (++nrange > rangehashmask) { + mask = rangehashmask; + k = S->krangehash; + rq = rq0 = rangehash; + rqe = rq + mask; + rangehashmask = mask = (mask << 1) | 1; + rp = rangehash = (range**)new_mblkzap(asl, S->krangehash = k + 1); + while(rq <= rqe) { + for(r2 = *rq++; r2; r2 = r3) { + r3 = r2->hunext; + rp1 = rp + (r2->chksum & mask); + r2->hunext = *rp1; + *rp1 = r2; + } + } + del_mblk(k, rq0); + } + return r1; + } + + static int +lacompar(const void *a, const void *b, void *v) +{ + ograd *oa, *ob; + int i; + real t; + + if (a == b) + return 0; + Not_Used(v); + oa = (*(linarg **)a)->nz; + ob = (*(linarg **)b)->nz; + for(;;) { + if (!oa) { + if (ob) + return -1; + break; + } + if (!ob) + return 1; + if ((i = oa->varno - ob->varno)) + return i; + if ((t = oa->coef - ob->coef)) + return t > 0. ? 1 : -1; + oa = oa->next; + ob = ob->next; + } + return 0; + } + + static int +ndiff(range *r, range *r1) +{ + /* Return 0 if r and r1 have the same linargs; else return 1. */ + + linarg **la, **la1, **la1e, **lae; + + la = r->lap; + lae = la + r->n; + la1 = r1->lap; + la1e = la1 + r1->n; + for(;;++la, ++la1) { + if (la >= lae) + return la1 < la1e; + if (la1 >= la1e) + return 1; + if (lacompar((char*)la, (char *)la1, NULL)) + return 1; + } + return 0; + } + + static range * +uhash(Static *S, range *r) +{ + ASLTYPE *asl; + int len, n, nv, *ui, *uie; + range *r1, **rp; + size_t L; + + asl = S->asl; + L = 0; + nv = r->nv; + ui = r->ui; + uie = ui + nv; + len = sizeof(int)*nv; + while (ui < uie) + L = 37*L + *ui++; + r->chksum = L; + ui = r->ui; + rp = rangehash + (L & rangehashmask); + n = r->n; + if (asl->P.merge) + while((r1 = *rp)) { + if (r1->nv == nv && r1->n == n && (len == 0 || !memcmp(ui, r1->ui, len))) { + if (n == 0 || !memcmp(r->lap, r1->lap, n*sizeof(linarg*))) + return *rp; + if (!ndiff(r, r1)) + return *rp; + } + rp = &r1->hunext; + } + return new_range(S, r, rp); + } + + static range * +rhash(Static *S, range *r, int addnew) +{ + ASLTYPE *asl; + int len, n; + linarg **lae, **lap; + ograd *og; + range *r1, **rp; + unsigned long L; + + asl = S->asl; + L = 0; + lap = r->lap; + n = r->n; + lae = lap + n; + len = n * sizeof(linarg*); + while(lap < lae) { + L *= 37; + for(og = (*lap++)->nz; og; og = og->next) + L = 101*L + og->varno; + } + r->chksum = L; + lap = r->lap; + rp = rangehash + (L & rangehashmask); + if (asl->P.merge) + while((r1 = *rp)) { + if (r1->n == n && !memcmp(lap, r1->lap, len)) + return r1; + rp = &r1->hnext; + } + if (addnew) + return new_range(S, r, rp); + return *rp; + } + + static int +compar(const void *a0, const void *b0, void *v) +{ + int a, b, c; + Static *S = (Static *)v; + + if ((a = *(int*)a0) >= max_var) { + a = zl[a-nv0x]; + if ((b = *(int*)b0) >= max_var) { + if ((c = a - zl[b-nv0x])) + return c; + return *(int*)a0 - b; + } + if (a == b) + return -1; + } + else if ((b = *(int*)b0) >= max_var) { + b = zl[b-nv0x]; + if (a == b) + return 1; + } + return a - b; + } + + static int +hscompar(const void *a0, const void *b0, void *v) +{ + int a, b, c; + Static *S = (Static *)v; + + if ((a = *(int*)a0) >= Ncom) { + a = zl[a]; + if ((b = *(int*)b0) >= Ncom) { + if ((c = a - zl[b])) + return c; + return *(int*)a0 - b; + } + a -= nv0x; + if (a == b) + return -1; + } + else if ((b = *(int*)b0) >= Ncom) { + b = zl[b] - nv0x; + if (a == b) + return 1; + } + return a - b; + } + + static void +zcsort(Static *S, int *c, int *ci, int i, int n, int p) +{ + int j; + + if (n < nzclim || p < 0) + qsortv(ci, n, sizeof(int), compar, S); + else for(j = 0; i < p; i++) + if (c[i]) + ci[j++] = i; + } + + static ograd* +compress(Static *S, ograd *og, real *cp, int *comvar) +{ + ograd *og1; + int i, ix, j, j1, k, nzc1, *zc1, *zci1; + real c, t; + + c = og->varno < 0 ? ogfree1(S,&og) : 0.; + ix = nzc1 = 0; + zc1 = S->zc1; + zci1 = S->zci1; + for(og1 = og; og1; og1 = og1->next) { + zc1[i = og1->varno] = 1; + zci1[nzc1++] = i; + rnz[i] = og1->coef; + if (ix < i) + ix = i; + } + if (ix < nv0x) { + *cp = c; + *comvar = 0; + for(og1 = og; og1; og1 = og1->next) + zc1[og1->varno] = 0; + return og; + } + *comvar = 1; + for(i = 0; i < nzc1; ) + if ((j = zci1[i]) < nv0x) + i++; + else { + if (!zc[j]++) + zci[nzc++] = j; + t = rnz[j]; + j1 = j - nv0x; + if ((og1 = (S->asl->P.dv + j1)->ll)) { + if (og1->varno < 0) { + c += t*og1->coef; + og1 = og1->next; + } + for(; og1; og1 = og1->next) + if (!zc1[k = og1->varno]++) { + zci1[nzc1++] = k; + rnz[k] = t*og1->coef; + } + else + rnz[k] += t*og1->coef; + } + zc1[j] = 0; + zci1[i] = zci1[--nzc1]; + } + *cp = c; + ogfree(S, og); + og = 0; + if (nzc1 > 0) { + zcsort(S, zc1, zci1, 0, nzc1, max_var); + while(nzc1 > 0) { + i = zci1[--nzc1]; + zc1[i] = 0; + if ((t = rnz[i])) { + og = new_ograd(S, og, i, t); + if (!zc[i]++) + zci[nzc++] = i; + } + } + } + return og; + } + + static linarg * +afree(Static *S, ograd *og, expr **ep) +{ + linarg *la, *la1, *rv; + real c, s, t; + ograd *og1, *ogx; + int comvar, nnz; + la_ref *refs; + ASLTYPE *asl = S->asl; + + rv = 0; + if (!og || !(og = compress(S,og,&c,&comvar))) + goto done; + + if ((la = ltfree)) + ltfree = la->hnext; + else { + la = (linarg *)mem(sizeof(linarg)); + la->refs = 0; + } + + s = og->coef; + if (s < 0) + s = -s; + la->nz = ogx = og1 = og; + nnz = 1; + while((og1 = og1->next)) { + nnz++; + t = og1->coef; + if (t < 0) + t = -t; + if (s < t) { + s = t; + ogx = og1; + } + } + + la->nnz = nnz; + if ((s = ogx->coef) != 1.) + for(og1 = og; og1; og1 = og1->next) + og1->coef /= s; + lt_scale = s; + la1 = lahash(S, la); + if (la1 == la) { + la->refs = 0; + la->v = 0; + la->termno = Termno; + la->tnext = tlist; + tlist = la; + la->lnext = asl->P.lalist; + asl->P.lalist = la; + la->hnext = 0; + } + else { + /* !!? Eventually create or adjust la1->v to save cycles. */ + /* This requires possible adjustment for differing */ + /* constants and scale factors. */ + if (la1->termno == Termno) + asl->P.ndupst++; /* research statistic */ + else { + free_laref(S, &la->refs); + la1->termno = Termno; + la1->tnext = tlist; + tlist = la1; + asl->P.ndupdt++; /* research statistic */ + } + ogfree(S, og); + la->hnext = ltfree; + ltfree = la; + la = la1; + } + if (ep && (nnz > 1 || comvar)) { + la->refs = refs = new_laref(S, la->refs); + refs->ep = ep; + refs->c = c; + refs->scale = s; + } + if (nnz > 1) + rv = la; + done: + return rv; + } + + static ograd * +af_sum(Static *S, ograd *Log, ograd *Rog) +{ + ograd *oL, *oR, *oR1, *og, **ogp; + + oL = Log; + oR = Rog; + ogp = &og; + for(;;) { + if (!oL) { + *ogp = oR; + break; + } + if (!oR) { + *ogp = oL; + break; + } + if (oL->varno > oR->varno) { + *ogp = oR; + ogp = &oR->next; + oR = *ogp; + } + else { + if (oL->varno == oR->varno) { + oL->coef += oR->coef; + oR1 = oR->next; + oR->next = 0; + ogfree(S, oR); + oR = oR1; + if (oL->coef == 0.) { + oR1 = oL->next; + oL->next = 0; + ogfree(S, oL); + oL = oR1; + continue; + } + } + *ogp = oL; + ogp = &oL->next; + oL = *ogp; + } + } + return og; + } + + static cgrad * +cf_sum(Static *S, cgrad *Lcg, ograd *Rog) +{ + cgrad *cg, **cgp, *oL; + ograd *oR, *oR1; + + oL = Lcg; + oR = Rog; + cgp = &cg; + cg = 0; + for(;;) { + if (!oL) { + assert(!oR); + break; + } + if (!oR) { + *cgp = oL; + break; + } + assert(oL->varno <= oR->varno); + if (oL->varno == oR->varno) { + oL->coef += oR->coef; + oR1 = oR->next; + oR->next = 0; + ogfree(S, oR); + oR = oR1; + } + *cgp = oL; + cgp = &oL->next; + oL = *cgp; + } + return cg; + } + + static void +sumlist_afree(Static *S, ograd *Laf, expr *e, expr **argso, expr **sls, expr **slscr) +{ + int nlin, nnl; + expr *e1, *e2, **ep, **ep1; + ASL *asl = S->a; + + nlin = sls - slscr; + ep = e->L.ep; + nnl = argso - ep; + switch(nlin) { + case 1: + e1 = slscr[0]; + break; + case 2: + e1 = new_expr(S, f_OPPLUS, slscr[0], slscr[1]); + e1->dL = e1->dR = 1.; + break; + default: + ep1 = (expr**)new_mblk(htcl(nlin*sizeof(expr*))); + e1 = new_expr(S, f_OPSUMLIST, (expr*)ep1, (expr*)(ep1+nlin)); + memcpy(ep1, slscr, nlin*sizeof(expr*)); + } + switch(nnl) { + case 1: + e2 = *ep; + break; + case 2: + e2 = new_expr(S, f_OPPLUS, ep[0], ep[1]); + e2->dL = e2->dR = 1.; + break; + default: + ep1 = (expr**)new_mblk(htcl(nnl*sizeof(expr*))); + e2 = new_expr(S, f_OPSUMLIST, (expr*)ep1, (expr*)(ep1+nnl)); + memcpy(ep1, ep, nnl*sizeof(expr*)); + } + del_mblk(htcl((e->R.ep - ep)*sizeof(expr*)), ep); + e->op = (efunc*)f_OPPLUS; + e->dL = e->dR = 1.; + e->L.e = e1; + e->R.e = e2; + afree(S, Laf, &e->L.e); + } + + static void +sdvcite(Static *S, int k) +{ + range *r; + split_ce *cs; + linarg *la, **lap, **lape; + + cs = S->asl->P.Split_ce + (k - max_var); + r = cs->r; + lap = r->lap; + lape = lap + r->n; + while(lap < lape) { + la = *lap++; + if (la->termno != Termno) { + free_laref(S, &la->refs); + la->termno = Termno; + la->tnext = tlist; + tlist = la; + } + } + } + + static ograd * +awalk(Static *S, expr *e) /* return 0 if e is not linear */ +{ + ASLTYPE *asl; + int k, k1, kscr; + expr *L, **args, **argse, **argso, **sls, **slscr; + expr_va *rva; + de *d; + expr_if *rvif; + expr_f *rvf; + ograd *Laf, *Raf, *rv, *taf; + linarg *la, **nl; + argpair *ap, *ape; + real t; + + switch(optypeb[k = Intcast e->op]) { + + case 1: /* unary */ + Laf = awalk(S, e->L.e); + if (k == f_OPUMINUS) { + if ((taf = Laf)) { + do taf->coef = -taf->coef; + while((taf = taf->next)); + return Laf; + } + } + if (Laf) + afree(S, Laf, &e->L.e); + break; + + case 2: /* binary */ + Laf = awalk(S, e->L.e); + Raf = awalk(S, e->R.e); + if (Laf && Raf) + switch(k) { + case f_OPMINUS: + taf = Raf; + do taf->coef = -taf->coef; + while((taf = taf->next)); + /* no break; */ + case f_OPPLUS: + return af_sum(S, Laf, Raf); + + case f_OPMULT: + if (Raf->varno < 0 && !Raf->next) { + swap: + taf = Laf; + Laf = Raf; + Raf = taf; + goto scale; + } + if (Laf->varno < 0 && !Laf->next) { + taf = Raf; + scale: t = Laf->coef; + if (t == 0.) { + ogfree(S, Raf); + return Laf; + } + do taf->coef *= t; + while((taf = taf->next)); + ogfree(S, Laf); + return Raf; + } + break; + + case f_OPDIV: + if (Raf->varno < 0 && !Raf->next) { + Raf->coef = 1. / Raf->coef; + goto swap; + } + } + afree(S, Laf, &e->L.e); + afree(S, Raf, &e->R.e); + break; + + case 3: /* vararg (min, max) */ + rva = (expr_va *)e; + for(d = rva->L.d; (L = d->e); d++) + afree(S, awalk(S,L), &d->e); + break; + + case 4: /* piece-wise linear */ + afree(S, awalk(S,e->R.e), &e->R.e); + break; + + case 5: /* if */ + rvif = (expr_if *)e; + afree(S, awalk(S,rvif->T), &rvif->T); + afree(S, awalk(S,rvif->F), &rvif->F); + break; + + case 6: /* sumlist */ + args = e->L.ep; + argse = e->R.ep; + k1 = argse - args; + while(!(Laf = awalk(S,*args++))) + if (args >= argse) + return 0; + kscr = -1; + asl = S->asl; + if (!slscratchlev++) + slscr = slscratch; + else { + kscr = htcl(k1*sizeof(expr*)); + slscr = (expr**)new_mblk(kscr); + } + sls = slscr; + argso = args - 1; + *sls++ = *argso; + while(args < argse) + if ((Raf = awalk(S,*args))) { + *sls++ = *args++; + Laf = af_sum(S, Laf, Raf); + } + else + *argso++ = *args++; + rv = 0; + if (argso == e->L.ep) { + rv = Laf; + goto delscratch; + } + sumlist_afree(S, Laf, e, argso, sls, slscr); + delscratch: + --slscratchlev; + if (kscr >= 0) + del_mblk(kscr, slscr); + return rv; + + case 7: /* function call */ + rvf = (expr_f *)e; + ap = rvf->ap; + for(ape = rvf->ape; ap < ape; ap++) + afree(S, awalk(S,ap->e), &ap->e); + + case 8: /* OPHOL (Hollerith) */ + break; + + case 9: /* OPNUM */ + return new_ograd(S, 0, -1, ((expr_n*)e)->v); + + case 10:/* OPVARVAL */ + asl = S->asl; + k = (expr_v *)e - var_e; + if ((k < 0 || k >= max_var) + && (k = ((expr_vx*)e)->a0) < 0) { + if ((la = ((expr_vx*)e)->la) + && la->termno != Termno) { + la->termno = Termno; + la->tnext = tlist; + tlist = la; + } + return 0; + } + if ((k1 = (k - nv0x)) < 0) { + if (!zc[k]++) + zci[nzc++] = k; + goto ogret; + } + if (k1 >= Ncom) { + if (!zc[k = ((expr_vx*)e)->a1]++) + zci[nzc++] = k; + sdvcite(S,k); + return 0; + } + if ((nl = (asl->P.dv + k1)->nl)) { + if (!zc[k]++) + zci[nzc++] = k; + while((la = *nl++)) + if (la->termno != Termno) { + la->termno = Termno; + la->tnext = tlist; + tlist = la; + } + return 0; + } + if (!zc[k]++) + zci[nzc++] = k; + ogret: + return new_ograd(S, 0, k, 1.); + + case 11: /* OPCOUNT */ + break; + + default: + scream(S->R, ASL_readerr_bug, + "awalk: unexpected optype[%d] = %d\n", + k, optype[k]); + } + return 0; + } + + static void +cexp_upgrade(Static *S, int t) +{ + int k, n1, n2; + cexp *ce; + int *z; + expr_v **vp; + split_ce *cs; + ASLTYPE *asl = S->asl; + + n2 = t - Ncom; + k = htcl(t*(sizeof(cexp) + sizeof(int) + sizeof(expr_v*)) + + n2*sizeof(split_ce)); + memset(ce = (cexp *)new_mblk(k), 0, n1 = sizeof(char*) << k); + n1 = (n1 + Ncom*sizeof(split_ce)) + / (sizeof(cexp) + sizeof(int) + sizeof(split_ce) + + sizeof(expr_v*)); + n2 = n1 - Ncom; + cs = (split_ce *)(ce + n1); + vp = (expr_v **)(cs + n2); + z = (int *)(vp + n1); + if (cexps) { + if (nsce) + memcpy(cs, asl->P.Split_ce, nsce*sizeof(split_ce)); + memcpy(ce, cexps, cexp_n*sizeof(cexp)); + memcpy(z, zl, cexp_n*sizeof(int)); + memcpy(vp, varp, cexp_n*sizeof(expr_v*)); + del_mblk(cexp_k, cexps); + } + nsce = n2; + asl->P.Split_ce = cs; + cexps = ce; + zl = z; + cexp_k = k; + cexp_n = n1; + varp = vp; + } + + static void +zc_upgrade(Static *S) +{ + int k, n, n0; + int *z; + ASLTYPE *asl = S->asl; + + k = htcl(sizeof(int)*(max_var1 + 1)) + 1; + z = (int*)new_mblk(k); + n = (sizeof(char*)/sizeof(int)) << (k-1); + memset(z + n, 0, n*sizeof(int)); + if (zci) { + n0 = (sizeof(char*)/sizeof(int)) << (kzc - 1); + memcpy(z, zci, n0*sizeof(int)); + memcpy(z+n, zci+n0, n0*sizeof(int)); + del_mblk(kzc, zci); + } + kzc = k; + zci = z; + zc = z + n + 1; + zc_lim = n; + } + + static void +la_replace(Static *S, linarg *la) +{ + la_ref *r; + expr_v *v; + expr_vx *vx; + expr *eout; + ASLTYPE *asl = S->asl; + + if (la->nnz > 1) { + if (!(v = la->v)) { + vx = (expr_vx *)mem(sizeof(expr_vx)); + vx->la = la; + vx->a0 = vx->a1 = -1; + la->v = v = (expr_v *)vx; + v->a = max_var + nndv++; + lasta00++; + v->op = (efunc *)f_OPVARVAL; + if (larvlist) { + *larvlist = v; + larvlist = (expr_v**)&v->v; + } + } + } + else + v = var_e + la->nz->varno; + for(r = la->refs; r; r = r->next) { + efree(S, *r->ep); + eout = (expr *)v; + if (r->scale != 1.) { + if (r->scale == -1.) { + eout = new_expr(S, f_OPUMINUS, eout, 0); + eout->dL = -1.; + } + else + eout = new_expr(S, f_OPMULT, eout, + (expr*)new_expr_n(S, r->scale)); + } + if (r->c) { + eout = new_expr(S, f_OPPLUS, eout, + (expr*)new_expr_n(S, r->c)); + eout->dL = 1.; + } + *r->ep = eout; + } + free_laref(S, &la->refs); + } + + static void +tlistgen(Static *S, ps_func *f) +{ + int *ci, *cie, i, t; + linarg *la, **lap, **lape, *tl; + ograd *og; + range *r; + psb_elem *b, *be; + + t = f->nb; + b = f->b; + be = b + t; + t = ++Termno; + tl = 0; + for(; b < be; b++) { + if ((ci = b->ce)) { + cie = ci + *ci; + do { + if (!zc[i = nv0x + *++ci]) + zci[nzc++] = i; + } + while(ci < cie); + } + r = b->U; + lap = r->lap; + lape = lap + r->n; + while(lap < lape) { + la = *lap++; + if (la->termno != t) { + la->termno = t; + la->tnext = tl; + tl = la; + for(og = la->nz; og; og = og->next) + if (!zc[og->varno]++) + zci[nzc++] = og->varno; + } + } + } + tlist = tl; + } + + static int +might_expand(Static *S, expr *e) +{ + top: + switch(Intcast e->op) { + case f_OPPLUS: + case f_OPMINUS: + case f_OPUMINUS: + case f_OPSUMLIST: + return 1; + case f_OPMULT: + if (Intcast e->R.e->op == f_OPNUM) { + e = e->L.e; + goto top; + } + if (Intcast e->L.e->op == f_OPNUM) { + e = e->R.e; + goto top; + } + break; + case f_OPVARVAL: + if (((expr_v*)e)->a >= nv0x) + return 1; + } + return 0; + } + + static void +ce_split(Static *S, int i, ps_func *f) +{ + int j, j1, je, k, n; + cexp *c, *ce; + expr **ep; + expr_v **vp, **vp0; + expr_vx *vx; + split_ce *cs; + psb_elem *b; + ASLTYPE *asl = S->asl; + + asl->P.ndvsplit++; + n = f->nb; + j1 = asl->P.ndvspout; + j = k = j1 + Ncom; + asl->P.ndvspout += n; + je = j + n; + if (je > cexp_n) + cexp_upgrade(S, je); + c = cexps + j; + b = f->b; + cs = asl->P.Split_ce + (j-Ncom); + for(ce = c + n; c < ce; c++, b++, cs++) { + c->e = b->D.e; + cs->r = b->U; + cs->ce = b->ce; + } + c = cexps + i; + vp0 = vp = varp + j; + j = max_var + nndv; + je = j + n; + nndv += n; + j1 += max_var; + while(j < je) { + vx = (expr_vx *)mem(sizeof(expr_vx)); + vx->la = 0; + *vp++ = (expr_v*)vx; + vx->v.a = vx->a0 = j++; + vx->a1 = j1++; + vx->v.op = (efunc *)f_OPVARVAL; + } + if (n == 2) + c->e = new_expr(S, f_OPPLUS, (expr*)vp0[0], (expr*)vp0[1]); + else { + ep = (expr**)new_mblk(htcl(n*sizeof(expr*))); + memcpy(ep, vp0, n*sizeof(expr*)); + c->e = new_expr(S, f_OPSUMLIST, (expr*)ep, (expr*)(ep+n)); + } + if ((max_var1 += n) >= zc_lim) + zc_upgrade(S); + + /* adjust zl for use in compar() */ + i += nv0x; + j = k + nv0x; + n += k; + do { + zl[k++] = i; + zci[nzc++] = j++; /* for crefs */ + } + while(k < n); + } + + static void +linpart_augment(Static *S, cexp *c, ograd *og, ps_func *f) +{ + linpart *L, *Le; + int n; + ograd *og1; + real t; + ASL *asl = S->a; + + if (og->varno < 0) { + if ((t = ogfree1(S, &og))) + f->b->D.e = new_expr(S, f_OPPLUS, f->b->D.e, + (expr*)new_expr_n(S, t)); + if (!og) + return; + } + if ((n = c->nlin)) { + L = c->L; + Le = L + n; + og1 = 0; + do { + --Le; + og1 = new_ograd(S, og1, Le->v.i, Le->fac); + } + while(Le > L); + del_mblk(htcl(htcl(n*sizeof(linpart))), c->L); + og = af_sum(S, og, og1); + } + og1 = og; + for(n = 0; og; og = og->next) + n++; + c->L = L = (linpart*)new_mblk(htcl((c->nlin = n)*sizeof(linpart))); + for(og = og1; og; og = og->next, L++) { + L->v.i = og->varno; + L->fac = og->coef; + } + ogfree(S, og1); + } + + static ograd * +linterms(Static *S, cexp *c, real scale) +{ + linpart *L = c->L; + linpart *Le = L + c->nlin; + ograd *og = 0; + while(Le > L) { + --Le; + og = new_ograd(S, og, Le->v.i, scale*Le->fac); + } + return og; + } + + static void +dvwalk(Static *S, int i) +{ + cexp *c; + dv_info *dvi; + expr *e; + int n, split; + linarg *tl, **tp; + ograd *og, *og0; + ps_func f; + ASLTYPE *asl = S->asl; + + Termno++; + tlist = 0; + c = cexps + i; + split = zl[i] & 2; + zl[i] = 0; + if (split && might_expand(S, c->e)) { + asl->P.ndvspcand++; + if ((og = cotermwalk(S, &c->e, &f, 0, 0)) && f.nb >= 1) + linpart_augment(S, c, og, &f); + if (f.nb > 1) { + zl[i] = f.nb; + ce_split(S, i, &f); + c = cexps + i; + og = 0; + } + else if (f.nb == 1) { + c->e = f.b->D.e; + og = 0; + } + tlistgen(S, &f); + del_Elemtemp(S, last_psb_elem); + } + else if ((e = c->e)) + og = awalk(S, e); + else + og = 0; + og0 = og; + if (c->nlin) + og = af_sum(S, og, linterms(S,c,1.)); + asl->P.dvsp0[i+1] = asl->P.ndvspout + Ncom; + dvi = asl->P.dv + i; + if (og0) { + c->cref = crefs(S); + dvi->ll = og; + dvi->nl = 0; + return; + } + if ((c->la = dvi->lt = afree(S, og, 0))) + dvi->scale = lt_scale; + for(n = 1, tl = tlist; tl; tl = tl->tnext) + n++; + dvi->ll = 0; + dvi->nl = tp = (linarg **)mem(n*sizeof(linarg*)); + for(tl = tlist; tl; tl = tl->tnext) + la_replace(S, *tp++ = tl); + *tp = 0; + c->cref = crefs(S); + } + + static int +colindvref(Static *S, expr *e, int ndv) +{ + expr **a, **ae; + int j, k, rv = 0; + + top: + switch(Intcast e->op) { + case f_OPPLUS: + case f_OPMINUS: + rv |= colindvref(S, e->R.e, ndv); + /* no break */ + case f_OPUMINUS: + e = e->L.e; + goto top; + case f_OPSUMLIST: + a = e->L.ep; + ae = e->R.ep; + while(a < ae) + rv |= colindvref(S, *a++, ndv); + break; + case f_OPVARVAL: + k = ((expr_v *)e)->a; + if ((k -= nv0x) < 0) + break; + if (zl[k]) { + rv |= zl[k]; + break; + } + zl[k] = 1; + if ((j = colindvref(S, (S->cexps + k)->e, k))) { + rv |= j; + zl[k] |= j; + } + break; + case f_OPMULT: + if (Intcast e->R.e->op == f_OPNUM) { + e = e->L.e; + goto top; + } + if (Intcast e->L.e->op == f_OPNUM) { + e = e->R.e; + goto top; + } + /* no break */ + default: + if (ndv >= 0) + rv = zl[ndv] |= 2; + } + return rv; + } + + static void +dv_walk(Static *S) +{ + int i, m, n; + expr_v *v; + ASLTYPE *asl = S->asl; + + m = asl->i.n_con0; + if ((n = Ncom)) { + for(i = 0; i < n_obj; i++) + colindvref(S, obj_de[i].e, -1); + for(i = 0; i < m; i++) + colindvref(S, con_de[i].e, -1); + larvlist = &v; + for(i = 0; i < n; i++) + dvwalk(S, i); + *larvlist = 0; + if ((i = asl->P.ndvspout)) + nv0b = nv0x + Ncom + i; + larvlist = 0; + } + } + + static int /* adjust zci, return k s.t. zci[i] < nv0 for i < k */ +nzcperm(Static *S) +{ + int i, j, k; + + k = nzc; + for(i = 0; i < k; ) + if (zci[i] >= nv0x) { + j = zci[--k]; + zci[k] = zci[i]; + zci[i] = j; + } + else + i++; + return k; + } + + static void +sumlist_adj(ASLTYPE *asl, expr *e, expr *e1) +{ + int k, n; + expr **ep, **ep0, **ep1; + + ep = e->R.ep; + k = htcl((n = ep - e->L.ep)*sizeof(expr*)); + if (n == (1 << k)) { + ep1 = (expr**)new_mblk(k+1); + memcpy(ep1, ep0 = e->L.ep, n*sizeof(expr*)); + del_mblk(k, ep0); + e->L.ep = ep1; + ep = ep1 + n; + } + *ep++ = e1; + e->R.ep = ep; + } + + static void +termwalk(Static *S, expr **ep, PSfind *p) +{ + ograd *og; + linarg **lap, **lap1; + linarg *tl; + int *cp, *cp0, *ce, *cee, i, j, k, kl, n, ncp, nzc2, *ui; + size_t len; + range *r; + list *L; + expr *e, *e1, *e2; + split_ce *cs; + psb_elem *b; + Elemtemp *bt; + ps_func *f; + ASLTYPE *asl = S->asl; + + Termno++; + tlist = 0; + afree(S, awalk(S, *ep), ep); + + i = k = nzcperm(S); + f = p->f; + if (!larvlist) + for(; i < nzc; i++) { + if ((j = zci[i]) < max_var) { + for(L = cexps[j-nv0x].cref; L; L = L->next) + if (!zc[L->item.i]++) + zci[nzc++] = L->item.i; + } + else { + cs = asl->P.Split_ce + (j-max_var); + if ((ce = cs->ce)) { + cee = ce + *ce; + while(ce++ < cee) + if (!zc[j = *ce + nv0x]++) + zci[nzc++] = j; + } + } + } + + r = (range *)rnz; /* scratch */ + if ((ncp = nzc - k)) { + zcsort(S, zc, zci+k, nv0x, ncp, -1 /*max_var*/); + cp = cp0 = (int*)(r+1); + i = k; + *cp = ncp; + do { + zc[j = zci[i++]] = 0; + *++cp = j - nv0x; + } while(i < nzc); + } + else + cp0 = 0; + + for(n = 0, tl = tlist; tl; tl = tl->tnext) { + n++; + for(og = tl->nz; og; og = og->next) + if (!zc[og->varno]++) + zci[k++] = og->varno; + } + nzc2 = k; + if (zc[-1]) + --nzc2; /* ignore constant */ + if (n <= 0 && Intcast (*ep)->op == f_OPNUM) + goto done; + + r->n = n; + r->nv = nzc2; + r->lap = lap1 = lap = (linarg**) + new_mblk(kl = htcl(n*sizeof(linarg*))); + for(tl = tlist; tl; tl = tl->tnext) + la_replace(S, *lap1++ = tl); + if (n > 1) + qsortv(lap, n, sizeof(linarg*), lacompar, NULL); + r->ui = ui = cp0 ? cp+1 : (int*)(r+1); + zcsort(S, zc, zci, 0, nzc2, nv0x); + for(j = 0; j < nzc2; j++) + *ui++ = zci[j]; + r = n >= nzc2 ? uhash(S,r) : rhash(S,r,1); + del_mblk(kl, lap); + + e1 = *ep; + if (!r || (i = r->lasttermno) == -1 || r->lastgroupno != Groupno) { + bt = p->b; + if ((i = f->nb++) >= bt->nmax) + upgrade_Elemtemp(S, bt); + b = f->b + i; + b->conno = Conno; + b->termno = i; + b->groupno = Groupno; + if ((b->U = r)) { + r->lasttermno = i; + r->lastgroupno = Groupno; + } + b->D.e = e1; + if (cp0) { + cp = (int*)mem(len = (ncp+1)*sizeof(int)); + memcpy(cp, cp0, len); + cp0 = cp; + } + b->ce = cp0; + while(k > 0) + zc[zci[--k]]= 0; + } + else { + b = f->b + i; + while(k > 0) + zc[zci[--k]] = 0; + if ((cp = b->ce) + && (*cp != ncp || memcmp(cp+1, cp0+1, ncp*sizeof(int)))) { + /* fix up f->ce */ + n = *cp; + while(k < n) + zc[zci[k++] = *++cp] = 1; + for(n = 0; n < ncp; ) + if (!zc[j = cp0[++n]]++) + zci[k++] = j; + qsortv(zci, k, sizeof(int), hscompar, S); + b->ce = cp = (int*)mem((k+1)*sizeof(int)); + *cp++ = k; + n = 0; + while(n < k) + zc[*cp++ = zci[n++]] = 0; + } + e = b->D.e; + switch(Intcast e->op) { + case f_OPPLUS: + ep = (expr**)new_mblk(2); + ep[0] = e->L.e; + ep[1] = e->R.e; + ep[2] = e1; + e->L.ep = ep; + e->R.ep = ep + 3; + e->op = (efunc *)f_OPSUMLIST; + break; + + case f_OPSUMLIST: + sumlist_adj(asl, e, e1); + break; + default: + b->D.e = e2 = new_expr(S, 0, e, e1); + e2->op = (efunc *)f_OPPLUS; + } + } + done: + nzc = 0; + } + + + static void +co_finish(ps_func *f) +{ + range *r; + psb_elem *b, *be; + psg_elem *g, *ge; + + b = f->b; + for(be = b + f->nb; b < be; b++) + if ((r = b->U)) + r->lasttermno = -1; + g = f->g; + for(ge = g + f->ng; g < ge; g++) + for(b = g->E, be = b + g->ns; b < be; b++) + if ((r = b->U)) + r->lasttermno = -1; + } + + static expr * +ecopy(Static *S, expr *e) +{ + expr **a, **a1, **ae; + int n, op; + + switch(op = Intcast e->op) { + case f_OPPLUS: + case f_OPMINUS: + e = new_expr(S, op, ecopy(S, e->L.e), ecopy(S, e->R.e)); + break; + case f_OPMULT: + if (Intcast e->L.e->op == f_OPNUM) { + e = new_expr(S, op, ecopy(S, e->R.e), (expr*) + new_expr_n(S, ((expr_n*)e->L.e)->v)); + break; + } + assert(Intcast e->R.e->op == f_OPNUM); + e = new_expr(S, op, ecopy(S, e->L.e), + (expr*)new_expr_n(S, ((expr_n*)e->R.e)->v)); + break; + case f_OPUMINUS: + e = new_expr(S, op, ecopy(S, e->L.e), 0); + break; + case f_OPSUMLIST: + a = e->L.ep; + ae = e->R.ep; + a1 = (expr**)new_mblk_ASL(S->a, htcl((n = ae-a)*sizeof(expr*))); + e = new_expr(S, op, (expr*)a1, (expr*)(a1+n)); + while(a < ae) + *a1++ = ecopy(S, *a++); + break; + case f_OPVARVAL: + break; +#ifdef DEBUG + default: + fprintf(Stderr, "Impossible case in ecopy!\n"); + exit(1); +#endif + } + return e; + } + + static int getgroup ANSI((Static *S, real scale, expr *e, PSfind *p)); + + static ograd * +ltermwalk(Static *S, real scale, expr **ep, PSfind *p) +{ + ASLTYPE *asl; + EU *en; + dv_info *dvi; + expr **a, **ae, *e; + int k, k0; + ograd *og, *rv; + cexp *ce; + + asl = S->asl; + rv = 0; + top: + e = *ep; + switch(Intcast e->op) { + case f_OPPLUS: + rv = af_sum(S, rv, ltermwalk(S, scale, &e->L.e, p)); + ep = &e->R.e; + goto top; + case f_OPMINUS: + rv = af_sum(S, rv, ltermwalk(S, scale, &e->L.e, p)); + ep = &e->R.e; + scale = -scale; + goto top; + case f_OPUMINUS: + ep = &e->L.e; + scale = -scale; + goto top; + case f_OPSUMLIST: + a = e->L.ep; + ae = e->R.ep; + while(a < ae) + rv = af_sum(S, rv, ltermwalk(S, scale, a++, p)); + break; + case f_OPVARVAL: + k = e->a; + if (k < nv0x) { + rv = af_sum(S, rv, new_ograd(S,0,k,scale)); + break; + } + k0 = k; + if ((k -= nv0x) >= Ncom) + goto dflt; + if (zl[k]) { + ce = cexps + k; + *ep = ecopy(S, ce->e); + if (ce->nlin) + rv = af_sum(S, rv, linterms(S,ce,scale)); + goto top; + } + if (!zc[e->a]++) + zci[nzc++] = e->a; + dvi = asl->P.dv + k; + if (dvi->nl) + goto dflt; + og = new_ograd(S, 0, k0, scale); + rv = af_sum(S, rv, og); + break; + case f_OPNUM: + rv = af_sum(S, rv, + new_ograd(S, 0, -1, scale*((expr_n *)e)->v)); + break; + case f_OPMULT: + en = (EU *)e->R.e; + if (Intcast en->op == f_OPNUM) { + *ep = e->L.e; + case_opnum: + if (en->u.v == 0.) { + efree(S,*ep); + *ep = (expr*)en; + } + else { + rv = af_sum(S, rv, ltermwalk(S, + scale*en->u.v, ep, p)); + en->u.p = expr_n_free; + expr_n_free = en; + } + e->L.e = expr_free; + expr_free = e; + break; + } + en = (EU *)e->L.e; + if (Intcast en->op == f_OPNUM) { + *ep = e->R.e; + goto case_opnum; + } + default: + dflt: + if (p->g && getgroup(S, scale, e, p)) + break; + if (scale != 1.) + *ep = scale == -1. + ? new_expr(S, f_OPUMINUS, e, 0) + : new_expr(S, f_OPMULT, e, + (expr *)new_expr_n(S, scale)); + termwalk(S, ep, p); + asl->P.ns0++; + } + return rv; + } + + static void +PSfind_init(Static *S, ps_func *f, PSfind *psf, int wantg) +{ + f->nxval = f->nb = f->ng = 0; + psf->f = f; + psf->b = new_Elemtemp(S, sizeof(psb_elem), (void**)&f->b); + if (wantg) + psf->g = new_Elemtemp(S, sizeof(psg_elem), (void**)&f->g); + else { + psf->g = 0; + f->g = 0; + last_psb_elem = psf->b; + } + } + + static void +psb_copy(psb_elem *b, psb_elem *b0, int n) +{ + range *r; + psb_elem *be; + + memcpy(b, b0, n*sizeof(psb_elem)); + for(be = b + n; b < be; b++) { + /* *b = *b0++; */ /* DEC Alpha gcc got this wrong */ + if (b->conno != -1 && (r = b->U)) { + b->next = r->refs; + r->refs = b; + } + } + } + + static int +getgroup(Static *S, real scale, expr *e, PSfind *p) +{ + ps_func *f, f1; + PSfind p1; + expr *e0, *e1, *e2; + ograd *og, *og1; + psb_elem *b, *be; + psg_elem *g; + Elemtemp *gt; + linarg **lap, **lape; + linpart *L; + range *U; + int i, nzc1, *zc1, *zci1; + ASL *asl = S->a; + + for(e1 = e; optype[Intcast e1->op] == 1; e1 = e1->L.e) + e0 = e1; + if (e1 == e || !might_expand(S, e1)) + return 0; + PSfind_init(S, &f1, &p1, 0); + f = p->f; + gt = p->g; + if ((i = f->ng++) >= gt->nmax) + upgrade_Elemtemp(S, gt); + Groupno = f->ng; + memset(g = f->g + i, 0, sizeof(psg_elem)); + g->g = e; + g->ge = e0; + g->scale = scale; + if ((og = ltermwalk(S, 1., &e0->L.e, &p1))) + og = compress(S, og, &g->g0, &i); + for(e1 = e; e1 != e0; e1 = e2) { + e2 = e1->L.e; + e2->R.e = e1; /* back pointer, used in psderprop */ + } + zc1 = S->zc1; + zci1 = S->zci1; + Groupno = nzc1 = 0; + if ((og1 = og)) { + for(i = 1; (og = og->next); i++); + g->nlin = i; + g->L = L = (linpart*)mem(i*sizeof(linpart)); + for(og = og1;; L++) { + zc1[zci1[nzc1++] = L->v.i = og->varno]= 1; + L->fac = og->coef; + if (!(og = og->next)) + break; + } + ogfree(S, og1); + } + g->esum.op = (efunc_n *)f_OPNUM; + g->ns = f1.nb; + g->E = b = (psb_elem*)mem(f1.nb * sizeof(psb_elem)); + psb_copy(b, f1.b, f1.nb); + del_Elemtemp(S, p1.b); + for(be = b + f1.nb; b < be; b++) { + if (!(U = b->U)) + continue; + lap = U->lap; + lape = lap + U->n; + while(lap < lape) + for(og = (*lap++)->nz; og; og = og->next) + if (!zc1[og->varno]++) + zci1[nzc1++] = og->varno; + } + zcsort(S, zc1, zci1, 0, nzc1, nv0x); + og = 0; + while(nzc1 > 0) { + og = new_ograd(S, og, i = zci1[--nzc1], 0.); + zc1[i] = 0; + } + g->og = og; + return 1; + } + + static ograd * +cotermwalk(Static *S, expr **ep, ps_func *f, int wantg, int omitdv) +{ + int comvar, n, x; + ograd *og; + real t; + PSfind psf; + psb_elem *b; + psg_elem *g, *g1, *ge; + + PSfind_init(S, f, &psf, wantg); + t = 0.; + og = ltermwalk(S, 1., ep, &psf); + if (omitdv && og) + og = compress(S, og, &t, &comvar); + b = f->b; + if (f->nb + f->ng == 0) + *ep = (expr*)new_expr_n(S, t); + else if (t) { + if (f->nb) + b->D.e = new_expr(S, f_OPPLUS, b->D.e, + (expr*)new_expr_n(S, t)); + else { + f->nb = 1; + memset(b, 0, sizeof(psb_elem)); + b->D.e = (expr *)new_expr_n(S, t); + } + } + co_finish(f); + if (!larvlist) { + n = f->nb; + if ((x = n * sizeof(psb_elem) + f->ng * sizeof(psg_elem))) { + ASL *asl = S->a; + g = (psg_elem*)(x >= 256 ? M1alloc(x) : mem(x)); + b = (psb_elem*)(g + f->ng); + if (n) + psb_copy(b, f->b, n); + else + b = 0; + if ((n = f->ng)) { + memcpy(g1 = g, f->g, n*sizeof(psg_elem)); + for(ge = g + n; g1 < ge; g1++) + g1->ge->L.e = (expr*)&g1->esum; + } + else + g = 0; + del_Elemtemp(S, psf.b); + if (wantg) + del_Elemtemp(S, psf.g); + f->b = b; + f->g = g; + } + } + return og; + } + + static void +psfind(Static *S) +{ + ASLTYPE *asl; + fint x; + int i, j, k, m, mx, nv, nx; + linarg *la, **lap; + ps_func *f, *f1; + range *r, *r0; + size_t L; +#ifdef PSHVREAD + Ihinfo *ihi, *ihi0, *ihi1; + int ihdlim, ihdlim1, ihdmax, ihdmin, n; +#endif + asl = S->asl; + m = asl->i.n_con0; + mx = m + (nx = asl->i.nsufext[ASL_Sufkind_con]); + x = (n_obj+mx)*sizeof(ps_func) + + slmax*sizeof(expr*) + + (Ncom+1)*sizeof(int); + asl->P.ops = f = (ps_func *)M1alloc(x); + asl->P.cps = f1 = f + n_obj; + if (nx) + memset(f1 + m, 0, nx*sizeof(ps_func)); + slscratch = (expr**)(f1 + mx); + asl->P.dvsp0 = (int*)(slscratch + slmax); + *asl->P.dvsp0 = Ncom; + + k = 3; + if ((nv = n_var) > 7) { + nv >>= 3; + do ++k; while(nv >>= 1); + } + S->klthash = S->krangehash = k; + lthash = (linarg**)new_mblkzap(asl, k); + rangehash = (range**)new_mblkzap(asl, k); + L = 1; + lthashmask = rangehashmask = (L << k) - 1; + nrange = 0; + + /* Separate allocation of zc1, zci1, rnz to allow freeing them later. */ + /* Allow room in rnz for use as scratch in termwalk. */ + i = max_var + (sizeof(range)+sizeof(int)+sizeof(real)-1)/sizeof(real); + rnz = (real *)M1alloc(i*(sizeof(real)+2*sizeof(int))); + S->zc1 = (int *)(rnz + i); + S->zci1 = S->zc1 + i; + memset(S->zc1, 0, max_var*sizeof(int)); + + Conno = -1; /* do not record refs for split defined vars */ + dv_walk(S); + asl->P.ndvspin = asl->P.ns0; + asl->P.ns0 = 0; + PSHV(n = 0); + for(i = 0; i < n_obj; i++, f++) { + Conno = -2 - i; + Ograd[i] = af_sum(S, Ograd[i], cotermwalk(S, &obj_de[i].e, f, + wantOgroups, 1)); + PSHV(n += f->ng); + } + PSHV(asl->P.nobjgroups = n); + PSHV(n = 0); + for(i = 0; i < m; i++, f1++) { + Conno = i; + Cgrad[i] = cf_sum(S, Cgrad[i], cotermwalk(S, &con_de[i].e, f1, + wantCgroups, 1)); + PSHV(n += f1->ng); + } + PSHV(asl->P.ncongroups = n); + + /* Compute nintv value for ranges */ + /* and set lasttermno to first var. */ + + r0 = (range*)&asl->P.rlist; + for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) { + i = 0; + if ((j = r->n) > 0) { + lap = r->lap; + r->lasttermno = (*lap)->nz->varno; + while(i < j) { + la = lap[i]; + if (la->v) { + i++; + la->termno = 0; /* for psgcomp */ + } + else { + lap[i] = lap[--j]; + lap[j] = la; + } + } + } + r->nintv = i; + } +#ifdef PSHVREAD + if ((ihdlim = ihd_limit) > 0) { + ihdmin = ihdlim1 = ihdlim + 1; + n = ihdlim1*sizeof(Ihinfo); + asl->P.ihi = ihi0 = (Ihinfo*)M1zapalloc(n); + ihdmax = 0; + for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) + if ((n = r->n) > 0) { + if (n > r->nv) + n = r->nv; + if (n > ihdlim) + n = ihdlim1; + else { + if (ihdmax < n) + ihdmax = n; + if (ihdmin > n) + ihdmin = n; + } + ihi = ihi0 + n - 1; + r->rlist.prev = ihi->r; + ihi->r = r; + ihi->nr++; + } + asl->P.ihdmax = ihdmax; + asl->P.ihdmin = asl->P.ndhmax = ihdmin; + ihi1 = ihi = ihi0 + ihdlim; + ihi->ihd = ihdlim1; /* sentinel */ + for(i = ihdlim; i > 0; --i) + if ((n = (--ihi)->nr)) { + ihi->next = ihi1; + ihi1 = ihi; + ihi->ihd = i; + ihi->k = htcl(((i*(i+1))>>1)*n*sizeof(real)); + } + asl->P.ihi1 = ihi1; + } +#endif + del_mblk(S->krangehash, rangehash); + del_mblk(S->klthash, lthash); + lthash = 0; + rangehash = 0; + } + +/*******/ + + static void +ewalk(Static *S, expr *e, int deriv) +{ + int a0, a1, i, j, j1, k, kf, numargs, op; + real *b, *b0, *ra; + unsigned int len; +#ifdef PSHVREAD + ASL *asl; + argpair *da; + int i1, j0, k0, k1, kn, *nn, *nn0; + real **fh; + expr **args1; +#endif + expr *L, *arg, **args, **argse; + expr_va *rva; + de *d; + derp *dsave; + expr_if *rvif; + expr_f *rvf; + arglist *al; + argpair *ap, *ap0, *ape; + char *dig; + + switch(optypeb[k = Intcast e->op]) { + + case 1: /* unary */ + e->dL = dvalue[k]; /* for UMINUS, FLOOR, CEIL */ + ewalk(S, e->L.e, deriv); + if (deriv) + dexpr(S, e, e->L.e, 0); + break; + + case 2: /* binary */ + e->dL = 1.; + e->dR = dvalue[k]; /* for PLUS, MINUS, REM */ + if (dvalue[k] == 11) + return; + ewalk(S, e->L.e, deriv); + ewalk(S, e->R.e, deriv); + if (deriv) + dexpr(S, e, e->L.e, e->R.e); + break; + + case 3: /* vararg (min, max) */ + rva = (expr_va *)e; + rva->next2 = varg2list; + varg2list = rva; + if (!last_d) { + new_derp(S, lasta, lasta, &edagread_one); + lasta++; + } + rva->d0 = dsave = last_d; +#ifdef PSHVREAD + rva->bak = last_e; +#endif + a0 = a1 = lasta; + j = 0; + for(d = rva->L.d; (L = d->e); d++) { + last_d = dsave; + op = Intcast L->op; + PSHV(last_e = 0;) + ewalk(S, L, deriv); + PSHV(d->ee = last_e;) + if (op == f_OPNUM || L->a == noa) { + d->d = dsave; + d->dv.i = noa; + } + else if (deriv) { + d->d = new_relo(S, L, dsave, &d->dv.i); + j++; + if (a1 < lasta) + a1 = lasta; + lasta = a0; + } + } + last_d = dsave; + if (j) { + rva->a = lasta = a1; + new_derp(S, 0, lasta++, &edagread_one); + /* f_MINLIST or f_MAXLIST will replace the 0 */ + rva->R.D = last_d; + nocopy = 1; +#ifdef PSHVREAD + last_e = (expr *)rva; + rva->dO.i = Hv_vararg; +#endif + } + else { + rva->a = noa; + rva->R.D = 0; + PSHV(last_e = rva->bak); + } + break; + + case 4: /* piece-wise linear */ + ewalk(S, e->R.e, deriv); /* should be a no-op */ + if (deriv) { + new_derp(S, e->R.e->a, e->a = lasta++, &e->dL); +#ifdef PSHVREAD + e->bak = last_e; + last_e = e; + e->dO.i = Hv_plterm; +#endif + } + break; + + case 5: /* if */ + rvif = (expr_if *)e; + rvif->next2 = if2list; + if2list = rvif; + PSHV(rvif->bak = last_e;) + ewalk(S, rvif->e, 0); + if (deriv && !last_d) { + new_derp(S, lasta, lasta, &edagread_one); + lasta++; + } + rvif->d0 = dsave = last_d; + a0 = lasta; + L = rvif->T; + op = Intcast L->op; + PSHV(last_e = 0;) + ewalk(S, L, deriv); + PSHV(rvif->Te = last_e;) + j = 0; + if (op == f_OPNUM || L->a == noa) { + rvif->dT = dsave; + rvif->Tv.i = noa; + } + else if ((j = deriv)) + rvif->dT = new_relo(S, L, dsave, &rvif->Tv.i); + a1 = lasta; + lasta = a0; + last_d = dsave; + L = rvif->F; + op = Intcast L->op; + PSHV(last_e = 0;) + ewalk(S, L, deriv); + PSHV(rvif->Fe = last_e;) + if (op == f_OPNUM || L->a == noa) { + rvif->dF = dsave; + rvif->Fv.i = noa; + } + else if ((j = deriv)) + rvif->dF = new_relo(S, L, dsave, &rvif->Fv.i); + if (lasta < a1) + lasta = a1; + last_d = dsave; + if (j) { + new_derp(S, 0, rvif->a = lasta++, &edagread_one); + rvif->D = last_d; + nocopy = 1; +#ifdef PSHVREAD + last_e = (expr *)rvif; + rvif->dO.i = Hv_if; +#endif + } + else { + rvif->a = noa; + rvif->D = 0; + PSHV(last_e = rvif->bak); + } + break; + + case 6: /* sumlist */ + a0 = lasta; + j = 0; + args = e->L.ep; + argse = e->R.ep; + e->a = lasta++; + while(args < argse) { + L = *args++; + op = Intcast L->op; + ewalk(S, L, deriv); + if (op != f_OPNUM && L->a != noa && deriv) { + new_derp(S, L->a, e->a, &edagread_one); + j++; + } + } + if (!j) { + e->a = noa; + lasta = a0; + } +#ifdef PSHVREAD + else { + e->bak = last_e; + last_e = e; + e->dO.i = Hv_sumlist; + j0 = e->R.ep - e->L.ep; + k0 = htcl(j0*sizeof(expr*)); + j1 = j0 + j + 1; + k1 = htcl(j1*sizeof(expr*)); + if (k1 > k0) { + ASL *asl = S->a; + args1 = (expr**)new_mblk(k1); + memcpy(args1, e->L.ep, + j0*sizeof(expr*)); + del_mblk(k0, e->L.ep); + e->L.ep = args1; + args = e->R.ep = args1 + j0; + } + argse = e->R.ep; + args1 = e->L.ep; + while(args1 < args) { + arg = *args1++; + if (arg->op != OPNUM && arg->a != noa) + *argse++ = arg; + } + *argse = 0; + } +#endif + break; + + case 7: /* function call */ + rvf = (expr_f *)e; + ap = rvf->sap; + for(ape = rvf->sape; ap < ape; ap++) + ewalk(S, ap->e,0); + ap = rvf->ap; + ape = rvf->ape; + kf = ape - ap; + al = rvf->al; + ra = al->ra; + numargs = al->nr; + for(i = 0; ap < ape; ap++) { + arg = ap->e; + op = Intcast arg->op; + ewalk(S, arg, deriv); + arg->op = (efunc *)(size_t)op; + if (arg->a != noa) + i += deriv; + } + rvf->a = i ? lasta++ : noa; + if (!numargs) { +#ifdef PSHVREAD + rvf->da = rvf->dae = 0; + rvf->fh= 0; /* for debugging */ +#endif + break; + } + dig = 0; + if (i) { + b = (real *)mem_ASL(S->a, len = +#ifdef PSHVREAD + (numargs*(numargs+3) >> 1) +#else + numargs +#endif + *sizeof(real)); + memset(b, 0, len); +#ifdef PSHVREAD + al->hes = b + numargs; +#endif + if (kf < numargs) { + ASLTYPE *asl = S->asl; + dig = (char*)mem(numargs); + for(j1 = 0; j1 < numargs; ++j1) + dig[j1] = 1; + } + } + else + b = 0; + al->derivs = b; + al->dig = dig; + b0 = b; +#ifdef PSHVREAD + asl = S->a; + if (i) { + da = (argpair*)mem(kf*sizeof(argpair) + + kf*kf*sizeof(real*)); + fh = (real**)(da + kf); + kn = htcl(2*sizeof(int)*numargs); + nn = (int*)new_mblk(kn); + } + else { + da = 0; + fh = 0; + kn = 12345; /* silence false "used uninitialized" warning */ + nn = 0; + } + rvf->da = da; + rvf->fh = fh; + nn0 = nn; +#endif + + for(ap = ap0 = rvf->ap; ap < ape; ap++) { + j1 = 0; + arg = ap->e; + arg->op = r_ops[op = Intcast arg->op]; + switch(op) { + case f_OPNUM: + goto loopend; + case f_OPVARVAL: + arg->op = (efunc *)(size_t)op; + } + b = b0 + (j = ap->u.v - ra); + *b = 0; + if (arg->a == noa) + goto loopend; + new_derp(S, arg->a, rvf->a, b); +#ifdef PSHVREAD + da->e = arg; + (da++)->u.v = b; + *nn++ = j; + j1 = 1; +#endif + loopend: + if (dig && j1) + dig[j] = 0; + } +#ifdef PSHVREAD + rvf->dae = da; + if (i) { + rvf->bak = last_e; + last_e = (expr *)rvf; + rvf->dO.i = Hv_func; + b = al->hes; + for(i1 = 0; i1 < kf; i1++) { + i = nn0[i1]; + for(j1 = 0; j1 < kf; j1++) { + j = nn0[j1]; + *fh++ = &b[i >= j + ? (i*(i+1)>>1)+j + : (j*(j+1)>>1)+i]; + } + } + } + if (nn0) + del_mblk(kn, nn0); +#endif + /* no break */ + + case 8: /* OPHOL (Hollerith) */ + case 9: /* OPNUM */ + break; + + case 10:/* OPVARVAL */ + k = (expr_v *)e - S->var_e; + if (k < 0 || k >= max_var) { + if ((k = ((expr_vx*)e)->a1) < 0) + return; + } + if (k >= 0 && deriv && !zc[k]++) + zci[nzc++] = k; + return; + + case 11: /* OPCOUNT */ + args = e->L.ep; + argse = e->R.ep; + while(args < argse) + ewalk(S, *args++, 0); + break; + + default: + scream(S->R, ASL_readerr_bug, + "unexpected optype[%d] = %d\n", k, optype[k]); + } +#ifdef DEBUG + if (e == opzork) + printf(""); +#endif + e->op = r_ops[k]; + } + + static list * +crefs(Static *S) +{ + int Nv0, Nzc, i, maxvar1 = S->max_var1; + list *rv = 0; + + Nv0 = nv0x; + for(Nzc = nzc; Nzc > 0; ) { + if ((i = zci[--Nzc]) >= Nv0 && i < maxvar1) { + rv = new_list(S->a, rv); + rv->item.i = i; + } + zc[i] = 0; + } + nzc = Nzc; + return rv; + } + + static funnel * +funnelfix(funnel *f) +{ + cexp *ce; + funnel *fnext, *fprev; + derp *d; + + for(fprev = 0; f; f = fnext) { + fnext = f->next; + f->next = fprev; + fprev = f; + ce = f->ce; + if ((d = ce->d)) + ce->z.i = d->b.i; + } + return fprev; + } + + static derp * +derpadjust(Static *S, derp *d0, int a, derp *dnext) +{ + derp *d, *d1; + int *r, *re; + relo *rl; + expr_if *il, *ile; + expr_va *vl, *vle; + de *de1; + ASL *asl; + + if (!(d = d0)) + return dnext; + asl = S->a; + r = imap + lasta0; + re = imap + lasta; + while(r < re) + *r++ = a++; + if (amax < a) + amax = a; + r = imap; + for(;; d = d1) { + d->a.i = r[d->a.i]; + d->b.i = r[d->b.i]; + if (!(d1 = d->next)) + break; + } + d->next = dnext; + if ((rl = relo2list)) { + relo2list = 0; + do { + d = rl->Dcond; + do { + d->a.i = r[d->a.i]; + d->b.i = r[d->b.i]; + } + while((d = d->next)); + } + while((rl = rl->next2)); + } + if (if2list != if2list_end) { + ile = if2list_end; + if2list_end = il = if2list; + do { + il->Tv.i = r[il->Tv.i]; + il->Fv.i = r[il->Fv.i]; + } + while((il = il->next2) != ile); + } + if (varg2list != varg2list_end) { + vle = varg2list_end; + varg2list_end = vl = varg2list; + do { + for(de1 = vl->L.d; de1->e; de1++) + de1->dv.i = r[de1->dv.i]; + } + while((vl = vl->next2) != vle); + } + return d0; + } + + static derp * +derpcopy(Static *S, cexp *ce, derp *dnext) +{ + derp *d, *dprev; + int *map; + derp d00; + + if (!(d = ce->d)) + return dnext; + map = imap; + for(dprev = &d00; d; d = d->next) { + new_derp(S, map[d->a.i], map[d->b.i], d->c.rp); + dprev = dprev->next = last_d; + } + dprev->next = dnext; + return d00.next; + } + + static derp * +derp_ogcopy(Static *S, ograd *og, derp *dnext, int k) +{ + last_d = dnext; + if (og->varno < 0) + og = og->next; + while(og) { + new_derp(S, og->varno, k, &og->coef); + og = og->next; + } + return last_d; + } + + static void +imap_alloc(Static *S) +{ + expr_v *v; + int i, *r; + linarg *tl; + ASLTYPE *asl = S->asl; + + if (imap) { + r = (int*)new_mblk(i = htcl(lasta*sizeof(int))); + memcpy(r, imap, imap_len*sizeof(int)); + del_mblk(kimap, imap); + imap = r; + kimap = i; + imap_len = (sizeof(char*)/sizeof(int)) << i; + return; + } + i = amax1 > lasta ? amax1 : lasta; + r = imap = (int*)new_mblk(kimap = htcl((i+100)*sizeof(int))); + imap_len = (sizeof(char*)/sizeof(int)) << kimap; + r += i = max_var1; + while(i > 0) + *--r = --i; + i = nv0x; + for(tl = asl->P.lalist; tl; tl = tl->lnext) + if ((v = tl->v)) + r[v->a] = i++; + r[noa] = i; + } + + static void +comsubs(Static *S, int alen, cde *d) +{ + list *L; + int a, i, j, jx, k, nzc1; + int *r, *re; + cexp *ce; + derp *D, *dnext; + dv_info *dvi; + relo *R; + expr *e; + split_ce *cs; + ograd *og; + ASLTYPE *asl = S->asl; + + D = last_d; + a = lasta00; + dnext = 0; + R = 0; + nzc1 = nzc; + nzc = 0; + for(i = j = 0; i < nzc1; i++) + if ((k = zci[i]) >= nv0x && k < max_var1) + zci[j++] = k; + else + zc[k] = 0; + if ((nzc1 = j)) { + for(i = 0; i < nzc1; i++) { + if ((j = zci[i] - nv0x) >= Ncom) { + cs = asl->P.Split_ce + (j-Ncom); + if ((r = cs->ce)) { + re = r + *r; + while(r++ < re) + if (!zc[j = *r + nv0x]++) + zci[nzc1++] = j; + } + } + else if (j >= 0 && (L = cexps[j].cref)) { + if (cexps[j].funneled) do { + if (zc[j = L->item.i] + || asl->P.dv[j-nv0x].ll) + continue; + zc[j] = 1; + zci[nzc1++] = j; + } + while((L = L->next)); + else do { + if (!zc[L->item.i]++) + zci[nzc1++] = L->item.i; + } + while((L = L->next)); + } + } + if (nzc1 > 1) + zcsort(S, zc, zci, 0, nzc1, -1); + R = new_relo1(S, dnext); + i = 0; + do { + j = zci[i]; + jx = j - nv0x; + zc[j] = 0; + ce = &cexps[jx]; + e = ce->e; + k = varp[jx]->a; + if (ce->funneled) { + if (j >= max_var) + imap[((expr_vx*)varp[jx])->a0] = a; + imap[k] = a++; + } + else { + r = imap + ce->z.i; + re = r + ce->zlen; + while(r < re) + *r++ = a++; + /* zlen == 1 if e->op == OPNUM */ + imap[k] = e->op == OPNUM ? a-1 : imap[e->a]; + if (!ce->d + && jx < Ncom + && (og = (dvi = &asl->P.dv[jx])->ll) + && (dvi->nl || dvi->lt)) { + dnext = R->D = + derp_ogcopy(S,og,R->D,imap[j]); + continue; + } + } + dnext = R->D = derpcopy(S, ce, R->D); + } + while(++i < nzc1); + } + if (D || R) { + if (!R) + R = new_relo1(S, dnext); + D = R->D = derpadjust(S, D, a, R->D); + if (Intcast d->e->op != f_OPVARVAL) + d->e->a = imap[d->e->a]; + } + d->d = D; + a += alen; + d->zaplen = (a > lasta00 ? a - nv1 : 0)*sizeof(real); + if (amax < a) + amax = a; + } + + static void +co_read(EdRead *R, cde *d, int k) +{ + d = (cde *)((char *)d + k*sizeof(cde)); + d->e = eread(R); + } + + static void +co_walk(Static *S, cde *d) +{ + int alen; + + if (amax1 < lasta) + amax1 = lasta; + lasta = lasta0; + last_d = 0; + PSHV(last_e = 0;) + ewalk(S, d->e, 1); + PSHV(d->ee = last_e;) + alen = lasta - lasta0; + if (imap_len < lasta) + imap_alloc(S); + comsubs(S, alen, d); + } + + static int +lpcompar(const void *a, const void *b, void *v) +{ + Not_Used(v); + return ((linpart *)a)->v.i - ((linpart *)b)->v.i; + } + + static linpart * +linpt_read(EdRead *R, int nlin) +{ + ASL *asl; + int i0, needsort; + int (*Xscanf)(EdRead*, const char*, ...); + linpart *L, *rv; + + asl = R->asl; + if (nlin <= 0) + return 0; + Xscanf = xscanf; + L = rv = (linpart *)new_mblk(htcl(nlin*sizeof(linpart))); + i0 = needsort = 0; + do { + if (Xscanf(R, "%d %lf", &L->v.i, &L->fac) != 2) + badline(R); + if (i0 > L->v.i) + needsort++; + i0 = L->v.i; + L++; + } + while(--nlin > 0); + if (needsort) + qsortv(rv, L-rv, sizeof(linpart), lpcompar, NULL); + return rv; + } + + static int +funnelkind(Static *S, int i0, int *ip) +{ + cexp *ce, *cej, *ce2; + dv_info *dvi; + int i, j, k, k1, k2, nzc0, rv; + int *vr, *vre; + linarg *la, **lap, **lape; + linpart *L, *Le; + list *tl; + range *r; + ASLTYPE *asl = S->asl; + + ce = cexps + i0; + ce->vref = 0; + if (!(nzc0 = nzc)) + return nocopy; + k1 = k2 = nzcperm(S); + dvi = asl->P.dv; + if (i0 >= Ncom || !dvi[i0].nl) + dvi = 0; + ce2 = 0; /* used to get funnelkind 2 right in rare cases */ + for(i = k = 0; i < nzc; i++) + if ((j = zci[i]) < nv0x) { + if (k >= maxfwd) { + k = 0; + break; + } + vrefx[k++] = j; + } + else { + cej = cexps + (j -= nv0x); + if (!(vr = cej->vref)) { + if (dvi && j < Ncom && (tl = cej->cref)) { + if (dvi[j].ll) { + for(; tl; tl = tl->next) + if (!zc[tl->item.i]++) + zci[nzc++] = tl->item.i; + } + else { + cej->vref = (int*)ce2; + ce2 = cej; + } + } + continue; + } + if (!*++vr) { + k = 0; + break; + } + vre = vr + *vr; + while(++vr <= vre) + if (!zc[*vr]++) + zci[nzc++] = *vr; + } + if (k2 < k) + k2 = k; + k2 += 2; + if (k2 > nvref) { + nvref = (maxfwd + 1)*(ncom_togo < vrefGulp + ? ncom_togo : vrefGulp); + if (nvref < k2) + nvref = 2*k2; + vrefnext = (int *)M1alloc(nvref*Sizeof(int)); + } + vr = ce->vref = vrefnext; + vrefnext += k2; + nvref -= k2; + vr[0] = k1; + vr[1] = k; + vr += 2; + i = 0; + if (k) + while(i < k) + *vr++ = vrefx[i++]; + else + while(i < k1) + *vr++ = zci[i++]; + while(nzc > nzc0) + zc[zci[--nzc]] = 0; + if (!nocopy) { + k1 = 0; + if (i0 < Ncom) { + if ((lap = (asl->P.dv + i0)->nl)) + while((la = *lap++)) + if (la->nnz > 1) + k1++; + } + else { + r = asl->P.Split_ce[i0-Ncom].r; + lap = r->lap; + lape = lap + r->nintv; + while(lap < lape) + if ((*lap++)->nnz > 1) + k1++; + } + if (k) { + if (nderp > 3*(k+k1)) { + if (!ce2) + goto ret2; + while((cej = ce2)) { + ce2 = (cexp*)cej->vref; + cej->vref = 0; + for(tl = cej->cref; tl; tl = tl->next) + if (!zc[tl->item.i]++) + zci[nzc++] = tl->item.i; + } + for(i = k = 0; i < nzc; i++) { + if ((j = zci[i]) < nv0x) { + vrefx[k++] = j; + continue; + } + cej = cexps + (j -= nv0x); + if ((vr = cej->vref)) { + if (*++vr) { + for(vre = vr + *vr; ++vr < vre; ) + if (!zc[*vr]++) + zci[nzc++] = *vr; + } + } + if (j < Ncom) { + for(tl = cej->cref; tl; tl = tl->next) + if (!zc[tl->item.i]++) + zci[nzc++] = tl->item.i; + } + if ((L = cej->L)) { + for(Le = L + cej->nlin; L < Le; ++L) + if (!zc[L->v.i]++) + zci[nzc++] = L->v.i; + } + } + while(nzc > nzc0) + zc[zci[--nzc]] = 0; + ret2: + *ip = k; + return 2; + } + } + } + while((cej = ce2)) { + ce2 = (cexp*)cej->vref; + cej->vref = 0; + } + rv = 0; + if (nocopy || nderp > 3*(nzc0+k1)) + rv = 1; + return rv; + } + + static void +cexp_read(EdRead *R, int k, int nlin) +{ + cexp *ce; + expr *e, *zero; + Static *S = (Static *)R->S; + ASLTYPE *asl = S->asl; + + ce = cexps + k - nv0x; + ce->L = linpt_read(R, ce->nlin = nlin); + e = eread(R); + if (e->op == (efunc *)f_OPVARVAL) { + /* avoid confusion with very simple defined variables */ + zero = (expr*)new_expr_n(S, 0.); + zero->op = (efunc*)f_OPNUM; + e = new_expr(S, 0, e, zero); + } + ce->e = e; + } + + static cplist * +funnelderp(Static *S, int a, cplist *cl0) +{ + cplist *cl; + ASL *asl = S->a; + + new_derp(S, a, lasta0, 0); + cl = (cplist *)mem(sizeof(cplist)); + cl->next = cl0; + cl->ca.i = imap[a]; + cl->cfa = last_d->c.rp = (real *)mem(sizeof(real)); + return cl; + } + + static void +cexp_walk(Static *S, int k) +{ + ASLTYPE *asl = S->asl; + cexp *ce; + cplist *cl; + dv_info *dvi; + expr *e; + funnel *f, **fp; + int a, *ap, fk, i, j, ka, la0, na, nlin; + linarg *la, **lap, **lape; + linpart *L, *Le; + ograd *og; + range *r; + + /* for now, assume all common exprs need to be differentiated */ + + ce = cexps + k; + nlin = ce->nlin; + L = ce->L; + nocopy = 0; + last_d = 0; + la0 = lasta; + nderps += nderp; + nderp = 0; + e = ce->e; + switch(Intcast e->op) { + case f_OPVARVAL: + case f_OPNUM: + ap = 0; + break; + default: + ap = &e->a; + } + PSHV(last_e = 0;) + ewalk(S, e, 1); + PSHV(ce->ee = last_e;) + ka = k + nv0x; + if ((na = lasta - la0)) + ce->z.i = la0; + else { + ce->z.i = k < Ncom ? ka : ((expr_vx*)varp[k-Ncom])->a0; + na = 1; + } + ce->zlen = na; + j = ap ? *ap : ka; + if (nlin) { + if (nlin == 1) + new_derp(S, L->v.i, j, &L->fac); + else if (k < Ncom) { + dvi = asl->P.dv + k; + if (dvi->lt) + new_derp(S, dvi->lt->v->a, j, &dvi->scale); + } + for(Le = L + nlin; L < Le; L++) + if (!zc[i = L->v.i]++) + zci[nzc++] = i; + } + i = 0; /* only needed to shut up an erroneous warning */ + if ((fk = funnelkind(S, k, &i))) { + /* arrange to funnel */ + fp = ka < nv0b ? &f_b : ka < nv0c ? &f_c : &f_o; + ce->funneled = f = (funnel *)mem(sizeof(funnel)); + f->next = *fp; + *fp = f; + f->ce = ce; + if (imap_len < lasta) + imap_alloc(S); + if (fk == 1) { + f->fulld = last_d; + a = lasta00; + for(i = nzc; --i >= 0; ) + if ((j = zci[i]) >= nv0x && j < max_var1) + imap[varp[j-nv0x]->a] = a++; + na += a - nv1; + f->fcde.zaplen = na*sizeof(real); + i = nzc; + derpadjust(S, last_d, a, 0); + } + else { + f->fulld = 0; + f->fcde.e = e; + comsubs(S, na, &f->fcde); + memcpy(zci, vrefx, i*sizeof(int)); + } + last_d = 0; + cl = 0; + while(i > 0) { + if ((a = zci[--i]) >= nv0x && a < max_var1) + a = varp[a-nv0x]->a; + if (a != noa) + cl = funnelderp(S, a, cl); + } + if (k < Ncom) { + if ((lap = (asl->P.dv + k)->nl)) + while((la = *lap++)) + if (la->v) + cl = funnelderp(S, la->v->a, cl); + } + else { + r = asl->P.Split_ce[k-Ncom].r; + lap = r->lap; + lape = lap + r->nintv; + while(lap < lape) + if ((la = *lap++)->nnz > 1) + cl = funnelderp(S, la->v->a, cl); + } + f->cl = cl; + varp[k]->a = lasta0++; + lasta = lasta0; + } + lasta0 = lasta; + if (!(ce->d = last_d) + && e->op == OPNUM + && (og = asl->P.dv[k].ll) + && og->varno < 0) + ((expr_n*)e)->v = 0; + while(nzc > 0) + zc[zci[--nzc]] = 0; + --ncom_togo; + } + +#ifdef PSHVREAD + static expr * +hvadjust(expr *e) +{ + expr *e0; + + for(e0 = 0; e; e = e->bak) { + e->fwd = e0; + e0 = e; + e->a = e->dO.i; + } + return e0; + } + + static void +co_adjust(ps_func *p, int n) +{ + ps_func *pe; + psb_elem *b, *be; + psg_elem *g, *ge; + + for(pe = p + n; p < pe; p++) { + for(b = p->b, be = b + p->nb; b < be; b++) + b->D.ef = hvadjust(b->D.ee); + for(g = p->g, ge = g + p->ng; g < ge; g++) { + g->esum.op = (efunc_n*)OPNUM; + for(b = g->E, be = b + g->ns; b < be; b++) + b->D.ef = hvadjust(b->D.ee); + } + } + } +#else /*PSHVREAD*/ +#define co_adjust(a,b) /*nothing*/ +#endif /*PSHVREAD*/ + + static void +ifadjust(ASLTYPE *asl, expr_if *e) +{ + real *Adjoints = asl->i.adjoints_; + for(; e; e = e->next) { + e->Tv.rp = &Adjoints[e->Tv.i]; + e->Fv.rp = &Adjoints[e->Fv.i]; +#ifdef PSHVREAD + e->Tf = hvadjust(e->Te); + e->Ff = hvadjust(e->Fe); +#endif + } + } + + static void +vargadjust(ASLTYPE *asl, expr_va *e) +{ + de *d; + real *Adjoints = asl->i.adjoints_; + + for(; e; e = e->next) { + for(d = e->L.d; d->e; d++) { + d->dv.rp = &Adjoints[d->dv.i]; +#ifdef PSHVREAD + d->ef = hvadjust(d->ee); +#endif + } + } + } + + static void +funneladj1(ASLTYPE *asl, funnel *f) +{ + real *a = adjoints; + derp *d; + cplist *cl; + funnel *fnext; + + for(; f; f = fnext) { + fnext = f->next; + f->next = 0; + if ((d = f->fulld)) { + f->fcde.d = d; + do { + d->a.rp = &a[d->a.i]; + d->b.rp = &a[d->b.i]; + } + while((d = d->next)); + } + for(cl = f->cl; cl; cl = cl->next) + cl->ca.rp = &a[cl->ca.i]; + } + } + + static void +funneladjust(Static *S) +{ + cexp *c, *ce; + linpart *L, *Le; + ASLTYPE *asl = S->asl; + + c = cexps; + for(ce = c + Ncom; c < ce; c++) { + if ((L = c->L)) + for(Le = L + c->nlin; L < Le; L++) + L->v.vp = (void*)&var_e[L->v.i]; +#ifdef PSHVREAD + c->ef = hvadjust(c->ee); +#endif + } +#ifdef PSHVREAD + for(ce += asl->P.ndvspout; c < ce; c++) + c->ef = hvadjust(c->ee); +#endif + + funneladj1(asl, f_b); + funneladj1(asl, f_c); + funneladj1(asl, f_o); + } + + static void +cg_zzap(ASLTYPE *asl) +{ + cgrad *cg, **cgp,**cgp1, **cgpe; + + cgp1 = Cgrad; + cgpe = cgp1 + asl->i.n_con0; + while(cgp1 < cgpe) { + cgp = cgp1++; + while((cg = *cgp)) + if (cg->coef) + cgp = &cg->next; + else + *cgp = cg->next; + } + } + + static void +adjust_compl_rhs(ASLTYPE *asl) +{ + cde *C; + expr *e; + int *Cvar, i, j, nc, stride; + real *L, *U, t, t1; + + L = LUrhs; + if ((U = Urhsx)) + stride = 1; + else { + U = L + 1; + stride = 2; + } + C = con_de; + Cvar = cvar; + nc = n_con; + for(i = nlc; i < nc; i++) + if (Cvar[i] && (e = C[i].e) && Intcast e->op == f_OPNUM + && (t = ((expr_n*)e)->v) != 0.) { + t1 = t; + if (L[j = stride*i] > negInfinity) { + L[j] -= t; + t1 = 0.; + } + if (U[j] < Infinity) { + U[j] -= t; + t1 = 0.; + } + ((expr_n*)e)->v = t1; + } + } + + static void +adjust(Static *S, int flags) +{ + derp *d, **dp; + ASLTYPE *asl = S->asl; + real *a = adjoints; + relo *r; + + for(r = relolist; r; r = r->next) { + dp = &r->D; + while((d = *dp)) { + d->a.rp = &a[d->a.i]; + d->b.rp = &a[d->b.i]; + dp = &d->next; + } + *dp = r->Dnext; + } + ifadjust(asl, iflist); + vargadjust(asl, varglist); + if (Ncom) + funneladjust(S); + co_adjust(asl->P.cps, asl->i.n_con0); + co_adjust(asl->P.ops, n_obj); + if (asl->i.n_con0 && !allJ) + cg_zzap(asl); + if (k_seen) { + if (!A_vals) + goff_comp_ASL((ASL*)asl); + else if (Fortran) + colstart_inc_ASL((ASL*)asl); + } + if (n_cc > nlcc && nlc < n_con + && !(flags & ASL_no_linear_cc_rhs_adjust)) + adjust_compl_rhs(asl); + } + + static void +br_read(EdRead *R, int nc, real *L, real *U, int *Cvar, int nv) +{ + ASL *asl; + int i, inc, j, k; + int (*Xscanf)(EdRead*, const char*, ...); + + asl = R->asl; + if (U) + inc = 1; + else { + U = L + 1; + inc = 2; + } + Xscanf = xscanf; + Xscanf(R, ""); /* purge line */ + for(i = 0; i < nc; i++, L += inc, U += inc) { + switch(edag_peek(R) - '0') { + case 0: + if (Xscanf(R,"%lf %lf",L,U)!= 2) + badline(R); + break; + case 1: + if (Xscanf(R, "%lf", U) != 1) + badline(R); + *L = negInfinity; + break; + case 2: + if (Xscanf(R, "%lf", L) != 1) + badline(R); + *U = Infinity; + break; + case 3: + *L = negInfinity; + *U = Infinity; + Xscanf(R, ""); /* purge line */ + break; + case 4: + if (Xscanf(R, "%lf", L) != 1) + badline(R); + *U = *L; + break; + case 5: + if (Cvar) { + if (Xscanf(R, "%d %d", &k, &j) == 2 + && j > 0 && j <= nv) { + Cvar[i] = j; + *L = k & 2 ? negInfinity : 0.; + *U = k & 1 ? Infinity : 0.; + break; + } + } + default: + badline(R); + } + } + } + + static expr * +aholread(EdRead *R) +{ + int i, k; + expr_h *rvh; + char *s1; + FILE *nl = R->nl; + Static *S = (Static *)R->S; + + k = getc(nl); + if (k < '1' || k > '9') + badline(R); + i = k - '0'; + while((k = getc(nl)) != ':') { + if (k < '0' || k > '9') + badline(R); + i = 10*i + k - '0'; + } + rvh = (expr_h *)mem_ASL(R->asl, memadj(sizeof(expr_h) + i)); + for(s1 = rvh->sym;;) { + if ((k = getc(nl)) < 0) { + fprintf(Stderr, + "Premature end of file in aholread, line %ld of %s\n", + R->Line, R->asl->i.filename_); + exit_ASL(R,1); + } + if (k == '\n') { + R->Line++; + if (!i) + break; + } + if (--i < 0) + badline(R); + *s1++ = k; + } + *s1 = 0; + rvh->op = (efunc *)f_OPHOL; + rvh->a = noa; + return (expr *)rvh; + } + + static expr * +bholread(EdRead *R) +{ + int i; + expr_h *rvh; + char *s; + FILE *nl = R->nl; + ASL *asl = R->asl; + Static *S = (Static *)R->S; + + if (xscanf(R, "%d", &i) != 1) + badline(R); + rvh = (expr_h *)mem(memadj(sizeof(expr_h) + i)); + s = rvh->sym; + if (fread(s, i, 1, nl) != 1) + badline(R); + s[i] = 0; + rvh->op = (efunc *)f_OPHOL; + rvh->a = noa; + for(;;) switch(*s++) { + case 0: goto break2; /* so we return at end of fcn */ + case '\n': R->Line++; + } + break2: + return (expr *)rvh; + } + + static int +qwalk(Static *S, expr *e) +{ + int i, j, k; + expr **args, **argse; + expr_f *rvf; + argpair *ap, *ape; + + top: + switch(optype[k = Intcast e->op]) { + + case 1: /* unary */ + switch(k) { + case f_OPCPOW: + if (qwalk(S, e->R.e)) + return 3; + return 0; + case f_OPUMINUS: + e = e->L.e; + goto top; + case f_OP2POW: + i = qwalk(S, e->L.e); + if (i < 2) + return i << 1; + } + return 3; + + case 2: /* binary */ + switch(k) { + case f_OPPLUS: + case f_OPMINUS: + i = qwalk(S, e->L.e); + if (i == 3) + break; + j = qwalk(S, e->R.e); + if (i < j) + i = j; + return i; + case f_OPDIV: + if (qwalk(S, e->R.e)) + return 3; + e = e->L.e; + goto top; + case f_OPMULT: + i = qwalk(S, e->L.e); + if (i < 3) { + i += qwalk(S, e->R.e); + if (i < 3) + return i; + } + } + return 3; + + case 6: /* sumlist */ + j = 0; + args = e->L.ep; + argse = e->R.ep; + while(args < argse) { + i = qwalk(S, *args++); + if (j < i) { + j = i; + if (j == 3) + break; + } + } + return j; + + case 7: /* function call */ + rvf = (expr_f *)e; + ap = rvf->ap; + ape = rvf->ape; + for(i = 0; ap < ape; ap++) { + if (qwalk(S, ap->e)) + return 3; + } + + case 9: /* OPNUM */ + return 0; + case 10:/* OPVARVAL */ + k = (expr_v *)e - S->var_e; + if (k < 0) + goto k_adj; + if (k < nv0x) + return 1; + if (k >= max_var) { + k_adj: + k = ((expr_vx*)e)->a1; + return k < 0 ? 1 : S->asl->I.v_class[k-nv0x]; + } + return S->asl->I.v_class[k - nv0x]; + } + return 3; + } + + static int +co_walkloop(Static *S, ps_func *f, int n, char *c, ograd **o) +{ + ps_func *fe; + psb_elem *b, *be; + psg_elem *g, *ge; + int k, k1, kx; + + kx = 0; + for(fe = f + n; f < fe; f++) { + if (c) { + k = *o++ != 0; + for(g = f->g, ge = g + f->ng; g < ge; g++) { + if (Intcast g->g->op != f_OP2POW) { + k = 3; + goto have_c; + } + if (g->nlin) + k = 2; + for(b = g->E, be = b + g->ns; b < be; b++) { + if ((k1 = qwalk(S, b->D.e)) > 1) { + k = 3; + goto have_c; + } + k = 2; + } + } + for(b = f->b, be = b + f->nb; b < be; b++) { + if ((k1 = qwalk(S, b->D.e)) > k) { + k = k1; + if (k == 3) + goto have_c; + } + } + have_c: *c++ = (char)k; + if (kx < k) + kx = k; + } + for(b = f->b, be = b + f->nb; b < be; b++) + co_walk(S, &b->D); + for(g = f->g, ge = g + f->ng; g < ge; g++) { + ewalk(S, g->g, 1); + for(b = g->E, be = b + g->ns; b < be; b++) + co_walk(S, &b->D); + } + } + return kx; + } + + static void +do_ewalk(Static *S) +{ + int i, *map, n, nc; + char *v, *ve; + expr_v *e, *ee, **vp, **vpe; + efunc *vv; + linarg *la, **lap; + ograd *og; + cexp *c; + ASLTYPE *asl = S->asl; + + psfind(S); + nc = max_var1 - nv0x; + noa = lasta00 + nc; + nv1 = lasta00++; + amax = lasta = lasta0 = max_var + nndv + 1; + + if (Ncom && asl->I.o_class) { + i = max_var1 - nv0x; + v = asl->I.v_class = (char*)M1alloc(i); + ve = v + i; + c = cexps; + while(v < ve) + *v++ = (char)((i = qwalk(S, (c++)->e)) ? i : 1); + } + + n = Ncom + asl->P.ndvspout; + while(nzc > 0) + zc[zci[--nzc]] = 0; + for(i = 0; i < n; i++) { +#ifdef DEBUG + if (i == izork) + printf(""); +#endif + cexp_walk(S, i); + } + + if (imap_len < lasta) + imap_alloc(S); + f_b = funnelfix(f_b); + f_c = funnelfix(f_c); + f_o = funnelfix(f_o); + + asl->I.c_class_max = co_walkloop(S, asl->P.cps, asl->i.n_con0, + asl->I.c_class, (ograd**)Cgrad); + asl->I.o_class_max = co_walkloop(S, asl->P.ops, n_obj, + asl->I.o_class, Ograd); + del_mblk(kzc, zci); + zci = zc = 0; /* for debugging */ + e = var_e; + vv = OPVARVAL; + for(ee = e + nv0x + Ncom; e < ee; e++) + e->op = vv; + lap = &asl->P.lalist; + map = imap; + while((la = *lap)) + if ((e = la->v)) { + e->op = vv; + e->a = map[e->a]; + lap = &la->lnext; + } + else { + og = la->nz; + assert(og != 0); + assert(og->next == 0); + assert(og->varno < nv0x); + la->v = var_e + og->varno; + *lap = la->lnext; + } + if (n) { + asl->P.vp = vp = varp; +#ifdef PSHVREAD + asl->P.zlsave = zl; +#endif + e = var_ex; + for(ee = e + Ncom; e < ee; e++) + *vp++ = e; + vpe = vp + asl->P.ndvspout; + while(vp < vpe) + (*vp++)->op = vv; + } + } + + int +pfg_read_ASL(ASL *a, FILE *nl, int flags) +{ + ASLTYPE *asl; + EdRead ER, *R; + Jmp_buf JB; + Static SS, *S; + cgrad *cg, **cgp; + char fname[128]; + expr_v *e; + func_info *fi; + int allG, i, i1, j, k, *ka, kseen, maxfwd1, nc, nc0, nco, nlin, nlogc, no; + int nv, nvc, nvo, nvr, nvextra, nvx, readall; + int (*Xscanf)(EdRead*, const char*, ...); + ograd *og, **ogp; + real t; + size_t *kaz, nz; + unsigned x; + + ASL_CHECK(a, asltype, who); + flagsave_ASL(a, flags); /* includes allocation of LUv, LUrhs, A_vals or Cgrad, etc. */ + asl = (ASLTYPE*)a; + S = S_init(&SS, asl); + ed_reset(asl); + SS.R = R = EdReadInit_ASL(&ER, a, nl, S); + if (flags & ASL_return_read_err) { + a->i.err_jmp_ = &JB; + i = setjmp(JB.jb); + if (i) { + a->i.err_jmp_ = 0; + return i; + } + } + if ((nlogc = a->i.n_lcon_) && !(flags & ASL_allow_CLP)) { + if (a->i.err_jmp_) + return ASL_readerr_CLP; + sorry_CLP(R, "logical constraints"); + } + if (!(flags & ASL_find_default_no_groups)) + flags |= ASL_findgroups; + Xscanf = xscanf; + readall = flags & ASL_keep_all_suffixes; + PSHV(asl->P.pshv_g1 = 1;) + k_Elemtemp = htcl(sizeof(Elemtemp)); + allJ = (flags & ASL_J_zerodrop) == 0; + allG = (flags & ASL_G_zerodrop) == 0; + wantCgroups = flags & ASL_findCgroups; + wantOgroups = flags & ASL_findOgroups; + OPNUM = r_ops[f_OPNUM]; + OPVARVAL = r_ops[f_OPVARVAL]; + if (!size_expr_n) + size_expr_n = sizeof(expr_n); + S->size_exprn = size_expr_n; + asl->P.rlist.next = asl->P.rlist.prev = (range*)&asl->P.rlist; + if (nfunc) + func_add(a); + if (binary_nl) + holread = bholread; + else + holread = aholread; + + asl->P.ncom = Ncom = comb + comc + como + comc1 + como1; + nc0 = n_con; + nc = nc0 + a->i.nsufext[ASL_Sufkind_con]; + no = n_obj; + nvc = c_vars; + nvo = o_vars; + nco = nc + no + nlogc; + if (no < 0 || nco <= 0) + scream(R, ASL_readerr_corrupt, + "pshvread: nc = %d, no = %d, nlogc = %d\n", + nc0, no, nlogc); + if (pi0) { + memset(pi0, 0, nc*sizeof(real)); + if (havepi0) + memset(havepi0, 0, nc); + } + nvextra = a->i.nsufext[ASL_Sufkind_var]; + nv0r = nvr = n_var; + asl->P.nv0_ = lasta00 = nv0x = nvx = n_var + nvextra; + max_var1 = max_var = nv = nvr + Ncom + nvextra; + combc = comb + comc; + ncom0 = ncom_togo = combc + como; + nzclim = max_var >> 3; + ncom1 = comc1 + como1; + nv0b = nvr + comb; + nv0c = nv0b + comc; + if ((maxfwd1 = maxfwd + 1) > 1) + nvref = maxfwd1*((ncom0 < vrefGulp ? ncom0 : vrefGulp) + 1); + x = nco*sizeof(cde) + no*sizeof(ograd *) + + nv*sizeof(expr_v) + + nfunc*sizeof(func_info *) + + nvref*sizeof(int) + + no; + SS.nvar0 = a->i.n_var0; + if (!(SS.nvinc = a->i.nvinc = a->i.n_var_ - SS.nvar0 + nvextra)) + SS.nvar0 += ncom0 + ncom1; + if (flags & ASL_find_co_class) + x += nco; + if (X0) + memset(X0, 0, nvr*sizeof(real)); + if (havex0) + memset(havex0, 0, nvr); + e = var_e = (expr_v *)M1zapalloc(x); + con_de = (cde *)(e + nv); + lcon_de = con_de + nc; + obj_de = lcon_de + nlogc; + Ograd = (ograd **)(obj_de + no); + var_ex = e + nvx; + var_ex1 = var_ex + ncom0; + for(k = 0; k < nv; e++) { + e->op = (efunc *)f_OPVARVAL; + e->a = k++; + } + funcs = (func_info **)(Ograd + no); + zci = 0; + zc_upgrade(S); /* initially allow nv+1 entries in zci and zc[-1] */ + vrefx = (int *)(funcs + nfunc); + if (Ncom) { + cexps = 0; + asl->P.ndvspout = 0; + memset(asl->P.dv = (dv_info*)mem(Ncom*sizeof(dv_info)), + 0, Ncom*sizeof(dv_info)); + cexp_upgrade(S, Ncom); + for(k = 0; k < Ncom; k++) + varp[k] = &var_ex[k]; + } + objtype = (char *)(vrefx + nvref); + if (flags & ASL_find_co_class) { + asl->I.o_class = objtype + no; + asl->I.c_class = asl->I.o_class + no; + } + if (nvref) { + vrefnext = vrefx + maxfwd1; + nvref -= maxfwd1; + } + if (n_cc && !cvar) + cvar = (int*)M1alloc(nc*sizeof(int)); + if (cvar) + memset(cvar, 0, nc*sizeof(int)); + last_d = 0; + ka = 0; + kaz = 0; + kseen = 0; + nz = 0; + for(;;) { + ER.can_end = 1; + i = edag_peek(R); + if (i == EOF) { + fclose(nl); + do_ewalk(S); + if (imap) + del_mblk(kimap, imap); + /* Make amax long enough for nlc to handle */ + /* var_e[i].a for common variables i. */ + if (ncom0) { + i = comb + como; + if (i < combc) + i = combc; + if ((i += noa + 1) > amax) + amax = i; + } + adjoints = (real *)M1zapalloc(amax*Sizeof(real)); + adjoints_nv1 = &adjoints[nv1]; + nderps += nderp; + adjust(S, flags); + nzjac = nz; + if (!Lastx) + Lastx = (real *)M1alloc(nvx*sizeof(real)); +#ifdef PSHVREAD + a->i.ncxval = (int*)M1zapalloc(nco*sizeof(int)); + a->i.noxval = a->i.ncxval + nc; + a->p.Conival = conpival_ASL; + a->p.Conival_nomap = conpival_nomap_ASL; + a->p.Congrd = conpgrd_ASL; + a->p.Congrd_nomap = conpgrd_nomap_ASL; + a->p.Objval = a->p.Objval_nomap = objpval_ASL; + a->p.Objgrd = a->p.Objgrd_nomap = objpgrd_ASL; + a->p.Conval = conpval_ASL; + a->p.Jacval = jacpval_ASL; + a->p.Lconval= lconpval_ASL; + a->p.Hvcomp = a->p.Hvcomp_nomap = hvpcomp_ASL; + a->p.Hvcompd = hvpcompd_ASL; + a->p.Hvcompde = hvpcompde_ASL; + a->p.Hvcomps = hvpcomps_ASL; + a->p.Hvcompse = hvpcompse_ASL; + a->p.Hvinit = a->p.Hvinit_nomap = hvpinit_ASL; + a->p.Hvinite = a->p.Hvinite_nomap = hvpinite_ASL; + a->p.Hesset = hes_setup; + a->p.Xknown = xp2known_ASL; + a->p.Duthes = a->p.Duthes_nomap = duthes_ASL; + a->p.Duthese = a->p.Duthese_nomap = duthese_ASL; + a->p.Fulhes = a->p.Fulhes_nomap = fullhes_ASL; + a->p.Fulhese = a->p.Fulhese_nomap = fullhese_ASL; + a->p.Sphes = a->p.Sphes_nomap = sphes_ASL; + a->p.Sphese = a->p.Sphese_nomap = sphese_ASL; + a->p.Sphset = a->p.Sphset_nomap = sphes_setup_ASL; +#else + a->p.Xknown = xp1known_ASL; +#endif + return prob_adj_ASL(a); + } + ER.can_end = 0; + k = -1; + switch(i) { + case 'C': + Xscanf(R, "%d", &k); + if (k < 0 || k >= nc0) + badline(R); + co_read(R, con_de, k); + break; + case 'F': + if (Xscanf(R, "%d %d %d %127s", + &i, &j, &k, fname) != 4 + || i < 0 || i >= nfunc) + badline(R); + if ((fi = func_lookup(a, fname,0))) { + if (fi->nargs != k && fi->nargs >= 0 + && (k >= 0 || fi->nargs < -(k+1))) + scream(R, ASL_readerr_argerr, + "function %s: disagreement of nargs: %d and %d\n", + fname,fi->nargs, k); + } + else { + fi = (func_info *)mem(sizeof(func_info)); + fi->ftype = j; + fi->nargs = k; + fi->funcp = 0; + fi->name = (Const char *)strcpy((char*) + mem(memadj(strlen(fname)+1)), + fname); + } + if (!fi->funcp && !(fi->funcp = dynlink(fname))) + scream(R, ASL_readerr_unavail, + "function %s not available\n", + fname); + funcs[i] = fi; + break; + case 'L': + Xscanf(R, "%d", &k); + if (k < 0 || k >= nlogc) + badline(R); + co_read(R, lcon_de, k); + ewalk(S, lcon_de[k].e, 0); + break; + case 'V': + if (Xscanf(R, "%d %d %d", &k, &nlin, &j) != 3) + badline(R); + if (k >= SS.nvar0) + k += SS.nvinc; + if (k < nvr || k >= nv) + badline(R); + cexp_read(R, k, nlin); + break; + case 'G': + if (Xscanf(R, "%d %d", &j, &k) != 2 + || j < 0 || j >= no || k <= 0 || k > nvo) + badline(R); + ogp = Ograd + j; + while(k--) { + if (Xscanf(R, "%d %lf", &i, &t) != 2) + badline(R); + if (allG || t) { + *ogp = og = (ograd *) + mem(sizeof(ograd)); + ogp = &og->next; + og->varno = i; + og->coef = t; + } + } + *ogp = 0; + break; + case 'J': + if (Xscanf(R, "%d %d", &j, &k) != 2 + || j < 0 || j >= nc || k <= 0 || k > nvc) + badline(R); + nz += k; + if (A_vals) { + j += Fortran; + if (ka) { + while(k--) { + if (Xscanf(R, "%d %lf", + &i, &t) != 2) + badline(R); + i1 = ka[i]++; + A_vals[i1] = t; + A_rownos[i1] = j; + } + } + else { + while(k--) { + if (Xscanf(R, "%d %lf", + &i, &t) != 2) + badline(R); + i1 = kaz[i]++; + A_vals[i1] = t; + A_rownos[i1] = j; + } + } + break; + } + cgp = Cgrad + j; + j = 0; + while(k--) { + if (kseen) { + if (Xscanf(R, "%d %lf", &i, &t) != 2) + badline(R); + } + else + if (Xscanf(R, "%d %d %lf", &i, &j, &t) != 3) + badline(R); + *cgp = cg = (cgrad *) + mem(sizeof(cgrad)); + cgp = &cg->next; + cg->varno = i; + cg->goff = j; + cg->coef = t; + } + *cgp = 0; + break; + case 'O': + if (Xscanf(R, "%d %d", &k, &j) != 2 + || k < 0 || k >= no) + badline(R); + objtype[k] = j; + co_read(R, obj_de, k); + break; + case 'S': + Suf_read_ASL(R, readall); + break; + case 'r': + br_read(R, asl->i.n_con0, LUrhs, Urhsx, cvar, nvr); + break; + case 'b': + br_read(R, asl->i.n_var0, LUv, Uvx, 0, 0); + break; + case 'K': + case 'k': + k_seen = ++kseen; + if (ka_read_ASL(a, R, i, &ka, &kaz)) + badline(R); + break; + case 'x': + if (!Xscanf(R,"%d",&k) + || k < 0 || k > nvr) + badline(R); + if (!X0 && want_xpi0 & 1) { + x = nvx*sizeof(real); + if (want_xpi0 & 4) + x += nvx; + X0 = (real *)M1zapalloc(x); + if (want_xpi0 & 4) + havex0 = (char*)(X0 + nvx); + } + while(k--) { + if (Xscanf(R, "%d %lf", &j, &t) != 2 + || j < 0 || j >= nvr) + badline(R); + if (X0) { + X0[j] = t; + if (havex0) + havex0[j] = 1; + } + } + break; + case 'd': + if (!Xscanf(R,"%d",&k) + || k < 0 || k > nc0) + badline(R); + if (!pi0 && want_xpi0 & 2) { + x = nc*sizeof(real); + if (want_xpi0 & 4) + x += nc; + pi0 = (real *)M1zapalloc(x); + if (want_xpi0 & 4) + havepi0 = (char*)(pi0 + nc); + } + while(k--) { + if (Xscanf(R, "%d %lf", &j, &t) != 2 + || j < 0 || j >= nc0) + badline(R); + if (pi0) { + pi0[j] = t; + if (havepi0) + havepi0[j] = 1; + } + } + break; + default: + badline(R); + } + } + } +#undef asl + +#ifdef PSHVREAD + + static int +heswork(expr *e) +{ + argpair *da, *dae; + int i, j, n; + expr *e1, **ep; + de *d; + + n = 0; + if (e && e->op != OPVARVAL && e->op != OPNUM) do { + switch(e->a) { + + case Hv_timesR: + case Hv_timesL: + case Hv_negate: + case Hv_plterm: + n += 4; + break; + + case Hv_binaryR: + case Hv_unary: + n += 6; + break; + + case Hv_timesLR: + n += 10; + break; + + case Hv_divLR: + n += 13; + break; + + case Hv_binaryLR: + n += 14; + break; + + case Hv_vararg: + d = ((expr_va*)e)->L.d; + i = heswork(d->e); + while((e1 = (++d)->e)) { + j = heswork(e1); + if (i < j) + i = j; + } + n = i + 2; + break; + + case Hv_if: + i = heswork(((expr_if *)e)->Tf); + j = heswork(((expr_if *)e)->Ff); + if (i < j) + i = j; + n += i + 2; + break; + + + case Hv_sumlist: + for(ep = e->R.ep; (e1 = *ep); ep++) + n++; + break; + + case Hv_func: + da = ((expr_f *)e)->da; + dae = ((expr_f *)e)->dae; + i = dae - da; + n += i*(i+4); + break; + + case Hv_plusR: + case Hv_plusL: + case Hv_minusR: + n += 2; + break; + + case Hv_plusLR: + case Hv_minusLR: + n += 3; + break; + + default:/*DEBUG*/ + fprintf(Stderr, "bad e->a = %d in heswork\n", e->a); + exit(1); + } + } while((e = e->fwd)); + return n; + } + + static cexp * +hesfunnel(Static *S, int *ip, int nrefs, ograd **ogp, linarg ***lapp, linarg ***lapep) +{ + cexp *c; + expr *e; + int hw, i, k, m, n; + linarg *la, **lap, **lape; + ograd *og; + dv_info *dvi; + /*list *L;*//*20170615*/ + range *r; + split_ce *cs; /*20170615*/ + ASLTYPE *asl = S->asl; + + i = *ip; + c = cexps + i; +#if 1 /*20170615*/ + if (c->cref) + return 0; +#endif + n = 0; + if (i >= Ncom) { + cs = &asl->P.Split_ce[i-Ncom]; +#if 1 /*20170615*/ + if (cs->ce) + return 0; +#endif + r = cs->r; + lap = r->lap; + lape = lap + (n = r->n); + } + else if (!(lap = (dvi = asl->P.dv + i)->nl)) { + asl->P.linmultr++; + og = dvi->ll; + if (og->varno < 0) + og = og->next; + for(*ogp = og; og; og = og->next) + n++; + if (n > 1 && asl->p.hffactor > 0) { + asl->P.linhesfun++; + *ip = n; + return c; + } + return 0; + } + else { + lape = lap; + while(*lape) + lape++; + n = lape - lap; + } + if (!(e = c->ef)) + return 0; + *lapp = lap; + *lapep = lape; + *ogp = 0; + asl->P.nlmultr++; + while(lap < lape) { + la = *lap++; + for(og = la->nz; og; og = og->next) + if (!zc[og->varno]++) + zci[nzc++] = og->varno; + } + m = nzc; + while(nzc > 0) + zc[zci[--nzc]] = 0; +#if 0 /*20170615*/ + for(L = c->cref; L; L = L->next) + n++; +#endif + if (m > n) + m = n; + if ((k = m*nrefs - n) <= 0) + return 0; + hw = heswork(e); + if (k*hw <= nrefs*n*(n+3)*asl->p.hffactor) + return 0; + *ip = n; + asl->P.nlhesfun++; + return c; + } + + ASL_pfgh * +pscheck_ASL(ASL *a, const char *who1) +{ + ASL_CHECK(a, ASL_read_pfgh, who1); + return (ASL_pfgh*)a; + } + + static void +hes_setup(ASL *a, int flags, int obj, int nobj, int con, int ncon) +{ + range *r, *r0; + int *c, *ce, *cei, i, k, n, n0, nmax, nvx, *z; + hes_fun *hf, *hfth; + cexp *C; + linarg *la, **lap, **lape; + ograd *og; + list *L; + expr_v **vp; + psb_elem *b; + ASL_pfgh *asl; + Static SS, *S; + + asl = pscheck_ASL(a, "hes_setup"); + + S = S_init(&SS, asl); + Ncom = asl->P.ncom; + nv0x = nvx = asl->P.nv0_; + max_var = nvx + Ncom; + max_var1 = max_var + asl->P.ndvspout; + nzclim = max_var1 >> 3; + varp = asl->P.vp; + zl = asl->P.zlsave; + + zc_upgrade(S); /* make zc and zci available */ + n = nmax = 0; + r0 = (range*)&asl->P.rlist; + for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) { + if (r->n <= 0) { + r->cei = 0; + continue; + } + if (nmax < r->n) + nmax = r->n; + n0 = 0; + cei = 0; + for(b = r->refs; b; b = b->next) { + i = b->conno; + if (i < 0) { + i = -i - 2; + if (i < obj || i-obj >= nobj) + continue; + } + else { + if (i < con || i-con >= ncon) + continue; + } + if ((c = b->ce)) { + if (!n0) { + cei = c; + n0 = *c; + } + ce = c + *c; + while(c < ce) + if (!zc[i = *++c]++) + zci[n++] = i; + } + } + if (n0 < n) + cei = 0; + if ((r->cei = cei)) { + while(n > 0) + zc[zci[--n]] = 0; + continue; + } + if (!n) + continue; + r->cei = cei = (int*)mem((n+1)*sizeof(int)); + *cei++ = n; + if (n > 1) + qsortv(zci, n, sizeof(int), hscompar, S); + cei += n; + do zc[*--cei = zci[--n]] = 0; + while(n > 0); + } + z = zc + nvx; + k = -1; + for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) { + if (!(cei = r->cei)) + continue; + c = cei + 1; + ce = c + *cei; + do { + if (z[i = *c++]++ && k < i) + k = i; + } while(c < ce); + } + hfth = 0; + asl->P.linmultr = asl->P.linhesfun = 0; + asl->P.nlmultr = asl->P.nlhesfun = 0; + lap = lape = 0; /* silence false "used uninitialized" warning */ + while (k >= 0) + if (z[i = k--] > 1 + && (C = hesfunnel(S,&i,z[i],&og,&lap,&lape))) { + hf = (hes_fun*)mem(sizeof(hes_fun)); + hf->hfthread = hfth; + C->hfun = hfth = hf; + hf->c = C; + hf->n = i; + if ((hf->og = og)) { + hf->vp = 0; + hf->grdhes = 0; + } + else { + hf->vp = vp = (expr_v **) + mem(i*sizeof(expr_v*)); + while(lap < lape) { + la = *lap++; + *vp++ = la->v; + } + for(L = C->cref; L; L = L->next) + *vp++ = varp[L->item.i - nvx]; + i += i*i; + hf->grdhes = (real*)mem(i*sizeof(real)); + } + } + if ((asl->I.hesthread = hfth)) { + asl->I.x0kind_init = ASL_need_funnel; + if (x0kind != ASL_first_x) + x0kind |= ASL_need_funnel; + } + else + asl->I.x0kind_init = 0; + del_mblk(kzc, zci); + asl->P.hes_setup_called = 1 | (flags << 2); + asl->P.khesoprod = 5; + asl->P.hop_free = 0; + asl->P.nmax = nmax; + if (!(flags & 1)) + return; + i = htcl(n = nvx*sizeof(range*)); + k = htcl(3*n); + if (k == i + 1) { + asl->P.rtodo = (range**)new_mblkzap(asl, k); + asl->P.utodo = (uHeswork **)(asl->P.rtodo + nvx); + asl->P.otodo = (Hesoprod**)(asl->P.utodo + nvx); + } + else { + asl->P.rtodo = (range**)new_mblkzap(asl, i); + asl->P.utodo = (uHeswork**)new_mblkzap(asl, i); + asl->P.otodo = (Hesoprod**)new_mblkzap(asl, i); + } + k = htcl(nmax*(sizeof(int)+sizeof(real))); + asl->P.dOscratch = (real *)new_mblkzap(asl, k); + asl->P.iOscratch = (int *)(asl->P.dOscratch + nmax); + } + +#endif /* PSHVREAD */ + +#ifdef __cplusplus +} +#endif diff --git a/src/solvers2/pfghread.c b/src/solvers2/pfghread.c index 11e963c..52b2c7c 100644 --- a/src/solvers2/pfghread.c +++ b/src/solvers2/pfghread.c @@ -3182,7 +3182,7 @@ ewalk(Static *S, expr *e, uint *deriv, uint atop) i1 = opg3[2*i+1]; for(j = 0; j < k; ++j) { j1 = opg3[2*j+1]; - *fh++ = i1 <= j1 ? (i1*(i1+1)>>1) + j1 + *fh++ = i1 >= j1 ? (i1*(i1+1)>>1) + j1 : (j1*(j1+1)>>1) + i1; } } diff --git a/src/solvers2/sphes.c b/src/solvers2/sphes.c index d7edd77..b13820f 100644 --- a/src/solvers2/sphes.c +++ b/src/solvers2/sphes.c @@ -1578,7 +1578,6 @@ sphes_ew_ASL(EvalWorkspace *ew, SputInfo **pspi, real *H, int nobj, real *ow, re } } *rtodo++ = 0; /* reset */ - V[i].dO = 1.; while((uhw = uhwi)) { uhwi = uhwi->next; r = uhw->r; @@ -1623,7 +1622,6 @@ sphes_ew_ASL(EvalWorkspace *ew, SputInfo **pspi, real *H, int nobj, real *ow, re *utodoj = uhw; } } - V[i].dO = 0.; hop1 = *otodoi; *otodoi++ = 0;