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