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