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