isl_map_piplib.c

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