polytope_scan.c

   1 #include <unistd.h>
   2 #include <stdlib.h>
   3 #include <strings.h>
   4 #include <polylib/polylibgmp.h>
   5 #include <barvinok/util.h>
   6 #include "basis_reduction.h"
   7 #include "config.h"
   8
   9 #ifdef HAVE_GROWING_CHERNIKOVA
  10 #define MAXRAYS    (POL_NO_DUAL | POL_INTEGER)
  11 #else
  12 #define MAXRAYS  600
  13 #endif
  14
  15 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  16
  17 #ifndef HAVE_GETOPT_H
  18 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
  19 #else
  20 #include <getopt.h>
  21 struct option options[] = {
  22     { "direct",    no_argument,  0,  'd' },
  23     { "version",   no_argument,  0,  'V' },
  24     { 0, 0, 0, 0 }
  25 };
  26 #endif
  27
  28 static Polyhedron *Polyhedron_Read()
  29 {
  30     int vertices = 0;
  31     unsigned NbRows, NbColumns;
  32     Matrix *M;
  33     Polyhedron *P;
  34     char s[128];
  35
  36     while (fgets(s, sizeof(s), stdin)) {
  37         if (*s == '#')
  38             continue;
  39         if (strncasecmp(s, "vertices", sizeof("vertices")-1) == 0)
  40             vertices = 1;
  41         if (sscanf(s, "%u %u", &NbRows, &NbColumns) == 2)
  42             break;
  43     }
  44     if (feof(stdin))
  45         return NULL;
  46     M = Matrix_Alloc(NbRows,NbColumns);
  47     Matrix_Read_Input(M);
  48     if (vertices)
  49         P = Rays2Polyhedron(M, MAXRAYS);
  50     else
  51         P = Constraints2Polyhedron(M, MAXRAYS);
  52     Matrix_Free(M);
  53     return P;
  54 }
  55
  56 static void scan_poly(Polyhedron *S, int pos, Value *z, Matrix *T)
  57 {
  58     if (!S) {
  59         int k;
  60         Vector *v;
  61
  62         v = Vector_Alloc(T->NbRows);
  63         Matrix_Vector_Product(T, z+1, v->p);
  64         value_print(stdout, VALUE_FMT, v->p[0]);
  65         for (k=1; k < pos; ++k) {
  66             printf(", ");
  67             value_print(stdout,VALUE_FMT, v->p[k]);
  68         }
  69         Vector_Free(v);
  70         printf("\n");
  71     } else {
  72         int ok;
  73         Value LB, UB, tmp;
  74         value_init(LB);
  75         value_init(UB);
  76         value_init(tmp);
  77         ok = !(lower_upper_bounds(1+pos, S, z, &LB, &UB));
  78         assert(ok);
  79         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
  80             value_assign(z[pos+1], tmp);
  81             scan_poly(S->next, pos+1, z, T);
  82         }
  83         value_set_si(z[pos+1], 0);
  84         value_clear(LB);
  85         value_clear(UB);
  86         value_clear(tmp);
  87     }
  88 }
  89
  90 int main(int argc, char **argv)
  91 {
  92     Polyhedron *A, *P, *U, *S;
  93     Value *p;
  94     int i, j, ok;
  95     Matrix *basis, *T, *inv;
  96     int c, ind = 0;
  97     int direct = 0;
  98
  99     while ((c = getopt_long(argc, argv, "dV", options, &ind)) != -1) {
 100         switch (c) {
 101         case 'd':
 102             direct = 1;
 103             break;
 104         case 'V':
 105             printf(barvinok_version());
 106             exit(0);
 107             break;
 108         }
 109     }
 110
 111     A = Polyhedron_Read();
 112
 113     if (direct) {
 114         inv = Identity(A->Dimension+1);
 115         P = A;
 116     } else {
 117         basis = reduced_basis(A);
 118
 119         T = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
 120         inv = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
 121         for (i = 0; i < A->Dimension; ++i)
 122             for (j = 0; j < A->Dimension; ++j)
 123                 value_assign(T->p[i][j], basis->p[i][j]);
 124         value_set_si(T->p[A->Dimension][A->Dimension], 1);
 125         Matrix_Free(basis);
 126
 127         ok = Matrix_Inverse(T, inv);
 128         assert(ok);
 129         Matrix_Free(T);
 130
 131         P = Polyhedron_Preimage(A, inv, MAXRAYS);
 132         Polyhedron_Free(A);
 133     }
 134
 135     U = Universe_Polyhedron(0);
 136     S = Polyhedron_Scan(P, U, MAXRAYS);
 137
 138     p = ALLOCN(Value, P->Dimension+2);
 139     for(i=0;i<=P->Dimension;i++) {
 140         value_init(p[i]);
 141         value_set_si(p[i],0);
 142     }
 143     value_init(p[i]);
 144     value_set_si(p[i],1);
 145
 146     scan_poly(S, 0, p, inv);
 147
 148     Matrix_Free(inv);
 149     for(i=0;i<=(P->Dimension+1);i++)
 150         value_clear(p[i]);
 151     free(p);
 152     Domain_Free(S);
 153     Polyhedron_Free(P);
 154     Polyhedron_Free(U);
 155     return 0;
 156 }