/**************************************************************** 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. ****************************************************************/ #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 #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; #define NFHASH 23 #undef f_OPNUM #include "r_opn0.hd" #undef nzc /*******/ 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; Char **mp; } Elemtemp; typedef struct PSfind { ps_func *f; Elemtemp *b, *g; } PSfind; static int nlthash = 1021; /* hash table length (fixed for now) */ static int nrangehash = 1021; /* range hash table length (fixed for now) */ 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; expr_n *_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; 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 _nv0; int _nv0b; int _nv0c; 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 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 nsce S->_nsce #define nv0 S->_nv0 #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 rangehash S->_rangehash #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 * #ifdef KR_headers S_init(S, asl) Static *S; ASLTYPE *asl; #else S_init(Static *S, ASLTYPE *asl) #endif { memset(S, 0, sizeof(Static)); S->asl = asl; S->a = (ASL*)asl; return S; } 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 void #ifdef KR_headers ed_reset(asl) ASLTYPE *asl; #else ed_reset(ASLTYPE *asl) #endif { asl->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 expr_n *enzork; static int izork = -1; #endif static Elemtemp * #ifdef KR_headers new_Elemtemp(S, esize, mp) Static *S; unsigned int esize; Char **mp; #else new_Elemtemp(Static *S, unsigned int esize, Char **mp) #endif { 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(Char*) << k) / esize; return e; } static void #ifdef KR_headers del_Elemtemp(S, e) Static *S; Elemtemp *e; #else del_Elemtemp(Static *S, Elemtemp *e) #endif { ASL *asl = S->a; del_mblk(e->k, *e->mp); del_mblk(k_Elemtemp, e); } static void #ifdef KR_headers upgrade_Elemtemp(S, e) Static *S; Elemtemp *e; #else upgrade_Elemtemp(Static *S, Elemtemp *e) #endif { Char *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(Char*) << k) / e->esize; } static void #ifdef KR_headers free_laref(S, L) Static *S; la_ref **L; #else free_laref(Static *S, la_ref **L) #endif { 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 * #ifdef KR_headers new_laref(S, nxt) Static *S; la_ref *nxt; #else new_laref(Static *S, la_ref *nxt) #endif { 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 * #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 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; int a, b; real *c; #else new_derp(Static *S, int a, int b, real *c) #endif { 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(""); #endif } 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); } if (!last_d) return 0; 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 void #ifdef KR_headers efree(S, e) Static *S; expr *e; #else efree(Static *S, expr *e) #endif { expr **args, **argse, *e1; top: switch(optype[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: *(expr_n **)&((expr_n*)e)->v = expr_n_free; expr_n_free = (expr_n*)e; } } static expr * #ifdef KR_headers new_expr(S, o, L, R) Static *S; int o; expr *L, *R; #else new_expr(Static *S, int o, expr *L, expr *R) #endif { expr *rv; 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; else if (Intcast L->op == f_OPNUM) o = f_OPCPOW; rv->op = (efunc *)(long)o; rv->L.e = L; rv->R.e = R; #ifdef DEBUG if (rv == ezork) printf(""); #endif return rv; } static void #ifdef KR_headers dexpr(S, e, L, R) Static *S; expr *e; expr *L, *R; #else dexpr(Static *S, expr *e, expr *L, expr *R) #endif { 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; 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; } e->dO.i = i; #endif } } static real dvalue[] = { #include "dvalue.hd" }; static expr_n * #ifdef KR_headers new_expr_n(S, t) Static *S; real t; #else new_expr_n(Static *S, real t) #endif { expr_n *en; if (en = expr_n_free) expr_n_free = *(expr_n **)&en->v; else en = (expr_n *)mem_ASL(S->a, S->size_exprn); #ifdef DEBUG if (en == enzork) printf(""); #endif en->v = t; en->op = (efunc_n *)f_OPNUM; return en; } static expr * #ifdef KR_headers eread(R) EdRead *R; #else eread(EdRead *R) #endif { short sh; fint L1; real r; char **sa; int *at, i, j, k, ks, numargs, op, symargs; real *b, *ra; expr *L, *arg, **args, **argse, *rv; expr_va *rva; plterm *p; de *d; expr_if *rvif; expr_f *rvf; func_info *fi; arglist *al; argpair *ap, *sap; Static *S = (Static *)R->S; ASLTYPE *asl = S->asl; 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", &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: 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, "%d", &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 *)(long)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",&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 = (efunc *)(long)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;) rvif->op = (efunc *)(long)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 #ifdef KR_headers ogfree(S, og) Static *S; ograd *og; #else ogfree(Static *S, ograd *og) #endif { 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 #ifdef KR_headers ogfree1(S, ogp) Static *S; ograd **ogp; #else ogfree1(Static *S, ograd **ogp) #endif { ograd *og = *ogp; *ogp = og->next; og->next = freeog; freeog = og; return og->coef; } static ograd * #ifdef KR_headers new_ograd(S, next, varno, coef) Static *S; ograd *next; int varno; real coef; #else new_ograd(Static *S, ograd *next, int varno, real coef) #endif { 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 * #ifdef KR_headers lahash(S, la) Static *S; linarg *la; #else lahash(Static *S, linarg *la) #endif { ograd *og, *og1; linarg *la1, **lap; int nnz = la->nnz; Ulong x = nnz; for(og = la->nz; og; og = og->next) x = (x << 1 | (x & 0x80000000) >> 31) ^ og->varno*101 + ((Ulong *)&og->coef)[0] + ((Ulong *)&og->coef)[1]; for(lap = <hash[x % nlthash]; 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; } } S->asl->P.nlttot++; return *lap = la; } static range * #ifdef KR_headers new_range(asl, r, linkup) ASLTYPE *asl; range *r; int linkup; #else new_range(ASLTYPE *asl, range *r, int linkup) #endif { int len, uilen; range *r1, *r2; uilen = r->nv*sizeof(int); len = sizeof(range) + uilen; if (r->n >= r->nv) r1 = (range*)(linkup ? mem(len) : new_mblk(htcl(len))); else r1 = (range*)mem(len); r1->nintv = 0; r1->n = r->n; r1->nv = r->nv; 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); if (linkup) { r2 = r1->rlist.next = asl->P.rlist.next; r1->rlist.prev = r2->rlist.prev; r2->rlist.prev = asl->P.rlist.next = r1; } else r1->rlist.prev = r1->rlist.next = 0; return r1; } static int #ifdef KR_headers lacompar(a, b, v) char *a, *b, *v; #else lacompar(const void *a, const void *b, void *v) #endif { 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 #ifdef KR_headers ndiff(r, r1, nn1p) range *r, *r1; int *nn1p; #else ndiff(range *r, range *r1, int *nn1p) #endif { /* Return the number of linargs in r but not r1; */ /* set *nn1p to the number of linargs in r1 but not r. */ int i, nn, nn1; linarg **la, **la1, **la1e, **lae; nn = nn1 = 0; la = r->lap; lae = la + r->n; la1 = r1->lap; la1e = la1 + r1->n; for(;;) { if (la >= lae) { nn1 += la1e - la1; break; } if (la1 >= la1e) { nn += lae - la; break; } if (!(i = lacompar((char*)la, (char *)la1, NULL))) { la++; la1++; } else if (i < 0) { nn++; la++; } else { nn1++; la1++; } } if (nn1p) *nn1p = nn1; return nn; } static void #ifdef KR_headers lap_merge(asl, nn, r, r1) ASL *asl; int nn; range *r, *r1; #else lap_merge(ASL *asl, int nn, range *r, range *r1) #endif { /* merge nn new linargs from r into r1 */ int i, n1; linarg **la, **la1, **la1e, **lae, **lap, **lap0; la = r->lap; lae = la + r->n; la1 = r1->lap; la1e = la1 + r1->n; n1 = r1->n + nn; lap = lap0 = (linarg**)new_mblk(htcl(n1*sizeof(linarg*))); for(;;) { if (la >= lae) { while(la1 < la1e) *lap++ = *la1++; break; } if (la1 >= la1e) { do *lap++ = *la++; while(la < lae); break; } if (!(i = lacompar((char*)la, (char *)la1, NULL))) { *lap++ = *la++; la1++; } else if (i < 0) *lap++ = *la++; else *lap++ = *la1++; } del_mblk(htcl(r1->n*sizeof(linarg*)), r1->lap); r1->lap = lap0; r1->n = n1; } static void #ifdef KR_headers ucopy(asl, r, rp) ASL *asl; range *r, **rp; #else ucopy(ASL *asl, range *r, range **rp) #endif { range *r1; int nn; r1 = *rp; /* see whether to update this copy */ if (r->n == r1->n && !memcmp(r->lap, r1->lap, r->n*sizeof(linarg*))) return; if (!(nn = ndiff(r, r1, 0))) return; lap_merge(asl, nn, r, r1); } static range ** #ifdef KR_headers uhash(S, r) Static *S; range *r; #else uhash(Static *S, range *r) #endif { range *r1, **rp; int len, n, nv, *ui, *uie; unsigned long L = 0; ASLTYPE *asl = S->asl; nv = r->nv; ui = r->ui; uie = ui + nv; len = sizeof(int)*nv; while (ui < uie) L = 37*L + *ui++; ui = r->ui; rp = rangehash + L%nrangehash; n = r->n; if (S->asl->P.merge) while(r1 = *rp) { if (r1->nv == nv && r1->n == n && !memcmp(ui, r1->ui, len)) { ucopy((ASL*)asl, r, rp); return rp; } rp = &r1->hunext; } *rp = new_range(asl, r, 1); return rp; } static range ** #ifdef KR_headers rhash(S, r, addnew) Static *S; range *r; int addnew; #else rhash(Static *S, range *r, int addnew) #endif { range *r1, **rp; ograd *og; linarg **lae, **lap; int len, n; ASLTYPE *asl = S->asl; unsigned long 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; } lap = r->lap; rp = rangehash + L%nrangehash; if (asl->P.merge) while(r1 = *rp) { if (r1->n == n && !memcmp(lap, r1->lap, len)) return rp; rp = &r1->hnext; } if (addnew) *rp = new_range(asl, r, 1); return rp; } static int #ifdef KR_headers compar(a0, b0, v) char *a0, *b0, *v; #else compar(const void *a0, const void *b0, void *v) #endif { int a, b, c; Static *S = (Static *)v; if ((a = *(int*)a0) >= max_var) { a = zl[a-nv0]; if ((b = *(int*)b0) >= max_var) { if (c = a - zl[b-nv0]) return c; return *(int*)a0 - b; } if (a == b) return -1; } else if ((b = *(int*)b0) >= max_var) { b = zl[b-nv0]; if (a == b) return 1; } return a - b; } static int #ifdef KR_headers hscompar(a0, b0, v) char *a0, *b0, *v; #else hscompar(const void *a0, const void *b0, void *v) #endif { 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 -= nv0; if (a == b) return -1; } else if ((b = *(int*)b0) >= Ncom) { b = zl[b] - nv0; if (a == b) return 1; } return a - b; } static void #ifdef KR_headers zcsort(S, c, ci, i, n, p) Static *S; int *c, *ci, i, n, p; #else zcsort(Static *S, int *c, int *ci, int i, int n, int p) #endif { 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* #ifdef KR_headers compress(S, og, cp, comvar) Static *S; ograd *og; real *cp; int *comvar; #else compress(Static *S, ograd *og, real *cp, int *comvar) #endif { 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 < nv0) { *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]) < nv0) i++; else { if (!zc[j]++) zci[nzc++] = j; t = rnz[j]; j1 = j - nv0; 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 * #ifdef KR_headers afree(S, og, ep) Static *S; ograd *og; expr **ep; #else afree(Static *S, ograd *og, expr **ep) #endif { 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 * #ifdef KR_headers af_sum(S, Log, Rog) Static *S; ograd *Log, *Rog; #else af_sum(Static *S, ograd *Log, ograd *Rog) #endif { 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 * #ifdef KR_headers cf_sum(S, Lcg, Rog) Static *S; cgrad *Lcg; ograd *Rog; #else cf_sum(Static *S, cgrad *Lcg, ograd *Rog) #endif { 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 #ifdef KR_headers sumlist_afree(S, Laf, e, argso, sls, slscr) Static *S; ograd *Laf; expr *e; expr **argso; expr **sls, **slscr; #else sumlist_afree(Static *S, ograd *Laf, expr *e, expr **argso, expr **sls, expr **slscr) #endif { 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 #ifdef KR_headers sdvcite(S, k) Static *S; int k; #else sdvcite(Static *S, int k) #endif { 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 * #ifdef KR_headers awalk(S, e) Static *S; expr *e; #else awalk(Static *S, expr *e) /* return 0 if e is not linear */ #endif { 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(optype[k = Intcast e->op]) { case 1: /* unary */ if (k == f_OPCPOW) { if (Raf = awalk(S, e->R.e)) afree(S, Raf, &e->R.e); break; } 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->e), &rvif->e); /* The above may introduce "range" rows that */ /* are not involved in any derivatives... */ 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; asl = 0; if (!slscratchlev++) slscr = slscratch; else { asl = S->asl; 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 (asl) 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 */ {ASLTYPE *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 - nv0)) < 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 */ args = e->L.ep; for(argse = e->R.ep; args < argse; args++) afree(S, awalk(S, *args), args); break; default: scream(S->R, ASL_readerr_bug, "awalk: unexpected optype[%d] = %d\n", k, optype[k]); } return 0; } static void #ifdef KR_headers cexp_upgrade(S, t) Static *S; int t; #else cexp_upgrade(Static *S, int t) #endif { 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 #ifdef KR_headers zc_upgrade(S) Static *S; #else zc_upgrade(Static *S) #endif { 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 #ifdef KR_headers la_replace(S, la) Static *S; linarg *la; #else la_replace(Static *S, linarg *la) #endif { 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 #ifdef KR_headers tlistgen(S, f) Static *S; ps_func *f; #else tlistgen(Static *S, ps_func *f) #endif { 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 = nv0 + *++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 #ifdef KR_headers might_expand(S, e) Static *S; expr *e; #else might_expand(Static *S, expr *e) #endif { 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 >= nv0) return 1; } return 0; } static void #ifdef KR_headers ce_split(S, i, f) Static *S; int i; ps_func *f; #else ce_split(Static *S, int i, ps_func *f) #endif { 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 += nv0; j = k + nv0; n += k; do { zl[k++] = i; zci[nzc++] = j++; /* for crefs */ } while(k < n); } static void #ifdef KR_headers linpart_augment(S, c, og, f) Static *S; cexp *c; ograd *og; ps_func *f; #else linpart_augment(Static *S, cexp *c, ograd *og, ps_func *f) #endif { 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 * #ifdef KR_headers linterms(S, c, scale) Static *S; cexp *c; real scale; #else linterms(Static *S, cexp *c, real scale) #endif { 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 #ifdef KR_headers dvwalk(S, i) Static *S; int i; #else dvwalk(Static *S, int i) #endif { int n; ograd *og, *og0; linarg *tl, **tp; ps_func f; dv_info *dvi; cexp *c; int split; 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 og = awalk(S, c->e); 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 (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 #ifdef KR_headers colindvref(S, e, ndv) Static *S; expr *e; int ndv; #else colindvref(Static *S, expr *e, int ndv) #endif { 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 -= nv0) < 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 #ifdef KR_headers dv_walk(S) Static *S; #else dv_walk(Static *S) #endif { int i, m; expr_v *v; ASLTYPE *asl = S->asl; m = asl->i.n_con0; if (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 < Ncom; i++) dvwalk(S, i); *larvlist = 0; if (i = asl->P.ndvspout) nv0b = nv0 + Ncom + i; larvlist = 0; } } static int /* adjust zci, return k s.t. zci[i] < nv0 for i < k */ #ifdef KR_headers nzcperm(S) Static *S; #else nzcperm(Static *S) #endif { int i, j, k; k = nzc; for(i = 0; i < k; ) if (zci[i] >= nv0) { j = zci[--k]; zci[k] = zci[i]; zci[i] = j; } else i++; return k; } static void #ifdef KR_headers sumlist_adj(asl, e, e1) ASLTYPE *asl; expr *e, *e1; #else sumlist_adj(ASLTYPE *asl, expr *e, expr *e1) #endif { 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 #ifdef KR_headers termwalk(S, ep, p) Static *S; expr **ep; PSfind *p; #else termwalk(Static *S, expr **ep, PSfind *p) #endif { 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-nv0].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 + nv0]++) zci[nzc++] = j; } } r = (range *)rnz; /* scratch */ if (ncp = nzc - k) { zcsort(S, zc, zci+k, nv0, ncp, -1 /*max_var*/); cp = cp0 = (int*)(r+1); i = k; *cp = ncp; do { zc[j = zci[i++]] = 0; *++cp = j - nv0; } 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, nv0); 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 #ifdef KR_headers co_finish(f) ps_func *f; #else co_finish(ps_func *f) #endif { 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 * #ifdef KR_headers ecopy(S, e) Static *S; expr *e; #else ecopy(Static *S, expr *e) #endif { 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 * #ifdef KR_headers ltermwalk(S, scale, ep, p) Static *S; real scale; expr **ep; PSfind *p; #else ltermwalk(Static *S, real scale, expr **ep, PSfind *p) #endif { dv_info *dvi; expr **a, **ae, *e; expr_n *en; int k, k0; ograd *og, *rv; cexp *ce; ASLTYPE *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 < nv0) { rv = af_sum(S, rv, new_ograd(S,0,k,scale)); break; } k0 = k; if ((k -= nv0) >= 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 = (expr_n *)e->R.e; if (Intcast en->op == f_OPNUM) { *ep = e->L.e; case_opnum: if (en->v == 0.) { efree(S,*ep); *ep = (expr*)en; } else { rv = af_sum(S, rv, ltermwalk(S, scale*en->v, ep, p)); *(expr_n **)&en->v = expr_n_free; expr_n_free = en; } e->L.e = expr_free; expr_free = e; break; } en = (expr_n *)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 #ifdef KR_headers PSfind_init(S, f, psf, wantg) Static *S; ps_func *f; PSfind *psf; int wantg; #else PSfind_init(Static *S, ps_func *f, PSfind *psf, int wantg) #endif { 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 #ifdef KR_headers psb_copy(b, b0, n) psb_elem *b; psb_elem *b0; int n; #else psb_copy(psb_elem *b, psb_elem *b0, int n) #endif { 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 #ifdef KR_headers getgroup(S, scale, e, p) Static *S; real scale; expr *e; PSfind *p; #else getgroup(Static *S, real scale, expr *e, PSfind *p) #endif { 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, nv0); og = 0; while(nzc1 > 0) { og = new_ograd(S, og, i = zci1[--nzc1], 0.); zc1[i] = 0; } g->og = og; return 1; } static ograd * #ifdef KR_headers cotermwalk(S, ep, f, wantg, omitdv) Static *S; expr **ep; ps_func *f; int wantg, omitdv; #else cotermwalk(Static *S, expr **ep, ps_func *f, int wantg, int omitdv) #endif { 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 #ifdef KR_headers psfind(S) Static *S; #else psfind(Static *S) #endif { int i, j, m; fint x; ps_func *f, *f1; range *r, *r0; linarg *la, **lap; #ifdef PSHVREAD int ihdlim, ihdlim1, ihdmax, ihdmin, n; Ihinfo *ihi, *ihi0, *ihi1; #endif ASLTYPE *asl = S->asl; m = asl->i.n_con0; x = (n_obj+m)*sizeof(ps_func) + nlthash*sizeof(linarg*) + nrangehash*sizeof(range*) + slmax*sizeof(expr*) + (Ncom+1)*sizeof(int); asl->P.ops = f = (ps_func *)M1alloc(x); asl->P.cps = f1 = f + n_obj; lthash = (linarg**)(f1 + m); rangehash = (range**)(lthash + nlthash); slscratch = (expr**)(rangehash + nrangehash); memset(lthash, 0, (nlthash+nrangehash)*sizeof(linarg*)); asl->P.dvsp0 = (int*)(slscratch + slmax); *asl->P.dvsp0 = Ncom; /* 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 = 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 } /*******/ static void #ifdef KR_headers ewalk(S, e, deriv) Static *S; expr *e; int deriv; #else ewalk(Static *S, expr *e, int deriv) #endif { int a0, a1, i, j, k, kf, numargs, op; real *b, *b0, *ra; unsigned int len; #ifdef PSHVREAD ASL *asl; argpair *da; int i1, j0, j1, 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(optype[k = Intcast e->op]) { case 1: /* unary */ e->dL = dvalue[k]; /* for UMINUS, FLOOR, CEIL */ if (k == f_OPCPOW) { ewalk(S, e->R.e, deriv); if (deriv) dexpr(S, e, 0, e->R.e); break; } 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 *)(long)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); } } 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; nn = 0; } rvf->da = da; rvf->fh = fh; nn0 = nn; #endif for(ap = ap0 = rvf->ap; ap < ape; ap++) { j = 1; arg = ap->e; arg->op = r_ops[op = Intcast arg->op]; switch(op) { case f_OPNUM: goto loopend; case f_OPVARVAL: arg->op = (efunc *)(long)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; j = 0; #endif loopend: if (dig) *dig++ = j; } #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 * #ifdef KR_headers crefs(S) Static *S; #else crefs(Static *S) #endif { int i, maxvar1 = S->max_var1; list *rv = 0; while(nzc > 0) { if ((i = zci[--nzc]) >= nv0 && i < maxvar1) { 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; 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 * #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 *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 * #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 derp * #ifdef KR_headers derp_ogcopy(S, og, dnext, k) Static *S; ograd *og; derp *dnext; int k; #else derp_ogcopy(Static *S, ograd *og, derp *dnext, int k) #endif { 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 #ifdef KR_headers imap_alloc(S) Static *S; #else imap_alloc(Static *S) #endif { 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 = nv0; while(i > 0) *--r = --i; i = nv0; for(tl = asl->P.lalist; tl; tl = tl->lnext) if (v = tl->v) r[v->a] = i++; r[noa] = i; } static void #ifdef KR_headers comsubs(S, alen, d) Static *S; cde *d; int alen; #else comsubs(Static *S, int alen, cde *d) #endif { list *L; int a, i, j, k; int *r, *re; cexp *ce; derp *D, *dnext; relo *R; expr *e; split_ce *cs; ograd *og; ASLTYPE *asl = S->asl; D = last_d; a = lasta00; dnext = 0; R = 0; for(i = j = 0; i < nzc; i++) if ((k = zci[i]) >= nv0 && k < max_var1) zci[j++] = k; else zc[k] = 0; if (nzc = j) { for(i = 0; i < nzc; i++) { if ((j = zci[i] - nv0) >= Ncom) { cs = asl->P.Split_ce + (j-Ncom); if (r = cs->ce) { re = r + *r; while(r++ < re) if (!zc[j = *r + nv0]++) zci[nzc++] = j; } } else if (j >= 0 && (L = cexps[j].cref)) { if (cexps[j].funneled) do { if (zc[j = L->item.i] || asl->P.dv[j-nv0].ll) continue; zc[j] = 1; zci[nzc++] = j; } while(L = L->next); else do { if (!zc[L->item.i]++) zci[nzc++] = L->item.i; } while(L = L->next); } } if (nzc > 1) zcsort(S, zc, zci, 0, nzc, -1); } if (nzc > 0) { R = new_relo1(S, dnext); i = 0; do { j = zci[i]; zc[j] = 0; ce = &cexps[j - nv0]; e = ce->e; k = varp[j-nv0]->a; if (ce->funneled) { if (j >= max_var) imap[((expr_vx*)varp[j-nv0])->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 && (k = j - nv0) < Ncom && (og = asl->P.dv[k].ll)) { dnext = R->D = derp_ogcopy(S,og,R->D,imap[j]); continue; } } 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 (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 #ifdef KR_headers co_read(R, d, k) EdRead *R; cde *d; int k; #else co_read(EdRead *R, cde *d, int k) #endif { d = (cde *)((char *)d + k*sizeof(cde)); d->e = eread(R); } static void #ifdef KR_headers co_walk(S, d) Static *S; cde *d; #else co_walk(Static *S, cde *d) #endif { 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 #ifdef KR_headers lpcompar(a, b, v) char *a, *b, *v; #else lpcompar(const void *a, const void *b, void *v) #endif { Not_Used(v); return ((linpart *)a)->v.i - ((linpart *)b)->v.i; } static linpart * #ifdef KR_headers linpt_read(R, nlin) EdRead *R; int nlin; #else linpt_read(EdRead *R, int nlin) #endif { linpart *L, *rv; int i0, needsort; ASL *asl = R->asl; if (nlin <= 0) return 0; 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 #ifdef KR_headers funnelkind(S, i0, ip) Static *S; int i0, *ip; #else funnelkind(Static *S, int i0, int *ip) #endif { cexp *ce, *cej; dv_info *dvi; int i, j, k, k1, k2, nzc0, rv; int *vr, *vre; linarg *la, **lap, **lape; 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; for(i = k = 0; i < nzc; i++) if ((j = zci[i]) < nv0) { if (k >= maxfwd) { k = 0; break; } vrefx[k++] = j; } else { cej = cexps + (j -= nv0); if (!(vr = cej->vref)) { if (dvi && j < Ncom && dvi[j].ll) for(tl = cej->cref; tl; tl = tl->next) if (!zc[tl->item.i]++) zci[nzc++] = tl->item.i; 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++]; 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)) { *ip = k; return 2; } } } rv = 0; if (nocopy || nderp > 3*(nzc0+k1)) 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 { cexp *ce; expr *e, *zero; Static *S = (Static *)R->S; ASLTYPE *asl = S->asl; ce = cexps + k - nv0; 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 * #ifdef KR_headers funnelderp(S, a, cl0) Static *S; int a; cplist *cl0; #else funnelderp(Static *S, int a, cplist *cl0) #endif { 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 #ifdef KR_headers cexp_walk(S, k) Static *S; int k; #else cexp_walk(Static *S, int k) #endif { 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 + nv0; if (ce->zlen = lasta - la0) ce->z.i = la0; else { ce->z.i = k < Ncom ? ka : ((expr_vx*)varp[k-Ncom])->a0; ce->zlen = 1; } 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; } 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]) >= nv0 && j < max_var1) imap[varp[j-nv0]->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); memcpy(zci, vrefx, i*sizeof(int)); } last_d = 0; cl = 0; while(i > 0) { if ((a = zci[--i]) >= nv0 && a < max_var1) a = varp[a-nv0]->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 * #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(p, n) ps_func *p; int n; #else co_adjust(ps_func *p, int n) #endif { 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 #ifdef KR_headers ifadjust(asl, e) ASLTYPE *asl; expr_if *e; #else ifadjust(ASLTYPE *asl, expr_if *e) #endif { 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 #ifdef KR_headers vargadjust(asl, e) ASLTYPE *asl; expr_va *e; #else vargadjust(ASLTYPE *asl, expr_va *e) #endif { 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 #ifdef KR_headers funneladj1(asl, f) ASLTYPE *asl; funnel *f; #else funneladj1(ASLTYPE *asl, funnel *f) #endif { 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 #ifdef KR_headers funneladjust(S) Static *S; #else funneladjust(Static *S) #endif { 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 = (Char*)&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 #ifdef KR_headers goff_comp(asl) ASLTYPE *asl; #else goff_comp(ASLTYPE *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 { ASLTYPE *asl = S->asl; int *ka, *kae; ka = asl->i.A_colstarts_; kae = ka + asl->i.n_var0; while(ka <= kae) ++*ka++; } static void #ifdef KR_headers zerograd_chk(S) Static *S; #else zerograd_chk(Static *S) #endif { int j, n, nv, *z, **zg; ograd *og, **ogp, **ogpe; ASLTYPE *asl = S->asl; 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 #ifdef KR_headers cg_zzap(asl) ASLTYPE *asl; #else cg_zzap(ASLTYPE *asl) #endif { 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 #ifdef KR_headers adjust_compl_rhs(asl) ASL_fg *asl; #else adjust_compl_rhs(ASLTYPE *asl) #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) && Intcast 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(S, flags) Static *S; int flags; #else adjust(Static *S, int flags) #endif { 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 (n_obj) zerograd_chk(S); if (asl->i.n_con0 && !allJ) cg_zzap(asl); 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); } 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 *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 * #ifdef KR_headers bholread(R) EdRead *R; #else bholread(EdRead *R) #endif { 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 #ifdef KR_headers qwalk(S, e) Static *S; expr *e; #else qwalk(Static *S, expr *e) #endif { 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 < nv0) return 1; if (k >= max_var) { k_adj: k = ((expr_vx*)e)->a1; return k < 0 ? 1 : S->asl->I.v_class[k-nv0]; } return S->asl->I.v_class[k - nv0]; } return 3; } static int #ifdef KR_headers co_walkloop(S, f, n, c, o) Static *S; ps_func *f; int n; char *c; ograd **o; #else co_walkloop(Static *S, ps_func *f, int n, char *c, ograd **o) #endif { 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 #ifdef KR_headers do_ewalk(S) Static *S; #else do_ewalk(Static *S) #endif { 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 - nv0; noa = lasta00 + nc; nv1 = lasta00++; amax = lasta = lasta0 = max_var + nndv + 1; if (Ncom && asl->I.o_class) { i = max_var1 - nv0; 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 = r_ops[f_OPVARVAL]; for(ee = e + nv0 + 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 < nv0); 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 #ifdef KR_headers pfg_read_ASL(a, nl, flags) ASL *a; FILE *nl; int flags; #else pfg_read_ASL(ASL *a, FILE *nl, int flags) #endif { 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, maxfwd1, nc, nco, nlin, nlogc, no; int nv, nvc, nvo, nz, readall; ograd *og, **ogp; real t; unsigned x; ASL_CHECK(a, asltype, who); 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; if ((readall = flags & ASL_keep_all_suffixes) && a->i.nsuffixes) readall |= 1; 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]; 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; xscanf = bscanf; } else { holread = aholread; xscanf = ascanf; } asl->P.ncom = Ncom = comb + comc + como + comc1 + como1; nc = n_con; no = n_obj; nvc = c_vars; nvo = o_vars; if (no < 0 || (nco = nc + no + nlogc) <= 0) scream(R, ASL_readerr_corrupt, "pshvread: nc = %d, no = %d, nlogc = %d\n", nc, no, nlogc); if (pi0) { memset(pi0, 0, nc*sizeof(real)); if (havepi0) memset(havepi0, 0, nc); } asl->P.nv0_ = lasta00 = nv0 = nvc > nvo ? nvc : nvo; max_var1 = max_var = nv = nv0 + Ncom; combc = comb + comc; ncom0 = ncom_togo = combc + como; nzclim = max_var >> 3; ncom1 = comc1 + como1; nv0b = nv0 + 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.n_var_ - SS.nvar0)) SS.nvar0 += ncom0 + ncom1; if (flags & ASL_find_co_class) x += nco; 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 + nlogc; Ograd = (ograd **)(obj_de + no); if (A_vals) { if (!A_rownos) A_rownos = (int *)M1alloc(a->i.nzc_*sizeof(int)); } else if (nc) Cgrad = (cgrad **)M1zapalloc(nc*sizeof(cgrad *)); var_ex = e + nv0; 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; 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(nv0*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.Congrd = conpgrd_ASL; a->p.Objval = objpval_ASL; a->p.Objgrd = objpgrd_ASL; a->p.Conval = conpval_ASL; a->p.Jacval = jacpval_ASL; a->p.Lconval = lconpval_ASL; a->p.Hvcomp = hvpcomp_ASL; a->p.Hvinit = hvpinit_ASL; a->p.Hesset = hes_setup; a->p.Xknown = xp2known_ASL; a->p.Duthes = duthes_ASL; a->p.Fulhes = fullhes_ASL; a->p.Sphes = sphes_ASL; a->p.Sphset = sphes_setup_ASL; #else a->p.Xknown = xp1known_ASL; #endif 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); 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 < nv0 || 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 (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--) { if (ka) { 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, 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); } } } #undef asl #ifdef PSHVREAD static Char * #ifdef KR_headers new_mblkzap(asl, k) ASL_pfgh *asl; int k; #else new_mblkzap(ASL_pfgh *asl, int k) #endif { Char *rv = new_mblk_ASL((ASL*)asl,k); memset(rv, 0, sizeof(Char*)<fwd) { 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_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); } } return n; } static cexp * #ifdef KR_headers hesfunnel(S, ip, nrefs, ogp, lapp, lapep) Static *S; int *ip, nrefs; ograd **ogp; linarg ***lapp, ***lapep; #else hesfunnel(Static *S, int *ip, int nrefs, ograd **ogp, linarg ***lapp, linarg ***lapep) #endif { cexp *c; expr *e; int hw, i, k, m, n; linarg *la, **lap, **lape; ograd *og; dv_info *dvi; list *L; range *r; ASLTYPE *asl = S->asl; i = *ip; c = cexps + i; n = 0; if (i >= Ncom) { r = asl->P.Split_ce[i-Ncom].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; for(L = c->cref; L; L = L->next) n++; 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 * #ifdef KR_headers pscheck_ASL(a, who1) ASL *a; char *who1; #else pscheck_ASL(ASL *a, char *who1) #endif { ASL_CHECK(a, ASL_read_pfgh, who1); return (ASL_pfgh*)a; } static void #ifdef KR_headers hes_setup(a, flags, obj, nobj, con, ncon) ASL *a; int flags; int obj; int nobj; int con; int ncon; #else hes_setup(ASL *a, int flags, int obj, int nobj, int con, int ncon) #endif { range *r, *r0; int *c, *ce, *cei, i, k, n, n0, nmax, *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; nv0 = asl->P.nv0_; max_var = nv0 + 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 + nv0; 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; 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 - nv0]; 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 = nv0*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 + nv0); asl->P.otodo = (Hesoprod**)(asl->P.utodo + nv0); } 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