basis_reduction_templ.c

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