3 #include "isl_map_private.h"
5 struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
6 unsigned n_row, unsigned n_col)
11 mat = isl_alloc_type(ctx, struct isl_mat);
16 mat->block = isl_blk_alloc(ctx, n_row * n_col);
17 if (isl_blk_is_error(mat->block))
19 mat->row = isl_alloc_array(ctx, isl_int *, n_row);
23 for (i = 0; i < n_row; ++i)
24 mat->row[i] = mat->block.data + i * n_col;
33 isl_blk_free(ctx, mat->block);
38 struct isl_mat *isl_mat_sub_alloc(struct isl_ctx *ctx, isl_int **row,
39 unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
44 mat = isl_alloc_type(ctx, struct isl_mat);
47 mat->row = isl_alloc_array(ctx, isl_int *, n_row);
50 for (i = 0; i < n_row; ++i)
51 mat->row[i] = row[first_row+i] + first_col;
55 mat->block = isl_blk_empty();
56 mat->flags = ISL_MAT_BORROWED;
63 void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
64 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
68 for (i = 0; i < n_row; ++i)
69 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
72 void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
73 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
77 for (i = 0; i < n_row; ++i)
78 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
81 struct isl_mat *isl_mat_copy(struct isl_ctx *ctx, struct isl_mat *mat)
90 struct isl_mat *isl_mat_dup(struct isl_ctx *ctx, struct isl_mat *mat)
97 mat2 = isl_mat_alloc(ctx, mat->n_row, mat->n_col);
100 for (i = 0; i < mat->n_row; ++i)
101 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
105 struct isl_mat *isl_mat_cow(struct isl_ctx *ctx, struct isl_mat *mat)
107 struct isl_mat *mat2;
111 if (mat->ref == 1 && !F_ISSET(mat, ISL_MAT_BORROWED))
114 mat2 = isl_mat_dup(ctx, mat);
115 isl_mat_free(ctx, mat);
119 void isl_mat_free(struct isl_ctx *ctx, struct isl_mat *mat)
127 if (!F_ISSET(mat, ISL_MAT_BORROWED))
128 isl_blk_free(ctx, mat->block);
133 struct isl_mat *isl_mat_identity(struct isl_ctx *ctx, unsigned n_row)
138 mat = isl_mat_alloc(ctx, n_row, n_row);
141 for (i = 0; i < n_row; ++i) {
142 isl_seq_clr(mat->row[i], i);
143 isl_int_set_si(mat->row[i][i], 1);
144 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
150 struct isl_vec *isl_mat_vec_product(struct isl_ctx *ctx,
151 struct isl_mat *mat, struct isl_vec *vec)
154 struct isl_vec *prod;
159 isl_assert(ctx, mat->n_col == vec->size, goto error);
161 prod = isl_vec_alloc(ctx, mat->n_row);
165 for (i = 0; i < prod->size; ++i)
166 isl_seq_inner_product(mat->row[i], vec->block.data, vec->size,
167 &prod->block.data[i]);
168 isl_mat_free(ctx, mat);
169 isl_vec_free(ctx, vec);
172 isl_mat_free(ctx, mat);
173 isl_vec_free(ctx, vec);
177 struct isl_mat *isl_mat_aff_direct_sum(struct isl_ctx *ctx,
178 struct isl_mat *left, struct isl_mat *right)
186 isl_assert(ctx, left->n_row == right->n_row, goto error);
187 isl_assert(ctx, left->n_row >= 1, goto error);
188 isl_assert(ctx, left->n_col >= 1, goto error);
189 isl_assert(ctx, right->n_col >= 1, goto error);
191 isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
194 isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
197 sum = isl_mat_alloc(ctx, left->n_row, left->n_col + right->n_col - 1);
200 isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
201 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
202 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
204 isl_seq_clr(sum->row[0]+1, sum->n_col-1);
205 for (i = 1; i < sum->n_row; ++i) {
206 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
207 isl_int_addmul(sum->row[i][0],
208 right->row[0][0], right->row[i][0]);
209 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
211 isl_seq_scale(sum->row[i]+left->n_col,
212 right->row[i]+1, right->row[0][0],
216 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
217 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
218 isl_mat_free(ctx, left);
219 isl_mat_free(ctx, right);
222 isl_mat_free(ctx, left);
223 isl_mat_free(ctx, right);
227 static void exchange(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
228 struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
231 for (r = row; r < M->n_row; ++r)
232 isl_int_swap(M->row[r][i], M->row[r][j]);
234 for (r = 0; r < (*U)->n_row; ++r)
235 isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
238 isl_mat_swap_rows(ctx, *Q, i, j);
241 static void subtract(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
242 struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
245 for (r = row; r < M->n_row; ++r)
246 isl_int_submul(M->row[r][j], m, M->row[r][i]);
248 for (r = 0; r < (*U)->n_row; ++r)
249 isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
252 for (r = 0; r < (*Q)->n_col; ++r)
253 isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
257 static void oppose(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
258 struct isl_mat **Q, unsigned row, unsigned col)
261 for (r = row; r < M->n_row; ++r)
262 isl_int_neg(M->row[r][col], M->row[r][col]);
264 for (r = 0; r < (*U)->n_row; ++r)
265 isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
268 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
271 /* Given matrix M, compute
276 * with U and Q unimodular matrices and H a matrix in column echelon form
277 * such that on each echelon row the entries in the non-echelon column
278 * are non-negative and stricly smaller than the entries in the echelon
279 * column. If U or Q are NULL, then these matrices are not computed.
281 struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
282 struct isl_mat *M, struct isl_mat **U, struct isl_mat **Q)
293 M = isl_mat_cow(ctx, M);
297 *U = isl_mat_identity(ctx, M->n_col);
302 *Q = isl_mat_identity(ctx, M->n_col);
309 for (row = 0; row < M->n_row; ++row) {
311 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
316 exchange(ctx, M, U, Q, row, first, col);
317 if (isl_int_is_neg(M->row[row][col]))
318 oppose(ctx, M, U, Q, row, col);
320 while ((off = isl_seq_first_non_zero(M->row[row]+first,
321 M->n_col-first)) != -1) {
323 isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
324 subtract(ctx, M, U, Q, row, col, first, c);
325 if (!isl_int_is_zero(M->row[row][first]))
326 exchange(ctx, M, U, Q, row, first, col);
330 for (i = 0; i < col; ++i) {
331 if (isl_int_is_zero(M->row[row][i]))
333 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
334 if (isl_int_is_zero(c))
336 subtract(ctx, M, U, Q, row, col, i, c);
345 isl_mat_free(ctx, *Q);
349 isl_mat_free(ctx, *U);
355 struct isl_mat *isl_mat_lin_to_aff(struct isl_ctx *ctx, struct isl_mat *mat)
358 struct isl_mat *mat2;
362 mat2 = isl_mat_alloc(ctx, 1+mat->n_row, 1+mat->n_col);
365 isl_int_set_si(mat2->row[0][0], 1);
366 isl_seq_clr(mat2->row[0]+1, mat->n_col);
367 for (i = 0; i < mat->n_row; ++i) {
368 isl_int_set_si(mat2->row[1+i][0], 0);
369 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
371 isl_mat_free(ctx, mat);
375 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
379 for (i = 0; i < n_row; ++i)
380 if (!isl_int_is_zero(row[i][col]))
385 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
387 int i, min = row_first_non_zero(row, n_row, col);
390 for (i = min + 1; i < n_row; ++i) {
391 if (isl_int_is_zero(row[i][col]))
393 if (isl_int_abs_lt(row[i][col], row[min][col]))
399 static void inv_exchange(struct isl_ctx *ctx,
400 struct isl_mat *left, struct isl_mat *right,
401 unsigned i, unsigned j)
403 left = isl_mat_swap_rows(ctx, left, i, j);
404 right = isl_mat_swap_rows(ctx, right, i, j);
407 static void inv_oppose(struct isl_ctx *ctx,
408 struct isl_mat *left, struct isl_mat *right, unsigned row)
410 isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
411 isl_seq_neg(right->row[row], right->row[row], right->n_col);
414 static void inv_subtract(struct isl_ctx *ctx,
415 struct isl_mat *left, struct isl_mat *right,
416 unsigned row, unsigned i, isl_int m)
419 isl_seq_combine(left->row[i]+row,
420 ctx->one, left->row[i]+row,
421 m, left->row[row]+row,
423 isl_seq_combine(right->row[i], ctx->one, right->row[i],
424 m, right->row[row], right->n_col);
427 /* Compute inv(left)*right
429 struct isl_mat *isl_mat_inverse_product(struct isl_ctx *ctx,
430 struct isl_mat *left, struct isl_mat *right)
438 isl_assert(ctx, left->n_row == left->n_col, goto error);
439 isl_assert(ctx, left->n_row == right->n_row, goto error);
441 if (left->n_row == 0) {
442 isl_mat_free(ctx, left);
446 left = isl_mat_cow(ctx, left);
447 right = isl_mat_cow(ctx, right);
453 for (row = 0; row < left->n_row; ++row) {
454 int pivot, first, i, off;
455 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
459 isl_assert(ctx, pivot >= 0, goto error);
463 inv_exchange(ctx, left, right, pivot, row);
464 if (isl_int_is_neg(left->row[row][row]))
465 inv_oppose(ctx, left, right, row);
467 while ((off = row_first_non_zero(left->row+first,
468 left->n_row-first, row)) != -1) {
470 isl_int_fdiv_q(a, left->row[first][row],
471 left->row[row][row]);
472 inv_subtract(ctx, left, right, row, first, a);
473 if (!isl_int_is_zero(left->row[first][row]))
474 inv_exchange(ctx, left, right, row, first);
478 for (i = 0; i < row; ++i) {
479 if (isl_int_is_zero(left->row[i][row]))
481 isl_int_gcd(a, left->row[row][row], left->row[i][row]);
482 isl_int_divexact(b, left->row[i][row], a);
483 isl_int_divexact(a, left->row[row][row], a);
485 isl_seq_combine(left->row[i]+row,
487 b, left->row[row]+row,
489 isl_seq_combine(right->row[i], a, right->row[i],
490 b, right->row[row], right->n_col);
495 isl_int_set(a, left->row[0][0]);
496 for (row = 1; row < left->n_row; ++row)
497 isl_int_lcm(a, a, left->row[row][row]);
498 if (isl_int_is_zero(a)){
500 isl_assert(ctx, 0, goto error);
502 for (row = 0; row < left->n_row; ++row) {
503 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
504 if (isl_int_is_one(left->row[row][row]))
506 isl_seq_scale(right->row[row], right->row[row],
507 left->row[row][row], right->n_col);
511 isl_mat_free(ctx, left);
514 isl_mat_free(ctx, left);
515 isl_mat_free(ctx, right);
519 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
520 struct isl_mat *mat, unsigned i, unsigned j)
526 mat = isl_mat_cow(ctx, mat);
530 mat->row[i] = mat->row[j];
535 struct isl_mat *isl_mat_product(struct isl_ctx *ctx,
536 struct isl_mat *left, struct isl_mat *right)
539 struct isl_mat *prod;
543 isl_assert(ctx, left->n_col == right->n_row, goto error);
544 prod = isl_mat_alloc(ctx, left->n_row, right->n_col);
547 if (left->n_col == 0) {
548 for (i = 0; i < prod->n_row; ++i)
549 isl_seq_clr(prod->row[i], prod->n_col);
552 for (i = 0; i < prod->n_row; ++i) {
553 for (j = 0; j < prod->n_col; ++j) {
554 isl_int_mul(prod->row[i][j],
555 left->row[i][0], right->row[0][j]);
556 for (k = 1; k < left->n_col; ++k)
557 isl_int_addmul(prod->row[i][j],
558 left->row[i][k], right->row[k][j]);
561 isl_mat_free(ctx, left);
562 isl_mat_free(ctx, right);
565 isl_mat_free(ctx, left);
566 isl_mat_free(ctx, right);
570 struct isl_basic_set *isl_basic_set_preimage(struct isl_ctx *ctx,
571 struct isl_basic_set *bset, struct isl_mat *mat)
579 bset = isl_basic_set_cow(ctx, bset);
583 isl_assert(ctx, bset->nparam == 0, goto error);
584 isl_assert(ctx, bset->n_div == 0, goto error);
585 isl_assert(ctx, 1+bset->dim == mat->n_row, goto error);
587 t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
588 t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
591 for (i = 0; i < bset->n_eq; ++i)
592 isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
593 isl_mat_free(ctx, t);
595 t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
596 t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
599 for (i = 0; i < bset->n_ineq; ++i)
600 isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
601 isl_mat_free(ctx, t);
603 bset->dim -= mat->n_row - mat->n_col;
604 bset = isl_basic_set_simplify(ctx, bset);
605 bset = isl_basic_set_finalize(ctx, bset);
607 isl_mat_free(ctx, mat);
610 isl_mat_free(ctx, mat);
611 isl_basic_set_free(ctx, bset);
615 void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
616 FILE *out, int indent)
621 fprintf(out, "%*s[]\n", indent, "");
623 for (i = 0; i < mat->n_row; ++i) {
625 fprintf(out, "%*s[[", indent, "");
627 fprintf(out, "%*s[", indent+1, "");
628 for (j = 0; j < mat->n_col; ++j) {
631 isl_int_print(out, mat->row[i][j]);
633 if (i == mat->n_row-1)
634 fprintf(out, "]]\n");
640 struct isl_mat *isl_mat_drop_col(struct isl_ctx *ctx, struct isl_mat *mat,
648 if (col != mat->n_col-1) {
649 for (r = 0; r < mat->n_row; ++r)
650 isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+1,
651 mat->n_col - col - 1);
657 struct isl_mat *isl_mat_drop_rows(struct isl_ctx *ctx, struct isl_mat *mat,
658 unsigned row, unsigned n)
665 for (r = row; r+n < mat->n_row; ++r)
666 mat->row[r] = mat->row[r+n];