return sample;
}
-/* Given a tableau that is known to represent a bounded set, find and return
+/* Return a matrix containing the equalities of the tableau
+ * in constraint form. The tableau is assumed to have
+ * an associated bset that has been kept up-to-date.
+ */
+static struct isl_mat *tab_equalities(struct isl_tab *tab)
+{
+ int i, j;
+ int n_eq;
+ struct isl_mat *eq;
+ struct isl_basic_set *bset;
+
+ if (!tab)
+ return NULL;
+
+ isl_assert(tab->mat->ctx, tab->bset, return NULL);
+ bset = tab->bset;
+
+ n_eq = tab->n_var - tab->n_col + tab->n_dead;
+ if (tab->empty || n_eq == 0)
+ return isl_mat_alloc(tab->mat->ctx, 0, tab->n_var);
+ if (n_eq == tab->n_var)
+ return isl_mat_identity(tab->mat->ctx, tab->n_var);
+
+ eq = isl_mat_alloc(tab->mat->ctx, n_eq, tab->n_var);
+ if (!eq)
+ return NULL;
+ for (i = 0, j = 0; i < tab->n_con; ++i) {
+ if (tab->con[i].is_row)
+ continue;
+ if (tab->con[i].index >= 0 && tab->con[i].index >= tab->n_dead)
+ continue;
+ if (i < bset->n_eq)
+ isl_seq_cpy(eq->row[j], bset->eq[i] + 1, tab->n_var);
+ else
+ isl_seq_cpy(eq->row[j],
+ bset->ineq[i - bset->n_eq] + 1, tab->n_var);
+ ++j;
+ }
+ isl_assert(bset->ctx, j == n_eq, goto error);
+ return eq;
+error:
+ isl_mat_free(eq);
+ return NULL;
+}
+
+/* Compute and return an initial basis for the bounded tableau "tab".
+ *
+ * If the tableau is either full-dimensional or zero-dimensional,
+ * the we simply return an identity matrix.
+ * Otherwise, we construct a basis whose first directions correspond
+ * to equalities.
+ */
+static struct isl_mat *initial_basis(struct isl_tab *tab)
+{
+ int n_eq;
+ struct isl_mat *eq;
+ struct isl_mat *Q;
+
+ n_eq = tab->n_var - tab->n_col + tab->n_dead;
+ if (tab->empty || n_eq == 0 || n_eq == tab->n_var)
+ return isl_mat_identity(tab->mat->ctx, 1 + tab->n_var);
+
+ eq = tab_equalities(tab);
+ eq = isl_mat_left_hermite(eq, 0, NULL, &Q);
+ if (!eq)
+ return NULL;
+ isl_mat_free(eq);
+
+ Q = isl_mat_lin_to_aff(Q);
+ return Q;
+}
+
+/* Given a tableau representing a set, find and return
* an integer point in the set, if there is any.
*
* We perform a depth first search
* for an integer point, by scanning all possible values in the range
- * attained by a basis vector, where the initial basis is assumed
- * to have been set by the calling function.
+ * attained by a basis vector, where an initial basis may have been set
+ * by the calling function. Otherwise an initial basis that exploits
+ * the equalities in the tableau is created.
* tab->n_zero is currently ignored and is clobbered by this function.
*
+ * The tableau is allowed to have unbounded direction, but then
+ * the calling function needs to set an initial basis, with the
+ * unbounded directions last and with tab->n_unbounded set
+ * to the number of unbounded directions.
+ * Furthermore, the calling functions needs to add shifted copies
+ * of all constraints involving unbounded directions to ensure
+ * that any feasible rational value in these directions can be rounded
+ * up to yield a feasible integer value.
+ * In particular, let B define the given basis x' = B x
+ * and let T be the inverse of B, i.e., X = T x'.
+ * Let a x + c >= 0 be a constraint of the set represented by the tableau,
+ * or a T x' + c >= 0 in terms of the given basis. Assume that
+ * the bounded directions have an integer value, then we can safely
+ * round up the values for the unbounded directions if we make sure
+ * that x' not only satisfies the original constraint, but also
+ * the constraint "a T x' + c + s >= 0" with s the sum of all
+ * negative values in the last n_unbounded entries of "a T".
+ * The calling function therefore needs to add the constraint
+ * a x + c + s >= 0. The current function then scans the first
+ * directions for an integer value and once those have been found,
+ * it can compute "T ceil(B x)" to yield an integer point in the set.
+ * Note that during the search, the first rows of B may be changed
+ * by a basis reduction, but the last n_unbounded rows of B remain
+ * unaltered and are also not mixed into the first rows.
+ *
* The search is implemented iteratively. "level" identifies the current
* basis vector. "init" is true if we want the first value at the current
* level and false if we want the next value.
if (tab->empty)
return isl_vec_alloc(tab->mat->ctx, 0);
+ if (!tab->basis)
+ tab->basis = initial_basis(tab);
+ if (!tab->basis)
+ return NULL;
+ isl_assert(tab->mat->ctx, tab->basis->n_row == tab->n_var + 1,
+ return NULL);
+ isl_assert(tab->mat->ctx, tab->basis->n_col == tab->n_var + 1,
+ return NULL);
+
ctx = tab->mat->ctx;
dim = tab->n_var;
gbr = ctx->gbr;
- isl_assert(ctx, tab->basis, return NULL);
+ if (tab->n_unbounded == tab->n_var) {
+ sample = isl_tab_get_sample_value(tab);
+ sample = isl_mat_vec_product(isl_mat_copy(tab->basis), sample);
+ sample = isl_vec_ceil(sample);
+ sample = isl_mat_vec_inverse_product(isl_mat_copy(tab->basis),
+ sample);
+ return sample;
+ }
if (isl_tab_extend_cons(tab, dim + 1) < 0)
return NULL;
isl_int_neg(tab->basis->row[1 + level][0], min->el[level]);
tab = isl_tab_add_valid_eq(tab, tab->basis->row[1 + level]);
isl_int_set_si(tab->basis->row[1 + level][0], 0);
- if (level < dim - 1) {
+ if (level + tab->n_unbounded < dim - 1) {
++level;
init = 1;
continue;
break;
}
- if (level >= 0)
+ if (level >= 0) {
sample = isl_tab_get_sample_value(tab);
- else
+ if (!sample)
+ goto error;
+ if (tab->n_unbounded && !isl_int_is_one(sample->el[0])) {
+ sample = isl_mat_vec_product(isl_mat_copy(tab->basis),
+ sample);
+ sample = isl_vec_ceil(sample);
+ sample = isl_mat_vec_inverse_product(
+ isl_mat_copy(tab->basis), sample);
+ }
+ } else
sample = isl_vec_alloc(ctx, 0);
ctx->gbr = gbr;
ctx = bset->ctx;
tab = isl_tab_from_basic_set(bset);
+ if (!ISL_F_ISSET(bset, ISL_BASIC_SET_NO_IMPLICIT))
+ tab = isl_tab_detect_implicit_equalities(tab);
if (!tab)
goto error;
- tab->basis = isl_mat_identity(bset->ctx, 1 + dim);
- if (!tab->basis)
- goto error;
+ tab->bset = isl_basic_set_copy(bset);
sample = isl_tab_sample(tab);
if (!sample)
isl_int_neg(bset->eq[k][0], vec->el[1 + i]);
isl_int_set(bset->eq[k][1 + i], vec->el[0]);
}
- isl_vec_free(vec);
+ bset->sample = vec;
return bset;
error:
projects / platform / upstream / isl.git / blobdiff