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     GBR_lp_solve(lp);
  70     GBR_lp_get_obj_val(lp, &F[0]);
  71
  72     do {
  73         if (use_saved) {
  74             row = GBR_lp_next_row(lp);
  75             GBR_set(F_new, F_saved);
  76             GBR_set(alpha, alpha_saved[i]);
  77         } else {
  78             row = GBR_lp_add_row(lp, basis->p[i], dim);
  79             GBR_lp_set_obj(lp, basis->p[i+1], dim);
  80             GBR_lp_solve(lp);
  81             GBR_lp_get_obj_val(lp, &F_new);
  82
  83             GBR_lp_get_alpha(lp, row, &alpha);
  84
  85             if (i > 0)
  86                 save_alpha(lp, row-i, i, alpha_saved);
  87
  88             GBR_lp_del_row(lp);
  89         }
  90         GBR_set(F[i+1], F_new);
  91
  92         GBR_floor(mu[0], alpha);
  93         GBR_ceil(mu[1], alpha);
  94
  95         if (value_eq(mu[0], mu[1]))
  96             value_assign(tmp, mu[0]);
  97         else {
  98             int j;
  99
 100             for (j = 0; j <= 1; ++j) {
 101                 value_assign(tmp, mu[j]);
 102                 Vector_Combine(basis->p[i+1], basis->p[i], b_tmp->p, one, tmp, dim);
 103                 GBR_lp_set_obj(lp, b_tmp->p, dim);
 104                 GBR_lp_solve(lp);
 105                 GBR_lp_get_obj_val(lp, &mu_F[j]);
 106                 if (i > 0)
 107                     save_alpha(lp, row-i, i, alpha_buffer[j]);
 108             }
 109
 110             if (GBR_lt(mu_F[0], mu_F[1]))
 111                 j = 0;
 112             else
 113                 j = 1;
 114
 115             value_assign(tmp, mu[j]);
 116             GBR_set(F_new, mu_F[j]);
 117             alpha_saved = alpha_buffer[j];
 118         }
 119         Vector_Combine(basis->p[i+1], basis->p[i], basis->p[i+1], one, tmp, dim);
 120
 121         GBR_set(F_old, F[i]);
 122
 123         use_saved = 0;
 124         /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
 125         GBR_set_ui(mu_F[0], 4);
 126         GBR_mul(mu_F[0], mu_F[0], F_new);
 127         GBR_set_ui(mu_F[1], 3);
 128         GBR_mul(mu_F[1], mu_F[1], F_old);
 129         if (GBR_lt(mu_F[0], mu_F[1])) {
 130             Vector_Exchange(basis->p[i], basis->p[i+1], dim);
 131             if (i > 0) {
 132                 use_saved = 1;
 133                 GBR_set(F_saved, F_new);
 134                 GBR_lp_del_row(lp);
 135                 --i;
 136             } else
 137                 GBR_set(F[0], F_new);
 138         } else {
 139             GBR_lp_add_row(lp, basis->p[i], dim);
 140             ++i;
 141         }
 142     } while (i < dim-1);
 143
 144     Vector_Free(b_tmp);
 145
 146     value_clear(one);
 147     value_clear(tmp);
 148     value_clear(mu[0]);
 149     value_clear(mu[1]);
 150     for (i = 0; i < dim; ++i) {
 151         GBR_clear(F[i]);
 152         GBR_clear(alpha_buffer[0][i]);
 153         GBR_clear(alpha_buffer[1][i]);
 154     }
 155     free(F);
 156     free(alpha_buffer[0]);
 157     free(alpha_buffer[1]);
 158     GBR_clear(alpha);
 159     GBR_clear(F_old);
 160     GBR_clear(F_new);
 161     GBR_clear(F_saved);
 162     GBR_clear(mu_F[0]);
 163     GBR_clear(mu_F[1]);
 164
 165     GBR_lp_delete(lp);
 166
 167     return basis;
 168 }