1/* 2 * Copyright 2008-2009 Katholieke Universiteit Leuven 3 * Copyright 2010 INRIA Saclay 4 * Copyright 2014 Ecole Normale Superieure 5 * Copyright 2017 Sven Verdoolaege 6 * 7 * Use of this software is governed by the MIT license 8 * 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 14 */ 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, compute 630 * 631 * M U = H 632 * M = H Q 633 * 634 * with U and Q unimodular matrices and H a matrix in column echelon form 635 * such that on each echelon row the entries in the non-echelon column 636 * are non-negative (if neg == 0) or non-positive (if neg == 1) 637 * and strictly smaller (in absolute value) than the entries in the echelon 638 * 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 else 692 ++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 else 700 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 start 754 * 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) and 760 * 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 that 799 * 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 first 828 * zero column? That is how many initial columns are non-zero? 829 * Start looking at column "first_col" and only consider 830 * 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 entry 834 * 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 matrix 917 * 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)*right 1007 */ 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 else 1164 ++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 and 1269 * 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 original 1315 * number of variables, then the rows q have already been 1316 * preextended. If the new number is smaller, then the coefficients 1317 * of the divs, which are not changed, need to be shifted down. 1318 * The row q may be the equalities, the inequalities or the 1319 * div expressions. In the latter case, has_div is true and 1320 * 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 else 1332 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', with 1351 * M the matrix mat. 1352 * 1353 * If there are fewer variables x' then there are x, then we perform 1354 * the transformation in place, which means that, in principle, 1355 * this frees up some extra variables as the number 1356 * of columns remains constant, but we would have to extend 1357 * the div array too as the number of rows in this array is assumed 1358 * 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 else 1487 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 else 1496 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" and 1679 * 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 and 2032 * 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 one 2040 * that covers both "mat1" and "mat2". 2041 * The Hermite normal form of the concatenation of the two matrices is 2042 * 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 rows 2048 * in Q1 and Q2. Q1 is a basis for M1, while Q2 extends this basis 2049 * 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 rows 2083 * that is not already present in either "mat1" or "mat2"? 2084 * In other words, is the rank of "mat1" and "mat2" combined equal 2085 * 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