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