| 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(bset: dom) : NULL; |
| 64 | |
| 65 | bmap_in = isl_basic_map_dim(bmap, type: isl_dim_in); |
| 66 | bmap_param = isl_basic_map_dim(bmap, type: isl_dim_param); |
| 67 | bmap_all = isl_basic_map_dim(bmap, type: 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, size: 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(dst: var->el, src: bmap->ineq[first] + 1 + n_in, len: n_out); |
| 82 | for (i = second + 1, n = 2; i < bmap->n_ineq; ++i) { |
| 83 | if (isl_seq_eq(p1: var->el, p2: bmap->ineq[i] + 1 + n_in, len: n_out) && |
| 84 | all_single_occurrence(bmap, ineq: i, n: n_in)) |
| 85 | list[n++] = i; |
| 86 | } |
| 87 | |
| 88 | cst = isl_mat_alloc(ctx, n_row: n, n_col: 1 + n_in); |
| 89 | if (!cst) |
| 90 | goto error; |
| 91 | |
| 92 | for (i = 0; i < n; ++i) |
| 93 | isl_seq_cpy(dst: cst->row[i], src: bmap->ineq[list[i]], len: 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, pos: list[i]) < 0) |
| 100 | goto error; |
| 101 | |
| 102 | bmap = isl_basic_map_add_dims(bmap, type: isl_dim_in, n: 1); |
| 103 | bmap = isl_basic_map_extend_constraints(base: bmap, n_eq: 0, n_ineq: 1); |
| 104 | k = isl_basic_map_alloc_inequality(bmap); |
| 105 | if (k < 0) |
| 106 | goto error; |
| 107 | isl_seq_clr(p: bmap->ineq[k], len: 1 + n_in); |
| 108 | isl_int_set_si(bmap->ineq[k][1 + n_in], 1); |
| 109 | isl_seq_cpy(dst: bmap->ineq[k] + 1 + n_in + 1, src: var->el, len: n_out); |
| 110 | bmap = isl_basic_map_finalize(bmap); |
| 111 | |
| 112 | n_div = isl_basic_set_dim(bset: dom, type: isl_dim_div); |
| 113 | dom = isl_basic_set_add_dims(bset: dom, type: isl_dim_set, n: 1); |
| 114 | dom = isl_basic_set_extend_constraints(base: dom, n_eq: 0, n_ineq: n); |
| 115 | for (i = 0; i < n; ++i) { |
| 116 | k = isl_basic_set_alloc_inequality(bset: dom); |
| 117 | if (k < 0) |
| 118 | goto error; |
| 119 | isl_seq_cpy(dst: dom->ineq[k], src: cst->row[i], len: 1 + n_in); |
| 120 | isl_int_set_si(dom->ineq[k][1 + n_in], -1); |
| 121 | isl_seq_clr(p: dom->ineq[k] + 1 + n_in + 1, len: n_div); |
| 122 | } |
| 123 | |
| 124 | isl_vec_free(vec: var); |
| 125 | free(ptr: 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(space: map_space); |
| 131 | isl_space_free(space: set_space); |
| 132 | isl_mat_free(mat: cst); |
| 133 | isl_vec_free(vec: var); |
| 134 | free(ptr: list); |
| 135 | isl_basic_set_free(bset: 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 | |
| 156 | if (!bmap) |
| 157 | goto error; |
| 158 | |
| 159 | if (bmap->ctx->opt->pip_symmetry) |
| 160 | par = parallel_constraints(bmap, first: &first, second: &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); |
| 169 | error: |
| 170 | isl_basic_set_free(bset: dom); |
| 171 | isl_basic_map_free(bmap); |
| 172 | return NULL; |
| 173 | } |
| 174 | |
| 175 | /* Compute the lexicographic minimum (or maximum if "flags" includes |
| 176 | * ISL_OPT_MAX) of "bmap" over the domain "dom" and return the result as |
| 177 | * either a map or a piecewise multi-affine expression depending on TYPE. |
| 178 | * If "empty" is not NULL, then *empty is assigned a set that |
| 179 | * contains those parts of the domain where there is no solution. |
| 180 | * If "flags" includes ISL_OPT_FULL, then "dom" is NULL and the optimum |
| 181 | * should be computed over the domain of "bmap". "empty" is also NULL |
| 182 | * in this case. |
| 183 | * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL), |
| 184 | * then we compute the rational optimum. Otherwise, we compute |
| 185 | * the integral optimum. |
| 186 | * |
| 187 | * We perform some preprocessing. As the PILP solver does not |
| 188 | * handle implicit equalities very well, we first make sure all |
| 189 | * the equalities are explicitly available. |
| 190 | * |
| 191 | * We also add context constraints to the basic map and remove |
| 192 | * redundant constraints. This is only needed because of the |
| 193 | * way we handle simple symmetries. In particular, we currently look |
| 194 | * 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 is |
| 197 | * 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, bset: 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(bset: dom, type: 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 | bset: isl_basic_set_copy(bset: 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); |
| 231 | error: |
| 232 | isl_basic_set_free(bset: dom); |
| 233 | isl_basic_map_free(bmap); |
| 234 | return NULL; |
| 235 | } |
| 236 | |