polytope_scan.c

   1 #include <assert.h>
   2 #include "isl_basis_reduction.h"
   3 #include "isl_equalities.h"
   4 #include "isl_tab.h"
   5 #include "isl_vec.h"
   6
   7 /* The input of this program is the same as that of the "polytope_scan"
   8  * program from the barvinok distribution.
   9  *
  10  * Constraints of set is PolyLib format.
  11  *
  12  * The input set is assumed to be bounded.
  13  */
  14
  15 static struct isl_mat *isl_basic_set_samples(struct isl_basic_set *bset);
  16
  17 static struct isl_mat *samples_eq(struct isl_basic_set *bset)
  18 {
  19         struct isl_mat *T;
  20         struct isl_mat *samples;
  21
  22         bset = isl_basic_set_remove_equalities(bset, &T, NULL);
  23         samples = isl_basic_set_samples(bset);
  24         return isl_mat_product(samples, isl_mat_transpose(T));
  25 }
  26
  27 /* Add the current sample value of the tableau "tab" to the list
  28  * in "samples".
  29  */
  30 static struct isl_mat *add_solution(struct isl_mat *samples,
  31         struct isl_tab *tab)
  32 {
  33         struct isl_vec *sample;
  34
  35         if (!samples || !tab)
  36                 goto error;
  37         samples = isl_mat_extend(samples, samples->n_row + 1, samples->n_col);
  38         if (!samples)
  39                 return NULL;
  40         sample = isl_tab_get_sample_value(tab);
  41         if (!sample)
  42                 goto error;
  43         isl_seq_cpy(samples->row[samples->n_row - 1], sample->el, sample->size);
  44         isl_vec_free(sample);
  45         return samples;
  46 error:
  47         isl_mat_free(samples);
  48         return NULL;
  49 }
  50
  51 /* Look for and return all integer points in "bset", which is assumed
  52  * to be unbounded.
  53  *
  54  * We first compute a reduced basis for the set and then scan
  55  * the set in the directions of this basis.
  56  * We basically perform a depth first search, where in each level i
  57  * we compute the range in the i-th basis vector direction, given
  58  * fixed values in the directions of the previous basis vector.
  59  * We then add an equality to the tableau fixing the value in the
  60  * direction of the current basis vector to each value in the range
  61  * in turn and then continue to the next level.
  62  *
  63  * The search is implemented iteratively.  "level" identifies the current
  64  * basis vector.  "init" is true if we want the first value at the current
  65  * level and false if we want the next value.
  66  * Solutions are added in the leaves of the search tree, i.e., after
  67  * we have fixed a value in each direction of the basis.
  68  */
  69 static struct isl_mat *isl_basic_set_samples(struct isl_basic_set *bset)
  70 {
  71         unsigned dim;
  72         struct isl_mat *B = NULL;
  73         struct isl_tab *tab = NULL;
  74         struct isl_vec *min;
  75         struct isl_vec *max;
  76         struct isl_mat *samples;
  77         struct isl_tab_undo **snap;
  78         int level;
  79         int init;
  80         enum isl_lp_result res;
  81
  82         if (bset->n_eq)
  83                 return samples_eq(bset);
  84
  85         dim = isl_basic_set_total_dim(bset);
  86
  87         min = isl_vec_alloc(bset->ctx, dim);
  88         max = isl_vec_alloc(bset->ctx, dim);
  89         samples = isl_mat_alloc(bset->ctx, 0, 1 + dim);
  90         snap = isl_alloc_array(bset->ctx, struct isl_tab_undo *, dim);
  91
  92         if (!min || !max || !samples || !snap)
  93                 goto error;
  94
  95         if (1)
  96                 B = isl_basic_set_reduced_basis(bset);
  97         else
  98                 B = isl_mat_identity(bset->ctx, dim);
  99         B = isl_mat_lin_to_aff(B);
 100         if (!B)
 101                 goto error;
 102
 103         tab = isl_tab_from_basic_set(bset);
 104
 105         level = 0;
 106         init = 1;
 107
 108         while (level >= 0) {
 109                 int empty = 0;
 110                 if (init) {
 111                         res = isl_tab_min(tab, B->row[1 + level],
 112                                     bset->ctx->one, &min->el[level], NULL, 0);
 113                         if (res == isl_lp_empty)
 114                                 empty = 1;
 115                         if (res == isl_lp_error || res == isl_lp_unbounded)
 116                                 goto error;
 117                         isl_seq_neg(B->row[1 + level] + 1,
 118                                     B->row[1 + level] + 1, dim);
 119                         res = isl_tab_min(tab, B->row[1 + level],
 120                                     bset->ctx->one, &max->el[level], NULL, 0);
 121                         isl_seq_neg(B->row[1 + level] + 1,
 122                                     B->row[1 + level] + 1, dim);
 123                         isl_int_neg(max->el[level], max->el[level]);
 124                         if (res == isl_lp_empty)
 125                                 empty = 1;
 126                         if (res == isl_lp_error || res == isl_lp_unbounded)
 127                                 goto error;
 128                         snap[level] = isl_tab_snap(tab);
 129                 } else
 130                         isl_int_add_ui(min->el[level], min->el[level], 1);
 131
 132                 if (empty || isl_int_gt(min->el[level], max->el[level])) {
 133                         level--;
 134                         init = 0;
 135                         if (level >= 0)
 136                                 isl_tab_rollback(tab, snap[level]);
 137                         continue;
 138                 }
 139                 isl_int_neg(B->row[1 + level][0], min->el[level]);
 140                 tab = isl_tab_add_valid_eq(tab, B->row[1 + level]);
 141                 isl_int_set_si(B->row[1 + level][0], 0);
 142                 if (level < dim - 1) {
 143                         ++level;
 144                         init = 1;
 145                         continue;
 146                 }
 147                 samples = add_solution(samples, tab);
 148                 init = 0;
 149                 isl_tab_rollback(tab, snap[level]);
 150         }
 151
 152         isl_tab_free(tab);
 153         free(snap);
 154         isl_vec_free(min);
 155         isl_vec_free(max);
 156         isl_basic_set_free(bset);
 157         isl_mat_free(B);
 158         return samples;
 159 error:
 160         isl_tab_free(tab);
 161         free(snap);
 162         isl_mat_free(samples);
 163         isl_vec_free(min);
 164         isl_vec_free(max);
 165         isl_basic_set_free(bset);
 166         isl_mat_free(B);
 167         return NULL;
 168 }
 169
 170 int main(int argc, char **argv)
 171 {
 172         struct isl_ctx *ctx = isl_ctx_alloc();
 173         struct isl_basic_set *bset;
 174         struct isl_mat *samples;
 175
 176         bset = isl_basic_set_read_from_file(ctx, stdin, 0, ISL_FORMAT_POLYLIB);
 177         samples = isl_basic_set_samples(bset);
 178         isl_mat_dump(samples, stdout, 0);
 179         isl_mat_free(samples);
 180         isl_ctx_free(ctx);
 181
 182         return 0;
 183 }