version 1.64, 2003/09/12 01:12:40 |
version 1.68, 2003/09/12 15:07:11 |
|
|
/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.63 2003/09/11 09:03:53 noro Exp $ */ |
/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.67 2003/09/12 14:51:27 noro Exp $ */ |
|
|
#include "ca.h" |
#include "ca.h" |
#include "inline.h" |
#include "inline.h" |
|
|
typedef unsigned int UINT; |
typedef unsigned int UINT; |
|
|
#define USE_GEOBUCKET 1 |
#define USE_GEOBUCKET 1 |
|
#define USE_UNROLL 1 |
|
|
#define REDTAB_LEN 32003 |
#define REDTAB_LEN 32003 |
|
|
Line 102 typedef struct oNM_ind_pair |
|
Line 103 typedef struct oNM_ind_pair |
|
int index; |
int index; |
} *NM_ind_pair; |
} *NM_ind_pair; |
|
|
|
typedef struct oIndArray |
|
{ |
|
char width; |
|
int head; |
|
union { |
|
unsigned char *c; |
|
unsigned short *s; |
|
unsigned int *i; |
|
} index; |
|
} *IndArray; |
|
|
int (*ndl_compare_function)(UINT *a1,UINT *a2); |
int (*ndl_compare_function)(UINT *a1,UINT *a2); |
|
|
Line 348 P ndvtop(int mod,VL vl,VL dvl,NDV p); |
|
Line 359 P ndvtop(int mod,VL vl,VL dvl,NDV p); |
|
NDV ndtondv(int mod,ND p); |
NDV ndtondv(int mod,ND p); |
ND ndvtond(int mod,NDV p); |
ND ndvtond(int mod,NDV p); |
int nm_ind_pair_to_vect(int m,UINT *s0,int n,NM_ind_pair pair,UINT *r); |
int nm_ind_pair_to_vect(int m,UINT *s0,int n,NM_ind_pair pair,UINT *r); |
void nm_ind_pair_to_vect_compress(int m,UINT *s0,int n,NM_ind_pair pair,UINT *r,int *ind); |
IndArray nm_ind_pair_to_vect_compress(int m,UINT *s0,int n,NM_ind_pair pair); |
int nd_to_vect(int mod,UINT *s0,int n,ND d,UINT *r); |
int nd_to_vect(int mod,UINT *s0,int n,ND d,UINT *r); |
|
|
void nd_free_private_storage() |
void nd_free_private_storage() |
Line 413 INLINE int ndl_reducible(UINT *d1,UINT *d2) |
|
Line 424 INLINE int ndl_reducible(UINT *d1,UINT *d2) |
|
int i,j; |
int i,j; |
|
|
if ( TD(d1) < TD(d2) ) return 0; |
if ( TD(d1) < TD(d2) ) return 0; |
|
#if USE_UNROLL |
switch ( nd_bpe ) { |
switch ( nd_bpe ) { |
#if 0 |
|
case 3: |
case 3: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u1 = d1[i]; u2 = d2[i]; |
u1 = d1[i]; u2 = d2[i]; |
Line 479 INLINE int ndl_reducible(UINT *d1,UINT *d2) |
|
Line 490 INLINE int ndl_reducible(UINT *d1,UINT *d2) |
|
if ( d1[i] < d2[i] ) return 0; |
if ( d1[i] < d2[i] ) return 0; |
return 1; |
return 1; |
break; |
break; |
#endif |
|
default: |
default: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u1 = d1[i]; u2 = d2[i]; |
u1 = d1[i]; u2 = d2[i]; |
Line 488 INLINE int ndl_reducible(UINT *d1,UINT *d2) |
|
Line 498 INLINE int ndl_reducible(UINT *d1,UINT *d2) |
|
} |
} |
return 1; |
return 1; |
} |
} |
|
#else |
|
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
|
u1 = d1[i]; u2 = d2[i]; |
|
for ( j = 0; j < nd_epw; j++ ) |
|
if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0; |
|
} |
|
return 1; |
|
#endif |
} |
} |
|
|
/* |
/* |
Line 555 void ndl_lcm(UINT *d1,unsigned *d2,UINT *d) |
|
Line 573 void ndl_lcm(UINT *d1,unsigned *d2,UINT *d) |
|
UINT t1,t2,u,u1,u2; |
UINT t1,t2,u,u1,u2; |
int i,j,l; |
int i,j,l; |
|
|
|
#if USE_UNROLL |
switch ( nd_bpe ) { |
switch ( nd_bpe ) { |
#if 0 |
|
case 3: |
case 3: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u1 = d1[i]; u2 = d2[i]; |
u1 = d1[i]; u2 = d2[i]; |
Line 622 void ndl_lcm(UINT *d1,unsigned *d2,UINT *d) |
|
Line 640 void ndl_lcm(UINT *d1,unsigned *d2,UINT *d) |
|
d[i] = u1>u2?u1:u2; |
d[i] = u1>u2?u1:u2; |
} |
} |
break; |
break; |
#endif |
|
default: |
default: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u1 = d1[i]; u2 = d2[i]; |
u1 = d1[i]; u2 = d2[i]; |
Line 633 void ndl_lcm(UINT *d1,unsigned *d2,UINT *d) |
|
Line 650 void ndl_lcm(UINT *d1,unsigned *d2,UINT *d) |
|
} |
} |
break; |
break; |
} |
} |
|
#else |
|
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
|
u1 = d1[i]; u2 = d2[i]; |
|
for ( j = 0, u = 0; j < nd_epw; j++ ) { |
|
t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2; |
|
} |
|
d[i] = u; |
|
} |
|
#endif |
TD(d) = ndl_weight(d); |
TD(d) = ndl_weight(d); |
if ( nd_blockmask ) ndl_weight_mask(d); |
if ( nd_blockmask ) ndl_weight_mask(d); |
} |
} |
Line 838 int ndl_disjoint(UINT *d1,UINT *d2) |
|
Line 864 int ndl_disjoint(UINT *d1,UINT *d2) |
|
UINT t1,t2,u,u1,u2; |
UINT t1,t2,u,u1,u2; |
int i,j; |
int i,j; |
|
|
|
#if USE_UNROLL |
switch ( nd_bpe ) { |
switch ( nd_bpe ) { |
#if 0 |
|
case 3: |
case 3: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u1 = d1[i]; u2 = d2[i]; |
u1 = d1[i]; u2 = d2[i]; |
Line 904 int ndl_disjoint(UINT *d1,UINT *d2) |
|
Line 930 int ndl_disjoint(UINT *d1,UINT *d2) |
|
if ( d1[i] && d2[i] ) return 0; |
if ( d1[i] && d2[i] ) return 0; |
return 1; |
return 1; |
break; |
break; |
#endif |
|
default: |
default: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u1 = d1[i]; u2 = d2[i]; |
u1 = d1[i]; u2 = d2[i]; |
Line 916 int ndl_disjoint(UINT *d1,UINT *d2) |
|
Line 941 int ndl_disjoint(UINT *d1,UINT *d2) |
|
return 1; |
return 1; |
break; |
break; |
} |
} |
|
#else |
|
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
|
u1 = d1[i]; u2 = d2[i]; |
|
for ( j = 0; j < nd_epw; j++ ) { |
|
if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0; |
|
u1 >>= nd_bpe; u2 >>= nd_bpe; |
|
} |
|
} |
|
return 1; |
|
#endif |
} |
} |
|
|
int ndl_check_bound2(int index,UINT *d2) |
int ndl_check_bound2(int index,UINT *d2) |
Line 926 int ndl_check_bound2(int index,UINT *d2) |
|
Line 961 int ndl_check_bound2(int index,UINT *d2) |
|
|
|
d1 = nd_bound[index]; |
d1 = nd_bound[index]; |
ind = 0; |
ind = 0; |
|
#if USE_UNROLL |
switch ( nd_bpe ) { |
switch ( nd_bpe ) { |
#if 0 |
|
case 3: |
case 3: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u2 = d2[i]; |
u2 = d2[i]; |
Line 992 int ndl_check_bound2(int index,UINT *d2) |
|
Line 1027 int ndl_check_bound2(int index,UINT *d2) |
|
if ( d1[i]+d2[i]<d1[i] ) return 1; |
if ( d1[i]+d2[i]<d1[i] ) return 1; |
return 0; |
return 0; |
break; |
break; |
#endif |
|
default: |
default: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u2 = d2[i]; |
u2 = d2[i]; |
Line 1003 int ndl_check_bound2(int index,UINT *d2) |
|
Line 1037 int ndl_check_bound2(int index,UINT *d2) |
|
return 0; |
return 0; |
break; |
break; |
} |
} |
|
#else |
|
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
|
u2 = d2[i]; |
|
k = (nd_epw-1)*nd_bpe; |
|
for ( j = 0; j < nd_epw; j++, k -= nd_bpe ) |
|
if ( d1[ind++]+((u2>>k)&nd_mask0) > nd_mask0 ) return 1; |
|
} |
|
return 0; |
|
#endif |
} |
} |
|
|
int ndl_check_bound2_direct(UINT *d1,UINT *d2) |
int ndl_check_bound2_direct(UINT *d1,UINT *d2) |
Line 1011 int ndl_check_bound2_direct(UINT *d1,UINT *d2) |
|
Line 1054 int ndl_check_bound2_direct(UINT *d1,UINT *d2) |
|
int i,j,ind,k; |
int i,j,ind,k; |
|
|
ind = 0; |
ind = 0; |
|
#if USE_UNROLL |
switch ( nd_bpe ) { |
switch ( nd_bpe ) { |
#if 0 |
|
case 3: |
case 3: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u2 = d2[i]; |
u2 = d2[i]; |
Line 1077 int ndl_check_bound2_direct(UINT *d1,UINT *d2) |
|
Line 1120 int ndl_check_bound2_direct(UINT *d1,UINT *d2) |
|
if ( d1[i]+d2[i]<d1[i] ) return 1; |
if ( d1[i]+d2[i]<d1[i] ) return 1; |
return 0; |
return 0; |
break; |
break; |
#endif |
|
default: |
default: |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
u2 = d2[i]; |
u2 = d2[i]; |
Line 1088 int ndl_check_bound2_direct(UINT *d1,UINT *d2) |
|
Line 1130 int ndl_check_bound2_direct(UINT *d1,UINT *d2) |
|
return 0; |
return 0; |
break; |
break; |
} |
} |
|
#else |
|
for ( i = nd_exporigin; i < nd_wpd; i++ ) { |
|
u2 = d2[i]; |
|
k = (nd_epw-1)*nd_bpe; |
|
for ( j = 0; j < nd_epw; j++, k -= nd_bpe ) |
|
if ( d1[ind++]+((u2>>k)&nd_mask0) > nd_mask0 ) return 1; |
|
} |
|
return 0; |
|
#endif |
} |
} |
|
|
INLINE int ndl_hash_value(UINT *d) |
INLINE int ndl_hash_value(UINT *d) |
Line 3746 int nm_ind_pair_to_vect(int mod,UINT *s0,int n,NM_ind_ |
|
Line 3797 int nm_ind_pair_to_vect(int mod,UINT *s0,int n,NM_ind_ |
|
return i; |
return i; |
} |
} |
|
|
void nm_ind_pair_to_vect_compress(int mod,UINT *s0,int n,NM_ind_pair pair,UINT *r,int *ind) |
IndArray nm_ind_pair_to_vect_compress(int mod,UINT *s0,int n,NM_ind_pair pair) |
{ |
{ |
NM m; |
NM m; |
NMV mr; |
NMV mr; |
UINT *d,*t,*s; |
UINT *d,*t,*s; |
NDV p; |
NDV p; |
int i,j,len; |
unsigned char *ivc; |
|
unsigned short *ivs; |
|
UINT *v,*ivi; |
|
int i,j,len,prev,diff,cdiff; |
|
IndArray r; |
|
|
m = pair->mul; |
m = pair->mul; |
d = DL(m); |
d = DL(m); |
p = nd_ps[pair->index]; |
p = nd_ps[pair->index]; |
len = LEN(p); |
len = LEN(p); |
t = (UINT *)ALLOCA(nd_wpd*sizeof(UINT)); |
t = (UINT *)ALLOCA(nd_wpd*sizeof(UINT)); |
|
r = (IndArray)MALLOC(sizeof(struct oIndArray)); |
|
v = (unsigned int *)ALLOCA(len*sizeof(unsigned int)); |
for ( i = j = 0, s = s0, mr = BDY(p); j < len; j++, NMV_ADV(mr) ) { |
for ( i = j = 0, s = s0, mr = BDY(p); j < len; j++, NMV_ADV(mr) ) { |
ndl_add(d,DL(mr),t); |
ndl_add(d,DL(mr),t); |
for ( ; !ndl_equal(t,s); s += nd_wpd, i++ ); |
for ( ; !ndl_equal(t,s); s += nd_wpd, i++ ); |
r[j] = CM(mr); |
v[j] = i; |
ind[j] = i; |
|
} |
} |
|
r->head = v[0]; |
|
diff = 0; |
|
for ( i = 1; i < len; i++ ) { |
|
cdiff = v[i]-v[i-1]; diff = MAX(cdiff,diff); |
|
} |
|
if ( diff < 256 ) { |
|
r->width = 1; |
|
ivc = (unsigned char *)MALLOC_ATOMIC(len*sizeof(unsigned char)); |
|
r->index.c = ivc; |
|
for ( i = 1, ivc[0] = 0; i < len; i++ ) ivc[i] = v[i]-v[i-1]; |
|
} else if ( diff < 65536 ) { |
|
r->width = 2; |
|
ivs = (unsigned short *)MALLOC_ATOMIC(len*sizeof(unsigned short)); |
|
r->index.s = ivs; |
|
for ( i = 1, ivs[0] = 0; i < len; i++ ) ivs[i] = v[i]-v[i-1]; |
|
} else { |
|
r->width = 4; |
|
ivi = (unsigned int *)MALLOC_ATOMIC(len*sizeof(unsigned int)); |
|
r->index.i = ivi; |
|
for ( i = 1, ivi[0] = 0; i < len; i++ ) ivi[i] = v[i]-v[i-1]; |
|
} |
|
return r; |
} |
} |
|
|
|
|
|
void ndv_reduce_vect(int m,UINT *svect,int col,IndArray *imat,NODE rp0) |
|
{ |
|
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 ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) { |
|
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]; |
|
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]; c1 = CM(mr); c2 = svect[pos]; |
|
prev = pos; |
|
DMA(c1,c,c2,up,lo); |
|
if ( up ) { DSAB(m,up,lo,dmy,c3); svect[pos] = c3; |
|
} else svect[pos] = lo; |
|
} |
|
break; |
|
case 2: |
|
ivs = ivect->index.s; |
|
for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { |
|
pos = prev+ivs[j]; c1 = CM(mr); c2 = svect[pos]; |
|
prev = pos; |
|
DMA(c1,c,c2,up,lo); |
|
if ( up ) { DSAB(m,up,lo,dmy,c3); svect[pos] = c3; |
|
} else svect[pos] = lo; |
|
} |
|
break; |
|
case 4: |
|
ivi = ivect->index.i; |
|
for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { |
|
pos = prev+ivi[j]; c1 = CM(mr); c2 = svect[pos]; |
|
prev = pos; |
|
DMA(c1,c,c2,up,lo); |
|
if ( up ) { DSAB(m,up,lo,dmy,c3); svect[pos] = c3; |
|
} else svect[pos] = lo; |
|
} |
|
break; |
|
} |
|
} |
|
} |
|
for ( i = 0; i < col; i++ ) |
|
if ( svect[i] >= (UINT)m ) svect[i] %= m; |
|
} |
|
|
|
NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhead,UINT *s0vect) |
|
{ |
|
int j,k,len; |
|
UINT *p; |
|
UINT c; |
|
NDV r; |
|
NMV mr0,mr; |
|
|
|
for ( j = 0, len = 0; j < spcol; j++ ) if ( vect[j] ) len++; |
|
if ( !len ) return 0; |
|
else { |
|
mr0 = (NMV)MALLOC_ATOMIC(nmv_adv*len); |
|
mr = mr0; |
|
p = s0vect; |
|
for ( j = k = 0; j < col; j++, p += nd_wpd ) |
|
if ( !rhead[j] ) { |
|
if ( c = vect[k++] ) { |
|
ndl_copy(p,DL(mr)); CM(mr) = c; NMV_ADV(mr); |
|
} |
|
} |
|
MKNDV(nd_nvar,mr0,len,r); |
|
return r; |
|
} |
|
} |
|
|
|
int nd_sp_f4(int m,ND_pairs l,PGeoBucket bucket,NODE *s) |
|
{ |
|
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; |
|
} |
|
|
|
int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vect,NODE *r) |
|
{ |
|
NODE rp0,rp; |
|
NM mul,head,s0,s; |
|
int index,col,i; |
|
RHist h; |
|
UINT *s0v,*p; |
|
NM_ind_pair pair; |
|
ND red; |
|
|
|
s0 = 0; rp0 = 0; col = 0; |
|
while ( 1 ) { |
|
head = remove_head_pbucket_symbolic(bucket); |
|
if ( !head ) break; |
|
if ( !s0 ) s0 = head; |
|
else NEXT(s) = head; |
|
s = head; |
|
index = ndl_find_reducer(DL(head)); |
|
if ( index >= 0 ) { |
|
h = nd_psh[index]; |
|
NEWNM(mul); |
|
ndl_sub(DL(head),DL(h),DL(mul)); |
|
if ( ndl_check_bound2(index,DL(mul)) ) return 0; |
|
MKNM_ind_pair(pair,mul,index); |
|
red = ndv_mul_nm_symbolic(mul,nd_ps[index]); |
|
add_pbucket_symbolic(bucket,nd_remove_head(red)); |
|
NEXTNODE(rp0,rp); BDY(rp) = (pointer)pair; |
|
} |
|
col++; |
|
} |
|
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); |
|
*s0vect = s0v; |
|
*r = rp0; |
|
return col; |
|
} |
|
|
NODE nd_f4(int m) |
NODE nd_f4(int m) |
{ |
{ |
int i,nh,stat,index; |
int i,nh,stat,index; |
NODE r,g; |
NODE r,g; |
ND_pairs d,l,t; |
ND_pairs d,l,t; |
ND spol,red; |
ND spol,red; |
NDV nf; |
NDV nf,redv; |
NM_ind_pair pair,pair1,pair2; |
|
NM s0,s; |
NM s0,s; |
NODE sp0,sp,rp0,rp; |
NODE sp0,sp,rp0,rp; |
RHist h; |
int nsp,nred,col,rank,len,k,j,a; |
int nsp,nred,col,rank,len,k,j; |
UINT c; |
NMV mr0,mr; |
|
UINT c,c1,c2,c3; |
|
NM mul,head; |
|
UINT **spmat; |
UINT **spmat; |
UINT *s0vect,*rvect,*svect,*p,*ivect; |
UINT *s0vect,*svect,*p,*v; |
int *colstat; |
int *colstat; |
int sugar,pos; |
IndArray *imat; |
|
int *rhead; |
|
int spcol,sprow; |
|
int sugar; |
PGeoBucket bucket; |
PGeoBucket bucket; |
int t_0,t_1,t_2,t_3,t_4,t_5,t_c,t_e,t_20,t_21; |
struct oEGT eg0,eg1,eg_f4; |
|
|
if ( !m ) |
if ( !m ) |
error("nd_f4 : not implemented"); |
error("nd_f4 : not implemented"); |
Line 3798 NODE nd_f4(int m) |
|
Line 4017 NODE nd_f4(int m) |
|
g = update_base(g,i); |
g = update_base(g,i); |
} |
} |
while ( d ) { |
while ( d ) { |
|
get_eg(&eg0); |
l = nd_minsugarp(d,&d); |
l = nd_minsugarp(d,&d); |
sugar = SG(l); |
sugar = SG(l); |
sp0 = 0; |
|
bucket = create_pbucket(); |
bucket = create_pbucket(); |
t_0 = clock(); |
stat = nd_sp_f4(m,l,bucket,&sp0); |
for ( t = l, nsp = 0; t; t = NEXT(t) ) { |
if ( !stat ) { |
stat = nd_sp(m,0,t,&spol); |
for ( t = l; NEXT(t); t = NEXT(t) ); |
if ( !stat ) goto again; |
NEXT(t) = d; d = l; |
if ( spol ) { |
d = nd_reconstruct(m,0,d); |
NEXTNODE(sp0,sp); BDY(sp) = (pointer)nd_dup(spol); |
continue; |
add_pbucket_symbolic(bucket,spol); |
|
nsp++; |
|
} |
|
} |
} |
if ( sp0 ) NEXT(sp) = 0; |
if ( !sp0 ) continue; |
s0 = 0; |
col = nd_symbolic_preproc(bucket,&s0vect,&rp0); |
rp0 = 0; |
if ( !col ) { |
while ( 1 ) { |
for ( t = l; NEXT(t); t = NEXT(t) ); |
head = remove_head_pbucket_symbolic(bucket); |
NEXT(t) = d; d = l; |
if ( !head ) break; |
d = nd_reconstruct(m,0,d); |
if ( !s0 ) s0 = head; |
continue; |
else NEXT(s) = head; |
|
s = head; |
|
index = ndl_find_reducer(DL(head)); |
|
if ( index >= 0 ) { |
|
h = nd_psh[index]; |
|
NEWNM(mul); |
|
ndl_sub(DL(head),DL(h),DL(mul)); |
|
if ( ndl_check_bound2(index,DL(mul)) ) goto again; |
|
MKNM_ind_pair(pair,mul,index); |
|
red = ndv_mul_nm_symbolic(mul,nd_ps[index]); |
|
add_pbucket_symbolic(bucket,nd_remove_head(red)); |
|
NEXTNODE(rp0,rp); BDY(rp) = (pointer)pair; |
|
nred++; |
|
} |
|
} |
} |
t_1 = clock(); |
|
if ( s0 ) NEXT(s) = 0; |
nsp = length(sp0); nred = length(rp0); spcol = col-nred; |
for ( i = 0, s = s0; s; s = NEXT(s), i++ ); |
imat = (IndArray *)MALLOC(nred*sizeof(IndArray)); |
col = i; |
rhead = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
s0vect = (UINT *)MALLOC(col*nd_wpd*sizeof(UINT)); |
for ( i = 0; i < col; i++ ) rhead[i] = 0; |
for ( i = 0, p = s0vect, s = s0; i < col; |
|
i++, p += nd_wpd, s = NEXT(s) ) ndl_copy(DL(s),p); |
/* construction of index arrays */ |
spmat = (UINT **)MALLOC(nsp*sizeof(UINT *)); |
for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) { |
for ( i = 0, sp = sp0; i < nsp; i++, sp = NEXT(sp) ) { |
imat[i] = nm_ind_pair_to_vect_compress(m,s0vect,col,(NM_ind_pair)BDY(rp)); |
spmat[i] = (UINT *)MALLOC(col*sizeof(UINT)); |
rhead[imat[i]->head] = 1; |
nd_to_vect(m,s0vect,col,BDY(sp),spmat[i]); |
|
} |
} |
t_2 = clock(); |
|
t_c = 0; |
/* elimination (1st step) */ |
t_e = 0; |
spmat = (UINT **)MALLOC(nsp*sizeof(UINT *)); |
rvect = (UINT *)ALLOCA(col*sizeof(UINT)); |
svect = (UINT *)MALLOC_ATOMIC(col*sizeof(UINT)); |
ivect = (int *)ALLOCA(col*sizeof(int)); |
for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) { |
for ( rp = rp0; rp; rp = NEXT(rp) ) { |
nd_to_vect(m,s0vect,col,BDY(sp),svect); |
t_20 = clock(); |
ndv_reduce_vect(m,svect,col,imat,rp0); |
nm_ind_pair_to_vect_compress(m,s0vect,col,(NM_ind_pair)BDY(rp),rvect,ivect); |
for ( i = 0; i < col; i++ ) if ( svect[i] ) break; |
t_21 = clock(); |
if ( i < col ) { |
t_c += t_21-t_20; |
spmat[sprow] = v = (UINT *)MALLOC_ATOMIC(spcol*sizeof(UINT)); |
k = ivect[0]; |
for ( j = k = 0; j < col; j++ ) |
len = LEN(nd_ps[((NM_ind_pair)BDY(rp))->index]); |
if ( !rhead[j] ) v[k++] = svect[j]; |
for ( i = 0; i < nsp; i++ ) { |
sprow++; |
svect = spmat[i]; |
|
if ( c = svect[k] ) |
|
for ( j = 0; j < len; j++ ) { |
|
pos = ivect[j]; |
|
c1 = m-rvect[j]; |
|
c2 = svect[pos]; DMAR(c1,c,c2,m,c3); |
|
svect[pos] = c3; |
|
} |
|
} |
} |
t_e += clock()-t_21; |
|
} |
} |
t_4 = clock(); |
/* free index arrays */ |
colstat = (int *)ALLOCA(col*sizeof(int)); |
for ( i = 0; i < nred; i++ ) GC_free(imat[i]->index.c); |
rank = generic_gauss_elim_mod(spmat,nsp,col,m,colstat); |
|
t_5 = clock(); |
|
fprintf(asir_out,"sugar=%d,rank=%d ",sugar,rank); |
|
fprintf(asir_out,"symb=%f,conv1=%f,conv2=%f,elim1=%f,elim2=%f\n", |
|
(t_1-t_0)/(double)CLOCKS_PER_SEC, |
|
(t_2-t_1)/(double)CLOCKS_PER_SEC, |
|
(t_c)/(double)CLOCKS_PER_SEC, |
|
(t_e)/(double)CLOCKS_PER_SEC, |
|
(t_5-t_4)/(double)CLOCKS_PER_SEC); |
|
if ( rank ) |
|
for ( j = 0, i = 0; j < col; j++ ) { |
|
if ( colstat[j] ) { |
|
for ( k = j, len = 0; k < col; k++ ) |
|
if ( spmat[i][k] ) len++; |
|
mr0 = (NMV)MALLOC_ATOMIC(nmv_adv*len); |
|
mr = mr0; |
|
p = s0vect+nd_wpd*j; |
|
for ( k = j; k < col; k++, p += nd_wpd ) |
|
if ( spmat[i][k] ) { |
|
ndl_copy(p,DL(mr)); |
|
CM(mr) = spmat[i][k]; |
|
NMV_ADV(mr); |
|
} |
|
MKNDV(nd_nvar,mr0,len,nf); |
|
SG(nf) = sugar; |
|
ndv_removecont(m,nf); |
|
nh = ndv_newps(nf,0); |
|
d = update_pairs(d,g,nh); |
|
g = update_base(g,nh); |
|
i++; |
|
} |
|
} |
|
continue; |
|
|
|
again: |
/* elimination (2nd step) */ |
for ( t = l; NEXT(t); t = NEXT(t) ); |
colstat = (int *)ALLOCA(spcol*sizeof(int)); |
NEXT(t) = d; d = l; |
rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat); |
d = nd_reconstruct(m,0,d); |
|
NEWNM(mul); |
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); |
|
|
|
/* 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]); |
|
} |
|
for ( ; i < sprow; i++ ) GC_free(spmat[i]); |
} |
} |
for ( r = g; r; r = NEXT(r) ) BDY(r) = (pointer)nd_ps[(int)BDY(r)]; |
for ( r = g; r; r = NEXT(r) ) BDY(r) = (pointer)nd_ps[(int)BDY(r)]; |
return g; |
return g; |