basis_reduction_templ.c

   1 #include <stdlib.h>
   2 #include <barvinok/basis_reduction.h>
   3 #include <barvinok/options.h>
   4
   5 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   6
   7 static void save_alpha(GBR_LP *lp, int first, int n, GBR_type *alpha)
   8 {
   9     int i;
  10
  11     for (i = 0; i < n; ++i)
  12         GBR_lp_get_alpha(lp, first+i, &alpha[i]);
  13 }
  14
  15 /* This function implements the algorithm described in
  16  * "An Implementation of the Generalized Basis Reduction Algorithm
  17  *  for Integer Programming" of Cook el al. to compute a reduced basis.
  18  * We use \epsilon = 1/4.
  19  *
  20  * If options->gbr_only_first is set, the user is only interested
  21  * in the first direction.  In this case we stop the basis reduction when
  22  *      - the width in the first direction becomes smaller than 2
  23  * or
  24  *      - we have moved forward all the way to the last direction
  25  *        and then back again all the way to the first direction.
  26  */
  27 Matrix *Polyhedron_Reduced_Basis(Polyhedron *P,
  28                                  struct barvinok_options *options)
  29 {
  30     int dim = P->Dimension;
  31     int i;
  32     Matrix *basis = Identity(dim);
  33     GBR_LP *lp;
  34     GBR_type F_old, alpha, F_new;
  35     int row;
  36     Value one, tmp;
  37     Vector *b_tmp;
  38     GBR_type *F;
  39     GBR_type *alpha_buffer[2];
  40     GBR_type *alpha_saved;
  41     GBR_type F_saved;
  42     int use_saved = 0;
  43     Value mu[2];
  44     GBR_type mu_F[2];
  45     GBR_type two;
  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             Vector_Exchange(basis->p[i], basis->p[i+1], dim);
 151             if (i > 0) {
 152                 use_saved = 1;
 153                 GBR_set(F_saved, F_new);
 154                 GBR_lp_del_row(lp);
 155                 --i;
 156             } else {
 157                 GBR_set(F[0], F_new);
 158                 if (options->gbr_only_first && GBR_lt(F[0], two))
 159                     break;
 160             }
 161         } else {
 162             GBR_lp_add_row(lp, basis->p[i], dim);
 163             ++i;
 164         }
 165     } while (i < dim-1);
 166
 167     if (0) {
 168 unbounded:
 169         Matrix_Free(basis);
 170         basis = NULL;
 171     }
 172     Vector_Free(b_tmp);
 173
 174     value_clear(one);
 175     value_clear(tmp);
 176     value_clear(mu[0]);
 177     value_clear(mu[1]);
 178     for (i = 0; i < dim; ++i) {
 179         GBR_clear(F[i]);
 180         GBR_clear(alpha_buffer[0][i]);
 181         GBR_clear(alpha_buffer[1][i]);
 182     }
 183     free(F);
 184     free(alpha_buffer[0]);
 185     free(alpha_buffer[1]);
 186     GBR_clear(alpha);
 187     GBR_clear(F_old);
 188     GBR_clear(F_new);
 189     GBR_clear(F_saved);
 190     GBR_clear(mu_F[0]);
 191     GBR_clear(mu_F[1]);
 192     GBR_clear(two);
 193
 194     GBR_lp_delete(lp);
 195
 196     return basis;
 197 }