isl_map_piplib.c

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