update pet for explicitly linking in gmp
[barvinok.git] / topcom.c
blob72d737c2c6d81218f4930c13672e13a617b6253f
1 #include <assert.h>
2 #include <barvinok/util.h>
3 #include <barvinok/options.h>
4 #include <unistd.h>
5 #include "normalization.h"
6 #include "topcom.h"
7 #include "config.h"
9 #define ALLOC(type) (type*)malloc(sizeof(type))
11 void run_points2triangs(pid_t *child, int *in, int *out)
13 int in_fd[2], out_fd[2];
15 if (pipe(in_fd))
16 assert(0);
17 if (pipe(out_fd))
18 assert(0);
19 *child = fork();
20 assert(*child >= 0);
21 if (!*child) {
22 int rc;
24 dup2(in_fd[0], 0);
25 dup2(out_fd[1], 1);
26 close(in_fd[0]);
27 close(out_fd[1]);
28 close(in_fd[1]);
29 close(out_fd[0]);
31 rc = execl(POINTS2TRIANGS_PATH, "points2triangs", "--regular", NULL);
32 assert(0);
34 close(in_fd[0]);
35 close(out_fd[1]);
36 *in = in_fd[1];
37 *out = out_fd[0];
40 struct domain {
41 Param_Domain domain;
42 int F_len;
45 static Param_Vertices *construct_vertex(unsigned *vertex_facets,
46 Matrix *Constraints,
47 int d, unsigned nparam, unsigned MaxRays)
49 unsigned nvar = Constraints->NbColumns-2 - nparam;
50 Matrix *A = Matrix_Alloc(nvar+1, nvar+1);
51 Matrix *inv = Matrix_Alloc(nvar+1, nvar+1);
52 Matrix *B = Matrix_Alloc(nvar, nparam+2);
53 Matrix *V = Matrix_Alloc(nvar, nparam+2);
54 Matrix *Domain = Matrix_Alloc(d-nvar, nparam+2);
55 Polyhedron *AD;
56 unsigned bx;
57 int i, j, ix;
58 int ok;
59 Param_Vertices *vertex;
61 for (j = 0, i = 0, ix = 0, bx = MSB; i < d; ++i) {
62 if ((vertex_facets[ix] & bx) == bx) {
63 Vector_Copy(Constraints->p[i]+1, A->p[j], nvar);
64 Vector_Oppose(Constraints->p[i]+1+nvar, B->p[j++], nparam+1);
66 NEXT(ix, bx);
68 assert(j == nvar);
69 value_set_si(A->p[nvar][nvar], 1);
70 ok = Matrix_Inverse(A, inv);
71 assert(ok);
72 Matrix_Free(A);
73 inv->NbRows = nvar;
74 inv->NbColumns = nvar;
75 Matrix_Product(inv, B, V);
76 Matrix_Free(B);
77 for (i = 0; i < nvar; ++i) {
78 value_assign(V->p[i][nparam+1], inv->p[nvar][nvar]);
79 Vector_Normalize(V->p[i], V->NbColumns);
81 Matrix_Free(inv);
82 for (j = 0, i = 0, ix = 0, bx = MSB; i < d; ++i) {
83 if ((vertex_facets[ix] & bx) == bx) {
84 NEXT(ix, bx);
85 continue;
87 Param_Inner_Product(Constraints->p[i], V, Domain->p[j]);
88 if (First_Non_Zero(Domain->p[j]+1, nparam+1) == -1)
89 vertex_facets[ix] |= bx;
90 else
91 value_set_si(Domain->p[j++][0], 1);
92 NEXT(ix, bx);
94 Domain->NbRows = j;
95 A = Matrix_Copy(Domain);
96 AD = Constraints2Polyhedron(A, MaxRays);
97 Matrix_Free(A);
98 POL_ENSURE_VERTICES(AD);
99 /* A vertex with a lower-dimensional activity domain
100 * saturates more facets than derived above and is actually
101 * the superimposition of two or more vertices.
102 * We therefore discard the domain and (ultimately)
103 * the chamber containing it.
104 * We keep the vertex, though, since it may reappear
105 * in other chambers, which will then likewise be discarded.
106 * The same holds if the activity domain is empty.
108 if (AD->NbEq > 0) {
109 Matrix_Free(Domain);
110 Domain = NULL;
112 Polyhedron_Free(AD);
113 vertex = calloc(1, sizeof(Param_Vertices));
114 vertex->Facets = vertex_facets;
115 vertex->Vertex = V;
116 vertex->Domain = Domain;
117 return vertex;
120 static int add_vertex_to_domain(Param_Vertices **vertices, int words,
121 unsigned *vertex_facets,
122 Matrix *Constraints, int d, unsigned nparam,
123 struct domain *domain,
124 unsigned MaxRays)
126 Param_Vertices *vertex;
127 unsigned vbx;
128 int vi, vix;
130 for (vi = 0, vix = 0, vbx = MSB;
131 *vertices;
132 vertices = &(*vertices)->next, ++vi) {
133 int i;
134 for (i = 0; i < words; ++i)
135 if (((*vertices)->Facets[i] & vertex_facets[i]) != vertex_facets[i])
136 break;
137 if (i >= words) {
138 if (!(*vertices)->Domain)
139 domain->F_len = 0;
140 else
141 domain->domain.F[vix] |= vbx;
142 free(vertex_facets);
143 return 0;
145 NEXT(vix, vbx);
147 if (domain->F_len <= vix) {
148 domain->F_len++;
149 domain->domain.F = realloc(domain->domain.F,
150 domain->F_len * sizeof(unsigned));
151 domain->domain.F[domain->F_len-1] = 0;
153 vertex = construct_vertex(vertex_facets, Constraints, d, nparam, MaxRays);
154 if (!vertex->Domain)
155 domain->F_len = 0;
156 else
157 domain->domain.F[vix] |= vbx;
158 vertex->next = *vertices;
159 *vertices = vertex;
160 return 1;
163 static void compute_domain(struct domain *domain, Param_Vertices *vertices,
164 Polyhedron *C, unsigned MaxRays)
166 unsigned bx;
167 int ix, j;
168 int nbV = bit_vector_count(domain->domain.F, domain->F_len);
169 unsigned cols;
170 unsigned rows;
171 Matrix *Constraints;
173 for (ix = 0, bx = MSB; !vertices->Domain; vertices = vertices->next)
174 NEXT(ix, bx);
176 cols = vertices->Domain->NbColumns;
177 rows = vertices->Domain->NbRows;
178 Constraints = Matrix_Alloc(nbV * rows + C->NbConstraints, cols);
180 for (j = 0; vertices; vertices = vertices->next) {
181 if ((domain->domain.F[ix] & bx) == bx)
182 Vector_Copy(vertices->Domain->p[0],
183 Constraints->p[(j++)*rows], rows * cols);
184 NEXT(ix, bx);
186 Vector_Copy(C->Constraint[0], Constraints->p[j*rows], C->NbConstraints * cols);
187 domain->domain.Domain = Constraints2Polyhedron(Constraints, MaxRays);
188 Matrix_Free(Constraints);
191 static void add_domain(struct domain **domains, struct domain *domain,
192 Param_Vertices *vertices, Polyhedron *C,
193 struct barvinok_options *options)
195 options->stats->topcom_chambers++;
197 for (; *domains; domains = (struct domain **)&(*domains)->domain.next) {
198 int i;
199 for (i = 0; i < (*domains)->F_len; ++i)
200 if (((*domains)->domain.F[i] & domain->domain.F[i])
201 != domain->domain.F[i])
202 break;
203 if (i < (*domains)->F_len)
204 continue;
205 for (; i < domain->F_len; ++i)
206 if (domain->domain.F[i])
207 break;
208 if (i >= domain->F_len) {
209 Param_Domain_Free(&domain->domain);
210 return;
213 options->stats->topcom_distinct_chambers++;
214 compute_domain(domain, vertices, C, options->MaxRays);
215 *domains = domain;
218 #define INT_BITS (sizeof(unsigned) * 8)
220 /* Remove any empty or lower-dimensional chamber. The latter
221 * lie on the boundary of the context and are facets of other chambers.
223 * While we are examining each chamber, also extend the F vector
224 * of each chamber to the maximum.
226 static void remove_empty_chambers(Param_Domain **PD, unsigned vertex_words)
228 while (*PD) {
229 int remove = 0;
230 int i;
232 if ((*PD)->Domain->NbEq > 0)
233 remove = 1;
234 else {
235 POL_ENSURE_FACETS((*PD)->Domain);
236 if ((*PD)->Domain->NbEq > 0)
237 remove = 1;
239 if (remove) {
240 Param_Domain *D = *PD;
241 *PD = (*PD)->next;
242 D->next = NULL;
243 Param_Domain_Free(D);
244 continue;
246 if ((i = ((struct domain*)(*PD))->F_len) < vertex_words)
247 (*PD)->F = realloc((*PD)->F, vertex_words * sizeof(unsigned));
248 for (; i < vertex_words; ++i)
249 (*PD)->F[i] = 0;
250 PD = &(*PD)->next;
254 static Param_Polyhedron *points2triangs(Matrix *K, Polyhedron *P, Polyhedron *C,
255 struct barvinok_options *options)
257 int in, out;
258 int i, j;
259 pid_t child;
260 FILE *fin, *fout;
261 int d = K->NbRows;
262 int words = (d+INT_BITS-1)/INT_BITS;
263 Param_Vertices *vertices = NULL;
264 struct domain *domains = NULL;
265 int vertex_words = 1;
266 Param_Polyhedron *PP = ALLOC(Param_Polyhedron);
267 unsigned MaxRays = options->MaxRays;
269 PP->Rays = NULL;
270 PP->nbV = 0;
271 PP->Constraints = Polyhedron2Constraints(P);
272 /* We need the exact facets, because we may make some of them open later */
273 POL_UNSET(options->MaxRays, POL_INTEGER);
275 run_points2triangs(&child, &in, &out);
277 fin = fdopen(in, "w");
278 fprintf(fin, "[\n");
279 for (i = 0; i < K->NbRows; ++i) {
280 fprintf(fin, "[");
281 for (j = 0; j < K->NbColumns; ++j)
282 value_print(fin, P_VALUE_FMT, K->p[i][j]);
283 fprintf(fin, "]");
285 fprintf(fin, "]\n");
286 fclose(fin);
288 fout = fdopen(out, "r");
289 while (fscanf(fout, "T[%d]:=", &i) == 1) {
290 int a, b, c;
291 struct domain *domain = ALLOC(struct domain);
292 memset(domain, 0, sizeof(*domain));
293 domain->domain.F = calloc(vertex_words, sizeof(unsigned));
294 domain->F_len = vertex_words;
296 c = fgetc(fout);
297 if (c == '[')
298 fscanf(fout, "%d,%d:{", &a, &b);
300 while (fgetc(fout) == '{') { /* '{' or closing '}' */
301 int c;
302 unsigned *F = calloc(words, sizeof(unsigned));
304 for (j = 0; j < K->NbColumns; ++j) {
305 unsigned v, shift;
306 fscanf(fout, "%d", &v);
307 shift = INT_BITS - (v % INT_BITS) - 1;
308 F[v / INT_BITS] |= 1u << shift;
309 fgetc(fout); /* , or } */
311 if (!domain->F_len)
312 free(F);
313 else if (add_vertex_to_domain(&vertices, words, F, PP->Constraints,
314 d, C->Dimension,
315 domain, options->MaxRays))
316 ++PP->nbV;
317 if ((c = fgetc(fout)) != ',') /* , or } */
318 ungetc(c, fout);
320 if (domain->F_len)
321 vertex_words = domain->F_len;
322 if (c == '[')
323 fgetc(fout); /* ] */
324 fgetc(fout); /* ; */
325 fgetc(fout); /* \n */
326 if (bit_vector_count(domain->domain.F, domain->F_len) > 0)
327 add_domain(&domains, domain, vertices, C, options);
328 else {
329 options->stats->topcom_empty_chambers++;
330 Param_Domain_Free(&domain->domain);
333 fclose(fout);
335 PP->V = vertices;
336 PP->D = &domains->domain;
338 remove_empty_chambers(&PP->D, vertex_words);
340 options->MaxRays = MaxRays;
342 return PP;
345 /* Assumes M is of full row rank */
346 static Matrix *null_space(Matrix *M)
348 Matrix *H, *Q, *U;
349 Matrix *N;
350 int i;
352 left_hermite(M, &H, &Q, &U);
353 N = Matrix_Alloc(M->NbColumns, M->NbColumns - M->NbRows);
354 for (i = 0; i < N->NbRows; ++i)
355 Vector_Copy(U->p[i] + M->NbRows, N->p[i], N->NbColumns);
356 Matrix_Free(H);
357 Matrix_Free(Q);
358 Matrix_Free(U);
359 return N;
362 static void SwapColumns(Value **V, int n, int i, int j)
364 int r;
366 for (r = 0; r < n; ++r)
367 value_swap(V[r][i], V[r][j]);
370 /* C is assumed to be the "true" context, i.e., it has been intersected
371 * with the projection of P onto the parameter space.
372 * Furthermore, P and C are assumed to be full-dimensional.
374 Param_Polyhedron *TC_P2PP(Polyhedron *P, Polyhedron *C,
375 struct barvinok_options *options)
377 unsigned nparam = C->Dimension;
378 unsigned nvar = P->Dimension - C->Dimension;
379 int i, j;
380 Matrix *H;
381 Matrix *A;
382 Matrix *K;
383 int rows;
384 Param_Polyhedron *PP;
385 Matrix M;
387 assert(C->NbEq == 0);
389 Polyhedron_Matrix_View(P, &M, P->NbConstraints);
390 H = standard_constraints(&M, nparam, &rows, NULL);
392 A = Matrix_Alloc(rows, nvar+rows);
393 for (i = nvar; i < H->NbRows; ++i) {
394 Vector_Oppose(H->p[i], A->p[i-nvar], H->NbColumns);
395 value_set_si(A->p[i-nvar][i], 1);
397 for (i = 0, j = H->NbRows-nvar; i < nvar; ++i) {
398 if (First_Non_Zero(H->p[i], i) == -1)
399 continue;
400 Vector_Oppose(H->p[i], A->p[j], H->NbColumns);
401 value_set_si(A->p[j][j+nvar], 1);
402 SwapColumns(A->p, A->NbRows, j+nvar, i);
403 ++j;
405 K = null_space(A);
406 Matrix_Free(A);
407 /* Ignore extra constraints */
408 K->NbRows = H->NbRows;
409 Matrix_Free(H);
410 PP = points2triangs(K, P, C, options);
411 Matrix_Free(K);
412 return PP;