basis_reduction.c

   1 #include <math.h>
   2 #include <glpk.h>
   3 #include <polylib/polylibgmp.h>
   4
   5 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   6
   7 static LPX *init_lp(Polyhedron *P)
   8 {
   9     LPX *lp;
  10     int *ind;
  11     double *val;
  12     int i, j, k, l;
  13
  14     ind = ALLOCN(int, 1+P->Dimension);
  15     val = ALLOCN(double, 1+P->Dimension);
  16     lp = lpx_create_prob();
  17     lpx_set_obj_dir(lp, LPX_MAX);
  18     lpx_add_rows(lp, 2*P->NbConstraints);
  19     lpx_add_cols(lp, 2*P->Dimension);
  20     for (i = 0; i < 2; ++i) {
  21         for (j = 0; j < P->NbConstraints; ++j) {
  22             for (k = 0, l = 0; k < P->Dimension; ++k) {
  23                 if (value_zero_p(P->Constraint[j][1+k]))
  24                     continue;
  25                 ind[1+l] = 1+i*P->Dimension+k;
  26                 val[1+l] = VALUE_TO_DOUBLE(P->Constraint[j][1+k]);
  27                 ++l;
  28             }
  29             lpx_set_mat_row(lp, 1+i*P->NbConstraints+j, l, ind, val);
  30             lpx_set_row_bnds(lp, 1+i*P->NbConstraints+j, LPX_LO,
  31                              -VALUE_TO_DOUBLE(P->Constraint[j][1+P->Dimension]), 0);
  32         }
  33         for (k = 0, l = 0; k < P->Dimension; ++k) {
  34             lpx_set_col_bnds(lp, 1+i*P->Dimension+k, LPX_FR, 0, 0);
  35         }
  36     }
  37     free(ind);
  38     free(val);
  39     return lp;
  40 }
  41
  42 static void set_lp_obj(LPX *lp, Value *row, int dim)
  43 {
  44     int j;
  45     for (j = 0; j < dim; ++j) {
  46         lpx_set_obj_coef(lp, 1+j, VALUE_TO_DOUBLE(row[j]));
  47         lpx_set_obj_coef(lp, 1+dim+j, -VALUE_TO_DOUBLE(row[j]));
  48     }
  49 }
  50
  51 static int add_lp_row(LPX *lp, Value *row, int dim)
  52 {
  53     int j, l;
  54     int nr = lpx_add_rows(lp, 1);
  55     int *ind;
  56     double *val;
  57
  58     ind = ALLOCN(int, 1+2*dim);
  59     val = ALLOCN(double, 1+2*dim);
  60     for (j = 0, l = 0; j < dim; ++j) {
  61         if (value_zero_p(row[j]))
  62             continue;
  63         ind[1+l] = 1+j;
  64         val[1+l] = VALUE_TO_DOUBLE(row[j]);
  65         ind[1+l+1] = 1+dim+j;
  66         val[1+l+1] = -VALUE_TO_DOUBLE(row[j]);
  67         l += 2;
  68     }
  69     lpx_set_mat_row(lp, nr, l, ind, val);
  70     lpx_set_row_bnds(lp, nr, LPX_FX, 0, 0);
  71     free(ind);
  72     free(val);
  73
  74     return nr;
  75 }
  76
  77 static void del_lp_row(LPX *lp)
  78 {
  79     int rows[2];
  80     rows[1] = lpx_get_num_rows(lp);
  81     lpx_del_rows(lp, 1, rows);
  82 }
  83
  84 static void save_alpha(LPX *lp, int first, int n, double *alpha)
  85 {
  86     int i;
  87
  88     for (i = 0; i < n; ++i)
  89         alpha[i] = -lpx_get_row_dual(lp, first+i);
  90 }
  91
  92 /* This function implements the algorithm described in
  93  * "An Implementation of the Generalized Basis Reduction Algorithm
  94  *  for Integer Programming" of Cook el al. to compute a reduced basis.
  95  * We use \epsilon = 1/4.
  96  */
  97 Matrix *reduced_basis(Polyhedron *P)
  98 {
  99     int dim = P->Dimension;
 100     int i;
 101     Matrix *basis = Identity(dim);
 102     LPX *lp;
 103     double F_old, alpha, F_new;
 104     int row;
 105     Value one, tmp;
 106     Vector *b_tmp;
 107     double *F = ALLOCN(double, dim);
 108     double *alpha_buffer[2];
 109     double *alpha_saved;
 110     double F_saved;
 111     int use_saved = 0;
 112
 113     value_init(one);
 114     value_init(tmp);
 115     value_set_si(one, 1);
 116
 117     b_tmp = Vector_Alloc(dim);
 118
 119     alpha_buffer[0] = ALLOCN(double, dim);
 120     alpha_buffer[1] = ALLOCN(double, dim);
 121     alpha_saved = alpha_buffer[0];
 122
 123     lp = init_lp(P);
 124
 125     lpx_set_int_parm(lp, LPX_K_MSGLEV, 0);
 126
 127     i = 0;
 128
 129     set_lp_obj(lp, basis->p[0], dim);
 130
 131     lpx_adv_basis(lp);
 132     lpx_simplex(lp);
 133     F[0] = lpx_get_obj_val(lp);
 134
 135     do {
 136         int mu[2];
 137
 138         if (use_saved) {
 139             row = lpx_get_num_rows(lp)+1;
 140             F_new = F_saved;
 141             alpha = alpha_saved[i];
 142         } else {
 143             row = add_lp_row(lp, basis->p[i], dim);
 144             set_lp_obj(lp, basis->p[i+1], dim);
 145             lpx_adv_basis(lp);
 146             lpx_simplex(lp);
 147             F_new = lpx_get_obj_val(lp);
 148
 149             alpha = -lpx_get_row_dual(lp, row);
 150
 151             if (i > 0)
 152                 save_alpha(lp, row-i, i, alpha_saved);
 153
 154             del_lp_row(lp);
 155         }
 156         F[i+1] = F_new;
 157
 158
 159         mu[0] = (int)floor(alpha+1e-10);
 160         mu[1] = (int)ceil(alpha-1e-10);
 161
 162         if (mu[0] == mu[1])
 163             value_set_si(tmp, mu[0]);
 164         else {
 165             double mu_F[2];
 166             int j;
 167
 168             for (j = 0; j <= 1; ++j) {
 169                 value_set_si(tmp, mu[j]);
 170                 Vector_Combine(basis->p[i+1], basis->p[i], b_tmp->p, one, tmp, dim);
 171                 set_lp_obj(lp, b_tmp->p, dim);
 172                 lpx_adv_basis(lp);
 173                 lpx_simplex(lp);
 174                 mu_F[j] = lpx_get_obj_val(lp);
 175                 if (i > 0)
 176                     save_alpha(lp, row-i, i, alpha_buffer[j]);
 177             }
 178
 179             if (mu_F[0] < mu_F[1])
 180                 j = 0;
 181             else
 182                 j = 1;
 183
 184             value_set_si(tmp, mu[j]);
 185             F_new = mu_F[j];
 186             alpha_saved = alpha_buffer[j];
 187         }
 188         Vector_Combine(basis->p[i+1], basis->p[i], basis->p[i+1], one, tmp, dim);
 189
 190         F_old = F[i];
 191
 192         use_saved = 0;
 193         if (F_new < 3*F_old/4) {
 194             Vector_Exchange(basis->p[i], basis->p[i+1], dim);
 195             if (i > 0) {
 196                 use_saved = 1;
 197                 F_saved = F_new;
 198                 del_lp_row(lp);
 199                 --i;
 200             } else
 201                 F[0] = F_new;
 202         } else {
 203             add_lp_row(lp, basis->p[i], dim);
 204             ++i;
 205         }
 206     } while (i < dim-1);
 207
 208     Vector_Free(b_tmp);
 209
 210     value_clear(one);
 211     value_clear(tmp);
 212     free(F);
 213     free(alpha_buffer[0]);
 214     free(alpha_buffer[1]);
 215
 216     return basis;
 217 }