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 void isl_mat_dump(struct isl_mat *mat, FILE *out, int indent)
1129 {
1130         int i, j;
1131
1132         if (!mat) {
1133                 fprintf(out, "%*snull mat\n", indent, "");
1134                 return;
1135         }
1136
1137         if (mat->n_row == 0)
1138                 fprintf(out, "%*s[]\n", indent, "");
1139
1140         for (i = 0; i < mat->n_row; ++i) {
1141                 if (!i)
1142                         fprintf(out, "%*s[[", indent, "");
1143                 else
1144                         fprintf(out, "%*s[", indent+1, "");
1145                 for (j = 0; j < mat->n_col; ++j) {
1146                         if (j)
1147                             fprintf(out, ",");
1148                         isl_int_print(out, mat->row[i][j], 0);
1149                 }
1150                 if (i == mat->n_row-1)
1151                         fprintf(out, "]]\n");
1152                 else
1153                         fprintf(out, "]\n");
1154         }
1155 }
1156
1157 struct isl_mat *isl_mat_drop_cols(struct isl_mat *mat, unsigned col, unsigned n)
1158 {
1159         int r;
1160
1161         mat = isl_mat_cow(mat);
1162         if (!mat)
1163                 return NULL;
1164
1165         if (col != mat->n_col-n) {
1166                 for (r = 0; r < mat->n_row; ++r)
1167                         isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1168                                         mat->n_col - col - n);
1169         }
1170         mat->n_col -= n;
1171         return mat;
1172 }
1173
1174 struct isl_mat *isl_mat_drop_rows(struct isl_mat *mat, unsigned row, unsigned n)
1175 {
1176         int r;
1177
1178         mat = isl_mat_cow(mat);
1179         if (!mat)
1180                 return NULL;
1181
1182         for (r = row; r+n < mat->n_row; ++r)
1183                 mat->row[r] = mat->row[r+n];
1184
1185         mat->n_row -= n;
1186         return mat;
1187 }
1188
1189 __isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1190                                 unsigned col, unsigned n)
1191 {
1192         isl_mat *ext;
1193
1194         if (!mat)
1195                 return NULL;
1196         if (n == 0)
1197                 return mat;
1198
1199         ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1200         if (!ext)
1201                 goto error;
1202
1203         isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1204         isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1205                                 col + n, col, mat->n_col - col);
1206
1207         isl_mat_free(mat);
1208         return ext;
1209 error:
1210         isl_mat_free(mat);
1211         return NULL;
1212 }
1213
1214 __isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1215         unsigned first, unsigned n)
1216 {
1217         int i;
1218
1219         if (!mat)
1220                 return NULL;
1221         mat = isl_mat_insert_cols(mat, first, n);
1222         if (!mat)
1223                 return NULL;
1224
1225         for (i = 0; i < mat->n_row; ++i)
1226                 isl_seq_clr(mat->row[i] + first, n);
1227
1228         return mat;
1229 }
1230
1231 __isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1232 {
1233         if (!mat)
1234                 return NULL;
1235
1236         return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1237 }
1238
1239 __isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1240                                 unsigned row, unsigned n)
1241 {
1242         isl_mat *ext;
1243
1244         if (!mat)
1245                 return NULL;
1246         if (n == 0)
1247                 return mat;
1248
1249         ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1250         if (!ext)
1251                 goto error;
1252
1253         isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1254         isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1255                                 mat->n_row - row, 0, 0, mat->n_col);
1256
1257         isl_mat_free(mat);
1258         return ext;
1259 error:
1260         isl_mat_free(mat);
1261         return NULL;
1262 }
1263
1264 __isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1265 {
1266         if (!mat)
1267                 return NULL;
1268
1269         return isl_mat_insert_rows(mat, mat->n_row, n);
1270 }
1271
1272 void isl_mat_col_submul(struct isl_mat *mat,
1273                         int dst_col, isl_int f, int src_col)
1274 {
1275         int i;
1276
1277         for (i = 0; i < mat->n_row; ++i)
1278                 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1279 }
1280
1281 void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)
1282 {
1283         int i;
1284
1285         if (!mat)
1286                 return;
1287
1288         for (i = 0; i < mat->n_row; ++i)
1289                 isl_int_add(mat->row[i][dst_col],
1290                             mat->row[i][dst_col], mat->row[i][src_col]);
1291 }
1292
1293 void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1294 {
1295         int i;
1296
1297         for (i = 0; i < mat->n_row; ++i)
1298                 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1299 }
1300
1301 struct isl_mat *isl_mat_unimodular_complete(struct isl_mat *M, int row)
1302 {
1303         int r;
1304         struct isl_mat *H = NULL, *Q = NULL;
1305
1306         if (!M)
1307                 return NULL;
1308
1309         isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1310         M->n_row = row;
1311         H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1312         M->n_row = M->n_col;
1313         if (!H)
1314                 goto error;
1315         for (r = 0; r < row; ++r)
1316                 isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1317         for (r = row; r < M->n_row; ++r)
1318                 isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1319         isl_mat_free(H);
1320         isl_mat_free(Q);
1321         return M;
1322 error:
1323         isl_mat_free(H);
1324         isl_mat_free(Q);
1325         isl_mat_free(M);
1326         return NULL;
1327 }
1328
1329 __isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1330         __isl_take isl_mat *bot)
1331 {
1332         struct isl_mat *mat;
1333
1334         if (!top || !bot)
1335                 goto error;
1336
1337         isl_assert(top->ctx, top->n_col == bot->n_col, goto error);
1338         if (top->n_row == 0) {
1339                 isl_mat_free(top);
1340                 return bot;
1341         }
1342         if (bot->n_row == 0) {
1343                 isl_mat_free(bot);
1344                 return top;
1345         }
1346
1347         mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1348         if (!mat)
1349                 goto error;
1350         isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1351                          0, 0, mat->n_col);
1352         isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1353                          0, 0, mat->n_col);
1354         isl_mat_free(top);
1355         isl_mat_free(bot);
1356         return mat;
1357 error:
1358         isl_mat_free(top);
1359         isl_mat_free(bot);
1360         return NULL;
1361 }
1362
1363 int isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1364 {
1365         int i;
1366
1367         if (!mat1 || !mat2)
1368                 return -1;
1369
1370         if (mat1->n_row != mat2->n_row)
1371                 return 0;
1372
1373         if (mat1->n_col != mat2->n_col)
1374                 return 0;
1375
1376         for (i = 0; i < mat1->n_row; ++i)
1377                 if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col))
1378                         return 0;
1379
1380         return 1;
1381 }
1382
1383 __isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)
1384 {
1385         struct isl_mat *mat;
1386
1387         if (!vec)
1388                 return NULL;
1389         mat = isl_mat_alloc(vec->ctx, 1, vec->size);
1390         if (!mat)
1391                 goto error;
1392
1393         isl_seq_cpy(mat->row[0], vec->el, vec->size);
1394
1395         isl_vec_free(vec);
1396         return mat;
1397 error:
1398         isl_vec_free(vec);
1399         return NULL;
1400 }
1401
1402 __isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,
1403         __isl_take isl_vec *bot)
1404 {
1405         return isl_mat_concat(top, isl_mat_from_row_vec(bot));
1406 }
1407
1408 __isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,
1409         unsigned dst_col, unsigned src_col, unsigned n)
1410 {
1411         isl_mat *res;
1412
1413         if (!mat)
1414                 return NULL;
1415         if (n == 0 || dst_col == src_col)
1416                 return mat;
1417
1418         res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1419         if (!res)
1420                 goto error;
1421
1422         if (dst_col < src_col) {
1423                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1424                                  0, 0, dst_col);
1425                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1426                                  dst_col, src_col, n);
1427                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1428                                  dst_col + n, dst_col, src_col - dst_col);
1429                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1430                                  src_col + n, src_col + n,
1431                                  res->n_col - src_col - n);
1432         } else {
1433                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1434                                  0, 0, src_col);
1435                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1436                                  src_col, src_col + n, dst_col - src_col);
1437                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1438                                  dst_col, src_col, n);
1439                 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1440                                  dst_col + n, dst_col + n,
1441                                  res->n_col - dst_col - n);
1442         }
1443         isl_mat_free(mat);
1444
1445         return res;
1446 error:
1447         isl_mat_free(mat);
1448         return NULL;
1449 }
1450
1451 void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1452 {
1453         int i;
1454         isl_int g;
1455
1456         isl_int_set_si(*gcd, 0);
1457         if (!mat)
1458                 return;
1459
1460         isl_int_init(g);
1461         for (i = 0; i < mat->n_row; ++i) {
1462                 isl_seq_gcd(mat->row[i], mat->n_col, &g);
1463                 isl_int_gcd(*gcd, *gcd, g);
1464         }
1465         isl_int_clear(g);
1466 }
1467
1468 __isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1469 {
1470         int i;
1471
1472         if (!mat)
1473                 return NULL;
1474
1475         for (i = 0; i < mat->n_row; ++i)
1476                 isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1477
1478         return mat;
1479 }
1480
1481 __isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1482 {
1483         isl_int gcd;
1484
1485         if (!mat)
1486                 return NULL;
1487
1488         isl_int_init(gcd);
1489         isl_mat_gcd(mat, &gcd);
1490         mat = isl_mat_scale_down(mat, gcd);
1491         isl_int_clear(gcd);
1492
1493         return mat;
1494 }
1495
1496 /* Number of initial non-zero columns.
1497  */
1498 int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
1499 {
1500         int i;
1501
1502         if (!mat)
1503                 return -1;
1504
1505         for (i = 0; i < mat->n_col; ++i)
1506                 if (row_first_non_zero(mat->row, mat->n_row, i) < 0)
1507                         break;
1508
1509         return i;
1510 }