2 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
13 #include "isl_map_private.h"
15 struct isl_mat
*isl_mat_alloc(struct isl_ctx
*ctx
,
16 unsigned n_row
, unsigned n_col
)
21 mat
= isl_alloc_type(ctx
, struct isl_mat
);
26 mat
->block
= isl_blk_alloc(ctx
, n_row
* n_col
);
27 if (isl_blk_is_error(mat
->block
))
29 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
33 for (i
= 0; i
< n_row
; ++i
)
34 mat
->row
[i
] = mat
->block
.data
+ i
* n_col
;
46 isl_blk_free(ctx
, mat
->block
);
51 struct isl_mat
*isl_mat_extend(struct isl_mat
*mat
,
52 unsigned n_row
, unsigned n_col
)
60 if (mat
->max_col
>= n_col
&& mat
->n_row
>= n_row
) {
61 if (mat
->n_col
< n_col
)
66 if (mat
->max_col
< n_col
) {
67 struct isl_mat
*new_mat
;
69 if (n_row
< mat
->n_row
)
71 new_mat
= isl_mat_alloc(mat
->ctx
, n_row
, n_col
);
74 for (i
= 0; i
< mat
->n_row
; ++i
)
75 isl_seq_cpy(new_mat
->row
[i
], mat
->row
[i
], mat
->n_col
);
80 mat
= isl_mat_cow(mat
);
84 assert(mat
->ref
== 1);
85 old
= mat
->block
.data
;
86 mat
->block
= isl_blk_extend(mat
->ctx
, mat
->block
, n_row
* mat
->max_col
);
87 if (isl_blk_is_error(mat
->block
))
89 mat
->row
= isl_realloc_array(mat
->ctx
, mat
->row
, isl_int
*, n_row
);
93 for (i
= 0; i
< mat
->n_row
; ++i
)
94 mat
->row
[i
] = mat
->block
.data
+ (mat
->row
[i
] - old
);
95 for (i
= mat
->n_row
; i
< n_row
; ++i
)
96 mat
->row
[i
] = mat
->block
.data
+ i
* mat
->max_col
;
98 if (mat
->n_col
< n_col
)
107 struct isl_mat
*isl_mat_sub_alloc(struct isl_ctx
*ctx
, isl_int
**row
,
108 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
113 mat
= isl_alloc_type(ctx
, struct isl_mat
);
116 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
119 for (i
= 0; i
< n_row
; ++i
)
120 mat
->row
[i
] = row
[first_row
+i
] + first_col
;
126 mat
->block
= isl_blk_empty();
127 mat
->flags
= ISL_MAT_BORROWED
;
134 void isl_mat_sub_copy(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
135 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
139 for (i
= 0; i
< n_row
; ++i
)
140 isl_seq_cpy(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
143 void isl_mat_sub_neg(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
144 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
148 for (i
= 0; i
< n_row
; ++i
)
149 isl_seq_neg(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
152 struct isl_mat
*isl_mat_copy(struct isl_mat
*mat
)
161 struct isl_mat
*isl_mat_dup(struct isl_mat
*mat
)
164 struct isl_mat
*mat2
;
168 mat2
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
171 for (i
= 0; i
< mat
->n_row
; ++i
)
172 isl_seq_cpy(mat2
->row
[i
], mat
->row
[i
], mat
->n_col
);
176 struct isl_mat
*isl_mat_cow(struct isl_mat
*mat
)
178 struct isl_mat
*mat2
;
182 if (mat
->ref
== 1 && !ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
185 mat2
= isl_mat_dup(mat
);
190 void isl_mat_free(struct isl_mat
*mat
)
198 if (!ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
199 isl_blk_free(mat
->ctx
, mat
->block
);
200 isl_ctx_deref(mat
->ctx
);
205 struct isl_mat
*isl_mat_identity(struct isl_ctx
*ctx
, unsigned n_row
)
210 mat
= isl_mat_alloc(ctx
, n_row
, n_row
);
213 for (i
= 0; i
< n_row
; ++i
) {
214 isl_seq_clr(mat
->row
[i
], i
);
215 isl_int_set_si(mat
->row
[i
][i
], 1);
216 isl_seq_clr(mat
->row
[i
]+i
+1, n_row
-(i
+1));
222 struct isl_vec
*isl_mat_vec_product(struct isl_mat
*mat
, struct isl_vec
*vec
)
225 struct isl_vec
*prod
;
230 isl_assert(mat
->ctx
, mat
->n_col
== vec
->size
, goto error
);
232 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_row
);
236 for (i
= 0; i
< prod
->size
; ++i
)
237 isl_seq_inner_product(mat
->row
[i
], vec
->el
, vec
->size
,
238 &prod
->block
.data
[i
]);
248 __isl_give isl_vec
*isl_mat_vec_inverse_product(__isl_take isl_mat
*mat
,
249 __isl_take isl_vec
*vec
)
251 struct isl_mat
*vec_mat
;
256 vec_mat
= isl_mat_alloc(vec
->ctx
, vec
->size
, 1);
259 for (i
= 0; i
< vec
->size
; ++i
)
260 isl_int_set(vec_mat
->row
[i
][0], vec
->el
[i
]);
261 vec_mat
= isl_mat_inverse_product(mat
, vec_mat
);
265 vec
= isl_vec_alloc(vec_mat
->ctx
, vec_mat
->n_row
);
267 for (i
= 0; i
< vec
->size
; ++i
)
268 isl_int_set(vec
->el
[i
], vec_mat
->row
[i
][0]);
269 isl_mat_free(vec_mat
);
277 struct isl_vec
*isl_vec_mat_product(struct isl_vec
*vec
, struct isl_mat
*mat
)
280 struct isl_vec
*prod
;
285 isl_assert(mat
->ctx
, mat
->n_row
== vec
->size
, goto error
);
287 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_col
);
291 for (i
= 0; i
< prod
->size
; ++i
) {
292 isl_int_set_si(prod
->el
[i
], 0);
293 for (j
= 0; j
< vec
->size
; ++j
)
294 isl_int_addmul(prod
->el
[i
], vec
->el
[j
], mat
->row
[j
][i
]);
305 struct isl_mat
*isl_mat_aff_direct_sum(struct isl_mat
*left
,
306 struct isl_mat
*right
)
314 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
315 isl_assert(left
->ctx
, left
->n_row
>= 1, goto error
);
316 isl_assert(left
->ctx
, left
->n_col
>= 1, goto error
);
317 isl_assert(left
->ctx
, right
->n_col
>= 1, goto error
);
318 isl_assert(left
->ctx
,
319 isl_seq_first_non_zero(left
->row
[0]+1, left
->n_col
-1) == -1,
321 isl_assert(left
->ctx
,
322 isl_seq_first_non_zero(right
->row
[0]+1, right
->n_col
-1) == -1,
325 sum
= isl_mat_alloc(left
->ctx
, left
->n_row
, left
->n_col
+ right
->n_col
- 1);
328 isl_int_lcm(sum
->row
[0][0], left
->row
[0][0], right
->row
[0][0]);
329 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
330 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
332 isl_seq_clr(sum
->row
[0]+1, sum
->n_col
-1);
333 for (i
= 1; i
< sum
->n_row
; ++i
) {
334 isl_int_mul(sum
->row
[i
][0], left
->row
[0][0], left
->row
[i
][0]);
335 isl_int_addmul(sum
->row
[i
][0],
336 right
->row
[0][0], right
->row
[i
][0]);
337 isl_seq_scale(sum
->row
[i
]+1, left
->row
[i
]+1, left
->row
[0][0],
339 isl_seq_scale(sum
->row
[i
]+left
->n_col
,
340 right
->row
[i
]+1, right
->row
[0][0],
344 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
345 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
355 static void exchange(struct isl_mat
*M
, struct isl_mat
**U
,
356 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
)
359 for (r
= row
; r
< M
->n_row
; ++r
)
360 isl_int_swap(M
->row
[r
][i
], M
->row
[r
][j
]);
362 for (r
= 0; r
< (*U
)->n_row
; ++r
)
363 isl_int_swap((*U
)->row
[r
][i
], (*U
)->row
[r
][j
]);
366 isl_mat_swap_rows(*Q
, i
, j
);
369 static void subtract(struct isl_mat
*M
, struct isl_mat
**U
,
370 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
, isl_int m
)
373 for (r
= row
; r
< M
->n_row
; ++r
)
374 isl_int_submul(M
->row
[r
][j
], m
, M
->row
[r
][i
]);
376 for (r
= 0; r
< (*U
)->n_row
; ++r
)
377 isl_int_submul((*U
)->row
[r
][j
], m
, (*U
)->row
[r
][i
]);
380 for (r
= 0; r
< (*Q
)->n_col
; ++r
)
381 isl_int_addmul((*Q
)->row
[i
][r
], m
, (*Q
)->row
[j
][r
]);
385 static void oppose(struct isl_mat
*M
, struct isl_mat
**U
,
386 struct isl_mat
**Q
, unsigned row
, unsigned col
)
389 for (r
= row
; r
< M
->n_row
; ++r
)
390 isl_int_neg(M
->row
[r
][col
], M
->row
[r
][col
]);
392 for (r
= 0; r
< (*U
)->n_row
; ++r
)
393 isl_int_neg((*U
)->row
[r
][col
], (*U
)->row
[r
][col
]);
396 isl_seq_neg((*Q
)->row
[col
], (*Q
)->row
[col
], (*Q
)->n_col
);
399 /* Given matrix M, compute
404 * with U and Q unimodular matrices and H a matrix in column echelon form
405 * such that on each echelon row the entries in the non-echelon column
406 * are non-negative (if neg == 0) or non-positive (if neg == 1)
407 * and stricly smaller (in absolute value) than the entries in the echelon
409 * If U or Q are NULL, then these matrices are not computed.
411 struct isl_mat
*isl_mat_left_hermite(struct isl_mat
*M
, int neg
,
412 struct isl_mat
**U
, struct isl_mat
**Q
)
427 *U
= isl_mat_identity(M
->ctx
, M
->n_col
);
432 *Q
= isl_mat_identity(M
->ctx
, M
->n_col
);
439 for (row
= 0; row
< M
->n_row
; ++row
) {
441 first
= isl_seq_abs_min_non_zero(M
->row
[row
]+col
, M
->n_col
-col
);
446 exchange(M
, U
, Q
, row
, first
, col
);
447 if (isl_int_is_neg(M
->row
[row
][col
]))
448 oppose(M
, U
, Q
, row
, col
);
450 while ((off
= isl_seq_first_non_zero(M
->row
[row
]+first
,
451 M
->n_col
-first
)) != -1) {
453 isl_int_fdiv_q(c
, M
->row
[row
][first
], M
->row
[row
][col
]);
454 subtract(M
, U
, Q
, row
, col
, first
, c
);
455 if (!isl_int_is_zero(M
->row
[row
][first
]))
456 exchange(M
, U
, Q
, row
, first
, col
);
460 for (i
= 0; i
< col
; ++i
) {
461 if (isl_int_is_zero(M
->row
[row
][i
]))
464 isl_int_cdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
466 isl_int_fdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
467 if (isl_int_is_zero(c
))
469 subtract(M
, U
, Q
, row
, col
, i
, c
);
488 struct isl_mat
*isl_mat_right_kernel(struct isl_mat
*mat
)
491 struct isl_mat
*U
= NULL
;
494 mat
= isl_mat_left_hermite(mat
, 0, &U
, NULL
);
498 for (i
= 0, rank
= 0; rank
< mat
->n_col
; ++rank
) {
499 while (i
< mat
->n_row
&& isl_int_is_zero(mat
->row
[i
][rank
]))
504 K
= isl_mat_alloc(U
->ctx
, U
->n_row
, U
->n_col
- rank
);
507 isl_mat_sub_copy(K
->ctx
, K
->row
, U
->row
, U
->n_row
, 0, rank
, U
->n_col
-rank
);
517 struct isl_mat
*isl_mat_lin_to_aff(struct isl_mat
*mat
)
520 struct isl_mat
*mat2
;
524 mat2
= isl_mat_alloc(mat
->ctx
, 1+mat
->n_row
, 1+mat
->n_col
);
527 isl_int_set_si(mat2
->row
[0][0], 1);
528 isl_seq_clr(mat2
->row
[0]+1, mat
->n_col
);
529 for (i
= 0; i
< mat
->n_row
; ++i
) {
530 isl_int_set_si(mat2
->row
[1+i
][0], 0);
531 isl_seq_cpy(mat2
->row
[1+i
]+1, mat
->row
[i
], mat
->n_col
);
537 static int row_first_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
541 for (i
= 0; i
< n_row
; ++i
)
542 if (!isl_int_is_zero(row
[i
][col
]))
547 static int row_abs_min_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
549 int i
, min
= row_first_non_zero(row
, n_row
, col
);
552 for (i
= min
+ 1; i
< n_row
; ++i
) {
553 if (isl_int_is_zero(row
[i
][col
]))
555 if (isl_int_abs_lt(row
[i
][col
], row
[min
][col
]))
561 static void inv_exchange(struct isl_mat
*left
, struct isl_mat
*right
,
562 unsigned i
, unsigned j
)
564 left
= isl_mat_swap_rows(left
, i
, j
);
565 right
= isl_mat_swap_rows(right
, i
, j
);
568 static void inv_oppose(
569 struct isl_mat
*left
, struct isl_mat
*right
, unsigned row
)
571 isl_seq_neg(left
->row
[row
]+row
, left
->row
[row
]+row
, left
->n_col
-row
);
572 isl_seq_neg(right
->row
[row
], right
->row
[row
], right
->n_col
);
575 static void inv_subtract(struct isl_mat
*left
, struct isl_mat
*right
,
576 unsigned row
, unsigned i
, isl_int m
)
579 isl_seq_combine(left
->row
[i
]+row
,
580 left
->ctx
->one
, left
->row
[i
]+row
,
581 m
, left
->row
[row
]+row
,
583 isl_seq_combine(right
->row
[i
], right
->ctx
->one
, right
->row
[i
],
584 m
, right
->row
[row
], right
->n_col
);
587 /* Compute inv(left)*right
589 struct isl_mat
*isl_mat_inverse_product(struct isl_mat
*left
,
590 struct isl_mat
*right
)
598 isl_assert(left
->ctx
, left
->n_row
== left
->n_col
, goto error
);
599 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
601 if (left
->n_row
== 0) {
606 left
= isl_mat_cow(left
);
607 right
= isl_mat_cow(right
);
613 for (row
= 0; row
< left
->n_row
; ++row
) {
614 int pivot
, first
, i
, off
;
615 pivot
= row_abs_min_non_zero(left
->row
+row
, left
->n_row
-row
, row
);
619 isl_assert(left
->ctx
, pivot
>= 0, goto error
);
623 inv_exchange(left
, right
, pivot
, row
);
624 if (isl_int_is_neg(left
->row
[row
][row
]))
625 inv_oppose(left
, right
, row
);
627 while ((off
= row_first_non_zero(left
->row
+first
,
628 left
->n_row
-first
, row
)) != -1) {
630 isl_int_fdiv_q(a
, left
->row
[first
][row
],
631 left
->row
[row
][row
]);
632 inv_subtract(left
, right
, row
, first
, a
);
633 if (!isl_int_is_zero(left
->row
[first
][row
]))
634 inv_exchange(left
, right
, row
, first
);
638 for (i
= 0; i
< row
; ++i
) {
639 if (isl_int_is_zero(left
->row
[i
][row
]))
641 isl_int_gcd(a
, left
->row
[row
][row
], left
->row
[i
][row
]);
642 isl_int_divexact(b
, left
->row
[i
][row
], a
);
643 isl_int_divexact(a
, left
->row
[row
][row
], a
);
645 isl_seq_combine(left
->row
[i
] + i
,
647 b
, left
->row
[row
] + i
,
649 isl_seq_combine(right
->row
[i
], a
, right
->row
[i
],
650 b
, right
->row
[row
], right
->n_col
);
655 isl_int_set(a
, left
->row
[0][0]);
656 for (row
= 1; row
< left
->n_row
; ++row
)
657 isl_int_lcm(a
, a
, left
->row
[row
][row
]);
658 if (isl_int_is_zero(a
)){
660 isl_assert(left
->ctx
, 0, goto error
);
662 for (row
= 0; row
< left
->n_row
; ++row
) {
663 isl_int_divexact(left
->row
[row
][row
], a
, left
->row
[row
][row
]);
664 if (isl_int_is_one(left
->row
[row
][row
]))
666 isl_seq_scale(right
->row
[row
], right
->row
[row
],
667 left
->row
[row
][row
], right
->n_col
);
679 void isl_mat_col_scale(struct isl_mat
*mat
, unsigned col
, isl_int m
)
683 for (i
= 0; i
< mat
->n_row
; ++i
)
684 isl_int_mul(mat
->row
[i
][col
], mat
->row
[i
][col
], m
);
687 void isl_mat_col_combine(struct isl_mat
*mat
, unsigned dst
,
688 isl_int m1
, unsigned src1
, isl_int m2
, unsigned src2
)
694 for (i
= 0; i
< mat
->n_row
; ++i
) {
695 isl_int_mul(tmp
, m1
, mat
->row
[i
][src1
]);
696 isl_int_addmul(tmp
, m2
, mat
->row
[i
][src2
]);
697 isl_int_set(mat
->row
[i
][dst
], tmp
);
702 struct isl_mat
*isl_mat_right_inverse(struct isl_mat
*mat
)
708 mat
= isl_mat_cow(mat
);
712 inv
= isl_mat_identity(mat
->ctx
, mat
->n_col
);
713 inv
= isl_mat_cow(inv
);
719 for (row
= 0; row
< mat
->n_row
; ++row
) {
720 int pivot
, first
, i
, off
;
721 pivot
= isl_seq_abs_min_non_zero(mat
->row
[row
]+row
, mat
->n_col
-row
);
725 isl_assert(mat
->ctx
, pivot
>= 0, goto error
);
729 exchange(mat
, &inv
, NULL
, row
, pivot
, row
);
730 if (isl_int_is_neg(mat
->row
[row
][row
]))
731 oppose(mat
, &inv
, NULL
, row
, row
);
733 while ((off
= isl_seq_first_non_zero(mat
->row
[row
]+first
,
734 mat
->n_col
-first
)) != -1) {
736 isl_int_fdiv_q(a
, mat
->row
[row
][first
],
738 subtract(mat
, &inv
, NULL
, row
, row
, first
, a
);
739 if (!isl_int_is_zero(mat
->row
[row
][first
]))
740 exchange(mat
, &inv
, NULL
, row
, row
, first
);
744 for (i
= 0; i
< row
; ++i
) {
745 if (isl_int_is_zero(mat
->row
[row
][i
]))
747 isl_int_gcd(a
, mat
->row
[row
][row
], mat
->row
[row
][i
]);
748 isl_int_divexact(b
, mat
->row
[row
][i
], a
);
749 isl_int_divexact(a
, mat
->row
[row
][row
], a
);
751 isl_mat_col_combine(mat
, i
, a
, i
, b
, row
);
752 isl_mat_col_combine(inv
, i
, a
, i
, b
, row
);
757 isl_int_set(a
, mat
->row
[0][0]);
758 for (row
= 1; row
< mat
->n_row
; ++row
)
759 isl_int_lcm(a
, a
, mat
->row
[row
][row
]);
760 if (isl_int_is_zero(a
)){
764 for (row
= 0; row
< mat
->n_row
; ++row
) {
765 isl_int_divexact(mat
->row
[row
][row
], a
, mat
->row
[row
][row
]);
766 if (isl_int_is_one(mat
->row
[row
][row
]))
768 isl_mat_col_scale(inv
, row
, mat
->row
[row
][row
]);
780 struct isl_mat
*isl_mat_transpose(struct isl_mat
*mat
)
782 struct isl_mat
*transpose
= NULL
;
785 if (mat
->n_col
== mat
->n_row
) {
786 mat
= isl_mat_cow(mat
);
789 for (i
= 0; i
< mat
->n_row
; ++i
)
790 for (j
= i
+ 1; j
< mat
->n_col
; ++j
)
791 isl_int_swap(mat
->row
[i
][j
], mat
->row
[j
][i
]);
794 transpose
= isl_mat_alloc(mat
->ctx
, mat
->n_col
, mat
->n_row
);
797 for (i
= 0; i
< mat
->n_row
; ++i
)
798 for (j
= 0; j
< mat
->n_col
; ++j
)
799 isl_int_set(transpose
->row
[j
][i
], mat
->row
[i
][j
]);
807 struct isl_mat
*isl_mat_swap_cols(struct isl_mat
*mat
, unsigned i
, unsigned j
)
811 mat
= isl_mat_cow(mat
);
814 isl_assert(mat
->ctx
, i
< mat
->n_col
, goto error
);
815 isl_assert(mat
->ctx
, j
< mat
->n_col
, goto error
);
817 for (r
= 0; r
< mat
->n_row
; ++r
)
818 isl_int_swap(mat
->row
[r
][i
], mat
->row
[r
][j
]);
825 struct isl_mat
*isl_mat_swap_rows(struct isl_mat
*mat
, unsigned i
, unsigned j
)
831 mat
= isl_mat_cow(mat
);
835 mat
->row
[i
] = mat
->row
[j
];
840 struct isl_mat
*isl_mat_product(struct isl_mat
*left
, struct isl_mat
*right
)
843 struct isl_mat
*prod
;
847 isl_assert(left
->ctx
, left
->n_col
== right
->n_row
, goto error
);
848 prod
= isl_mat_alloc(left
->ctx
, left
->n_row
, right
->n_col
);
851 if (left
->n_col
== 0) {
852 for (i
= 0; i
< prod
->n_row
; ++i
)
853 isl_seq_clr(prod
->row
[i
], prod
->n_col
);
856 for (i
= 0; i
< prod
->n_row
; ++i
) {
857 for (j
= 0; j
< prod
->n_col
; ++j
) {
858 isl_int_mul(prod
->row
[i
][j
],
859 left
->row
[i
][0], right
->row
[0][j
]);
860 for (k
= 1; k
< left
->n_col
; ++k
)
861 isl_int_addmul(prod
->row
[i
][j
],
862 left
->row
[i
][k
], right
->row
[k
][j
]);
874 /* Replace the variables x in the rows q by x' given by x = M x',
875 * with M the matrix mat.
877 * If the number of new variables is greater than the original
878 * number of variables, then the rows q have already been
879 * preextended. If the new number is smaller, then the coefficients
880 * of the divs, which are not changed, need to be shifted down.
881 * The row q may be the equalities, the inequalities or the
882 * div expressions. In the latter case, has_div is true and
883 * we need to take into account the extra denominator column.
885 static int preimage(struct isl_ctx
*ctx
, isl_int
**q
, unsigned n
,
886 unsigned n_div
, int has_div
, struct isl_mat
*mat
)
892 if (mat
->n_col
>= mat
->n_row
)
895 e
= mat
->n_row
- mat
->n_col
;
897 for (i
= 0; i
< n
; ++i
)
898 isl_int_mul(q
[i
][0], q
[i
][0], mat
->row
[0][0]);
899 t
= isl_mat_sub_alloc(mat
->ctx
, q
, 0, n
, has_div
, mat
->n_row
);
900 t
= isl_mat_product(t
, mat
);
903 for (i
= 0; i
< n
; ++i
) {
904 isl_seq_swp_or_cpy(q
[i
] + has_div
, t
->row
[i
], t
->n_col
);
905 isl_seq_cpy(q
[i
] + has_div
+ t
->n_col
,
906 q
[i
] + has_div
+ t
->n_col
+ e
, n_div
);
907 isl_seq_clr(q
[i
] + has_div
+ t
->n_col
+ n_div
, e
);
913 /* Replace the variables x in bset by x' given by x = M x', with
916 * If there are fewer variables x' then there are x, then we perform
917 * the transformation in place, which that, in principle,
918 * this frees up some extra variables as the number
919 * of columns remains constant, but we would have to extend
920 * the div array too as the number of rows in this array is assumed
921 * to be equal to extra.
923 struct isl_basic_set
*isl_basic_set_preimage(struct isl_basic_set
*bset
,
932 bset
= isl_basic_set_cow(bset
);
936 isl_assert(ctx
, bset
->dim
->nparam
== 0, goto error
);
937 isl_assert(ctx
, 1+bset
->dim
->n_out
== mat
->n_row
, goto error
);
939 if (mat
->n_col
> mat
->n_row
)
940 bset
= isl_basic_set_extend(bset
, 0, mat
->n_col
-1, 0,
942 else if (mat
->n_col
< mat
->n_row
) {
943 bset
->dim
= isl_dim_cow(bset
->dim
);
946 bset
->dim
->n_out
-= mat
->n_row
- mat
->n_col
;
949 if (preimage(ctx
, bset
->eq
, bset
->n_eq
, bset
->n_div
, 0,
950 isl_mat_copy(mat
)) < 0)
953 if (preimage(ctx
, bset
->ineq
, bset
->n_ineq
, bset
->n_div
, 0,
954 isl_mat_copy(mat
)) < 0)
957 if (preimage(ctx
, bset
->div
, bset
->n_div
, bset
->n_div
, 1, mat
) < 0)
960 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
961 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_REDUNDANT
);
962 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED
);
963 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED_DIVS
);
964 ISL_F_CLR(bset
, ISL_BASIC_SET_ALL_EQUALITIES
);
966 bset
= isl_basic_set_simplify(bset
);
967 bset
= isl_basic_set_finalize(bset
);
973 isl_basic_set_free(bset
);
977 struct isl_set
*isl_set_preimage(struct isl_set
*set
, struct isl_mat
*mat
)
982 set
= isl_set_cow(set
);
987 for (i
= 0; i
< set
->n
; ++i
) {
988 set
->p
[i
] = isl_basic_set_preimage(set
->p
[i
],
993 if (mat
->n_col
!= mat
->n_row
) {
994 set
->dim
= isl_dim_cow(set
->dim
);
997 set
->dim
->n_out
+= mat
->n_col
;
998 set
->dim
->n_out
-= mat
->n_row
;
1001 ISL_F_CLR(set
, ISL_SET_NORMALIZED
);
1009 void isl_mat_dump(struct isl_mat
*mat
, FILE *out
, int indent
)
1014 fprintf(out
, "%*snull mat\n", indent
, "");
1018 if (mat
->n_row
== 0)
1019 fprintf(out
, "%*s[]\n", indent
, "");
1021 for (i
= 0; i
< mat
->n_row
; ++i
) {
1023 fprintf(out
, "%*s[[", indent
, "");
1025 fprintf(out
, "%*s[", indent
+1, "");
1026 for (j
= 0; j
< mat
->n_col
; ++j
) {
1029 isl_int_print(out
, mat
->row
[i
][j
], 0);
1031 if (i
== mat
->n_row
-1)
1032 fprintf(out
, "]]\n");
1034 fprintf(out
, "]\n");
1038 struct isl_mat
*isl_mat_drop_cols(struct isl_mat
*mat
, unsigned col
, unsigned n
)
1042 mat
= isl_mat_cow(mat
);
1046 if (col
!= mat
->n_col
-n
) {
1047 for (r
= 0; r
< mat
->n_row
; ++r
)
1048 isl_seq_cpy(mat
->row
[r
]+col
, mat
->row
[r
]+col
+n
,
1049 mat
->n_col
- col
- n
);
1055 struct isl_mat
*isl_mat_drop_rows(struct isl_mat
*mat
, unsigned row
, unsigned n
)
1059 mat
= isl_mat_cow(mat
);
1063 for (r
= row
; r
+n
< mat
->n_row
; ++r
)
1064 mat
->row
[r
] = mat
->row
[r
+n
];
1070 __isl_give isl_mat
*isl_mat_insert_cols(__isl_take isl_mat
*mat
,
1071 unsigned col
, unsigned n
)
1080 ext
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
+ n
);
1084 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
, 0, 0, col
);
1085 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
,
1086 col
+ n
, col
, mat
->n_col
- col
);
1095 void isl_mat_col_submul(struct isl_mat
*mat
,
1096 int dst_col
, isl_int f
, int src_col
)
1100 for (i
= 0; i
< mat
->n_row
; ++i
)
1101 isl_int_submul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1104 void isl_mat_col_mul(struct isl_mat
*mat
, int dst_col
, isl_int f
, int src_col
)
1108 for (i
= 0; i
< mat
->n_row
; ++i
)
1109 isl_int_mul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1112 struct isl_mat
*isl_mat_unimodular_complete(struct isl_mat
*M
, int row
)
1115 struct isl_mat
*H
= NULL
, *Q
= NULL
;
1120 isl_assert(M
->ctx
, M
->n_row
== M
->n_col
, goto error
);
1122 H
= isl_mat_left_hermite(isl_mat_copy(M
), 0, NULL
, &Q
);
1123 M
->n_row
= M
->n_col
;
1126 for (r
= 0; r
< row
; ++r
)
1127 isl_assert(M
->ctx
, isl_int_is_one(H
->row
[r
][r
]), goto error
);
1128 for (r
= row
; r
< M
->n_row
; ++r
)
1129 isl_seq_cpy(M
->row
[r
], Q
->row
[r
], M
->n_col
);
1140 __isl_give isl_mat
*isl_mat_concat(__isl_take isl_mat
*top
,
1141 __isl_take isl_mat
*bot
)
1143 struct isl_mat
*mat
;
1148 isl_assert(top
->ctx
, top
->n_col
== bot
->n_col
, goto error
);
1149 if (top
->n_row
== 0) {
1153 if (bot
->n_row
== 0) {
1158 mat
= isl_mat_alloc(top
->ctx
, top
->n_row
+ bot
->n_row
, top
->n_col
);
1161 isl_mat_sub_copy(mat
->ctx
, mat
->row
, top
->row
, top
->n_row
,
1163 isl_mat_sub_copy(mat
->ctx
, mat
->row
+ top
->n_row
, bot
->row
, bot
->n_row
,
1174 int isl_mat_is_equal(__isl_keep isl_mat
*mat1
, __isl_keep isl_mat
*mat2
)
1181 if (mat1
->n_row
!= mat2
->n_row
)
1184 if (mat1
->n_col
!= mat2
->n_col
)
1187 for (i
= 0; i
< mat1
->n_row
; ++i
)
1188 if (!isl_seq_eq(mat1
->row
[i
], mat2
->row
[i
], mat1
->n_col
))
1194 __isl_give isl_mat
*isl_mat_from_row_vec(__isl_take isl_vec
*vec
)
1196 struct isl_mat
*mat
;
1200 mat
= isl_mat_alloc(vec
->ctx
, 1, vec
->size
);
1204 isl_seq_cpy(mat
->row
[0], vec
->el
, vec
->size
);
1213 __isl_give isl_mat
*isl_mat_vec_concat(__isl_take isl_mat
*top
,
1214 __isl_take isl_vec
*bot
)
1216 return isl_mat_concat(top
, isl_mat_from_row_vec(bot
));