* Added matprune executable

* Added auto-maxrelsinFF (TM) feature :)
* Modified factLat.pl accordingly


git-svn-id: https://ggnfs.svn.sourceforge.net/svnroot/ggnfs/trunk@236 80ef69ca-9341-0410-9d60-8acc7774fa31

Changelog

@@ -1,4 +1,11 @@
-07/07/21 (relf)
+07/22/06 (asl)
+  * Separated matrix pruning into separate program.
+  * factLat.pl was changed in order to track changes in flow.
+  * Added smart adjustment of maxRelsInFF: if supplied value is
+    zero, program will select it as small as possible.
+  * Set default for maxRelsInFF for factLat.pl to zero.
+
+07/21/06 (relf)
   * GGNFS_VERBOSE introduced (on by default)
   * some auxiliary output in rels.c wrapped with #ifdef GGNFS_VERBOSE .. #endif
 

Makefile

@@ -1,4 +1,4 @@
-VERSION=0.77.1-20060721
+VERSION=0.77.1-20060722
 THISDIR=branch_0
 
 TOOLSPREFIX=

include/ggnfs.h

@@ -65,7 +65,7 @@ typedef long ssize_t;
 	#endif
 
 	#ifndef M_PI
-		#define M_PI		3.14159265358979323846f
+		#define M_PI        3.14159265358979323846f
 	#endif
 
 	#ifndef M_SQRT2
@@ -276,8 +276,8 @@ typedef struct {
   u32 *relIndex; /* relIndex[i] is the index of the start of relation 'i' in
                      the following array: */
   s32 *relData;  /* This is where all the relation data is held, in a very
-                     specific format: see rels.c for the format description.
-                   */
+                    specific format: see rels.c for the format description.
+                 */
 } rel_list;
 
 /* This is a more generic structure, (which will be) used
@@ -290,13 +290,10 @@ typedef struct {
   s32  numFields;
 } llist_t;
 
-
 typedef struct {
   s32 x, y;
 } lpair_t;
 
-
-
 typedef struct {
   char       name[MAX_NAMESIZE];
   mpz_poly   f;
@@ -334,7 +331,6 @@ typedef struct {
   int        MFB_a;      /* Bits in largest value to factor into large primes */
 } nfs_fb_t;
 
