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
Param_Vertex_Common_Denominator(Param_Vertices * V)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 */
Param_Inner_Product(Value * constraint,Matrix * Vertex,Value * row)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
PL_P2PP(Polyhedron * Din,Polyhedron * Cin,struct barvinok_options * options)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 */
Param_Polyhedron_Facet(Param_Polyhedron * PP,Param_Domain * D,Value * constraint)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 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
TC_P2PP(Polyhedron * P,Polyhedron * C,struct barvinok_options * options)115 Param_Polyhedron *TC_P2PP(Polyhedron *P, Polyhedron *C,
116 struct barvinok_options *options)
117 {
118 return NULL;
119 }
120 #endif
121
Polyhedron2Param_Polyhedron(Polyhedron * P,Polyhedron * C,struct barvinok_options * options)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 case BV_CHAMBERS_ISL:
133 return ISL_P2PP(P, C, options);
134 break;
135 default:
136 assert(0);
137 }
138 }
139
140 #define INT_BITS (sizeof(unsigned) * 8)
141
142 /* Wegner's method for counting the number of ones in a bit vector */
bit_vector_count(unsigned * F,int F_len)143 int bit_vector_count(unsigned *F, int F_len)
144 {
145 int i;
146 int count = 0;
147
148 for (i = 0; i < F_len; ++i) {
149 unsigned v = F[i];
150 while (v) {
151 v &= v-1;
152 ++count;
153 }
154 }
155 return count;
156 }
157
Param_Vertex_Set_Facets(Param_Polyhedron * PP,Param_Vertices * V)158 void Param_Vertex_Set_Facets(Param_Polyhedron *PP, Param_Vertices *V)
159 {
160 if (!V->Facets) {
161 unsigned nparam = V->Vertex->NbColumns - 2;
162 int len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
163 int i, ix;
164 unsigned bx;
165 Vector *row = Vector_Alloc(1 + nparam + 1);
166
167 V->Facets = (unsigned *)calloc(len, sizeof(unsigned));
168 for (i = 0, ix = 0, bx = MSB; i < PP->Constraints->NbRows; ++i) {
169 Param_Inner_Product(PP->Constraints->p[i], V->Vertex, row->p);
170 if (First_Non_Zero(row->p+1, nparam+1) == -1)
171 V->Facets[ix] |= bx;
172 NEXT(ix, bx);
173 }
174 Vector_Free(row);
175 }
176 }
177
Param_Vertex_Cone(Param_Polyhedron * PP,Param_Vertices * V,struct barvinok_options * options)178 Polyhedron *Param_Vertex_Cone(Param_Polyhedron *PP, Param_Vertices *V,
179 struct barvinok_options *options)
180 {
181 int i, j, ix;
182 unsigned bx;
183 int len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
184 int n;
185 Matrix *M;
186 Polyhedron *C;
187 unsigned nvar = V->Vertex->NbRows;
188
189 if (!V->Facets)
190 Param_Vertex_Set_Facets(PP, V);
191 n = bit_vector_count(V->Facets, len);
192
193 M = Matrix_Alloc(n, 1+nvar+1);
194 assert(M);
195 for (i = 0, j = 0, ix = 0, bx = MSB; i < PP->Constraints->NbRows; ++i) {
196 if (V->Facets[ix] & bx)
197 Vector_Copy(PP->Constraints->p[i], M->p[j++], 1+nvar);
198 NEXT(ix, bx);
199 }
200 C = Constraints2Polyhedron(M, options->MaxRays);
201 assert(C);
202 Matrix_Free(M);
203 return C;
204 }
205