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 #include "isl_equalities.h"
   9
  10 static void copy_values_from(isl_int *dst, Entier *src, unsigned n)
  11 {
  12         int i;
  13
  14         for (i = 0; i < n; ++i)
  15                 entier_assign(dst[i], src[i]);
  16 }
  17
  18 static void add_value(isl_int *dst, Entier *src)
  19 {
  20         mpz_add(*dst, *dst, *src);
  21 }
  22
  23 static void copy_constraint_from(isl_int *dst, PipVector *src,
  24                 unsigned nparam, unsigned n_in, unsigned n_out,
  25                 unsigned extra, int *pos)
  26 {
  27         int i;
  28
  29         copy_values_from(dst, src->the_vector+src->nb_elements-1, 1);
  30         copy_values_from(dst+1, src->the_vector, nparam+n_in);
  31         isl_seq_clr(dst+1+nparam+n_in, n_out);
  32         isl_seq_clr(dst+1+nparam+n_in+n_out, extra);
  33         for (i = 0; i + n_in + nparam < src->nb_elements-1; ++i) {
  34                 int p = pos[i];
  35                 add_value(&dst[1+nparam+n_in+n_out+p],
  36                           &src->the_vector[n_in+nparam+i]);
  37         }
  38 }
  39
  40 static int add_inequality(struct isl_ctx *ctx,
  41                    struct isl_basic_map *bmap, int *pos, PipVector *vec)
  42 {
  43         unsigned nparam = isl_basic_map_n_param(bmap);
  44         unsigned n_in = isl_basic_map_n_in(bmap);
  45         unsigned n_out = isl_basic_map_n_out(bmap);
  46         unsigned n_div = isl_basic_map_n_div(bmap);
  47         int i = isl_basic_map_alloc_inequality(bmap);
  48         if (i < 0)
  49                 return -1;
  50         copy_constraint_from(bmap->ineq[i], vec,
  51             nparam, n_in, n_out, n_div, pos);
  52
  53         return i;
  54 }
  55
  56 /* For a div d = floor(f/m), add the constraints
  57  *
  58  *              f - m d >= 0
  59  *              -(f-(n-1)) + m d >= 0
  60  *
  61  * Note that the second constraint is the negation of
  62  *
  63  *              f - m d >= n
  64  */
  65 static int add_div_constraints(struct isl_ctx *ctx,
  66         struct isl_basic_map *bmap, int *pos, PipNewparm *p, unsigned div)
  67 {
  68         int i, j;
  69         unsigned total = isl_basic_map_total_dim(bmap);
  70         unsigned div_pos = 1 + total - bmap->n_div + div;
  71
  72         i = add_inequality(ctx, bmap, pos, p->vector);
  73         if (i < 0)
  74                 return -1;
  75         copy_values_from(&bmap->ineq[i][div_pos], &p->deno, 1);
  76         isl_int_neg(bmap->ineq[i][div_pos], bmap->ineq[i][div_pos]);
  77
  78         j = isl_basic_map_alloc_inequality(bmap);
  79         if (j < 0)
  80                 return -1;
  81         isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
  82         isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
  83         isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
  84         return j;
  85 }
  86
  87 static int add_equality(struct isl_ctx *ctx,
  88                    struct isl_basic_map *bmap, int *pos,
  89                    unsigned var, PipVector *vec)
  90 {
  91         int i;
  92         unsigned nparam = isl_basic_map_n_param(bmap);
  93         unsigned n_in = isl_basic_map_n_in(bmap);
  94         unsigned n_out = isl_basic_map_n_out(bmap);
  95
  96         isl_assert(ctx, var < n_out, return -1);
  97
  98         i = isl_basic_map_alloc_equality(bmap);
  99         if (i < 0)
 100                 return -1;
 101         copy_constraint_from(bmap->eq[i], vec,
 102             nparam, n_in, n_out, bmap->extra, pos);
 103         isl_int_set_si(bmap->eq[i][1+nparam+n_in+var], -1);
 104
 105         return i;
 106 }
 107
 108 static int find_div(struct isl_ctx *ctx,
 109                    struct isl_basic_map *bmap, int *pos, PipNewparm *p)
 110 {
 111         int i, j;
 112         unsigned nparam = isl_basic_map_n_param(bmap);
 113         unsigned n_in = isl_basic_map_n_in(bmap);
 114         unsigned n_out = isl_basic_map_n_out(bmap);
 115
 116         i = isl_basic_map_alloc_div(bmap);
 117         if (i < 0)
 118                 return -1;
 119
 120         copy_constraint_from(bmap->div[i]+1, p->vector,
 121             nparam, n_in, n_out, bmap->extra, pos);
 122
 123         copy_values_from(bmap->div[i], &p->deno, 1);
 124         for (j = 0; j < i; ++j)
 125                 if (isl_seq_eq(bmap->div[i], bmap->div[j],
 126                                 1+1+isl_basic_map_total_dim(bmap)+j)) {
 127                         isl_basic_map_free_div(bmap, 1);
 128                         return j;
 129                 }
 130
 131         if (add_div_constraints(ctx, bmap, pos, p, i) < 0)
 132                 return -1;
 133
 134         return i;
 135 }
 136
 137 /* Count some properties of a quast
 138  * - maximal number of new parameters
 139  * - maximal depth
 140  * - total number of solutions
 141  * - total number of empty branches
 142  */
 143 static void quast_count(PipQuast *q, int *maxnew, int depth, int *maxdepth,
 144                         int *sol, int *nosol)
 145 {
 146         PipNewparm *p;
 147
 148         for (p = q->newparm; p; p = p->next)
 149                 if (p->rank > *maxnew)
 150                         *maxnew = p->rank;
 151         if (q->condition) {
 152                 if (++depth > *maxdepth)
 153                         *maxdepth = depth;
 154                 quast_count(q->next_else, maxnew, depth, maxdepth, sol, nosol);
 155                 quast_count(q->next_then, maxnew, depth, maxdepth, sol, nosol);
 156         } else {
 157                 if (q->list)
 158                         ++(*sol);
 159                 else
 160                         ++(*nosol);
 161         }
 162 }
 163
 164 /*
 165  * pos: array of length bmap->set.extra, mapping each of the existential
 166  *              variables PIP proposes to an existential variable in bmap
 167  * bmap: collects the currently active constraints
 168  * rest: collects the empty leaves of the quast (if not NULL)
 169  */
 170 struct scan_data {
 171         struct isl_ctx                  *ctx;
 172         struct isl_basic_map            *bmap;
 173         struct isl_set                  **rest;
 174         int        *pos;
 175 };
 176
 177 /*
 178  * New existentially quantified variables are places after the existing ones.
 179  */
 180 static struct isl_map *scan_quast_r(struct scan_data *data, PipQuast *q,
 181                                     struct isl_map *map)
 182 {
 183         PipNewparm *p;
 184         struct isl_basic_map *bmap = data->bmap;
 185         unsigned old_n_div = bmap->n_div;
 186         unsigned nparam = isl_basic_map_n_param(bmap);
 187         unsigned n_in = isl_basic_map_n_in(bmap);
 188         unsigned n_out = isl_basic_map_n_out(bmap);
 189
 190         if (!map)
 191                 goto error;
 192
 193         for (p = q->newparm; p; p = p->next) {
 194                 int pos;
 195                 unsigned pip_param = nparam + n_in;
 196
 197                 pos = find_div(data->ctx, bmap, data->pos, p);
 198                 if (pos < 0)
 199                         goto error;
 200                 data->pos[p->rank - pip_param] = pos;
 201         }
 202
 203         if (q->condition) {
 204                 int pos = add_inequality(data->ctx, bmap, data->pos,
 205                                          q->condition);
 206                 if (pos < 0)
 207                         goto error;
 208                 map = scan_quast_r(data, q->next_then, map);
 209
 210                 if (isl_inequality_negate(bmap, pos))
 211                         goto error;
 212                 map = scan_quast_r(data, q->next_else, map);
 213
 214                 if (isl_basic_map_free_inequality(bmap, 1))
 215                         goto error;
 216         } else if (q->list) {
 217                 PipList *l;
 218                 int j;
 219                 /* if bmap->n_out is zero, we are only interested in the domains
 220                  * where a solution exists and not in the actual solution
 221                  */
 222                 for (j = 0, l = q->list; j < n_out && l; ++j, l = l->next)
 223                         if (add_equality(data->ctx, bmap, data->pos, j,
 224                                                 l->vector) < 0)
 225                                 goto error;
 226                 map = isl_map_add(map, isl_basic_map_copy(bmap));
 227                 if (isl_basic_map_free_equality(bmap, n_out))
 228                         goto error;
 229         } else if (data->rest) {
 230                 struct isl_basic_set *bset;
 231                 bset = isl_basic_set_from_basic_map(isl_basic_map_copy(bmap));
 232                 bset = isl_basic_set_drop_dims(bset, n_in, n_out);
 233                 if (!bset)
 234                         goto error;
 235                 *data->rest = isl_set_add(*data->rest, bset);
 236         }
 237
 238         if (isl_basic_map_free_inequality(bmap, 2*(bmap->n_div - old_n_div)))
 239                 goto error;
 240         if (isl_basic_map_free_div(bmap, bmap->n_div - old_n_div))
 241                 goto error;
 242         return map;
 243 error:
 244         isl_map_free(map);
 245         return NULL;
 246 }
 247
 248 /*
 249  * Returns a map of dimension "keep_dim" with "context" as domain and
 250  * as range the first "isl_dim_size(keep_dim, isl_dim_out)" variables
 251  * in the quast lists.
 252  */
 253 static struct isl_map *isl_map_from_quast(struct isl_ctx *ctx, PipQuast *q,
 254                 struct isl_dim *keep_dim,
 255                 struct isl_basic_set *context,
 256                 struct isl_set **rest)
 257 {
 258         int             pip_param;
 259         int             nexist;
 260         int             max_depth;
 261         int             n_sol, n_nosol;
 262         struct scan_data        data;
 263         struct isl_map          *map = NULL;
 264         struct isl_dim          *dims;
 265         unsigned                nparam;
 266         unsigned                dim;
 267         unsigned                keep;
 268
 269         data.ctx = ctx;
 270         data.rest = rest;
 271         data.bmap = NULL;
 272         data.pos = NULL;
 273
 274         if (!context || !keep_dim)
 275                 goto error;
 276
 277         dim = isl_basic_set_n_dim(context);
 278         nparam = isl_basic_set_n_param(context);
 279         keep = isl_dim_size(keep_dim, isl_dim_out);
 280         pip_param = nparam + dim;
 281
 282         max_depth = 0;
 283         n_sol = 0;
 284         n_nosol = 0;
 285         nexist = pip_param-1;
 286         quast_count(q, &nexist, 0, &max_depth, &n_sol, &n_nosol);
 287         nexist -= pip_param-1;
 288
 289         if (rest) {
 290                 *rest = isl_set_alloc_dim(isl_dim_copy(context->dim), n_nosol,
 291                                         ISL_MAP_DISJOINT);
 292                 if (!*rest)
 293                         goto error;
 294         }
 295         map = isl_map_alloc_dim(isl_dim_copy(keep_dim), n_sol,
 296                                 ISL_MAP_DISJOINT);
 297         if (!map)
 298                 goto error;
 299
 300         dims = isl_dim_reverse(isl_dim_copy(context->dim));
 301         data.bmap = isl_basic_map_from_basic_set(context, dims);
 302         data.bmap = isl_basic_map_extend_dim(data.bmap,
 303                 keep_dim, nexist, keep, max_depth+2*nexist);
 304         if (!data.bmap)
 305                 goto error2;
 306
 307         if (data.bmap->extra) {
 308                 int i;
 309                 data.pos = isl_alloc_array(ctx, int, data.bmap->extra);
 310                 if (!data.pos)
 311                         goto error;
 312                 for (i = 0; i < data.bmap->n_div; ++i)
 313                         data.pos[i] = i;
 314         }
 315
 316         map = scan_quast_r(&data, q, map);
 317         map = isl_map_finalize(map);
 318         if (!map)
 319                 goto error2;
 320         if (rest) {
 321                 *rest = isl_set_finalize(*rest);
 322                 if (!*rest)
 323                         goto error2;
 324         }
 325         isl_basic_map_free(data.bmap);
 326         if (data.pos)
 327                 free(data.pos);
 328         return map;
 329 error:
 330         isl_basic_set_free(context);
 331         isl_dim_free(keep_dim);
 332 error2:
 333         if (data.pos)
 334                 free(data.pos);
 335         isl_basic_map_free(data.bmap);
 336         isl_map_free(map);
 337         if (rest) {
 338                 isl_set_free(*rest);
 339                 *rest = NULL;
 340         }
 341         return NULL;
 342 }
 343
 344 static void copy_values_to(Entier *dst, isl_int *src, unsigned n)
 345 {
 346         int i;
 347
 348         for (i = 0; i < n; ++i)
 349                 entier_assign(dst[i], src[i]);
 350 }
 351
 352 static void copy_constraint_to(Entier *dst, isl_int *src,
 353                 unsigned pip_param, unsigned pip_var,
 354                 unsigned extra_front, unsigned extra_back)
 355 {
 356         copy_values_to(dst+1+extra_front+pip_var+pip_param+extra_back, src, 1);
 357         copy_values_to(dst+1+extra_front+pip_var, src+1, pip_param);
 358         copy_values_to(dst+1+extra_front, src+1+pip_param, pip_var);
 359 }
 360
 361 PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
 362                          unsigned extra_front, unsigned extra_back)
 363 {
 364         int i;
 365         unsigned nrow;
 366         unsigned ncol;
 367         PipMatrix *M;
 368         unsigned off;
 369         unsigned pip_var = isl_basic_map_total_dim(bmap) - pip_param;
 370
 371         nrow = extra_front + bmap->n_eq + bmap->n_ineq;
 372         ncol = 1 + extra_front + pip_var + pip_param + extra_back + 1;
 373         M = pip_matrix_alloc(nrow, ncol);
 374         if (!M)
 375                 return NULL;
 376
 377         off = extra_front;
 378         for (i = 0; i < bmap->n_eq; ++i) {
 379                 entier_set_si(M->p[off+i][0], 0);
 380                 copy_constraint_to(M->p[off+i], bmap->eq[i],
 381                                    pip_param, pip_var, extra_front, extra_back);
 382         }
 383         off += bmap->n_eq;
 384         for (i = 0; i < bmap->n_ineq; ++i) {
 385                 entier_set_si(M->p[off+i][0], 1);
 386                 copy_constraint_to(M->p[off+i], bmap->ineq[i],
 387                                    pip_param, pip_var, extra_front, extra_back);
 388         }
 389         return M;
 390 }
 391
 392 PipMatrix *isl_basic_set_to_pip(struct isl_basic_set *bset, unsigned pip_param,
 393                          unsigned extra_front, unsigned extra_back)
 394 {
 395         return isl_basic_map_to_pip((struct isl_basic_map *)bset,
 396                                         pip_param, extra_front, extra_back);
 397 }
 398
 399 static struct isl_map *extremum_on(
 400                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
 401                 struct isl_set **empty, int max)
 402 {
 403         PipOptions      *options;
 404         PipQuast        *sol;
 405         struct isl_map  *map;
 406         struct isl_ctx  *ctx;
 407         PipMatrix *domain = NULL, *context = NULL;
 408         unsigned         nparam, n_in, n_out;
 409
 410         bmap = isl_basic_map_detect_equalities(bmap);
 411         if (!bmap || !dom)
 412                 goto error;
 413
 414         ctx = bmap->ctx;
 415         isl_assert(ctx, isl_basic_map_compatible_domain(bmap, dom), goto error);
 416         nparam = isl_basic_map_n_param(bmap);
 417         n_in = isl_basic_map_n_in(bmap);
 418         n_out = isl_basic_map_n_out(bmap);
 419
 420         domain = isl_basic_map_to_pip(bmap, nparam + n_in, 0, dom->n_div);
 421         if (!domain)
 422                 goto error;
 423         context = isl_basic_map_to_pip((struct isl_basic_map *)dom, 0, 0, 0);
 424         if (!context)
 425                 goto error;
 426
 427         options = pip_options_init();
 428         options->Simplify = 1;
 429         options->Maximize = max;
 430         options->Urs_unknowns = -1;
 431         options->Urs_parms = -1;
 432         sol = pip_solve(domain, context, -1, options);
 433
 434         if (sol) {
 435                 struct isl_basic_set *copy;
 436                 copy = isl_basic_set_copy(dom);
 437                 map = isl_map_from_quast(ctx, sol,
 438                                 isl_dim_copy(bmap->dim), copy, empty);
 439         } else {
 440                 map = isl_map_empty_like_basic_map(bmap);
 441                 if (empty)
 442                         *empty = NULL;
 443         }
 444         if (!map)
 445                 goto error;
 446         if (map->n == 0 && empty) {
 447                 isl_set_free(*empty);
 448                 *empty = isl_set_from_basic_set(dom);
 449         } else
 450                 isl_basic_set_free(dom);
 451         isl_basic_map_free(bmap);
 452
 453         pip_quast_free(sol);
 454         pip_options_free(options);
 455         pip_matrix_free(domain);
 456         pip_matrix_free(context);
 457
 458         return map;
 459 error:
 460         if (domain)
 461                 pip_matrix_free(domain);
 462         if (context)
 463                 pip_matrix_free(context);
 464         isl_basic_map_free(bmap);
 465         isl_basic_set_free(dom);
 466         return NULL;
 467 }
 468
 469 struct isl_map *isl_pip_basic_map_lexmax(
 470                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
 471                 struct isl_set **empty)
 472 {
 473         return extremum_on(bmap, dom, empty, 1);
 474 }
 475
 476 struct isl_map *isl_pip_basic_map_lexmin(
 477                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
 478                 struct isl_set **empty)
 479 {
 480         return extremum_on(bmap, dom, empty, 0);
 481 }
 482
 483 /* Project the given basic set onto its parameter domain, possibly introducing
 484  * new, explicit, existential variables in the constraints.
 485  * The input has parameters and output variables.
 486  * The output has the same parameters, but no output variables, only
 487  * explicit existentially quantified variables.
 488  */
 489 static struct isl_set *compute_divs_no_eq(struct isl_basic_set *bset)
 490 {
 491         PipMatrix *domain = NULL, *context = NULL;
 492         PipOptions      *options;
 493         PipQuast        *sol;
 494         struct isl_ctx  *ctx;
 495         struct isl_dim  *dim;
 496         struct isl_map  *map;
 497         struct isl_set  *set;
 498         struct isl_basic_set    *dom;
 499         unsigned         nparam;
 500
 501         if (!bset)
 502                 goto error;
 503
 504         ctx = bset->ctx;
 505         nparam = isl_basic_set_n_param(bset);
 506
 507         domain = isl_basic_set_to_pip(bset, nparam, 0, 0);
 508         if (!domain)
 509                 goto error;
 510         context = pip_matrix_alloc(0, nparam + 2);
 511         if (!context)
 512                 goto error;
 513
 514         options = pip_options_init();
 515         options->Simplify = 1;
 516         options->Urs_unknowns = -1;
 517         options->Urs_parms = -1;
 518         sol = pip_solve(domain, context, -1, options);
 519
 520         dom = isl_basic_set_alloc(ctx, nparam, 0, 0, 0, 0);
 521         map = isl_map_from_quast(ctx, sol, isl_dim_copy(dom->dim), dom, NULL);
 522         set = (struct isl_set *)map;
 523
 524         pip_quast_free(sol);
 525         pip_options_free(options);
 526         pip_matrix_free(domain);
 527         pip_matrix_free(context);
 528
 529         isl_basic_set_free(bset);
 530
 531         return set;
 532 error:
 533         if (domain)
 534                 pip_matrix_free(domain);
 535         if (context)
 536                 pip_matrix_free(context);
 537         isl_basic_set_free(dom);
 538         isl_basic_set_free(bset);
 539         return NULL;
 540 }
 541
 542 static struct isl_map *isl_map_reset_dim(struct isl_map *map,
 543         struct isl_dim *dim)
 544 {
 545         int i;
 546
 547         if (!map || !dim)
 548                 goto error;
 549
 550         for (i = 0; i < map->n; ++i) {
 551                 isl_dim_free(map->p[i]->dim);
 552                 map->p[i]->dim = isl_dim_copy(dim);
 553         }
 554         isl_dim_free(map->dim);
 555         map->dim = dim;
 556
 557         return map;
 558 error:
 559         isl_map_free(map);
 560         isl_dim_free(dim);
 561         return NULL;
 562 }
 563
 564 static struct isl_set *isl_set_reset_dim(struct isl_set *set,
 565         struct isl_dim *dim)
 566 {
 567         return (struct isl_set *) isl_map_reset_dim((struct isl_map *)set, dim);
 568 }
 569
 570 /* Given a matrix M (mat) and a size n (size), replace mat
 571  * by the matrix
 572  *
 573  *              [ M 0 ]
 574  *              [ 0 I ]
 575  *
 576  * where I is an n x n identity matrix.
 577  */
 578 static struct isl_mat *append_identity(struct isl_ctx *ctx,
 579         struct isl_mat *mat, unsigned size)
 580 {
 581         int i;
 582         unsigned n_row, n_col;
 583
 584         n_row = mat->n_row;
 585         n_col = mat->n_col;
 586         mat = isl_mat_extend(ctx, mat, n_row + size, n_col + size);
 587         if (!mat)
 588                 return NULL;
 589         for (i = 0; i < n_row; ++i)
 590                 isl_seq_clr(mat->row[i] + n_col, size);
 591         for (i = 0; i < size; ++i) {
 592                 isl_seq_clr(mat->row[n_row + i], n_col + size);
 593                 isl_int_set_si(mat->row[n_row + i][n_col + i], 1);
 594         }
 595         return mat;
 596 }
 597
 598 /* Apply a preimage specified by "mat" on the parameters of "bset".
 599  */
 600 static struct isl_basic_set *basic_set_parameter_preimage(
 601         struct isl_basic_set *bset, struct isl_mat *mat)
 602 {
 603         unsigned nparam, n_out;
 604
 605         if (!bset || !mat)
 606                 goto error;
 607
 608         bset->dim = isl_dim_cow(bset->dim);
 609         if (!bset->dim)
 610                 goto error;
 611
 612         nparam = isl_basic_set_dim(bset, isl_dim_param);
 613         n_out = isl_basic_set_dim(bset, isl_dim_set);
 614
 615         isl_assert(bset->ctx, mat->n_row == 1 + nparam, goto error);
 616
 617         mat = append_identity(bset->ctx, mat, n_out);
 618         if (!mat)
 619                 goto error;
 620
 621         bset->dim->nparam = 0;
 622         bset->dim->n_out += nparam;
 623         bset = isl_basic_set_preimage(bset, mat);
 624         if (bset) {
 625                 bset->dim->nparam = bset->dim->n_out - n_out;
 626                 bset->dim->n_out = n_out;
 627         }
 628         return bset;
 629 error:
 630         isl_mat_free(bset ? bset->ctx : NULL, mat);
 631         isl_basic_set_free(bset);
 632         return NULL;
 633 }
 634
 635 /* Apply a preimage specified by "mat" on the parameters of "set".
 636  */
 637 static struct isl_set *set_parameter_preimage(
 638         struct isl_set *set, struct isl_mat *mat)
 639 {
 640         struct isl_dim *dim = NULL;
 641         unsigned nparam, n_out;
 642
 643         if (!set || !mat)
 644                 goto error;
 645
 646         dim = isl_dim_copy(set->dim);
 647         dim = isl_dim_cow(dim);
 648         if (!dim)
 649                 goto error;
 650
 651         nparam = isl_set_dim(set, isl_dim_param);
 652         n_out = isl_set_dim(set, isl_dim_set);
 653
 654         isl_assert(set->ctx, mat->n_row == 1 + nparam, goto error);
 655
 656         mat = append_identity(set->ctx, mat, n_out);
 657         if (!mat)
 658                 goto error;
 659
 660         dim->nparam = 0;
 661         dim->n_out += nparam;
 662         isl_set_reset_dim(set, dim);
 663         set = isl_set_preimage(set, mat);
 664         if (!set)
 665                 goto error2;
 666         dim = isl_dim_copy(set->dim);
 667         dim = isl_dim_cow(dim);
 668         if (!dim)
 669                 goto error2;
 670         dim->nparam = dim->n_out - n_out;
 671         dim->n_out = n_out;
 672         isl_set_reset_dim(set, dim);
 673         return set;
 674 error:
 675         isl_dim_free(dim);
 676         isl_mat_free(set ? set->ctx : NULL, mat);
 677 error2:
 678         isl_set_free(set);
 679         return NULL;
 680 }
 681
 682 /* Intersect the basic set "bset" with the affine space specified by the
 683  * equalities in "eq".
 684  */
 685 static struct isl_basic_set *basic_set_append_equalities(
 686         struct isl_basic_set *bset, struct isl_mat *eq)
 687 {
 688         int i, k;
 689         unsigned len;
 690
 691         if (!bset || !eq)
 692                 goto error;
 693
 694         bset = isl_basic_set_extend_dim(bset, isl_dim_copy(bset->dim), 0,
 695                                         eq->n_row, 0);
 696         if (!bset)
 697                 goto error;
 698
 699         len = 1 + isl_dim_total(bset->dim) + bset->extra;
 700         for (i = 0; i < eq->n_row; ++i) {
 701                 k = isl_basic_set_alloc_equality(bset);
 702                 if (k < 0)
 703                         goto error;
 704                 isl_seq_cpy(bset->eq[k], eq->row[i], eq->n_col);
 705                 isl_seq_clr(bset->eq[k] + eq->n_col, len - eq->n_col);
 706         }
 707         isl_mat_free(bset->ctx, eq);
 708
 709         return bset;
 710 error:
 711         isl_mat_free(bset ? bset->ctx : NULL, eq);
 712         isl_basic_set_free(bset);
 713         return NULL;
 714 }
 715
 716 /* Intersect the set "set" with the affine space specified by the
 717  * equalities in "eq".
 718  */
 719 static struct isl_set *set_append_equalities(struct isl_set *set,
 720         struct isl_mat *eq)
 721 {
 722         int i;
 723
 724         if (!set || !eq)
 725                 goto error;
 726
 727         for (i = 0; i < set->n; ++i) {
 728                 set->p[i] = basic_set_append_equalities(set->p[i],
 729                                         isl_mat_copy(set->ctx, eq));
 730                 if (!set->p[i])
 731                         goto error;
 732         }
 733         isl_mat_free(set->ctx, eq);
 734         return set;
 735 error:
 736         isl_mat_free(set ? set->ctx : NULL, eq);
 737         isl_set_free(set);
 738         return NULL;
 739 }
 740
 741 /* Project the given basic set onto its parameter domain, possibly introducing
 742  * new, explicit, existential variables in the constraints.
 743  * The input has parameters and output variables.
 744  * The output has the same parameters, but no output variables, only
 745  * explicit existentially quantified variables.
 746  *
 747  * The actual projection is performed by pip, but pip doesn't seem
 748  * to like equalities very much, so we first remove the equalities
 749  * among the parameters by performing a variable compression on
 750  * the parameters.  Afterward, an inverse transformation is performed
 751  * and the equalities among the parameters are inserted back in.
 752  */
 753 static struct isl_set *compute_divs(struct isl_basic_set *bset)
 754 {
 755         int i, j;
 756         struct isl_mat *eq;
 757         struct isl_mat *T, *T2;
 758         struct isl_set *set;
 759         unsigned nparam, n_out;
 760
 761         bset = isl_basic_set_cow(bset);
 762         if (!bset)
 763                 return NULL;
 764
 765         if (bset->n_eq == 0)
 766                 return compute_divs_no_eq(bset);
 767
 768         isl_basic_set_gauss(bset, NULL);
 769
 770         nparam = isl_basic_set_dim(bset, isl_dim_param);
 771         n_out = isl_basic_set_dim(bset, isl_dim_out);
 772
 773         for (i = 0, j = n_out - 1; i < bset->n_eq && j >= 0; --j) {
 774                 if (!isl_int_is_zero(bset->eq[i][1 + nparam + j]))
 775                         ++i;
 776         }
 777         if (i == bset->n_eq)
 778                 return compute_divs_no_eq(bset);
 779
 780         eq = isl_mat_sub_alloc(bset->ctx, bset->eq, i, bset->n_eq - i,
 781                 0, 1 + nparam);
 782         eq = isl_mat_cow(bset->ctx, eq);
 783         T = isl_mat_variable_compression(bset->ctx,
 784                 isl_mat_copy(bset->ctx, eq), &T2);
 785         bset = basic_set_parameter_preimage(bset, T);
 786
 787         set = compute_divs_no_eq(bset);
 788         set = set_parameter_preimage(set, T2);
 789         set = set_append_equalities(set, eq);
 790         return set;
 791 }
 792
 793 /* Compute an explicit representation for all the existentially
 794  * quantified variables.
 795  * The input and output dimensions are first turned into parameters
 796  * and the existential variables into output dimensions.
 797  * compute_divs then returns a map with the same parameters and
 798  * no input or output dimensions and the dimension specification
 799  * is reset to that of the input.
 800  */
 801 struct isl_map *isl_pip_basic_map_compute_divs(struct isl_basic_map *bmap)
 802 {
 803         struct isl_basic_set *bset;
 804         struct isl_set *set;
 805         struct isl_map *map;
 806         struct isl_dim *dim, *orig_dim = NULL;
 807         unsigned         nparam;
 808         unsigned         n_in;
 809         unsigned         n_out;
 810
 811         bmap = isl_basic_map_cow(bmap);
 812         if (!bmap)
 813                 return NULL;
 814
 815         nparam = isl_basic_map_dim(bmap, isl_dim_param);
 816         n_in = isl_basic_map_dim(bmap, isl_dim_in);
 817         n_out = isl_basic_map_dim(bmap, isl_dim_out);
 818         dim = isl_dim_set_alloc(bmap->ctx, nparam + n_in + n_out, bmap->n_div);
 819         if (!dim)
 820                 goto error;
 821
 822         orig_dim = bmap->dim;
 823         bmap->dim = dim;
 824         bmap->extra -= bmap->n_div;
 825         bmap->n_div = 0;
 826         bset = (struct isl_basic_set *)bmap;
 827
 828         set = compute_divs(bset);
 829         map = (struct isl_map *)set;
 830         map = isl_map_reset_dim(map, orig_dim);
 831
 832         return map;
 833 error:
 834         isl_basic_map_free(bmap);
 835         return NULL;
 836 }