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);
17 mat->block = isl_blk_alloc(ctx, n_row * n_col);
20 mat->row = isl_alloc_array(ctx, isl_int *, n_row);
24 for (i = 0; i < n_row; ++i)
25 mat->row[i] = mat->block.data + i * n_col;
34 isl_blk_free(ctx, mat->block);
39 struct isl_mat *isl_mat_sub_alloc(struct isl_ctx *ctx, isl_int **row,
40 unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
45 mat = isl_alloc_type(ctx, struct isl_mat);
48 mat->row = isl_alloc_array(ctx, isl_int *, n_row);
51 for (i = 0; i < n_row; ++i)
52 mat->row[i] = row[first_row+i] + first_col;
56 mat->block = isl_blk_empty();
57 mat->flags = ISL_MAT_BORROWED;
64 void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
65 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
69 for (i = 0; i < n_row; ++i)
70 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
73 void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
74 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
78 for (i = 0; i < n_row; ++i)
79 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
82 struct isl_mat *isl_mat_copy(struct isl_ctx *ctx, struct isl_mat *mat)
91 struct isl_mat *isl_mat_dup(struct isl_ctx *ctx, struct isl_mat *mat)
98 mat2 = isl_mat_alloc(ctx, mat->n_row, mat->n_col);
101 for (i = 0; i < mat->n_row; ++i)
102 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
106 struct isl_mat *isl_mat_cow(struct isl_ctx *ctx, struct isl_mat *mat)
108 struct isl_mat *mat2;
112 if (mat->ref == 1 && !F_ISSET(mat, ISL_MAT_BORROWED))
115 mat2 = isl_mat_dup(ctx, mat);
116 isl_mat_free(ctx, mat);
120 void isl_mat_free(struct isl_ctx *ctx, struct isl_mat *mat)
128 if (!F_ISSET(mat, ISL_MAT_BORROWED))
129 isl_blk_free(ctx, mat->block);
134 struct isl_mat *isl_mat_identity(struct isl_ctx *ctx, unsigned n_row)
139 mat = isl_mat_alloc(ctx, n_row, n_row);
142 for (i = 0; i < n_row; ++i) {
143 isl_seq_clr(mat->row[i], i);
144 isl_int_set_si(mat->row[i][i], 1);
145 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
151 struct isl_vec *isl_mat_vec_product(struct isl_ctx *ctx,
152 struct isl_mat *mat, struct isl_vec *vec)
155 struct isl_vec *prod;
160 isl_assert(ctx, mat->n_col == vec->size, goto error);
162 prod = isl_vec_alloc(ctx, mat->n_row);
166 for (i = 0; i < prod->size; ++i)
167 isl_seq_inner_product(mat->row[i], vec->block.data, vec->size,
168 &prod->block.data[i]);
169 isl_mat_free(ctx, mat);
170 isl_vec_free(ctx, vec);
173 isl_mat_free(ctx, mat);
174 isl_vec_free(ctx, vec);
178 struct isl_mat *isl_mat_aff_direct_sum(struct isl_ctx *ctx,
179 struct isl_mat *left, struct isl_mat *right)
187 isl_assert(ctx, left->n_row == right->n_row, goto error);
188 isl_assert(ctx, left->n_row >= 1, goto error);
189 isl_assert(ctx, left->n_col >= 1, goto error);
190 isl_assert(ctx, right->n_col >= 1, goto error);
192 isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
195 isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
198 sum = isl_mat_alloc(ctx, left->n_row, left->n_col + right->n_col - 1);
201 isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
202 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
203 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
205 isl_seq_clr(sum->row[0]+1, sum->n_col-1);
206 for (i = 1; i < sum->n_row; ++i) {
207 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
208 isl_int_addmul(sum->row[i][0],
209 right->row[0][0], right->row[i][0]);
210 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
212 isl_seq_scale(sum->row[i]+left->n_col,
213 right->row[i]+1, right->row[0][0],
217 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
218 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
219 isl_mat_free(ctx, left);
220 isl_mat_free(ctx, right);
223 isl_mat_free(ctx, left);
224 isl_mat_free(ctx, right);
228 static void exchange(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
229 struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
232 for (r = row; r < M->n_row; ++r)
233 isl_int_swap(M->row[r][i], M->row[r][j]);
235 for (r = 0; r < (*U)->n_row; ++r)
236 isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
239 isl_mat_swap_rows(ctx, *Q, i, j);
242 static void subtract(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
243 struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
246 for (r = row; r < M->n_row; ++r)
247 isl_int_submul(M->row[r][j], m, M->row[r][i]);
249 for (r = 0; r < (*U)->n_row; ++r)
250 isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
253 for (r = 0; r < (*Q)->n_col; ++r)
254 isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
258 static void oppose(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
259 struct isl_mat **Q, unsigned row, unsigned col)
262 for (r = row; r < M->n_row; ++r)
263 isl_int_neg(M->row[r][col], M->row[r][col]);
265 for (r = 0; r < (*U)->n_row; ++r)
266 isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
269 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
272 /* Given matrix M, compute
277 * with U and Q unimodular matrices and H a matrix in column echelon form
278 * such that on each echelon row the entries in the non-echelon column
279 * are non-negative and stricly smaller than the entries in the echelon
280 * column. If U or Q are NULL, then these matrices are not computed.
282 struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
283 struct isl_mat *M, struct isl_mat **U, struct isl_mat **Q)
294 M = isl_mat_cow(ctx, M);
298 *U = isl_mat_identity(ctx, M->n_col);
303 *Q = isl_mat_identity(ctx, M->n_col);
310 for (row = 0; row < M->n_row; ++row) {
312 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
317 exchange(ctx, M, U, Q, row, first, col);
318 if (isl_int_is_neg(M->row[row][col]))
319 oppose(ctx, M, U, Q, row, col);
321 while ((off = isl_seq_first_non_zero(M->row[row]+first,
322 M->n_col-first)) != -1) {
324 isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
325 subtract(ctx, M, U, Q, row, col, first, c);
326 if (!isl_int_is_zero(M->row[row][first]))
327 exchange(ctx, M, U, Q, row, first, col);
331 for (i = 0; i < col; ++i) {
332 if (isl_int_is_zero(M->row[row][i]))
334 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
335 if (isl_int_is_zero(c))
337 subtract(ctx, M, U, Q, row, col, i, c);
346 isl_mat_free(ctx, *Q);
350 isl_mat_free(ctx, *U);
356 struct isl_mat *isl_mat_lin_to_aff(struct isl_ctx *ctx, struct isl_mat *mat)
359 struct isl_mat *mat2;
363 mat2 = isl_mat_alloc(ctx, 1+mat->n_row, 1+mat->n_col);
366 isl_int_set_si(mat2->row[0][0], 1);
367 isl_seq_clr(mat2->row[0]+1, mat->n_col);
368 for (i = 0; i < mat->n_row; ++i) {
369 isl_int_set_si(mat2->row[1+i][0], 0);
370 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
372 isl_mat_free(ctx, mat);
376 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
380 for (i = 0; i < n_row; ++i)
381 if (!isl_int_is_zero(row[i][col]))
386 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
388 int i, min = row_first_non_zero(row, n_row, col);
391 for (i = min + 1; i < n_row; ++i) {
392 if (isl_int_is_zero(row[i][col]))
394 if (isl_int_abs_lt(row[i][col], row[min][col]))
400 static void inv_exchange(struct isl_ctx *ctx,
401 struct isl_mat *left, struct isl_mat *right,
402 unsigned i, unsigned j)
404 left = isl_mat_swap_rows(ctx, left, i, j);
405 right = isl_mat_swap_rows(ctx, right, i, j);
408 static void inv_oppose(struct isl_ctx *ctx,
409 struct isl_mat *left, struct isl_mat *right, unsigned row)
411 isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
412 isl_seq_neg(right->row[row], right->row[row], right->n_col);
415 static void inv_subtract(struct isl_ctx *ctx,
416 struct isl_mat *left, struct isl_mat *right,
417 unsigned row, unsigned i, isl_int m)
420 isl_seq_combine(left->row[i]+row,
421 ctx->one, left->row[i]+row,
422 m, left->row[row]+row,
424 isl_seq_combine(right->row[i], ctx->one, right->row[i],
425 m, right->row[row], right->n_col);
428 /* Compute inv(left)*right
430 struct isl_mat *isl_mat_inverse_product(struct isl_ctx *ctx,
431 struct isl_mat *left, struct isl_mat *right)
439 isl_assert(ctx, left->n_row == left->n_col, goto error);
440 isl_assert(ctx, left->n_row == right->n_row, goto error);
442 if (left->n_row == 0) {
443 isl_mat_free(ctx, left);
447 left = isl_mat_cow(ctx, left);
448 right = isl_mat_cow(ctx, right);
454 for (row = 0; row < left->n_row; ++row) {
455 int pivot, first, i, off;
456 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
460 isl_assert(ctx, pivot >= 0, goto error);
464 inv_exchange(ctx, left, right, pivot, row);
465 if (isl_int_is_neg(left->row[row][row]))
466 inv_oppose(ctx, left, right, row);
468 while ((off = row_first_non_zero(left->row+first,
469 left->n_row-first, row)) != -1) {
471 isl_int_fdiv_q(a, left->row[first][row],
472 left->row[row][row]);
473 inv_subtract(ctx, left, right, row, first, a);
474 if (!isl_int_is_zero(left->row[first][row]))
475 inv_exchange(ctx, left, right, row, first);
479 for (i = 0; i < row; ++i) {
480 if (isl_int_is_zero(left->row[i][row]))
482 isl_int_gcd(a, left->row[row][row], left->row[i][row]);
483 isl_int_divexact(b, left->row[i][row], a);
484 isl_int_divexact(a, left->row[row][row], a);
486 isl_seq_combine(left->row[i]+row,
488 b, left->row[row]+row,
490 isl_seq_combine(right->row[i], a, right->row[i],
491 b, right->row[row], right->n_col);
496 isl_int_set(a, left->row[0][0]);
497 for (row = 1; row < left->n_row; ++row)
498 isl_int_lcm(a, a, left->row[row][row]);
499 if (isl_int_is_zero(a)){
501 isl_assert(ctx, 0, goto error);
503 for (row = 0; row < left->n_row; ++row) {
504 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
505 if (isl_int_is_one(left->row[row][row]))
507 isl_seq_scale(right->row[row], right->row[row],
508 left->row[row][row], right->n_col);
512 isl_mat_free(ctx, left);
515 isl_mat_free(ctx, left);
516 isl_mat_free(ctx, right);
520 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
521 struct isl_mat *mat, unsigned i, unsigned j)
527 mat = isl_mat_cow(ctx, mat);
531 mat->row[i] = mat->row[j];
536 struct isl_mat *isl_mat_product(struct isl_ctx *ctx,
537 struct isl_mat *left, struct isl_mat *right)
540 struct isl_mat *prod;
544 isl_assert(ctx, left->n_col == right->n_row, goto error);
545 prod = isl_mat_alloc(ctx, left->n_row, right->n_col);
548 if (left->n_col == 0) {
549 for (i = 0; i < prod->n_row; ++i)
550 isl_seq_clr(prod->row[i], prod->n_col);
553 for (i = 0; i < prod->n_row; ++i) {
554 for (j = 0; j < prod->n_col; ++j) {
555 isl_int_mul(prod->row[i][j],
556 left->row[i][0], right->row[0][j]);
557 for (k = 1; k < left->n_col; ++k)
558 isl_int_addmul(prod->row[i][j],
559 left->row[i][k], right->row[k][j]);
562 isl_mat_free(ctx, left);
563 isl_mat_free(ctx, right);
566 isl_mat_free(ctx, left);
567 isl_mat_free(ctx, right);
571 struct isl_basic_set *isl_basic_set_preimage(struct isl_ctx *ctx,
572 struct isl_basic_set *bset, struct isl_mat *mat)
580 bset = isl_basic_set_cow(ctx, bset);
584 isl_assert(ctx, bset->nparam == 0, goto error);
585 isl_assert(ctx, bset->n_div == 0, goto error);
586 isl_assert(ctx, 1+bset->dim == mat->n_row, goto error);
588 t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
589 t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
592 for (i = 0; i < bset->n_eq; ++i)
593 isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
594 isl_mat_free(ctx, t);
596 t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
597 t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
600 for (i = 0; i < bset->n_ineq; ++i)
601 isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
602 isl_mat_free(ctx, t);
604 bset->dim -= mat->n_row - mat->n_col;
605 bset = isl_basic_set_simplify(ctx, bset);
606 bset = isl_basic_set_finalize(ctx, bset);
608 isl_mat_free(ctx, mat);
611 isl_mat_free(ctx, mat);
612 isl_basic_set_free(ctx, bset);
616 void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
617 FILE *out, int indent)
622 fprintf(out, "%*s[]\n", indent, "");
624 for (i = 0; i < mat->n_row; ++i) {
626 fprintf(out, "%*s[[", indent, "");
628 fprintf(out, "%*s[", indent+1, "");
629 for (j = 0; j < mat->n_col; ++j) {
632 isl_int_print(out, mat->row[i][j]);
634 if (i == mat->n_row-1)
635 fprintf(out, "]]\n");
641 struct isl_mat *isl_mat_drop_col(struct isl_ctx *ctx, struct isl_mat *mat,
649 if (col != mat->n_col-1) {
650 for (r = 0; r < mat->n_row; ++r)
651 isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+1,
652 mat->n_col - col - 1);
658 struct isl_mat *isl_mat_drop_rows(struct isl_ctx *ctx, struct isl_mat *mat,
659 unsigned row, unsigned n)
666 for (r = row; r+n < mat->n_row; ++r)
667 mat->row[r] = mat->row[r+n];