summate.c: sum_with_equalities: extract out sum_over_polytope_0D
[barvinok.git] / test.c
blob62d108357f6dcf2cc30c7ce82c0e8a1385dfd3f4
1 #include <assert.h>
2 #include <stdio.h>
3 #include <string.h>
4 #include <unistd.h>
5 #include <isl_set_polylib.h>
6 #include <barvinok/util.h>
7 #include <barvinok/barvinok.h>
8 #include <barvinok/sample.h>
9 #include "config.h"
10 #include "evalue_read.h"
11 #include "lattice_width.h"
13 #ifdef HAVE_SYS_TIMES_H
15 #include <sys/times.h>
17 static void time_diff(struct tms *before, struct tms *after)
19 long ticks = sysconf(_SC_CLK_TCK);
20 printf("User: %g; Sys: %g\n",
21 (0.0 + after->tms_utime - before->tms_utime) / ticks,
22 (0.0 + after->tms_stime - before->tms_stime) / ticks);
25 #else
27 struct tms {};
28 static void times(struct tms* time)
31 static void time_diff(struct tms *before, struct tms *after)
35 #endif
37 int main(int argc, char **argv)
39 isl_ctx *ctx;
40 int i, nbPol, nbVec, nbMat, func, j, n;
41 Polyhedron *A, *B, *C, *D, *E, *F, *G;
42 char s[128];
43 struct barvinok_options *options = barvinok_options_new_with_defaults();
45 argc = barvinok_options_parse(options, argc, argv, ISL_ARG_ALL);
46 ctx = isl_ctx_alloc_with_options(&barvinok_options_args, options);
48 nbPol = nbVec = nbMat = 0;
49 fgets(s, 128, stdin);
50 while ((*s=='#') ||
51 ((sscanf(s, "D %d", &nbPol) < 1) &&
52 (sscanf(s, "V %d", &nbVec) < 1) &&
53 (sscanf(s, "M %d", &nbMat) < 1)))
54 fgets(s, 128, stdin);
56 for (i = 0; i < nbPol; ++i) {
57 Matrix *M = Matrix_Read();
58 A = Constraints2Polyhedron(M, options->MaxRays);
59 Matrix_Free(M);
60 fgets(s, 128, stdin);
61 while ((*s=='#') || (sscanf(s, "F %d", &func)<1))
62 fgets(s, 128, stdin);
64 switch(func) {
65 case 0: {
66 Value cb, ck;
67 value_init(cb);
68 value_init(ck);
69 fgets(s, 128, stdin);
70 /* workaround for apparent bug in older gmps */
71 *strchr(s, '\n') = '\0';
72 while ((*s=='#') || (value_read(ck, s) != 0)) {
73 fgets(s, 128, stdin);
74 /* workaround for apparent bug in older gmps */
75 *strchr(s, '\n') = '\0';
77 barvinok_count_with_options(A, &cb, options);
78 if (value_ne(cb, ck))
79 return -1;
80 value_clear(cb);
81 value_clear(ck);
82 break;
84 case 1:
85 Polyhedron_Print(stdout, P_VALUE_FMT, A);
86 B = Polyhedron_Polar(A, options->MaxRays);
87 Polyhedron_Print(stdout, P_VALUE_FMT, B);
88 C = Polyhedron_Polar(B, options->MaxRays);
89 Polyhedron_Print(stdout, P_VALUE_FMT, C);
90 Polyhedron_Free(C);
91 Polyhedron_Free(B);
92 break;
93 case 2:
94 Polyhedron_Print(stdout, P_VALUE_FMT, A);
95 for (j = 0; j < A->NbRays; ++j) {
96 B = supporting_cone(A, j);
97 Polyhedron_Print(stdout, P_VALUE_FMT, B);
98 Polyhedron_Free(B);
100 break;
101 case 3:
102 Polyhedron_Print(stdout, P_VALUE_FMT, A);
103 C = B = NULL;
104 barvinok_decompose(A,&B,&C);
105 puts("Pos:");
106 Polyhedron_Print(stdout, P_VALUE_FMT, B);
107 puts("Neg:");
108 Polyhedron_Print(stdout, P_VALUE_FMT, C);
109 Domain_Free(B);
110 Domain_Free(C);
111 break;
112 case 4: {
113 Value cm, cb;
114 struct tms tms_before, tms_between, tms_after;
115 value_init(cm);
116 value_init(cb);
117 Polyhedron_Print(stdout, P_VALUE_FMT, A);
118 times(&tms_before);
119 manual_count(A, &cm);
120 times(&tms_between);
121 barvinok_count(A, &cb, 100);
122 times(&tms_after);
123 printf("manual: ");
124 value_print(stdout, P_VALUE_FMT, cm);
125 puts("");
126 time_diff(&tms_before, &tms_between);
127 printf("Barvinok: ");
128 value_print(stdout, P_VALUE_FMT, cb);
129 puts("");
130 time_diff(&tms_between, &tms_after);
131 value_clear(cm);
132 value_clear(cb);
133 break;
135 case 5:
136 Polyhedron_Print(stdout, P_VALUE_FMT, A);
137 B = triangulate_cone(A, 100);
138 Polyhedron_Print(stdout, P_VALUE_FMT, B);
139 check_triangulization(A, B);
140 Domain_Free(B);
141 break;
142 case 6:
143 Polyhedron_Print(stdout, P_VALUE_FMT, A);
144 B = remove_equalities(A, options->MaxRays);
145 Polyhedron_Print(stdout, P_VALUE_FMT, B);
146 Polyhedron_Free(B);
147 break;
148 case 8: {
149 evalue *EP;
150 Matrix *M = Matrix_Read();
151 const char **param_name;
152 C = Constraints2Polyhedron(M, options->MaxRays);
153 Matrix_Free(M);
154 Polyhedron_Print(stdout, P_VALUE_FMT, A);
155 Polyhedron_Print(stdout, P_VALUE_FMT, C);
156 EP = barvinok_enumerate_with_options(A, C, options);
157 param_name = Read_ParamNames(stdin, C->Dimension);
158 print_evalue(stdout, EP, (const char**)param_name);
159 evalue_free(EP);
160 Polyhedron_Free(C);
162 case 9:
163 Polyhedron_Print(stdout, P_VALUE_FMT, A);
164 Polyhedron_Polarize(A);
165 C = B = NULL;
166 barvinok_decompose(A,&B,&C);
167 for (D = B; D; D = D->next)
168 Polyhedron_Polarize(D);
169 for (D = C; D; D = D->next)
170 Polyhedron_Polarize(D);
171 puts("Pos:");
172 Polyhedron_Print(stdout, P_VALUE_FMT, B);
173 puts("Neg:");
174 Polyhedron_Print(stdout, P_VALUE_FMT, C);
175 Domain_Free(B);
176 Domain_Free(C);
177 break;
178 case 10: {
179 evalue *EP;
180 Value cb, ck;
182 value_init(cb);
183 value_init(ck);
184 fgets(s, 128, stdin);
185 sscanf(s, "%d", &n);
186 for (j = 0; j < n; ++j) {
187 Polyhedron *P;
188 M = Matrix_Read();
189 P = Constraints2Polyhedron(M, options->MaxRays);
190 Matrix_Free(M);
191 A = DomainConcat(P, A);
193 fgets(s, 128, stdin);
194 /* workaround for apparent bug in older gmps */
195 *strchr(s, '\n') = '\0';
196 while ((*s=='#') || (value_read(ck, s) != 0)) {
197 fgets(s, 128, stdin);
198 /* workaround for apparent bug in older gmps */
199 *strchr(s, '\n') = '\0';
201 C = Universe_Polyhedron(0);
202 EP = barvinok_enumerate_union(A, C, options->MaxRays);
203 value_set_double(cb, compute_evalue(EP, &ck)+.25);
204 if (value_ne(cb, ck))
205 return -1;
206 Domain_Free(C);
207 value_clear(cb);
208 value_clear(ck);
209 evalue_free(EP);
210 break;
212 case 11: {
213 isl_space *dim;
214 isl_basic_set *bset;
215 isl_pw_qpolynomial *expected, *computed;
216 unsigned nparam;
218 expected = isl_pw_qpolynomial_read_from_file(ctx, stdin);
219 nparam = isl_pw_qpolynomial_dim(expected, isl_dim_param);
220 dim = isl_space_set_alloc(ctx, nparam, A->Dimension - nparam);
221 bset = isl_basic_set_new_from_polylib(A, dim);
222 computed = isl_basic_set_lattice_width(bset);
223 computed = isl_pw_qpolynomial_sub(computed, expected);
224 if (!isl_pw_qpolynomial_is_zero(computed))
225 return -1;
226 isl_pw_qpolynomial_free(computed);
227 break;
229 case 12: {
230 Vector *sample;
231 int has_sample;
232 fgets(s, 128, stdin);
233 sscanf(s, "%d", &has_sample);
235 sample = Polyhedron_Sample(A, options);
236 if (!sample && has_sample)
237 return -1;
238 if (sample && !has_sample)
239 return -1;
240 if (sample && !in_domain(A, sample->p))
241 return -1;
242 Vector_Free(sample);
245 Domain_Free(A);
247 for (i = 0; i < nbVec; ++i) {
248 int ok;
249 Vector *V = Vector_Read();
250 Matrix *M = Matrix_Alloc(V->Size, V->Size);
251 Vector_Copy(V->p, M->p[0], V->Size);
252 ok = unimodular_complete(M, 1);
253 assert(ok);
254 Matrix_Print(stdout, P_VALUE_FMT, M);
255 Matrix_Free(M);
256 Vector_Free(V);
258 for (i = 0; i < nbMat; ++i) {
259 Matrix *U, *V, *S;
260 Matrix *M = Matrix_Read();
261 Smith(M, &U, &V, &S);
262 Matrix_Print(stdout, P_VALUE_FMT, U);
263 Matrix_Print(stdout, P_VALUE_FMT, V);
264 Matrix_Print(stdout, P_VALUE_FMT, S);
265 Matrix_Free(M);
266 Matrix_Free(U);
267 Matrix_Free(V);
268 Matrix_Free(S);
271 isl_ctx_free(ctx);
272 return 0;