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;
 108     double *alpha_buffer[2];
 109     double *alpha_saved;
 110     double F_saved;
 111     int use_saved = 0;
 112
 113     if (P->Dimension == 1)
 114         return basis;
 115
 116     value_init(one);
 117     value_init(tmp);
 118     value_set_si(one, 1);
 119
 120     b_tmp = Vector_Alloc(dim);
 121
 122     F = ALLOCN(double, dim);
 123     alpha_buffer[0] = ALLOCN(double, dim);
 124     alpha_buffer[1] = ALLOCN(double, dim);
 125     alpha_saved = alpha_buffer[0];
 126
 127     lp = init_lp(P);
 128
 129     lpx_set_int_parm(lp, LPX_K_MSGLEV, 0);
 130
 131     i = 0;
 132
 133     set_lp_obj(lp, basis->p[0], dim);
 134
 135     lpx_adv_basis(lp);
 136     lpx_simplex(lp);
 137     F[0] = lpx_get_obj_val(lp);
 138     assert(F[0] > -1e-10);
 139
 140     do {
 141         int mu[2];
 142
 143         if (use_saved) {
 144             row = lpx_get_num_rows(lp)+1;
 145             F_new = F_saved;
 146             alpha = alpha_saved[i];
 147         } else {
 148             row = add_lp_row(lp, basis->p[i], dim);
 149             set_lp_obj(lp, basis->p[i+1], dim);
 150             lpx_adv_basis(lp);
 151             lpx_simplex(lp);
 152             F_new = lpx_get_obj_val(lp);
 153
 154             alpha = -lpx_get_row_dual(lp, row);
 155
 156             if (i > 0)
 157                 save_alpha(lp, row-i, i, alpha_saved);
 158
 159             del_lp_row(lp);
 160         }
 161         F[i+1] = F_new;
 162         assert(F[i+1] > -1e-10);
 163
 164
 165         mu[0] = (int)floor(alpha+1e-10);
 166         mu[1] = (int)ceil(alpha-1e-10);
 167
 168         if (mu[0] == mu[1])
 169             value_set_si(tmp, mu[0]);
 170         else {
 171             double mu_F[2];
 172             int j;
 173
 174             for (j = 0; j <= 1; ++j) {
 175                 value_set_si(tmp, mu[j]);
 176                 Vector_Combine(basis->p[i+1], basis->p[i], b_tmp->p, one, tmp, dim);
 177                 set_lp_obj(lp, b_tmp->p, dim);
 178                 lpx_adv_basis(lp);
 179                 lpx_simplex(lp);
 180                 mu_F[j] = lpx_get_obj_val(lp);
 181                 if (i > 0)
 182                     save_alpha(lp, row-i, i, alpha_buffer[j]);
 183             }
 184
 185             if (mu_F[0] < mu_F[1])
 186                 j = 0;
 187             else
 188                 j = 1;
 189
 190             value_set_si(tmp, mu[j]);
 191             F_new = mu_F[j];
 192             alpha_saved = alpha_buffer[j];
 193         }
 194         Vector_Combine(basis->p[i+1], basis->p[i], basis->p[i+1], one, tmp, dim);
 195
 196         assert(F_new > -1e-10);
 197         F_old = F[i];
 198
 199         use_saved = 0;
 200         if (F_new < 3*F_old/4) {
 201             Vector_Exchange(basis->p[i], basis->p[i+1], dim);
 202             if (i > 0) {
 203                 use_saved = 1;
 204                 F_saved = F_new;
 205                 del_lp_row(lp);
 206                 --i;
 207             } else
 208                 F[0] = F_new;
 209         } else {
 210             add_lp_row(lp, basis->p[i], dim);
 211             ++i;
 212         }
 213     } while (i < dim-1);
 214
 215     Vector_Free(b_tmp);
 216
 217     value_clear(one);
 218     value_clear(tmp);
 219     free(F);
 220     free(alpha_buffer[0]);
 221     free(alpha_buffer[1]);
 222
 223     return basis;
 224 }