[BACK]Return to nd.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/nd.c between version 1.65 and 1.68

version 1.65, 2003/09/12 08:01:51 version 1.68, 2003/09/12 15:07:11
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.64 2003/09/12 01:12:40 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"
Line 103  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 349  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,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 3787  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,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++ );
                 ind[j] = i;                  v[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 **imat,NODE rp0)  void ndv_reduce_vect(int m,UINT *svect,int col,IndArray *imat,NODE rp0)
 {  {
         int i,j,k,len,pos;          int i,j,k,len,pos,prev;
         UINT c,c1,c2,c3;          UINT c,c1,c2,c3,up,lo,dmy;
         UINT *ivect;          IndArray ivect;
           unsigned char *ivc;
           unsigned short *ivs;
           unsigned int *ivi;
         NDV redv;          NDV redv;
         NMV mr0,mr;          NMV mr;
         NODE rp;          NODE rp;
   
         for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) {          for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) {
                 ivect = imat[i];                  ivect = imat[i];
                 k = ivect[0];                  k = ivect->head; svect[k] %= m;
                 if ( c = svect[k] ) {                  if ( c = svect[k] ) {
                         c = m-c;                          c = m-c; redv = nd_ps[((NM_ind_pair)BDY(rp))->index];
                         redv = nd_ps[((NM_ind_pair)BDY(rp))->index];                          len = LEN(redv); mr = BDY(redv);
                         len = LEN(redv);                          svect[k] = 0; prev = k;
                         mr0 = BDY(redv);                          switch ( ivect->width ) {
                         for ( j = 0, mr = mr0; j < len; j++, NMV_ADV(mr) ) {                                  case 1:
                                 pos = ivect[j];                                          ivc = ivect->index.c;
                                 c1 = CM(mr);                                          for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {
                                 c2 = svect[pos]; DMAR(c1,c,c2,m,c3);                                                  pos = prev+ivc[j]; c1 = CM(mr); c2 = svect[pos];
                                 svect[pos] = c3;                                                  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)  NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhead,UINT *s0vect)
Line 3860  NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea
Line 3928  NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea
         }          }
 }  }
   
 NODE nd_sp_f4(int m,ND_pairs l,PGeoBucket bucket)  int nd_sp_f4(int m,ND_pairs l,PGeoBucket bucket,NODE *s)
 {  {
         ND_pairs t;          ND_pairs t;
         NODE sp0,sp;          NODE sp0,sp;
Line 3876  NODE nd_sp_f4(int m,ND_pairs l,PGeoBucket bucket)
Line 3944  NODE nd_sp_f4(int m,ND_pairs l,PGeoBucket bucket)
                         add_pbucket_symbolic(bucket,spol);                          add_pbucket_symbolic(bucket,spol);
                 }                  }
         }          }
         return sp0;          *s = sp0;
           return 1;
 }  }
   
 int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vect,NODE *r)  int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vect,NODE *r)
Line 3932  NODE nd_f4(int m)
Line 4001  NODE nd_f4(int m)
         UINT **spmat;          UINT **spmat;
         UINT *s0vect,*svect,*p,*v;          UINT *s0vect,*svect,*p,*v;
         int *colstat;          int *colstat;
         int **imat;          IndArray *imat;
         int *rhead;          int *rhead;
         int spcol,sprow;          int spcol,sprow;
         int sugar;          int sugar;
         PGeoBucket bucket;          PGeoBucket bucket;
           struct oEGT eg0,eg1,eg_f4;
   
         if ( !m )          if ( !m )
                 error("nd_f4 : not implemented");                  error("nd_f4 : not implemented");
