lexmin.cc: move some code around to prepare for new version
[barvinok.git] / polytope_scan.c
blob5b3cf5bcfdc887aae5159dc950c4eac9dcce665f
1 #include <unistd.h>
2 #include <stdlib.h>
3 #include <strings.h>
4 #include <polylib/polylibgmp.h>
5 #include <barvinok/util.h>
6 #include "basis_reduction.h"
7 #include "config.h"
9 #ifdef HAVE_GROWING_CHERNIKOVA
10 #define MAXRAYS (POL_NO_DUAL | POL_INTEGER)
11 #else
12 #define MAXRAYS 600
13 #endif
15 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
17 #ifndef HAVE_GETOPT_H
18 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
19 #else
20 #include <getopt.h>
21 struct option options[] = {
22 { "direct", no_argument, 0, 'd' },
23 { "version", no_argument, 0, 'V' },
24 { 0, 0, 0, 0 }
26 #endif
28 static Polyhedron *Polyhedron_Read()
30 int vertices = 0;
31 unsigned NbRows, NbColumns;
32 Matrix *M;
33 Polyhedron *P;
34 char s[128];
36 while (fgets(s, sizeof(s), stdin)) {
37 if (*s == '#')
38 continue;
39 if (strncasecmp(s, "vertices", sizeof("vertices")-1) == 0)
40 vertices = 1;
41 if (sscanf(s, "%u %u", &NbRows, &NbColumns) == 2)
42 break;
44 if (feof(stdin))
45 return NULL;
46 M = Matrix_Alloc(NbRows,NbColumns);
47 Matrix_Read_Input(M);
48 if (vertices)
49 P = Rays2Polyhedron(M, MAXRAYS);
50 else
51 P = Constraints2Polyhedron(M, MAXRAYS);
52 Matrix_Free(M);
53 return P;
56 static void scan_poly(Polyhedron *S, int pos, Value *z, Matrix *T)
58 if (!S) {
59 int k;
60 Vector *v;
62 v = Vector_Alloc(T->NbRows);
63 Matrix_Vector_Product(T, z+1, v->p);
64 value_print(stdout, VALUE_FMT, v->p[0]);
65 for (k=1; k < pos; ++k) {
66 printf(", ");
67 value_print(stdout,VALUE_FMT, v->p[k]);
69 Vector_Free(v);
70 printf("\n");
71 } else {
72 int ok;
73 Value LB, UB, tmp;
74 value_init(LB);
75 value_init(UB);
76 value_init(tmp);
77 ok = !(lower_upper_bounds(1+pos, S, z, &LB, &UB));
78 assert(ok);
79 for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
80 value_assign(z[pos+1], tmp);
81 scan_poly(S->next, pos+1, z, T);
83 value_set_si(z[pos+1], 0);
84 value_clear(LB);
85 value_clear(UB);
86 value_clear(tmp);
90 int main(int argc, char **argv)
92 Polyhedron *A, *P, *U, *S;
93 Value *p;
94 int i, j, ok;
95 Matrix *basis, *T, *inv;
96 int c, ind = 0;
97 int direct = 0;
99 while ((c = getopt_long(argc, argv, "dV", options, &ind)) != -1) {
100 switch (c) {
101 case 'd':
102 direct = 1;
103 break;
104 case 'V':
105 printf(barvinok_version());
106 exit(0);
107 break;
111 A = Polyhedron_Read();
113 if (direct) {
114 inv = Identity(A->Dimension+1);
115 P = A;
116 } else {
117 basis = reduced_basis(A);
119 T = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
120 inv = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
121 for (i = 0; i < A->Dimension; ++i)
122 for (j = 0; j < A->Dimension; ++j)
123 value_assign(T->p[i][j], basis->p[i][j]);
124 value_set_si(T->p[A->Dimension][A->Dimension], 1);
125 Matrix_Free(basis);
127 ok = Matrix_Inverse(T, inv);
128 assert(ok);
129 Matrix_Free(T);
131 P = Polyhedron_Preimage(A, inv, MAXRAYS);
132 Polyhedron_Free(A);
135 U = Universe_Polyhedron(0);
136 S = Polyhedron_Scan(P, U, MAXRAYS);
138 p = ALLOCN(Value, P->Dimension+2);
139 for(i=0;i<=P->Dimension;i++) {
140 value_init(p[i]);
141 value_set_si(p[i],0);
143 value_init(p[i]);
144 value_set_si(p[i],1);
146 scan_poly(S, 0, p, inv);
148 Matrix_Free(inv);
149 for(i=0;i<=(P->Dimension+1);i++)
150 value_clear(p[i]);
151 free(p);
152 Domain_Free(S);
153 Polyhedron_Free(P);
154 Polyhedron_Free(U);
155 return 0;