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->extra, 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         int error;
 173         struct isl_basic_map *bmap = data->bmap;
 174         unsigned old_n_div = bmap->n_div;
 175
 176         if (!map)
 177                 goto error;
 178
 179         for (p = q->newparm; p; p = p->next) {
 180                 int pos;
 181                 unsigned pip_param = bmap->nparam + bmap->n_in;
 182
 183                 pos = find_div(data->ctx, bmap, data->pos, p);
 184                 if (pos < 0)
 185                         goto error;
 186                 data->pos[p->rank - pip_param] = pos;
 187         }
 188
 189         if (q->condition) {
 190                 int j;
 191                 int pos = add_inequality(data->ctx, bmap, data->pos,
 192                                          q->condition);
 193                 if (pos < 0)
 194                         goto error;
 195                 map = scan_quast_r(data, q->next_then, map);
 196
 197                 if (isl_inequality_negate(bmap, pos))
 198                         goto error;
 199                 map = scan_quast_r(data, q->next_else, map);
 200
 201                 if (isl_basic_map_free_inequality(bmap, 1))
 202                         goto error;
 203         } else if (q->list) {
 204                 PipList *l;
 205                 int j;
 206                 /* if bmap->n_out is zero, we are only interested in the domains
 207                  * where a solution exists and not in the actual solution
 208                  */
 209                 for (j = 0, l = q->list; j < bmap->n_out && l; ++j, l = l->next)
 210                         if (add_equality(data->ctx, bmap, data->pos, j,
 211                                                 l->vector) < 0)
 212                                 goto error;
 213                 map = isl_map_add(map, isl_basic_map_copy(bmap));
 214                 if (isl_basic_map_free_equality(bmap, bmap->n_out))
 215                         goto error;
 216         } else if (map->n && data->rest) {
 217                 /* not interested in rest if no sol */
 218                 struct isl_basic_set *bset;
 219                 bset = isl_basic_set_from_basic_map(isl_basic_map_copy(bmap));
 220                 bset = isl_basic_set_drop_vars(bset, bmap->n_in, bmap->n_out);
 221                 if (!bset)
 222                         goto error;
 223                 *data->rest = isl_set_add(*data->rest, bset);
 224         }
 225
 226         if (isl_basic_map_free_inequality(bmap, 2*(bmap->n_div - old_n_div)))
 227                 goto error;
 228         if (isl_basic_map_free_div(bmap, bmap->n_div - old_n_div))
 229                 goto error;
 230         return map;
 231 error:
 232         isl_map_free(map);
 233         return NULL;
 234 }
 235
 236 /*
 237  * Returns a map with "context" as domain and as range the first
 238  * "keep" variables in the quast lists.
 239  */
 240 static struct isl_map *isl_map_from_quast(struct isl_ctx *ctx, PipQuast *q,
 241                 unsigned keep,
 242                 struct isl_basic_set *context,
 243                 struct isl_set **rest)
 244 {
 245         int             pip_param;
 246         int             nexist;
 247         int             max_depth;
 248         int             n_sol, n_nosol;
 249         struct scan_data        data;
 250         struct isl_map          *map;
 251         unsigned                nparam;
 252
 253         data.ctx = ctx;
 254         data.rest = rest;
 255         data.bmap = NULL;
 256         data.pos = NULL;
 257
 258         if (!context)
 259                 goto error;
 260
 261         nparam = context->nparam;
 262         pip_param = nparam + context->dim;
 263
 264         max_depth = 0;
 265         n_sol = 0;
 266         n_nosol = 0;
 267         nexist = pip_param-1;
 268         quast_count(q, &nexist, 0, &max_depth, &n_sol, &n_nosol);
 269         nexist -= pip_param-1;
 270
 271         if (rest) {
 272                 *rest = isl_set_alloc(data.ctx, nparam, context->dim, n_nosol,
 273                                         ISL_MAP_DISJOINT);
 274                 if (!*rest)
 275                         goto error;
 276         }
 277         map = isl_map_alloc(data.ctx, nparam, context->dim, keep, n_sol,
 278                                 ISL_MAP_DISJOINT);
 279         if (!map)
 280                 goto error;
 281
 282         data.bmap = isl_basic_map_from_basic_set(context,
 283                                 context->dim, 0);
 284         data.bmap = isl_basic_map_extend(data.bmap,
 285             nparam, data.bmap->n_in, keep, nexist, keep, max_depth+2*nexist);
 286         if (!data.bmap)
 287                 goto error;
 288
 289         if (data.bmap->extra) {
 290                 int i;
 291                 data.pos = isl_alloc_array(ctx, int, data.bmap->extra);
 292                 if (!data.pos)
 293                         goto error;
 294                 for (i = 0; i < data.bmap->n_div; ++i)
 295                         data.pos[i] = i;
 296         }
 297
 298         map = scan_quast_r(&data, q, map);
 299         map = isl_map_finalize(map);
 300         if (!map)
 301                 goto error;
 302         if (rest) {
 303                 *rest = isl_set_finalize(*rest);
 304                 if (!*rest)
 305                         goto error;
 306         }
 307         isl_basic_map_free(data.bmap);
 308         if (data.pos)
 309                 free(data.pos);
 310         return map;
 311 error:
 312         if (data.pos)
 313                 free(data.pos);
 314         isl_basic_map_free(data.bmap);
 315         isl_map_free(map);
 316         if (rest) {
 317                 isl_set_free(*rest);
 318                 *rest = NULL;
 319         }
 320         return NULL;
 321 }
 322
 323 static void copy_values_to(Entier *dst, isl_int *src, unsigned n)
 324 {
 325         int i;
 326
 327         for (i = 0; i < n; ++i)
 328                 entier_assign(dst[i], src[i]);
 329 }
 330
 331 static void copy_constraint_to(Entier *dst, isl_int *src,
 332                 unsigned pip_param, unsigned pip_var,
 333                 unsigned extra_front, unsigned extra_back)
 334 {
 335         copy_values_to(dst+1+extra_front+pip_var+pip_param+extra_back, src, 1);
 336         copy_values_to(dst+1+extra_front+pip_var, src+1, pip_param);
 337         copy_values_to(dst+1+extra_front, src+1+pip_param, pip_var);
 338 }
 339
 340 PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
 341                          unsigned extra_front, unsigned extra_back)
 342 {
 343         int i;
 344         unsigned nrow;
 345         unsigned ncol;
 346         PipMatrix *M;
 347         unsigned off;
 348         unsigned pip_var = bmap->nparam + bmap->n_in + bmap->n_out
 349                                 + bmap->n_div - pip_param;
 350         unsigned pip_dim = pip_var - bmap->n_div;
 351
 352         nrow = extra_front + bmap->n_eq + bmap->n_ineq;
 353         ncol = 1 + extra_front + pip_var + pip_param + extra_back + 1;
 354         M = pip_matrix_alloc(nrow, ncol);
 355         if (!M)
 356                 return NULL;
 357
 358         off = extra_front;
 359         for (i = 0; i < bmap->n_eq; ++i) {
 360                 entier_set_si(M->p[off+i][0], 0);
 361                 copy_constraint_to(M->p[off+i], bmap->eq[i],
 362                                    pip_param, pip_var, extra_front, extra_back);
 363         }
 364         off += bmap->n_eq;
 365         for (i = 0; i < bmap->n_ineq; ++i) {
 366                 entier_set_si(M->p[off+i][0], 1);
 367                 copy_constraint_to(M->p[off+i], bmap->ineq[i],
 368                                    pip_param, pip_var, extra_front, extra_back);
 369         }
 370         return M;
 371 }
 372
 373 static struct isl_map *extremum_on(
 374                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
 375                 struct isl_set **empty, int max)
 376 {
 377         PipOptions      *options;
 378         PipQuast        *sol;
 379         struct isl_map  *map;
 380         struct isl_ctx  *ctx;
 381         PipMatrix *domain = NULL, *context = NULL;
 382
 383         if (!bmap || !dom)
 384                 goto error;
 385
 386         ctx = bmap->ctx;
 387         isl_assert(ctx, bmap->nparam == dom->nparam, goto error);
 388         isl_assert(ctx, bmap->n_in == dom->dim, goto error);
 389
 390         domain = isl_basic_map_to_pip(bmap, bmap->nparam + bmap->n_in,
 391                                         0, dom->n_div);
 392         if (!domain)
 393                 goto error;
 394         context = isl_basic_map_to_pip((struct isl_basic_map *)dom, 0, 0, 0);
 395         if (!context)
 396                 goto error;
 397
 398         options = pip_options_init();
 399         options->Simplify = 1;
 400         options->Maximize = max;
 401         options->Urs_unknowns = -1;
 402         options->Urs_parms = -1;
 403         sol = pip_solve(domain, context, -1, options);
 404
 405         if (sol) {
 406                 struct isl_basic_set *copy;
 407                 copy = isl_basic_set_copy(dom);
 408                 map = isl_map_from_quast(ctx, sol, bmap->n_out, copy, empty);
 409         } else {
 410                 map = isl_map_empty(ctx, bmap->nparam, bmap->n_in, bmap->n_out);
 411                 if (empty)
 412                         *empty = NULL;
 413         }
 414         if (!map)
 415                 goto error;
 416         if (map->n == 0 && empty) {
 417                 isl_set_free(*empty);
 418                 *empty = isl_set_from_basic_set(dom);
 419         } else
 420                 isl_basic_set_free(dom);
 421         isl_basic_map_free(bmap);
 422
 423         pip_quast_free(sol);
 424         pip_options_free(options);
 425         pip_matrix_free(domain);
 426         pip_matrix_free(context);
 427
 428         return map;
 429 error:
 430         if (domain)
 431                 pip_matrix_free(domain);
 432         if (context)
 433                 pip_matrix_free(context);
 434         isl_basic_map_free(bmap);
 435         isl_basic_set_free(dom);
 436         return NULL;
 437 }
 438
 439 struct isl_map *isl_pip_basic_map_lexmax(
 440                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
 441                 struct isl_set **empty)
 442 {
 443         return extremum_on(bmap, dom, empty, 1);
 444 }
 445
 446 struct isl_map *isl_pip_basic_map_lexmin(
 447                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
 448                 struct isl_set **empty)
 449 {
 450         return extremum_on(bmap, dom, empty, 0);
 451 }
 452
 453 struct isl_map *isl_pip_basic_map_compute_divs(struct isl_basic_map *bmap)
 454 {
 455         PipMatrix *domain = NULL, *context = NULL;
 456         PipOptions      *options;
 457         PipQuast        *sol;
 458         struct isl_ctx  *ctx;
 459         struct isl_map  *map;
 460         struct isl_set  *set;
 461         struct isl_basic_set    *dom;
 462         unsigned         n_in;
 463         unsigned         n_out;
 464
 465         if (!bmap)
 466                 goto error;
 467
 468         ctx = bmap->ctx;
 469         n_in = bmap->n_in;
 470         n_out = bmap->n_out;
 471
 472         domain = isl_basic_map_to_pip(bmap, bmap->nparam + n_in + n_out, 0, 0);
 473         if (!domain)
 474                 goto error;
 475         context = pip_matrix_alloc(0, bmap->nparam + n_in + n_out + 2);
 476         if (!context)
 477                 goto error;
 478
 479         options = pip_options_init();
 480         options->Simplify = 1;
 481         options->Urs_unknowns = -1;
 482         options->Urs_parms = -1;
 483         sol = pip_solve(domain, context, -1, options);
 484
 485         dom = isl_basic_set_alloc(ctx, bmap->nparam, n_in + n_out, 0, 0, 0);
 486         map = isl_map_from_quast(ctx, sol, 0, dom, NULL);
 487
 488         pip_quast_free(sol);
 489         pip_options_free(options);
 490         pip_matrix_free(domain);
 491         pip_matrix_free(context);
 492
 493         isl_basic_map_free(bmap);
 494
 495         set = isl_map_domain(map);
 496
 497         return isl_map_from_set(set, n_in, n_out);
 498 error:
 499         if (domain)
 500                 pip_matrix_free(domain);
 501         if (context)
 502                 pip_matrix_free(context);
 503         isl_basic_set_free(dom);
 504         isl_basic_map_free(bmap);
 505         return NULL;
 506 }