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_mat *M, struct isl_mat **U,
 242         struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
 243 {
 244         int r;
 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 (if neg == 0) or non-positive (if neg == 1)
 279  * and stricly smaller (in absolute value) than the entries in the echelon
 280  * column.
 281  * If U or Q are NULL, then these matrices are not computed.
 282  */
 283 struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
 284         struct isl_mat *M, int neg, struct isl_mat **U, struct isl_mat **Q)
 285 {
 286         isl_int c;
 287         int row, col;
 288
 289         if (U)
 290                 *U = NULL;
 291         if (Q)
 292                 *Q = NULL;
 293         if (!M)
 294                 goto error;
 295         M = isl_mat_cow(ctx, M);
 296         if (!M)
 297                 goto error;
 298         if (U) {
 299                 *U = isl_mat_identity(ctx, M->n_col);
 300                 if (!*U)
 301                         goto error;
 302         }
 303         if (Q) {
 304                 *Q = isl_mat_identity(ctx, M->n_col);
 305                 if (!*Q)
 306                         goto error;
 307         }
 308
 309         col = 0;
 310         isl_int_init(c);
 311         for (row = 0; row < M->n_row; ++row) {
 312                 int first, i, off;
 313                 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
 314                 if (first == -1)
 315                         continue;
 316                 first += col;
 317                 if (first != col)
 318                         exchange(ctx, M, U, Q, row, first, col);
 319                 if (isl_int_is_neg(M->row[row][col]))
 320                         oppose(ctx, M, U, Q, row, col);
 321                 first = col+1;
 322                 while ((off = isl_seq_first_non_zero(M->row[row]+first,
 323                                                        M->n_col-first)) != -1) {
 324                         first += off;
 325                         isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
 326                         subtract(M, U, Q, row, col, first, c);
 327                         if (!isl_int_is_zero(M->row[row][first]))
 328                                 exchange(ctx, M, U, Q, row, first, col);
 329                         else
 330                                 ++first;
 331                 }
 332                 for (i = 0; i < col; ++i) {
 333                         if (isl_int_is_zero(M->row[row][i]))
 334                                 continue;
 335                         if (neg)
 336                                 isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
 337                         else
 338                                 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
 339                         if (isl_int_is_zero(c))
 340                                 continue;
 341                         subtract(M, U, Q, row, col, i, c);
 342                 }
 343                 ++col;
 344         }
 345         isl_int_clear(c);
 346
 347         return M;
 348 error:
 349         if (Q) {
 350                 isl_mat_free(ctx, *Q);
 351                 *Q = NULL;
 352         }
 353         if (U) {
 354                 isl_mat_free(ctx, *U);
 355                 *U = NULL;
 356         }
 357         return NULL;
 358 }
 359
 360 struct isl_mat *isl_mat_lin_to_aff(struct isl_ctx *ctx, struct isl_mat *mat)
 361 {
 362         int i;
 363         struct isl_mat *mat2;
 364
 365         if (!mat)
 366                 return NULL;
 367         mat2 = isl_mat_alloc(ctx, 1+mat->n_row, 1+mat->n_col);
 368         if (!mat2)
 369                 return NULL;
 370         isl_int_set_si(mat2->row[0][0], 1);
 371         isl_seq_clr(mat2->row[0]+1, mat->n_col);
 372         for (i = 0; i < mat->n_row; ++i) {
 373                 isl_int_set_si(mat2->row[1+i][0], 0);
 374                 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
 375         }
 376         isl_mat_free(ctx, mat);
 377         return mat2;
 378 }
 379
 380 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
 381 {
 382         int i;
 383
 384         for (i = 0; i < n_row; ++i)
 385                 if (!isl_int_is_zero(row[i][col]))
 386                         return i;
 387         return -1;
 388 }
 389
 390 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
 391 {
 392         int i, min = row_first_non_zero(row, n_row, col);
 393         if (min < 0)
 394                 return -1;
 395         for (i = min + 1; i < n_row; ++i) {
 396                 if (isl_int_is_zero(row[i][col]))
 397                         continue;
 398                 if (isl_int_abs_lt(row[i][col], row[min][col]))
 399                         min = i;
 400         }
 401         return min;
 402 }
 403
 404 static void inv_exchange(struct isl_ctx *ctx,
 405         struct isl_mat *left, struct isl_mat *right,
 406         unsigned i, unsigned j)
 407 {
 408         left = isl_mat_swap_rows(ctx, left, i, j);
 409         right = isl_mat_swap_rows(ctx, right, i, j);
 410 }
 411
 412 static void inv_oppose(struct isl_ctx *ctx,
 413         struct isl_mat *left, struct isl_mat *right, unsigned row)
 414 {
 415         isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
 416         isl_seq_neg(right->row[row], right->row[row], right->n_col);
 417 }
 418
 419 static void inv_subtract(struct isl_ctx *ctx,
 420         struct isl_mat *left, struct isl_mat *right,
 421         unsigned row, unsigned i, isl_int m)
 422 {
 423         isl_int_neg(m, m);
 424         isl_seq_combine(left->row[i]+row,
 425                         ctx->one, left->row[i]+row,
 426                         m, left->row[row]+row,
 427                         left->n_col-row);
 428         isl_seq_combine(right->row[i], ctx->one, right->row[i],
 429                         m, right->row[row], right->n_col);
 430 }
 431
 432 /* Compute inv(left)*right
 433  */
 434 struct isl_mat *isl_mat_inverse_product(struct isl_ctx *ctx,
 435         struct isl_mat *left, struct isl_mat *right)
 436 {
 437         int row;
 438         isl_int a, b;
 439
 440         if (!left || !right)
 441                 goto error;
 442
 443         isl_assert(ctx, left->n_row == left->n_col, goto error);
 444         isl_assert(ctx, left->n_row == right->n_row, goto error);
 445
 446         if (left->n_row == 0) {
 447                 isl_mat_free(ctx, left);
 448                 return right;
 449         }
 450
 451         left = isl_mat_cow(ctx, left);
 452         right = isl_mat_cow(ctx, right);
 453         if (!left || !right)
 454                 goto error;
 455
 456         isl_int_init(a);
 457         isl_int_init(b);
 458         for (row = 0; row < left->n_row; ++row) {
 459                 int pivot, first, i, off;
 460                 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
 461                 if (pivot < 0) {
 462                         isl_int_clear(a);
 463                         isl_int_clear(b);
 464                         isl_assert(ctx, pivot >= 0, goto error);
 465                 }
 466                 pivot += row;
 467                 if (pivot != row)
 468                         inv_exchange(ctx, left, right, pivot, row);
 469                 if (isl_int_is_neg(left->row[row][row]))
 470                         inv_oppose(ctx, left, right, row);
 471                 first = row+1;
 472                 while ((off = row_first_non_zero(left->row+first,
 473                                         left->n_row-first, row)) != -1) {
 474                         first += off;
 475                         isl_int_fdiv_q(a, left->row[first][row],
 476                                         left->row[row][row]);
 477                         inv_subtract(ctx, left, right, row, first, a);
 478                         if (!isl_int_is_zero(left->row[first][row]))
 479                                 inv_exchange(ctx, left, right, row, first);
 480                         else
 481                                 ++first;
 482                 }
 483                 for (i = 0; i < row; ++i) {
 484                         if (isl_int_is_zero(left->row[i][row]))
 485                                 continue;
 486                         isl_int_gcd(a, left->row[row][row], left->row[i][row]);
 487                         isl_int_divexact(b, left->row[i][row], a);
 488                         isl_int_divexact(a, left->row[row][row], a);
 489                         isl_int_neg(a, a);
 490                         isl_seq_combine(left->row[i]+row,
 491                                         a, left->row[i]+row,
 492                                         b, left->row[row]+row,
 493                                         left->n_col-row);
 494                         isl_seq_combine(right->row[i], a, right->row[i],
 495                                         b, right->row[row], right->n_col);
 496                 }
 497         }
 498         isl_int_clear(b);
 499
 500         isl_int_set(a, left->row[0][0]);
 501         for (row = 1; row < left->n_row; ++row)
 502                 isl_int_lcm(a, a, left->row[row][row]);
 503         if (isl_int_is_zero(a)){
 504                 isl_int_clear(a);
 505                 isl_assert(ctx, 0, goto error);
 506         }
 507         for (row = 0; row < left->n_row; ++row) {
 508                 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
 509                 if (isl_int_is_one(left->row[row][row]))
 510                         continue;
 511                 isl_seq_scale(right->row[row], right->row[row],
 512                                 left->row[row][row], right->n_col);
 513         }
 514         isl_int_clear(a);
 515
 516         isl_mat_free(ctx, left);
 517         return right;
 518 error:
 519         isl_mat_free(ctx, left);
 520         isl_mat_free(ctx, right);
 521         return NULL;
 522 }
 523
 524 void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
 525 {
 526         int i;
 527
 528         for (i = 0; i < mat->n_row; ++i)
 529                 isl_int_mul(mat->row[i][col], mat->row[i][col], m);
 530 }
 531
 532 void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
 533         isl_int m1, unsigned src1, isl_int m2, unsigned src2)
 534 {
 535         int i;
 536         isl_int tmp;
 537
 538         isl_int_init(tmp);
 539         for (i = 0; i < mat->n_row; ++i) {
 540                 isl_int_mul(tmp, m1, mat->row[i][src1]);
 541                 isl_int_addmul(tmp, m2, mat->row[i][src2]);
 542                 isl_int_set(mat->row[i][dst], tmp);
 543         }
 544         isl_int_clear(tmp);
 545 }
 546
 547 struct isl_mat *isl_mat_right_inverse(struct isl_ctx *ctx,
 548         struct isl_mat *mat)
 549 {
 550         struct isl_mat *inv;
 551         int row;
 552         isl_int a, b;
 553
 554         mat = isl_mat_cow(ctx, mat);
 555         if (!mat)
 556                 return NULL;
 557
 558         inv = isl_mat_identity(ctx, mat->n_col);
 559         inv = isl_mat_cow(ctx, inv);
 560         if (!inv)
 561                 goto error;
 562
 563         isl_int_init(a);
 564         isl_int_init(b);
 565         for (row = 0; row < mat->n_row; ++row) {
 566                 int pivot, first, i, off;
 567                 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
 568                 if (pivot < 0) {
 569                         isl_int_clear(a);
 570                         isl_int_clear(b);
 571                         goto error;
 572                 }
 573                 pivot += row;
 574                 if (pivot != row)
 575                         exchange(ctx, mat, &inv, NULL, row, pivot, row);
 576                 if (isl_int_is_neg(mat->row[row][row]))
 577                         oppose(ctx, mat, &inv, NULL, row, row);
 578                 first = row+1;
 579                 while ((off = isl_seq_first_non_zero(mat->row[row]+first,
 580                                                     mat->n_col-first)) != -1) {
 581                         first += off;
 582                         isl_int_fdiv_q(a, mat->row[row][first],
 583                                                     mat->row[row][row]);
 584                         subtract(mat, &inv, NULL, row, row, first, a);
 585                         if (!isl_int_is_zero(mat->row[row][first]))
 586                                 exchange(ctx, mat, &inv, NULL, row, row, first);
 587                         else
 588                                 ++first;
 589                 }
 590                 for (i = 0; i < row; ++i) {
 591                         if (isl_int_is_zero(mat->row[row][i]))
 592                                 continue;
 593                         isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
 594                         isl_int_divexact(b, mat->row[row][i], a);
 595                         isl_int_divexact(a, mat->row[row][row], a);
 596                         isl_int_neg(a, a);
 597                         isl_mat_col_combine(mat, i, a, i, b, row);
 598                         isl_mat_col_combine(inv, i, a, i, b, row);
 599                 }
 600         }
 601         isl_int_clear(b);
 602
 603         isl_int_set(a, mat->row[0][0]);
 604         for (row = 1; row < mat->n_row; ++row)
 605                 isl_int_lcm(a, a, mat->row[row][row]);
 606         if (isl_int_is_zero(a)){
 607                 isl_int_clear(a);
 608                 goto error;
 609         }
 610         for (row = 0; row < mat->n_row; ++row) {
 611                 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
 612                 if (isl_int_is_one(mat->row[row][row]))
 613                         continue;
 614                 isl_mat_col_scale(inv, row, mat->row[row][row]);
 615         }
 616         isl_int_clear(a);
 617
 618         isl_mat_free(ctx, mat);
 619
 620         return inv;
 621 error:
 622         isl_mat_free(ctx, mat);
 623         return NULL;
 624 }
 625
 626 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
 627         struct isl_mat *mat, unsigned i, unsigned j)
 628 {
 629         isl_int *t;
 630
 631         if (!mat)
 632                 return NULL;
 633         mat = isl_mat_cow(ctx, mat);
 634         if (!mat)
 635                 return NULL;
 636         t = mat->row[i];
 637         mat->row[i] = mat->row[j];
 638         mat->row[j] = t;
 639         return mat;
 640 }
 641
 642 struct isl_mat *isl_mat_product(struct isl_ctx *ctx,
 643         struct isl_mat *left, struct isl_mat *right)
 644 {
 645         int i, j, k;
 646         struct isl_mat *prod;
 647
 648         if (!left || !right)
 649                 goto error;
 650         isl_assert(ctx, left->n_col == right->n_row, goto error);
 651         prod = isl_mat_alloc(ctx, left->n_row, right->n_col);
 652         if (!prod)
 653                 goto error;
 654         if (left->n_col == 0) {
 655                 for (i = 0; i < prod->n_row; ++i)
 656                         isl_seq_clr(prod->row[i], prod->n_col);
 657                 return prod;
 658         }
 659         for (i = 0; i < prod->n_row; ++i) {
 660                 for (j = 0; j < prod->n_col; ++j) {
 661                         isl_int_mul(prod->row[i][j],
 662                                     left->row[i][0], right->row[0][j]);
 663                         for (k = 1; k < left->n_col; ++k)
 664                                 isl_int_addmul(prod->row[i][j],
 665                                             left->row[i][k], right->row[k][j]);
 666                 }
 667         }
 668         isl_mat_free(ctx, left);
 669         isl_mat_free(ctx, right);
 670         return prod;
 671 error:
 672         isl_mat_free(ctx, left);
 673         isl_mat_free(ctx, right);
 674         return NULL;
 675 }
 676
 677 struct isl_basic_set *isl_basic_set_preimage(struct isl_ctx *ctx,
 678         struct isl_basic_set *bset, struct isl_mat *mat)
 679 {
 680         struct isl_mat *t;
 681         int i;
 682
 683         if (!bset || !mat)
 684                 goto error;
 685
 686         bset = isl_basic_set_cow(bset);
 687         if (!bset)
 688                 goto error;
 689
 690         isl_assert(ctx, bset->nparam == 0, goto error);
 691         isl_assert(ctx, bset->n_div == 0, goto error);
 692         isl_assert(ctx, 1+bset->dim == mat->n_row, goto error);
 693
 694         if (mat->n_col > mat->n_row)
 695                 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0,
 696                                                 0, 0);
 697         else {
 698                 bset->dim -= mat->n_row - mat->n_col;
 699                 bset->extra += mat->n_row - mat->n_col;
 700         }
 701
 702         t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
 703         t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
 704         if (!t)
 705                 goto error;
 706         for (i = 0; i < bset->n_eq; ++i) {
 707                 isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
 708                 isl_seq_clr(bset->eq[i]+t->n_col, bset->extra);
 709         }
 710         isl_mat_free(ctx, t);
 711
 712         t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
 713         t = isl_mat_product(ctx, t, mat);
 714         if (!t)
 715                 goto error2;
 716         for (i = 0; i < bset->n_ineq; ++i) {
 717                 isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
 718                 isl_seq_clr(bset->ineq[i]+t->n_col, bset->extra);
 719         }
 720         isl_mat_free(ctx, t);
 721
 722         bset = isl_basic_set_simplify(bset);
 723         bset = isl_basic_set_finalize(bset);
 724
 725         return bset;
 726 error:
 727         isl_mat_free(ctx, mat);
 728 error2:
 729         isl_basic_set_free(bset);
 730         return NULL;
 731 }
 732
 733 struct isl_set *isl_set_preimage(struct isl_ctx *ctx,
 734         struct isl_set *set, struct isl_mat *mat)
 735 {
 736         int i;
 737
 738         set = isl_set_cow(set);
 739         if (!set)
 740                 return NULL;
 741
 742         ctx = set->ctx;
 743         for (i = 0; i < set->n; ++i) {
 744                 set->p[i] = isl_basic_set_preimage(ctx, set->p[i],
 745                                                     isl_mat_copy(ctx, mat));
 746                 if (!set->p[i])
 747                         goto error;
 748         }
 749         set->dim += mat->n_col;
 750         set->dim -= mat->n_row;
 751         isl_mat_free(ctx, mat);
 752         return set;
 753 error:
 754         isl_set_free(set);
 755         isl_mat_free(ctx, mat);
 756         return NULL;
 757 }
 758
 759 void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
 760                                 FILE *out, int indent)
 761 {
 762         int i, j;
 763
 764         if (!mat) {
 765                 fprintf(out, "null mat\n");
 766                 return;
 767         }
 768
 769         if (mat->n_row == 0)
 770                 fprintf(out, "%*s[]\n", indent, "");
 771
 772         for (i = 0; i < mat->n_row; ++i) {
 773                 if (!i)
 774                         fprintf(out, "%*s[[", indent, "");
 775                 else
 776                         fprintf(out, "%*s[", indent+1, "");
 777                 for (j = 0; j < mat->n_col; ++j) {
 778                         if (j)
 779                             fprintf(out, ",");
 780                         isl_int_print(out, mat->row[i][j]);
 781                 }
 782                 if (i == mat->n_row-1)
 783                         fprintf(out, "]]\n");
 784                 else
 785                         fprintf(out, "]\n");
 786         }
 787 }
 788
 789 struct isl_mat *isl_mat_drop_cols(struct isl_ctx *ctx, struct isl_mat *mat,
 790                                 unsigned col, unsigned n)
 791 {
 792         int r;
 793
 794         if (!mat)
 795                 return NULL;
 796
 797         if (col != mat->n_col-n) {
 798                 for (r = 0; r < mat->n_row; ++r)
 799                         isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
 800                                         mat->n_col - col - n);
 801         }
 802         mat->n_col -= n;
 803         return mat;
 804 }
 805
 806 struct isl_mat *isl_mat_drop_rows(struct isl_ctx *ctx, struct isl_mat *mat,
 807                                 unsigned row, unsigned n)
 808 {
 809         int r;
 810
 811         if (!mat)
 812                 return NULL;
 813
 814         for (r = row; r+n < mat->n_row; ++r)
 815                 mat->row[r] = mat->row[r+n];
 816
 817         mat->n_row -= n;
 818         return mat;
 819 }