isl_mat.c

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