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