isl_map_piplib.c

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