[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.75 and 1.76

version 1.75, 2003/09/19 10:09:42 version 1.76, 2003/09/28 09:18:57
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.74 2003/09/19 02:33:12 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.75 2003/09/19 10:09:42 noro Exp $ */
   
 #include "ca.h"  #include "ca.h"
 #include "parse.h"  #include "parse.h"
Line 100  typedef struct oBaseSet {
Line 100  typedef struct oBaseSet {
 typedef struct oNM_ind_pair  typedef struct oNM_ind_pair
 {  {
         NM mul;          NM mul;
         int index;          int index,sugar;
 } *NM_ind_pair;  } *NM_ind_pair;
   
 typedef struct oIndArray  typedef struct oIndArray
Line 213  if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEX
Line 213  if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEX
 if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}  if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
 #define NEXTND_pairs(r,c) \  #define NEXTND_pairs(r,c) \
 if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);}  if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);}
 #define MKNM_ind_pair(p,m,i) (NEWNM_ind_pair(p),(p)->mul=(m),(p)->index=(i))  #define MKNM_ind_pair(p,m,i,s) (NEWNM_ind_pair(p),(p)->mul=(m),(p)->index=(i),(p)->sugar = (s))
   
 /* deallocators */  /* deallocators */
 #define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m)  #define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m)
