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