isl_mat.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   7  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
   8  */
   9
  10 #include <isl/dim.h>
  11 #include <isl/seq.h>
  12 #include <isl_mat_private.h>
  13 #include "isl_map_private.h"
  14 #include <isl_dim_private.h>
  15
  16 struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
  17         unsigned n_row, unsigned n_col)
  18 {
  19         int i;
  20         struct isl_mat *mat;
  21
  22         mat = isl_alloc_type(ctx, struct isl_mat);
  23         if (!mat)
  24                 return NULL;
  25
  26         mat->row = NULL;
  27         mat->block = isl_blk_alloc(ctx, n_row * n_col);
  28         if (isl_blk_is_error(mat->block))
  29                 goto error;
  30         mat->row = isl_alloc_array(ctx, isl_int *, n_row);
  31         if (!mat->row)
  32                 goto error;
  33
  34         for (i = 0; i < n_row; ++i)
  35                 mat->row[i] = mat->block.data + i * n_col;
  36
  37         mat->ctx = ctx;
  38         isl_ctx_ref(ctx);
  39         mat->ref = 1;
  40         mat->n_row = n_row;
  41         mat->n_col = n_col;
  42         mat->max_col = n_col;
  43         mat->flags = 0;
  44
  45         return mat;
  46 error:
  47         isl_blk_free(ctx, mat->block);
  48         free(mat);
  49         return NULL;
  50 }
  51
  52 struct isl_mat *isl_mat_extend(struct isl_mat *mat,
  53         unsigned n_row, unsigned n_col)
  54 {
  55         int i;
  56         isl_int *old;
  57
  58         if (!mat)
  59                 return NULL;
  60
  61         if (mat->max_col >= n_col && mat->n_row >= n_row) {
  62                 if (mat->n_col < n_col)
  63                         mat->n_col = n_col;
  64                 return mat;
  65         }
  66
  67         if (mat->max_col < n_col) {
  68                 struct isl_mat *new_mat;
  69
  70                 if (n_row < mat->n_row)
  71                         n_row = mat->n_row;
  72                 new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
  73                 if (!new_mat)
  74                         goto error;
  75                 for (i = 0; i < mat->n_row; ++i)
  76                         isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
  77                 isl_mat_free(mat);
  78                 return new_mat;
  79         }
  80
  81         mat = isl_mat_cow(mat);
  82         if (!mat)
  83                 goto error;
  84
  85         old = mat->block.data;
  86         mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
  87         if (isl_blk_is_error(mat->block))
  88                 goto error;
  89         mat->row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);
  90         if (!mat->row)
  91                 goto error;
  92
  93         for (i = 0; i < mat->n_row; ++i)
  94                 mat->row[i] = mat->block.data + (mat->row[i] - old);
  95         for (i = mat->n_row; i < n_row; ++i)
  96                 mat->row[i] = mat->block.data + i * mat->max_col;
  97         mat->n_row = n_row;
  98         if (mat->n_col < n_col)
  99                 mat->n_col = n_col;
 100
 101         return mat;
 102 error:
 103         isl_mat_free(mat);
 104         return NULL;
 105 }
 106
 107 struct isl_mat *isl_mat_sub_alloc(struct isl_ctx *ctx, isl_int **row,
 108         unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
 109 {
 110         int i;
 111         struct isl_mat *mat;
 112
 113         mat = isl_alloc_type(ctx, struct isl_mat);
 114         if (!mat)
 115                 return NULL;
 116         mat->row = isl_alloc_array(ctx, isl_int *, n_row);
 117         if (!mat->row)
 118                 goto error;
 119         for (i = 0; i < n_row; ++i)
 120                 mat->row[i] = row[first_row+i] + first_col;
 121         mat->ctx = ctx;
 122         isl_ctx_ref(ctx);
 123         mat->ref = 1;
 124         mat->n_row = n_row;
 125         mat->n_col = n_col;
 126         mat->block = isl_blk_empty();
 127         mat->flags = ISL_MAT_BORROWED;
 128         return mat;
 129 error:
 130         free(mat);
 131         return NULL;
 132 }
 133
 134 void isl_mat_sub_copy(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_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
 141 }
 142
 143 void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
 144         unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
 145 {
 146         int i;
 147
 148         for (i = 0; i < n_row; ++i)
 149                 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
 150 }
 151
 152 struct isl_mat *isl_mat_copy(struct isl_mat *mat)
 153 {
 154         if (!mat)
 155                 return NULL;
 156
 157         mat->ref++;
 158         return mat;
 159 }
 160
 161 struct isl_mat *isl_mat_dup(struct isl_mat *mat)
 162 {
 163         int i;
 164         struct isl_mat *mat2;
 165
 166         if (!mat)
 167                 return NULL;
 168         mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
 169         if (!mat2)
 170                 return NULL;
 171         for (i = 0; i < mat->n_row; ++i)
 172                 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
 173         return mat2;
 174 }
 175
 176 struct isl_mat *isl_mat_cow(struct isl_mat *mat)
 177 {
 178         struct isl_mat *mat2;
 179         if (!mat)
 180                 return NULL;
 181
 182         if (mat->ref == 1 && !ISL_F_ISSET(mat, ISL_MAT_BORROWED))
 183                 return mat;
 184
 185         mat2 = isl_mat_dup(mat);
 186         isl_mat_free(mat);
 187         return mat2;
 188 }
 189
 190 void isl_mat_free(struct isl_mat *mat)
 191 {
 192         if (!mat)
 193                 return;
 194
 195         if (--mat->ref > 0)
 196                 return;
 197
 198         if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
 199                 isl_blk_free(mat->ctx, mat->block);
 200         isl_ctx_deref(mat->ctx);
 201         free(mat->row);
 202         free(mat);
 203 }
 204
 205 int isl_mat_rows(__isl_keep isl_mat *mat)
 206 {
 207         return mat ? mat->n_row : -1;
 208 }
 209
 210 int isl_mat_cols(__isl_keep isl_mat *mat)
 211 {
 212         return mat ? mat->n_col : -1;
 213 }
 214
 215 int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
 216 {
 217         if (!mat)
 218                 return -1;
 219         if (row < 0 || row >= mat->n_row)
 220                 isl_die(mat->ctx, isl_error_invalid, "row out of range",
 221                         return -1);
 222         if (col < 0 || col >= mat->n_col)
 223                 isl_die(mat->ctx, isl_error_invalid, "column out of range",
 224                         return -1);
 225         isl_int_set(*v, mat->row[row][col]);
 226         return 0;
 227 }
 228
 229 __isl_give isl_mat *isl_mat_set_element(__isl_take isl_mat *mat,
 230         int row, int col, isl_int v)
 231 {
 232         mat = isl_mat_cow(mat);
 233         if (!mat)
 234                 return NULL;
 235         if (row < 0 || row >= mat->n_row)
 236                 isl_die(mat->ctx, isl_error_invalid, "row out of range",
 237                         goto error);
 238         if (col < 0 || col >= mat->n_col)
 239                 isl_die(mat->ctx, isl_error_invalid, "column out of range",
 240                         goto error);
 241         isl_int_set(mat->row[row][col], v);
 242         return mat;
 243 error:
 244         isl_mat_free(mat);
 245         return NULL;
 246 }
 247
 248 struct isl_mat *isl_mat_identity(struct isl_ctx *ctx, unsigned n_row)
 249 {
 250         int i;
 251         struct isl_mat *mat;
 252
 253         mat = isl_mat_alloc(ctx, n_row, n_row);
 254         if (!mat)
 255                 return NULL;
 256         for (i = 0; i < n_row; ++i) {
 257                 isl_seq_clr(mat->row[i], i);
 258                 isl_int_set_si(mat->row[i][i], 1);
 259                 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
 260         }
 261
 262         return mat;
 263 }
 264
 265 struct isl_vec *isl_mat_vec_product(struct isl_mat *mat, struct isl_vec *vec)
 266 {
 267         int i;
 268         struct isl_vec *prod;
 269
 270         if (!mat || !vec)
 271                 goto error;
 272
 273         isl_assert(mat->ctx, mat->n_col == vec->size, goto error);
 274
 275         prod = isl_vec_alloc(mat->ctx, mat->n_row);
 276         if (!prod)
 277                 goto error;
 278
 279         for (i = 0; i < prod->size; ++i)
 280                 isl_seq_inner_product(mat->row[i], vec->el, vec->size,
 281                                         &prod->block.data[i]);
 282         isl_mat_free(mat);
 283         isl_vec_free(vec);
 284         return prod;
 285 error:
 286         isl_mat_free(mat);
 287         isl_vec_free(vec);
 288         return NULL;
 289 }
 290
 291 __isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
 292         __isl_take isl_vec *vec)
 293 {
 294         struct isl_mat *vec_mat;
 295         int i;
 296
 297         if (!mat || !vec)
 298                 goto error;
 299         vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
 300         if (!vec_mat)
 301                 goto error;
 302         for (i = 0; i < vec->size; ++i)
 303                 isl_int_set(vec_mat->row[i][0], vec->el[i]);
 304         vec_mat = isl_mat_inverse_product(mat, vec_mat);
 305         isl_vec_free(vec);
 306         if (!vec_mat)
 307                 return NULL;
 308         vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
 309         if (vec)
 310                 for (i = 0; i < vec->size; ++i)
 311                         isl_int_set(vec->el[i], vec_mat->row[i][0]);
 312         isl_mat_free(vec_mat);
 313         return vec;
 314 error:
 315         isl_mat_free(mat);
 316         isl_vec_free(vec);
 317         return NULL;
 318 }
 319
 320 struct isl_vec *isl_vec_mat_product(struct isl_vec *vec, struct isl_mat *mat)
 321 {
 322         int i, j;
 323         struct isl_vec *prod;
 324
 325         if (!mat || !vec)
 326                 goto error;
 327
 328         isl_assert(mat->ctx, mat->n_row == vec->size, goto error);
 329
 330         prod = isl_vec_alloc(mat->ctx, mat->n_col);
 331         if (!prod)
 332                 goto error;
 333
 334         for (i = 0; i < prod->size; ++i) {
 335                 isl_int_set_si(prod->el[i], 0);
 336                 for (j = 0; j < vec->size; ++j)
 337                         isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
 338         }
 339         isl_mat_free(mat);
 340         isl_vec_free(vec);
 341         return prod;
 342 error:
 343         isl_mat_free(mat);
 344         isl_vec_free(vec);
 345         return NULL;
 346 }
 347
 348 struct isl_mat *isl_mat_aff_direct_sum(struct isl_mat *left,
 349         struct isl_mat *right)
 350 {
 351         int i;
 352         struct isl_mat *sum;
 353
 354         if (!left || !right)
 355                 goto error;
 356
 357         isl_assert(left->ctx, left->n_row == right->n_row, goto error);
 358         isl_assert(left->ctx, left->n_row >= 1, goto error);
 359         isl_assert(left->ctx, left->n_col >= 1, goto error);
 360         isl_assert(left->ctx, right->n_col >= 1, goto error);
 361         isl_assert(left->ctx,
 362             isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
 363             goto error);
 364         isl_assert(left->ctx,
 365             isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
 366             goto error);
 367
 368         sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
 369         if (!sum)
 370                 goto error;
 371         isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
 372         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
 373         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
 374
 375         isl_seq_clr(sum->row[0]+1, sum->n_col-1);
 376         for (i = 1; i < sum->n_row; ++i) {
 377                 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
 378                 isl_int_addmul(sum->row[i][0],
 379                                 right->row[0][0], right->row[i][0]);
 380                 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
 381                                 left->n_col-1);
 382                 isl_seq_scale(sum->row[i]+left->n_col,
 383                                 right->row[i]+1, right->row[0][0],
 384                                 right->n_col-1);
 385         }
 386
 387         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
 388         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
 389         isl_mat_free(left);
 390         isl_mat_free(right);
 391         return sum;
 392 error:
 393         isl_mat_free(left);
 394         isl_mat_free(right);
 395         return NULL;
 396 }
 397
 398 static void exchange(struct isl_mat *M, struct isl_mat **U,
 399         struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
 400 {
 401         int r;
 402         for (r = row; r < M->n_row; ++r)
 403                 isl_int_swap(M->row[r][i], M->row[r][j]);
 404         if (U) {
 405                 for (r = 0; r < (*U)->n_row; ++r)
 406                         isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
 407         }
 408         if (Q)
 409                 isl_mat_swap_rows(*Q, i, j);
 410 }
 411
 412 static void subtract(struct isl_mat *M, struct isl_mat **U,
 413         struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
 414 {
 415         int r;
 416         for (r = row; r < M->n_row; ++r)
 417                 isl_int_submul(M->row[r][j], m, M->row[r][i]);
 418         if (U) {
 419                 for (r = 0; r < (*U)->n_row; ++r)
 420                         isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
 421         }
 422         if (Q) {
 423                 for (r = 0; r < (*Q)->n_col; ++r)
 424                         isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
 425         }
 426 }
 427
 428 static void oppose(struct isl_mat *M, struct isl_mat **U,
 429         struct isl_mat **Q, unsigned row, unsigned col)
 430 {
 431         int r;
 432         for (r = row; r < M->n_row; ++r)
 433                 isl_int_neg(M->row[r][col], M->row[r][col]);
 434         if (U) {
 435                 for (r = 0; r < (*U)->n_row; ++r)
 436                         isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
 437         }
 438         if (Q)
 439                 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
 440 }
 441
 442 /* Given matrix M, compute
 443  *
 444  *              M U = H
 445  *              M   = H Q
 446  *
 447  * with U and Q unimodular matrices and H a matrix in column echelon form
 448  * such that on each echelon row the entries in the non-echelon column
 449  * are non-negative (if neg == 0) or non-positive (if neg == 1)
 450  * and stricly smaller (in absolute value) than the entries in the echelon
 451  * column.
 452  * If U or Q are NULL, then these matrices are not computed.
 453  */
 454 struct isl_mat *isl_mat_left_hermite(struct isl_mat *M, int neg,
 455         struct isl_mat **U, struct isl_mat **Q)
 456 {
 457         isl_int c;
 458         int row, col;
 459
 460         if (U)
 461                 *U = NULL;
 462         if (Q)
 463                 *Q = NULL;
 464         if (!M)
 465                 goto error;
 466         M = isl_mat_cow(M);
 467         if (!M)
 468                 goto error;
 469         if (U) {
 470                 *U = isl_mat_identity(M->ctx, M->n_col);
 471                 if (!*U)
 472                         goto error;
 473         }
 474         if (Q) {
 475                 *Q = isl_mat_identity(M->ctx, M->n_col);
 476                 if (!*Q)
 477                         goto error;
 478         }
 479
 480         col = 0;
 481         isl_int_init(c);
 482         for (row = 0; row < M->n_row; ++row) {
 483                 int first, i, off;
 484                 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
 485                 if (first == -1)
 486                         continue;
 487                 first += col;
 488                 if (first != col)
 489                         exchange(M, U, Q, row, first, col);
 490                 if (isl_int_is_neg(M->row[row][col]))
 491                         oppose(M, U, Q, row, col);
 492                 first = col+1;
 493                 while ((off = isl_seq_first_non_zero(M->row[row]+first,
 494                                                        M->n_col-first)) != -1) {
 495                         first += off;
 496                         isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
 497                         subtract(M, U, Q, row, col, first, c);
 498                         if (!isl_int_is_zero(M->row[row][first]))
 499                                 exchange(M, U, Q, row, first, col);
 500                         else
 501                                 ++first;
 502                 }
 503                 for (i = 0; i < col; ++i) {
 504                         if (isl_int_is_zero(M->row[row][i]))
 505                                 continue;
 506                         if (neg)
 507                                 isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
 508                         else
 509                                 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
 510                         if (isl_int_is_zero(c))
 511                                 continue;
 512                         subtract(M, U, Q, row, col, i, c);
 513                 }
 514                 ++col;
 515         }
 516         isl_int_clear(c);
 517
 518         return M;
 519 error:
 520         if (Q) {
 521                 isl_mat_free(*Q);
 522                 *Q = NULL;
 523         }
 524         if (U) {
 525                 isl_mat_free(*U);
 526                 *U = NULL;
 527         }
 528         return NULL;
 529 }
 530
 531 struct isl_mat *isl_mat_right_kernel(struct isl_mat *mat)
 532 {
 533         int i, rank;
 534         struct isl_mat *U = NULL;
 535         struct isl_mat *K;
 536
 537         mat = isl_mat_left_hermite(mat, 0, &U, NULL);
 538         if (!mat || !U)
 539                 goto error;
 540
 541         for (i = 0, rank = 0; rank < mat->n_col; ++rank) {
 542                 while (i < mat->n_row && isl_int_is_zero(mat->row[i][rank]))
 543                         ++i;
 544                 if (i >= mat->n_row)
 545                         break;
 546         }
 547         K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank);
 548         if (!K)
 549                 goto error;
 550         isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
 551         isl_mat_free(mat);
 552         isl_mat_free(U);
 553         return K;
 554 error:
 555         isl_mat_free(mat);
 556         isl_mat_free(U);
 557         return NULL;
 558 }
 559
 560 struct isl_mat *isl_mat_lin_to_aff(struct isl_mat *mat)
 561 {
 562         int i;
 563         struct isl_mat *mat2;
 564
 565         if (!mat)
 566                 return NULL;
 567         mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
 568         if (!mat2)
 569                 goto error;
 570         isl_int_set_si(mat2->row[0][0], 1);
 571         isl_seq_clr(mat2->row[0]+1, mat->n_col);
 572         for (i = 0; i < mat->n_row; ++i) {
 573                 isl_int_set_si(mat2->row[1+i][0], 0);
 574                 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
 575         }
 576         isl_mat_free(mat);
 577         return mat2;
 578 error:
 579         isl_mat_free(mat);
 580         return NULL;
 581 }
 582
 583 /* Given two matrices M1 and M2, return the block matrix
 584  *
 585  *      [ M1  0  ]
 586  *      [ 0   M2 ]
 587  */
 588 __isl_give isl_mat *isl_mat_diagonal(__isl_take isl_mat *mat1,
 589         __isl_take isl_mat *mat2)
 590 {
 591         int i;
 592         isl_mat *mat;
 593
 594         if (!mat1 || !mat2)
 595                 goto error;
 596
 597         mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
 598                                        mat1->n_col + mat2->n_col);
 599         if (!mat)
 600                 goto error;
 601         for (i = 0; i < mat1->n_row; ++i) {
 602                 isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
 603                 isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
 604         }
 605         for (i = 0; i < mat2->n_row; ++i) {
 606                 isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
 607                 isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
 608                                                     mat2->row[i], mat2->n_col);
 609         }
 610         isl_mat_free(mat1);
 611         isl_mat_free(mat2);
 612         return mat;
 613 error:
 614         isl_mat_free(mat1);
 615         isl_mat_free(mat2);
 616         return NULL;
 617 }
 618
 619 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
 620 {
 621         int i;
 622
 623         for (i = 0; i < n_row; ++i)
 624                 if (!isl_int_is_zero(row[i][col]))
 625                         return i;
 626         return -1;
 627 }
 628
 629 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
 630 {
 631         int i, min = row_first_non_zero(row, n_row, col);
 632         if (min < 0)
 633                 return -1;
 634         for (i = min + 1; i < n_row; ++i) {
 635                 if (isl_int_is_zero(row[i][col]))
 636                         continue;
 637                 if (isl_int_abs_lt(row[i][col], row[min][col]))
 638                         min = i;
 639         }
 640         return min;
 641 }
 642
 643 static void inv_exchange(struct isl_mat *left, struct isl_mat *right,
 644         unsigned i, unsigned j)
 645 {
 646         left = isl_mat_swap_rows(left, i, j);
 647         right = isl_mat_swap_rows(right, i, j);
 648 }
 649
 650 static void inv_oppose(
 651         struct isl_mat *left, struct isl_mat *right, unsigned row)
 652 {
 653         isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
 654         isl_seq_neg(right->row[row], right->row[row], right->n_col);
 655 }
 656
 657 static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
 658         unsigned row, unsigned i, isl_int m)
 659 {
 660         isl_int_neg(m, m);
 661         isl_seq_combine(left->row[i]+row,
 662                         left->ctx->one, left->row[i]+row,
 663                         m, left->row[row]+row,
 664                         left->n_col-row);
 665         isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
 666                         m, right->row[row], right->n_col);
 667 }
 668
 669 /* Compute inv(left)*right
 670  */
 671 struct isl_mat *isl_mat_inverse_product(struct isl_mat *left,
 672         struct isl_mat *right)
 673 {
 674         int row;
 675         isl_int a, b;
 676
 677         if (!left || !right)
 678                 goto error;
 679
 680         isl_assert(left->ctx, left->n_row == left->n_col, goto error);
 681         isl_assert(left->ctx, left->n_row == right->n_row, goto error);
 682
 683         if (left->n_row == 0) {
 684                 isl_mat_free(left);
 685                 return right;
 686         }
 687
 688         left = isl_mat_cow(left);
 689         right = isl_mat_cow(right);
 690         if (!left || !right)
 691                 goto error;
 692
 693         isl_int_init(a);
 694         isl_int_init(b);
 695         for (row = 0; row < left->n_row; ++row) {
 696                 int pivot, first, i, off;
 697                 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
 698                 if (pivot < 0) {
 699                         isl_int_clear(a);
 700                         isl_int_clear(b);
 701                         isl_assert(left->ctx, pivot >= 0, goto error);
 702                 }
 703                 pivot += row;
 704                 if (pivot != row)
 705                         inv_exchange(left, right, pivot, row);
 706                 if (isl_int_is_neg(left->row[row][row]))
 707                         inv_oppose(left, right, row);
 708                 first = row+1;
 709                 while ((off = row_first_non_zero(left->row+first,
 710                                         left->n_row-first, row)) != -1) {
 711                         first += off;
 712                         isl_int_fdiv_q(a, left->row[first][row],
 713                                         left->row[row][row]);
 714                         inv_subtract(left, right, row, first, a);
 715                         if (!isl_int_is_zero(left->row[first][row]))
 716                                 inv_exchange(left, right, row, first);
 717                         else
 718                                 ++first;
 719                 }
 720                 for (i = 0; i < row; ++i) {
 721                         if (isl_int_is_zero(left->row[i][row]))
 722                                 continue;
 723                         isl_int_gcd(a, left->row[row][row], left->row[i][row]);
 724                         isl_int_divexact(b, left->row[i][row], a);
 725                         isl_int_divexact(a, left->row[row][row], a);
 726                         isl_int_neg(b, b);
 727                         isl_seq_combine(left->row[i] + i,
 728                                         a, left->row[i] + i,
 729                                         b, left->row[row] + i,
 730                                         left->n_col - i);
 731                         isl_seq_combine(right->row[i], a, right->row[i],
 732                                         b, right->row[row], right->n_col);
 733                 }
 734         }
 735         isl_int_clear(b);
 736
 737         isl_int_set(a, left->row[0][0]);
 738         for (row = 1; row < left->n_row; ++row)
 739                 isl_int_lcm(a, a, left->row[row][row]);
 740         if (isl_int_is_zero(a)){
 741                 isl_int_clear(a);
 742                 isl_assert(left->ctx, 0, goto error);
 743         }
 744         for (row = 0; row < left->n_row; ++row) {
 745                 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
 746                 if (isl_int_is_one(left->row[row][row]))
 747                         continue;
 748                 isl_seq_scale(right->row[row], right->row[row],
 749                                 left->row[row][row], right->n_col);
 750         }
 751         isl_int_clear(a);
 752
 753         isl_mat_free(left);
 754         return right;
 755 error:
 756         isl_mat_free(left);
 757         isl_mat_free(right);
 758         return NULL;
 759 }
 760
 761 void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
 762 {
 763         int i;
 764
 765         for (i = 0; i < mat->n_row; ++i)
 766                 isl_int_mul(mat->row[i][col], mat->row[i][col], m);
 767 }
 768
 769 void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
 770         isl_int m1, unsigned src1, isl_int m2, unsigned src2)
 771 {
 772         int i;
 773         isl_int tmp;
 774
 775         isl_int_init(tmp);
 776         for (i = 0; i < mat->n_row; ++i) {
 777                 isl_int_mul(tmp, m1, mat->row[i][src1]);
 778                 isl_int_addmul(tmp, m2, mat->row[i][src2]);
 779                 isl_int_set(mat->row[i][dst], tmp);
 780         }
 781         isl_int_clear(tmp);
 782 }
 783
 784 struct isl_mat *isl_mat_right_inverse(struct isl_mat *mat)
 785 {
 786         struct isl_mat *inv;
 787         int row;
 788         isl_int a, b;
 789
 790         mat = isl_mat_cow(mat);
 791         if (!mat)
 792                 return NULL;
 793
 794         inv = isl_mat_identity(mat->ctx, mat->n_col);
 795         inv = isl_mat_cow(inv);
 796         if (!inv)
 797                 goto error;
 798
 799         isl_int_init(a);
 800         isl_int_init(b);
 801         for (row = 0; row < mat->n_row; ++row) {
 802                 int pivot, first, i, off;
 803                 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
 804                 if (pivot < 0) {
 805                         isl_int_clear(a);
 806                         isl_int_clear(b);
 807                         isl_assert(mat->ctx, pivot >= 0, goto error);
 808                 }
 809                 pivot += row;
 810                 if (pivot != row)
 811                         exchange(mat, &inv, NULL, row, pivot, row);
 812                 if (isl_int_is_neg(mat->row[row][row]))
 813                         oppose(mat, &inv, NULL, row, row);
 814                 first = row+1;
 815                 while ((off = isl_seq_first_non_zero(mat->row[row]+first,
 816                                                     mat->n_col-first)) != -1) {
 817                         first += off;
 818                         isl_int_fdiv_q(a, mat->row[row][first],
 819                                                     mat->row[row][row]);
 820                         subtract(mat, &inv, NULL, row, row, first, a);
 821                         if (!isl_int_is_zero(mat->row[row][first]))
 822                                 exchange(mat, &inv, NULL, row, row, first);
 823                         else
 824                                 ++first;
 825                 }
 826                 for (i = 0; i < row; ++i) {
 827                         if (isl_int_is_zero(mat->row[row][i]))
 828                                 continue;
 829                         isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
 830                         isl_int_divexact(b, mat->row[row][i], a);
 831                         isl_int_divexact(a, mat->row[row][row], a);
 832                         isl_int_neg(a, a);
 833                         isl_mat_col_combine(mat, i, a, i, b, row);
 834                         isl_mat_col_combine(inv, i, a, i, b, row);
 835                 }
 836         }
 837         isl_int_clear(b);
 838
 839         isl_int_set(a, mat->row[0][0]);
 840         for (row = 1; row < mat->n_row; ++row)
 841                 isl_int_lcm(a, a, mat->row[row][row]);
 842         if (isl_int_is_zero(a)){
 843                 isl_int_clear(a);
 844                 goto error;
 845         }
 846         for (row = 0; row < mat->n_row; ++row) {
 847                 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
 848                 if (isl_int_is_one(mat->row[row][row]))
 849                         continue;
 850                 isl_mat_col_scale(inv, row, mat->row[row][row]);
 851         }
 852         isl_int_clear(a);
 853
 854         isl_mat_free(mat);
 855
 856         return inv;
 857 error:
 858         isl_mat_free(mat);
 859         isl_mat_free(inv);
 860         return NULL;
 861 }
 862
 863 struct isl_mat *isl_mat_transpose(struct isl_mat *mat)
 864 {
 865         struct isl_mat *transpose = NULL;
 866         int i, j;
 867
 868         if (mat->n_col == mat->n_row) {
 869                 mat = isl_mat_cow(mat);
 870                 if (!mat)
 871                         return NULL;
 872                 for (i = 0; i < mat->n_row; ++i)
 873                         for (j = i + 1; j < mat->n_col; ++j)
 874                                 isl_int_swap(mat->row[i][j], mat->row[j][i]);
 875                 return mat;
 876         }
 877         transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);
 878         if (!transpose)
 879                 goto error;
 880         for (i = 0; i < mat->n_row; ++i)
 881                 for (j = 0; j < mat->n_col; ++j)
 882                         isl_int_set(transpose->row[j][i], mat->row[i][j]);
 883         isl_mat_free(mat);
 884         return transpose;
 885 error:
 886         isl_mat_free(mat);
 887         return NULL;
 888 }
 889
 890 struct isl_mat *isl_mat_swap_cols(struct isl_mat *mat, unsigned i, unsigned j)
 891 {
 892         int r;
 893
 894         mat = isl_mat_cow(mat);
 895         if (!mat)
 896                 return NULL;
 897         isl_assert(mat->ctx, i < mat->n_col, goto error);
 898         isl_assert(mat->ctx, j < mat->n_col, goto error);
 899
 900         for (r = 0; r < mat->n_row; ++r)
 901                 isl_int_swap(mat->row[r][i], mat->row[r][j]);
 902         return mat;
 903 error:
 904         isl_mat_free(mat);
 905         return NULL;
 906 }
 907
 908 struct isl_mat *isl_mat_swap_rows(struct isl_mat *mat, unsigned i, unsigned j)
 909 {
 910         isl_int *t;
 911
 912         if (!mat)
 913                 return NULL;
 914         mat = isl_mat_cow(mat);
 915         if (!mat)
 916                 return NULL;
 917         t = mat->row[i];
 918         mat->row[i] = mat->row[j];
 919         mat->row[j] = t;
 920         return mat;
 921 }
 922
 923 struct isl_mat *isl_mat_product(struct isl_mat *left, struct isl_mat *right)
 924 {
 925         int i, j, k;
 926         struct isl_mat *prod;
 927
 928         if (!left || !right)
 929                 goto error;
 930         isl_assert(left->ctx, left->n_col == right->n_row, goto error);
 931         prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
 932         if (!prod)
 933                 goto error;
 934         if (left->n_col == 0) {
 935                 for (i = 0; i < prod->n_row; ++i)
 936                         isl_seq_clr(prod->row[i], prod->n_col);
 937                 return prod;
 938         }
 939         for (i = 0; i < prod->n_row; ++i) {
 940                 for (j = 0; j < prod->n_col; ++j) {
 941                         isl_int_mul(prod->row[i][j],
 942                                     left->row[i][0], right->row[0][j]);
 943                         for (k = 1; k < left->n_col; ++k)
 944                                 isl_int_addmul(prod->row[i][j],
 945                                             left->row[i][k], right->row[k][j]);
 946                 }
 947         }
 948         isl_mat_free(left);
 949         isl_mat_free(right);
 950         return prod;
 951 error:
 952         isl_mat_free(left);
 953         isl_mat_free(right);
 954         return NULL;
 955 }
 956
 957 /* Replace the variables x in the rows q by x' given by x = M x',
 958  * with M the matrix mat.
 959  *
 960  * If the number of new variables is greater than the original
 961  * number of variables, then the rows q have already been
 962  * preextended.  If the new number is smaller, then the coefficients
 963  * of the divs, which are not changed, need to be shifted down.
 964  * The row q may be the equalities, the inequalities or the
 965  * div expressions.  In the latter case, has_div is true and
 966  * we need to take into account the extra denominator column.
 967  */
 968 static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,
 969         unsigned n_div, int has_div, struct isl_mat *mat)
 970 {
 971         int i;
 972         struct isl_mat *t;
 973         int e;
 974
 975         if (mat->n_col >= mat->n_row)
 976                 e = 0;
 977         else
 978                 e = mat->n_row - mat->n_col;
 979         if (has_div)
 980                 for (i = 0; i < n; ++i)
 981                         isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);
 982         t = isl_mat_sub_alloc(mat->ctx, q, 0, n, has_div, mat->n_row);
 983         t = isl_mat_product(t, mat);
 984         if (!t)
 985                 return -1;
 986         for (i = 0; i < n; ++i) {
 987                 isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
 988                 isl_seq_cpy(q[i] + has_div + t->n_col,
 989                             q[i] + has_div + t->n_col + e, n_div);
 990                 isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
 991         }
 992         isl_mat_free(t);
 993         return 0;
 994 }
 995
 996 /* Replace the variables x in bset by x' given by x = M x', with
 997  * M the matrix mat.
 998  *
 999  * If there are fewer variables x' then there are x, then we perform
1000  * the transformation in place, which that, in principle,
1001  * this frees up some extra variables as the number
1002  * of columns remains constant, but we would have to extend
1003  * the div array too as the number of rows in this array is assumed
1004  * to be equal to extra.
1005  */
1006 struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset,
1007         struct isl_mat *mat)
1008 {
1009         struct isl_ctx *ctx;
1010
1011         if (!bset || !mat)
1012                 goto error;
1013
1014         ctx = bset->ctx;
1015         bset = isl_basic_set_cow(bset);
1016         if (!bset)
1017                 goto error;
1018
1019         isl_assert(ctx, bset->dim->nparam == 0, goto error);
1020         isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1021         isl_assert(ctx, mat->n_col > 0, goto error);
1022
1023         if (mat->n_col > mat->n_row) {
1024                 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1025                 if (!bset)
1026                         goto error;
1027         } else if (mat->n_col < mat->n_row) {
1028                 bset->dim = isl_dim_cow(bset->dim);
1029                 if (!bset->dim)
1030                         goto error;
1031                 bset->dim->n_out -= mat->n_row - mat->n_col;
1032         }
1033
1034         if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
1035                         isl_mat_copy(mat)) < 0)
1036                 goto error;
1037
1038         if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
1039                         isl_mat_copy(mat)) < 0)
1040                 goto error;
1041
1042         if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)
1043                 goto error2;
1044
1045         ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
1046         ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1047         ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1048         ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1049         ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1050
1051         bset = isl_basic_set_simplify(bset);
1052         bset = isl_basic_set_finalize(bset);
1053
1054         return bset;
1055 error:
1056         isl_mat_free(mat);
1057 error2:
1058         isl_basic_set_free(bset);
1059         return NULL;
1060 }
1061
1062 struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat)
1063 {
1064         struct isl_ctx *ctx;
1065         int i;
1066
1067         set = isl_set_cow(set);
1068         if (!set)
1069                 return NULL;
1070
1071         ctx = set->ctx;
1072         for (i = 0; i < set->n; ++i) {
1073                 set->p[i] = isl_basic_set_preimage(set->p[i],
1074                                                     isl_mat_copy(mat));
1075                 if (!set->p[i])
1076                         goto error;
1077         }
1078         if (mat->n_col != mat->n_row) {
1079                 set->dim = isl_dim_cow(set->dim);
1080                 if (!set->dim)
1081                         goto error;
1082                 set->dim->n_out += mat->n_col;
1083                 set->dim->n_out -= mat->n_row;
1084         }
1085         isl_mat_free(mat);
1086         ISL_F_CLR(set, ISL_SET_NORMALIZED);
1087         return set;
1088 error:
1089         isl_set_free(set);
1090         isl_mat_free(mat);
1091         return NULL;
1092 }
1093
1094 void isl_mat_dump(struct isl_mat *mat, FILE *out, int indent)
1095 {
1096         int i, j;
1097
1098         if (!mat) {
1099                 fprintf(out, "%*snull mat\n", indent, "");
1100                 return;
1101         }
1102
1103         if (mat->n_row == 0)
1104                 fprintf(out, "%*s[]\n", indent, "");
1105
1106         for (i = 0; i < mat->n_row; ++i) {
1107                 if (!i)
1108                         fprintf(out, "%*s[[", indent, "");
1109                 else
1110                         fprintf(out, "%*s[", indent+1, "");
1111                 for (j = 0; j < mat->n_col; ++j) {
1112                         if (j)
1113                             fprintf(out, ",");
1114                         isl_int_print(out, mat->row[i][j], 0);
1115                 }
1116                 if (i == mat->n_row-1)
1117                         fprintf(out, "]]\n");
1118                 else
1119                         fprintf(out, "]\n");
1120         }
1121 }
1122
1123 struct isl_mat *isl_mat_drop_cols(struct isl_mat *mat, unsigned col, unsigned n)
1124 {
1125         int r;
1126
1127         mat = isl_mat_cow(mat);
1128         if (!mat)
1129                 return NULL;
1130
1131         if (col != mat->n_col-n) {
1132                 for (r = 0; r < mat->n_row; ++r)
1133                         isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1134                                         mat->n_col - col - n);
1135         }
1136         mat->n_col -= n;
1137         return mat;
1138 }
1139
1140 struct isl_mat *isl_mat_drop_rows(struct isl_mat *mat, unsigned row, unsigned n)
1141 {
1142         int r;
1143
1144         mat = isl_mat_cow(mat);
1145         if (!mat)
1146                 return NULL;
1147
1148         for (r = row; r+n < mat->n_row; ++r)
1149                 mat->row[r] = mat->row[r+n];
1150
1151         mat->n_row -= n;
1152         return mat;
1153 }
1154
1155 __isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1156                                 unsigned col, unsigned n)
1157 {
1158         isl_mat *ext;
1159
1160         if (!mat)
1161                 return NULL;
1162         if (n == 0)
1163                 return mat;
1164
1165         ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1166         if (!ext)
1167                 goto error;
1168
1169         isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1170         isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1171                                 col + n, col, mat->n_col - col);
1172
1173         isl_mat_free(mat);
1174         return ext;
1175 error:
1176         isl_mat_free(mat);
1177         return NULL;
1178 }
1179
1180 __isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1181         unsigned first, unsigned n)
1182 {
1183         int i;
1184
1185         if (!mat)
1186                 return NULL;
1187         mat = isl_mat_insert_cols(mat, first, n);
1188         if (!mat)
1189                 return NULL;
1190
1191         for (i = 0; i < mat->n_row; ++i)
1192                 isl_seq_clr(mat->row[i] + first, n);
1193
1194         return mat;
1195 }
1196
1197 __isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1198 {
1199         if (!mat)
1200                 return NULL;
1201
1202         return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1203 }
1204
1205 __isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1206                                 unsigned row, unsigned n)
1207 {
1208         isl_mat *ext;
1209
1210         if (!mat)
1211                 return NULL;
1212         if (n == 0)
1213                 return mat;
1214
1215         ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1216         if (!ext)
1217                 goto error;
1218
1219         isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1220         isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1221                                 mat->n_row - row, 0, 0, mat->n_col);
1222
1223         isl_mat_free(mat);
1224         return ext;
1225 error:
1226         isl_mat_free(mat);
1227         return NULL;
1228 }
1229
1230 __isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1231 {
1232         if (!mat)
1233                 return NULL;
1234
1235         return isl_mat_insert_rows(mat, mat->n_row, n);
1236 }
1237
1238 void isl_mat_col_submul(struct isl_mat *mat,
1239                         int dst_col, isl_int f, int src_col)
1240 {
1241         int i;
1242
1243         for (i = 0; i < mat->n_row; ++i)
1244                 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1245 }
1246
1247 void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)
1248 {
1249         int i;
1250
1251         if (!mat)
1252                 return;
1253
1254         for (i = 0; i < mat->n_row; ++i)
1255                 isl_int_add(mat->row[i][dst_col],
1256                             mat->row[i][dst_col], mat->row[i][src_col]);
1257 }
1258
1259 void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1260 {
1261         int i;
1262
1263         for (i = 0; i < mat->n_row; ++i)
1264                 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1265 }
1266
1267 struct isl_mat *isl_mat_unimodular_complete(struct isl_mat *M, int row)
1268 {
1269         int r;
1270         struct isl_mat *H = NULL, *Q = NULL;
1271
1272         if (!M)
1273                 return NULL;
1274
1275         isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1276         M->n_row = row;
1277         H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1278         M->n_row = M->n_col;
1279         if (!H)
1280                 goto error;
1281         for (r = 0; r < row; ++r)
1282                 isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1283         for (r = row; r < M->n_row; ++r)
1284                 isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1285         isl_mat_free(H);
1286         isl_mat_free(Q);
1287         return M;
1288 error:
1289         isl_mat_free(H);
1290         isl_mat_free(Q);
1291         isl_mat_free(M);
1292         return NULL;
1293 }
1294
1295 __isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1296         __isl_take isl_mat *bot)
1297 {
1298         struct isl_mat *mat;
1299
1300         if (!top || !bot)
1301                 goto error;
1302
1303         isl_assert(top->ctx, top->n_col == bot->n_col, goto error);
1304         if (top->n_row == 0) {
1305                 isl_mat_free(top);
1306                 return bot;
1307         }
1308         if (bot->n_row == 0) {
1309                 isl_mat_free(bot);
1310                 return top;
1311         }
1312
1313         mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1314         if (!mat)
1315                 goto error;
1316         isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1317                          0, 0, mat->n_col);
1318         isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1319                          0, 0, mat->n_col);
1320         isl_mat_free(top);
1321         isl_mat_free(bot);
1322         return mat;
1323 error:
1324         isl_mat_free(top);
1325         isl_mat_free(bot);
1326         return NULL;
1327 }
1328
1329 int isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1330 {
1331         int i;
1332
1333         if (!mat1 || !mat2)
1334                 return -1;
1335
1336         if (mat1->n_row != mat2->n_row)
1337                 return 0;
1338
1339         if (mat1->n_col != mat2->n_col)
1340                 return 0;
1341
1342         for (i = 0; i < mat1->n_row; ++i)
1343                 if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col))
1344                         return 0;
1345
1346         return 1;
1347 }
1348
1349 __isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)
1350 {
1351         struct isl_mat *mat;
1352
1353         if (!vec)
1354                 return NULL;
1355         mat = isl_mat_alloc(vec->ctx, 1, vec->size);
1356         if (!mat)
1357                 goto error;
1358
1359         isl_seq_cpy(mat->row[0], vec->el, vec->size);
1360
1361         isl_vec_free(vec);
1362         return mat;
1363 error:
1364         isl_vec_free(vec);
1365         return NULL;
1366 }
1367
1368 __isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,
1369         __isl_take isl_vec *bot)
1370 {
1371         return isl_mat_concat(top, isl_mat_from_row_vec(bot));
1372 }
1373
1374 __isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,
1375         unsigned dst_col, unsigned src_col, unsigned n)
1376 {
1377         isl_mat *res;
1378
1379         if (!mat)
1380                 return NULL;
1381         if (n == 0 || dst_col == src_col)
1382                 return mat;
1383
1384         res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1385         if (!res)
1386                 goto error;
1387
1388         if (dst_col < src_col) {
1389                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1390                                  0, 0, dst_col);
1391                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1392                                  dst_col, src_col, n);
1393                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1394                                  dst_col + n, dst_col, src_col - dst_col);
1395                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1396                                  src_col + n, src_col + n,
1397                                  res->n_col - src_col - n);
1398         } else {
1399                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1400                                  0, 0, src_col);
1401                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1402                                  src_col, src_col + n, dst_col - src_col);
1403                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1404                                  dst_col, src_col, n);
1405                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1406                                  dst_col + n, dst_col + n,
1407                                  res->n_col - dst_col - n);
1408         }
1409         isl_mat_free(mat);
1410
1411         return res;
1412 error:
1413         isl_mat_free(mat);
1414         return NULL;
1415 }
1416
1417 void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1418 {
1419         int i;
1420         isl_int g;
1421
1422         isl_int_set_si(*gcd, 0);
1423         if (!mat)
1424                 return;
1425
1426         isl_int_init(g);
1427         for (i = 0; i < mat->n_row; ++i) {
1428                 isl_seq_gcd(mat->row[i], mat->n_col, &g);
1429                 isl_int_gcd(*gcd, *gcd, g);
1430         }
1431         isl_int_clear(g);
1432 }
1433
1434 __isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1435 {
1436         int i;
1437
1438         if (!mat)
1439                 return NULL;
1440
1441         for (i = 0; i < mat->n_row; ++i)
1442                 isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1443
1444         return mat;
1445 }
1446
1447 __isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1448 {
1449         isl_int gcd;
1450
1451         if (!mat)
1452                 return NULL;
1453
1454         isl_int_init(gcd);
1455         isl_mat_gcd(mat, &gcd);
1456         mat = isl_mat_scale_down(mat, gcd);
1457         isl_int_clear(gcd);
1458
1459         return mat;
1460 }