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