isl_tab_lexopt_templ.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  * Copyright 2010      INRIA Saclay
   4  * Copyright 2011      Sven Verdoolaege
   5  *
   6  * Use of this software is governed by the MIT license
   7  *
   8  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   9  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
  10  * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
  11  * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
  12  */
  13
  14 #define xSF(TYPE,SUFFIX) TYPE ## SUFFIX
  15 #define SF(TYPE,SUFFIX) xSF(TYPE,SUFFIX)
  16
  17 /* Given a basic map with at least two parallel constraints (as found
  18  * by the function parallel_constraints), first look for more constraints
  19  * parallel to the two constraint and replace the found list of parallel
  20  * constraints by a single constraint with as "input" part the minimum
  21  * of the input parts of the list of constraints.  Then, recursively call
  22  * basic_map_partial_lexopt (possibly finding more parallel constraints)
  23  * and plug in the definition of the minimum in the result.
  24  *
  25  * As in parallel_constraints, only inequality constraints that only
  26  * involve input variables that do not occur in any other inequality
  27  * constraints are considered.
  28  *
  29  * More specifically, given a set of constraints
  30  *
  31  *      a x + b_i(p) >= 0
  32  *
  33  * Replace this set by a single constraint
  34  *
  35  *      a x + u >= 0
  36  *
  37  * with u a new parameter with constraints
  38  *
  39  *      u <= b_i(p)
  40  *
  41  * Any solution to the new system is also a solution for the original system
  42  * since
  43  *
  44  *      a x >= -u >= -b_i(p)
  45  *
  46  * Moreover, m = min_i(b_i(p)) satisfies the constraints on u and can
  47  * therefore be plugged into the solution.
  48  */
  49 static TYPE *SF(basic_map_partial_lexopt_symm,SUFFIX)(
  50         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
  51         __isl_give isl_set **empty, int max, int first, int second)
  52 {
  53         int i, n, k;
  54         int *list = NULL;
  55         isl_size bmap_in, bmap_param, bmap_all;
  56         unsigned n_in, n_out, n_div;
  57         isl_ctx *ctx;
  58         isl_vec *var = NULL;
  59         isl_mat *cst = NULL;
  60         isl_space *map_space, *set_space;
  61
  62         map_space = isl_basic_map_get_space(bmap);
  63         set_space = empty ? isl_basic_set_get_space(dom) : NULL;
  64
  65         bmap_in = isl_basic_map_dim(bmap, isl_dim_in);
  66         bmap_param = isl_basic_map_dim(bmap, isl_dim_param);
  67         bmap_all = isl_basic_map_dim(bmap, isl_dim_all);
  68         if (bmap_in < 0 || bmap_param < 0 || bmap_all < 0)
  69                 goto error;
  70         n_in = bmap_param + bmap_in;
  71         n_out = bmap_all - n_in;
  72
  73         ctx = isl_basic_map_get_ctx(bmap);
  74         list = isl_alloc_array(ctx, int, bmap->n_ineq);
  75         var = isl_vec_alloc(ctx, n_out);
  76         if ((bmap->n_ineq && !list) || (n_out && !var))
  77                 goto error;
  78
  79         list[0] = first;
  80         list[1] = second;
  81         isl_seq_cpy(var->el, bmap->ineq[first] + 1 + n_in, n_out);
  82         for (i = second + 1, n = 2; i < bmap->n_ineq; ++i) {
  83                 if (isl_seq_eq(var->el, bmap->ineq[i] + 1 + n_in, n_out) &&
  84                     all_single_occurrence(bmap, i, n_in))
  85                         list[n++] = i;
  86         }
  87
  88         cst = isl_mat_alloc(ctx, n, 1 + n_in);
  89         if (!cst)
  90                 goto error;
  91
  92         for (i = 0; i < n; ++i)
  93                 isl_seq_cpy(cst->row[i], bmap->ineq[list[i]], 1 + n_in);
  94
  95         bmap = isl_basic_map_cow(bmap);
  96         if (!bmap)
  97                 goto error;
  98         for (i = n - 1; i >= 0; --i)
  99                 if (isl_basic_map_drop_inequality(bmap, list[i]) < 0)
 100                         goto error;
 101
 102         bmap = isl_basic_map_add_dims(bmap, isl_dim_in, 1);
 103         bmap = isl_basic_map_extend_constraints(bmap, 0, 1);
 104         k = isl_basic_map_alloc_inequality(bmap);
 105         if (k < 0)
 106                 goto error;
 107         isl_seq_clr(bmap->ineq[k], 1 + n_in);
 108         isl_int_set_si(bmap->ineq[k][1 + n_in], 1);
 109         isl_seq_cpy(bmap->ineq[k] + 1 + n_in + 1, var->el, n_out);
 110         bmap = isl_basic_map_finalize(bmap);
 111
 112         n_div = isl_basic_set_dim(dom, isl_dim_div);
 113         dom = isl_basic_set_add_dims(dom, isl_dim_set, 1);
 114         dom = isl_basic_set_extend_constraints(dom, 0, n);
 115         for (i = 0; i < n; ++i) {
 116                 k = isl_basic_set_alloc_inequality(dom);
 117                 if (k < 0)
 118                         goto error;
 119                 isl_seq_cpy(dom->ineq[k], cst->row[i], 1 + n_in);
 120                 isl_int_set_si(dom->ineq[k][1 + n_in], -1);
 121                 isl_seq_clr(dom->ineq[k] + 1 + n_in + 1, n_div);
 122         }
 123
 124         isl_vec_free(var);
 125         free(list);
 126
 127         return SF(basic_map_partial_lexopt_symm_core,SUFFIX)(bmap, dom, empty,
 128                                                 max, cst, map_space, set_space);
 129 error:
 130         isl_space_free(map_space);
 131         isl_space_free(set_space);
 132         isl_mat_free(cst);
 133         isl_vec_free(var);
 134         free(list);
 135         isl_basic_set_free(dom);
 136         isl_basic_map_free(bmap);
 137         return NULL;
 138 }
 139
 140 /* Recursive part of isl_tab_basic_map_partial_lexopt*, after detecting
 141  * equalities and removing redundant constraints.
 142  *
 143  * We first check if there are any parallel constraints (left).
 144  * If not, we are in the base case.
 145  * If there are parallel constraints, we replace them by a single
 146  * constraint in basic_map_partial_lexopt_symm_pma and then call
 147  * this function recursively to look for more parallel constraints.
 148  */
 149 static __isl_give TYPE *SF(basic_map_partial_lexopt,SUFFIX)(
 150         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
 151         __isl_give isl_set **empty, int max)
 152 {
 153         isl_bool par = isl_bool_false;
 154         int first, second;
 155         isl_ctx *ctx;
 156
 157         if (!bmap)
 158                 goto error;
 159
 160         ctx = isl_basic_map_get_ctx(bmap);
 161         if (ctx->opt->pip_symmetry)
 162                 par = parallel_constraints(bmap, &first, &second);
 163         if (par < 0)
 164                 goto error;
 165         if (!par)
 166                 return SF(basic_map_partial_lexopt_base,SUFFIX)(bmap, dom,
 167                                                                 empty, max);
 168
 169         return SF(basic_map_partial_lexopt_symm,SUFFIX)(bmap, dom, empty, max,
 170                                                          first, second);
 171 error:
 172         isl_basic_set_free(dom);
 173         isl_basic_map_free(bmap);
 174         return NULL;
 175 }
 176
 177 /* Compute the lexicographic minimum (or maximum if "flags" includes
 178  * ISL_OPT_MAX) of "bmap" over the domain "dom" and return the result as
 179  * either a map or a piecewise multi-affine expression depending on TYPE.
 180  * If "empty" is not NULL, then *empty is assigned a set that
 181  * contains those parts of the domain where there is no solution.
 182  * If "flags" includes ISL_OPT_FULL, then "dom" is NULL and the optimum
 183  * should be computed over the domain of "bmap".  "empty" is also NULL
 184  * in this case.
 185  * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
 186  * then we compute the rational optimum.  Otherwise, we compute
 187  * the integral optimum.
 188  *
 189  * We perform some preprocessing.  As the PILP solver does not
 190  * handle implicit equalities very well, we first make sure all
 191  * the equalities are explicitly available.
 192  *
 193  * We also add context constraints to the basic map and remove
 194  * redundant constraints.  This is only needed because of the
 195  * way we handle simple symmetries.  In particular, we currently look
 196  * for symmetries on the constraints, before we set up the main tableau.
 197  * It is then no good to look for symmetries on possibly redundant constraints.
 198  * If the domain was extracted from the basic map, then there is
 199  * no need to add back those constraints again.
 200  */
 201 __isl_give TYPE *SF(isl_tab_basic_map_partial_lexopt,SUFFIX)(
 202         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
 203         __isl_give isl_set **empty, unsigned flags)
 204 {
 205         int max, full;
 206         isl_bool compatible;
 207
 208         if (empty)
 209                 *empty = NULL;
 210
 211         full = ISL_FL_ISSET(flags, ISL_OPT_FULL);
 212         if (full)
 213                 dom = extract_domain(bmap, flags);
 214         compatible = isl_basic_map_compatible_domain(bmap, dom);
 215         if (compatible < 0)
 216                 goto error;
 217         if (!compatible)
 218                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
 219                         "domain does not match input", goto error);
 220
 221         max = ISL_FL_ISSET(flags, ISL_OPT_MAX);
 222         if (isl_basic_set_dim(dom, isl_dim_all) == 0)
 223                 return SF(basic_map_partial_lexopt,SUFFIX)(bmap, dom, empty,
 224                                                             max);
 225
 226         if (!full)
 227                 bmap = isl_basic_map_intersect_domain(bmap,
 228                                                     isl_basic_set_copy(dom));
 229         bmap = isl_basic_map_detect_equalities(bmap);
 230         bmap = isl_basic_map_remove_redundancies(bmap);
 231
 232         return SF(basic_map_partial_lexopt,SUFFIX)(bmap, dom, empty, max);
 233 error:
 234         isl_basic_set_free(dom);
 235         isl_basic_map_free(bmap);
 236         return NULL;
 237 }