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