add isl_basic_set_interval and isl_basic_set_product
[platform/upstream/isl.git] / isl_sample_piplib.c
1 #include "isl_mat.h"
2 #include "isl_vec.h"
3 #include "isl_seq.h"
4 #include "isl_piplib.h"
5 #include "isl_sample_piplib.h"
6
7 static void swap_inequality(struct isl_basic_set *bset, int a, int b)
8 {
9         isl_int *t = bset->ineq[a];
10         bset->ineq[a] = bset->ineq[b];
11         bset->ineq[b] = t;
12 }
13
14 static struct isl_mat *independent_bounds(struct isl_ctx *ctx,
15         struct isl_basic_set *bset)
16 {
17         int i, j, n;
18         struct isl_mat *dirs = NULL;
19         struct isl_mat *bounds = NULL;
20
21         if (!bset)
22                 return NULL;
23
24         bounds = isl_mat_alloc(ctx, 1+bset->dim, 1+bset->dim);
25         if (!bounds)
26                 return NULL;
27
28         isl_int_set_si(bounds->row[0][0], 1);
29         isl_seq_clr(bounds->row[0]+1, bset->dim);
30         bounds->n_row = 1;
31
32         if (bset->n_ineq == 0)
33                 return bounds;
34
35         dirs = isl_mat_alloc(ctx, bset->dim, bset->dim);
36         if (!dirs) {
37                 isl_mat_free(ctx, bounds);
38                 return NULL;
39         }
40         isl_seq_cpy(dirs->row[0], bset->ineq[0]+1, dirs->n_col);
41         isl_seq_cpy(bounds->row[1], bset->ineq[0], bounds->n_col);
42         for (j = 1, n = 1; n < bset->dim && j < bset->n_ineq; ++j) {
43                 int pos;
44
45                 isl_seq_cpy(dirs->row[n], bset->ineq[j]+1, dirs->n_col);
46
47                 pos = isl_seq_first_non_zero(dirs->row[n], dirs->n_col);
48                 if (pos < 0)
49                         continue;
50                 for (i = 0; i < n; ++i) {
51                         int pos_i;
52                         pos_i = isl_seq_first_non_zero(dirs->row[i], dirs->n_col);
53                         if (pos_i < pos)
54                                 continue;
55                         if (pos_i > pos)
56                                 break;
57                         isl_seq_elim(dirs->row[n], dirs->row[i], pos,
58                                         dirs->n_col, NULL);
59                         pos = isl_seq_first_non_zero(dirs->row[n], dirs->n_col);
60                         if (pos < 0)
61                                 break;
62                 }
63                 if (pos < 0)
64                         continue;
65                 if (i < n) {
66                         int k;
67                         isl_int *t = dirs->row[n];
68                         for (k = n; k > i; --k)
69                                 dirs->row[k] = dirs->row[k-1];
70                         dirs->row[i] = t;
71                 }
72                 ++n;
73                 isl_seq_cpy(bounds->row[n], bset->ineq[j], bounds->n_col);
74         }
75         isl_mat_free(ctx, dirs);
76         bounds->n_row = 1+n;
77         return bounds;
78 }
79
80 /* Skew into positive orthant and project out lineality space */
81 static struct isl_basic_set *isl_basic_set_skew_to_positive_orthant(
82         struct isl_basic_set *bset, struct isl_mat **T)
83 {
84         struct isl_mat *U = NULL;
85         struct isl_mat *bounds = NULL;
86         int i, j;
87         unsigned old_dim, new_dim;
88         struct isl_ctx *ctx;
89
90         *T = NULL;
91         if (!bset)
92                 return NULL;
93
94         ctx = bset->ctx;
95         isl_assert(ctx, bset->nparam == 0, goto error);
96         isl_assert(ctx, bset->n_div == 0, goto error);
97         isl_assert(ctx, bset->n_eq == 0, goto error);
98         
99         /* Try to move (multiples of) unit rows up. */
100         for (i = 0, j = 0; i < bset->n_ineq; ++i) {
101                 int pos = isl_seq_first_non_zero(bset->ineq[i]+1,
102                                                             bset->dim);
103                 if (pos < 0)
104                         continue;
105                 if (isl_seq_first_non_zero(bset->ineq[i]+1+pos+1,
106                                                 bset->dim-pos-1) >= 0)
107                         continue;
108                 if (i != j)
109                         swap_inequality(bset, i, j);
110                 ++j;
111         }
112         bounds = independent_bounds(ctx, bset);
113         if (!bounds)
114                 goto error;
115         old_dim = bset->dim;
116         new_dim = bounds->n_row - 1;
117         bounds = isl_mat_left_hermite(ctx, bounds, 1, &U, NULL);
118         if (!bounds)
119                 goto error;
120         U = isl_mat_drop_cols(ctx, U, 1 + new_dim, old_dim - new_dim);
121         bset = isl_basic_set_preimage(ctx, bset, isl_mat_copy(ctx, U));
122         if (!bset)
123                 goto error;
124         *T = U;
125         isl_mat_free(ctx, bounds);
126         return bset;
127 error:
128         isl_mat_free(ctx, bounds);
129         isl_mat_free(ctx, U);
130         isl_basic_set_free(bset);
131         return NULL;
132 }
133
134 struct isl_vec *isl_pip_basic_set_sample(struct isl_basic_set *bset)
135 {
136         PipOptions      *options = NULL;
137         PipMatrix       *domain = NULL;
138         PipQuast        *sol = NULL;
139         struct isl_vec *vec = NULL;
140         unsigned        dim;
141         struct isl_mat *T;
142         struct isl_ctx *ctx;
143
144         if (!bset)
145                 goto error;
146         ctx = bset->ctx;
147         isl_assert(ctx, bset->nparam == 0, goto error);
148         isl_assert(ctx, bset->n_div == 0, goto error);
149         bset = isl_basic_set_skew_to_positive_orthant(bset, &T);
150         if (!bset)
151                 goto error;
152         dim = bset->dim;
153         domain = isl_basic_map_to_pip((struct isl_basic_map *)bset, 0, 0, 0);
154         if (!domain)
155                 goto error;
156
157         options = pip_options_init();
158         if (!options)
159                 goto error;
160         sol = pip_solve(domain, NULL, -1, options);
161         if (!sol)
162                 goto error;
163         if (!sol->list) {
164                 vec = isl_vec_alloc(ctx, 0);
165                 isl_mat_free(ctx, T);
166         } else {
167                 PipList *l;
168                 int i;
169                 vec = isl_vec_alloc(ctx, 1 + dim);
170                 if (!vec)
171                         goto error;
172                 isl_int_set_si(vec->block.data[0], 1);
173                 for (i = 0, l = sol->list; l && i < dim; ++i, l = l->next) {
174                         isl_seq_cpy_from_pip(&vec->block.data[1+i],
175                                         &l->vector->the_vector[0], 1);
176                         isl_assert(ctx, !entier_zero_p(l->vector->the_deno[0]),
177                                         goto error);
178                 }
179                 isl_assert(ctx, i == dim, goto error);
180                 vec = isl_mat_vec_product(ctx, T, vec);
181         }
182
183         pip_quast_free(sol);
184         pip_options_free(options);
185         pip_matrix_free(domain);
186
187         isl_basic_set_free(bset);
188         return vec;
189 error:
190         isl_vec_free(ctx, vec);
191         isl_basic_set_free(bset);
192         if (sol)
193                 pip_quast_free(sol);
194         if (domain)
195                 pip_matrix_free(domain);
196         return NULL;
197 }