isl_morph.c

   1 /*
   2  * Copyright 2010      INRIA Saclay
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
   7  * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
   8  * 91893 Orsay, France
   9  */
  10
  11 #include <isl_map_private.h>
  12 #include <isl_morph.h>
  13 #include <isl/seq.h>
  14 #include <isl_mat_private.h>
  15 #include <isl_space_private.h>
  16 #include <isl_equalities.h>
  17
  18 __isl_give isl_morph *isl_morph_alloc(
  19         __isl_take isl_basic_set *dom, __isl_take isl_basic_set *ran,
  20         __isl_take isl_mat *map, __isl_take isl_mat *inv)
  21 {
  22         isl_morph *morph;
  23
  24         if (!dom || !ran || !map || !inv)
  25                 goto error;
  26
  27         morph = isl_alloc_type(dom->ctx, struct isl_morph);
  28         if (!morph)
  29                 goto error;
  30
  31         morph->ref = 1;
  32         morph->dom = dom;
  33         morph->ran = ran;
  34         morph->map = map;
  35         morph->inv = inv;
  36
  37         return morph;
  38 error:
  39         isl_basic_set_free(dom);
  40         isl_basic_set_free(ran);
  41         isl_mat_free(map);
  42         isl_mat_free(inv);
  43         return NULL;
  44 }
  45
  46 __isl_give isl_morph *isl_morph_copy(__isl_keep isl_morph *morph)
  47 {
  48         if (!morph)
  49                 return NULL;
  50
  51         morph->ref++;
  52         return morph;
  53 }
  54
  55 __isl_give isl_morph *isl_morph_dup(__isl_keep isl_morph *morph)
  56 {
  57         if (!morph)
  58                 return NULL;
  59
  60         return isl_morph_alloc(isl_basic_set_copy(morph->dom),
  61                 isl_basic_set_copy(morph->ran),
  62                 isl_mat_copy(morph->map), isl_mat_copy(morph->inv));
  63 }
  64
  65 __isl_give isl_morph *isl_morph_cow(__isl_take isl_morph *morph)
  66 {
  67         if (!morph)
  68                 return NULL;
  69
  70         if (morph->ref == 1)
  71                 return morph;
  72         morph->ref--;
  73         return isl_morph_dup(morph);
  74 }
  75
  76 void isl_morph_free(__isl_take isl_morph *morph)
  77 {
  78         if (!morph)
  79                 return;
  80
  81         if (--morph->ref > 0)
  82                 return;
  83
  84         isl_basic_set_free(morph->dom);
  85         isl_basic_set_free(morph->ran);
  86         isl_mat_free(morph->map);
  87         isl_mat_free(morph->inv);
  88         free(morph);
  89 }
  90
  91 __isl_give isl_space *isl_morph_get_ran_space(__isl_keep isl_morph *morph)
  92 {
  93         if (!morph)
  94                 return NULL;
  95
  96         return isl_space_copy(morph->ran->dim);
  97 }
  98
  99 unsigned isl_morph_dom_dim(__isl_keep isl_morph *morph, enum isl_dim_type type)
 100 {
 101         if (!morph)
 102                 return 0;
 103
 104         return isl_basic_set_dim(morph->dom, type);
 105 }
 106
 107 unsigned isl_morph_ran_dim(__isl_keep isl_morph *morph, enum isl_dim_type type)
 108 {
 109         if (!morph)
 110                 return 0;
 111
 112         return isl_basic_set_dim(morph->ran, type);
 113 }
 114
 115 __isl_give isl_morph *isl_morph_remove_dom_dims(__isl_take isl_morph *morph,
 116         enum isl_dim_type type, unsigned first, unsigned n)
 117 {
 118         unsigned dom_offset;
 119
 120         if (n == 0)
 121                 return morph;
 122
 123         morph = isl_morph_cow(morph);
 124         if (!morph)
 125                 return NULL;
 126
 127         dom_offset = 1 + isl_space_offset(morph->dom->dim, type);
 128
 129         morph->dom = isl_basic_set_remove_dims(morph->dom, type, first, n);
 130
 131         morph->map = isl_mat_drop_cols(morph->map, dom_offset + first, n);
 132
 133         morph->inv = isl_mat_drop_rows(morph->inv, dom_offset + first, n);
 134
 135         if (morph->dom && morph->ran && morph->map && morph->inv)
 136                 return morph;
 137
 138         isl_morph_free(morph);
 139         return NULL;
 140 }
 141
 142 __isl_give isl_morph *isl_morph_remove_ran_dims(__isl_take isl_morph *morph,
 143         enum isl_dim_type type, unsigned first, unsigned n)
 144 {
 145         unsigned ran_offset;
 146
 147         if (n == 0)
 148                 return morph;
 149
 150         morph = isl_morph_cow(morph);
 151         if (!morph)
 152                 return NULL;
 153
 154         ran_offset = 1 + isl_space_offset(morph->ran->dim, type);
 155
 156         morph->ran = isl_basic_set_remove_dims(morph->ran, type, first, n);
 157
 158         morph->map = isl_mat_drop_rows(morph->map, ran_offset + first, n);
 159
 160         morph->inv = isl_mat_drop_cols(morph->inv, ran_offset + first, n);
 161
 162         if (morph->dom && morph->ran && morph->map && morph->inv)
 163                 return morph;
 164
 165         isl_morph_free(morph);
 166         return NULL;
 167 }
 168
 169 /* Project domain of morph onto its parameter domain.
 170  */
 171 __isl_give isl_morph *isl_morph_dom_params(__isl_take isl_morph *morph)
 172 {
 173         unsigned n;
 174
 175         if (!morph)
 176                 return NULL;
 177         n = isl_basic_set_dim(morph->dom, isl_dim_set);
 178         morph = isl_morph_remove_dom_dims(morph, isl_dim_set, 0, n);
 179         if (!morph)
 180                 return NULL;
 181         morph->dom = isl_basic_set_params(morph->dom);
 182         if (morph->dom)
 183                 return morph;
 184
 185         isl_morph_free(morph);
 186         return NULL;
 187 }
 188
 189 /* Project range of morph onto its parameter domain.
 190  */
 191 __isl_give isl_morph *isl_morph_ran_params(__isl_take isl_morph *morph)
 192 {
 193         unsigned n;
 194
 195         if (!morph)
 196                 return NULL;
 197         n = isl_basic_set_dim(morph->ran, isl_dim_set);
 198         morph = isl_morph_remove_ran_dims(morph, isl_dim_set, 0, n);
 199         if (!morph)
 200                 return NULL;
 201         morph->ran = isl_basic_set_params(morph->ran);
 202         if (morph->ran)
 203                 return morph;
 204
 205         isl_morph_free(morph);
 206         return NULL;
 207 }
 208
 209 void isl_morph_dump(__isl_take isl_morph *morph, FILE *out)
 210 {
 211         if (!morph)
 212                 return;
 213
 214         isl_basic_set_print(morph->dom, out, 0, "", "", ISL_FORMAT_ISL);
 215         isl_basic_set_print(morph->ran, out, 0, "", "", ISL_FORMAT_ISL);
 216         isl_mat_print_internal(morph->map, out, 4);
 217         isl_mat_print_internal(morph->inv, out, 4);
 218 }
 219
 220 __isl_give isl_morph *isl_morph_identity(__isl_keep isl_basic_set *bset)
 221 {
 222         isl_mat *id;
 223         isl_basic_set *universe;
 224         unsigned total;
 225
 226         if (!bset)
 227                 return NULL;
 228
 229         total = isl_basic_set_total_dim(bset);
 230         id = isl_mat_identity(bset->ctx, 1 + total);
 231         universe = isl_basic_set_universe(isl_space_copy(bset->dim));
 232
 233         return isl_morph_alloc(universe, isl_basic_set_copy(universe),
 234                 id, isl_mat_copy(id));
 235 }
 236
 237 /* Create a(n identity) morphism between empty sets of the same dimension
 238  * a "bset".
 239  */
 240 __isl_give isl_morph *isl_morph_empty(__isl_keep isl_basic_set *bset)
 241 {
 242         isl_mat *id;
 243         isl_basic_set *empty;
 244         unsigned total;
 245
 246         if (!bset)
 247                 return NULL;
 248
 249         total = isl_basic_set_total_dim(bset);
 250         id = isl_mat_identity(bset->ctx, 1 + total);
 251         empty = isl_basic_set_empty(isl_space_copy(bset->dim));
 252
 253         return isl_morph_alloc(empty, isl_basic_set_copy(empty),
 254                 id, isl_mat_copy(id));
 255 }
 256
 257 /* Given a matrix that maps a (possibly) parametric domain to
 258  * a parametric domain, add in rows that map the "nparam" parameters onto
 259  * themselves.
 260  */
 261 static __isl_give isl_mat *insert_parameter_rows(__isl_take isl_mat *mat,
 262         unsigned nparam)
 263 {
 264         int i;
 265
 266         if (nparam == 0)
 267                 return mat;
 268         if (!mat)
 269                 return NULL;
 270
 271         mat = isl_mat_insert_rows(mat, 1, nparam);
 272         if (!mat)
 273                 return NULL;
 274
 275         for (i = 0; i < nparam; ++i) {
 276                 isl_seq_clr(mat->row[1 + i], mat->n_col);
 277                 isl_int_set(mat->row[1 + i][1 + i], mat->row[0][0]);
 278         }
 279
 280         return mat;
 281 }
 282
 283 /* Construct a basic set described by the "n" equalities of "bset" starting
 284  * at "first".
 285  */
 286 static __isl_give isl_basic_set *copy_equalities(__isl_keep isl_basic_set *bset,
 287         unsigned first, unsigned n)
 288 {
 289         int i, k;
 290         isl_basic_set *eq;
 291         unsigned total;
 292
 293         isl_assert(bset->ctx, bset->n_div == 0, return NULL);
 294
 295         total = isl_basic_set_total_dim(bset);
 296         eq = isl_basic_set_alloc_space(isl_space_copy(bset->dim), 0, n, 0);
 297         if (!eq)
 298                 return NULL;
 299         for (i = 0; i < n; ++i) {
 300                 k = isl_basic_set_alloc_equality(eq);
 301                 if (k < 0)
 302                         goto error;
 303                 isl_seq_cpy(eq->eq[k], bset->eq[first + k], 1 + total);
 304         }
 305
 306         return eq;
 307 error:
 308         isl_basic_set_free(eq);
 309         return NULL;
 310 }
 311
 312 /* Given a basic set, exploit the equalties in the a basic set to construct
 313  * a morphishm that maps the basic set to a lower-dimensional space.
 314  * Specifically, the morphism reduces the number of dimensions of type "type".
 315  *
 316  * This function is a slight generalization of isl_mat_variable_compression
 317  * in that it allows the input to be parametric and that it allows for the
 318  * compression of either parameters or set variables.
 319  *
 320  * We first select the equalities of interest, that is those that involve
 321  * variables of type "type" and no later variables.
 322  * Denote those equalities as
 323  *
 324  *              -C(p) + M x = 0
 325  *
 326  * where C(p) depends on the parameters if type == isl_dim_set and
 327  * is a constant if type == isl_dim_param.
 328  *
 329  * First compute the (left) Hermite normal form of M,
 330  *
 331  *              M [U1 U2] = M U = H = [H1 0]
 332  * or
 333  *                            M = H Q = [H1 0] [Q1]
 334  *                                             [Q2]
 335  *
 336  * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
 337  * Define the transformed variables as
 338  *
 339  *              x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
 340  *                          [ x2' ]           [Q2]
 341  *
 342  * The equalities then become
 343  *
 344  *              -C(p) + H1 x1' = 0   or   x1' = H1^{-1} C(p) = C'(p)
 345  *
 346  * If the denominator of the constant term does not divide the
 347  * the common denominator of the parametric terms, then every
 348  * integer point is mapped to a non-integer point and then the original set has no
 349  * integer solutions (since the x' are a unimodular transformation
 350  * of the x).  In this case, an empty morphism is returned.
 351  * Otherwise, the transformation is given by
 352  *
 353  *              x = U1 H1^{-1} C(p) + U2 x2'
 354  *
 355  * The inverse transformation is simply
 356  *
 357  *              x2' = Q2 x
 358  *
 359  * Both matrices are extended to map the full original space to the full
 360  * compressed space.
 361  */
 362 __isl_give isl_morph *isl_basic_set_variable_compression(
 363         __isl_keep isl_basic_set *bset, enum isl_dim_type type)
 364 {
 365         unsigned otype;
 366         unsigned ntype;
 367         unsigned orest;
 368         unsigned nrest;
 369         int f_eq, n_eq;
 370         isl_space *dim;
 371         isl_mat *H, *U, *Q, *C = NULL, *H1, *U1, *U2;
 372         isl_basic_set *dom, *ran;
 373
 374         if (!bset)
 375                 return NULL;
 376
 377         if (isl_basic_set_plain_is_empty(bset))
 378                 return isl_morph_empty(bset);
 379
 380         isl_assert(bset->ctx, bset->n_div == 0, return NULL);
 381
 382         otype = 1 + isl_space_offset(bset->dim, type);
 383         ntype = isl_basic_set_dim(bset, type);
 384         orest = otype + ntype;
 385         nrest = isl_basic_set_total_dim(bset) - (orest - 1);
 386
 387         for (f_eq = 0; f_eq < bset->n_eq; ++f_eq)
 388                 if (isl_seq_first_non_zero(bset->eq[f_eq] + orest, nrest) == -1)
 389                         break;
 390         for (n_eq = 0; f_eq + n_eq < bset->n_eq; ++n_eq)
 391                 if (isl_seq_first_non_zero(bset->eq[f_eq + n_eq] + otype, ntype) == -1)
 392                         break;
 393         if (n_eq == 0)
 394                 return isl_morph_identity(bset);
 395
 396         H = isl_mat_sub_alloc6(bset->ctx, bset->eq, f_eq, n_eq, otype, ntype);
 397         H = isl_mat_left_hermite(H, 0, &U, &Q);
 398         if (!H || !U || !Q)
 399                 goto error;
 400         Q = isl_mat_drop_rows(Q, 0, n_eq);
 401         Q = isl_mat_diagonal(isl_mat_identity(bset->ctx, otype), Q);
 402         Q = isl_mat_diagonal(Q, isl_mat_identity(bset->ctx, nrest));
 403         C = isl_mat_alloc(bset->ctx, 1 + n_eq, otype);
 404         if (!C)
 405                 goto error;
 406         isl_int_set_si(C->row[0][0], 1);
 407         isl_seq_clr(C->row[0] + 1, otype - 1);
 408         isl_mat_sub_neg(C->ctx, C->row + 1, bset->eq + f_eq, n_eq, 0, 0, otype);
 409         H1 = isl_mat_sub_alloc(H, 0, H->n_row, 0, H->n_row);
 410         H1 = isl_mat_lin_to_aff(H1);
 411         C = isl_mat_inverse_product(H1, C);
 412         if (!C)
 413                 goto error;
 414         isl_mat_free(H);
 415
 416         if (!isl_int_is_one(C->row[0][0])) {
 417                 int i;
 418                 isl_int g;
 419
 420                 isl_int_init(g);
 421                 for (i = 0; i < n_eq; ++i) {
 422                         isl_seq_gcd(C->row[1 + i] + 1, otype - 1, &g);
 423                         isl_int_gcd(g, g, C->row[0][0]);
 424                         if (!isl_int_is_divisible_by(C->row[1 + i][0], g))
 425                                 break;
 426                 }
 427                 isl_int_clear(g);
 428
 429                 if (i < n_eq) {
 430                         isl_mat_free(C);
 431                         isl_mat_free(U);
 432                         isl_mat_free(Q);
 433                         return isl_morph_empty(bset);
 434                 }
 435
 436                 C = isl_mat_normalize(C);
 437         }
 438
 439         U1 = isl_mat_sub_alloc(U, 0, U->n_row, 0, n_eq);
 440         U1 = isl_mat_lin_to_aff(U1);
 441         U2 = isl_mat_sub_alloc(U, 0, U->n_row, n_eq, U->n_row - n_eq);
 442         U2 = isl_mat_lin_to_aff(U2);
 443         isl_mat_free(U);
 444
 445         C = isl_mat_product(U1, C);
 446         C = isl_mat_aff_direct_sum(C, U2);
 447         C = insert_parameter_rows(C, otype - 1);
 448         C = isl_mat_diagonal(C, isl_mat_identity(bset->ctx, nrest));
 449
 450         dim = isl_space_copy(bset->dim);
 451         dim = isl_space_drop_dims(dim, type, 0, ntype);
 452         dim = isl_space_add_dims(dim, type, ntype - n_eq);
 453         ran = isl_basic_set_universe(dim);
 454         dom = copy_equalities(bset, f_eq, n_eq);
 455
 456         return isl_morph_alloc(dom, ran, Q, C);
 457 error:
 458         isl_mat_free(C);
 459         isl_mat_free(H);
 460         isl_mat_free(U);
 461         isl_mat_free(Q);
 462         return NULL;
 463 }
 464
 465 /* Construct a parameter compression for "bset".
 466  * We basically just call isl_mat_parameter_compression with the right input
 467  * and then extend the resulting matrix to include the variables.
 468  *
 469  * Let the equalities be given as
 470  *
 471  *      B(p) + A x = 0
 472  *
 473  * and let [H 0] be the Hermite Normal Form of A, then
 474  *
 475  *      H^-1 B(p)
 476  *
 477  * needs to be integer, so we impose that each row is divisible by
 478  * the denominator.
 479  */
 480 __isl_give isl_morph *isl_basic_set_parameter_compression(
 481         __isl_keep isl_basic_set *bset)
 482 {
 483         unsigned nparam;
 484         unsigned nvar;
 485         int n_eq;
 486         isl_mat *H, *B;
 487         isl_vec *d;
 488         isl_mat *map, *inv;
 489         isl_basic_set *dom, *ran;
 490
 491         if (!bset)
 492                 return NULL;
 493
 494         if (isl_basic_set_plain_is_empty(bset))
 495                 return isl_morph_empty(bset);
 496         if (bset->n_eq == 0)
 497                 return isl_morph_identity(bset);
 498
 499         isl_assert(bset->ctx, bset->n_div == 0, return NULL);
 500
 501         n_eq = bset->n_eq;
 502         nparam = isl_basic_set_dim(bset, isl_dim_param);
 503         nvar = isl_basic_set_dim(bset, isl_dim_set);
 504
 505         isl_assert(bset->ctx, n_eq <= nvar, return NULL);
 506
 507         d = isl_vec_alloc(bset->ctx, n_eq);
 508         B = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, n_eq, 0, 1 + nparam);
 509         H = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, n_eq, 1 + nparam, nvar);
 510         H = isl_mat_left_hermite(H, 0, NULL, NULL);
 511         H = isl_mat_drop_cols(H, n_eq, nvar - n_eq);
 512         H = isl_mat_lin_to_aff(H);
 513         H = isl_mat_right_inverse(H);
 514         if (!H || !d)
 515                 goto error;
 516         isl_seq_set(d->el, H->row[0][0], d->size);
 517         H = isl_mat_drop_rows(H, 0, 1);
 518         H = isl_mat_drop_cols(H, 0, 1);
 519         B = isl_mat_product(H, B);
 520         inv = isl_mat_parameter_compression(B, d);
 521         inv = isl_mat_diagonal(inv, isl_mat_identity(bset->ctx, nvar));
 522         map = isl_mat_right_inverse(isl_mat_copy(inv));
 523
 524         dom = isl_basic_set_universe(isl_space_copy(bset->dim));
 525         ran = isl_basic_set_universe(isl_space_copy(bset->dim));
 526
 527         return isl_morph_alloc(dom, ran, map, inv);
 528 error:
 529         isl_mat_free(H);
 530         isl_mat_free(B);
 531         isl_vec_free(d);
 532         return NULL;
 533 }
 534
 535 /* Add stride constraints to "bset" based on the inverse mapping
 536  * that was plugged in.  In particular, if morph maps x' to x,
 537  * the the constraints of the original input
 538  *
 539  *      A x' + b >= 0
 540  *
 541  * have been rewritten to
 542  *
 543  *      A inv x + b >= 0
 544  *
 545  * However, this substitution may loose information on the integrality of x',
 546  * so we need to impose that
 547  *
 548  *      inv x
 549  *
 550  * is integral.  If inv = B/d, this means that we need to impose that
 551  *
 552  *      B x = 0         mod d
 553  *
 554  * or
 555  *
 556  *      exists alpha in Z^m: B x = d alpha
 557  *
 558  */
 559 static __isl_give isl_basic_set *add_strides(__isl_take isl_basic_set *bset,
 560         __isl_keep isl_morph *morph)
 561 {
 562         int i, div, k;
 563         isl_int gcd;
 564
 565         if (isl_int_is_one(morph->inv->row[0][0]))
 566                 return bset;
 567
 568         isl_int_init(gcd);
 569
 570         for (i = 0; 1 + i < morph->inv->n_row; ++i) {
 571                 isl_seq_gcd(morph->inv->row[1 + i], morph->inv->n_col, &gcd);
 572                 if (isl_int_is_divisible_by(gcd, morph->inv->row[0][0]))
 573                         continue;
 574                 div = isl_basic_set_alloc_div(bset);
 575                 if (div < 0)
 576                         goto error;
 577                 k = isl_basic_set_alloc_equality(bset);
 578                 if (k < 0)
 579                         goto error;
 580                 isl_seq_cpy(bset->eq[k], morph->inv->row[1 + i],
 581                             morph->inv->n_col);
 582                 isl_seq_clr(bset->eq[k] + morph->inv->n_col, bset->n_div);
 583                 isl_int_set(bset->eq[k][morph->inv->n_col + div],
 584                             morph->inv->row[0][0]);
 585         }
 586
 587         isl_int_clear(gcd);
 588
 589         return bset;
 590 error:
 591         isl_int_clear(gcd);
 592         isl_basic_set_free(bset);
 593         return NULL;
 594 }
 595
 596 /* Apply the morphism to the basic set.
 597  * We basically just compute the preimage of "bset" under the inverse mapping
 598  * in morph, add in stride constraints and intersect with the range
 599  * of the morphism.
 600  */
 601 __isl_give isl_basic_set *isl_morph_basic_set(__isl_take isl_morph *morph,
 602         __isl_take isl_basic_set *bset)
 603 {
 604         isl_basic_set *res = NULL;
 605         isl_mat *mat = NULL;
 606         int i, k;
 607         int max_stride;
 608
 609         if (!morph || !bset)
 610                 goto error;
 611
 612         isl_assert(bset->ctx, isl_space_is_equal(bset->dim, morph->dom->dim),
 613                     goto error);
 614
 615         max_stride = morph->inv->n_row - 1;
 616         if (isl_int_is_one(morph->inv->row[0][0]))
 617                 max_stride = 0;
 618         res = isl_basic_set_alloc_space(isl_space_copy(morph->ran->dim),
 619                 bset->n_div + max_stride, bset->n_eq + max_stride, bset->n_ineq);
 620
 621         for (i = 0; i < bset->n_div; ++i)
 622                 if (isl_basic_set_alloc_div(res) < 0)
 623                         goto error;
 624
 625         mat = isl_mat_sub_alloc6(bset->ctx, bset->eq, 0, bset->n_eq,
 626                                         0, morph->inv->n_row);
 627         mat = isl_mat_product(mat, isl_mat_copy(morph->inv));
 628         if (!mat)
 629                 goto error;
 630         for (i = 0; i < bset->n_eq; ++i) {
 631                 k = isl_basic_set_alloc_equality(res);
 632                 if (k < 0)
 633                         goto error;
 634                 isl_seq_cpy(res->eq[k], mat->row[i], mat->n_col);
 635                 isl_seq_scale(res->eq[k] + mat->n_col, bset->eq[i] + mat->n_col,
 636                                 morph->inv->row[0][0], bset->n_div);
 637         }
 638         isl_mat_free(mat);
 639
 640         mat = isl_mat_sub_alloc6(bset->ctx, bset->ineq, 0, bset->n_ineq,
 641                                         0, morph->inv->n_row);
 642         mat = isl_mat_product(mat, isl_mat_copy(morph->inv));
 643         if (!mat)
 644                 goto error;
 645         for (i = 0; i < bset->n_ineq; ++i) {
 646                 k = isl_basic_set_alloc_inequality(res);
 647                 if (k < 0)
 648                         goto error;
 649                 isl_seq_cpy(res->ineq[k], mat->row[i], mat->n_col);
 650                 isl_seq_scale(res->ineq[k] + mat->n_col,
 651                                 bset->ineq[i] + mat->n_col,
 652                                 morph->inv->row[0][0], bset->n_div);
 653         }
 654         isl_mat_free(mat);
 655
 656         mat = isl_mat_sub_alloc6(bset->ctx, bset->div, 0, bset->n_div,
 657                                         1, morph->inv->n_row);
 658         mat = isl_mat_product(mat, isl_mat_copy(morph->inv));
 659         if (!mat)
 660                 goto error;
 661         for (i = 0; i < bset->n_div; ++i) {
 662                 isl_int_mul(res->div[i][0],
 663                                 morph->inv->row[0][0], bset->div[i][0]);
 664                 isl_seq_cpy(res->div[i] + 1, mat->row[i], mat->n_col);
 665                 isl_seq_scale(res->div[i] + 1 + mat->n_col,
 666                                 bset->div[i] + 1 + mat->n_col,
 667                                 morph->inv->row[0][0], bset->n_div);
 668         }
 669         isl_mat_free(mat);
 670
 671         res = add_strides(res, morph);
 672
 673         if (isl_basic_set_is_rational(bset))
 674                 res = isl_basic_set_set_rational(res);
 675
 676         res = isl_basic_set_simplify(res);
 677         res = isl_basic_set_finalize(res);
 678
 679         res = isl_basic_set_intersect(res, isl_basic_set_copy(morph->ran));
 680
 681         isl_morph_free(morph);
 682         isl_basic_set_free(bset);
 683         return res;
 684 error:
 685         isl_mat_free(mat);
 686         isl_morph_free(morph);
 687         isl_basic_set_free(bset);
 688         isl_basic_set_free(res);
 689         return NULL;
 690 }
 691
 692 /* Apply the morphism to the set.
 693  */
 694 __isl_give isl_set *isl_morph_set(__isl_take isl_morph *morph,
 695         __isl_take isl_set *set)
 696 {
 697         int i;
 698
 699         if (!morph || !set)
 700                 goto error;
 701
 702         isl_assert(set->ctx, isl_space_is_equal(set->dim, morph->dom->dim), goto error);
 703
 704         set = isl_set_cow(set);
 705         if (!set)
 706                 goto error;
 707
 708         isl_space_free(set->dim);
 709         set->dim = isl_space_copy(morph->ran->dim);
 710         if (!set->dim)
 711                 goto error;
 712
 713         for (i = 0; i < set->n; ++i) {
 714                 set->p[i] = isl_morph_basic_set(isl_morph_copy(morph), set->p[i]);
 715                 if (!set->p[i])
 716                         goto error;
 717         }
 718
 719         isl_morph_free(morph);
 720
 721         ISL_F_CLR(set, ISL_SET_NORMALIZED);
 722
 723         return set;
 724 error:
 725         isl_set_free(set);
 726         isl_morph_free(morph);
 727         return NULL;
 728 }
 729
 730 /* Construct a morphism that first does morph2 and then morph1.
 731  */
 732 __isl_give isl_morph *isl_morph_compose(__isl_take isl_morph *morph1,
 733         __isl_take isl_morph *morph2)
 734 {
 735         isl_mat *map, *inv;
 736         isl_basic_set *dom, *ran;
 737
 738         if (!morph1 || !morph2)
 739                 goto error;
 740
 741         map = isl_mat_product(isl_mat_copy(morph1->map), isl_mat_copy(morph2->map));
 742         inv = isl_mat_product(isl_mat_copy(morph2->inv), isl_mat_copy(morph1->inv));
 743         dom = isl_morph_basic_set(isl_morph_inverse(isl_morph_copy(morph2)),
 744                                   isl_basic_set_copy(morph1->dom));
 745         dom = isl_basic_set_intersect(dom, isl_basic_set_copy(morph2->dom));
 746         ran = isl_morph_basic_set(isl_morph_copy(morph1),
 747                                   isl_basic_set_copy(morph2->ran));
 748         ran = isl_basic_set_intersect(ran, isl_basic_set_copy(morph1->ran));
 749
 750         isl_morph_free(morph1);
 751         isl_morph_free(morph2);
 752
 753         return isl_morph_alloc(dom, ran, map, inv);
 754 error:
 755         isl_morph_free(morph1);
 756         isl_morph_free(morph2);
 757         return NULL;
 758 }
 759
 760 __isl_give isl_morph *isl_morph_inverse(__isl_take isl_morph *morph)
 761 {
 762         isl_basic_set *bset;
 763         isl_mat *mat;
 764
 765         morph = isl_morph_cow(morph);
 766         if (!morph)
 767                 return NULL;
 768
 769         bset = morph->dom;
 770         morph->dom = morph->ran;
 771         morph->ran = bset;
 772
 773         mat = morph->map;
 774         morph->map = morph->inv;
 775         morph->inv = mat;
 776
 777         return morph;
 778 }
 779
 780 __isl_give isl_morph *isl_basic_set_full_compression(
 781         __isl_keep isl_basic_set *bset)
 782 {
 783         isl_morph *morph, *morph2;
 784
 785         bset = isl_basic_set_copy(bset);
 786
 787         morph = isl_basic_set_variable_compression(bset, isl_dim_param);
 788         bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
 789
 790         morph2 = isl_basic_set_parameter_compression(bset);
 791         bset = isl_morph_basic_set(isl_morph_copy(morph2), bset);
 792
 793         morph = isl_morph_compose(morph2, morph);
 794
 795         morph2 = isl_basic_set_variable_compression(bset, isl_dim_set);
 796         isl_basic_set_free(bset);
 797
 798         morph = isl_morph_compose(morph2, morph);
 799
 800         return morph;
 801 }
 802
 803 __isl_give isl_vec *isl_morph_vec(__isl_take isl_morph *morph,
 804         __isl_take isl_vec *vec)
 805 {
 806         if (!morph)
 807                 goto error;
 808
 809         vec = isl_mat_vec_product(isl_mat_copy(morph->map), vec);
 810
 811         isl_morph_free(morph);
 812         return vec;
 813 error:
 814         isl_morph_free(morph);
 815         isl_vec_free(vec);
 816         return NULL;
 817 }