diff --git a/CMakeLists.txt b/CMakeLists.txt index 63ecb1c..c6b62a0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,7 +14,7 @@ # BUILD_MT_LIBS to build the multithreaded libraries # BUILD_F2C to build the fortran to c converter (https://www.netlib.org/f2c/f2c.pdf) # BUILD_EXAMPLES to build the examples - +cmake_minimum_required(VERSION 3.5) if (${CMAKE_VERSION} VERSION_GREATER "3.13.0") cmake_policy(SET CMP0077 NEW) endif() @@ -27,7 +27,7 @@ if (NOT ${CMAKE_VERSION} VERSION_LESS "3.9.0") endif() project(ASL) -cmake_minimum_required(VERSION 3.0) + set_property(GLOBAL PROPERTY USE_FOLDERS ON) # Set the path to CMake modules. set(AMPL_CMAKE_MODULE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/support/cmake) @@ -376,7 +376,7 @@ createSingleASL(asl ${ASL_SOURCE_DIR} ASL_SOURCES) # Create ASL 2 set(ADDITIONALDEFS "") if(WIN32) - set(ADDITIONALDEFS NO_MBLK_LOCK NO_PTHREADS) + set(ADDITIONALDEFS NO_MBLK_LOCK ) endif() createSingleASL(asl2 ${ASL2_SOURCE_DIR} ASL2_SOURCES DEFINITIONS ${ADDITIONALDEFS}) diff --git a/src/solvers/asldate.c b/src/solvers/asldate.c index 6a530b8..63f5a1c 100644 --- a/src/solvers/asldate.c +++ b/src/solvers/asldate.c @@ -1 +1 @@ -long ASLdate_ASL = 20241111; +long ASLdate_ASL = 20241202; diff --git a/src/solvers/changes b/src/solvers/changes index 6bea407..885c2c5 100644 --- a/src/solvers/changes +++ b/src/solvers/changes @@ -3433,3 +3433,31 @@ Current sph_opts bits: on a large example. Note that parallel Hessian evaluations are not available in solvers.tgz. + +20241111 + fpinit.c in solvers.tgz and solvers2.tgz: insert + #define _GNU_SOURCE +before + include "fenv.h" +to make fedisableexcept() visible on some systems (as someone +requested). It is "fenv.h" rather than because it was +sometimes necessary to supply a custom fenv.h, as indicated in +the comment + +/* Some Intel Linux systems (e.g., S.u.S.E. 5.2) come with a suitable fenv.h, */ +/* but some (e.g., S.U.S.E. 6.1) do not. This is for the latter systems. */ + +20241115 + pfghread.c in solvers2.tgz: fix an allocation bug introduced 20240618. + xp2known.c in solvers2.tgz: fix a glitch that appeared when not +compiled with -DALLOW_OPENMP . + +20241121 + solvers.c in solvers2.tgz: bypass a bug seen in the preprocessor of +at least one Microsoft compiler. On most systems, the change should +be invisible, except that under a debugger, some line numbers in +sphes.c will have changed. + +20241122 + asl.h in solvers2.tgz: only require pthreads.h when compiled with +-DMULTIPLE_THREADS (and not with -DALLOW_OPENMP). diff --git a/src/solvers/fgh_read.c b/src/solvers/fgh_read.c index d7d82af..8c83216 100644 --- a/src/solvers/fgh_read.c +++ b/src/solvers/fgh_read.c @@ -98,7 +98,7 @@ sorry_CLP(EdRead *R, const char *what) #define memadj(x) (((x) + (sizeof(long)-1)) & ~(sizeof(long)-1)) #endif - extern char* f_OPHOL; + extern char* f_OPHOL ANSI((expr* A_ASL)); extern efunc f_OPPLTERM, f_OPVARVAL, f_OPFUNCALL; extern sfunc f_OPIFSYM; diff --git a/src/solvers/makefile.u b/src/solvers/makefile.u index ef6c3fa..12a1948 100644 --- a/src/solvers/makefile.u +++ b/src/solvers/makefile.u @@ -305,7 +305,6 @@ xs0 = \ amplsolv.sy \ arith.ibm \ arith.h0 \ - arith.h1 \ arithchk.c \ asl.h \ asl_pfg.h \ diff --git a/src/solvers/xsum0.out b/src/solvers/xsum0.out index a271d54..0b527ae 100644 --- a/src/solvers/xsum0.out +++ b/src/solvers/xsum0.out @@ -4,12 +4,11 @@ amplsolv.lbc e92fe2c2 1075 amplsolv.sy 5451119 1343 arith.ibm f398be2b 1391 arith.h0 e206209c 2309 -arith.h1 f053bd25 177 arithchk.c 63b0185 6532 asl.h 1fe3527a 42925 asl_pfg.h 6483336 1268 asl_pfgh.h 49b5a53 1288 -asldate.c 8721c18 29 +asldate.c e6a239aa 29 atof.c 1fabc7d3 1747 auxinfo.c 1a3c8690 1488 avltree.c 10aeaa58 11346 @@ -41,7 +40,7 @@ errchk.h f88aec0 1695 f_read.c 12b32457 1227 fg_read.c f77acdfc 36818 fg_write.c 198882e0 17423 -fgh_read.c ed246ea7 34450 +fgh_read.c fb1c0aad 34470 float.h0 11d95cda 1822 fpecatch.c e1d8a914 1574 fpinit.c 66e5c26 5685 @@ -75,7 +74,7 @@ mach.c f480fbf5 2662 mainexit.c e4761b9b 1713 makefile.lc e5454c81 6036 makefile.sy 1e488eda 5581 -makefile.u 131b6278 11191 +makefile.u 1df62d7f 11179 makefile.vc 1bb69704 7500 makefile.wat c465856 5860 mip_pri.c f1dd1d47 7887 diff --git a/src/solvers2/asl.h b/src/solvers2/asl.h index f58a973..70dc209 100644 --- a/src/solvers2/asl.h +++ b/src/solvers2/asl.h @@ -378,6 +378,7 @@ Exitcall { void *v; }; +#define NO_PTHREADS #ifdef ALLOW_OPENMP #undef MULTIPLE_THREADS #define MULTIPLE_THREADS @@ -387,12 +388,13 @@ Exitcall { #define ACQUIRE_MBLK_LOCK(I,L) if ((I)->rd_M1z_bytes) omp_set_lock(&((I)->L)) #define FREE_MBLK_LOCK(I,L) if ((I)->rd_M1z_bytes) omp_unset_lock(&((I)->L)) #else -#ifdef NO_MBLK_LOCK +#if !defined(MULTIPLE_THREADS) || defined(NO_MBLK_LOCK) #undef MBLK_LOCK #undef MULTIPLE_THREADS #define ACQUIRE_MBLK_LOCK(I,L) /*nothing*/ #define FREE_MBLK_LOCK(I,L) /*nothing*/ #else +#undef NO_PTHREADS #include #define MBLK_LOCK pthread_mutex_t #ifndef MULTIPLE_THREADS diff --git a/src/solvers2/asldate.c b/src/solvers2/asldate.c index 6a530b8..5c7de60 100644 --- a/src/solvers2/asldate.c +++ b/src/solvers2/asldate.c @@ -1 +1 @@ -long ASLdate_ASL = 20241111; +long ASLdate_ASL = 20241122; diff --git a/src/solvers2/pfg_read.c b/src/solvers2/pfg_read.c deleted file mode 100644 index 4624c86..0000000 --- a/src/solvers2/pfg_read.c +++ /dev/null @@ -1,5007 +0,0 @@ -/******************************************************************* -Copyright (C) 2017, 2020 AMPL Optimization, Inc.; written by David M. Gay. - -Permission to use, copy, modify, and distribute this software and its -documentation for any purpose and without fee is hereby granted, -provided that the above copyright notice appear in all copies and that -both that the copyright notice and this permission notice and warranty -disclaimer appear in supporting documentation. - -The author and AMPL Optimization, Inc. disclaim all warranties with -regard to this software, including all implied warranties of -merchantability and fitness. In no event shall the author be liable -for any special, indirect or consequential damages or any damages -whatsoever resulting from loss of use, data or profits, whether in an -action of contract, negligence or other tortious action, arising out -of or in connection with the use or performance of this software. -*******************************************************************/ - -#ifdef DEBUG -#include "assert.h" -#else -#define assert(x) /*nothing*/ -#endif - -#ifdef PSHVREAD -#include "jacpdim.h" -#define efunc efunc2 -#include "opnos.hd" -#define PSHV(x) x -#define asltype ASL_read_pfgh -#define who "pfgh_read" -#define ASLTYPE ASL_pfgh -#else -#define PSHV(x) /*nothing*/ -#include "asl_pfg.h" -#define asltype ASL_read_pfg -#define who "pfg_read" -#define ASLTYPE ASL_pfg -#endif - -#ifdef __cplusplus -extern "C" { -#endif - -extern real conpival_nomap_ASL(ASL*, int, real*, fint*); -extern void conpgrd_nomap_ASL(ASL*, int, real *, real*, fint*); - -#define GAP_MAX 10 - -#ifdef PSHVREAD - static void hes_setup ANSI((ASL*, int, int, int, int, int)); -/* #define Static to avoid confusion in pi. */ -#define Static Static_pshv -#else -#define Static Static_psed -#endif -typedef struct Static Static; - static ograd *cotermwalk ANSI((Static*, expr **ep, ps_func *f, int wantg, int omitdv)); - static efunc *OPNUM, *OPVARVAL; - -#define NFHASH 23 - -#undef f_OPNUM -#include "r_opn0.hd" - -#undef nzc - -#if defined(IEEE_MC68k) || defined(IBM) - enum { W0 = 0, W1 = 1 }; -#else - enum { W0 = 1, W1 = 0 }; -#endif - -/*******/ - - typedef struct -expr_vx { - expr_v v; - linarg *la; - int a0; - int a1; - } expr_vx; - - typedef struct -Elemtemp { - unsigned int esize; - int nmax; - int k; - void **mp; - } Elemtemp; - - typedef struct -PSfind { - ps_func *f; - Elemtemp *b, *g; - } PSfind; - - typedef struct EU EU; - struct EU { efunc_n *op; union {real v; EU *p; } u; }; - - struct -Static { - ASLTYPE *asl; - ASL *a; - Elemtemp *_last_psb_elem; - derp *_last_d; - expr *(*_holread) ANSI((EdRead*)); - expr *_expr_free, **_slscratch; /* used in awalk(sumlist) */ - expr *_last_e; - expr_if *_iflist, *_if2list, *_if2list_end; - EU *_expr_n_free; - expr_v **_larvlist; - expr_v **_varp; - expr_va *_varglist, *_varg2list, *_varg2list_end; - int *_imap, *_vrefnext, *_vrefx; - int *_zc, *zc1, *_zci, *zci1, *_zl; - size_t _lthashmask, _nrange, _rangehashmask; - int klthash, krangehash; - int _allJ; - int _amax1; - int _cexp_k; - int _cexp_n; - int _conno; - int _groupno; - int _imap_len; - int _k_Elemtemp; - int _k_seen; - int _kimap; - int _kzc; - int _lasta; - int _lasta0; - int _lasta00; - int _max_var; - int _ncom; - int _ncom_togo; - int _nderp; - int _noa; - int _nocopy; - int _nndv; - int _nsce; - int _nv0b; - int _nv0c; - int _nv0r; - int _nv0x; - int _nv1; - int _nvref; - int _nzc; - int _nzclim; - int _slmax; - int _slscratchlev; - int _termno; /* current term number */ - int _wantCgroups; - int _wantOgroups; - int _zc_lim; - int nvar0, nvinc; - int size_exprn; - la_ref *_laref_free; - linarg **_lthash; /* hash table */ - linarg *_ltfree; /* free list */ - linarg *_tlist; /* linargs in this term */ - ograd *_freeog; /* free list */ - range **_rangehash; - real *_rnz; - real _lt_scale; - relo *_relolist, *_relo2list; - EdRead *R; - }; - -#define allJ S->_allJ -#define amax1 S->_amax1 -#define cexp_k S->_cexp_k -#define cexp_n S->_cexp_n -#define Conno S->_conno -#define expr_free S->_expr_free -#define expr_n_free S->_expr_n_free -#define freeog S->_freeog -#define Groupno S->_groupno -#define holread S->_holread -#define if2list S->_if2list -#define if2list_end S->_if2list_end -#define iflist S->_iflist -#define imap S->_imap -#define imap_len S->_imap_len -#define k_Elemtemp S->_k_Elemtemp -#define k_seen S->_k_seen -#define kimap S->_kimap -#define kzc S->_kzc -#define laref_free S->_laref_free -#define larvlist S->_larvlist -#define last_d S->_last_d -#define last_e S->_last_e -#define last_psb_elem S->_last_psb_elem -#define lasta S->_lasta -#define lasta0 S->_lasta0 -#define lasta00 S->_lasta00 -#define lt_scale S->_lt_scale -#define ltfree S->_ltfree -#define lthash S->_lthash -#define lthashmask S->_lthashmask -#define max_var S->_max_var -#define max_var1 asl->P.max_var1_ -#define Ncom S->_ncom -#define ncom_togo S->_ncom_togo -#define nderp S->_nderp -#define noa S->_noa -#define nocopy S->_nocopy -#define nndv S->_nndv -#define nrange S->_nrange -#define nsce S->_nsce -#define nv0b S->_nv0b -#define nv0c S->_nv0c -#define nv0r S->_nv0r -#define nv0x S->_nv0x -#define nv1 S->_nv1 -#define nvref S->_nvref -#define nzc S->_nzc -#define nzclim S->_nzclim -#define rangehash S->_rangehash -#define rangehashmask S->_rangehashmask -#define relo2list S->_relo2list -#define relolist S->_relolist -#define rnz S->_rnz -#define slmax S->_slmax -#define slscratch S->_slscratch -#define slscratchlev S->_slscratchlev -#define Termno S->_termno -#define tlist S->_tlist -#define varg2list S->_varg2list -#define varg2list_end S->_varg2list_end -#define varglist S->_varglist -#define varp S->_varp -#define vrefnext S->_vrefnext -#define vrefx S->_vrefx -#define wantCgroups S->_wantCgroups -#define wantOgroups S->_wantOgroups -#define zc S->_zc -#define zc_lim S->_zc_lim -#define zci S->_zci -#define zl S->_zl - -static list *crefs ANSI((Static*)); - - static Static * -S_init(Static *S, ASLTYPE *asl) -{ - memset(S, 0, sizeof(Static)); - S->asl = asl; - S->a = (ASL*)asl; - return S; - } - - static void * -new_mblkzap(ASLTYPE *asl, int k) -{ - void *rv = new_mblk_ASL((ASL*)asl,k); - memset(rv, 0, sizeof(void*)<i.memLast = asl->i.memNext = 0; - memset(asl->mblk_free, 0, MBLK_KMAX*sizeof(char*)); - memset(&asl->P.merge + 1, 0, sizeof(ps_info) - sizeof(long)); -#ifdef PSHVREAD - asl->P.hop_free = 0; - asl->P.hes_setup_called = 0; -#endif - } - -#ifdef Double_Align -#define memadj(x) x -#else -#define memadj(x) (((x) + (sizeof(long)-1)) & ~(sizeof(long)-1)) -#endif -#ifdef DEBUG - static expr *opzork; - static derp *dzork; - static ograd *ogzork; - static expr *ezork; - static EU *enzork; - static int dzork1, dzork2, exprzork, exprzork1, izork = -1; -#endif - - static Elemtemp * -new_Elemtemp(Static *S, unsigned int esize, void **mp) -{ - Elemtemp *e; - int k; - ASL *asl = S->a; - - e = (Elemtemp *)new_mblk(k_Elemtemp); - e->esize = esize; - e->mp = mp; - e->k = k = htcl(8*esize); - *mp = new_mblk(k); - e->nmax = (sizeof(void*) << k) / esize; - return e; - } - - static void -del_Elemtemp(Static *S, Elemtemp *e) -{ - ASL *asl = S->a; - del_mblk(e->k, *e->mp); - del_mblk(k_Elemtemp, e); - } - - static void -upgrade_Elemtemp(Static *S, Elemtemp *e) -{ - void *m, *m0; - int k; - ASL *asl = S->a; - - k = e->k++; - memcpy(m = new_mblk(e->k), m0 = *e->mp, e->esize * e->nmax); - del_mblk(k++, m0); - *e->mp = m; - e->nmax = (sizeof(void*) << k) / e->esize; - } - - static void -free_laref(Static *S, la_ref **L) -{ - la_ref *L1, *L2; - - if ((L1 = *L)) { - while((L2 = L1->next)) - L1 = L2; - L1->next = laref_free; - laref_free = *L; - *L = 0; - } - } - - static la_ref * -new_laref(Static *S, la_ref *nxt) -{ - la_ref *rv; - - if ((rv = laref_free)) - laref_free = rv->next; - else - rv = (la_ref *)mem_ASL(S->a, sizeof(la_ref)); - rv->next = nxt; - return rv; - } - - static list * -new_list(ASL *asl, list *nxt) -{ - list *rv = (list *)mem(sizeof(list)); - rv->next = nxt; - return rv; - } - - static void -fscream(EdRead *R, const char *name, int nargs, const char *kind) -{ - scream(R, ASL_readerr_argerr, - "line %ld: attempt to call %s with %d %sargs\n", - R->Line, name, nargs, kind); - } - - static void -new_derp(Static *S, int a, int b, real *c) -{ - derp *d; - nderp++; - d = (derp *)mem_ASL(S->a, sizeof(derp)); - d->next = last_d; - last_d = d; - d->a.i = a; - d->b.i = b; - d->c.rp = c; -#ifdef DEBUG - if (d == dzork) - printf(""); - if (++dzork1 == dzork2) - printf(""); -#endif - } - - static derp * -new_relo(Static *S, expr *e, derp *Dnext, int *ap) -{ - relo *r; - derp *d; - - r = (relo *)mem_ASL(S->a, sizeof(relo)); - r->next = relolist; - r->next2 = relo2list; - relo2list = relolist = r; - if (last_d != Dnext) { - *ap = e->a; - for(d = last_d; d->next != Dnext; d = d->next); - d->next = 0; - } - else { - last_d = 0; - new_derp(S, e->a, *ap = lasta++, &edagread_one); - } - if (!last_d) - return 0; - r->D = r->Dcond = last_d; - r->Dnext = Dnext; - return r->D; - } - - static relo * -new_relo1(Static *S, derp *Dnext) -{ - relo *r; - - r = (relo *)mem_ASL(S->a, sizeof(relo)); - r->next = relolist; - relolist = r; - r->D = 0; - r->Dnext = Dnext; - return r; - } - - static void -efree(Static *S, expr *e) -{ - EU *en; - expr **args, **argse, *e1; - top: - switch(optypeb[Intcast e->op]) { - - case 2: /* binary */ - efree(S, e->R.e); - /* no break */ - - case 1: /* unary */ - e1 = e->L.e; - e->L.e = expr_free; - expr_free = e; - e = e1; - goto top; - - case 6: /* sumlist */ - args = e->L.ep; - argse = e->R.ep; - while(args < argse) - efree(S, *args++); - e->L.e = expr_free; - expr_free = e; - break; - - case 9: - en = (EU*)e; - en->u.p = expr_n_free; - expr_n_free = en; - } - } - - static expr * -new_expr(Static *S, int o, expr *L, expr *R) -{ - expr *rv; - -#ifdef DEBUG - if (++exprzork1 == exprzork) - printf("exprzork1 = %d\n", exprzork1);; -#endif - if ((rv = expr_free)) - expr_free = rv->L.e; - else - rv = (expr *)mem_ASL(S->a, sizeof(expr)); - PSHV(rv->dL2 = 0); - if (o == f_OPPOW) { - if (Intcast R->op == f_OPNUM) { - if (((expr_n *)R)->v == 2.) { - o = f_OP2POW; - R = 0; - PSHV(rv->dL2 = 2); - } - else { - o = f_OP1POW; - rv->dR = ((expr_n *)R)->v; - /* rv->R.e may be used for a back pointer */ - } - } - else if (Intcast L->op == f_OPNUM) - o = f_OPCPOW; - } - rv->op = (efunc *)(size_t)o; - rv->L.e = L; - rv->R.e = R; - rv->a = 0; -#ifdef DEBUG - if (rv == ezork) - printf(""); -#endif - return rv; - } - - static void -dexpr(Static *S, expr *e, expr *L, expr *R) -{ - int L1, R1; -#ifdef PSHVREAD - int i, opcode; -#endif - - e->a = noa; - L1 = L && L->op != OPNUM && L->a != noa; - R1 = R && R->op != OPNUM && R->a != noa; - if (L1 | R1) { - if (L1) - new_derp(S, L->a, lasta, &e->dL); - if (R1) - new_derp(S, R->a, lasta, &e->dR); - e->a = lasta++; -#ifdef PSHVREAD - e->bak = last_e; - last_e = e; - opcode = Intcast e->op; - if (R) - e->dLR = e->dR2 = 0; - if (L1) - if (R1) - switch(opcode) { - case PLUS: i = Hv_plusLR; break; - case MINUS: i = Hv_minusLR; break; - case MULT: i = Hv_timesLR; break; - case DIV: i = Hv_divLR; break; - default: i = Hv_binaryLR; - } - else switch(opcode) { - case PLUS: - case MINUS: i = Hv_plusL; break; - case DIV: - case MULT: i = Hv_timesL; break; - case UMINUS: i = Hv_negate; break; - default: i = Hv_unary; - } - else - switch(opcode) { - case PLUS: i = Hv_plusR; break; - case MINUS: i = Hv_minusR; break; - case MULT: i = Hv_timesR; break; - default: i = Hv_binaryR; - } - e->dO.i = i; -#endif - } - } - - static real dvalue[] = { -#include "dvalue.hd" - }; - - static expr_n * -new_expr_n(Static *S, real t) -{ - EU *en; - - if ((en = expr_n_free)) - expr_n_free = en->u.p; - else - en = (EU *)mem_ASL(S->a, S->size_exprn); -#ifdef DEBUG - if (en == enzork) - printf(""); -#endif - en->u.v = t; - en->op = (efunc_n *)f_OPNUM; - return (expr_n*)en; - } - - static expr * -eread(EdRead *R) -{ - ASLTYPE *asl; - Static *S; - arglist *al; - argpair *ap, *sap; - char **sa; - de *d; - expr *L, *arg, **args, **argse, *rv; - expr_f *rvf; - expr_if *rvif; - expr_va *rva; - fint L1; - func_info *fi; - int *at, i, j, k, ks, numargs, op, symargs; - int (*Xscanf)(EdRead*, const char*, ...); - plterm *p; - real *b, *ra; - real r; - - S = (Static *)R->S; - asl = S->asl; - Xscanf = xscanf; - switch(edag_peek(R)) { - case 'f': - if (Xscanf(R, "%d %d", &i, &j) != 2) - badline(R); - fi = funcs[i]; - if (fi->nargs >= 0) { - if (fi->nargs != j) { - bad_nargs: - fscream(R, fi->name, j, ""); - } - } - else if (-(1+fi->nargs) > j) - goto bad_nargs; - rvf = (expr_f *)mem(sizeof(expr_f) - + (j-1)*sizeof(expr *)); - rvf->op = (efunc *)f_OPFUNCALL; - rvf->fi = fi; - PSHV(rvf->dO.i = Hv_func); - args = rvf->args; - argse = args + j; - k = ks = symargs = numargs = 0; - while(args < argse) { - arg = *args++ = eread(R); - if ((op = Intcast arg->op) == f_OPHOL) - symargs++; - else if (op == f_OPIFSYM) - ks++; - else { - numargs++; - if (op != f_OPNUM) - k++; - } - } - symargs += ks; - if (symargs && !(fi->ftype & 1)) - fscream(R, fi->name, symargs, "symbolic "); - ra = (real *)mem(sizeof(arglist) - + (k+ks)*sizeof(argpair) - + numargs*sizeof(real) - + symargs*sizeof(char *) - + j*sizeof(int)); - al = rvf->al = (arglist *)(ra + numargs); - al->n = numargs + symargs; - al->nr = numargs; - al->ra = ra; - al->derivs = al->hes = 0; - al->funcinfo = fi->funcinfo; - al->AE = asl->i.ae; - al->sa = (Const char**)(sa = (char **)(al + 1)); - ap = rvf->ap = (argpair *)(sa + symargs); - sap = rvf->sap = ap + k; - at = al->at = (int *)(sap + ks); - symargs = numargs = 0; - for(args = rvf->args; args < argse; at++) { - arg = *args++; - if ((op = Intcast arg->op) == f_OPHOL) { - *at = --symargs; - *sa++ = ((expr_h *)arg)->sym; - } - else if (op == f_OPIFSYM) { - *at = --symargs; - sap->e = arg; - (sap++)->u.s = sa++; - } - else { - *at = numargs++; - if (op == f_OPNUM) - *ra = ((expr_n *)arg)->v; - else { - ap->e = arg; - (ap++)->u.v = ra; - } - ra++; - } - } - rvf->ape = ap; - rvf->sape = sap; - return (expr *)rvf; - - case 'h': - return holread(R); - case 's': - if (Xscanf(R, "%hd", &L1) != 1) - badline(R); - r = L1; - goto have_r; - - case 'l': - if (Xscanf(R, "%ld", &L1) != 1) - badline(R); - r = L1; - goto have_r; - - case 'n': - if (Xscanf(R, "%lf", &r) != 1) - badline(R); - have_r: - return (expr *)new_expr_n(S, r); - - case 'o': - break; - - case 'v': - if (Xscanf(R,"%d",&k) != 1 || k < 0) - badline(R); - if (k >= S->nvar0) - k += S->nvinc; - if (k > max_var) - badline(R); - return (expr *)(var_e + k); - - default: - badline(R); - } - - if (Xscanf(R, asl->i.opfmt, &k) != 1 || k < 0 || k >= N_OPS) - badline(R); - switch(optype[k]) { - - case 1: /* unary */ - rv = new_expr(S, k, eread(R), 0); - return rv; - - case 2: /* binary */ - L = eread(R); - rv = new_expr(S, k, L, eread(R)); - return rv; - - case 3: /* vararg (min, max) */ - i = -1; - Xscanf(R, "%d", &i); - if (i <= 0) - badline(R); - rva = (expr_va *)mem(sizeof(expr_va)); - PSHV(rva->val = 0;) - rva->op = (efunc *)(size_t)k; - rva->L.d = d = (de *)mem(i*sizeof(de) + sizeof(expr *)); - rva->next = varglist; - varglist = rva; - for(j = 0; i > 0; i--, d++) - d->e = eread(R); - d->e = 0; /* sentinnel expr * */ - return (expr *)rva; - - case 4: /* piece-wise linear */ - i = -1; - Xscanf(R, "%d", &i); - if (i <= 1) - badline(R); - plterms++; - j = 2*i - 1; - p = (plterm *)mem(sizeof(plterm) + (j-1)*sizeof(real)); - p->n = i; - b = p->bs; - do { - switch(edag_peek(R)) { - case 's': - if (Xscanf(R,"%hd",&L1) != 1) - badline(R); - r = L1; - break; - case 'l': - if (Xscanf(R,"%ld",&L1) != 1) - badline(R); - r = L1; - break; - case 'n': - if (Xscanf(R,"%lf",&r) == 1) - break; - default: - badline(R); - } - *b++ = r; - } - while(--j > 0); - if (b[-2] <= 0.) - p->z = 2*i - 2; - else { - b = p->bs + 1; - while(*b <= 0.) - b += 2; - p->z = (b - p->bs) - 1; - } - rv = (expr *)mem(sizeof(expr)); - rv->op = (efunc *)(size_t)k; - rv->L.p = p; - rv->R.e = eread(R); - return rv; - - case 5: /* if */ - rvif = (expr_if *)mem(sizeof(expr_if)); - PSHV(rvif->val = 0;) - PSHV(rvif->Te = rvif->Fe = 0;) /* for bogus .nl files */ - rvif->op = (efunc *)(size_t)k; - rvif->next = iflist; - iflist = rvif; - rvif->e = eread(R); - rvif->T = L = eread(R); - rvif->F = L = eread(R); - return (expr *)rvif; - - case 11: /* OPCOUNT */ - case 6: /* sumlist */ - i = 0; - Xscanf(R, "%d", &i); - if (i <= 2 && (optype[k] == 6 || i < 1)) - badline(R); - if (slmax < i) - slmax = i; - args = (expr**)new_mblk(htcl(i*sizeof(expr*))); - rv = new_expr(S, k, (expr*)args, (expr*)(args+i)); - do - *args++ = eread(R); - while(--i > 0); - return rv; - } - badline(R); - return 0; - } - - static void -ogfree(Static *S, ograd *og) -{ - ograd *og1; - - if ((og1 = og)) { -#ifdef DEBUG - if (ogzork) { - do if (og == ogzork) - printf(""); - while((og = og->next)); - og = og1; - } -#endif - while(og1->next) - og1 = og1->next; - og1->next = freeog; - freeog = og; - } - } - - static real -ogfree1(Static *S, ograd **ogp) -{ - ograd *og = *ogp; - *ogp = og->next; - og->next = freeog; - freeog = og; - return og->coef; - } - - static ograd * -new_ograd(Static *S, ograd *next, int varno, real coef) -{ - ograd *rv; - - if ((rv = freeog)) - freeog = rv->next; - else - rv = (ograd *)mem_ASL(S->a, sizeof(ograd)); - rv->next = next; - rv->varno = varno; - rv->coef = coef; -#ifdef DEBUG - if (rv == ogzork) - printf(""); -#endif - return rv; - } - - static linarg * -lahash(Static *S, linarg *la) -{ - ASLTYPE *asl; - Ulong x; - int k, nnz; - linarg *la1, *la2, **lap, **lap1, **q, **q0, **qe; - ograd *og, *og1; - size_t mask; - union {real r; Ulong L[2];} u; - - asl = S->asl; - x = nnz = la->nnz; - for(og = la->nz; og; og = og->next) { - u.r = og->coef; - x = (x << 1 | (x & 0x80000000) >> 31) ^ - (og->varno*101 + u.L[W0] + 257*u.L[W1]); - } - for(lap = <hash[x & lthashmask]; (la1 = *lap); lap = &la1->hnext) - if (la1->nnz == nnz) { - og = la->nz; - og1 = la1->nz; - for(;;) { - if (!og) { - if (!og1) - return la1; - break; - } - if (!og1) - break; - if (og->varno != og1->varno - || og->coef != og1->coef) - break; - og = og->next; - og1 = og1->next; - } - } - *lap = la; - la->hnext = 0; - if (++S->asl->P.nlttot > lthashmask) { - mask = lthashmask; - k = S->klthash; - q = q0 = lthash; - qe = q + mask; - lthashmask = mask = (mask << 1) | 1; - lap = lthash = (linarg**)new_mblkzap(asl, S->klthash = k + 1); - while(q <= qe) { - for(la2 = *q++; la2; la2 = la1) { - la1 = la2->hnext; - x = la2->nnz; - for(og = la2->nz; og; og = og->next) { - u.r = og->coef; - x = (x << 1 | (x & 0x80000000) >> 31) ^ - (og->varno*101 + u.L[W0] + 257*u.L[W1]); - } - lap1 = lap + (x & mask); - la2->hnext = *lap1; - *lap1 = la2; - } - } - del_mblk(k, q0); - } - return la; - } - - static range * -new_range(Static *S, range *r, range **rp) -{ - ASLTYPE *asl; - int k, len, uilen; - range *r1, *r2, *r3, **rp1, **rq, **rq0, **rqe; - size_t mask; - - asl = S->asl; - uilen = r->nv*sizeof(int); - len = sizeof(range) + uilen; - r1 = (range*)mem(len); - r1->nintv = 0; - r1->n = r->n; - r1->nv = r->nv; - r1->chksum = r->chksum; - r1->refs = 0; - r1->lasttermno = -1; - r1->hnext = r1->hunext = 0; - if (uilen) - memcpy(r1->ui = (int*)(r1+1), r->ui, uilen); - r1->lap = (linarg**)new_mblk(htcl(len = r->n*sizeof(linarg*))); - memcpy(r1->lap, r->lap, len); - r2 = r1->rlist.next = asl->P.rlist.next; - r1->rlist.prev = r2->rlist.prev; - r2->rlist.prev = asl->P.rlist.next = r1; - *rp = r1; - if (++nrange > rangehashmask) { - mask = rangehashmask; - k = S->krangehash; - rq = rq0 = rangehash; - rqe = rq + mask; - rangehashmask = mask = (mask << 1) | 1; - rp = rangehash = (range**)new_mblkzap(asl, S->krangehash = k + 1); - while(rq <= rqe) { - for(r2 = *rq++; r2; r2 = r3) { - r3 = r2->hunext; - rp1 = rp + (r2->chksum & mask); - r2->hunext = *rp1; - *rp1 = r2; - } - } - del_mblk(k, rq0); - } - return r1; - } - - static int -lacompar(const void *a, const void *b, void *v) -{ - ograd *oa, *ob; - int i; - real t; - - if (a == b) - return 0; - Not_Used(v); - oa = (*(linarg **)a)->nz; - ob = (*(linarg **)b)->nz; - for(;;) { - if (!oa) { - if (ob) - return -1; - break; - } - if (!ob) - return 1; - if ((i = oa->varno - ob->varno)) - return i; - if ((t = oa->coef - ob->coef)) - return t > 0. ? 1 : -1; - oa = oa->next; - ob = ob->next; - } - return 0; - } - - static int -ndiff(range *r, range *r1) -{ - /* Return 0 if r and r1 have the same linargs; else return 1. */ - - linarg **la, **la1, **la1e, **lae; - - la = r->lap; - lae = la + r->n; - la1 = r1->lap; - la1e = la1 + r1->n; - for(;;++la, ++la1) { - if (la >= lae) - return la1 < la1e; - if (la1 >= la1e) - return 1; - if (lacompar((char*)la, (char *)la1, NULL)) - return 1; - } - return 0; - } - - static range * -uhash(Static *S, range *r) -{ - ASLTYPE *asl; - int len, n, nv, *ui, *uie; - range *r1, **rp; - size_t L; - - asl = S->asl; - L = 0; - nv = r->nv; - ui = r->ui; - uie = ui + nv; - len = sizeof(int)*nv; - while (ui < uie) - L = 37*L + *ui++; - r->chksum = L; - ui = r->ui; - rp = rangehash + (L & rangehashmask); - n = r->n; - if (asl->P.merge) - while((r1 = *rp)) { - if (r1->nv == nv && r1->n == n && (len == 0 || !memcmp(ui, r1->ui, len))) { - if (n == 0 || !memcmp(r->lap, r1->lap, n*sizeof(linarg*))) - return *rp; - if (!ndiff(r, r1)) - return *rp; - } - rp = &r1->hunext; - } - return new_range(S, r, rp); - } - - static range * -rhash(Static *S, range *r, int addnew) -{ - ASLTYPE *asl; - int len, n; - linarg **lae, **lap; - ograd *og; - range *r1, **rp; - unsigned long L; - - asl = S->asl; - L = 0; - lap = r->lap; - n = r->n; - lae = lap + n; - len = n * sizeof(linarg*); - while(lap < lae) { - L *= 37; - for(og = (*lap++)->nz; og; og = og->next) - L = 101*L + og->varno; - } - r->chksum = L; - lap = r->lap; - rp = rangehash + (L & rangehashmask); - if (asl->P.merge) - while((r1 = *rp)) { - if (r1->n == n && !memcmp(lap, r1->lap, len)) - return r1; - rp = &r1->hnext; - } - if (addnew) - return new_range(S, r, rp); - return *rp; - } - - static int -compar(const void *a0, const void *b0, void *v) -{ - int a, b, c; - Static *S = (Static *)v; - - if ((a = *(int*)a0) >= max_var) { - a = zl[a-nv0x]; - if ((b = *(int*)b0) >= max_var) { - if ((c = a - zl[b-nv0x])) - return c; - return *(int*)a0 - b; - } - if (a == b) - return -1; - } - else if ((b = *(int*)b0) >= max_var) { - b = zl[b-nv0x]; - if (a == b) - return 1; - } - return a - b; - } - - static int -hscompar(const void *a0, const void *b0, void *v) -{ - int a, b, c; - Static *S = (Static *)v; - - if ((a = *(int*)a0) >= Ncom) { - a = zl[a]; - if ((b = *(int*)b0) >= Ncom) { - if ((c = a - zl[b])) - return c; - return *(int*)a0 - b; - } - a -= nv0x; - if (a == b) - return -1; - } - else if ((b = *(int*)b0) >= Ncom) { - b = zl[b] - nv0x; - if (a == b) - return 1; - } - return a - b; - } - - static void -zcsort(Static *S, int *c, int *ci, int i, int n, int p) -{ - int j; - - if (n < nzclim || p < 0) - qsortv(ci, n, sizeof(int), compar, S); - else for(j = 0; i < p; i++) - if (c[i]) - ci[j++] = i; - } - - static ograd* -compress(Static *S, ograd *og, real *cp, int *comvar) -{ - ograd *og1; - int i, ix, j, j1, k, nzc1, *zc1, *zci1; - real c, t; - - c = og->varno < 0 ? ogfree1(S,&og) : 0.; - ix = nzc1 = 0; - zc1 = S->zc1; - zci1 = S->zci1; - for(og1 = og; og1; og1 = og1->next) { - zc1[i = og1->varno] = 1; - zci1[nzc1++] = i; - rnz[i] = og1->coef; - if (ix < i) - ix = i; - } - if (ix < nv0x) { - *cp = c; - *comvar = 0; - for(og1 = og; og1; og1 = og1->next) - zc1[og1->varno] = 0; - return og; - } - *comvar = 1; - for(i = 0; i < nzc1; ) - if ((j = zci1[i]) < nv0x) - i++; - else { - if (!zc[j]++) - zci[nzc++] = j; - t = rnz[j]; - j1 = j - nv0x; - if ((og1 = (S->asl->P.dv + j1)->ll)) { - if (og1->varno < 0) { - c += t*og1->coef; - og1 = og1->next; - } - for(; og1; og1 = og1->next) - if (!zc1[k = og1->varno]++) { - zci1[nzc1++] = k; - rnz[k] = t*og1->coef; - } - else - rnz[k] += t*og1->coef; - } - zc1[j] = 0; - zci1[i] = zci1[--nzc1]; - } - *cp = c; - ogfree(S, og); - og = 0; - if (nzc1 > 0) { - zcsort(S, zc1, zci1, 0, nzc1, max_var); - while(nzc1 > 0) { - i = zci1[--nzc1]; - zc1[i] = 0; - if ((t = rnz[i])) { - og = new_ograd(S, og, i, t); - if (!zc[i]++) - zci[nzc++] = i; - } - } - } - return og; - } - - static linarg * -afree(Static *S, ograd *og, expr **ep) -{ - linarg *la, *la1, *rv; - real c, s, t; - ograd *og1, *ogx; - int comvar, nnz; - la_ref *refs; - ASLTYPE *asl = S->asl; - - rv = 0; - if (!og || !(og = compress(S,og,&c,&comvar))) - goto done; - - if ((la = ltfree)) - ltfree = la->hnext; - else { - la = (linarg *)mem(sizeof(linarg)); - la->refs = 0; - } - - s = og->coef; - if (s < 0) - s = -s; - la->nz = ogx = og1 = og; - nnz = 1; - while((og1 = og1->next)) { - nnz++; - t = og1->coef; - if (t < 0) - t = -t; - if (s < t) { - s = t; - ogx = og1; - } - } - - la->nnz = nnz; - if ((s = ogx->coef) != 1.) - for(og1 = og; og1; og1 = og1->next) - og1->coef /= s; - lt_scale = s; - la1 = lahash(S, la); - if (la1 == la) { - la->refs = 0; - la->v = 0; - la->termno = Termno; - la->tnext = tlist; - tlist = la; - la->lnext = asl->P.lalist; - asl->P.lalist = la; - la->hnext = 0; - } - else { - /* !!? Eventually create or adjust la1->v to save cycles. */ - /* This requires possible adjustment for differing */ - /* constants and scale factors. */ - if (la1->termno == Termno) - asl->P.ndupst++; /* research statistic */ - else { - free_laref(S, &la->refs); - la1->termno = Termno; - la1->tnext = tlist; - tlist = la1; - asl->P.ndupdt++; /* research statistic */ - } - ogfree(S, og); - la->hnext = ltfree; - ltfree = la; - la = la1; - } - if (ep && (nnz > 1 || comvar)) { - la->refs = refs = new_laref(S, la->refs); - refs->ep = ep; - refs->c = c; - refs->scale = s; - } - if (nnz > 1) - rv = la; - done: - return rv; - } - - static ograd * -af_sum(Static *S, ograd *Log, ograd *Rog) -{ - ograd *oL, *oR, *oR1, *og, **ogp; - - oL = Log; - oR = Rog; - ogp = &og; - for(;;) { - if (!oL) { - *ogp = oR; - break; - } - if (!oR) { - *ogp = oL; - break; - } - if (oL->varno > oR->varno) { - *ogp = oR; - ogp = &oR->next; - oR = *ogp; - } - else { - if (oL->varno == oR->varno) { - oL->coef += oR->coef; - oR1 = oR->next; - oR->next = 0; - ogfree(S, oR); - oR = oR1; - if (oL->coef == 0.) { - oR1 = oL->next; - oL->next = 0; - ogfree(S, oL); - oL = oR1; - continue; - } - } - *ogp = oL; - ogp = &oL->next; - oL = *ogp; - } - } - return og; - } - - static cgrad * -cf_sum(Static *S, cgrad *Lcg, ograd *Rog) -{ - cgrad *cg, **cgp, *oL; - ograd *oR, *oR1; - - oL = Lcg; - oR = Rog; - cgp = &cg; - cg = 0; - for(;;) { - if (!oL) { - assert(!oR); - break; - } - if (!oR) { - *cgp = oL; - break; - } - assert(oL->varno <= oR->varno); - if (oL->varno == oR->varno) { - oL->coef += oR->coef; - oR1 = oR->next; - oR->next = 0; - ogfree(S, oR); - oR = oR1; - } - *cgp = oL; - cgp = &oL->next; - oL = *cgp; - } - return cg; - } - - static void -sumlist_afree(Static *S, ograd *Laf, expr *e, expr **argso, expr **sls, expr **slscr) -{ - int nlin, nnl; - expr *e1, *e2, **ep, **ep1; - ASL *asl = S->a; - - nlin = sls - slscr; - ep = e->L.ep; - nnl = argso - ep; - switch(nlin) { - case 1: - e1 = slscr[0]; - break; - case 2: - e1 = new_expr(S, f_OPPLUS, slscr[0], slscr[1]); - e1->dL = e1->dR = 1.; - break; - default: - ep1 = (expr**)new_mblk(htcl(nlin*sizeof(expr*))); - e1 = new_expr(S, f_OPSUMLIST, (expr*)ep1, (expr*)(ep1+nlin)); - memcpy(ep1, slscr, nlin*sizeof(expr*)); - } - switch(nnl) { - case 1: - e2 = *ep; - break; - case 2: - e2 = new_expr(S, f_OPPLUS, ep[0], ep[1]); - e2->dL = e2->dR = 1.; - break; - default: - ep1 = (expr**)new_mblk(htcl(nnl*sizeof(expr*))); - e2 = new_expr(S, f_OPSUMLIST, (expr*)ep1, (expr*)(ep1+nnl)); - memcpy(ep1, ep, nnl*sizeof(expr*)); - } - del_mblk(htcl((e->R.ep - ep)*sizeof(expr*)), ep); - e->op = (efunc*)f_OPPLUS; - e->dL = e->dR = 1.; - e->L.e = e1; - e->R.e = e2; - afree(S, Laf, &e->L.e); - } - - static void -sdvcite(Static *S, int k) -{ - range *r; - split_ce *cs; - linarg *la, **lap, **lape; - - cs = S->asl->P.Split_ce + (k - max_var); - r = cs->r; - lap = r->lap; - lape = lap + r->n; - while(lap < lape) { - la = *lap++; - if (la->termno != Termno) { - free_laref(S, &la->refs); - la->termno = Termno; - la->tnext = tlist; - tlist = la; - } - } - } - - static ograd * -awalk(Static *S, expr *e) /* return 0 if e is not linear */ -{ - ASLTYPE *asl; - int k, k1, kscr; - expr *L, **args, **argse, **argso, **sls, **slscr; - expr_va *rva; - de *d; - expr_if *rvif; - expr_f *rvf; - ograd *Laf, *Raf, *rv, *taf; - linarg *la, **nl; - argpair *ap, *ape; - real t; - - switch(optypeb[k = Intcast e->op]) { - - case 1: /* unary */ - Laf = awalk(S, e->L.e); - if (k == f_OPUMINUS) { - if ((taf = Laf)) { - do taf->coef = -taf->coef; - while((taf = taf->next)); - return Laf; - } - } - if (Laf) - afree(S, Laf, &e->L.e); - break; - - case 2: /* binary */ - Laf = awalk(S, e->L.e); - Raf = awalk(S, e->R.e); - if (Laf && Raf) - switch(k) { - case f_OPMINUS: - taf = Raf; - do taf->coef = -taf->coef; - while((taf = taf->next)); - /* no break; */ - case f_OPPLUS: - return af_sum(S, Laf, Raf); - - case f_OPMULT: - if (Raf->varno < 0 && !Raf->next) { - swap: - taf = Laf; - Laf = Raf; - Raf = taf; - goto scale; - } - if (Laf->varno < 0 && !Laf->next) { - taf = Raf; - scale: t = Laf->coef; - if (t == 0.) { - ogfree(S, Raf); - return Laf; - } - do taf->coef *= t; - while((taf = taf->next)); - ogfree(S, Laf); - return Raf; - } - break; - - case f_OPDIV: - if (Raf->varno < 0 && !Raf->next) { - Raf->coef = 1. / Raf->coef; - goto swap; - } - } - afree(S, Laf, &e->L.e); - afree(S, Raf, &e->R.e); - break; - - case 3: /* vararg (min, max) */ - rva = (expr_va *)e; - for(d = rva->L.d; (L = d->e); d++) - afree(S, awalk(S,L), &d->e); - break; - - case 4: /* piece-wise linear */ - afree(S, awalk(S,e->R.e), &e->R.e); - break; - - case 5: /* if */ - rvif = (expr_if *)e; - afree(S, awalk(S,rvif->T), &rvif->T); - afree(S, awalk(S,rvif->F), &rvif->F); - break; - - case 6: /* sumlist */ - args = e->L.ep; - argse = e->R.ep; - k1 = argse - args; - while(!(Laf = awalk(S,*args++))) - if (args >= argse) - return 0; - kscr = -1; - asl = S->asl; - if (!slscratchlev++) - slscr = slscratch; - else { - kscr = htcl(k1*sizeof(expr*)); - slscr = (expr**)new_mblk(kscr); - } - sls = slscr; - argso = args - 1; - *sls++ = *argso; - while(args < argse) - if ((Raf = awalk(S,*args))) { - *sls++ = *args++; - Laf = af_sum(S, Laf, Raf); - } - else - *argso++ = *args++; - rv = 0; - if (argso == e->L.ep) { - rv = Laf; - goto delscratch; - } - sumlist_afree(S, Laf, e, argso, sls, slscr); - delscratch: - --slscratchlev; - if (kscr >= 0) - del_mblk(kscr, slscr); - return rv; - - case 7: /* function call */ - rvf = (expr_f *)e; - ap = rvf->ap; - for(ape = rvf->ape; ap < ape; ap++) - afree(S, awalk(S,ap->e), &ap->e); - - case 8: /* OPHOL (Hollerith) */ - break; - - case 9: /* OPNUM */ - return new_ograd(S, 0, -1, ((expr_n*)e)->v); - - case 10:/* OPVARVAL */ - asl = S->asl; - k = (expr_v *)e - var_e; - if ((k < 0 || k >= max_var) - && (k = ((expr_vx*)e)->a0) < 0) { - if ((la = ((expr_vx*)e)->la) - && la->termno != Termno) { - la->termno = Termno; - la->tnext = tlist; - tlist = la; - } - return 0; - } - if ((k1 = (k - nv0x)) < 0) { - if (!zc[k]++) - zci[nzc++] = k; - goto ogret; - } - if (k1 >= Ncom) { - if (!zc[k = ((expr_vx*)e)->a1]++) - zci[nzc++] = k; - sdvcite(S,k); - return 0; - } - if ((nl = (asl->P.dv + k1)->nl)) { - if (!zc[k]++) - zci[nzc++] = k; - while((la = *nl++)) - if (la->termno != Termno) { - la->termno = Termno; - la->tnext = tlist; - tlist = la; - } - return 0; - } - if (!zc[k]++) - zci[nzc++] = k; - ogret: - return new_ograd(S, 0, k, 1.); - - case 11: /* OPCOUNT */ - break; - - default: - scream(S->R, ASL_readerr_bug, - "awalk: unexpected optype[%d] = %d\n", - k, optype[k]); - } - return 0; - } - - static void -cexp_upgrade(Static *S, int t) -{ - int k, n1, n2; - cexp *ce; - int *z; - expr_v **vp; - split_ce *cs; - ASLTYPE *asl = S->asl; - - n2 = t - Ncom; - k = htcl(t*(sizeof(cexp) + sizeof(int) + sizeof(expr_v*)) - + n2*sizeof(split_ce)); - memset(ce = (cexp *)new_mblk(k), 0, n1 = sizeof(char*) << k); - n1 = (n1 + Ncom*sizeof(split_ce)) - / (sizeof(cexp) + sizeof(int) + sizeof(split_ce) - + sizeof(expr_v*)); - n2 = n1 - Ncom; - cs = (split_ce *)(ce + n1); - vp = (expr_v **)(cs + n2); - z = (int *)(vp + n1); - if (cexps) { - if (nsce) - memcpy(cs, asl->P.Split_ce, nsce*sizeof(split_ce)); - memcpy(ce, cexps, cexp_n*sizeof(cexp)); - memcpy(z, zl, cexp_n*sizeof(int)); - memcpy(vp, varp, cexp_n*sizeof(expr_v*)); - del_mblk(cexp_k, cexps); - } - nsce = n2; - asl->P.Split_ce = cs; - cexps = ce; - zl = z; - cexp_k = k; - cexp_n = n1; - varp = vp; - } - - static void -zc_upgrade(Static *S) -{ - int k, n, n0; - int *z; - ASLTYPE *asl = S->asl; - - k = htcl(sizeof(int)*(max_var1 + 1)) + 1; - z = (int*)new_mblk(k); - n = (sizeof(char*)/sizeof(int)) << (k-1); - memset(z + n, 0, n*sizeof(int)); - if (zci) { - n0 = (sizeof(char*)/sizeof(int)) << (kzc - 1); - memcpy(z, zci, n0*sizeof(int)); - memcpy(z+n, zci+n0, n0*sizeof(int)); - del_mblk(kzc, zci); - } - kzc = k; - zci = z; - zc = z + n + 1; - zc_lim = n; - } - - static void -la_replace(Static *S, linarg *la) -{ - la_ref *r; - expr_v *v; - expr_vx *vx; - expr *eout; - ASLTYPE *asl = S->asl; - - if (la->nnz > 1) { - if (!(v = la->v)) { - vx = (expr_vx *)mem(sizeof(expr_vx)); - vx->la = la; - vx->a0 = vx->a1 = -1; - la->v = v = (expr_v *)vx; - v->a = max_var + nndv++; - lasta00++; - v->op = (efunc *)f_OPVARVAL; - if (larvlist) { - *larvlist = v; - larvlist = (expr_v**)&v->v; - } - } - } - else - v = var_e + la->nz->varno; - for(r = la->refs; r; r = r->next) { - efree(S, *r->ep); - eout = (expr *)v; - if (r->scale != 1.) { - if (r->scale == -1.) { - eout = new_expr(S, f_OPUMINUS, eout, 0); - eout->dL = -1.; - } - else - eout = new_expr(S, f_OPMULT, eout, - (expr*)new_expr_n(S, r->scale)); - } - if (r->c) { - eout = new_expr(S, f_OPPLUS, eout, - (expr*)new_expr_n(S, r->c)); - eout->dL = 1.; - } - *r->ep = eout; - } - free_laref(S, &la->refs); - } - - static void -tlistgen(Static *S, ps_func *f) -{ - int *ci, *cie, i, t; - linarg *la, **lap, **lape, *tl; - ograd *og; - range *r; - psb_elem *b, *be; - - t = f->nb; - b = f->b; - be = b + t; - t = ++Termno; - tl = 0; - for(; b < be; b++) { - if ((ci = b->ce)) { - cie = ci + *ci; - do { - if (!zc[i = nv0x + *++ci]) - zci[nzc++] = i; - } - while(ci < cie); - } - r = b->U; - lap = r->lap; - lape = lap + r->n; - while(lap < lape) { - la = *lap++; - if (la->termno != t) { - la->termno = t; - la->tnext = tl; - tl = la; - for(og = la->nz; og; og = og->next) - if (!zc[og->varno]++) - zci[nzc++] = og->varno; - } - } - } - tlist = tl; - } - - static int -might_expand(Static *S, expr *e) -{ - top: - switch(Intcast e->op) { - case f_OPPLUS: - case f_OPMINUS: - case f_OPUMINUS: - case f_OPSUMLIST: - return 1; - case f_OPMULT: - if (Intcast e->R.e->op == f_OPNUM) { - e = e->L.e; - goto top; - } - if (Intcast e->L.e->op == f_OPNUM) { - e = e->R.e; - goto top; - } - break; - case f_OPVARVAL: - if (((expr_v*)e)->a >= nv0x) - return 1; - } - return 0; - } - - static void -ce_split(Static *S, int i, ps_func *f) -{ - int j, j1, je, k, n; - cexp *c, *ce; - expr **ep; - expr_v **vp, **vp0; - expr_vx *vx; - split_ce *cs; - psb_elem *b; - ASLTYPE *asl = S->asl; - - asl->P.ndvsplit++; - n = f->nb; - j1 = asl->P.ndvspout; - j = k = j1 + Ncom; - asl->P.ndvspout += n; - je = j + n; - if (je > cexp_n) - cexp_upgrade(S, je); - c = cexps + j; - b = f->b; - cs = asl->P.Split_ce + (j-Ncom); - for(ce = c + n; c < ce; c++, b++, cs++) { - c->e = b->D.e; - cs->r = b->U; - cs->ce = b->ce; - } - c = cexps + i; - vp0 = vp = varp + j; - j = max_var + nndv; - je = j + n; - nndv += n; - j1 += max_var; - while(j < je) { - vx = (expr_vx *)mem(sizeof(expr_vx)); - vx->la = 0; - *vp++ = (expr_v*)vx; - vx->v.a = vx->a0 = j++; - vx->a1 = j1++; - vx->v.op = (efunc *)f_OPVARVAL; - } - if (n == 2) - c->e = new_expr(S, f_OPPLUS, (expr*)vp0[0], (expr*)vp0[1]); - else { - ep = (expr**)new_mblk(htcl(n*sizeof(expr*))); - memcpy(ep, vp0, n*sizeof(expr*)); - c->e = new_expr(S, f_OPSUMLIST, (expr*)ep, (expr*)(ep+n)); - } - if ((max_var1 += n) >= zc_lim) - zc_upgrade(S); - - /* adjust zl for use in compar() */ - i += nv0x; - j = k + nv0x; - n += k; - do { - zl[k++] = i; - zci[nzc++] = j++; /* for crefs */ - } - while(k < n); - } - - static void -linpart_augment(Static *S, cexp *c, ograd *og, ps_func *f) -{ - linpart *L, *Le; - int n; - ograd *og1; - real t; - ASL *asl = S->a; - - if (og->varno < 0) { - if ((t = ogfree1(S, &og))) - f->b->D.e = new_expr(S, f_OPPLUS, f->b->D.e, - (expr*)new_expr_n(S, t)); - if (!og) - return; - } - if ((n = c->nlin)) { - L = c->L; - Le = L + n; - og1 = 0; - do { - --Le; - og1 = new_ograd(S, og1, Le->v.i, Le->fac); - } - while(Le > L); - del_mblk(htcl(htcl(n*sizeof(linpart))), c->L); - og = af_sum(S, og, og1); - } - og1 = og; - for(n = 0; og; og = og->next) - n++; - c->L = L = (linpart*)new_mblk(htcl((c->nlin = n)*sizeof(linpart))); - for(og = og1; og; og = og->next, L++) { - L->v.i = og->varno; - L->fac = og->coef; - } - ogfree(S, og1); - } - - static ograd * -linterms(Static *S, cexp *c, real scale) -{ - linpart *L = c->L; - linpart *Le = L + c->nlin; - ograd *og = 0; - while(Le > L) { - --Le; - og = new_ograd(S, og, Le->v.i, scale*Le->fac); - } - return og; - } - - static void -dvwalk(Static *S, int i) -{ - cexp *c; - dv_info *dvi; - expr *e; - int n, split; - linarg *tl, **tp; - ograd *og, *og0; - ps_func f; - ASLTYPE *asl = S->asl; - - Termno++; - tlist = 0; - c = cexps + i; - split = zl[i] & 2; - zl[i] = 0; - if (split && might_expand(S, c->e)) { - asl->P.ndvspcand++; - if ((og = cotermwalk(S, &c->e, &f, 0, 0)) && f.nb >= 1) - linpart_augment(S, c, og, &f); - if (f.nb > 1) { - zl[i] = f.nb; - ce_split(S, i, &f); - c = cexps + i; - og = 0; - } - else if (f.nb == 1) { - c->e = f.b->D.e; - og = 0; - } - tlistgen(S, &f); - del_Elemtemp(S, last_psb_elem); - } - else if ((e = c->e)) - og = awalk(S, e); - else - og = 0; - og0 = og; - if (c->nlin) - og = af_sum(S, og, linterms(S,c,1.)); - asl->P.dvsp0[i+1] = asl->P.ndvspout + Ncom; - dvi = asl->P.dv + i; - if (og0) { - c->cref = crefs(S); - dvi->ll = og; - dvi->nl = 0; - return; - } - if ((c->la = dvi->lt = afree(S, og, 0))) - dvi->scale = lt_scale; - for(n = 1, tl = tlist; tl; tl = tl->tnext) - n++; - dvi->ll = 0; - dvi->nl = tp = (linarg **)mem(n*sizeof(linarg*)); - for(tl = tlist; tl; tl = tl->tnext) - la_replace(S, *tp++ = tl); - *tp = 0; - c->cref = crefs(S); - } - - static int -colindvref(Static *S, expr *e, int ndv) -{ - expr **a, **ae; - int j, k, rv = 0; - - top: - switch(Intcast e->op) { - case f_OPPLUS: - case f_OPMINUS: - rv |= colindvref(S, e->R.e, ndv); - /* no break */ - case f_OPUMINUS: - e = e->L.e; - goto top; - case f_OPSUMLIST: - a = e->L.ep; - ae = e->R.ep; - while(a < ae) - rv |= colindvref(S, *a++, ndv); - break; - case f_OPVARVAL: - k = ((expr_v *)e)->a; - if ((k -= nv0x) < 0) - break; - if (zl[k]) { - rv |= zl[k]; - break; - } - zl[k] = 1; - if ((j = colindvref(S, (S->cexps + k)->e, k))) { - rv |= j; - zl[k] |= j; - } - break; - case f_OPMULT: - if (Intcast e->R.e->op == f_OPNUM) { - e = e->L.e; - goto top; - } - if (Intcast e->L.e->op == f_OPNUM) { - e = e->R.e; - goto top; - } - /* no break */ - default: - if (ndv >= 0) - rv = zl[ndv] |= 2; - } - return rv; - } - - static void -dv_walk(Static *S) -{ - int i, m, n; - expr_v *v; - ASLTYPE *asl = S->asl; - - m = asl->i.n_con0; - if ((n = Ncom)) { - for(i = 0; i < n_obj; i++) - colindvref(S, obj_de[i].e, -1); - for(i = 0; i < m; i++) - colindvref(S, con_de[i].e, -1); - larvlist = &v; - for(i = 0; i < n; i++) - dvwalk(S, i); - *larvlist = 0; - if ((i = asl->P.ndvspout)) - nv0b = nv0x + Ncom + i; - larvlist = 0; - } - } - - static int /* adjust zci, return k s.t. zci[i] < nv0 for i < k */ -nzcperm(Static *S) -{ - int i, j, k; - - k = nzc; - for(i = 0; i < k; ) - if (zci[i] >= nv0x) { - j = zci[--k]; - zci[k] = zci[i]; - zci[i] = j; - } - else - i++; - return k; - } - - static void -sumlist_adj(ASLTYPE *asl, expr *e, expr *e1) -{ - int k, n; - expr **ep, **ep0, **ep1; - - ep = e->R.ep; - k = htcl((n = ep - e->L.ep)*sizeof(expr*)); - if (n == (1 << k)) { - ep1 = (expr**)new_mblk(k+1); - memcpy(ep1, ep0 = e->L.ep, n*sizeof(expr*)); - del_mblk(k, ep0); - e->L.ep = ep1; - ep = ep1 + n; - } - *ep++ = e1; - e->R.ep = ep; - } - - static void -termwalk(Static *S, expr **ep, PSfind *p) -{ - ograd *og; - linarg **lap, **lap1; - linarg *tl; - int *cp, *cp0, *ce, *cee, i, j, k, kl, n, ncp, nzc2, *ui; - size_t len; - range *r; - list *L; - expr *e, *e1, *e2; - split_ce *cs; - psb_elem *b; - Elemtemp *bt; - ps_func *f; - ASLTYPE *asl = S->asl; - - Termno++; - tlist = 0; - afree(S, awalk(S, *ep), ep); - - i = k = nzcperm(S); - f = p->f; - if (!larvlist) - for(; i < nzc; i++) { - if ((j = zci[i]) < max_var) { - for(L = cexps[j-nv0x].cref; L; L = L->next) - if (!zc[L->item.i]++) - zci[nzc++] = L->item.i; - } - else { - cs = asl->P.Split_ce + (j-max_var); - if ((ce = cs->ce)) { - cee = ce + *ce; - while(ce++ < cee) - if (!zc[j = *ce + nv0x]++) - zci[nzc++] = j; - } - } - } - - r = (range *)rnz; /* scratch */ - if ((ncp = nzc - k)) { - zcsort(S, zc, zci+k, nv0x, ncp, -1 /*max_var*/); - cp = cp0 = (int*)(r+1); - i = k; - *cp = ncp; - do { - zc[j = zci[i++]] = 0; - *++cp = j - nv0x; - } while(i < nzc); - } - else - cp0 = 0; - - for(n = 0, tl = tlist; tl; tl = tl->tnext) { - n++; - for(og = tl->nz; og; og = og->next) - if (!zc[og->varno]++) - zci[k++] = og->varno; - } - nzc2 = k; - if (zc[-1]) - --nzc2; /* ignore constant */ - if (n <= 0 && Intcast (*ep)->op == f_OPNUM) - goto done; - - r->n = n; - r->nv = nzc2; - r->lap = lap1 = lap = (linarg**) - new_mblk(kl = htcl(n*sizeof(linarg*))); - for(tl = tlist; tl; tl = tl->tnext) - la_replace(S, *lap1++ = tl); - if (n > 1) - qsortv(lap, n, sizeof(linarg*), lacompar, NULL); - r->ui = ui = cp0 ? cp+1 : (int*)(r+1); - zcsort(S, zc, zci, 0, nzc2, nv0x); - for(j = 0; j < nzc2; j++) - *ui++ = zci[j]; - r = n >= nzc2 ? uhash(S,r) : rhash(S,r,1); - del_mblk(kl, lap); - - e1 = *ep; - if (!r || (i = r->lasttermno) == -1 || r->lastgroupno != Groupno) { - bt = p->b; - if ((i = f->nb++) >= bt->nmax) - upgrade_Elemtemp(S, bt); - b = f->b + i; - b->conno = Conno; - b->termno = i; - b->groupno = Groupno; - if ((b->U = r)) { - r->lasttermno = i; - r->lastgroupno = Groupno; - } - b->D.e = e1; - if (cp0) { - cp = (int*)mem(len = (ncp+1)*sizeof(int)); - memcpy(cp, cp0, len); - cp0 = cp; - } - b->ce = cp0; - while(k > 0) - zc[zci[--k]]= 0; - } - else { - b = f->b + i; - while(k > 0) - zc[zci[--k]] = 0; - if ((cp = b->ce) - && (*cp != ncp || memcmp(cp+1, cp0+1, ncp*sizeof(int)))) { - /* fix up f->ce */ - n = *cp; - while(k < n) - zc[zci[k++] = *++cp] = 1; - for(n = 0; n < ncp; ) - if (!zc[j = cp0[++n]]++) - zci[k++] = j; - qsortv(zci, k, sizeof(int), hscompar, S); - b->ce = cp = (int*)mem((k+1)*sizeof(int)); - *cp++ = k; - n = 0; - while(n < k) - zc[*cp++ = zci[n++]] = 0; - } - e = b->D.e; - switch(Intcast e->op) { - case f_OPPLUS: - ep = (expr**)new_mblk(2); - ep[0] = e->L.e; - ep[1] = e->R.e; - ep[2] = e1; - e->L.ep = ep; - e->R.ep = ep + 3; - e->op = (efunc *)f_OPSUMLIST; - break; - - case f_OPSUMLIST: - sumlist_adj(asl, e, e1); - break; - default: - b->D.e = e2 = new_expr(S, 0, e, e1); - e2->op = (efunc *)f_OPPLUS; - } - } - done: - nzc = 0; - } - - - static void -co_finish(ps_func *f) -{ - range *r; - psb_elem *b, *be; - psg_elem *g, *ge; - - b = f->b; - for(be = b + f->nb; b < be; b++) - if ((r = b->U)) - r->lasttermno = -1; - g = f->g; - for(ge = g + f->ng; g < ge; g++) - for(b = g->E, be = b + g->ns; b < be; b++) - if ((r = b->U)) - r->lasttermno = -1; - } - - static expr * -ecopy(Static *S, expr *e) -{ - expr **a, **a1, **ae; - int n, op; - - switch(op = Intcast e->op) { - case f_OPPLUS: - case f_OPMINUS: - e = new_expr(S, op, ecopy(S, e->L.e), ecopy(S, e->R.e)); - break; - case f_OPMULT: - if (Intcast e->L.e->op == f_OPNUM) { - e = new_expr(S, op, ecopy(S, e->R.e), (expr*) - new_expr_n(S, ((expr_n*)e->L.e)->v)); - break; - } - assert(Intcast e->R.e->op == f_OPNUM); - e = new_expr(S, op, ecopy(S, e->L.e), - (expr*)new_expr_n(S, ((expr_n*)e->R.e)->v)); - break; - case f_OPUMINUS: - e = new_expr(S, op, ecopy(S, e->L.e), 0); - break; - case f_OPSUMLIST: - a = e->L.ep; - ae = e->R.ep; - a1 = (expr**)new_mblk_ASL(S->a, htcl((n = ae-a)*sizeof(expr*))); - e = new_expr(S, op, (expr*)a1, (expr*)(a1+n)); - while(a < ae) - *a1++ = ecopy(S, *a++); - break; - case f_OPVARVAL: - break; -#ifdef DEBUG - default: - fprintf(Stderr, "Impossible case in ecopy!\n"); - exit(1); -#endif - } - return e; - } - - static int getgroup ANSI((Static *S, real scale, expr *e, PSfind *p)); - - static ograd * -ltermwalk(Static *S, real scale, expr **ep, PSfind *p) -{ - ASLTYPE *asl; - EU *en; - dv_info *dvi; - expr **a, **ae, *e; - int k, k0; - ograd *og, *rv; - cexp *ce; - - asl = S->asl; - rv = 0; - top: - e = *ep; - switch(Intcast e->op) { - case f_OPPLUS: - rv = af_sum(S, rv, ltermwalk(S, scale, &e->L.e, p)); - ep = &e->R.e; - goto top; - case f_OPMINUS: - rv = af_sum(S, rv, ltermwalk(S, scale, &e->L.e, p)); - ep = &e->R.e; - scale = -scale; - goto top; - case f_OPUMINUS: - ep = &e->L.e; - scale = -scale; - goto top; - case f_OPSUMLIST: - a = e->L.ep; - ae = e->R.ep; - while(a < ae) - rv = af_sum(S, rv, ltermwalk(S, scale, a++, p)); - break; - case f_OPVARVAL: - k = e->a; - if (k < nv0x) { - rv = af_sum(S, rv, new_ograd(S,0,k,scale)); - break; - } - k0 = k; - if ((k -= nv0x) >= Ncom) - goto dflt; - if (zl[k]) { - ce = cexps + k; - *ep = ecopy(S, ce->e); - if (ce->nlin) - rv = af_sum(S, rv, linterms(S,ce,scale)); - goto top; - } - if (!zc[e->a]++) - zci[nzc++] = e->a; - dvi = asl->P.dv + k; - if (dvi->nl) - goto dflt; - og = new_ograd(S, 0, k0, scale); - rv = af_sum(S, rv, og); - break; - case f_OPNUM: - rv = af_sum(S, rv, - new_ograd(S, 0, -1, scale*((expr_n *)e)->v)); - break; - case f_OPMULT: - en = (EU *)e->R.e; - if (Intcast en->op == f_OPNUM) { - *ep = e->L.e; - case_opnum: - if (en->u.v == 0.) { - efree(S,*ep); - *ep = (expr*)en; - } - else { - rv = af_sum(S, rv, ltermwalk(S, - scale*en->u.v, ep, p)); - en->u.p = expr_n_free; - expr_n_free = en; - } - e->L.e = expr_free; - expr_free = e; - break; - } - en = (EU *)e->L.e; - if (Intcast en->op == f_OPNUM) { - *ep = e->R.e; - goto case_opnum; - } - default: - dflt: - if (p->g && getgroup(S, scale, e, p)) - break; - if (scale != 1.) - *ep = scale == -1. - ? new_expr(S, f_OPUMINUS, e, 0) - : new_expr(S, f_OPMULT, e, - (expr *)new_expr_n(S, scale)); - termwalk(S, ep, p); - asl->P.ns0++; - } - return rv; - } - - static void -PSfind_init(Static *S, ps_func *f, PSfind *psf, int wantg) -{ - f->nxval = f->nb = f->ng = 0; - psf->f = f; - psf->b = new_Elemtemp(S, sizeof(psb_elem), (void**)&f->b); - if (wantg) - psf->g = new_Elemtemp(S, sizeof(psg_elem), (void**)&f->g); - else { - psf->g = 0; - f->g = 0; - last_psb_elem = psf->b; - } - } - - static void -psb_copy(psb_elem *b, psb_elem *b0, int n) -{ - range *r; - psb_elem *be; - - memcpy(b, b0, n*sizeof(psb_elem)); - for(be = b + n; b < be; b++) { - /* *b = *b0++; */ /* DEC Alpha gcc got this wrong */ - if (b->conno != -1 && (r = b->U)) { - b->next = r->refs; - r->refs = b; - } - } - } - - static int -getgroup(Static *S, real scale, expr *e, PSfind *p) -{ - ps_func *f, f1; - PSfind p1; - expr *e0, *e1, *e2; - ograd *og, *og1; - psb_elem *b, *be; - psg_elem *g; - Elemtemp *gt; - linarg **lap, **lape; - linpart *L; - range *U; - int i, nzc1, *zc1, *zci1; - ASL *asl = S->a; - - for(e1 = e; optype[Intcast e1->op] == 1; e1 = e1->L.e) - e0 = e1; - if (e1 == e || !might_expand(S, e1)) - return 0; - PSfind_init(S, &f1, &p1, 0); - f = p->f; - gt = p->g; - if ((i = f->ng++) >= gt->nmax) - upgrade_Elemtemp(S, gt); - Groupno = f->ng; - memset(g = f->g + i, 0, sizeof(psg_elem)); - g->g = e; - g->ge = e0; - g->scale = scale; - if ((og = ltermwalk(S, 1., &e0->L.e, &p1))) - og = compress(S, og, &g->g0, &i); - for(e1 = e; e1 != e0; e1 = e2) { - e2 = e1->L.e; - e2->R.e = e1; /* back pointer, used in psderprop */ - } - zc1 = S->zc1; - zci1 = S->zci1; - Groupno = nzc1 = 0; - if ((og1 = og)) { - for(i = 1; (og = og->next); i++); - g->nlin = i; - g->L = L = (linpart*)mem(i*sizeof(linpart)); - for(og = og1;; L++) { - zc1[zci1[nzc1++] = L->v.i = og->varno]= 1; - L->fac = og->coef; - if (!(og = og->next)) - break; - } - ogfree(S, og1); - } - g->esum.op = (efunc_n *)f_OPNUM; - g->ns = f1.nb; - g->E = b = (psb_elem*)mem(f1.nb * sizeof(psb_elem)); - psb_copy(b, f1.b, f1.nb); - del_Elemtemp(S, p1.b); - for(be = b + f1.nb; b < be; b++) { - if (!(U = b->U)) - continue; - lap = U->lap; - lape = lap + U->n; - while(lap < lape) - for(og = (*lap++)->nz; og; og = og->next) - if (!zc1[og->varno]++) - zci1[nzc1++] = og->varno; - } - zcsort(S, zc1, zci1, 0, nzc1, nv0x); - og = 0; - while(nzc1 > 0) { - og = new_ograd(S, og, i = zci1[--nzc1], 0.); - zc1[i] = 0; - } - g->og = og; - return 1; - } - - static ograd * -cotermwalk(Static *S, expr **ep, ps_func *f, int wantg, int omitdv) -{ - int comvar, n, x; - ograd *og; - real t; - PSfind psf; - psb_elem *b; - psg_elem *g, *g1, *ge; - - PSfind_init(S, f, &psf, wantg); - t = 0.; - og = ltermwalk(S, 1., ep, &psf); - if (omitdv && og) - og = compress(S, og, &t, &comvar); - b = f->b; - if (f->nb + f->ng == 0) - *ep = (expr*)new_expr_n(S, t); - else if (t) { - if (f->nb) - b->D.e = new_expr(S, f_OPPLUS, b->D.e, - (expr*)new_expr_n(S, t)); - else { - f->nb = 1; - memset(b, 0, sizeof(psb_elem)); - b->D.e = (expr *)new_expr_n(S, t); - } - } - co_finish(f); - if (!larvlist) { - n = f->nb; - if ((x = n * sizeof(psb_elem) + f->ng * sizeof(psg_elem))) { - ASL *asl = S->a; - g = (psg_elem*)(x >= 256 ? M1alloc(x) : mem(x)); - b = (psb_elem*)(g + f->ng); - if (n) - psb_copy(b, f->b, n); - else - b = 0; - if ((n = f->ng)) { - memcpy(g1 = g, f->g, n*sizeof(psg_elem)); - for(ge = g + n; g1 < ge; g1++) - g1->ge->L.e = (expr*)&g1->esum; - } - else - g = 0; - del_Elemtemp(S, psf.b); - if (wantg) - del_Elemtemp(S, psf.g); - f->b = b; - f->g = g; - } - } - return og; - } - - static void -psfind(Static *S) -{ - ASLTYPE *asl; - fint x; - int i, j, k, m, mx, nv, nx; - linarg *la, **lap; - ps_func *f, *f1; - range *r, *r0; - size_t L; -#ifdef PSHVREAD - Ihinfo *ihi, *ihi0, *ihi1; - int ihdlim, ihdlim1, ihdmax, ihdmin, n; -#endif - asl = S->asl; - m = asl->i.n_con0; - mx = m + (nx = asl->i.nsufext[ASL_Sufkind_con]); - x = (n_obj+mx)*sizeof(ps_func) - + slmax*sizeof(expr*) - + (Ncom+1)*sizeof(int); - asl->P.ops = f = (ps_func *)M1alloc(x); - asl->P.cps = f1 = f + n_obj; - if (nx) - memset(f1 + m, 0, nx*sizeof(ps_func)); - slscratch = (expr**)(f1 + mx); - asl->P.dvsp0 = (int*)(slscratch + slmax); - *asl->P.dvsp0 = Ncom; - - k = 3; - if ((nv = n_var) > 7) { - nv >>= 3; - do ++k; while(nv >>= 1); - } - S->klthash = S->krangehash = k; - lthash = (linarg**)new_mblkzap(asl, k); - rangehash = (range**)new_mblkzap(asl, k); - L = 1; - lthashmask = rangehashmask = (L << k) - 1; - nrange = 0; - - /* Separate allocation of zc1, zci1, rnz to allow freeing them later. */ - /* Allow room in rnz for use as scratch in termwalk. */ - i = max_var + (sizeof(range)+sizeof(int)+sizeof(real)-1)/sizeof(real); - rnz = (real *)M1alloc(i*(sizeof(real)+2*sizeof(int))); - S->zc1 = (int *)(rnz + i); - S->zci1 = S->zc1 + i; - memset(S->zc1, 0, max_var*sizeof(int)); - - Conno = -1; /* do not record refs for split defined vars */ - dv_walk(S); - asl->P.ndvspin = asl->P.ns0; - asl->P.ns0 = 0; - PSHV(n = 0); - for(i = 0; i < n_obj; i++, f++) { - Conno = -2 - i; - Ograd[i] = af_sum(S, Ograd[i], cotermwalk(S, &obj_de[i].e, f, - wantOgroups, 1)); - PSHV(n += f->ng); - } - PSHV(asl->P.nobjgroups = n); - PSHV(n = 0); - for(i = 0; i < m; i++, f1++) { - Conno = i; - Cgrad[i] = cf_sum(S, Cgrad[i], cotermwalk(S, &con_de[i].e, f1, - wantCgroups, 1)); - PSHV(n += f1->ng); - } - PSHV(asl->P.ncongroups = n); - - /* Compute nintv value for ranges */ - /* and set lasttermno to first var. */ - - r0 = (range*)&asl->P.rlist; - for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) { - i = 0; - if ((j = r->n) > 0) { - lap = r->lap; - r->lasttermno = (*lap)->nz->varno; - while(i < j) { - la = lap[i]; - if (la->v) { - i++; - la->termno = 0; /* for psgcomp */ - } - else { - lap[i] = lap[--j]; - lap[j] = la; - } - } - } - r->nintv = i; - } -#ifdef PSHVREAD - if ((ihdlim = ihd_limit) > 0) { - ihdmin = ihdlim1 = ihdlim + 1; - n = ihdlim1*sizeof(Ihinfo); - asl->P.ihi = ihi0 = (Ihinfo*)M1zapalloc(n); - ihdmax = 0; - for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) - if ((n = r->n) > 0) { - if (n > r->nv) - n = r->nv; - if (n > ihdlim) - n = ihdlim1; - else { - if (ihdmax < n) - ihdmax = n; - if (ihdmin > n) - ihdmin = n; - } - ihi = ihi0 + n - 1; - r->rlist.prev = ihi->r; - ihi->r = r; - ihi->nr++; - } - asl->P.ihdmax = ihdmax; - asl->P.ihdmin = asl->P.ndhmax = ihdmin; - ihi1 = ihi = ihi0 + ihdlim; - ihi->ihd = ihdlim1; /* sentinel */ - for(i = ihdlim; i > 0; --i) - if ((n = (--ihi)->nr)) { - ihi->next = ihi1; - ihi1 = ihi; - ihi->ihd = i; - ihi->k = htcl(((i*(i+1))>>1)*n*sizeof(real)); - } - asl->P.ihi1 = ihi1; - } -#endif - del_mblk(S->krangehash, rangehash); - del_mblk(S->klthash, lthash); - lthash = 0; - rangehash = 0; - } - -/*******/ - - static void -ewalk(Static *S, expr *e, int deriv) -{ - int a0, a1, i, j, j1, k, kf, numargs, op; - real *b, *b0, *ra; - unsigned int len; -#ifdef PSHVREAD - ASL *asl; - argpair *da; - int i1, j0, k0, k1, kn, *nn, *nn0; - real **fh; - expr **args1; -#endif - expr *L, *arg, **args, **argse; - expr_va *rva; - de *d; - derp *dsave; - expr_if *rvif; - expr_f *rvf; - arglist *al; - argpair *ap, *ap0, *ape; - char *dig; - - switch(optypeb[k = Intcast e->op]) { - - case 1: /* unary */ - e->dL = dvalue[k]; /* for UMINUS, FLOOR, CEIL */ - ewalk(S, e->L.e, deriv); - if (deriv) - dexpr(S, e, e->L.e, 0); - break; - - case 2: /* binary */ - e->dL = 1.; - e->dR = dvalue[k]; /* for PLUS, MINUS, REM */ - if (dvalue[k] == 11) - return; - ewalk(S, e->L.e, deriv); - ewalk(S, e->R.e, deriv); - if (deriv) - dexpr(S, e, e->L.e, e->R.e); - break; - - case 3: /* vararg (min, max) */ - rva = (expr_va *)e; - rva->next2 = varg2list; - varg2list = rva; - if (!last_d) { - new_derp(S, lasta, lasta, &edagread_one); - lasta++; - } - rva->d0 = dsave = last_d; -#ifdef PSHVREAD - rva->bak = last_e; -#endif - a0 = a1 = lasta; - j = 0; - for(d = rva->L.d; (L = d->e); d++) { - last_d = dsave; - op = Intcast L->op; - PSHV(last_e = 0;) - ewalk(S, L, deriv); - PSHV(d->ee = last_e;) - if (op == f_OPNUM || L->a == noa) { - d->d = dsave; - d->dv.i = noa; - } - else if (deriv) { - d->d = new_relo(S, L, dsave, &d->dv.i); - j++; - if (a1 < lasta) - a1 = lasta; - lasta = a0; - } - } - last_d = dsave; - if (j) { - rva->a = lasta = a1; - new_derp(S, 0, lasta++, &edagread_one); - /* f_MINLIST or f_MAXLIST will replace the 0 */ - rva->R.D = last_d; - nocopy = 1; -#ifdef PSHVREAD - last_e = (expr *)rva; - rva->dO.i = Hv_vararg; -#endif - } - else { - rva->a = noa; - rva->R.D = 0; - PSHV(last_e = rva->bak); - } - break; - - case 4: /* piece-wise linear */ - ewalk(S, e->R.e, deriv); /* should be a no-op */ - if (deriv) { - new_derp(S, e->R.e->a, e->a = lasta++, &e->dL); -#ifdef PSHVREAD - e->bak = last_e; - last_e = e; - e->dO.i = Hv_plterm; -#endif - } - break; - - case 5: /* if */ - rvif = (expr_if *)e; - rvif->next2 = if2list; - if2list = rvif; - PSHV(rvif->bak = last_e;) - ewalk(S, rvif->e, 0); - if (deriv && !last_d) { - new_derp(S, lasta, lasta, &edagread_one); - lasta++; - } - rvif->d0 = dsave = last_d; - a0 = lasta; - L = rvif->T; - op = Intcast L->op; - PSHV(last_e = 0;) - ewalk(S, L, deriv); - PSHV(rvif->Te = last_e;) - j = 0; - if (op == f_OPNUM || L->a == noa) { - rvif->dT = dsave; - rvif->Tv.i = noa; - } - else if ((j = deriv)) - rvif->dT = new_relo(S, L, dsave, &rvif->Tv.i); - a1 = lasta; - lasta = a0; - last_d = dsave; - L = rvif->F; - op = Intcast L->op; - PSHV(last_e = 0;) - ewalk(S, L, deriv); - PSHV(rvif->Fe = last_e;) - if (op == f_OPNUM || L->a == noa) { - rvif->dF = dsave; - rvif->Fv.i = noa; - } - else if ((j = deriv)) - rvif->dF = new_relo(S, L, dsave, &rvif->Fv.i); - if (lasta < a1) - lasta = a1; - last_d = dsave; - if (j) { - new_derp(S, 0, rvif->a = lasta++, &edagread_one); - rvif->D = last_d; - nocopy = 1; -#ifdef PSHVREAD - last_e = (expr *)rvif; - rvif->dO.i = Hv_if; -#endif - } - else { - rvif->a = noa; - rvif->D = 0; - PSHV(last_e = rvif->bak); - } - break; - - case 6: /* sumlist */ - a0 = lasta; - j = 0; - args = e->L.ep; - argse = e->R.ep; - e->a = lasta++; - while(args < argse) { - L = *args++; - op = Intcast L->op; - ewalk(S, L, deriv); - if (op != f_OPNUM && L->a != noa && deriv) { - new_derp(S, L->a, e->a, &edagread_one); - j++; - } - } - if (!j) { - e->a = noa; - lasta = a0; - } -#ifdef PSHVREAD - else { - e->bak = last_e; - last_e = e; - e->dO.i = Hv_sumlist; - j0 = e->R.ep - e->L.ep; - k0 = htcl(j0*sizeof(expr*)); - j1 = j0 + j + 1; - k1 = htcl(j1*sizeof(expr*)); - if (k1 > k0) { - ASL *asl = S->a; - args1 = (expr**)new_mblk(k1); - memcpy(args1, e->L.ep, - j0*sizeof(expr*)); - del_mblk(k0, e->L.ep); - e->L.ep = args1; - args = e->R.ep = args1 + j0; - } - argse = e->R.ep; - args1 = e->L.ep; - while(args1 < args) { - arg = *args1++; - if (arg->op != OPNUM && arg->a != noa) - *argse++ = arg; - } - *argse = 0; - } -#endif - break; - - case 7: /* function call */ - rvf = (expr_f *)e; - ap = rvf->sap; - for(ape = rvf->sape; ap < ape; ap++) - ewalk(S, ap->e,0); - ap = rvf->ap; - ape = rvf->ape; - kf = ape - ap; - al = rvf->al; - ra = al->ra; - numargs = al->nr; - for(i = 0; ap < ape; ap++) { - arg = ap->e; - op = Intcast arg->op; - ewalk(S, arg, deriv); - arg->op = (efunc *)(size_t)op; - if (arg->a != noa) - i += deriv; - } - rvf->a = i ? lasta++ : noa; - if (!numargs) { -#ifdef PSHVREAD - rvf->da = rvf->dae = 0; - rvf->fh= 0; /* for debugging */ -#endif - break; - } - dig = 0; - if (i) { - b = (real *)mem_ASL(S->a, len = -#ifdef PSHVREAD - (numargs*(numargs+3) >> 1) -#else - numargs -#endif - *sizeof(real)); - memset(b, 0, len); -#ifdef PSHVREAD - al->hes = b + numargs; -#endif - if (kf < numargs) { - ASLTYPE *asl = S->asl; - dig = (char*)mem(numargs); - for(j1 = 0; j1 < numargs; ++j1) - dig[j1] = 1; - } - } - else - b = 0; - al->derivs = b; - al->dig = dig; - b0 = b; -#ifdef PSHVREAD - asl = S->a; - if (i) { - da = (argpair*)mem(kf*sizeof(argpair) - + kf*kf*sizeof(real*)); - fh = (real**)(da + kf); - kn = htcl(2*sizeof(int)*numargs); - nn = (int*)new_mblk(kn); - } - else { - da = 0; - fh = 0; - kn = 12345; /* silence false "used uninitialized" warning */ - nn = 0; - } - rvf->da = da; - rvf->fh = fh; - nn0 = nn; -#endif - - for(ap = ap0 = rvf->ap; ap < ape; ap++) { - j1 = 0; - arg = ap->e; - arg->op = r_ops[op = Intcast arg->op]; - switch(op) { - case f_OPNUM: - goto loopend; - case f_OPVARVAL: - arg->op = (efunc *)(size_t)op; - } - b = b0 + (j = ap->u.v - ra); - *b = 0; - if (arg->a == noa) - goto loopend; - new_derp(S, arg->a, rvf->a, b); -#ifdef PSHVREAD - da->e = arg; - (da++)->u.v = b; - *nn++ = j; - j1 = 1; -#endif - loopend: - if (dig && j1) - dig[j] = 0; - } -#ifdef PSHVREAD - rvf->dae = da; - if (i) { - rvf->bak = last_e; - last_e = (expr *)rvf; - rvf->dO.i = Hv_func; - b = al->hes; - for(i1 = 0; i1 < kf; i1++) { - i = nn0[i1]; - for(j1 = 0; j1 < kf; j1++) { - j = nn0[j1]; - *fh++ = &b[i >= j - ? (i*(i+1)>>1)+j - : (j*(j+1)>>1)+i]; - } - } - } - if (nn0) - del_mblk(kn, nn0); -#endif - /* no break */ - - case 8: /* OPHOL (Hollerith) */ - case 9: /* OPNUM */ - break; - - case 10:/* OPVARVAL */ - k = (expr_v *)e - S->var_e; - if (k < 0 || k >= max_var) { - if ((k = ((expr_vx*)e)->a1) < 0) - return; - } - if (k >= 0 && deriv && !zc[k]++) - zci[nzc++] = k; - return; - - case 11: /* OPCOUNT */ - args = e->L.ep; - argse = e->R.ep; - while(args < argse) - ewalk(S, *args++, 0); - break; - - default: - scream(S->R, ASL_readerr_bug, - "unexpected optype[%d] = %d\n", k, optype[k]); - } -#ifdef DEBUG - if (e == opzork) - printf(""); -#endif - e->op = r_ops[k]; - } - - static list * -crefs(Static *S) -{ - int Nv0, Nzc, i, maxvar1 = S->max_var1; - list *rv = 0; - - Nv0 = nv0x; - for(Nzc = nzc; Nzc > 0; ) { - if ((i = zci[--Nzc]) >= Nv0 && i < maxvar1) { - rv = new_list(S->a, rv); - rv->item.i = i; - } - zc[i] = 0; - } - nzc = Nzc; - return rv; - } - - static funnel * -funnelfix(funnel *f) -{ - cexp *ce; - funnel *fnext, *fprev; - derp *d; - - for(fprev = 0; f; f = fnext) { - fnext = f->next; - f->next = fprev; - fprev = f; - ce = f->ce; - if ((d = ce->d)) - ce->z.i = d->b.i; - } - return fprev; - } - - static derp * -derpadjust(Static *S, derp *d0, int a, derp *dnext) -{ - derp *d, *d1; - int *r, *re; - relo *rl; - expr_if *il, *ile; - expr_va *vl, *vle; - de *de1; - ASL *asl; - - if (!(d = d0)) - return dnext; - asl = S->a; - r = imap + lasta0; - re = imap + lasta; - while(r < re) - *r++ = a++; - if (amax < a) - amax = a; - r = imap; - for(;; d = d1) { - d->a.i = r[d->a.i]; - d->b.i = r[d->b.i]; - if (!(d1 = d->next)) - break; - } - d->next = dnext; - if ((rl = relo2list)) { - relo2list = 0; - do { - d = rl->Dcond; - do { - d->a.i = r[d->a.i]; - d->b.i = r[d->b.i]; - } - while((d = d->next)); - } - while((rl = rl->next2)); - } - if (if2list != if2list_end) { - ile = if2list_end; - if2list_end = il = if2list; - do { - il->Tv.i = r[il->Tv.i]; - il->Fv.i = r[il->Fv.i]; - } - while((il = il->next2) != ile); - } - if (varg2list != varg2list_end) { - vle = varg2list_end; - varg2list_end = vl = varg2list; - do { - for(de1 = vl->L.d; de1->e; de1++) - de1->dv.i = r[de1->dv.i]; - } - while((vl = vl->next2) != vle); - } - return d0; - } - - static derp * -derpcopy(Static *S, cexp *ce, derp *dnext) -{ - derp *d, *dprev; - int *map; - derp d00; - - if (!(d = ce->d)) - return dnext; - map = imap; - for(dprev = &d00; d; d = d->next) { - new_derp(S, map[d->a.i], map[d->b.i], d->c.rp); - dprev = dprev->next = last_d; - } - dprev->next = dnext; - return d00.next; - } - - static derp * -derp_ogcopy(Static *S, ograd *og, derp *dnext, int k) -{ - last_d = dnext; - if (og->varno < 0) - og = og->next; - while(og) { - new_derp(S, og->varno, k, &og->coef); - og = og->next; - } - return last_d; - } - - static void -imap_alloc(Static *S) -{ - expr_v *v; - int i, *r; - linarg *tl; - ASLTYPE *asl = S->asl; - - if (imap) { - r = (int*)new_mblk(i = htcl(lasta*sizeof(int))); - memcpy(r, imap, imap_len*sizeof(int)); - del_mblk(kimap, imap); - imap = r; - kimap = i; - imap_len = (sizeof(char*)/sizeof(int)) << i; - return; - } - i = amax1 > lasta ? amax1 : lasta; - r = imap = (int*)new_mblk(kimap = htcl((i+100)*sizeof(int))); - imap_len = (sizeof(char*)/sizeof(int)) << kimap; - r += i = max_var1; - while(i > 0) - *--r = --i; - i = nv0x; - for(tl = asl->P.lalist; tl; tl = tl->lnext) - if ((v = tl->v)) - r[v->a] = i++; - r[noa] = i; - } - - static void -comsubs(Static *S, int alen, cde *d) -{ - list *L; - int a, i, j, jx, k, nzc1; - int *r, *re; - cexp *ce; - derp *D, *dnext; - dv_info *dvi; - relo *R; - expr *e; - split_ce *cs; - ograd *og; - ASLTYPE *asl = S->asl; - - D = last_d; - a = lasta00; - dnext = 0; - R = 0; - nzc1 = nzc; - nzc = 0; - for(i = j = 0; i < nzc1; i++) - if ((k = zci[i]) >= nv0x && k < max_var1) - zci[j++] = k; - else - zc[k] = 0; - if ((nzc1 = j)) { - for(i = 0; i < nzc1; i++) { - if ((j = zci[i] - nv0x) >= Ncom) { - cs = asl->P.Split_ce + (j-Ncom); - if ((r = cs->ce)) { - re = r + *r; - while(r++ < re) - if (!zc[j = *r + nv0x]++) - zci[nzc1++] = j; - } - } - else if (j >= 0 && (L = cexps[j].cref)) { - if (cexps[j].funneled) do { - if (zc[j = L->item.i] - || asl->P.dv[j-nv0x].ll) - continue; - zc[j] = 1; - zci[nzc1++] = j; - } - while((L = L->next)); - else do { - if (!zc[L->item.i]++) - zci[nzc1++] = L->item.i; - } - while((L = L->next)); - } - } - if (nzc1 > 1) - zcsort(S, zc, zci, 0, nzc1, -1); - R = new_relo1(S, dnext); - i = 0; - do { - j = zci[i]; - jx = j - nv0x; - zc[j] = 0; - ce = &cexps[jx]; - e = ce->e; - k = varp[jx]->a; - if (ce->funneled) { - if (j >= max_var) - imap[((expr_vx*)varp[jx])->a0] = a; - imap[k] = a++; - } - else { - r = imap + ce->z.i; - re = r + ce->zlen; - while(r < re) - *r++ = a++; - /* zlen == 1 if e->op == OPNUM */ - imap[k] = e->op == OPNUM ? a-1 : imap[e->a]; - if (!ce->d - && jx < Ncom - && (og = (dvi = &asl->P.dv[jx])->ll) - && (dvi->nl || dvi->lt)) { - dnext = R->D = - derp_ogcopy(S,og,R->D,imap[j]); - continue; - } - } - dnext = R->D = derpcopy(S, ce, R->D); - } - while(++i < nzc1); - } - if (D || R) { - if (!R) - R = new_relo1(S, dnext); - D = R->D = derpadjust(S, D, a, R->D); - if (Intcast d->e->op != f_OPVARVAL) - d->e->a = imap[d->e->a]; - } - d->d = D; - a += alen; - d->zaplen = (a > lasta00 ? a - nv1 : 0)*sizeof(real); - if (amax < a) - amax = a; - } - - static void -co_read(EdRead *R, cde *d, int k) -{ - d = (cde *)((char *)d + k*sizeof(cde)); - d->e = eread(R); - } - - static void -co_walk(Static *S, cde *d) -{ - int alen; - - if (amax1 < lasta) - amax1 = lasta; - lasta = lasta0; - last_d = 0; - PSHV(last_e = 0;) - ewalk(S, d->e, 1); - PSHV(d->ee = last_e;) - alen = lasta - lasta0; - if (imap_len < lasta) - imap_alloc(S); - comsubs(S, alen, d); - } - - static int -lpcompar(const void *a, const void *b, void *v) -{ - Not_Used(v); - return ((linpart *)a)->v.i - ((linpart *)b)->v.i; - } - - static linpart * -linpt_read(EdRead *R, int nlin) -{ - ASL *asl; - int i0, needsort; - int (*Xscanf)(EdRead*, const char*, ...); - linpart *L, *rv; - - asl = R->asl; - if (nlin <= 0) - return 0; - Xscanf = xscanf; - L = rv = (linpart *)new_mblk(htcl(nlin*sizeof(linpart))); - i0 = needsort = 0; - do { - if (Xscanf(R, "%d %lf", &L->v.i, &L->fac) != 2) - badline(R); - if (i0 > L->v.i) - needsort++; - i0 = L->v.i; - L++; - } - while(--nlin > 0); - if (needsort) - qsortv(rv, L-rv, sizeof(linpart), lpcompar, NULL); - return rv; - } - - static int -funnelkind(Static *S, int i0, int *ip) -{ - cexp *ce, *cej, *ce2; - dv_info *dvi; - int i, j, k, k1, k2, nzc0, rv; - int *vr, *vre; - linarg *la, **lap, **lape; - linpart *L, *Le; - list *tl; - range *r; - ASLTYPE *asl = S->asl; - - ce = cexps + i0; - ce->vref = 0; - if (!(nzc0 = nzc)) - return nocopy; - k1 = k2 = nzcperm(S); - dvi = asl->P.dv; - if (i0 >= Ncom || !dvi[i0].nl) - dvi = 0; - ce2 = 0; /* used to get funnelkind 2 right in rare cases */ - for(i = k = 0; i < nzc; i++) - if ((j = zci[i]) < nv0x) { - if (k >= maxfwd) { - k = 0; - break; - } - vrefx[k++] = j; - } - else { - cej = cexps + (j -= nv0x); - if (!(vr = cej->vref)) { - if (dvi && j < Ncom && (tl = cej->cref)) { - if (dvi[j].ll) { - for(; tl; tl = tl->next) - if (!zc[tl->item.i]++) - zci[nzc++] = tl->item.i; - } - else { - cej->vref = (int*)ce2; - ce2 = cej; - } - } - continue; - } - if (!*++vr) { - k = 0; - break; - } - vre = vr + *vr; - while(++vr <= vre) - if (!zc[*vr]++) - zci[nzc++] = *vr; - } - if (k2 < k) - k2 = k; - k2 += 2; - if (k2 > nvref) { - nvref = (maxfwd + 1)*(ncom_togo < vrefGulp - ? ncom_togo : vrefGulp); - if (nvref < k2) - nvref = 2*k2; - vrefnext = (int *)M1alloc(nvref*Sizeof(int)); - } - vr = ce->vref = vrefnext; - vrefnext += k2; - nvref -= k2; - vr[0] = k1; - vr[1] = k; - vr += 2; - i = 0; - if (k) - while(i < k) - *vr++ = vrefx[i++]; - else - while(i < k1) - *vr++ = zci[i++]; - while(nzc > nzc0) - zc[zci[--nzc]] = 0; - if (!nocopy) { - k1 = 0; - if (i0 < Ncom) { - if ((lap = (asl->P.dv + i0)->nl)) - while((la = *lap++)) - if (la->nnz > 1) - k1++; - } - else { - r = asl->P.Split_ce[i0-Ncom].r; - lap = r->lap; - lape = lap + r->nintv; - while(lap < lape) - if ((*lap++)->nnz > 1) - k1++; - } - if (k) { - if (nderp > 3*(k+k1)) { - if (!ce2) - goto ret2; - while((cej = ce2)) { - ce2 = (cexp*)cej->vref; - cej->vref = 0; - for(tl = cej->cref; tl; tl = tl->next) - if (!zc[tl->item.i]++) - zci[nzc++] = tl->item.i; - } - for(i = k = 0; i < nzc; i++) { - if ((j = zci[i]) < nv0x) { - vrefx[k++] = j; - continue; - } - cej = cexps + (j -= nv0x); - if ((vr = cej->vref)) { - if (*++vr) { - for(vre = vr + *vr; ++vr < vre; ) - if (!zc[*vr]++) - zci[nzc++] = *vr; - } - } - if (j < Ncom) { - for(tl = cej->cref; tl; tl = tl->next) - if (!zc[tl->item.i]++) - zci[nzc++] = tl->item.i; - } - if ((L = cej->L)) { - for(Le = L + cej->nlin; L < Le; ++L) - if (!zc[L->v.i]++) - zci[nzc++] = L->v.i; - } - } - while(nzc > nzc0) - zc[zci[--nzc]] = 0; - ret2: - *ip = k; - return 2; - } - } - } - while((cej = ce2)) { - ce2 = (cexp*)cej->vref; - cej->vref = 0; - } - rv = 0; - if (nocopy || nderp > 3*(nzc0+k1)) - rv = 1; - return rv; - } - - static void -cexp_read(EdRead *R, int k, int nlin) -{ - cexp *ce; - expr *e, *zero; - Static *S = (Static *)R->S; - ASLTYPE *asl = S->asl; - - ce = cexps + k - nv0x; - ce->L = linpt_read(R, ce->nlin = nlin); - e = eread(R); - if (e->op == (efunc *)f_OPVARVAL) { - /* avoid confusion with very simple defined variables */ - zero = (expr*)new_expr_n(S, 0.); - zero->op = (efunc*)f_OPNUM; - e = new_expr(S, 0, e, zero); - } - ce->e = e; - } - - static cplist * -funnelderp(Static *S, int a, cplist *cl0) -{ - cplist *cl; - ASL *asl = S->a; - - new_derp(S, a, lasta0, 0); - cl = (cplist *)mem(sizeof(cplist)); - cl->next = cl0; - cl->ca.i = imap[a]; - cl->cfa = last_d->c.rp = (real *)mem(sizeof(real)); - return cl; - } - - static void -cexp_walk(Static *S, int k) -{ - ASLTYPE *asl = S->asl; - cexp *ce; - cplist *cl; - dv_info *dvi; - expr *e; - funnel *f, **fp; - int a, *ap, fk, i, j, ka, la0, na, nlin; - linarg *la, **lap, **lape; - linpart *L, *Le; - ograd *og; - range *r; - - /* for now, assume all common exprs need to be differentiated */ - - ce = cexps + k; - nlin = ce->nlin; - L = ce->L; - nocopy = 0; - last_d = 0; - la0 = lasta; - nderps += nderp; - nderp = 0; - e = ce->e; - switch(Intcast e->op) { - case f_OPVARVAL: - case f_OPNUM: - ap = 0; - break; - default: - ap = &e->a; - } - PSHV(last_e = 0;) - ewalk(S, e, 1); - PSHV(ce->ee = last_e;) - ka = k + nv0x; - if ((na = lasta - la0)) - ce->z.i = la0; - else { - ce->z.i = k < Ncom ? ka : ((expr_vx*)varp[k-Ncom])->a0; - na = 1; - } - ce->zlen = na; - j = ap ? *ap : ka; - if (nlin) { - if (nlin == 1) - new_derp(S, L->v.i, j, &L->fac); - else if (k < Ncom) { - dvi = asl->P.dv + k; - if (dvi->lt) - new_derp(S, dvi->lt->v->a, j, &dvi->scale); - } - for(Le = L + nlin; L < Le; L++) - if (!zc[i = L->v.i]++) - zci[nzc++] = i; - } - i = 0; /* only needed to shut up an erroneous warning */ - if ((fk = funnelkind(S, k, &i))) { - /* arrange to funnel */ - fp = ka < nv0b ? &f_b : ka < nv0c ? &f_c : &f_o; - ce->funneled = f = (funnel *)mem(sizeof(funnel)); - f->next = *fp; - *fp = f; - f->ce = ce; - if (imap_len < lasta) - imap_alloc(S); - if (fk == 1) { - f->fulld = last_d; - a = lasta00; - for(i = nzc; --i >= 0; ) - if ((j = zci[i]) >= nv0x && j < max_var1) - imap[varp[j-nv0x]->a] = a++; - na += a - nv1; - f->fcde.zaplen = na*sizeof(real); - i = nzc; - derpadjust(S, last_d, a, 0); - } - else { - f->fulld = 0; - f->fcde.e = e; - comsubs(S, na, &f->fcde); - memcpy(zci, vrefx, i*sizeof(int)); - } - last_d = 0; - cl = 0; - while(i > 0) { - if ((a = zci[--i]) >= nv0x && a < max_var1) - a = varp[a-nv0x]->a; - if (a != noa) - cl = funnelderp(S, a, cl); - } - if (k < Ncom) { - if ((lap = (asl->P.dv + k)->nl)) - while((la = *lap++)) - if (la->v) - cl = funnelderp(S, la->v->a, cl); - } - else { - r = asl->P.Split_ce[k-Ncom].r; - lap = r->lap; - lape = lap + r->nintv; - while(lap < lape) - if ((la = *lap++)->nnz > 1) - cl = funnelderp(S, la->v->a, cl); - } - f->cl = cl; - varp[k]->a = lasta0++; - lasta = lasta0; - } - lasta0 = lasta; - if (!(ce->d = last_d) - && e->op == OPNUM - && (og = asl->P.dv[k].ll) - && og->varno < 0) - ((expr_n*)e)->v = 0; - while(nzc > 0) - zc[zci[--nzc]] = 0; - --ncom_togo; - } - -#ifdef PSHVREAD - static expr * -hvadjust(expr *e) -{ - expr *e0; - - for(e0 = 0; e; e = e->bak) { - e->fwd = e0; - e0 = e; - e->a = e->dO.i; - } - return e0; - } - - static void -co_adjust(ps_func *p, int n) -{ - ps_func *pe; - psb_elem *b, *be; - psg_elem *g, *ge; - - for(pe = p + n; p < pe; p++) { - for(b = p->b, be = b + p->nb; b < be; b++) - b->D.ef = hvadjust(b->D.ee); - for(g = p->g, ge = g + p->ng; g < ge; g++) { - g->esum.op = (efunc_n*)OPNUM; - for(b = g->E, be = b + g->ns; b < be; b++) - b->D.ef = hvadjust(b->D.ee); - } - } - } -#else /*PSHVREAD*/ -#define co_adjust(a,b) /*nothing*/ -#endif /*PSHVREAD*/ - - static void -ifadjust(ASLTYPE *asl, expr_if *e) -{ - real *Adjoints = asl->i.adjoints_; - for(; e; e = e->next) { - e->Tv.rp = &Adjoints[e->Tv.i]; - e->Fv.rp = &Adjoints[e->Fv.i]; -#ifdef PSHVREAD - e->Tf = hvadjust(e->Te); - e->Ff = hvadjust(e->Fe); -#endif - } - } - - static void -vargadjust(ASLTYPE *asl, expr_va *e) -{ - de *d; - real *Adjoints = asl->i.adjoints_; - - for(; e; e = e->next) { - for(d = e->L.d; d->e; d++) { - d->dv.rp = &Adjoints[d->dv.i]; -#ifdef PSHVREAD - d->ef = hvadjust(d->ee); -#endif - } - } - } - - static void -funneladj1(ASLTYPE *asl, funnel *f) -{ - real *a = adjoints; - derp *d; - cplist *cl; - funnel *fnext; - - for(; f; f = fnext) { - fnext = f->next; - f->next = 0; - if ((d = f->fulld)) { - f->fcde.d = d; - do { - d->a.rp = &a[d->a.i]; - d->b.rp = &a[d->b.i]; - } - while((d = d->next)); - } - for(cl = f->cl; cl; cl = cl->next) - cl->ca.rp = &a[cl->ca.i]; - } - } - - static void -funneladjust(Static *S) -{ - cexp *c, *ce; - linpart *L, *Le; - ASLTYPE *asl = S->asl; - - c = cexps; - for(ce = c + Ncom; c < ce; c++) { - if ((L = c->L)) - for(Le = L + c->nlin; L < Le; L++) - L->v.vp = (void*)&var_e[L->v.i]; -#ifdef PSHVREAD - c->ef = hvadjust(c->ee); -#endif - } -#ifdef PSHVREAD - for(ce += asl->P.ndvspout; c < ce; c++) - c->ef = hvadjust(c->ee); -#endif - - funneladj1(asl, f_b); - funneladj1(asl, f_c); - funneladj1(asl, f_o); - } - - static void -cg_zzap(ASLTYPE *asl) -{ - cgrad *cg, **cgp,**cgp1, **cgpe; - - cgp1 = Cgrad; - cgpe = cgp1 + asl->i.n_con0; - while(cgp1 < cgpe) { - cgp = cgp1++; - while((cg = *cgp)) - if (cg->coef) - cgp = &cg->next; - else - *cgp = cg->next; - } - } - - static void -adjust_compl_rhs(ASLTYPE *asl) -{ - cde *C; - expr *e; - int *Cvar, i, j, nc, stride; - real *L, *U, t, t1; - - L = LUrhs; - if ((U = Urhsx)) - stride = 1; - else { - U = L + 1; - stride = 2; - } - C = con_de; - Cvar = cvar; - nc = n_con; - for(i = nlc; i < nc; i++) - if (Cvar[i] && (e = C[i].e) && Intcast e->op == f_OPNUM - && (t = ((expr_n*)e)->v) != 0.) { - t1 = t; - if (L[j = stride*i] > negInfinity) { - L[j] -= t; - t1 = 0.; - } - if (U[j] < Infinity) { - U[j] -= t; - t1 = 0.; - } - ((expr_n*)e)->v = t1; - } - } - - static void -adjust(Static *S, int flags) -{ - derp *d, **dp; - ASLTYPE *asl = S->asl; - real *a = adjoints; - relo *r; - - for(r = relolist; r; r = r->next) { - dp = &r->D; - while((d = *dp)) { - d->a.rp = &a[d->a.i]; - d->b.rp = &a[d->b.i]; - dp = &d->next; - } - *dp = r->Dnext; - } - ifadjust(asl, iflist); - vargadjust(asl, varglist); - if (Ncom) - funneladjust(S); - co_adjust(asl->P.cps, asl->i.n_con0); - co_adjust(asl->P.ops, n_obj); - if (asl->i.n_con0 && !allJ) - cg_zzap(asl); - if (k_seen) { - if (!A_vals) - goff_comp_ASL((ASL*)asl); - else if (Fortran) - colstart_inc_ASL((ASL*)asl); - } - if (n_cc > nlcc && nlc < n_con - && !(flags & ASL_no_linear_cc_rhs_adjust)) - adjust_compl_rhs(asl); - } - - static void -br_read(EdRead *R, int nc, real *L, real *U, int *Cvar, int nv) -{ - ASL *asl; - int i, inc, j, k; - int (*Xscanf)(EdRead*, const char*, ...); - - asl = R->asl; - if (U) - inc = 1; - else { - U = L + 1; - inc = 2; - } - Xscanf = xscanf; - Xscanf(R, ""); /* purge line */ - for(i = 0; i < nc; i++, L += inc, U += inc) { - switch(edag_peek(R) - '0') { - case 0: - if (Xscanf(R,"%lf %lf",L,U)!= 2) - badline(R); - break; - case 1: - if (Xscanf(R, "%lf", U) != 1) - badline(R); - *L = negInfinity; - break; - case 2: - if (Xscanf(R, "%lf", L) != 1) - badline(R); - *U = Infinity; - break; - case 3: - *L = negInfinity; - *U = Infinity; - Xscanf(R, ""); /* purge line */ - break; - case 4: - if (Xscanf(R, "%lf", L) != 1) - badline(R); - *U = *L; - break; - case 5: - if (Cvar) { - if (Xscanf(R, "%d %d", &k, &j) == 2 - && j > 0 && j <= nv) { - Cvar[i] = j; - *L = k & 2 ? negInfinity : 0.; - *U = k & 1 ? Infinity : 0.; - break; - } - } - default: - badline(R); - } - } - } - - static expr * -aholread(EdRead *R) -{ - int i, k; - expr_h *rvh; - char *s1; - FILE *nl = R->nl; - Static *S = (Static *)R->S; - - k = getc(nl); - if (k < '1' || k > '9') - badline(R); - i = k - '0'; - while((k = getc(nl)) != ':') { - if (k < '0' || k > '9') - badline(R); - i = 10*i + k - '0'; - } - rvh = (expr_h *)mem_ASL(R->asl, memadj(sizeof(expr_h) + i)); - for(s1 = rvh->sym;;) { - if ((k = getc(nl)) < 0) { - fprintf(Stderr, - "Premature end of file in aholread, line %ld of %s\n", - R->Line, R->asl->i.filename_); - exit_ASL(R,1); - } - if (k == '\n') { - R->Line++; - if (!i) - break; - } - if (--i < 0) - badline(R); - *s1++ = k; - } - *s1 = 0; - rvh->op = (efunc *)f_OPHOL; - rvh->a = noa; - return (expr *)rvh; - } - - static expr * -bholread(EdRead *R) -{ - int i; - expr_h *rvh; - char *s; - FILE *nl = R->nl; - ASL *asl = R->asl; - Static *S = (Static *)R->S; - - if (xscanf(R, "%d", &i) != 1) - badline(R); - rvh = (expr_h *)mem(memadj(sizeof(expr_h) + i)); - s = rvh->sym; - if (fread(s, i, 1, nl) != 1) - badline(R); - s[i] = 0; - rvh->op = (efunc *)f_OPHOL; - rvh->a = noa; - for(;;) switch(*s++) { - case 0: goto break2; /* so we return at end of fcn */ - case '\n': R->Line++; - } - break2: - return (expr *)rvh; - } - - static int -qwalk(Static *S, expr *e) -{ - int i, j, k; - expr **args, **argse; - expr_f *rvf; - argpair *ap, *ape; - - top: - switch(optype[k = Intcast e->op]) { - - case 1: /* unary */ - switch(k) { - case f_OPCPOW: - if (qwalk(S, e->R.e)) - return 3; - return 0; - case f_OPUMINUS: - e = e->L.e; - goto top; - case f_OP2POW: - i = qwalk(S, e->L.e); - if (i < 2) - return i << 1; - } - return 3; - - case 2: /* binary */ - switch(k) { - case f_OPPLUS: - case f_OPMINUS: - i = qwalk(S, e->L.e); - if (i == 3) - break; - j = qwalk(S, e->R.e); - if (i < j) - i = j; - return i; - case f_OPDIV: - if (qwalk(S, e->R.e)) - return 3; - e = e->L.e; - goto top; - case f_OPMULT: - i = qwalk(S, e->L.e); - if (i < 3) { - i += qwalk(S, e->R.e); - if (i < 3) - return i; - } - } - return 3; - - case 6: /* sumlist */ - j = 0; - args = e->L.ep; - argse = e->R.ep; - while(args < argse) { - i = qwalk(S, *args++); - if (j < i) { - j = i; - if (j == 3) - break; - } - } - return j; - - case 7: /* function call */ - rvf = (expr_f *)e; - ap = rvf->ap; - ape = rvf->ape; - for(i = 0; ap < ape; ap++) { - if (qwalk(S, ap->e)) - return 3; - } - - case 9: /* OPNUM */ - return 0; - case 10:/* OPVARVAL */ - k = (expr_v *)e - S->var_e; - if (k < 0) - goto k_adj; - if (k < nv0x) - return 1; - if (k >= max_var) { - k_adj: - k = ((expr_vx*)e)->a1; - return k < 0 ? 1 : S->asl->I.v_class[k-nv0x]; - } - return S->asl->I.v_class[k - nv0x]; - } - return 3; - } - - static int -co_walkloop(Static *S, ps_func *f, int n, char *c, ograd **o) -{ - ps_func *fe; - psb_elem *b, *be; - psg_elem *g, *ge; - int k, k1, kx; - - kx = 0; - for(fe = f + n; f < fe; f++) { - if (c) { - k = *o++ != 0; - for(g = f->g, ge = g + f->ng; g < ge; g++) { - if (Intcast g->g->op != f_OP2POW) { - k = 3; - goto have_c; - } - if (g->nlin) - k = 2; - for(b = g->E, be = b + g->ns; b < be; b++) { - if ((k1 = qwalk(S, b->D.e)) > 1) { - k = 3; - goto have_c; - } - k = 2; - } - } - for(b = f->b, be = b + f->nb; b < be; b++) { - if ((k1 = qwalk(S, b->D.e)) > k) { - k = k1; - if (k == 3) - goto have_c; - } - } - have_c: *c++ = (char)k; - if (kx < k) - kx = k; - } - for(b = f->b, be = b + f->nb; b < be; b++) - co_walk(S, &b->D); - for(g = f->g, ge = g + f->ng; g < ge; g++) { - ewalk(S, g->g, 1); - for(b = g->E, be = b + g->ns; b < be; b++) - co_walk(S, &b->D); - } - } - return kx; - } - - static void -do_ewalk(Static *S) -{ - int i, *map, n, nc; - char *v, *ve; - expr_v *e, *ee, **vp, **vpe; - efunc *vv; - linarg *la, **lap; - ograd *og; - cexp *c; - ASLTYPE *asl = S->asl; - - psfind(S); - nc = max_var1 - nv0x; - noa = lasta00 + nc; - nv1 = lasta00++; - amax = lasta = lasta0 = max_var + nndv + 1; - - if (Ncom && asl->I.o_class) { - i = max_var1 - nv0x; - v = asl->I.v_class = (char*)M1alloc(i); - ve = v + i; - c = cexps; - while(v < ve) - *v++ = (char)((i = qwalk(S, (c++)->e)) ? i : 1); - } - - n = Ncom + asl->P.ndvspout; - while(nzc > 0) - zc[zci[--nzc]] = 0; - for(i = 0; i < n; i++) { -#ifdef DEBUG - if (i == izork) - printf(""); -#endif - cexp_walk(S, i); - } - - if (imap_len < lasta) - imap_alloc(S); - f_b = funnelfix(f_b); - f_c = funnelfix(f_c); - f_o = funnelfix(f_o); - - asl->I.c_class_max = co_walkloop(S, asl->P.cps, asl->i.n_con0, - asl->I.c_class, (ograd**)Cgrad); - asl->I.o_class_max = co_walkloop(S, asl->P.ops, n_obj, - asl->I.o_class, Ograd); - del_mblk(kzc, zci); - zci = zc = 0; /* for debugging */ - e = var_e; - vv = OPVARVAL; - for(ee = e + nv0x + Ncom; e < ee; e++) - e->op = vv; - lap = &asl->P.lalist; - map = imap; - while((la = *lap)) - if ((e = la->v)) { - e->op = vv; - e->a = map[e->a]; - lap = &la->lnext; - } - else { - og = la->nz; - assert(og != 0); - assert(og->next == 0); - assert(og->varno < nv0x); - la->v = var_e + og->varno; - *lap = la->lnext; - } - if (n) { - asl->P.vp = vp = varp; -#ifdef PSHVREAD - asl->P.zlsave = zl; -#endif - e = var_ex; - for(ee = e + Ncom; e < ee; e++) - *vp++ = e; - vpe = vp + asl->P.ndvspout; - while(vp < vpe) - (*vp++)->op = vv; - } - } - - int -pfg_read_ASL(ASL *a, FILE *nl, int flags) -{ - ASLTYPE *asl; - EdRead ER, *R; - Jmp_buf JB; - Static SS, *S; - cgrad *cg, **cgp; - char fname[128]; - expr_v *e; - func_info *fi; - int allG, i, i1, j, k, *ka, kseen, maxfwd1, nc, nc0, nco, nlin, nlogc, no; - int nv, nvc, nvo, nvr, nvextra, nvx, readall; - int (*Xscanf)(EdRead*, const char*, ...); - ograd *og, **ogp; - real t; - size_t *kaz, nz; - unsigned x; - - ASL_CHECK(a, asltype, who); - flagsave_ASL(a, flags); /* includes allocation of LUv, LUrhs, A_vals or Cgrad, etc. */ - asl = (ASLTYPE*)a; - S = S_init(&SS, asl); - ed_reset(asl); - SS.R = R = EdReadInit_ASL(&ER, a, nl, S); - if (flags & ASL_return_read_err) { - a->i.err_jmp_ = &JB; - i = setjmp(JB.jb); - if (i) { - a->i.err_jmp_ = 0; - return i; - } - } - if ((nlogc = a->i.n_lcon_) && !(flags & ASL_allow_CLP)) { - if (a->i.err_jmp_) - return ASL_readerr_CLP; - sorry_CLP(R, "logical constraints"); - } - if (!(flags & ASL_find_default_no_groups)) - flags |= ASL_findgroups; - Xscanf = xscanf; - readall = flags & ASL_keep_all_suffixes; - PSHV(asl->P.pshv_g1 = 1;) - k_Elemtemp = htcl(sizeof(Elemtemp)); - allJ = (flags & ASL_J_zerodrop) == 0; - allG = (flags & ASL_G_zerodrop) == 0; - wantCgroups = flags & ASL_findCgroups; - wantOgroups = flags & ASL_findOgroups; - OPNUM = r_ops[f_OPNUM]; - OPVARVAL = r_ops[f_OPVARVAL]; - if (!size_expr_n) - size_expr_n = sizeof(expr_n); - S->size_exprn = size_expr_n; - asl->P.rlist.next = asl->P.rlist.prev = (range*)&asl->P.rlist; - if (nfunc) - func_add(a); - if (binary_nl) - holread = bholread; - else - holread = aholread; - - asl->P.ncom = Ncom = comb + comc + como + comc1 + como1; - nc0 = n_con; - nc = nc0 + a->i.nsufext[ASL_Sufkind_con]; - no = n_obj; - nvc = c_vars; - nvo = o_vars; - nco = nc + no + nlogc; - if (no < 0 || nco <= 0) - scream(R, ASL_readerr_corrupt, - "pshvread: nc = %d, no = %d, nlogc = %d\n", - nc0, no, nlogc); - if (pi0) { - memset(pi0, 0, nc*sizeof(real)); - if (havepi0) - memset(havepi0, 0, nc); - } - nvextra = a->i.nsufext[ASL_Sufkind_var]; - nv0r = nvr = n_var; - asl->P.nv0_ = lasta00 = nv0x = nvx = n_var + nvextra; - max_var1 = max_var = nv = nvr + Ncom + nvextra; - combc = comb + comc; - ncom0 = ncom_togo = combc + como; - nzclim = max_var >> 3; - ncom1 = comc1 + como1; - nv0b = nvr + comb; - nv0c = nv0b + comc; - if ((maxfwd1 = maxfwd + 1) > 1) - nvref = maxfwd1*((ncom0 < vrefGulp ? ncom0 : vrefGulp) + 1); - x = nco*sizeof(cde) + no*sizeof(ograd *) - + nv*sizeof(expr_v) - + nfunc*sizeof(func_info *) - + nvref*sizeof(int) - + no; - SS.nvar0 = a->i.n_var0; - if (!(SS.nvinc = a->i.nvinc = a->i.n_var_ - SS.nvar0 + nvextra)) - SS.nvar0 += ncom0 + ncom1; - if (flags & ASL_find_co_class) - x += nco; - if (X0) - memset(X0, 0, nvr*sizeof(real)); - if (havex0) - memset(havex0, 0, nvr); - e = var_e = (expr_v *)M1zapalloc(x); - con_de = (cde *)(e + nv); - lcon_de = con_de + nc; - obj_de = lcon_de + nlogc; - Ograd = (ograd **)(obj_de + no); - var_ex = e + nvx; - var_ex1 = var_ex + ncom0; - for(k = 0; k < nv; e++) { - e->op = (efunc *)f_OPVARVAL; - e->a = k++; - } - funcs = (func_info **)(Ograd + no); - zci = 0; - zc_upgrade(S); /* initially allow nv+1 entries in zci and zc[-1] */ - vrefx = (int *)(funcs + nfunc); - if (Ncom) { - cexps = 0; - asl->P.ndvspout = 0; - memset(asl->P.dv = (dv_info*)mem(Ncom*sizeof(dv_info)), - 0, Ncom*sizeof(dv_info)); - cexp_upgrade(S, Ncom); - for(k = 0; k < Ncom; k++) - varp[k] = &var_ex[k]; - } - objtype = (char *)(vrefx + nvref); - if (flags & ASL_find_co_class) { - asl->I.o_class = objtype + no; - asl->I.c_class = asl->I.o_class + no; - } - if (nvref) { - vrefnext = vrefx + maxfwd1; - nvref -= maxfwd1; - } - if (n_cc && !cvar) - cvar = (int*)M1alloc(nc*sizeof(int)); - if (cvar) - memset(cvar, 0, nc*sizeof(int)); - last_d = 0; - ka = 0; - kaz = 0; - kseen = 0; - nz = 0; - for(;;) { - ER.can_end = 1; - i = edag_peek(R); - if (i == EOF) { - fclose(nl); - do_ewalk(S); - if (imap) - del_mblk(kimap, imap); - /* Make amax long enough for nlc to handle */ - /* var_e[i].a for common variables i. */ - if (ncom0) { - i = comb + como; - if (i < combc) - i = combc; - if ((i += noa + 1) > amax) - amax = i; - } - adjoints = (real *)M1zapalloc(amax*Sizeof(real)); - adjoints_nv1 = &adjoints[nv1]; - nderps += nderp; - adjust(S, flags); - nzjac = nz; - if (!Lastx) - Lastx = (real *)M1alloc(nvx*sizeof(real)); -#ifdef PSHVREAD - a->i.ncxval = (int*)M1zapalloc(nco*sizeof(int)); - a->i.noxval = a->i.ncxval + nc; - a->p.Conival = conpival_ASL; - a->p.Conival_nomap = conpival_nomap_ASL; - a->p.Congrd = conpgrd_ASL; - a->p.Congrd_nomap = conpgrd_nomap_ASL; - a->p.Objval = a->p.Objval_nomap = objpval_ASL; - a->p.Objgrd = a->p.Objgrd_nomap = objpgrd_ASL; - a->p.Conval = conpval_ASL; - a->p.Jacval = jacpval_ASL; - a->p.Lconval= lconpval_ASL; - a->p.Hvcomp = a->p.Hvcomp_nomap = hvpcomp_ASL; - a->p.Hvcompd = hvpcompd_ASL; - a->p.Hvcompde = hvpcompde_ASL; - a->p.Hvcomps = hvpcomps_ASL; - a->p.Hvcompse = hvpcompse_ASL; - a->p.Hvinit = a->p.Hvinit_nomap = hvpinit_ASL; - a->p.Hvinite = a->p.Hvinite_nomap = hvpinite_ASL; - a->p.Hesset = hes_setup; - a->p.Xknown = xp2known_ASL; - a->p.Duthes = a->p.Duthes_nomap = duthes_ASL; - a->p.Duthese = a->p.Duthese_nomap = duthese_ASL; - a->p.Fulhes = a->p.Fulhes_nomap = fullhes_ASL; - a->p.Fulhese = a->p.Fulhese_nomap = fullhese_ASL; - a->p.Sphes = a->p.Sphes_nomap = sphes_ASL; - a->p.Sphese = a->p.Sphese_nomap = sphese_ASL; - a->p.Sphset = a->p.Sphset_nomap = sphes_setup_ASL; -#else - a->p.Xknown = xp1known_ASL; -#endif - return prob_adj_ASL(a); - } - ER.can_end = 0; - k = -1; - switch(i) { - case 'C': - Xscanf(R, "%d", &k); - if (k < 0 || k >= nc0) - badline(R); - co_read(R, con_de, k); - break; - case 'F': - if (Xscanf(R, "%d %d %d %127s", - &i, &j, &k, fname) != 4 - || i < 0 || i >= nfunc) - badline(R); - if ((fi = func_lookup(a, fname,0))) { - if (fi->nargs != k && fi->nargs >= 0 - && (k >= 0 || fi->nargs < -(k+1))) - scream(R, ASL_readerr_argerr, - "function %s: disagreement of nargs: %d and %d\n", - fname,fi->nargs, k); - } - else { - fi = (func_info *)mem(sizeof(func_info)); - fi->ftype = j; - fi->nargs = k; - fi->funcp = 0; - fi->name = (Const char *)strcpy((char*) - mem(memadj(strlen(fname)+1)), - fname); - } - if (!fi->funcp && !(fi->funcp = dynlink(fname))) - scream(R, ASL_readerr_unavail, - "function %s not available\n", - fname); - funcs[i] = fi; - break; - case 'L': - Xscanf(R, "%d", &k); - if (k < 0 || k >= nlogc) - badline(R); - co_read(R, lcon_de, k); - ewalk(S, lcon_de[k].e, 0); - break; - case 'V': - if (Xscanf(R, "%d %d %d", &k, &nlin, &j) != 3) - badline(R); - if (k >= SS.nvar0) - k += SS.nvinc; - if (k < nvr || k >= nv) - badline(R); - cexp_read(R, k, nlin); - break; - case 'G': - if (Xscanf(R, "%d %d", &j, &k) != 2 - || j < 0 || j >= no || k <= 0 || k > nvo) - badline(R); - ogp = Ograd + j; - while(k--) { - if (Xscanf(R, "%d %lf", &i, &t) != 2) - badline(R); - if (allG || t) { - *ogp = og = (ograd *) - mem(sizeof(ograd)); - ogp = &og->next; - og->varno = i; - og->coef = t; - } - } - *ogp = 0; - break; - case 'J': - if (Xscanf(R, "%d %d", &j, &k) != 2 - || j < 0 || j >= nc || k <= 0 || k > nvc) - badline(R); - nz += k; - if (A_vals) { - j += Fortran; - if (ka) { - while(k--) { - if (Xscanf(R, "%d %lf", - &i, &t) != 2) - badline(R); - i1 = ka[i]++; - A_vals[i1] = t; - A_rownos[i1] = j; - } - } - else { - while(k--) { - if (Xscanf(R, "%d %lf", - &i, &t) != 2) - badline(R); - i1 = kaz[i]++; - A_vals[i1] = t; - A_rownos[i1] = j; - } - } - break; - } - cgp = Cgrad + j; - j = 0; - while(k--) { - if (kseen) { - if (Xscanf(R, "%d %lf", &i, &t) != 2) - badline(R); - } - else - if (Xscanf(R, "%d %d %lf", &i, &j, &t) != 3) - badline(R); - *cgp = cg = (cgrad *) - mem(sizeof(cgrad)); - cgp = &cg->next; - cg->varno = i; - cg->goff = j; - cg->coef = t; - } - *cgp = 0; - break; - case 'O': - if (Xscanf(R, "%d %d", &k, &j) != 2 - || k < 0 || k >= no) - badline(R); - objtype[k] = j; - co_read(R, obj_de, k); - break; - case 'S': - Suf_read_ASL(R, readall); - break; - case 'r': - br_read(R, asl->i.n_con0, LUrhs, Urhsx, cvar, nvr); - break; - case 'b': - br_read(R, asl->i.n_var0, LUv, Uvx, 0, 0); - break; - case 'K': - case 'k': - k_seen = ++kseen; - if (ka_read_ASL(a, R, i, &ka, &kaz)) - badline(R); - break; - case 'x': - if (!Xscanf(R,"%d",&k) - || k < 0 || k > nvr) - badline(R); - if (!X0 && want_xpi0 & 1) { - x = nvx*sizeof(real); - if (want_xpi0 & 4) - x += nvx; - X0 = (real *)M1zapalloc(x); - if (want_xpi0 & 4) - havex0 = (char*)(X0 + nvx); - } - while(k--) { - if (Xscanf(R, "%d %lf", &j, &t) != 2 - || j < 0 || j >= nvr) - badline(R); - if (X0) { - X0[j] = t; - if (havex0) - havex0[j] = 1; - } - } - break; - case 'd': - if (!Xscanf(R,"%d",&k) - || k < 0 || k > nc0) - badline(R); - if (!pi0 && want_xpi0 & 2) { - x = nc*sizeof(real); - if (want_xpi0 & 4) - x += nc; - pi0 = (real *)M1zapalloc(x); - if (want_xpi0 & 4) - havepi0 = (char*)(pi0 + nc); - } - while(k--) { - if (Xscanf(R, "%d %lf", &j, &t) != 2 - || j < 0 || j >= nc0) - badline(R); - if (pi0) { - pi0[j] = t; - if (havepi0) - havepi0[j] = 1; - } - } - break; - default: - badline(R); - } - } - } -#undef asl - -#ifdef PSHVREAD - - static int -heswork(expr *e) -{ - argpair *da, *dae; - int i, j, n; - expr *e1, **ep; - de *d; - - n = 0; - if (e && e->op != OPVARVAL && e->op != OPNUM) do { - switch(e->a) { - - case Hv_timesR: - case Hv_timesL: - case Hv_negate: - case Hv_plterm: - n += 4; - break; - - case Hv_binaryR: - case Hv_unary: - n += 6; - break; - - case Hv_timesLR: - n += 10; - break; - - case Hv_divLR: - n += 13; - break; - - case Hv_binaryLR: - n += 14; - break; - - case Hv_vararg: - d = ((expr_va*)e)->L.d; - i = heswork(d->e); - while((e1 = (++d)->e)) { - j = heswork(e1); - if (i < j) - i = j; - } - n = i + 2; - break; - - case Hv_if: - i = heswork(((expr_if *)e)->Tf); - j = heswork(((expr_if *)e)->Ff); - if (i < j) - i = j; - n += i + 2; - break; - - - case Hv_sumlist: - for(ep = e->R.ep; (e1 = *ep); ep++) - n++; - break; - - case Hv_func: - da = ((expr_f *)e)->da; - dae = ((expr_f *)e)->dae; - i = dae - da; - n += i*(i+4); - break; - - case Hv_plusR: - case Hv_plusL: - case Hv_minusR: - n += 2; - break; - - case Hv_plusLR: - case Hv_minusLR: - n += 3; - break; - - default:/*DEBUG*/ - fprintf(Stderr, "bad e->a = %d in heswork\n", e->a); - exit(1); - } - } while((e = e->fwd)); - return n; - } - - static cexp * -hesfunnel(Static *S, int *ip, int nrefs, ograd **ogp, linarg ***lapp, linarg ***lapep) -{ - cexp *c; - expr *e; - int hw, i, k, m, n; - linarg *la, **lap, **lape; - ograd *og; - dv_info *dvi; - /*list *L;*//*20170615*/ - range *r; - split_ce *cs; /*20170615*/ - ASLTYPE *asl = S->asl; - - i = *ip; - c = cexps + i; -#if 1 /*20170615*/ - if (c->cref) - return 0; -#endif - n = 0; - if (i >= Ncom) { - cs = &asl->P.Split_ce[i-Ncom]; -#if 1 /*20170615*/ - if (cs->ce) - return 0; -#endif - r = cs->r; - lap = r->lap; - lape = lap + (n = r->n); - } - else if (!(lap = (dvi = asl->P.dv + i)->nl)) { - asl->P.linmultr++; - og = dvi->ll; - if (og->varno < 0) - og = og->next; - for(*ogp = og; og; og = og->next) - n++; - if (n > 1 && asl->p.hffactor > 0) { - asl->P.linhesfun++; - *ip = n; - return c; - } - return 0; - } - else { - lape = lap; - while(*lape) - lape++; - n = lape - lap; - } - if (!(e = c->ef)) - return 0; - *lapp = lap; - *lapep = lape; - *ogp = 0; - asl->P.nlmultr++; - while(lap < lape) { - la = *lap++; - for(og = la->nz; og; og = og->next) - if (!zc[og->varno]++) - zci[nzc++] = og->varno; - } - m = nzc; - while(nzc > 0) - zc[zci[--nzc]] = 0; -#if 0 /*20170615*/ - for(L = c->cref; L; L = L->next) - n++; -#endif - if (m > n) - m = n; - if ((k = m*nrefs - n) <= 0) - return 0; - hw = heswork(e); - if (k*hw <= nrefs*n*(n+3)*asl->p.hffactor) - return 0; - *ip = n; - asl->P.nlhesfun++; - return c; - } - - ASL_pfgh * -pscheck_ASL(ASL *a, const char *who1) -{ - ASL_CHECK(a, ASL_read_pfgh, who1); - return (ASL_pfgh*)a; - } - - static void -hes_setup(ASL *a, int flags, int obj, int nobj, int con, int ncon) -{ - range *r, *r0; - int *c, *ce, *cei, i, k, n, n0, nmax, nvx, *z; - hes_fun *hf, *hfth; - cexp *C; - linarg *la, **lap, **lape; - ograd *og; - list *L; - expr_v **vp; - psb_elem *b; - ASL_pfgh *asl; - Static SS, *S; - - asl = pscheck_ASL(a, "hes_setup"); - - S = S_init(&SS, asl); - Ncom = asl->P.ncom; - nv0x = nvx = asl->P.nv0_; - max_var = nvx + Ncom; - max_var1 = max_var + asl->P.ndvspout; - nzclim = max_var1 >> 3; - varp = asl->P.vp; - zl = asl->P.zlsave; - - zc_upgrade(S); /* make zc and zci available */ - n = nmax = 0; - r0 = (range*)&asl->P.rlist; - for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) { - if (r->n <= 0) { - r->cei = 0; - continue; - } - if (nmax < r->n) - nmax = r->n; - n0 = 0; - cei = 0; - for(b = r->refs; b; b = b->next) { - i = b->conno; - if (i < 0) { - i = -i - 2; - if (i < obj || i-obj >= nobj) - continue; - } - else { - if (i < con || i-con >= ncon) - continue; - } - if ((c = b->ce)) { - if (!n0) { - cei = c; - n0 = *c; - } - ce = c + *c; - while(c < ce) - if (!zc[i = *++c]++) - zci[n++] = i; - } - } - if (n0 < n) - cei = 0; - if ((r->cei = cei)) { - while(n > 0) - zc[zci[--n]] = 0; - continue; - } - if (!n) - continue; - r->cei = cei = (int*)mem((n+1)*sizeof(int)); - *cei++ = n; - if (n > 1) - qsortv(zci, n, sizeof(int), hscompar, S); - cei += n; - do zc[*--cei = zci[--n]] = 0; - while(n > 0); - } - z = zc + nvx; - k = -1; - for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) { - if (!(cei = r->cei)) - continue; - c = cei + 1; - ce = c + *cei; - do { - if (z[i = *c++]++ && k < i) - k = i; - } while(c < ce); - } - hfth = 0; - asl->P.linmultr = asl->P.linhesfun = 0; - asl->P.nlmultr = asl->P.nlhesfun = 0; - lap = lape = 0; /* silence false "used uninitialized" warning */ - while (k >= 0) - if (z[i = k--] > 1 - && (C = hesfunnel(S,&i,z[i],&og,&lap,&lape))) { - hf = (hes_fun*)mem(sizeof(hes_fun)); - hf->hfthread = hfth; - C->hfun = hfth = hf; - hf->c = C; - hf->n = i; - if ((hf->og = og)) { - hf->vp = 0; - hf->grdhes = 0; - } - else { - hf->vp = vp = (expr_v **) - mem(i*sizeof(expr_v*)); - while(lap < lape) { - la = *lap++; - *vp++ = la->v; - } - for(L = C->cref; L; L = L->next) - *vp++ = varp[L->item.i - nvx]; - i += i*i; - hf->grdhes = (real*)mem(i*sizeof(real)); - } - } - if ((asl->I.hesthread = hfth)) { - asl->I.x0kind_init = ASL_need_funnel; - if (x0kind != ASL_first_x) - x0kind |= ASL_need_funnel; - } - else - asl->I.x0kind_init = 0; - del_mblk(kzc, zci); - asl->P.hes_setup_called = 1 | (flags << 2); - asl->P.khesoprod = 5; - asl->P.hop_free = 0; - asl->P.nmax = nmax; - if (!(flags & 1)) - return; - i = htcl(n = nvx*sizeof(range*)); - k = htcl(3*n); - if (k == i + 1) { - asl->P.rtodo = (range**)new_mblkzap(asl, k); - asl->P.utodo = (uHeswork **)(asl->P.rtodo + nvx); - asl->P.otodo = (Hesoprod**)(asl->P.utodo + nvx); - } - else { - asl->P.rtodo = (range**)new_mblkzap(asl, i); - asl->P.utodo = (uHeswork**)new_mblkzap(asl, i); - asl->P.otodo = (Hesoprod**)new_mblkzap(asl, i); - } - k = htcl(nmax*(sizeof(int)+sizeof(real))); - asl->P.dOscratch = (real *)new_mblkzap(asl, k); - asl->P.iOscratch = (int *)(asl->P.dOscratch + nmax); - } - -#endif /* PSHVREAD */ - -#ifdef __cplusplus -} -#endif diff --git a/src/solvers2/pfghread.c b/src/solvers2/pfghread.c index f905922..a7ee2dd 100644 --- a/src/solvers2/pfghread.c +++ b/src/solvers2/pfghread.c @@ -6647,7 +6647,7 @@ hes_setup(Static *S) ASL_pfgh *asl; cexp *C; hes_fun *hf, *hfth; - int *c, *ce, *cei, *dvsp0, h, h0, i, i0, k, k1, k2, n, n0, ncom; + int *c, *ce, *cei, *dvsp0, h, h0, i, i0, j, j1, k, k1, k2, n, n0, ncom; int *ndvsp, nmax, ns, nvx, *vp, *vr, *vre, *z; linarg *la, **lap, **lape; ograd *og; @@ -6764,8 +6764,9 @@ hes_setup(Static *S) asl->P.khesoprod = 5; asl->P.nmax = nmax; asl->P.rtodo = h0 = h; - /* 20240616 changing asl->P.nran to asl->i.defvar0 in the next line */ - h += n = (asl->i.defvar0*sizeof(range*) + sizeof(real) - 1) / sizeof(real); + if ((j = asl->P.nran*sizeof(range*)) < (j1 = asl->i.defvar0*sizeof(real))) + j = j1; + h += n = (j + sizeof(real) - 1) / sizeof(real); asl->P.utodo = h; h += n; asl->P.otodo = h; diff --git a/src/solvers2/pshvprod.c b/src/solvers2/pshvprod.c index 7977ad6..e744412 100644 --- a/src/solvers2/pshvprod.c +++ b/src/solvers2/pshvprod.c @@ -28,22 +28,12 @@ extern "C" { #define alignarg(x) #endif -#ifdef ALLOW_OPENMP -#include -#define NO_PTHREADS -#endif #ifdef NO_PTHREADS #define PTHREADS(x) #define pthread_mutex_lock(t) /*nothing*/ #define pthread_mutex_unlock(t) /*nothing*/ -#ifndef ALLOW_OPENMP -#undef MULTIPLE_THREADS -#endif #else #define PTHREADS(x) x -#undef MULTIPLE_THREADS -#define MULTIPLE_THREADS -#include #endif #ifdef MULTIPLE_THREADS @@ -1639,7 +1629,6 @@ hvtodo(ParHvInfo *tp, Thparshv1 *tp1) { Ihinfo *ihi; range *r; - int k; if (!(ihi = tp->ihi)) return 0; diff --git a/src/solvers2/sphes.c b/src/solvers2/sphes.c index f7bf39a..1517c19 100644 --- a/src/solvers2/sphes.c +++ b/src/solvers2/sphes.c @@ -24,21 +24,12 @@ of or in connection with the use or performance of this software. #else #define alignarg(x) /*nothing*/ #endif -#ifdef ALLOW_OPENMP -/* asl.h has #defined MULTIPLE_THREADS */ -#include -#define NO_PTHREADS -#elif !defined(NO_PTHREADS) -#include -#ifndef MULTIPLE_THREADS -#define MULTIPLE_THREADS -#endif -#define PTHREADS(x) x -#endif #ifdef NO_PTHREADS #define pthread_mutex_lock(t) /*nothing*/ #define pthread_mutex_unlock(t) /*nothing*/ #define PTHREADS(x) /*nothing*/ +#else +#define PTHREADS(x) x #endif #undef exit extern void wk_init_ASL(real*, int*, real); @@ -1674,10 +1665,7 @@ sphes_setup_ew_ASL(EvalWorkspace *ew, SputInfo **pspi, int nobj, int ow, int y, tp.iow = ow; tp.iy = y; tp.utodo = utodo; - hs0 = hs = (real*)new_mblk(htcl(m*( -#ifndef ALLOW_OPENMP - sizeof(pthread_t*) + -#endif + hs0 = hs = (real*)new_mblk(htcl(m*(PTHREADS(sizeof(pthread_t*) +) sizeof(Thpars1)) + asl->P.thlen*sizeof(real) + nrng*(sizeof(RangeInfo)))); @@ -2259,10 +2247,7 @@ sphes_ew_ASL(EvalWorkspace *ew, SputInfo **pspi, real *H, int nobj, real *ow, re if (!(hs0 = hs = asl->P.hesparwk)) { alloc_anew: asl->P.hesparwk = - hs = (real*)M1alloc(m*( -#ifndef ALLOW_OPENMP - sizeof(pthread_t*) + -#endif + hs = (real*)M1alloc(m*(PTHREADS(sizeof(pthread_t*) +) sizeof(Thpars1)) + asl->P.thlen*sizeof(real) + nrng*(sizeof(RangeInfo))); diff --git a/src/solvers2/xp2known.c b/src/solvers2/xp2known.c index 9f037de..42cb0a4 100644 --- a/src/solvers2/xp2known.c +++ b/src/solvers2/xp2known.c @@ -363,6 +363,7 @@ hvpinit_nc_ASL(EvalWorkspace *ew, int ndhmax, int nobj, real *ow, real *y) PTHREADS(pthread_t *T;) real *W; #ifndef ALLOW_OPENMP + Thparshv2 *tpi; int it; void *rc; #endif @@ -424,12 +425,11 @@ hvpinit_nc_ASL(EvalWorkspace *ew, int ndhmax, int nobj, real *ow, real *y) phvi->ihi = asl->P.ihi1; tp1 = phvi->tpv1; phvi->r = 0; - Thparshv2* tpi; #ifdef ALLOW_OPENMP #pragma omp parallel num_threads(m) { int it, jt; - + Thparshv2 *tpi; it = omp_get_thread_num(); tpi = tp1 + it; diff --git a/src/solvers2/xsum0.out b/src/solvers2/xsum0.out index e2df939..fad1667 100644 --- a/src/solvers2/xsum0.out +++ b/src/solvers2/xsum0.out @@ -5,10 +5,10 @@ arith.h0 e206209c 2309 arith.h1 12e5e42c 1871 arith.ibm f398be2b 1391 arithchk.c 63b0185 6532 -asl.h ff744e06 48918 +asl.h 1d93ebd2 48993 asl_pfg.h 6483336 1268 asl_pfgh.h 49b5a53 1288 -asldate.c 8721c18 29 +asldate.c f4afa32d 29 atof.c 1fabc7d3 1747 auxinfo.c 1a3c8690 1488 avltree.c 10aeaa58 11346 @@ -91,9 +91,9 @@ opcode.hd 12bc8de3 1348 opno.hd f8e46748 2123 opno2.h f82baf62 4720 opnos.hd 1be96111 120 -pfghread.c f3b93e8b 131968 +pfghread.c f898c96c 131961 printf.c 147de217 24293 -pshvprod.c f6c5f139 61858 +pshvprod.c 1d09aeff 61664 psinfo.h e05960fc 13888 punknown.c fa46cf10 1704 qpcheck.c 1712f297 1588 @@ -107,7 +107,7 @@ rnd_prod.s f4996bed 1471 sigcatch.c fe35ba8b 2374 sjac0dim.c ff686adb 1546 sos_add.c f59f41c0 22398 -sphes.c 354e74f 54173 +sphes.c 187471ef 53914 sprintf.c e6ce1d45 3050 sscanf.c e8e79b09 3302 stderr.c fc99b4f4 2096 @@ -121,4 +121,4 @@ wrtsol_.c 15a9faa4 2431 ws_desc.c fceb8f65 214 wsu_desc.c fbb1c30c 170 xectim.c e56314e8 3543 -xp2known.c 3d83d8f 11944 +xp2known.c 1ca179c 11961