4 #include "isl_map_private.h"
6 struct isl_mat
*isl_mat_alloc(struct isl_ctx
*ctx
,
7 unsigned n_row
, unsigned n_col
)
12 mat
= isl_alloc_type(ctx
, struct isl_mat
);
17 mat
->block
= isl_blk_alloc(ctx
, n_row
* n_col
);
18 if (isl_blk_is_error(mat
->block
))
20 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
24 for (i
= 0; i
< n_row
; ++i
)
25 mat
->row
[i
] = mat
->block
.data
+ i
* n_col
;
37 isl_blk_free(ctx
, mat
->block
);
42 struct isl_mat
*isl_mat_extend(struct isl_mat
*mat
,
43 unsigned n_row
, unsigned n_col
)
51 if (mat
->max_col
>= n_col
&& mat
->n_row
>= n_row
) {
52 if (mat
->n_col
< n_col
)
57 if (mat
->max_col
< n_col
) {
58 struct isl_mat
*new_mat
;
60 if (n_row
< mat
->n_row
)
62 new_mat
= isl_mat_alloc(mat
->ctx
, n_row
, n_col
);
65 for (i
= 0; i
< mat
->n_row
; ++i
)
66 isl_seq_cpy(new_mat
->row
[i
], mat
->row
[i
], mat
->n_col
);
71 mat
= isl_mat_cow(mat
);
75 assert(mat
->ref
== 1);
76 old
= mat
->block
.data
;
77 mat
->block
= isl_blk_extend(mat
->ctx
, mat
->block
, n_row
* mat
->max_col
);
78 if (isl_blk_is_error(mat
->block
))
80 mat
->row
= isl_realloc_array(mat
->ctx
, mat
->row
, isl_int
*, n_row
);
84 for (i
= 0; i
< mat
->n_row
; ++i
)
85 mat
->row
[i
] = mat
->block
.data
+ (mat
->row
[i
] - old
);
86 for (i
= mat
->n_row
; i
< n_row
; ++i
)
87 mat
->row
[i
] = mat
->block
.data
+ i
* mat
->max_col
;
89 if (mat
->n_col
< n_col
)
98 struct isl_mat
*isl_mat_sub_alloc(struct isl_ctx
*ctx
, isl_int
**row
,
99 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
104 mat
= isl_alloc_type(ctx
, struct isl_mat
);
107 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
110 for (i
= 0; i
< n_row
; ++i
)
111 mat
->row
[i
] = row
[first_row
+i
] + first_col
;
117 mat
->block
= isl_blk_empty();
118 mat
->flags
= ISL_MAT_BORROWED
;
125 void isl_mat_sub_copy(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
126 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
130 for (i
= 0; i
< n_row
; ++i
)
131 isl_seq_cpy(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
134 void isl_mat_sub_neg(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_neg(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
143 struct isl_mat
*isl_mat_copy(struct isl_mat
*mat
)
152 struct isl_mat
*isl_mat_dup(struct isl_mat
*mat
)
155 struct isl_mat
*mat2
;
159 mat2
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
162 for (i
= 0; i
< mat
->n_row
; ++i
)
163 isl_seq_cpy(mat2
->row
[i
], mat
->row
[i
], mat
->n_col
);
167 struct isl_mat
*isl_mat_cow(struct isl_mat
*mat
)
169 struct isl_mat
*mat2
;
173 if (mat
->ref
== 1 && !ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
176 mat2
= isl_mat_dup(mat
);
181 void isl_mat_free(struct isl_mat
*mat
)
189 if (!ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
190 isl_blk_free(mat
->ctx
, mat
->block
);
191 isl_ctx_deref(mat
->ctx
);
196 struct isl_mat
*isl_mat_identity(struct isl_ctx
*ctx
, unsigned n_row
)
201 mat
= isl_mat_alloc(ctx
, n_row
, n_row
);
204 for (i
= 0; i
< n_row
; ++i
) {
205 isl_seq_clr(mat
->row
[i
], i
);
206 isl_int_set_si(mat
->row
[i
][i
], 1);
207 isl_seq_clr(mat
->row
[i
]+i
+1, n_row
-(i
+1));
213 struct isl_vec
*isl_mat_vec_product(struct isl_mat
*mat
, struct isl_vec
*vec
)
216 struct isl_vec
*prod
;
221 isl_assert(mat
->ctx
, mat
->n_col
== vec
->size
, goto error
);
223 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_row
);
227 for (i
= 0; i
< prod
->size
; ++i
)
228 isl_seq_inner_product(mat
->row
[i
], vec
->el
, vec
->size
,
229 &prod
->block
.data
[i
]);
239 __isl_give isl_vec
*isl_mat_vec_inverse_product(__isl_take isl_mat
*mat
,
240 __isl_take isl_vec
*vec
)
242 struct isl_mat
*vec_mat
;
247 vec_mat
= isl_mat_alloc(vec
->ctx
, vec
->size
, 1);
250 for (i
= 0; i
< vec
->size
; ++i
)
251 isl_int_set(vec_mat
->row
[i
][0], vec
->el
[i
]);
252 vec_mat
= isl_mat_inverse_product(mat
, vec_mat
);
256 vec
= isl_vec_alloc(vec_mat
->ctx
, vec_mat
->n_row
);
258 for (i
= 0; i
< vec
->size
; ++i
)
259 isl_int_set(vec
->el
[i
], vec_mat
->row
[i
][0]);
260 isl_mat_free(vec_mat
);
268 struct isl_vec
*isl_vec_mat_product(struct isl_vec
*vec
, struct isl_mat
*mat
)
271 struct isl_vec
*prod
;
276 isl_assert(mat
->ctx
, mat
->n_row
== vec
->size
, goto error
);
278 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_col
);
282 for (i
= 0; i
< prod
->size
; ++i
) {
283 isl_int_set_si(prod
->el
[i
], 0);
284 for (j
= 0; j
< vec
->size
; ++j
)
285 isl_int_addmul(prod
->el
[i
], vec
->el
[j
], mat
->row
[j
][i
]);
296 struct isl_mat
*isl_mat_aff_direct_sum(struct isl_mat
*left
,
297 struct isl_mat
*right
)
305 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
306 isl_assert(left
->ctx
, left
->n_row
>= 1, goto error
);
307 isl_assert(left
->ctx
, left
->n_col
>= 1, goto error
);
308 isl_assert(left
->ctx
, right
->n_col
>= 1, goto error
);
309 isl_assert(left
->ctx
,
310 isl_seq_first_non_zero(left
->row
[0]+1, left
->n_col
-1) == -1,
312 isl_assert(left
->ctx
,
313 isl_seq_first_non_zero(right
->row
[0]+1, right
->n_col
-1) == -1,
316 sum
= isl_mat_alloc(left
->ctx
, left
->n_row
, left
->n_col
+ right
->n_col
- 1);
319 isl_int_lcm(sum
->row
[0][0], left
->row
[0][0], right
->row
[0][0]);
320 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
321 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
323 isl_seq_clr(sum
->row
[0]+1, sum
->n_col
-1);
324 for (i
= 1; i
< sum
->n_row
; ++i
) {
325 isl_int_mul(sum
->row
[i
][0], left
->row
[0][0], left
->row
[i
][0]);
326 isl_int_addmul(sum
->row
[i
][0],
327 right
->row
[0][0], right
->row
[i
][0]);
328 isl_seq_scale(sum
->row
[i
]+1, left
->row
[i
]+1, left
->row
[0][0],
330 isl_seq_scale(sum
->row
[i
]+left
->n_col
,
331 right
->row
[i
]+1, right
->row
[0][0],
335 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
336 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
346 static void exchange(struct isl_mat
*M
, struct isl_mat
**U
,
347 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
)
350 for (r
= row
; r
< M
->n_row
; ++r
)
351 isl_int_swap(M
->row
[r
][i
], M
->row
[r
][j
]);
353 for (r
= 0; r
< (*U
)->n_row
; ++r
)
354 isl_int_swap((*U
)->row
[r
][i
], (*U
)->row
[r
][j
]);
357 isl_mat_swap_rows(*Q
, i
, j
);
360 static void subtract(struct isl_mat
*M
, struct isl_mat
**U
,
361 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
, isl_int m
)
364 for (r
= row
; r
< M
->n_row
; ++r
)
365 isl_int_submul(M
->row
[r
][j
], m
, M
->row
[r
][i
]);
367 for (r
= 0; r
< (*U
)->n_row
; ++r
)
368 isl_int_submul((*U
)->row
[r
][j
], m
, (*U
)->row
[r
][i
]);
371 for (r
= 0; r
< (*Q
)->n_col
; ++r
)
372 isl_int_addmul((*Q
)->row
[i
][r
], m
, (*Q
)->row
[j
][r
]);
376 static void oppose(struct isl_mat
*M
, struct isl_mat
**U
,
377 struct isl_mat
**Q
, unsigned row
, unsigned col
)
380 for (r
= row
; r
< M
->n_row
; ++r
)
381 isl_int_neg(M
->row
[r
][col
], M
->row
[r
][col
]);
383 for (r
= 0; r
< (*U
)->n_row
; ++r
)
384 isl_int_neg((*U
)->row
[r
][col
], (*U
)->row
[r
][col
]);
387 isl_seq_neg((*Q
)->row
[col
], (*Q
)->row
[col
], (*Q
)->n_col
);
390 /* Given matrix M, compute
395 * with U and Q unimodular matrices and H a matrix in column echelon form
396 * such that on each echelon row the entries in the non-echelon column
397 * are non-negative (if neg == 0) or non-positive (if neg == 1)
398 * and stricly smaller (in absolute value) than the entries in the echelon
400 * If U or Q are NULL, then these matrices are not computed.
402 struct isl_mat
*isl_mat_left_hermite(struct isl_mat
*M
, int neg
,
403 struct isl_mat
**U
, struct isl_mat
**Q
)
418 *U
= isl_mat_identity(M
->ctx
, M
->n_col
);
423 *Q
= isl_mat_identity(M
->ctx
, M
->n_col
);
430 for (row
= 0; row
< M
->n_row
; ++row
) {
432 first
= isl_seq_abs_min_non_zero(M
->row
[row
]+col
, M
->n_col
-col
);
437 exchange(M
, U
, Q
, row
, first
, col
);
438 if (isl_int_is_neg(M
->row
[row
][col
]))
439 oppose(M
, U
, Q
, row
, col
);
441 while ((off
= isl_seq_first_non_zero(M
->row
[row
]+first
,
442 M
->n_col
-first
)) != -1) {
444 isl_int_fdiv_q(c
, M
->row
[row
][first
], M
->row
[row
][col
]);
445 subtract(M
, U
, Q
, row
, col
, first
, c
);
446 if (!isl_int_is_zero(M
->row
[row
][first
]))
447 exchange(M
, U
, Q
, row
, first
, col
);
451 for (i
= 0; i
< col
; ++i
) {
452 if (isl_int_is_zero(M
->row
[row
][i
]))
455 isl_int_cdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
457 isl_int_fdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
458 if (isl_int_is_zero(c
))
460 subtract(M
, U
, Q
, row
, col
, i
, c
);
479 struct isl_mat
*isl_mat_right_kernel(struct isl_mat
*mat
)
482 struct isl_mat
*U
= NULL
;
485 mat
= isl_mat_left_hermite(mat
, 0, &U
, NULL
);
489 for (i
= 0, rank
= 0; rank
< mat
->n_col
; ++rank
) {
490 while (i
< mat
->n_row
&& isl_int_is_zero(mat
->row
[i
][rank
]))
495 K
= isl_mat_alloc(U
->ctx
, U
->n_row
, U
->n_col
- rank
);
498 isl_mat_sub_copy(K
->ctx
, K
->row
, U
->row
, U
->n_row
, 0, rank
, U
->n_col
-rank
);
508 struct isl_mat
*isl_mat_lin_to_aff(struct isl_mat
*mat
)
511 struct isl_mat
*mat2
;
515 mat2
= isl_mat_alloc(mat
->ctx
, 1+mat
->n_row
, 1+mat
->n_col
);
518 isl_int_set_si(mat2
->row
[0][0], 1);
519 isl_seq_clr(mat2
->row
[0]+1, mat
->n_col
);
520 for (i
= 0; i
< mat
->n_row
; ++i
) {
521 isl_int_set_si(mat2
->row
[1+i
][0], 0);
522 isl_seq_cpy(mat2
->row
[1+i
]+1, mat
->row
[i
], mat
->n_col
);
528 static int row_first_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
532 for (i
= 0; i
< n_row
; ++i
)
533 if (!isl_int_is_zero(row
[i
][col
]))
538 static int row_abs_min_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
540 int i
, min
= row_first_non_zero(row
, n_row
, col
);
543 for (i
= min
+ 1; i
< n_row
; ++i
) {
544 if (isl_int_is_zero(row
[i
][col
]))
546 if (isl_int_abs_lt(row
[i
][col
], row
[min
][col
]))
552 static void inv_exchange(struct isl_mat
*left
, struct isl_mat
*right
,
553 unsigned i
, unsigned j
)
555 left
= isl_mat_swap_rows(left
, i
, j
);
556 right
= isl_mat_swap_rows(right
, i
, j
);
559 static void inv_oppose(
560 struct isl_mat
*left
, struct isl_mat
*right
, unsigned row
)
562 isl_seq_neg(left
->row
[row
]+row
, left
->row
[row
]+row
, left
->n_col
-row
);
563 isl_seq_neg(right
->row
[row
], right
->row
[row
], right
->n_col
);
566 static void inv_subtract(struct isl_mat
*left
, struct isl_mat
*right
,
567 unsigned row
, unsigned i
, isl_int m
)
570 isl_seq_combine(left
->row
[i
]+row
,
571 left
->ctx
->one
, left
->row
[i
]+row
,
572 m
, left
->row
[row
]+row
,
574 isl_seq_combine(right
->row
[i
], right
->ctx
->one
, right
->row
[i
],
575 m
, right
->row
[row
], right
->n_col
);
578 /* Compute inv(left)*right
580 struct isl_mat
*isl_mat_inverse_product(struct isl_mat
*left
,
581 struct isl_mat
*right
)
589 isl_assert(left
->ctx
, left
->n_row
== left
->n_col
, goto error
);
590 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
592 if (left
->n_row
== 0) {
597 left
= isl_mat_cow(left
);
598 right
= isl_mat_cow(right
);
604 for (row
= 0; row
< left
->n_row
; ++row
) {
605 int pivot
, first
, i
, off
;
606 pivot
= row_abs_min_non_zero(left
->row
+row
, left
->n_row
-row
, row
);
610 isl_assert(left
->ctx
, pivot
>= 0, goto error
);
614 inv_exchange(left
, right
, pivot
, row
);
615 if (isl_int_is_neg(left
->row
[row
][row
]))
616 inv_oppose(left
, right
, row
);
618 while ((off
= row_first_non_zero(left
->row
+first
,
619 left
->n_row
-first
, row
)) != -1) {
621 isl_int_fdiv_q(a
, left
->row
[first
][row
],
622 left
->row
[row
][row
]);
623 inv_subtract(left
, right
, row
, first
, a
);
624 if (!isl_int_is_zero(left
->row
[first
][row
]))
625 inv_exchange(left
, right
, row
, first
);
629 for (i
= 0; i
< row
; ++i
) {
630 if (isl_int_is_zero(left
->row
[i
][row
]))
632 isl_int_gcd(a
, left
->row
[row
][row
], left
->row
[i
][row
]);
633 isl_int_divexact(b
, left
->row
[i
][row
], a
);
634 isl_int_divexact(a
, left
->row
[row
][row
], a
);
636 isl_seq_combine(left
->row
[i
] + i
,
638 b
, left
->row
[row
] + i
,
640 isl_seq_combine(right
->row
[i
], a
, right
->row
[i
],
641 b
, right
->row
[row
], right
->n_col
);
646 isl_int_set(a
, left
->row
[0][0]);
647 for (row
= 1; row
< left
->n_row
; ++row
)
648 isl_int_lcm(a
, a
, left
->row
[row
][row
]);
649 if (isl_int_is_zero(a
)){
651 isl_assert(left
->ctx
, 0, goto error
);
653 for (row
= 0; row
< left
->n_row
; ++row
) {
654 isl_int_divexact(left
->row
[row
][row
], a
, left
->row
[row
][row
]);
655 if (isl_int_is_one(left
->row
[row
][row
]))
657 isl_seq_scale(right
->row
[row
], right
->row
[row
],
658 left
->row
[row
][row
], right
->n_col
);
670 void isl_mat_col_scale(struct isl_mat
*mat
, unsigned col
, isl_int m
)
674 for (i
= 0; i
< mat
->n_row
; ++i
)
675 isl_int_mul(mat
->row
[i
][col
], mat
->row
[i
][col
], m
);
678 void isl_mat_col_combine(struct isl_mat
*mat
, unsigned dst
,
679 isl_int m1
, unsigned src1
, isl_int m2
, unsigned src2
)
685 for (i
= 0; i
< mat
->n_row
; ++i
) {
686 isl_int_mul(tmp
, m1
, mat
->row
[i
][src1
]);
687 isl_int_addmul(tmp
, m2
, mat
->row
[i
][src2
]);
688 isl_int_set(mat
->row
[i
][dst
], tmp
);
693 struct isl_mat
*isl_mat_right_inverse(struct isl_mat
*mat
)
699 mat
= isl_mat_cow(mat
);
703 inv
= isl_mat_identity(mat
->ctx
, mat
->n_col
);
704 inv
= isl_mat_cow(inv
);
710 for (row
= 0; row
< mat
->n_row
; ++row
) {
711 int pivot
, first
, i
, off
;
712 pivot
= isl_seq_abs_min_non_zero(mat
->row
[row
]+row
, mat
->n_col
-row
);
720 exchange(mat
, &inv
, NULL
, row
, pivot
, row
);
721 if (isl_int_is_neg(mat
->row
[row
][row
]))
722 oppose(mat
, &inv
, NULL
, row
, row
);
724 while ((off
= isl_seq_first_non_zero(mat
->row
[row
]+first
,
725 mat
->n_col
-first
)) != -1) {
727 isl_int_fdiv_q(a
, mat
->row
[row
][first
],
729 subtract(mat
, &inv
, NULL
, row
, row
, first
, a
);
730 if (!isl_int_is_zero(mat
->row
[row
][first
]))
731 exchange(mat
, &inv
, NULL
, row
, row
, first
);
735 for (i
= 0; i
< row
; ++i
) {
736 if (isl_int_is_zero(mat
->row
[row
][i
]))
738 isl_int_gcd(a
, mat
->row
[row
][row
], mat
->row
[row
][i
]);
739 isl_int_divexact(b
, mat
->row
[row
][i
], a
);
740 isl_int_divexact(a
, mat
->row
[row
][row
], a
);
742 isl_mat_col_combine(mat
, i
, a
, i
, b
, row
);
743 isl_mat_col_combine(inv
, i
, a
, i
, b
, row
);
748 isl_int_set(a
, mat
->row
[0][0]);
749 for (row
= 1; row
< mat
->n_row
; ++row
)
750 isl_int_lcm(a
, a
, mat
->row
[row
][row
]);
751 if (isl_int_is_zero(a
)){
755 for (row
= 0; row
< mat
->n_row
; ++row
) {
756 isl_int_divexact(mat
->row
[row
][row
], a
, mat
->row
[row
][row
]);
757 if (isl_int_is_one(mat
->row
[row
][row
]))
759 isl_mat_col_scale(inv
, row
, mat
->row
[row
][row
]);
771 struct isl_mat
*isl_mat_transpose(struct isl_mat
*mat
)
773 struct isl_mat
*transpose
= NULL
;
776 if (mat
->n_col
== mat
->n_row
) {
777 mat
= isl_mat_cow(mat
);
780 for (i
= 0; i
< mat
->n_row
; ++i
)
781 for (j
= i
+ 1; j
< mat
->n_col
; ++j
)
782 isl_int_swap(mat
->row
[i
][j
], mat
->row
[j
][i
]);
785 transpose
= isl_mat_alloc(mat
->ctx
, mat
->n_col
, mat
->n_row
);
788 for (i
= 0; i
< mat
->n_row
; ++i
)
789 for (j
= 0; j
< mat
->n_col
; ++j
)
790 isl_int_set(transpose
->row
[j
][i
], mat
->row
[i
][j
]);
798 struct isl_mat
*isl_mat_swap_cols(struct isl_mat
*mat
, unsigned i
, unsigned j
)
802 mat
= isl_mat_cow(mat
);
805 isl_assert(mat
->ctx
, i
< mat
->n_col
, goto error
);
806 isl_assert(mat
->ctx
, j
< mat
->n_col
, goto error
);
808 for (r
= 0; r
< mat
->n_row
; ++r
)
809 isl_int_swap(mat
->row
[r
][i
], mat
->row
[r
][j
]);
816 struct isl_mat
*isl_mat_swap_rows(struct isl_mat
*mat
, unsigned i
, unsigned j
)
822 mat
= isl_mat_cow(mat
);
826 mat
->row
[i
] = mat
->row
[j
];
831 struct isl_mat
*isl_mat_product(struct isl_mat
*left
, struct isl_mat
*right
)
834 struct isl_mat
*prod
;
838 isl_assert(left
->ctx
, left
->n_col
== right
->n_row
, goto error
);
839 prod
= isl_mat_alloc(left
->ctx
, left
->n_row
, right
->n_col
);
842 if (left
->n_col
== 0) {
843 for (i
= 0; i
< prod
->n_row
; ++i
)
844 isl_seq_clr(prod
->row
[i
], prod
->n_col
);
847 for (i
= 0; i
< prod
->n_row
; ++i
) {
848 for (j
= 0; j
< prod
->n_col
; ++j
) {
849 isl_int_mul(prod
->row
[i
][j
],
850 left
->row
[i
][0], right
->row
[0][j
]);
851 for (k
= 1; k
< left
->n_col
; ++k
)
852 isl_int_addmul(prod
->row
[i
][j
],
853 left
->row
[i
][k
], right
->row
[k
][j
]);
865 /* Replace the variables x in the rows q by x' given by x = M x',
866 * with M the matrix mat.
868 * If the number of new variables is greater than the original
869 * number of variables, then the rows q have already been
870 * preextended. If the new number is smaller, then the coefficients
871 * of the divs, which are not changed, need to be shifted down.
872 * The row q may be the equalities, the inequalities or the
873 * div expressions. In the latter case, has_div is true and
874 * we need to take into account the extra denominator column.
876 static int preimage(struct isl_ctx
*ctx
, isl_int
**q
, unsigned n
,
877 unsigned n_div
, int has_div
, struct isl_mat
*mat
)
883 if (mat
->n_col
>= mat
->n_row
)
886 e
= mat
->n_row
- mat
->n_col
;
888 for (i
= 0; i
< n
; ++i
)
889 isl_int_mul(q
[i
][0], q
[i
][0], mat
->row
[0][0]);
890 t
= isl_mat_sub_alloc(mat
->ctx
, q
, 0, n
, has_div
, mat
->n_row
);
891 t
= isl_mat_product(t
, mat
);
894 for (i
= 0; i
< n
; ++i
) {
895 isl_seq_swp_or_cpy(q
[i
] + has_div
, t
->row
[i
], t
->n_col
);
896 isl_seq_cpy(q
[i
] + has_div
+ t
->n_col
,
897 q
[i
] + has_div
+ t
->n_col
+ e
, n_div
);
898 isl_seq_clr(q
[i
] + has_div
+ t
->n_col
+ n_div
, e
);
904 /* Replace the variables x in bset by x' given by x = M x', with
907 * If there are fewer variables x' then there are x, then we perform
908 * the transformation in place, which that, in principle,
909 * this frees up some extra variables as the number
910 * of columns remains constant, but we would have to extend
911 * the div array too as the number of rows in this array is assumed
912 * to be equal to extra.
914 struct isl_basic_set
*isl_basic_set_preimage(struct isl_basic_set
*bset
,
923 bset
= isl_basic_set_cow(bset
);
927 isl_assert(ctx
, bset
->dim
->nparam
== 0, goto error
);
928 isl_assert(ctx
, 1+bset
->dim
->n_out
== mat
->n_row
, goto error
);
930 if (mat
->n_col
> mat
->n_row
)
931 bset
= isl_basic_set_extend(bset
, 0, mat
->n_col
-1, 0,
933 else if (mat
->n_col
< mat
->n_row
) {
934 bset
->dim
= isl_dim_cow(bset
->dim
);
937 bset
->dim
->n_out
-= mat
->n_row
- mat
->n_col
;
940 if (preimage(ctx
, bset
->eq
, bset
->n_eq
, bset
->n_div
, 0,
941 isl_mat_copy(mat
)) < 0)
944 if (preimage(ctx
, bset
->ineq
, bset
->n_ineq
, bset
->n_div
, 0,
945 isl_mat_copy(mat
)) < 0)
948 if (preimage(ctx
, bset
->div
, bset
->n_div
, bset
->n_div
, 1, mat
) < 0)
951 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
952 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_REDUNDANT
);
953 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED
);
954 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED_DIVS
);
955 ISL_F_CLR(bset
, ISL_BASIC_SET_ALL_EQUALITIES
);
957 bset
= isl_basic_set_simplify(bset
);
958 bset
= isl_basic_set_finalize(bset
);
964 isl_basic_set_free(bset
);
968 struct isl_set
*isl_set_preimage(struct isl_set
*set
, struct isl_mat
*mat
)
973 set
= isl_set_cow(set
);
978 for (i
= 0; i
< set
->n
; ++i
) {
979 set
->p
[i
] = isl_basic_set_preimage(set
->p
[i
],
984 if (mat
->n_col
!= mat
->n_row
) {
985 set
->dim
= isl_dim_cow(set
->dim
);
988 set
->dim
->n_out
+= mat
->n_col
;
989 set
->dim
->n_out
-= mat
->n_row
;
992 ISL_F_CLR(set
, ISL_SET_NORMALIZED
);
1000 void isl_mat_dump(struct isl_mat
*mat
, FILE *out
, int indent
)
1005 fprintf(out
, "%*snull mat\n", indent
, "");
1009 if (mat
->n_row
== 0)
1010 fprintf(out
, "%*s[]\n", indent
, "");
1012 for (i
= 0; i
< mat
->n_row
; ++i
) {
1014 fprintf(out
, "%*s[[", indent
, "");
1016 fprintf(out
, "%*s[", indent
+1, "");
1017 for (j
= 0; j
< mat
->n_col
; ++j
) {
1020 isl_int_print(out
, mat
->row
[i
][j
], 0);
1022 if (i
== mat
->n_row
-1)
1023 fprintf(out
, "]]\n");
1025 fprintf(out
, "]\n");
1029 struct isl_mat
*isl_mat_drop_cols(struct isl_mat
*mat
, unsigned col
, unsigned n
)
1033 mat
= isl_mat_cow(mat
);
1037 if (col
!= mat
->n_col
-n
) {
1038 for (r
= 0; r
< mat
->n_row
; ++r
)
1039 isl_seq_cpy(mat
->row
[r
]+col
, mat
->row
[r
]+col
+n
,
1040 mat
->n_col
- col
- n
);
1046 struct isl_mat
*isl_mat_drop_rows(struct isl_mat
*mat
, unsigned row
, unsigned n
)
1050 mat
= isl_mat_cow(mat
);
1054 for (r
= row
; r
+n
< mat
->n_row
; ++r
)
1055 mat
->row
[r
] = mat
->row
[r
+n
];
1061 void isl_mat_col_submul(struct isl_mat
*mat
,
1062 int dst_col
, isl_int f
, int src_col
)
1066 for (i
= 0; i
< mat
->n_row
; ++i
)
1067 isl_int_submul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1070 void isl_mat_col_mul(struct isl_mat
*mat
, int dst_col
, isl_int f
, int src_col
)
1074 for (i
= 0; i
< mat
->n_row
; ++i
)
1075 isl_int_mul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1078 struct isl_mat
*isl_mat_unimodular_complete(struct isl_mat
*M
, int row
)
1081 struct isl_mat
*H
= NULL
, *Q
= NULL
;
1086 isl_assert(M
->ctx
, M
->n_row
== M
->n_col
, goto error
);
1088 H
= isl_mat_left_hermite(isl_mat_copy(M
), 0, NULL
, &Q
);
1089 M
->n_row
= M
->n_col
;
1092 for (r
= 0; r
< row
; ++r
)
1093 isl_assert(M
->ctx
, isl_int_is_one(H
->row
[r
][r
]), goto error
);
1094 for (r
= row
; r
< M
->n_row
; ++r
)
1095 isl_seq_cpy(M
->row
[r
], Q
->row
[r
], M
->n_col
);
1106 __isl_give isl_mat
*isl_mat_concat(__isl_take isl_mat
*top
,
1107 __isl_take isl_mat
*bot
)
1109 struct isl_mat
*mat
;
1114 isl_assert(top
->ctx
, top
->n_col
== bot
->n_col
, goto error
);
1115 if (top
->n_row
== 0) {
1119 if (bot
->n_row
== 0) {
1124 mat
= isl_mat_alloc(top
->ctx
, top
->n_row
+ bot
->n_row
, top
->n_col
);
1127 isl_mat_sub_copy(mat
->ctx
, mat
->row
, top
->row
, top
->n_row
,
1129 isl_mat_sub_copy(mat
->ctx
, mat
->row
+ top
->n_row
, bot
->row
, bot
->n_row
,
1140 int isl_mat_is_equal(__isl_keep isl_mat
*mat1
, __isl_keep isl_mat
*mat2
)
1147 if (mat1
->n_row
!= mat2
->n_row
)
1150 if (mat1
->n_col
!= mat2
->n_col
)
1153 for (i
= 0; i
< mat1
->n_row
; ++i
)
1154 if (!isl_seq_eq(mat1
->row
[i
], mat2
->row
[i
], mat1
->n_col
))