isl_tab_compute_reduced_basis: allow incremental computation

diff --git a/basis_reduction_templ.c b/basis_reduction_templ.c
index 56aac95..60dcead 100644 (file)
--- a/basis_reduction_templ.c
+++ b/basis_reduction_templ.c
@@ -9,7 +9,15 @@ static void save_alpha(GBR_LP *lp, int first, int n, GBR_type *alpha)
                GBR_lp_get_alpha(lp, first + i, &alpha[i]);
 }
 
-/* This function implements the algorithm described in
+/* Compute a reduced basis for the set represented by the tableau "tab".
+ * tab->basis, must be initialized by the calling function to a unimodular
+ * basis, is updated to reflect the reduced basis.
+ * The first tab->n_zero rows of the basis are assumed to correspond
+ * to equalities and are left untouched.
+ * tab->n_zero is updated to reflect any additional equalities that
+ * have been detected in the first rows of the new basis.
+ *
+ * This function implements the algorithm described in
  * "An Implementation of the Generalized Basis Reduction Algorithm
  *  for Integer Programming" of Cook el al. to compute a reduced basis.
  * We use \epsilon = 1/4.
@@ -18,7 +26,7 @@ static void save_alpha(GBR_LP *lp, int first, int n, GBR_type *alpha)
  * in the first direction.  In this case we stop the basis reduction when
  * the width in the first direction becomes smaller than 2.
  */
