isl_mat.c

   1 #include "isl_seq.h"
   2 #include "isl_mat.h"
   3 #include "isl_map_private.h"
   4
   5 struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
   6         unsigned n_row, unsigned n_col)
   7 {
   8         int i;
   9         struct isl_mat *mat;
  10
  11         mat = isl_alloc_type(ctx, struct isl_mat);
  12         if (!mat)
  13                 return NULL;
  14
  15         mat->row = NULL;
  16         mat->block = isl_blk_alloc(ctx, n_row * n_col);
  17         if (isl_blk_is_error(mat->block))
  18                 goto error;
  19         mat->row = isl_alloc_array(ctx, isl_int *, n_row);
  20         if (!mat->row)
  21                 goto error;
  22
  23         for (i = 0; i < n_row; ++i)
  24                 mat->row[i] = mat->block.data + i * n_col;
  25
  26         mat->ref = 1;
  27         mat->n_row = n_row;
  28         mat->n_col = n_col;
  29         mat->flags = 0;
  30
  31         return mat;
  32 error:
  33         isl_blk_free(ctx, mat->block);
  34         free(mat);
  35         return NULL;
  36 }
  37
  38 struct isl_mat *isl_mat_sub_alloc(struct isl_ctx *ctx, isl_int **row,
  39         unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
  40 {
  41         int i;
  42         struct isl_mat *mat;
  43
  44         mat = isl_alloc_type(ctx, struct isl_mat);
  45         if (!mat)
  46                 return NULL;
  47         mat->row = isl_alloc_array(ctx, isl_int *, n_row);
  48         if (!mat->row)
  49                 goto error;
  50         for (i = 0; i < n_row; ++i)
  51                 mat->row[i] = row[first_row+i] + first_col;
  52         mat->ref = 1;
  53         mat->n_row = n_row;
  54         mat->n_col = n_col;
  55         mat->block = isl_blk_empty();
  56         mat->flags = ISL_MAT_BORROWED;
  57         return mat;
  58 error:
  59         free(mat);
  60         return NULL;
  61 }
  62
  63 void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
  64         unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
  65 {
  66         int i;
  67
  68         for (i = 0; i < n_row; ++i)
  69                 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
  70 }
  71
  72 void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
  73         unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
  74 {
  75         int i;
  76
  77         for (i = 0; i < n_row; ++i)
  78                 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
  79 }
  80
  81 struct isl_mat *isl_mat_copy(struct isl_ctx *ctx, struct isl_mat *mat)
  82 {
  83         if (!mat)
  84                 return NULL;
  85
  86         mat->ref++;
  87         return mat;
  88 }
  89
  90 struct isl_mat *isl_mat_dup(struct isl_ctx *ctx, struct isl_mat *mat)
  91 {
  92         int i;
  93         struct isl_mat *mat2;
  94
  95         if (!mat)
  96                 return NULL;
  97         mat2 = isl_mat_alloc(ctx, mat->n_row, mat->n_col);
  98         if (!mat2)
  99                 return NULL;
 100         for (i = 0; i < mat->n_row; ++i)
 101                 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
 102         return mat2;
 103 }
 104
 105 struct isl_mat *isl_mat_cow(struct isl_ctx *ctx, struct isl_mat *mat)
 106 {
 107         struct isl_mat *mat2;
 108         if (!mat)
 109                 return NULL;
 110
 111         if (mat->ref == 1 && !F_ISSET(mat, ISL_MAT_BORROWED))
 112                 return mat;
 113
 114         mat2 = isl_mat_dup(ctx, mat);
 115         isl_mat_free(ctx, mat);
 116         return mat2;
 117 }
 118
 119 void isl_mat_free(struct isl_ctx *ctx, struct isl_mat *mat)
 120 {
 121         if (!mat)
 122                 return;
 123
 124         if (--mat->ref > 0)
 125                 return;
 126
 127         if (!F_ISSET(mat, ISL_MAT_BORROWED))
 128                 isl_blk_free(ctx, mat->block);
 129         free(mat->row);
 130         free(mat);
 131 }
 132
 133 struct isl_mat *isl_mat_identity(struct isl_ctx *ctx, unsigned n_row)
 134 {
 135         int i;
 136         struct isl_mat *mat;
 137
 138         mat = isl_mat_alloc(ctx, n_row, n_row);
 139         if (!mat)
 140                 return NULL;
 141         for (i = 0; i < n_row; ++i) {
 142                 isl_seq_clr(mat->row[i], i);
 143                 isl_int_set_si(mat->row[i][i], 1);
 144                 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
 145         }
 146
 147         return mat;
 148 }
 149
 150 struct isl_vec *isl_mat_vec_product(struct isl_ctx *ctx,
 151         struct isl_mat *mat, struct isl_vec *vec)
 152 {
 153         int i;
 154         struct isl_vec *prod;
 155
 156         if (!mat || !vec)
 157                 goto error;
 158
 159         isl_assert(ctx, mat->n_col == vec->size, goto error);
 160
 161         prod = isl_vec_alloc(ctx, mat->n_row);
 162         if (!prod)
 163                 goto error;
 164
 165         for (i = 0; i < prod->size; ++i)
 166                 isl_seq_inner_product(mat->row[i], vec->block.data, vec->size,
 167                                         &prod->block.data[i]);
 168         isl_mat_free(ctx, mat);
 169         isl_vec_free(ctx, vec);
 170         return prod;
 171 error:
 172         isl_mat_free(ctx, mat);
 173         isl_vec_free(ctx, vec);
 174         return NULL;
 175 }
 176
 177 struct isl_mat *isl_mat_aff_direct_sum(struct isl_ctx *ctx,
 178         struct isl_mat *left, struct isl_mat *right)
 179 {
 180         int i;
 181         struct isl_mat *sum;
 182
 183         if (!left || !right)
 184                 goto error;
 185
 186         isl_assert(ctx, left->n_row == right->n_row, goto error);
 187         isl_assert(ctx, left->n_row >= 1, goto error);
 188         isl_assert(ctx, left->n_col >= 1, goto error);
 189         isl_assert(ctx, right->n_col >= 1, goto error);
 190         isl_assert(ctx,
 191             isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
 192             goto error);
 193         isl_assert(ctx,
 194             isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
 195             goto error);
 196
 197         sum = isl_mat_alloc(ctx, left->n_row, left->n_col + right->n_col - 1);
 198         if (!sum)
 199                 goto error;
 200         isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
 201         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
 202         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
 203
 204         isl_seq_clr(sum->row[0]+1, sum->n_col-1);
 205         for (i = 1; i < sum->n_row; ++i) {
 206                 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
 207                 isl_int_addmul(sum->row[i][0],
 208                                 right->row[0][0], right->row[i][0]);
 209                 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
 210                                 left->n_col-1);
 211                 isl_seq_scale(sum->row[i]+left->n_col,
 212                                 right->row[i]+1, right->row[0][0],
 213                                 right->n_col-1);
 214         }
 215
 216         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
 217         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
 218         isl_mat_free(ctx, left);
 219         isl_mat_free(ctx, right);
 220         return sum;
 221 error:
 222         isl_mat_free(ctx, left);
 223         isl_mat_free(ctx, right);
 224         return NULL;
 225 }
 226
 227 static void exchange(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
 228         struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
 229 {
 230         int r;
 231         for (r = row; r < M->n_row; ++r)
 232                 isl_int_swap(M->row[r][i], M->row[r][j]);
 233         if (U) {
 234                 for (r = 0; r < (*U)->n_row; ++r)
 235                         isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
 236         }
 237         if (Q)
 238                 isl_mat_swap_rows(ctx, *Q, i, j);
 239 }
 240
 241 static void subtract(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
 242         struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
 243 {
 244         int r, c;
 245         for (r = row; r < M->n_row; ++r)
 246                 isl_int_submul(M->row[r][j], m, M->row[r][i]);
 247         if (U) {
 248                 for (r = 0; r < (*U)->n_row; ++r)
 249                         isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
 250         }
 251         if (Q) {
 252                 for (r = 0; r < (*Q)->n_col; ++r)
 253                         isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
 254         }
 255 }
 256
 257 static void oppose(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
 258         struct isl_mat **Q, unsigned row, unsigned col)
 259 {
 260         int r;
 261         for (r = row; r < M->n_row; ++r)
 262                 isl_int_neg(M->row[r][col], M->row[r][col]);
 263         if (U) {
 264                 for (r = 0; r < (*U)->n_row; ++r)
 265                         isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
 266         }
 267         if (Q)
 268                 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
 269 }
 270
 271 /* Given matrix M, compute
 272  *
 273  *              M U = H
 274  *              M   = H Q
 275  *
 276  * with U and Q unimodular matrices and H a matrix in column echelon form
 277  * such that on each echelon row the entries in the non-echelon column
 278  * are non-negative and stricly smaller than the entries in the echelon
 279  * column.  If U or Q are NULL, then these matrices are not computed.
 280  */
 281 struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
 282         struct isl_mat *M, struct isl_mat **U, struct isl_mat **Q)
 283 {
 284         isl_int c;
 285         int row, col;
 286
 287         if (U)
 288                 *U = NULL;
 289         if (Q)
 290                 *Q = NULL;
 291         if (!M)
 292                 goto error;
 293         M = isl_mat_cow(ctx, M);
 294         if (!M)
 295                 goto error;
 296         if (U) {
 297                 *U = isl_mat_identity(ctx, M->n_col);
 298                 if (!*U)
 299                         goto error;
 300         }
 301         if (Q) {
 302                 *Q = isl_mat_identity(ctx, M->n_col);
 303                 if (!*Q)
 304                         goto error;
 305         }
 306
 307         col = 0;
 308         isl_int_init(c);
 309         for (row = 0; row < M->n_row; ++row) {
 310                 int first, i, off;
 311                 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
 312                 if (first == -1)
 313                         continue;
 314                 first += col;
 315                 if (first != col)
 316                         exchange(ctx, M, U, Q, row, first, col);
 317                 if (isl_int_is_neg(M->row[row][col]))
 318                         oppose(ctx, M, U, Q, row, col);
 319                 first = col+1;
 320                 while ((off = isl_seq_first_non_zero(M->row[row]+first,
 321                                                        M->n_col-first)) != -1) {
 322                         first += off;
 323                         isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
 324                         subtract(ctx, M, U, Q, row, col, first, c);
 325                         if (!isl_int_is_zero(M->row[row][first]))
 326                                 exchange(ctx, M, U, Q, row, first, col);
 327                         else
 328                                 ++first;
 329                 }
 330                 for (i = 0; i < col; ++i) {
 331                         if (isl_int_is_zero(M->row[row][i]))
 332                                 continue;
 333                         isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
 334                         if (isl_int_is_zero(c))
 335                                 continue;
 336                         subtract(ctx, M, U, Q, row, col, i, c);
 337                 }
 338                 ++col;
 339         }
 340         isl_int_clear(c);
 341
 342         return M;
 343 error:
 344         if (Q) {
 345                 isl_mat_free(ctx, *Q);
 346                 *Q = NULL;
 347         }
 348         if (U) {
 349                 isl_mat_free(ctx, *U);
 350                 *U = NULL;
 351         }
 352         return NULL;
 353 }
 354
 355 struct isl_mat *isl_mat_lin_to_aff(struct isl_ctx *ctx, struct isl_mat *mat)
 356 {
 357         int i;
 358         struct isl_mat *mat2;
 359
 360         if (!mat)
 361                 return NULL;
 362         mat2 = isl_mat_alloc(ctx, 1+mat->n_row, 1+mat->n_col);
 363         if (!mat2)
 364                 return NULL;
 365         isl_int_set_si(mat2->row[0][0], 1);
 366         isl_seq_clr(mat2->row[0]+1, mat->n_col);
 367         for (i = 0; i < mat->n_row; ++i) {
 368                 isl_int_set_si(mat2->row[1+i][0], 0);
 369                 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
 370         }
 371         isl_mat_free(ctx, mat);
 372         return mat2;
 373 }
 374
 375 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
 376 {
 377         int i;
 378
 379         for (i = 0; i < n_row; ++i)
 380                 if (!isl_int_is_zero(row[i][col]))
 381                         return i;
 382         return -1;
 383 }
 384
 385 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
 386 {
 387         int i, min = row_first_non_zero(row, n_row, col);
 388         if (min < 0)
 389                 return -1;
 390         for (i = min + 1; i < n_row; ++i) {
 391                 if (isl_int_is_zero(row[i][col]))
 392                         continue;
 393                 if (isl_int_abs_lt(row[i][col], row[min][col]))
 394                         min = i;
 395         }
 396         return min;
 397 }
 398
 399 static void inv_exchange(struct isl_ctx *ctx,
 400         struct isl_mat *left, struct isl_mat *right,
 401         unsigned i, unsigned j)
 402 {
 403         left = isl_mat_swap_rows(ctx, left, i, j);
 404         right = isl_mat_swap_rows(ctx, right, i, j);
 405 }
 406
 407 static void inv_oppose(struct isl_ctx *ctx,
 408         struct isl_mat *left, struct isl_mat *right, unsigned row)
 409 {
 410         isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
 411         isl_seq_neg(right->row[row], right->row[row], right->n_col);
 412 }
 413
 414 static void inv_subtract(struct isl_ctx *ctx,
 415         struct isl_mat *left, struct isl_mat *right,
 416         unsigned row, unsigned i, isl_int m)
 417 {
 418         isl_int_neg(m, m);
 419         isl_seq_combine(left->row[i]+row,
 420                         ctx->one, left->row[i]+row,
 421                         m, left->row[row]+row,
 422                         left->n_col-row);
 423         isl_seq_combine(right->row[i], ctx->one, right->row[i],
 424                         m, right->row[row], right->n_col);
 425 }
 426
 427 /* Compute inv(left)*right
 428  */
 429 struct isl_mat *isl_mat_inverse_product(struct isl_ctx *ctx,
 430         struct isl_mat *left, struct isl_mat *right)
 431 {
 432         int row;
 433         isl_int a, b;
 434
 435         if (!left || !right)
 436                 goto error;
 437
 438         isl_assert(ctx, left->n_row == left->n_col, goto error);
 439         isl_assert(ctx, left->n_row == right->n_row, goto error);
 440
 441         if (left->n_row == 0) {
 442                 isl_mat_free(ctx, left);
 443                 return right;
 444         }
 445
 446         left = isl_mat_cow(ctx, left);
 447         right = isl_mat_cow(ctx, right);
 448         if (!left || !right)
 449                 goto error;
 450
 451         isl_int_init(a);
 452         isl_int_init(b);
 453         for (row = 0; row < left->n_row; ++row) {
 454                 int pivot, first, i, off;
 455                 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
 456                 if (pivot < 0) {
 457                         isl_int_clear(a);
 458                         isl_int_clear(b);
 459                         isl_assert(ctx, pivot >= 0, goto error);
 460                 }
 461                 pivot += row;
 462                 if (pivot != row)
 463                         inv_exchange(ctx, left, right, pivot, row);
 464                 if (isl_int_is_neg(left->row[row][row]))
 465                         inv_oppose(ctx, left, right, row);
 466                 first = row+1;
 467                 while ((off = row_first_non_zero(left->row+first,
 468                                         left->n_row-first, row)) != -1) {
 469                         first += off;
 470                         isl_int_fdiv_q(a, left->row[first][row],
 471                                         left->row[row][row]);
 472                         inv_subtract(ctx, left, right, row, first, a);
 473                         if (!isl_int_is_zero(left->row[first][row]))
 474                                 inv_exchange(ctx, left, right, row, first);
 475                         else
 476                                 ++first;
 477                 }
 478                 for (i = 0; i < row; ++i) {
 479                         if (isl_int_is_zero(left->row[i][row]))
 480                                 continue;
 481                         isl_int_gcd(a, left->row[row][row], left->row[i][row]);
 482                         isl_int_divexact(b, left->row[i][row], a);
 483                         isl_int_divexact(a, left->row[row][row], a);
 484                         isl_int_neg(a, a);
 485                         isl_seq_combine(left->row[i]+row,
 486                                         a, left->row[i]+row,
 487                                         b, left->row[row]+row,
 488                                         left->n_col-row);
 489                         isl_seq_combine(right->row[i], a, right->row[i],
 490                                         b, right->row[row], right->n_col);
 491                 }
 492         }
 493         isl_int_clear(b);
 494
 495         isl_int_set(a, left->row[0][0]);
 496         for (row = 1; row < left->n_row; ++row)
 497                 isl_int_lcm(a, a, left->row[row][row]);
 498         if (isl_int_is_zero(a)){
 499                 isl_int_clear(a);
 500                 isl_assert(ctx, 0, goto error);
 501         }
 502         for (row = 0; row < left->n_row; ++row) {
 503                 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
 504                 if (isl_int_is_one(left->row[row][row]))
 505                         continue;
 506                 isl_seq_scale(right->row[row], right->row[row],
 507                                 left->row[row][row], right->n_col);
 508         }
 509         isl_int_clear(a);
 510
 511         isl_mat_free(ctx, left);
 512         return right;
 513 error:
 514         isl_mat_free(ctx, left);
 515         isl_mat_free(ctx, right);
 516         return NULL;
 517 }
 518
 519 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
 520         struct isl_mat *mat, unsigned i, unsigned j)
 521 {
 522         isl_int *t;
 523
 524         if (!mat)
 525                 return NULL;
 526         mat = isl_mat_cow(ctx, mat);
 527         if (!mat)
 528                 return NULL;
 529         t = mat->row[i];
 530         mat->row[i] = mat->row[j];
 531         mat->row[j] = t;
 532         return mat;
 533 }
 534
 535 struct isl_mat *isl_mat_product(struct isl_ctx *ctx,
 536         struct isl_mat *left, struct isl_mat *right)
 537 {
 538         int i, j, k;
 539         struct isl_mat *prod;
 540
 541         if (!left || !right)
 542                 goto error;
 543         isl_assert(ctx, left->n_col == right->n_row, goto error);
 544         prod = isl_mat_alloc(ctx, left->n_row, right->n_col);
 545         if (!prod)
 546                 goto error;
 547         if (left->n_col == 0) {
 548                 for (i = 0; i < prod->n_row; ++i)
 549                         isl_seq_clr(prod->row[i], prod->n_col);
 550                 return prod;
 551         }
 552         for (i = 0; i < prod->n_row; ++i) {
 553                 for (j = 0; j < prod->n_col; ++j) {
 554                         isl_int_mul(prod->row[i][j],
 555                                     left->row[i][0], right->row[0][j]);
 556                         for (k = 1; k < left->n_col; ++k)
 557                                 isl_int_addmul(prod->row[i][j],
 558                                             left->row[i][k], right->row[k][j]);
 559                 }
 560         }
 561         isl_mat_free(ctx, left);
 562         isl_mat_free(ctx, right);
 563         return prod;
 564 error:
 565         isl_mat_free(ctx, left);
 566         isl_mat_free(ctx, right);
 567         return NULL;
 568 }
 569
 570 struct isl_basic_set *isl_basic_set_preimage(struct isl_ctx *ctx,
 571         struct isl_basic_set *bset, struct isl_mat *mat)
 572 {
 573         struct isl_mat *t;
 574         int i;
 575
 576         if (!bset || !mat)
 577                 goto error;
 578
 579         bset = isl_basic_set_cow(ctx, bset);
 580         if (!bset)
 581                 goto error;
 582
 583         isl_assert(ctx, bset->nparam == 0, goto error);
 584         isl_assert(ctx, bset->n_div == 0, goto error);
 585         isl_assert(ctx, 1+bset->dim == mat->n_row, goto error);
 586
 587         t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
 588         t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
 589         if (!t)
 590                 goto error;
 591         for (i = 0; i < bset->n_eq; ++i)
 592                 isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
 593         isl_mat_free(ctx, t);
 594
 595         t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
 596         t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
 597         if (!t)
 598                 goto error;
 599         for (i = 0; i < bset->n_ineq; ++i)
 600                 isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
 601         isl_mat_free(ctx, t);
 602
 603         bset->dim -= mat->n_row - mat->n_col;
 604         bset = isl_basic_set_simplify(ctx, bset);
 605         bset = isl_basic_set_finalize(ctx, bset);
 606
 607         isl_mat_free(ctx, mat);
 608         return bset;
 609 error:
 610         isl_mat_free(ctx, mat);
 611         isl_basic_set_free(ctx, bset);
 612         return NULL;
 613 }
 614
 615 void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
 616                                 FILE *out, int indent)
 617 {
 618         int i, j;
 619
 620         if (mat->n_row == 0)
 621                 fprintf(out, "%*s[]\n", indent, "");
 622
 623         for (i = 0; i < mat->n_row; ++i) {
 624                 if (!i)
 625                         fprintf(out, "%*s[[", indent, "");
 626                 else
 627                         fprintf(out, "%*s[", indent+1, "");
 628                 for (j = 0; j < mat->n_col; ++j) {
 629                         if (j)
 630                             fprintf(out, ",");
 631                         isl_int_print(out, mat->row[i][j]);
 632                 }
 633                 if (i == mat->n_row-1)
 634                         fprintf(out, "]]\n");
 635                 else
 636                         fprintf(out, "]\n");
 637         }
 638 }
 639
 640 struct isl_mat *isl_mat_drop_col(struct isl_ctx *ctx, struct isl_mat *mat,
 641                                 unsigned col)
 642 {
 643         int r;
 644
 645         if (!mat)
 646                 return NULL;
 647
 648         if (col != mat->n_col-1) {
 649                 for (r = 0; r < mat->n_row; ++r)
 650                         isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+1,
 651                                         mat->n_col - col - 1);
 652         }
 653         mat->n_col--;
 654         return mat;
 655 }
 656
 657 struct isl_mat *isl_mat_drop_rows(struct isl_ctx *ctx, struct isl_mat *mat,
 658                                 unsigned row, unsigned n)
 659 {
 660         int r;
 661
 662         if (!mat)
 663                 return NULL;
 664
 665         for (r = row; r+n < mat->n_row; ++r)
 666                 mat->row[r] = mat->row[r+n];
 667
 668         mat->n_row -= n;
 669         return mat;
 670 }