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