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