isl_map_piplib.c

   1 #include "isl_set.h"
   2 #include "isl_map.h"
   3 #include "isl_mat.h"
   4 #include "isl_seq.h"
   5 #include "isl_piplib.h"
   6 #include "isl_map_piplib.h"
   7 #include "isl_map_private.h"
   8
   9 static void copy_values_from(isl_int *dst, Entier *src, unsigned n)
  10 {
  11         int i;
  12
  13         for (i = 0; i < n; ++i)
  14                 entier_assign(dst[i], src[i]);
  15 }
  16
  17 static void add_value(isl_int *dst, Entier *src)
  18 {
  19         mpz_add(*dst, *dst, *src);
  20 }
  21
  22 static void copy_constraint_from(isl_int *dst, PipVector *src,
  23                 unsigned nparam, unsigned n_in, unsigned n_out,
  24                 unsigned extra, int *pos)
  25 {
  26         int i;
  27
  28         copy_values_from(dst, src->the_vector+src->nb_elements-1, 1);
  29         copy_values_from(dst+1, src->the_vector, nparam+n_in);
  30         isl_seq_clr(dst+1+nparam+n_in, n_out);
  31         isl_seq_clr(dst+1+nparam+n_in+n_out, extra);
  32         for (i = 0; i + n_in + nparam < src->nb_elements-1; ++i) {
  33                 int p = pos[i];
  34                 add_value(&dst[1+nparam+n_in+n_out+p],
  35                           &src->the_vector[n_in+nparam+i]);
  36         }
  37 }
  38
  39 static int add_inequality(struct isl_ctx *ctx,
  40                    struct isl_basic_map *bmap, int *pos, PipVector *vec)
  41 {
  42         unsigned nparam = isl_basic_map_n_param(bmap);
  43         unsigned n_in = isl_basic_map_n_in(bmap);
  44         unsigned n_out = isl_basic_map_n_out(bmap);
  45         unsigned n_div = isl_basic_map_n_div(bmap);
  46         int i = isl_basic_map_alloc_inequality(bmap);
  47         if (i < 0)
  48                 return -1;
  49         copy_constraint_from(bmap->ineq[i], vec,
  50             nparam, n_in, n_out, n_div, pos);
  51
  52         return i;
  53 }
  54
  55 /* For a div d = floor(f/m), add the constraints
  56  *
  57  *              f - m d >= 0
  58  *              -(f-(n-1)) + m d >= 0
  59  *
  60  * Note that the second constraint is the negation of
  61  *
  62  *              f - m d >= n
  63  */
  64 static int add_div_constraints(struct isl_ctx *ctx,
  65         struct isl_basic_map *bmap, int *pos, PipNewparm *p, unsigned div)
  66 {
  67         int i, j;
  68         unsigned total = isl_basic_map_total_dim(bmap);
  69         unsigned div_pos = 1 + total - bmap->n_div + div;
  70
  71         i = add_inequality(ctx, bmap, pos, p->vector);
  72         if (i < 0)
  73                 return -1;
  74         copy_values_from(&bmap->ineq[i][div_pos], &p->deno, 1);
  75         isl_int_neg(bmap->ineq[i][div_pos], bmap->ineq[i][div_pos]);
  76
  77         j = isl_basic_map_alloc_inequality(bmap);
  78         if (j < 0)
  79                 return -1;
  80         isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
  81         isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
  82         isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
  83         return j;
  84 }
  85
  86 static int add_equality(struct isl_ctx *ctx,
  87                    struct isl_basic_map *bmap, int *pos,
  88                    unsigned var, PipVector *vec)
  89 {
  90         int i;
  91         unsigned nparam = isl_basic_map_n_param(bmap);
  92         unsigned n_in = isl_basic_map_n_in(bmap);
  93         unsigned n_out = isl_basic_map_n_out(bmap);
  94
  95         isl_assert(ctx, var < n_out, return -1);
  96
  97         i = isl_basic_map_alloc_equality(bmap);
  98         if (i < 0)
  99                 return -1;
 100         copy_constraint_from(bmap->eq[i], vec,
 101             nparam, n_in, n_out, bmap->extra, pos);
 102         isl_int_set_si(bmap->eq[i][1+nparam+n_in+var], -1);
 103
 104         return i;
 105 }
 106
 107 static int find_div(struct isl_ctx *ctx,
 108                    struct isl_basic_map *bmap, int *pos, PipNewparm *p)
 109 {
 110         int i, j;
 111         unsigned nparam = isl_basic_map_n_param(bmap);
 112         unsigned n_in = isl_basic_map_n_in(bmap);
 113         unsigned n_out = isl_basic_map_n_out(bmap);
 114
 115         i = isl_basic_map_alloc_div(bmap);
 116         if (i < 0)
 117                 return -1;
 118
 119         copy_constraint_from(bmap->div[i]+1, p->vector,
 120             nparam, n_in, n_out, bmap->extra, pos);
 121
 122         copy_values_from(bmap->div[i], &p->deno, 1);
 123         for (j = 0; j < i; ++j)
 124                 if (isl_seq_eq(bmap->div[i], bmap->div[j],
 125                                 1+1+isl_basic_map_total_dim(bmap)+j)) {
 126                         isl_basic_map_free_div(bmap, 1);
 127                         return j;
 128                 }
 129
 130         if (add_div_constraints(ctx, bmap, pos, p, i) < 0)
 131                 return -1;
 132
 133         return i;
 134 }
 135
 136 /* Count some properties of a quast
 137  * - maximal number of new parameters
 138  * - maximal depth
 139  * - total number of solutions
 140  * - total number of empty branches
 141  */
 142 static void quast_count(PipQuast *q, int *maxnew, int depth, int *maxdepth,
 143                         int *sol, int *nosol)
 144 {
 145         PipNewparm *p;
 146
 147         for (p = q->newparm; p; p = p->next)
 148                 if (p->rank > *maxnew)
 149                         *maxnew = p->rank;
 150         if (q->condition) {
 151                 if (++depth > *maxdepth)
 152                         *maxdepth = depth;
 153                 quast_count(q->next_else, maxnew, depth, maxdepth, sol, nosol);
 154                 quast_count(q->next_then, maxnew, depth, maxdepth, sol, nosol);
 155         } else {
 156                 if (q->list)
 157                         ++(*sol);
 158                 else
 159                         ++(*nosol);
 160         }
 161 }
 162
 163 /*
 164  * pos: array of length bmap->set.extra, mapping each of the existential
 165  *              variables PIP proposes to an existential variable in bmap
 166  * bmap: collects the currently active constraints
 167  * rest: collects the empty leaves of the quast (if not NULL)
 168  */
 169 struct scan_data {
 170         struct isl_ctx                  *ctx;
 171         struct isl_basic_map            *bmap;
 172         struct isl_set                  **rest;
 173         int        *pos;
 174 };
 175
 176 /*
 177  * New existentially quantified variables are places after the existing ones.
 178  */
 179 static struct isl_map *scan_quast_r(struct scan_data *data, PipQuast *q,
 180                                     struct isl_map *map)
 181 {
 182         PipNewparm *p;
 183         struct isl_basic_map *bmap = data->bmap;
 184         unsigned old_n_div = bmap->n_div;
 185         unsigned nparam = isl_basic_map_n_param(bmap);
 186         unsigned n_in = isl_basic_map_n_in(bmap);
 187         unsigned n_out = isl_basic_map_n_out(bmap);
 188
 189         if (!map)
 190                 goto error;
 191
 192         for (p = q->newparm; p; p = p->next) {
 193                 int pos;
 194                 unsigned pip_param = nparam + n_in;
 195
 196                 pos = find_div(data->ctx, bmap, data->pos, p);
 197                 if (pos < 0)
 198                         goto error;
 199                 data->pos[p->rank - pip_param] = pos;
 200         }
 201
 202         if (q->condition) {
 203                 int pos = add_inequality(data->ctx, bmap, data->pos,
 204                                          q->condition);
 205                 if (pos < 0)
 206                         goto error;
 207                 map = scan_quast_r(data, q->next_then, map);
 208
 209                 if (isl_inequality_negate(bmap, pos))
 210                         goto error;
 211                 map = scan_quast_r(data, q->next_else, map);
 212
 213                 if (isl_basic_map_free_inequality(bmap, 1))
 214                         goto error;
 215         } else if (q->list) {
 216                 PipList *l;
 217                 int j;
 218                 /* if bmap->n_out is zero, we are only interested in the domains
 219                  * where a solution exists and not in the actual solution
 220                  */
 221                 for (j = 0, l = q->list; j < n_out && l; ++j, l = l->next)
 222                         if (add_equality(data->ctx, bmap, data->pos, j,
 223                                                 l->vector) < 0)
 224                                 goto error;
 225                 map = isl_map_add(map, isl_basic_map_copy(bmap));
 226                 if (isl_basic_map_free_equality(bmap, n_out))
 227                         goto error;
 228         } else if (data->rest) {
 229                 struct isl_basic_set *bset;
 230                 bset = isl_basic_set_from_basic_map(isl_basic_map_copy(bmap));
 231                 bset = isl_basic_set_drop_dims(bset, n_in, n_out);
 232                 if (!bset)
 233                         goto error;
 234                 *data->rest = isl_set_add(*data->rest, bset);
 235         }
 236
 237         if (isl_basic_map_free_inequality(bmap, 2*(bmap->n_div - old_n_div)))
 238                 goto error;
 239         if (isl_basic_map_free_div(bmap, bmap->n_div - old_n_div))
 240                 goto error;
 241         return map;
 242 error:
 243         isl_map_free(map);
 244         return NULL;
 245 }
 246
 247 /*
 248  * Returns a map of dimension "keep_dim" with "context" as domain and
 249  * as range the first "isl_dim_size(keep_dim, isl_dim_out)" variables
 250  * in the quast lists.
 251  */
 252 static struct isl_map *isl_map_from_quast(struct isl_ctx *ctx, PipQuast *q,
 253                 struct isl_dim *keep_dim,
 254                 struct isl_basic_set *context,
 255                 struct isl_set **rest)
 256 {
 257         int             pip_param;
 258         int             nexist;
 259         int             max_depth;
 260         int             n_sol, n_nosol;
 261         struct scan_data        data;
 262         struct isl_map          *map = NULL;
 263         struct isl_dim          *dims;
 264         unsigned                nparam;
 265         unsigned                dim;
 266         unsigned                keep;
 267
 268         data.ctx = ctx;
 269         data.rest = rest;
 270         data.bmap = NULL;
 271         data.pos = NULL;
 272
 273         if (!context || !keep_dim)
 274                 goto error;
 275
 276         dim = isl_basic_set_n_dim(context);
 277         nparam = isl_basic_set_n_param(context);
 278         keep = isl_dim_size(keep_dim, isl_dim_out);
 279         pip_param = nparam + dim;
 280
 281         max_depth = 0;
 282         n_sol = 0;
 283         n_nosol = 0;
 284         nexist = pip_param-1;
 285         quast_count(q, &nexist, 0, &max_depth, &n_sol, &n_nosol);
 286         nexist -= pip_param-1;
 287
 288         if (rest) {
 289                 *rest = isl_set_alloc_dim(isl_dim_copy(context->dim), n_nosol,
 290                                         ISL_MAP_DISJOINT);
 291                 if (!*rest)
 292                         goto error;
 293         }
 294         map = isl_map_alloc_dim(isl_dim_copy(keep_dim), n_sol,
 295                                 ISL_MAP_DISJOINT);
 296         if (!map)
 297                 goto error;
 298
 299         dims = isl_dim_reverse(isl_dim_copy(context->dim));
 300         data.bmap = isl_basic_map_from_basic_set(context, dims);
 301         data.bmap = isl_basic_map_extend_dim(data.bmap,
 302                 keep_dim, nexist, keep, max_depth+2*nexist);
 303         if (!data.bmap)
 304                 goto error2;
 305
 306         if (data.bmap->extra) {
 307                 int i;
 308                 data.pos = isl_alloc_array(ctx, int, data.bmap->extra);
 309                 if (!data.pos)
 310                         goto error;
 311                 for (i = 0; i < data.bmap->n_div; ++i)
 312                         data.pos[i] = i;
 313         }
 314
 315         map = scan_quast_r(&data, q, map);
 316         map = isl_map_finalize(map);
 317         if (!map)
 318                 goto error2;
 319         if (rest) {
 320                 *rest = isl_set_finalize(*rest);
 321                 if (!*rest)
 322                         goto error2;
 323         }
 324         isl_basic_map_free(data.bmap);
 325         if (data.pos)
 326                 free(data.pos);
 327         return map;
 328 error:
 329         isl_basic_set_free(context);
 330         isl_dim_free(keep_dim);
 331 error2:
 332         if (data.pos)
 333                 free(data.pos);
 334         isl_basic_map_free(data.bmap);
 335         isl_map_free(map);
 336         if (rest) {
 337                 isl_set_free(*rest);
 338                 *rest = NULL;
 339         }
 340         return NULL;
 341 }
 342
 343 static void copy_values_to(Entier *dst, isl_int *src, unsigned n)
 344 {
 345         int i;
 346
 347         for (i = 0; i < n; ++i)
 348                 entier_assign(dst[i], src[i]);
 349 }
 350
 351 static void copy_constraint_to(Entier *dst, isl_int *src,
 352                 unsigned pip_param, unsigned pip_var,
 353                 unsigned extra_front, unsigned extra_back)
 354 {
 355         copy_values_to(dst+1+extra_front+pip_var+pip_param+extra_back, src, 1);
 356         copy_values_to(dst+1+extra_front+pip_var, src+1, pip_param);
 357         copy_values_to(dst+1+extra_front, src+1+pip_param, pip_var);
 358 }
 359
 360 PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
 361                          unsigned extra_front, unsigned extra_back)
 362 {
 363         int i;
 364         unsigned nrow;
 365         unsigned ncol;
 366         PipMatrix *M;
 367         unsigned off;
 368         unsigned pip_var = isl_basic_map_total_dim(bmap) - pip_param;
 369
 370         nrow = extra_front + bmap->n_eq + bmap->n_ineq;
 371         ncol = 1 + extra_front + pip_var + pip_param + extra_back + 1;
 372         M = pip_matrix_alloc(nrow, ncol);
 373         if (!M)
 374                 return NULL;
 375
 376         off = extra_front;
 377         for (i = 0; i < bmap->n_eq; ++i) {
 378                 entier_set_si(M->p[off+i][0], 0);
 379                 copy_constraint_to(M->p[off+i], bmap->eq[i],
 380                                    pip_param, pip_var, extra_front, extra_back);
 381         }
 382         off += bmap->n_eq;
 383         for (i = 0; i < bmap->n_ineq; ++i) {
 384                 entier_set_si(M->p[off+i][0], 1);
 385                 copy_constraint_to(M->p[off+i], bmap->ineq[i],
 386                                    pip_param, pip_var, extra_front, extra_back);
 387         }
 388         return M;
 389 }
 390
 391 PipMatrix *isl_basic_set_to_pip(struct isl_basic_set *bset, unsigned pip_param,
 392                          unsigned extra_front, unsigned extra_back)
 393 {
 394         return isl_basic_map_to_pip((struct isl_basic_map *)bset,
 395                                         pip_param, extra_front, extra_back);
 396 }
 397
 398 struct isl_map *isl_pip_basic_map_lexopt(
 399                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
 400                 struct isl_set **empty, int max)
 401 {
 402         PipOptions      *options;
 403         PipQuast        *sol;
 404         struct isl_map  *map;
 405         struct isl_ctx  *ctx;
 406         PipMatrix *domain = NULL, *context = NULL;
 407         unsigned         nparam, n_in, n_out;
 408
 409         bmap = isl_basic_map_detect_equalities(bmap);
 410         if (!bmap || !dom)
 411                 goto error;
 412
 413         ctx = bmap->ctx;
 414         isl_assert(ctx, isl_basic_map_compatible_domain(bmap, dom), goto error);
 415         nparam = isl_basic_map_n_param(bmap);
 416         n_in = isl_basic_map_n_in(bmap);
 417         n_out = isl_basic_map_n_out(bmap);
 418
 419         domain = isl_basic_map_to_pip(bmap, nparam + n_in, 0, dom->n_div);
 420         if (!domain)
 421                 goto error;
 422         context = isl_basic_map_to_pip((struct isl_basic_map *)dom, 0, 0, 0);
 423         if (!context)
 424                 goto error;
 425
 426         options = pip_options_init();
 427         options->Simplify = 1;
 428         options->Maximize = max;
 429         options->Urs_unknowns = -1;
 430         options->Urs_parms = -1;
 431         sol = pip_solve(domain, context, -1, options);
 432
 433         if (sol) {
 434                 struct isl_basic_set *copy;
 435                 copy = isl_basic_set_copy(dom);
 436                 map = isl_map_from_quast(ctx, sol,
 437                                 isl_dim_copy(bmap->dim), copy, empty);
 438         } else {
 439                 map = isl_map_empty_like_basic_map(bmap);
 440                 if (empty)
 441                         *empty = NULL;
 442         }
 443         if (!map)
 444                 goto error;
 445         if (map->n == 0 && empty) {
 446                 isl_set_free(*empty);
 447                 *empty = isl_set_from_basic_set(dom);
 448         } else
 449                 isl_basic_set_free(dom);
 450         isl_basic_map_free(bmap);
 451
 452         pip_quast_free(sol);
 453         pip_options_free(options);
 454         pip_matrix_free(domain);
 455         pip_matrix_free(context);
 456
 457         return map;
 458 error:
 459         if (domain)
 460                 pip_matrix_free(domain);
 461         if (context)
 462                 pip_matrix_free(context);
 463         isl_basic_map_free(bmap);
 464         isl_basic_set_free(dom);
 465         return NULL;
 466 }