+/*
+ * Copyright 2008-2009 Katholieke Universiteit Leuven
+ *
+ * Use of this software is governed by the GNU LGPLv2.1 license
+ *
+ * Written by Sven Verdoolaege, K.U.Leuven, Departement
+ * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
+ */
+
#include <assert.h>
#include <string.h>
#include "isl_set.h"
#include "isl_tab.h"
#include "isl_map_private.h"
+#include "isl_sample.h"
+#include "isl_scan.h"
+#include "isl_seq.h"
+#include "isl_ilp.h"
+#include <isl_point_private.h>
/* The input of this program is the same as that of the "example" program
* from the PipLib distribution, except that the "big parameter column"
* Urs_unknowns don't assume unknowns are non-negative
*/
+struct options {
+ struct isl_options *isl;
+ unsigned verify;
+};
+
+struct isl_arg options_arg[] = {
+ISL_ARG_CHILD(struct options, isl, "isl", isl_options_arg, "isl options")
+ISL_ARG_BOOL(struct options, verify, 'T', "verify", 0, NULL)
+ISL_ARG_END
+};
+
+ISL_ARG_DEF(options, struct options, options_arg)
+
static struct isl_basic_set *to_parameter_domain(struct isl_basic_set *context)
{
struct isl_dim *param_dim;
return context;
}
+static isl_basic_set *construct_cube(isl_basic_set *context)
+{
+ int i;
+ unsigned dim;
+ int range;
+ isl_int m, M;
+ struct isl_basic_set *cube;
+ struct isl_basic_set *interval;
+ struct isl_basic_set_list *list;
+
+ dim = isl_basic_set_total_dim(context);
+ if (dim == 0) {
+ struct isl_dim *dims;
+ dims = isl_dim_set_alloc(context->ctx, 0, 0);
+ return isl_basic_set_universe(dims);
+ }
+
+ if (dim >= 8)
+ range = 4;
+ else if (dim >= 5)
+ range = 6;
+ else
+ range = 30;
+
+ isl_int_init(m);
+ isl_int_init(M);
+
+ isl_int_set_si(m, -range);
+ isl_int_set_si(M, range);
+
+ interval = isl_basic_set_interval(context->ctx, m, M);
+ list = isl_basic_set_list_alloc(context->ctx, dim);
+ for (i = 0; i < dim; ++i)
+ list = isl_basic_set_list_add(list, isl_basic_set_copy(interval));
+ isl_basic_set_free(interval);
+ cube = isl_basic_set_product(list);
+
+ isl_int_clear(m);
+ isl_int_clear(M);
+
+ return cube;
+}
+
+isl_basic_set *plug_in_parameters(isl_basic_set *bset, struct isl_vec *params)
+{
+ int i;
+
+ for (i = 0; i < params->size - 1; ++i)
+ bset = isl_basic_set_fix(bset,
+ isl_dim_param, i, params->el[1 + i]);
+
+ bset = isl_basic_set_remove(bset, isl_dim_param, 0, params->size - 1);
+
+ isl_vec_free(params);
+
+ return bset;
+}
+
+isl_set *set_plug_in_parameters(isl_set *set, struct isl_vec *params)
+{
+ int i;
+
+ for (i = 0; i < params->size - 1; ++i)
+ set = isl_set_fix(set, isl_dim_param, i, params->el[1 + i]);
+
+ set = isl_set_remove(set, isl_dim_param, 0, params->size - 1);
+
+ isl_vec_free(params);
+
+ return set;
+}
+
+/* Compute the lexicographically minimal (or maximal if max is set)
+ * element of bset for the given values of the parameters, by
+ * successively solving an ilp problem in each direction.
+ */
+struct isl_vec *opt_at(struct isl_basic_set *bset,
+ struct isl_vec *params, int max)
+{
+ unsigned dim;
+ struct isl_vec *opt;
+ struct isl_vec *obj;
+ int i;
+
+ dim = isl_basic_set_dim(bset, isl_dim_set);
+
+ bset = plug_in_parameters(bset, params);
+
+ if (isl_basic_set_fast_is_empty(bset)) {
+ opt = isl_vec_alloc(bset->ctx, 0);
+ isl_basic_set_free(bset);
+ return opt;
+ }
+
+ opt = isl_vec_alloc(bset->ctx, 1 + dim);
+ assert(opt);
+
+ obj = isl_vec_alloc(bset->ctx, 1 + dim);
+ assert(obj);
+
+ isl_int_set_si(opt->el[0], 1);
+ isl_int_set_si(obj->el[0], 0);
+
+ for (i = 0; i < dim; ++i) {
+ enum isl_lp_result res;
+
+ isl_seq_clr(obj->el + 1, dim);
+ isl_int_set_si(obj->el[1 + i], 1);
+ res = isl_basic_set_solve_ilp(bset, max, obj->el,
+ &opt->el[1 + i], NULL);
+ if (res == isl_lp_empty)
+ goto empty;
+ assert(res == isl_lp_ok);
+ bset = isl_basic_set_fix(bset, isl_dim_set, i, opt->el[1 + i]);
+ }
+
+ isl_basic_set_free(bset);
+ isl_vec_free(obj);
+
+ return opt;
+empty:
+ isl_vec_free(opt);
+ opt = isl_vec_alloc(bset->ctx, 0);
+ isl_basic_set_free(bset);
+ isl_vec_free(obj);
+
+ return opt;
+}
+
+struct isl_scan_pip {
+ struct isl_scan_callback callback;
+ isl_basic_set *bset;
+ isl_set *sol;
+ isl_set *empty;
+ int stride;
+ int n;
+ int max;
+};
+
+/* Check if the "manually" computed optimum of bset at the "sample"
+ * values of the parameters agrees with the solution of pilp problem
+ * represented by the pair (sol, empty).
+ * In particular, if there is no solution for this value of the parameters,
+ * then it should be an element of the parameter domain "empty".
+ * Otherwise, the optimal solution, should be equal to the result of
+ * plugging in the value of the parameters in "sol".
+ */
+static int scan_one(struct isl_scan_callback *callback,
+ __isl_take isl_vec *sample)
+{
+ struct isl_scan_pip *sp = (struct isl_scan_pip *)callback;
+ struct isl_vec *opt;
+
+ opt = opt_at(isl_basic_set_copy(sp->bset), isl_vec_copy(sample), sp->max);
+ assert(opt);
+
+ if (opt->size == 0) {
+ isl_point *sample_pnt;
+ sample_pnt = isl_point_alloc(isl_set_get_dim(sp->empty), sample);
+ assert(isl_set_contains_point(sp->empty, sample_pnt));
+ isl_point_free(sample_pnt);
+ isl_vec_free(opt);
+ } else {
+ isl_set *sol;
+ isl_set *opt_set;
+ opt_set = isl_set_from_basic_set(isl_basic_set_from_vec(opt));
+ sol = set_plug_in_parameters(isl_set_copy(sp->sol), sample);
+ assert(isl_set_is_equal(opt_set, sol));
+ isl_set_free(sol);
+ isl_set_free(opt_set);
+ }
+
+ if (!(sp->n++ % sp->stride)) {
+ printf("o");
+ fflush(stdout);
+ }
+
+ return 0;
+}
+
+struct isl_scan_count {
+ struct isl_scan_callback callback;
+ int n;
+};
+
+static int count_one(struct isl_scan_callback *callback,
+ __isl_take isl_vec *sample)
+{
+ struct isl_scan_count *sc = (struct isl_scan_count *)callback;
+ isl_vec_free(sample);
+
+ sc->n++;
+
+ return 0;
+}
+
+static void check_solution(isl_basic_set *bset, isl_basic_set *context,
+ isl_set *sol, isl_set *empty, int max)
+{
+ isl_basic_set *cube;
+ struct isl_scan_count sc;
+ struct isl_scan_pip sp;
+ int i;
+ int r;
+
+ context = isl_basic_set_underlying_set(context);
+
+ cube = construct_cube(context);
+ context = isl_basic_set_intersect(context, cube);
+
+ sc.callback.add = count_one;
+ sc.n = 0;
+
+ r = isl_basic_set_scan(isl_basic_set_copy(context), &sc.callback);
+ assert (r == 0);
+
+ sp.callback.add = scan_one;
+ sp.bset = bset;
+ sp.sol = sol;
+ sp.empty = empty;
+ sp.n = 0;
+ sp.stride = sc.n > 70 ? 1 + (sc.n + 1)/70 : 1;
+ sp.max = max;
+
+ for (i = 0; i < sc.n; i += sp.stride)
+ printf(".");
+ printf("\r");
+ fflush(stdout);
+
+ r = isl_basic_set_scan(context, &sp.callback);
+ assert(r == 0);
+
+ printf("\n");
+
+ isl_basic_set_free(bset);
+}
+
int main(int argc, char **argv)
{
- struct isl_ctx *ctx = isl_ctx_alloc();
- struct isl_basic_set *context, *bset;
+ struct isl_ctx *ctx;
+ struct isl_basic_set *context, *bset, *copy, *context_copy;
struct isl_set *set;
struct isl_set *empty;
int neg_one;
int urs_parms = 0;
int urs_unknowns = 0;
int max = 0;
+ int rational = 0;
+ int n;
+ struct options *options;
- context = isl_basic_set_read_from_file(ctx, stdin, 0, ISL_FORMAT_POLYLIB);
+ options = options_new_with_defaults();
+ assert(options);
+ argc = options_parse(options, argc, argv, ISL_ARG_ALL);
+
+ ctx = isl_ctx_alloc_with_options(options_arg, options);
+
+ context = isl_basic_set_read_from_file(ctx, stdin, 0);
assert(context);
- fscanf(stdin, "%d", &neg_one);
+ n = fscanf(stdin, "%d", &neg_one);
+ assert(n == 1);
assert(neg_one == -1);
bset = isl_basic_set_read_from_file(ctx, stdin,
- isl_basic_set_dim(context, isl_dim_set), ISL_FORMAT_POLYLIB);
+ isl_basic_set_dim(context, isl_dim_set));
while (fgets(s, sizeof(s), stdin)) {
if (strncasecmp(s, "Maximize", 8) == 0)
max = 1;
- if (strncasecmp(s, "Rational", 8) == 0)
+ if (strncasecmp(s, "Rational", 8) == 0) {
+ rational = 1;
bset = isl_basic_set_set_rational(bset);
+ }
if (strncasecmp(s, "Urs_parms", 9) == 0)
urs_parms = 1;
if (strncasecmp(s, "Urs_unknowns", 12) == 0)
bset = isl_basic_set_intersect(bset,
isl_basic_set_positive_orthant(isl_basic_set_get_dim(bset)));
+ if (options->verify) {
+ copy = isl_basic_set_copy(bset);
+ context_copy = isl_basic_set_copy(context);
+ }
+
if (max)
set = isl_basic_set_partial_lexmax(bset, context, &empty);
else
set = isl_basic_set_partial_lexmin(bset, context, &empty);
- isl_set_dump(set, stdout, 0);
- fprintf(stdout, "no solution:\n");
- isl_set_dump(empty, stdout, 4);
+ if (options->verify) {
+ assert(!rational);
+ check_solution(copy, context_copy, set, empty, max);
+ } else {
+ isl_set_print(set, stdout, 0, ISL_FORMAT_ISL);
+ fprintf(stdout, "\n");
+ fprintf(stdout, "no solution: ");
+ isl_set_print(empty, stdout, 0, ISL_FORMAT_ISL);
+ fprintf(stdout, "\n");
+ }
isl_set_free(set);
isl_set_free(empty);