This is the mail archive of the gcc-patches@gcc.gnu.org mailing list for the GCC project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[autovect] Add entry point for omega and polyhedral data dep analyzers


The new data dependence analyzers (omega and polyhedral) are not yet
finished, and thus they are completely disabled by default.

Bootstrapped on amd64-linux.

Index: ChangeLog.autovect
===================================================================
RCS file: /cvs/gcc/gcc/gcc/Attic/ChangeLog.autovect,v
retrieving revision 1.1.2.59
diff -d -u -p -r1.1.2.59 ChangeLog.autovect
--- ChangeLog.autovect	20 Apr 2005 17:34:39 -0000	1.1.2.59
+++ ChangeLog.autovect	21 Apr 2005 15:05:49 -0000
@@ -1,3 +1,39 @@
+2005-04-21  Sebastian Pop  <pop@cri.ensmp.fr>
+
+	* Makefile.in (POLYHEDRON_H): New variable.
+	(TREE_DATA_REF_H): Depends on POLYHEDRON_H.
+	(OBJS-common): Depends on polyhedron.o.
+	(polyhedron.o): New rule.
+	* lambda-code.c (gcd): Moved to lambda.h.
+	(gcd_vector): Renamed lambda_vector_gcd, and moved to lambda.h.
+	(lambda_compute_target_space): Rename uses of gcd_vector.
+	* lambda-mat.c (lambda_matrix_row_exchange): Moved to lambda.h.
+	* lambda.h (lambda_matrix_row_exchange): Make it static inline.
+	(lambda_vector_div_const, lambda_vector_gcd, lambda_vector_normalize,
+	lambda_vector_scalar_product, lambda_vector_linear_combine_1, 
+	lambda_matrix_extend, lambda_matrix_exchange_columns, 
+	lambda_matrix_copy_columns, lambda_matrix_copy_c2v,
+	lambda_matrix_copy_v2c, lambda_vector_exchange): New functions.
+	* omega.c (wild_name): Allocate more elements.
+	(omega_eqn_is_red): Check bounds of printed variable, avoid the use
+	of omega_current_get_var_name.
+	(omega_delete_geq, omega_delete_geq_extra): Make static.
+	(omega_initialize): Initialize current variables names for
+	pretty printing.
+	* omega.h (omega_current_get_var_name): Remove declaration.
+	(omega_delete_geq, omega_delete_geq_extra): Moved to omega.c.
+	* polyhedron.c, polyhedron.h: New files.
+	* tree-data-ref.c (gcd): Moved to lambda.h.
+	(init_csys_eq_with_af, find_loops_surrounding, init_csys_for_ddr,
+	polyhedra_dependence_tester, omega_compute_classic_representations,
+	omega_dependence_tester): New.
+	(compute_affine_dependence): Insert entry point to the new 
+	dependence analyzers, but they are disabled for now.
+	* tree-data-ref.h: Depend on omega.h and polyhedron.h.
+	(data_dependence_relation): Add two fields:
+	dependence_constraint_system, and dependence_polyhedron.
+	(DDR_CSYS, DDR_POLYHEDRON): New.
+
 2005-04-20  Sebastian Pop  <pop@cri.ensmp.fr>
 
 	* omega.h, omega.c: Write some documentation for functions, 
Index: Makefile.in
===================================================================
RCS file: /cvs/gcc/gcc/gcc/Makefile.in,v
retrieving revision 1.1419.2.9
diff -d -u -p -r1.1419.2.9 Makefile.in
--- Makefile.in	20 Apr 2005 11:12:54 -0000	1.1419.2.9
+++ Makefile.in	21 Apr 2005 15:05:50 -0000
@@ -763,7 +763,8 @@ DIAGNOSTIC_H = diagnostic.h diagnostic.d
 C_PRETTY_PRINT_H = c-pretty-print.h $(PRETTY_PRINT_H) $(C_COMMON_H) $(TREE_H)
 SCEV_H = tree-scalar-evolution.h $(GGC_H) tree-chrec.h
 LAMBDA_H = lambda.h tree.h vec.h $(GGC_H)
-TREE_DATA_REF_H = tree-data-ref.h $(LAMBDA_H)
+POLYHEDRON_H = polyhedron.h $(LAMBDA_H)
+TREE_DATA_REF_H = tree-data-ref.h $(LAMBDA_H) $(POLYHEDRON_H)
 
 #
 # Now figure out from those variables how to compile and link.
