2111 lines · c
1/*2 * Copyright 2008-2009 Katholieke Universiteit Leuven3 * Copyright 2010 INRIA Saclay4 * Copyright 2014 Ecole Normale Superieure5 * Copyright 2017 Sven Verdoolaege6 *7 * Use of this software is governed by the MIT license8 *9 * Written by Sven Verdoolaege, K.U.Leuven, Departement10 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium11 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,12 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France13 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France14 */15 16#include <isl_ctx_private.h>17#include <isl_map_private.h>18#include <isl/space.h>19#include <isl_seq.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 25isl_ctx *isl_mat_get_ctx(__isl_keep isl_mat *mat)26{27 return mat ? mat->ctx : NULL;28}29 30/* Return a hash value that digests "mat".31 */32uint32_t isl_mat_get_hash(__isl_keep isl_mat *mat)33{34 int i;35 uint32_t hash;36 37 if (!mat)38 return 0;39 40 hash = isl_hash_init();41 isl_hash_byte(hash, mat->n_row & 0xFF);42 isl_hash_byte(hash, mat->n_col & 0xFF);43 for (i = 0; i < mat->n_row; ++i) {44 uint32_t row_hash;45 46 row_hash = isl_seq_get_hash(mat->row[i], mat->n_col);47 isl_hash_hash(hash, row_hash);48 }49 50 return hash;51}52 53__isl_give isl_mat *isl_mat_alloc(isl_ctx *ctx,54 unsigned n_row, unsigned n_col)55{56 int i;57 struct isl_mat *mat;58 59 mat = isl_alloc_type(ctx, struct isl_mat);60 if (!mat)61 return NULL;62 63 mat->row = NULL;64 mat->block = isl_blk_alloc(ctx, n_row * n_col);65 if (isl_blk_is_error(mat->block))66 goto error;67 mat->row = isl_calloc_array(ctx, isl_int *, n_row);68 if (n_row && !mat->row)69 goto error;70 71 if (n_col != 0) {72 for (i = 0; i < n_row; ++i)73 mat->row[i] = mat->block.data + i * n_col;74 }75 76 mat->ctx = ctx;77 isl_ctx_ref(ctx);78 mat->ref = 1;79 mat->n_row = n_row;80 mat->n_col = n_col;81 mat->max_col = n_col;82 mat->flags = 0;83 84 return mat;85error:86 isl_blk_free(ctx, mat->block);87 free(mat);88 return NULL;89}90 91__isl_give isl_mat *isl_mat_extend(__isl_take isl_mat *mat,92 unsigned n_row, unsigned n_col)93{94 int i;95 isl_int *old;96 isl_int **row;97 98 if (!mat)99 return NULL;100 101 if (mat->max_col >= n_col && mat->n_row >= n_row) {102 if (mat->n_col < n_col)103 mat->n_col = n_col;104 return mat;105 }106 107 if (mat->max_col < n_col) {108 struct isl_mat *new_mat;109 110 if (n_row < mat->n_row)111 n_row = mat->n_row;112 new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);113 if (!new_mat)114 goto error;115 for (i = 0; i < mat->n_row; ++i)116 isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);117 isl_mat_free(mat);118 return new_mat;119 }120 121 mat = isl_mat_cow(mat);122 if (!mat)123 goto error;124 125 old = mat->block.data;126 mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);127 if (isl_blk_is_error(mat->block))128 goto error;129 row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);130 if (n_row && !row)131 goto error;132 mat->row = row;133 134 for (i = 0; i < mat->n_row; ++i)135 mat->row[i] = mat->block.data + (mat->row[i] - old);136 for (i = mat->n_row; i < n_row; ++i)137 mat->row[i] = mat->block.data + i * mat->max_col;138 mat->n_row = n_row;139 if (mat->n_col < n_col)140 mat->n_col = n_col;141 142 return mat;143error:144 isl_mat_free(mat);145 return NULL;146}147 148__isl_give isl_mat *isl_mat_sub_alloc6(isl_ctx *ctx, isl_int **row,149 unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)150{151 int i;152 struct isl_mat *mat;153 154 mat = isl_alloc_type(ctx, struct isl_mat);155 if (!mat)156 return NULL;157 mat->row = isl_alloc_array(ctx, isl_int *, n_row);158 if (n_row && !mat->row)159 goto error;160 for (i = 0; i < n_row; ++i)161 mat->row[i] = row[first_row+i] + first_col;162 mat->ctx = ctx;163 isl_ctx_ref(ctx);164 mat->ref = 1;165 mat->n_row = n_row;166 mat->n_col = n_col;167 mat->block = isl_blk_empty();168 mat->flags = ISL_MAT_BORROWED;169 return mat;170error:171 free(mat);172 return NULL;173}174 175__isl_give isl_mat *isl_mat_sub_alloc(__isl_keep isl_mat *mat,176 unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)177{178 if (!mat)179 return NULL;180 return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,181 first_col, n_col);182}183 184void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,185 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)186{187 int i;188 189 for (i = 0; i < n_row; ++i)190 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);191}192 193void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,194 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)195{196 int i;197 198 for (i = 0; i < n_row; ++i)199 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);200}201 202__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)203{204 if (!mat)205 return NULL;206 207 mat->ref++;208 return mat;209}210 211__isl_give isl_mat *isl_mat_dup(__isl_keep isl_mat *mat)212{213 int i;214 struct isl_mat *mat2;215 216 if (!mat)217 return NULL;218 mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);219 if (!mat2)220 return NULL;221 for (i = 0; i < mat->n_row; ++i)222 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);223 return mat2;224}225 226__isl_give isl_mat *isl_mat_cow(__isl_take isl_mat *mat)227{228 struct isl_mat *mat2;229 if (!mat)230 return NULL;231 232 if (mat->ref == 1 && !ISL_F_ISSET(mat, ISL_MAT_BORROWED))233 return mat;234 235 mat2 = isl_mat_dup(mat);236 isl_mat_free(mat);237 return mat2;238}239 240__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)241{242 if (!mat)243 return NULL;244 245 if (--mat->ref > 0)246 return NULL;247 248 if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))249 isl_blk_free(mat->ctx, mat->block);250 isl_ctx_deref(mat->ctx);251 free(mat->row);252 free(mat);253 254 return NULL;255}256 257isl_size isl_mat_rows(__isl_keep isl_mat *mat)258{259 return mat ? mat->n_row : isl_size_error;260}261 262isl_size isl_mat_cols(__isl_keep isl_mat *mat)263{264 return mat ? mat->n_col : isl_size_error;265}266 267/* Check that "col" is a valid column position for "mat".268 */269static isl_stat check_col(__isl_keep isl_mat *mat, int col)270{271 if (!mat)272 return isl_stat_error;273 if (col < 0 || col >= mat->n_col)274 isl_die(isl_mat_get_ctx(mat), isl_error_invalid,275 "column out of range", return isl_stat_error);276 return isl_stat_ok;277}278 279/* Check that "row" is a valid row position for "mat".280 */281static isl_stat check_row(__isl_keep isl_mat *mat, int row)282{283 if (!mat)284 return isl_stat_error;285 if (row < 0 || row >= mat->n_row)286 isl_die(isl_mat_get_ctx(mat), isl_error_invalid,287 "row out of range", return isl_stat_error);288 return isl_stat_ok;289}290 291/* Check that there are "n" columns starting at position "first" in "mat".292 */293static isl_stat check_col_range(__isl_keep isl_mat *mat, unsigned first,294 unsigned n)295{296 if (!mat)297 return isl_stat_error;298 if (first + n > mat->n_col || first + n < first)299 isl_die(isl_mat_get_ctx(mat), isl_error_invalid,300 "column position or range out of bounds",301 return isl_stat_error);302 return isl_stat_ok;303}304 305/* Check that there are "n" rows starting at position "first" in "mat".306 */307static isl_stat check_row_range(__isl_keep isl_mat *mat, unsigned first,308 unsigned n)309{310 if (!mat)311 return isl_stat_error;312 if (first + n > mat->n_row || first + n < first)313 isl_die(isl_mat_get_ctx(mat), isl_error_invalid,314 "row position or range out of bounds",315 return isl_stat_error);316 return isl_stat_ok;317}318 319int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)320{321 if (check_row(mat, row) < 0)322 return -1;323 if (check_col(mat, col) < 0)324 return -1;325 isl_int_set(*v, mat->row[row][col]);326 return 0;327}328 329/* Extract the element at row "row", oolumn "col" of "mat".330 */331__isl_give isl_val *isl_mat_get_element_val(__isl_keep isl_mat *mat,332 int row, int col)333{334 isl_ctx *ctx;335 336 if (check_row(mat, row) < 0)337 return NULL;338 if (check_col(mat, col) < 0)339 return NULL;340 ctx = isl_mat_get_ctx(mat);341 return isl_val_int_from_isl_int(ctx, mat->row[row][col]);342}343 344__isl_give isl_mat *isl_mat_set_element(__isl_take isl_mat *mat,345 int row, int col, isl_int v)346{347 mat = isl_mat_cow(mat);348 if (check_row(mat, row) < 0)349 return isl_mat_free(mat);350 if (check_col(mat, col) < 0)351 return isl_mat_free(mat);352 isl_int_set(mat->row[row][col], v);353 return mat;354}355 356__isl_give isl_mat *isl_mat_set_element_si(__isl_take isl_mat *mat,357 int row, int col, int v)358{359 mat = isl_mat_cow(mat);360 if (check_row(mat, row) < 0)361 return isl_mat_free(mat);362 if (check_col(mat, col) < 0)363 return isl_mat_free(mat);364 isl_int_set_si(mat->row[row][col], v);365 return mat;366}367 368/* Replace the element at row "row", column "col" of "mat" by "v".369 */370__isl_give isl_mat *isl_mat_set_element_val(__isl_take isl_mat *mat,371 int row, int col, __isl_take isl_val *v)372{373 if (!v)374 return isl_mat_free(mat);375 if (!isl_val_is_int(v))376 isl_die(isl_val_get_ctx(v), isl_error_invalid,377 "expecting integer value", goto error);378 mat = isl_mat_set_element(mat, row, col, v->n);379 isl_val_free(v);380 return mat;381error:382 isl_val_free(v);383 return isl_mat_free(mat);384}385 386__isl_give isl_mat *isl_mat_diag(isl_ctx *ctx, unsigned n_row, isl_int d)387{388 int i;389 struct isl_mat *mat;390 391 mat = isl_mat_alloc(ctx, n_row, n_row);392 if (!mat)393 return NULL;394 for (i = 0; i < n_row; ++i) {395 isl_seq_clr(mat->row[i], i);396 isl_int_set(mat->row[i][i], d);397 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));398 }399 400 return mat;401}402 403/* Create an "n_row" by "n_col" matrix with zero elements.404 */405__isl_give isl_mat *isl_mat_zero(isl_ctx *ctx, unsigned n_row, unsigned n_col)406{407 int i;408 isl_mat *mat;409 410 mat = isl_mat_alloc(ctx, n_row, n_col);411 if (!mat)412 return NULL;413 for (i = 0; i < n_row; ++i)414 isl_seq_clr(mat->row[i], n_col);415 416 return mat;417}418 419__isl_give isl_mat *isl_mat_identity(isl_ctx *ctx, unsigned n_row)420{421 if (!ctx)422 return NULL;423 return isl_mat_diag(ctx, n_row, ctx->one);424}425 426/* Is "mat" a (possibly scaled) identity matrix?427 */428isl_bool isl_mat_is_scaled_identity(__isl_keep isl_mat *mat)429{430 int i;431 432 if (!mat)433 return isl_bool_error;434 if (mat->n_row != mat->n_col)435 return isl_bool_false;436 437 for (i = 0; i < mat->n_row; ++i) {438 if (isl_seq_first_non_zero(mat->row[i], i) != -1)439 return isl_bool_false;440 if (isl_int_ne(mat->row[0][0], mat->row[i][i]))441 return isl_bool_false;442 if (isl_seq_first_non_zero(mat->row[i] + i + 1,443 mat->n_col - (i + 1)) != -1)444 return isl_bool_false;445 }446 447 return isl_bool_true;448}449 450__isl_give isl_vec *isl_mat_vec_product(__isl_take isl_mat *mat,451 __isl_take isl_vec *vec)452{453 int i;454 struct isl_vec *prod;455 456 if (!mat || !vec)457 goto error;458 459 isl_assert(mat->ctx, mat->n_col == vec->size, goto error);460 461 prod = isl_vec_alloc(mat->ctx, mat->n_row);462 if (!prod)463 goto error;464 465 for (i = 0; i < prod->size; ++i)466 isl_seq_inner_product(mat->row[i], vec->el, vec->size,467 &prod->block.data[i]);468 isl_mat_free(mat);469 isl_vec_free(vec);470 return prod;471error:472 isl_mat_free(mat);473 isl_vec_free(vec);474 return NULL;475}476 477__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,478 __isl_take isl_vec *vec)479{480 struct isl_mat *vec_mat;481 int i;482 483 if (!mat || !vec)484 goto error;485 vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);486 if (!vec_mat)487 goto error;488 for (i = 0; i < vec->size; ++i)489 isl_int_set(vec_mat->row[i][0], vec->el[i]);490 vec_mat = isl_mat_inverse_product(mat, vec_mat);491 isl_vec_free(vec);492 if (!vec_mat)493 return NULL;494 vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);495 if (vec)496 for (i = 0; i < vec->size; ++i)497 isl_int_set(vec->el[i], vec_mat->row[i][0]);498 isl_mat_free(vec_mat);499 return vec;500error:501 isl_mat_free(mat);502 isl_vec_free(vec);503 return NULL;504}505 506__isl_give isl_vec *isl_vec_mat_product(__isl_take isl_vec *vec,507 __isl_take isl_mat *mat)508{509 int i, j;510 struct isl_vec *prod;511 512 if (!mat || !vec)513 goto error;514 515 isl_assert(mat->ctx, mat->n_row == vec->size, goto error);516 517 prod = isl_vec_alloc(mat->ctx, mat->n_col);518 if (!prod)519 goto error;520 521 for (i = 0; i < prod->size; ++i) {522 isl_int_set_si(prod->el[i], 0);523 for (j = 0; j < vec->size; ++j)524 isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);525 }526 isl_mat_free(mat);527 isl_vec_free(vec);528 return prod;529error:530 isl_mat_free(mat);531 isl_vec_free(vec);532 return NULL;533}534 535__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,536 __isl_take isl_mat *right)537{538 int i;539 struct isl_mat *sum;540 541 if (!left || !right)542 goto error;543 544 isl_assert(left->ctx, left->n_row == right->n_row, goto error);545 isl_assert(left->ctx, left->n_row >= 1, goto error);546 isl_assert(left->ctx, left->n_col >= 1, goto error);547 isl_assert(left->ctx, right->n_col >= 1, goto error);548 isl_assert(left->ctx,549 isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,550 goto error);551 isl_assert(left->ctx,552 isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,553 goto error);554 555 sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);556 if (!sum)557 goto error;558 isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);559 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);560 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);561 562 isl_seq_clr(sum->row[0]+1, sum->n_col-1);563 for (i = 1; i < sum->n_row; ++i) {564 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);565 isl_int_addmul(sum->row[i][0],566 right->row[0][0], right->row[i][0]);567 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],568 left->n_col-1);569 isl_seq_scale(sum->row[i]+left->n_col,570 right->row[i]+1, right->row[0][0],571 right->n_col-1);572 }573 574 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);575 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);576 isl_mat_free(left);577 isl_mat_free(right);578 return sum;579error:580 isl_mat_free(left);581 isl_mat_free(right);582 return NULL;583}584 585static void exchange(__isl_keep isl_mat *M, __isl_keep isl_mat **U,586 __isl_keep isl_mat **Q, unsigned row, unsigned i, unsigned j)587{588 int r;589 for (r = row; r < M->n_row; ++r)590 isl_int_swap(M->row[r][i], M->row[r][j]);591 if (U) {592 for (r = 0; r < (*U)->n_row; ++r)593 isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);594 }595 if (Q)596 isl_mat_swap_rows(*Q, i, j);597}598 599static void subtract(__isl_keep isl_mat *M, __isl_keep isl_mat **U,600 __isl_keep isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)601{602 int r;603 for (r = row; r < M->n_row; ++r)604 isl_int_submul(M->row[r][j], m, M->row[r][i]);605 if (U) {606 for (r = 0; r < (*U)->n_row; ++r)607 isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);608 }609 if (Q) {610 for (r = 0; r < (*Q)->n_col; ++r)611 isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);612 }613}614 615static void oppose(__isl_keep isl_mat *M, __isl_keep isl_mat **U,616 __isl_keep isl_mat **Q, unsigned row, unsigned col)617{618 int r;619 for (r = row; r < M->n_row; ++r)620 isl_int_neg(M->row[r][col], M->row[r][col]);621 if (U) {622 for (r = 0; r < (*U)->n_row; ++r)623 isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);624 }625 if (Q)626 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);627}628 629/* Given matrix M, compute630 *631 * M U = H632 * M = H Q633 *634 * with U and Q unimodular matrices and H a matrix in column echelon form635 * such that on each echelon row the entries in the non-echelon column636 * are non-negative (if neg == 0) or non-positive (if neg == 1)637 * and strictly smaller (in absolute value) than the entries in the echelon638 * column.639 * If U or Q are NULL, then these matrices are not computed.640 */641__isl_give isl_mat *isl_mat_left_hermite(__isl_take isl_mat *M, int neg,642 __isl_give isl_mat **U, __isl_give isl_mat **Q)643{644 isl_int c;645 int row, col;646 647 if (U)648 *U = NULL;649 if (Q)650 *Q = NULL;651 if (!M)652 goto error;653 if (U) {654 *U = isl_mat_identity(M->ctx, M->n_col);655 if (!*U)656 goto error;657 }658 if (Q) {659 *Q = isl_mat_identity(M->ctx, M->n_col);660 if (!*Q)661 goto error;662 }663 664 if (M->n_col == 0)665 return M;666 667 M = isl_mat_cow(M);668 if (!M)669 goto error;670 671 col = 0;672 isl_int_init(c);673 for (row = 0; row < M->n_row; ++row) {674 int first, i, off;675 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);676 if (first == -1)677 continue;678 first += col;679 if (first != col)680 exchange(M, U, Q, row, first, col);681 if (isl_int_is_neg(M->row[row][col]))682 oppose(M, U, Q, row, col);683 first = col+1;684 while ((off = isl_seq_first_non_zero(M->row[row]+first,685 M->n_col-first)) != -1) {686 first += off;687 isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);688 subtract(M, U, Q, row, col, first, c);689 if (!isl_int_is_zero(M->row[row][first]))690 exchange(M, U, Q, row, first, col);691 else692 ++first;693 }694 for (i = 0; i < col; ++i) {695 if (isl_int_is_zero(M->row[row][i]))696 continue;697 if (neg)698 isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);699 else700 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);701 if (isl_int_is_zero(c))702 continue;703 subtract(M, U, Q, row, col, i, c);704 }705 ++col;706 }707 isl_int_clear(c);708 709 return M;710error:711 if (Q) {712 isl_mat_free(*Q);713 *Q = NULL;714 }715 if (U) {716 isl_mat_free(*U);717 *U = NULL;718 }719 isl_mat_free(M);720 return NULL;721}722 723/* Use row "row" of "mat" to eliminate column "col" from all other rows.724 */725static __isl_give isl_mat *eliminate(__isl_take isl_mat *mat, int row, int col)726{727 int k;728 isl_size nr, nc;729 isl_ctx *ctx;730 731 nr = isl_mat_rows(mat);732 nc = isl_mat_cols(mat);733 if (nr < 0 || nc < 0)734 return isl_mat_free(mat);735 736 ctx = isl_mat_get_ctx(mat);737 738 for (k = 0; k < nr; ++k) {739 if (k == row)740 continue;741 if (isl_int_is_zero(mat->row[k][col]))742 continue;743 mat = isl_mat_cow(mat);744 if (!mat)745 return NULL;746 isl_seq_elim(mat->row[k], mat->row[row], col, nc, NULL);747 isl_seq_normalize(ctx, mat->row[k], nc);748 }749 750 return mat;751}752 753/* Perform Gaussian elimination on the rows of "mat", but start754 * from the final row and the final column.755 * Any zero rows that result from the elimination are removed.756 *757 * In particular, for each column from last to first,758 * look for the last row with a non-zero coefficient in that column,759 * move it last (but before other rows moved last in previous steps) and760 * use it to eliminate the column from the other rows.761 */762__isl_give isl_mat *isl_mat_reverse_gauss(__isl_take isl_mat *mat)763{764 int k, row, last;765 isl_size nr, nc;766 767 nr = isl_mat_rows(mat);768 nc = isl_mat_cols(mat);769 if (nr < 0 || nc < 0)770 return isl_mat_free(mat);771 772 last = nc - 1;773 for (row = nr - 1; row >= 0; --row) {774 for (; last >= 0; --last) {775 for (k = row; k >= 0; --k)776 if (!isl_int_is_zero(mat->row[k][last]))777 break;778 if (k >= 0)779 break;780 }781 if (last < 0)782 break;783 if (k != row)784 mat = isl_mat_swap_rows(mat, k, row);785 if (!mat)786 return NULL;787 if (isl_int_is_neg(mat->row[row][last]))788 mat = isl_mat_row_neg(mat, row);789 mat = eliminate(mat, row, last);790 if (!mat)791 return NULL;792 }793 mat = isl_mat_drop_rows(mat, 0, row + 1);794 795 return mat;796}797 798/* Negate the lexicographically negative rows of "mat" such that799 * all rows in the result are lexicographically non-negative.800 */801__isl_give isl_mat *isl_mat_lexnonneg_rows(__isl_take isl_mat *mat)802{803 int i;804 isl_size nr, nc;805 806 nr = isl_mat_rows(mat);807 nc = isl_mat_cols(mat);808 if (nr < 0 || nc < 0)809 return isl_mat_free(mat);810 811 for (i = 0; i < nr; ++i) {812 int pos;813 814 pos = isl_seq_first_non_zero(mat->row[i], nc);815 if (pos < 0)816 continue;817 if (isl_int_is_nonneg(mat->row[i][pos]))818 continue;819 mat = isl_mat_row_neg(mat, i);820 if (!mat)821 return NULL;822 }823 824 return mat;825}826 827/* Given a matrix "H" is column echelon form, what is the first828 * zero column? That is how many initial columns are non-zero?829 * Start looking at column "first_col" and only consider830 * the columns to be of size "n_row".831 * "H" is assumed to be non-NULL.832 *833 * Since "H" is in column echelon form, the first non-zero entry834 * in a column is always in a later position compared to the previous column.835 */836static int hermite_first_zero_col(__isl_keep isl_mat *H, int first_col,837 int n_row)838{839 int row, col;840 841 for (col = first_col, row = 0; col < H->n_col; ++col) {842 for (; row < n_row; ++row)843 if (!isl_int_is_zero(H->row[row][col]))844 break;845 if (row == n_row)846 return col;847 }848 849 return H->n_col;850}851 852/* Return the rank of "mat", or isl_size_error in case of error.853 */854isl_size isl_mat_rank(__isl_keep isl_mat *mat)855{856 int rank;857 isl_mat *H;858 859 H = isl_mat_left_hermite(isl_mat_copy(mat), 0, NULL, NULL);860 if (!H)861 return isl_size_error;862 863 rank = hermite_first_zero_col(H, 0, H->n_row);864 isl_mat_free(H);865 866 return rank;867}868 869__isl_give isl_mat *isl_mat_right_kernel(__isl_take isl_mat *mat)870{871 int rank;872 struct isl_mat *U = NULL;873 struct isl_mat *K;874 875 mat = isl_mat_left_hermite(mat, 0, &U, NULL);876 if (!mat || !U)877 goto error;878 879 rank = hermite_first_zero_col(mat, 0, mat->n_row);880 K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank);881 if (!K)882 goto error;883 isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);884 isl_mat_free(mat);885 isl_mat_free(U);886 return K;887error:888 isl_mat_free(mat);889 isl_mat_free(U);890 return NULL;891}892 893__isl_give isl_mat *isl_mat_lin_to_aff(__isl_take isl_mat *mat)894{895 int i;896 struct isl_mat *mat2;897 898 if (!mat)899 return NULL;900 mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);901 if (!mat2)902 goto error;903 isl_int_set_si(mat2->row[0][0], 1);904 isl_seq_clr(mat2->row[0]+1, mat->n_col);905 for (i = 0; i < mat->n_row; ++i) {906 isl_int_set_si(mat2->row[1+i][0], 0);907 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);908 }909 isl_mat_free(mat);910 return mat2;911error:912 isl_mat_free(mat);913 return NULL;914}915 916/* Given two matrices M1 and M2, return the block matrix917 *918 * [ M1 0 ]919 * [ 0 M2 ]920 */921__isl_give isl_mat *isl_mat_diagonal(__isl_take isl_mat *mat1,922 __isl_take isl_mat *mat2)923{924 int i;925 isl_mat *mat;926 927 if (!mat1 || !mat2)928 goto error;929 930 mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,931 mat1->n_col + mat2->n_col);932 if (!mat)933 goto error;934 for (i = 0; i < mat1->n_row; ++i) {935 isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);936 isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);937 }938 for (i = 0; i < mat2->n_row; ++i) {939 isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);940 isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,941 mat2->row[i], mat2->n_col);942 }943 isl_mat_free(mat1);944 isl_mat_free(mat2);945 return mat;946error:947 isl_mat_free(mat1);948 isl_mat_free(mat2);949 return NULL;950}951 952static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)953{954 int i;955 956 for (i = 0; i < n_row; ++i)957 if (!isl_int_is_zero(row[i][col]))958 return i;959 return -1;960}961 962static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)963{964 int i, min = row_first_non_zero(row, n_row, col);965 if (min < 0)966 return -1;967 for (i = min + 1; i < n_row; ++i) {968 if (isl_int_is_zero(row[i][col]))969 continue;970 if (isl_int_abs_lt(row[i][col], row[min][col]))971 min = i;972 }973 return min;974}975 976static isl_stat inv_exchange(__isl_keep isl_mat **left,977 __isl_keep isl_mat **right, unsigned i, unsigned j)978{979 *left = isl_mat_swap_rows(*left, i, j);980 *right = isl_mat_swap_rows(*right, i, j);981 982 if (!*left || !*right)983 return isl_stat_error;984 return isl_stat_ok;985}986 987static void inv_oppose(988 __isl_keep isl_mat *left, __isl_keep isl_mat *right, unsigned row)989{990 isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);991 isl_seq_neg(right->row[row], right->row[row], right->n_col);992}993 994static void inv_subtract(__isl_keep isl_mat *left, __isl_keep isl_mat *right,995 unsigned row, unsigned i, isl_int m)996{997 isl_int_neg(m, m);998 isl_seq_combine(left->row[i]+row,999 left->ctx->one, left->row[i]+row,1000 m, left->row[row]+row,1001 left->n_col-row);1002 isl_seq_combine(right->row[i], right->ctx->one, right->row[i],1003 m, right->row[row], right->n_col);1004}1005 1006/* Compute inv(left)*right1007 */1008__isl_give isl_mat *isl_mat_inverse_product(__isl_take isl_mat *left,1009 __isl_take isl_mat *right)1010{1011 int row;1012 isl_int a, b;1013 1014 if (!left || !right)1015 goto error;1016 1017 isl_assert(left->ctx, left->n_row == left->n_col, goto error);1018 isl_assert(left->ctx, left->n_row == right->n_row, goto error);1019 1020 if (left->n_row == 0) {1021 isl_mat_free(left);1022 return right;1023 }1024 1025 left = isl_mat_cow(left);1026 right = isl_mat_cow(right);1027 if (!left || !right)1028 goto error;1029 1030 isl_int_init(a);1031 isl_int_init(b);1032 for (row = 0; row < left->n_row; ++row) {1033 int pivot, first, i, off;1034 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);1035 if (pivot < 0) {1036 isl_int_clear(a);1037 isl_int_clear(b);1038 isl_assert(left->ctx, pivot >= 0, goto error);1039 }1040 pivot += row;1041 if (pivot != row)1042 if (inv_exchange(&left, &right, pivot, row) < 0)1043 goto error;1044 if (isl_int_is_neg(left->row[row][row]))1045 inv_oppose(left, right, row);1046 first = row+1;1047 while ((off = row_first_non_zero(left->row+first,1048 left->n_row-first, row)) != -1) {1049 first += off;1050 isl_int_fdiv_q(a, left->row[first][row],1051 left->row[row][row]);1052 inv_subtract(left, right, row, first, a);1053 if (!isl_int_is_zero(left->row[first][row])) {1054 if (inv_exchange(&left, &right, row, first) < 0)1055 goto error;1056 } else {1057 ++first;1058 }1059 }1060 for (i = 0; i < row; ++i) {1061 if (isl_int_is_zero(left->row[i][row]))1062 continue;1063 isl_int_gcd(a, left->row[row][row], left->row[i][row]);1064 isl_int_divexact(b, left->row[i][row], a);1065 isl_int_divexact(a, left->row[row][row], a);1066 isl_int_neg(b, b);1067 isl_seq_combine(left->row[i] + i,1068 a, left->row[i] + i,1069 b, left->row[row] + i,1070 left->n_col - i);1071 isl_seq_combine(right->row[i], a, right->row[i],1072 b, right->row[row], right->n_col);1073 }1074 }1075 isl_int_clear(b);1076 1077 isl_int_set(a, left->row[0][0]);1078 for (row = 1; row < left->n_row; ++row)1079 isl_int_lcm(a, a, left->row[row][row]);1080 if (isl_int_is_zero(a)){1081 isl_int_clear(a);1082 isl_assert(left->ctx, 0, goto error);1083 }1084 for (row = 0; row < left->n_row; ++row) {1085 isl_int_divexact(left->row[row][row], a, left->row[row][row]);1086 if (isl_int_is_one(left->row[row][row]))1087 continue;1088 isl_seq_scale(right->row[row], right->row[row],1089 left->row[row][row], right->n_col);1090 }1091 isl_int_clear(a);1092 1093 isl_mat_free(left);1094 return right;1095error:1096 isl_mat_free(left);1097 isl_mat_free(right);1098 return NULL;1099}1100 1101void isl_mat_col_scale(__isl_keep isl_mat *mat, unsigned col, isl_int m)1102{1103 int i;1104 1105 for (i = 0; i < mat->n_row; ++i)1106 isl_int_mul(mat->row[i][col], mat->row[i][col], m);1107}1108 1109void isl_mat_col_combine(__isl_keep isl_mat *mat, unsigned dst,1110 isl_int m1, unsigned src1, isl_int m2, unsigned src2)1111{1112 int i;1113 isl_int tmp;1114 1115 isl_int_init(tmp);1116 for (i = 0; i < mat->n_row; ++i) {1117 isl_int_mul(tmp, m1, mat->row[i][src1]);1118 isl_int_addmul(tmp, m2, mat->row[i][src2]);1119 isl_int_set(mat->row[i][dst], tmp);1120 }1121 isl_int_clear(tmp);1122}1123 1124__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)1125{1126 struct isl_mat *inv;1127 int row;1128 isl_int a, b;1129 1130 mat = isl_mat_cow(mat);1131 if (!mat)1132 return NULL;1133 1134 inv = isl_mat_identity(mat->ctx, mat->n_col);1135 inv = isl_mat_cow(inv);1136 if (!inv)1137 goto error;1138 1139 isl_int_init(a);1140 isl_int_init(b);1141 for (row = 0; row < mat->n_row; ++row) {1142 int pivot, first, i, off;1143 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);1144 if (pivot < 0) {1145 isl_int_clear(a);1146 isl_int_clear(b);1147 isl_assert(mat->ctx, pivot >= 0, goto error);1148 }1149 pivot += row;1150 if (pivot != row)1151 exchange(mat, &inv, NULL, row, pivot, row);1152 if (isl_int_is_neg(mat->row[row][row]))1153 oppose(mat, &inv, NULL, row, row);1154 first = row+1;1155 while ((off = isl_seq_first_non_zero(mat->row[row]+first,1156 mat->n_col-first)) != -1) {1157 first += off;1158 isl_int_fdiv_q(a, mat->row[row][first],1159 mat->row[row][row]);1160 subtract(mat, &inv, NULL, row, row, first, a);1161 if (!isl_int_is_zero(mat->row[row][first]))1162 exchange(mat, &inv, NULL, row, row, first);1163 else1164 ++first;1165 }1166 for (i = 0; i < row; ++i) {1167 if (isl_int_is_zero(mat->row[row][i]))1168 continue;1169 isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);1170 isl_int_divexact(b, mat->row[row][i], a);1171 isl_int_divexact(a, mat->row[row][row], a);1172 isl_int_neg(a, a);1173 isl_mat_col_combine(mat, i, a, i, b, row);1174 isl_mat_col_combine(inv, i, a, i, b, row);1175 }1176 }1177 isl_int_clear(b);1178 1179 isl_int_set(a, mat->row[0][0]);1180 for (row = 1; row < mat->n_row; ++row)1181 isl_int_lcm(a, a, mat->row[row][row]);1182 if (isl_int_is_zero(a)){1183 isl_int_clear(a);1184 goto error;1185 }1186 for (row = 0; row < mat->n_row; ++row) {1187 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);1188 if (isl_int_is_one(mat->row[row][row]))1189 continue;1190 isl_mat_col_scale(inv, row, mat->row[row][row]);1191 }1192 isl_int_clear(a);1193 1194 isl_mat_free(mat);1195 1196 return inv;1197error:1198 isl_mat_free(mat);1199 isl_mat_free(inv);1200 return NULL;1201}1202 1203__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)1204{1205 struct isl_mat *transpose = NULL;1206 int i, j;1207 1208 if (!mat)1209 return NULL;1210 1211 if (mat->n_col == mat->n_row) {1212 mat = isl_mat_cow(mat);1213 if (!mat)1214 return NULL;1215 for (i = 0; i < mat->n_row; ++i)1216 for (j = i + 1; j < mat->n_col; ++j)1217 isl_int_swap(mat->row[i][j], mat->row[j][i]);1218 return mat;1219 }1220 transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);1221 if (!transpose)1222 goto error;1223 for (i = 0; i < mat->n_row; ++i)1224 for (j = 0; j < mat->n_col; ++j)1225 isl_int_set(transpose->row[j][i], mat->row[i][j]);1226 isl_mat_free(mat);1227 return transpose;1228error:1229 isl_mat_free(mat);1230 return NULL;1231}1232 1233__isl_give isl_mat *isl_mat_swap_cols(__isl_take isl_mat *mat,1234 unsigned i, unsigned j)1235{1236 int r;1237 1238 mat = isl_mat_cow(mat);1239 if (check_col_range(mat, i, 1) < 0 ||1240 check_col_range(mat, j, 1) < 0)1241 return isl_mat_free(mat);1242 1243 for (r = 0; r < mat->n_row; ++r)1244 isl_int_swap(mat->row[r][i], mat->row[r][j]);1245 return mat;1246}1247 1248__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,1249 unsigned i, unsigned j)1250{1251 isl_int *t;1252 1253 if (!mat)1254 return NULL;1255 mat = isl_mat_cow(mat);1256 if (check_row_range(mat, i, 1) < 0 ||1257 check_row_range(mat, j, 1) < 0)1258 return isl_mat_free(mat);1259 1260 t = mat->row[i];1261 mat->row[i] = mat->row[j];1262 mat->row[j] = t;1263 return mat;1264}1265 1266/* Calculate the product of two matrices.1267 *1268 * This function is optimized for operand matrices that contain many zeros and1269 * skips multiplications where we know one of the operands is zero.1270 */1271__isl_give isl_mat *isl_mat_product(__isl_take isl_mat *left,1272 __isl_take isl_mat *right)1273{1274 int i, j, k;1275 struct isl_mat *prod;1276 1277 if (!left || !right)1278 goto error;1279 isl_assert(left->ctx, left->n_col == right->n_row, goto error);1280 prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);1281 if (!prod)1282 goto error;1283 if (left->n_col == 0) {1284 for (i = 0; i < prod->n_row; ++i)1285 isl_seq_clr(prod->row[i], prod->n_col);1286 isl_mat_free(left);1287 isl_mat_free(right);1288 return prod;1289 }1290 for (i = 0; i < prod->n_row; ++i) {1291 for (j = 0; j < prod->n_col; ++j)1292 isl_int_mul(prod->row[i][j],1293 left->row[i][0], right->row[0][j]);1294 for (k = 1; k < left->n_col; ++k) {1295 if (isl_int_is_zero(left->row[i][k]))1296 continue;1297 for (j = 0; j < prod->n_col; ++j)1298 isl_int_addmul(prod->row[i][j],1299 left->row[i][k], right->row[k][j]);1300 }1301 }1302 isl_mat_free(left);1303 isl_mat_free(right);1304 return prod;1305error:1306 isl_mat_free(left);1307 isl_mat_free(right);1308 return NULL;1309}1310 1311/* Replace the variables x in the rows q by x' given by x = M x',1312 * with M the matrix mat.1313 *1314 * If the number of new variables is greater than the original1315 * number of variables, then the rows q have already been1316 * preextended. If the new number is smaller, then the coefficients1317 * of the divs, which are not changed, need to be shifted down.1318 * The row q may be the equalities, the inequalities or the1319 * div expressions. In the latter case, has_div is true and1320 * we need to take into account the extra denominator column.1321 */1322static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,1323 unsigned n_div, int has_div, struct isl_mat *mat)1324{1325 int i;1326 struct isl_mat *t;1327 int e;1328 1329 if (mat->n_col >= mat->n_row)1330 e = 0;1331 else1332 e = mat->n_row - mat->n_col;1333 if (has_div)1334 for (i = 0; i < n; ++i)1335 isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);1336 t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);1337 t = isl_mat_product(t, mat);1338 if (!t)1339 return -1;1340 for (i = 0; i < n; ++i) {1341 isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);1342 isl_seq_cpy(q[i] + has_div + t->n_col,1343 q[i] + has_div + t->n_col + e, n_div);1344 isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);1345 }1346 isl_mat_free(t);1347 return 0;1348}1349 1350/* Replace the variables x in bset by x' given by x = M x', with1351 * M the matrix mat.1352 *1353 * If there are fewer variables x' then there are x, then we perform1354 * the transformation in place, which means that, in principle,1355 * this frees up some extra variables as the number1356 * of columns remains constant, but we would have to extend1357 * the div array too as the number of rows in this array is assumed1358 * to be equal to extra.1359 */1360__isl_give isl_basic_set *isl_basic_set_preimage(1361 __isl_take isl_basic_set *bset, __isl_take isl_mat *mat)1362{1363 struct isl_ctx *ctx;1364 1365 if (!bset || !mat)1366 goto error;1367 1368 ctx = bset->ctx;1369 bset = isl_basic_set_cow(bset);1370 if (isl_basic_set_check_no_params(bset) < 0)1371 goto error;1372 1373 isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);1374 isl_assert(ctx, mat->n_col > 0, goto error);1375 1376 if (mat->n_col > mat->n_row) {1377 bset = isl_basic_set_add_dims(bset, isl_dim_set,1378 mat->n_col - mat->n_row);1379 if (!bset)1380 goto error;1381 } else if (mat->n_col < mat->n_row) {1382 bset->dim = isl_space_cow(bset->dim);1383 if (!bset->dim)1384 goto error;1385 bset->dim->n_out -= mat->n_row - mat->n_col;1386 }1387 1388 if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,1389 isl_mat_copy(mat)) < 0)1390 goto error;1391 1392 if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,1393 isl_mat_copy(mat)) < 0)1394 goto error;1395 1396 if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)1397 goto error2;1398 1399 ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);1400 ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);1401 ISL_F_CLR(bset, ISL_BASIC_SET_SORTED);1402 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);1403 ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);1404 1405 bset = isl_basic_set_simplify(bset);1406 bset = isl_basic_set_finalize(bset);1407 1408 return bset;1409error:1410 isl_mat_free(mat);1411error2:1412 isl_basic_set_free(bset);1413 return NULL;1414}1415 1416__isl_give isl_set *isl_set_preimage(1417 __isl_take isl_set *set, __isl_take isl_mat *mat)1418{1419 int i;1420 1421 set = isl_set_cow(set);1422 if (!set)1423 goto error;1424 1425 for (i = 0; i < set->n; ++i) {1426 set->p[i] = isl_basic_set_preimage(set->p[i],1427 isl_mat_copy(mat));1428 if (!set->p[i])1429 goto error;1430 }1431 if (mat->n_col != mat->n_row) {1432 set->dim = isl_space_cow(set->dim);1433 if (!set->dim)1434 goto error;1435 set->dim->n_out += mat->n_col;1436 set->dim->n_out -= mat->n_row;1437 }1438 isl_mat_free(mat);1439 ISL_F_CLR(set, ISL_SET_NORMALIZED);1440 return set;1441error:1442 isl_set_free(set);1443 isl_mat_free(mat);1444 return NULL;1445}1446 1447/* Replace the variables x starting at "first_col" in the rows "rows"1448 * of some coefficient matrix by x' with x = M x' with M the matrix mat.1449 * That is, replace the corresponding coefficients c by c M.1450 */1451isl_stat isl_mat_sub_transform(isl_int **row, unsigned n_row,1452 unsigned first_col, __isl_take isl_mat *mat)1453{1454 int i;1455 isl_ctx *ctx;1456 isl_mat *t;1457 1458 if (!mat)1459 return isl_stat_error;1460 ctx = isl_mat_get_ctx(mat);1461 t = isl_mat_sub_alloc6(ctx, row, 0, n_row, first_col, mat->n_row);1462 t = isl_mat_product(t, mat);1463 if (!t)1464 return isl_stat_error;1465 for (i = 0; i < n_row; ++i)1466 isl_seq_swp_or_cpy(row[i] + first_col, t->row[i], t->n_col);1467 isl_mat_free(t);1468 return isl_stat_ok;1469}1470 1471void isl_mat_print_internal(__isl_keep isl_mat *mat, FILE *out, int indent)1472{1473 int i, j;1474 1475 if (!mat) {1476 fprintf(out, "%*snull mat\n", indent, "");1477 return;1478 }1479 1480 if (mat->n_row == 0)1481 fprintf(out, "%*s[]\n", indent, "");1482 1483 for (i = 0; i < mat->n_row; ++i) {1484 if (!i)1485 fprintf(out, "%*s[[", indent, "");1486 else1487 fprintf(out, "%*s[", indent+1, "");1488 for (j = 0; j < mat->n_col; ++j) {1489 if (j)1490 fprintf(out, ",");1491 isl_int_print(out, mat->row[i][j], 0);1492 }1493 if (i == mat->n_row-1)1494 fprintf(out, "]]\n");1495 else1496 fprintf(out, "]\n");1497 }1498}1499 1500void isl_mat_dump(__isl_keep isl_mat *mat)1501{1502 isl_mat_print_internal(mat, stderr, 0);1503}1504 1505__isl_give isl_mat *isl_mat_drop_cols(__isl_take isl_mat *mat,1506 unsigned col, unsigned n)1507{1508 int r;1509 1510 if (n == 0)1511 return mat;1512 1513 mat = isl_mat_cow(mat);1514 if (check_col_range(mat, col, n) < 0)1515 return isl_mat_free(mat);1516 1517 if (col != mat->n_col-n) {1518 for (r = 0; r < mat->n_row; ++r)1519 isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,1520 mat->n_col - col - n);1521 }1522 mat->n_col -= n;1523 return mat;1524}1525 1526__isl_give isl_mat *isl_mat_drop_rows(__isl_take isl_mat *mat,1527 unsigned row, unsigned n)1528{1529 int r;1530 1531 mat = isl_mat_cow(mat);1532 if (check_row_range(mat, row, n) < 0)1533 return isl_mat_free(mat);1534 1535 for (r = row; r+n < mat->n_row; ++r)1536 mat->row[r] = mat->row[r+n];1537 1538 mat->n_row -= n;1539 return mat;1540}1541 1542__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,1543 unsigned col, unsigned n)1544{1545 isl_mat *ext;1546 1547 if (check_col_range(mat, col, 0) < 0)1548 return isl_mat_free(mat);1549 if (n == 0)1550 return mat;1551 1552 ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);1553 if (!ext)1554 goto error;1555 1556 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);1557 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,1558 col + n, col, mat->n_col - col);1559 1560 isl_mat_free(mat);1561 return ext;1562error:1563 isl_mat_free(mat);1564 return NULL;1565}1566 1567__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,1568 unsigned first, unsigned n)1569{1570 int i;1571 1572 if (!mat)1573 return NULL;1574 mat = isl_mat_insert_cols(mat, first, n);1575 if (!mat)1576 return NULL;1577 1578 for (i = 0; i < mat->n_row; ++i)1579 isl_seq_clr(mat->row[i] + first, n);1580 1581 return mat;1582}1583 1584__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)1585{1586 if (!mat)1587 return NULL;1588 1589 return isl_mat_insert_zero_cols(mat, mat->n_col, n);1590}1591 1592__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,1593 unsigned row, unsigned n)1594{1595 isl_mat *ext;1596 1597 if (check_row_range(mat, row, 0) < 0)1598 return isl_mat_free(mat);1599 if (n == 0)1600 return mat;1601 1602 ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);1603 if (!ext)1604 goto error;1605 1606 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);1607 isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,1608 mat->n_row - row, 0, 0, mat->n_col);1609 1610 isl_mat_free(mat);1611 return ext;1612error:1613 isl_mat_free(mat);1614 return NULL;1615}1616 1617__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)1618{1619 if (!mat)1620 return NULL;1621 1622 return isl_mat_insert_rows(mat, mat->n_row, n);1623}1624 1625__isl_give isl_mat *isl_mat_insert_zero_rows(__isl_take isl_mat *mat,1626 unsigned row, unsigned n)1627{1628 int i;1629 1630 mat = isl_mat_insert_rows(mat, row, n);1631 if (!mat)1632 return NULL;1633 1634 for (i = 0; i < n; ++i)1635 isl_seq_clr(mat->row[row + i], mat->n_col);1636 1637 return mat;1638}1639 1640__isl_give isl_mat *isl_mat_add_zero_rows(__isl_take isl_mat *mat, unsigned n)1641{1642 if (!mat)1643 return NULL;1644 1645 return isl_mat_insert_zero_rows(mat, mat->n_row, n);1646}1647 1648void isl_mat_col_submul(__isl_keep isl_mat *mat,1649 int dst_col, isl_int f, int src_col)1650{1651 int i;1652 1653 for (i = 0; i < mat->n_row; ++i)1654 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);1655}1656 1657void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)1658{1659 int i;1660 1661 if (!mat)1662 return;1663 1664 for (i = 0; i < mat->n_row; ++i)1665 isl_int_add(mat->row[i][dst_col],1666 mat->row[i][dst_col], mat->row[i][src_col]);1667}1668 1669void isl_mat_col_mul(__isl_keep isl_mat *mat, int dst_col, isl_int f,1670 int src_col)1671{1672 int i;1673 1674 for (i = 0; i < mat->n_row; ++i)1675 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);1676}1677 1678/* Add "f" times column "src_col" to column "dst_col" of "mat" and1679 * return the result.1680 */1681__isl_give isl_mat *isl_mat_col_addmul(__isl_take isl_mat *mat, int dst_col,1682 isl_int f, int src_col)1683{1684 int i;1685 1686 if (check_col(mat, dst_col) < 0 || check_col(mat, src_col) < 0)1687 return isl_mat_free(mat);1688 1689 for (i = 0; i < mat->n_row; ++i) {1690 if (isl_int_is_zero(mat->row[i][src_col]))1691 continue;1692 mat = isl_mat_cow(mat);1693 if (!mat)1694 return NULL;1695 isl_int_addmul(mat->row[i][dst_col], f, mat->row[i][src_col]);1696 }1697 1698 return mat;1699}1700 1701/* Negate column "col" of "mat" and return the result.1702 */1703__isl_give isl_mat *isl_mat_col_neg(__isl_take isl_mat *mat, int col)1704{1705 int i;1706 1707 if (check_col(mat, col) < 0)1708 return isl_mat_free(mat);1709 1710 for (i = 0; i < mat->n_row; ++i) {1711 if (isl_int_is_zero(mat->row[i][col]))1712 continue;1713 mat = isl_mat_cow(mat);1714 if (!mat)1715 return NULL;1716 isl_int_neg(mat->row[i][col], mat->row[i][col]);1717 }1718 1719 return mat;1720}1721 1722/* Negate row "row" of "mat" and return the result.1723 */1724__isl_give isl_mat *isl_mat_row_neg(__isl_take isl_mat *mat, int row)1725{1726 if (check_row(mat, row) < 0)1727 return isl_mat_free(mat);1728 if (isl_seq_first_non_zero(mat->row[row], mat->n_col) == -1)1729 return mat;1730 mat = isl_mat_cow(mat);1731 if (!mat)1732 return NULL;1733 isl_seq_neg(mat->row[row], mat->row[row], mat->n_col);1734 return mat;1735}1736 1737__isl_give isl_mat *isl_mat_unimodular_complete(__isl_take isl_mat *M, int row)1738{1739 int r;1740 struct isl_mat *H = NULL, *Q = NULL;1741 1742 if (!M)1743 return NULL;1744 1745 isl_assert(M->ctx, M->n_row == M->n_col, goto error);1746 M->n_row = row;1747 H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);1748 M->n_row = M->n_col;1749 if (!H)1750 goto error;1751 for (r = 0; r < row; ++r)1752 isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);1753 for (r = row; r < M->n_row; ++r)1754 isl_seq_cpy(M->row[r], Q->row[r], M->n_col);1755 isl_mat_free(H);1756 isl_mat_free(Q);1757 return M;1758error:1759 isl_mat_free(H);1760 isl_mat_free(Q);1761 isl_mat_free(M);1762 return NULL;1763}1764 1765__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,1766 __isl_take isl_mat *bot)1767{1768 struct isl_mat *mat;1769 1770 if (!top || !bot)1771 goto error;1772 1773 isl_assert(top->ctx, top->n_col == bot->n_col, goto error);1774 if (top->n_row == 0) {1775 isl_mat_free(top);1776 return bot;1777 }1778 if (bot->n_row == 0) {1779 isl_mat_free(bot);1780 return top;1781 }1782 1783 mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);1784 if (!mat)1785 goto error;1786 isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,1787 0, 0, mat->n_col);1788 isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,1789 0, 0, mat->n_col);1790 isl_mat_free(top);1791 isl_mat_free(bot);1792 return mat;1793error:1794 isl_mat_free(top);1795 isl_mat_free(bot);1796 return NULL;1797}1798 1799isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)1800{1801 int i;1802 1803 if (!mat1 || !mat2)1804 return isl_bool_error;1805 1806 if (mat1->n_row != mat2->n_row)1807 return isl_bool_false;1808 1809 if (mat1->n_col != mat2->n_col)1810 return isl_bool_false;1811 1812 for (i = 0; i < mat1->n_row; ++i)1813 if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col))1814 return isl_bool_false;1815 1816 return isl_bool_true;1817}1818 1819__isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)1820{1821 struct isl_mat *mat;1822 1823 if (!vec)1824 return NULL;1825 mat = isl_mat_alloc(vec->ctx, 1, vec->size);1826 if (!mat)1827 goto error;1828 1829 isl_seq_cpy(mat->row[0], vec->el, vec->size);1830 1831 isl_vec_free(vec);1832 return mat;1833error:1834 isl_vec_free(vec);1835 return NULL;1836}1837 1838/* Return a copy of row "row" of "mat" as an isl_vec.1839 */1840__isl_give isl_vec *isl_mat_get_row(__isl_keep isl_mat *mat, unsigned row)1841{1842 isl_vec *v;1843 1844 if (!mat)1845 return NULL;1846 if (row >= mat->n_row)1847 isl_die(mat->ctx, isl_error_invalid, "row out of range",1848 return NULL);1849 1850 v = isl_vec_alloc(isl_mat_get_ctx(mat), mat->n_col);1851 if (!v)1852 return NULL;1853 isl_seq_cpy(v->el, mat->row[row], mat->n_col);1854 1855 return v;1856}1857 1858__isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,1859 __isl_take isl_vec *bot)1860{1861 return isl_mat_concat(top, isl_mat_from_row_vec(bot));1862}1863 1864__isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,1865 unsigned dst_col, unsigned src_col, unsigned n)1866{1867 isl_mat *res;1868 1869 if (!mat)1870 return NULL;1871 if (n == 0 || dst_col == src_col)1872 return mat;1873 1874 res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);1875 if (!res)1876 goto error;1877 1878 if (dst_col < src_col) {1879 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,1880 0, 0, dst_col);1881 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,1882 dst_col, src_col, n);1883 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,1884 dst_col + n, dst_col, src_col - dst_col);1885 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,1886 src_col + n, src_col + n,1887 res->n_col - src_col - n);1888 } else {1889 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,1890 0, 0, src_col);1891 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,1892 src_col, src_col + n, dst_col - src_col);1893 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,1894 dst_col, src_col, n);1895 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,1896 dst_col + n, dst_col + n,1897 res->n_col - dst_col - n);1898 }1899 isl_mat_free(mat);1900 1901 return res;1902error:1903 isl_mat_free(mat);1904 return NULL;1905}1906 1907/* Return the gcd of the elements in row "row" of "mat" in *gcd.1908 * Return isl_stat_ok on success and isl_stat_error on failure.1909 */1910isl_stat isl_mat_row_gcd(__isl_keep isl_mat *mat, int row, isl_int *gcd)1911{1912 if (check_row(mat, row) < 0)1913 return isl_stat_error;1914 1915 isl_seq_gcd(mat->row[row], mat->n_col, gcd);1916 1917 return isl_stat_ok;1918}1919 1920void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)1921{1922 int i;1923 isl_int g;1924 1925 isl_int_set_si(*gcd, 0);1926 if (!mat)1927 return;1928 1929 isl_int_init(g);1930 for (i = 0; i < mat->n_row; ++i) {1931 isl_seq_gcd(mat->row[i], mat->n_col, &g);1932 isl_int_gcd(*gcd, *gcd, g);1933 }1934 isl_int_clear(g);1935}1936 1937/* Return the result of scaling "mat" by a factor of "m".1938 */1939__isl_give isl_mat *isl_mat_scale(__isl_take isl_mat *mat, isl_int m)1940{1941 int i;1942 1943 if (isl_int_is_one(m))1944 return mat;1945 1946 mat = isl_mat_cow(mat);1947 if (!mat)1948 return NULL;1949 1950 for (i = 0; i < mat->n_row; ++i)1951 isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);1952 1953 return mat;1954}1955 1956__isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)1957{1958 int i;1959 1960 if (isl_int_is_one(m))1961 return mat;1962 1963 mat = isl_mat_cow(mat);1964 if (!mat)1965 return NULL;1966 1967 for (i = 0; i < mat->n_row; ++i)1968 isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);1969 1970 return mat;1971}1972 1973__isl_give isl_mat *isl_mat_scale_down_row(__isl_take isl_mat *mat, int row,1974 isl_int m)1975{1976 if (isl_int_is_one(m))1977 return mat;1978 1979 mat = isl_mat_cow(mat);1980 if (!mat)1981 return NULL;1982 1983 isl_seq_scale_down(mat->row[row], mat->row[row], m, mat->n_col);1984 1985 return mat;1986}1987 1988__isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)1989{1990 isl_int gcd;1991 1992 if (!mat)1993 return NULL;1994 1995 isl_int_init(gcd);1996 isl_mat_gcd(mat, &gcd);1997 mat = isl_mat_scale_down(mat, gcd);1998 isl_int_clear(gcd);1999 2000 return mat;2001}2002 2003__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)2004{2005 mat = isl_mat_cow(mat);2006 if (!mat)2007 return NULL;2008 2009 isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);2010 2011 return mat;2012}2013 2014/* Number of initial non-zero columns.2015 */2016int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)2017{2018 int i;2019 2020 if (!mat)2021 return -1;2022 2023 for (i = 0; i < mat->n_col; ++i)2024 if (row_first_non_zero(mat->row, mat->n_row, i) < 0)2025 break;2026 2027 return i;2028}2029 2030/* Return a basis for the space spanned by the rows of "mat".2031 * Any basis will do, so simply perform Gaussian elimination and2032 * remove the empty rows.2033 */2034__isl_give isl_mat *isl_mat_row_basis(__isl_take isl_mat *mat)2035{2036 return isl_mat_reverse_gauss(mat);2037}2038 2039/* Return rows that extend a basis of "mat1" to one2040 * that covers both "mat1" and "mat2".2041 * The Hermite normal form of the concatenation of the two matrices is2042 *2043 * [ Q1 ]2044 * [ M1 ] = [ H1 0 0 ] [ Q2 ]2045 * [ M2 ] = [ H2 H3 0 ] [ Q3 ]2046 *2047 * The number of columns in H1 and H3 determine the number of rows2048 * in Q1 and Q2. Q1 is a basis for M1, while Q2 extends this basis2049 * to also cover M2.2050 */2051__isl_give isl_mat *isl_mat_row_basis_extension(2052 __isl_take isl_mat *mat1, __isl_take isl_mat *mat2)2053{2054 isl_size n_row;2055 int r1, r;2056 isl_size n1;2057 isl_mat *H, *Q;2058 2059 n1 = isl_mat_rows(mat1);2060 H = isl_mat_concat(mat1, mat2);2061 H = isl_mat_left_hermite(H, 0, NULL, &Q);2062 if (n1 < 0 || !H || !Q)2063 goto error;2064 2065 r1 = hermite_first_zero_col(H, 0, n1);2066 r = hermite_first_zero_col(H, r1, H->n_row);2067 n_row = isl_mat_rows(Q);2068 if (n_row < 0)2069 goto error;2070 Q = isl_mat_drop_rows(Q, r, n_row - r);2071 Q = isl_mat_drop_rows(Q, 0, r1);2072 2073 isl_mat_free(H);2074 return Q;2075error:2076 isl_mat_free(H);2077 isl_mat_free(Q);2078 return NULL;2079}2080 2081/* Are the rows of "mat1" linearly independent of those of "mat2"?2082 * That is, is there no linear dependence among the combined rows2083 * that is not already present in either "mat1" or "mat2"?2084 * In other words, is the rank of "mat1" and "mat2" combined equal2085 * to the sum of the ranks of "mat1" and "mat2"?2086 */2087isl_bool isl_mat_has_linearly_independent_rows(__isl_keep isl_mat *mat1,2088 __isl_keep isl_mat *mat2)2089{2090 isl_size r1, r2, r;2091 isl_mat *mat;2092 2093 r1 = isl_mat_rank(mat1);2094 if (r1 < 0)2095 return isl_bool_error;2096 if (r1 == 0)2097 return isl_bool_true;2098 r2 = isl_mat_rank(mat2);2099 if (r2 < 0)2100 return isl_bool_error;2101 if (r2 == 0)2102 return isl_bool_true;2103 2104 mat = isl_mat_concat(isl_mat_copy(mat1), isl_mat_copy(mat2));2105 r = isl_mat_rank(mat);2106 isl_mat_free(mat);2107 if (r < 0)2108 return isl_bool_error;2109 return isl_bool_ok(r == r1 + r2);2110}2111