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 "config.h"
   8
   9 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  10
  11 #ifndef HAVE_GETOPT_H
  12 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
  13 #else
  14 #include <getopt.h>
  15 struct option options[] = {
  16     { "direct",    no_argument,  0,  'd' },
  17     { "version",   no_argument,  0,  'V' },
  18     { 0, 0, 0, 0 }
  19 };
  20 #endif
  21
  22 static void scan_poly(Polyhedron *S, int pos, Value *z, Matrix *T)
  23 {
  24     if (!S) {
  25         int k;
  26         Vector *v;
  27
  28         v = Vector_Alloc(T->NbRows);
  29         Matrix_Vector_Product(T, z+1, v->p);
  30         value_print(stdout, VALUE_FMT, v->p[0]);
  31         for (k=1; k < pos; ++k) {
  32             printf(", ");
  33             value_print(stdout,VALUE_FMT, v->p[k]);
  34         }
  35         Vector_Free(v);
  36         printf("\n");
  37     } else {
  38         int ok;
  39         Value LB, UB, tmp;
  40         value_init(LB);
  41         value_init(UB);
  42         value_init(tmp);
  43         ok = !(lower_upper_bounds(1+pos, S, z, &LB, &UB));
  44         assert(ok);
  45         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
  46             value_assign(z[pos+1], tmp);
  47             scan_poly(S->next, pos+1, z, T);
  48         }
  49         value_set_si(z[pos+1], 0);
  50         value_clear(LB);
  51         value_clear(UB);
  52         value_clear(tmp);
  53     }
  54 }
  55
  56 int main(int argc, char **argv)
  57 {
  58     Polyhedron *A, *P, *U, *S;
  59     Value *p;
  60     int i, j, ok;
  61     Matrix *basis, *T, *inv;
  62     int c, ind = 0;
  63     int direct = 0;
  64     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
  65
  66     while ((c = getopt_long(argc, argv, "dV", options, &ind)) != -1) {
  67         switch (c) {
  68         case 'd':
  69             direct = 1;
  70             break;
  71         case 'V':
  72             printf(barvinok_version());
  73             exit(0);
  74             break;
  75         }
  76     }
  77
  78     A = Polyhedron_Read(bv_options->MaxRays);
  79
  80     if (direct) {
  81         inv = Identity(A->Dimension+1);
  82         P = A;
  83     } else {
  84         basis = Polyhedron_Reduced_Basis(A, bv_options);
  85
  86         T = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
  87         inv = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
  88         for (i = 0; i < A->Dimension; ++i)
  89             for (j = 0; j < A->Dimension; ++j)
  90                 value_assign(T->p[i][j], basis->p[i][j]);
  91         value_set_si(T->p[A->Dimension][A->Dimension], 1);
  92         Matrix_Free(basis);
  93
  94         ok = Matrix_Inverse(T, inv);
  95         assert(ok);
  96         Matrix_Free(T);
  97
  98         P = Polyhedron_Preimage(A, inv, bv_options->MaxRays);
  99         Polyhedron_Free(A);
 100     }
 101
 102     U = Universe_Polyhedron(0);
 103     S = Polyhedron_Scan(P, U, bv_options->MaxRays);
 104
 105     p = ALLOCN(Value, P->Dimension+2);
 106     for(i=0;i<=P->Dimension;i++) {
 107         value_init(p[i]);
 108         value_set_si(p[i],0);
 109     }
 110     value_init(p[i]);
 111     value_set_si(p[i],1);
 112
 113     scan_poly(S, 0, p, inv);
 114
 115     Matrix_Free(inv);
 116     for(i=0;i<=(P->Dimension+1);i++)
 117         value_clear(p[i]);
 118     free(p);
 119     Domain_Free(S);
 120     Polyhedron_Free(P);
 121     Polyhedron_Free(U);
 122     barvinok_options_free(bv_options);
 123     return 0;
 124 }