ilp.c

   1 #include <barvinok/options.h>
   2 #include <barvinok/sample.h>
   3 #include "ilp.h"
   4
   5 static int matrix_add(Matrix *M, int n, Value *v)
   6 {
   7     if (n >= M->NbRows)
   8         Matrix_Extend(M, 3*(M->NbRows+10)/2);
   9     Vector_Copy(v, M->p[n], M->NbColumns);
  10     return n+1;
  11 }
  12
  13 /* Computes and returns an integer point of P attaining the minimum
  14  * value in [min, max] with respect to the linear objective function obj.
  15  * If there is no such point, then a NULL pointer is returned.
  16  * If misses is not a NULL pointer, then any other points found along
  17  * the way are stored in this matrix, updating *n_misses.
  18  */
  19 Vector *Polyhedron_Integer_Minimum(Polyhedron *P, Value *obj,
  20                                     Value min, Value max,
  21                                     Matrix *misses, int *n_misses,
  22                                     struct barvinok_options *options)
  23 {
  24     int divide = 0;
  25     Vector *opt = NULL;
  26     Value l, u, tmp;
  27     value_init(l);
  28     value_init(u);
  29     value_init(tmp);
  30     value_assign(l, min);
  31     value_assign(u, max);
  32
  33     while (value_le(l, u)) {
  34         int i;
  35         Matrix *M = Matrix_Alloc(P->NbConstraints + 2, 2 + P->Dimension);
  36         Polyhedron *Q;
  37         Vector *sample;
  38
  39         if (divide) {
  40             value_subtract(tmp, u, l);
  41             mpz_fdiv_q_ui(tmp, tmp, 2);
  42             value_addto(tmp, tmp, l);
  43         } else {
  44             value_assign(tmp, u);
  45         }
  46
  47         for (i = 0; i < P->NbConstraints; ++i)
  48             Vector_Copy(P->Constraint[i], M->p[i+2], P->Dimension+2);
  49         value_set_si(M->p[0][0], 1);
  50         Vector_Copy(obj, M->p[0]+1, P->Dimension);
  51         value_oppose(M->p[0][1+P->Dimension], l);
  52         value_set_si(M->p[1][0], 1);
  53         Vector_Oppose(obj, M->p[1]+1, P->Dimension);
  54         value_assign(M->p[1][1+P->Dimension], tmp);
  55
  56         Q = Constraints2Polyhedron(M, options->MaxRays);
  57         Matrix_Free(M);
  58
  59         sample = Polyhedron_Sample(Q, options);
  60         Polyhedron_Free(Q);
  61
  62         if (sample) {
  63             Inner_Product(obj, sample->p, P->Dimension, &u);
  64             value_decrement(u, u);
  65             divide = 1;
  66             if (opt) {
  67                 if (misses)
  68                     *n_misses = matrix_add(misses, *n_misses, opt->p);
  69                 Vector_Free(opt);
  70             }
  71             opt = sample;
  72         } else {
  73             if (!divide)
  74                 break;
  75             value_increment(l, tmp);
  76             divide = 0;
  77         }
  78     }
  79
  80     value_clear(l);
  81     value_clear(u);
  82     value_clear(tmp);
  83
  84     return opt;
  85 }