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
);
39 Param_Vertices vertex
;
48 static struct vertex
*construct_vertex(unsigned *vertex_facets
, Polyhedron
*P
,
49 int d
, unsigned nparam
, unsigned MaxRays
)
51 unsigned nvar
= P
->Dimension
- nparam
;
52 Matrix
*A
= Matrix_Alloc(nvar
+1, nvar
+1);
53 Matrix
*inv
= Matrix_Alloc(nvar
+1, nvar
+1);
54 Matrix
*B
= Matrix_Alloc(nvar
, nparam
+2);
55 Matrix
*V
= Matrix_Alloc(nvar
, nparam
+2);
56 Matrix
*Domain
= Matrix_Alloc(d
-nvar
, nparam
+2);
61 struct vertex
*vertex
;
63 for (j
= 0, i
= 0, ix
= 0, bx
= MSB
; i
< d
; ++i
) {
64 if ((vertex_facets
[ix
] & bx
) == bx
) {
65 Vector_Copy(P
->Constraint
[i
]+1, A
->p
[j
], nvar
);
66 Vector_Oppose(P
->Constraint
[i
]+1+nvar
, B
->p
[j
++], nparam
+1);
71 value_set_si(A
->p
[nvar
][nvar
], 1);
72 ok
= Matrix_Inverse(A
, inv
);
76 inv
->NbColumns
= nvar
;
77 Matrix_Product(inv
, B
, V
);
79 for (i
= 0; i
< nvar
; ++i
) {
80 value_assign(V
->p
[i
][nparam
+1], inv
->p
[nvar
][nvar
]);
81 Vector_Normalize(V
->p
[i
], V
->NbColumns
);
84 for (j
= 0, i
= 0, ix
= 0, bx
= MSB
; i
< d
; ++i
) {
85 if ((vertex_facets
[ix
] & bx
) == bx
) {
89 Param_Inner_Product(P
->Constraint
[i
], V
, Domain
->p
[j
]);
90 if (First_Non_Zero(Domain
->p
[j
]+1, nparam
+1) == -1)
91 vertex_facets
[ix
] |= bx
;
93 value_set_si(Domain
->p
[j
++][0], 1);
97 A
= Matrix_Copy(Domain
);
98 AD
= Constraints2Polyhedron(A
, MaxRays
);
100 POL_ENSURE_VERTICES(AD
);
101 /* A vertex with a lower-dimensional activity domain
102 * saturates more facets than derived above and is actually
103 * the superimposition of two or more vertices.
104 * We therefore discard the domain and (ultimately)
105 * the chamber containing it.
106 * We keep the vertex, though, since it may reappear
107 * in other chambers, which will then likewise be discarded.
108 * The same holds if the activity domain is empty.
115 vertex
= calloc(1, sizeof(struct vertex
));
116 vertex
->facets
= vertex_facets
;
117 vertex
->vertex
.Vertex
= V
;
118 vertex
->vertex
.Domain
= Domain
;
122 static int add_vertex_to_domain(struct vertex
**vertices
, int words
,
123 unsigned *vertex_facets
,
124 Polyhedron
*P
, int d
, unsigned nparam
,
125 struct domain
*domain
,
128 struct vertex
*vertex
;
132 for (vi
= 0, vix
= 0, vbx
= MSB
;
134 vertices
= (struct vertex
**)&(*vertices
)->vertex
.next
, ++vi
) {
136 for (i
= 0; i
< words
; ++i
)
137 if (((*vertices
)->facets
[i
] & vertex_facets
[i
]) != vertex_facets
[i
])
140 if (!(*vertices
)->vertex
.Domain
)
143 domain
->domain
.F
[vix
] |= vbx
;
149 if (domain
->F_len
<= vix
) {
151 domain
->domain
.F
= realloc(domain
->domain
.F
,
152 domain
->F_len
* sizeof(unsigned));
153 domain
->domain
.F
[domain
->F_len
-1] = 0;
155 vertex
= construct_vertex(vertex_facets
, P
, d
, nparam
, MaxRays
);
156 if (!vertex
->vertex
.Domain
)
159 domain
->domain
.F
[vix
] |= vbx
;
160 vertex
->vertex
.next
= &(*vertices
)->vertex
;
165 static int bit_count(unsigned *F
, int F_len
)
170 for (i
= 0; i
< F_len
; ++i
) {
180 static void compute_domain(struct domain
*domain
, struct vertex
*vertices
,
181 Polyhedron
*C
, unsigned MaxRays
)
185 int nbV
= bit_count(domain
->domain
.F
, domain
->F_len
);
186 unsigned cols
= vertices
->vertex
.Domain
->NbColumns
;
187 unsigned rows
= vertices
->vertex
.Domain
->NbRows
;
188 Matrix
*Constraints
= Matrix_Alloc(nbV
* rows
+ C
->NbConstraints
, cols
);
190 for (i
= 0, j
= 0, ix
= 0, bx
= MSB
;
192 vertices
= (struct vertex
*)vertices
->vertex
.next
, ++i
) {
193 if ((domain
->domain
.F
[ix
] & bx
) == bx
)
194 Vector_Copy(vertices
->vertex
.Domain
->p
[0],
195 Constraints
->p
[(j
++)*rows
], rows
* cols
);
198 Vector_Copy(C
->Constraint
[0], Constraints
->p
[j
*rows
], C
->NbConstraints
* cols
);
199 domain
->domain
.Domain
= Constraints2Polyhedron(Constraints
, MaxRays
);
200 Matrix_Free(Constraints
);
203 static void add_domain(struct domain
**domains
, struct domain
*domain
,
204 struct vertex
*vertices
, Polyhedron
*C
,
205 struct barvinok_options
*options
)
207 options
->stats
->topcom_chambers
++;
209 for (; *domains
; domains
= (struct domain
**)&(*domains
)->domain
.next
) {
211 for (i
= 0; i
< (*domains
)->F_len
; ++i
)
212 if (((*domains
)->domain
.F
[i
] & domain
->domain
.F
[i
])
213 != domain
->domain
.F
[i
])
215 if (i
< (*domains
)->F_len
)
217 for (; i
< domain
->F_len
; ++i
)
218 if (domain
->domain
.F
[i
])
220 if (i
>= domain
->F_len
) {
221 Param_Domain_Free(&domain
->domain
);
225 options
->stats
->topcom_distinct_chambers
++;
226 compute_domain(domain
, vertices
, C
, options
->MaxRays
);
230 #define INT_BITS (sizeof(unsigned) * 8)
232 /* Remove any empty or lower-dimensional chamber. The latter
233 * lie on the boundary of the context and are facets of other chambers.
235 * While we are examining each chamber, also extend the F vector
236 * of each chamber to the maximum.
238 static void remove_empty_chambers(Param_Domain
**PD
, unsigned vertex_words
)
244 if ((*PD
)->Domain
->NbEq
> 0)
247 POL_ENSURE_FACETS((*PD
)->Domain
);
248 if ((*PD
)->Domain
->NbEq
> 0)
252 Param_Domain
*D
= *PD
;
255 Param_Domain_Free(D
);
258 if ((i
= ((struct domain
*)(*PD
))->F_len
) < vertex_words
)
259 (*PD
)->F
= realloc((*PD
)->F
, vertex_words
* sizeof(unsigned));
260 for (; i
< vertex_words
; ++i
)
266 /* Clean up memory in struct vertex not in Param_Vertices.
267 * We could also remove the vertices we don't need, but
268 * then we'd have to fix up the vertex masks in the domains.
270 static void clean_up_vertices(Param_Vertices
*V
)
272 for (; V
; V
= V
->next
)
273 free(((struct vertex
*)V
)->facets
);
276 static Param_Polyhedron
*points2triangs(Matrix
*K
, Polyhedron
*P
, Polyhedron
*C
,
277 struct barvinok_options
*options
)
284 int words
= (d
+INT_BITS
-1)/INT_BITS
;
285 struct vertex
*vertices
= NULL
;
286 struct domain
*domains
= NULL
;
287 int vertex_words
= 1;
288 Param_Polyhedron
*PP
= ALLOC(Param_Polyhedron
);
289 unsigned MaxRays
= options
->MaxRays
;
292 /* We need the exact facets, because we may make some of them open later */
293 POL_UNSET(options
->MaxRays
, POL_INTEGER
);
295 run_points2triangs(&child
, &in
, &out
);
297 fin
= fdopen(in
, "w");
299 for (i
= 0; i
< K
->NbRows
; ++i
) {
301 for (j
= 0; j
< K
->NbColumns
; ++j
)
302 value_print(fin
, P_VALUE_FMT
, K
->p
[i
][j
]);
308 fout
= fdopen(out
, "r");
309 while (fscanf(fout
, "T[%d]:={", &i
) == 1) {
310 struct domain
*domain
= ALLOC(struct domain
);
311 memset(domain
, 0, sizeof(*domain
));
312 domain
->domain
.F
= calloc(vertex_words
, sizeof(unsigned));
313 domain
->F_len
= vertex_words
;
315 while (fgetc(fout
) == '{') { /* '{' or closing '}' */
317 unsigned *F
= calloc(words
, sizeof(unsigned));
319 for (j
= 0; j
< K
->NbColumns
; ++j
) {
321 fscanf(fout
, "%d", &v
);
322 shift
= INT_BITS
- (v
% INT_BITS
) - 1;
323 F
[v
/ INT_BITS
] |= 1u << shift
;
324 fgetc(fout
); /* , or } */
328 else if (add_vertex_to_domain(&vertices
, words
, F
, P
, d
, C
->Dimension
,
329 domain
, options
->MaxRays
))
331 if ((c
= fgetc(fout
)) != ',') /* , or } */
335 vertex_words
= domain
->F_len
;
337 fgetc(fout
); /* \n */
338 if (bit_count(domain
->domain
.F
, domain
->F_len
) > 0)
339 add_domain(&domains
, domain
, vertices
, C
, options
);
341 options
->stats
->topcom_empty_chambers
++;
342 Param_Domain_Free(&domain
->domain
);
347 PP
->V
= &vertices
->vertex
;
348 PP
->D
= &domains
->domain
;
350 remove_empty_chambers(&PP
->D
, vertex_words
);
351 clean_up_vertices(PP
->V
);
353 options
->MaxRays
= MaxRays
;
358 /* Assumes M is of full row rank */
359 static Matrix
*null_space(Matrix
*M
)
365 left_hermite(M
, &H
, &Q
, &U
);
366 N
= Matrix_Alloc(M
->NbColumns
, M
->NbColumns
- M
->NbRows
);
367 for (i
= 0; i
< N
->NbRows
; ++i
)
368 Vector_Copy(U
->p
[i
] + M
->NbRows
, N
->p
[i
], N
->NbColumns
);
376 * left_hermite may leave positive entries below the main diagonal in H.
377 * This function postprocesses the output of left_hermite to make
378 * the non-zero entries below the main diagonal negative.
380 static void neg_left_hermite(Matrix
*A
, Matrix
**H_p
, Matrix
**Q_p
, Matrix
**U_p
)
385 left_hermite(A
, &H
, &Q
, &U
);
390 for (row
= 0, col
= 0; col
< H
->NbColumns
; ++col
, ++row
) {
391 while (value_zero_p(H
->p
[row
][col
]))
393 for (i
= 0; i
< col
; ++i
) {
394 if (value_negz_p(H
->p
[row
][i
]))
397 /* subtract column col from column i in H and U */
398 for (j
= 0; j
< H
->NbRows
; ++j
)
399 value_subtract(H
->p
[j
][i
], H
->p
[j
][i
], H
->p
[j
][col
]);
400 for (j
= 0; j
< U
->NbRows
; ++j
)
401 value_subtract(U
->p
[j
][i
], U
->p
[j
][i
], U
->p
[j
][col
]);
403 /* add row i to row col in Q */
404 for (j
= 0; j
< Q
->NbColumns
; ++j
)
405 value_addto(Q
->p
[col
][j
], Q
->p
[col
][j
], Q
->p
[i
][j
]);
410 static void SwapColumns(Value
**V
, int n
, int i
, int j
)
414 for (r
= 0; r
< n
; ++r
)
415 value_swap(V
[r
][i
], V
[r
][j
]);
418 static int is_unit_row(Value
*row
, int pos
, int len
)
420 if (!value_one_p(row
[pos
]) && !value_mone_p(row
[pos
]))
422 return First_Non_Zero(row
+pos
+1, len
-(pos
+1)) == -1;
425 /* C is assumed to be the "true" context, i.e., it has been intersected
426 * with the projection of P onto the parameter space.
427 * Furthermore, P and C are assumed to be full-dimensional.
429 Param_Polyhedron
*TC_P2PP(Polyhedron
*P
, Polyhedron
*C
,
430 struct barvinok_options
*options
)
432 unsigned nparam
= C
->Dimension
;
433 unsigned nvar
= P
->Dimension
- C
->Dimension
;
440 Param_Polyhedron
*PP
;
442 assert(P
->NbEq
== 0);
443 assert(C
->NbEq
== 0);
444 /* move constraints only involving parameters down
445 * and move unit vectors (if there are any) to the right place.
447 for (d
= 0, j
= P
->NbConstraints
; d
< j
; ++d
) {
449 index
= First_Non_Zero(P
->Constraint
[d
]+1, nvar
);
452 is_unit_row(P
->Constraint
[d
]+1, index
, nvar
)) {
453 Vector_Exchange(P
->Constraint
[d
], P
->Constraint
[index
],
460 while (d
< --j
&& First_Non_Zero(P
->Constraint
[j
]+1, nvar
) == -1)
464 Vector_Exchange(P
->Constraint
[d
], P
->Constraint
[j
], P
->Dimension
+2);
466 M
= Matrix_Alloc(d
+nvar
, nvar
);
467 for (j
= 0; j
< d
; ++j
)
468 Vector_Copy(P
->Constraint
[j
]+1, M
->p
[j
], nvar
);
470 neg_left_hermite(M
, &H
, &Q
, &U
);
476 /* Rearrange rows such that top of H is lower diagonal and
477 * add extra unit vector rows that are positive linear
478 * combinations of two or more of the top rows.
480 for (i
= 0; i
< H
->NbColumns
; ++i
) {
481 for (j
= i
; j
< H
->NbRows
; ++j
)
482 if (value_notzero_p(H
->p
[j
][i
]))
485 Vector_Exchange(P
->Constraint
[i
], P
->Constraint
[j
], P
->Dimension
+2);
486 Vector_Exchange(H
->p
[i
], H
->p
[j
], H
->NbColumns
);
488 if (First_Non_Zero(H
->p
[i
], i
) != -1)
489 value_set_si(H
->p
[H
->NbRows
++][i
], 1);
492 rows
= H
->NbRows
-nvar
;
493 A
= Matrix_Alloc(rows
, nvar
+rows
);
494 for (i
= nvar
; i
< d
; ++i
) {
495 Vector_Oppose(H
->p
[i
], A
->p
[i
-nvar
], H
->NbColumns
);
496 value_set_si(A
->p
[i
-nvar
][i
], 1);
498 for (; i
< H
->NbRows
; ++i
) {
499 j
= First_Non_Zero(H
->p
[i
], nvar
);
501 Vector_Oppose(H
->p
[j
], A
->p
[i
-nvar
], H
->NbColumns
);
502 value_set_si(A
->p
[i
-nvar
][i
], 1);
503 SwapColumns(A
->p
, A
->NbRows
, i
, j
);
508 /* Ignore extra constraints */
510 PP
= points2triangs(K
, P
, C
, options
);