=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v retrieving revision 1.70 retrieving revision 1.75 diff -u -p -r1.70 -r1.75 --- OpenXM_contrib2/asir2000/engine/nd.c 2003/09/15 10:51:45 1.70 +++ OpenXM_contrib2/asir2000/engine/nd.c 2003/09/19 10:09:42 1.75 @@ -1,6 +1,8 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.69 2003/09/15 09:49:44 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.74 2003/09/19 02:33:12 noro Exp $ */ #include "ca.h" +#include "parse.h" +#include "ox.h" #include "inline.h" #include @@ -114,6 +116,7 @@ typedef struct oIndArray int (*ndl_compare_function)(UINT *a1,UINT *a2); +static int ndv_alloc; static int nd_f4_nsp=0x7fffffff; static double nd_scale=2; static UINT **nd_bound; @@ -142,8 +145,9 @@ static int nm_adv; static int nmv_adv; static int nd_dcomp; +extern struct order_spec dp_current_spec; extern VL CO; -extern int Top,Reverse,dp_nelim,do_weyl; +extern int Top,Reverse,DP_Print,dp_nelim,do_weyl; extern int *current_weyl_weight_vector; /* fundamental macros */ @@ -224,7 +228,9 @@ if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT #define NMV_OPREV(m) (m = (NMV)(((char *)m)-oadv)) /* external functions */ +#if 1 void GC_gcollect(); +#endif NODE append_one(NODE,int); /* manipulation of coefficients */ @@ -331,7 +337,10 @@ ND nd_copy(ND p); ND nd_merge(ND p1,ND p2); ND nd_add(int mod,ND p1,ND p2); ND nd_add_q(ND p1,ND p2); +ND nd_add_sf(ND p1,ND p2); INLINE int nd_length(ND p); +NODE nd_f4_red(int m,ND_pairs sp0,UINT *s0vect,int col,NODE rp0); +NODE nd_f4_red_dist(int m,ND_pairs sp0,UINT *s0vect,int col,NODE rp0); /* NDV functions */ ND weyl_ndv_mul_nm(int mod,NM m0,NDV p); @@ -361,11 +370,11 @@ int nd_to_vect(int mod,UINT *s0,int n,ND d,UINT *r); void nd_free_private_storage() { - _nd_free_list = 0; _nm_free_list = 0; _ndp_free_list = 0; - bzero(nd_red,sizeof(REDTAB_LEN*sizeof(RHist))); +#if 0 GC_gcollect(); +#endif } void _NM_alloc() @@ -1141,6 +1150,7 @@ ND nd_add(int mod,ND p1,ND p2) if ( !p1 ) return p2; else if ( !p2 ) return p1; + else if ( mod == -1 ) return nd_add_sf(p1,p2); else if ( !mod ) return nd_add_q(p1,p2); else { can = 0; @@ -1228,6 +1238,53 @@ ND nd_add_q(ND p1,ND p2) } } +ND nd_add_sf(ND p1,ND p2) +{ + int n,c,can; + ND r; + NM m1,m2,mr0,mr,s; + int t; + + if ( !p1 ) return p2; + else if ( !p2 ) return p1; + else { + can = 0; + for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) { + c = DL_COMPARE(DL(m1),DL(m2)); + switch ( c ) { + case 0: + t = _addsf(CM(m1),CM(m2)); + s = m1; m1 = NEXT(m1); + if ( t ) { + can++; NEXTNM2(mr0,mr,s); CM(mr) = (t); + } else { + can += 2; FREENM(s); + } + s = m2; m2 = NEXT(m2); FREENM(s); + break; + case 1: + s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s); + break; + case -1: + s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s); + break; + } + } + if ( !mr0 ) + if ( m1 ) mr0 = m1; + else if ( m2 ) mr0 = m2; + else return 0; + else if ( m1 ) NEXT(mr) = m1; + else if ( m2 ) NEXT(mr) = m2; + else NEXT(mr) = 0; + BDY(p1) = mr0; + SG(p1) = MAX(SG(p1),SG(p2)); + LEN(p1) = LEN(p1)+LEN(p2)-can; + FREEND(p2); + return p1; + } +} + /* ret=1 : success, ret=0 : overflow */ int nd_nf(int mod,ND g,NDV *ps,int full,NDC dn,ND *rp) { @@ -1266,7 +1323,9 @@ int nd_nf(int mod,ND g,NDV *ps,int full,NDC dn,ND *rp) return 0; } p = ps[index]; - if ( mod ) { + if ( mod == -1 ) + CM(mul) = _mulsf(_invsf(HCM(p)),_chsgnsf(HCM(g))); + else if ( mod ) { c1 = invm(HCM(p),mod); c2 = mod-HCM(g); DMAR(c1,c2,0,mod,c); CM(mul) = c; } else { @@ -1353,7 +1412,9 @@ int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp return 0; } p = ps[index]; - if ( mod ) { + if ( mod == -1 ) + CM(mul) = _mulsf(_invsf(HCM(p)),_chsgnsf(HCM(g))); + else if ( mod ) { c1 = invm(HCM(p),mod); c2 = mod-HCM(g); DMAR(c1,c2,0,mod,c); CM(mul) = c; } else { @@ -1430,9 +1491,9 @@ again: nd_reconstruct(0,0,0); goto again; } else if ( nf ) return 0; - printf("."); fflush(stdout); + if ( DP_Print ) { printf("."); fflush(stdout); } } - printf("\n"); + if ( DP_Print ) { printf("\n"); } /* gbcheck : cand is a GB of Id(cand) ? */ if ( !nd_gb(0,1) ) return 0; /* XXX */ @@ -1594,8 +1655,12 @@ int head_pbucket(int mod,PGeoBucket g) dj = HDL(gj); sum = HCM(gj); } else if ( c == 0 ) { - sum = sum+HCM(gi)-mod; - if ( sum < 0 ) sum += mod; + if ( mod == -1 ) + sum = _addsf(sum,HCM(gi)); + else { + sum = sum+HCM(gi)-mod; + if ( sum < 0 ) sum += mod; + } g->body[i] = nd_remove_head(gi); } } @@ -1686,7 +1751,7 @@ again: l = nd_minp(d,&d); if ( SG(l) != sugar ) { sugar = SG(l); - fprintf(asir_out,"%d",sugar); + if ( DP_Print ) fprintf(asir_out,"%d",sugar); } stat = nd_sp(m,0,l,&h); if ( !stat ) { @@ -1705,7 +1770,7 @@ again: goto again; } else if ( nf ) { if ( checkonly ) return 0; - printf("+"); fflush(stdout); + if ( DP_Print ) { printf("+"); fflush(stdout); } nd_removecont(m,nf); nfv = ndtondv(m,nf); nd_free(nf); nh = ndv_newps(nfv,0); @@ -1713,7 +1778,7 @@ again: g = update_base(g,nh); FREENDP(l); } else { - printf("."); fflush(stdout); + if ( DP_Print ) { printf("."); fflush(stdout); } FREENDP(l); } } @@ -1742,7 +1807,7 @@ again: l = nd_minp(d,&d); if ( SG(l) != sugar ) { sugar = SG(l); - fprintf(asir_out,"%d",sugar); + if ( DP_Print ) fprintf(asir_out,"%d",sugar); } stat = nd_sp(m,0,l,&h); if ( !stat ) { @@ -1769,17 +1834,17 @@ again: /* m|HC(nfq) => failure */ if ( !rem(NM(HCQ(nfq)),m) ) return 0; - printf("+"); fflush(stdout); + if ( DP_Print ) { printf("+"); fflush(stdout); } nd_removecont(m,nf); nfv = ndtondv(m,nf); nd_free(nf); nd_removecont(0,nfq); nfqv = ndtondv(0,nfq); nd_free(nfq); nh = ndv_newps(nfv,nfqv); d = update_pairs(d,g,nh); g = update_base(g,nh); } else { - printf("*"); fflush(stdout); + if ( DP_Print ) { printf("*"); fflush(stdout); } } } else { - printf("."); fflush(stdout); + if ( DP_Print ) { printf("."); fflush(stdout); } } FREENDP(l); } @@ -1823,7 +1888,7 @@ NODE ndv_reduceall(int m,NODE f) if ( !stat ) nd_reconstruct(m,0,0); else { - printf("."); fflush(stdout); + if ( DP_Print ) { printf("."); fflush(stdout); } if ( !m ) { mulq(HCQ(head),dn.z,&q); HCQ(head) = q; } nf = nd_add(m,head,nf); ndv_free(nd_ps[i]); @@ -1833,7 +1898,7 @@ NODE ndv_reduceall(int m,NODE f) i++; } } - printf("\n"); + if ( DP_Print ) { printf("\n"); } for ( a0 = 0, i = 0; i < n; i++ ) { NEXTNODE(a0,a); BDY(a) = (pointer)nd_ps[i]; @@ -2156,8 +2221,7 @@ void ndv_setup(int mod,int trace,NODE f) if ( !nd_red ) nd_red = (RHist *)MALLOC(REDTAB_LEN*sizeof(RHist)); - bzero(nd_red,REDTAB_LEN*sizeof(RHist)); - nd_free_private_storage(); + for ( i = 0; i < REDTAB_LEN; i++ ) nd_red[i] = 0; for ( i = 0; i < nd_psn; i++ ) { if ( trace ) { a = nd_ps_trace[i] = ndv_dup(0,w[i]); @@ -2182,8 +2246,9 @@ void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe int e,max,nvar; NDV b; + ndv_alloc = 0; get_vars((Obj)f,&fv); pltovl(v,&vv); - nvar = length(vv); + for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); nd_init_ord(ord); for ( t = BDY(f), max = 0; t; t = NEXT(t) ) for ( tv = vv; tv; tv = NEXT(tv) ) { @@ -2207,6 +2272,7 @@ void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe } if ( r0 ) NEXT(r) = 0; MKLIST(*rp,r0); + fprintf(asir_out,"ndv_alloc=%d\n",ndv_alloc); } void nd_gr_trace(LIST f,LIST v,int trace,int homo,struct order_spec *ord,LIST *rp) @@ -2222,7 +2288,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru int obpe,oadv,wmax,i,len,cbpe; get_vars((Obj)f,&fv); pltovl(v,&vv); - nvar = length(vv); + for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); nocheck = 0; mindex = 0; @@ -2302,7 +2368,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru } /* dp->p */ nd_bpe = cbpe; - nd_setup_parameters(0,0); + nd_setup_parameters(nd_nvar,0); for ( r = cand; r; r = NEXT(r) ) BDY(r) = (pointer)ndvtop(0,CO,vv,BDY(r)); MKLIST(*rp,cand); } @@ -2383,7 +2449,8 @@ void nd_print(ND p) printf("0\n"); else { for ( m = BDY(p); m; m = NEXT(m) ) { - printf("+%d*",CM(m)); + if ( CM(m) & 0x80000000 ) printf("+@_%d*",IFTOF(CM(m))); + else printf("+%d*",CM(m)); ndl_print(DL(m)); } printf("\n"); @@ -2399,7 +2466,7 @@ void nd_print_q(ND p) else { for ( m = BDY(p); m; m = NEXT(m) ) { printf("+"); - printexpr(CO,CQ(m)); + printexpr(CO,(Obj)CQ(m)); printf("*"); ndl_print(DL(m)); } @@ -2424,7 +2491,8 @@ void nd_removecont(int mod,ND p) struct oVECT v; N q,r; - if ( mod ) nd_mul_c(mod,p,invm(HCM(p),mod)); + if ( mod == -1 ) nd_mul_c(mod,p,_invsf(HCM(p))); + else if ( mod ) nd_mul_c(mod,p,invm(HCM(p),mod)); else { for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ ); w = (Q *)ALLOCA(n*sizeof(Q)); @@ -2470,7 +2538,9 @@ void ndv_removecont(int mod,NDV p) Q dvr,t; NMV m; - if ( mod ) + if ( mod == -1 ) + ndv_mul_c(mod,p,_invsf(HCM(p))); + else if ( mod ) ndv_mul_c(mod,p,invm(HCM(p),mod)); else { len = p->len; @@ -2571,11 +2641,13 @@ void nd_mul_c(int mod,ND p,int mul) int c,c1; if ( !p ) return; - for ( m = BDY(p); m; m = NEXT(m) ) { - c1 = CM(m); - DMAR(c1,mul,0,mod,c); - CM(m) = c; - } + if ( mod == -1 ) + for ( m = BDY(p); m; m = NEXT(m) ) + CM(m) = _mulsf(CM(m),mul); + else + for ( m = BDY(p); m; m = NEXT(m) ) { + c1 = CM(m); DMAR(c1,mul,0,mod,c); CM(m) = c; + } } void nd_mul_c_q(ND p,Q mul) @@ -2678,23 +2750,23 @@ int nd_get_exporigin(struct order_spec *ord) void nd_setup_parameters(int nvar,int max) { int i,j,n,elen,ord_o,ord_l,l,s; struct order_pair *op; + int bpe; - /* if max == 0, don't touch nd_bpe */ - if ( max > 0 ) { - if ( max < 2 ) nd_bpe = 1; - if ( max < 4 ) nd_bpe = 2; - else if ( max < 8 ) nd_bpe = 3; - else if ( max < 16 ) nd_bpe = 4; - else if ( max < 32 ) nd_bpe = 5; - else if ( max < 64 ) nd_bpe = 6; - else if ( max < 256 ) nd_bpe = 8; - else if ( max < 1024 ) nd_bpe = 10; - else if ( max < 65536 ) nd_bpe = 16; - else nd_bpe = 32; - } - /* nvar == 0, don't touch nd_nvar */ - if ( nvar > 0 ) nd_nvar = nvar; - + if ( !max ) bpe = nd_bpe; + else if ( max < 2 ) bpe = 1; + else if ( max < 4 ) bpe = 2; + else if ( max < 8 ) bpe = 3; + else if ( max < 16 ) bpe = 4; + else if ( max < 32 ) bpe = 5; + else if ( max < 64 ) bpe = 6; + else if ( max < 256 ) bpe = 8; + else if ( max < 1024 ) bpe = 10; + else if ( max < 65536 ) bpe = 16; + else bpe = 32; + if ( bpe != nd_bpe || nvar != nd_nvar ) + nd_free_private_storage(); + nd_bpe = bpe; + nd_nvar = nvar; nd_epw = (sizeof(UINT)*8)/nd_bpe; elen = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0); @@ -2721,11 +2793,12 @@ void nd_setup_parameters(int nvar,int max) { ND_pairs nd_reconstruct(int mod,int trace,ND_pairs d) { int i,obpe,oadv,h; - NM prev_nm_free_list; + static NM prev_nm_free_list; + static ND_pairs prev_ndp_free_list; RHist mr0,mr; RHist r; RHist *old_red; - ND_pairs s0,s,t,prev_ndp_free_list; + ND_pairs s0,s,t; EPOS oepos; obpe = nd_bpe; @@ -2742,7 +2815,7 @@ ND_pairs nd_reconstruct(int mod,int trace,ND_pairs d) else if ( obpe < 32 ) nd_bpe = 32; else error("nd_reconstruct : exponent too large"); - nd_setup_parameters(0,0); + nd_setup_parameters(nd_nvar,0); prev_nm_free_list = _nm_free_list; prev_ndp_free_list = _ndp_free_list; _nm_free_list = 0; @@ -2785,7 +2858,9 @@ ND_pairs nd_reconstruct(int mod,int trace,ND_pairs d) if ( s0 ) NEXT(s) = 0; prev_nm_free_list = 0; prev_ndp_free_list = 0; +#if 0 GC_gcollect(); +#endif return s0; } @@ -2860,7 +2935,8 @@ int nd_sp(int mod,int trace,ND_pairs p,ND *rp) if ( ndl_check_bound2(p->i1,DL(m)) ) return 0; t1 = ndv_mul_nm(mod,m,p1); - if ( mod ) CM(m) = mod-HCM(p1); + if ( mod == -1 ) CM(m) = _chsgnsf(HCM(p1)); + else if ( mod ) CM(m) = mod-HCM(p1); else chsgnq(HCQ(p1),&CQ(m)); ndl_sub(lcm,HDL(p2),DL(m)); if ( ndl_check_bound2(p->i2,DL(m)) ) { @@ -2880,9 +2956,13 @@ void ndv_mul_c(int mod,NDV p,int mul) if ( !p ) return; len = LEN(p); - for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) { - c1 = CM(m); DMAR(c1,mul,0,mod,c); CM(m) = c; - } + if ( mod == -1 ) + for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) + CM(m) = _mulsf(CM(m),mul); + else + for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) { + c1 = CM(m); DMAR(c1,mul,0,mod,c); CM(m) = c; + } } void ndv_mul_c_q(NDV p,Q mul) @@ -3072,17 +3152,27 @@ ND ndv_mul_nm(int mod,NM m0,NDV p) if ( !p ) return 0; else if ( do_weyl ) - return weyl_ndv_mul_nm(mod,m0,p); + if ( mod == -1 ) + error("ndv_mul_nm : not implemented (weyl)"); + else + return weyl_ndv_mul_nm(mod,m0,p); else { n = NV(p); m = BDY(p); d = DL(m0); len = LEN(p); mr0 = 0; td = TD(d); - if ( mod ) { + if ( mod == -1 ) { c = CM(m0); for ( i = 0; i < len; i++, NMV_ADV(m) ) { NEXTNM(mr0,mr); + CM(mr) = _mulsf(CM(m),c); + ndl_add(DL(m),d,DL(mr)); + } + } else if ( mod ) { + c = CM(m0); + for ( i = 0; i < len; i++, NMV_ADV(m) ) { + NEXTNM(mr0,mr); c1 = CM(m); DMAR(c1,c,0,mod,c2); CM(mr) = c2; @@ -3187,21 +3277,32 @@ void ndv_mod(int mod,NDV p) NMV t,d; int r; int i,len,dlen; + Obj gfs; if ( !p ) return; len = LEN(p); dlen = 0; - for ( t = d = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) { - r = rem(NM(CQ(t)),mod); - if ( r ) { - if ( SGN(CQ(t)) < 0 ) - r = mod-r; + if ( mod == -1 ) + for ( t = d = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) { + simp_ff((Obj)CP(t),&gfs); + r = FTOIF(CONT((GFS)gfs)); CM(d) = r; ndl_copy(DL(t),DL(d)); NMV_ADV(d); dlen++; } - } + else + for ( t = d = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) { + r = rem(NM(CQ(t)),mod); + if ( r ) { + if ( SGN(CQ(t)) < 0 ) + r = mod-r; + CM(d) = r; + ndl_copy(DL(t),DL(d)); + NMV_ADV(d); + dlen++; + } + } LEN(p) = dlen; } @@ -3276,6 +3377,7 @@ P ndvtop(int mod,VL vl,VL dvl,NDV p) P c; UINT *d; P s,r,u,t,w; + GFS gfs; if ( !p ) return 0; else { @@ -3283,7 +3385,9 @@ P ndvtop(int mod,VL vl,VL dvl,NDV p) n = NV(p); m = (NMV)(((char *)BDY(p))+nmv_adv*(len-1)); for ( j = len-1, s = 0; j >= 0; j--, NMV_PREV(m) ) { - if ( mod ) { + if ( mod == -1 ) { + e = IFTOF(CM(m)); MKGFS(e,gfs); c = (P)gfs; + } else if ( mod ) { STOQ(CM(m),q); c = (P)q; } else c = CP(m); @@ -3307,7 +3411,11 @@ NDV ndtondv(int mod,ND p) if ( !p ) return 0; len = LEN(p); - m0 = m = (NMV)(mod?MALLOC_ATOMIC(len*nmv_adv):MALLOC(len*nmv_adv)); + if ( mod ) + m0 = m = (NMV)GC_malloc_atomic_ignore_off_page(len*nmv_adv); + else + m0 = m = MALLOC(len*nmv_adv); + ndv_alloc += nmv_adv*len; for ( t = BDY(p), i = 0; t; t = NEXT(t), i++, NMV_ADV(m) ) { ndl_copy(DL(t),DL(m)); CQ(m) = CQ(t); @@ -3347,7 +3455,8 @@ void ndv_print(NDV p) else { len = LEN(p); for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) { - printf("+%d*",CM(m)); + if ( CM(m) & 0x80000000 ) printf("+@_%d*",IFTOF(CM(m))); + else printf("+%d*",CM(m)); ndl_print(DL(m)); } printf("\n"); @@ -3364,7 +3473,7 @@ void ndv_print_q(NDV p) len = LEN(p); for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) { printf("+"); - printexpr(CO,CQ(m)); + printexpr(CO,(Obj)CQ(m)); printf("*"); ndl_print(DL(m)); } @@ -3519,7 +3628,7 @@ void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec union oNDC dn; pltovl(v,&vv); - nvar = length(vv); + for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); /* get the degree bound */ for ( t = BDY(g), max = 0; t; t = NEXT(t) ) @@ -3578,6 +3687,24 @@ int nd_to_vect(int mod,UINT *s0,int n,ND d,UINT *r) return i; } +int ndv_to_vect(int mod,UINT *s0,int n,NDV d,UINT *r) +{ + NMV m; + UINT *t,*s; + int i,j,len; + + for ( i = 0; i < n; i++ ) r[i] = 0; + m = BDY(d); + len = LEN(d); + for ( i = j = 0, s = s0; j < len; j++, NMV_ADV(m)) { + t = DL(m); + for ( ; !ndl_equal(t,s); s += nd_wpd, i++ ); + r[i] = CM(m); + } + for ( i = 0; !r[i]; i++ ); + return i; +} + int nm_ind_pair_to_vect(int mod,UINT *s0,int n,NM_ind_pair pair,UINT *r) { NM m; @@ -3650,7 +3777,7 @@ IndArray nm_ind_pair_to_vect_compress(int mod,UINT *s0 } -void ndv_reduce_vect(int m,UINT *svect,int col,IndArray *imat,NODE rp0) +void ndv_reduce_vect(int m,UINT *svect,int col,IndArray *imat,NM_ind_pair *rp0,int nred) { int i,j,k,len,pos,prev; UINT c,c1,c2,c3,up,lo,dmy; @@ -3662,11 +3789,11 @@ void ndv_reduce_vect(int m,UINT *svect,int col,IndArra NMV mr; NODE rp; - for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) { + for ( i = 0; i < nred; i++ ) { ivect = imat[i]; k = ivect->head; svect[k] %= m; if ( c = svect[k] ) { - c = m-c; redv = nd_ps[((NM_ind_pair)BDY(rp))->index]; + c = m-c; redv = nd_ps[rp0[i]->index]; len = LEN(redv); mr = BDY(redv); svect[k] = 0; prev = k; switch ( ivect->width ) { @@ -3707,6 +3834,52 @@ void ndv_reduce_vect(int m,UINT *svect,int col,IndArra if ( svect[i] >= (UINT)m ) svect[i] %= m; } +void ndv_reduce_vect_sf(int m,UINT *svect,int col,IndArray *imat,NM_ind_pair *rp0,int nred) +{ + int i,j,k,len,pos,prev; + UINT c,c1,c2,c3,up,lo,dmy; + IndArray ivect; + unsigned char *ivc; + unsigned short *ivs; + unsigned int *ivi; + NDV redv; + NMV mr; + NODE rp; + + for ( i = 0; i < nred; i++ ) { + ivect = imat[i]; + k = ivect->head; svect[k] %= m; + if ( c = svect[k] ) { + c = _chsgnsf(c); redv = nd_ps[rp0[i]->index]; + len = LEN(redv); mr = BDY(redv); + svect[k] = 0; prev = k; + switch ( ivect->width ) { + case 1: + ivc = ivect->index.c; + for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { + pos = prev+ivc[j]; prev = pos; + svect[pos] = _addsf(_mulsf(CM(mr),c),svect[pos]); + } + break; + case 2: + ivs = ivect->index.s; + for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { + pos = prev+ivs[j]; prev = pos; + svect[pos] = _addsf(_mulsf(CM(mr),c),svect[pos]); + } + break; + case 4: + ivi = ivect->index.i; + for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { + pos = prev+ivi[j]; prev = pos; + svect[pos] = _addsf(_mulsf(CM(mr),c),svect[pos]); + } + break; + } + } + } +} + NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhead,UINT *s0vect) { int j,k,len; @@ -3718,7 +3891,8 @@ NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea for ( j = 0, len = 0; j < spcol; j++ ) if ( vect[j] ) len++; if ( !len ) return 0; else { - mr0 = (NMV)MALLOC_ATOMIC(nmv_adv*len); + mr0 = (NMV)GC_malloc_atomic_ignore_off_page(nmv_adv*len); + ndv_alloc += nmv_adv*len; mr = mr0; p = s0vect; for ( j = k = 0; j < col; j++, p += nd_wpd ) @@ -3732,23 +3906,20 @@ NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea } } -int nd_sp_f4(int m,ND_pairs l,PGeoBucket bucket,NODE *s) +int nd_sp_f4(int m,ND_pairs l,PGeoBucket bucket) { ND_pairs t; NODE sp0,sp; int stat; ND spol; - sp0 = 0; for ( t = l; t; t = NEXT(t) ) { stat = nd_sp(m,0,t,&spol); if ( !stat ) return 0; if ( spol ) { - NEXTNODE(sp0,sp); BDY(sp) = (pointer)nd_dup(spol); add_pbucket_symbolic(bucket,spol); } } - *s = sp0; return 1; } @@ -3782,7 +3953,8 @@ int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vec } col++; } - NEXT(rp) = 0; NEXT(s) = 0; + if ( rp0 ) NEXT(rp) = 0; + NEXT(s) = 0; s0v = (UINT *)MALLOC_ATOMIC(col*nd_wpd*sizeof(UINT)); for ( i = 0, p = s0v, s = s0; i < col; i++, p += nd_wpd, s = NEXT(s) ) ndl_copy(DL(s),p); @@ -3791,7 +3963,6 @@ int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vec return col; } -#if 0 NODE nd_f4(int m) { int i,nh,stat,index; @@ -3800,9 +3971,10 @@ NODE nd_f4(int m) ND spol,red; NDV nf,redv; NM s0,s; - NODE sp0,sp,rp0,rp; + NODE rp0,sp0,srp0,nflist; int nsp,nred,col,rank,len,k,j,a; UINT c; + UINT **spmat; UINT *s0vect,*svect,*p,*v; int *colstat; IndArray *imat; @@ -3811,29 +3983,21 @@ NODE nd_f4(int m) int sugar; PGeoBucket bucket; struct oEGT eg0,eg1,eg_f4; - static UINT **spmat; - static UINT *spb; - static int spblen; if ( !m ) error("nd_f4 : not implemented"); - + ndv_alloc = 0; g = 0; d = 0; for ( i = 0; i < nd_psn; i++ ) { d = update_pairs(d,g,i); g = update_base(g,i); } - if ( !spmat ) { - spmat = (UINT **)MALLOC(nd_f4_nsp*sizeof(UINT *)); - spblen = 10000; - spb = (UINT *)MALLOC_ATOMIC(nd_f4_nsp*spblen*sizeof(UINT)); - } while ( d ) { get_eg(&eg0); l = nd_minsugarp(d,&d); sugar = SG(l); bucket = create_pbucket(); - stat = nd_sp_f4(m,l,bucket,&sp0); + stat = nd_sp_f4(m,l,bucket); if ( !stat ) { for ( t = l; NEXT(t); t = NEXT(t) ); NEXT(t) = d; d = l; @@ -3848,53 +4012,17 @@ NODE nd_f4(int m) d = nd_reconstruct(m,0,d); continue; } - - nsp = length(sp0); nred = length(rp0); spcol = col-nred; - imat = (IndArray *)MALLOC(nred*sizeof(IndArray)); - rhead = (int *)MALLOC_ATOMIC(col*sizeof(int)); - for ( i = 0; i < col; i++ ) rhead[i] = 0; - - /* construction of index arrays */ - for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) { - imat[i] = nm_ind_pair_to_vect_compress(m,s0vect,col,(NM_ind_pair)BDY(rp)); - rhead[imat[i]->head] = 1; - } - - /* elimination (1st step) */ - svect = (UINT *)MALLOC_ATOMIC(col*sizeof(UINT)); - if ( spcol > spblen ) { - spblen = spcol; - spb = REALLOC(spb,spblen*nd_f4_nsp*sizeof(UINT)); - } - - for ( a = sprow = 0, sp = sp0, p = spb; a < nsp; a++, sp = NEXT(sp) ) { - nd_to_vect(m,s0vect,col,BDY(sp),svect); - ndv_reduce_vect(m,svect,col,imat,rp0); - for ( i = 0; i < col; i++ ) if ( svect[i] ) break; - if ( i < col ) { - spmat[sprow] = p; - v = spmat[sprow]; - for ( j = k = 0; j < col; j++ ) - if ( !rhead[j] ) v[k++] = svect[j]; - sprow++; - p += k; - } - } - /* free index arrays */ - for ( i = 0; i < nred; i++ ) GC_free(imat[i]->index.c); - - /* elimination (2nd step) */ - colstat = (int *)ALLOCA(spcol*sizeof(int)); - rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat); - get_eg(&eg1); init_eg(&eg_f4); add_eg(&eg_f4,&eg0,&eg1); - fprintf(asir_out,"sugar=%d,nsp=%d,nred=%d,spmat=(%d,%d),rank=%d ", - sugar,nsp,nred,sprow,spcol,rank); - fprintf(asir_out,"%fsec\n",eg_f4.exectime+eg_f4.gctime); - + if ( DP_Print ) + fprintf(asir_out,"sugar=%d,symb=%fsec,", + sugar,eg_f4.exectime+eg_f4.gctime); + if ( 1 ) + nflist = nd_f4_red(m,l,s0vect,col,rp0); + else + nflist = nd_f4_red_dist(m,l,s0vect,col,rp0); /* adding new bases */ - for ( i = 0; i < rank; i++ ) { - nf = vect_to_ndv(spmat[i],spcol,col,rhead,s0vect); + for ( r = nflist; r; r = NEXT(r) ) { + nf = (NDV)BDY(r); SG(nf) = sugar; ndv_removecont(m,nf); nh = ndv_newps(nf,0); @@ -3903,109 +4031,327 @@ NODE nd_f4(int m) } } for ( r = g; r; r = NEXT(r) ) BDY(r) = (pointer)nd_ps[(int)BDY(r)]; + fprintf(asir_out,"ndv_alloc=%d\n",ndv_alloc); return g; } -#else -NODE nd_f4(int m) + +NODE nd_f4_red(int m,ND_pairs sp0,UINT *s0vect,int col,NODE rp0) { - int i,nh,stat,index; - NODE r,g; - ND_pairs d,l,t; - ND spol,red; - NDV nf,redv; - NM s0,s; - NODE sp0,sp,rp0,rp; - int nsp,nred,col,rank,len,k,j,a; - UINT c; - UINT **spmat; - UINT *s0vect,*svect,*p,*v; - int *colstat; IndArray *imat; + int nsp,nred,spcol,sprow,a; int *rhead; - int spcol,sprow; - int sugar; - PGeoBucket bucket; + int i,j,k,l,rank; + NODE rp,r0,r; + ND_pairs sp; + ND spol; + int **spmat; + UINT *svect,*v; + int *colstat; struct oEGT eg0,eg1,eg_f4; + NM_ind_pair *rvect; - if ( !m ) - error("nd_f4 : not implemented"); + get_eg(&eg0); + for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ ); + nred = length(rp0); spcol = col-nred; + imat = (IndArray *)ALLOCA(nred*sizeof(IndArray)); + rhead = (int *)ALLOCA(col*sizeof(int)); + for ( i = 0; i < col; i++ ) rhead[i] = 0; - g = 0; d = 0; - for ( i = 0; i < nd_psn; i++ ) { - d = update_pairs(d,g,i); - g = update_base(g,i); + /* construction of index arrays */ + rvect = (NM_ind_pair *)ALLOCA(nred*sizeof(NM_ind_pair)); + for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) { + rvect[i] = (NM_ind_pair)BDY(rp); + imat[i] = nm_ind_pair_to_vect_compress(m,s0vect,col,rvect[i]); + rhead[imat[i]->head] = 1; } - while ( d ) { - get_eg(&eg0); - l = nd_minsugarp(d,&d); - sugar = SG(l); - bucket = create_pbucket(); - stat = nd_sp_f4(m,l,bucket,&sp0); - if ( !stat ) { - for ( t = l; NEXT(t); t = NEXT(t) ); - NEXT(t) = d; d = l; - d = nd_reconstruct(m,0,d); - continue; + + /* elimination (1st step) */ + spmat = (int **)ALLOCA(nsp*sizeof(UINT *)); + svect = (UINT *)ALLOCA(col*sizeof(UINT)); + for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) { + nd_sp(m,0,sp,&spol); + if ( !spol ) continue; + nd_to_vect(m,s0vect,col,spol,svect); + nd_free(spol); + if ( m == -1 ) ndv_reduce_vect_sf(m,svect,col,imat,rvect,nred); + else ndv_reduce_vect(m,svect,col,imat,rvect,nred); + for ( i = 0; i < col; i++ ) if ( svect[i] ) break; + if ( i < col ) { + spmat[sprow] = v = (UINT *)MALLOC_ATOMIC(spcol*sizeof(UINT)); + for ( j = k = 0; j < col; j++ ) + if ( !rhead[j] ) v[k++] = svect[j]; + sprow++; } - if ( !sp0 ) continue; - col = nd_symbolic_preproc(bucket,&s0vect,&rp0); - if ( !col ) { - for ( t = l; NEXT(t); t = NEXT(t) ); - NEXT(t) = d; d = l; - d = nd_reconstruct(m,0,d); - continue; + } + /* free index arrays */ + for ( i = 0; i < nred; i++ ) GC_free(imat[i]->index.c); + + /* elimination (2nd step) */ + colstat = (int *)ALLOCA(spcol*sizeof(int)); + if ( m == -1 ) + rank = generic_gauss_elim_sf(spmat,sprow,spcol,m,colstat); + else + rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat); + r0 = 0; + for ( i = 0; i < rank; i++ ) { + NEXTNODE(r0,r); BDY(r) = + (pointer)vect_to_ndv(spmat[i],spcol,col,rhead,s0vect); + GC_free(spmat[i]); + } + for ( ; i < sprow; i++ ) GC_free(spmat[i]); + get_eg(&eg1); init_eg(&eg_f4); add_eg(&eg_f4,&eg0,&eg1); + if ( DP_Print ) { + fprintf(asir_out,"nsp=%d,nred=%d,spmat=(%d,%d),rank=%d ", + nsp,nred,sprow,spcol,rank); + fprintf(asir_out,"%fsec\n",eg_f4.exectime+eg_f4.gctime); + } + return r0; +} + +FILE *nd_write,*nd_read; + +void nd_send_int(int a) { + write_int(nd_write,&a); +} + +void nd_send_intarray(int *p,int len) { + write_intarray(nd_write,p,len); +} + +int nd_recv_int() { + int a; + + read_int(nd_read,&a); + return a; +} + +void nd_recv_intarray(int *p,int len) { + read_intarray(nd_read,p,len); +} + +void nd_send_ndv(NDV p) { + int len,i; + NMV m; + + if ( !p ) nd_send_int(0); + else { + len = LEN(p); + nd_send_int(len); + m = BDY(p); + for ( i = 0; i < len; i++, NMV_ADV(m) ) { + nd_send_int(CM(m)); + nd_send_intarray(DL(m),nd_wpd); } + } +} - nsp = length(sp0); nred = length(rp0); spcol = col-nred; - imat = (IndArray *)MALLOC(nred*sizeof(IndArray)); - rhead = (int *)MALLOC_ATOMIC(col*sizeof(int)); - for ( i = 0; i < col; i++ ) rhead[i] = 0; +void nd_send_nd(ND p) { + int len,i; + NM m; - /* construction of index arrays */ - for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) { - imat[i] = nm_ind_pair_to_vect_compress(m,s0vect,col,(NM_ind_pair)BDY(rp)); - rhead[imat[i]->head] = 1; + if ( !p ) nd_send_int(0); + else { + len = LEN(p); + nd_send_int(len); + m = BDY(p); + for ( i = 0; i < len; i++, m = NEXT(m) ) { + nd_send_int(CM(m)); + nd_send_intarray(DL(m),nd_wpd); } + } +} - /* elimination (1st step) */ - spmat = (UINT **)MALLOC(nsp*sizeof(UINT *)); - svect = (UINT *)MALLOC_ATOMIC(col*sizeof(UINT)); - for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) { - nd_to_vect(m,s0vect,col,BDY(sp),svect); - ndv_reduce_vect(m,svect,col,imat,rp0); - for ( i = 0; i < col; i++ ) if ( svect[i] ) break; - if ( i < col ) { - spmat[sprow] = v = (UINT *)MALLOC_ATOMIC(spcol*sizeof(UINT)); - for ( j = k = 0; j < col; j++ ) - if ( !rhead[j] ) v[k++] = svect[j]; - sprow++; - } +NDV nd_recv_ndv() +{ + int len,i; + NMV m,m0; + NDV r; + + len = nd_recv_int(); + if ( !len ) return 0; + else { + m0 = m = (NMV)GC_malloc_atomic_ignore_off_page(nmv_adv*len); + ndv_alloc += len*nmv_adv; + for ( i = 0; i < len; i++, NMV_ADV(m) ) { + CM(m) = nd_recv_int(); + nd_recv_intarray(DL(m),nd_wpd); } - /* free index arrays */ - for ( i = 0; i < nred; i++ ) GC_free(imat[i]->index.c); + MKNDV(nd_nvar,m0,len,r); + return r; + } +} - /* elimination (2nd step) */ - colstat = (int *)ALLOCA(spcol*sizeof(int)); - rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat); +int ox_exec_f4_red(Q proc) +{ + Obj obj; + STRING fname; + NODE arg; + int s; + extern int ox_need_conv,ox_file_io; - get_eg(&eg1); init_eg(&eg_f4); add_eg(&eg_f4,&eg0,&eg1); - fprintf(asir_out,"sugar=%d,nsp=%d,nred=%d,spmat=(%d,%d),rank=%d ", - sugar,nsp,nred,sprow,spcol,rank); - fprintf(asir_out,"%fsec\n",eg_f4.exectime+eg_f4.gctime); + MKSTR(fname,"nd_exec_f4_red"); + arg = mknode(2,proc,fname); + Pox_cmo_rpc(arg,&obj); + s = get_ox_server_id(QTOS(proc)); + nd_write = iofp[s].out; + nd_read = iofp[s].in; + ox_need_conv = ox_file_io = 0; + return s; +} - /* adding new bases */ - for ( i = 0; i < rank; i++ ) { - nf = vect_to_ndv(spmat[i],spcol,col,rhead,s0vect); - SG(nf) = sugar; - ndv_removecont(m,nf); - nh = ndv_newps(nf,0); - d = update_pairs(d,g,nh); - g = update_base(g,nh); - GC_free(spmat[i]); +NODE nd_f4_red_dist(int m,ND_pairs sp0,UINT *s0vect,int col,NODE rp0) +{ + int nsp,nred; + int i,rank,s; + NODE rp,r0,r; + ND_pairs sp; + NM_ind_pair pair; + NMV nmv; + NM nm; + NDV nf; + Obj proc,dmy; + + ox_launch_main(0,0,&proc); + s = ox_exec_f4_red((Q)proc); + + nd_send_int(m); + nd_send_int(nd_nvar); + nd_send_int(nd_bpe); + nd_send_int(nd_wpd); + nd_send_int(nmv_adv); + + saveobj(nd_write,dp_current_spec.obj); fflush(nd_write); + + nd_send_int(nd_psn); + for ( i = 0; i < nd_psn; i++ ) nd_send_ndv(nd_ps[i]); + + for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ ); + nd_send_int(nsp); + for ( i = 0, sp = sp0; i < nsp; i++, sp = NEXT(sp) ) { + nd_send_int(sp->i1); nd_send_int(sp->i2); + } + + nd_send_int(col); nd_send_intarray(s0vect,col*nd_wpd); + + nred = length(rp0); nd_send_int(nred); + for ( i = 0, rp = rp0; i < nred; i++, rp = NEXT(rp) ) { + pair = (NM_ind_pair)BDY(rp); + nd_send_int(pair->index); + nd_send_intarray(pair->mul->dl,nd_wpd); + } + fflush(nd_write); + rank = nd_recv_int(); + fprintf(asir_out,"rank=%d\n",rank); + r0 = 0; + for ( i = 0; i < rank; i++ ) { + nf = nd_recv_ndv(); + NEXTNODE(r0,r); BDY(r) = (pointer)nf; + } + Pox_shutdown(mknode(1,proc),&dmy); + return r0; +} + +/* server side */ + +void nd_exec_f4_red_dist() +{ + int m,i,nsp,col,s0size,nred,spcol,j,k; + NM_ind_pair *rp0; + NDV nf; + UINT *s0vect; + IndArray *imat; + int *rhead; + int **spmat; + UINT *svect,*v; + ND_pairs *sp0; + int *colstat; + int a,sprow,rank; + struct order_spec ord; + Obj ordspec; + ND spol; + + nd_read = iofp[0].in; + nd_write = iofp[0].out; + m = nd_recv_int(); + nd_nvar = nd_recv_int(); + nd_bpe = nd_recv_int(); + nd_wpd = nd_recv_int(); + nmv_adv = nd_recv_int(); + + loadobj(nd_read,&ordspec); + create_order_spec(ordspec,&ord); + nd_init_ord(&ord); + nd_setup_parameters(nd_nvar,0); + + nd_psn = nd_recv_int(); + nd_ps = (NDV *)MALLOC(nd_psn*sizeof(NDV)); + nd_bound = (UINT **)MALLOC(nd_psn*sizeof(UINT *)); + for ( i = 0; i < nd_psn; i++ ) { + nd_ps[i] = nd_recv_ndv(); + nd_bound[i] = ndv_compute_bound(nd_ps[i]); + } + + nsp = nd_recv_int(); + sp0 = (ND_pairs *)MALLOC(nsp*sizeof(ND_pairs)); + for ( i = 0; i < nsp; i++ ) { + NEWND_pairs(sp0[i]); + sp0[i]->i1 = nd_recv_int(); sp0[i]->i2 = nd_recv_int(); + ndl_lcm(HDL(nd_ps[sp0[i]->i1]),HDL(nd_ps[sp0[i]->i2]),LCM(sp0[i])); + } + + col = nd_recv_int(); + s0size = col*nd_wpd; + s0vect = (UINT *)MALLOC(s0size*sizeof(UINT)); + nd_recv_intarray(s0vect,s0size); + + nred = nd_recv_int(); + rp0 = (NM_ind_pair *)MALLOC(nred*sizeof(NM_ind_pair)); + for ( i = 0; i < nred; i++ ) { + rp0[i] = (NM_ind_pair)MALLOC(sizeof(struct oNM_ind_pair)); + rp0[i]->index = nd_recv_int(); + rp0[i]->mul = (NM)MALLOC(sizeof(struct oNM)+(nd_wpd-1)*sizeof(UINT)); + nd_recv_intarray(rp0[i]->mul->dl,nd_wpd); + } + + spcol = col-nred; + imat = (IndArray *)MALLOC(nred*sizeof(IndArray)); + rhead = (int *)MALLOC(col*sizeof(int)); + for ( i = 0; i < col; i++ ) rhead[i] = 0; + + /* construction of index arrays */ + for ( i = 0; i < nred; i++ ) { + imat[i] = nm_ind_pair_to_vect_compress(m,s0vect,col,rp0[i]); + rhead[imat[i]->head] = 1; + } + + /* elimination (1st step) */ + spmat = (int **)MALLOC(nsp*sizeof(UINT *)); + svect = (UINT *)MALLOC(col*sizeof(UINT)); + for ( a = sprow = 0; a < nsp; a++ ) { + nd_sp(m,0,sp0[a],&spol); + if ( !spol ) continue; + nd_to_vect(m,s0vect,col,spol,svect); + nd_free(spol); + if ( m == -1 ) ndv_reduce_vect_sf(m,svect,col,imat,rp0,nred); + else ndv_reduce_vect(m,svect,col,imat,rp0,nred); + for ( i = 0; i < col; i++ ) if ( svect[i] ) break; + if ( i < col ) { + spmat[sprow] = v = (UINT *)MALLOC(spcol*sizeof(UINT)); + for ( j = k = 0; j < col; j++ ) + if ( !rhead[j] ) v[k++] = svect[j]; + sprow++; } - for ( ; i < sprow; i++ ) GC_free(spmat[i]); } - for ( r = g; r; r = NEXT(r) ) BDY(r) = (pointer)nd_ps[(int)BDY(r)]; - return g; + /* elimination (2nd step) */ + colstat = (int *)ALLOCA(spcol*sizeof(int)); + if ( m == -1 ) + rank = generic_gauss_elim_sf(spmat,sprow,spcol,m,colstat); + else + rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat); + nd_send_int(rank); + for ( i = 0; i < rank; i++ ) { + nf = vect_to_ndv(spmat[i],spcol,col,rhead,s0vect); + nd_send_ndv(nf); + } + fflush(nd_write); } -#endif