/**************************************************************** Copyright (C) 1997-2001 Lucent Technologies All Rights Reserved Permission to use, copy, modify, and distribute this software and its documentation for any purpose and without fee is hereby granted, provided that the above copyright notice appear in all copies and that both that the copyright notice and this permission notice and warranty disclaimer appear in supporting documentation, and that the name of Lucent or any of its entities not be used in advertising or publicity pertaining to distribution of the software without specific, written prior permission. LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. ****************************************************************/ #include "nlp.h" #define Egulp 400 #define GAP_MAX 10 #ifdef Just_Linear #define who "f_read" #define fg_read_ASL f_read_ASL #else #define who "fg_read" #ifdef __cplusplus extern "C" { static real Missing_func(arglist*); static int compar(const void*, const void*, void*); } #endif /* __cplusplus */ #endif /* Just_Linear */ #undef nzc #undef r_ops typedef struct Static { int _k_seen, _nv0; ASL *a; ASL_fg *asl; efunc **_r_ops; #ifndef Just_Linear derp *_last_d; expr *(*_holread) ANSI((EdRead*)); expr_if *_iflist, *_if2list, *_if2list_end; expr_va *_varglist, *_varg2list, *_varg2list_end; relo *_relolist, *_relo2list; int *_imap, *_vrefnext, *_vrefx, *_zc, *_zci; int _amax1, _co_first, _firstc1, _imap_len; int _last_cex, _lasta, _lasta0, _lasta00, _lastc1, _lastj; int _max_var, _ncom_togo, _nderp, _nocopy; int _nv01, _nv011, _nv0b, _nv0c, _nv1, _nvref, _nzc, _nzclim; int nvar0, nvinc; #endif /* Just_Linear */ } Static; #define amax1 S->_amax1 #define co_first S->_co_first #define firstc1 S->_firstc1 #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_seen S->_k_seen #define last_cex S->_last_cex #define last_d S->_last_d #define lasta S->_lasta #define lasta0 S->_lasta0 #define lasta00 S->_lasta00 #define lastc1 S->_lastc1 #define lastj S->_lastj #define max_var S->_max_var #define ncom_togo S->_ncom_togo #define nderp S->_nderp #define nocopy S->_nocopy #define nv0 S->_nv0 #define nv01 S->_nv01 #define nv011 S->_nv011 #define nv0b S->_nv0b #define nv0c S->_nv0c #define nv1 S->_nv1 #define nvref S->_nvref #define nzc S->_nzc #define nzclim S->_nzclim #define r_ops S->_r_ops #define relo2list S->_relo2list #define relolist S->_relolist #define varg2list S->_varg2list #define varg2list_end S->_varg2list_end #define varglist S->_varglist #define vrefnext S->_vrefnext #define vrefx S->_vrefx #define zc S->_zc #define zci S->_zci #ifdef Just_Linear static void #ifdef KR_headers sorry_nonlin(R) EdRead *R; #else sorry_nonlin(EdRead *R) #endif { fprintf(Stderr, "Sorry, %s can't handle nonlinearities.\n", progname ? progname : ""); exit_ASL(R,ASL_readerr_nonlin); } #else /* Just_Linear */ #include "r_opn.hd" #endif /* Just_Linear */ 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); } static Static * #ifdef KR_headers ed_reset(S, a) Static *S; ASL *a; #else ed_reset(Static *S, ASL *a) #endif { memset(S, 0, sizeof(Static)); S->asl = (ASL_fg*)a; S->a = a; a->i.memLast = a->i.memNext = 0; #ifndef Just_Linear co_first = 1; #endif return S; } #ifdef Double_Align #define memadj(x) x #else #define memadj(x) (((x) + (sizeof(long)-1)) & ~(sizeof(long)-1)) #endif #ifndef Just_Linear static void #ifdef KR_headers fscream(R, name, nargs, kind) EdRead *R; char *name; 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 new_derp(S, a, b, c) Static *S; real *c; #else new_derp(Static *S, int a, int b, real *c) #endif { derp *d; if (a == nv1) return; 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; } static derp * #ifdef KR_headers new_relo(S, e, Dnext, ap) Static *S; expr *e; derp *Dnext; int *ap; #else new_relo(Static *S, expr *e, derp *Dnext, int *ap) #endif { 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); } r->D = r->Dcond = last_d; r->Dnext = Dnext; return r->D; } static relo * #ifdef KR_headers new_relo1(S, Dnext) Static *S; derp *Dnext; #else new_relo1(Static *S, derp *Dnext) #endif { relo *r; r = (relo *)mem_ASL(S->a, sizeof(relo)); r->next = relolist; relolist = r; r->D = 0; r->Dnext = Dnext; return r; } static expr * #ifdef KR_headers new_expr(S, opcode, L, R, deriv) Static *S; expr *L, *R; #else new_expr(Static *S, int opcode, expr *L, expr *R, int deriv) #endif { expr *rv; efunc *o; int L1, R1; 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_ASL(S->a, sizeof(expr)); rv->op = o; rv->L.e = L; rv->R.e = R; rv->a = nv1; if (deriv) { L1 = L && L->op != f_OPNUM && L->a != nv1; R1 = R && R->op != f_OPNUM && R->a != nv1; if (L1 | R1) { rv->a = lasta++; if (L1) new_derp(S, L->a, rv->a, &rv->dL); if (R1) new_derp(S, R->a, rv->a, &rv->dR); } } return rv; } #endif /* Just_Linear */ static expr * #ifdef KR_headers eread(R, deriv) EdRead *R; int deriv; #else eread(EdRead *R, int deriv) #endif { expr_n *rvn; short sh; fint L1; real r; Static *S = (Static *)R->S; ASL_fg *asl = S->asl; #ifndef Just_Linear char **sa; int a0, a1, *at, i, j, k, kd, ks, numargs, symargs; real *b, *ra; expr *L, *arg, **args, **argse, *rv; expr_va *rva; plterm *p; de *d; derp *dsave; efunc *op; expr_if *rvif; expr_f *rvf; func_info *fi; arglist *al; argpair *ap, *sap; char *dig; static real dvalue[] = { #include "dvalue.hd" }; #else /* Just_Linear */ Not_Used(deriv); #endif /* Just_Linear */ switch(edag_peek(R)) { #ifdef Just_Linear case 'f': case 'h': case 'o': case 'v': sorry_nonlin(R); #else 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; 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 == 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 && k) { kd = numargs; rvf->a = lasta++; } else { kd = 0; rvf->a = nv1; } ra = (real *)mem(sizeof(arglist) + (k+ks)*sizeof(argpair) + (numargs+kd)*sizeof(real) + symargs*sizeof(char *) + j*sizeof(int)); dig = 0; if (k < kd) dig = (char*)mem(numargs); b = ra + kd; al = rvf->al = (arglist *)(b + numargs); al->n = numargs + symargs; al->nr = numargs; al->ra = ra; if (kd) memset(al->derivs = b, 0, kd*sizeof(real)); else al->derivs = 0; 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); 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 = 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++; 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(S, arg->a, rvf->a, b); *b = 0; } } if (dig) *dig++ = j; b++; ra++; } } rvf->ape = ap; rvf->sape = sap; return (expr *)rvf; case 'h': return holread(R); #endif /* Just_Linear */ 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(size_expr_n); rvn->op = (efunc_n*)f_OPNUM; rvn->v = r; return (expr *)rvn; #ifndef Just_Linear 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); if (k < nv01 && deriv && !zc[k]++) zci[nzc++] = k; return (expr *)(var_e + k); #endif /* Just_Linear */ default: badline(R); } #ifndef Just_Linear if (xscanf(R, "%d", &k) != 1 || k < 0 || k >= N_OPS) badline(R); switch(optype[k]) { case 1: /* unary */ rv = new_expr(S, 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(S, 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 && deriv) { new_derp(S, lasta, lasta, &edagread_one); lasta++; } rva->d0 = dsave = last_d; a0 = a1 = lasta; for(j = 0; i > 0; i--, d++) { last_d = dsave; d->e = L = eread(R, deriv); if (L->op == f_OPNUM || L->a == nv1 || !deriv) { d->d = dsave; d->dv.i = nv1; } else { d->d = new_relo(S, L, dsave, &d->dv.i); j++; if (a1 < lasta) a1 = lasta; lasta = a0; } } d->e = 0; /* sentinnel expr * */ 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; } else { rva->a = nv1; rva->R.D = 0; } 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",&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); rv = (expr *)mem(sizeof(expr)); rv->op = f_OPPLTERM; rv->L.p = p; rv->R.e = L = eread(R, deriv); if (deriv) new_derp(S, L->a, rv->a = lasta++, &rv->dL); 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 && deriv) { new_derp(S, lasta, lasta, &edagread_one); lasta++; } rvif->d0 = dsave = last_d; rvif->e = eread(R, 0); a0 = lasta; rvif->T = L = eread(R, deriv); j = 0; if (L->op == f_OPNUM || L->a == nv1 || !(j = deriv)) { rvif->dT = dsave; rvif->Tv.i = nv1; } else rvif->dT = new_relo(S, L, dsave, &rvif->Tv.i); a1 = lasta; lasta = a0; last_d = dsave; rvif->F = L = eread(R, deriv); if (L->op == f_OPNUM || L->a == nv1 || !(j = deriv)) { rvif->dF = dsave; rvif->Fv.i = nv1; } else 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; } else { rvif->a = nv1; rvif->D = 0; } return (expr *)rvif; case 11: /* OPCOUNT */ deriv = 0; /* no break */ case 6: /* sumlist */ i = j = 0; xscanf(R, "%d", &i); if (i <= 2 && (optype[k] == 6 || i < 1)) badline(R); rv = (expr *)mem(sizeof(expr) - sizeof(real) + i*sizeof(expr *)); rv->op = r_ops[k]; a0 = lasta; rv->L.ep = args = (expr **)&rv->dR; if (deriv) { rv->a = lasta++; do { *args++ = L = eread(R, deriv); if (L->op != f_OPNUM && L->a != nv1) { new_derp(S, L->a, rv->a, &edagread_one); j++; } } while(--i > 0); } else do *args++ = eread(R, deriv); while(--i > 0); rv->R.ep = args; if (!j) { rv->a = nv1; lasta = a0; } return rv; } #endif /* Just_Linear */ badline(R); return 0; } #ifndef Just_Linear static list * #ifdef KR_headers new_list(asl, nxt) ASL *asl; list *nxt; #else new_list(ASL *asl, list *nxt) #endif { list *rv = (list *)mem(sizeof(list)); rv->next = nxt; return rv; } static list * #ifdef KR_headers crefs(S) Static *S; #else crefs(Static *S) #endif { int i; list *rv = 0; while(nzc > 0) { if ((i = zci[--nzc]) >= nv0) { rv = new_list(S->a, 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(S, d0, a, dnext) Static *S; derp *d0; int a; derp *dnext; #else derpadjust(Static *S, 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; ASL_fg *asl; if (!(d = d0)) return dnext; asl = S->asl; 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(S, ce, dnext) Static *S; cexp *ce; derp *dnext; #else derpcopy(Static *S, 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(S, map[d->a.i], map[d->b.i], d->c.rp); dprev = dprev->next = last_d; } dprev->next = dnext; return d00.next; } static void #ifdef KR_headers imap_alloc(S) Static *S; #else imap_alloc(Static *S) #endif { 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(S, alen, d, z) Static *S; cde *d; int alen, **z; #else comsubs(Static *S, int alen, cde *d, int **z) #endif { list *L; int a, i, j, k; int *r, *re, *z1; cexp *ce; derp *D, *dnext; relo *R; ASL_fg *asl = S->asl; 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 (z && (k = lastc1 - firstc1 + nzc)) { i = (2*k + 1)*sizeof(int); *z = z1 = k > 20 ? (int *)M1alloc(i) : (int *)mem(i); *z1++ = k; } if (nzc > 0) { R = new_relo1(S, 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++; } if (z) { *z1++ = j; *z1++ = a - 1; } dnext = R->D = derpcopy(S, ce, R->D); } while(++i < nzc); nzc = 0; } if (D || R) { if (!R) R = new_relo1(S, dnext); D = R->D = derpadjust(S, D, a, R->D); if (d->e->op != f_OPVARVAL) d->e->a = imap[d->e->a]; } if (z) for(i = firstc1 + nv01, j = lastc1 + nv01; i < j; i++) { *z1++ = i; *z1++ = imap[var_e[i].a]; } d->d = D; a += alen; d->zaplen = (a > lasta00 ? a - nv1 : 0)*sizeof(real); if (amax < a) amax = a; } #endif /* Just_Linear */ static void #ifdef KR_headers co_read(R, d, cexp1_end, k, z, wd) EdRead *R; cde *d; int *cexp1_end, k, wd, **z; #else co_read(EdRead *R, cde *d, int *cexp1_end, int k, int **z, int wd) #endif { #ifdef Just_Linear Not_Used(cexp1_end); Not_Used(z); #else Static *S = (Static *)R->S; ASL_fg *asl = S->asl; lastc1 = last_cex - nv011; if (cexp1_end) cexp1_end[k+1] = lastc1; if (amax1 < lasta) amax1 = lasta; if (co_first) { co_first = 0; if (imap_len < lasta) imap_alloc(S); f_b = funnelfix(f_b); f_c = funnelfix(f_c); f_o = funnelfix(f_o); } if (!lastj) { lasta = lasta0; last_d = 0; } lastj = 0; #endif /* Just_Linear */ d += k; d->e = eread(R, wd); #ifndef Just_Linear { int alen; alen = lasta - lasta0; if (imap_len < lasta) imap_alloc(S); if (z) { z += k; *z = 0; } comsubs(S, alen, d, z); firstc1 = lastc1; } #endif /* Just_Linear */ } #ifndef Just_Linear static linpart * #ifdef KR_headers linpt_read(R, nlin) EdRead *R; #else linpt_read(EdRead *R, int nlin) #endif { linpart *L, *rv; ASL *asl = R->asl; 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(S, ce, ip) Static *S; cexp *ce; int *ip; #else funnelkind(Static *S, cexp *ce, int *ip) #endif { int i, j, k, nzc0, rv; int *vr, *vre; ASL_fg *asl = S->asl; ce->vref = 0; if (!(nzc0 = nzc) || maxfwd <= 0) 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, nlin; #else cexp_read(EdRead *R, int k, int nlin) #endif { int a, fk, i, j, la0, na; int *z1, **zp; funnel *f, **fp; linpart *L, *Le; expr *e; cplist *cl, *cl0; cexp *ce; Static *S = (Static *)R->S; ASL_fg *asl = S->asl; ce = cexps + k - nv0; L = ce->L = linpt_read(R, ce->nlin = nlin); nocopy = 0; last_d = 0; ce->z.i = la0 = lasta; nderps += nderp; nderp = 0; e = ce->e = eread(R, want_derivs); if (la0 == lasta) { a = lasta++; if (e->op != f_OPNUM) new_derp(S, 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(S, i = L->v.i, a, &L->fac); if (!zc[i]++) zci[nzc++] = i; } if (zp = zaC) *(zp += k - nv0) = 0; if (fk = funnelkind(S, 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(S); if (fk == 1) { f->fulld = last_d; a = lasta00; z1 = 0; if (zp) { for(i = j = 0; i < nzc; i++) if (zci[i] >= nv0) j++; if (j) { i = (2*j + 1)*sizeof(int); *zp = z1 = j > 20 ? (int *)M1alloc(i) : (int *)mem(i); *z1++ = j; } } for(i = nzc; --i >= 0; ) if ((j = zci[i]) >= nv0) { if (z1) { *z1++ = j; *z1++ = a; } imap[var_e[j].a] = a++; } if ((na = ce->zlen) || a > lasta00) 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, ce->zlen, &f->fcde, zp); memcpy(zci, vrefx, i*sizeof(int)); } last_d = 0; cl0 = 0; while(i > 0) if ((a = var_e[zci[--i]].a) != nv1) { new_derp(S, 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->cref = crefs(S); --ncom_togo; } static void #ifdef KR_headers cexp1_read(R, j, k, nlin) EdRead *R; #else cexp1_read(EdRead *R, int j, int k, int nlin) #endif { Static *S = (Static *)R->S; ASL_fg *asl = S->asl; linpart *L, *Le; cexp1 *ce = cexps1 + (k - nv01); expr *e; int la0; L = ce->L = linpt_read(R, ce->nlin = nlin); if (!lastj) { last_d = 0; if (amax1 < lasta) amax1 = lasta; lasta = lasta0; lastj = j; } la0 = lasta; e = ce->e = eread(R, want_derivs); if (la0 == lasta) { j = lasta++; if (e->op != f_OPNUM) new_derp(S, e->a, j, &edagread_one); } else j = e->a; var_e[k].a = j; for(Le = L + nlin; L < Le; L++) new_derp(S, L->v.i, j, &L->fac); last_cex = k; } static void #ifdef KR_headers ifadjust(a, e) real *a; expr_if *e; #else ifadjust(real *a, expr_if *e) #endif { for(; e; e = e->next) { e->Tv.rp = &a[e->Tv.i]; e->Fv.rp = &a[e->Fv.i]; } } static void #ifdef KR_headers vargadjust(a, e) real *a; expr_va *e; #else vargadjust(real *a, expr_va *e) #endif { de *d; for(; e; e = e->next) { for(d = e->L.d; d->e; d++) d->dv.rp = &a[d->dv.i]; } } static void #ifdef KR_headers funneladj1(a, f) real *a; funnel *f; #else funneladj1(real *a, funnel *f) #endif { derp *d; cplist *cl; for(; 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 #ifdef KR_headers funneladjust(asl) ASL_fg *asl; #else funneladjust(ASL_fg *asl) #endif { cexp *c, *ce; linpart *L, *Le; real *a = adjoints; c = cexps; for(ce = c + ncom0; c < ce; c++) if (L = c->L) for(Le = L + c->nlin; L < Le; L++) L->v.rp = &var_e[L->v.i].v; funneladj1(a, f_b); funneladj1(a, f_c); funneladj1(a, f_o); } static void #ifdef KR_headers com1adjust(asl) ASL_fg *asl; #else com1adjust(ASL_fg *asl) #endif { cexp1 *c, *ce; linpart *L, *Le; expr_v *v = var_e; for(c = cexps1, ce = c + ncom1; c < ce; c++) for(L = c->L, Le = L + c->nlin; L < Le; L++) L->v.rp = &v[L->v.i].v; } #endif /* Just_Linear */ static void #ifdef KR_headers goff_comp(asl) ASL_fg *asl; #else goff_comp(ASL_fg *asl) #endif { 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 #ifdef KR_headers colstart_inc(S) Static *S; #else colstart_inc(Static *S) #endif { int *ka, *kae; ASL_fg *asl = S->asl; ka = A_colstarts; kae = ka + asl->i.n_var0; while(ka <= kae) ++*ka++; } #ifndef Just_Linear static void #ifdef KR_headers zerograd_chk(S) Static *S; #else zerograd_chk(Static *S) #endif { ASL_fg *asl = S->asl; 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; } } #endif /* Just_Linear */ static void #ifdef KR_headers adjust_compl_rhs(asl, opnum) ASL_fg *asl; efunc *opnum; #else adjust_compl_rhs(ASL_fg *asl, efunc *opnum) #endif { 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 == 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(S, flags) Static *S; int flags; #else adjust(Static *S, int flags) #endif { ASL_fg *asl = S->asl; #ifndef Just_Linear 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(adjoints, iflist); vargadjust(adjoints, varglist); if (ncom0) funneladjust(asl); com1adjust(asl); if (n_obj) zerograd_chk(S); #endif /* Just_Linear */ if (k_seen) { if (!A_vals) goff_comp(asl); else if (Fortran) colstart_inc(S); } if (n_cc > nlcc && nlc < n_con && !(flags & ASL_no_linear_cc_rhs_adjust)) adjust_compl_rhs(asl, f_OPNUM); } 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); } } } #ifndef Just_Linear static expr * #ifdef KR_headers aholread(R) EdRead *R; #else aholread(EdRead *R) #endif { 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 = 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; Static *S = (Static *)R->S; ASL_fg *asl = S->asl; 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; } static void #ifdef KR_headers nlvzap(S, i, j) Static *S; int i, j; #else nlvzap(Static *S, int i, int j) #endif { int n = nv1; expr_v *v = S->var_e; i -= j; while(--j >= 0) v[i+j].a = n; } static real #ifdef KR_headers Missing_func(al) arglist *al; #else Missing_func(arglist *al) #endif { AmplExports *ae = al->AE; func_info *fi = (func_info*)al->funcinfo; char *s = (char*)(*ae->Tempmem)(al->TMI, strlen(fi->name) + 64); (*ae->SprintF)(al->Errmsg = s, "Attempt to call unavailable function %s.", fi->name); return 0.; } #endif /* Just_Linear */ int #ifdef KR_headers fg_read_ASL(a, nl, flags) ASL *a; FILE *nl; int flags; #else fg_read_ASL(ASL *a, FILE *nl, int flags) #endif { cgrad *cg, **cgp; expr_v *e; int *ka; int i, i1, j, k, nc, nco, nlcon, no, nv, nvc, nvo, nz, readall; ograd *og, **ogp; real t; unsigned x; #ifdef Just_Linear #define ASL_readtype ASL_read_f #else /* Just_Linear */ #define ASL_readtype ASL_read_fg int maxfwd1, ncom, nlin; func_info *fi; char fname[128]; #endif /* Just_Linear */ EdRead ER, *R; Static SS, *S; Jmp_buf JB; ASL_CHECK(a, ASL_readtype, who); S = ed_reset(&SS, a); 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; } } 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; #define asl ((ASL_fg*)a) ER.lineinc = 1; if (!size_expr_n) size_expr_n = sizeof(expr_n); #ifdef Just_Linear /* Tried "xscanf = binary ? bscanf : ascanf;", but this * encounters a bug in SGI's C compiler. */ if (binary_nl) xscanf = bscanf; else xscanf = ascanf; #else if (!(SS._r_ops = asl->I.r_ops_)) SS._r_ops = r_ops_ASL; if (c_cexp1st) *c_cexp1st = 0; if (o_cexp1st) *o_cexp1st = comc1; if (nfunc) func_add(a); if (binary_nl) { holread = bholread; xscanf = bscanf; } else { holread = aholread; xscanf = ascanf; } ncom = comb + comc + como + comc1 + como1; #endif /* Just_Linear */ 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, "edagread: nc = %d, no = %d, nlcon = %d\n", nc, no, nlcon); if (pi0) { memset(pi0, 0, nc*sizeof(real)); if (havepi0) memset(havepi0, 0, nc); } #ifdef Just_Linear nv = nv0 = nvc > nvo ? nvc : nvo; x = nco*sizeof(cde) + no*sizeof(ograd *) + nv*sizeof(expr_v) + no; #else 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; last_cex = nv011 = nv01 - 1; amax = lasta = lasta0 = lasta00 = nv1 + 1; 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; SS.nvar0 = a->i.n_var0; if (!(SS.nvinc = a->i.n_var_ - SS.nvar0)) SS.nvar0 += ncom0 + ncom1; #endif /* Just_Linear */ if (X0) memset(X0, 0, nv0*sizeof(real)); if (havex0) memset(havex0, 0, nv0); e = var_e = (expr_v *)M1zapalloc(x); con_de = (cde *)(e + nv); lcon_de = con_de + nc; obj_de = lcon_de + nlcon; Ograd = (ograd **)(obj_de + no); if (A_vals) { if (!A_rownos) A_rownos = (int *)M1alloc(asl->i.nzc_*sizeof(int)); } else if (nc) Cgrad = (cgrad **)M1zapalloc(nc*sizeof(cgrad *)); #ifdef Just_Linear objtype = (char *)(Ograd + no); #else var_ex = e + nv0; var_ex1 = var_ex + ncom0; for(k = 0; k < nv; e++) { e->op = f_OPVARVAL; e->a = k++; } if (skip_int_derivs) { if (nlvbi) nlvzap(S, nlvb, nlvbi); if (nlvci) nlvzap(S, nlvb+nlvc, nlvci); if (nlvoi) nlvzap(S, nlvb+nlvc+nlvo, nlvoi); } cexps = (cexp *)(Ograd + no); cexps1 = (cexp1 *)(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; } last_d = 0; #endif if (n_cc && !cvar) cvar = (int*)M1alloc(nc*sizeof(int)); if (cvar) memset(cvar, 0, nc*sizeof(int)); ka = 0; nz = 0; for(;;) { ER.can_end = 1; i = edag_peek(R); if (i == EOF) { #ifndef Just_Linear free(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 += nv1 + 1) > amax) amax = i; } adjoints = (real *)M1zapalloc(amax*Sizeof(real)); adjoints_nv1 = &adjoints[nv1]; nderps += nderp; #endif /* Just_Linear */ adjust(S, flags); nzjac = nz; if (!Lastx) Lastx = (real *)M1alloc(nv0*sizeof(real)); fclose(nl); #ifndef Just_Linear a->p.Objval = obj1val_ASL; a->p.Objgrd = obj1grd_ASL; a->p.Conval = con1val_ASL; a->p.Jacval = jac1val_ASL; a->p.Conival = con1ival_ASL; a->p.Congrd = con1grd_ASL; a->p.Lconval = lcon1val_ASL; a->p.Xknown = x1known_ASL; #endif /* Just_Linear */ 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,c_cexp1st,k,zac,want_derivs); break; #ifdef Just_Linear case 'F': case 'V': sorry_nonlin(R); #else 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))){ if (!(flags & ASL_allow_missing_funcs)) scream(R, ASL_readerr_unavail, "function %s not available\n", fname); fi->funcp = Missing_func; fi->funcinfo = (Char*)fi; } funcs[i] = fi; break; case 'L': xscanf(R, "%d", &k); if (k < 0 || k >= nlcon) badline(R); co_read(R, lcon_de, 0, k, 0, 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 < nv0 || k >= nv) badline(R); if (j) cexp1_read(R, j, k, nlin); else cexp_read(R, k, nlin); break; #endif /* Just_Linear */ 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); i1 = ka[i]++; A_vals[i1] = t; A_rownos[i1] = 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 'O': if (xscanf(R, "%d %d", &k, &j) != 2 || k < 0 || k >= no) badline(R); objtype[k] = j; co_read(R,obj_de,o_cexp1st,k,zao,want_derivs); 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 > nv0) 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 || j < 0 || j >= nv) 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); } } }