isl_tab_lexopt_templ.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  * Copyright 2010      INRIA Saclay
   4  * Copyright 2011      Sven Verdoolaege
   5  * Copyright 2023      Cerebras Systems
   6  *
   7  * Use of this software is governed by the MIT license
   8  *
   9  * Written by Sven Verdoolaege, K.U.Leuven, Departement
  10  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
  11  * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
  12  * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
  13  * and Cerebras Systems, 1237 E Arques Ave, Sunnyvale, CA, USA
  14  */
  15
  16 #define xSF(TYPE,SUFFIX) TYPE ## SUFFIX
  17 #define SF(TYPE,SUFFIX) xSF(TYPE,SUFFIX)
  18
  19 /* Given a basic map with at least two parallel constraints (as found
  20  * by the function parallel_constraints), first look for more constraints
  21  * parallel to the two constraint and replace the found list of parallel
  22  * constraints by a single constraint with as "input" part the minimum
  23  * of the input parts of the list of constraints.  Then, recursively call
  24  * basic_map_partial_lexopt (possibly finding more parallel constraints)
  25  * and plug in the definition of the minimum in the result.
  26  *
  27  * As in parallel_constraints, only inequality constraints that only
  28  * involve input variables that do not occur in any other inequality
  29  * constraints are considered.
  30  *
  31  * More specifically, given a set of constraints
  32  *
  33  *      a x + b_i(p) >= 0
  34  *
  35  * Replace this set by a single constraint
  36  *
  37  *      a x + u >= 0
  38  *
  39  * with u a new parameter with constraints
  40  *
  41  *      u <= b_i(p)
  42  *
  43  * Any solution to the new system is also a solution for the original system
  44  * since
  45  *
  46  *      a x >= -u >= -b_i(p)
  47  *
  48  * Moreover, m = min_i(b_i(p)) satisfies the constraints on u and can
  49  * therefore be plugged into the solution.
  50  */
  51 static TYPE *SF(basic_map_partial_lexopt_symm,SUFFIX)(
  52         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
  53         __isl_give isl_set **empty, int max, int first, int second)
  54 {
  55         int i, n, k;
  56         int *list = NULL;
  57         isl_size bmap_in, bmap_param, bmap_all;
  58         unsigned n_in, n_out, n_div;
  59         isl_ctx *ctx;
  60         isl_vec *var = NULL;
  61         isl_mat *cst = NULL;
  62         isl_space *map_space, *set_space;
  63
  64         map_space = isl_basic_map_get_space(bmap);
  65         set_space = empty ? isl_basic_set_get_space(dom) : NULL;
  66
  67         bmap_in = isl_basic_map_dim(bmap, isl_dim_in);
  68         bmap_param = isl_basic_map_dim(bmap, isl_dim_param);
  69         bmap_all = isl_basic_map_dim(bmap, isl_dim_all);
  70         if (bmap_in < 0 || bmap_param < 0 || bmap_all < 0)
  71                 goto error;
  72         n_in = bmap_param + bmap_in;
  73         n_out = bmap_all - n_in;
  74
  75         ctx = isl_basic_map_get_ctx(bmap);
  76         list = isl_alloc_array(ctx, int, bmap->n_ineq);
  77         var = isl_vec_alloc(ctx, n_out);
  78         if ((bmap->n_ineq && !list) || (n_out && !var))
  79                 goto error;
  80
  81         list[0] = first;
  82         list[1] = second;
  83         isl_seq_cpy(var->el, bmap->ineq[first] + 1 + n_in, n_out);
  84         for (i = second + 1, n = 2; i < bmap->n_ineq; ++i) {
  85                 if (isl_seq_eq(var->el, bmap->ineq[i] + 1 + n_in, n_out) &&
  86                     all_single_occurrence(bmap, i, n_in))
  87                         list[n++] = i;
  88         }
  89
  90         cst = isl_mat_alloc(ctx, n, 1 + n_in);
  91         if (!cst)
  92                 goto error;
  93
  94         for (i = 0; i < n; ++i)
  95                 isl_seq_cpy(cst->row[i], bmap->ineq[list[i]], 1 + n_in);
  96
  97         bmap = isl_basic_map_cow(bmap);
  98         if (!bmap)
  99                 goto error;
 100         for (i = n - 1; i >= 0; --i)
 101                 if (isl_basic_map_drop_inequality(bmap, list[i]) < 0)
 102                         goto error;
 103
 104         bmap = isl_basic_map_add_dims(bmap, isl_dim_in, 1);
 105         bmap = isl_basic_map_extend_constraints(bmap, 0, 1);
 106         k = isl_basic_map_alloc_inequality(bmap);
 107         if (k < 0)
 108                 goto error;
 109         isl_seq_clr(bmap->ineq[k], 1 + n_in);
 110         isl_int_set_si(bmap->ineq[k][1 + n_in], 1);
 111         isl_seq_cpy(bmap->ineq[k] + 1 + n_in + 1, var->el, n_out);
 112         bmap = isl_basic_map_finalize(bmap);
 113
 114         n_div = isl_basic_set_dim(dom, isl_dim_div);
 115         dom = isl_basic_set_add_dims(dom, isl_dim_set, 1);
 116         dom = isl_basic_set_extend_constraints(dom, 0, n);
 117         for (i = 0; i < n; ++i) {
 118                 k = isl_basic_set_alloc_inequality(dom);
 119                 if (k < 0)
 120                         goto error;
 121                 isl_seq_cpy(dom->ineq[k], cst->row[i], 1 + n_in);
 122                 isl_int_set_si(dom->ineq[k][1 + n_in], -1);
 123                 isl_seq_clr(dom->ineq[k] + 1 + n_in + 1, n_div);
 124         }
 125
 126         isl_vec_free(var);
 127         free(list);
 128
 129         return SF(basic_map_partial_lexopt_symm_core,SUFFIX)(bmap, dom, empty,
 130                                                 max, cst, map_space, set_space);
 131 error:
 132         isl_space_free(map_space);
 133         isl_space_free(set_space);
 134         isl_mat_free(cst);
 135         isl_vec_free(var);
 136         free(list);
 137         isl_basic_set_free(dom);
 138         isl_basic_map_free(bmap);
 139         return NULL;
 140 }
 141
 142 static __isl_give TYPE *SF(basic_map_partial_lexopt_intersected,SUFFIX)(
 143         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
 144         __isl_give isl_set **empty, unsigned flags);
 145
 146 /* Given that the output dimension of "bmap" at position "d" is equal to "aff",
 147  * exploit this information to reduce the effective dimensionality of "bmap" and
 148  * then call basic_map_partial_lexopt_intersected recursively.
 149  * "flags" is simply passed along to the recursive call.
 150  * If "flags" includes ISL_OPT_FULL, then "dom" is NULL and
 151  * then also a NULL domain is passed to the recursive call.
 152  *
 153  * In particular, introduce a dimension in the context "dom" (and the domain
 154  * of "bmap") that is equal to "aff" and equate output dimension "d"
 155  * to this new input dimension.
 156  * This essentially moves the output dimension to the input, but
 157  * leaves a placeholder so that the value "aff" can easily be plugged
 158  * into the result of the recursive call.
 159  */
 160 static __isl_give TYPE *SF(basic_map_partial_lexopt_plugin,SUFFIX)(
 161         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
 162         __isl_give isl_set **empty, unsigned flags, int d,
 163         __isl_take isl_aff *aff)
 164 {
 165         isl_size n_in;
 166         isl_multi_aff *ma;
 167         isl_basic_map *insert;
 168         TYPE *res;
 169
 170         n_in = isl_aff_dim(aff, isl_dim_in);
 171         if (n_in < 0)
 172                 bmap = isl_basic_map_free(bmap);
 173
 174         ma = isl_aff_as_domain_extension(aff);
 175         insert = isl_basic_map_from_multi_aff2(isl_multi_aff_copy(ma), 0);
 176
 177         bmap = isl_basic_map_apply_domain(bmap, isl_basic_map_copy(insert));
 178         dom = isl_basic_set_apply(dom, insert);
 179         bmap = isl_basic_map_equate(bmap, isl_dim_in, n_in, isl_dim_out, d);
 180
 181         res = SF(basic_map_partial_lexopt_intersected,SUFFIX)(bmap, dom, empty,
 182                                                                 flags);
 183         if (empty)
 184                 *empty = isl_set_preimage_multi_aff(*empty,
 185                                                 isl_multi_aff_copy(ma));
 186         res = FN(TYPE,pullback_multi_aff)(res, ma);
 187
 188         return res;
 189 }
 190
 191 /* Recursive part of isl_tab_basic_map_partial_lexopt*, after detecting
 192  * equalities and removing redundant constraints.
 193  *
 194  * Check if there are any parallel constraints (left).
 195  * If not, we are in the base case.
 196  * If there are parallel constraints, we replace them by a single
 197  * constraint in basic_map_partial_lexopt_symm_pma and then call
 198  * this function recursively to look for more parallel constraints.
 199  */
 200 static __isl_give TYPE *SF(basic_map_partial_lexopt,SUFFIX)(
 201         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
 202         __isl_give isl_set **empty, int max)
 203 {
 204         isl_bool par = isl_bool_false;
 205         int first, second;
 206         isl_ctx *ctx;
 207
 208         if (!bmap)
 209                 goto error;
 210
 211         ctx = isl_basic_map_get_ctx(bmap);
 212         if (ctx->opt->pip_symmetry)
 213                 par = parallel_constraints(bmap, &first, &second);
 214         if (par < 0)
 215                 goto error;
 216         if (!par)
 217                 return SF(basic_map_partial_lexopt_base,SUFFIX)(bmap, dom,
 218                                                                 empty, max);
 219
 220         return SF(basic_map_partial_lexopt_symm,SUFFIX)(bmap, dom, empty, max,
 221                                                          first, second);
 222 error:
 223         isl_basic_set_free(dom);
 224         isl_basic_map_free(bmap);
 225         return NULL;
 226 }
 227
 228 /* Compute the lexicographic minimum (or maximum if "flags" includes
 229  * ISL_OPT_MAX) of "bmap" over the domain "dom" and return the result as
 230  * either a map or a piecewise multi-affine expression depending on TYPE.
 231  * If "empty" is not NULL, then *empty is assigned a set that
 232  * contains those parts of the domain where there is no solution.
 233  * If "flags" includes ISL_OPT_FULL, then "dom" is NULL and the optimum
 234  * should be computed over the domain of "bmap".  "empty" is also NULL
 235  * in this case.
 236  * All information in "dom" (if any) is assumed to be available in "bmap"
 237  * as well.
 238  * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
 239  * then we compute the rational optimum.  Otherwise, we compute
 240  * the integral optimum.
 241  *
 242  * First check if some combination of constraints can be found that force
 243  * a given dimension to be equal to the floor or modulo
 244  * of some affine combination of the input dimensions.
 245  * If so, plug in this expression and continue.
 246  *
 247  * Otherwise, perform some preprocessing.
 248  * As the PILP solver does not
 249  * handle implicit equalities very well, we first make sure all
 250  * the equalities are explicitly available.
 251  *
 252  * We also remove redundant constraints.  This is only needed because of the
 253  * way we handle simple symmetries.  In particular, we currently look
 254  * for symmetries on the constraints, before we set up the main tableau.
 255  * It is then no good to look for symmetries on possibly redundant constraints.
 256  */
 257 static __isl_give TYPE *SF(basic_map_partial_lexopt_intersected,SUFFIX)(
 258         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
 259         __isl_give isl_set **empty, unsigned flags)
 260 {
 261         int d;
 262         int max;
 263         isl_maybe_isl_aff div_mod;
 264
 265         div_mod = isl_basic_map_try_find_any_output_div_mod(bmap, &d);
 266         if (div_mod.valid < 0)
 267                 bmap = isl_basic_map_free(bmap);
 268         else if (div_mod.valid)
 269                 return SF(basic_map_partial_lexopt_plugin,SUFFIX)(bmap, dom,
 270                                                 empty, flags, d, div_mod.value);
 271
 272         if (empty)
 273                 *empty = NULL;
 274
 275         if (ISL_FL_ISSET(flags, ISL_OPT_FULL))
 276                 dom = extract_domain(bmap, flags);
 277
 278         max = ISL_FL_ISSET(flags, ISL_OPT_MAX);
 279         if (isl_basic_set_dim(dom, isl_dim_all) == 0)
 280                 return SF(basic_map_partial_lexopt,SUFFIX)(bmap, dom, empty,
 281                                                             max);
 282
 283         bmap = isl_basic_map_detect_equalities(bmap);
 284         bmap = isl_basic_map_remove_redundancies(bmap);
 285
 286         return SF(basic_map_partial_lexopt,SUFFIX)(bmap, dom, empty, max);
 287 }
 288
 289 /* Compute the lexicographic minimum (or maximum if "flags" includes
 290  * ISL_OPT_MAX) of "bmap" over the domain "dom" and return the result as
 291  * either a map or a piecewise multi-affine expression depending on TYPE.
 292  * If "empty" is not NULL, then *empty is assigned a set that
 293  * contains those parts of the domain where there is no solution.
 294  * If "flags" includes ISL_OPT_FULL, then "dom" is NULL and the optimum
 295  * should be computed over the domain of "bmap".  "empty" is also NULL
 296  * in this case.
 297  * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
 298  * then we compute the rational optimum.  Otherwise, we compute
 299  * the integral optimum.
 300  *
 301  * Intersect the domain of "bmap" with "dom" (if any)
 302  * to make all information available to "bmap" and
 303  * continue with further processing.
 304  */
 305 __isl_give TYPE *SF(isl_tab_basic_map_partial_lexopt,SUFFIX)(
 306         __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
 307         __isl_give isl_set **empty, unsigned flags)
 308 {
 309         if (!ISL_FL_ISSET(flags, ISL_OPT_FULL))
 310                 bmap = isl_basic_map_intersect_domain(bmap,
 311                                                     isl_basic_set_copy(dom));
 312         return SF(basic_map_partial_lexopt_intersected,SUFFIX)(bmap, dom,
 313                                                                 empty, flags);
 314 }