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