polytope_scan.c

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