2 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Use of this software is governed by the MIT license
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
10 #include <isl_ctx_private.h>
11 #include <isl_map_private.h>
12 #include <isl/space.h>
14 #include <isl_mat_private.h>
15 #include <isl_space_private.h>
16 #include <isl_val_private.h>
18 isl_ctx
*isl_mat_get_ctx(__isl_keep isl_mat
*mat
)
20 return mat
? mat
->ctx
: NULL
;
23 struct isl_mat
*isl_mat_alloc(struct isl_ctx
*ctx
,
24 unsigned n_row
, unsigned n_col
)
29 mat
= isl_alloc_type(ctx
, struct isl_mat
);
34 mat
->block
= isl_blk_alloc(ctx
, n_row
* n_col
);
35 if (isl_blk_is_error(mat
->block
))
37 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
38 if (n_row
&& !mat
->row
)
41 for (i
= 0; i
< n_row
; ++i
)
42 mat
->row
[i
] = mat
->block
.data
+ i
* n_col
;
54 isl_blk_free(ctx
, mat
->block
);
59 struct isl_mat
*isl_mat_extend(struct isl_mat
*mat
,
60 unsigned n_row
, unsigned n_col
)
69 if (mat
->max_col
>= n_col
&& mat
->n_row
>= n_row
) {
70 if (mat
->n_col
< n_col
)
75 if (mat
->max_col
< n_col
) {
76 struct isl_mat
*new_mat
;
78 if (n_row
< mat
->n_row
)
80 new_mat
= isl_mat_alloc(mat
->ctx
, n_row
, n_col
);
83 for (i
= 0; i
< mat
->n_row
; ++i
)
84 isl_seq_cpy(new_mat
->row
[i
], mat
->row
[i
], mat
->n_col
);
89 mat
= isl_mat_cow(mat
);
93 old
= mat
->block
.data
;
94 mat
->block
= isl_blk_extend(mat
->ctx
, mat
->block
, n_row
* mat
->max_col
);
95 if (isl_blk_is_error(mat
->block
))
97 row
= isl_realloc_array(mat
->ctx
, mat
->row
, isl_int
*, n_row
);
102 for (i
= 0; i
< mat
->n_row
; ++i
)
103 mat
->row
[i
] = mat
->block
.data
+ (mat
->row
[i
] - old
);
104 for (i
= mat
->n_row
; i
< n_row
; ++i
)
105 mat
->row
[i
] = mat
->block
.data
+ i
* mat
->max_col
;
107 if (mat
->n_col
< n_col
)
116 __isl_give isl_mat
*isl_mat_sub_alloc6(isl_ctx
*ctx
, isl_int
**row
,
117 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
122 mat
= isl_alloc_type(ctx
, struct isl_mat
);
125 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
126 if (n_row
&& !mat
->row
)
128 for (i
= 0; i
< n_row
; ++i
)
129 mat
->row
[i
] = row
[first_row
+i
] + first_col
;
135 mat
->block
= isl_blk_empty();
136 mat
->flags
= ISL_MAT_BORROWED
;
143 __isl_give isl_mat
*isl_mat_sub_alloc(__isl_keep isl_mat
*mat
,
144 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
148 return isl_mat_sub_alloc6(mat
->ctx
, mat
->row
, first_row
, n_row
,
152 void isl_mat_sub_copy(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
153 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
157 for (i
= 0; i
< n_row
; ++i
)
158 isl_seq_cpy(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
161 void isl_mat_sub_neg(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
162 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
166 for (i
= 0; i
< n_row
; ++i
)
167 isl_seq_neg(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
170 struct isl_mat
*isl_mat_copy(struct isl_mat
*mat
)
179 struct isl_mat
*isl_mat_dup(struct isl_mat
*mat
)
182 struct isl_mat
*mat2
;
186 mat2
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
189 for (i
= 0; i
< mat
->n_row
; ++i
)
190 isl_seq_cpy(mat2
->row
[i
], mat
->row
[i
], mat
->n_col
);
194 struct isl_mat
*isl_mat_cow(struct isl_mat
*mat
)
196 struct isl_mat
*mat2
;
200 if (mat
->ref
== 1 && !ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
203 mat2
= isl_mat_dup(mat
);
208 void *isl_mat_free(struct isl_mat
*mat
)
216 if (!ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
217 isl_blk_free(mat
->ctx
, mat
->block
);
218 isl_ctx_deref(mat
->ctx
);
225 int isl_mat_rows(__isl_keep isl_mat
*mat
)
227 return mat
? mat
->n_row
: -1;
230 int isl_mat_cols(__isl_keep isl_mat
*mat
)
232 return mat
? mat
->n_col
: -1;
235 int isl_mat_get_element(__isl_keep isl_mat
*mat
, int row
, int col
, isl_int
*v
)
239 if (row
< 0 || row
>= mat
->n_row
)
240 isl_die(mat
->ctx
, isl_error_invalid
, "row out of range",
242 if (col
< 0 || col
>= mat
->n_col
)
243 isl_die(mat
->ctx
, isl_error_invalid
, "column out of range",
245 isl_int_set(*v
, mat
->row
[row
][col
]);
249 /* Extract the element at row "row", oolumn "col" of "mat".
251 __isl_give isl_val
*isl_mat_get_element_val(__isl_keep isl_mat
*mat
,
258 ctx
= isl_mat_get_ctx(mat
);
259 if (row
< 0 || row
>= mat
->n_row
)
260 isl_die(ctx
, isl_error_invalid
, "row out of range",
262 if (col
< 0 || col
>= mat
->n_col
)
263 isl_die(ctx
, isl_error_invalid
, "column out of range",
265 return isl_val_int_from_isl_int(ctx
, mat
->row
[row
][col
]);
268 __isl_give isl_mat
*isl_mat_set_element(__isl_take isl_mat
*mat
,
269 int row
, int col
, isl_int v
)
271 mat
= isl_mat_cow(mat
);
274 if (row
< 0 || row
>= mat
->n_row
)
275 isl_die(mat
->ctx
, isl_error_invalid
, "row out of range",
277 if (col
< 0 || col
>= mat
->n_col
)
278 isl_die(mat
->ctx
, isl_error_invalid
, "column out of range",
280 isl_int_set(mat
->row
[row
][col
], v
);
287 __isl_give isl_mat
*isl_mat_set_element_si(__isl_take isl_mat
*mat
,
288 int row
, int col
, int v
)
290 mat
= isl_mat_cow(mat
);
293 if (row
< 0 || row
>= mat
->n_row
)
294 isl_die(mat
->ctx
, isl_error_invalid
, "row out of range",
296 if (col
< 0 || col
>= mat
->n_col
)
297 isl_die(mat
->ctx
, isl_error_invalid
, "column out of range",
299 isl_int_set_si(mat
->row
[row
][col
], v
);
306 /* Replace the element at row "row", column "col" of "mat" by "v".
308 __isl_give isl_mat
*isl_mat_set_element_val(__isl_take isl_mat
*mat
,
309 int row
, int col
, __isl_take isl_val
*v
)
312 return isl_mat_free(mat
);
313 if (!isl_val_is_int(v
))
314 isl_die(isl_val_get_ctx(v
), isl_error_invalid
,
315 "expecting integer value", goto error
);
316 mat
= isl_mat_set_element(mat
, row
, col
, v
->n
);
321 return isl_mat_free(mat
);
324 __isl_give isl_mat
*isl_mat_diag(isl_ctx
*ctx
, unsigned n_row
, isl_int d
)
329 mat
= isl_mat_alloc(ctx
, n_row
, n_row
);
332 for (i
= 0; i
< n_row
; ++i
) {
333 isl_seq_clr(mat
->row
[i
], i
);
334 isl_int_set(mat
->row
[i
][i
], d
);
335 isl_seq_clr(mat
->row
[i
]+i
+1, n_row
-(i
+1));
341 __isl_give isl_mat
*isl_mat_identity(isl_ctx
*ctx
, unsigned n_row
)
345 return isl_mat_diag(ctx
, n_row
, ctx
->one
);
348 struct isl_vec
*isl_mat_vec_product(struct isl_mat
*mat
, struct isl_vec
*vec
)
351 struct isl_vec
*prod
;
356 isl_assert(mat
->ctx
, mat
->n_col
== vec
->size
, goto error
);
358 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_row
);
362 for (i
= 0; i
< prod
->size
; ++i
)
363 isl_seq_inner_product(mat
->row
[i
], vec
->el
, vec
->size
,
364 &prod
->block
.data
[i
]);
374 __isl_give isl_vec
*isl_mat_vec_inverse_product(__isl_take isl_mat
*mat
,
375 __isl_take isl_vec
*vec
)
377 struct isl_mat
*vec_mat
;
382 vec_mat
= isl_mat_alloc(vec
->ctx
, vec
->size
, 1);
385 for (i
= 0; i
< vec
->size
; ++i
)
386 isl_int_set(vec_mat
->row
[i
][0], vec
->el
[i
]);
387 vec_mat
= isl_mat_inverse_product(mat
, vec_mat
);
391 vec
= isl_vec_alloc(vec_mat
->ctx
, vec_mat
->n_row
);
393 for (i
= 0; i
< vec
->size
; ++i
)
394 isl_int_set(vec
->el
[i
], vec_mat
->row
[i
][0]);
395 isl_mat_free(vec_mat
);
403 struct isl_vec
*isl_vec_mat_product(struct isl_vec
*vec
, struct isl_mat
*mat
)
406 struct isl_vec
*prod
;
411 isl_assert(mat
->ctx
, mat
->n_row
== vec
->size
, goto error
);
413 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_col
);
417 for (i
= 0; i
< prod
->size
; ++i
) {
418 isl_int_set_si(prod
->el
[i
], 0);
419 for (j
= 0; j
< vec
->size
; ++j
)
420 isl_int_addmul(prod
->el
[i
], vec
->el
[j
], mat
->row
[j
][i
]);
431 struct isl_mat
*isl_mat_aff_direct_sum(struct isl_mat
*left
,
432 struct isl_mat
*right
)
440 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
441 isl_assert(left
->ctx
, left
->n_row
>= 1, goto error
);
442 isl_assert(left
->ctx
, left
->n_col
>= 1, goto error
);
443 isl_assert(left
->ctx
, right
->n_col
>= 1, goto error
);
444 isl_assert(left
->ctx
,
445 isl_seq_first_non_zero(left
->row
[0]+1, left
->n_col
-1) == -1,
447 isl_assert(left
->ctx
,
448 isl_seq_first_non_zero(right
->row
[0]+1, right
->n_col
-1) == -1,
451 sum
= isl_mat_alloc(left
->ctx
, left
->n_row
, left
->n_col
+ right
->n_col
- 1);
454 isl_int_lcm(sum
->row
[0][0], left
->row
[0][0], right
->row
[0][0]);
455 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
456 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
458 isl_seq_clr(sum
->row
[0]+1, sum
->n_col
-1);
459 for (i
= 1; i
< sum
->n_row
; ++i
) {
460 isl_int_mul(sum
->row
[i
][0], left
->row
[0][0], left
->row
[i
][0]);
461 isl_int_addmul(sum
->row
[i
][0],
462 right
->row
[0][0], right
->row
[i
][0]);
463 isl_seq_scale(sum
->row
[i
]+1, left
->row
[i
]+1, left
->row
[0][0],
465 isl_seq_scale(sum
->row
[i
]+left
->n_col
,
466 right
->row
[i
]+1, right
->row
[0][0],
470 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
471 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
481 static void exchange(struct isl_mat
*M
, struct isl_mat
**U
,
482 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
)
485 for (r
= row
; r
< M
->n_row
; ++r
)
486 isl_int_swap(M
->row
[r
][i
], M
->row
[r
][j
]);
488 for (r
= 0; r
< (*U
)->n_row
; ++r
)
489 isl_int_swap((*U
)->row
[r
][i
], (*U
)->row
[r
][j
]);
492 isl_mat_swap_rows(*Q
, i
, j
);
495 static void subtract(struct isl_mat
*M
, struct isl_mat
**U
,
496 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
, isl_int m
)
499 for (r
= row
; r
< M
->n_row
; ++r
)
500 isl_int_submul(M
->row
[r
][j
], m
, M
->row
[r
][i
]);
502 for (r
= 0; r
< (*U
)->n_row
; ++r
)
503 isl_int_submul((*U
)->row
[r
][j
], m
, (*U
)->row
[r
][i
]);
506 for (r
= 0; r
< (*Q
)->n_col
; ++r
)
507 isl_int_addmul((*Q
)->row
[i
][r
], m
, (*Q
)->row
[j
][r
]);
511 static void oppose(struct isl_mat
*M
, struct isl_mat
**U
,
512 struct isl_mat
**Q
, unsigned row
, unsigned col
)
515 for (r
= row
; r
< M
->n_row
; ++r
)
516 isl_int_neg(M
->row
[r
][col
], M
->row
[r
][col
]);
518 for (r
= 0; r
< (*U
)->n_row
; ++r
)
519 isl_int_neg((*U
)->row
[r
][col
], (*U
)->row
[r
][col
]);
522 isl_seq_neg((*Q
)->row
[col
], (*Q
)->row
[col
], (*Q
)->n_col
);
525 /* Given matrix M, compute
530 * with U and Q unimodular matrices and H a matrix in column echelon form
531 * such that on each echelon row the entries in the non-echelon column
532 * are non-negative (if neg == 0) or non-positive (if neg == 1)
533 * and stricly smaller (in absolute value) than the entries in the echelon
535 * If U or Q are NULL, then these matrices are not computed.
537 struct isl_mat
*isl_mat_left_hermite(struct isl_mat
*M
, int neg
,
538 struct isl_mat
**U
, struct isl_mat
**Q
)
553 *U
= isl_mat_identity(M
->ctx
, M
->n_col
);
558 *Q
= isl_mat_identity(M
->ctx
, M
->n_col
);
565 for (row
= 0; row
< M
->n_row
; ++row
) {
567 first
= isl_seq_abs_min_non_zero(M
->row
[row
]+col
, M
->n_col
-col
);
572 exchange(M
, U
, Q
, row
, first
, col
);
573 if (isl_int_is_neg(M
->row
[row
][col
]))
574 oppose(M
, U
, Q
, row
, col
);
576 while ((off
= isl_seq_first_non_zero(M
->row
[row
]+first
,
577 M
->n_col
-first
)) != -1) {
579 isl_int_fdiv_q(c
, M
->row
[row
][first
], M
->row
[row
][col
]);
580 subtract(M
, U
, Q
, row
, col
, first
, c
);
581 if (!isl_int_is_zero(M
->row
[row
][first
]))
582 exchange(M
, U
, Q
, row
, first
, col
);
586 for (i
= 0; i
< col
; ++i
) {
587 if (isl_int_is_zero(M
->row
[row
][i
]))
590 isl_int_cdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
592 isl_int_fdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
593 if (isl_int_is_zero(c
))
595 subtract(M
, U
, Q
, row
, col
, i
, c
);
615 struct isl_mat
*isl_mat_right_kernel(struct isl_mat
*mat
)
618 struct isl_mat
*U
= NULL
;
621 mat
= isl_mat_left_hermite(mat
, 0, &U
, NULL
);
625 for (i
= 0, rank
= 0; rank
< mat
->n_col
; ++rank
) {
626 while (i
< mat
->n_row
&& isl_int_is_zero(mat
->row
[i
][rank
]))
631 K
= isl_mat_alloc(U
->ctx
, U
->n_row
, U
->n_col
- rank
);
634 isl_mat_sub_copy(K
->ctx
, K
->row
, U
->row
, U
->n_row
, 0, rank
, U
->n_col
-rank
);
644 struct isl_mat
*isl_mat_lin_to_aff(struct isl_mat
*mat
)
647 struct isl_mat
*mat2
;
651 mat2
= isl_mat_alloc(mat
->ctx
, 1+mat
->n_row
, 1+mat
->n_col
);
654 isl_int_set_si(mat2
->row
[0][0], 1);
655 isl_seq_clr(mat2
->row
[0]+1, mat
->n_col
);
656 for (i
= 0; i
< mat
->n_row
; ++i
) {
657 isl_int_set_si(mat2
->row
[1+i
][0], 0);
658 isl_seq_cpy(mat2
->row
[1+i
]+1, mat
->row
[i
], mat
->n_col
);
667 /* Given two matrices M1 and M2, return the block matrix
672 __isl_give isl_mat
*isl_mat_diagonal(__isl_take isl_mat
*mat1
,
673 __isl_take isl_mat
*mat2
)
681 mat
= isl_mat_alloc(mat1
->ctx
, mat1
->n_row
+ mat2
->n_row
,
682 mat1
->n_col
+ mat2
->n_col
);
685 for (i
= 0; i
< mat1
->n_row
; ++i
) {
686 isl_seq_cpy(mat
->row
[i
], mat1
->row
[i
], mat1
->n_col
);
687 isl_seq_clr(mat
->row
[i
] + mat1
->n_col
, mat2
->n_col
);
689 for (i
= 0; i
< mat2
->n_row
; ++i
) {
690 isl_seq_clr(mat
->row
[mat1
->n_row
+ i
], mat1
->n_col
);
691 isl_seq_cpy(mat
->row
[mat1
->n_row
+ i
] + mat1
->n_col
,
692 mat2
->row
[i
], mat2
->n_col
);
703 static int row_first_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
707 for (i
= 0; i
< n_row
; ++i
)
708 if (!isl_int_is_zero(row
[i
][col
]))
713 static int row_abs_min_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
715 int i
, min
= row_first_non_zero(row
, n_row
, col
);
718 for (i
= min
+ 1; i
< n_row
; ++i
) {
719 if (isl_int_is_zero(row
[i
][col
]))
721 if (isl_int_abs_lt(row
[i
][col
], row
[min
][col
]))
727 static void inv_exchange(struct isl_mat
*left
, struct isl_mat
*right
,
728 unsigned i
, unsigned j
)
730 left
= isl_mat_swap_rows(left
, i
, j
);
731 right
= isl_mat_swap_rows(right
, i
, j
);
734 static void inv_oppose(
735 struct isl_mat
*left
, struct isl_mat
*right
, unsigned row
)
737 isl_seq_neg(left
->row
[row
]+row
, left
->row
[row
]+row
, left
->n_col
-row
);
738 isl_seq_neg(right
->row
[row
], right
->row
[row
], right
->n_col
);
741 static void inv_subtract(struct isl_mat
*left
, struct isl_mat
*right
,
742 unsigned row
, unsigned i
, isl_int m
)
745 isl_seq_combine(left
->row
[i
]+row
,
746 left
->ctx
->one
, left
->row
[i
]+row
,
747 m
, left
->row
[row
]+row
,
749 isl_seq_combine(right
->row
[i
], right
->ctx
->one
, right
->row
[i
],
750 m
, right
->row
[row
], right
->n_col
);
753 /* Compute inv(left)*right
755 struct isl_mat
*isl_mat_inverse_product(struct isl_mat
*left
,
756 struct isl_mat
*right
)
764 isl_assert(left
->ctx
, left
->n_row
== left
->n_col
, goto error
);
765 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
767 if (left
->n_row
== 0) {
772 left
= isl_mat_cow(left
);
773 right
= isl_mat_cow(right
);
779 for (row
= 0; row
< left
->n_row
; ++row
) {
780 int pivot
, first
, i
, off
;
781 pivot
= row_abs_min_non_zero(left
->row
+row
, left
->n_row
-row
, row
);
785 isl_assert(left
->ctx
, pivot
>= 0, goto error
);
789 inv_exchange(left
, right
, pivot
, row
);
790 if (isl_int_is_neg(left
->row
[row
][row
]))
791 inv_oppose(left
, right
, row
);
793 while ((off
= row_first_non_zero(left
->row
+first
,
794 left
->n_row
-first
, row
)) != -1) {
796 isl_int_fdiv_q(a
, left
->row
[first
][row
],
797 left
->row
[row
][row
]);
798 inv_subtract(left
, right
, row
, first
, a
);
799 if (!isl_int_is_zero(left
->row
[first
][row
]))
800 inv_exchange(left
, right
, row
, first
);
804 for (i
= 0; i
< row
; ++i
) {
805 if (isl_int_is_zero(left
->row
[i
][row
]))
807 isl_int_gcd(a
, left
->row
[row
][row
], left
->row
[i
][row
]);
808 isl_int_divexact(b
, left
->row
[i
][row
], a
);
809 isl_int_divexact(a
, left
->row
[row
][row
], a
);
811 isl_seq_combine(left
->row
[i
] + i
,
813 b
, left
->row
[row
] + i
,
815 isl_seq_combine(right
->row
[i
], a
, right
->row
[i
],
816 b
, right
->row
[row
], right
->n_col
);
821 isl_int_set(a
, left
->row
[0][0]);
822 for (row
= 1; row
< left
->n_row
; ++row
)
823 isl_int_lcm(a
, a
, left
->row
[row
][row
]);
824 if (isl_int_is_zero(a
)){
826 isl_assert(left
->ctx
, 0, goto error
);
828 for (row
= 0; row
< left
->n_row
; ++row
) {
829 isl_int_divexact(left
->row
[row
][row
], a
, left
->row
[row
][row
]);
830 if (isl_int_is_one(left
->row
[row
][row
]))
832 isl_seq_scale(right
->row
[row
], right
->row
[row
],
833 left
->row
[row
][row
], right
->n_col
);
845 void isl_mat_col_scale(struct isl_mat
*mat
, unsigned col
, isl_int m
)
849 for (i
= 0; i
< mat
->n_row
; ++i
)
850 isl_int_mul(mat
->row
[i
][col
], mat
->row
[i
][col
], m
);
853 void isl_mat_col_combine(struct isl_mat
*mat
, unsigned dst
,
854 isl_int m1
, unsigned src1
, isl_int m2
, unsigned src2
)
860 for (i
= 0; i
< mat
->n_row
; ++i
) {
861 isl_int_mul(tmp
, m1
, mat
->row
[i
][src1
]);
862 isl_int_addmul(tmp
, m2
, mat
->row
[i
][src2
]);
863 isl_int_set(mat
->row
[i
][dst
], tmp
);
868 struct isl_mat
*isl_mat_right_inverse(struct isl_mat
*mat
)
874 mat
= isl_mat_cow(mat
);
878 inv
= isl_mat_identity(mat
->ctx
, mat
->n_col
);
879 inv
= isl_mat_cow(inv
);
885 for (row
= 0; row
< mat
->n_row
; ++row
) {
886 int pivot
, first
, i
, off
;
887 pivot
= isl_seq_abs_min_non_zero(mat
->row
[row
]+row
, mat
->n_col
-row
);
891 isl_assert(mat
->ctx
, pivot
>= 0, goto error
);
895 exchange(mat
, &inv
, NULL
, row
, pivot
, row
);
896 if (isl_int_is_neg(mat
->row
[row
][row
]))
897 oppose(mat
, &inv
, NULL
, row
, row
);
899 while ((off
= isl_seq_first_non_zero(mat
->row
[row
]+first
,
900 mat
->n_col
-first
)) != -1) {
902 isl_int_fdiv_q(a
, mat
->row
[row
][first
],
904 subtract(mat
, &inv
, NULL
, row
, row
, first
, a
);
905 if (!isl_int_is_zero(mat
->row
[row
][first
]))
906 exchange(mat
, &inv
, NULL
, row
, row
, first
);
910 for (i
= 0; i
< row
; ++i
) {
911 if (isl_int_is_zero(mat
->row
[row
][i
]))
913 isl_int_gcd(a
, mat
->row
[row
][row
], mat
->row
[row
][i
]);
914 isl_int_divexact(b
, mat
->row
[row
][i
], a
);
915 isl_int_divexact(a
, mat
->row
[row
][row
], a
);
917 isl_mat_col_combine(mat
, i
, a
, i
, b
, row
);
918 isl_mat_col_combine(inv
, i
, a
, i
, b
, row
);
923 isl_int_set(a
, mat
->row
[0][0]);
924 for (row
= 1; row
< mat
->n_row
; ++row
)
925 isl_int_lcm(a
, a
, mat
->row
[row
][row
]);
926 if (isl_int_is_zero(a
)){
930 for (row
= 0; row
< mat
->n_row
; ++row
) {
931 isl_int_divexact(mat
->row
[row
][row
], a
, mat
->row
[row
][row
]);
932 if (isl_int_is_one(mat
->row
[row
][row
]))
934 isl_mat_col_scale(inv
, row
, mat
->row
[row
][row
]);
947 struct isl_mat
*isl_mat_transpose(struct isl_mat
*mat
)
949 struct isl_mat
*transpose
= NULL
;
955 if (mat
->n_col
== mat
->n_row
) {
956 mat
= isl_mat_cow(mat
);
959 for (i
= 0; i
< mat
->n_row
; ++i
)
960 for (j
= i
+ 1; j
< mat
->n_col
; ++j
)
961 isl_int_swap(mat
->row
[i
][j
], mat
->row
[j
][i
]);
964 transpose
= isl_mat_alloc(mat
->ctx
, mat
->n_col
, mat
->n_row
);
967 for (i
= 0; i
< mat
->n_row
; ++i
)
968 for (j
= 0; j
< mat
->n_col
; ++j
)
969 isl_int_set(transpose
->row
[j
][i
], mat
->row
[i
][j
]);
977 struct isl_mat
*isl_mat_swap_cols(struct isl_mat
*mat
, unsigned i
, unsigned j
)
981 mat
= isl_mat_cow(mat
);
984 isl_assert(mat
->ctx
, i
< mat
->n_col
, goto error
);
985 isl_assert(mat
->ctx
, j
< mat
->n_col
, goto error
);
987 for (r
= 0; r
< mat
->n_row
; ++r
)
988 isl_int_swap(mat
->row
[r
][i
], mat
->row
[r
][j
]);
995 struct isl_mat
*isl_mat_swap_rows(struct isl_mat
*mat
, unsigned i
, unsigned j
)
1001 mat
= isl_mat_cow(mat
);
1005 mat
->row
[i
] = mat
->row
[j
];
1010 __isl_give isl_mat
*isl_mat_product(__isl_take isl_mat
*left
,
1011 __isl_take isl_mat
*right
)
1014 struct isl_mat
*prod
;
1016 if (!left
|| !right
)
1018 isl_assert(left
->ctx
, left
->n_col
== right
->n_row
, goto error
);
1019 prod
= isl_mat_alloc(left
->ctx
, left
->n_row
, right
->n_col
);
1022 if (left
->n_col
== 0) {
1023 for (i
= 0; i
< prod
->n_row
; ++i
)
1024 isl_seq_clr(prod
->row
[i
], prod
->n_col
);
1026 isl_mat_free(right
);
1029 for (i
= 0; i
< prod
->n_row
; ++i
) {
1030 for (j
= 0; j
< prod
->n_col
; ++j
) {
1031 isl_int_mul(prod
->row
[i
][j
],
1032 left
->row
[i
][0], right
->row
[0][j
]);
1033 for (k
= 1; k
< left
->n_col
; ++k
)
1034 isl_int_addmul(prod
->row
[i
][j
],
1035 left
->row
[i
][k
], right
->row
[k
][j
]);
1039 isl_mat_free(right
);
1043 isl_mat_free(right
);
1047 /* Replace the variables x in the rows q by x' given by x = M x',
1048 * with M the matrix mat.
1050 * If the number of new variables is greater than the original
1051 * number of variables, then the rows q have already been
1052 * preextended. If the new number is smaller, then the coefficients
1053 * of the divs, which are not changed, need to be shifted down.
1054 * The row q may be the equalities, the inequalities or the
1055 * div expressions. In the latter case, has_div is true and
1056 * we need to take into account the extra denominator column.
1058 static int preimage(struct isl_ctx
*ctx
, isl_int
**q
, unsigned n
,
1059 unsigned n_div
, int has_div
, struct isl_mat
*mat
)
1065 if (mat
->n_col
>= mat
->n_row
)
1068 e
= mat
->n_row
- mat
->n_col
;
1070 for (i
= 0; i
< n
; ++i
)
1071 isl_int_mul(q
[i
][0], q
[i
][0], mat
->row
[0][0]);
1072 t
= isl_mat_sub_alloc6(mat
->ctx
, q
, 0, n
, has_div
, mat
->n_row
);
1073 t
= isl_mat_product(t
, mat
);
1076 for (i
= 0; i
< n
; ++i
) {
1077 isl_seq_swp_or_cpy(q
[i
] + has_div
, t
->row
[i
], t
->n_col
);
1078 isl_seq_cpy(q
[i
] + has_div
+ t
->n_col
,
1079 q
[i
] + has_div
+ t
->n_col
+ e
, n_div
);
1080 isl_seq_clr(q
[i
] + has_div
+ t
->n_col
+ n_div
, e
);
1086 /* Replace the variables x in bset by x' given by x = M x', with
1089 * If there are fewer variables x' then there are x, then we perform
1090 * the transformation in place, which that, in principle,
1091 * this frees up some extra variables as the number
1092 * of columns remains constant, but we would have to extend
1093 * the div array too as the number of rows in this array is assumed
1094 * to be equal to extra.
1096 struct isl_basic_set
*isl_basic_set_preimage(struct isl_basic_set
*bset
,
1097 struct isl_mat
*mat
)
1099 struct isl_ctx
*ctx
;
1105 bset
= isl_basic_set_cow(bset
);
1109 isl_assert(ctx
, bset
->dim
->nparam
== 0, goto error
);
1110 isl_assert(ctx
, 1+bset
->dim
->n_out
== mat
->n_row
, goto error
);
1111 isl_assert(ctx
, mat
->n_col
> 0, goto error
);
1113 if (mat
->n_col
> mat
->n_row
) {
1114 bset
= isl_basic_set_extend(bset
, 0, mat
->n_col
-1, 0, 0, 0);
1117 } else if (mat
->n_col
< mat
->n_row
) {
1118 bset
->dim
= isl_space_cow(bset
->dim
);
1121 bset
->dim
->n_out
-= mat
->n_row
- mat
->n_col
;
1124 if (preimage(ctx
, bset
->eq
, bset
->n_eq
, bset
->n_div
, 0,
1125 isl_mat_copy(mat
)) < 0)
1128 if (preimage(ctx
, bset
->ineq
, bset
->n_ineq
, bset
->n_div
, 0,
1129 isl_mat_copy(mat
)) < 0)
1132 if (preimage(ctx
, bset
->div
, bset
->n_div
, bset
->n_div
, 1, mat
) < 0)
1135 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
1136 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_REDUNDANT
);
1137 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED
);
1138 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED_DIVS
);
1139 ISL_F_CLR(bset
, ISL_BASIC_SET_ALL_EQUALITIES
);
1141 bset
= isl_basic_set_simplify(bset
);
1142 bset
= isl_basic_set_finalize(bset
);
1148 isl_basic_set_free(bset
);
1152 struct isl_set
*isl_set_preimage(struct isl_set
*set
, struct isl_mat
*mat
)
1154 struct isl_ctx
*ctx
;
1157 set
= isl_set_cow(set
);
1162 for (i
= 0; i
< set
->n
; ++i
) {
1163 set
->p
[i
] = isl_basic_set_preimage(set
->p
[i
],
1168 if (mat
->n_col
!= mat
->n_row
) {
1169 set
->dim
= isl_space_cow(set
->dim
);
1172 set
->dim
->n_out
+= mat
->n_col
;
1173 set
->dim
->n_out
-= mat
->n_row
;
1176 ISL_F_CLR(set
, ISL_SET_NORMALIZED
);
1184 /* Replace the variables x starting at pos in the rows q
1185 * by x' with x = M x' with M the matrix mat.
1186 * That is, replace the corresponding coefficients c by c M.
1188 static int transform(isl_ctx
*ctx
, isl_int
**q
, unsigned n
,
1189 unsigned pos
, __isl_take isl_mat
*mat
)
1194 t
= isl_mat_sub_alloc6(ctx
, q
, 0, n
, pos
, mat
->n_row
);
1195 t
= isl_mat_product(t
, mat
);
1198 for (i
= 0; i
< n
; ++i
)
1199 isl_seq_swp_or_cpy(q
[i
] + pos
, t
->row
[i
], t
->n_col
);
1204 /* Replace the variables x of type "type" starting at "first" in "bset"
1205 * by x' with x = M x' with M the matrix trans.
1206 * That is, replace the corresponding coefficients c by c M.
1208 * The transformation matrix should be a square matrix.
1210 __isl_give isl_basic_set
*isl_basic_set_transform_dims(
1211 __isl_take isl_basic_set
*bset
, enum isl_dim_type type
, unsigned first
,
1212 __isl_take isl_mat
*trans
)
1217 bset
= isl_basic_set_cow(bset
);
1218 if (!bset
|| !trans
)
1221 ctx
= isl_basic_set_get_ctx(bset
);
1222 if (trans
->n_row
!= trans
->n_col
)
1223 isl_die(trans
->ctx
, isl_error_invalid
,
1224 "expecting square transformation matrix", goto error
);
1225 if (first
+ trans
->n_row
> isl_basic_set_dim(bset
, type
))
1226 isl_die(trans
->ctx
, isl_error_invalid
,
1227 "oversized transformation matrix", goto error
);
1229 pos
= isl_basic_set_offset(bset
, type
) + first
;
1231 if (transform(ctx
, bset
->eq
, bset
->n_eq
, pos
, isl_mat_copy(trans
)) < 0)
1233 if (transform(ctx
, bset
->ineq
, bset
->n_ineq
, pos
,
1234 isl_mat_copy(trans
)) < 0)
1236 if (transform(ctx
, bset
->div
, bset
->n_div
, 1 + pos
,
1237 isl_mat_copy(trans
)) < 0)
1240 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED
);
1241 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED_DIVS
);
1243 isl_mat_free(trans
);
1246 isl_mat_free(trans
);
1247 isl_basic_set_free(bset
);
1251 void isl_mat_print_internal(__isl_keep isl_mat
*mat
, FILE *out
, int indent
)
1256 fprintf(out
, "%*snull mat\n", indent
, "");
1260 if (mat
->n_row
== 0)
1261 fprintf(out
, "%*s[]\n", indent
, "");
1263 for (i
= 0; i
< mat
->n_row
; ++i
) {
1265 fprintf(out
, "%*s[[", indent
, "");
1267 fprintf(out
, "%*s[", indent
+1, "");
1268 for (j
= 0; j
< mat
->n_col
; ++j
) {
1271 isl_int_print(out
, mat
->row
[i
][j
], 0);
1273 if (i
== mat
->n_row
-1)
1274 fprintf(out
, "]]\n");
1276 fprintf(out
, "]\n");
1280 void isl_mat_dump(__isl_keep isl_mat
*mat
)
1282 isl_mat_print_internal(mat
, stderr
, 0);
1285 struct isl_mat
*isl_mat_drop_cols(struct isl_mat
*mat
, unsigned col
, unsigned n
)
1292 mat
= isl_mat_cow(mat
);
1296 if (col
!= mat
->n_col
-n
) {
1297 for (r
= 0; r
< mat
->n_row
; ++r
)
1298 isl_seq_cpy(mat
->row
[r
]+col
, mat
->row
[r
]+col
+n
,
1299 mat
->n_col
- col
- n
);
1305 struct isl_mat
*isl_mat_drop_rows(struct isl_mat
*mat
, unsigned row
, unsigned n
)
1309 mat
= isl_mat_cow(mat
);
1313 for (r
= row
; r
+n
< mat
->n_row
; ++r
)
1314 mat
->row
[r
] = mat
->row
[r
+n
];
1320 __isl_give isl_mat
*isl_mat_insert_cols(__isl_take isl_mat
*mat
,
1321 unsigned col
, unsigned n
)
1330 ext
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
+ n
);
1334 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
, 0, 0, col
);
1335 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
,
1336 col
+ n
, col
, mat
->n_col
- col
);
1345 __isl_give isl_mat
*isl_mat_insert_zero_cols(__isl_take isl_mat
*mat
,
1346 unsigned first
, unsigned n
)
1352 mat
= isl_mat_insert_cols(mat
, first
, n
);
1356 for (i
= 0; i
< mat
->n_row
; ++i
)
1357 isl_seq_clr(mat
->row
[i
] + first
, n
);
1362 __isl_give isl_mat
*isl_mat_add_zero_cols(__isl_take isl_mat
*mat
, unsigned n
)
1367 return isl_mat_insert_zero_cols(mat
, mat
->n_col
, n
);
1370 __isl_give isl_mat
*isl_mat_insert_rows(__isl_take isl_mat
*mat
,
1371 unsigned row
, unsigned n
)
1380 ext
= isl_mat_alloc(mat
->ctx
, mat
->n_row
+ n
, mat
->n_col
);
1384 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, row
, 0, 0, mat
->n_col
);
1385 isl_mat_sub_copy(mat
->ctx
, ext
->row
+ row
+ n
, mat
->row
+ row
,
1386 mat
->n_row
- row
, 0, 0, mat
->n_col
);
1395 __isl_give isl_mat
*isl_mat_add_rows(__isl_take isl_mat
*mat
, unsigned n
)
1400 return isl_mat_insert_rows(mat
, mat
->n_row
, n
);
1403 __isl_give isl_mat
*isl_mat_insert_zero_rows(__isl_take isl_mat
*mat
,
1404 unsigned row
, unsigned n
)
1408 mat
= isl_mat_insert_rows(mat
, row
, n
);
1412 for (i
= 0; i
< n
; ++i
)
1413 isl_seq_clr(mat
->row
[row
+ i
], mat
->n_col
);
1418 __isl_give isl_mat
*isl_mat_add_zero_rows(__isl_take isl_mat
*mat
, unsigned n
)
1423 return isl_mat_insert_zero_rows(mat
, mat
->n_row
, n
);
1426 void isl_mat_col_submul(struct isl_mat
*mat
,
1427 int dst_col
, isl_int f
, int src_col
)
1431 for (i
= 0; i
< mat
->n_row
; ++i
)
1432 isl_int_submul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1435 void isl_mat_col_add(__isl_keep isl_mat
*mat
, int dst_col
, int src_col
)
1442 for (i
= 0; i
< mat
->n_row
; ++i
)
1443 isl_int_add(mat
->row
[i
][dst_col
],
1444 mat
->row
[i
][dst_col
], mat
->row
[i
][src_col
]);
1447 void isl_mat_col_mul(struct isl_mat
*mat
, int dst_col
, isl_int f
, int src_col
)
1451 for (i
= 0; i
< mat
->n_row
; ++i
)
1452 isl_int_mul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1455 struct isl_mat
*isl_mat_unimodular_complete(struct isl_mat
*M
, int row
)
1458 struct isl_mat
*H
= NULL
, *Q
= NULL
;
1463 isl_assert(M
->ctx
, M
->n_row
== M
->n_col
, goto error
);
1465 H
= isl_mat_left_hermite(isl_mat_copy(M
), 0, NULL
, &Q
);
1466 M
->n_row
= M
->n_col
;
1469 for (r
= 0; r
< row
; ++r
)
1470 isl_assert(M
->ctx
, isl_int_is_one(H
->row
[r
][r
]), goto error
);
1471 for (r
= row
; r
< M
->n_row
; ++r
)
1472 isl_seq_cpy(M
->row
[r
], Q
->row
[r
], M
->n_col
);
1483 __isl_give isl_mat
*isl_mat_concat(__isl_take isl_mat
*top
,
1484 __isl_take isl_mat
*bot
)
1486 struct isl_mat
*mat
;
1491 isl_assert(top
->ctx
, top
->n_col
== bot
->n_col
, goto error
);
1492 if (top
->n_row
== 0) {
1496 if (bot
->n_row
== 0) {
1501 mat
= isl_mat_alloc(top
->ctx
, top
->n_row
+ bot
->n_row
, top
->n_col
);
1504 isl_mat_sub_copy(mat
->ctx
, mat
->row
, top
->row
, top
->n_row
,
1506 isl_mat_sub_copy(mat
->ctx
, mat
->row
+ top
->n_row
, bot
->row
, bot
->n_row
,
1517 int isl_mat_is_equal(__isl_keep isl_mat
*mat1
, __isl_keep isl_mat
*mat2
)
1524 if (mat1
->n_row
!= mat2
->n_row
)
1527 if (mat1
->n_col
!= mat2
->n_col
)
1530 for (i
= 0; i
< mat1
->n_row
; ++i
)
1531 if (!isl_seq_eq(mat1
->row
[i
], mat2
->row
[i
], mat1
->n_col
))
1537 __isl_give isl_mat
*isl_mat_from_row_vec(__isl_take isl_vec
*vec
)
1539 struct isl_mat
*mat
;
1543 mat
= isl_mat_alloc(vec
->ctx
, 1, vec
->size
);
1547 isl_seq_cpy(mat
->row
[0], vec
->el
, vec
->size
);
1556 __isl_give isl_mat
*isl_mat_vec_concat(__isl_take isl_mat
*top
,
1557 __isl_take isl_vec
*bot
)
1559 return isl_mat_concat(top
, isl_mat_from_row_vec(bot
));
1562 __isl_give isl_mat
*isl_mat_move_cols(__isl_take isl_mat
*mat
,
1563 unsigned dst_col
, unsigned src_col
, unsigned n
)
1569 if (n
== 0 || dst_col
== src_col
)
1572 res
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
1576 if (dst_col
< src_col
) {
1577 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1579 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1580 dst_col
, src_col
, n
);
1581 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1582 dst_col
+ n
, dst_col
, src_col
- dst_col
);
1583 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1584 src_col
+ n
, src_col
+ n
,
1585 res
->n_col
- src_col
- n
);
1587 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1589 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1590 src_col
, src_col
+ n
, dst_col
- src_col
);
1591 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1592 dst_col
, src_col
, n
);
1593 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1594 dst_col
+ n
, dst_col
+ n
,
1595 res
->n_col
- dst_col
- n
);
1605 void isl_mat_gcd(__isl_keep isl_mat
*mat
, isl_int
*gcd
)
1610 isl_int_set_si(*gcd
, 0);
1615 for (i
= 0; i
< mat
->n_row
; ++i
) {
1616 isl_seq_gcd(mat
->row
[i
], mat
->n_col
, &g
);
1617 isl_int_gcd(*gcd
, *gcd
, g
);
1622 __isl_give isl_mat
*isl_mat_scale_down(__isl_take isl_mat
*mat
, isl_int m
)
1626 if (isl_int_is_one(m
))
1629 mat
= isl_mat_cow(mat
);
1633 for (i
= 0; i
< mat
->n_row
; ++i
)
1634 isl_seq_scale_down(mat
->row
[i
], mat
->row
[i
], m
, mat
->n_col
);
1639 __isl_give isl_mat
*isl_mat_scale_down_row(__isl_take isl_mat
*mat
, int row
,
1642 if (isl_int_is_one(m
))
1645 mat
= isl_mat_cow(mat
);
1649 isl_seq_scale_down(mat
->row
[row
], mat
->row
[row
], m
, mat
->n_col
);
1654 __isl_give isl_mat
*isl_mat_normalize(__isl_take isl_mat
*mat
)
1662 isl_mat_gcd(mat
, &gcd
);
1663 mat
= isl_mat_scale_down(mat
, gcd
);
1669 __isl_give isl_mat
*isl_mat_normalize_row(__isl_take isl_mat
*mat
, int row
)
1671 mat
= isl_mat_cow(mat
);
1675 isl_seq_normalize(mat
->ctx
, mat
->row
[row
], mat
->n_col
);
1680 /* Number of initial non-zero columns.
1682 int isl_mat_initial_non_zero_cols(__isl_keep isl_mat
*mat
)
1689 for (i
= 0; i
< mat
->n_col
; ++i
)
1690 if (row_first_non_zero(mat
->row
, mat
->n_row
, i
) < 0)