2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
11 #include <isl_morph.h>
13 #include <isl_map_private.h>
14 #include <isl_dim_private.h>
15 #include <isl_equalities.h>
17 __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)
23 if (!dom || !ran || !map || !inv)
26 morph = isl_alloc_type(dom->ctx, struct isl_morph);
38 isl_basic_set_free(dom);
39 isl_basic_set_free(ran);
45 __isl_give isl_morph *isl_morph_copy(__isl_keep isl_morph *morph)
54 __isl_give isl_morph *isl_morph_dup(__isl_keep isl_morph *morph)
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));
64 __isl_give isl_morph *isl_morph_cow(__isl_take isl_morph *morph)
72 return isl_morph_dup(morph);
75 void isl_morph_free(__isl_take isl_morph *morph)
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);
90 __isl_give isl_dim *isl_morph_get_ran_dim(__isl_keep isl_morph *morph)
95 return isl_dim_copy(morph->ran->dim);
98 unsigned isl_morph_dom_dim(__isl_keep isl_morph *morph, enum isl_dim_type type)
103 return isl_basic_set_dim(morph->dom, type);
106 unsigned isl_morph_ran_dim(__isl_keep isl_morph *morph, enum isl_dim_type type)
111 return isl_basic_set_dim(morph->ran, type);
114 __isl_give isl_morph *isl_morph_remove_dom_dims(__isl_take isl_morph *morph,
115 enum isl_dim_type type, unsigned first, unsigned n)
122 morph = isl_morph_cow(morph);
126 dom_offset = 1 + isl_dim_offset(morph->dom->dim, type);
128 morph->dom = isl_basic_set_remove(morph->dom, type, first, n);
130 morph->map = isl_mat_drop_cols(morph->map, dom_offset + first, n);
132 morph->inv = isl_mat_drop_rows(morph->inv, dom_offset + first, n);
134 if (morph->dom && morph->ran && morph->map && morph->inv)
137 isl_morph_free(morph);
141 __isl_give isl_morph *isl_morph_remove_ran_dims(__isl_take isl_morph *morph,
142 enum isl_dim_type type, unsigned first, unsigned n)
149 morph = isl_morph_cow(morph);
153 ran_offset = 1 + isl_dim_offset(morph->ran->dim, type);
155 morph->ran = isl_basic_set_remove(morph->ran, type, first, n);
157 morph->map = isl_mat_drop_rows(morph->map, ran_offset + first, n);
159 morph->inv = isl_mat_drop_cols(morph->inv, ran_offset + first, n);
161 if (morph->dom && morph->ran && morph->map && morph->inv)
164 isl_morph_free(morph);
168 void isl_morph_dump(__isl_take isl_morph *morph, FILE *out)
173 isl_basic_set_print(morph->dom, out, 0, "", "", ISL_FORMAT_ISL);
174 isl_basic_set_print(morph->ran, out, 0, "", "", ISL_FORMAT_ISL);
175 isl_mat_dump(morph->map, out, 4);
176 isl_mat_dump(morph->inv, out, 4);
179 __isl_give isl_morph *isl_morph_identity(__isl_keep isl_basic_set *bset)
182 isl_basic_set *universe;
188 total = isl_basic_set_total_dim(bset);
189 id = isl_mat_identity(bset->ctx, 1 + total);
190 universe = isl_basic_set_universe(isl_dim_copy(bset->dim));
192 return isl_morph_alloc(universe, isl_basic_set_copy(universe),
193 id, isl_mat_copy(id));
196 /* Create a(n identity) morphism between empty sets of the same dimension
199 __isl_give isl_morph *isl_morph_empty(__isl_keep isl_basic_set *bset)
202 isl_basic_set *empty;
208 total = isl_basic_set_total_dim(bset);
209 id = isl_mat_identity(bset->ctx, 1 + total);
210 empty = isl_basic_set_empty(isl_dim_copy(bset->dim));
212 return isl_morph_alloc(empty, isl_basic_set_copy(empty),
213 id, isl_mat_copy(id));
216 /* Given a matrix that maps a (possibly) parametric domain to
217 * a parametric domain, add in rows that map the "nparam" parameters onto
220 static __isl_give isl_mat *insert_parameter_rows(__isl_take isl_mat *mat,
230 mat = isl_mat_insert_rows(mat, 1, nparam);
234 for (i = 0; i < nparam; ++i) {
235 isl_seq_clr(mat->row[1 + i], mat->n_col);
236 isl_int_set(mat->row[1 + i][1 + i], mat->row[0][0]);
242 /* Construct a basic set described by the "n" equalities of "bset" starting
245 static __isl_give isl_basic_set *copy_equalities(__isl_keep isl_basic_set *bset,
246 unsigned first, unsigned n)
252 isl_assert(bset->ctx, bset->n_div == 0, return NULL);
254 total = isl_basic_set_total_dim(bset);
255 eq = isl_basic_set_alloc_dim(isl_dim_copy(bset->dim), 0, n, 0);
258 for (i = 0; i < n; ++i) {
259 k = isl_basic_set_alloc_equality(eq);
262 isl_seq_cpy(eq->eq[k], bset->eq[first + k], 1 + total);
267 isl_basic_set_free(eq);
271 /* Given a basic set, exploit the equalties in the a basic set to construct
272 * a morphishm that maps the basic set to a lower-dimensional space.
273 * Specifically, the morphism reduces the number of dimensions of type "type".
275 * This function is a slight generalization of isl_mat_variable_compression
276 * in that it allows the input to be parametric and that it allows for the
277 * compression of either parameters or set variables.
279 * We first select the equalities of interest, that is those that involve
280 * variables of type "type" and no later variables.
281 * Denote those equalities as
285 * where C(p) depends on the parameters if type == isl_dim_set and
286 * is a constant if type == isl_dim_param.
288 * First compute the (left) Hermite normal form of M,
290 * M [U1 U2] = M U = H = [H1 0]
292 * M = H Q = [H1 0] [Q1]
295 * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
296 * Define the transformed variables as
298 * x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
301 * The equalities then become
303 * -C(p) + H1 x1' = 0 or x1' = H1^{-1} C(p) = C'(p)
305 * If the denominator of the constant term does not divide the
306 * the common denominator of the parametric terms, then every
307 * integer point is mapped to a non-integer point and then the original set has no
308 * integer solutions (since the x' are a unimodular transformation
309 * of the x). In this case, an empty morphism is returned.
310 * Otherwise, the transformation is given by
312 * x = U1 H1^{-1} C(p) + U2 x2'
314 * The inverse transformation is simply
318 * Both matrices are extended to map the full original space to the full
321 __isl_give isl_morph *isl_basic_set_variable_compression(
322 __isl_keep isl_basic_set *bset, enum isl_dim_type type)
331 isl_mat *H, *U, *Q, *C = NULL, *H1, *U1, *U2;
332 isl_basic_set *dom, *ran;
337 if (isl_basic_set_fast_is_empty(bset))
338 return isl_morph_empty(bset);
340 isl_assert(bset->ctx, bset->n_div == 0, return NULL);
342 otype = 1 + isl_dim_offset(bset->dim, type);
343 ntype = isl_basic_set_dim(bset, type);
344 orest = otype + ntype;
345 nrest = isl_basic_set_total_dim(bset) - (orest - 1);
347 for (f_eq = 0; f_eq < bset->n_eq; ++f_eq)
348 if (isl_seq_first_non_zero(bset->eq[f_eq] + orest, nrest) == -1)
350 for (n_eq = 0; f_eq + n_eq < bset->n_eq; ++n_eq)
351 if (isl_seq_first_non_zero(bset->eq[f_eq + n_eq] + otype, ntype) == -1)
354 return isl_morph_identity(bset);
356 H = isl_mat_sub_alloc(bset->ctx, bset->eq, f_eq, n_eq, otype, ntype);
357 H = isl_mat_left_hermite(H, 0, &U, &Q);
360 Q = isl_mat_drop_rows(Q, 0, n_eq);
361 Q = isl_mat_diagonal(isl_mat_identity(bset->ctx, otype), Q);
362 Q = isl_mat_diagonal(Q, isl_mat_identity(bset->ctx, nrest));
363 C = isl_mat_alloc(bset->ctx, 1 + n_eq, otype);
366 isl_int_set_si(C->row[0][0], 1);
367 isl_seq_clr(C->row[0] + 1, otype - 1);
368 isl_mat_sub_neg(C->ctx, C->row + 1, bset->eq + f_eq, n_eq, 0, 0, otype);
369 H1 = isl_mat_sub_alloc(H->ctx, H->row, 0, H->n_row, 0, H->n_row);
370 H1 = isl_mat_lin_to_aff(H1);
371 C = isl_mat_inverse_product(H1, C);
376 if (!isl_int_is_one(C->row[0][0])) {
381 for (i = 0; i < n_eq; ++i) {
382 isl_seq_gcd(C->row[1 + i] + 1, otype - 1, &g);
383 isl_int_gcd(g, g, C->row[0][0]);
384 if (!isl_int_is_divisible_by(C->row[1 + i][0], g))
393 return isl_morph_empty(bset);
396 C = isl_mat_normalize(C);
399 U1 = isl_mat_sub_alloc(U->ctx, U->row, 0, U->n_row, 0, n_eq);
400 U1 = isl_mat_lin_to_aff(U1);
401 U2 = isl_mat_sub_alloc(U->ctx, U->row, 0, U->n_row, n_eq, U->n_row - n_eq);
402 U2 = isl_mat_lin_to_aff(U2);
405 C = isl_mat_product(U1, C);
406 C = isl_mat_aff_direct_sum(C, U2);
407 C = insert_parameter_rows(C, otype - 1);
408 C = isl_mat_diagonal(C, isl_mat_identity(bset->ctx, nrest));
410 dim = isl_dim_copy(bset->dim);
411 dim = isl_dim_drop(dim, type, 0, ntype);
412 dim = isl_dim_add(dim, type, ntype - n_eq);
413 ran = isl_basic_set_universe(dim);
414 dom = copy_equalities(bset, f_eq, n_eq);
416 return isl_morph_alloc(dom, ran, Q, C);
425 /* Construct a parameter compression for "bset".
426 * We basically just call isl_mat_parameter_compression with the right input
427 * and then extend the resulting matrix to include the variables.
429 * Let the equalities be given as
433 * and let [H 0] be the Hermite Normal Form of A, then
437 * needs to be integer, so we impose that each row is divisible by
440 __isl_give isl_morph *isl_basic_set_parameter_compression(
441 __isl_keep isl_basic_set *bset)
449 isl_basic_set *dom, *ran;
454 if (isl_basic_set_fast_is_empty(bset))
455 return isl_morph_empty(bset);
457 return isl_morph_identity(bset);
459 isl_assert(bset->ctx, bset->n_div == 0, return NULL);
462 nparam = isl_basic_set_dim(bset, isl_dim_param);
463 nvar = isl_basic_set_dim(bset, isl_dim_set);
465 isl_assert(bset->ctx, n_eq <= nvar, return NULL);
467 d = isl_vec_alloc(bset->ctx, n_eq);
468 B = isl_mat_sub_alloc(bset->ctx, bset->eq, 0, n_eq, 0, 1 + nparam);
469 H = isl_mat_sub_alloc(bset->ctx, bset->eq, 0, n_eq, 1 + nparam, nvar);
470 H = isl_mat_left_hermite(H, 0, NULL, NULL);
471 H = isl_mat_drop_cols(H, n_eq, nvar - n_eq);
472 H = isl_mat_lin_to_aff(H);
473 H = isl_mat_right_inverse(H);
476 isl_seq_set(d->el, H->row[0][0], d->size);
477 H = isl_mat_drop_rows(H, 0, 1);
478 H = isl_mat_drop_cols(H, 0, 1);
479 B = isl_mat_product(H, B);
480 inv = isl_mat_parameter_compression(B, d);
481 inv = isl_mat_diagonal(inv, isl_mat_identity(bset->ctx, nvar));
482 map = isl_mat_right_inverse(isl_mat_copy(inv));
484 dom = isl_basic_set_universe(isl_dim_copy(bset->dim));
485 ran = isl_basic_set_universe(isl_dim_copy(bset->dim));
487 return isl_morph_alloc(dom, ran, map, inv);
495 /* Add stride constraints to "bset" based on the inverse mapping
496 * that was plugged in. In particular, if morph maps x' to x,
497 * the the constraints of the original input
501 * have been rewritten to
505 * However, this substitution may loose information on the integrality of x',
506 * so we need to impose that
510 * is integral. If inv = B/d, this means that we need to impose that
516 * exists alpha in Z^m: B x = d alpha
519 static __isl_give isl_basic_set *add_strides(__isl_take isl_basic_set *bset,
520 __isl_keep isl_morph *morph)
525 if (isl_int_is_one(morph->inv->row[0][0]))
530 for (i = 0; 1 + i < morph->inv->n_row; ++i) {
531 isl_seq_gcd(morph->inv->row[1 + i], morph->inv->n_col, &gcd);
532 if (isl_int_is_divisible_by(gcd, morph->inv->row[0][0]))
534 div = isl_basic_set_alloc_div(bset);
537 k = isl_basic_set_alloc_equality(bset);
540 isl_seq_cpy(bset->eq[k], morph->inv->row[1 + i],
542 isl_seq_clr(bset->eq[k] + morph->inv->n_col, bset->n_div);
543 isl_int_set(bset->eq[k][morph->inv->n_col + div],
544 morph->inv->row[0][0]);
552 isl_basic_set_free(bset);
556 /* Apply the morphism to the basic set.
557 * We basically just compute the preimage of "bset" under the inverse mapping
558 * in morph, add in stride constraints and intersect with the range
561 __isl_give isl_basic_set *isl_morph_basic_set(__isl_take isl_morph *morph,
562 __isl_take isl_basic_set *bset)
564 isl_basic_set *res = NULL;
572 isl_assert(bset->ctx, isl_dim_equal(bset->dim, morph->dom->dim),
575 max_stride = morph->inv->n_row - 1;
576 if (isl_int_is_one(morph->inv->row[0][0]))
578 res = isl_basic_set_alloc_dim(isl_dim_copy(morph->ran->dim),
579 bset->n_div + max_stride, bset->n_eq + max_stride, bset->n_ineq);
581 for (i = 0; i < bset->n_div; ++i)
582 if (isl_basic_set_alloc_div(res) < 0)
585 mat = isl_mat_sub_alloc(bset->ctx, bset->eq, 0, bset->n_eq,
586 0, morph->inv->n_row);
587 mat = isl_mat_product(mat, isl_mat_copy(morph->inv));
590 for (i = 0; i < bset->n_eq; ++i) {
591 k = isl_basic_set_alloc_equality(res);
594 isl_seq_cpy(res->eq[k], mat->row[i], mat->n_col);
595 isl_seq_scale(res->eq[k] + mat->n_col, bset->eq[i] + mat->n_col,
596 morph->inv->row[0][0], bset->n_div);
600 mat = isl_mat_sub_alloc(bset->ctx, bset->ineq, 0, bset->n_ineq,
601 0, morph->inv->n_row);
602 mat = isl_mat_product(mat, isl_mat_copy(morph->inv));
605 for (i = 0; i < bset->n_ineq; ++i) {
606 k = isl_basic_set_alloc_inequality(res);
609 isl_seq_cpy(res->ineq[k], mat->row[i], mat->n_col);
610 isl_seq_scale(res->ineq[k] + mat->n_col,
611 bset->ineq[i] + mat->n_col,
612 morph->inv->row[0][0], bset->n_div);
616 mat = isl_mat_sub_alloc(bset->ctx, bset->div, 0, bset->n_div,
617 1, morph->inv->n_row);
618 mat = isl_mat_product(mat, isl_mat_copy(morph->inv));
621 for (i = 0; i < bset->n_div; ++i) {
622 isl_int_mul(res->div[i][0],
623 morph->inv->row[0][0], bset->div[i][0]);
624 isl_seq_cpy(res->div[i] + 1, mat->row[i], mat->n_col);
625 isl_seq_scale(res->div[i] + 1 + mat->n_col,
626 bset->div[i] + 1 + mat->n_col,
627 morph->inv->row[0][0], bset->n_div);
631 res = add_strides(res, morph);
633 res = isl_basic_set_simplify(res);
634 res = isl_basic_set_finalize(res);
636 res = isl_basic_set_intersect(res, isl_basic_set_copy(morph->ran));
638 isl_morph_free(morph);
639 isl_basic_set_free(bset);
643 isl_morph_free(morph);
644 isl_basic_set_free(bset);
645 isl_basic_set_free(res);
649 /* Apply the morphism to the set.
651 __isl_give isl_set *isl_morph_set(__isl_take isl_morph *morph,
652 __isl_take isl_set *set)
659 isl_assert(set->ctx, isl_dim_equal(set->dim, morph->dom->dim), goto error);
661 set = isl_set_cow(set);
665 isl_dim_free(set->dim);
666 set->dim = isl_dim_copy(morph->ran->dim);
670 for (i = 0; i < set->n; ++i) {
671 set->p[i] = isl_morph_basic_set(isl_morph_copy(morph), set->p[i]);
676 isl_morph_free(morph);
678 ISL_F_CLR(set, ISL_SET_NORMALIZED);
683 isl_morph_free(morph);
687 /* Construct a morphism that first does morph2 and then morph1.
689 __isl_give isl_morph *isl_morph_compose(__isl_take isl_morph *morph1,
690 __isl_take isl_morph *morph2)
693 isl_basic_set *dom, *ran;
695 if (!morph1 || !morph2)
698 map = isl_mat_product(isl_mat_copy(morph1->map), isl_mat_copy(morph2->map));
699 inv = isl_mat_product(isl_mat_copy(morph2->inv), isl_mat_copy(morph1->inv));
700 dom = isl_morph_basic_set(isl_morph_inverse(isl_morph_copy(morph2)),
701 isl_basic_set_copy(morph1->dom));
702 dom = isl_basic_set_intersect(dom, isl_basic_set_copy(morph2->dom));
703 ran = isl_morph_basic_set(isl_morph_copy(morph1),
704 isl_basic_set_copy(morph2->ran));
705 ran = isl_basic_set_intersect(ran, isl_basic_set_copy(morph1->ran));
707 isl_morph_free(morph1);
708 isl_morph_free(morph2);
710 return isl_morph_alloc(dom, ran, map, inv);
712 isl_morph_free(morph1);
713 isl_morph_free(morph2);
717 __isl_give isl_morph *isl_morph_inverse(__isl_take isl_morph *morph)
722 morph = isl_morph_cow(morph);
727 morph->dom = morph->ran;
731 morph->map = morph->inv;
737 __isl_give isl_morph *isl_basic_set_full_compression(
738 __isl_keep isl_basic_set *bset)
740 isl_morph *morph, *morph2;
742 bset = isl_basic_set_copy(bset);
744 morph = isl_basic_set_variable_compression(bset, isl_dim_param);
745 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
747 morph2 = isl_basic_set_parameter_compression(bset);
748 bset = isl_morph_basic_set(isl_morph_copy(morph2), bset);
750 morph = isl_morph_compose(morph2, morph);
752 morph2 = isl_basic_set_variable_compression(bset, isl_dim_set);
753 isl_basic_set_free(bset);
755 morph = isl_morph_compose(morph2, morph);