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
  46     if (P->Dimension == 1)
  47         return basis;
  48
  49     value_init(one);
  50     value_init(tmp);
  51     value_set_si(one, 1);
  52     value_init(mu[0]);
  53     value_init(mu[1]);
  54
  55     b_tmp = Vector_Alloc(dim);
  56
  57     F = ALLOCN(GBR_type, dim);
  58     alpha_buffer[0] = ALLOCN(GBR_type, dim);
  59     alpha_buffer[1] = ALLOCN(GBR_type, dim);
  60     alpha_saved = alpha_buffer[0];
  61
  62     for (i = 0; i < dim; ++i) {
  63         GBR_init(F[i]);
  64         GBR_init(alpha_buffer[0][i]);
  65         GBR_init(alpha_buffer[1][i]);
  66     }
  67     GBR_init(alpha);
  68     GBR_init(F_old);
  69     GBR_init(F_new);
  70     GBR_init(F_saved);
  71     GBR_init(mu_F[0]);
  72     GBR_init(mu_F[1]);
  73     GBR_init(two);
  74
  75     GBR_set_ui(two, 2);
  76
  77     lp = GBR_lp_init(P);
  78
  79     i = 0;
  80
  81     GBR_lp_set_obj(lp, basis->p[0], dim);
  82     options->stats->gbr_solved_lps++;
  83     if (GBR_lp_solve(lp))
  84         goto unbounded;
  85     GBR_lp_get_obj_val(lp, &F[0]);
  86
  87     do {
  88         if (use_saved) {
  89             row = GBR_lp_next_row(lp);
  90             GBR_set(F_new, F_saved);
  91             GBR_set(alpha, alpha_saved[i]);
  92         } else {
  93             row = GBR_lp_add_row(lp, basis->p[i], dim);
  94             GBR_lp_set_obj(lp, basis->p[i+1], dim);
  95             options->stats->gbr_solved_lps++;
  96             if (GBR_lp_solve(lp))
  97                 goto unbounded;
  98             GBR_lp_get_obj_val(lp, &F_new);
  99
 100             GBR_lp_get_alpha(lp, row, &alpha);
 101
 102             if (i > 0)
 103                 save_alpha(lp, row-i, i, alpha_saved);
 104
 105             GBR_lp_del_row(lp);
 106         }
 107         GBR_set(F[i+1], F_new);
 108
 109         GBR_floor(mu[0], alpha);
 110         GBR_ceil(mu[1], alpha);
 111
 112         if (value_eq(mu[0], mu[1]))
 113             value_assign(tmp, mu[0]);
 114         else {
 115             int j;
 116
 117             for (j = 0; j <= 1; ++j) {
 118                 value_assign(tmp, mu[j]);
 119                 Vector_Combine(basis->p[i+1], basis->p[i], b_tmp->p, one, tmp, dim);
 120                 GBR_lp_set_obj(lp, b_tmp->p, dim);
 121                 options->stats->gbr_solved_lps++;
 122                 if (GBR_lp_solve(lp))
 123                     goto unbounded;
 124                 GBR_lp_get_obj_val(lp, &mu_F[j]);
 125                 if (i > 0)
 126                     save_alpha(lp, row-i, i, alpha_buffer[j]);
 127             }
 128
 129             if (GBR_lt(mu_F[0], mu_F[1]))
 130                 j = 0;
 131             else
 132                 j = 1;
 133
 134             value_assign(tmp, mu[j]);
 135             GBR_set(F_new, mu_F[j]);
 136             alpha_saved = alpha_buffer[j];
 137         }
 138         Vector_Combine(basis->p[i+1], basis->p[i], basis->p[i+1], one, tmp, dim);
 139
 140         GBR_set(F_old, F[i]);
 141
 142         use_saved = 0;
 143         /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
 144         GBR_set_ui(mu_F[0], 4);
 145         GBR_mul(mu_F[0], mu_F[0], F_new);
 146         GBR_set_ui(mu_F[1], 3);
 147         GBR_mul(mu_F[1], mu_F[1], F_old);
 148         if (GBR_lt(mu_F[0], mu_F[1])) {
 149             Vector_Exchange(basis->p[i], basis->p[i+1], dim);
 150             if (i > 0) {
 151                 use_saved = 1;
 152                 GBR_set(F_saved, F_new);
 153                 GBR_lp_del_row(lp);
 154                 --i;
 155             } else {
 156                 GBR_set(F[0], F_new);
 157                 if (options->gbr_only_first && GBR_lt(F[0], two))
 158                     break;
 159             }
 160         } else {
 161             GBR_lp_add_row(lp, basis->p[i], dim);
 162             ++i;
 163         }
 164     } while (i < dim-1);
 165
 166     if (0) {
 167 unbounded:
 168         Matrix_Free(basis);
 169         basis = NULL;
 170     }
 171     Vector_Free(b_tmp);
 172
 173     value_clear(one);
 174     value_clear(tmp);
 175     value_clear(mu[0]);
 176     value_clear(mu[1]);
 177     for (i = 0; i < dim; ++i) {
 178         GBR_clear(F[i]);
 179         GBR_clear(alpha_buffer[0][i]);
 180         GBR_clear(alpha_buffer[1][i]);
 181     }
 182     free(F);
 183     free(alpha_buffer[0]);
 184     free(alpha_buffer[1]);
 185     GBR_clear(alpha);
 186     GBR_clear(F_old);
 187     GBR_clear(F_new);
 188     GBR_clear(F_saved);
 189     GBR_clear(mu_F[0]);
 190     GBR_clear(mu_F[1]);
 191     GBR_clear(two);
 192
 193     GBR_lp_delete(lp);
 194
 195     return basis;
 196 }