isl_basic_set_compute_vertices: explicitly turn chambers into parameter sets
[platform/upstream/isl.git] / 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 }