isl_mat.c

   1 #include "isl_dim.h"
   2 #include "isl_seq.h"
   3 #include "isl_mat.h"
   4 #include "isl_map_private.h"
   5
   6 struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
   7         unsigned n_row, unsigned n_col)
   8 {
   9         int i;
  10         struct isl_mat *mat;
  11
  12         mat = isl_alloc_type(ctx, struct isl_mat);
  13         if (!mat)
  14                 return NULL;
  15
  16         mat->row = NULL;
  17         mat->block = isl_blk_alloc(ctx, n_row * n_col);
  18         if (isl_blk_is_error(mat->block))
  19                 goto error;
  20         mat->row = isl_alloc_array(ctx, isl_int *, n_row);
  21         if (!mat->row)
  22                 goto error;
  23
  24         for (i = 0; i < n_row; ++i)
  25                 mat->row[i] = mat->block.data + i * n_col;
  26
  27         mat->ref = 1;
  28         mat->n_row = n_row;
  29         mat->n_col = n_col;
  30         mat->flags = 0;
  31
  32         return mat;
  33 error:
  34         isl_blk_free(ctx, mat->block);
  35         free(mat);
  36         return NULL;
  37 }
  38
  39 struct isl_mat *isl_mat_extend(struct isl_ctx *ctx, struct isl_mat *mat,
  40         unsigned n_row, unsigned n_col)
  41 {
  42         int i;
  43
  44         if (!mat)
  45                 return NULL;
  46
  47         if (mat->n_col >= n_col && mat->n_row >= n_row)
  48                 return mat;
  49
  50         if (mat->n_col < n_col) {
  51                 struct isl_mat *new_mat;
  52
  53                 new_mat = isl_mat_alloc(ctx, n_row, n_col);
  54                 if (!new_mat)
  55                         goto error;
  56                 for (i = 0; i < mat->n_row; ++i)
  57                         isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
  58                 isl_mat_free(ctx, mat);
  59                 return new_mat;
  60         }
  61
  62         mat = isl_mat_cow(ctx, mat);
  63         if (!mat)
  64                 goto error;
  65
  66         assert(mat->ref == 1);
  67         mat->block = isl_blk_extend(ctx, mat->block, n_row * mat->n_col);
  68         if (isl_blk_is_error(mat->block))
  69                 goto error;
  70         mat->row = isl_realloc_array(ctx, mat->row, isl_int *, n_row);
  71         if (!mat->row)
  72                 goto error;
  73
  74         for (i = 0; i < n_row; ++i)
  75                 mat->row[i] = mat->block.data + i * mat->n_col;
  76         mat->n_row = n_row;
  77
  78         return mat;
  79 error:
  80         isl_mat_free(ctx, mat);
  81         return NULL;
  82 }
  83
  84 struct isl_mat *isl_mat_sub_alloc(struct isl_ctx *ctx, isl_int **row,
  85         unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
  86 {
  87         int i;
  88         struct isl_mat *mat;
  89
  90         mat = isl_alloc_type(ctx, struct isl_mat);
  91         if (!mat)
  92                 return NULL;
  93         mat->row = isl_alloc_array(ctx, isl_int *, n_row);
  94         if (!mat->row)
  95                 goto error;
  96         for (i = 0; i < n_row; ++i)
  97                 mat->row[i] = row[first_row+i] + first_col;
  98         mat->ref = 1;
  99         mat->n_row = n_row;
 100         mat->n_col = n_col;
 101         mat->block = isl_blk_empty();
 102         mat->flags = ISL_MAT_BORROWED;
 103         return mat;
 104 error:
 105         free(mat);
 106         return NULL;
 107 }
 108
 109 void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
 110         unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
 111 {
 112         int i;
 113
 114         for (i = 0; i < n_row; ++i)
 115                 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
 116 }
 117
 118 void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
 119         unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
 120 {
 121         int i;
 122
 123         for (i = 0; i < n_row; ++i)
 124                 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
 125 }
 126
 127 struct isl_mat *isl_mat_copy(struct isl_ctx *ctx, struct isl_mat *mat)
 128 {
 129         if (!mat)
 130                 return NULL;
 131
 132         mat->ref++;
 133         return mat;
 134 }
 135
 136 struct isl_mat *isl_mat_dup(struct isl_ctx *ctx, struct isl_mat *mat)
 137 {
 138         int i;
 139         struct isl_mat *mat2;
 140
 141         if (!mat)
 142                 return NULL;
 143         mat2 = isl_mat_alloc(ctx, mat->n_row, mat->n_col);
 144         if (!mat2)
 145                 return NULL;
 146         for (i = 0; i < mat->n_row; ++i)
 147                 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
 148         return mat2;
 149 }
 150
 151 struct isl_mat *isl_mat_cow(struct isl_ctx *ctx, struct isl_mat *mat)
 152 {
 153         struct isl_mat *mat2;
 154         if (!mat)
 155                 return NULL;
 156
 157         if (mat->ref == 1 && !ISL_F_ISSET(mat, ISL_MAT_BORROWED))
 158                 return mat;
 159
 160         mat2 = isl_mat_dup(ctx, mat);
 161         isl_mat_free(ctx, mat);
 162         return mat2;
 163 }
 164
 165 void isl_mat_free(struct isl_ctx *ctx, struct isl_mat *mat)
 166 {
 167         if (!mat)
 168                 return;
 169
 170         if (--mat->ref > 0)
 171                 return;
 172
 173         if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
 174                 isl_blk_free(ctx, mat->block);
 175         free(mat->row);
 176         free(mat);
 177 }
 178
 179 struct isl_mat *isl_mat_identity(struct isl_ctx *ctx, unsigned n_row)
 180 {
 181         int i;
 182         struct isl_mat *mat;
 183
 184         mat = isl_mat_alloc(ctx, n_row, n_row);
 185         if (!mat)
 186                 return NULL;
 187         for (i = 0; i < n_row; ++i) {
 188                 isl_seq_clr(mat->row[i], i);
 189                 isl_int_set_si(mat->row[i][i], 1);
 190                 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
 191         }
 192
 193         return mat;
 194 }
 195
 196 struct isl_vec *isl_mat_vec_product(struct isl_ctx *ctx,
 197         struct isl_mat *mat, struct isl_vec *vec)
 198 {
 199         int i;
 200         struct isl_vec *prod;
 201
 202         if (!mat || !vec)
 203                 goto error;
 204
 205         isl_assert(ctx, mat->n_col == vec->size, goto error);
 206
 207         prod = isl_vec_alloc(ctx, mat->n_row);
 208         if (!prod)
 209                 goto error;
 210
 211         for (i = 0; i < prod->size; ++i)
 212                 isl_seq_inner_product(mat->row[i], vec->block.data, vec->size,
 213                                         &prod->block.data[i]);
 214         isl_mat_free(ctx, mat);
 215         isl_vec_free(ctx, vec);
 216         return prod;
 217 error:
 218         isl_mat_free(ctx, mat);
 219         isl_vec_free(ctx, vec);
 220         return NULL;
 221 }
 222
 223 struct isl_mat *isl_mat_aff_direct_sum(struct isl_ctx *ctx,
 224         struct isl_mat *left, struct isl_mat *right)
 225 {
 226         int i;
 227         struct isl_mat *sum;
 228
 229         if (!left || !right)
 230                 goto error;
 231
 232         isl_assert(ctx, left->n_row == right->n_row, goto error);
 233         isl_assert(ctx, left->n_row >= 1, goto error);
 234         isl_assert(ctx, left->n_col >= 1, goto error);
 235         isl_assert(ctx, right->n_col >= 1, goto error);
 236         isl_assert(ctx,
 237             isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
 238             goto error);
 239         isl_assert(ctx,
 240             isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
 241             goto error);
 242
 243         sum = isl_mat_alloc(ctx, left->n_row, left->n_col + right->n_col - 1);
 244         if (!sum)
 245                 goto error;
 246         isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
 247         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
 248         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
 249
 250         isl_seq_clr(sum->row[0]+1, sum->n_col-1);
 251         for (i = 1; i < sum->n_row; ++i) {
 252                 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
 253                 isl_int_addmul(sum->row[i][0],
 254                                 right->row[0][0], right->row[i][0]);
 255                 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
 256                                 left->n_col-1);
 257                 isl_seq_scale(sum->row[i]+left->n_col,
 258                                 right->row[i]+1, right->row[0][0],
 259                                 right->n_col-1);
 260         }
 261
 262         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
 263         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
 264         isl_mat_free(ctx, left);
 265         isl_mat_free(ctx, right);
 266         return sum;
 267 error:
 268         isl_mat_free(ctx, left);
 269         isl_mat_free(ctx, right);
 270         return NULL;
 271 }
 272
 273 static void exchange(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
 274         struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
 275 {
 276         int r;
 277         for (r = row; r < M->n_row; ++r)
 278                 isl_int_swap(M->row[r][i], M->row[r][j]);
 279         if (U) {
 280                 for (r = 0; r < (*U)->n_row; ++r)
 281                         isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
 282         }
 283         if (Q)
 284                 isl_mat_swap_rows(ctx, *Q, i, j);
 285 }
 286
 287 static void subtract(struct isl_mat *M, struct isl_mat **U,
 288         struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
 289 {
 290         int r;
 291         for (r = row; r < M->n_row; ++r)
 292                 isl_int_submul(M->row[r][j], m, M->row[r][i]);
 293         if (U) {
 294                 for (r = 0; r < (*U)->n_row; ++r)
 295                         isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
 296         }
 297         if (Q) {
 298                 for (r = 0; r < (*Q)->n_col; ++r)
 299                         isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
 300         }
 301 }
 302
 303 static void oppose(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
 304         struct isl_mat **Q, unsigned row, unsigned col)
 305 {
 306         int r;
 307         for (r = row; r < M->n_row; ++r)
 308                 isl_int_neg(M->row[r][col], M->row[r][col]);
 309         if (U) {
 310                 for (r = 0; r < (*U)->n_row; ++r)
 311                         isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
 312         }
 313         if (Q)
 314                 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
 315 }
 316
 317 /* Given matrix M, compute
 318  *
 319  *              M U = H
 320  *              M   = H Q
 321  *
 322  * with U and Q unimodular matrices and H a matrix in column echelon form
 323  * such that on each echelon row the entries in the non-echelon column
 324  * are non-negative (if neg == 0) or non-positive (if neg == 1)
 325  * and stricly smaller (in absolute value) than the entries in the echelon
 326  * column.
 327  * If U or Q are NULL, then these matrices are not computed.
 328  */
 329 struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
 330         struct isl_mat *M, int neg, struct isl_mat **U, struct isl_mat **Q)
 331 {
 332         isl_int c;
 333         int row, col;
 334
 335         if (U)
 336                 *U = NULL;
 337         if (Q)
 338                 *Q = NULL;
 339         if (!M)
 340                 goto error;
 341         M = isl_mat_cow(ctx, M);
 342         if (!M)
 343                 goto error;
 344         if (U) {
 345                 *U = isl_mat_identity(ctx, M->n_col);
 346                 if (!*U)
 347                         goto error;
 348         }
 349         if (Q) {
 350                 *Q = isl_mat_identity(ctx, M->n_col);
 351                 if (!*Q)
 352                         goto error;
 353         }
 354
 355         col = 0;
 356         isl_int_init(c);
 357         for (row = 0; row < M->n_row; ++row) {
 358                 int first, i, off;
 359                 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
 360                 if (first == -1)
 361                         continue;
 362                 first += col;
 363                 if (first != col)
 364                         exchange(ctx, M, U, Q, row, first, col);
 365                 if (isl_int_is_neg(M->row[row][col]))
 366                         oppose(ctx, M, U, Q, row, col);
 367                 first = col+1;
 368                 while ((off = isl_seq_first_non_zero(M->row[row]+first,
 369                                                        M->n_col-first)) != -1) {
 370                         first += off;
 371                         isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
 372                         subtract(M, U, Q, row, col, first, c);
 373                         if (!isl_int_is_zero(M->row[row][first]))
 374                                 exchange(ctx, M, U, Q, row, first, col);
 375                         else
 376                                 ++first;
 377                 }
 378                 for (i = 0; i < col; ++i) {
 379                         if (isl_int_is_zero(M->row[row][i]))
 380                                 continue;
 381                         if (neg)
 382                                 isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
 383                         else
 384                                 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
 385                         if (isl_int_is_zero(c))
 386                                 continue;
 387                         subtract(M, U, Q, row, col, i, c);
 388                 }
 389                 ++col;
 390         }
 391         isl_int_clear(c);
 392
 393         return M;
 394 error:
 395         if (Q) {
 396                 isl_mat_free(ctx, *Q);
 397                 *Q = NULL;
 398         }
 399         if (U) {
 400                 isl_mat_free(ctx, *U);
 401                 *U = NULL;
 402         }
 403         return NULL;
 404 }
 405
 406 struct isl_mat *isl_mat_right_kernel(struct isl_ctx *ctx, struct isl_mat *mat)
 407 {
 408         int i, rank;
 409         struct isl_mat *U = NULL;
 410         struct isl_mat *K;
 411
 412         mat = isl_mat_left_hermite(ctx, mat, 0, &U, NULL);
 413         if (!mat || !U)
 414                 goto error;
 415
 416         for (i = 0, rank = 0; rank < mat->n_col; ++rank) {
 417                 while (i < mat->n_row && isl_int_is_zero(mat->row[i][rank]))
 418                         ++i;
 419                 if (i >= mat->n_row)
 420                         break;
 421         }
 422         K = isl_mat_alloc(ctx, U->n_row, U->n_col - rank);
 423         if (!K)
 424                 goto error;
 425         isl_mat_sub_copy(ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
 426         isl_mat_free(ctx, mat);
 427         isl_mat_free(ctx, U);
 428         return K;
 429 error:
 430         isl_mat_free(ctx, mat);
 431         isl_mat_free(ctx, U);
 432         return NULL;
 433 }
 434
 435 struct isl_mat *isl_mat_lin_to_aff(struct isl_ctx *ctx, struct isl_mat *mat)
 436 {
 437         int i;
 438         struct isl_mat *mat2;
 439
 440         if (!mat)
 441                 return NULL;
 442         mat2 = isl_mat_alloc(ctx, 1+mat->n_row, 1+mat->n_col);
 443         if (!mat2)
 444                 return NULL;
 445         isl_int_set_si(mat2->row[0][0], 1);
 446         isl_seq_clr(mat2->row[0]+1, mat->n_col);
 447         for (i = 0; i < mat->n_row; ++i) {
 448                 isl_int_set_si(mat2->row[1+i][0], 0);
 449                 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
 450         }
 451         isl_mat_free(ctx, mat);
 452         return mat2;
 453 }
 454
 455 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
 456 {
 457         int i;
 458
 459         for (i = 0; i < n_row; ++i)
 460                 if (!isl_int_is_zero(row[i][col]))
 461                         return i;
 462         return -1;
 463 }
 464
 465 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
 466 {
 467         int i, min = row_first_non_zero(row, n_row, col);
 468         if (min < 0)
 469                 return -1;
 470         for (i = min + 1; i < n_row; ++i) {
 471                 if (isl_int_is_zero(row[i][col]))
 472                         continue;
 473                 if (isl_int_abs_lt(row[i][col], row[min][col]))
 474                         min = i;
 475         }
 476         return min;
 477 }
 478
 479 static void inv_exchange(struct isl_ctx *ctx,
 480         struct isl_mat *left, struct isl_mat *right,
 481         unsigned i, unsigned j)
 482 {
 483         left = isl_mat_swap_rows(ctx, left, i, j);
 484         right = isl_mat_swap_rows(ctx, right, i, j);
 485 }
 486
 487 static void inv_oppose(struct isl_ctx *ctx,
 488         struct isl_mat *left, struct isl_mat *right, unsigned row)
 489 {
 490         isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
 491         isl_seq_neg(right->row[row], right->row[row], right->n_col);
 492 }
 493
 494 static void inv_subtract(struct isl_ctx *ctx,
 495         struct isl_mat *left, struct isl_mat *right,
 496         unsigned row, unsigned i, isl_int m)
 497 {
 498         isl_int_neg(m, m);
 499         isl_seq_combine(left->row[i]+row,
 500                         ctx->one, left->row[i]+row,
 501                         m, left->row[row]+row,
 502                         left->n_col-row);
 503         isl_seq_combine(right->row[i], ctx->one, right->row[i],
 504                         m, right->row[row], right->n_col);
 505 }
 506
 507 /* Compute inv(left)*right
 508  */
 509 struct isl_mat *isl_mat_inverse_product(struct isl_ctx *ctx,
 510         struct isl_mat *left, struct isl_mat *right)
 511 {
 512         int row;
 513         isl_int a, b;
 514
 515         if (!left || !right)
 516                 goto error;
 517
 518         isl_assert(ctx, left->n_row == left->n_col, goto error);
 519         isl_assert(ctx, left->n_row == right->n_row, goto error);
 520
 521         if (left->n_row == 0) {
 522                 isl_mat_free(ctx, left);
 523                 return right;
 524         }
 525
 526         left = isl_mat_cow(ctx, left);
 527         right = isl_mat_cow(ctx, right);
 528         if (!left || !right)
 529                 goto error;
 530
 531         isl_int_init(a);
 532         isl_int_init(b);
 533         for (row = 0; row < left->n_row; ++row) {
 534                 int pivot, first, i, off;
 535                 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
 536                 if (pivot < 0) {
 537                         isl_int_clear(a);
 538                         isl_int_clear(b);
 539                         isl_assert(ctx, pivot >= 0, goto error);
 540                 }
 541                 pivot += row;
 542                 if (pivot != row)
 543                         inv_exchange(ctx, left, right, pivot, row);
 544                 if (isl_int_is_neg(left->row[row][row]))
 545                         inv_oppose(ctx, left, right, row);
 546                 first = row+1;
 547                 while ((off = row_first_non_zero(left->row+first,
 548                                         left->n_row-first, row)) != -1) {
 549                         first += off;
 550                         isl_int_fdiv_q(a, left->row[first][row],
 551                                         left->row[row][row]);
 552                         inv_subtract(ctx, left, right, row, first, a);
 553                         if (!isl_int_is_zero(left->row[first][row]))
 554                                 inv_exchange(ctx, left, right, row, first);
 555                         else
 556                                 ++first;
 557                 }
 558                 for (i = 0; i < row; ++i) {
 559                         if (isl_int_is_zero(left->row[i][row]))
 560                                 continue;
 561                         isl_int_gcd(a, left->row[row][row], left->row[i][row]);
 562                         isl_int_divexact(b, left->row[i][row], a);
 563                         isl_int_divexact(a, left->row[row][row], a);
 564                         isl_int_neg(a, a);
 565                         isl_seq_combine(left->row[i]+row,
 566                                         a, left->row[i]+row,
 567                                         b, left->row[row]+row,
 568                                         left->n_col-row);
 569                         isl_seq_combine(right->row[i], a, right->row[i],
 570                                         b, right->row[row], right->n_col);
 571                 }
 572         }
 573         isl_int_clear(b);
 574
 575         isl_int_set(a, left->row[0][0]);
 576         for (row = 1; row < left->n_row; ++row)
 577                 isl_int_lcm(a, a, left->row[row][row]);
 578         if (isl_int_is_zero(a)){
 579                 isl_int_clear(a);
 580                 isl_assert(ctx, 0, goto error);
 581         }
 582         for (row = 0; row < left->n_row; ++row) {
 583                 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
 584                 if (isl_int_is_one(left->row[row][row]))
 585                         continue;
 586                 isl_seq_scale(right->row[row], right->row[row],
 587                                 left->row[row][row], right->n_col);
 588         }
 589         isl_int_clear(a);
 590
 591         isl_mat_free(ctx, left);
 592         return right;
 593 error:
 594         isl_mat_free(ctx, left);
 595         isl_mat_free(ctx, right);
 596         return NULL;
 597 }
 598
 599 void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
 600 {
 601         int i;
 602
 603         for (i = 0; i < mat->n_row; ++i)
 604                 isl_int_mul(mat->row[i][col], mat->row[i][col], m);
 605 }
 606
 607 void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
 608         isl_int m1, unsigned src1, isl_int m2, unsigned src2)
 609 {
 610         int i;
 611         isl_int tmp;
 612
 613         isl_int_init(tmp);
 614         for (i = 0; i < mat->n_row; ++i) {
 615                 isl_int_mul(tmp, m1, mat->row[i][src1]);
 616                 isl_int_addmul(tmp, m2, mat->row[i][src2]);
 617                 isl_int_set(mat->row[i][dst], tmp);
 618         }
 619         isl_int_clear(tmp);
 620 }
 621
 622 struct isl_mat *isl_mat_right_inverse(struct isl_ctx *ctx,
 623         struct isl_mat *mat)
 624 {
 625         struct isl_mat *inv;
 626         int row;
 627         isl_int a, b;
 628
 629         mat = isl_mat_cow(ctx, mat);
 630         if (!mat)
 631                 return NULL;
 632
 633         inv = isl_mat_identity(ctx, mat->n_col);
 634         inv = isl_mat_cow(ctx, inv);
 635         if (!inv)
 636                 goto error;
 637
 638         isl_int_init(a);
 639         isl_int_init(b);
 640         for (row = 0; row < mat->n_row; ++row) {
 641                 int pivot, first, i, off;
 642                 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
 643                 if (pivot < 0) {
 644                         isl_int_clear(a);
 645                         isl_int_clear(b);
 646                         goto error;
 647                 }
 648                 pivot += row;
 649                 if (pivot != row)
 650                         exchange(ctx, mat, &inv, NULL, row, pivot, row);
 651                 if (isl_int_is_neg(mat->row[row][row]))
 652                         oppose(ctx, mat, &inv, NULL, row, row);
 653                 first = row+1;
 654                 while ((off = isl_seq_first_non_zero(mat->row[row]+first,
 655                                                     mat->n_col-first)) != -1) {
 656                         first += off;
 657                         isl_int_fdiv_q(a, mat->row[row][first],
 658                                                     mat->row[row][row]);
 659                         subtract(mat, &inv, NULL, row, row, first, a);
 660                         if (!isl_int_is_zero(mat->row[row][first]))
 661                                 exchange(ctx, mat, &inv, NULL, row, row, first);
 662                         else
 663                                 ++first;
 664                 }
 665                 for (i = 0; i < row; ++i) {
 666                         if (isl_int_is_zero(mat->row[row][i]))
 667                                 continue;
 668                         isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
 669                         isl_int_divexact(b, mat->row[row][i], a);
 670                         isl_int_divexact(a, mat->row[row][row], a);
 671                         isl_int_neg(a, a);
 672                         isl_mat_col_combine(mat, i, a, i, b, row);
 673                         isl_mat_col_combine(inv, i, a, i, b, row);
 674                 }
 675         }
 676         isl_int_clear(b);
 677
 678         isl_int_set(a, mat->row[0][0]);
 679         for (row = 1; row < mat->n_row; ++row)
 680                 isl_int_lcm(a, a, mat->row[row][row]);
 681         if (isl_int_is_zero(a)){
 682                 isl_int_clear(a);
 683                 goto error;
 684         }
 685         for (row = 0; row < mat->n_row; ++row) {
 686                 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
 687                 if (isl_int_is_one(mat->row[row][row]))
 688                         continue;
 689                 isl_mat_col_scale(inv, row, mat->row[row][row]);
 690         }
 691         isl_int_clear(a);
 692
 693         isl_mat_free(ctx, mat);
 694
 695         return inv;
 696 error:
 697         isl_mat_free(ctx, mat);
 698         return NULL;
 699 }
 700
 701 struct isl_mat *isl_mat_transpose(struct isl_ctx *ctx, struct isl_mat *mat)
 702 {
 703         struct isl_mat *transpose = NULL;
 704         int i, j;
 705
 706         if (mat->n_col == mat->n_row) {
 707                 mat = isl_mat_cow(ctx, mat);
 708                 if (!mat)
 709                         return NULL;
 710                 for (i = 0; i < mat->n_row; ++i)
 711                         for (j = i + 1; j < mat->n_col; ++j)
 712                                 isl_int_swap(mat->row[i][j], mat->row[j][i]);
 713                 return mat;
 714         }
 715         transpose = isl_mat_alloc(ctx, mat->n_col, mat->n_row);
 716         if (!transpose)
 717                 goto error;
 718         for (i = 0; i < mat->n_row; ++i)
 719                 for (j = 0; j < mat->n_col; ++j)
 720                         isl_int_set(transpose->row[j][i], mat->row[i][j]);
 721         isl_mat_free(ctx, mat);
 722         return transpose;
 723 error:
 724         isl_mat_free(ctx, mat);
 725         return NULL;
 726 }
 727
 728 struct isl_mat *isl_mat_swap_cols(struct isl_ctx *ctx,
 729         struct isl_mat *mat, unsigned i, unsigned j)
 730 {
 731         int r;
 732
 733         mat = isl_mat_cow(ctx, mat);
 734         if (!mat)
 735                 return NULL;
 736         isl_assert(ctx, i < mat->n_col, goto error);
 737         isl_assert(ctx, j < mat->n_col, goto error);
 738
 739         for (r = 0; r < mat->n_row; ++r)
 740                 isl_int_swap(mat->row[r][i], mat->row[r][j]);
 741         return mat;
 742 error:
 743         isl_mat_free(ctx, mat);
 744         return NULL;
 745 }
 746
 747 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
 748         struct isl_mat *mat, unsigned i, unsigned j)
 749 {
 750         isl_int *t;
 751
 752         if (!mat)
 753                 return NULL;
 754         mat = isl_mat_cow(ctx, mat);
 755         if (!mat)
 756                 return NULL;
 757         t = mat->row[i];
 758         mat->row[i] = mat->row[j];
 759         mat->row[j] = t;
 760         return mat;
 761 }
 762
 763 struct isl_mat *isl_mat_product(struct isl_ctx *ctx,
 764         struct isl_mat *left, struct isl_mat *right)
 765 {
 766         int i, j, k;
 767         struct isl_mat *prod;
 768
 769         if (!left || !right)
 770                 goto error;
 771         isl_assert(ctx, left->n_col == right->n_row, goto error);
 772         prod = isl_mat_alloc(ctx, left->n_row, right->n_col);
 773         if (!prod)
 774                 goto error;
 775         if (left->n_col == 0) {
 776                 for (i = 0; i < prod->n_row; ++i)
 777                         isl_seq_clr(prod->row[i], prod->n_col);
 778                 return prod;
 779         }
 780         for (i = 0; i < prod->n_row; ++i) {
 781                 for (j = 0; j < prod->n_col; ++j) {
 782                         isl_int_mul(prod->row[i][j],
 783                                     left->row[i][0], right->row[0][j]);
 784                         for (k = 1; k < left->n_col; ++k)
 785                                 isl_int_addmul(prod->row[i][j],
 786                                             left->row[i][k], right->row[k][j]);
 787                 }
 788         }
 789         isl_mat_free(ctx, left);
 790         isl_mat_free(ctx, right);
 791         return prod;
 792 error:
 793         isl_mat_free(ctx, left);
 794         isl_mat_free(ctx, right);
 795         return NULL;
 796 }
 797
 798 /* Replace the variables x in bset by x' given by x = M x', with
 799  * M the matrix mat.
 800  *
 801  * If there are fewer variables x' then there are x, then we perform
 802  * the transformation in place, which that, in principle,
 803  * this frees up some extra variables as the number
 804  * of columns remains constant, but we would have to extend
 805  * the div array too as the number of rows in this array is assumed
 806  * to be equal to extra.
 807  */
 808 struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset,
 809         struct isl_mat *mat)
 810 {
 811         struct isl_ctx *ctx;
 812         struct isl_mat *t;
 813         int i;
 814
 815         if (!bset || !mat)
 816                 goto error;
 817
 818         ctx = bset->ctx;
 819         bset = isl_basic_set_cow(bset);
 820         if (!bset)
 821                 goto error;
 822
 823         isl_assert(ctx, bset->dim->nparam == 0, goto error);
 824         isl_assert(ctx, bset->n_div == 0, goto error);
 825         isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
 826
 827         if (mat->n_col > mat->n_row)
 828                 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0,
 829                                                 0, 0);
 830         else if (mat->n_col < mat->n_row) {
 831                 bset->dim = isl_dim_cow(bset->dim);
 832                 if (!bset->dim)
 833                         goto error;
 834                 bset->dim->n_out -= mat->n_row - mat->n_col;
 835         }
 836
 837         t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
 838         t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
 839         if (!t)
 840                 goto error;
 841         for (i = 0; i < bset->n_eq; ++i) {
 842                 isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
 843                 isl_seq_clr(bset->eq[i]+t->n_col, bset->extra);
 844         }
 845         isl_mat_free(ctx, t);
 846
 847         t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
 848         t = isl_mat_product(ctx, t, mat);
 849         if (!t)
 850                 goto error2;
 851         for (i = 0; i < bset->n_ineq; ++i) {
 852                 isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
 853                 isl_seq_clr(bset->ineq[i]+t->n_col, bset->extra);
 854         }
 855         isl_mat_free(ctx, t);
 856
 857         ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
 858         ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
 859         ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
 860         ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
 861         ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
 862
 863         bset = isl_basic_set_simplify(bset);
 864         bset = isl_basic_set_finalize(bset);
 865
 866         return bset;
 867 error:
 868         isl_mat_free(ctx, mat);
 869 error2:
 870         isl_basic_set_free(bset);
 871         return NULL;
 872 }
 873
 874 struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat)
 875 {
 876         struct isl_ctx *ctx;
 877         int i;
 878
 879         set = isl_set_cow(set);
 880         if (!set)
 881                 return NULL;
 882
 883         ctx = set->ctx;
 884         for (i = 0; i < set->n; ++i) {
 885                 set->p[i] = isl_basic_set_preimage(set->p[i],
 886                                                     isl_mat_copy(ctx, mat));
 887                 if (!set->p[i])
 888                         goto error;
 889         }
 890         if (mat->n_col != mat->n_row) {
 891                 set->dim = isl_dim_cow(set->dim);
 892                 if (!set->dim)
 893                         goto error;
 894                 set->dim->n_out += mat->n_col;
 895                 set->dim->n_out -= mat->n_row;
 896         }
 897         isl_mat_free(ctx, mat);
 898         ISL_F_CLR(set, ISL_SET_NORMALIZED);
 899         return set;
 900 error:
 901         isl_set_free(set);
 902         isl_mat_free(ctx, mat);
 903         return NULL;
 904 }
 905
 906 void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
 907                                 FILE *out, int indent)
 908 {
 909         int i, j;
 910
 911         if (!mat) {
 912                 fprintf(out, "%*snull mat\n", indent, "");
 913                 return;
 914         }
 915
 916         if (mat->n_row == 0)
 917                 fprintf(out, "%*s[]\n", indent, "");
 918
 919         for (i = 0; i < mat->n_row; ++i) {
 920                 if (!i)
 921                         fprintf(out, "%*s[[", indent, "");
 922                 else
 923                         fprintf(out, "%*s[", indent+1, "");
 924                 for (j = 0; j < mat->n_col; ++j) {
 925                         if (j)
 926                             fprintf(out, ",");
 927                         isl_int_print(out, mat->row[i][j], 0);
 928                 }
 929                 if (i == mat->n_row-1)
 930                         fprintf(out, "]]\n");
 931                 else
 932                         fprintf(out, "]\n");
 933         }
 934 }
 935
 936 struct isl_mat *isl_mat_drop_cols(struct isl_ctx *ctx, struct isl_mat *mat,
 937                                 unsigned col, unsigned n)
 938 {
 939         int r;
 940
 941         mat = isl_mat_cow(ctx, mat);
 942         if (!mat)
 943                 return NULL;
 944
 945         if (col != mat->n_col-n) {
 946                 for (r = 0; r < mat->n_row; ++r)
 947                         isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
 948                                         mat->n_col - col - n);
 949         }
 950         mat->n_col -= n;
 951         return mat;
 952 }
 953
 954 struct isl_mat *isl_mat_drop_rows(struct isl_ctx *ctx, struct isl_mat *mat,
 955                                 unsigned row, unsigned n)
 956 {
 957         int r;
 958
 959         mat = isl_mat_cow(ctx, mat);
 960         if (!mat)
 961                 return NULL;
 962
 963         for (r = row; r+n < mat->n_row; ++r)
 964                 mat->row[r] = mat->row[r+n];
 965
 966         mat->n_row -= n;
 967         return mat;
 968 }
 969
 970 void isl_mat_col_submul(struct isl_mat *mat,
 971                         int dst_col, isl_int f, int src_col)
 972 {
 973         int i;
 974
 975         for (i = 0; i < mat->n_row; ++i)
 976                 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
 977 }
 978
 979 void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
 980 {
 981         int i;
 982
 983         for (i = 0; i < mat->n_row; ++i)
 984                 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
 985 }