reduce_domain.c

   1 #include <assert.h>
   2 #include <barvinok/options.h>
   3 #include <barvinok/util.h>
   4 #include "reduce_domain.h"
   5
   6 Polyhedron *true_context(Polyhedron *P, Polyhedron *C, unsigned MaxRays)
   7 {
   8     unsigned nparam = C->Dimension;
   9     Polyhedron *tmp = Polyhedron_Project(P, nparam);
  10     C = DomainIntersection(tmp, C, MaxRays);
  11     Polyhedron_Free(tmp);
  12     return C;
  13 }
  14
  15 Vector *inner_point(Polyhedron *P)
  16 {
  17     Vector *average;
  18     int nv = 0;
  19     Value fc, fv;
  20     unsigned dim = P->Dimension;
  21     int i;
  22
  23     average = Vector_Alloc(1+dim+1);
  24
  25     POL_ENSURE_VERTICES(P);
  26     for (i = 0; i < P->NbRays; ++i) {
  27         if (value_zero_p(P->Ray[i][1+dim]))
  28             continue;
  29         ++nv;
  30         if (nv == 1) {
  31             Vector_Copy(P->Ray[i]+1, average->p+1, dim+1);
  32             continue;
  33         }
  34         if (nv == 2) {
  35             value_init(fc);
  36             value_init(fv);
  37         }
  38         value_assign(fc, average->p[1+dim]);
  39         value_lcm(average->p[1+dim], fc, P->Ray[i][1+dim]);
  40         value_division(fc, average->p[1+dim], fc);
  41         value_division(fv, average->p[1+dim], P->Ray[i][1+dim]);
  42         Vector_Combine(average->p+1, P->Ray[i]+1, average->p+1, fc, fv, dim);
  43     }
  44     if (nv > 1) {
  45         value_set_si(fc, nv);
  46         value_multiply(average->p[1+dim], average->p[1+dim], fc);
  47         Vector_Normalize(average->p+1, dim+1);
  48         value_clear(fc);
  49         value_clear(fv);
  50     }
  51     for (i = 0; i < P->NbRays; ++i) {
  52         if (value_notzero_p(P->Ray[i][1+dim]))
  53             continue;
  54         Vector_Add(average->p+1, P->Ray[i]+1, average->p+1, dim);
  55     }
  56     return average;
  57 }
  58
  59 int is_internal(Vector *point, Value *constraint)
  60 {
  61     int p;
  62     unsigned dim = point->Size-2;
  63
  64     Inner_Product(constraint+1, point->p+1, dim+1, &point->p[0]);
  65     if (value_notzero_p(point->p[0]))
  66         return value_pos_p(point->p[0]);
  67
  68     p = First_Non_Zero(constraint+1, dim);
  69     return value_pos_p(constraint[1+p]);
  70 }
  71
  72 Polyhedron *reduce_domain(Polyhedron *D, int nd,
  73                           Vector *inner, struct barvinok_options *options)
  74 {
  75     Polyhedron *rVD;
  76     Value c;
  77     int i;
  78     Matrix *constraints;
  79     int changed;
  80
  81     rVD = D->next ? DomainConvex(D, options->MaxRays) : Polyhedron_Copy(D);
  82
  83     /* If there is only one chamber, then we don't need to take care
  84      * of possible overlaps.
  85      * Plus, this decomposition may be the result of a recursive call
  86      * and then some of the assumptions used in determining whether
  87      * the domain is too small in geometric dimension no longer apply.
  88      */
  89     if (nd == 1) {
  90         return rVD;
  91     }
  92
  93     assert(rVD->Dimension == inner->Size-2);
  94     constraints = Polyhedron2Constraints(rVD);
  95     changed = 0;
  96     for (i = 0; i < constraints->NbRows; ++i) {
  97         if (!is_internal(inner, constraints->p[i])) {
  98             value_decrement(constraints->p[i][1+rVD->Dimension],
  99                             constraints->p[i][1+rVD->Dimension]);
 100             changed = 1;
 101         }
 102     }
 103     if (changed) {
 104         Polyhedron_Free(rVD);
 105         rVD = Constraints2Polyhedron(constraints, options->MaxRays);
 106     }
 107     Matrix_Free(constraints);
 108
 109     rVD = DomainConstraintSimplify(rVD, options->MaxRays);
 110     POL_ENSURE_FACETS(rVD);
 111     if (emptyQ(rVD)) {
 112         Domain_Free(rVD);
 113         return 0;
 114     }
 115
 116     value_init(c);
 117     barvinok_count_with_options(rVD, &c, options);
 118     if (value_zero_p(c)) {
 119         Domain_Free(rVD);
 120         rVD = 0;
 121     }
 122     value_clear(c);
 123
 124     return rVD;
 125 }