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
  98         if (1)
  99                 B = isl_tab_reduced_basis(tab);
 100         else
 101                 B = isl_mat_identity(bset->ctx, dim);
 102         B = isl_mat_lin_to_aff(B);
 103         if (!B)
 104                 goto error;
 105
 106         level = 0;
 107         init = 1;
 108
 109         while (level >= 0) {
 110                 int empty = 0;
 111                 if (init) {
 112                         res = isl_tab_min(tab, B->row[1 + level],
 113                                     bset->ctx->one, &min->el[level], NULL, 0);
 114                         if (res == isl_lp_empty)
 115                                 empty = 1;
 116                         if (res == isl_lp_error || res == isl_lp_unbounded)
 117                                 goto error;
 118                         isl_seq_neg(B->row[1 + level] + 1,
 119                                     B->row[1 + level] + 1, dim);
 120                         res = isl_tab_min(tab, B->row[1 + level],
 121                                     bset->ctx->one, &max->el[level], NULL, 0);
 122                         isl_seq_neg(B->row[1 + level] + 1,
 123                                     B->row[1 + level] + 1, dim);
 124                         isl_int_neg(max->el[level], max->el[level]);
 125                         if (res == isl_lp_empty)
 126                                 empty = 1;
 127                         if (res == isl_lp_error || res == isl_lp_unbounded)
 128                                 goto error;
 129                         snap[level] = isl_tab_snap(tab);
 130                 } else
 131                         isl_int_add_ui(min->el[level], min->el[level], 1);
 132
 133                 if (empty || isl_int_gt(min->el[level], max->el[level])) {
 134                         level--;
 135                         init = 0;
 136                         if (level >= 0)
 137                                 isl_tab_rollback(tab, snap[level]);
 138                         continue;
 139                 }
 140                 isl_int_neg(B->row[1 + level][0], min->el[level]);
 141                 tab = isl_tab_add_valid_eq(tab, B->row[1 + level]);
 142                 isl_int_set_si(B->row[1 + level][0], 0);
 143                 if (level < dim - 1) {
 144                         ++level;
 145                         init = 1;
 146                         continue;
 147                 }
 148                 samples = add_solution(samples, tab);
 149                 init = 0;
 150                 isl_tab_rollback(tab, snap[level]);
 151         }
 152
 153         isl_tab_free(tab);
 154         free(snap);
 155         isl_vec_free(min);
 156         isl_vec_free(max);
 157         isl_basic_set_free(bset);
 158         isl_mat_free(B);
 159         return samples;
 160 error:
 161         isl_tab_free(tab);
 162         free(snap);
 163         isl_mat_free(samples);
 164         isl_vec_free(min);
 165         isl_vec_free(max);
 166         isl_basic_set_free(bset);
 167         isl_mat_free(B);
 168         return NULL;
 169 }
 170
 171 int main(int argc, char **argv)
 172 {
 173         struct isl_ctx *ctx = isl_ctx_alloc();
 174         struct isl_basic_set *bset;
 175         struct isl_mat *samples;
 176
 177         bset = isl_basic_set_read_from_file(ctx, stdin, 0, ISL_FORMAT_POLYLIB);
 178         samples = isl_basic_set_samples(bset);
 179         isl_mat_dump(samples, stdout, 0);
 180         isl_mat_free(samples);
 181         isl_ctx_free(ctx);
 182
 183         return 0;
 184 }