param_util.c

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