-struct isl_mat *isl_tab_reduced_basis(struct isl_tab *tab)
+struct isl_tab *isl_tab_compute_reduced_basis(struct isl_tab *tab)
 {
        unsigned dim;
        struct isl_ctx *ctx;
@@ -39,7 +47,6 @@ struct isl_mat *isl_tab_reduced_basis(struct isl_tab *tab)
        GBR_type mu_F[2];
        GBR_type two;
        GBR_type one;
-       int n_zero = 0;
        int empty = 0;
        int fixed = 0;
        int fixed_saved = 0;
@@ -50,12 +57,12 @@ struct isl_mat *isl_tab_reduced_basis(struct isl_tab *tab)
 
        ctx = tab->mat->ctx;
        dim = tab->n_var;
-       basis = isl_mat_identity(ctx, dim);
+       basis = tab->basis;
        if (!basis)
-               return NULL;
+               return tab;
 
-       if (dim == 1)
-               return basis;
+       if (dim <= tab->n_zero + 1)
+               return tab;
 
        isl_int_init(tmp);
        isl_int_init(mu[0]);
@@ -95,26 +102,26 @@ struct isl_mat *isl_tab_reduced_basis(struct isl_tab *tab)
        if (!lp)
                goto error;
 
-       i = 0;
+       i = tab->n_zero;
 
-       GBR_lp_set_obj(lp, basis->row[0], dim);
+       GBR_lp_set_obj(lp, basis->row[i], dim);
        ctx->stats->gbr_solved_lps++;
        unbounded = GBR_lp_solve(lp);
        isl_assert(ctx, !unbounded, goto error);
-       GBR_lp_get_obj_val(lp, &F[0]);
+       GBR_lp_get_obj_val(lp, &F[i]);
 
-       if (GBR_lt(F[0], one)) {
-               if (!GBR_is_zero(F[0])) {
-                       empty = GBR_lp_cut(lp, basis->row[0]);
+       if (GBR_lt(F[i], one)) {
+               if (!GBR_is_zero(F[i])) {
+                       empty = GBR_lp_cut(lp, basis->row[i]);
                        if (empty)
                                goto done;
-                       GBR_set_ui(F[0], 0);
+                       GBR_set_ui(F[i], 0);
                }
-               n_zero++;
+               tab->n_zero++;
        }
 
        do {
-               if (i+1 == n_zero) {
+               if (i+1 == tab->n_zero) {
                        GBR_lp_set_obj(lp, basis->row[i+1], dim);
                        ctx->stats->gbr_solved_lps++;
                        unbounded = GBR_lp_solve(lp);
@@ -182,14 +189,14 @@ struct isl_mat *isl_tab_reduced_basis(struct isl_tab *tab)
                isl_seq_combine(basis->row[i+1], ctx->one, basis->row[i+1],
                                tmp, basis->row[i], dim);
 
-               if (i+1 == n_zero && fixed) {
+               if (i+1 == tab->n_zero && fixed) {
                        if (!GBR_is_zero(F[i+1])) {
                                empty = GBR_lp_cut(lp, basis->row[i+1]);
                                if (empty)
                                        goto done;
                                GBR_set_ui(F[i+1], 0);
                        }
-                       n_zero++;
+                       tab->n_zero++;
                }
 
                GBR_set(F_old, F[i]);
@@ -202,25 +209,25 @@ struct isl_mat *isl_tab_reduced_basis(struct isl_tab *tab)
                GBR_mul(mu_F[1], mu_F[1], F_old);
                if (GBR_lt(mu_F[0], mu_F[1])) {
                        basis = isl_mat_swap_rows(basis, i, i + 1);
-                       if (i > 0) {
+                       if (i > tab->n_zero) {
                                use_saved = 1;
                                GBR_set(F_saved, F_new);
                                fixed_saved = fixed;
                                GBR_lp_del_row(lp);
                                --i;
                        } else {
-                               GBR_set(F[0], F_new);
-                               if (ctx->gbr_only_first && GBR_lt(F[0], two))
+                               GBR_set(F[tab->n_zero], F_new);
+                               if (ctx->gbr_only_first && GBR_lt(F[tab->n_zero], two))
                                        break;
 
                                if (fixed) {
-                                       if (!GBR_is_zero(F[0])) {
-                                               empty = GBR_lp_cut(lp, basis->row[0]);
+                                       if (!GBR_is_zero(F[tab->n_zero])) {
+                                               empty = GBR_lp_cut(lp, basis->row[tab->n_zero]);
                                                if (empty)
                                                        goto done;
-                                               GBR_set_ui(F[0], 0);
+                                               GBR_set_ui(F[tab->n_zero], 0);
                                        }
-                                       n_zero++;
+                                       tab->n_zero++;
                                }
                        }
                } else {
@@ -265,7 +272,9 @@ error:
        isl_int_clear(mu[0]);
        isl_int_clear(mu[1]);
 
-       return basis;
+       tab->basis = basis;
+
+       return tab;
 }
 
 struct isl_mat *isl_basic_set_reduced_basis(struct isl_basic_set *bset)
@@ -276,7 +285,12 @@ struct isl_mat *isl_basic_set_reduced_basis(struct isl_basic_set *bset)
        isl_assert(bset->ctx, bset->n_eq == 0, return NULL);
 
        tab = isl_tab_from_basic_set(bset);
-       basis = isl_tab_reduced_basis(tab);
+       tab->basis = isl_mat_identity(bset->ctx, tab->n_var);
+       tab = isl_tab_compute_reduced_basis(tab);
+       if (!tab)
+               return NULL;
+
+       basis = isl_mat_copy(tab->basis);
 
        isl_tab_free(tab);
 

diff --git a/isl_basis_reduction.h b/isl_basis_reduction.h
index dffb724..4b1187e 100644 (file)
--- a/isl_basis_reduction.h
+++ b/isl_basis_reduction.h
@@ -9,7 +9,7 @@
 extern "C" {
 #endif
 
-struct isl_mat *isl_tab_reduced_basis(struct isl_tab *tab);
+struct isl_tab *isl_tab_compute_reduced_basis(struct isl_tab *tab);
 struct isl_mat *isl_basic_set_reduced_basis(struct isl_basic_set *bset);
 
 #if defined(__cplusplus)

diff --git a/isl_tab.c b/isl_tab.c
index b697024..84345e0 100644 (file)
--- a/isl_tab.c
+++ b/isl_tab.c
@@ -63,6 +63,10 @@ struct isl_tab *isl_tab_alloc(struct isl_ctx *ctx,
        tab->bottom.type = isl_tab_undo_bottom;
        tab->bottom.next = NULL;
        tab->top = &tab->bottom;
+
+       tab->n_zero = 0;
+       tab->basis = NULL;
+
        return tab;
 error:
        isl_tab_free(tab);
@@ -178,6 +182,7 @@ void isl_tab_free(struct isl_tab *tab)
        free(tab->col_var);
        free(tab->row_sign);
        isl_mat_free(tab->samples);
+       isl_mat_free(tab->basis);
        free(tab);
 }
 
@@ -251,6 +256,10 @@ struct isl_tab *isl_tab_dup(struct isl_tab *tab)
        dup->bottom.type = isl_tab_undo_bottom;
        dup->bottom.next = NULL;
        dup->top = &dup->bottom;
+
+       dup->n_zero = tab->n_zero;
+       dup->basis = isl_mat_dup(tab->basis);
+
        return dup;
 error:
        isl_tab_free(dup);
@@ -493,6 +502,10 @@ struct isl_tab *isl_tab_product(struct isl_tab *tab1, struct isl_tab *tab2)
        prod->bottom.type = isl_tab_undo_bottom;
        prod->bottom.next = NULL;
        prod->top = &prod->bottom;
+
+       prod->n_zero = 0;
+       prod->basis = NULL;
+
        return prod;
 error:
        isl_tab_free(prod);

diff --git a/isl_tab.h b/isl_tab.h
index a0cf789..1bfaae7 100644 (file)
--- a/isl_tab.h
+++ b/isl_tab.h
@@ -138,6 +138,9 @@ struct isl_tab {
        unsigned n_outside;
        struct isl_mat *samples;
 
+       int n_zero;
+       struct isl_mat *basis;
+
        unsigned need_undo : 1;
        unsigned rational : 1;
        unsigned empty : 1;

diff --git a/polytope_scan.c b/polytope_scan.c
index a69ba9c..190d44f 100644 (file)
--- a/polytope_scan.c
+++ b/polytope_scan.c
@@ -94,12 +94,15 @@ static struct isl_mat *isl_basic_set_samples(struct isl_basic_set *bset)
                goto error;
 
        tab = isl_tab_from_basic_set(bset);
+       if (!tab)
+               goto error;
 
+       tab->basis = isl_mat_identity(bset->ctx, dim);
        if (1)
-               B = isl_tab_reduced_basis(tab);
-       else
-               B = isl_mat_identity(bset->ctx, dim);
-       B = isl_mat_lin_to_aff(B);
+               tab = isl_tab_compute_reduced_basis(tab);
+       if (!tab)
+               goto error;
+       B = isl_mat_lin_to_aff(isl_mat_copy(tab->basis));
        if (!B)
                goto error;