-
 typedef struct {
   nfs_fb_t *FB;
   relation_t *R;
@@ -381,7 +377,6 @@ typedef struct {
   char     fbName[MAXFNAMESIZE+1];
   int      nDigits;
 } nfs_sieve_job_t;
-
 
 typedef struct {
   s32 p;
@@ -408,7 +403,6 @@ typedef struct {
   int      f;     /* This ideal has norm p^f. */
 } prime_id_t;
 
-
 typedef struct {
   mpz_poly      f;      /* NFS polynomial, f=c_0 + c_1x + ... + c_dx^d. */
   mpz_poly      T;      /* Defining monic polynomial: T(x) = c_d^{d-1}f(x/c_d). */
@@ -464,18 +458,12 @@ typedef struct {
   s32 rows[MAX_ROWS_IN_COL];
 } column_t;
 
-
-
 /* matstuff.c */
 int getDependencies(nfs_sparse_mat_t *M, llist_t *C, s32 *deps, s32 origC, long testMode);
 int writeSparseMat(char *fname, nfs_sparse_mat_t *M);
 int readSparseMat(nfs_sparse_mat_t *M, char *fname);
-int cm_init(llist_t *C, nfs_sparse_mat_t *M);
-int cm_removeCols(llist_t *C, s32 *cols, s32 numCols);
-int removeCols(nfs_sparse_mat_t *M, llist_t *C, s32 *cols, s32 m);
-int pruneMatrix(nfs_sparse_mat_t *M, s32 minExtraCols, double wtFactor, llist_t *C);
 int checkMat(nfs_sparse_mat_t *M, s32 *delCols, s32 *numDel);
-
+s32 matrixWeight(nfs_sparse_mat_t *M);
 
 /* mpz_poly.c */
 void   mpz_poly_init(mpz_poly f);
@@ -502,7 +490,7 @@ int    mpz_poly_mod_pp(mpz_poly r, mpz_poly a, mpz_poly b, mpz_t p);
 int    mpz_poly_gcd(mpz_poly g, mpz_poly h, mpz_t p);
 int    mpz_poly_modp(poly_t res, mpz_poly q, s32 p);
 int    mpz_polyCoprime(mpz_poly f, mpz_poly g, mpz_t p);
-s32   mpz_poly_evalModp(mpz_poly h, s32 p, s32 x);
+s32    mpz_poly_evalModp(mpz_poly h, s32 p, s32 x);
 void   mpz_poly_mul(mpz_poly res, mpz_poly op1, mpz_poly op2);
 void   mpz_poly_add(mpz_poly res, mpz_poly op1, mpz_poly op2);
 void   mpz_poly_mulmod(mpz_poly res, mpz_poly op1, mpz_poly op2, mpz_poly mod);
@@ -514,7 +502,6 @@ int    mpz_poly_div(mpz_poly q, mpz_poly x, mpz_poly y, mpz_t p);
 int    mpz_poly_fact(mpz_poly *Pi, int *exps, mpz_poly A, mpz_t p);
 int    mpz_poly_inv(mpz_poly inv, mpz_poly f, mpz_poly m, mpz_t p);
 
-
 int    Zalpha_sqrt(mpz_poly res, mpz_poly a, mpz_poly _f, mpz_t N, mpz_t p);
 
 /* poly.c */
@@ -711,10 +698,7 @@ int  ll_read(llist_t *C, char *fname);
 
 
 /* combparts.c */
-s32 combParts(llist_t *R, llist_t *P, int maxRelsInFF, s32 minFF);
-
-
-
+u32 combParts(llist_t *R, llist_t *P, u32 maxRelsInFF, u32 minFF, u32 minFull);
 
 /* rels.c */
 #define GETNUMRFB(_s) ((_s)>>24)

src/Makefile

@@ -12,7 +12,7 @@ BINDIR=../bin
 LIBFLAGS=$(LOCALLIB)
 
 MATBUILD_TPIE=0
-TPIE_DIR=F:/tmp/tpie/tpie_091905
+TPIE_DIR=../../tpie
 
 # Turn on\off workaround for GMP 4.2 bug
 # Consult http://swox.com/list-archives/gmp-bugs/2006-May/000475.html
@@ -46,7 +46,7 @@ OBJS=getprimes.o fbmisc.o squfof.o rels.o $(LANCZOS).o poly.o mpz_poly.o \
      matstuff.o dickman.o fbgen.o llist.o if.o rellist.o intutils.o lasieve4/mpz-ull.o
 
 BINS=$(BINDIR)/sieve $(BINDIR)/procrels $(BINDIR)/sqrt $(BINDIR)/polyselect \
-     $(BINDIR)/makefb $(BINDIR)/matsolve $(BINDIR)/matbuild
+     $(BINDIR)/makefb $(BINDIR)/matsolve $(BINDIR)/matbuild $(BINDIR)/matprune
 
 LSBINS=latsiever polsel
 
@@ -115,6 +115,9 @@ $(BINDIR)/matbuild-tpie : matbuild.c $(OBJS)
 	  experimental/combparts_tpie.cpp experimental/llist_tpie.cpp \
 	  $(OBJS) $(LIBS)
 
+$(BINDIR)/matprune : matprune.c $(OBJS)
+	$(CC) $(INC) $(CFLAGS) $(LIBFLAGS) -o $@ matprune.c $(OBJS) $(LIBS)
+
 $(BINDIR)/matsolve : matsolve.c $(OBJS)
 	$(CC) $(INC) $(CFLAGS) $(LIBFLAGS) -o $@ matsolve.c $(OBJS) $(LIBS)
 

src/combparts.c

@@ -60,7 +60,6 @@ void init_bitcount()
   }
 }
 
-
 /*************************************************************/
 void sortByNumLP(s32 *fieldsBySize, llist_t *P, char *str)
 /*************************************************************/
@@ -425,7 +424,6 @@ int merge(llist_t *R, llist_t *P, llist_t *revP, s32 P0, s32 P1, int level)
   return 0;
 }
 
