2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the MIT 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_map_private.h>
12 #include <isl_morph.h>
14 #include <isl_mat_private.h>
15 #include <isl_space_private.h>
16 #include <isl_equalities.h>
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
)
24 if (!dom
|| !ran
|| !map
|| !inv
)
27 morph
= isl_alloc_type(dom
->ctx
, struct isl_morph
);
39 isl_basic_set_free(dom
);
40 isl_basic_set_free(ran
);
46 __isl_give isl_morph
*isl_morph_copy(__isl_keep isl_morph
*morph
)
55 __isl_give isl_morph
*isl_morph_dup(__isl_keep isl_morph
*morph
)
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
));
65 __isl_give isl_morph
*isl_morph_cow(__isl_take isl_morph
*morph
)
73 return isl_morph_dup(morph
);
76 void isl_morph_free(__isl_take isl_morph
*morph
)
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
);
91 __isl_give isl_space
*isl_morph_get_ran_space(__isl_keep isl_morph
*morph
)
96 return isl_space_copy(morph
->ran
->dim
);
99 unsigned isl_morph_dom_dim(__isl_keep isl_morph
*morph
, enum isl_dim_type type
)
104 return isl_basic_set_dim(morph
->dom
, type
);
107 unsigned isl_morph_ran_dim(__isl_keep isl_morph
*morph
, enum isl_dim_type type
)
112 return isl_basic_set_dim(morph
->ran
, type
);
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
)
123 morph
= isl_morph_cow(morph
);
127 dom_offset
= 1 + isl_space_offset(morph
->dom
->dim
, type
);
129 morph
->dom
= isl_basic_set_remove_dims(morph
->dom
, type
, first
, n
);
131 morph
->map
= isl_mat_drop_cols(morph
->map
, dom_offset
+ first
, n
);
133 morph
->inv
= isl_mat_drop_rows(morph
->inv
, dom_offset
+ first
, n
);
135 if (morph
->dom
&& morph
->ran
&& morph
->map
&& morph
->inv
)
138 isl_morph_free(morph
);
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
)
150 morph
= isl_morph_cow(morph
);
154 ran_offset
= 1 + isl_space_offset(morph
->ran
->dim
, type
);
156 morph
->ran
= isl_basic_set_remove_dims(morph
->ran
, type
, first
, n
);
158 morph
->map
= isl_mat_drop_rows(morph
->map
, ran_offset
+ first
, n
);
160 morph
->inv
= isl_mat_drop_cols(morph
->inv
, ran_offset
+ first
, n
);
162 if (morph
->dom
&& morph
->ran
&& morph
->map
&& morph
->inv
)
165 isl_morph_free(morph
);
169 /* Project domain of morph onto its parameter domain.
171 __isl_give isl_morph
*isl_morph_dom_params(__isl_take isl_morph
*morph
)
175 morph
= isl_morph_cow(morph
);
178 n
= isl_basic_set_dim(morph
->dom
, isl_dim_set
);
179 morph
= isl_morph_remove_dom_dims(morph
, isl_dim_set
, 0, n
);
182 morph
->dom
= isl_basic_set_params(morph
->dom
);
186 isl_morph_free(morph
);
190 /* Project range of morph onto its parameter domain.
192 __isl_give isl_morph
*isl_morph_ran_params(__isl_take isl_morph
*morph
)
196 morph
= isl_morph_cow(morph
);
199 n
= isl_basic_set_dim(morph
->ran
, isl_dim_set
);
200 morph
= isl_morph_remove_ran_dims(morph
, isl_dim_set
, 0, n
);
203 morph
->ran
= isl_basic_set_params(morph
->ran
);
207 isl_morph_free(morph
);
211 void isl_morph_print_internal(__isl_take isl_morph
*morph
, FILE *out
)
216 isl_basic_set_print(morph
->dom
, out
, 0, "", "", ISL_FORMAT_ISL
);
217 isl_basic_set_print(morph
->ran
, out
, 0, "", "", ISL_FORMAT_ISL
);
218 isl_mat_print_internal(morph
->map
, out
, 4);
219 isl_mat_print_internal(morph
->inv
, out
, 4);
222 void isl_morph_dump(__isl_take isl_morph
*morph
)
224 isl_morph_print_internal(morph
, stderr
);
227 __isl_give isl_morph
*isl_morph_identity(__isl_keep isl_basic_set
*bset
)
230 isl_basic_set
*universe
;
236 total
= isl_basic_set_total_dim(bset
);
237 id
= isl_mat_identity(bset
->ctx
, 1 + total
);
238 universe
= isl_basic_set_universe(isl_space_copy(bset
->dim
));
240 return isl_morph_alloc(universe
, isl_basic_set_copy(universe
),
241 id
, isl_mat_copy(id
));
244 /* Create a(n identity) morphism between empty sets of the same dimension
247 __isl_give isl_morph
*isl_morph_empty(__isl_keep isl_basic_set
*bset
)
250 isl_basic_set
*empty
;
256 total
= isl_basic_set_total_dim(bset
);
257 id
= isl_mat_identity(bset
->ctx
, 1 + total
);
258 empty
= isl_basic_set_empty(isl_space_copy(bset
->dim
));
260 return isl_morph_alloc(empty
, isl_basic_set_copy(empty
),
261 id
, isl_mat_copy(id
));
264 /* Given a matrix that maps a (possibly) parametric domain to
265 * a parametric domain, add in rows that map the "nparam" parameters onto
268 static __isl_give isl_mat
*insert_parameter_rows(__isl_take isl_mat
*mat
,
278 mat
= isl_mat_insert_rows(mat
, 1, nparam
);
282 for (i
= 0; i
< nparam
; ++i
) {
283 isl_seq_clr(mat
->row
[1 + i
], mat
->n_col
);
284 isl_int_set(mat
->row
[1 + i
][1 + i
], mat
->row
[0][0]);
290 /* Construct a basic set described by the "n" equalities of "bset" starting
293 static __isl_give isl_basic_set
*copy_equalities(__isl_keep isl_basic_set
*bset
,
294 unsigned first
, unsigned n
)
300 isl_assert(bset
->ctx
, bset
->n_div
== 0, return NULL
);
302 total
= isl_basic_set_total_dim(bset
);
303 eq
= isl_basic_set_alloc_space(isl_space_copy(bset
->dim
), 0, n
, 0);
306 for (i
= 0; i
< n
; ++i
) {
307 k
= isl_basic_set_alloc_equality(eq
);
310 isl_seq_cpy(eq
->eq
[k
], bset
->eq
[first
+ k
], 1 + total
);
315 isl_basic_set_free(eq
);
319 /* Given a basic set, exploit the equalties in the basic set to construct
320 * a morphishm that maps the basic set to a lower-dimensional space.
321 * Specifically, the morphism reduces the number of dimensions of type "type".
323 * This function is a slight generalization of isl_mat_variable_compression
324 * in that it allows the input to be parametric and that it allows for the
325 * compression of either parameters or set variables.
327 * We first select the equalities of interest, that is those that involve
328 * variables of type "type" and no later variables.
329 * Denote those equalities as
333 * where C(p) depends on the parameters if type == isl_dim_set and
334 * is a constant if type == isl_dim_param.
336 * First compute the (left) Hermite normal form of M,
338 * M [U1 U2] = M U = H = [H1 0]
340 * M = H Q = [H1 0] [Q1]
343 * with U, Q unimodular, Q = U^{-1} (and H lower triangular).
344 * Define the transformed variables as
346 * x = [U1 U2] [ x1' ] = [U1 U2] [Q1] x
349 * The equalities then become
351 * -C(p) + H1 x1' = 0 or x1' = H1^{-1} C(p) = C'(p)
353 * If the denominator of the constant term does not divide the
354 * the common denominator of the parametric terms, then every
355 * integer point is mapped to a non-integer point and then the original set has no
356 * integer solutions (since the x' are a unimodular transformation
357 * of the x). In this case, an empty morphism is returned.
358 * Otherwise, the transformation is given by
360 * x = U1 H1^{-1} C(p) + U2 x2'
362 * The inverse transformation is simply
366 * Both matrices are extended to map the full original space to the full
369 __isl_give isl_morph
*isl_basic_set_variable_compression(
370 __isl_keep isl_basic_set
*bset
, enum isl_dim_type type
)
378 isl_mat
*H
, *U
, *Q
, *C
= NULL
, *H1
, *U1
, *U2
;
379 isl_basic_set
*dom
, *ran
;
384 if (isl_basic_set_plain_is_empty(bset
))
385 return isl_morph_empty(bset
);
387 isl_assert(bset
->ctx
, bset
->n_div
== 0, return NULL
);
389 otype
= 1 + isl_space_offset(bset
->dim
, type
);
390 ntype
= isl_basic_set_dim(bset
, type
);
391 orest
= otype
+ ntype
;
392 nrest
= isl_basic_set_total_dim(bset
) - (orest
- 1);
394 for (f_eq
= 0; f_eq
< bset
->n_eq
; ++f_eq
)
395 if (isl_seq_first_non_zero(bset
->eq
[f_eq
] + orest
, nrest
) == -1)
397 for (n_eq
= 0; f_eq
+ n_eq
< bset
->n_eq
; ++n_eq
)
398 if (isl_seq_first_non_zero(bset
->eq
[f_eq
+ n_eq
] + otype
, ntype
) == -1)
401 return isl_morph_identity(bset
);
403 H
= isl_mat_sub_alloc6(bset
->ctx
, bset
->eq
, f_eq
, n_eq
, otype
, ntype
);
404 H
= isl_mat_left_hermite(H
, 0, &U
, &Q
);
407 Q
= isl_mat_drop_rows(Q
, 0, n_eq
);
408 Q
= isl_mat_diagonal(isl_mat_identity(bset
->ctx
, otype
), Q
);
409 Q
= isl_mat_diagonal(Q
, isl_mat_identity(bset
->ctx
, nrest
));
410 C
= isl_mat_alloc(bset
->ctx
, 1 + n_eq
, otype
);
413 isl_int_set_si(C
->row
[0][0], 1);
414 isl_seq_clr(C
->row
[0] + 1, otype
- 1);
415 isl_mat_sub_neg(C
->ctx
, C
->row
+ 1, bset
->eq
+ f_eq
, n_eq
, 0, 0, otype
);
416 H1
= isl_mat_sub_alloc(H
, 0, H
->n_row
, 0, H
->n_row
);
417 H1
= isl_mat_lin_to_aff(H1
);
418 C
= isl_mat_inverse_product(H1
, C
);
423 if (!isl_int_is_one(C
->row
[0][0])) {
428 for (i
= 0; i
< n_eq
; ++i
) {
429 isl_seq_gcd(C
->row
[1 + i
] + 1, otype
- 1, &g
);
430 isl_int_gcd(g
, g
, C
->row
[0][0]);
431 if (!isl_int_is_divisible_by(C
->row
[1 + i
][0], g
))
440 return isl_morph_empty(bset
);
443 C
= isl_mat_normalize(C
);
446 U1
= isl_mat_sub_alloc(U
, 0, U
->n_row
, 0, n_eq
);
447 U1
= isl_mat_lin_to_aff(U1
);
448 U2
= isl_mat_sub_alloc(U
, 0, U
->n_row
, n_eq
, U
->n_row
- n_eq
);
449 U2
= isl_mat_lin_to_aff(U2
);
452 C
= isl_mat_product(U1
, C
);
453 C
= isl_mat_aff_direct_sum(C
, U2
);
454 C
= insert_parameter_rows(C
, otype
- 1);
455 C
= isl_mat_diagonal(C
, isl_mat_identity(bset
->ctx
, nrest
));
457 dim
= isl_space_copy(bset
->dim
);
458 dim
= isl_space_drop_dims(dim
, type
, 0, ntype
);
459 dim
= isl_space_add_dims(dim
, type
, ntype
- n_eq
);
460 ran
= isl_basic_set_universe(dim
);
461 dom
= copy_equalities(bset
, f_eq
, n_eq
);
463 return isl_morph_alloc(dom
, ran
, Q
, C
);
472 /* Construct a parameter compression for "bset".
473 * We basically just call isl_mat_parameter_compression with the right input
474 * and then extend the resulting matrix to include the variables.
476 * The implementation assumes that "bset" does not have any equalities
477 * that only involve the parameters and that isl_basic_set_gauss has
478 * been applied to "bset".
480 * Let the equalities be given as
484 * and let [H 0] be the Hermite Normal Form of A, then
488 * needs to be integer, so we impose that each row is divisible by
491 __isl_give isl_morph
*isl_basic_set_parameter_compression(
492 __isl_keep isl_basic_set
*bset
)
501 isl_basic_set
*dom
, *ran
;
506 if (isl_basic_set_plain_is_empty(bset
))
507 return isl_morph_empty(bset
);
509 return isl_morph_identity(bset
);
512 nparam
= isl_basic_set_dim(bset
, isl_dim_param
);
513 nvar
= isl_basic_set_dim(bset
, isl_dim_set
);
514 n_div
= isl_basic_set_dim(bset
, isl_dim_div
);
516 if (isl_seq_first_non_zero(bset
->eq
[bset
->n_eq
- 1] + 1 + nparam
,
518 isl_die(isl_basic_set_get_ctx(bset
), isl_error_invalid
,
519 "input not allowed to have parameter equalities",
521 if (n_eq
> nvar
+ n_div
)
522 isl_die(isl_basic_set_get_ctx(bset
), isl_error_invalid
,
523 "input not gaussed", return NULL
);
525 d
= isl_vec_alloc(bset
->ctx
, n_eq
);
526 B
= isl_mat_sub_alloc6(bset
->ctx
, bset
->eq
, 0, n_eq
, 0, 1 + nparam
);
527 H
= isl_mat_sub_alloc6(bset
->ctx
, bset
->eq
,
528 0, n_eq
, 1 + nparam
, nvar
+ n_div
);
529 H
= isl_mat_left_hermite(H
, 0, NULL
, NULL
);
530 H
= isl_mat_drop_cols(H
, n_eq
, (nvar
+ n_div
) - n_eq
);
531 H
= isl_mat_lin_to_aff(H
);
532 H
= isl_mat_right_inverse(H
);
535 d
= isl_vec_set(d
, H
->row
[0][0]);
536 H
= isl_mat_drop_rows(H
, 0, 1);
537 H
= isl_mat_drop_cols(H
, 0, 1);
538 B
= isl_mat_product(H
, B
);
539 inv
= isl_mat_parameter_compression(B
, d
);
540 inv
= isl_mat_diagonal(inv
, isl_mat_identity(bset
->ctx
, nvar
));
541 map
= isl_mat_right_inverse(isl_mat_copy(inv
));
543 dom
= isl_basic_set_universe(isl_space_copy(bset
->dim
));
544 ran
= isl_basic_set_universe(isl_space_copy(bset
->dim
));
546 return isl_morph_alloc(dom
, ran
, map
, inv
);
554 /* Add stride constraints to "bset" based on the inverse mapping
555 * that was plugged in. In particular, if morph maps x' to x,
556 * the the constraints of the original input
560 * have been rewritten to
564 * However, this substitution may loose information on the integrality of x',
565 * so we need to impose that
569 * is integral. If inv = B/d, this means that we need to impose that
575 * exists alpha in Z^m: B x = d alpha
578 static __isl_give isl_basic_set
*add_strides(__isl_take isl_basic_set
*bset
,
579 __isl_keep isl_morph
*morph
)
584 if (isl_int_is_one(morph
->inv
->row
[0][0]))
589 for (i
= 0; 1 + i
< morph
->inv
->n_row
; ++i
) {
590 isl_seq_gcd(morph
->inv
->row
[1 + i
], morph
->inv
->n_col
, &gcd
);
591 if (isl_int_is_divisible_by(gcd
, morph
->inv
->row
[0][0]))
593 div
= isl_basic_set_alloc_div(bset
);
596 isl_int_set_si(bset
->div
[div
][0], 0);
597 k
= isl_basic_set_alloc_equality(bset
);
600 isl_seq_cpy(bset
->eq
[k
], morph
->inv
->row
[1 + i
],
602 isl_seq_clr(bset
->eq
[k
] + morph
->inv
->n_col
, bset
->n_div
);
603 isl_int_set(bset
->eq
[k
][morph
->inv
->n_col
+ div
],
604 morph
->inv
->row
[0][0]);
612 isl_basic_set_free(bset
);
616 /* Apply the morphism to the basic set.
617 * We basically just compute the preimage of "bset" under the inverse mapping
618 * in morph, add in stride constraints and intersect with the range
621 __isl_give isl_basic_set
*isl_morph_basic_set(__isl_take isl_morph
*morph
,
622 __isl_take isl_basic_set
*bset
)
624 isl_basic_set
*res
= NULL
;
632 isl_assert(bset
->ctx
, isl_space_is_equal(bset
->dim
, morph
->dom
->dim
),
635 max_stride
= morph
->inv
->n_row
- 1;
636 if (isl_int_is_one(morph
->inv
->row
[0][0]))
638 res
= isl_basic_set_alloc_space(isl_space_copy(morph
->ran
->dim
),
639 bset
->n_div
+ max_stride
, bset
->n_eq
+ max_stride
, bset
->n_ineq
);
641 for (i
= 0; i
< bset
->n_div
; ++i
)
642 if (isl_basic_set_alloc_div(res
) < 0)
645 mat
= isl_mat_sub_alloc6(bset
->ctx
, bset
->eq
, 0, bset
->n_eq
,
646 0, morph
->inv
->n_row
);
647 mat
= isl_mat_product(mat
, isl_mat_copy(morph
->inv
));
650 for (i
= 0; i
< bset
->n_eq
; ++i
) {
651 k
= isl_basic_set_alloc_equality(res
);
654 isl_seq_cpy(res
->eq
[k
], mat
->row
[i
], mat
->n_col
);
655 isl_seq_scale(res
->eq
[k
] + mat
->n_col
, bset
->eq
[i
] + mat
->n_col
,
656 morph
->inv
->row
[0][0], bset
->n_div
);
660 mat
= isl_mat_sub_alloc6(bset
->ctx
, bset
->ineq
, 0, bset
->n_ineq
,
661 0, morph
->inv
->n_row
);
662 mat
= isl_mat_product(mat
, isl_mat_copy(morph
->inv
));
665 for (i
= 0; i
< bset
->n_ineq
; ++i
) {
666 k
= isl_basic_set_alloc_inequality(res
);
669 isl_seq_cpy(res
->ineq
[k
], mat
->row
[i
], mat
->n_col
);
670 isl_seq_scale(res
->ineq
[k
] + mat
->n_col
,
671 bset
->ineq
[i
] + mat
->n_col
,
672 morph
->inv
->row
[0][0], bset
->n_div
);
676 mat
= isl_mat_sub_alloc6(bset
->ctx
, bset
->div
, 0, bset
->n_div
,
677 1, morph
->inv
->n_row
);
678 mat
= isl_mat_product(mat
, isl_mat_copy(morph
->inv
));
681 for (i
= 0; i
< bset
->n_div
; ++i
) {
682 isl_int_mul(res
->div
[i
][0],
683 morph
->inv
->row
[0][0], bset
->div
[i
][0]);
684 isl_seq_cpy(res
->div
[i
] + 1, mat
->row
[i
], mat
->n_col
);
685 isl_seq_scale(res
->div
[i
] + 1 + mat
->n_col
,
686 bset
->div
[i
] + 1 + mat
->n_col
,
687 morph
->inv
->row
[0][0], bset
->n_div
);
691 res
= add_strides(res
, morph
);
693 if (isl_basic_set_is_rational(bset
))
694 res
= isl_basic_set_set_rational(res
);
696 res
= isl_basic_set_simplify(res
);
697 res
= isl_basic_set_finalize(res
);
699 res
= isl_basic_set_intersect(res
, isl_basic_set_copy(morph
->ran
));
701 isl_morph_free(morph
);
702 isl_basic_set_free(bset
);
706 isl_morph_free(morph
);
707 isl_basic_set_free(bset
);
708 isl_basic_set_free(res
);
712 /* Apply the morphism to the set.
714 __isl_give isl_set
*isl_morph_set(__isl_take isl_morph
*morph
,
715 __isl_take isl_set
*set
)
722 isl_assert(set
->ctx
, isl_space_is_equal(set
->dim
, morph
->dom
->dim
), goto error
);
724 set
= isl_set_cow(set
);
728 isl_space_free(set
->dim
);
729 set
->dim
= isl_space_copy(morph
->ran
->dim
);
733 for (i
= 0; i
< set
->n
; ++i
) {
734 set
->p
[i
] = isl_morph_basic_set(isl_morph_copy(morph
), set
->p
[i
]);
739 isl_morph_free(morph
);
741 ISL_F_CLR(set
, ISL_SET_NORMALIZED
);
746 isl_morph_free(morph
);
750 /* Construct a morphism that first does morph2 and then morph1.
752 __isl_give isl_morph
*isl_morph_compose(__isl_take isl_morph
*morph1
,
753 __isl_take isl_morph
*morph2
)
756 isl_basic_set
*dom
, *ran
;
758 if (!morph1
|| !morph2
)
761 map
= isl_mat_product(isl_mat_copy(morph1
->map
), isl_mat_copy(morph2
->map
));
762 inv
= isl_mat_product(isl_mat_copy(morph2
->inv
), isl_mat_copy(morph1
->inv
));
763 dom
= isl_morph_basic_set(isl_morph_inverse(isl_morph_copy(morph2
)),
764 isl_basic_set_copy(morph1
->dom
));
765 dom
= isl_basic_set_intersect(dom
, isl_basic_set_copy(morph2
->dom
));
766 ran
= isl_morph_basic_set(isl_morph_copy(morph1
),
767 isl_basic_set_copy(morph2
->ran
));
768 ran
= isl_basic_set_intersect(ran
, isl_basic_set_copy(morph1
->ran
));
770 isl_morph_free(morph1
);
771 isl_morph_free(morph2
);
773 return isl_morph_alloc(dom
, ran
, map
, inv
);
775 isl_morph_free(morph1
);
776 isl_morph_free(morph2
);
780 __isl_give isl_morph
*isl_morph_inverse(__isl_take isl_morph
*morph
)
785 morph
= isl_morph_cow(morph
);
790 morph
->dom
= morph
->ran
;
794 morph
->map
= morph
->inv
;
800 /* We detect all the equalities first to avoid implicit equalties
801 * being discovered during the computations. In particular,
802 * the compression on the variables could expose additional stride
803 * constraints on the parameters. This would result in existentially
804 * quantified variables after applying the resulting morph, which
805 * in turn could break invariants of the calling functions.
807 __isl_give isl_morph
*isl_basic_set_full_compression(
808 __isl_keep isl_basic_set
*bset
)
810 isl_morph
*morph
, *morph2
;
812 bset
= isl_basic_set_copy(bset
);
813 bset
= isl_basic_set_detect_equalities(bset
);
815 morph
= isl_basic_set_variable_compression(bset
, isl_dim_param
);
816 bset
= isl_morph_basic_set(isl_morph_copy(morph
), bset
);
818 morph2
= isl_basic_set_parameter_compression(bset
);
819 bset
= isl_morph_basic_set(isl_morph_copy(morph2
), bset
);
821 morph
= isl_morph_compose(morph2
, morph
);
823 morph2
= isl_basic_set_variable_compression(bset
, isl_dim_set
);
824 isl_basic_set_free(bset
);
826 morph
= isl_morph_compose(morph2
, morph
);
831 __isl_give isl_vec
*isl_morph_vec(__isl_take isl_morph
*morph
,
832 __isl_take isl_vec
*vec
)
837 vec
= isl_mat_vec_product(isl_mat_copy(morph
->map
), vec
);
839 isl_morph_free(morph
);
842 isl_morph_free(morph
);