basis_reduction_templ.c

   1 #include <barvinok/basis_reduction.h>
   2
   3 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   4
   5 static void save_alpha(GBR_LP *lp, int first, int n, GBR_type *alpha)
   6 {
   7     int i;
   8
   9     for (i = 0; i < n; ++i)
  10         GBR_lp_get_alpha(lp, first+i, &alpha[i]);
  11 }
  12
  13 /* This function implements the algorithm described in
  14  * "An Implementation of the Generalized Basis Reduction Algorithm
  15  *  for Integer Programming" of Cook el al. to compute a reduced basis.
  16  * We use \epsilon = 1/4.
  17  */
  18 Matrix *Polyhedron_Reduced_Basis(Polyhedron *P)
  19 {
  20     int dim = P->Dimension;
  21     int i;
  22     Matrix *basis = Identity(dim);
  23     GBR_LP *lp;
  24     GBR_type F_old, alpha, F_new;
  25     int row;
  26     Value one, tmp;
  27     Vector *b_tmp;
  28     GBR_type *F;
  29     GBR_type *alpha_buffer[2];
  30     GBR_type *alpha_saved;
  31     GBR_type F_saved;
  32     int use_saved = 0;
  33     Value mu[2];
  34     GBR_type mu_F[2];
  35
  36     if (P->Dimension == 1)
  37         return basis;
  38
  39     value_init(one);
  40     value_init(tmp);
  41     value_set_si(one, 1);
  42     value_init(mu[0]);
  43     value_init(mu[1]);
  44
  45     b_tmp = Vector_Alloc(dim);
  46
  47     F = ALLOCN(GBR_type, dim);
  48     alpha_buffer[0] = ALLOCN(GBR_type, dim);
  49     alpha_buffer[1] = ALLOCN(GBR_type, dim);
  50     alpha_saved = alpha_buffer[0];
  51
  52     for (i = 0; i < dim; ++i) {
  53         GBR_init(F[i]);
  54         GBR_init(alpha_buffer[0][i]);
  55         GBR_init(alpha_buffer[1][i]);
  56     }
  57     GBR_init(alpha);
  58     GBR_init(F_old);
  59     GBR_init(F_new);
  60     GBR_init(F_saved);
  61     GBR_init(mu_F[0]);
  62     GBR_init(mu_F[1]);
  63
  64     lp = GBR_lp_init(P);
  65
  66     i = 0;
  67
  68     GBR_lp_set_obj(lp, basis->p[0], dim);
  69     if (GBR_lp_solve(lp))
  70         goto unbounded;
  71     GBR_lp_get_obj_val(lp, &F[0]);
  72
  73     do {
  74         if (use_saved) {
  75             row = GBR_lp_next_row(lp);
  76             GBR_set(F_new, F_saved);
  77             GBR_set(alpha, alpha_saved[i]);
  78         } else {
  79             row = GBR_lp_add_row(lp, basis->p[i], dim);
  80             GBR_lp_set_obj(lp, basis->p[i+1], dim);
  81             if (GBR_lp_solve(lp))
  82                 goto unbounded;
  83             GBR_lp_get_obj_val(lp, &F_new);
  84
  85             GBR_lp_get_alpha(lp, row, &alpha);
  86
  87             if (i > 0)
  88                 save_alpha(lp, row-i, i, alpha_saved);
  89
  90             GBR_lp_del_row(lp);
  91         }
  92         GBR_set(F[i+1], F_new);
  93
  94         GBR_floor(mu[0], alpha);
  95         GBR_ceil(mu[1], alpha);
  96
  97         if (value_eq(mu[0], mu[1]))
  98             value_assign(tmp, mu[0]);
  99         else {
 100             int j;
 101
 102             for (j = 0; j <= 1; ++j) {
 103                 value_assign(tmp, mu[j]);
 104                 Vector_Combine(basis->p[i+1], basis->p[i], b_tmp->p, one, tmp, dim);
 105                 GBR_lp_set_obj(lp, b_tmp->p, dim);
 106                 if (GBR_lp_solve(lp))
 107                     goto unbounded;
 108                 GBR_lp_get_obj_val(lp, &mu_F[j]);
 109                 if (i > 0)
 110                     save_alpha(lp, row-i, i, alpha_buffer[j]);
 111             }
 112
 113             if (GBR_lt(mu_F[0], mu_F[1]))
 114                 j = 0;
 115             else
 116                 j = 1;
 117
 118             value_assign(tmp, mu[j]);
 119             GBR_set(F_new, mu_F[j]);
 120             alpha_saved = alpha_buffer[j];
 121         }
 122         Vector_Combine(basis->p[i+1], basis->p[i], basis->p[i+1], one, tmp, dim);
 123
 124         GBR_set(F_old, F[i]);
 125
 126         use_saved = 0;
 127         /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
 128         GBR_set_ui(mu_F[0], 4);
 129         GBR_mul(mu_F[0], mu_F[0], F_new);
 130         GBR_set_ui(mu_F[1], 3);
 131         GBR_mul(mu_F[1], mu_F[1], F_old);
 132         if (GBR_lt(mu_F[0], mu_F[1])) {
 133             Vector_Exchange(basis->p[i], basis->p[i+1], dim);
 134             if (i > 0) {
 135                 use_saved = 1;
 136                 GBR_set(F_saved, F_new);
 137                 GBR_lp_del_row(lp);
 138                 --i;
 139             } else
 140                 GBR_set(F[0], F_new);
 141         } else {
 142             GBR_lp_add_row(lp, basis->p[i], dim);
 143             ++i;
 144         }
 145     } while (i < dim-1);
 146
 147     if (0) {
 148 unbounded:
 149         Matrix_Free(basis);
 150         basis = NULL;
 151     }
 152     Vector_Free(b_tmp);
 153
 154     value_clear(one);
 155     value_clear(tmp);
 156     value_clear(mu[0]);
 157     value_clear(mu[1]);
 158     for (i = 0; i < dim; ++i) {
 159         GBR_clear(F[i]);
 160         GBR_clear(alpha_buffer[0][i]);
 161         GBR_clear(alpha_buffer[1][i]);
 162     }
 163     free(F);
 164     free(alpha_buffer[0]);
 165     free(alpha_buffer[1]);
 166     GBR_clear(alpha);
 167     GBR_clear(F_old);
 168     GBR_clear(F_new);
 169     GBR_clear(F_saved);
 170     GBR_clear(mu_F[0]);
 171     GBR_clear(mu_F[1]);
 172
 173     GBR_lp_delete(lp);
 174
 175     return basis;
 176 }