polytope_scan.c

   1 #include <assert.h>
   2 #include <unistd.h>
   3 #include <stdlib.h>
   4 #include <strings.h>
   5 #include <barvinok/util.h>
   6 #include <barvinok/options.h>
   7 #include <barvinok/basis_reduction.h>
   8
   9 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  10
  11 struct options {
  12     struct barvinok_options *barvinok;
  13     int direct;
  14 };
  15
  16 struct isl_arg options_arg[] = {
  17 ISL_ARG_CHILD(struct options, barvinok, NULL, barvinok_options_arg, NULL)
  18 ISL_ARG_BOOL(struct options, direct, 'd', "direct", 0,
  19         "don't apply generalized basis reduction first")
  20 ISL_ARG_END
  21 };
  22
  23 ISL_ARG_DEF(options, struct options, options_arg)
  24
  25 static void scan_poly(Polyhedron *S, int pos, Value *z, Matrix *T)
  26 {
  27     if (!S) {
  28         int k;
  29         Vector *v;
  30
  31         v = Vector_Alloc(T->NbRows);
  32         Matrix_Vector_Product(T, z+1, v->p);
  33         value_print(stdout, VALUE_FMT, v->p[0]);
  34         for (k=1; k < pos; ++k) {
  35             printf(", ");
  36             value_print(stdout,VALUE_FMT, v->p[k]);
  37         }
  38         Vector_Free(v);
  39         printf("\n");
  40     } else {
  41         int ok;
  42         Value LB, UB, tmp;
  43         value_init(LB);
  44         value_init(UB);
  45         value_init(tmp);
  46         ok = !(lower_upper_bounds(1+pos, S, z, &LB, &UB));
  47         assert(ok);
  48         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
  49             value_assign(z[pos+1], tmp);
  50             scan_poly(S->next, pos+1, z, T);
  51         }
  52         value_set_si(z[pos+1], 0);
  53         value_clear(LB);
  54         value_clear(UB);
  55         value_clear(tmp);
  56     }
  57 }
  58
  59 int main(int argc, char **argv)
  60 {
  61     Polyhedron *A, *P, *U, *S;
  62     Value *p;
  63     int i, j, ok;
  64     Matrix *basis, *T, *inv;
  65     int c, ind = 0;
  66     struct options *options = options_new_with_defaults();
  67
  68     argc = options_parse(options, argc, argv, ISL_ARG_ALL);
  69
  70     A = Polyhedron_Read(options->barvinok->MaxRays);
  71
  72     if (options->direct) {
  73         inv = Identity(A->Dimension+1);
  74         P = A;
  75     } else {
  76         basis = Polyhedron_Reduced_Basis(A, options->barvinok);
  77
  78         T = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
  79         inv = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
  80         for (i = 0; i < A->Dimension; ++i)
  81             for (j = 0; j < A->Dimension; ++j)
  82                 value_assign(T->p[i][j], basis->p[i][j]);
  83         value_set_si(T->p[A->Dimension][A->Dimension], 1);
  84         Matrix_Free(basis);
  85
  86         ok = Matrix_Inverse(T, inv);
  87         assert(ok);
  88         Matrix_Free(T);
  89
  90         P = Polyhedron_Preimage(A, inv, options->barvinok->MaxRays);
  91         Polyhedron_Free(A);
  92     }
  93
  94     U = Universe_Polyhedron(0);
  95     S = Polyhedron_Scan(P, U, options->barvinok->MaxRays);
  96
  97     p = ALLOCN(Value, P->Dimension+2);
  98     for(i=0;i<=P->Dimension;i++) {
  99         value_init(p[i]);
 100         value_set_si(p[i],0);
 101     }
 102     value_init(p[i]);
 103     value_set_si(p[i],1);
 104
 105     scan_poly(S, 0, p, inv);
 106
 107     Matrix_Free(inv);
 108     for(i=0;i<=(P->Dimension+1);i++)
 109         value_clear(p[i]);
 110     free(p);
 111     Domain_Free(S);
 112     Polyhedron_Free(P);
 113     Polyhedron_Free(U);
 114     options_free(options);
 115     return 0;
 116 }