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