isl_map_piplib.c

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