Line 368  int nm_ind_pair_to_vect(int m,UINT *s0,int n,NM_ind_pa
Line 368  int nm_ind_pair_to_vect(int m,UINT *s0,int n,NM_ind_pa
 IndArray nm_ind_pair_to_vect_compress(int m,UINT *s0,int n,NM_ind_pair pair);  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);
   
   /* elimination */
   int nd_gauss_elim_mod(int **mat0,int *sugar,int row,int col,int md,int *colstat);
   int nd_gauss_elim_sf(int **mat0,int *sugar,int row,int col,int md,int *colstat);
   
 void nd_free_private_storage()  void nd_free_private_storage()
 {  {
         _nm_free_list = 0;          _nm_free_list = 0;
Line 3777  IndArray nm_ind_pair_to_vect_compress(int mod,UINT *s0
Line 3781  IndArray nm_ind_pair_to_vect_compress(int mod,UINT *s0
 }  }
   
   
 void ndv_reduce_vect(int m,UINT *svect,int col,IndArray *imat,NM_ind_pair *rp0,int nred)  int ndv_reduce_vect(int m,UINT *svect,int col,IndArray *imat,NM_ind_pair *rp0,int nred)
 {  {
         int i,j,k,len,pos,prev;          int i,j,k,len,pos,prev;
         UINT c,c1,c2,c3,up,lo,dmy;          UINT c,c1,c2,c3,up,lo,dmy;
Line 3788  void ndv_reduce_vect(int m,UINT *svect,int col,IndArra
Line 3792  void ndv_reduce_vect(int m,UINT *svect,int col,IndArra
         NDV redv;          NDV redv;
         NMV mr;          NMV mr;
         NODE rp;          NODE rp;
           int maxrs;
   
           maxrs = 0;
         for ( i = 0; i < nred; i++ ) {          for ( i = 0; i < nred; i++ ) {
                 ivect = imat[i];                  ivect = imat[i];
                 k = ivect->head; svect[k] %= m;                  k = ivect->head; svect[k] %= m;
                 if ( c = svect[k] ) {                  if ( c = svect[k] ) {
                           maxrs = MAX(maxrs,rp0[i]->sugar);
                         c = m-c; redv = nd_ps[rp0[i]->index];                          c = m-c; redv = nd_ps[rp0[i]->index];
                         len = LEN(redv); mr = BDY(redv);                          len = LEN(redv); mr = BDY(redv);
                         svect[k] = 0; prev = k;                          svect[k] = 0; prev = k;
Line 3832  void ndv_reduce_vect(int m,UINT *svect,int col,IndArra
Line 3839  void ndv_reduce_vect(int m,UINT *svect,int col,IndArra
         }          }
         for ( i = 0; i < col; i++ )          for ( i = 0; i < col; i++ )
                 if ( svect[i] >= (UINT)m ) svect[i] %= m;                  if ( svect[i] >= (UINT)m ) svect[i] %= m;
           return maxrs;
 }  }
   
 void ndv_reduce_vect_sf(int m,UINT *svect,int col,IndArray *imat,NM_ind_pair *rp0,int nred)  int 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;          int i,j,k,len,pos,prev;
         UINT c,c1,c2,c3,up,lo,dmy;          UINT c,c1,c2,c3,up,lo,dmy;
Line 3845  void ndv_reduce_vect_sf(int m,UINT *svect,int col,IndA
Line 3853  void ndv_reduce_vect_sf(int m,UINT *svect,int col,IndA
         NDV redv;          NDV redv;
         NMV mr;          NMV mr;
         NODE rp;          NODE rp;
           int maxrs;
   
           maxrs = 0;
         for ( i = 0; i < nred; i++ ) {          for ( i = 0; i < nred; i++ ) {
                 ivect = imat[i];                  ivect = imat[i];
                 k = ivect->head; svect[k] %= m;                  k = ivect->head; svect[k] %= m;
                 if ( c = svect[k] ) {                  if ( c = svect[k] ) {
                           maxrs = MAX(maxrs,rp0[i]->sugar);
                         c = _chsgnsf(c); redv = nd_ps[rp0[i]->index];                          c = _chsgnsf(c); redv = nd_ps[rp0[i]->index];
                         len = LEN(redv); mr = BDY(redv);                          len = LEN(redv); mr = BDY(redv);
                         svect[k] = 0; prev = k;                          svect[k] = 0; prev = k;
Line 3878  void ndv_reduce_vect_sf(int m,UINT *svect,int col,IndA
Line 3889  void ndv_reduce_vect_sf(int m,UINT *svect,int col,IndA
                         }                          }
                 }                  }
         }          }
           return maxrs;
 }  }
   
 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 3927  int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vec
Line 3939  int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vec
 {  {
         NODE rp0,rp;          NODE rp0,rp;
         NM mul,head,s0,s;          NM mul,head,s0,s;
         int index,col,i;          int index,col,i,sugar;
         RHist h;          RHist h;
         UINT *s0v,*p;          UINT *s0v,*p;
         NM_ind_pair pair;          NM_ind_pair pair;
Line 3946  int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vec
Line 3958  int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vec
                         NEWNM(mul);                          NEWNM(mul);
                         ndl_sub(DL(head),DL(h),DL(mul));                          ndl_sub(DL(head),DL(h),DL(mul));
                         if ( ndl_check_bound2(index,DL(mul)) ) return 0;                          if ( ndl_check_bound2(index,DL(mul)) ) return 0;
                         MKNM_ind_pair(pair,mul,index);                          sugar = TD(DL(mul))+SG(nd_ps[index]);
                           MKNM_ind_pair(pair,mul,index,sugar);
                         red = ndv_mul_nm_symbolic(mul,nd_ps[index]);                          red = ndv_mul_nm_symbolic(mul,nd_ps[index]);
                         add_pbucket_symbolic(bucket,nd_remove_head(red));                          add_pbucket_symbolic(bucket,nd_remove_head(red));
                         NEXTNODE(rp0,rp); BDY(rp) = (pointer)pair;                          NEXTNODE(rp0,rp); BDY(rp) = (pointer)pair;
Line 4023  NODE nd_f4(int m)
Line 4036  NODE nd_f4(int m)
                 /* adding new bases */                  /* adding new bases */
                 for ( r = nflist; r; r = NEXT(r) ) {                  for ( r = nflist; r; r = NEXT(r) ) {
                         nf = (NDV)BDY(r);                          nf = (NDV)BDY(r);
                         SG(nf) = sugar;  
                         ndv_removecont(m,nf);                          ndv_removecont(m,nf);
                         nh = ndv_newps(nf,0);                          nh = ndv_newps(nf,0);
                         d = update_pairs(d,g,nh);                          d = update_pairs(d,g,nh);
Line 4049  NODE nd_f4_red(int m,ND_pairs sp0,UINT *s0vect,int col
Line 4061  NODE nd_f4_red(int m,ND_pairs sp0,UINT *s0vect,int col
         int *colstat;          int *colstat;
         struct oEGT eg0,eg1,eg_f4;          struct oEGT eg0,eg1,eg_f4;
         NM_ind_pair *rvect;          NM_ind_pair *rvect;
           int maxrs;
           int *spsugar;
   
         get_eg(&eg0);          get_eg(&eg0);
         for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ );          for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ );
Line 4068  NODE nd_f4_red(int m,ND_pairs sp0,UINT *s0vect,int col
Line 4082  NODE nd_f4_red(int m,ND_pairs sp0,UINT *s0vect,int col
         /* elimination (1st step) */          /* elimination (1st step) */
         spmat = (int **)ALLOCA(nsp*sizeof(UINT *));          spmat = (int **)ALLOCA(nsp*sizeof(UINT *));
         svect = (UINT *)ALLOCA(col*sizeof(UINT));          svect = (UINT *)ALLOCA(col*sizeof(UINT));
           spsugar = (int *)ALLOCA(nsp*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_sp(m,0,sp,&spol);                  nd_sp(m,0,sp,&spol);
                 if ( !spol ) continue;                  if ( !spol ) continue;
                 nd_to_vect(m,s0vect,col,spol,svect);                  nd_to_vect(m,s0vect,col,spol,svect);
                 nd_free(spol);                  if ( m == -1 )
                 if ( m == -1 ) ndv_reduce_vect_sf(m,svect,col,imat,rvect,nred);                          maxrs = ndv_reduce_vect_sf(m,svect,col,imat,rvect,nred);
                 else ndv_reduce_vect(m,svect,col,imat,rvect,nred);                  else
                           maxrs = ndv_reduce_vect(m,svect,col,imat,rvect,nred);
                 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));
                         for ( j = k = 0; j < col; j++ )                          for ( j = k = 0; j < col; j++ )
                                 if ( !rhead[j] ) v[k++] = svect[j];                                  if ( !rhead[j] ) v[k++] = svect[j];
                           spsugar[sprow] = MAX(maxrs,SG(spol));
                         sprow++;                          sprow++;
                 }                  }
                   nd_free(spol);
         }          }
         /* free index arrays */          /* free index arrays */
         for ( i = 0; i < nred; i++ ) GC_free(imat[i]->index.c);          for ( i = 0; i < nred; i++ ) GC_free(imat[i]->index.c);
Line 4089  NODE nd_f4_red(int m,ND_pairs sp0,UINT *s0vect,int col
Line 4107  NODE nd_f4_red(int m,ND_pairs sp0,UINT *s0vect,int col
         /* elimination (2nd step) */          /* elimination (2nd step) */
         colstat = (int *)ALLOCA(spcol*sizeof(int));          colstat = (int *)ALLOCA(spcol*sizeof(int));
         if ( m == -1 )          if ( m == -1 )
                 rank = generic_gauss_elim_sf(spmat,sprow,spcol,m,colstat);                  rank = nd_gauss_elim_sf(spmat,spsugar,sprow,spcol,m,colstat);
         else          else
                 rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat);                  rank = nd_gauss_elim_mod(spmat,spsugar,sprow,spcol,m,colstat);
         r0 = 0;          r0 = 0;
         for ( i = 0; i < rank; i++ ) {          for ( i = 0; i < rank; i++ ) {
                 NEXTNODE(r0,r); BDY(r) =                  NEXTNODE(r0,r); BDY(r) =
                         (pointer)vect_to_ndv(spmat[i],spcol,col,rhead,s0vect);                          (pointer)vect_to_ndv(spmat[i],spcol,col,rhead,s0vect);
                   SG((NDV)BDY(r)) = spsugar[i];
                 GC_free(spmat[i]);                  GC_free(spmat[i]);
         }          }
         for ( ; i < sprow; i++ ) GC_free(spmat[i]);          for ( ; i < sprow; i++ ) GC_free(spmat[i]);
Line 4269  void nd_exec_f4_red_dist()
Line 4288  void nd_exec_f4_red_dist()
         struct order_spec ord;          struct order_spec ord;
         Obj ordspec;          Obj ordspec;
         ND spol;          ND spol;
           int maxrs;
           int *spsugar;
   
         nd_read = iofp[0].in;          nd_read = iofp[0].in;
         nd_write = iofp[0].out;          nd_write = iofp[0].out;
Line 4327  void nd_exec_f4_red_dist()
Line 4348  void nd_exec_f4_red_dist()
         /* elimination (1st step) */          /* elimination (1st step) */
         spmat = (int **)MALLOC(nsp*sizeof(UINT *));          spmat = (int **)MALLOC(nsp*sizeof(UINT *));
         svect = (UINT *)MALLOC(col*sizeof(UINT));          svect = (UINT *)MALLOC(col*sizeof(UINT));
           spsugar = (int *)ALLOCA(nsp*sizeof(UINT));
         for ( a = sprow = 0; a < nsp; a++ ) {          for ( a = sprow = 0; a < nsp; a++ ) {
                 nd_sp(m,0,sp0[a],&spol);                  nd_sp(m,0,sp0[a],&spol);
                 if ( !spol ) continue;                  if ( !spol ) continue;
                 nd_to_vect(m,s0vect,col,spol,svect);                  nd_to_vect(m,s0vect,col,spol,svect);
                 nd_free(spol);                  if ( m == -1 )
                 if ( m == -1 ) ndv_reduce_vect_sf(m,svect,col,imat,rp0,nred);                          maxrs = ndv_reduce_vect_sf(m,svect,col,imat,rp0,nred);
                 else ndv_reduce_vect(m,svect,col,imat,rp0,nred);                  else
                           maxrs = ndv_reduce_vect(m,svect,col,imat,rp0,nred);
                 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(spcol*sizeof(UINT));                          spmat[sprow] = v = (UINT *)MALLOC(spcol*sizeof(UINT));
                         for ( j = k = 0; j < col; j++ )                          for ( j = k = 0; j < col; j++ )
                                 if ( !rhead[j] ) v[k++] = svect[j];                                  if ( !rhead[j] ) v[k++] = svect[j];
                           spsugar[sprow] = MAX(maxrs,SG(spol));
                         sprow++;                          sprow++;
                 }                  }
                   nd_free(spol);
         }          }
         /* elimination (2nd step) */          /* elimination (2nd step) */
         colstat = (int *)ALLOCA(spcol*sizeof(int));          colstat = (int *)ALLOCA(spcol*sizeof(int));
         if ( m == -1 )          if ( m == -1 )
                 rank = generic_gauss_elim_sf(spmat,sprow,spcol,m,colstat);                  rank = nd_gauss_elim_sf(spmat,spsugar,sprow,spcol,m,colstat);
         else          else
                 rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat);                  rank = nd_gauss_elim_mod(spmat,spsugar,sprow,spcol,m,colstat);
         nd_send_int(rank);          nd_send_int(rank);
         for ( i = 0; i < rank; i++ ) {          for ( i = 0; i < rank; i++ ) {
                 nf = vect_to_ndv(spmat[i],spcol,col,rhead,s0vect);                  nf = vect_to_ndv(spmat[i],spcol,col,rhead,s0vect);
                 nd_send_ndv(nf);                  nd_send_ndv(nf);
         }          }
         fflush(nd_write);          fflush(nd_write);
   }
   
   int nd_gauss_elim_mod(int **mat0,int *sugar,int row,int col,int md,int *colstat)
   {
           int i,j,k,l,inv,a,rank,s;
           unsigned int *t,*pivot,*pk;
           unsigned int **mat;
   
           mat = (unsigned int **)mat0;
           for ( rank = 0, j = 0; j < col; j++ ) {
                   for ( i = rank; i < row; i++ )
                           mat[i][j] %= md;
                   for ( i = rank; i < row; i++ )
                           if ( mat[i][j] )
                                   break;
                   if ( i == row ) {
                           colstat[j] = 0;
                           continue;
                   } else
                           colstat[j] = 1;
                   if ( i != rank ) {
                           t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                           s = sugar[i]; sugar[i] = sugar[rank]; sugar[rank] = s;
                   }
                   pivot = mat[rank];
                   s = sugar[rank];
                   inv = invm(pivot[j],md);
                   for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                           if ( *pk ) {
                                   if ( *pk >= (unsigned int)md )
                                           *pk %= md;
                                   DMAR(*pk,inv,0,md,*pk)
                           }
                   for ( i = rank+1; i < row; i++ ) {
                           t = mat[i];
                           if ( a = t[j] ) {
                                   sugar[i] = MAX(sugar[i],s);
                                   red_by_vect(md,t+j,pivot+j,md-a,col-j);
                           }
                   }
                   rank++;
           }
           for ( j = col-1, l = rank-1; j >= 0; j-- )
                   if ( colstat[j] ) {
                           pivot = mat[l];
                           s = sugar[l];
                           for ( i = 0; i < l; i++ ) {
                                   t = mat[i];
                                   t[j] %= md;
                                   if ( a = t[j] ) {
                                           sugar[i] = MAX(sugar[i],s);
                                           red_by_vect(md,t+j,pivot+j,md-a,col-j);
                                   }
                           }
                           l--;
                   }
           for ( j = 0, l = 0; l < rank; j++ )
                   if ( colstat[j] ) {
                           t = mat[l];
                           for ( k = j; k < col; k++ )
                                   if ( t[k] >= (unsigned int)md )
                                           t[k] %= md;
                           l++;
                   }
           return rank;
   }
   
   int nd_gauss_elim_sf(int **mat0,int *sugar,int row,int col,int md,int *colstat)
   {
           int i,j,k,l,inv,a,rank,s;
           unsigned int *t,*pivot,*pk;
           unsigned int **mat;
   
           mat = (unsigned int **)mat0;
           for ( rank = 0, j = 0; j < col; j++ ) {
                   for ( i = rank; i < row; i++ )
                           if ( mat[i][j] )
                                   break;
                   if ( i == row ) {
                           colstat[j] = 0;
                           continue;
                   } else
                           colstat[j] = 1;
                   if ( i != rank ) {
                           t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                           s = sugar[i]; sugar[i] = sugar[rank]; sugar[rank] = s;
                   }
                   pivot = mat[rank];
                   s = sugar[rank];
                   inv = _invsf(pivot[j]);
                   for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                           if ( *pk )
                                   *pk = _mulsf(*pk,inv);
                   for ( i = rank+1; i < row; i++ ) {
                           t = mat[i];
                           if ( a = t[j] ) {
                                   sugar[i] = MAX(sugar[i],s);
                                   red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
                           }
                   }
                   rank++;
           }
           for ( j = col-1, l = rank-1; j >= 0; j-- )
                   if ( colstat[j] ) {
                           pivot = mat[l];
                           s = sugar[l];
                           for ( i = 0; i < l; i++ ) {
                                   t = mat[i];
                                   if ( a = t[j] ) {
                                           sugar[i] = MAX(sugar[i],s);
                                           red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
                                   }
                           }
                           l--;
                   }
           return rank;
 }  }
Legend:
Removed from v.1.75  
changed lines
  Added in v.1.76
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>