236 lines · c
1/*2 * Copyright 2008-2009 Katholieke Universiteit Leuven3 * Copyright 2010 INRIA Saclay4 * Copyright 2011 Sven Verdoolaege5 *6 * Use of this software is governed by the MIT license7 *8 * Written by Sven Verdoolaege, K.U.Leuven, Departement9 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium10 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,11 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France12 */13 14#define xSF(TYPE,SUFFIX) TYPE ## SUFFIX15#define SF(TYPE,SUFFIX) xSF(TYPE,SUFFIX)16 17/* Given a basic map with at least two parallel constraints (as found18 * by the function parallel_constraints), first look for more constraints19 * parallel to the two constraint and replace the found list of parallel20 * constraints by a single constraint with as "input" part the minimum21 * of the input parts of the list of constraints. Then, recursively call22 * 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 only26 * involve input variables that do not occur in any other inequality27 * constraints are considered.28 *29 * More specifically, given a set of constraints30 *31 * a x + b_i(p) >= 032 *33 * Replace this set by a single constraint34 *35 * a x + u >= 036 *37 * with u a new parameter with constraints38 *39 * u <= b_i(p)40 *41 * Any solution to the new system is also a solution for the original system42 * since43 *44 * a x >= -u >= -b_i(p)45 *46 * Moreover, m = min_i(b_i(p)) satisfies the constraints on u and can47 * therefore be plugged into the solution.48 */49static 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);129error: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 detecting141 * 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 single146 * constraint in basic_map_partial_lexopt_symm_pma and then call147 * this function recursively to look for more parallel constraints.148 */149static __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 156 if (!bmap)157 goto error;158 159 if (bmap->ctx->opt->pip_symmetry)160 par = parallel_constraints(bmap, &first, &second);161 if (par < 0)162 goto error;163 if (!par)164 return SF(basic_map_partial_lexopt_base,SUFFIX)(bmap, dom,165 empty, max);166 167 return SF(basic_map_partial_lexopt_symm,SUFFIX)(bmap, dom, empty, max,168 first, second);169error:170 isl_basic_set_free(dom);171 isl_basic_map_free(bmap);172 return NULL;173}174 175/* Compute the lexicographic minimum (or maximum if "flags" includes176 * ISL_OPT_MAX) of "bmap" over the domain "dom" and return the result as177 * either a map or a piecewise multi-affine expression depending on TYPE.178 * If "empty" is not NULL, then *empty is assigned a set that179 * contains those parts of the domain where there is no solution.180 * If "flags" includes ISL_OPT_FULL, then "dom" is NULL and the optimum181 * should be computed over the domain of "bmap". "empty" is also NULL182 * in this case.183 * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),184 * then we compute the rational optimum. Otherwise, we compute185 * the integral optimum.186 *187 * We perform some preprocessing. As the PILP solver does not188 * handle implicit equalities very well, we first make sure all189 * the equalities are explicitly available.190 *191 * We also add context constraints to the basic map and remove192 * redundant constraints. This is only needed because of the193 * way we handle simple symmetries. In particular, we currently look194 * for symmetries on the constraints, before we set up the main tableau.195 * It is then no good to look for symmetries on possibly redundant constraints.196 * If the domain was extracted from the basic map, then there is197 * no need to add back those constraints again.198 */199__isl_give TYPE *SF(isl_tab_basic_map_partial_lexopt,SUFFIX)(200 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,201 __isl_give isl_set **empty, unsigned flags)202{203 int max, full;204 isl_bool compatible;205 206 if (empty)207 *empty = NULL;208 209 full = ISL_FL_ISSET(flags, ISL_OPT_FULL);210 if (full)211 dom = extract_domain(bmap, flags);212 compatible = isl_basic_map_compatible_domain(bmap, dom);213 if (compatible < 0)214 goto error;215 if (!compatible)216 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,217 "domain does not match input", goto error);218 219 max = ISL_FL_ISSET(flags, ISL_OPT_MAX);220 if (isl_basic_set_dim(dom, isl_dim_all) == 0)221 return SF(basic_map_partial_lexopt,SUFFIX)(bmap, dom, empty,222 max);223 224 if (!full)225 bmap = isl_basic_map_intersect_domain(bmap,226 isl_basic_set_copy(dom));227 bmap = isl_basic_map_detect_equalities(bmap);228 bmap = isl_basic_map_remove_redundancies(bmap);229 230 return SF(basic_map_partial_lexopt,SUFFIX)(bmap, dom, empty, max);231error:232 isl_basic_set_free(dom);233 isl_basic_map_free(bmap);234 return NULL;235}236