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_vec_private.h>
16 #include <isl_space_private.h>
17 #include <isl_val_private.h>
19 isl_ctx
*isl_mat_get_ctx(__isl_keep isl_mat
*mat
)
21 return mat
? mat
->ctx
: NULL
;
24 struct isl_mat
*isl_mat_alloc(struct isl_ctx
*ctx
,
25 unsigned n_row
, unsigned n_col
)
30 mat
= isl_alloc_type(ctx
, struct isl_mat
);
35 mat
->block
= isl_blk_alloc(ctx
, n_row
* n_col
);
36 if (isl_blk_is_error(mat
->block
))
38 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
42 for (i
= 0; i
< n_row
; ++i
)
43 mat
->row
[i
] = mat
->block
.data
+ i
* n_col
;
55 isl_blk_free(ctx
, mat
->block
);
60 struct isl_mat
*isl_mat_extend(struct isl_mat
*mat
,
61 unsigned n_row
, unsigned n_col
)
70 if (mat
->max_col
>= n_col
&& mat
->n_row
>= n_row
) {
71 if (mat
->n_col
< n_col
)
76 if (mat
->max_col
< n_col
) {
77 struct isl_mat
*new_mat
;
79 if (n_row
< mat
->n_row
)
81 new_mat
= isl_mat_alloc(mat
->ctx
, n_row
, n_col
);
84 for (i
= 0; i
< mat
->n_row
; ++i
)
85 isl_seq_cpy(new_mat
->row
[i
], mat
->row
[i
], mat
->n_col
);
90 mat
= isl_mat_cow(mat
);
94 old
= mat
->block
.data
;
95 mat
->block
= isl_blk_extend(mat
->ctx
, mat
->block
, n_row
* mat
->max_col
);
96 if (isl_blk_is_error(mat
->block
))
98 row
= isl_realloc_array(mat
->ctx
, mat
->row
, isl_int
*, n_row
);
103 for (i
= 0; i
< mat
->n_row
; ++i
)
104 mat
->row
[i
] = mat
->block
.data
+ (mat
->row
[i
] - old
);
105 for (i
= mat
->n_row
; i
< n_row
; ++i
)
106 mat
->row
[i
] = mat
->block
.data
+ i
* mat
->max_col
;
108 if (mat
->n_col
< n_col
)
117 __isl_give isl_mat
*isl_mat_sub_alloc6(isl_ctx
*ctx
, isl_int
**row
,
118 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
123 mat
= isl_alloc_type(ctx
, struct isl_mat
);
126 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
129 for (i
= 0; i
< n_row
; ++i
)
130 mat
->row
[i
] = row
[first_row
+i
] + first_col
;
136 mat
->block
= isl_blk_empty();
137 mat
->flags
= ISL_MAT_BORROWED
;
144 __isl_give isl_mat
*isl_mat_sub_alloc(__isl_keep isl_mat
*mat
,
145 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
149 return isl_mat_sub_alloc6(mat
->ctx
, mat
->row
, first_row
, n_row
,
153 void isl_mat_sub_copy(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
154 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
158 for (i
= 0; i
< n_row
; ++i
)
159 isl_seq_cpy(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
162 void isl_mat_sub_neg(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
163 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
167 for (i
= 0; i
< n_row
; ++i
)
168 isl_seq_neg(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
171 struct isl_mat
*isl_mat_copy(struct isl_mat
*mat
)
180 struct isl_mat
*isl_mat_dup(struct isl_mat
*mat
)
183 struct isl_mat
*mat2
;
187 mat2
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
190 for (i
= 0; i
< mat
->n_row
; ++i
)
191 isl_seq_cpy(mat2
->row
[i
], mat
->row
[i
], mat
->n_col
);
195 struct isl_mat
*isl_mat_cow(struct isl_mat
*mat
)
197 struct isl_mat
*mat2
;
201 if (mat
->ref
== 1 && !ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
204 mat2
= isl_mat_dup(mat
);
209 void *isl_mat_free(struct isl_mat
*mat
)
217 if (!ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
218 isl_blk_free(mat
->ctx
, mat
->block
);
219 isl_ctx_deref(mat
->ctx
);
226 int isl_mat_rows(__isl_keep isl_mat
*mat
)
228 return mat
? mat
->n_row
: -1;
231 int isl_mat_cols(__isl_keep isl_mat
*mat
)
233 return mat
? mat
->n_col
: -1;
236 int isl_mat_get_element(__isl_keep isl_mat
*mat
, int row
, int col
, isl_int
*v
)
240 if (row
< 0 || row
>= mat
->n_row
)
241 isl_die(mat
->ctx
, isl_error_invalid
, "row out of range",
243 if (col
< 0 || col
>= mat
->n_col
)
244 isl_die(mat
->ctx
, isl_error_invalid
, "column out of range",
246 isl_int_set(*v
, mat
->row
[row
][col
]);
250 /* Extract the element at row "row", oolumn "col" of "mat".
252 __isl_give isl_val
*isl_mat_get_element_val(__isl_keep isl_mat
*mat
,
259 ctx
= isl_mat_get_ctx(mat
);
260 if (row
< 0 || row
>= mat
->n_row
)
261 isl_die(ctx
, isl_error_invalid
, "row out of range",
263 if (col
< 0 || col
>= mat
->n_col
)
264 isl_die(ctx
, isl_error_invalid
, "column out of range",
266 return isl_val_int_from_isl_int(ctx
, mat
->row
[row
][col
]);
269 __isl_give isl_mat
*isl_mat_set_element(__isl_take isl_mat
*mat
,
270 int row
, int col
, isl_int v
)
272 mat
= isl_mat_cow(mat
);
275 if (row
< 0 || row
>= mat
->n_row
)
276 isl_die(mat
->ctx
, isl_error_invalid
, "row out of range",
278 if (col
< 0 || col
>= mat
->n_col
)
279 isl_die(mat
->ctx
, isl_error_invalid
, "column out of range",
281 isl_int_set(mat
->row
[row
][col
], v
);
288 __isl_give isl_mat
*isl_mat_set_element_si(__isl_take isl_mat
*mat
,
289 int row
, int col
, int v
)
291 mat
= isl_mat_cow(mat
);
294 if (row
< 0 || row
>= mat
->n_row
)
295 isl_die(mat
->ctx
, isl_error_invalid
, "row out of range",
297 if (col
< 0 || col
>= mat
->n_col
)
298 isl_die(mat
->ctx
, isl_error_invalid
, "column out of range",
300 isl_int_set_si(mat
->row
[row
][col
], v
);
307 /* Replace the element at row "row", column "col" of "mat" by "v".
309 __isl_give isl_mat
*isl_mat_set_element_val(__isl_take isl_mat
*mat
,
310 int row
, int col
, __isl_take isl_val
*v
)
313 return isl_mat_free(mat
);
314 if (!isl_val_is_int(v
))
315 isl_die(isl_val_get_ctx(v
), isl_error_invalid
,
316 "expecting integer value", goto error
);
317 mat
= isl_mat_set_element(mat
, row
, col
, v
->n
);
322 return isl_mat_free(mat
);
325 __isl_give isl_mat
*isl_mat_diag(isl_ctx
*ctx
, unsigned n_row
, isl_int d
)
330 mat
= isl_mat_alloc(ctx
, n_row
, n_row
);
333 for (i
= 0; i
< n_row
; ++i
) {
334 isl_seq_clr(mat
->row
[i
], i
);
335 isl_int_set(mat
->row
[i
][i
], d
);
336 isl_seq_clr(mat
->row
[i
]+i
+1, n_row
-(i
+1));
342 __isl_give isl_mat
*isl_mat_identity(isl_ctx
*ctx
, unsigned n_row
)
346 return isl_mat_diag(ctx
, n_row
, ctx
->one
);
349 struct isl_vec
*isl_mat_vec_product(struct isl_mat
*mat
, struct isl_vec
*vec
)
352 struct isl_vec
*prod
;
357 isl_assert(mat
->ctx
, mat
->n_col
== vec
->size
, goto error
);
359 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_row
);
363 for (i
= 0; i
< prod
->size
; ++i
)
364 isl_seq_inner_product(mat
->row
[i
], vec
->el
, vec
->size
,
365 &prod
->block
.data
[i
]);
375 __isl_give isl_vec
*isl_mat_vec_inverse_product(__isl_take isl_mat
*mat
,
376 __isl_take isl_vec
*vec
)
378 struct isl_mat
*vec_mat
;
383 vec_mat
= isl_mat_alloc(vec
->ctx
, vec
->size
, 1);
386 for (i
= 0; i
< vec
->size
; ++i
)
387 isl_int_set(vec_mat
->row
[i
][0], vec
->el
[i
]);
388 vec_mat
= isl_mat_inverse_product(mat
, vec_mat
);
392 vec
= isl_vec_alloc(vec_mat
->ctx
, vec_mat
->n_row
);
394 for (i
= 0; i
< vec
->size
; ++i
)
395 isl_int_set(vec
->el
[i
], vec_mat
->row
[i
][0]);
396 isl_mat_free(vec_mat
);
404 struct isl_vec
*isl_vec_mat_product(struct isl_vec
*vec
, struct isl_mat
*mat
)
407 struct isl_vec
*prod
;
412 isl_assert(mat
->ctx
, mat
->n_row
== vec
->size
, goto error
);
414 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_col
);
418 for (i
= 0; i
< prod
->size
; ++i
) {
419 isl_int_set_si(prod
->el
[i
], 0);
420 for (j
= 0; j
< vec
->size
; ++j
)
421 isl_int_addmul(prod
->el
[i
], vec
->el
[j
], mat
->row
[j
][i
]);
432 struct isl_mat
*isl_mat_aff_direct_sum(struct isl_mat
*left
,
433 struct isl_mat
*right
)
441 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
442 isl_assert(left
->ctx
, left
->n_row
>= 1, goto error
);
443 isl_assert(left
->ctx
, left
->n_col
>= 1, goto error
);
444 isl_assert(left
->ctx
, right
->n_col
>= 1, goto error
);
445 isl_assert(left
->ctx
,
446 isl_seq_first_non_zero(left
->row
[0]+1, left
->n_col
-1) == -1,
448 isl_assert(left
->ctx
,
449 isl_seq_first_non_zero(right
->row
[0]+1, right
->n_col
-1) == -1,
452 sum
= isl_mat_alloc(left
->ctx
, left
->n_row
, left
->n_col
+ right
->n_col
- 1);
455 isl_int_lcm(sum
->row
[0][0], left
->row
[0][0], right
->row
[0][0]);
456 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
457 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
459 isl_seq_clr(sum
->row
[0]+1, sum
->n_col
-1);
460 for (i
= 1; i
< sum
->n_row
; ++i
) {
461 isl_int_mul(sum
->row
[i
][0], left
->row
[0][0], left
->row
[i
][0]);
462 isl_int_addmul(sum
->row
[i
][0],
463 right
->row
[0][0], right
->row
[i
][0]);
464 isl_seq_scale(sum
->row
[i
]+1, left
->row
[i
]+1, left
->row
[0][0],
466 isl_seq_scale(sum
->row
[i
]+left
->n_col
,
467 right
->row
[i
]+1, right
->row
[0][0],
471 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
472 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
482 static void exchange(struct isl_mat
*M
, struct isl_mat
**U
,
483 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
)
486 for (r
= row
; r
< M
->n_row
; ++r
)
487 isl_int_swap(M
->row
[r
][i
], M
->row
[r
][j
]);
489 for (r
= 0; r
< (*U
)->n_row
; ++r
)
490 isl_int_swap((*U
)->row
[r
][i
], (*U
)->row
[r
][j
]);
493 isl_mat_swap_rows(*Q
, i
, j
);
496 static void subtract(struct isl_mat
*M
, struct isl_mat
**U
,
497 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
, isl_int m
)
500 for (r
= row
; r
< M
->n_row
; ++r
)
501 isl_int_submul(M
->row
[r
][j
], m
, M
->row
[r
][i
]);
503 for (r
= 0; r
< (*U
)->n_row
; ++r
)
504 isl_int_submul((*U
)->row
[r
][j
], m
, (*U
)->row
[r
][i
]);
507 for (r
= 0; r
< (*Q
)->n_col
; ++r
)
508 isl_int_addmul((*Q
)->row
[i
][r
], m
, (*Q
)->row
[j
][r
]);
512 static void oppose(struct isl_mat
*M
, struct isl_mat
**U
,
513 struct isl_mat
**Q
, unsigned row
, unsigned col
)
516 for (r
= row
; r
< M
->n_row
; ++r
)
517 isl_int_neg(M
->row
[r
][col
], M
->row
[r
][col
]);
519 for (r
= 0; r
< (*U
)->n_row
; ++r
)
520 isl_int_neg((*U
)->row
[r
][col
], (*U
)->row
[r
][col
]);
523 isl_seq_neg((*Q
)->row
[col
], (*Q
)->row
[col
], (*Q
)->n_col
);
526 /* Given matrix M, compute
531 * with U and Q unimodular matrices and H a matrix in column echelon form
532 * such that on each echelon row the entries in the non-echelon column
533 * are non-negative (if neg == 0) or non-positive (if neg == 1)
534 * and stricly smaller (in absolute value) than the entries in the echelon
536 * If U or Q are NULL, then these matrices are not computed.
538 struct isl_mat
*isl_mat_left_hermite(struct isl_mat
*M
, int neg
,
539 struct isl_mat
**U
, struct isl_mat
**Q
)
554 *U
= isl_mat_identity(M
->ctx
, M
->n_col
);
559 *Q
= isl_mat_identity(M
->ctx
, M
->n_col
);
566 for (row
= 0; row
< M
->n_row
; ++row
) {
568 first
= isl_seq_abs_min_non_zero(M
->row
[row
]+col
, M
->n_col
-col
);
573 exchange(M
, U
, Q
, row
, first
, col
);
574 if (isl_int_is_neg(M
->row
[row
][col
]))
575 oppose(M
, U
, Q
, row
, col
);
577 while ((off
= isl_seq_first_non_zero(M
->row
[row
]+first
,
578 M
->n_col
-first
)) != -1) {
580 isl_int_fdiv_q(c
, M
->row
[row
][first
], M
->row
[row
][col
]);
581 subtract(M
, U
, Q
, row
, col
, first
, c
);
582 if (!isl_int_is_zero(M
->row
[row
][first
]))
583 exchange(M
, U
, Q
, row
, first
, col
);
587 for (i
= 0; i
< col
; ++i
) {
588 if (isl_int_is_zero(M
->row
[row
][i
]))
591 isl_int_cdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
593 isl_int_fdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
594 if (isl_int_is_zero(c
))
596 subtract(M
, U
, Q
, row
, col
, i
, c
);
616 struct isl_mat
*isl_mat_right_kernel(struct isl_mat
*mat
)
619 struct isl_mat
*U
= NULL
;
622 mat
= isl_mat_left_hermite(mat
, 0, &U
, NULL
);
626 for (i
= 0, rank
= 0; rank
< mat
->n_col
; ++rank
) {
627 while (i
< mat
->n_row
&& isl_int_is_zero(mat
->row
[i
][rank
]))
632 K
= isl_mat_alloc(U
->ctx
, U
->n_row
, U
->n_col
- rank
);
635 isl_mat_sub_copy(K
->ctx
, K
->row
, U
->row
, U
->n_row
, 0, rank
, U
->n_col
-rank
);
645 struct isl_mat
*isl_mat_lin_to_aff(struct isl_mat
*mat
)
648 struct isl_mat
*mat2
;
652 mat2
= isl_mat_alloc(mat
->ctx
, 1+mat
->n_row
, 1+mat
->n_col
);
655 isl_int_set_si(mat2
->row
[0][0], 1);
656 isl_seq_clr(mat2
->row
[0]+1, mat
->n_col
);
657 for (i
= 0; i
< mat
->n_row
; ++i
) {
658 isl_int_set_si(mat2
->row
[1+i
][0], 0);
659 isl_seq_cpy(mat2
->row
[1+i
]+1, mat
->row
[i
], mat
->n_col
);
668 /* Given two matrices M1 and M2, return the block matrix
673 __isl_give isl_mat
*isl_mat_diagonal(__isl_take isl_mat
*mat1
,
674 __isl_take isl_mat
*mat2
)
682 mat
= isl_mat_alloc(mat1
->ctx
, mat1
->n_row
+ mat2
->n_row
,
683 mat1
->n_col
+ mat2
->n_col
);
686 for (i
= 0; i
< mat1
->n_row
; ++i
) {
687 isl_seq_cpy(mat
->row
[i
], mat1
->row
[i
], mat1
->n_col
);
688 isl_seq_clr(mat
->row
[i
] + mat1
->n_col
, mat2
->n_col
);
690 for (i
= 0; i
< mat2
->n_row
; ++i
) {
691 isl_seq_clr(mat
->row
[mat1
->n_row
+ i
], mat1
->n_col
);
692 isl_seq_cpy(mat
->row
[mat1
->n_row
+ i
] + mat1
->n_col
,
693 mat2
->row
[i
], mat2
->n_col
);
704 static int row_first_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
708 for (i
= 0; i
< n_row
; ++i
)
709 if (!isl_int_is_zero(row
[i
][col
]))
714 static int row_abs_min_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
716 int i
, min
= row_first_non_zero(row
, n_row
, col
);
719 for (i
= min
+ 1; i
< n_row
; ++i
) {
720 if (isl_int_is_zero(row
[i
][col
]))
722 if (isl_int_abs_lt(row
[i
][col
], row
[min
][col
]))
728 static void inv_exchange(struct isl_mat
*left
, struct isl_mat
*right
,
729 unsigned i
, unsigned j
)
731 left
= isl_mat_swap_rows(left
, i
, j
);
732 right
= isl_mat_swap_rows(right
, i
, j
);
735 static void inv_oppose(
736 struct isl_mat
*left
, struct isl_mat
*right
, unsigned row
)
738 isl_seq_neg(left
->row
[row
]+row
, left
->row
[row
]+row
, left
->n_col
-row
);
739 isl_seq_neg(right
->row
[row
], right
->row
[row
], right
->n_col
);
742 static void inv_subtract(struct isl_mat
*left
, struct isl_mat
*right
,
743 unsigned row
, unsigned i
, isl_int m
)
746 isl_seq_combine(left
->row
[i
]+row
,
747 left
->ctx
->one
, left
->row
[i
]+row
,
748 m
, left
->row
[row
]+row
,
750 isl_seq_combine(right
->row
[i
], right
->ctx
->one
, right
->row
[i
],
751 m
, right
->row
[row
], right
->n_col
);
754 /* Compute inv(left)*right
756 struct isl_mat
*isl_mat_inverse_product(struct isl_mat
*left
,
757 struct isl_mat
*right
)
765 isl_assert(left
->ctx
, left
->n_row
== left
->n_col
, goto error
);
766 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
768 if (left
->n_row
== 0) {
773 left
= isl_mat_cow(left
);
774 right
= isl_mat_cow(right
);
780 for (row
= 0; row
< left
->n_row
; ++row
) {
781 int pivot
, first
, i
, off
;
782 pivot
= row_abs_min_non_zero(left
->row
+row
, left
->n_row
-row
, row
);
786 isl_assert(left
->ctx
, pivot
>= 0, goto error
);
790 inv_exchange(left
, right
, pivot
, row
);
791 if (isl_int_is_neg(left
->row
[row
][row
]))
792 inv_oppose(left
, right
, row
);
794 while ((off
= row_first_non_zero(left
->row
+first
,
795 left
->n_row
-first
, row
)) != -1) {
797 isl_int_fdiv_q(a
, left
->row
[first
][row
],
798 left
->row
[row
][row
]);
799 inv_subtract(left
, right
, row
, first
, a
);
800 if (!isl_int_is_zero(left
->row
[first
][row
]))
801 inv_exchange(left
, right
, row
, first
);
805 for (i
= 0; i
< row
; ++i
) {
806 if (isl_int_is_zero(left
->row
[i
][row
]))
808 isl_int_gcd(a
, left
->row
[row
][row
], left
->row
[i
][row
]);
809 isl_int_divexact(b
, left
->row
[i
][row
], a
);
810 isl_int_divexact(a
, left
->row
[row
][row
], a
);
812 isl_seq_combine(left
->row
[i
] + i
,
814 b
, left
->row
[row
] + i
,
816 isl_seq_combine(right
->row
[i
], a
, right
->row
[i
],
817 b
, right
->row
[row
], right
->n_col
);
822 isl_int_set(a
, left
->row
[0][0]);
823 for (row
= 1; row
< left
->n_row
; ++row
)
824 isl_int_lcm(a
, a
, left
->row
[row
][row
]);
825 if (isl_int_is_zero(a
)){
827 isl_assert(left
->ctx
, 0, goto error
);
829 for (row
= 0; row
< left
->n_row
; ++row
) {
830 isl_int_divexact(left
->row
[row
][row
], a
, left
->row
[row
][row
]);
831 if (isl_int_is_one(left
->row
[row
][row
]))
833 isl_seq_scale(right
->row
[row
], right
->row
[row
],
834 left
->row
[row
][row
], right
->n_col
);
846 void isl_mat_col_scale(struct isl_mat
*mat
, unsigned col
, isl_int m
)
850 for (i
= 0; i
< mat
->n_row
; ++i
)
851 isl_int_mul(mat
->row
[i
][col
], mat
->row
[i
][col
], m
);
854 void isl_mat_col_combine(struct isl_mat
*mat
, unsigned dst
,
855 isl_int m1
, unsigned src1
, isl_int m2
, unsigned src2
)
861 for (i
= 0; i
< mat
->n_row
; ++i
) {
862 isl_int_mul(tmp
, m1
, mat
->row
[i
][src1
]);
863 isl_int_addmul(tmp
, m2
, mat
->row
[i
][src2
]);
864 isl_int_set(mat
->row
[i
][dst
], tmp
);
869 struct isl_mat
*isl_mat_right_inverse(struct isl_mat
*mat
)
875 mat
= isl_mat_cow(mat
);
879 inv
= isl_mat_identity(mat
->ctx
, mat
->n_col
);
880 inv
= isl_mat_cow(inv
);
886 for (row
= 0; row
< mat
->n_row
; ++row
) {
887 int pivot
, first
, i
, off
;
888 pivot
= isl_seq_abs_min_non_zero(mat
->row
[row
]+row
, mat
->n_col
-row
);
892 isl_assert(mat
->ctx
, pivot
>= 0, goto error
);
896 exchange(mat
, &inv
, NULL
, row
, pivot
, row
);
897 if (isl_int_is_neg(mat
->row
[row
][row
]))
898 oppose(mat
, &inv
, NULL
, row
, row
);
900 while ((off
= isl_seq_first_non_zero(mat
->row
[row
]+first
,
901 mat
->n_col
-first
)) != -1) {
903 isl_int_fdiv_q(a
, mat
->row
[row
][first
],
905 subtract(mat
, &inv
, NULL
, row
, row
, first
, a
);
906 if (!isl_int_is_zero(mat
->row
[row
][first
]))
907 exchange(mat
, &inv
, NULL
, row
, row
, first
);
911 for (i
= 0; i
< row
; ++i
) {
912 if (isl_int_is_zero(mat
->row
[row
][i
]))
914 isl_int_gcd(a
, mat
->row
[row
][row
], mat
->row
[row
][i
]);
915 isl_int_divexact(b
, mat
->row
[row
][i
], a
);
916 isl_int_divexact(a
, mat
->row
[row
][row
], a
);
918 isl_mat_col_combine(mat
, i
, a
, i
, b
, row
);
919 isl_mat_col_combine(inv
, i
, a
, i
, b
, row
);
924 isl_int_set(a
, mat
->row
[0][0]);
925 for (row
= 1; row
< mat
->n_row
; ++row
)
926 isl_int_lcm(a
, a
, mat
->row
[row
][row
]);
927 if (isl_int_is_zero(a
)){
931 for (row
= 0; row
< mat
->n_row
; ++row
) {
932 isl_int_divexact(mat
->row
[row
][row
], a
, mat
->row
[row
][row
]);
933 if (isl_int_is_one(mat
->row
[row
][row
]))
935 isl_mat_col_scale(inv
, row
, mat
->row
[row
][row
]);
948 struct isl_mat
*isl_mat_transpose(struct isl_mat
*mat
)
950 struct isl_mat
*transpose
= NULL
;
956 if (mat
->n_col
== mat
->n_row
) {
957 mat
= isl_mat_cow(mat
);
960 for (i
= 0; i
< mat
->n_row
; ++i
)
961 for (j
= i
+ 1; j
< mat
->n_col
; ++j
)
962 isl_int_swap(mat
->row
[i
][j
], mat
->row
[j
][i
]);
965 transpose
= isl_mat_alloc(mat
->ctx
, mat
->n_col
, mat
->n_row
);
968 for (i
= 0; i
< mat
->n_row
; ++i
)
969 for (j
= 0; j
< mat
->n_col
; ++j
)
970 isl_int_set(transpose
->row
[j
][i
], mat
->row
[i
][j
]);
978 struct isl_mat
*isl_mat_swap_cols(struct isl_mat
*mat
, unsigned i
, unsigned j
)
982 mat
= isl_mat_cow(mat
);
985 isl_assert(mat
->ctx
, i
< mat
->n_col
, goto error
);
986 isl_assert(mat
->ctx
, j
< mat
->n_col
, goto error
);
988 for (r
= 0; r
< mat
->n_row
; ++r
)
989 isl_int_swap(mat
->row
[r
][i
], mat
->row
[r
][j
]);
996 struct isl_mat
*isl_mat_swap_rows(struct isl_mat
*mat
, unsigned i
, unsigned j
)
1002 mat
= isl_mat_cow(mat
);
1006 mat
->row
[i
] = mat
->row
[j
];
1011 __isl_give isl_mat
*isl_mat_product(__isl_take isl_mat
*left
,
1012 __isl_take isl_mat
*right
)
1015 struct isl_mat
*prod
;
1017 if (!left
|| !right
)
1019 isl_assert(left
->ctx
, left
->n_col
== right
->n_row
, goto error
);
1020 prod
= isl_mat_alloc(left
->ctx
, left
->n_row
, right
->n_col
);
1023 if (left
->n_col
== 0) {
1024 for (i
= 0; i
< prod
->n_row
; ++i
)
1025 isl_seq_clr(prod
->row
[i
], prod
->n_col
);
1027 isl_mat_free(right
);
1030 for (i
= 0; i
< prod
->n_row
; ++i
) {
1031 for (j
= 0; j
< prod
->n_col
; ++j
) {
1032 isl_int_mul(prod
->row
[i
][j
],
1033 left
->row
[i
][0], right
->row
[0][j
]);
1034 for (k
= 1; k
< left
->n_col
; ++k
)
1035 isl_int_addmul(prod
->row
[i
][j
],
1036 left
->row
[i
][k
], right
->row
[k
][j
]);
1040 isl_mat_free(right
);
1044 isl_mat_free(right
);
1048 /* Replace the variables x in the rows q by x' given by x = M x',
1049 * with M the matrix mat.
1051 * If the number of new variables is greater than the original
1052 * number of variables, then the rows q have already been
1053 * preextended. If the new number is smaller, then the coefficients
1054 * of the divs, which are not changed, need to be shifted down.
1055 * The row q may be the equalities, the inequalities or the
1056 * div expressions. In the latter case, has_div is true and
1057 * we need to take into account the extra denominator column.
1059 static int preimage(struct isl_ctx
*ctx
, isl_int
**q
, unsigned n
,
1060 unsigned n_div
, int has_div
, struct isl_mat
*mat
)
1066 if (mat
->n_col
>= mat
->n_row
)
1069 e
= mat
->n_row
- mat
->n_col
;
1071 for (i
= 0; i
< n
; ++i
)
1072 isl_int_mul(q
[i
][0], q
[i
][0], mat
->row
[0][0]);
1073 t
= isl_mat_sub_alloc6(mat
->ctx
, q
, 0, n
, has_div
, mat
->n_row
);
1074 t
= isl_mat_product(t
, mat
);
1077 for (i
= 0; i
< n
; ++i
) {
1078 isl_seq_swp_or_cpy(q
[i
] + has_div
, t
->row
[i
], t
->n_col
);
1079 isl_seq_cpy(q
[i
] + has_div
+ t
->n_col
,
1080 q
[i
] + has_div
+ t
->n_col
+ e
, n_div
);
1081 isl_seq_clr(q
[i
] + has_div
+ t
->n_col
+ n_div
, e
);
1087 /* Replace the variables x in bset by x' given by x = M x', with
1090 * If there are fewer variables x' then there are x, then we perform
1091 * the transformation in place, which that, in principle,
1092 * this frees up some extra variables as the number
1093 * of columns remains constant, but we would have to extend
1094 * the div array too as the number of rows in this array is assumed
1095 * to be equal to extra.
1097 struct isl_basic_set
*isl_basic_set_preimage(struct isl_basic_set
*bset
,
1098 struct isl_mat
*mat
)
1100 struct isl_ctx
*ctx
;
1106 bset
= isl_basic_set_cow(bset
);
1110 isl_assert(ctx
, bset
->dim
->nparam
== 0, goto error
);
1111 isl_assert(ctx
, 1+bset
->dim
->n_out
== mat
->n_row
, goto error
);
1112 isl_assert(ctx
, mat
->n_col
> 0, goto error
);
1114 if (mat
->n_col
> mat
->n_row
) {
1115 bset
= isl_basic_set_extend(bset
, 0, mat
->n_col
-1, 0, 0, 0);
1118 } else if (mat
->n_col
< mat
->n_row
) {
1119 bset
->dim
= isl_space_cow(bset
->dim
);
1122 bset
->dim
->n_out
-= mat
->n_row
- mat
->n_col
;
1125 if (preimage(ctx
, bset
->eq
, bset
->n_eq
, bset
->n_div
, 0,
1126 isl_mat_copy(mat
)) < 0)
1129 if (preimage(ctx
, bset
->ineq
, bset
->n_ineq
, bset
->n_div
, 0,
1130 isl_mat_copy(mat
)) < 0)
1133 if (preimage(ctx
, bset
->div
, bset
->n_div
, bset
->n_div
, 1, mat
) < 0)
1136 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
1137 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_REDUNDANT
);
1138 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED
);
1139 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED_DIVS
);
1140 ISL_F_CLR(bset
, ISL_BASIC_SET_ALL_EQUALITIES
);
1142 bset
= isl_basic_set_simplify(bset
);
1143 bset
= isl_basic_set_finalize(bset
);
1149 isl_basic_set_free(bset
);
1153 struct isl_set
*isl_set_preimage(struct isl_set
*set
, struct isl_mat
*mat
)
1155 struct isl_ctx
*ctx
;
1158 set
= isl_set_cow(set
);
1163 for (i
= 0; i
< set
->n
; ++i
) {
1164 set
->p
[i
] = isl_basic_set_preimage(set
->p
[i
],
1169 if (mat
->n_col
!= mat
->n_row
) {
1170 set
->dim
= isl_space_cow(set
->dim
);
1173 set
->dim
->n_out
+= mat
->n_col
;
1174 set
->dim
->n_out
-= mat
->n_row
;
1177 ISL_F_CLR(set
, ISL_SET_NORMALIZED
);
1185 /* Replace the variables x starting at pos in the rows q
1186 * by x' with x = M x' with M the matrix mat.
1187 * That is, replace the corresponding coefficients c by c M.
1189 static int transform(isl_ctx
*ctx
, isl_int
**q
, unsigned n
,
1190 unsigned pos
, __isl_take isl_mat
*mat
)
1195 t
= isl_mat_sub_alloc6(ctx
, q
, 0, n
, pos
, mat
->n_row
);
1196 t
= isl_mat_product(t
, mat
);
1199 for (i
= 0; i
< n
; ++i
)
1200 isl_seq_swp_or_cpy(q
[i
] + pos
, t
->row
[i
], t
->n_col
);
1205 /* Replace the variables x of type "type" starting at "first" in "bset"
1206 * by x' with x = M x' with M the matrix trans.
1207 * That is, replace the corresponding coefficients c by c M.
1209 * The transformation matrix should be a square matrix.
1211 __isl_give isl_basic_set
*isl_basic_set_transform_dims(
1212 __isl_take isl_basic_set
*bset
, enum isl_dim_type type
, unsigned first
,
1213 __isl_take isl_mat
*trans
)
1218 bset
= isl_basic_set_cow(bset
);
1219 if (!bset
|| !trans
)
1222 ctx
= isl_basic_set_get_ctx(bset
);
1223 if (trans
->n_row
!= trans
->n_col
)
1224 isl_die(trans
->ctx
, isl_error_invalid
,
1225 "expecting square transformation matrix", goto error
);
1226 if (first
+ trans
->n_row
> isl_basic_set_dim(bset
, type
))
1227 isl_die(trans
->ctx
, isl_error_invalid
,
1228 "oversized transformation matrix", goto error
);
1230 pos
= isl_basic_set_offset(bset
, type
) + first
;
1232 if (transform(ctx
, bset
->eq
, bset
->n_eq
, pos
, isl_mat_copy(trans
)) < 0)
1234 if (transform(ctx
, bset
->ineq
, bset
->n_ineq
, pos
,
1235 isl_mat_copy(trans
)) < 0)
1237 if (transform(ctx
, bset
->div
, bset
->n_div
, 1 + pos
,
1238 isl_mat_copy(trans
)) < 0)
1241 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED
);
1242 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED_DIVS
);
1244 isl_mat_free(trans
);
1247 isl_mat_free(trans
);
1248 isl_basic_set_free(bset
);
1252 void isl_mat_print_internal(__isl_keep isl_mat
*mat
, FILE *out
, int indent
)
1257 fprintf(out
, "%*snull mat\n", indent
, "");
1261 if (mat
->n_row
== 0)
1262 fprintf(out
, "%*s[]\n", indent
, "");
1264 for (i
= 0; i
< mat
->n_row
; ++i
) {
1266 fprintf(out
, "%*s[[", indent
, "");
1268 fprintf(out
, "%*s[", indent
+1, "");
1269 for (j
= 0; j
< mat
->n_col
; ++j
) {
1272 isl_int_print(out
, mat
->row
[i
][j
], 0);
1274 if (i
== mat
->n_row
-1)
1275 fprintf(out
, "]]\n");
1277 fprintf(out
, "]\n");
1281 void isl_mat_dump(__isl_keep isl_mat
*mat
)
1283 isl_mat_print_internal(mat
, stderr
, 0);
1286 struct isl_mat
*isl_mat_drop_cols(struct isl_mat
*mat
, unsigned col
, unsigned n
)
1293 mat
= isl_mat_cow(mat
);
1297 if (col
!= mat
->n_col
-n
) {
1298 for (r
= 0; r
< mat
->n_row
; ++r
)
1299 isl_seq_cpy(mat
->row
[r
]+col
, mat
->row
[r
]+col
+n
,
1300 mat
->n_col
- col
- n
);
1306 struct isl_mat
*isl_mat_drop_rows(struct isl_mat
*mat
, unsigned row
, unsigned n
)
1310 mat
= isl_mat_cow(mat
);
1314 for (r
= row
; r
+n
< mat
->n_row
; ++r
)
1315 mat
->row
[r
] = mat
->row
[r
+n
];
1321 __isl_give isl_mat
*isl_mat_insert_cols(__isl_take isl_mat
*mat
,
1322 unsigned col
, unsigned n
)
1331 ext
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
+ n
);
1335 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
, 0, 0, col
);
1336 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
,
1337 col
+ n
, col
, mat
->n_col
- col
);
1346 __isl_give isl_mat
*isl_mat_insert_zero_cols(__isl_take isl_mat
*mat
,
1347 unsigned first
, unsigned n
)
1353 mat
= isl_mat_insert_cols(mat
, first
, n
);
1357 for (i
= 0; i
< mat
->n_row
; ++i
)
1358 isl_seq_clr(mat
->row
[i
] + first
, n
);
1363 __isl_give isl_mat
*isl_mat_add_zero_cols(__isl_take isl_mat
*mat
, unsigned n
)
1368 return isl_mat_insert_zero_cols(mat
, mat
->n_col
, n
);
1371 __isl_give isl_mat
*isl_mat_insert_rows(__isl_take isl_mat
*mat
,
1372 unsigned row
, unsigned n
)
1381 ext
= isl_mat_alloc(mat
->ctx
, mat
->n_row
+ n
, mat
->n_col
);
1385 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, row
, 0, 0, mat
->n_col
);
1386 isl_mat_sub_copy(mat
->ctx
, ext
->row
+ row
+ n
, mat
->row
+ row
,
1387 mat
->n_row
- row
, 0, 0, mat
->n_col
);
1396 __isl_give isl_mat
*isl_mat_add_rows(__isl_take isl_mat
*mat
, unsigned n
)
1401 return isl_mat_insert_rows(mat
, mat
->n_row
, n
);
1404 __isl_give isl_mat
*isl_mat_insert_zero_rows(__isl_take isl_mat
*mat
,
1405 unsigned row
, unsigned n
)
1409 mat
= isl_mat_insert_rows(mat
, row
, n
);
1413 for (i
= 0; i
< n
; ++i
)
1414 isl_seq_clr(mat
->row
[row
+ i
], mat
->n_col
);
1419 __isl_give isl_mat
*isl_mat_add_zero_rows(__isl_take isl_mat
*mat
, unsigned n
)
1424 return isl_mat_insert_zero_rows(mat
, mat
->n_row
, n
);
1427 void isl_mat_col_submul(struct isl_mat
*mat
,
1428 int dst_col
, isl_int f
, int src_col
)
1432 for (i
= 0; i
< mat
->n_row
; ++i
)
1433 isl_int_submul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1436 void isl_mat_col_add(__isl_keep isl_mat
*mat
, int dst_col
, int src_col
)
1443 for (i
= 0; i
< mat
->n_row
; ++i
)
1444 isl_int_add(mat
->row
[i
][dst_col
],
1445 mat
->row
[i
][dst_col
], mat
->row
[i
][src_col
]);
1448 void isl_mat_col_mul(struct isl_mat
*mat
, int dst_col
, isl_int f
, int src_col
)
1452 for (i
= 0; i
< mat
->n_row
; ++i
)
1453 isl_int_mul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1456 struct isl_mat
*isl_mat_unimodular_complete(struct isl_mat
*M
, int row
)
1459 struct isl_mat
*H
= NULL
, *Q
= NULL
;
1464 isl_assert(M
->ctx
, M
->n_row
== M
->n_col
, goto error
);
1466 H
= isl_mat_left_hermite(isl_mat_copy(M
), 0, NULL
, &Q
);
1467 M
->n_row
= M
->n_col
;
1470 for (r
= 0; r
< row
; ++r
)
1471 isl_assert(M
->ctx
, isl_int_is_one(H
->row
[r
][r
]), goto error
);
1472 for (r
= row
; r
< M
->n_row
; ++r
)
1473 isl_seq_cpy(M
->row
[r
], Q
->row
[r
], M
->n_col
);
1484 __isl_give isl_mat
*isl_mat_concat(__isl_take isl_mat
*top
,
1485 __isl_take isl_mat
*bot
)
1487 struct isl_mat
*mat
;
1492 isl_assert(top
->ctx
, top
->n_col
== bot
->n_col
, goto error
);
1493 if (top
->n_row
== 0) {
1497 if (bot
->n_row
== 0) {
1502 mat
= isl_mat_alloc(top
->ctx
, top
->n_row
+ bot
->n_row
, top
->n_col
);
1505 isl_mat_sub_copy(mat
->ctx
, mat
->row
, top
->row
, top
->n_row
,
1507 isl_mat_sub_copy(mat
->ctx
, mat
->row
+ top
->n_row
, bot
->row
, bot
->n_row
,
1518 int isl_mat_is_equal(__isl_keep isl_mat
*mat1
, __isl_keep isl_mat
*mat2
)
1525 if (mat1
->n_row
!= mat2
->n_row
)
1528 if (mat1
->n_col
!= mat2
->n_col
)
1531 for (i
= 0; i
< mat1
->n_row
; ++i
)
1532 if (!isl_seq_eq(mat1
->row
[i
], mat2
->row
[i
], mat1
->n_col
))
1538 __isl_give isl_mat
*isl_mat_from_row_vec(__isl_take isl_vec
*vec
)
1540 struct isl_mat
*mat
;
1544 mat
= isl_mat_alloc(vec
->ctx
, 1, vec
->size
);
1548 isl_seq_cpy(mat
->row
[0], vec
->el
, vec
->size
);
1557 __isl_give isl_mat
*isl_mat_vec_concat(__isl_take isl_mat
*top
,
1558 __isl_take isl_vec
*bot
)
1560 return isl_mat_concat(top
, isl_mat_from_row_vec(bot
));
1563 __isl_give isl_mat
*isl_mat_move_cols(__isl_take isl_mat
*mat
,
1564 unsigned dst_col
, unsigned src_col
, unsigned n
)
1570 if (n
== 0 || dst_col
== src_col
)
1573 res
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
1577 if (dst_col
< src_col
) {
1578 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1580 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1581 dst_col
, src_col
, n
);
1582 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1583 dst_col
+ n
, dst_col
, src_col
- dst_col
);
1584 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1585 src_col
+ n
, src_col
+ n
,
1586 res
->n_col
- src_col
- n
);
1588 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1590 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1591 src_col
, src_col
+ n
, dst_col
- src_col
);
1592 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1593 dst_col
, src_col
, n
);
1594 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1595 dst_col
+ n
, dst_col
+ n
,
1596 res
->n_col
- dst_col
- n
);
1606 void isl_mat_gcd(__isl_keep isl_mat
*mat
, isl_int
*gcd
)
1611 isl_int_set_si(*gcd
, 0);
1616 for (i
= 0; i
< mat
->n_row
; ++i
) {
1617 isl_seq_gcd(mat
->row
[i
], mat
->n_col
, &g
);
1618 isl_int_gcd(*gcd
, *gcd
, g
);
1623 __isl_give isl_mat
*isl_mat_scale_down(__isl_take isl_mat
*mat
, isl_int m
)
1627 if (isl_int_is_one(m
))
1630 mat
= isl_mat_cow(mat
);
1634 for (i
= 0; i
< mat
->n_row
; ++i
)
1635 isl_seq_scale_down(mat
->row
[i
], mat
->row
[i
], m
, mat
->n_col
);
1640 __isl_give isl_mat
*isl_mat_scale_down_row(__isl_take isl_mat
*mat
, int row
,
1643 if (isl_int_is_one(m
))
1646 mat
= isl_mat_cow(mat
);
1650 isl_seq_scale_down(mat
->row
[row
], mat
->row
[row
], m
, mat
->n_col
);
1655 __isl_give isl_mat
*isl_mat_normalize(__isl_take isl_mat
*mat
)
1663 isl_mat_gcd(mat
, &gcd
);
1664 mat
= isl_mat_scale_down(mat
, gcd
);
1670 __isl_give isl_mat
*isl_mat_normalize_row(__isl_take isl_mat
*mat
, int row
)
1672 mat
= isl_mat_cow(mat
);
1676 isl_seq_normalize(mat
->ctx
, mat
->row
[row
], mat
->n_col
);
1681 /* Number of initial non-zero columns.
1683 int isl_mat_initial_non_zero_cols(__isl_keep isl_mat
*mat
)
1690 for (i
= 0; i
< mat
->n_col
; ++i
)
1691 if (row_first_non_zero(mat
->row
, mat
->n_row
, i
) < 0)