1/* 2 * Copyright 2008-2009 Katholieke Universiteit Leuven 3 * 4 * Use of this software is governed by the MIT license 5 * 6 * Written by Sven Verdoolaege, K.U.Leuven, Departement 7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium 8 */ 9 10#include <isl_ctx_private.h> 11#include <isl_map_private.h> 12#include <isl/space.h> 13#include <isl/seq.h> 14#include <isl_mat_private.h> 15#include <isl_space_private.h> 16#include <isl_val_private.h> 17 18isl_ctx *isl_mat_get_ctx(__isl_keep isl_mat *mat) 19{ 20 return mat ? mat->ctx : NULL; 21} 22 23struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx, 24 unsigned n_row, unsigned n_col) 25{ 26 int i; 27 struct isl_mat *mat; 28 29 mat = isl_alloc_type(ctx, struct isl_mat); 30 if (!mat) 31 return NULL; 32 33 mat->row = NULL; 34 mat->block = isl_blk_alloc(ctx, n_row * n_col); 35 if (isl_blk_is_error(mat->block)) 36 goto error; 37 mat->row = isl_alloc_array(ctx, isl_int *, n_row); 38 if (n_row && !mat->row) 39 goto error; 40 41 for (i = 0; i < n_row; ++i) 42 mat->row[i] = mat->block.data + i * n_col; 43 44 mat->ctx = ctx; 45 isl_ctx_ref(ctx); 46 mat->ref = 1; 47 mat->n_row = n_row; 48 mat->n_col = n_col; 49 mat->max_col = n_col; 50 mat->flags = 0; 51 52 return mat; 53error: 54 isl_blk_free(ctx, mat->block); 55 free(mat); 56 return NULL; 57} 58 59struct isl_mat *isl_mat_extend(struct isl_mat *mat, 60 unsigned n_row, unsigned n_col) 61{ 62 int i; 63 isl_int *old; 64 isl_int **row; 65 66 if (!mat) 67 return NULL; 68 69 if (mat->max_col >= n_col && mat->n_row >= n_row) { 70 if (mat->n_col < n_col) 71 mat->n_col = n_col; 72 return mat; 73 } 74 75 if (mat->max_col < n_col) { 76 struct isl_mat *new_mat; 77 78 if (n_row < mat->n_row) 79 n_row = mat->n_row; 80 new_mat = isl_mat_alloc(mat->ctx, n_row, n_col); 81 if (!new_mat) 82 goto error; 83 for (i = 0; i < mat->n_row; ++i) 84 isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col); 85 isl_mat_free(mat); 86 return new_mat; 87 } 88 89 mat = isl_mat_cow(mat); 90 if (!mat) 91 goto error; 92 93 old = mat->block.data; 94 mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col); 95 if (isl_blk_is_error(mat->block)) 96 goto error; 97 row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row); 98 if (n_row && !row) 99 goto error; 100 mat->row = row; 101 102 for (i = 0; i < mat->n_row; ++i) 103 mat->row[i] = mat->block.data + (mat->row[i] - old); 104 for (i = mat->n_row; i < n_row; ++i) 105 mat->row[i] = mat->block.data + i * mat->max_col; 106 mat->n_row = n_row; 107 if (mat->n_col < n_col) 108 mat->n_col = n_col; 109 110 return mat; 111error: 112 isl_mat_free(mat); 113 return NULL; 114} 115 116__isl_give isl_mat *isl_mat_sub_alloc6(isl_ctx *ctx, isl_int **row, 117 unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col) 118{ 119 int i; 120 struct isl_mat *mat; 121 122 mat = isl_alloc_type(ctx, struct isl_mat); 123 if (!mat) 124 return NULL; 125 mat->row = isl_alloc_array(ctx, isl_int *, n_row); 126 if (n_row && !mat->row) 127 goto error; 128 for (i = 0; i < n_row; ++i) 129 mat->row[i] = row[first_row+i] + first_col; 130 mat->ctx = ctx; 131 isl_ctx_ref(ctx); 132 mat->ref = 1; 133 mat->n_row = n_row; 134 mat->n_col = n_col; 135 mat->block = isl_blk_empty(); 136 mat->flags = ISL_MAT_BORROWED; 137 return mat; 138error: 139 free(mat); 140 return NULL; 141} 142 143__isl_give isl_mat *isl_mat_sub_alloc(__isl_keep isl_mat *mat, 144 unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col) 145{ 146 if (!mat) 147 return NULL; 148 return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row, 149 first_col, n_col); 150} 151 152void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src, 153 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col) 154{ 155 int i; 156 157 for (i = 0; i < n_row; ++i) 158 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col); 159} 160 161void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src, 162 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col) 163{ 164 int i; 165 166 for (i = 0; i < n_row; ++i) 167 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col); 168} 169 170struct isl_mat *isl_mat_copy(struct isl_mat *mat) 171{ 172 if (!mat) 173 return NULL; 174 175 mat->ref++; 176 return mat; 177} 178 179struct isl_mat *isl_mat_dup(struct isl_mat *mat) 180{ 181 int i; 182 struct isl_mat *mat2; 183 184 if (!mat) 185 return NULL; 186 mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col); 187 if (!mat2) 188 return NULL; 189 for (i = 0; i < mat->n_row; ++i) 190 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col); 191 return mat2; 192} 193 194struct isl_mat *isl_mat_cow(struct isl_mat *mat) 195{ 196 struct isl_mat *mat2; 197 if (!mat) 198 return NULL; 199 200 if (mat->ref == 1 && !ISL_F_ISSET(mat, ISL_MAT_BORROWED)) 201 return mat; 202 203 mat2 = isl_mat_dup(mat); 204 isl_mat_free(mat); 205 return mat2; 206} 207 208void *isl_mat_free(struct isl_mat *mat) 209{ 210 if (!mat) 211 return NULL; 212 213 if (--mat->ref > 0) 214 return NULL; 215 216 if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED)) 217 isl_blk_free(mat->ctx, mat->block); 218 isl_ctx_deref(mat->ctx); 219 free(mat->row); 220 free(mat); 221 222 return NULL; 223} 224 225int isl_mat_rows(__isl_keep isl_mat *mat) 226{ 227 return mat ? mat->n_row : -1; 228} 229 230int isl_mat_cols(__isl_keep isl_mat *mat) 231{ 232 return mat ? mat->n_col : -1; 233} 234 235int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v) 236{ 237 if (!mat) 238 return -1; 239 if (row < 0 || row >= mat->n_row) 240 isl_die(mat->ctx, isl_error_invalid, "row out of range", 241 return -1); 242 if (col < 0 || col >= mat->n_col) 243 isl_die(mat->ctx, isl_error_invalid, "column out of range", 244 return -1); 245 isl_int_set(*v, mat->row[row][col]); 246 return 0; 247} 248 249/* Extract the element at row "row", oolumn "col" of "mat". 250 */ 251__isl_give isl_val *isl_mat_get_element_val(__isl_keep isl_mat *mat, 252 int row, int col) 253{ 254 isl_ctx *ctx; 255 256 if (!mat) 257 return NULL; 258 ctx = isl_mat_get_ctx(mat); 259 if (row < 0 || row >= mat->n_row) 260 isl_die(ctx, isl_error_invalid, "row out of range", 261 return NULL); 262 if (col < 0 || col >= mat->n_col) 263 isl_die(ctx, isl_error_invalid, "column out of range", 264 return NULL); 265 return isl_val_int_from_isl_int(ctx, mat->row[row][col]); 266} 267 268__isl_give isl_mat *isl_mat_set_element(__isl_take isl_mat *mat, 269 int row, int col, isl_int v) 270{ 271 mat = isl_mat_cow(mat); 272 if (!mat) 273 return NULL; 274 if (row < 0 || row >= mat->n_row) 275 isl_die(mat->ctx, isl_error_invalid, "row out of range", 276 goto error); 277 if (col < 0 || col >= mat->n_col) 278 isl_die(mat->ctx, isl_error_invalid, "column out of range", 279 goto error); 280 isl_int_set(mat->row[row][col], v); 281 return mat; 282error: 283 isl_mat_free(mat); 284 return NULL; 285} 286 287__isl_give isl_mat *isl_mat_set_element_si(__isl_take isl_mat *mat, 288 int row, int col, int v) 289{ 290 mat = isl_mat_cow(mat); 291 if (!mat) 292 return NULL; 293 if (row < 0 || row >= mat->n_row) 294 isl_die(mat->ctx, isl_error_invalid, "row out of range", 295 goto error); 296 if (col < 0 || col >= mat->n_col) 297 isl_die(mat->ctx, isl_error_invalid, "column out of range", 298 goto error); 299 isl_int_set_si(mat->row[row][col], v); 300 return mat; 301error: 302 isl_mat_free(mat); 303 return NULL; 304} 305 306/* Replace the element at row "row", column "col" of "mat" by "v". 307 */ 308__isl_give isl_mat *isl_mat_set_element_val(__isl_take isl_mat *mat, 309 int row, int col, __isl_take isl_val *v) 310{ 311 if (!v) 312 return isl_mat_free(mat); 313 if (!isl_val_is_int(v)) 314 isl_die(isl_val_get_ctx(v), isl_error_invalid, 315 "expecting integer value", goto error); 316 mat = isl_mat_set_element(mat, row, col, v->n); 317 isl_val_free(v); 318 return mat; 319error: 320 isl_val_free(v); 321 return isl_mat_free(mat); 322} 323 324__isl_give isl_mat *isl_mat_diag(isl_ctx *ctx, unsigned n_row, isl_int d) 325{ 326 int i; 327 struct isl_mat *mat; 328 329 mat = isl_mat_alloc(ctx, n_row, n_row); 330 if (!mat) 331 return NULL; 332 for (i = 0; i < n_row; ++i) { 333 isl_seq_clr(mat->row[i], i); 334 isl_int_set(mat->row[i][i], d); 335 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1)); 336 } 337 338 return mat; 339} 340 341__isl_give isl_mat *isl_mat_identity(isl_ctx *ctx, unsigned n_row) 342{ 343 if (!ctx) 344 return NULL; 345 return isl_mat_diag(ctx, n_row, ctx->one); 346} 347 348struct isl_vec *isl_mat_vec_product(struct isl_mat *mat, struct isl_vec *vec) 349{ 350 int i; 351 struct isl_vec *prod; 352 353 if (!mat || !vec) 354 goto error; 355 356 isl_assert(mat->ctx, mat->n_col == vec->size, goto error); 357 358 prod = isl_vec_alloc(mat->ctx, mat->n_row); 359 if (!prod) 360 goto error; 361 362 for (i = 0; i < prod->size; ++i) 363 isl_seq_inner_product(mat->row[i], vec->el, vec->size, 364 &prod->block.data[i]); 365 isl_mat_free(mat); 366 isl_vec_free(vec); 367 return prod; 368error: 369 isl_mat_free(mat); 370 isl_vec_free(vec); 371 return NULL; 372} 373 374__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat, 375 __isl_take isl_vec *vec) 376{ 377 struct isl_mat *vec_mat; 378 int i; 379 380 if (!mat || !vec) 381 goto error; 382 vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1); 383 if (!vec_mat) 384 goto error; 385 for (i = 0; i < vec->size; ++i) 386 isl_int_set(vec_mat->row[i][0], vec->el[i]); 387 vec_mat = isl_mat_inverse_product(mat, vec_mat); 388 isl_vec_free(vec); 389 if (!vec_mat) 390 return NULL; 391 vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row); 392 if (vec) 393 for (i = 0; i < vec->size; ++i) 394 isl_int_set(vec->el[i], vec_mat->row[i][0]); 395 isl_mat_free(vec_mat); 396 return vec; 397error: 398 isl_mat_free(mat); 399 isl_vec_free(vec); 400 return NULL; 401} 402 403struct isl_vec *isl_vec_mat_product(struct isl_vec *vec, struct isl_mat *mat) 404{ 405 int i, j; 406 struct isl_vec *prod; 407 408 if (!mat || !vec) 409 goto error; 410 411 isl_assert(mat->ctx, mat->n_row == vec->size, goto error); 412 413 prod = isl_vec_alloc(mat->ctx, mat->n_col); 414 if (!prod) 415 goto error; 416 417 for (i = 0; i < prod->size; ++i) { 418 isl_int_set_si(prod->el[i], 0); 419 for (j = 0; j < vec->size; ++j) 420 isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]); 421 } 422 isl_mat_free(mat); 423 isl_vec_free(vec); 424 return prod; 425error: 426 isl_mat_free(mat); 427 isl_vec_free(vec); 428 return NULL; 429} 430 431struct isl_mat *isl_mat_aff_direct_sum(struct isl_mat *left, 432 struct isl_mat *right) 433{ 434 int i; 435 struct isl_mat *sum; 436 437 if (!left || !right) 438 goto error; 439 440 isl_assert(left->ctx, left->n_row == right->n_row, goto error); 441 isl_assert(left->ctx, left->n_row >= 1, goto error); 442 isl_assert(left->ctx, left->n_col >= 1, goto error); 443 isl_assert(left->ctx, right->n_col >= 1, goto error); 444 isl_assert(left->ctx, 445 isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1, 446 goto error); 447 isl_assert(left->ctx, 448 isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1, 449 goto error); 450 451 sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1); 452 if (!sum) 453 goto error; 454 isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]); 455 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]); 456 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]); 457 458 isl_seq_clr(sum->row[0]+1, sum->n_col-1); 459 for (i = 1; i < sum->n_row; ++i) { 460 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]); 461 isl_int_addmul(sum->row[i][0], 462 right->row[0][0], right->row[i][0]); 463 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0], 464 left->n_col-1); 465 isl_seq_scale(sum->row[i]+left->n_col, 466 right->row[i]+1, right->row[0][0], 467 right->n_col-1); 468 } 469 470 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]); 471 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]); 472 isl_mat_free(left); 473 isl_mat_free(right); 474 return sum; 475error: 476 isl_mat_free(left); 477 isl_mat_free(right); 478 return NULL; 479} 480 481static void exchange(struct isl_mat *M, struct isl_mat **U, 482 struct isl_mat **Q, unsigned row, unsigned i, unsigned j) 483{ 484 int r; 485 for (r = row; r < M->n_row; ++r) 486 isl_int_swap(M->row[r][i], M->row[r][j]); 487 if (U) { 488 for (r = 0; r < (*U)->n_row; ++r) 489 isl_int_swap((*U)->row[r][i], (*U)->row[r][j]); 490 } 491 if (Q) 492 isl_mat_swap_rows(*Q, i, j); 493} 494 495static void subtract(struct isl_mat *M, struct isl_mat **U, 496 struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m) 497{ 498 int r; 499 for (r = row; r < M->n_row; ++r) 500 isl_int_submul(M->row[r][j], m, M->row[r][i]); 501 if (U) { 502 for (r = 0; r < (*U)->n_row; ++r) 503 isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]); 504 } 505 if (Q) { 506 for (r = 0; r < (*Q)->n_col; ++r) 507 isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]); 508 } 509} 510 511static void oppose(struct isl_mat *M, struct isl_mat **U, 512 struct isl_mat **Q, unsigned row, unsigned col) 513{ 514 int r; 515 for (r = row; r < M->n_row; ++r) 516 isl_int_neg(M->row[r][col], M->row[r][col]); 517 if (U) { 518 for (r = 0; r < (*U)->n_row; ++r) 519 isl_int_neg((*U)->row[r][col], (*U)->row[r][col]); 520 } 521 if (Q) 522 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col); 523} 524 525/* Given matrix M, compute 526 * 527 * M U = H 528 * M = H Q 529 * 530 * with U and Q unimodular matrices and H a matrix in column echelon form 531 * such that on each echelon row the entries in the non-echelon column 532 * are non-negative (if neg == 0) or non-positive (if neg == 1) 533 * and stricly smaller (in absolute value) than the entries in the echelon 534 * column. 535 * If U or Q are NULL, then these matrices are not computed. 536 */ 537struct isl_mat *isl_mat_left_hermite(struct isl_mat *M, int neg, 538 struct isl_mat **U, struct isl_mat **Q) 539{ 540 isl_int c; 541 int row, col; 542 543 if (U) 544 *U = NULL; 545 if (Q) 546 *Q = NULL; 547 if (!M) 548 goto error; 549 M = isl_mat_cow(M); 550 if (!M) 551 goto error; 552 if (U) { 553 *U = isl_mat_identity(M->ctx, M->n_col); 554 if (!*U) 555 goto error; 556 } 557 if (Q) { 558 *Q = isl_mat_identity(M->ctx, M->n_col); 559 if (!*Q) 560 goto error; 561 } 562 563 col = 0; 564 isl_int_init(c); 565 for (row = 0; row < M->n_row; ++row) { 566 int first, i, off; 567 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col); 568 if (first == -1) 569 continue; 570 first += col; 571 if (first != col) 572 exchange(M, U, Q, row, first, col); 573 if (isl_int_is_neg(M->row[row][col])) 574 oppose(M, U, Q, row, col); 575 first = col+1; 576 while ((off = isl_seq_first_non_zero(M->row[row]+first, 577 M->n_col-first)) != -1) { 578 first += off; 579 isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]); 580 subtract(M, U, Q, row, col, first, c); 581 if (!isl_int_is_zero(M->row[row][first])) 582 exchange(M, U, Q, row, first, col); 583 else 584 ++first; 585 } 586 for (i = 0; i < col; ++i) { 587 if (isl_int_is_zero(M->row[row][i])) 588 continue; 589 if (neg) 590 isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]); 591 else 592 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]); 593 if (isl_int_is_zero(c)) 594 continue; 595 subtract(M, U, Q, row, col, i, c); 596 } 597 ++col; 598 } 599 isl_int_clear(c); 600 601 return M; 602error: 603 if (Q) { 604 isl_mat_free(*Q); 605 *Q = NULL; 606 } 607 if (U) { 608 isl_mat_free(*U); 609 *U = NULL; 610 } 611 isl_mat_free(M); 612 return NULL; 613} 614 615struct isl_mat *isl_mat_right_kernel(struct isl_mat *mat) 616{ 617 int i, rank; 618 struct isl_mat *U = NULL; 619 struct isl_mat *K; 620 621 mat = isl_mat_left_hermite(mat, 0, &U, NULL); 622 if (!mat || !U) 623 goto error; 624 625 for (i = 0, rank = 0; rank < mat->n_col; ++rank) { 626 while (i < mat->n_row && isl_int_is_zero(mat->row[i][rank])) 627 ++i; 628 if (i >= mat->n_row) 629 break; 630 } 631 K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank); 632 if (!K) 633 goto error; 634 isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank); 635 isl_mat_free(mat); 636 isl_mat_free(U); 637 return K; 638error: 639 isl_mat_free(mat); 640 isl_mat_free(U); 641 return NULL; 642} 643 644struct isl_mat *isl_mat_lin_to_aff(struct isl_mat *mat) 645{ 646 int i; 647 struct isl_mat *mat2; 648 649 if (!mat) 650 return NULL; 651 mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col); 652 if (!mat2) 653 goto error; 654 isl_int_set_si(mat2->row[0][0], 1); 655 isl_seq_clr(mat2->row[0]+1, mat->n_col); 656 for (i = 0; i < mat->n_row; ++i) { 657 isl_int_set_si(mat2->row[1+i][0], 0); 658 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col); 659 } 660 isl_mat_free(mat); 661 return mat2; 662error: 663 isl_mat_free(mat); 664 return NULL; 665} 666 667/* Given two matrices M1 and M2, return the block matrix 668 * 669 * [ M1 0 ] 670 * [ 0 M2 ] 671 */ 672__isl_give isl_mat *isl_mat_diagonal(__isl_take isl_mat *mat1, 673 __isl_take isl_mat *mat2) 674{ 675 int i; 676 isl_mat *mat; 677 678 if (!mat1 || !mat2) 679 goto error; 680 681 mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row, 682 mat1->n_col + mat2->n_col); 683 if (!mat) 684 goto error; 685 for (i = 0; i < mat1->n_row; ++i) { 686 isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col); 687 isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col); 688 } 689 for (i = 0; i < mat2->n_row; ++i) { 690 isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col); 691 isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col, 692 mat2->row[i], mat2->n_col); 693 } 694 isl_mat_free(mat1); 695 isl_mat_free(mat2); 696 return mat; 697error: 698 isl_mat_free(mat1); 699 isl_mat_free(mat2); 700 return NULL; 701} 702 703static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col) 704{ 705 int i; 706 707 for (i = 0; i < n_row; ++i) 708 if (!isl_int_is_zero(row[i][col])) 709 return i; 710 return -1; 711} 712 713static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col) 714{ 715 int i, min = row_first_non_zero(row, n_row, col); 716 if (min < 0) 717 return -1; 718 for (i = min + 1; i < n_row; ++i) { 719 if (isl_int_is_zero(row[i][col])) 720 continue; 721 if (isl_int_abs_lt(row[i][col], row[min][col])) 722 min = i; 723 } 724 return min; 725} 726 727static void inv_exchange(struct isl_mat *left, struct isl_mat *right, 728 unsigned i, unsigned j) 729{ 730 left = isl_mat_swap_rows(left, i, j); 731 right = isl_mat_swap_rows(right, i, j); 732} 733 734static void inv_oppose( 735 struct isl_mat *left, struct isl_mat *right, unsigned row) 736{ 737 isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row); 738 isl_seq_neg(right->row[row], right->row[row], right->n_col); 739} 740 741static void inv_subtract(struct isl_mat *left, struct isl_mat *right, 742 unsigned row, unsigned i, isl_int m) 743{ 744 isl_int_neg(m, m); 745 isl_seq_combine(left->row[i]+row, 746 left->ctx->one, left->row[i]+row, 747 m, left->row[row]+row, 748 left->n_col-row); 749 isl_seq_combine(right->row[i], right->ctx->one, right->row[i], 750 m, right->row[row], right->n_col); 751} 752 753/* Compute inv(left)*right 754 */ 755struct isl_mat *isl_mat_inverse_product(struct isl_mat *left, 756 struct isl_mat *right) 757{ 758 int row; 759 isl_int a, b; 760 761 if (!left || !right) 762 goto error; 763 764 isl_assert(left->ctx, left->n_row == left->n_col, goto error); 765 isl_assert(left->ctx, left->n_row == right->n_row, goto error); 766 767 if (left->n_row == 0) { 768 isl_mat_free(left); 769 return right; 770 } 771 772 left = isl_mat_cow(left); 773 right = isl_mat_cow(right); 774 if (!left || !right) 775 goto error; 776 777 isl_int_init(a); 778 isl_int_init(b); 779 for (row = 0; row < left->n_row; ++row) { 780 int pivot, first, i, off; 781 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row); 782 if (pivot < 0) { 783 isl_int_clear(a); 784 isl_int_clear(b); 785 isl_assert(left->ctx, pivot >= 0, goto error); 786 } 787 pivot += row; 788 if (pivot != row) 789 inv_exchange(left, right, pivot, row); 790 if (isl_int_is_neg(left->row[row][row])) 791 inv_oppose(left, right, row); 792 first = row+1; 793 while ((off = row_first_non_zero(left->row+first, 794 left->n_row-first, row)) != -1) { 795 first += off; 796 isl_int_fdiv_q(a, left->row[first][row], 797 left->row[row][row]); 798 inv_subtract(left, right, row, first, a); 799 if (!isl_int_is_zero(left->row[first][row])) 800 inv_exchange(left, right, row, first); 801 else 802 ++first; 803 } 804 for (i = 0; i < row; ++i) { 805 if (isl_int_is_zero(left->row[i][row])) 806 continue; 807 isl_int_gcd(a, left->row[row][row], left->row[i][row]); 808 isl_int_divexact(b, left->row[i][row], a); 809 isl_int_divexact(a, left->row[row][row], a); 810 isl_int_neg(b, b); 811 isl_seq_combine(left->row[i] + i, 812 a, left->row[i] + i, 813 b, left->row[row] + i, 814 left->n_col - i); 815 isl_seq_combine(right->row[i], a, right->row[i], 816 b, right->row[row], right->n_col); 817 } 818 } 819 isl_int_clear(b); 820 821 isl_int_set(a, left->row[0][0]); 822 for (row = 1; row < left->n_row; ++row) 823 isl_int_lcm(a, a, left->row[row][row]); 824 if (isl_int_is_zero(a)){ 825 isl_int_clear(a); 826 isl_assert(left->ctx, 0, goto error); 827 } 828 for (row = 0; row < left->n_row; ++row) { 829 isl_int_divexact(left->row[row][row], a, left->row[row][row]); 830 if (isl_int_is_one(left->row[row][row])) 831 continue; 832 isl_seq_scale(right->row[row], right->row[row], 833 left->row[row][row], right->n_col); 834 } 835 isl_int_clear(a); 836 837 isl_mat_free(left); 838 return right; 839error: 840 isl_mat_free(left); 841 isl_mat_free(right); 842 return NULL; 843} 844 845void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m) 846{ 847 int i; 848 849 for (i = 0; i < mat->n_row; ++i) 850 isl_int_mul(mat->row[i][col], mat->row[i][col], m); 851} 852 853void isl_mat_col_combine(struct isl_mat *mat, unsigned dst, 854 isl_int m1, unsigned src1, isl_int m2, unsigned src2) 855{ 856 int i; 857 isl_int tmp; 858 859 isl_int_init(tmp); 860 for (i = 0; i < mat->n_row; ++i) { 861 isl_int_mul(tmp, m1, mat->row[i][src1]); 862 isl_int_addmul(tmp, m2, mat->row[i][src2]); 863 isl_int_set(mat->row[i][dst], tmp); 864 } 865 isl_int_clear(tmp); 866} 867 868struct isl_mat *isl_mat_right_inverse(struct isl_mat *mat) 869{ 870 struct isl_mat *inv; 871 int row; 872 isl_int a, b; 873 874 mat = isl_mat_cow(mat); 875 if (!mat) 876 return NULL; 877 878 inv = isl_mat_identity(mat->ctx, mat->n_col); 879 inv = isl_mat_cow(inv); 880 if (!inv) 881 goto error; 882 883 isl_int_init(a); 884 isl_int_init(b); 885 for (row = 0; row < mat->n_row; ++row) { 886 int pivot, first, i, off; 887 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row); 888 if (pivot < 0) { 889 isl_int_clear(a); 890 isl_int_clear(b); 891 isl_assert(mat->ctx, pivot >= 0, goto error); 892 } 893 pivot += row; 894 if (pivot != row) 895 exchange(mat, &inv, NULL, row, pivot, row); 896 if (isl_int_is_neg(mat->row[row][row])) 897 oppose(mat, &inv, NULL, row, row); 898 first = row+1; 899 while ((off = isl_seq_first_non_zero(mat->row[row]+first, 900 mat->n_col-first)) != -1) { 901 first += off; 902 isl_int_fdiv_q(a, mat->row[row][first], 903 mat->row[row][row]); 904 subtract(mat, &inv, NULL, row, row, first, a); 905 if (!isl_int_is_zero(mat->row[row][first])) 906 exchange(mat, &inv, NULL, row, row, first); 907 else 908 ++first; 909 } 910 for (i = 0; i < row; ++i) { 911 if (isl_int_is_zero(mat->row[row][i])) 912 continue; 913 isl_int_gcd(a, mat->row[row][row], mat->row[row][i]); 914 isl_int_divexact(b, mat->row[row][i], a); 915 isl_int_divexact(a, mat->row[row][row], a); 916 isl_int_neg(a, a); 917 isl_mat_col_combine(mat, i, a, i, b, row); 918 isl_mat_col_combine(inv, i, a, i, b, row); 919 } 920 } 921 isl_int_clear(b); 922 923 isl_int_set(a, mat->row[0][0]); 924 for (row = 1; row < mat->n_row; ++row) 925 isl_int_lcm(a, a, mat->row[row][row]); 926 if (isl_int_is_zero(a)){ 927 isl_int_clear(a); 928 goto error; 929 } 930 for (row = 0; row < mat->n_row; ++row) { 931 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]); 932 if (isl_int_is_one(mat->row[row][row])) 933 continue; 934 isl_mat_col_scale(inv, row, mat->row[row][row]); 935 } 936 isl_int_clear(a); 937 938 isl_mat_free(mat); 939 940 return inv; 941error: 942 isl_mat_free(mat); 943 isl_mat_free(inv); 944 return NULL; 945} 946 947struct isl_mat *isl_mat_transpose(struct isl_mat *mat) 948{ 949 struct isl_mat *transpose = NULL; 950 int i, j; 951 952 if (!mat) 953 return NULL; 954 955 if (mat->n_col == mat->n_row) { 956 mat = isl_mat_cow(mat); 957 if (!mat) 958 return NULL; 959 for (i = 0; i < mat->n_row; ++i) 960 for (j = i + 1; j < mat->n_col; ++j) 961 isl_int_swap(mat->row[i][j], mat->row[j][i]); 962 return mat; 963 } 964 transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row); 965 if (!transpose) 966 goto error; 967 for (i = 0; i < mat->n_row; ++i) 968 for (j = 0; j < mat->n_col; ++j) 969 isl_int_set(transpose->row[j][i], mat->row[i][j]); 970 isl_mat_free(mat); 971 return transpose; 972error: 973 isl_mat_free(mat); 974 return NULL; 975} 976 977struct isl_mat *isl_mat_swap_cols(struct isl_mat *mat, unsigned i, unsigned j) 978{ 979 int r; 980 981 mat = isl_mat_cow(mat); 982 if (!mat) 983 return NULL; 984 isl_assert(mat->ctx, i < mat->n_col, goto error); 985 isl_assert(mat->ctx, j < mat->n_col, goto error); 986 987 for (r = 0; r < mat->n_row; ++r) 988 isl_int_swap(mat->row[r][i], mat->row[r][j]); 989 return mat; 990error: 991 isl_mat_free(mat); 992 return NULL; 993} 994 995struct isl_mat *isl_mat_swap_rows(struct isl_mat *mat, unsigned i, unsigned j) 996{ 997 isl_int *t; 998 999 if (!mat) 1000 return NULL; 1001 mat = isl_mat_cow(mat); 1002 if (!mat) 1003 return NULL; 1004 t = mat->row[i]; 1005 mat->row[i] = mat->row[j]; 1006 mat->row[j] = t; 1007 return mat; 1008} 1009 1010__isl_give isl_mat *isl_mat_product(__isl_take isl_mat *left, 1011 __isl_take isl_mat *right) 1012{ 1013 int i, j, k; 1014 struct isl_mat *prod; 1015 1016 if (!left || !right) 1017 goto error; 1018 isl_assert(left->ctx, left->n_col == right->n_row, goto error); 1019 prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col); 1020 if (!prod) 1021 goto error; 1022 if (left->n_col == 0) { 1023 for (i = 0; i < prod->n_row; ++i) 1024 isl_seq_clr(prod->row[i], prod->n_col); 1025 isl_mat_free(left); 1026 isl_mat_free(right); 1027 return prod; 1028 } 1029 for (i = 0; i < prod->n_row; ++i) { 1030 for (j = 0; j < prod->n_col; ++j) { 1031 isl_int_mul(prod->row[i][j], 1032 left->row[i][0], right->row[0][j]); 1033 for (k = 1; k < left->n_col; ++k) 1034 isl_int_addmul(prod->row[i][j], 1035 left->row[i][k], right->row[k][j]); 1036 } 1037 } 1038 isl_mat_free(left); 1039 isl_mat_free(right); 1040 return prod; 1041error: 1042 isl_mat_free(left); 1043 isl_mat_free(right); 1044 return NULL; 1045} 1046 1047/* Replace the variables x in the rows q by x' given by x = M x', 1048 * with M the matrix mat. 1049 * 1050 * If the number of new variables is greater than the original 1051 * number of variables, then the rows q have already been 1052 * preextended. If the new number is smaller, then the coefficients 1053 * of the divs, which are not changed, need to be shifted down. 1054 * The row q may be the equalities, the inequalities or the 1055 * div expressions. In the latter case, has_div is true and 1056 * we need to take into account the extra denominator column. 1057 */ 1058static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n, 1059 unsigned n_div, int has_div, struct isl_mat *mat) 1060{ 1061 int i; 1062 struct isl_mat *t; 1063 int e; 1064 1065 if (mat->n_col >= mat->n_row) 1066 e = 0; 1067 else 1068 e = mat->n_row - mat->n_col; 1069 if (has_div) 1070 for (i = 0; i < n; ++i) 1071 isl_int_mul(q[i][0], q[i][0], mat->row[0][0]); 1072 t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row); 1073 t = isl_mat_product(t, mat); 1074 if (!t) 1075 return -1; 1076 for (i = 0; i < n; ++i) { 1077 isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col); 1078 isl_seq_cpy(q[i] + has_div + t->n_col, 1079 q[i] + has_div + t->n_col + e, n_div); 1080 isl_seq_clr(q[i] + has_div + t->n_col + n_div, e); 1081 } 1082 isl_mat_free(t); 1083 return 0; 1084} 1085 1086/* Replace the variables x in bset by x' given by x = M x', with 1087 * M the matrix mat. 1088 * 1089 * If there are fewer variables x' then there are x, then we perform 1090 * the transformation in place, which that, in principle, 1091 * this frees up some extra variables as the number 1092 * of columns remains constant, but we would have to extend 1093 * the div array too as the number of rows in this array is assumed 1094 * to be equal to extra. 1095 */ 1096struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset, 1097 struct isl_mat *mat) 1098{ 1099 struct isl_ctx *ctx; 1100 1101 if (!bset || !mat) 1102 goto error; 1103 1104 ctx = bset->ctx; 1105 bset = isl_basic_set_cow(bset); 1106 if (!bset) 1107 goto error; 1108 1109 isl_assert(ctx, bset->dim->nparam == 0, goto error); 1110 isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error); 1111 isl_assert(ctx, mat->n_col > 0, goto error); 1112 1113 if (mat->n_col > mat->n_row) { 1114 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0); 1115 if (!bset) 1116 goto error; 1117 } else if (mat->n_col < mat->n_row) { 1118 bset->dim = isl_space_cow(bset->dim); 1119 if (!bset->dim) 1120 goto error; 1121 bset->dim->n_out -= mat->n_row - mat->n_col; 1122 } 1123 1124 if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0, 1125 isl_mat_copy(mat)) < 0) 1126 goto error; 1127 1128 if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0, 1129 isl_mat_copy(mat)) < 0) 1130 goto error; 1131 1132 if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0) 1133 goto error2; 1134 1135 ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT); 1136 ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT); 1137 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED); 1138 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS); 1139 ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES); 1140 1141 bset = isl_basic_set_simplify(bset); 1142 bset = isl_basic_set_finalize(bset); 1143 1144 return bset; 1145error: 1146 isl_mat_free(mat); 1147error2: 1148 isl_basic_set_free(bset); 1149 return NULL; 1150} 1151 1152struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat) 1153{ 1154 struct isl_ctx *ctx; 1155 int i; 1156 1157 set = isl_set_cow(set); 1158 if (!set) 1159 return NULL; 1160 1161 ctx = set->ctx; 1162 for (i = 0; i < set->n; ++i) { 1163 set->p[i] = isl_basic_set_preimage(set->p[i], 1164 isl_mat_copy(mat)); 1165 if (!set->p[i]) 1166 goto error; 1167 } 1168 if (mat->n_col != mat->n_row) { 1169 set->dim = isl_space_cow(set->dim); 1170 if (!set->dim) 1171 goto error; 1172 set->dim->n_out += mat->n_col; 1173 set->dim->n_out -= mat->n_row; 1174 } 1175 isl_mat_free(mat); 1176 ISL_F_CLR(set, ISL_SET_NORMALIZED); 1177 return set; 1178error: 1179 isl_set_free(set); 1180 isl_mat_free(mat); 1181 return NULL; 1182} 1183 1184/* Replace the variables x starting at pos in the rows q 1185 * by x' with x = M x' with M the matrix mat. 1186 * That is, replace the corresponding coefficients c by c M. 1187 */ 1188static int transform(isl_ctx *ctx, isl_int **q, unsigned n, 1189 unsigned pos, __isl_take isl_mat *mat) 1190{ 1191 int i; 1192 isl_mat *t; 1193 1194 t = isl_mat_sub_alloc6(ctx, q, 0, n, pos, mat->n_row); 1195 t = isl_mat_product(t, mat); 1196 if (!t) 1197 return -1; 1198 for (i = 0; i < n; ++i) 1199 isl_seq_swp_or_cpy(q[i] + pos, t->row[i], t->n_col); 1200 isl_mat_free(t); 1201 return 0; 1202} 1203 1204/* Replace the variables x of type "type" starting at "first" in "bset" 1205 * by x' with x = M x' with M the matrix trans. 1206 * That is, replace the corresponding coefficients c by c M. 1207 * 1208 * The transformation matrix should be a square matrix. 1209 */ 1210__isl_give isl_basic_set *isl_basic_set_transform_dims( 1211 __isl_take isl_basic_set *bset, enum isl_dim_type type, unsigned first, 1212 __isl_take isl_mat *trans) 1213{ 1214 isl_ctx *ctx; 1215 unsigned pos; 1216 1217 bset = isl_basic_set_cow(bset); 1218 if (!bset || !trans) 1219 goto error; 1220 1221 ctx = isl_basic_set_get_ctx(bset); 1222 if (trans->n_row != trans->n_col) 1223 isl_die(trans->ctx, isl_error_invalid, 1224 "expecting square transformation matrix", goto error); 1225 if (first + trans->n_row > isl_basic_set_dim(bset, type)) 1226 isl_die(trans->ctx, isl_error_invalid, 1227 "oversized transformation matrix", goto error); 1228 1229 pos = isl_basic_set_offset(bset, type) + first; 1230 1231 if (transform(ctx, bset->eq, bset->n_eq, pos, isl_mat_copy(trans)) < 0) 1232 goto error; 1233 if (transform(ctx, bset->ineq, bset->n_ineq, pos, 1234 isl_mat_copy(trans)) < 0) 1235 goto error; 1236 if (transform(ctx, bset->div, bset->n_div, 1 + pos, 1237 isl_mat_copy(trans)) < 0) 1238 goto error; 1239 1240 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED); 1241 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS); 1242 1243 isl_mat_free(trans); 1244 return bset; 1245error: 1246 isl_mat_free(trans); 1247 isl_basic_set_free(bset); 1248 return NULL; 1249} 1250 1251void isl_mat_print_internal(__isl_keep isl_mat *mat, FILE *out, int indent) 1252{ 1253 int i, j; 1254 1255 if (!mat) { 1256 fprintf(out, "%*snull mat\n", indent, ""); 1257 return; 1258 } 1259 1260 if (mat->n_row == 0) 1261 fprintf(out, "%*s[]\n", indent, ""); 1262 1263 for (i = 0; i < mat->n_row; ++i) { 1264 if (!i) 1265 fprintf(out, "%*s[[", indent, ""); 1266 else 1267 fprintf(out, "%*s[", indent+1, ""); 1268 for (j = 0; j < mat->n_col; ++j) { 1269 if (j) 1270 fprintf(out, ","); 1271 isl_int_print(out, mat->row[i][j], 0); 1272 } 1273 if (i == mat->n_row-1) 1274 fprintf(out, "]]\n"); 1275 else 1276 fprintf(out, "]\n"); 1277 } 1278} 1279 1280void isl_mat_dump(__isl_keep isl_mat *mat) 1281{ 1282 isl_mat_print_internal(mat, stderr, 0); 1283} 1284 1285struct isl_mat *isl_mat_drop_cols(struct isl_mat *mat, unsigned col, unsigned n) 1286{ 1287 int r; 1288 1289 if (n == 0) 1290 return mat; 1291 1292 mat = isl_mat_cow(mat); 1293 if (!mat) 1294 return NULL; 1295 1296 if (col != mat->n_col-n) { 1297 for (r = 0; r < mat->n_row; ++r) 1298 isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n, 1299 mat->n_col - col - n); 1300 } 1301 mat->n_col -= n; 1302 return mat; 1303} 1304 1305struct isl_mat *isl_mat_drop_rows(struct isl_mat *mat, unsigned row, unsigned n) 1306{ 1307 int r; 1308 1309 mat = isl_mat_cow(mat); 1310 if (!mat) 1311 return NULL; 1312 1313 for (r = row; r+n < mat->n_row; ++r) 1314 mat->row[r] = mat->row[r+n]; 1315 1316 mat->n_row -= n; 1317 return mat; 1318} 1319 1320__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat, 1321 unsigned col, unsigned n) 1322{ 1323 isl_mat *ext; 1324 1325 if (!mat) 1326 return NULL; 1327 if (n == 0) 1328 return mat; 1329 1330 ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n); 1331 if (!ext) 1332 goto error; 1333 1334 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col); 1335 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 1336 col + n, col, mat->n_col - col); 1337 1338 isl_mat_free(mat); 1339 return ext; 1340error: 1341 isl_mat_free(mat); 1342 return NULL; 1343} 1344 1345__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat, 1346 unsigned first, unsigned n) 1347{ 1348 int i; 1349 1350 if (!mat) 1351 return NULL; 1352 mat = isl_mat_insert_cols(mat, first, n); 1353 if (!mat) 1354 return NULL; 1355 1356 for (i = 0; i < mat->n_row; ++i) 1357 isl_seq_clr(mat->row[i] + first, n); 1358 1359 return mat; 1360} 1361 1362__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n) 1363{ 1364 if (!mat) 1365 return NULL; 1366 1367 return isl_mat_insert_zero_cols(mat, mat->n_col, n); 1368} 1369 1370__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat, 1371 unsigned row, unsigned n) 1372{ 1373 isl_mat *ext; 1374 1375 if (!mat) 1376 return NULL; 1377 if (n == 0) 1378 return mat; 1379 1380 ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col); 1381 if (!ext) 1382 goto error; 1383 1384 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col); 1385 isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row, 1386 mat->n_row - row, 0, 0, mat->n_col); 1387 1388 isl_mat_free(mat); 1389 return ext; 1390error: 1391 isl_mat_free(mat); 1392 return NULL; 1393} 1394 1395__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n) 1396{ 1397 if (!mat) 1398 return NULL; 1399 1400 return isl_mat_insert_rows(mat, mat->n_row, n); 1401} 1402 1403__isl_give isl_mat *isl_mat_insert_zero_rows(__isl_take isl_mat *mat, 1404 unsigned row, unsigned n) 1405{ 1406 int i; 1407 1408 mat = isl_mat_insert_rows(mat, row, n); 1409 if (!mat) 1410 return NULL; 1411 1412 for (i = 0; i < n; ++i) 1413 isl_seq_clr(mat->row[row + i], mat->n_col); 1414 1415 return mat; 1416} 1417 1418__isl_give isl_mat *isl_mat_add_zero_rows(__isl_take isl_mat *mat, unsigned n) 1419{ 1420 if (!mat) 1421 return NULL; 1422 1423 return isl_mat_insert_zero_rows(mat, mat->n_row, n); 1424} 1425 1426void isl_mat_col_submul(struct isl_mat *mat, 1427 int dst_col, isl_int f, int src_col) 1428{ 1429 int i; 1430 1431 for (i = 0; i < mat->n_row; ++i) 1432 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]); 1433} 1434 1435void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col) 1436{ 1437 int i; 1438 1439 if (!mat) 1440 return; 1441 1442 for (i = 0; i < mat->n_row; ++i) 1443 isl_int_add(mat->row[i][dst_col], 1444 mat->row[i][dst_col], mat->row[i][src_col]); 1445} 1446 1447void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col) 1448{ 1449 int i; 1450 1451 for (i = 0; i < mat->n_row; ++i) 1452 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]); 1453} 1454 1455struct isl_mat *isl_mat_unimodular_complete(struct isl_mat *M, int row) 1456{ 1457 int r; 1458 struct isl_mat *H = NULL, *Q = NULL; 1459 1460 if (!M) 1461 return NULL; 1462 1463 isl_assert(M->ctx, M->n_row == M->n_col, goto error); 1464 M->n_row = row; 1465 H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q); 1466 M->n_row = M->n_col; 1467 if (!H) 1468 goto error; 1469 for (r = 0; r < row; ++r) 1470 isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error); 1471 for (r = row; r < M->n_row; ++r) 1472 isl_seq_cpy(M->row[r], Q->row[r], M->n_col); 1473 isl_mat_free(H); 1474 isl_mat_free(Q); 1475 return M; 1476error: 1477 isl_mat_free(H); 1478 isl_mat_free(Q); 1479 isl_mat_free(M); 1480 return NULL; 1481} 1482 1483__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top, 1484 __isl_take isl_mat *bot) 1485{ 1486 struct isl_mat *mat; 1487 1488 if (!top || !bot) 1489 goto error; 1490 1491 isl_assert(top->ctx, top->n_col == bot->n_col, goto error); 1492 if (top->n_row == 0) { 1493 isl_mat_free(top); 1494 return bot; 1495 } 1496 if (bot->n_row == 0) { 1497 isl_mat_free(bot); 1498 return top; 1499 } 1500 1501 mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col); 1502 if (!mat) 1503 goto error; 1504 isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row, 1505 0, 0, mat->n_col); 1506 isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row, 1507 0, 0, mat->n_col); 1508 isl_mat_free(top); 1509 isl_mat_free(bot); 1510 return mat; 1511error: 1512 isl_mat_free(top); 1513 isl_mat_free(bot); 1514 return NULL; 1515} 1516 1517int isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2) 1518{ 1519 int i; 1520 1521 if (!mat1 || !mat2) 1522 return -1; 1523 1524 if (mat1->n_row != mat2->n_row) 1525 return 0; 1526 1527 if (mat1->n_col != mat2->n_col) 1528 return 0; 1529 1530 for (i = 0; i < mat1->n_row; ++i) 1531 if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col)) 1532 return 0; 1533 1534 return 1; 1535} 1536 1537__isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec) 1538{ 1539 struct isl_mat *mat; 1540 1541 if (!vec) 1542 return NULL; 1543 mat = isl_mat_alloc(vec->ctx, 1, vec->size); 1544 if (!mat) 1545 goto error; 1546 1547 isl_seq_cpy(mat->row[0], vec->el, vec->size); 1548 1549 isl_vec_free(vec); 1550 return mat; 1551error: 1552 isl_vec_free(vec); 1553 return NULL; 1554} 1555 1556__isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top, 1557 __isl_take isl_vec *bot) 1558{ 1559 return isl_mat_concat(top, isl_mat_from_row_vec(bot)); 1560} 1561 1562__isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat, 1563 unsigned dst_col, unsigned src_col, unsigned n) 1564{ 1565 isl_mat *res; 1566 1567 if (!mat) 1568 return NULL; 1569 if (n == 0 || dst_col == src_col) 1570 return mat; 1571 1572 res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col); 1573 if (!res) 1574 goto error; 1575 1576 if (dst_col < src_col) { 1577 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row, 1578 0, 0, dst_col); 1579 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row, 1580 dst_col, src_col, n); 1581 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row, 1582 dst_col + n, dst_col, src_col - dst_col); 1583 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row, 1584 src_col + n, src_col + n, 1585 res->n_col - src_col - n); 1586 } else { 1587 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row, 1588 0, 0, src_col); 1589 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row, 1590 src_col, src_col + n, dst_col - src_col); 1591 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row, 1592 dst_col, src_col, n); 1593 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row, 1594 dst_col + n, dst_col + n, 1595 res->n_col - dst_col - n); 1596 } 1597 isl_mat_free(mat); 1598 1599 return res; 1600error: 1601 isl_mat_free(mat); 1602 return NULL; 1603} 1604 1605void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd) 1606{ 1607 int i; 1608 isl_int g; 1609 1610 isl_int_set_si(*gcd, 0); 1611 if (!mat) 1612 return; 1613 1614 isl_int_init(g); 1615 for (i = 0; i < mat->n_row; ++i) { 1616 isl_seq_gcd(mat->row[i], mat->n_col, &g); 1617 isl_int_gcd(*gcd, *gcd, g); 1618 } 1619 isl_int_clear(g); 1620} 1621 1622__isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m) 1623{ 1624 int i; 1625 1626 if (isl_int_is_one(m)) 1627 return mat; 1628 1629 mat = isl_mat_cow(mat); 1630 if (!mat) 1631 return NULL; 1632 1633 for (i = 0; i < mat->n_row; ++i) 1634 isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col); 1635 1636 return mat; 1637} 1638 1639__isl_give isl_mat *isl_mat_scale_down_row(__isl_take isl_mat *mat, int row, 1640 isl_int m) 1641{ 1642 if (isl_int_is_one(m)) 1643 return mat; 1644 1645 mat = isl_mat_cow(mat); 1646 if (!mat) 1647 return NULL; 1648 1649 isl_seq_scale_down(mat->row[row], mat->row[row], m, mat->n_col); 1650 1651 return mat; 1652} 1653 1654__isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat) 1655{ 1656 isl_int gcd; 1657 1658 if (!mat) 1659 return NULL; 1660 1661 isl_int_init(gcd); 1662 isl_mat_gcd(mat, &gcd); 1663 mat = isl_mat_scale_down(mat, gcd); 1664 isl_int_clear(gcd); 1665 1666 return mat; 1667} 1668 1669__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row) 1670{ 1671 mat = isl_mat_cow(mat); 1672 if (!mat) 1673 return NULL; 1674 1675 isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col); 1676 1677 return mat; 1678} 1679 1680/* Number of initial non-zero columns. 1681 */ 1682int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat) 1683{ 1684 int i; 1685 1686 if (!mat) 1687 return -1; 1688 1689 for (i = 0; i < mat->n_col; ++i) 1690 if (row_first_non_zero(mat->row, mat->n_row, i) < 0) 1691 break; 1692 1693 return i; 1694} 1695