Line 3947  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);
                 bucket = create_pbucket();                  bucket = create_pbucket();
                 if ( !(sp0 = nd_sp_f4(m,l,bucket))                  stat = nd_sp_f4(m,l,bucket,&sp0);
                         || !(col = nd_symbolic_preproc(bucket,&s0vect,&rp0)) ) {                  if ( !stat ) {
                         for ( t = l; NEXT(t); t = NEXT(t) );                          for ( t = l; NEXT(t); t = NEXT(t) );
                         NEXT(t) = d; d = l;                          NEXT(t) = d; d = l;
                         d = nd_reconstruct(m,0,d);                          d = nd_reconstruct(m,0,d);
                         continue;                          continue;
                 }                  }
                   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;
                   }
   
                 nsp = length(sp0); nred = length(rp0); spcol = col-nred;                  nsp = length(sp0); nred = length(rp0); spcol = col-nred;
                 imat = (int **)MALLOC(nred*sizeof(int *));                  imat = (IndArray *)MALLOC(nred*sizeof(IndArray));
                 rhead = (int *)MALLOC_ATOMIC(col*sizeof(int));                  rhead = (int *)MALLOC_ATOMIC(col*sizeof(int));
                 for ( i = 0; i < col; i++ ) rhead[i] = 0;                  for ( i = 0; i < col; i++ ) rhead[i] = 0;
   
                 /* construction of index arrays */                  /* construction of index arrays */
                 for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) {                  for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) {
                         redv = nd_ps[((NM_ind_pair)BDY(rp))->index];                          imat[i] = nm_ind_pair_to_vect_compress(m,s0vect,col,(NM_ind_pair)BDY(rp));
                         imat[i] = (int *)MALLOC_ATOMIC(LEN(redv)*sizeof(int));                          rhead[imat[i]->head] = 1;
                         nm_ind_pair_to_vect_compress(m,s0vect,col,  
                                 (NM_ind_pair)BDY(rp),imat[i]);  
                         rhead[imat[i][0]] = 1;  
                 }                  }
   
                 /* elimination (1st step) */                  /* elimination (1st step) */
Line 3977  NODE nd_f4(int m)
Line 4053  NODE nd_f4(int m)
                 svect = (UINT *)MALLOC_ATOMIC(col*sizeof(UINT));                  svect = (UINT *)MALLOC_ATOMIC(col*sizeof(UINT));
                 for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) {                  for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) {
                         nd_to_vect(m,s0vect,col,BDY(sp),svect);                          nd_to_vect(m,s0vect,col,BDY(sp),svect);
                         ndv_reduce_vect(m,svect,imat,rp0);                          ndv_reduce_vect(m,svect,col,imat,rp0);
                         for ( i = 0; i < col; i++ ) if ( svect[i] ) break;                          for ( i = 0; i < col; i++ ) if ( svect[i] ) break;
                         if ( i < col ) {                          if ( i < col ) {
                                 spmat[sprow] = v = (UINT *)MALLOC_ATOMIC(spcol*sizeof(UINT));                                  spmat[sprow] = v = (UINT *)MALLOC_ATOMIC(spcol*sizeof(UINT));
Line 3987  NODE nd_f4(int m)
Line 4063  NODE nd_f4(int m)
                         }                          }
                 }                  }
                 /* free index arrays */                  /* free index arrays */
                 for ( i = 0; i < nred; i++ ) GC_free(imat[i]);                  for ( i = 0; i < nred; i++ ) GC_free(imat[i]->index.c);
   
                 /* elimination (2nd step) */                  /* elimination (2nd step) */
                 colstat = (int *)ALLOCA(spcol*sizeof(int));                  colstat = (int *)ALLOCA(spcol*sizeof(int));
                 rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat);                  rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat);
   
                 fprintf(asir_out,"sugar=%d,rank=%d\n",sugar,rank); fflush(asir_out);                  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 */                  /* adding new bases */
                 for ( i = sprow-1; i >= rank; i-- ) GC_free(spmat[i]);                  for ( i = 0; i < rank; i++ ) {
                 for ( ; i >= 0; i-- ) {  
                         nf = vect_to_ndv(spmat[i],spcol,col,rhead,s0vect);                          nf = vect_to_ndv(spmat[i],spcol,col,rhead,s0vect);
                         SG(nf) = sugar;                          SG(nf) = sugar;
                         ndv_removecont(m,nf);                          ndv_removecont(m,nf);
Line 4006  NODE nd_f4(int m)
Line 4084  NODE nd_f4(int m)
                         g = update_base(g,nh);                          g = update_base(g,nh);
                         GC_free(spmat[i]);                          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;
Legend:
Removed from v.1.65  
changed lines
  Added in v.1.68
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>