/**************************************************************** Copyright (C) 1997-2001 Lucent Technologies All Rights Reserved Permission to use, copy, modify, and distribute this software and its documentation for any purpose and without fee is hereby granted, provided that the above copyright notice appear in all copies and that both that the copyright notice and this permission notice and warranty disclaimer appear in supporting documentation, and that the name of Lucent or any of its entities not be used in advertising or publicity pertaining to distribution of the software without specific, written prior permission. LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. ****************************************************************/ #include "jac2dim.h" #include "opnos.hd" #ifdef __cplusplus extern "C" { #endif #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 #ifdef KR_headers ed_reset(a) ASL *a; #else ed_reset(ASL *a) #endif { 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 #ifdef KR_headers fscream(R, name, nargs, kind) EdRead *R; char *name; int nargs; char *kind; #else fscream(EdRead *R, const char *name, int nargs, char *kind) #endif { scream(R, ASL_readerr_argerr, "line %ld: attempt to call %s with %d %sargs\n", R->Line, name, nargs, kind); } static void #ifdef KR_headers sorry_CLP(R, what) EdRead *R; char *what; #else sorry_CLP(EdRead *R, char *what) #endif { 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 efunc f_OPPLTERM, f_OPHOL, f_OPVARVAL, f_OPFUNCALL; extern sfunc f_OPIFSYM; static void #ifdef KR_headers new_derp(a, b, c) int a; int b; real *c; #else new_derp(int a, int b, real *c) #endif { 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 * #ifdef KR_headers new_relo(e, Dnext, ap) expr *e; derp *Dnext; int *ap; #else new_relo(expr *e, derp *Dnext, int *ap) #endif { 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 * #ifdef KR_headers new_relo1(Dnext) derp *Dnext; #else new_relo1(derp *Dnext) #endif { relo *r; r = (relo *)mem(sizeof(relo)); r->next = relolist; relolist = r; r->D = 0; r->Dnext = Dnext; return r; } static expr * #ifdef KR_headers new_expr(opcode, L, R, deriv) int opcode; expr *L; expr *R; int deriv; #else new_expr(int opcode, expr *L, expr *R, int deriv) #endif { 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; 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 * #ifdef KR_headers eread(R, deriv) EdRead *R; int deriv; #else eread(EdRead *R, int deriv) #endif { char **sa; int a0, a1, i, i1, j, j1, k, kd, kd2, ks, numargs, symargs; int *at, *nn, *nn0; real *b, **fh, *hes, r, *ra; expr *L, *arg, **args, **args1, **argse, *rv; expr_n *rvn; expr_va *rva; plterm *p; de *d; derp *dsave; efunc *op; expr_if *rvif; expr_f *rvf; func_info *fi; arglist *al; argpair *ap, *da, *sap; char *dig; short sh; long L1; static real dvalue[] = { #include "dvalue.hd" }; 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) == 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; } ra = (real *)mem(sizeof(arglist) + (k + kd + ks)*sizeof(argpair) + kd*kd*sizeof(real *) + (numargs+kd+kd2)*sizeof(real) + symargs*sizeof(char *) + j*sizeof(int)); 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) == 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", &sh) != 1) badline(R); r = sh; 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, "%d", &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",&sh) != 1) badline(R); r = sh; 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 (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 = deriv) 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 * #ifdef KR_headers new_list(nxt) list *nxt; #else new_list(list *nxt) #endif { 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 * #ifdef KR_headers funnelfix(f) funnel *f; #else funnelfix(funnel *f) #endif { 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 * #ifdef KR_headers derpadjust(d0, a, dnext) derp *d0; int a; derp *dnext; #else derpadjust(derp *d0, int a, derp *dnext) #endif { 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 * #ifdef KR_headers derpcopy(ce, dnext) cexp *ce; derp *dnext; #else derpcopy(cexp *ce, derp *dnext) #endif { 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 #ifdef KR_headers compar(a, b, v) char *a, *b, *v; #else compar(const void *a, const void *b, void *v) #endif { Not_Used(v); return *(int*)a - *(int*)b; } static void #ifdef KR_headers comsubs(alen, d) int alen; cde *d; #else comsubs(int alen, cde *d) #endif { 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 #ifdef KR_headers co_read(R, d, wd) EdRead *R; cde *d; int wd; #else co_read(EdRead *R, cde *d, int wd) #endif { 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 * #ifdef KR_headers linpt_read(R, nlin) EdRead *R; int nlin; #else linpt_read(EdRead *R, int nlin) #endif { linpart *L, *rv; if (nlin <= 0) return 0; 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 #ifdef KR_headers funnelkind(ce, ip) cexp *ce; int *ip; #else funnelkind(cexp *ce, int *ip) #endif { 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 #ifdef KR_headers cexp_read(R, k, nlin) EdRead *R; int k; int nlin; #else cexp_read(EdRead *R, int k, int nlin) #endif { 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; } 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 #ifdef KR_headers cexp1_read(R, j, k, nlin) EdRead *R; int j; int k; int nlin; #else cexp1_read(EdRead *R, int j, int k, int nlin) #endif { 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 * #ifdef KR_headers hvadjust(e) expr *e; #else hvadjust(expr *e) #endif { expr *e0; for(e0 = 0; e; e = e->bak) { e->fwd = e0; e0 = e; e->a = e->dO.i; } return e0; } static void #ifdef KR_headers co_adjust(d, n) cde *d; int n; #else co_adjust(cde *d, int n) #endif { cde *de1; for(de1 = d + n; d < de1; d++) d->ef = hvadjust(d->ee); } static void #ifdef KR_headers ifadjust(e) expr_if *e; #else ifadjust(expr_if *e) #endif { 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 #ifdef KR_headers vargadjust(e) expr_va *e; #else vargadjust(expr_va *e) #endif { 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 #ifdef KR_headers funneladj1(f) funnel *f; #else funneladj1(funnel *f) #endif { 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 goff_comp(VOID) { int *ka = A_colstarts + 1; cgrad **cgx, **cgxe; cgrad *cg; cgx = Cgrad; cgxe = cgx + asl->i.n_con0; while(cgx < cgxe) for(cg = *cgx++; cg; cg = cg->next) cg->goff = ka[cg->varno]++; } static void colstart_inc(VOID) { int *ka, *kae; ka = A_colstarts; kae = ka + asl->i.n_var0; while(ka <= kae) ++*ka++; } static void zerograd_chk(VOID) { int j, n, nv, *z, **zg; ograd *og, **ogp, **ogpe; if (!(nv = asl->i.nlvog)) nv = nv0; zerograds = 0; ogp = Ograd; ogpe = ogp + (j = n_obj); while(ogp < ogpe) { og = *ogp++; n = 0; while(og) { j += og->varno - n; n = og->varno + 1; if (n >= nv) break; og = og->next; } if (n < nv) j += nv - n; } if (j == n_obj) return; zerograds = zg = (int **)mem(n_obj*sizeof(int*)+j*sizeof(int)); z = (int*)(zg + n_obj); ogp = Ograd; while(ogp < ogpe) { *zg++ = z; og = *ogp++; n = 0; while(og) { while(n < og->varno) *z++ = n++; og = og->next; if (++n >= nv) break; } while(n < nv) *z++ = n++; *z++ = -1; } } static void adjust_compl_rhs(VOID) { cde *C; expr *e; int *Cvar, i, j, nc, stride; real *L, *U, t; 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.) { ((expr_n*)e)->v = 0.; if (L[j = stride*i] > negInfinity) L[j] -= t; if (U[j] < Infinity) U[j] -= t; } } static void #ifdef KR_headers adjust(flags) int flags; #else adjust(int flags) #endif { 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 (n_obj) zerograd_chk(); if (k_seen) if (!A_vals) goff_comp(); else if (Fortran) colstart_inc(); if (n_cc > nlcc && nlc < n_con && !(flags & ASL_no_linear_cc_rhs_adjust)) adjust_compl_rhs(); } static void #ifdef KR_headers br_read(R, nc, nc1, Lp, U, Cvar, nv) EdRead *R; real **Lp, *U; int nc, nc1, *Cvar, nv; #else br_read(EdRead *R, int nc, int nc1, real **Lp, real *U, int *Cvar, int nv) #endif { int i, inc, j, k; real *L; ASL *asl = R->asl; if (!(L = *Lp)) { if (nc1 < nc) nc1 = nc; L = *Lp = (real *)M1alloc(2*sizeof(real)*nc1); } if (U) inc = 1; else { U = L + 1; inc = 2; } 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 * #ifdef KR_headers aholread(R) EdRead *R; #else aholread(EdRead *R) #endif { int i, k; expr_h *rvh; char *s, *s1; FILE *nl = R->nl; s = read_line(R); k = *s; if (k < '1' || k > '9') badline(R); i = k - '0'; while((k = *++s) != ':') { 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; *s1 = *++s; s1++) if (--i < 0) badline(R); if (i) { for(*s1++ = '\n'; --i > 0; ) { k = getc(nl); if (k <= 0) badline(R); if ((*s1++ = k) == '\n') R->Line++; } if (getc(nl) != '\n') badline(R); } *s1 = 0; rvh->op = f_OPHOL; rvh->a = nv1; return (expr *)rvh; } static expr * #ifdef KR_headers bholread(R) EdRead *R; #else bholread(EdRead *R) #endif { 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 = 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; } int #ifdef KR_headers fgh_read_ASL(a, nl, flags) ASL *a; FILE *nl; int flags; #else fgh_read_ASL(ASL *a, FILE *nl, int flags) #endif { int i, j, k, maxfwd1, nc, nco, ncom, nlcon, nlin, no, nv, nvc, nvo, nz; int *ka, readall; unsigned x; expr_v *e; cgrad *cg, **cgp; ograd *og, **ogp; char fname[128]; func_info *fi; real t; EdRead ER, *R; Jmp_buf JB; ASL_CHECK(a, ASL_read_fgh, "fgh_read"); 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"); } if ((readall = flags & ASL_keep_all_suffixes) && a->i.nsuffixes) readall |= 1; func_add(a); if (binary_nl) { holread = bholread; xscanf = bscanf; } else { holread = aholread; xscanf = ascanf; } ncom = comb + comc + como + comc1 + como1; nc = n_con; no = n_obj; nvc = c_vars; nvo = o_vars; if (no < 0 || (nco = nc + no + nlcon) <= 0) scream(R, ASL_readerr_corrupt, "ed2read: nc = %d, no = %d, nlcon = %d\n", nc, no, nlcon); nv1 = nv0 = nvc > nvo ? nvc : nvo; max_var = nv = nv0 + ncom; combc = comb + comc; ncom0 = ncom_togo = combc + como; nzclim = ncom0 >> 3; ncom1 = comc1 + como1; nv0b = nv0 + comb; nv0c = nv0b + comc; nv01 = nv0 + ncom0; amax = lasta = lasta0 = lasta00 = nv1 + 1; ka = 0; if ((maxfwd1 = maxfwd + 1) > 1) nvref = maxfwd1*((ncom0 < vrefGulp ? ncom0 : vrefGulp) + 1); x = nco*sizeof(cde) + nc*sizeof(cgrad *) + 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.n_var_ - nvar0)) nvar0 += ncom0 + ncom1; if (pi0) { memset(pi0, 0, nc*sizeof(real)); if (havepi0) memset(havepi0, 0, nc); } if (X0) memset(X0, 0, nv0*sizeof(real)); if (havex0) memset(havex0, 0, nv0); 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; Cgrad = (cgrad **)(obj_de + no); Ograd = (ograd **)(Cgrad + nc); 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; 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(flags); nzjac = nz; if (!Lastx) Lastx = (real *)M1alloc(nv0*sizeof(real)); fclose(nl); nderps += nderp; a->p.Objval = obj2val_ASL; a->p.Objgrd = obj2grd_ASL; a->p.Conval = con2val_ASL; a->p.Jacval = jac2val_ASL; a->p.Hvcomp = hv2comp_ASL; a->p.Conival = con2ival_ASL; a->p.Congrd = con2grd_ASL; a->p.Lconval = lcon2val_ASL; a->p.Xknown = x2known_ASL; a->i.err_jmp_ = 0; return 0; } ER.can_end = 0; k = -1; switch(i) { case 'C': xscanf(R, "%d", &k); if (k < 0 || k >= nc) 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", &og->varno, &og->coef) != 2) badline(R); } *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 (ka) { if (!A_vals) goto cg_read; j += Fortran; while(k--) { if (xscanf(R, "%d %lf", &i, &t) != 2) badline(R); i = ka[i]++; A_vals[i] = t; A_rownos[i] = j; } break; } cg_read: cgp = Cgrad + j; j = 0; while(k--) { *cgp = cg = (cgrad *)mem(sizeof(cgrad)); cgp = &cg->next; if (ka) { if (xscanf(R, "%d %lf", &cg->varno, &cg->coef) != 2) badline(R); } else if (xscanf(R, "%d %d %lf", &cg->varno, &j, &cg->coef) != 3) badline(R); 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, nc, &LUrhs, Urhsx, cvar, nv0); break; case 'b': br_read(R, asl->i.n_var0, nv0, &LUv, Uvx, 0, 0); break; case 'k': k_seen++; k = asl->i.n_var0; if (!xscanf(R,"%d",&j) || j != k - 1) badline(R); if (!(ka = A_colstarts)) { if ((i = k) < n_var) i = n_var; ka = A_colstarts = (int *) M1alloc((i+1)*Sizeof(int)); } *ka++ = 0; *ka++ = 0; /* sic */ while(--k > 0) if (!xscanf(R, "%d", ka++)) badline(R); ka = A_colstarts + 1; break; case 'x': if (!xscanf(R,"%d",&k) || k < 0 || k > nv) badline(R); if (!X0 && want_xpi0 & 1) { x = nv0*sizeof(real); if (want_xpi0 & 4) x += nv0; X0 = (real *)M1zapalloc(x); if (want_xpi0 & 4) havex0 = (char*)(X0 + nv0); } 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 > nc) 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 >= nc) badline(R); if (pi0) { pi0[j] = t; if (havepi0) havepi0[j] = 1; } } break; default: badline(R); } } } #ifdef __cplusplus } #endif