1 #include <barvinok/util.h>
2 #include <barvinok/options.h>
7 #define ALLOC(type) (type*)malloc(sizeof(type))
9 void run_points2triangs(pid_t
*child
, int *in
, int *out
)
11 int in_fd
[2], out_fd
[2];
29 rc
= execl(POINTS2TRIANGS_PATH
, "points2triangs", "--regular", NULL
);
43 static Param_Vertices
*construct_vertex(unsigned *vertex_facets
,
45 int d
, unsigned nparam
, unsigned MaxRays
)
47 unsigned nvar
= Constraints
->NbColumns
-2 - nparam
;
48 Matrix
*A
= Matrix_Alloc(nvar
+1, nvar
+1);
49 Matrix
*inv
= Matrix_Alloc(nvar
+1, nvar
+1);
50 Matrix
*B
= Matrix_Alloc(nvar
, nparam
+2);
51 Matrix
*V
= Matrix_Alloc(nvar
, nparam
+2);
52 Matrix
*Domain
= Matrix_Alloc(d
-nvar
, nparam
+2);
57 Param_Vertices
*vertex
;
59 for (j
= 0, i
= 0, ix
= 0, bx
= MSB
; i
< d
; ++i
) {
60 if ((vertex_facets
[ix
] & bx
) == bx
) {
61 Vector_Copy(Constraints
->p
[i
]+1, A
->p
[j
], nvar
);
62 Vector_Oppose(Constraints
->p
[i
]+1+nvar
, B
->p
[j
++], nparam
+1);
67 value_set_si(A
->p
[nvar
][nvar
], 1);
68 ok
= Matrix_Inverse(A
, inv
);
72 inv
->NbColumns
= nvar
;
73 Matrix_Product(inv
, B
, V
);
75 for (i
= 0; i
< nvar
; ++i
) {
76 value_assign(V
->p
[i
][nparam
+1], inv
->p
[nvar
][nvar
]);
77 Vector_Normalize(V
->p
[i
], V
->NbColumns
);
80 for (j
= 0, i
= 0, ix
= 0, bx
= MSB
; i
< d
; ++i
) {
81 if ((vertex_facets
[ix
] & bx
) == bx
) {
85 Param_Inner_Product(Constraints
->p
[i
], V
, Domain
->p
[j
]);
86 if (First_Non_Zero(Domain
->p
[j
]+1, nparam
+1) == -1)
87 vertex_facets
[ix
] |= bx
;
89 value_set_si(Domain
->p
[j
++][0], 1);
93 A
= Matrix_Copy(Domain
);
94 AD
= Constraints2Polyhedron(A
, MaxRays
);
96 POL_ENSURE_VERTICES(AD
);
97 /* A vertex with a lower-dimensional activity domain
98 * saturates more facets than derived above and is actually
99 * the superimposition of two or more vertices.
100 * We therefore discard the domain and (ultimately)
101 * the chamber containing it.
102 * We keep the vertex, though, since it may reappear
103 * in other chambers, which will then likewise be discarded.
104 * The same holds if the activity domain is empty.
111 vertex
= calloc(1, sizeof(Param_Vertices
));
112 vertex
->Facets
= vertex_facets
;
114 vertex
->Domain
= Domain
;
118 static int add_vertex_to_domain(Param_Vertices
**vertices
, int words
,
119 unsigned *vertex_facets
,
120 Matrix
*Constraints
, int d
, unsigned nparam
,
121 struct domain
*domain
,
124 Param_Vertices
*vertex
;
128 for (vi
= 0, vix
= 0, vbx
= MSB
;
130 vertices
= &(*vertices
)->next
, ++vi
) {
132 for (i
= 0; i
< words
; ++i
)
133 if (((*vertices
)->Facets
[i
] & vertex_facets
[i
]) != vertex_facets
[i
])
136 if (!(*vertices
)->Domain
)
139 domain
->domain
.F
[vix
] |= vbx
;
145 if (domain
->F_len
<= vix
) {
147 domain
->domain
.F
= realloc(domain
->domain
.F
,
148 domain
->F_len
* sizeof(unsigned));
149 domain
->domain
.F
[domain
->F_len
-1] = 0;
151 vertex
= construct_vertex(vertex_facets
, Constraints
, d
, nparam
, MaxRays
);
155 domain
->domain
.F
[vix
] |= vbx
;
156 vertex
->next
= *vertices
;
161 static int bit_count(unsigned *F
, int F_len
)
166 for (i
= 0; i
< F_len
; ++i
) {
176 static void compute_domain(struct domain
*domain
, Param_Vertices
*vertices
,
177 Polyhedron
*C
, unsigned MaxRays
)
181 int nbV
= bit_count(domain
->domain
.F
, domain
->F_len
);
182 unsigned cols
= vertices
->Domain
->NbColumns
;
183 unsigned rows
= vertices
->Domain
->NbRows
;
184 Matrix
*Constraints
= Matrix_Alloc(nbV
* rows
+ C
->NbConstraints
, cols
);
186 for (i
= 0, j
= 0, ix
= 0, bx
= MSB
;
188 vertices
= vertices
->next
, ++i
) {
189 if ((domain
->domain
.F
[ix
] & bx
) == bx
)
190 Vector_Copy(vertices
->Domain
->p
[0],
191 Constraints
->p
[(j
++)*rows
], rows
* cols
);
194 Vector_Copy(C
->Constraint
[0], Constraints
->p
[j
*rows
], C
->NbConstraints
* cols
);
195 domain
->domain
.Domain
= Constraints2Polyhedron(Constraints
, MaxRays
);
196 Matrix_Free(Constraints
);
199 static void add_domain(struct domain
**domains
, struct domain
*domain
,
200 Param_Vertices
*vertices
, Polyhedron
*C
,
201 struct barvinok_options
*options
)
203 options
->stats
->topcom_chambers
++;
205 for (; *domains
; domains
= (struct domain
**)&(*domains
)->domain
.next
) {
207 for (i
= 0; i
< (*domains
)->F_len
; ++i
)
208 if (((*domains
)->domain
.F
[i
] & domain
->domain
.F
[i
])
209 != domain
->domain
.F
[i
])
211 if (i
< (*domains
)->F_len
)
213 for (; i
< domain
->F_len
; ++i
)
214 if (domain
->domain
.F
[i
])
216 if (i
>= domain
->F_len
) {
217 Param_Domain_Free(&domain
->domain
);
221 options
->stats
->topcom_distinct_chambers
++;
222 compute_domain(domain
, vertices
, C
, options
->MaxRays
);
226 #define INT_BITS (sizeof(unsigned) * 8)
228 /* Remove any empty or lower-dimensional chamber. The latter
229 * lie on the boundary of the context and are facets of other chambers.
231 * While we are examining each chamber, also extend the F vector
232 * of each chamber to the maximum.
234 static void remove_empty_chambers(Param_Domain
**PD
, unsigned vertex_words
)
240 if ((*PD
)->Domain
->NbEq
> 0)
243 POL_ENSURE_FACETS((*PD
)->Domain
);
244 if ((*PD
)->Domain
->NbEq
> 0)
248 Param_Domain
*D
= *PD
;
251 Param_Domain_Free(D
);
254 if ((i
= ((struct domain
*)(*PD
))->F_len
) < vertex_words
)
255 (*PD
)->F
= realloc((*PD
)->F
, vertex_words
* sizeof(unsigned));
256 for (; i
< vertex_words
; ++i
)
262 static Param_Polyhedron
*points2triangs(Matrix
*K
, Polyhedron
*P
, Polyhedron
*C
,
263 struct barvinok_options
*options
)
270 int words
= (d
+INT_BITS
-1)/INT_BITS
;
271 Param_Vertices
*vertices
= NULL
;
272 struct domain
*domains
= NULL
;
273 int vertex_words
= 1;
274 Param_Polyhedron
*PP
= ALLOC(Param_Polyhedron
);
275 unsigned MaxRays
= options
->MaxRays
;
278 PP
->Constraints
= Polyhedron2Constraints(P
);
279 /* We need the exact facets, because we may make some of them open later */
280 POL_UNSET(options
->MaxRays
, POL_INTEGER
);
282 run_points2triangs(&child
, &in
, &out
);
284 fin
= fdopen(in
, "w");
286 for (i
= 0; i
< K
->NbRows
; ++i
) {
288 for (j
= 0; j
< K
->NbColumns
; ++j
)
289 value_print(fin
, P_VALUE_FMT
, K
->p
[i
][j
]);
295 fout
= fdopen(out
, "r");
296 while (fscanf(fout
, "T[%d]:={", &i
) == 1) {
297 struct domain
*domain
= ALLOC(struct domain
);
298 memset(domain
, 0, sizeof(*domain
));
299 domain
->domain
.F
= calloc(vertex_words
, sizeof(unsigned));
300 domain
->F_len
= vertex_words
;
302 while (fgetc(fout
) == '{') { /* '{' or closing '}' */
304 unsigned *F
= calloc(words
, sizeof(unsigned));
306 for (j
= 0; j
< K
->NbColumns
; ++j
) {
308 fscanf(fout
, "%d", &v
);
309 shift
= INT_BITS
- (v
% INT_BITS
) - 1;
310 F
[v
/ INT_BITS
] |= 1u << shift
;
311 fgetc(fout
); /* , or } */
315 else if (add_vertex_to_domain(&vertices
, words
, F
, PP
->Constraints
,
317 domain
, options
->MaxRays
))
319 if ((c
= fgetc(fout
)) != ',') /* , or } */
323 vertex_words
= domain
->F_len
;
325 fgetc(fout
); /* \n */
326 if (bit_count(domain
->domain
.F
, domain
->F_len
) > 0)
327 add_domain(&domains
, domain
, vertices
, C
, options
);
329 options
->stats
->topcom_empty_chambers
++;
330 Param_Domain_Free(&domain
->domain
);
336 PP
->D
= &domains
->domain
;
338 remove_empty_chambers(&PP
->D
, vertex_words
);
340 options
->MaxRays
= MaxRays
;
345 /* Assumes M is of full row rank */
346 static Matrix
*null_space(Matrix
*M
)
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
);
362 static void SwapColumns(Value
**V
, int n
, int i
, int j
)
366 for (r
= 0; r
< n
; ++r
)
367 value_swap(V
[r
][i
], V
[r
][j
]);
370 static int is_unit_row(Value
*row
, int pos
, int len
)
372 if (!value_one_p(row
[pos
]) && !value_mone_p(row
[pos
]))
374 return First_Non_Zero(row
+pos
+1, len
-(pos
+1)) == -1;
377 /* C is assumed to be the "true" context, i.e., it has been intersected
378 * with the projection of P onto the parameter space.
379 * Furthermore, P and C are assumed to be full-dimensional.
381 Param_Polyhedron
*TC_P2PP(Polyhedron
*P
, Polyhedron
*C
,
382 struct barvinok_options
*options
)
384 unsigned nparam
= C
->Dimension
;
385 unsigned nvar
= P
->Dimension
- C
->Dimension
;
392 Param_Polyhedron
*PP
;
394 assert(P
->NbEq
== 0);
395 assert(C
->NbEq
== 0);
396 /* move constraints only involving parameters down
397 * and move unit vectors (if there are any) to the right place.
399 for (d
= 0, j
= P
->NbConstraints
; d
< j
; ++d
) {
401 index
= First_Non_Zero(P
->Constraint
[d
]+1, nvar
);
404 is_unit_row(P
->Constraint
[d
]+1, index
, nvar
)) {
405 Vector_Exchange(P
->Constraint
[d
], P
->Constraint
[index
],
412 while (d
< --j
&& First_Non_Zero(P
->Constraint
[j
]+1, nvar
) == -1)
416 Vector_Exchange(P
->Constraint
[d
], P
->Constraint
[j
], P
->Dimension
+2);
418 M
= Matrix_Alloc(d
+nvar
, nvar
);
419 for (j
= 0; j
< d
; ++j
)
420 Vector_Copy(P
->Constraint
[j
]+1, M
->p
[j
], nvar
);
422 neg_left_hermite(M
, &H
, &Q
, &U
);
428 /* Rearrange rows such that top of H is lower diagonal and
429 * add extra unit vector rows that are positive linear
430 * combinations of two or more of the top rows.
432 for (i
= 0; i
< H
->NbColumns
; ++i
) {
433 for (j
= i
; j
< H
->NbRows
; ++j
)
434 if (value_notzero_p(H
->p
[j
][i
]))
437 Vector_Exchange(P
->Constraint
[i
], P
->Constraint
[j
], P
->Dimension
+2);
438 Vector_Exchange(H
->p
[i
], H
->p
[j
], H
->NbColumns
);
440 if (First_Non_Zero(H
->p
[i
], i
) != -1)
441 value_set_si(H
->p
[H
->NbRows
++][i
], 1);
444 rows
= H
->NbRows
-nvar
;
445 A
= Matrix_Alloc(rows
, nvar
+rows
);
446 for (i
= nvar
; i
< d
; ++i
) {
447 Vector_Oppose(H
->p
[i
], A
->p
[i
-nvar
], H
->NbColumns
);
448 value_set_si(A
->p
[i
-nvar
][i
], 1);
450 for (; i
< H
->NbRows
; ++i
) {
451 j
= First_Non_Zero(H
->p
[i
], nvar
);
453 Vector_Oppose(H
->p
[j
], A
->p
[i
-nvar
], H
->NbColumns
);
454 value_set_si(A
->p
[i
-nvar
][i
], 1);
455 SwapColumns(A
->p
, A
->NbRows
, i
, j
);
460 /* Ignore extra constraints */
462 PP
= points2triangs(K
, P
, C
, options
);