-
 /* The maximum # of `s32s' we can allocate for a reverse-lookup
    table. The larger this is, the more memory makePass will use,
    but it will also be faster. When I'm less tired, I'll figure
@@ -748,26 +746,41 @@ s32 reduceRelSets(llist_t *R, llist_t *P)
   return R->index[R->numFields] - initW;
 }
 
-/*************************************************************/
-s32 removeHeavyRelSets(llist_t *R, llist_t *P, int maxRelsInRS)
-/*************************************************************/
-/* Remove any relation-set having more than maxRelsInRS      */
-/* relations.                                                */
-/* And dump a list of the number of relation-sets with each  */
-/* weight to stdout. This is to help the user if he/she needs*/
-/* to re-run with a different maxRelsInFF.                   */
-/*************************************************************/
+/**************************************************************************/
+s32 removeHeavyRelSets(llist_t *R, llist_t *P, s32 maxRelsInRS, u32 minFull)
+/**************************************************************************
+ * Remove any relation-set having more than maxRelsInRS relations.        *
+ * if maxRelsInRS == 0 do automatic adjustment of it based on minFull     *
+ * And dump a list of the number of relation-sets with each weight to     *
+ * stdout. This is to help the user if he/she needs to re-run with a      *
+ * different maxRelsInFF.                                                 *
+ *************************************************************************/
 { s32 i, *remove=NULL, numR, maxR, fSize;
-  u32 byWt[10*MAX_RELS_IN_FF], cum, cwt;
+  u32 byWt[10*MAX_RELS_IN_FF], cum=0, cwt;
 
   numR=0;
   maxR=8192;
   remove = (s32 *)lxmalloc(maxR*sizeof(s32), 1);
   for (i=0; i<10*MAX_RELS_IN_FF; i++)
     byWt[i]=0;
-  for (i=0; i<R->numFields; i++) {
+
+  /* Do one quick pass to measure weights */
+  for (i=0; i<R->numFields; ++i) {
     fSize = FIELDSIZE(R,i);
     byWt[fSize]+=1;
+  }
+
+  /* Adjust maxRelsInRS here, if automatic mode */
+  if (maxRelsInRS == 0) {
+    for (i=0, cum=0; (cum<minFull) && i<10*MAX_RELS_IN_FF; ++i) {
+      cum += byWt[i];
+    }
+
+    maxRelsInRS = i-1;
+  }
+
+  for (i=0; i<R->numFields; i++) {
+    fSize = FIELDSIZE(R,i);
     if (fSize>maxRelsInRS) {
       if (numR >= maxR) {
         maxR += 4*8192;
@@ -776,7 +789,14 @@ s32 removeHeavyRelSets(llist_t *R, llist_t *P, int maxRelsInRS)
       remove[numR++] = i;
     }
   }
-  printf("Before deleting relation sets heavier than wt %d, there were:\n", maxRelsInRS);
+
+  if (cum>0) {
+    /* Automatic mode */
+    printf("Before deleting relation sets heavier than wt %d [auto], there were:\n", maxRelsInRS);
+  } else {
+    /* User-specified mode */
+    printf("Before deleting relation sets heavier than wt %d, there were:\n", maxRelsInRS);
+  }    
   printf("Wt  |  # R-S   | Cum. R-S | Cum. wt.\n");
   printf("---------------------------------\n");
   for (i=0,cum=0,cwt=0; i<10*MAX_RELS_IN_FF; i++) {
@@ -800,13 +820,9 @@ s32 removeHeavyRelSets(llist_t *R, llist_t *P, int maxRelsInRS)
   return R->numFields;
 }
 
-
-
-
-
-/*************************************************************/
-s32 combParts(llist_t *R, llist_t *P, int maxRelsInFF, s32 minFF)
-/*************************************************************/
+/****************************************************************************/
+u32 combParts(llist_t *R, llist_t *P, u32 maxRelsInFF, u32 minFF, u32 minFull)
+/****************************************************************************/
 /*Input:
      'P' is a list containing the large primes. It should not
      actually contain the primes, but rather a unique integer associated
@@ -836,7 +852,7 @@ s32 combParts(llist_t *R, llist_t *P, int maxRelsInFF, s32 minFF)
 */
 /**************************************************************/ 
 { s32 i, wt0, wt1;
-  s32 lastFull, full;
+  u32 lastFull, full;
   int  pass=0;
   double shrink;
 
@@ -863,14 +879,16 @@ s32 combParts(llist_t *R, llist_t *P, int maxRelsInFF, s32 minFF)
     ll_appendField(R, &i, 1);
   }
 
-  full = -1; /* Force at least two passes. */
+  /* Force at least two passes. */
+
+  full = 0;  
   do {
     printf("  pass %d...\n", ++pass);
     lastFull = full;
     full = makePass(R, P);
     checkR(R);
     printf("* There are now %" PRId32 " full relations.\n", full);
-  } while (lastFull < full);
+  } while ((lastFull < full) || pass < 1);
 
   /* Drop any relation-sets still containing a large prime: */
   keepFulls(R, P); 
@@ -894,7 +912,7 @@ s32 combParts(llist_t *R, llist_t *P, int maxRelsInFF, s32 minFF)
     shrink = (double)(wt0-wt1)/wt0;
   }  while (shrink > 0.15);
 #endif
-  full = removeHeavyRelSets(R, P, maxRelsInFF);
+  full = removeHeavyRelSets(R, P, maxRelsInFF, minFull);
   msgLog("", "After removing heavy rel-sets, weight is %" PRId32 ".", R->index[R->numFields]);
   printf("After removing heavy rel-sets, weight is %" PRId32 ".\n", R->index[R->numFields]);
   if (ll_verify(R)) {

src/experimental/combparts_tpie.cpp

@@ -412,21 +412,39 @@ int combinedSize(LSList& R, s32 a, s32 b) {
  * weight to stdout. This is to help the user if he/she needs
  * to re-run with a different maxRelsInFF.
  **************************************************************/
-s32 removeHeavyRelSets(LSList& R, LSList& P, int maxRelsInRS) {
+u32 removeHeavyRelSets(LSList& R, LSList& P, u32 maxRelsInRS, u32 minFull) {
 	s32 *remove=NULL, numR, maxR, fSize;
 	u32 byWt[10*MAX_RELS_IN_FF];
+  u32 i, cum, cwt;
 
-	numR=0;
-	maxR=8192;
+  cum = 0;
+	numR = 0;
+	maxR = 8192;
+
 	if (!(remove = (s32 *)malloc(maxR*sizeof(s32)))) {
 		printf("removeHeavyRelSets(): Memory allocation error for remove!\n");
 		exit(-1);
 	}
-	for (u32 i=0; i<10*MAX_RELS_IN_FF; i++)
+	for (i=0; i<10*MAX_RELS_IN_FF; i++)
 		byWt[i]=0;
-	for (u32 i=0; i<R.numFields(); i++) {
+
+  /* Do one quick pass to measure weights */
+  for (i=0; i<R.numFields(); ++i) {
+    fSize = R.fieldSize(i);
+    byWt[fSize]+=1;
+  }
+
+  /* Adjust maxRelsInRS here, if automatic mode */
+  if (maxRelsInRS == 0) {
+    for (i=0, cum=0; (cum<minFull) && i<10*MAX_RELS_IN_FF; ++i ) {
+      cum += byWt[i];
+    }
+
+    maxRelsInRS = i;
+  }
+
+	for (i=0; i<R.numFields(); i++) {
 		fSize = R.fieldSize(i);
-		byWt[fSize]+=1;
 		if (fSize>maxRelsInRS) {
 			if (numR >= maxR) {
 				maxR += 8192;
@@ -439,10 +457,18 @@ s32 removeHeavyRelSets(LSList& R, LSList& P, int maxRelsInRS) {
 			remove[numR++] = i;
 		}
 	}
-	printf("Before deleting relation sets heavier than wt %d, there were:\n", maxRelsInRS);
+
+  if (cum>0) {
+    /* Automatic mode */
+    printf("Before deleting relation sets heavier than wt %d [auto], there were:\n", maxRelsInRS);
+  } else {
+    /* User-specified mode */
+    printf("Before deleting relation sets heavier than wt %d, there were:\n", maxRelsInRS);
+  }    
+
 	printf("Wt	|	# R-S	 | Cum. R-S | Cum. wt.\n");
 	printf("---------------------------------\n");
-	for (u32 i=0,cum=0,cwt=0; i<10*MAX_RELS_IN_FF; i++) {
+	for (i=0,cum=0,cwt=0; i<10*MAX_RELS_IN_FF; i++) {
 		cum += byWt[i];
 		cwt += i*byWt[i];
 		if (byWt[i]>0) {
@@ -736,9 +762,9 @@ s32 reduceRelSets(LSList& R, LSList& P) {
  * 	'P' will be modified in the process.
  * 	Return value: The number of full relations built.
  **************************************************************/ 
-s32 combParts_new(LSList& P, LSList& R, int maxRelsInFF, s32 minFF) {
+u32 combParts_new(LSList& P, LSList& R, u32 maxRelsInFF, u32 minFF, u32 minFull) {
 	s32 wt0, wt1;
-	s32 lastFull, full;
+	u32 lastFull, full;
 	int	pass=0;
 	double shrink;
 
@@ -783,7 +809,7 @@ s32 combParts_new(LSList& P, LSList& R, int maxRelsInFF, s32 minFF) {
 				wt0, wt1);
 		shrink = (double)(wt0-wt1)/wt0;
 	} while (shrink > 0.10);
-	full = removeHeavyRelSets(R, P, maxRelsInFF);
+	full = removeHeavyRelSets(R, P, maxRelsInFF, minFull);
 	msgLog("", "After removing heavy rel-sets, weight is %" PRId32 ".", R.Weight());
 	//printf("After removing heavy rel-sets, weight is %" PRId32 ".\n", R.Weight());
 
@@ -814,15 +840,15 @@ bool checkRelSets(LSList&R, LSList& dataList) {
 	return true;
 }
 
-s32 combParts_tpie(llist_t *lR, llist_t *lP, int maxRelsInFF, s32 minFF) {
+u32 combParts_tpie(llist_t *lR, llist_t *lP, u32 maxRelsInFF, u32 minFF, u32 minFull) {
 	MM_manager.ignore_memory_limit();
 
 	LSList P(lP);
 	LSList R(lR);
 	LSList cmpP(lP);
 
 	double t1=sTime();
-	s32 numFulls=combParts_new(P,R,maxRelsInFF,minFF);
+	s32 numFulls=combParts_new(P, R ,maxRelsInFF, minFF, minFull);
 	if (checkRelSets(R,cmpP)) { 
 		printf("checkRelSets() reports, that R has good chances to be fine ;)\n");
 	} else {