2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2014 Ecole Normale Superieure
5 * Use of this software is governed by the MIT license
7 * Written by Sven Verdoolaege, K.U.Leuven, Departement
8 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
9 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
12 #include <isl_ctx_private.h>
13 #include <isl_map_private.h>
14 #include <isl/space.h>
16 #include <isl_mat_private.h>
17 #include <isl_vec_private.h>
18 #include <isl_space_private.h>
19 #include <isl_val_private.h>
20 #include <isl/deprecated/mat_int.h>
22 isl_ctx
*isl_mat_get_ctx(__isl_keep isl_mat
*mat
)
24 return mat
? mat
->ctx
: NULL
;
27 /* Return a hash value that digests "mat".
29 uint32_t isl_mat_get_hash(__isl_keep isl_mat
*mat
)
37 hash
= isl_hash_init();
38 isl_hash_byte(hash
, mat
->n_row
& 0xFF);
39 isl_hash_byte(hash
, mat
->n_col
& 0xFF);
40 for (i
= 0; i
< mat
->n_row
; ++i
) {
43 row_hash
= isl_seq_get_hash(mat
->row
[i
], mat
->n_col
);
44 isl_hash_hash(hash
, row_hash
);
50 struct isl_mat
*isl_mat_alloc(struct isl_ctx
*ctx
,
51 unsigned n_row
, unsigned n_col
)
56 mat
= isl_alloc_type(ctx
, struct isl_mat
);
61 mat
->block
= isl_blk_alloc(ctx
, n_row
* n_col
);
62 if (isl_blk_is_error(mat
->block
))
64 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
65 if (n_row
&& !mat
->row
)
68 for (i
= 0; i
< n_row
; ++i
)
69 mat
->row
[i
] = mat
->block
.data
+ i
* n_col
;
81 isl_blk_free(ctx
, mat
->block
);
86 struct isl_mat
*isl_mat_extend(struct isl_mat
*mat
,
87 unsigned n_row
, unsigned n_col
)
96 if (mat
->max_col
>= n_col
&& mat
->n_row
>= n_row
) {
97 if (mat
->n_col
< n_col
)
102 if (mat
->max_col
< n_col
) {
103 struct isl_mat
*new_mat
;
105 if (n_row
< mat
->n_row
)
107 new_mat
= isl_mat_alloc(mat
->ctx
, n_row
, n_col
);
110 for (i
= 0; i
< mat
->n_row
; ++i
)
111 isl_seq_cpy(new_mat
->row
[i
], mat
->row
[i
], mat
->n_col
);
116 mat
= isl_mat_cow(mat
);
120 old
= mat
->block
.data
;
121 mat
->block
= isl_blk_extend(mat
->ctx
, mat
->block
, n_row
* mat
->max_col
);
122 if (isl_blk_is_error(mat
->block
))
124 row
= isl_realloc_array(mat
->ctx
, mat
->row
, isl_int
*, n_row
);
129 for (i
= 0; i
< mat
->n_row
; ++i
)
130 mat
->row
[i
] = mat
->block
.data
+ (mat
->row
[i
] - old
);
131 for (i
= mat
->n_row
; i
< n_row
; ++i
)
132 mat
->row
[i
] = mat
->block
.data
+ i
* mat
->max_col
;
134 if (mat
->n_col
< n_col
)
143 __isl_give isl_mat
*isl_mat_sub_alloc6(isl_ctx
*ctx
, isl_int
**row
,
144 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
149 mat
= isl_alloc_type(ctx
, struct isl_mat
);
152 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
153 if (n_row
&& !mat
->row
)
155 for (i
= 0; i
< n_row
; ++i
)
156 mat
->row
[i
] = row
[first_row
+i
] + first_col
;
162 mat
->block
= isl_blk_empty();
163 mat
->flags
= ISL_MAT_BORROWED
;
170 __isl_give isl_mat
*isl_mat_sub_alloc(__isl_keep isl_mat
*mat
,
171 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
175 return isl_mat_sub_alloc6(mat
->ctx
, mat
->row
, first_row
, n_row
,
179 void isl_mat_sub_copy(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
180 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
184 for (i
= 0; i
< n_row
; ++i
)
185 isl_seq_cpy(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
188 void isl_mat_sub_neg(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
189 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
193 for (i
= 0; i
< n_row
; ++i
)
194 isl_seq_neg(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
197 __isl_give isl_mat
*isl_mat_copy(__isl_keep isl_mat
*mat
)
206 __isl_give isl_mat
*isl_mat_dup(__isl_keep isl_mat
*mat
)
209 struct isl_mat
*mat2
;
213 mat2
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
216 for (i
= 0; i
< mat
->n_row
; ++i
)
217 isl_seq_cpy(mat2
->row
[i
], mat
->row
[i
], mat
->n_col
);
221 __isl_give isl_mat
*isl_mat_cow(__isl_take isl_mat
*mat
)
223 struct isl_mat
*mat2
;
227 if (mat
->ref
== 1 && !ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
230 mat2
= isl_mat_dup(mat
);
235 __isl_null isl_mat
*isl_mat_free(__isl_take isl_mat
*mat
)
243 if (!ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
244 isl_blk_free(mat
->ctx
, mat
->block
);
245 isl_ctx_deref(mat
->ctx
);
252 int isl_mat_rows(__isl_keep isl_mat
*mat
)
254 return mat
? mat
->n_row
: -1;
257 int isl_mat_cols(__isl_keep isl_mat
*mat
)
259 return mat
? mat
->n_col
: -1;
262 /* Check that "col" is a valid column position for "mat".
264 static isl_stat
check_col(__isl_keep isl_mat
*mat
, int col
)
267 return isl_stat_error
;
268 if (col
< 0 || col
>= mat
->n_col
)
269 isl_die(isl_mat_get_ctx(mat
), isl_error_invalid
,
270 "column out of range", return isl_stat_error
);
274 /* Check that "row" is a valid row position for "mat".
276 static isl_stat
check_row(__isl_keep isl_mat
*mat
, int row
)
279 return isl_stat_error
;
280 if (row
< 0 || row
>= mat
->n_row
)
281 isl_die(isl_mat_get_ctx(mat
), isl_error_invalid
,
282 "row out of range", return isl_stat_error
);
286 /* Check that there are "n" columns starting at position "first" in "mat".
288 static isl_stat
check_col_range(__isl_keep isl_mat
*mat
, unsigned first
,
292 return isl_stat_error
;
293 if (first
+ n
> mat
->n_col
|| first
+ n
< first
)
294 isl_die(isl_mat_get_ctx(mat
), isl_error_invalid
,
295 "column position or range out of bounds",
296 return isl_stat_error
);
300 /* Check that there are "n" rows starting at position "first" in "mat".
302 static isl_stat
check_row_range(__isl_keep isl_mat
*mat
, unsigned first
,
306 return isl_stat_error
;
307 if (first
+ n
> mat
->n_row
|| first
+ n
< first
)
308 isl_die(isl_mat_get_ctx(mat
), isl_error_invalid
,
309 "row position or range out of bounds",
310 return isl_stat_error
);
314 int isl_mat_get_element(__isl_keep isl_mat
*mat
, int row
, int col
, isl_int
*v
)
316 if (check_row(mat
, row
) < 0)
318 if (check_col(mat
, col
) < 0)
320 isl_int_set(*v
, mat
->row
[row
][col
]);
324 /* Extract the element at row "row", oolumn "col" of "mat".
326 __isl_give isl_val
*isl_mat_get_element_val(__isl_keep isl_mat
*mat
,
331 if (check_row(mat
, row
) < 0)
333 if (check_col(mat
, col
) < 0)
335 ctx
= isl_mat_get_ctx(mat
);
336 return isl_val_int_from_isl_int(ctx
, mat
->row
[row
][col
]);
339 __isl_give isl_mat
*isl_mat_set_element(__isl_take isl_mat
*mat
,
340 int row
, int col
, isl_int v
)
342 mat
= isl_mat_cow(mat
);
343 if (check_row(mat
, row
) < 0)
344 return isl_mat_free(mat
);
345 if (check_col(mat
, col
) < 0)
346 return isl_mat_free(mat
);
347 isl_int_set(mat
->row
[row
][col
], v
);
351 __isl_give isl_mat
*isl_mat_set_element_si(__isl_take isl_mat
*mat
,
352 int row
, int col
, int v
)
354 mat
= isl_mat_cow(mat
);
355 if (check_row(mat
, row
) < 0)
356 return isl_mat_free(mat
);
357 if (check_col(mat
, col
) < 0)
358 return isl_mat_free(mat
);
359 isl_int_set_si(mat
->row
[row
][col
], v
);
363 /* Replace the element at row "row", column "col" of "mat" by "v".
365 __isl_give isl_mat
*isl_mat_set_element_val(__isl_take isl_mat
*mat
,
366 int row
, int col
, __isl_take isl_val
*v
)
369 return isl_mat_free(mat
);
370 if (!isl_val_is_int(v
))
371 isl_die(isl_val_get_ctx(v
), isl_error_invalid
,
372 "expecting integer value", goto error
);
373 mat
= isl_mat_set_element(mat
, row
, col
, v
->n
);
378 return isl_mat_free(mat
);
381 __isl_give isl_mat
*isl_mat_diag(isl_ctx
*ctx
, unsigned n_row
, isl_int d
)
386 mat
= isl_mat_alloc(ctx
, n_row
, n_row
);
389 for (i
= 0; i
< n_row
; ++i
) {
390 isl_seq_clr(mat
->row
[i
], i
);
391 isl_int_set(mat
->row
[i
][i
], d
);
392 isl_seq_clr(mat
->row
[i
]+i
+1, n_row
-(i
+1));
398 /* Create an "n_row" by "n_col" matrix with zero elements.
400 __isl_give isl_mat
*isl_mat_zero(isl_ctx
*ctx
, unsigned n_row
, unsigned n_col
)
405 mat
= isl_mat_alloc(ctx
, n_row
, n_col
);
408 for (i
= 0; i
< n_row
; ++i
)
409 isl_seq_clr(mat
->row
[i
], n_col
);
414 __isl_give isl_mat
*isl_mat_identity(isl_ctx
*ctx
, unsigned n_row
)
418 return isl_mat_diag(ctx
, n_row
, ctx
->one
);
421 /* Is "mat" a (possibly scaled) identity matrix?
423 int isl_mat_is_scaled_identity(__isl_keep isl_mat
*mat
)
429 if (mat
->n_row
!= mat
->n_col
)
432 for (i
= 0; i
< mat
->n_row
; ++i
) {
433 if (isl_seq_first_non_zero(mat
->row
[i
], i
) != -1)
435 if (isl_int_ne(mat
->row
[0][0], mat
->row
[i
][i
]))
437 if (isl_seq_first_non_zero(mat
->row
[i
] + i
+ 1,
438 mat
->n_col
- (i
+ 1)) != -1)
445 __isl_give isl_vec
*isl_mat_vec_product(__isl_take isl_mat
*mat
,
446 __isl_take isl_vec
*vec
)
449 struct isl_vec
*prod
;
454 isl_assert(mat
->ctx
, mat
->n_col
== vec
->size
, goto error
);
456 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_row
);
460 for (i
= 0; i
< prod
->size
; ++i
)
461 isl_seq_inner_product(mat
->row
[i
], vec
->el
, vec
->size
,
462 &prod
->block
.data
[i
]);
472 __isl_give isl_vec
*isl_mat_vec_inverse_product(__isl_take isl_mat
*mat
,
473 __isl_take isl_vec
*vec
)
475 struct isl_mat
*vec_mat
;
480 vec_mat
= isl_mat_alloc(vec
->ctx
, vec
->size
, 1);
483 for (i
= 0; i
< vec
->size
; ++i
)
484 isl_int_set(vec_mat
->row
[i
][0], vec
->el
[i
]);
485 vec_mat
= isl_mat_inverse_product(mat
, vec_mat
);
489 vec
= isl_vec_alloc(vec_mat
->ctx
, vec_mat
->n_row
);
491 for (i
= 0; i
< vec
->size
; ++i
)
492 isl_int_set(vec
->el
[i
], vec_mat
->row
[i
][0]);
493 isl_mat_free(vec_mat
);
501 __isl_give isl_vec
*isl_vec_mat_product(__isl_take isl_vec
*vec
,
502 __isl_take isl_mat
*mat
)
505 struct isl_vec
*prod
;
510 isl_assert(mat
->ctx
, mat
->n_row
== vec
->size
, goto error
);
512 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_col
);
516 for (i
= 0; i
< prod
->size
; ++i
) {
517 isl_int_set_si(prod
->el
[i
], 0);
518 for (j
= 0; j
< vec
->size
; ++j
)
519 isl_int_addmul(prod
->el
[i
], vec
->el
[j
], mat
->row
[j
][i
]);
530 __isl_give isl_mat
*isl_mat_aff_direct_sum(__isl_take isl_mat
*left
,
531 __isl_take isl_mat
*right
)
539 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
540 isl_assert(left
->ctx
, left
->n_row
>= 1, goto error
);
541 isl_assert(left
->ctx
, left
->n_col
>= 1, goto error
);
542 isl_assert(left
->ctx
, right
->n_col
>= 1, goto error
);
543 isl_assert(left
->ctx
,
544 isl_seq_first_non_zero(left
->row
[0]+1, left
->n_col
-1) == -1,
546 isl_assert(left
->ctx
,
547 isl_seq_first_non_zero(right
->row
[0]+1, right
->n_col
-1) == -1,
550 sum
= isl_mat_alloc(left
->ctx
, left
->n_row
, left
->n_col
+ right
->n_col
- 1);
553 isl_int_lcm(sum
->row
[0][0], left
->row
[0][0], right
->row
[0][0]);
554 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
555 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
557 isl_seq_clr(sum
->row
[0]+1, sum
->n_col
-1);
558 for (i
= 1; i
< sum
->n_row
; ++i
) {
559 isl_int_mul(sum
->row
[i
][0], left
->row
[0][0], left
->row
[i
][0]);
560 isl_int_addmul(sum
->row
[i
][0],
561 right
->row
[0][0], right
->row
[i
][0]);
562 isl_seq_scale(sum
->row
[i
]+1, left
->row
[i
]+1, left
->row
[0][0],
564 isl_seq_scale(sum
->row
[i
]+left
->n_col
,
565 right
->row
[i
]+1, right
->row
[0][0],
569 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
570 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
580 static void exchange(struct isl_mat
*M
, struct isl_mat
**U
,
581 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
)
584 for (r
= row
; r
< M
->n_row
; ++r
)
585 isl_int_swap(M
->row
[r
][i
], M
->row
[r
][j
]);
587 for (r
= 0; r
< (*U
)->n_row
; ++r
)
588 isl_int_swap((*U
)->row
[r
][i
], (*U
)->row
[r
][j
]);
591 isl_mat_swap_rows(*Q
, i
, j
);
594 static void subtract(struct isl_mat
*M
, struct isl_mat
**U
,
595 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
, isl_int m
)
598 for (r
= row
; r
< M
->n_row
; ++r
)
599 isl_int_submul(M
->row
[r
][j
], m
, M
->row
[r
][i
]);
601 for (r
= 0; r
< (*U
)->n_row
; ++r
)
602 isl_int_submul((*U
)->row
[r
][j
], m
, (*U
)->row
[r
][i
]);
605 for (r
= 0; r
< (*Q
)->n_col
; ++r
)
606 isl_int_addmul((*Q
)->row
[i
][r
], m
, (*Q
)->row
[j
][r
]);
610 static void oppose(struct isl_mat
*M
, struct isl_mat
**U
,
611 struct isl_mat
**Q
, unsigned row
, unsigned col
)
614 for (r
= row
; r
< M
->n_row
; ++r
)
615 isl_int_neg(M
->row
[r
][col
], M
->row
[r
][col
]);
617 for (r
= 0; r
< (*U
)->n_row
; ++r
)
618 isl_int_neg((*U
)->row
[r
][col
], (*U
)->row
[r
][col
]);
621 isl_seq_neg((*Q
)->row
[col
], (*Q
)->row
[col
], (*Q
)->n_col
);
624 /* Given matrix M, compute
629 * with U and Q unimodular matrices and H a matrix in column echelon form
630 * such that on each echelon row the entries in the non-echelon column
631 * are non-negative (if neg == 0) or non-positive (if neg == 1)
632 * and strictly smaller (in absolute value) than the entries in the echelon
634 * If U or Q are NULL, then these matrices are not computed.
636 __isl_give isl_mat
*isl_mat_left_hermite(__isl_take isl_mat
*M
, int neg
,
637 __isl_give isl_mat
**U
, __isl_give isl_mat
**Q
)
652 *U
= isl_mat_identity(M
->ctx
, M
->n_col
);
657 *Q
= isl_mat_identity(M
->ctx
, M
->n_col
);
664 for (row
= 0; row
< M
->n_row
; ++row
) {
666 first
= isl_seq_abs_min_non_zero(M
->row
[row
]+col
, M
->n_col
-col
);
671 exchange(M
, U
, Q
, row
, first
, col
);
672 if (isl_int_is_neg(M
->row
[row
][col
]))
673 oppose(M
, U
, Q
, row
, col
);
675 while ((off
= isl_seq_first_non_zero(M
->row
[row
]+first
,
676 M
->n_col
-first
)) != -1) {
678 isl_int_fdiv_q(c
, M
->row
[row
][first
], M
->row
[row
][col
]);
679 subtract(M
, U
, Q
, row
, col
, first
, c
);
680 if (!isl_int_is_zero(M
->row
[row
][first
]))
681 exchange(M
, U
, Q
, row
, first
, col
);
685 for (i
= 0; i
< col
; ++i
) {
686 if (isl_int_is_zero(M
->row
[row
][i
]))
689 isl_int_cdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
691 isl_int_fdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
692 if (isl_int_is_zero(c
))
694 subtract(M
, U
, Q
, row
, col
, i
, c
);
714 /* Use row "row" of "mat" to eliminate column "col" from all other rows.
716 static __isl_give isl_mat
*eliminate(__isl_take isl_mat
*mat
, int row
, int col
)
724 ctx
= isl_mat_get_ctx(mat
);
725 nr
= isl_mat_rows(mat
);
726 nc
= isl_mat_cols(mat
);
728 for (k
= 0; k
< nr
; ++k
) {
731 if (isl_int_is_zero(mat
->row
[k
][col
]))
733 mat
= isl_mat_cow(mat
);
736 isl_seq_elim(mat
->row
[k
], mat
->row
[row
], col
, nc
, NULL
);
737 isl_seq_normalize(ctx
, mat
->row
[k
], nc
);
743 /* Perform Gaussian elimination on the rows of "mat", but start
744 * from the final row and the final column.
745 * Any zero rows that result from the elimination are removed.
747 * In particular, for each column from last to first,
748 * look for the last row with a non-zero coefficient in that column,
749 * move it last (but before other rows moved last in previous steps) and
750 * use it to eliminate the column from the other rows.
752 __isl_give isl_mat
*isl_mat_reverse_gauss(__isl_take isl_mat
*mat
)
754 int k
, row
, last
, nr
, nc
;
759 nr
= isl_mat_rows(mat
);
760 nc
= isl_mat_cols(mat
);
763 for (row
= nr
- 1; row
>= 0; --row
) {
764 for (; last
>= 0; --last
) {
765 for (k
= row
; k
>= 0; --k
)
766 if (!isl_int_is_zero(mat
->row
[k
][last
]))
774 mat
= isl_mat_swap_rows(mat
, k
, row
);
777 if (isl_int_is_neg(mat
->row
[row
][last
]))
778 mat
= isl_mat_row_neg(mat
, row
);
779 mat
= eliminate(mat
, row
, last
);
783 mat
= isl_mat_drop_rows(mat
, 0, row
+ 1);
788 /* Negate the lexicographically negative rows of "mat" such that
789 * all rows in the result are lexicographically non-negative.
791 __isl_give isl_mat
*isl_mat_lexnonneg_rows(__isl_take isl_mat
*mat
)
798 nr
= isl_mat_rows(mat
);
799 nc
= isl_mat_cols(mat
);
801 for (i
= 0; i
< nr
; ++i
) {
804 pos
= isl_seq_first_non_zero(mat
->row
[i
], nc
);
807 if (isl_int_is_nonneg(mat
->row
[i
][pos
]))
809 mat
= isl_mat_row_neg(mat
, i
);
817 struct isl_mat
*isl_mat_right_kernel(struct isl_mat
*mat
)
820 struct isl_mat
*U
= NULL
;
823 mat
= isl_mat_left_hermite(mat
, 0, &U
, NULL
);
827 for (i
= 0, rank
= 0; rank
< mat
->n_col
; ++rank
) {
828 while (i
< mat
->n_row
&& isl_int_is_zero(mat
->row
[i
][rank
]))
833 K
= isl_mat_alloc(U
->ctx
, U
->n_row
, U
->n_col
- rank
);
836 isl_mat_sub_copy(K
->ctx
, K
->row
, U
->row
, U
->n_row
, 0, rank
, U
->n_col
-rank
);
846 __isl_give isl_mat
*isl_mat_lin_to_aff(__isl_take isl_mat
*mat
)
849 struct isl_mat
*mat2
;
853 mat2
= isl_mat_alloc(mat
->ctx
, 1+mat
->n_row
, 1+mat
->n_col
);
856 isl_int_set_si(mat2
->row
[0][0], 1);
857 isl_seq_clr(mat2
->row
[0]+1, mat
->n_col
);
858 for (i
= 0; i
< mat
->n_row
; ++i
) {
859 isl_int_set_si(mat2
->row
[1+i
][0], 0);
860 isl_seq_cpy(mat2
->row
[1+i
]+1, mat
->row
[i
], mat
->n_col
);
869 /* Given two matrices M1 and M2, return the block matrix
874 __isl_give isl_mat
*isl_mat_diagonal(__isl_take isl_mat
*mat1
,
875 __isl_take isl_mat
*mat2
)
883 mat
= isl_mat_alloc(mat1
->ctx
, mat1
->n_row
+ mat2
->n_row
,
884 mat1
->n_col
+ mat2
->n_col
);
887 for (i
= 0; i
< mat1
->n_row
; ++i
) {
888 isl_seq_cpy(mat
->row
[i
], mat1
->row
[i
], mat1
->n_col
);
889 isl_seq_clr(mat
->row
[i
] + mat1
->n_col
, mat2
->n_col
);
891 for (i
= 0; i
< mat2
->n_row
; ++i
) {
892 isl_seq_clr(mat
->row
[mat1
->n_row
+ i
], mat1
->n_col
);
893 isl_seq_cpy(mat
->row
[mat1
->n_row
+ i
] + mat1
->n_col
,
894 mat2
->row
[i
], mat2
->n_col
);
905 static int row_first_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
909 for (i
= 0; i
< n_row
; ++i
)
910 if (!isl_int_is_zero(row
[i
][col
]))
915 static int row_abs_min_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
917 int i
, min
= row_first_non_zero(row
, n_row
, col
);
920 for (i
= min
+ 1; i
< n_row
; ++i
) {
921 if (isl_int_is_zero(row
[i
][col
]))
923 if (isl_int_abs_lt(row
[i
][col
], row
[min
][col
]))
929 static isl_stat
inv_exchange(__isl_keep isl_mat
**left
,
930 __isl_keep isl_mat
**right
, unsigned i
, unsigned j
)
932 *left
= isl_mat_swap_rows(*left
, i
, j
);
933 *right
= isl_mat_swap_rows(*right
, i
, j
);
935 if (!*left
|| !*right
)
936 return isl_stat_error
;
940 static void inv_oppose(
941 struct isl_mat
*left
, struct isl_mat
*right
, unsigned row
)
943 isl_seq_neg(left
->row
[row
]+row
, left
->row
[row
]+row
, left
->n_col
-row
);
944 isl_seq_neg(right
->row
[row
], right
->row
[row
], right
->n_col
);
947 static void inv_subtract(struct isl_mat
*left
, struct isl_mat
*right
,
948 unsigned row
, unsigned i
, isl_int m
)
951 isl_seq_combine(left
->row
[i
]+row
,
952 left
->ctx
->one
, left
->row
[i
]+row
,
953 m
, left
->row
[row
]+row
,
955 isl_seq_combine(right
->row
[i
], right
->ctx
->one
, right
->row
[i
],
956 m
, right
->row
[row
], right
->n_col
);
959 /* Compute inv(left)*right
961 __isl_give isl_mat
*isl_mat_inverse_product(__isl_take isl_mat
*left
,
962 __isl_take isl_mat
*right
)
970 isl_assert(left
->ctx
, left
->n_row
== left
->n_col
, goto error
);
971 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
973 if (left
->n_row
== 0) {
978 left
= isl_mat_cow(left
);
979 right
= isl_mat_cow(right
);
985 for (row
= 0; row
< left
->n_row
; ++row
) {
986 int pivot
, first
, i
, off
;
987 pivot
= row_abs_min_non_zero(left
->row
+row
, left
->n_row
-row
, row
);
991 isl_assert(left
->ctx
, pivot
>= 0, goto error
);
995 if (inv_exchange(&left
, &right
, pivot
, row
) < 0)
997 if (isl_int_is_neg(left
->row
[row
][row
]))
998 inv_oppose(left
, right
, row
);
1000 while ((off
= row_first_non_zero(left
->row
+first
,
1001 left
->n_row
-first
, row
)) != -1) {
1003 isl_int_fdiv_q(a
, left
->row
[first
][row
],
1004 left
->row
[row
][row
]);
1005 inv_subtract(left
, right
, row
, first
, a
);
1006 if (!isl_int_is_zero(left
->row
[first
][row
])) {
1007 if (inv_exchange(&left
, &right
, row
, first
) < 0)
1013 for (i
= 0; i
< row
; ++i
) {
1014 if (isl_int_is_zero(left
->row
[i
][row
]))
1016 isl_int_gcd(a
, left
->row
[row
][row
], left
->row
[i
][row
]);
1017 isl_int_divexact(b
, left
->row
[i
][row
], a
);
1018 isl_int_divexact(a
, left
->row
[row
][row
], a
);
1020 isl_seq_combine(left
->row
[i
] + i
,
1021 a
, left
->row
[i
] + i
,
1022 b
, left
->row
[row
] + i
,
1024 isl_seq_combine(right
->row
[i
], a
, right
->row
[i
],
1025 b
, right
->row
[row
], right
->n_col
);
1030 isl_int_set(a
, left
->row
[0][0]);
1031 for (row
= 1; row
< left
->n_row
; ++row
)
1032 isl_int_lcm(a
, a
, left
->row
[row
][row
]);
1033 if (isl_int_is_zero(a
)){
1035 isl_assert(left
->ctx
, 0, goto error
);
1037 for (row
= 0; row
< left
->n_row
; ++row
) {
1038 isl_int_divexact(left
->row
[row
][row
], a
, left
->row
[row
][row
]);
1039 if (isl_int_is_one(left
->row
[row
][row
]))
1041 isl_seq_scale(right
->row
[row
], right
->row
[row
],
1042 left
->row
[row
][row
], right
->n_col
);
1050 isl_mat_free(right
);
1054 void isl_mat_col_scale(struct isl_mat
*mat
, unsigned col
, isl_int m
)
1058 for (i
= 0; i
< mat
->n_row
; ++i
)
1059 isl_int_mul(mat
->row
[i
][col
], mat
->row
[i
][col
], m
);
1062 void isl_mat_col_combine(struct isl_mat
*mat
, unsigned dst
,
1063 isl_int m1
, unsigned src1
, isl_int m2
, unsigned src2
)
1069 for (i
= 0; i
< mat
->n_row
; ++i
) {
1070 isl_int_mul(tmp
, m1
, mat
->row
[i
][src1
]);
1071 isl_int_addmul(tmp
, m2
, mat
->row
[i
][src2
]);
1072 isl_int_set(mat
->row
[i
][dst
], tmp
);
1077 __isl_give isl_mat
*isl_mat_right_inverse(__isl_take isl_mat
*mat
)
1079 struct isl_mat
*inv
;
1083 mat
= isl_mat_cow(mat
);
1087 inv
= isl_mat_identity(mat
->ctx
, mat
->n_col
);
1088 inv
= isl_mat_cow(inv
);
1094 for (row
= 0; row
< mat
->n_row
; ++row
) {
1095 int pivot
, first
, i
, off
;
1096 pivot
= isl_seq_abs_min_non_zero(mat
->row
[row
]+row
, mat
->n_col
-row
);
1100 isl_assert(mat
->ctx
, pivot
>= 0, goto error
);
1104 exchange(mat
, &inv
, NULL
, row
, pivot
, row
);
1105 if (isl_int_is_neg(mat
->row
[row
][row
]))
1106 oppose(mat
, &inv
, NULL
, row
, row
);
1108 while ((off
= isl_seq_first_non_zero(mat
->row
[row
]+first
,
1109 mat
->n_col
-first
)) != -1) {
1111 isl_int_fdiv_q(a
, mat
->row
[row
][first
],
1112 mat
->row
[row
][row
]);
1113 subtract(mat
, &inv
, NULL
, row
, row
, first
, a
);
1114 if (!isl_int_is_zero(mat
->row
[row
][first
]))
1115 exchange(mat
, &inv
, NULL
, row
, row
, first
);
1119 for (i
= 0; i
< row
; ++i
) {
1120 if (isl_int_is_zero(mat
->row
[row
][i
]))
1122 isl_int_gcd(a
, mat
->row
[row
][row
], mat
->row
[row
][i
]);
1123 isl_int_divexact(b
, mat
->row
[row
][i
], a
);
1124 isl_int_divexact(a
, mat
->row
[row
][row
], a
);
1126 isl_mat_col_combine(mat
, i
, a
, i
, b
, row
);
1127 isl_mat_col_combine(inv
, i
, a
, i
, b
, row
);
1132 isl_int_set(a
, mat
->row
[0][0]);
1133 for (row
= 1; row
< mat
->n_row
; ++row
)
1134 isl_int_lcm(a
, a
, mat
->row
[row
][row
]);
1135 if (isl_int_is_zero(a
)){
1139 for (row
= 0; row
< mat
->n_row
; ++row
) {
1140 isl_int_divexact(mat
->row
[row
][row
], a
, mat
->row
[row
][row
]);
1141 if (isl_int_is_one(mat
->row
[row
][row
]))
1143 isl_mat_col_scale(inv
, row
, mat
->row
[row
][row
]);
1156 __isl_give isl_mat
*isl_mat_transpose(__isl_take isl_mat
*mat
)
1158 struct isl_mat
*transpose
= NULL
;
1164 if (mat
->n_col
== mat
->n_row
) {
1165 mat
= isl_mat_cow(mat
);
1168 for (i
= 0; i
< mat
->n_row
; ++i
)
1169 for (j
= i
+ 1; j
< mat
->n_col
; ++j
)
1170 isl_int_swap(mat
->row
[i
][j
], mat
->row
[j
][i
]);
1173 transpose
= isl_mat_alloc(mat
->ctx
, mat
->n_col
, mat
->n_row
);
1176 for (i
= 0; i
< mat
->n_row
; ++i
)
1177 for (j
= 0; j
< mat
->n_col
; ++j
)
1178 isl_int_set(transpose
->row
[j
][i
], mat
->row
[i
][j
]);
1186 __isl_give isl_mat
*isl_mat_swap_cols(__isl_take isl_mat
*mat
,
1187 unsigned i
, unsigned j
)
1191 mat
= isl_mat_cow(mat
);
1192 if (check_col_range(mat
, i
, 1) < 0 ||
1193 check_col_range(mat
, j
, 1) < 0)
1194 return isl_mat_free(mat
);
1196 for (r
= 0; r
< mat
->n_row
; ++r
)
1197 isl_int_swap(mat
->row
[r
][i
], mat
->row
[r
][j
]);
1201 __isl_give isl_mat
*isl_mat_swap_rows(__isl_take isl_mat
*mat
,
1202 unsigned i
, unsigned j
)
1208 mat
= isl_mat_cow(mat
);
1209 if (check_row_range(mat
, i
, 1) < 0 ||
1210 check_row_range(mat
, j
, 1) < 0)
1211 return isl_mat_free(mat
);
1214 mat
->row
[i
] = mat
->row
[j
];
1219 /* Calculate the product of two matrices.
1221 * This function is optimized for operand matrices that contain many zeros and
1222 * skips multiplications where we know one of the operands is zero.
1224 __isl_give isl_mat
*isl_mat_product(__isl_take isl_mat
*left
,
1225 __isl_take isl_mat
*right
)
1228 struct isl_mat
*prod
;
1230 if (!left
|| !right
)
1232 isl_assert(left
->ctx
, left
->n_col
== right
->n_row
, goto error
);
1233 prod
= isl_mat_alloc(left
->ctx
, left
->n_row
, right
->n_col
);
1236 if (left
->n_col
== 0) {
1237 for (i
= 0; i
< prod
->n_row
; ++i
)
1238 isl_seq_clr(prod
->row
[i
], prod
->n_col
);
1240 isl_mat_free(right
);
1243 for (i
= 0; i
< prod
->n_row
; ++i
) {
1244 for (j
= 0; j
< prod
->n_col
; ++j
)
1245 isl_int_mul(prod
->row
[i
][j
],
1246 left
->row
[i
][0], right
->row
[0][j
]);
1247 for (k
= 1; k
< left
->n_col
; ++k
) {
1248 if (isl_int_is_zero(left
->row
[i
][k
]))
1250 for (j
= 0; j
< prod
->n_col
; ++j
)
1251 isl_int_addmul(prod
->row
[i
][j
],
1252 left
->row
[i
][k
], right
->row
[k
][j
]);
1256 isl_mat_free(right
);
1260 isl_mat_free(right
);
1264 /* Replace the variables x in the rows q by x' given by x = M x',
1265 * with M the matrix mat.
1267 * If the number of new variables is greater than the original
1268 * number of variables, then the rows q have already been
1269 * preextended. If the new number is smaller, then the coefficients
1270 * of the divs, which are not changed, need to be shifted down.
1271 * The row q may be the equalities, the inequalities or the
1272 * div expressions. In the latter case, has_div is true and
1273 * we need to take into account the extra denominator column.
1275 static int preimage(struct isl_ctx
*ctx
, isl_int
**q
, unsigned n
,
1276 unsigned n_div
, int has_div
, struct isl_mat
*mat
)
1282 if (mat
->n_col
>= mat
->n_row
)
1285 e
= mat
->n_row
- mat
->n_col
;
1287 for (i
= 0; i
< n
; ++i
)
1288 isl_int_mul(q
[i
][0], q
[i
][0], mat
->row
[0][0]);
1289 t
= isl_mat_sub_alloc6(mat
->ctx
, q
, 0, n
, has_div
, mat
->n_row
);
1290 t
= isl_mat_product(t
, mat
);
1293 for (i
= 0; i
< n
; ++i
) {
1294 isl_seq_swp_or_cpy(q
[i
] + has_div
, t
->row
[i
], t
->n_col
);
1295 isl_seq_cpy(q
[i
] + has_div
+ t
->n_col
,
1296 q
[i
] + has_div
+ t
->n_col
+ e
, n_div
);
1297 isl_seq_clr(q
[i
] + has_div
+ t
->n_col
+ n_div
, e
);
1303 /* Replace the variables x in bset by x' given by x = M x', with
1306 * If there are fewer variables x' then there are x, then we perform
1307 * the transformation in place, which means that, in principle,
1308 * this frees up some extra variables as the number
1309 * of columns remains constant, but we would have to extend
1310 * the div array too as the number of rows in this array is assumed
1311 * to be equal to extra.
1313 __isl_give isl_basic_set
*isl_basic_set_preimage(
1314 __isl_take isl_basic_set
*bset
, __isl_take isl_mat
*mat
)
1316 struct isl_ctx
*ctx
;
1322 bset
= isl_basic_set_cow(bset
);
1326 isl_assert(ctx
, bset
->dim
->nparam
== 0, goto error
);
1327 isl_assert(ctx
, 1+bset
->dim
->n_out
== mat
->n_row
, goto error
);
1328 isl_assert(ctx
, mat
->n_col
> 0, goto error
);
1330 if (mat
->n_col
> mat
->n_row
) {
1331 bset
= isl_basic_set_extend(bset
, 0, mat
->n_col
-1, 0, 0, 0);
1334 } else if (mat
->n_col
< mat
->n_row
) {
1335 bset
->dim
= isl_space_cow(bset
->dim
);
1338 bset
->dim
->n_out
-= mat
->n_row
- mat
->n_col
;
1341 if (preimage(ctx
, bset
->eq
, bset
->n_eq
, bset
->n_div
, 0,
1342 isl_mat_copy(mat
)) < 0)
1345 if (preimage(ctx
, bset
->ineq
, bset
->n_ineq
, bset
->n_div
, 0,
1346 isl_mat_copy(mat
)) < 0)
1349 if (preimage(ctx
, bset
->div
, bset
->n_div
, bset
->n_div
, 1, mat
) < 0)
1352 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
1353 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_REDUNDANT
);
1354 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED
);
1355 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED_DIVS
);
1356 ISL_F_CLR(bset
, ISL_BASIC_SET_ALL_EQUALITIES
);
1358 bset
= isl_basic_set_simplify(bset
);
1359 bset
= isl_basic_set_finalize(bset
);
1365 isl_basic_set_free(bset
);
1369 __isl_give isl_set
*isl_set_preimage(
1370 __isl_take isl_set
*set
, __isl_take isl_mat
*mat
)
1374 set
= isl_set_cow(set
);
1378 for (i
= 0; i
< set
->n
; ++i
) {
1379 set
->p
[i
] = isl_basic_set_preimage(set
->p
[i
],
1384 if (mat
->n_col
!= mat
->n_row
) {
1385 set
->dim
= isl_space_cow(set
->dim
);
1388 set
->dim
->n_out
+= mat
->n_col
;
1389 set
->dim
->n_out
-= mat
->n_row
;
1392 ISL_F_CLR(set
, ISL_SET_NORMALIZED
);
1400 /* Replace the variables x starting at "first_col" in the rows "rows"
1401 * of some coefficient matrix by x' with x = M x' with M the matrix mat.
1402 * That is, replace the corresponding coefficients c by c M.
1404 isl_stat
isl_mat_sub_transform(isl_int
**row
, unsigned n_row
,
1405 unsigned first_col
, __isl_take isl_mat
*mat
)
1412 return isl_stat_error
;
1413 ctx
= isl_mat_get_ctx(mat
);
1414 t
= isl_mat_sub_alloc6(ctx
, row
, 0, n_row
, first_col
, mat
->n_row
);
1415 t
= isl_mat_product(t
, mat
);
1417 return isl_stat_error
;
1418 for (i
= 0; i
< n_row
; ++i
)
1419 isl_seq_swp_or_cpy(row
[i
] + first_col
, t
->row
[i
], t
->n_col
);
1424 void isl_mat_print_internal(__isl_keep isl_mat
*mat
, FILE *out
, int indent
)
1429 fprintf(out
, "%*snull mat\n", indent
, "");
1433 if (mat
->n_row
== 0)
1434 fprintf(out
, "%*s[]\n", indent
, "");
1436 for (i
= 0; i
< mat
->n_row
; ++i
) {
1438 fprintf(out
, "%*s[[", indent
, "");
1440 fprintf(out
, "%*s[", indent
+1, "");
1441 for (j
= 0; j
< mat
->n_col
; ++j
) {
1444 isl_int_print(out
, mat
->row
[i
][j
], 0);
1446 if (i
== mat
->n_row
-1)
1447 fprintf(out
, "]]\n");
1449 fprintf(out
, "]\n");
1453 void isl_mat_dump(__isl_keep isl_mat
*mat
)
1455 isl_mat_print_internal(mat
, stderr
, 0);
1458 __isl_give isl_mat
*isl_mat_drop_cols(__isl_take isl_mat
*mat
,
1459 unsigned col
, unsigned n
)
1466 mat
= isl_mat_cow(mat
);
1467 if (check_col_range(mat
, col
, n
) < 0)
1468 return isl_mat_free(mat
);
1470 if (col
!= mat
->n_col
-n
) {
1471 for (r
= 0; r
< mat
->n_row
; ++r
)
1472 isl_seq_cpy(mat
->row
[r
]+col
, mat
->row
[r
]+col
+n
,
1473 mat
->n_col
- col
- n
);
1479 __isl_give isl_mat
*isl_mat_drop_rows(__isl_take isl_mat
*mat
,
1480 unsigned row
, unsigned n
)
1484 mat
= isl_mat_cow(mat
);
1485 if (check_row_range(mat
, row
, n
) < 0)
1486 return isl_mat_free(mat
);
1488 for (r
= row
; r
+n
< mat
->n_row
; ++r
)
1489 mat
->row
[r
] = mat
->row
[r
+n
];
1495 __isl_give isl_mat
*isl_mat_insert_cols(__isl_take isl_mat
*mat
,
1496 unsigned col
, unsigned n
)
1500 if (check_col_range(mat
, col
, 0) < 0)
1501 return isl_mat_free(mat
);
1505 ext
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
+ n
);
1509 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
, 0, 0, col
);
1510 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
,
1511 col
+ n
, col
, mat
->n_col
- col
);
1520 __isl_give isl_mat
*isl_mat_insert_zero_cols(__isl_take isl_mat
*mat
,
1521 unsigned first
, unsigned n
)
1527 mat
= isl_mat_insert_cols(mat
, first
, n
);
1531 for (i
= 0; i
< mat
->n_row
; ++i
)
1532 isl_seq_clr(mat
->row
[i
] + first
, n
);
1537 __isl_give isl_mat
*isl_mat_add_zero_cols(__isl_take isl_mat
*mat
, unsigned n
)
1542 return isl_mat_insert_zero_cols(mat
, mat
->n_col
, n
);
1545 __isl_give isl_mat
*isl_mat_insert_rows(__isl_take isl_mat
*mat
,
1546 unsigned row
, unsigned n
)
1550 if (check_row_range(mat
, row
, 0) < 0)
1551 return isl_mat_free(mat
);
1555 ext
= isl_mat_alloc(mat
->ctx
, mat
->n_row
+ n
, mat
->n_col
);
1559 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, row
, 0, 0, mat
->n_col
);
1560 isl_mat_sub_copy(mat
->ctx
, ext
->row
+ row
+ n
, mat
->row
+ row
,
1561 mat
->n_row
- row
, 0, 0, mat
->n_col
);
1570 __isl_give isl_mat
*isl_mat_add_rows(__isl_take isl_mat
*mat
, unsigned n
)
1575 return isl_mat_insert_rows(mat
, mat
->n_row
, n
);
1578 __isl_give isl_mat
*isl_mat_insert_zero_rows(__isl_take isl_mat
*mat
,
1579 unsigned row
, unsigned n
)
1583 mat
= isl_mat_insert_rows(mat
, row
, n
);
1587 for (i
= 0; i
< n
; ++i
)
1588 isl_seq_clr(mat
->row
[row
+ i
], mat
->n_col
);
1593 __isl_give isl_mat
*isl_mat_add_zero_rows(__isl_take isl_mat
*mat
, unsigned n
)
1598 return isl_mat_insert_zero_rows(mat
, mat
->n_row
, n
);
1601 void isl_mat_col_submul(struct isl_mat
*mat
,
1602 int dst_col
, isl_int f
, int src_col
)
1606 for (i
= 0; i
< mat
->n_row
; ++i
)
1607 isl_int_submul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1610 void isl_mat_col_add(__isl_keep isl_mat
*mat
, int dst_col
, int src_col
)
1617 for (i
= 0; i
< mat
->n_row
; ++i
)
1618 isl_int_add(mat
->row
[i
][dst_col
],
1619 mat
->row
[i
][dst_col
], mat
->row
[i
][src_col
]);
1622 void isl_mat_col_mul(struct isl_mat
*mat
, int dst_col
, isl_int f
, int src_col
)
1626 for (i
= 0; i
< mat
->n_row
; ++i
)
1627 isl_int_mul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1630 /* Add "f" times column "src_col" to column "dst_col" of "mat" and
1631 * return the result.
1633 __isl_give isl_mat
*isl_mat_col_addmul(__isl_take isl_mat
*mat
, int dst_col
,
1634 isl_int f
, int src_col
)
1638 if (check_col(mat
, dst_col
) < 0 || check_col(mat
, src_col
) < 0)
1639 return isl_mat_free(mat
);
1641 for (i
= 0; i
< mat
->n_row
; ++i
) {
1642 if (isl_int_is_zero(mat
->row
[i
][src_col
]))
1644 mat
= isl_mat_cow(mat
);
1647 isl_int_addmul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1653 /* Negate column "col" of "mat" and return the result.
1655 __isl_give isl_mat
*isl_mat_col_neg(__isl_take isl_mat
*mat
, int col
)
1659 if (check_col(mat
, col
) < 0)
1660 return isl_mat_free(mat
);
1662 for (i
= 0; i
< mat
->n_row
; ++i
) {
1663 if (isl_int_is_zero(mat
->row
[i
][col
]))
1665 mat
= isl_mat_cow(mat
);
1668 isl_int_neg(mat
->row
[i
][col
], mat
->row
[i
][col
]);
1674 /* Negate row "row" of "mat" and return the result.
1676 __isl_give isl_mat
*isl_mat_row_neg(__isl_take isl_mat
*mat
, int row
)
1678 if (check_row(mat
, row
) < 0)
1679 return isl_mat_free(mat
);
1680 if (isl_seq_first_non_zero(mat
->row
[row
], mat
->n_col
) == -1)
1682 mat
= isl_mat_cow(mat
);
1685 isl_seq_neg(mat
->row
[row
], mat
->row
[row
], mat
->n_col
);
1689 __isl_give isl_mat
*isl_mat_unimodular_complete(__isl_take isl_mat
*M
, int row
)
1692 struct isl_mat
*H
= NULL
, *Q
= NULL
;
1697 isl_assert(M
->ctx
, M
->n_row
== M
->n_col
, goto error
);
1699 H
= isl_mat_left_hermite(isl_mat_copy(M
), 0, NULL
, &Q
);
1700 M
->n_row
= M
->n_col
;
1703 for (r
= 0; r
< row
; ++r
)
1704 isl_assert(M
->ctx
, isl_int_is_one(H
->row
[r
][r
]), goto error
);
1705 for (r
= row
; r
< M
->n_row
; ++r
)
1706 isl_seq_cpy(M
->row
[r
], Q
->row
[r
], M
->n_col
);
1717 __isl_give isl_mat
*isl_mat_concat(__isl_take isl_mat
*top
,
1718 __isl_take isl_mat
*bot
)
1720 struct isl_mat
*mat
;
1725 isl_assert(top
->ctx
, top
->n_col
== bot
->n_col
, goto error
);
1726 if (top
->n_row
== 0) {
1730 if (bot
->n_row
== 0) {
1735 mat
= isl_mat_alloc(top
->ctx
, top
->n_row
+ bot
->n_row
, top
->n_col
);
1738 isl_mat_sub_copy(mat
->ctx
, mat
->row
, top
->row
, top
->n_row
,
1740 isl_mat_sub_copy(mat
->ctx
, mat
->row
+ top
->n_row
, bot
->row
, bot
->n_row
,
1751 isl_bool
isl_mat_is_equal(__isl_keep isl_mat
*mat1
, __isl_keep isl_mat
*mat2
)
1756 return isl_bool_error
;
1758 if (mat1
->n_row
!= mat2
->n_row
)
1759 return isl_bool_false
;
1761 if (mat1
->n_col
!= mat2
->n_col
)
1762 return isl_bool_false
;
1764 for (i
= 0; i
< mat1
->n_row
; ++i
)
1765 if (!isl_seq_eq(mat1
->row
[i
], mat2
->row
[i
], mat1
->n_col
))
1766 return isl_bool_false
;
1768 return isl_bool_true
;
1771 __isl_give isl_mat
*isl_mat_from_row_vec(__isl_take isl_vec
*vec
)
1773 struct isl_mat
*mat
;
1777 mat
= isl_mat_alloc(vec
->ctx
, 1, vec
->size
);
1781 isl_seq_cpy(mat
->row
[0], vec
->el
, vec
->size
);
1790 /* Return a copy of row "row" of "mat" as an isl_vec.
1792 __isl_give isl_vec
*isl_mat_get_row(__isl_keep isl_mat
*mat
, unsigned row
)
1798 if (row
>= mat
->n_row
)
1799 isl_die(mat
->ctx
, isl_error_invalid
, "row out of range",
1802 v
= isl_vec_alloc(isl_mat_get_ctx(mat
), mat
->n_col
);
1805 isl_seq_cpy(v
->el
, mat
->row
[row
], mat
->n_col
);
1810 __isl_give isl_mat
*isl_mat_vec_concat(__isl_take isl_mat
*top
,
1811 __isl_take isl_vec
*bot
)
1813 return isl_mat_concat(top
, isl_mat_from_row_vec(bot
));
1816 __isl_give isl_mat
*isl_mat_move_cols(__isl_take isl_mat
*mat
,
1817 unsigned dst_col
, unsigned src_col
, unsigned n
)
1823 if (n
== 0 || dst_col
== src_col
)
1826 res
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
1830 if (dst_col
< src_col
) {
1831 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1833 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1834 dst_col
, src_col
, n
);
1835 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1836 dst_col
+ n
, dst_col
, src_col
- dst_col
);
1837 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1838 src_col
+ n
, src_col
+ n
,
1839 res
->n_col
- src_col
- n
);
1841 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1843 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1844 src_col
, src_col
+ n
, dst_col
- src_col
);
1845 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1846 dst_col
, src_col
, n
);
1847 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1848 dst_col
+ n
, dst_col
+ n
,
1849 res
->n_col
- dst_col
- n
);
1859 /* Return the gcd of the elements in row "row" of "mat" in *gcd.
1860 * Return isl_stat_ok on success and isl_stat_error on failure.
1862 isl_stat
isl_mat_row_gcd(__isl_keep isl_mat
*mat
, int row
, isl_int
*gcd
)
1864 if (check_row(mat
, row
) < 0)
1865 return isl_stat_error
;
1867 isl_seq_gcd(mat
->row
[row
], mat
->n_col
, gcd
);
1872 void isl_mat_gcd(__isl_keep isl_mat
*mat
, isl_int
*gcd
)
1877 isl_int_set_si(*gcd
, 0);
1882 for (i
= 0; i
< mat
->n_row
; ++i
) {
1883 isl_seq_gcd(mat
->row
[i
], mat
->n_col
, &g
);
1884 isl_int_gcd(*gcd
, *gcd
, g
);
1889 /* Return the result of scaling "mat" by a factor of "m".
1891 __isl_give isl_mat
*isl_mat_scale(__isl_take isl_mat
*mat
, isl_int m
)
1895 if (isl_int_is_one(m
))
1898 mat
= isl_mat_cow(mat
);
1902 for (i
= 0; i
< mat
->n_row
; ++i
)
1903 isl_seq_scale(mat
->row
[i
], mat
->row
[i
], m
, mat
->n_col
);
1908 __isl_give isl_mat
*isl_mat_scale_down(__isl_take isl_mat
*mat
, isl_int m
)
1912 if (isl_int_is_one(m
))
1915 mat
= isl_mat_cow(mat
);
1919 for (i
= 0; i
< mat
->n_row
; ++i
)
1920 isl_seq_scale_down(mat
->row
[i
], mat
->row
[i
], m
, mat
->n_col
);
1925 __isl_give isl_mat
*isl_mat_scale_down_row(__isl_take isl_mat
*mat
, int row
,
1928 if (isl_int_is_one(m
))
1931 mat
= isl_mat_cow(mat
);
1935 isl_seq_scale_down(mat
->row
[row
], mat
->row
[row
], m
, mat
->n_col
);
1940 __isl_give isl_mat
*isl_mat_normalize(__isl_take isl_mat
*mat
)
1948 isl_mat_gcd(mat
, &gcd
);
1949 mat
= isl_mat_scale_down(mat
, gcd
);
1955 __isl_give isl_mat
*isl_mat_normalize_row(__isl_take isl_mat
*mat
, int row
)
1957 mat
= isl_mat_cow(mat
);
1961 isl_seq_normalize(mat
->ctx
, mat
->row
[row
], mat
->n_col
);
1966 /* Number of initial non-zero columns.
1968 int isl_mat_initial_non_zero_cols(__isl_keep isl_mat
*mat
)
1975 for (i
= 0; i
< mat
->n_col
; ++i
)
1976 if (row_first_non_zero(mat
->row
, mat
->n_row
, i
) < 0)