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