reduce_domain.c

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