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