2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2010 INRIA Saclay
4 * Copyright 2014 Ecole Normale Superieure
5 * Copyright 2017 Sven Verdoolaege
7 * Use of this software is governed by the MIT license
9 * Written by Sven Verdoolaege, K.U.Leuven, Departement
10 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
11 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
12 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
13 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
16 #include <isl_ctx_private.h>
17 #include <isl_map_private.h>
18 #include <isl/space.h>
20 #include <isl_mat_private.h>
21 #include <isl_vec_private.h>
22 #include <isl_space_private.h>
23 #include <isl_val_private.h>
24 #include <isl/deprecated/mat_int.h>
26 isl_ctx
*isl_mat_get_ctx(__isl_keep isl_mat
*mat
)
28 return mat
? mat
->ctx
: NULL
;
31 /* Return a hash value that digests "mat".
33 uint32_t isl_mat_get_hash(__isl_keep isl_mat
*mat
)
41 hash
= isl_hash_init();
42 isl_hash_byte(hash
, mat
->n_row
& 0xFF);
43 isl_hash_byte(hash
, mat
->n_col
& 0xFF);
44 for (i
= 0; i
< mat
->n_row
; ++i
) {
47 row_hash
= isl_seq_get_hash(mat
->row
[i
], mat
->n_col
);
48 isl_hash_hash(hash
, row_hash
);
54 struct isl_mat
*isl_mat_alloc(struct isl_ctx
*ctx
,
55 unsigned n_row
, unsigned n_col
)
60 mat
= isl_alloc_type(ctx
, struct isl_mat
);
65 mat
->block
= isl_blk_alloc(ctx
, n_row
* n_col
);
66 if (isl_blk_is_error(mat
->block
))
68 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
69 if (n_row
&& !mat
->row
)
72 for (i
= 0; i
< n_row
; ++i
)
73 mat
->row
[i
] = mat
->block
.data
+ i
* n_col
;
85 isl_blk_free(ctx
, mat
->block
);
90 struct isl_mat
*isl_mat_extend(struct isl_mat
*mat
,
91 unsigned n_row
, unsigned n_col
)
100 if (mat
->max_col
>= n_col
&& mat
->n_row
>= n_row
) {
101 if (mat
->n_col
< n_col
)
106 if (mat
->max_col
< n_col
) {
107 struct isl_mat
*new_mat
;
109 if (n_row
< mat
->n_row
)
111 new_mat
= isl_mat_alloc(mat
->ctx
, n_row
, n_col
);
114 for (i
= 0; i
< mat
->n_row
; ++i
)
115 isl_seq_cpy(new_mat
->row
[i
], mat
->row
[i
], mat
->n_col
);
120 mat
= isl_mat_cow(mat
);
124 old
= mat
->block
.data
;
125 mat
->block
= isl_blk_extend(mat
->ctx
, mat
->block
, n_row
* mat
->max_col
);
126 if (isl_blk_is_error(mat
->block
))
128 row
= isl_realloc_array(mat
->ctx
, mat
->row
, isl_int
*, n_row
);
133 for (i
= 0; i
< mat
->n_row
; ++i
)
134 mat
->row
[i
] = mat
->block
.data
+ (mat
->row
[i
] - old
);
135 for (i
= mat
->n_row
; i
< n_row
; ++i
)
136 mat
->row
[i
] = mat
->block
.data
+ i
* mat
->max_col
;
138 if (mat
->n_col
< n_col
)
147 __isl_give isl_mat
*isl_mat_sub_alloc6(isl_ctx
*ctx
, isl_int
**row
,
148 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
153 mat
= isl_alloc_type(ctx
, struct isl_mat
);
156 mat
->row
= isl_alloc_array(ctx
, isl_int
*, n_row
);
157 if (n_row
&& !mat
->row
)
159 for (i
= 0; i
< n_row
; ++i
)
160 mat
->row
[i
] = row
[first_row
+i
] + first_col
;
166 mat
->block
= isl_blk_empty();
167 mat
->flags
= ISL_MAT_BORROWED
;
174 __isl_give isl_mat
*isl_mat_sub_alloc(__isl_keep isl_mat
*mat
,
175 unsigned first_row
, unsigned n_row
, unsigned first_col
, unsigned n_col
)
179 return isl_mat_sub_alloc6(mat
->ctx
, mat
->row
, first_row
, n_row
,
183 void isl_mat_sub_copy(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
184 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
188 for (i
= 0; i
< n_row
; ++i
)
189 isl_seq_cpy(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
192 void isl_mat_sub_neg(struct isl_ctx
*ctx
, isl_int
**dst
, isl_int
**src
,
193 unsigned n_row
, unsigned dst_col
, unsigned src_col
, unsigned n_col
)
197 for (i
= 0; i
< n_row
; ++i
)
198 isl_seq_neg(dst
[i
]+dst_col
, src
[i
]+src_col
, n_col
);
201 __isl_give isl_mat
*isl_mat_copy(__isl_keep isl_mat
*mat
)
210 __isl_give isl_mat
*isl_mat_dup(__isl_keep isl_mat
*mat
)
213 struct isl_mat
*mat2
;
217 mat2
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
220 for (i
= 0; i
< mat
->n_row
; ++i
)
221 isl_seq_cpy(mat2
->row
[i
], mat
->row
[i
], mat
->n_col
);
225 __isl_give isl_mat
*isl_mat_cow(__isl_take isl_mat
*mat
)
227 struct isl_mat
*mat2
;
231 if (mat
->ref
== 1 && !ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
234 mat2
= isl_mat_dup(mat
);
239 __isl_null isl_mat
*isl_mat_free(__isl_take isl_mat
*mat
)
247 if (!ISL_F_ISSET(mat
, ISL_MAT_BORROWED
))
248 isl_blk_free(mat
->ctx
, mat
->block
);
249 isl_ctx_deref(mat
->ctx
);
256 int isl_mat_rows(__isl_keep isl_mat
*mat
)
258 return mat
? mat
->n_row
: -1;
261 int isl_mat_cols(__isl_keep isl_mat
*mat
)
263 return mat
? mat
->n_col
: -1;
266 /* Check that "col" is a valid column position for "mat".
268 static isl_stat
check_col(__isl_keep isl_mat
*mat
, int col
)
271 return isl_stat_error
;
272 if (col
< 0 || col
>= mat
->n_col
)
273 isl_die(isl_mat_get_ctx(mat
), isl_error_invalid
,
274 "column out of range", return isl_stat_error
);
278 /* Check that "row" is a valid row position for "mat".
280 static isl_stat
check_row(__isl_keep isl_mat
*mat
, int row
)
283 return isl_stat_error
;
284 if (row
< 0 || row
>= mat
->n_row
)
285 isl_die(isl_mat_get_ctx(mat
), isl_error_invalid
,
286 "row out of range", return isl_stat_error
);
290 /* Check that there are "n" columns starting at position "first" in "mat".
292 static isl_stat
check_col_range(__isl_keep isl_mat
*mat
, unsigned first
,
296 return isl_stat_error
;
297 if (first
+ n
> mat
->n_col
|| first
+ n
< first
)
298 isl_die(isl_mat_get_ctx(mat
), isl_error_invalid
,
299 "column position or range out of bounds",
300 return isl_stat_error
);
304 /* Check that there are "n" rows starting at position "first" in "mat".
306 static isl_stat
check_row_range(__isl_keep isl_mat
*mat
, unsigned first
,
310 return isl_stat_error
;
311 if (first
+ n
> mat
->n_row
|| first
+ n
< first
)
312 isl_die(isl_mat_get_ctx(mat
), isl_error_invalid
,
313 "row position or range out of bounds",
314 return isl_stat_error
);
318 int isl_mat_get_element(__isl_keep isl_mat
*mat
, int row
, int col
, isl_int
*v
)
320 if (check_row(mat
, row
) < 0)
322 if (check_col(mat
, col
) < 0)
324 isl_int_set(*v
, mat
->row
[row
][col
]);
328 /* Extract the element at row "row", oolumn "col" of "mat".
330 __isl_give isl_val
*isl_mat_get_element_val(__isl_keep isl_mat
*mat
,
335 if (check_row(mat
, row
) < 0)
337 if (check_col(mat
, col
) < 0)
339 ctx
= isl_mat_get_ctx(mat
);
340 return isl_val_int_from_isl_int(ctx
, mat
->row
[row
][col
]);
343 __isl_give isl_mat
*isl_mat_set_element(__isl_take isl_mat
*mat
,
344 int row
, int col
, isl_int v
)
346 mat
= isl_mat_cow(mat
);
347 if (check_row(mat
, row
) < 0)
348 return isl_mat_free(mat
);
349 if (check_col(mat
, col
) < 0)
350 return isl_mat_free(mat
);
351 isl_int_set(mat
->row
[row
][col
], v
);
355 __isl_give isl_mat
*isl_mat_set_element_si(__isl_take isl_mat
*mat
,
356 int row
, int col
, int v
)
358 mat
= isl_mat_cow(mat
);
359 if (check_row(mat
, row
) < 0)
360 return isl_mat_free(mat
);
361 if (check_col(mat
, col
) < 0)
362 return isl_mat_free(mat
);
363 isl_int_set_si(mat
->row
[row
][col
], v
);
367 /* Replace the element at row "row", column "col" of "mat" by "v".
369 __isl_give isl_mat
*isl_mat_set_element_val(__isl_take isl_mat
*mat
,
370 int row
, int col
, __isl_take isl_val
*v
)
373 return isl_mat_free(mat
);
374 if (!isl_val_is_int(v
))
375 isl_die(isl_val_get_ctx(v
), isl_error_invalid
,
376 "expecting integer value", goto error
);
377 mat
= isl_mat_set_element(mat
, row
, col
, v
->n
);
382 return isl_mat_free(mat
);
385 __isl_give isl_mat
*isl_mat_diag(isl_ctx
*ctx
, unsigned n_row
, isl_int d
)
390 mat
= isl_mat_alloc(ctx
, n_row
, n_row
);
393 for (i
= 0; i
< n_row
; ++i
) {
394 isl_seq_clr(mat
->row
[i
], i
);
395 isl_int_set(mat
->row
[i
][i
], d
);
396 isl_seq_clr(mat
->row
[i
]+i
+1, n_row
-(i
+1));
402 /* Create an "n_row" by "n_col" matrix with zero elements.
404 __isl_give isl_mat
*isl_mat_zero(isl_ctx
*ctx
, unsigned n_row
, unsigned n_col
)
409 mat
= isl_mat_alloc(ctx
, n_row
, n_col
);
412 for (i
= 0; i
< n_row
; ++i
)
413 isl_seq_clr(mat
->row
[i
], n_col
);
418 __isl_give isl_mat
*isl_mat_identity(isl_ctx
*ctx
, unsigned n_row
)
422 return isl_mat_diag(ctx
, n_row
, ctx
->one
);
425 /* Is "mat" a (possibly scaled) identity matrix?
427 int isl_mat_is_scaled_identity(__isl_keep isl_mat
*mat
)
433 if (mat
->n_row
!= mat
->n_col
)
436 for (i
= 0; i
< mat
->n_row
; ++i
) {
437 if (isl_seq_first_non_zero(mat
->row
[i
], i
) != -1)
439 if (isl_int_ne(mat
->row
[0][0], mat
->row
[i
][i
]))
441 if (isl_seq_first_non_zero(mat
->row
[i
] + i
+ 1,
442 mat
->n_col
- (i
+ 1)) != -1)
449 __isl_give isl_vec
*isl_mat_vec_product(__isl_take isl_mat
*mat
,
450 __isl_take isl_vec
*vec
)
453 struct isl_vec
*prod
;
458 isl_assert(mat
->ctx
, mat
->n_col
== vec
->size
, goto error
);
460 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_row
);
464 for (i
= 0; i
< prod
->size
; ++i
)
465 isl_seq_inner_product(mat
->row
[i
], vec
->el
, vec
->size
,
466 &prod
->block
.data
[i
]);
476 __isl_give isl_vec
*isl_mat_vec_inverse_product(__isl_take isl_mat
*mat
,
477 __isl_take isl_vec
*vec
)
479 struct isl_mat
*vec_mat
;
484 vec_mat
= isl_mat_alloc(vec
->ctx
, vec
->size
, 1);
487 for (i
= 0; i
< vec
->size
; ++i
)
488 isl_int_set(vec_mat
->row
[i
][0], vec
->el
[i
]);
489 vec_mat
= isl_mat_inverse_product(mat
, vec_mat
);
493 vec
= isl_vec_alloc(vec_mat
->ctx
, vec_mat
->n_row
);
495 for (i
= 0; i
< vec
->size
; ++i
)
496 isl_int_set(vec
->el
[i
], vec_mat
->row
[i
][0]);
497 isl_mat_free(vec_mat
);
505 __isl_give isl_vec
*isl_vec_mat_product(__isl_take isl_vec
*vec
,
506 __isl_take isl_mat
*mat
)
509 struct isl_vec
*prod
;
514 isl_assert(mat
->ctx
, mat
->n_row
== vec
->size
, goto error
);
516 prod
= isl_vec_alloc(mat
->ctx
, mat
->n_col
);
520 for (i
= 0; i
< prod
->size
; ++i
) {
521 isl_int_set_si(prod
->el
[i
], 0);
522 for (j
= 0; j
< vec
->size
; ++j
)
523 isl_int_addmul(prod
->el
[i
], vec
->el
[j
], mat
->row
[j
][i
]);
534 __isl_give isl_mat
*isl_mat_aff_direct_sum(__isl_take isl_mat
*left
,
535 __isl_take isl_mat
*right
)
543 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
544 isl_assert(left
->ctx
, left
->n_row
>= 1, goto error
);
545 isl_assert(left
->ctx
, left
->n_col
>= 1, goto error
);
546 isl_assert(left
->ctx
, right
->n_col
>= 1, goto error
);
547 isl_assert(left
->ctx
,
548 isl_seq_first_non_zero(left
->row
[0]+1, left
->n_col
-1) == -1,
550 isl_assert(left
->ctx
,
551 isl_seq_first_non_zero(right
->row
[0]+1, right
->n_col
-1) == -1,
554 sum
= isl_mat_alloc(left
->ctx
, left
->n_row
, left
->n_col
+ right
->n_col
- 1);
557 isl_int_lcm(sum
->row
[0][0], left
->row
[0][0], right
->row
[0][0]);
558 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
559 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
561 isl_seq_clr(sum
->row
[0]+1, sum
->n_col
-1);
562 for (i
= 1; i
< sum
->n_row
; ++i
) {
563 isl_int_mul(sum
->row
[i
][0], left
->row
[0][0], left
->row
[i
][0]);
564 isl_int_addmul(sum
->row
[i
][0],
565 right
->row
[0][0], right
->row
[i
][0]);
566 isl_seq_scale(sum
->row
[i
]+1, left
->row
[i
]+1, left
->row
[0][0],
568 isl_seq_scale(sum
->row
[i
]+left
->n_col
,
569 right
->row
[i
]+1, right
->row
[0][0],
573 isl_int_divexact(left
->row
[0][0], sum
->row
[0][0], left
->row
[0][0]);
574 isl_int_divexact(right
->row
[0][0], sum
->row
[0][0], right
->row
[0][0]);
584 static void exchange(struct isl_mat
*M
, struct isl_mat
**U
,
585 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
)
588 for (r
= row
; r
< M
->n_row
; ++r
)
589 isl_int_swap(M
->row
[r
][i
], M
->row
[r
][j
]);
591 for (r
= 0; r
< (*U
)->n_row
; ++r
)
592 isl_int_swap((*U
)->row
[r
][i
], (*U
)->row
[r
][j
]);
595 isl_mat_swap_rows(*Q
, i
, j
);
598 static void subtract(struct isl_mat
*M
, struct isl_mat
**U
,
599 struct isl_mat
**Q
, unsigned row
, unsigned i
, unsigned j
, isl_int m
)
602 for (r
= row
; r
< M
->n_row
; ++r
)
603 isl_int_submul(M
->row
[r
][j
], m
, M
->row
[r
][i
]);
605 for (r
= 0; r
< (*U
)->n_row
; ++r
)
606 isl_int_submul((*U
)->row
[r
][j
], m
, (*U
)->row
[r
][i
]);
609 for (r
= 0; r
< (*Q
)->n_col
; ++r
)
610 isl_int_addmul((*Q
)->row
[i
][r
], m
, (*Q
)->row
[j
][r
]);
614 static void oppose(struct isl_mat
*M
, struct isl_mat
**U
,
615 struct isl_mat
**Q
, unsigned row
, unsigned col
)
618 for (r
= row
; r
< M
->n_row
; ++r
)
619 isl_int_neg(M
->row
[r
][col
], M
->row
[r
][col
]);
621 for (r
= 0; r
< (*U
)->n_row
; ++r
)
622 isl_int_neg((*U
)->row
[r
][col
], (*U
)->row
[r
][col
]);
625 isl_seq_neg((*Q
)->row
[col
], (*Q
)->row
[col
], (*Q
)->n_col
);
628 /* Given matrix M, compute
633 * with U and Q unimodular matrices and H a matrix in column echelon form
634 * such that on each echelon row the entries in the non-echelon column
635 * are non-negative (if neg == 0) or non-positive (if neg == 1)
636 * and strictly smaller (in absolute value) than the entries in the echelon
638 * If U or Q are NULL, then these matrices are not computed.
640 __isl_give isl_mat
*isl_mat_left_hermite(__isl_take isl_mat
*M
, int neg
,
641 __isl_give isl_mat
**U
, __isl_give isl_mat
**Q
)
656 *U
= isl_mat_identity(M
->ctx
, M
->n_col
);
661 *Q
= isl_mat_identity(M
->ctx
, M
->n_col
);
668 for (row
= 0; row
< M
->n_row
; ++row
) {
670 first
= isl_seq_abs_min_non_zero(M
->row
[row
]+col
, M
->n_col
-col
);
675 exchange(M
, U
, Q
, row
, first
, col
);
676 if (isl_int_is_neg(M
->row
[row
][col
]))
677 oppose(M
, U
, Q
, row
, col
);
679 while ((off
= isl_seq_first_non_zero(M
->row
[row
]+first
,
680 M
->n_col
-first
)) != -1) {
682 isl_int_fdiv_q(c
, M
->row
[row
][first
], M
->row
[row
][col
]);
683 subtract(M
, U
, Q
, row
, col
, first
, c
);
684 if (!isl_int_is_zero(M
->row
[row
][first
]))
685 exchange(M
, U
, Q
, row
, first
, col
);
689 for (i
= 0; i
< col
; ++i
) {
690 if (isl_int_is_zero(M
->row
[row
][i
]))
693 isl_int_cdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
695 isl_int_fdiv_q(c
, M
->row
[row
][i
], M
->row
[row
][col
]);
696 if (isl_int_is_zero(c
))
698 subtract(M
, U
, Q
, row
, col
, i
, c
);
718 /* Use row "row" of "mat" to eliminate column "col" from all other rows.
720 static __isl_give isl_mat
*eliminate(__isl_take isl_mat
*mat
, int row
, int col
)
728 ctx
= isl_mat_get_ctx(mat
);
729 nr
= isl_mat_rows(mat
);
730 nc
= isl_mat_cols(mat
);
732 for (k
= 0; k
< nr
; ++k
) {
735 if (isl_int_is_zero(mat
->row
[k
][col
]))
737 mat
= isl_mat_cow(mat
);
740 isl_seq_elim(mat
->row
[k
], mat
->row
[row
], col
, nc
, NULL
);
741 isl_seq_normalize(ctx
, mat
->row
[k
], nc
);
747 /* Perform Gaussian elimination on the rows of "mat", but start
748 * from the final row and the final column.
749 * Any zero rows that result from the elimination are removed.
751 * In particular, for each column from last to first,
752 * look for the last row with a non-zero coefficient in that column,
753 * move it last (but before other rows moved last in previous steps) and
754 * use it to eliminate the column from the other rows.
756 __isl_give isl_mat
*isl_mat_reverse_gauss(__isl_take isl_mat
*mat
)
758 int k
, row
, last
, nr
, nc
;
763 nr
= isl_mat_rows(mat
);
764 nc
= isl_mat_cols(mat
);
767 for (row
= nr
- 1; row
>= 0; --row
) {
768 for (; last
>= 0; --last
) {
769 for (k
= row
; k
>= 0; --k
)
770 if (!isl_int_is_zero(mat
->row
[k
][last
]))
778 mat
= isl_mat_swap_rows(mat
, k
, row
);
781 if (isl_int_is_neg(mat
->row
[row
][last
]))
782 mat
= isl_mat_row_neg(mat
, row
);
783 mat
= eliminate(mat
, row
, last
);
787 mat
= isl_mat_drop_rows(mat
, 0, row
+ 1);
792 /* Negate the lexicographically negative rows of "mat" such that
793 * all rows in the result are lexicographically non-negative.
795 __isl_give isl_mat
*isl_mat_lexnonneg_rows(__isl_take isl_mat
*mat
)
802 nr
= isl_mat_rows(mat
);
803 nc
= isl_mat_cols(mat
);
805 for (i
= 0; i
< nr
; ++i
) {
808 pos
= isl_seq_first_non_zero(mat
->row
[i
], nc
);
811 if (isl_int_is_nonneg(mat
->row
[i
][pos
]))
813 mat
= isl_mat_row_neg(mat
, i
);
821 /* Given a matrix "H" is column echelon form, what is the first
822 * zero column? That is how many initial columns are non-zero?
823 * Start looking at column "first_col" and only consider
824 * the columns to be of size "n_row".
825 * "H" is assumed to be non-NULL.
827 * Since "H" is in column echelon form, the first non-zero entry
828 * in a column is always in a later position compared to the previous column.
830 static int hermite_first_zero_col(__isl_keep isl_mat
*H
, int first_col
,
835 for (col
= first_col
, row
= 0; col
< H
->n_col
; ++col
) {
836 for (; row
< n_row
; ++row
)
837 if (!isl_int_is_zero(H
->row
[row
][col
]))
846 /* Return the rank of "mat", or -1 in case of error.
848 int isl_mat_rank(__isl_keep isl_mat
*mat
)
853 H
= isl_mat_left_hermite(isl_mat_copy(mat
), 0, NULL
, NULL
);
857 rank
= hermite_first_zero_col(H
, 0, H
->n_row
);
863 __isl_give isl_mat
*isl_mat_right_kernel(__isl_take isl_mat
*mat
)
866 struct isl_mat
*U
= NULL
;
869 mat
= isl_mat_left_hermite(mat
, 0, &U
, NULL
);
873 rank
= hermite_first_zero_col(mat
, 0, mat
->n_row
);
874 K
= isl_mat_alloc(U
->ctx
, U
->n_row
, U
->n_col
- rank
);
877 isl_mat_sub_copy(K
->ctx
, K
->row
, U
->row
, U
->n_row
, 0, rank
, U
->n_col
-rank
);
887 __isl_give isl_mat
*isl_mat_lin_to_aff(__isl_take isl_mat
*mat
)
890 struct isl_mat
*mat2
;
894 mat2
= isl_mat_alloc(mat
->ctx
, 1+mat
->n_row
, 1+mat
->n_col
);
897 isl_int_set_si(mat2
->row
[0][0], 1);
898 isl_seq_clr(mat2
->row
[0]+1, mat
->n_col
);
899 for (i
= 0; i
< mat
->n_row
; ++i
) {
900 isl_int_set_si(mat2
->row
[1+i
][0], 0);
901 isl_seq_cpy(mat2
->row
[1+i
]+1, mat
->row
[i
], mat
->n_col
);
910 /* Given two matrices M1 and M2, return the block matrix
915 __isl_give isl_mat
*isl_mat_diagonal(__isl_take isl_mat
*mat1
,
916 __isl_take isl_mat
*mat2
)
924 mat
= isl_mat_alloc(mat1
->ctx
, mat1
->n_row
+ mat2
->n_row
,
925 mat1
->n_col
+ mat2
->n_col
);
928 for (i
= 0; i
< mat1
->n_row
; ++i
) {
929 isl_seq_cpy(mat
->row
[i
], mat1
->row
[i
], mat1
->n_col
);
930 isl_seq_clr(mat
->row
[i
] + mat1
->n_col
, mat2
->n_col
);
932 for (i
= 0; i
< mat2
->n_row
; ++i
) {
933 isl_seq_clr(mat
->row
[mat1
->n_row
+ i
], mat1
->n_col
);
934 isl_seq_cpy(mat
->row
[mat1
->n_row
+ i
] + mat1
->n_col
,
935 mat2
->row
[i
], mat2
->n_col
);
946 static int row_first_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
950 for (i
= 0; i
< n_row
; ++i
)
951 if (!isl_int_is_zero(row
[i
][col
]))
956 static int row_abs_min_non_zero(isl_int
**row
, unsigned n_row
, unsigned col
)
958 int i
, min
= row_first_non_zero(row
, n_row
, col
);
961 for (i
= min
+ 1; i
< n_row
; ++i
) {
962 if (isl_int_is_zero(row
[i
][col
]))
964 if (isl_int_abs_lt(row
[i
][col
], row
[min
][col
]))
970 static isl_stat
inv_exchange(__isl_keep isl_mat
**left
,
971 __isl_keep isl_mat
**right
, unsigned i
, unsigned j
)
973 *left
= isl_mat_swap_rows(*left
, i
, j
);
974 *right
= isl_mat_swap_rows(*right
, i
, j
);
976 if (!*left
|| !*right
)
977 return isl_stat_error
;
981 static void inv_oppose(
982 struct isl_mat
*left
, struct isl_mat
*right
, unsigned row
)
984 isl_seq_neg(left
->row
[row
]+row
, left
->row
[row
]+row
, left
->n_col
-row
);
985 isl_seq_neg(right
->row
[row
], right
->row
[row
], right
->n_col
);
988 static void inv_subtract(struct isl_mat
*left
, struct isl_mat
*right
,
989 unsigned row
, unsigned i
, isl_int m
)
992 isl_seq_combine(left
->row
[i
]+row
,
993 left
->ctx
->one
, left
->row
[i
]+row
,
994 m
, left
->row
[row
]+row
,
996 isl_seq_combine(right
->row
[i
], right
->ctx
->one
, right
->row
[i
],
997 m
, right
->row
[row
], right
->n_col
);
1000 /* Compute inv(left)*right
1002 __isl_give isl_mat
*isl_mat_inverse_product(__isl_take isl_mat
*left
,
1003 __isl_take isl_mat
*right
)
1008 if (!left
|| !right
)
1011 isl_assert(left
->ctx
, left
->n_row
== left
->n_col
, goto error
);
1012 isl_assert(left
->ctx
, left
->n_row
== right
->n_row
, goto error
);
1014 if (left
->n_row
== 0) {
1019 left
= isl_mat_cow(left
);
1020 right
= isl_mat_cow(right
);
1021 if (!left
|| !right
)
1026 for (row
= 0; row
< left
->n_row
; ++row
) {
1027 int pivot
, first
, i
, off
;
1028 pivot
= row_abs_min_non_zero(left
->row
+row
, left
->n_row
-row
, row
);
1032 isl_assert(left
->ctx
, pivot
>= 0, goto error
);
1036 if (inv_exchange(&left
, &right
, pivot
, row
) < 0)
1038 if (isl_int_is_neg(left
->row
[row
][row
]))
1039 inv_oppose(left
, right
, row
);
1041 while ((off
= row_first_non_zero(left
->row
+first
,
1042 left
->n_row
-first
, row
)) != -1) {
1044 isl_int_fdiv_q(a
, left
->row
[first
][row
],
1045 left
->row
[row
][row
]);
1046 inv_subtract(left
, right
, row
, first
, a
);
1047 if (!isl_int_is_zero(left
->row
[first
][row
])) {
1048 if (inv_exchange(&left
, &right
, row
, first
) < 0)
1054 for (i
= 0; i
< row
; ++i
) {
1055 if (isl_int_is_zero(left
->row
[i
][row
]))
1057 isl_int_gcd(a
, left
->row
[row
][row
], left
->row
[i
][row
]);
1058 isl_int_divexact(b
, left
->row
[i
][row
], a
);
1059 isl_int_divexact(a
, left
->row
[row
][row
], a
);
1061 isl_seq_combine(left
->row
[i
] + i
,
1062 a
, left
->row
[i
] + i
,
1063 b
, left
->row
[row
] + i
,
1065 isl_seq_combine(right
->row
[i
], a
, right
->row
[i
],
1066 b
, right
->row
[row
], right
->n_col
);
1071 isl_int_set(a
, left
->row
[0][0]);
1072 for (row
= 1; row
< left
->n_row
; ++row
)
1073 isl_int_lcm(a
, a
, left
->row
[row
][row
]);
1074 if (isl_int_is_zero(a
)){
1076 isl_assert(left
->ctx
, 0, goto error
);
1078 for (row
= 0; row
< left
->n_row
; ++row
) {
1079 isl_int_divexact(left
->row
[row
][row
], a
, left
->row
[row
][row
]);
1080 if (isl_int_is_one(left
->row
[row
][row
]))
1082 isl_seq_scale(right
->row
[row
], right
->row
[row
],
1083 left
->row
[row
][row
], right
->n_col
);
1091 isl_mat_free(right
);
1095 void isl_mat_col_scale(struct isl_mat
*mat
, unsigned col
, isl_int m
)
1099 for (i
= 0; i
< mat
->n_row
; ++i
)
1100 isl_int_mul(mat
->row
[i
][col
], mat
->row
[i
][col
], m
);
1103 void isl_mat_col_combine(struct isl_mat
*mat
, unsigned dst
,
1104 isl_int m1
, unsigned src1
, isl_int m2
, unsigned src2
)
1110 for (i
= 0; i
< mat
->n_row
; ++i
) {
1111 isl_int_mul(tmp
, m1
, mat
->row
[i
][src1
]);
1112 isl_int_addmul(tmp
, m2
, mat
->row
[i
][src2
]);
1113 isl_int_set(mat
->row
[i
][dst
], tmp
);
1118 __isl_give isl_mat
*isl_mat_right_inverse(__isl_take isl_mat
*mat
)
1120 struct isl_mat
*inv
;
1124 mat
= isl_mat_cow(mat
);
1128 inv
= isl_mat_identity(mat
->ctx
, mat
->n_col
);
1129 inv
= isl_mat_cow(inv
);
1135 for (row
= 0; row
< mat
->n_row
; ++row
) {
1136 int pivot
, first
, i
, off
;
1137 pivot
= isl_seq_abs_min_non_zero(mat
->row
[row
]+row
, mat
->n_col
-row
);
1141 isl_assert(mat
->ctx
, pivot
>= 0, goto error
);
1145 exchange(mat
, &inv
, NULL
, row
, pivot
, row
);
1146 if (isl_int_is_neg(mat
->row
[row
][row
]))
1147 oppose(mat
, &inv
, NULL
, row
, row
);
1149 while ((off
= isl_seq_first_non_zero(mat
->row
[row
]+first
,
1150 mat
->n_col
-first
)) != -1) {
1152 isl_int_fdiv_q(a
, mat
->row
[row
][first
],
1153 mat
->row
[row
][row
]);
1154 subtract(mat
, &inv
, NULL
, row
, row
, first
, a
);
1155 if (!isl_int_is_zero(mat
->row
[row
][first
]))
1156 exchange(mat
, &inv
, NULL
, row
, row
, first
);
1160 for (i
= 0; i
< row
; ++i
) {
1161 if (isl_int_is_zero(mat
->row
[row
][i
]))
1163 isl_int_gcd(a
, mat
->row
[row
][row
], mat
->row
[row
][i
]);
1164 isl_int_divexact(b
, mat
->row
[row
][i
], a
);
1165 isl_int_divexact(a
, mat
->row
[row
][row
], a
);
1167 isl_mat_col_combine(mat
, i
, a
, i
, b
, row
);
1168 isl_mat_col_combine(inv
, i
, a
, i
, b
, row
);
1173 isl_int_set(a
, mat
->row
[0][0]);
1174 for (row
= 1; row
< mat
->n_row
; ++row
)
1175 isl_int_lcm(a
, a
, mat
->row
[row
][row
]);
1176 if (isl_int_is_zero(a
)){
1180 for (row
= 0; row
< mat
->n_row
; ++row
) {
1181 isl_int_divexact(mat
->row
[row
][row
], a
, mat
->row
[row
][row
]);
1182 if (isl_int_is_one(mat
->row
[row
][row
]))
1184 isl_mat_col_scale(inv
, row
, mat
->row
[row
][row
]);
1197 __isl_give isl_mat
*isl_mat_transpose(__isl_take isl_mat
*mat
)
1199 struct isl_mat
*transpose
= NULL
;
1205 if (mat
->n_col
== mat
->n_row
) {
1206 mat
= isl_mat_cow(mat
);
1209 for (i
= 0; i
< mat
->n_row
; ++i
)
1210 for (j
= i
+ 1; j
< mat
->n_col
; ++j
)
1211 isl_int_swap(mat
->row
[i
][j
], mat
->row
[j
][i
]);
1214 transpose
= isl_mat_alloc(mat
->ctx
, mat
->n_col
, mat
->n_row
);
1217 for (i
= 0; i
< mat
->n_row
; ++i
)
1218 for (j
= 0; j
< mat
->n_col
; ++j
)
1219 isl_int_set(transpose
->row
[j
][i
], mat
->row
[i
][j
]);
1227 __isl_give isl_mat
*isl_mat_swap_cols(__isl_take isl_mat
*mat
,
1228 unsigned i
, unsigned j
)
1232 mat
= isl_mat_cow(mat
);
1233 if (check_col_range(mat
, i
, 1) < 0 ||
1234 check_col_range(mat
, j
, 1) < 0)
1235 return isl_mat_free(mat
);
1237 for (r
= 0; r
< mat
->n_row
; ++r
)
1238 isl_int_swap(mat
->row
[r
][i
], mat
->row
[r
][j
]);
1242 __isl_give isl_mat
*isl_mat_swap_rows(__isl_take isl_mat
*mat
,
1243 unsigned i
, unsigned j
)
1249 mat
= isl_mat_cow(mat
);
1250 if (check_row_range(mat
, i
, 1) < 0 ||
1251 check_row_range(mat
, j
, 1) < 0)
1252 return isl_mat_free(mat
);
1255 mat
->row
[i
] = mat
->row
[j
];
1260 /* Calculate the product of two matrices.
1262 * This function is optimized for operand matrices that contain many zeros and
1263 * skips multiplications where we know one of the operands is zero.
1265 __isl_give isl_mat
*isl_mat_product(__isl_take isl_mat
*left
,
1266 __isl_take isl_mat
*right
)
1269 struct isl_mat
*prod
;
1271 if (!left
|| !right
)
1273 isl_assert(left
->ctx
, left
->n_col
== right
->n_row
, goto error
);
1274 prod
= isl_mat_alloc(left
->ctx
, left
->n_row
, right
->n_col
);
1277 if (left
->n_col
== 0) {
1278 for (i
= 0; i
< prod
->n_row
; ++i
)
1279 isl_seq_clr(prod
->row
[i
], prod
->n_col
);
1281 isl_mat_free(right
);
1284 for (i
= 0; i
< prod
->n_row
; ++i
) {
1285 for (j
= 0; j
< prod
->n_col
; ++j
)
1286 isl_int_mul(prod
->row
[i
][j
],
1287 left
->row
[i
][0], right
->row
[0][j
]);
1288 for (k
= 1; k
< left
->n_col
; ++k
) {
1289 if (isl_int_is_zero(left
->row
[i
][k
]))
1291 for (j
= 0; j
< prod
->n_col
; ++j
)
1292 isl_int_addmul(prod
->row
[i
][j
],
1293 left
->row
[i
][k
], right
->row
[k
][j
]);
1297 isl_mat_free(right
);
1301 isl_mat_free(right
);
1305 /* Replace the variables x in the rows q by x' given by x = M x',
1306 * with M the matrix mat.
1308 * If the number of new variables is greater than the original
1309 * number of variables, then the rows q have already been
1310 * preextended. If the new number is smaller, then the coefficients
1311 * of the divs, which are not changed, need to be shifted down.
1312 * The row q may be the equalities, the inequalities or the
1313 * div expressions. In the latter case, has_div is true and
1314 * we need to take into account the extra denominator column.
1316 static int preimage(struct isl_ctx
*ctx
, isl_int
**q
, unsigned n
,
1317 unsigned n_div
, int has_div
, struct isl_mat
*mat
)
1323 if (mat
->n_col
>= mat
->n_row
)
1326 e
= mat
->n_row
- mat
->n_col
;
1328 for (i
= 0; i
< n
; ++i
)
1329 isl_int_mul(q
[i
][0], q
[i
][0], mat
->row
[0][0]);
1330 t
= isl_mat_sub_alloc6(mat
->ctx
, q
, 0, n
, has_div
, mat
->n_row
);
1331 t
= isl_mat_product(t
, mat
);
1334 for (i
= 0; i
< n
; ++i
) {
1335 isl_seq_swp_or_cpy(q
[i
] + has_div
, t
->row
[i
], t
->n_col
);
1336 isl_seq_cpy(q
[i
] + has_div
+ t
->n_col
,
1337 q
[i
] + has_div
+ t
->n_col
+ e
, n_div
);
1338 isl_seq_clr(q
[i
] + has_div
+ t
->n_col
+ n_div
, e
);
1344 /* Replace the variables x in bset by x' given by x = M x', with
1347 * If there are fewer variables x' then there are x, then we perform
1348 * the transformation in place, which means that, in principle,
1349 * this frees up some extra variables as the number
1350 * of columns remains constant, but we would have to extend
1351 * the div array too as the number of rows in this array is assumed
1352 * to be equal to extra.
1354 __isl_give isl_basic_set
*isl_basic_set_preimage(
1355 __isl_take isl_basic_set
*bset
, __isl_take isl_mat
*mat
)
1357 struct isl_ctx
*ctx
;
1363 bset
= isl_basic_set_cow(bset
);
1367 isl_assert(ctx
, bset
->dim
->nparam
== 0, goto error
);
1368 isl_assert(ctx
, 1+bset
->dim
->n_out
== mat
->n_row
, goto error
);
1369 isl_assert(ctx
, mat
->n_col
> 0, goto error
);
1371 if (mat
->n_col
> mat
->n_row
) {
1372 bset
= isl_basic_set_extend(bset
, 0, mat
->n_col
-1, 0, 0, 0);
1375 } else if (mat
->n_col
< mat
->n_row
) {
1376 bset
->dim
= isl_space_cow(bset
->dim
);
1379 bset
->dim
->n_out
-= mat
->n_row
- mat
->n_col
;
1382 if (preimage(ctx
, bset
->eq
, bset
->n_eq
, bset
->n_div
, 0,
1383 isl_mat_copy(mat
)) < 0)
1386 if (preimage(ctx
, bset
->ineq
, bset
->n_ineq
, bset
->n_div
, 0,
1387 isl_mat_copy(mat
)) < 0)
1390 if (preimage(ctx
, bset
->div
, bset
->n_div
, bset
->n_div
, 1, mat
) < 0)
1393 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
1394 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_REDUNDANT
);
1395 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED
);
1396 ISL_F_CLR(bset
, ISL_BASIC_SET_NORMALIZED_DIVS
);
1397 ISL_F_CLR(bset
, ISL_BASIC_SET_ALL_EQUALITIES
);
1399 bset
= isl_basic_set_simplify(bset
);
1400 bset
= isl_basic_set_finalize(bset
);
1406 isl_basic_set_free(bset
);
1410 __isl_give isl_set
*isl_set_preimage(
1411 __isl_take isl_set
*set
, __isl_take isl_mat
*mat
)
1415 set
= isl_set_cow(set
);
1419 for (i
= 0; i
< set
->n
; ++i
) {
1420 set
->p
[i
] = isl_basic_set_preimage(set
->p
[i
],
1425 if (mat
->n_col
!= mat
->n_row
) {
1426 set
->dim
= isl_space_cow(set
->dim
);
1429 set
->dim
->n_out
+= mat
->n_col
;
1430 set
->dim
->n_out
-= mat
->n_row
;
1433 ISL_F_CLR(set
, ISL_SET_NORMALIZED
);
1441 /* Replace the variables x starting at "first_col" in the rows "rows"
1442 * of some coefficient matrix by x' with x = M x' with M the matrix mat.
1443 * That is, replace the corresponding coefficients c by c M.
1445 isl_stat
isl_mat_sub_transform(isl_int
**row
, unsigned n_row
,
1446 unsigned first_col
, __isl_take isl_mat
*mat
)
1453 return isl_stat_error
;
1454 ctx
= isl_mat_get_ctx(mat
);
1455 t
= isl_mat_sub_alloc6(ctx
, row
, 0, n_row
, first_col
, mat
->n_row
);
1456 t
= isl_mat_product(t
, mat
);
1458 return isl_stat_error
;
1459 for (i
= 0; i
< n_row
; ++i
)
1460 isl_seq_swp_or_cpy(row
[i
] + first_col
, t
->row
[i
], t
->n_col
);
1465 void isl_mat_print_internal(__isl_keep isl_mat
*mat
, FILE *out
, int indent
)
1470 fprintf(out
, "%*snull mat\n", indent
, "");
1474 if (mat
->n_row
== 0)
1475 fprintf(out
, "%*s[]\n", indent
, "");
1477 for (i
= 0; i
< mat
->n_row
; ++i
) {
1479 fprintf(out
, "%*s[[", indent
, "");
1481 fprintf(out
, "%*s[", indent
+1, "");
1482 for (j
= 0; j
< mat
->n_col
; ++j
) {
1485 isl_int_print(out
, mat
->row
[i
][j
], 0);
1487 if (i
== mat
->n_row
-1)
1488 fprintf(out
, "]]\n");
1490 fprintf(out
, "]\n");
1494 void isl_mat_dump(__isl_keep isl_mat
*mat
)
1496 isl_mat_print_internal(mat
, stderr
, 0);
1499 __isl_give isl_mat
*isl_mat_drop_cols(__isl_take isl_mat
*mat
,
1500 unsigned col
, unsigned n
)
1507 mat
= isl_mat_cow(mat
);
1508 if (check_col_range(mat
, col
, n
) < 0)
1509 return isl_mat_free(mat
);
1511 if (col
!= mat
->n_col
-n
) {
1512 for (r
= 0; r
< mat
->n_row
; ++r
)
1513 isl_seq_cpy(mat
->row
[r
]+col
, mat
->row
[r
]+col
+n
,
1514 mat
->n_col
- col
- n
);
1520 __isl_give isl_mat
*isl_mat_drop_rows(__isl_take isl_mat
*mat
,
1521 unsigned row
, unsigned n
)
1525 mat
= isl_mat_cow(mat
);
1526 if (check_row_range(mat
, row
, n
) < 0)
1527 return isl_mat_free(mat
);
1529 for (r
= row
; r
+n
< mat
->n_row
; ++r
)
1530 mat
->row
[r
] = mat
->row
[r
+n
];
1536 __isl_give isl_mat
*isl_mat_insert_cols(__isl_take isl_mat
*mat
,
1537 unsigned col
, unsigned n
)
1541 if (check_col_range(mat
, col
, 0) < 0)
1542 return isl_mat_free(mat
);
1546 ext
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
+ n
);
1550 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
, 0, 0, col
);
1551 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, mat
->n_row
,
1552 col
+ n
, col
, mat
->n_col
- col
);
1561 __isl_give isl_mat
*isl_mat_insert_zero_cols(__isl_take isl_mat
*mat
,
1562 unsigned first
, unsigned n
)
1568 mat
= isl_mat_insert_cols(mat
, first
, n
);
1572 for (i
= 0; i
< mat
->n_row
; ++i
)
1573 isl_seq_clr(mat
->row
[i
] + first
, n
);
1578 __isl_give isl_mat
*isl_mat_add_zero_cols(__isl_take isl_mat
*mat
, unsigned n
)
1583 return isl_mat_insert_zero_cols(mat
, mat
->n_col
, n
);
1586 __isl_give isl_mat
*isl_mat_insert_rows(__isl_take isl_mat
*mat
,
1587 unsigned row
, unsigned n
)
1591 if (check_row_range(mat
, row
, 0) < 0)
1592 return isl_mat_free(mat
);
1596 ext
= isl_mat_alloc(mat
->ctx
, mat
->n_row
+ n
, mat
->n_col
);
1600 isl_mat_sub_copy(mat
->ctx
, ext
->row
, mat
->row
, row
, 0, 0, mat
->n_col
);
1601 isl_mat_sub_copy(mat
->ctx
, ext
->row
+ row
+ n
, mat
->row
+ row
,
1602 mat
->n_row
- row
, 0, 0, mat
->n_col
);
1611 __isl_give isl_mat
*isl_mat_add_rows(__isl_take isl_mat
*mat
, unsigned n
)
1616 return isl_mat_insert_rows(mat
, mat
->n_row
, n
);
1619 __isl_give isl_mat
*isl_mat_insert_zero_rows(__isl_take isl_mat
*mat
,
1620 unsigned row
, unsigned n
)
1624 mat
= isl_mat_insert_rows(mat
, row
, n
);
1628 for (i
= 0; i
< n
; ++i
)
1629 isl_seq_clr(mat
->row
[row
+ i
], mat
->n_col
);
1634 __isl_give isl_mat
*isl_mat_add_zero_rows(__isl_take isl_mat
*mat
, unsigned n
)
1639 return isl_mat_insert_zero_rows(mat
, mat
->n_row
, n
);
1642 void isl_mat_col_submul(struct isl_mat
*mat
,
1643 int dst_col
, isl_int f
, int src_col
)
1647 for (i
= 0; i
< mat
->n_row
; ++i
)
1648 isl_int_submul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1651 void isl_mat_col_add(__isl_keep isl_mat
*mat
, int dst_col
, int src_col
)
1658 for (i
= 0; i
< mat
->n_row
; ++i
)
1659 isl_int_add(mat
->row
[i
][dst_col
],
1660 mat
->row
[i
][dst_col
], mat
->row
[i
][src_col
]);
1663 void isl_mat_col_mul(struct isl_mat
*mat
, int dst_col
, isl_int f
, int src_col
)
1667 for (i
= 0; i
< mat
->n_row
; ++i
)
1668 isl_int_mul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1671 /* Add "f" times column "src_col" to column "dst_col" of "mat" and
1672 * return the result.
1674 __isl_give isl_mat
*isl_mat_col_addmul(__isl_take isl_mat
*mat
, int dst_col
,
1675 isl_int f
, int src_col
)
1679 if (check_col(mat
, dst_col
) < 0 || check_col(mat
, src_col
) < 0)
1680 return isl_mat_free(mat
);
1682 for (i
= 0; i
< mat
->n_row
; ++i
) {
1683 if (isl_int_is_zero(mat
->row
[i
][src_col
]))
1685 mat
= isl_mat_cow(mat
);
1688 isl_int_addmul(mat
->row
[i
][dst_col
], f
, mat
->row
[i
][src_col
]);
1694 /* Negate column "col" of "mat" and return the result.
1696 __isl_give isl_mat
*isl_mat_col_neg(__isl_take isl_mat
*mat
, int col
)
1700 if (check_col(mat
, col
) < 0)
1701 return isl_mat_free(mat
);
1703 for (i
= 0; i
< mat
->n_row
; ++i
) {
1704 if (isl_int_is_zero(mat
->row
[i
][col
]))
1706 mat
= isl_mat_cow(mat
);
1709 isl_int_neg(mat
->row
[i
][col
], mat
->row
[i
][col
]);
1715 /* Negate row "row" of "mat" and return the result.
1717 __isl_give isl_mat
*isl_mat_row_neg(__isl_take isl_mat
*mat
, int row
)
1719 if (check_row(mat
, row
) < 0)
1720 return isl_mat_free(mat
);
1721 if (isl_seq_first_non_zero(mat
->row
[row
], mat
->n_col
) == -1)
1723 mat
= isl_mat_cow(mat
);
1726 isl_seq_neg(mat
->row
[row
], mat
->row
[row
], mat
->n_col
);
1730 __isl_give isl_mat
*isl_mat_unimodular_complete(__isl_take isl_mat
*M
, int row
)
1733 struct isl_mat
*H
= NULL
, *Q
= NULL
;
1738 isl_assert(M
->ctx
, M
->n_row
== M
->n_col
, goto error
);
1740 H
= isl_mat_left_hermite(isl_mat_copy(M
), 0, NULL
, &Q
);
1741 M
->n_row
= M
->n_col
;
1744 for (r
= 0; r
< row
; ++r
)
1745 isl_assert(M
->ctx
, isl_int_is_one(H
->row
[r
][r
]), goto error
);
1746 for (r
= row
; r
< M
->n_row
; ++r
)
1747 isl_seq_cpy(M
->row
[r
], Q
->row
[r
], M
->n_col
);
1758 __isl_give isl_mat
*isl_mat_concat(__isl_take isl_mat
*top
,
1759 __isl_take isl_mat
*bot
)
1761 struct isl_mat
*mat
;
1766 isl_assert(top
->ctx
, top
->n_col
== bot
->n_col
, goto error
);
1767 if (top
->n_row
== 0) {
1771 if (bot
->n_row
== 0) {
1776 mat
= isl_mat_alloc(top
->ctx
, top
->n_row
+ bot
->n_row
, top
->n_col
);
1779 isl_mat_sub_copy(mat
->ctx
, mat
->row
, top
->row
, top
->n_row
,
1781 isl_mat_sub_copy(mat
->ctx
, mat
->row
+ top
->n_row
, bot
->row
, bot
->n_row
,
1792 isl_bool
isl_mat_is_equal(__isl_keep isl_mat
*mat1
, __isl_keep isl_mat
*mat2
)
1797 return isl_bool_error
;
1799 if (mat1
->n_row
!= mat2
->n_row
)
1800 return isl_bool_false
;
1802 if (mat1
->n_col
!= mat2
->n_col
)
1803 return isl_bool_false
;
1805 for (i
= 0; i
< mat1
->n_row
; ++i
)
1806 if (!isl_seq_eq(mat1
->row
[i
], mat2
->row
[i
], mat1
->n_col
))
1807 return isl_bool_false
;
1809 return isl_bool_true
;
1812 __isl_give isl_mat
*isl_mat_from_row_vec(__isl_take isl_vec
*vec
)
1814 struct isl_mat
*mat
;
1818 mat
= isl_mat_alloc(vec
->ctx
, 1, vec
->size
);
1822 isl_seq_cpy(mat
->row
[0], vec
->el
, vec
->size
);
1831 /* Return a copy of row "row" of "mat" as an isl_vec.
1833 __isl_give isl_vec
*isl_mat_get_row(__isl_keep isl_mat
*mat
, unsigned row
)
1839 if (row
>= mat
->n_row
)
1840 isl_die(mat
->ctx
, isl_error_invalid
, "row out of range",
1843 v
= isl_vec_alloc(isl_mat_get_ctx(mat
), mat
->n_col
);
1846 isl_seq_cpy(v
->el
, mat
->row
[row
], mat
->n_col
);
1851 __isl_give isl_mat
*isl_mat_vec_concat(__isl_take isl_mat
*top
,
1852 __isl_take isl_vec
*bot
)
1854 return isl_mat_concat(top
, isl_mat_from_row_vec(bot
));
1857 __isl_give isl_mat
*isl_mat_move_cols(__isl_take isl_mat
*mat
,
1858 unsigned dst_col
, unsigned src_col
, unsigned n
)
1864 if (n
== 0 || dst_col
== src_col
)
1867 res
= isl_mat_alloc(mat
->ctx
, mat
->n_row
, mat
->n_col
);
1871 if (dst_col
< src_col
) {
1872 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1874 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1875 dst_col
, src_col
, n
);
1876 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1877 dst_col
+ n
, dst_col
, src_col
- dst_col
);
1878 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1879 src_col
+ n
, src_col
+ n
,
1880 res
->n_col
- src_col
- n
);
1882 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1884 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1885 src_col
, src_col
+ n
, dst_col
- src_col
);
1886 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1887 dst_col
, src_col
, n
);
1888 isl_mat_sub_copy(res
->ctx
, res
->row
, mat
->row
, mat
->n_row
,
1889 dst_col
+ n
, dst_col
+ n
,
1890 res
->n_col
- dst_col
- n
);
1900 /* Return the gcd of the elements in row "row" of "mat" in *gcd.
1901 * Return isl_stat_ok on success and isl_stat_error on failure.
1903 isl_stat
isl_mat_row_gcd(__isl_keep isl_mat
*mat
, int row
, isl_int
*gcd
)
1905 if (check_row(mat
, row
) < 0)
1906 return isl_stat_error
;
1908 isl_seq_gcd(mat
->row
[row
], mat
->n_col
, gcd
);
1913 void isl_mat_gcd(__isl_keep isl_mat
*mat
, isl_int
*gcd
)
1918 isl_int_set_si(*gcd
, 0);
1923 for (i
= 0; i
< mat
->n_row
; ++i
) {
1924 isl_seq_gcd(mat
->row
[i
], mat
->n_col
, &g
);
1925 isl_int_gcd(*gcd
, *gcd
, g
);
1930 /* Return the result of scaling "mat" by a factor of "m".
1932 __isl_give isl_mat
*isl_mat_scale(__isl_take isl_mat
*mat
, isl_int m
)
1936 if (isl_int_is_one(m
))
1939 mat
= isl_mat_cow(mat
);
1943 for (i
= 0; i
< mat
->n_row
; ++i
)
1944 isl_seq_scale(mat
->row
[i
], mat
->row
[i
], m
, mat
->n_col
);
1949 __isl_give isl_mat
*isl_mat_scale_down(__isl_take isl_mat
*mat
, isl_int m
)
1953 if (isl_int_is_one(m
))
1956 mat
= isl_mat_cow(mat
);
1960 for (i
= 0; i
< mat
->n_row
; ++i
)
1961 isl_seq_scale_down(mat
->row
[i
], mat
->row
[i
], m
, mat
->n_col
);
1966 __isl_give isl_mat
*isl_mat_scale_down_row(__isl_take isl_mat
*mat
, int row
,
1969 if (isl_int_is_one(m
))
1972 mat
= isl_mat_cow(mat
);
1976 isl_seq_scale_down(mat
->row
[row
], mat
->row
[row
], m
, mat
->n_col
);
1981 __isl_give isl_mat
*isl_mat_normalize(__isl_take isl_mat
*mat
)
1989 isl_mat_gcd(mat
, &gcd
);
1990 mat
= isl_mat_scale_down(mat
, gcd
);
1996 __isl_give isl_mat
*isl_mat_normalize_row(__isl_take isl_mat
*mat
, int row
)
1998 mat
= isl_mat_cow(mat
);
2002 isl_seq_normalize(mat
->ctx
, mat
->row
[row
], mat
->n_col
);
2007 /* Number of initial non-zero columns.
2009 int isl_mat_initial_non_zero_cols(__isl_keep isl_mat
*mat
)
2016 for (i
= 0; i
< mat
->n_col
; ++i
)
2017 if (row_first_non_zero(mat
->row
, mat
->n_row
, i
) < 0)
2023 /* Return a basis for the space spanned by the rows of "mat".
2024 * Any basis will do, so simply perform Gaussian elimination and
2025 * remove the empty rows.
2027 __isl_give isl_mat
*isl_mat_row_basis(__isl_take isl_mat
*mat
)
2029 return isl_mat_reverse_gauss(mat
);
2032 /* Return rows that extend a basis of "mat1" to one
2033 * that covers both "mat1" and "mat2".
2034 * The Hermite normal form of the concatenation of the two matrices is
2037 * [ M1 ] = [ H1 0 0 ] [ Q2 ]
2038 * [ M2 ] = [ H2 H3 0 ] [ Q3 ]
2040 * The number of columns in H1 and H3 determine the number of rows
2041 * in Q1 and Q2. Q1 is a basis for M1, while Q2 extends this basis
2044 __isl_give isl_mat
*isl_mat_row_basis_extension(
2045 __isl_take isl_mat
*mat1
, __isl_take isl_mat
*mat2
)
2051 n1
= isl_mat_rows(mat1
);
2052 H
= isl_mat_concat(mat1
, mat2
);
2053 H
= isl_mat_left_hermite(H
, 0, NULL
, &Q
);
2057 r1
= hermite_first_zero_col(H
, 0, n1
);
2058 r
= hermite_first_zero_col(H
, r1
, H
->n_row
);
2059 n_row
= isl_mat_rows(Q
);
2060 Q
= isl_mat_drop_rows(Q
, r
, n_row
- r
);
2061 Q
= isl_mat_drop_rows(Q
, 0, r1
);