configure.in: only compile occ if code_gen library is available
[barvinok.git] / piputil.c
blobb74a3151cf29a3286b31edc1f223534f4eabee0b
1 #include <stdio.h> /* needed for piplib ! */
2 #include <piplib/piplibMP.h>
3 #include <assert.h>
4 #include "piputil.h"
6 #define MAXRAYS POL_NO_DUAL
8 static PipMatrix *poly2pip(Polyhedron *P, int pos, int n, int nparam)
10 PipMatrix * matrix;
11 int i, j;
12 unsigned int extra = P->Dimension - pos - n - nparam;
14 matrix = pip_matrix_alloc(P->NbConstraints, P->Dimension+2);
15 for (i = 0; i < P->NbConstraints; ++i) {
16 value_assign(matrix->p[i][0], P->Constraint[i][0]);
17 for (j = 0; j < n; ++j)
18 value_assign(matrix->p[i][1+j], P->Constraint[i][1+pos+j]);
19 for (j = 0; j < pos; ++j)
20 value_assign(matrix->p[i][1+n+j], P->Constraint[i][1+j]);
21 for (j = 0; j < extra; ++j)
22 value_assign(matrix->p[i][1+n+pos+j], P->Constraint[i][1+n+pos+j]);
23 for (j = 1+pos+n+extra; j < P->Dimension+2; ++j)
24 value_assign(matrix->p[i][j], P->Constraint[i][j]);
27 return matrix;
30 static int max_new(PipQuast *q, int max, int d, int *maxd)
32 PipNewparm *p;
34 for (p = q->newparm; p; p = p->next)
35 if (p->rank > max)
36 max = p->rank;
37 if (q->condition) {
38 if (++d > *maxd)
39 *maxd = d;
40 max = max_new(q->next_else, max, d, maxd);
41 max = max_new(q->next_then, max, d, maxd);
43 return max;
47 * v: total number of variables in original problem
48 * pos: position of minimized variables
49 * n: dimension of space in which minimum was taken
50 * i.e., number of minimized vars in original problem
51 * e: number of extra existential variables
52 * d: total dimension
54 struct scan_data {
55 Matrix *M;
56 int v;
57 int pos;
58 int n;
59 int e;
60 int d;
61 void (*f)(struct scan_data *data, PipQuast *q, int i);
64 static int vectorpos2cpos(struct scan_data *data, int j)
66 int np = data->d - data->v - data->e;
67 int l;
69 if (j < data->pos)
70 l = 1+j;
71 else if (j < data->v - data->n)
72 l = 1+data->n+j;
73 else if (j < data->v-data->n+np)
74 l = 1+data->n+data->e+j;
75 else
76 l = 1+data->n+j-np;
78 return l;
82 * returns true if some solution exist in the whole quast q
84 static int full_quast(PipQuast *q)
86 return q->condition ? full_quast(q->next_then) && full_quast(q->next_else)
87 : q->list != NULL;
91 * i: next row
93 * Dimensions in the resulting domain are:
94 * v n e (d-v-n-e)
96 static void scan_quast_r(struct scan_data *data, PipQuast *q, int i)
98 PipNewparm *p;
99 int skip;
101 /* Skip to a (random) leaf if we are only interested in whether
102 * a solution exists.
104 skip = data->n == 0 && full_quast(q);
105 while (skip && q->condition)
106 q = q->next_then;
108 for (p = q->newparm; !skip && p; p = p->next) {
109 int j, l;
110 PipVector *vec = p->vector;
112 value_set_si(data->M->p[i][0], 1);
113 Vector_Set(data->M->p[i]+1, 0, data->d);
115 for (j = 0; j < vec->nb_elements-1; ++j) {
116 l = vectorpos2cpos(data, j);
117 value_assign(data->M->p[i][l], vec->the_vector[j]);
119 l = vectorpos2cpos(data, p->rank);
120 value_oppose(data->M->p[i][l], p->deno);
121 value_assign(data->M->p[i][1+data->d], vec->the_vector[j]);
122 ++i;
124 value_set_si(data->M->p[i][0], 1);
125 Vector_Set(data->M->p[i]+1, 0, data->d);
127 for (j = 0; j < vec->nb_elements-1; ++j) {
128 l = vectorpos2cpos(data, j);
129 value_oppose(data->M->p[i][l], vec->the_vector[j]);
131 l = vectorpos2cpos(data, p->rank);
132 value_assign(data->M->p[i][l], p->deno);
133 value_oppose(data->M->p[i][1+data->d], vec->the_vector[j]);
134 value_addto(data->M->p[i][1+data->d], data->M->p[i][1+data->d], p->deno);
135 value_decrement(data->M->p[i][1+data->d], data->M->p[i][1+data->d]);
136 ++i;
139 if (q->condition) {
140 int j;
141 value_set_si(data->M->p[i][0], 1);
142 Vector_Set(data->M->p[i]+1, 0, data->d);
144 for (j = 0; j < q->condition->nb_elements-1; ++j) {
145 int l = vectorpos2cpos(data, j);
146 value_assign(data->M->p[i][l], q->condition->the_vector[j]);
148 value_assign(data->M->p[i][1+data->d], q->condition->the_vector[j]);
149 scan_quast_r(data, q->next_then, i+1);
151 for (j = 0; j < q->condition->nb_elements-1; ++j) {
152 int l = vectorpos2cpos(data, j);
153 value_oppose(data->M->p[i][l], q->condition->the_vector[j]);
155 value_oppose(data->M->p[i][1+data->d], q->condition->the_vector[j]);
156 value_decrement(data->M->p[i][1+data->d], data->M->p[i][1+data->d]);
157 scan_quast_r(data, q->next_else, i+1);
159 return;
162 /* if data->n is zero, we are only interested in the domains
163 * where a solution exists and not in the actual solution
165 if (q->list && data->n) {
166 PipList *l;
167 int j, k;
168 for (j = 0, l = q->list; l; ++j, l = l->next) {
169 Vector_Set(data->M->p[i], 0, data->d+1);
170 value_set_si(data->M->p[i][1+data->pos+j], -1);
172 for (k = 0; k < l->vector->nb_elements-1; ++k) {
173 int ll = vectorpos2cpos(data, k);
174 value_assign(data->M->p[i][ll], l->vector->the_vector[k]);
176 value_assign(data->M->p[i][1+data->d], l->vector->the_vector[k]);
178 ++i;
181 data->f(data, q, i);
184 static void scan_quast(struct scan_data *data, PipQuast *q)
186 int nparam, nexist, d, i, rows;
188 nparam = data->d - data->n;
190 d = 0;
191 nexist = max_new(q, nparam-1, 0, &d) - nparam+1;
192 rows = 2 * nexist + d + data->n;
194 /* nparam now refers to the number of parameters in the original polyhedron */
195 nparam -= data->v - data->n;
197 data->e = nexist;
198 data->d = data->v + nexist + nparam;
200 data->M = Matrix_Alloc(rows, 1+data->d+1);
202 scan_quast_r(data, q, 0);
203 Matrix_Free(data->M);
206 struct quast2poly_data {
207 struct scan_data scan;
208 Polyhedron* D;
211 static void add_quast_base(struct scan_data *data, PipQuast *q, int i)
213 struct quast2poly_data *qdata = (struct quast2poly_data *)data;
214 if (q->list) {
215 Matrix *C;
216 Polyhedron *P;
217 C = Matrix_Alloc(i, 1+data->d+1);
218 Vector_Copy(data->M->p[0], C->p[0], i * (1+data->d+1));
219 P = Constraints2Polyhedron(C, MAXRAYS);
220 Matrix_Free(C);
221 P->next = qdata->D;
222 qdata->D = P;
227 * nvar: total number of variables, including the minimized variables,
228 * but excluding the parameters
229 * nparam: the number of parameters in the PIP problem, excluding the "big parameter"
230 * for maximization problems, i.e., the
231 * total number of variables minus the minimized variables
232 * pos: position of the first minimized variable
233 * n: number of minimized variables
235 static Polyhedron *quast2poly(PipQuast *q, int nvar, int nparam, int pos, int n)
237 struct quast2poly_data data;
239 data.scan.v = nvar;
240 data.scan.pos = pos;
241 data.scan.n = n;
242 data.scan.d = n+nparam;
244 data.D = 0;
245 data.scan.f = add_quast_base;
247 scan_quast(&data.scan, q);
249 return data.D;
252 Polyhedron *pip_projectout(Polyhedron *P, int pos, int n, int nparam)
254 PipOptions *options;
255 PipMatrix * domain;
256 unsigned int nvar = P->Dimension - nparam;
257 PipMatrix *context = pip_matrix_alloc(0, P->Dimension - n + 2);
258 PipQuast *sol;
259 Polyhedron *min;
261 POL_ENSURE_INEQUALITIES(P);
263 domain = poly2pip(P, pos, n, nparam);
264 #ifdef PIP_DEBUG
265 pip_matrix_print(stderr, domain);
266 pip_matrix_print(stderr, context);
267 #endif
269 options = pip_options_init();
270 options->Urs_unknowns = -1;
271 options->Urs_parms = -1;
272 options->Simplify = 1;
273 sol = pip_solve(domain, context, -1, options);
275 #ifdef PIP_DEBUG
276 pip_quast_print(stderr, sol, 0);
277 #endif
279 min = quast2poly(sol, nvar - n, P->Dimension - n, 0, 0);
281 pip_quast_free(sol);
282 pip_matrix_free(context);
283 pip_matrix_free(domain);
284 pip_options_free(options);
286 return min;