reduce_domain.c

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