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