param_util.c

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