diff options
author | Sven Verdoolaege <skimo@kotnet.org> | 2009-10-04 16:38:28 +0200 |
---|---|---|
committer | Sven Verdoolaege <skimo@kotnet.org> | 2009-10-10 13:16:33 +0200 |
commit | ad855c5588377065aca217b8d5bec4c34fc2d53c (patch) | |
tree | 5f01c314449ecb3db25857869a34cc0797b91185 /isl_sample.c | |
parent | 18933fc04fa4fb91a2aac1875c040c9bb14049fd (diff) | |
download | isl-ad855c5588377065aca217b8d5bec4c34fc2d53c.tar.gz isl-ad855c5588377065aca217b8d5bec4c34fc2d53c.tar.bz2 isl-ad855c5588377065aca217b8d5bec4c34fc2d53c.zip |
add isl_tab_set_initial_basis_with_cone
Diffstat (limited to 'isl_sample.c')
-rw-r--r-- | isl_sample.c | 126 |
1 files changed, 126 insertions, 0 deletions
diff --git a/isl_sample.c b/isl_sample.c index 35fc95e7..b89a3a2e 100644 --- a/isl_sample.c +++ b/isl_sample.c @@ -910,6 +910,132 @@ error: return NULL; } +static void vec_sum_of_neg(struct isl_vec *v, isl_int *s) +{ + int i; + + isl_int_set_si(*s, 0); + + for (i = 0; i < v->size; ++i) + if (isl_int_is_neg(v->el[i])) + isl_int_add(*s, *s, v->el[i]); +} + +/* Given a tableau "tab", a tableau "tab_cone" that corresponds + * to the recession cone and the inverse of a new basis U = inv(B), + * with the unbounded directions in B last, + * add constraints to "tab" that ensure any rational value + * in the unbounded directions can be rounded up to an integer value. + * + * The new basis is given by x' = B x, i.e., x = U x'. + * For any rational value of the last tab->n_unbounded coordinates + * in the update tableau, the value that is obtained by rounding + * up this value should be contained in the original tableau. + * For any constraint "a x + c >= 0", we therefore need to add + * a constraint "a x + c + s >= 0", with s the sum of all negative + * entries in the last elements of "a U". + * + * Since we are not interested in the first entries of any of the "a U", + * we first drop the columns of U that correpond to bounded directions. + */ +static int tab_shift_cone(struct isl_tab *tab, + struct isl_tab *tab_cone, struct isl_mat *U) +{ + int i; + isl_int v; + struct isl_basic_set *bset = NULL; + + if (tab && tab->n_unbounded == 0) { + isl_mat_free(U); + return 0; + } + isl_int_init(v); + if (!tab || !tab_cone || !U) + goto error; + bset = tab_cone->bset; + U = isl_mat_drop_cols(U, 0, tab->n_var - tab->n_unbounded); + for (i = 0; i < bset->n_ineq; ++i) { + struct isl_vec *row = NULL; + if (isl_tab_is_equality(tab_cone, tab_cone->n_eq + i)) + continue; + row = isl_vec_alloc(bset->ctx, tab_cone->n_var); + if (!row) + goto error; + isl_seq_cpy(row->el, bset->ineq[i] + 1, tab_cone->n_var); + row = isl_vec_mat_product(row, isl_mat_copy(U)); + if (!row) + goto error; + vec_sum_of_neg(row, &v); + isl_vec_free(row); + if (isl_int_is_zero(v)) + continue; + tab = isl_tab_extend(tab, 1); + isl_int_add(bset->ineq[i][0], bset->ineq[i][0], v); + tab = isl_tab_add_ineq(tab, bset->ineq[i]); + isl_int_sub(bset->ineq[i][0], bset->ineq[i][0], v); + if (!tab) + goto error; + } + + isl_mat_free(U); + isl_int_clear(v); + return 0; +error: + isl_mat_free(U); + isl_int_clear(v); + return -1; +} + +/* Compute and return an initial basis for the possibly + * unbounded tableau "tab". "tab_cone" is a tableau + * for the corresponding recession cone. + * Additionally, add constraints to "tab" that ensure + * that any rational value for the unbounded directions + * can be rounded up to an integer value. + * + * If the tableau is bounded, i.e., if the recession cone + * is zero-dimensional, then we just use inital_basis. + * Otherwise, we construct a basis whose first directions + * correspond to equalities, followed by bounded directions, + * i.e., equalities in the recession cone. + * The remaining directions are then unbounded. + */ +int isl_tab_set_initial_basis_with_cone(struct isl_tab *tab, + struct isl_tab *tab_cone) +{ + struct isl_mat *eq; + struct isl_mat *cone_eq; + struct isl_mat *U, *Q; + + if (!tab || !tab_cone) + return -1; + + if (tab_cone->n_col == tab_cone->n_dead) { + tab->basis = initial_basis(tab); + return tab->basis ? 0 : -1; + } + + eq = tab_equalities(tab); + if (!eq) + return -1; + tab->n_zero = eq->n_row; + cone_eq = tab_equalities(tab_cone); + eq = isl_mat_concat(eq, cone_eq); + if (!eq) + return -1; + tab->n_unbounded = tab->n_var - (eq->n_row - tab->n_zero); + eq = isl_mat_left_hermite(eq, 0, &U, &Q); + if (!eq) + return -1; + isl_mat_free(eq); + tab->basis = isl_mat_lin_to_aff(Q); + if (tab_shift_cone(tab, tab_cone, U) < 0) + return -1; + if (!tab->basis) + return -1; + return 0; +} + /* Compute and return a sample point in bset using generalized basis * reduction. We first check if the input set has a non-trivial * recession cone. If so, we perform some extra preprocessing in |