@@ -921,7 +922,7 @@ C_OBJS = c-lang.o stub-objc.o $(C_AND_OB
 
 # Language-independent object files.
 OBJS-common = \
- omega.o \
+ omega.o polyhedron.o \
  tree-chrec.o tree-scalar-evolution.o tree-data-ref.o tree-elim-check.o	   \
  tree-cfg.o tree-dfa.o tree-eh.o tree-ssa.o tree-optimize.o tree-gimple.o  \
  gimplify.o tree-pretty-print.o tree-into-ssa.o          \
@@ -1786,6 +1787,7 @@ tree-browser.o : tree-browser.c tree-bro
    $(TM_H) coretypes.h
 omega.o : omega.c omega.h $(CONFIG_H) $(SYSTEM_H) coretypes.h $(TM_H) \
    errors.h $(GGC_H) $(TREE_H) diagnostic.h varray.h tree-pass.h 
+polyhedron.o : polyhedron.c $(POLYHEDRON_H) $(CONFIG_H) $(SYSTEM_H)
 tree-chrec.o: tree-chrec.c $(CONFIG_H) $(SYSTEM_H) coretypes.h $(TM_H) \
    errors.h $(GGC_H) $(TREE_H) tree-chrec.h tree-pass.h cfgloop.h $(TREE_FLOW.h)
 tree-scalar-evolution.o: tree-scalar-evolution.c $(CONFIG_H) $(SYSTEM_H) \
Index: lambda-code.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/lambda-code.c,v
retrieving revision 2.18.2.3
diff -d -u -p -r2.18.2.3 lambda-code.c
--- lambda-code.c	10 Apr 2005 21:09:26 -0000	2.18.2.3
+++ lambda-code.c	21 Apr 2005 15:05:50 -0000
@@ -441,45 +441,6 @@ lambda_lattice_compute_base (lambda_loop
   return ret;
 }
 
-/* Compute the greatest common denominator of two numbers (A and B) using
-   Euclid's algorithm.  */
-
-static int
-gcd (int a, int b)
-{
-
-  int x, y, z;
-
-  x = abs (a);
-  y = abs (b);
-
-  while (x > 0)
-    {
-      z = y % x;
-      y = x;
-      x = z;
-    }
-
-  return (y);
-}
-
-/* Compute the greatest common denominator of a VECTOR of SIZE numbers.  */
-
-static int
-gcd_vector (lambda_vector vector, int size)
-{
-  int i;
-  int gcd1 = 0;
-
-  if (size > 0)
-    {
-      gcd1 = vector[0];
-      for (i = 1; i < size; i++)
-	gcd1 = gcd (gcd1, vector[i]);
-    }
-  return gcd1;
-}
-
 /* Compute the least common multiple of two numbers A and B .  */
 
 static int
@@ -848,7 +809,7 @@ lambda_compute_target_space (lambda_loop
       LN_LOOPS (target_nest)[i] = target_loop;
 
       /* Computes the gcd of the coefficients of the linear part.  */
-      gcd1 = gcd_vector (target[i], i);
+      gcd1 = lambda_vector_gcd (target[i], i);
 
       /* Include the denominator in the GCD.  */
       gcd1 = gcd (gcd1, determinant);
@@ -911,8 +872,8 @@ lambda_compute_target_space (lambda_loop
 	    }
 	  /* Find the gcd and divide by it here, rather than doing it
 	     at the tree level.  */
-	  gcd1 = gcd_vector (LLE_COEFFICIENTS (target_expr), depth);
-	  gcd2 = gcd_vector (LLE_INVARIANT_COEFFICIENTS (target_expr),
+	  gcd1 = lambda_vector_gcd (LLE_COEFFICIENTS (target_expr), depth);
+	  gcd2 = lambda_vector_gcd (LLE_INVARIANT_COEFFICIENTS (target_expr),
 			     invariants);
 	  gcd1 = gcd (gcd1, gcd2);
 	  gcd1 = gcd (gcd1, LLE_CONSTANT (target_expr));
@@ -967,8 +928,8 @@ lambda_compute_target_space (lambda_loop
 	    }
 	  /* Find the gcd and divide by it here, instead of at the
 	     tree level.  */
-	  gcd1 = gcd_vector (LLE_COEFFICIENTS (target_expr), depth);
-	  gcd2 = gcd_vector (LLE_INVARIANT_COEFFICIENTS (target_expr),
+	  gcd1 = lambda_vector_gcd (LLE_COEFFICIENTS (target_expr), depth);
+	  gcd2 = lambda_vector_gcd (LLE_INVARIANT_COEFFICIENTS (target_expr),
 			     invariants);
 	  gcd1 = gcd (gcd1, gcd2);
 	  gcd1 = gcd (gcd1, LLE_CONSTANT (target_expr));
Index: lambda-mat.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/lambda-mat.c,v
retrieving revision 2.8
diff -d -u -p -r2.8 lambda-mat.c
--- lambda-mat.c	25 Sep 2004 14:36:37 -0000	2.8
+++ lambda-mat.c	21 Apr 2005 15:05:50 -0000
@@ -194,18 +194,6 @@ lambda_matrix_delete_rows (lambda_matrix
     mat[i] = NULL;
 }
 
-/* Swap rows R1 and R2 in matrix MAT.  */
-
-void
-lambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2)
-{
-  lambda_vector row;
-
-  row = mat[r1];
-  mat[r1] = mat[r2];
-  mat[r2] = row;
-}
-
 /* Add a multiple of row R1 of matrix MAT with N columns to row R2:
    R2 = R2 + CONST1 * R1.  */
 
Index: lambda.h
===================================================================
RCS file: /cvs/gcc/gcc/gcc/lambda.h,v
retrieving revision 2.8
diff -d -u -p -r2.8 lambda.h
--- lambda.h	3 Nov 2004 17:32:34 -0000	2.8
+++ lambda.h	21 Apr 2005 15:05:50 -0000
@@ -161,7 +161,6 @@ void lambda_matrix_add_mc (lambda_matrix
 void lambda_matrix_mult (lambda_matrix, lambda_matrix, lambda_matrix,
 			 int, int, int);
 void lambda_matrix_delete_rows (lambda_matrix, int, int, int);
-void lambda_matrix_row_exchange (lambda_matrix, int, int);
 void lambda_matrix_row_add (lambda_matrix, int, int, int, int);
 void lambda_matrix_row_negate (lambda_matrix mat, int, int);
 void lambda_matrix_row_mc (lambda_matrix, int, int, int);
@@ -227,7 +226,18 @@ lambda_vector_new (int size)
   return ggc_alloc_cleared (size * sizeof(int));
 }
 
+/* Divide vector VEC1 of length SIZE by a constant CONST1,
+   and store the result in VEC2.  */
+ 
+static inline void
+lambda_vector_div_const (lambda_vector vec1, lambda_vector vec2,
+			 int size, int const1)
+{
+  int i;
 
+  for (i = 0; i < size; i++)
+    vec2[i] = vec1[i] / const1;
+}
 
 /* Multiply vector VEC1 of length SIZE by a constant CONST1,
    and store the result in VEC2.  */
@@ -378,8 +388,192 @@ print_lambda_vector (FILE * outfile, lam
   int i;
 
   for (i = 0; i < n; i++)
-    fprintf (outfile, "%3d ", vector[i]);
+    fprintf (outfile, "%4d ", vector[i]);
   fprintf (outfile, "\n");
 }
+
+/* Compute the greatest common divisor of two numbers using
+   Euclid's algorithm.  */
+
+static inline int 
+gcd (int a, int b)
+{
+  int x, y, z;
+
+  x = abs (a);
+  y = abs (b);
+
+  while (x > 0)
+    {
+      z = y % x;
+      y = x;
+      x = z;
+    }
+
+  return y;
+}
+
+/* Compute the greatest common divisor of a VECTOR of SIZE numbers.  */
+
+static inline int
+lambda_vector_gcd (lambda_vector vector, int size)
+{
+  int i;
+  int gcd1 = 0;
+
+  if (size > 0)
+    {
+      gcd1 = vector[0];
+      for (i = 1; i < size; i++)
+	gcd1 = gcd (gcd1, vector[i]);
+    }
+  return gcd1;
+}
+
+/* Divides all the coefficients of VECTOR by their greatest common
+   divisor.  VECTOR has SIZE elements.  */
+
+static inline void
+lambda_vector_normalize (lambda_vector vector, int size)
+{
+  int g = lambda_vector_gcd (vector, size);
+
+  if (g > 1)
+    lambda_vector_div_const (vector, vector, size, g);
+}
+
+/* Compute the scalar product VEC1 by VEC2.  Both vectors should have
+   the same SIZE.  */
+
+static inline int
+lambda_vector_scalar_product (lambda_vector vec1, lambda_vector vec2,
+			      int size)
+{
+  int i, res = 0;
+
+  for (i = 0; i < size; i++)
+    res += vec1[i] * vec2[i];
+
+  return res;
+}
+
+
+/* Compute RES as a linear combination of R1 and R2, such that
+   lambda_vector_scalar_product (A, RES, SIZE) == 0.  All vectors have
+   length SIZE.  */
+
+static inline void 
+lambda_vector_linear_combine_1 (lambda_vector r1,
+				lambda_vector r2,
+				lambda_vector a,
+				lambda_vector res,
+				int size)
+{
+  int lambda1 = lambda_vector_scalar_product (r2, a, size);
+  int lambda2 = lambda_vector_scalar_product (r1, a, size);
+  int g = gcd (lambda1, lambda2);
+
+  lambda_vector_add_mc (r1, lambda1 / g,
+			r2, - (lambda2 / g),
+			res, size);
+}
+
+/* Allocate a larger matrix NEW_M with NM lines and NN columns and
+   copy the content of the OM lines by ON columns OLD_M matrix.  */
+
+static inline void
+lambda_matrix_extend (lambda_matrix old_m, int om, int on, 
+		      lambda_matrix *new_m, int nm, int nn)
+{
+  if (om > nm || on > nn)
+    *new_m = old_m;
+
+  *new_m = lambda_matrix_new (nm, nn);
+  lambda_matrix_copy (old_m, *new_m, om, on);
+}
+
+/* Swap rows R1 and R2 in matrix MAT.  */
+
+static inline void
+lambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2)
+{
+  lambda_vector row;
+
+  row = mat[r1];
+  mat[r1] = mat[r2];
+  mat[r2] = row;
+}
+
+/* Swap values contained in columns C1 and C2 of size SIZE in matrix
+   MAT.  */
+
+static inline void
+lambda_matrix_exchange_columns (lambda_matrix mat, int c1, int c2,
+				int size)
+{
+  int i, val;
+
+  for (i = 0; i < size; i++)
+    {
+      val = mat[i][c1];
+      mat[i][c1] = mat[i][c2];
+      mat[i][c2] = val;
+    }
+}
+
+/* Copy values contained in column C1 to column C2 of size SIZE in
+   matrix MAT.  */
+
+static inline void
+lambda_matrix_copy_columns (lambda_matrix mat, int c1, int c2,
+			    int size)
+{
+  int i;
+
+  for (i = 0; i < size; i++)
+    mat[i][c2] = mat[i][c1];
+}
+
+/* Copy values contained in column C of size SIZE from matrix MAT to
+   vector V.  */
+
+static inline void
+lambda_matrix_copy_c2v (lambda_matrix mat, int c, lambda_vector v,
+			int size)
+{
+  int i;
+
+  for (i = 0; i < size; i++)
+    v[i] = mat[i][c];
+}
+
+/* Copy values contained in vector V of size SIZE to column C of
+   matrix MAT.  */
+
+static inline void
+lambda_matrix_copy_v2c (lambda_vector v, lambda_matrix mat, int c,
+			int size)
+{
+  int i;
+
+  for (i = 0; i < size; i++)
+    mat[i][c] = v[i];
+}
+
+/* Swap the values contained in vectors V1 and V2 of size SIZE.  */
+
+static inline void
+lambda_vector_exchange (lambda_vector v1, lambda_vector v2, int size)
+{
+  int i, val;
+
+  for (i = 0; i < size; i++)
+    {
+      val = v1[i];
+      v1[i] = v2[i];
+      v2[i] = val;
+    }
+}
+
 #endif /* LAMBDA_H  */
 
Index: omega.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/Attic/omega.c,v
retrieving revision 1.1.2.3
diff -d -u -p -r1.1.2.3 omega.c
--- omega.c	20 Apr 2005 17:34:40 -0000	1.1.2.3
+++ omega.c	21 Apr 2005 15:05:50 -0000
@@ -84,7 +84,7 @@ static int doItAgain;
 static int conservative = 0;
 
 static int next_key;
-static char wild_name[200][20];
+static char wild_name[200][40];
 static int next_wild_card = 0;
 static int trace = 1;
 static int depth = 0;
@@ -188,10 +188,13 @@ omega_eqn_is_red (eqn e, int desired_res
 static inline char *
 omega_variable_to_str (int i)
 {
-  if (i <= 0)
-    return wild_name[-i];
+  if (0 <= i && i < 20)
+    return wild_name[i];
 
-  return (*omega_current_get_var_name)(i);
+  if (-20 < i && i < 0)
+    return wild_name[40 - i];
+
+  gcc_unreachable ();
 }
 
 /* Do nothing function: used for default initializations.  */
@@ -951,6 +954,42 @@ omega_add_new_wild_card (omega_pb pb)
   return i;
 }
 
+/* Delete inequality E from problem PB that has NV variables.  */
+
+static void
+omega_delete_geq (omega_pb pb, int e, int nv)
+{
+  if (dump_file && (dump_flags & TDF_DETAILS))
+    {
+      fprintf (dump_file, "Deleting %d (last:%d): ", e, pb->num_geqs-1);
+      omega_print_geq (dump_file, pb, &pb->geqs[e]);
+      fprintf (dump_file, "\n");
+    }
+
+  if (e < pb->num_geqs-1)
+    omega_copy_eqn (&pb->geqs[e], &pb->geqs[pb->num_geqs - 1], nv);
+
+  pb->num_geqs--;
+}
+
+/* Delete extra inequality E from problem PB that has NV variables.  */
+
+static void
+omega_delete_geq_extra (omega_pb pb, int e, int nV)
+{
+  if (dump_file && (dump_flags & TDF_DETAILS))
+    {
+      fprintf (dump_file, "Deleting %d: ",e);
+      omega_print_geq_extra (dump_file, pb, &pb->geqs[e]);
+      fprintf (dump_file, "\n");
+    }
+
+  if (e < pb->num_geqs-1)
+    omega_copy_eqn (&pb->geqs[e], &pb->geqs[pb->num_geqs - 1],(nV));
+  pb->num_geqs--;
+}
+
+
 /* Remove variable I from problem PB.  */
 
 static void
@@ -5080,7 +5119,6 @@ omega_constrain_variable_value (omega_pb
       omega_copy_eqn (&pb->eqs[e], &pb->subs[k],
 		      pb->num_vars);
       pb->eqs[e].coef[0] -= value;
-
     }
   else
     {
@@ -5342,27 +5380,46 @@ omega_initialize (void)
     hash_master[i].touched = -1;
 
   sprintf (wild_name[0], "1");
-  sprintf (wild_name[1], "alpha");
-  sprintf (wild_name[2], "beta");
-  sprintf (wild_name[3], "gamma");
-  sprintf (wild_name[4], "delta");
-  sprintf (wild_name[5], "tau");
-  sprintf (wild_name[6], "sigma");
-  sprintf (wild_name[7], "chi");
-  sprintf (wild_name[8], "omega");
-  sprintf (wild_name[9], "pi");
-  sprintf (wild_name[10], "ni");
-  sprintf (wild_name[11], "Alpha");
-  sprintf (wild_name[12], "Beta");
-  sprintf (wild_name[13], "Gamma");
-  sprintf (wild_name[14], "Delta");
-  sprintf (wild_name[15], "Tau");
-  sprintf (wild_name[16], "Sigma");
-  sprintf (wild_name[17], "Chi");
-  sprintf (wild_name[18], "Omega");
-  sprintf (wild_name[19], "Pi");
+  sprintf (wild_name[1], "a");
+  sprintf (wild_name[2], "b");
+  sprintf (wild_name[3], "c");
+  sprintf (wild_name[4], "d");
+  sprintf (wild_name[5], "e");
+  sprintf (wild_name[6], "f");
+  sprintf (wild_name[7], "g");
+  sprintf (wild_name[8], "h");
+  sprintf (wild_name[9], "i");
+  sprintf (wild_name[10], "j");
+  sprintf (wild_name[11], "k");
+  sprintf (wild_name[12], "l");
+  sprintf (wild_name[13], "m");
+  sprintf (wild_name[14], "n");
+  sprintf (wild_name[15], "o");
+  sprintf (wild_name[16], "p");
+  sprintf (wild_name[17], "q");
+  sprintf (wild_name[18], "r");
+  sprintf (wild_name[19], "s");
+  sprintf (wild_name[20], "t");
+  sprintf (wild_name[40 - 1], "alpha");
+  sprintf (wild_name[40 - 2], "beta");
+  sprintf (wild_name[40 - 3], "gamma");
+  sprintf (wild_name[40 - 4], "delta");
+  sprintf (wild_name[40 - 5], "tau");
+  sprintf (wild_name[40 - 6], "sigma");
+  sprintf (wild_name[40 - 7], "chi");
+  sprintf (wild_name[40 - 8], "omega");
+  sprintf (wild_name[40 - 9], "pi");
+  sprintf (wild_name[40 - 10], "ni");
+  sprintf (wild_name[40 - 11], "Alpha");
+  sprintf (wild_name[40 - 12], "Beta");
+  sprintf (wild_name[40 - 13], "Gamma");
+  sprintf (wild_name[40 - 14], "Delta");
+  sprintf (wild_name[40 - 15], "Tau");
+  sprintf (wild_name[40 - 16], "Sigma");
+  sprintf (wild_name[40 - 17], "Chi");
+  sprintf (wild_name[40 - 18], "Omega");
+  sprintf (wild_name[40 - 19], "Pi");
 
   omega_initialized = true;
 }
 
-
Index: omega.h
===================================================================
RCS file: /cvs/gcc/gcc/gcc/Attic/omega.h,v
retrieving revision 1.1.2.3
diff -d -u -p -r1.1.2.3 omega.h
--- omega.h	20 Apr 2005 17:34:40 -0000	1.1.2.3
+++ omega.h	21 Apr 2005 15:05:50 -0000
@@ -71,7 +71,6 @@ typedef struct omega_pb
   struct eqn subs[max_vars + 1];
 } *omega_pb;
 
-extern char *(*omega_current_get_var_name) (unsigned int);
 extern bool omega_reduce_with_subs;
 extern bool omega_verify_simplification;
 
@@ -170,45 +169,10 @@ single_var_geq (struct eqn e, int nv ATT
   return (e.key != 0 && -max_vars <= e.key && e.key <= max_vars);
 }
 
-/* Delete inequality E from problem PB that has NV variables.  */
-
-static inline void
-omega_delete_geq (omega_pb pb, int e, int nv)
-{
-  if (dump_file && (dump_flags & TDF_DETAILS))
-    {
-      fprintf (dump_file, "Deleting %d (last:%d): ", e, pb->num_geqs-1);
-      omega_print_geq (dump_file, pb, &pb->geqs[e]);
-      fprintf (dump_file, "\n");
-    }
-
-  if (e < pb->num_geqs-1)
-    omega_copy_eqn (&pb->geqs[e], &pb->geqs[pb->num_geqs - 1], nv);
-
-  pb->num_geqs--;
-}
-
-/* Delete extra inequality E from problem PB that has NV variables.  */
-
-static inline void
-omega_delete_geq_extra (omega_pb pb, int e, int nV)
-{
-  if (dump_file && (dump_flags & TDF_DETAILS))
-    {
-      fprintf (dump_file, "Deleting %d: ",e);
-      omega_print_geq_extra (dump_file, pb, &pb->geqs[e]);
-      fprintf (dump_file, "\n");
-    }
-
-  if (e < pb->num_geqs-1)
-    omega_copy_eqn (&pb->geqs[e], &pb->geqs[pb->num_geqs - 1],(nV));
-  pb->num_geqs--;
-}
-
-
 #define cant_do_omega abort
 
-/* Initialize an Omega problem.  */
+/* Initialize P as an Omega problem with NVARS variables and
+   NPROT safe variables.  */
 
 static inline void
 init_problem (omega_pb p, unsigned nvars, unsigned nprot)
Index: polyhedron.c
===================================================================
RCS file: polyhedron.c
diff -N polyhedron.c
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ polyhedron.c	21 Apr 2005 15:05:50 -0000
@@ -0,0 +1,965 @@
+/* Linear constraint systems.
+   Copyright (C) 2005 Free Software Foundation, Inc.
+   Contributed by Sebastian Pop <sebastian.pop@cri.ensmp.fr>
+
+This file is part of GCC.
+
+GCC is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License as published by the Free
+Software Foundation; either version 2, or (at your option) any later
+version.
+
+GCC is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with GCC; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+02111-1307, USA.  */
+
+/* This file defines a structure for polyhedra and operations on
+   polyhedra.  Polyhedra have two possible representations: 
+   
+   - polyhedra as a set of linear equalities and inequalities,
+   aka. linear constraint systems, defined in this file,
+
+   - polyhedra as a set of vertices, rays, and lines, aka. generating
+     systems.
+
+   References:
+
+   - IRISA Publication Interne 785: "A library for doing polyhedral
+   operations" Doran Wilde.
+
+*/
+
+#include "config.h"
+#include "system.h"
+#include "coretypes.h"
+#include "tm.h"
+#include "errors.h"
+#include "ggc.h"
+#include "tree.h"
+#include "varray.h"
+#include "sbitmap.h"
+#include "vec.h"
+#include "polyhedron.h"
+
+/* Return a new constraint system.  */
+
+csys
+csys_new (int dimension, int nb_eqs, int nb_constraints)
+{
+  csys res = ggc_alloc_cleared (sizeof (struct csys));
+
+  CSYS_DIMENSION (res) = dimension;
+  CSYS_NB_EQS (res) = nb_eqs;
+  CSYS_NB_CONSTRAINTS (res) = nb_constraints;
+  CSYS_M (res) = lambda_matrix_new (nb_constraints, dimension + 2);
+  return res;
+}
+
+/* Return a copy of CY.  */
+
+csys 
+csys_copy (csys cy)
+{
+  csys res = ggc_alloc_cleared (sizeof (csys));
+
+  CSYS_DIMENSION (res) = CSYS_DIMENSION (cy);
+  CSYS_NB_EQS (res) = CSYS_NB_EQS (cy);
+  CSYS_NB_CONSTRAINTS (res) = CSYS_NB_CONSTRAINTS (cy);
+  lambda_matrix_copy (CSYS_M (cy), CSYS_M (res),
+		      CSYS_NB_CONSTRAINTS (cy), CSYS_DIMENSION (cy) + 2);
+  return res;
+}
+
+/* Print to FILE the constraint system CY.  */
+
+void 
+csys_print (FILE *file, csys cy)
+{
+  int i, j;
+
+  if (cy == NULL)
+    {
+      fprintf (file, "null\n");
+      return;
+    }
+
+  fprintf (file, "(constraint_system \n");
+  fprintf (file, "  dim = %d\n", CSYS_DIMENSION (cy));
+  fprintf (file, "  nb_eqs = %d\n", CSYS_NB_EQS (cy));
+  fprintf (file, "  nb_constraints = %d\n", CSYS_NB_CONSTRAINTS (cy));
+  fprintf (file, "  constraints_matrix = (\n");
+  fprintf (file, "%d %d\n", CSYS_NB_CONSTRAINTS (cy), CSYS_DIMENSION (cy) + 2);
+  print_lambda_matrix (file, CSYS_M (cy), CSYS_NB_CONSTRAINTS (cy),
+		       CSYS_DIMENSION (cy) + 2);
+  fprintf (file, "  )\n");
+  fprintf (file, "  constraints = (\n");
+
+  for (i = 0; i < CSYS_NB_CONSTRAINTS (cy); i++)
+    {
+      if (CSYS_IS_EQ (cy, i))
+	fprintf (file, "  0  = ");
+      else
+	fprintf (file, "  0 <= ");
+
+      for (j = 0; j < CSYS_DIMENSION (cy); j++)
+	if (CSYS_VEC (cy, i, j))
+	  fprintf (file, "%3d.x%d + ", CSYS_VEC (cy, i, j), j);
+      fprintf (file, "%3d\n", CSYS_CST (cy, i));
+    }
+  fprintf (file, "  )\n");
+  fprintf (file, ")\n");
+}
+
+/* Debugging function for gdb calls.  */
+
+void
+debug_csys (csys cy)
+{
+  csys_print (stderr, cy);
+}
+
+
+
+/* Return a new generating system.  */
+
+gsys
+gsys_new (int dimension, int nb_lines, int nb_rays)
+{
+  gsys res = ggc_alloc_cleared (sizeof (struct gsys));
+
+  GSYS_DIMENSION (res) = dimension;
+  GSYS_NB_LINES (res) = nb_lines;
+  GSYS_NB_RAYS (res) = nb_rays;
+  GSYS_M (res) = lambda_matrix_new (nb_rays, dimension + 2);
+  return res;
+}
+
+/* Return a copy of GY.  */
+
+gsys
+gsys_copy (gsys gy)
+{
+  gsys res = ggc_alloc_cleared (sizeof (gsys));
+
+  GSYS_DIMENSION (res) = GSYS_DIMENSION (gy);
+  GSYS_NB_LINES (res) = GSYS_NB_LINES (gy);
+  GSYS_NB_RAYS (res) = GSYS_NB_RAYS (gy);
+  lambda_matrix_copy (GSYS_M (gy), GSYS_M (res),
+		      GSYS_NB_RAYS (gy), GSYS_DIMENSION (gy) + 2);
+  return res;
+}
+
+/* Print to FILE the generating system GY.  */
+
+void
+gsys_print (FILE *file, gsys gy)
+{
+  if (gy == NULL)
+    {
+      fprintf (file, "null\n");
+      return;
+    }
+
+  fprintf (file, "(gsys \n");
+  fprintf (file, "  dim = %d\n", GSYS_DIMENSION (gy));
+  fprintf (file, "  nb_lines = %d\n", GSYS_NB_LINES (gy));
+  fprintf (file, "  nb_rays = %d\n", GSYS_NB_RAYS (gy));
+  fprintf (file, "  rays_matrix = (\n");
+  fprintf (file, "%d %d\n", GSYS_NB_RAYS (gy), GSYS_DIMENSION (gy) + 2);
+  print_lambda_matrix (file, GSYS_M (gy), GSYS_NB_RAYS (gy), 
+		       GSYS_DIMENSION (gy) + 2);
+  fprintf (file, "  )\n");
+  fprintf (file, ")\n");
+}
+
+/* Debugging function for gdb calls.  */
+
+void debug_gsys (gsys gy)
+{
+  gsys_print (stderr, gy);
+}
+
+
+
+/* Return a new polyhedron from a constraint system CY and a
+   generating system GY.  */
+
+polyhedron
+polyhedron_new (csys cy, gsys gy, lambda_matrix sat)
+{
+  polyhedron pol = ggc_alloc_cleared (sizeof (struct polyhedron));
+
+  POLYH_CSYS (pol) = cy;
+  POLYH_GSYS (pol) = gy;
+  POLYH_SAT (pol) = sat;
+  return pol;
+}
+
+
+
+/* Translate the constraint system CY to the omega problem RES.  RES
+   should be allocated before calling this function.  */
+
+void 
+csys_to_omega (csys cy, omega_pb res)
+{
+  int i, j;
+  omega_initialize ();
+  init_problem (res, CSYS_DIMENSION (cy), 0);
+
+  for (i = 0; i < CSYS_NB_CONSTRAINTS (cy); i++)
+    {
+      if (CSYS_IS_EQ (cy, i))
+	{
+	  unsigned idx = prob_add_zero_eq (res, omega_black);
+
+	  for (j = 1; j <= CSYS_DIMENSION (cy); j++)
+	    res->eqs[idx].coef[j] = CSYS_ELT (cy, i, j);
+	  res->eqs[idx].coef[0] = CSYS_CST (cy, i);
+	}
+      else
+	{
+	  unsigned idx = prob_add_zero_geq (res, omega_black);
+
+	  for (j = 1; j <= CSYS_DIMENSION (cy); j++)
+	    res->geqs[idx].coef[j] = CSYS_ELT (cy, i, j);
+	  res->geqs[idx].coef[0] = CSYS_CST (cy, i);
+	}
+    }
+}
+
+
+
+/* To generating systems.  */
+
+/* Compute RES as a linear combination of R1 and R2 such that 
+   RES [POS] == 0.  All vectors have length SIZE.  */
+
+static inline void 
+linear_combine (lambda_vector r1,
+		lambda_vector r2,
+		lambda_vector res,
+		int pos, int size)
+{
+  int lambda1 = r2[pos];
+  int lambda2 = r1[pos];
+  int g = gcd (lambda1, lambda2);
+
+  lambda_vector_add_mc (r1 + 1, lambda1 / g,
+			r2 + 1, - (lambda2 / g),
+			res + 1, size);
+  lambda_vector_normalize (res + 1, size);
+}
+
+/* Compute a minimal system of equations using Gausian elimination
+   method.  MAT is a matrix of constraints in which the first NB_EQ
+   constraints are equations.  The dimension of the homogenous system
+   is DIMENSION.  The function returns the rank of the matrix.  */
+
+static int
+csys_Gauss (lambda_matrix mat,
+	    int nb_eq,
+	    int dimension,
+	    int nb_rows)
+{
+  int i, j, k, rank;
+  lambda_vector column_index = lambda_vector_new (dimension);
+
+  rank=0;
+  for (j = 1; j <= dimension; j++)
+    {
+      for (i = rank; i < nb_eq; i++)
+	if (mat[i][j])
+	  break;
+
+      if (i != nb_eq)
+	{
+	  if (i != rank)
+	    lambda_vector_exchange (mat[rank] + 1, mat[i] + 1, dimension);
+
+	  lambda_vector_normalize (mat[rank] + 1, dimension);
+
+	  if (mat[rank][j] < 0)
+	    lambda_vector_negate (mat[rank] + 1, mat[rank] + 1, dimension);
+
+	  for (i = i + 1; i < nb_eq; i++)
+	    if (mat[i][j])
+	      linear_combine (mat[i], mat[rank], mat[i], j, dimension);
+
+	  column_index[rank] = j;
+	  rank++;
+	}
+    }
+
+  for (k = rank - 1; k >= 0; k--)
+    {
+      j = column_index[k];
+      
+      for (i = 0; i < k; i++)
+	if (mat[i][j])
+	  linear_combine (mat[i], mat[k], mat[i], j, dimension);
+      
+      for (i = nb_eq; i < nb_rows; i++)
+	if (mat[i][j])
+	  linear_combine (mat[i], mat[k], mat[i], j, dimension);
+    }
+
+  return rank;
+}
+
+/* Initialize CONS_P, RES_P, NB_CONS_P and NB_RAYS_P to the empty
+   polyhedron in dimension DIM.  */
+
+static polyhedron
+empty_polyhedron (polyhedron pol)
+{
+  int i;
+  int dim = POLYH_DIM (pol);
+
+  POLYH_CONS (pol) = lambda_matrix_new (dim + 1, dim + 2);
+  for (i = 0; i <= dim; i++)
+    POLYH_CONS (pol)[i][i+1] = 1;
+
+  POLYH_NB_CONS (pol) = dim + 1;
+  POLYH_NB_RAYS (pol) = 0;
+  POLYH_RAYS (pol) = NULL;
+
+  return pol;
+}
+
+#define MARK_SATURATED(S, R, C) (S[R][C] = 0)
+#define SATURATED_P(S, R, C) (S[R][C] == 0)
+#define MARK_NOT_SATURATED(S, R, C) (S[R][C] = 1)
+#define NOT_SATURATED_P(S, R, C) (S[R][C] == 1)
+
+/* R_TO_C should be true when translating from rays to constraints.
+
+   Compute the dual matrix representation RAYS_P starting from a an
+   initial representation START.  Because the rays and constraints
+   representations are dual, this algorithm can translate from one to
+   the other and vice versa.  The algorithm has been proposed by
+   Motzkin, then improved by Doran Wilde.
+
+   Notes: 
+   
+   - a ray "I" does not saturate a constraint "K" iff the scalar
+   product is not null: "I dot K != 0".
+
+   - matrix layouts: 
+   START has the following structure:
+   
+   |IS_INEQ  VECTORS  CST
+   |   0      1 0 5    9
+   |   1      1 2 3   -8
+
+   RES has the following structure:
+   
+   |IS_RAY   VECTORS  CST
+   |   0      1 0 5    0
+   |   1      1 2 3    3
+   
+   When IS_RAY is 0 VECTORS represent a line and constant is 0.
+   
+   When IS_RAY is 1 VECTORS represent: a ray if CST is 0, otherwise it
+   is a vertex represented in the homogeneous space.
+
+   Reference: IRISA Publication Interne 785: "A library for doing
+   polyhedral operations" Doran Wilde, pages 29 to 34.
+*/
+
+static polyhedron
+csys_gsys_compute_dual (polyhedron pol, bool r_to_c)
+{
+  lambda_matrix cons = POLYH_CONS (pol);
+  lambda_matrix rays = POLYH_RAYS (pol);
+  lambda_matrix sat = POLYH_SAT (pol);
+
+  int nb_constraints = POLYH_NB_CONS (pol);
+  int nb_rays = POLYH_NB_RAYS (pol);
+  int max_ray = nb_rays;
+  int nb_lines = POLYH_NB_LINES (pol);
+  int dim = POLYH_DIM (pol);
+
+  int i, j, k, l;
+
+  /* From now on, working in the homogeneous space.  */
+  dim++;
+
+  /* Initialize a basis for the space.  */
+  for (i = 0; i < nb_lines; i++)
+    rays[i][i + 1] = 1;
+
+  /* Include the point at the origin.  */
+  rays[nb_lines][0] = 1;
+  rays[nb_lines][dim] = 1;
+
+  for (k = 0; k < nb_constraints; k++) 
+    {
+      int first_non_saturating = nb_rays;
+
+      for (i = 0; i < nb_rays; i++)
+	{
+	  rays[i][0] = lambda_vector_scalar_product (rays[i] + 1, cons[k] + 1,
+						    dim);
+	  if (rays[i][0] && i < first_non_saturating)
+	    first_non_saturating = i;
+	}
+
+      if (first_non_saturating < nb_lines)
+	{
+	  nb_lines--;
+	  lambda_vector_exchange (rays[nb_lines], rays[first_non_saturating],
+				  dim + 1);
+
+	  for (i = 0; i < nb_rays; i++)
+	    {
+	      if (i == nb_lines)
+		{
+		  if (rays[nb_lines][0] < 0)
+		    lambda_vector_negate (rays[nb_lines], rays[nb_lines],
+					  dim + 1);
+		}
+	      else if (rays[i][0] != 0)
+		linear_combine (rays[i], rays[nb_lines], rays[i], 0, dim);
+	    }
+
+	  if (cons[k][0] != 0)
+	    {
+	      for (i = 0; i < nb_constraints; i++)
+		MARK_SATURATED (sat, nb_lines, i);
+	      MARK_NOT_SATURATED (sat, nb_lines, k);
+	    }
+	  else
+	    {
+	      nb_rays--;
+	      lambda_vector_copy (rays[nb_rays], rays[nb_lines], dim + 1);
+	      lambda_vector_copy (sat[nb_rays], sat[nb_lines], nb_constraints);
+	    }
+	}
+
+      else 
+	{
+	  int saved_nb_rays = nb_rays;
+	  int inf_bound = nb_rays;
+	  int sup_bound = nb_lines;
+	  int equal_bound = nb_lines;
+
+	  /* Sort the rays matrix.  */
+	  while (inf_bound > sup_bound)
+	    {
+	      if (rays[sup_bound][0] == 0)
+		{
+		  /* R0 = rays that saturate the constraint.  */
+		  lambda_vector_exchange (rays[equal_bound], rays[sup_bound],
+					  dim + 1);
+		  lambda_vector_exchange (sat[equal_bound], sat[sup_bound],
+					  nb_constraints);
+		  equal_bound++;
+		  sup_bound++;
+		}
+	      else
+		{
+		  MARK_NOT_SATURATED (sat, sup_bound, k);
+
+		  if (rays[sup_bound][0] < 0)
+		    {
+		      /* R- = rays that don't verify the constraint.  */
+		      inf_bound--;
+		      if (inf_bound != sup_bound)
+			{
+			  lambda_vector_exchange (rays[inf_bound], 
+						  rays[sup_bound],
+						  dim + 1);
+			  lambda_vector_exchange (sat[inf_bound], 
+						  sat[sup_bound],
+						  nb_constraints);
+			}
+		    }
+		  else 
+		    /* R+ = rays that verify the constraint.  */
+		    sup_bound++;
+		}
+	    }
+
+	  for (i = equal_bound; i < sup_bound; i++)
+	    for (j = sup_bound; j < saved_nb_rays; j++)
+	      {
+		int common_constraints = 0;
+		bool ray_only;
+
+		for (l = 0; l <= nb_constraints; l++)
+		  if (SATURATED_P (sat, i, l)
+		      && SATURATED_P (sat, j, l))
+		    common_constraints++;
+
+		ray_only = (rays[i][dim] == 0 && rays[j][dim] == 0 
+			    && r_to_c == 0);
+
+		if (ray_only)
+		  common_constraints++;
+
+		if (common_constraints + nb_lines >= dim - 2)
+		  {
+		    int iter;
+
+		    for (iter = nb_lines; iter < saved_nb_rays; iter++)
+		      if ((iter != i) && (iter != j) 
+			  && !(ray_only && rays[iter][dim]))
+			{
+			  for (l = 0; l < nb_constraints; l++)
+			    if (NOT_SATURATED_P (sat, iter, l)
+				&& SATURATED_P (sat, i, l)
+				&& SATURATED_P (sat, j, l))
+			      break;
+
+			  if (l == nb_constraints)
+			    /* The ray is redundant.  */
+			    goto next_ray;
+			}
+
+		    if (nb_rays == max_ray)
+		      {
+			max_ray *= 2;
+			lambda_matrix_extend (rays, nb_rays, dim + 1,
+					      &rays, max_ray, dim + 1);
+			lambda_matrix_extend (sat, nb_rays, nb_constraints,
+					      &sat, max_ray, nb_constraints);
+		      }
+
+		    linear_combine (rays[j], rays[i], rays[nb_rays], 0, dim);
+
+		    for (l = 0; l < nb_constraints; l++)
+		      if (NOT_SATURATED_P (sat, i, l)
+			  || NOT_SATURATED_P (sat, j, l))
+			MARK_NOT_SATURATED (sat, nb_rays, l);
+		      else
+			MARK_SATURATED (sat, nb_rays, l);
+		    MARK_SATURATED (sat, nb_rays, k);
+
+		    nb_rays++;
+		  }
+
+	      next_ray:;
+	      }
+
+	  /* Eliminates all non extremal rays.  Motzkin showed that a
+	     convex combination of a pair (r+ from R+, r- from R-)
+	     results in an extreme ray iff the minimal face that
+	     contains them:
+	     - is dimension one greater than r- and r+, 
+	     - and only contains the two rays r- and r+.
+	  */
+	  j = (cons[k][0] != 0) ? sup_bound : equal_bound;
+	  i = nb_rays;
+
+	  while ((j < saved_nb_rays) && (i > saved_nb_rays)) 
+	    {
+	      i--;
+	      lambda_vector_copy (rays[i], rays[j], dim + 1);
+	      lambda_vector_copy (sat[i], sat[j], nb_constraints);
+	      j++;
+	    }
+
+	  if (j == saved_nb_rays) 
+	    nb_rays = i;
+	  else 
+	    nb_rays = j;
+	}
+    }
+
+  POLYH_RAYS (pol) = rays;
+  POLYH_NB_RAYS (pol) = nb_rays;
+  POLYH_SAT (pol) = sat;
+  return pol;
+}
+
+/* Returns a polyhedron after having reduced the constraints system
+   and the rays of polyhedron POL.  For reference, see: IRISA
+   Publication Interne 785: "A library for doing polyhedral
+   operations" Doran Wilde.  */
+
+static polyhedron
+reduce_system (polyhedron pol)
+{
+  int i, j, k;
+  int dim_rayspace;
+  int NbBid, NbUni, NbEq, NbIneq;
+  int NbBid2, NbUni2, NbEq2, NbIneq2;
+  int aux;
+
+  lambda_vector trace;
+  lambda_matrix tmp_rays, tmp_cons;
+
+  /* Step0.  */
+  for (i = 0, aux = 0; i < POLYH_NB_RAYS (pol); i++)
+    {  
+      POLYH_RAYS (pol)[i][0] = 0;
+      if (POLYH_RAYS (pol)[i][POLYH_DIM (pol) + 1])
+	aux++;              
+    }
+  if (aux == 0)
+    return empty_polyhedron (pol);
+
+  /* Step1.  */
+  for (j = 0, NbEq = 0; j < POLYH_NB_CONS (pol); j++)
+    {
+      POLYH_CONS (pol)[j][0] = 0;
+
+      for (i = 1; i < POLYH_DIM (pol) + 1; i++)
+	if (POLYH_CONS (pol)[j][i])
+	  break;
+
+      if (i == POLYH_DIM (pol) + 1)
+	{
+	  for (i = 0; i < POLYH_NB_RAYS (pol); i++)
+	    if (SATURATED_P (POLYH_SAT (pol), i, j))
+	      POLYH_CONS (pol)[j][0]++;
+
+	  if (POLYH_CONS (pol)[j][0] == POLYH_NB_RAYS (pol)
+	      && POLYH_CONS (pol)[j][POLYH_DIM (pol) + 1])
+	    return empty_polyhedron (pol);
+	
+	  POLYH_NB_CONS (pol)--;
+	  if (j == POLYH_NB_CONS (pol))
+	    continue;
+
+	  lambda_vector_exchange (POLYH_CONS (pol)[j],
+				  POLYH_CONS (pol)[POLYH_NB_CONS (pol)],
+				  POLYH_DIM (pol) + 2);
+	  lambda_matrix_exchange_columns (POLYH_SAT (pol), j,
+					  POLYH_NB_CONS (pol), 
+					  POLYH_NB_RAYS (pol));
+	  j--;
+	  continue;
+	}
+      
+      for (i = 0; i < POLYH_NB_RAYS (pol); i++)
+	if (SATURATED_P (POLYH_SAT (pol), i, j))
+	  {
+	    POLYH_CONS (pol)[j][0]++;
+	    POLYH_RAYS (pol)[i][0]++;
+	  }
+
+      if (POLYH_CONS (pol)[j][0] == POLYH_NB_RAYS (pol))
+	NbEq++;
+    }
+
+  for (i = 0, NbBid = 0; i < POLYH_NB_RAYS (pol); i++)
+    {
+      if (POLYH_RAYS (pol)[i][POLYH_DIM (pol) + 1] == 0)
+	POLYH_RAYS (pol)[i][0]++;
+
+      if (POLYH_RAYS (pol)[i][0] == POLYH_NB_CONS (pol) + 1)
+	NbBid++;
+    }
+
+  /* Step2.  */
+  for (i = 0; i < NbEq; i++)
+    if (POLYH_CONS (pol)[i][0] != POLYH_NB_RAYS (pol))
+      {
+	lambda_vector temp1 = lambda_vector_new (POLYH_DIM (pol) + 2);
+
+	k = i + 1;
+	while (POLYH_CONS (pol)[k][0] != POLYH_NB_RAYS (pol) 
+	       && k < POLYH_NB_CONS (pol))
+	  k++;
+
+	if (k == POLYH_NB_CONS (pol))
+	  break;
+
+	lambda_vector_copy (POLYH_CONS (pol)[k], temp1, POLYH_DIM (pol) + 2);
+	lambda_matrix_copy_c2v (POLYH_SAT (pol), k, temp1, POLYH_NB_RAYS (pol));
+	for (; k > i; k--)
+	  {
+	    lambda_vector_copy (POLYH_CONS (pol)[k-1],
+				POLYH_CONS (pol)[k], POLYH_DIM (pol) + 2);
+	    lambda_matrix_copy_columns (POLYH_SAT (pol), k-1, k,
+					POLYH_NB_RAYS (pol));
+	  }
+	lambda_vector_copy (temp1, POLYH_CONS (pol)[i], POLYH_DIM (pol) + 2);
+	lambda_matrix_copy_v2c (temp1, POLYH_SAT (pol), i, POLYH_NB_RAYS (pol));
+      }
+
+  /* Step3.  */
+  NbEq2 = csys_Gauss (POLYH_CONS (pol), NbEq, POLYH_DIM (pol) + 1,
+		      POLYH_NB_CONS (pol));
+  if (NbEq2 >= POLYH_DIM (pol) + 1)
+    return empty_polyhedron (pol);
+
+  /* Step4.  */
+  for (i = 0, k = POLYH_NB_RAYS (pol); i < NbBid && k > i; i++)
+    if (POLYH_RAYS (pol)[i][0] != POLYH_NB_CONS (pol) + 1)
+      {
+	while (--k > i && POLYH_RAYS (pol)[k][0] != POLYH_NB_CONS (pol) + 1)
+	  ;
+
+	lambda_vector_exchange (POLYH_RAYS (pol)[i],
+				POLYH_RAYS (pol)[k], POLYH_DIM (pol) + 2);
+	lambda_vector_exchange (POLYH_SAT (pol)[i],
+				POLYH_SAT (pol)[k], POLYH_NB_CONS (pol));
+      }
+
+  /* Step5.  */
+  NbBid2 = csys_Gauss (POLYH_RAYS (pol), NbBid, POLYH_DIM (pol) + 1,
+		       POLYH_NB_RAYS (pol));
+  if (NbBid2 >= POLYH_DIM (pol) + 1)
+    return empty_polyhedron (pol);
+
+  /* Dimension of non-homogenous ray space.  */
+  dim_rayspace = POLYH_DIM (pol) - NbEq2 - NbBid2;
+
+  /* Step6.  */
+  for (j = 0, NbIneq = 0; j < POLYH_NB_CONS (pol); j++)
+    {
+      for (i = 1; i < POLYH_DIM (pol) + 1; i++)
+	if (POLYH_CONS (pol)[j][i] != 0)
+	  break;
+      
+      if (i == POLYH_DIM (pol) + 1)
+	{  
+	  if (POLYH_CONS (pol)[j][0] == POLYH_NB_RAYS (pol) 
+	      && POLYH_CONS (pol)[j][POLYH_DIM (pol) + 1] != 0)
+	    return empty_polyhedron (pol);
+	
+	  /* Constraint is redundant.  */
+	  POLYH_CONS (pol)[j][0] = 2;
+	  continue;
+	}
+
+      if (POLYH_CONS (pol)[j][0] == 0
+	  || POLYH_CONS (pol)[j][0] < dim_rayspace)
+	/* Constraint is redundant.  */
+	POLYH_CONS (pol)[j][0] = 2;
+
+      else if (POLYH_CONS (pol)[j][0] == POLYH_NB_RAYS (pol))
+	/* Constraint is an equality.  */
+	POLYH_CONS (pol)[j][0] = 0;
+
+      else
+	{
+	  /* Constraint is an inequality.  */
+	  NbIneq++; 	
+	  POLYH_CONS (pol)[j][0] = 1;
+	}
+    }
+    
+  /* Step7.  */
+  for (j = 0, NbUni = 0; j < POLYH_NB_RAYS (pol); j++)
+    if (POLYH_RAYS (pol)[j][0] < dim_rayspace)
+      POLYH_RAYS (pol)[j][0] = 2;
+    else if (POLYH_RAYS (pol)[j][0] == POLYH_NB_CONS (pol) + 1)
+      POLYH_RAYS (pol)[j][0] = 0;
+    else
+      {
+	NbUni++;
+	POLYH_RAYS (pol)[j][0] = 1;
+      }
+
+  /* Step8.  */
+  tmp_rays = lambda_matrix_new (NbUni + NbBid2, POLYH_DIM (pol) + 2);
+  tmp_cons = lambda_matrix_new (NbIneq + NbEq2 + 1, POLYH_DIM (pol) + 2);
+
+  if (NbBid2)
+    lambda_matrix_copy (POLYH_RAYS (pol), tmp_rays, NbBid2,
+			POLYH_DIM (pol) + 2);
+  if (NbEq2)
+    lambda_matrix_copy (POLYH_CONS (pol), tmp_cons, NbEq2,
+			POLYH_DIM (pol) + 2);
+
+  /* Step9.  */
+  trace = lambda_vector_new (POLYH_NB_CONS (pol));
+
+  for (j = NbEq, NbIneq2 = 0; j < POLYH_NB_CONS (pol); j++)
+    if (POLYH_CONS (pol)[j][0] == 1)
+      { 
+	bool redundant;
+	for (k = 0; k < POLYH_NB_CONS (pol); k++)
+	  trace[k]=0;
+
+	for (i = NbBid; i < POLYH_NB_RAYS (pol); i++) 
+	  if (POLYH_RAYS (pol)[i][0] == 1)
+	    if (SATURATED_P (POLYH_SAT (pol), i, j))
+	      for (k = 0; k < POLYH_NB_CONS (pol); k++)
+		trace[k] |= POLYH_SAT (pol)[i][k];
+	
+	for (i = NbEq, redundant = false; i < POLYH_NB_CONS (pol); i++)
+	  if (POLYH_CONS (pol)[i][0] == 1 && i != j && !trace[i])
+	    {
+	      redundant = true;
+	      break;
+	    }
+
+	if (redundant)
+	  POLYH_CONS (pol)[j][0] = 2;
+	else
+	  {
+	    lambda_vector_copy (POLYH_CONS (pol)[j],
+				tmp_cons[NbEq2 + NbIneq2],
+				POLYH_DIM (pol) + 2);
+	    NbIneq2++;
+	  }
+      }
+
+  /* Step10.  */
+  trace = lambda_vector_new (POLYH_NB_RAYS (pol));
+
+  for (i = NbBid, NbUni2 = 0, aux = 0; i < POLYH_NB_RAYS (pol); i++)
+    if (POLYH_RAYS (pol)[i][0] == 1)
+      { 
+	bool redundant;
+	if (POLYH_RAYS (pol)[i][POLYH_DIM (pol) + 1] != 0)
+	  for (k = NbBid; k < POLYH_NB_RAYS (pol); k++)
+	    trace[k]=0;
+	else
+	  for (k = NbBid; k < POLYH_NB_RAYS (pol); k++)
+	    trace[k] = (POLYH_RAYS (pol)[k][POLYH_DIM (pol) + 1] != 0);
+
+	for (j = NbEq; j < POLYH_NB_CONS (pol); j++)
+	  if (POLYH_CONS (pol)[j][0] == 1
+	      && SATURATED_P (POLYH_SAT (pol), i, j))
+	    for (k = NbBid; k < POLYH_NB_RAYS (pol); k++)
+	      trace[k] |= POLYH_SAT (pol)[k][j];
+
+	for (j = NbBid, redundant = false; j < POLYH_NB_RAYS (pol); j++)
+	  if (POLYH_RAYS (pol)[j][0] == 1 && i != j && !trace[j])
+	    { 
+	      redundant = true;
+	      break;
+	    }
+
+	if (redundant)
+	  POLYH_RAYS (pol)[i][0] = 2;	
+	else
+	  {
+	    lambda_vector_copy (POLYH_RAYS (pol)[i],
+				tmp_rays[NbBid2 + NbUni2],
+				POLYH_DIM (pol) + 2);
+	    NbUni2++;
+	    if (POLYH_RAYS (pol)[i][POLYH_DIM (pol) + 1] == 0)
+	      aux++;
+	  }
+      }
+
+  if (aux >= dim_rayspace)
+    {
+      lambda_vector_clear (tmp_cons[NbEq2 + NbIneq2], POLYH_DIM (pol) + 2);
+      tmp_cons[NbEq2 + NbIneq2][0] = 1;
+      tmp_cons[NbEq2 + NbIneq2][POLYH_DIM (pol) + 1] = 1;
+      NbIneq2++;
+    }
+
+  POLYH_CONS (pol) = tmp_cons;
+  POLYH_RAYS (pol) = tmp_rays;
+  POLYH_NB_RAYS (pol) = NbBid2 + NbUni2;
+  POLYH_NB_CONS (pol) = NbEq2 + NbIneq2;
+  return pol;
+}
+
+/* Returns a polyhedron built from the RAYS matrix.  */
+
+polyhedron
+rays_to_polyhedron (lambda_matrix rays ATTRIBUTE_UNUSED)
+{
+  return NULL;
+}
+
+/* Returns a polyhedron built from the CONSTRAINT matrix.  */
+
+polyhedron
+cons_to_polyhedron (lambda_matrix constraints ATTRIBUTE_UNUSED)
+{
+  return NULL;
+}
+
+/* Returns a polyhedron built from the constraint system CY.  */
+
+polyhedron
+polyhedron_from_csys (csys cy)
+{
+  int dimension = CSYS_DIMENSION (cy);
+  int nb_lines = dimension;
+  int nb_rays = dimension + 1;
+  gsys gy = gsys_new (dimension, nb_lines, nb_rays);
+  lambda_matrix sat = lambda_matrix_new (nb_rays, CSYS_NB_CONSTRAINTS (cy));
+  polyhedron pol = polyhedron_new (cy, gy, sat);
+
+  pol = csys_gsys_compute_dual (pol, false);
+  pol = reduce_system (pol);
+
+  return pol;
+}
+
+
+/* Transform a polyhedron POL into another polyhedron according to a
+   given affine mapping function FUN.  FUN is a FUN_ROWS by FUN_COLS
+   matrix.  */
+
+polyhedron
+polyhedron_image (polyhedron pol, lambda_matrix fun,
+		  int fun_rows, int fun_cols)
+{
+  int i, j, k;
+  lambda_matrix rays;
+
+  if (POLYH_DIM (pol) + 1 != fun_cols)
+    {
+      POLYH_DIM (pol) = fun_rows;
+      return empty_polyhedron (pol);
+    }
+
+  rays = lambda_matrix_new (POLYH_NB_RAYS (pol), fun_rows + 1);
+
+  for (i = 0; i < POLYH_NB_RAYS (pol); i++)
+    {
+      rays[i][0] = POLYH_RAYS (pol)[i][0];
+      for (j = 0; j < fun_rows; j++)
+	{
+	  int sum = 0;
+	  for (k = 0; k < POLYH_DIM (pol) + 1; k++)
+	    sum += POLYH_RAYS (pol)[i][k+1] * fun[j][k];
+	  rays[i][j+1] = sum;
+	}
+    }
+
+  return rays_to_polyhedron (rays);
+}
+
+/* Given a polyhedron POL and a transformation matrix FUN, return the
+   polyhedron which when transformed by mapping function FUN gives
+   POL.  FUN is a FUN_ROWS by FUN_COLS matrix.  */
+
+polyhedron
+polyhedron_preimage (polyhedron pol, lambda_matrix fun,
+		     int fun_rows, int fun_cols)
+{
+  int i, j, k;
+  lambda_matrix constraints;
+  
+  if (POLYH_DIM (pol) + 1 != fun_rows)
+    {
+      POLYH_DIM (pol) = fun_cols;
+      return empty_polyhedron (pol);
+    }
+
+  constraints = lambda_matrix_new (POLYH_NB_CONS (pol), fun_cols + 1);
+
+  for (i = 0; i < POLYH_NB_CONS (pol); i++)
+    {
+      constraints[i][0] = POLYH_CONS (pol)[i][0];
+      for (j = 0; j < fun_cols; j++)
+	{
+	  int sum = 0;
+	  for (k = 0; k < POLYH_DIM (pol) + 1; k++)
+	    sum += POLYH_CONS (pol)[i][k+1] * fun[k][j];
+	  constraints[i][j+1] = sum;
+	}
+    }
+
+  return cons_to_polyhedron (constraints);
+}
Index: polyhedron.h
===================================================================
RCS file: polyhedron.h
diff -N polyhedron.h
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ polyhedron.h	21 Apr 2005 15:05:50 -0000
@@ -0,0 +1,132 @@
+/* Structure of polyhedra.
+   Copyright (C) 2005 Free Software Foundation, Inc.
+   Contributed by Sebastian Pop <sebastian.pop@cri.ensmp.fr>
+
+This file is part of GCC.
+
+GCC is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License as published by the Free
+Software Foundation; either version 2, or (at your option) any later
+version.
+
+GCC is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with GCC; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+02111-1307, USA.  */
+
+#ifndef GCC_POLYHEDRON_H
+#define GCC_POLYHEDRON_H
+
+#include "lambda.h"
+#include "omega.h"
+
+/* A constraint system is composed of a system of linear equalities
+   and inequalities.  Example: the following constraints system
+   
+   0  = x + 0y + 5z + 9
+   0 <= x + 2y + 3z - 8
+
+   is encoded by the matrix: 
+
+   |IS_INEQ  VECTORS  CST
+   |   0      1 0 5    9
+   |   1      1 2 3   -8
+*/
+
+typedef struct csys
+{
+  /* NB_CONSTRAINTS = NB_EQS + nb_ineqs.  */
+  int dimension, nb_eqs, nb_constraints;
+
+  /* CONSTRAINTS is NB_CONSTRAINTS lines by DIMENSION + 2 columns.  */
+  lambda_matrix constraints;
+
+} *csys;
+
+#define CSYS_DIMENSION(C) ((C)->dimension)
+#define CSYS_NB_EQS(C) ((C)->nb_eqs)
+#define CSYS_NB_CONSTRAINTS(C) ((C)->nb_constraints)
+
+/* The raw matrix.  */
+#define CSYS_M(C) ((C)->constraints)
+/* An element.  */
+#define CSYS_ELT(C, I, J) (CSYS_M (C)[I][J])
+/* The constant.  */
+#define CSYS_CST(C, I) (CSYS_ELT (C, I, CSYS_DIMENSION (C) + 1))
+/* The vector.  */
+#define CSYS_VECTOR(C, I) (CSYS_M (C)[I] + 1)
+/* An element of the vector.  */
+#define CSYS_VEC(C, I, J) (CSYS_VECTOR (C, I)[J])
+/* Test for an equality constraint.  */
+#define CSYS_IS_EQ(C, I) ((CSYS_ELT (C, I, 0)) == 0)
+/* Test for an inequality constraint.  */
+#define CSYS_IS_INEQ(C, I) ((CSYS_ELT (C, I, 0)) != 0)
+
+/* A generating system is composed of a set of rays, and lines.  */
+
+typedef struct gsys
+{
+  int dimension, nb_lines, nb_rays;
+
+  /* RAYS is NB_RAYS lines by DIMENSION + 2 columns.  */
+  lambda_matrix rays;
+
+} *gsys;
+
+#define GSYS_DIMENSION(P) ((P)->dimension)
+#define GSYS_NB_LINES(P) ((P)->nb_lines)
+#define GSYS_NB_RAYS(P) ((P)->nb_rays)
+#define GSYS_M(P) ((P)->rays)
+
+/* A polyhedron is described by two dual representations:
+   - a set of linear constraints, aka. constraint system,
+   - a set of rays and lines, aka. generating system.
+*/
+
+typedef struct polyhedron
+{
+  csys cy;
+  gsys gy;
+  lambda_matrix saturation;
+} *polyhedron;
+
+#define POLYH_CSYS(P) (P->cy)
+#define POLYH_GSYS(P) (P->gy)
+#define POLYH_SAT(P) (P->saturation)
+#define POLYH_CONS(P) CSYS_M (POLYH_CSYS (P))
+#define POLYH_RAYS(P) GSYS_M (POLYH_GSYS (P))
+#define POLYH_NB_RAYS(P) GSYS_NB_RAYS (POLYH_GSYS (P))
+#define POLYH_NB_LINES(P) GSYS_NB_LINES (POLYH_GSYS (P))
+#define POLYH_NB_CONS(P) CSYS_NB_CONSTRAINTS (POLYH_CSYS (P))
+#define POLYH_DIM(P) CSYS_DIMENSION (POLYH_CSYS (P))
+
+extern csys csys_new (int, int, int);
+extern csys csys_copy (csys);
+extern void csys_print (FILE *, csys);
+extern void debug_csys (csys);
+extern void csys_to_omega (csys, omega_pb);
+
+extern gsys gsys_new (int, int, int);
+extern gsys gsys_copy (gsys);
+extern void gsys_print (FILE *, gsys);
+extern void debug_gsys (gsys);
+
+extern polyhedron polyhedron_new (csys, gsys, lambda_matrix);
+extern polyhedron polyhedron_from_csys (csys);
+extern polyhedron polyhedron_from_gsys (gsys);
+extern polyhedron polyhedron_copy (polyhedron);
+extern void polyhedron_print (FILE *, polyhedron);
+extern void debug_polyhedron (polyhedron);
+
+extern polyhedron polyhedron_image (polyhedron, lambda_matrix, int, int);
+extern polyhedron polyhedron_preimage (polyhedron, lambda_matrix, int, int);
+extern polyhedron rays_to_polyhedron (lambda_matrix);
+extern polyhedron cons_to_polyhedron (lambda_matrix);
+
+
+#endif
Index: tree-data-ref.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/tree-data-ref.c,v
retrieving revision 2.15.4.20
diff -d -u -p -r2.15.4.20 tree-data-ref.c
--- tree-data-ref.c	19 Apr 2005 16:36:50 -0000	2.15.4.20
+++ tree-data-ref.c	21 Apr 2005 15:05:51 -0000
@@ -350,28 +350,6 @@ tree_fold_divides_p (tree type, 
     (fold (build (MINUS_EXPR, type, a, tree_fold_gcd (a, b))));
 }
 
-/* Compute the greatest common denominator of two numbers using
-   Euclid's algorithm.  */
-
-static int 
-gcd (int a, int b)
-{
-  
-  int x, y, z;
-  
-  x = abs (a);
-  y = abs (b);
-
-  while (x>0)
-    {
-      z = y % x;
-      y = x;
-      x = z;
-    }
-
-  return (y);
-}
-
 /* Returns true iff A divides B.  */
 
 static inline bool 
@@ -3499,6 +3477,338 @@ access_functions_are_affine_or_constant_
   return true;
 }
 
