param_util.c

   1 #include <assert.h>
   2 #include <barvinok/options.h>
   3 #include "param_util.h"
   4 #include "topcom.h"
   5 #include "isl_param_util.h"
   6 #include "config.h"
   7
   8 #define ALLOC(type) (type*)malloc(sizeof(type))
   9 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  10
  11 void Param_Vertex_Common_Denominator(Param_Vertices *V)
  12 {
  13     unsigned dim;
  14     Value lcm;
  15     int i;
  16
  17     assert(V->Vertex->NbRows > 0);
  18     dim = V->Vertex->NbColumns-2;
  19
  20     value_init(lcm);
  21
  22     value_assign(lcm, V->Vertex->p[0][dim+1]);
  23     for (i = 1; i < V->Vertex->NbRows; ++i)
  24         value_lcm(lcm, V->Vertex->p[i][dim+1], lcm);
  25
  26     for (i = 0; i < V->Vertex->NbRows; ++i) {
  27         if (value_eq(V->Vertex->p[i][dim+1], lcm))
  28             continue;
  29         value_division(V->Vertex->p[i][dim+1], lcm, V->Vertex->p[i][dim+1]);
  30         Vector_Scale(V->Vertex->p[i], V->Vertex->p[i],
  31                      V->Vertex->p[i][dim+1], dim+1);
  32         value_assign(V->Vertex->p[i][dim+1], lcm);
  33     }
  34
  35     value_clear(lcm);
  36 }
  37
  38 /* Plug in the parametric vertex Vertex (nvar x (nparam + 2))
  39  * in the constraint constraint (1 + nvar + nparam + 1).
  40  * The result is stored in row (1 + nparam + 1),
  41  * with the denominator in position 0.
  42  */
  43 void Param_Inner_Product(Value *constraint, Matrix *Vertex, Value *row)
  44 {
  45     unsigned nparam = Vertex->NbColumns - 2;
  46     unsigned nvar = Vertex->NbRows;
  47     int j;
  48     Value tmp, tmp2;
  49
  50     value_set_si(row[0], 1);
  51     Vector_Set(row+1, 0, nparam+1);
  52
  53     value_init(tmp);
  54     value_init(tmp2);
  55
  56     for (j = 0 ; j < nvar; ++j) {
  57         value_set_si(tmp, 1);
  58         value_assign(tmp2,  constraint[1+j]);
  59         if (value_ne(row[0], Vertex->p[j][nparam+1])) {
  60             value_assign(tmp, row[0]);
  61             value_lcm(row[0], row[0], Vertex->p[j][nparam+1]);
  62             value_division(tmp, row[0], tmp);
  63             value_multiply(tmp2, tmp2, row[0]);
  64             value_division(tmp2, tmp2, Vertex->p[j][nparam+1]);
  65         }
  66         Vector_Combine(row+1, Vertex->p[j], row+1, tmp, tmp2, nparam+1);
  67     }
  68     value_set_si(tmp, 1);
  69     Vector_Combine(row+1, constraint+1+nvar, row+1, tmp, row[0], nparam+1);
  70
  71     value_clear(tmp);
  72     value_clear(tmp2);
  73 }
  74
  75 static Param_Polyhedron *PL_P2PP(Polyhedron *Din, Polyhedron *Cin,
  76                                               struct barvinok_options *options)
  77 {
  78     unsigned MaxRays = options->MaxRays;
  79     if (MaxRays & (POL_NO_DUAL | POL_INTEGER))
  80         MaxRays = 0;
  81     return Polyhedron2Param_Domain(Din, Cin, MaxRays);
  82 }
  83
  84 /* Compute a dummy Param_Domain that contains all vertices of Param_Domain D
  85  * (which contains the vertices of P) that lie on the facet obtained by
  86  * saturating constraint c of P
  87  */
  88 Param_Domain *Param_Polyhedron_Facet(Param_Polyhedron *PP, Param_Domain *D,
  89                                      Value *constraint)
  90 {
  91     int nv;
  92     Param_Vertices *V;
  93     unsigned nparam = PP->V->Vertex->NbColumns-2;
  94     Vector *row = Vector_Alloc(1+nparam+1);
  95     Param_Domain *FD = ALLOC(Param_Domain);
  96     FD->Domain = 0;
  97     FD->next = 0;
  98
  99     nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
 100     FD->F = ALLOCN(unsigned, nv);
 101     memset(FD->F, 0, nv * sizeof(unsigned));
 102
 103     FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _i, _ix, _bx internal counters */
 104         int n;
 105         Param_Inner_Product(constraint, V->Vertex, row->p);
 106         if (First_Non_Zero(row->p+1, nparam+1) == -1)
 107             FD->F[_ix] |= _bx;
 108     END_FORALL_PVertex_in_ParamPolyhedron;
 109
 110     Vector_Free(row);
 111
 112     return FD;
 113 }
 114
 115 #ifndef POINTS2TRIANGS_PATH
 116 Param_Polyhedron *TC_P2PP(Polyhedron *P, Polyhedron *C,
 117                           struct barvinok_options *options)
 118 {
 119     return NULL;
 120 }
 121 #endif
 122
 123 Param_Polyhedron *Polyhedron2Param_Polyhedron(Polyhedron *P, Polyhedron *C,
 124                                               struct barvinok_options *options)
 125 {
 126     switch(options->chambers) {
 127     case BV_CHAMBERS_POLYLIB:
 128         return PL_P2PP(P, C, options);
 129         break;
 130     case BV_CHAMBERS_TOPCOM:
 131         return TC_P2PP(P, C, options);
 132         break;
 133     case BV_CHAMBERS_ISL:
 134         return ISL_P2PP(P, C, options);
 135         break;
 136     default:
 137         assert(0);
 138     }
 139 }
 140
 141 #define INT_BITS (sizeof(unsigned) * 8)
 142
 143 /* Wegner's method for counting the number of ones in a bit vector */
 144 int bit_vector_count(unsigned *F, int F_len)
 145 {
 146     int i;
 147     int count = 0;
 148
 149     for (i = 0; i < F_len; ++i) {
 150         unsigned v = F[i];
 151         while (v) {
 152             v &= v-1;
 153             ++count;
 154         }
 155     }
 156     return count;
 157 }
 158
 159 void Param_Vertex_Set_Facets(Param_Polyhedron *PP, Param_Vertices *V)
 160 {
 161     if (!V->Facets) {
 162         unsigned nparam = V->Vertex->NbColumns - 2;
 163         int len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
 164         int i, ix;
 165         unsigned bx;
 166         Vector *row = Vector_Alloc(1 + nparam + 1);
 167
 168         V->Facets = (unsigned *)calloc(len, sizeof(unsigned));
 169         for (i = 0, ix = 0, bx = MSB; i < PP->Constraints->NbRows; ++i) {
 170             Param_Inner_Product(PP->Constraints->p[i], V->Vertex, row->p);
 171             if (First_Non_Zero(row->p+1, nparam+1) == -1)
 172                 V->Facets[ix] |= bx;
 173             NEXT(ix, bx);
 174         }
 175         Vector_Free(row);
 176     }
 177 }
 178
 179 Polyhedron *Param_Vertex_Cone(Param_Polyhedron *PP, Param_Vertices *V,
 180                               struct barvinok_options *options)
 181 {
 182     int i, j, ix;
 183     unsigned bx;
 184     int len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
 185     int n;
 186     Matrix *M;
 187     Polyhedron *C;
 188     unsigned nvar = V->Vertex->NbRows;
 189     unsigned nparam = V->Vertex->NbColumns - 2;
 190
 191     if (!V->Facets)
 192         Param_Vertex_Set_Facets(PP, V);
 193     n = bit_vector_count(V->Facets, len);
 194
 195     M = Matrix_Alloc(n, 1+nvar+1);
 196     assert(M);
 197     for (i = 0, j = 0, ix = 0, bx = MSB; i < PP->Constraints->NbRows; ++i) {
 198         if (V->Facets[ix] & bx)
 199             Vector_Copy(PP->Constraints->p[i], M->p[j++], 1+nvar);
 200         NEXT(ix, bx);
 201     }
 202     C = Constraints2Polyhedron(M, options->MaxRays);
 203     assert(C);
 204     Matrix_Free(M);
 205     return C;
 206 }
 207
 208 /* Compute activity domain of a parametric vertex. */
 209 void Param_Vertex_Domain(Param_Vertices *V, Polyhedron *P)
 210 {
 211     int i;
 212     unsigned nparam = V->Vertex->NbColumns-2;
 213
 214     if (V->Domain)
 215         return;
 216
 217     V->Domain = Matrix_Alloc(P->NbConstraints, 1+nparam+1);
 218     for (i = 0; i < P->NbConstraints; ++i) {
 219         Param_Inner_Product(P->Constraint[i], V->Vertex, V->Domain->p[i]);
 220         value_assign(V->Domain->p[i][0], P->Constraint[i][0]);
 221     }
 222 }