polytope_scan.c

   1 #include <unistd.h>
   2 #include <stdlib.h>
   3 #include <strings.h>
   4 #include <barvinok/util.h>
   5 #include <barvinok/options.h>
   6 #include <barvinok/basis_reduction.h>
   7 #include "argp.h"
   8
   9 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  10
  11 struct argp_option argp_options[] = {
  12 #if defined(HAVE_LIBCDDGMP) || defined(HAVE_LIBGLPK)
  13     { "direct",         'd',  0,        0,
  14         "don't apply generalized basis reduction first" },
  15 #endif
  16     { 0 }
  17 };
  18
  19 struct options {
  20     struct barvinok_options *barvinok;
  21     int direct;
  22 };
  23
  24 error_t parse_opt(int key, char *arg, struct argp_state *state)
  25 {
  26     struct options *options = state->input;
  27
  28     switch (key) {
  29     case ARGP_KEY_INIT:
  30         state->child_inputs[0] = options->barvinok;
  31 #if defined(HAVE_LIBCDDGMP) || defined(HAVE_LIBGLPK)
  32         options->direct = 0;
  33 #else
  34         options->direct = 1;
  35 #endif
  36         break;
  37     case 'd':
  38         options->direct = 1;
  39         break;
  40     case 'V':
  41         printf(barvinok_version());
  42         exit(0);
  43     default:
  44         return ARGP_ERR_UNKNOWN;
  45     }
  46     return 0;
  47 }
  48
  49 static void scan_poly(Polyhedron *S, int pos, Value *z, Matrix *T)
  50 {
  51     if (!S) {
  52         int k;
  53         Vector *v;
  54
  55         v = Vector_Alloc(T->NbRows);
  56         Matrix_Vector_Product(T, z+1, v->p);
  57         value_print(stdout, VALUE_FMT, v->p[0]);
  58         for (k=1; k < pos; ++k) {
  59             printf(", ");
  60             value_print(stdout,VALUE_FMT, v->p[k]);
  61         }
  62         Vector_Free(v);
  63         printf("\n");
  64     } else {
  65         int ok;
  66         Value LB, UB, tmp;
  67         value_init(LB);
  68         value_init(UB);
  69         value_init(tmp);
  70         ok = !(lower_upper_bounds(1+pos, S, z, &LB, &UB));
  71         assert(ok);
  72         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
  73             value_assign(z[pos+1], tmp);
  74             scan_poly(S->next, pos+1, z, T);
  75         }
  76         value_set_si(z[pos+1], 0);
  77         value_clear(LB);
  78         value_clear(UB);
  79         value_clear(tmp);
  80     }
  81 }
  82
  83 int main(int argc, char **argv)
  84 {
  85     Polyhedron *A, *P, *U, *S;
  86     Value *p;
  87     int i, j, ok;
  88     Matrix *basis, *T, *inv;
  89     int c, ind = 0;
  90     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
  91     static struct argp_child argp_children[] = {
  92         { &barvinok_argp,    0, 0,  0 },
  93         { 0 }
  94     };
  95     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
  96     struct options options;
  97
  98     options.barvinok = bv_options;
  99
 100     set_program_name(argv[0]);
 101     argp_parse(&argp, argc, argv, 0, 0, &options);
 102
 103     A = Polyhedron_Read(bv_options->MaxRays);
 104
 105     if (options.direct) {
 106         inv = Identity(A->Dimension+1);
 107         P = A;
 108     } else {
 109         basis = Polyhedron_Reduced_Basis(A, bv_options);
 110
 111         T = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
 112         inv = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
 113         for (i = 0; i < A->Dimension; ++i)
 114             for (j = 0; j < A->Dimension; ++j)
 115                 value_assign(T->p[i][j], basis->p[i][j]);
 116         value_set_si(T->p[A->Dimension][A->Dimension], 1);
 117         Matrix_Free(basis);
 118
 119         ok = Matrix_Inverse(T, inv);
 120         assert(ok);
 121         Matrix_Free(T);
 122
 123         P = Polyhedron_Preimage(A, inv, bv_options->MaxRays);
 124         Polyhedron_Free(A);
 125     }
 126
 127     U = Universe_Polyhedron(0);
 128     S = Polyhedron_Scan(P, U, bv_options->MaxRays);
 129
 130     p = ALLOCN(Value, P->Dimension+2);
 131     for(i=0;i<=P->Dimension;i++) {
 132         value_init(p[i]);
 133         value_set_si(p[i],0);
 134     }
 135     value_init(p[i]);
 136     value_set_si(p[i],1);
 137
 138     scan_poly(S, 0, p, inv);
 139
 140     Matrix_Free(inv);
 141     for(i=0;i<=(P->Dimension+1);i++)
 142         value_clear(p[i]);
 143     free(p);
 144     Domain_Free(S);
 145     Polyhedron_Free(P);
 146     Polyhedron_Free(U);
 147     barvinok_options_free(bv_options);
 148     return 0;
 149 }