+/* Initializes an equation using the information contained in the ACCESS_FUN.
+   Returns true when the operation succeeded.
+
+   CY is the constraint system.
+   EQ is the number of the equation to be initialized.
+   OFFSET is used for shifting the variables names in the constraints.
+   ACCESS_FUN is expected to be an affine chrec.  */
+
+static bool
+init_csys_eq_with_af (csys cy, unsigned eq, 
+		      unsigned int offset, tree access_fun, 
+		      varray_type vloops)
+{
+  switch (TREE_CODE (access_fun))
+    {
+    case POLYNOMIAL_CHREC:
+      {
+	tree left = CHREC_LEFT (access_fun);
+	tree right = CHREC_RIGHT (access_fun);
+	int var = CHREC_VARIABLE (access_fun);
+	unsigned var_idx;
+
+	if (TREE_CODE (right) != INTEGER_CST)
+	  return false;
+
+	/* Find the index of the current variable VAR_IDX in the
+	   VLOOPS array.  */
+	for (var_idx = 0; var_idx < VARRAY_ACTIVE_SIZE (vloops); var_idx++)
+	  {
+	    struct loop *loop = VARRAY_GENERIC_PTR (vloops, var_idx);
+	    if (loop->num == var)
+	      break;
+	  }
+
+	CSYS_VEC (cy, eq, offset + var_idx) = int_cst_value (right);
+	if (offset != 0)
+	  CSYS_VEC (cy, eq, var_idx) += int_cst_value (right);
+
+	switch (TREE_CODE (left))
+	  {
+	  case POLYNOMIAL_CHREC:
+	    return init_csys_eq_with_af (cy, eq, offset, left, vloops);
+
+	  case INTEGER_CST:
+	    CSYS_CST (cy, eq) += int_cst_value (left);
+	    return true;
+
+	  default:
+	    return false;
+	  }
+      }
+
+    case INTEGER_CST:
+      CSYS_CST (cy, eq) += int_cst_value (access_fun);
+      return true;
+
+    default:
+      return false;
+    }
+}
+
+/* Initialize VLOOPS with all the loops surrounding LOOP and inner to
+   STOP.  */
+
+static void
+find_loops_surrounding (struct loop *loop, struct loop *stop, 
+			varray_type *vloops)
+{
+  if (loop == stop
+      || loop == NULL
+      || loop->outer == NULL)
+    return;
+
+  VARRAY_PUSH_GENERIC_PTR (*vloops, loop);
+  find_loops_surrounding (loop->outer, stop, vloops);
+}
+
+/* Sets up the dependence constraint system for the data dependence
+   relation DDR.  Returns false when the constraint system cannot be
+   built, ie. when the test answers "don't know".  Returns true
+   otherwise, and when independence has been proved (using one of the
+   trivial dependence test), set MAYBE_DEPENDENT to false and the
+   DDR_CSYS is not initialized, otherwise set MAYBE_DEPENDENT to true.
+
+   Example: for setting up the dependence system corresponding to the
+   conflicting accesses 
+
+   loop_x
+     A[i] = ...
+     ... A[i+M]
+   endloop_x
+   
+   the following constraints come from the iteration domain: 
+
+   0 <= i <= N
+   0 <= i + di <= N
+
+   where di is the distance variable.  The conflicting elements
+   constraint inserted in the constraint system is:
+
+   i = i + di + M  
+
+   that gets simplified into 
+   
+   di + M = 0
+
+   Finally the constraint system initialized by the following function
+   looks like:
+
+   di + M = 0
+   0 <= i <= N
+   0 <= i + di <= N
+*/
+
+static bool
+init_csys_for_ddr (struct data_dependence_relation *ddr, bool *maybe_dependent)
+{
+  unsigned i;
+  int si;
+  struct data_reference *dra = DDR_A (ddr);
+  struct data_reference *drb = DDR_B (ddr);
+  struct loop *loop_a, *loop_b, *common_loop;
+  varray_type vloops;
+  unsigned nb_eqs, nb_loops, dimension;
+  unsigned eq, ineq;
+  csys cy;
+
+  *maybe_dependent = true;
+
+  /* Compute the size of the constraint system.
+     nb_loops_in_nest = number of loops surrounding both references.
+     dimension = 2 * nb_loops_in_nest.
+     nb_eqs = nb_subscripts
+     nb_ineqs = nb_loops_in_nest
+  */
+  VARRAY_GENERIC_PTR_INIT (vloops, 3, "vloops");
+  loop_a = bb_for_stmt (DR_STMT (dra))->loop_father;
+  loop_b = bb_for_stmt (DR_STMT (dra))->loop_father;
+  common_loop = find_common_loop (loop_a, loop_b);
+
+  if (common_loop != loop_a || common_loop != loop_b)
+    {
+      find_loops_surrounding (loop_a, common_loop, &vloops);
+      find_loops_surrounding (loop_b, common_loop, &vloops);
+    }
+  find_loops_surrounding (common_loop, NULL, &vloops);
+
+  nb_eqs = DDR_NUM_SUBSCRIPTS (ddr);
+  nb_loops = VARRAY_ACTIVE_SIZE (vloops);
+  dimension = 2 * nb_loops;
+  cy = csys_new (dimension, nb_eqs, 4 * nb_loops + nb_eqs);
+
+  /* For each subscript, insert an equality for representing the
+     conflicts.  */
+  for (i = 0, eq = 0; i < DDR_NUM_SUBSCRIPTS (ddr); i++)
+    {
+      tree access_fun_a = DR_ACCESS_FN (dra, i);
+      tree access_fun_b = DR_ACCESS_FN (drb, i);
+
+      /* ZIV test.  */
+      if (ziv_subscript_p (access_fun_a, access_fun_b))
+	{
+	  tree difference = chrec_fold_minus (integer_type_node, access_fun_a,
+					      access_fun_b);
+	  if (TREE_CODE (difference) == INTEGER_CST
+	      && !integer_zerop (difference))
+	    {
+	      /* There is no dependence.  */
+	      DDR_ARE_DEPENDENT (ddr) = chrec_known;
+	      *maybe_dependent = false;
+	      varray_clear (vloops);
+	      return true;
+	    }
+	}
+
+      access_fun_b = chrec_fold_multiply (chrec_type (access_fun_b), 
+					  access_fun_b, integer_minus_one_node);
+      
+      if (!init_csys_eq_with_af (cy, eq, 0, access_fun_a, vloops) 
+	  || !init_csys_eq_with_af (cy, eq, nb_loops, access_fun_b, vloops))
+	{
+	  /* There is probably a dependence, but the system of
+	     constraints cannot be built: answer "don't know".  */
+	  varray_clear (vloops);
+	  return false;
+	}
+
+      /* GCD test.  */
+      if (!int_divides_p (lambda_vector_gcd (CSYS_VECTOR (cy, eq), dimension),
+			  CSYS_CST (cy, eq)))
+	{
+	  /* There is no dependence.  */
+	  DDR_ARE_DEPENDENT (ddr) = chrec_known;
+	  *maybe_dependent = false;
+	  varray_clear (vloops);
+	  return true;
+	}
+      
+      eq++;
+    }
+
+  /* The rest are inequalities.  */
+  for (si = eq; si < CSYS_NB_CONSTRAINTS (cy); si++)
+    CSYS_ELT (cy, si, 0) = 1;
+
+  /* Insert the constraints corresponding to the iteration domain:
+     i.e. the loops surrounding the references "loop_x" and the
+     distance variables "dx".  */
+  for (i = 0, ineq = eq; i < VARRAY_ACTIVE_SIZE (vloops); i++)
+    {
+      struct loop *loop = VARRAY_GENERIC_PTR (vloops, i);
+      tree nb_iters = get_number_of_iters_for_loop (loop->num);
+
+      /* 0 <= loop_x */
+      CSYS_VEC (cy, ineq, i) = 1;
+      ineq++;
+
+      /* 0 <= loop_x + dx */
+      CSYS_VEC (cy, ineq, i) = 1;
+      CSYS_VEC (cy, ineq, i + nb_loops) = 1;
+      ineq++;
+
+      if (nb_iters != NULL_TREE
+	  && TREE_CODE (nb_iters) == INTEGER_CST)
+	{
+	  int nbi = int_cst_value (nb_iters);
+
+	  /* loop_x <= nb_iters */
+	  CSYS_VEC (cy, ineq, i) = -1;
+	  CSYS_CST (cy, ineq) = nbi;
+	  ineq++;
+
+	  /* loop_x + dx <= nb_iters */
+	  CSYS_VEC (cy, ineq, i) = -1;
+	  CSYS_VEC (cy, ineq, i + nb_loops) = -1;
+	  CSYS_CST (cy, ineq) = nbi;
+	  ineq++;
+	}
+    }
+
+  DDR_CSYS (ddr) = cy;
+
+  varray_clear (vloops);
+  return true;
+}
+
+/* Construct the constraint system for DDR, then solve it by polyhedra
+   solver.  */
+
+static void
+polyhedra_dependence_tester (struct data_dependence_relation *ddr)
+{
+  /* Translate to generating system (gs) representation, then detect
+     dep/indep.  */
+
+  debug_csys (DDR_CSYS (ddr));
+  DDR_POLYHEDRON (ddr) = polyhedron_from_csys (DDR_CSYS (ddr));
+  debug_gsys (POLYH_GSYS (DDR_POLYHEDRON (ddr)));
+
+  dependence_stats.num_dependence_undetermined++;
+  finalize_ddr_dependent (ddr, chrec_dont_know);
+}
+
+/* */
+
+static void
+omega_compute_classic_representations (struct data_dependence_relation *ddr ATTRIBUTE_UNUSED)
+{
+  
+}
+
+/* Construct the constraint system for DDR, then solve it using the
+   OMEGA solver.  */
+
+static void
+omega_dependence_tester (struct data_dependence_relation *ddr)
+{
+  int res = 0;
+  omega_pb pb = (omega_pb) xmalloc (sizeof (struct omega_pb));
+
+  if (dump_file && (dump_flags & TDF_DETAILS))
+    {
+      fprintf (dump_file, "(omega_dependence_tester \n");
+      dump_data_dependence_relation (dump_file, ddr);
+    }
+
+  csys_to_omega (DDR_CSYS (ddr), pb);
+  if (dump_file && (dump_flags & TDF_DETAILS))
+    omega_pretty_print_problem (dump_file, pb);
+
+  res = omega_solve_problem (pb, omega_unknown);
+
+  {
+    int i = 0, lower_bound, upper_bound, dist;
+    bool dist_known = false;
+    int dd_lt = 2;
+    int dd_eq = 4;
+    int dd_gt = 8;
+    int dir;
+
+    dir = omega_query_variable_signs (pb, i, dd_lt, dd_eq, dd_gt,
+				      lower_bound, upper_bound,
+				      &dist_known, &dist);
+    fprintf (stderr, "dir = %d\n", dir);
+  }
+
+  if (res == 0)
+    {
+      /* When there is no solution to the dependence problem, there is
+	 no dependence.  */
+      finalize_ddr_dependent (ddr, chrec_known);
+      dependence_stats.num_dependence_independent++;
+    }
+  else if (res == 1)
+    {
+      /* FIXME: Extract the classic representations: distances and
+	 directions. */
+      omega_compute_classic_representations (ddr);
+      dependence_stats.num_dependence_dependent++;
+    }
+  else
+    {
+      dependence_stats.num_dependence_undetermined++;
+      finalize_ddr_dependent (ddr, chrec_dont_know);
+    }
+  
+  if (dump_file && (dump_flags & TDF_DETAILS))
+    fprintf (dump_file, ")\n");
+
+  free (pb);
+}
+
 /* This computes the affine dependence relation between A and B.
    CHREC_KNOWN is used for representing the independence between two
    accesses, while CHREC_DONT_KNOW is used for representing the unknown
@@ -3531,13 +3841,32 @@ compute_affine_dependence (struct data_d
 
       if (access_functions_are_affine_or_constant_p (dra)
 	  && access_functions_are_affine_or_constant_p (drb))
-	subscript_dependence_tester (ddr);
+	{
+	  if (0)
+	    {
+	      bool maybe_dependent;
+
+	      if (!init_csys_for_ddr (ddr, &maybe_dependent))
+		goto csys_dont_know;
+
+	      if (maybe_dependent)
+		{
+  		  if (0)
+		    polyhedra_dependence_tester (ddr);
+		  else
+		    omega_dependence_tester (ddr);
+		}
+	    }
+	  else
+	    subscript_dependence_tester (ddr);
+	}
       
       /* As a last case, if the dependence cannot be determined, or if
 	 the dependence is considered too difficult to determine, answer
 	 "don't know".  */
       else
 	{
+	csys_dont_know:;
 	  dependence_stats.num_dependence_undetermined++;
 	  
 	  if (dump_file && (dump_flags & TDF_DETAILS))
Index: tree-data-ref.h
===================================================================
RCS file: /cvs/gcc/gcc/gcc/tree-data-ref.h,v
retrieving revision 2.5.4.8
diff -d -u -p -r2.5.4.8 tree-data-ref.h
--- tree-data-ref.h	19 Apr 2005 16:36:51 -0000	2.5.4.8
+++ tree-data-ref.h	21 Apr 2005 15:05:51 -0000
@@ -23,6 +23,8 @@ Software Foundation, 59 Temple Place - S
 #define GCC_TREE_DATA_REF_H
 
 #include "lambda.h"
+#include "omega.h"
+#include "polyhedron.h"
 
 struct data_reference
 {
@@ -177,6 +179,12 @@ struct data_dependence_relation
 
   /* The classic distance vector.  */
   lambda_vector dist_vect;
+
+  /* Representation of the dependence as a constraint system.  */
+  csys dependence_constraint_system;
+
+  /* Representation of the dependence as polyhedra.  */
+  polyhedron dependence_polyhedron;
 };
 
 #define DDR_A(DDR) DDR->a
@@ -191,6 +199,8 @@ struct data_dependence_relation
 #define DDR_SIZE_VECT(DDR) DDR->size_vect
 #define DDR_DIR_VECT(DDR) DDR->dir_vect
 #define DDR_DIST_VECT(DDR) DDR->dist_vect
+#define DDR_CSYS(DDR) DDR->dependence_constraint_system
+#define DDR_POLYHEDRON(DDR) DDR->dependence_polyhedron
 
 
 


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]