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 void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
 520 {
 521         int i;
 522
 523         for (i = 0; i < mat->n_row; ++i)
 524                 isl_int_mul(mat->row[i][col], mat->row[i][col], m);
 525 }
 526
 527 isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
 528         isl_int m1, unsigned src1, isl_int m2, unsigned src2)
 529 {
 530         int i;
 531         isl_int tmp;
 532
 533         isl_int_init(tmp);
 534         for (i = 0; i < mat->n_row; ++i) {
 535                 isl_int_mul(tmp, m1, mat->row[i][src1]);
 536                 isl_int_addmul(tmp, m2, mat->row[i][src2]);
 537                 isl_int_set(mat->row[i][dst], tmp);
 538         }
 539         isl_int_clear(tmp);
 540 }
 541
 542 struct isl_mat *isl_mat_right_inverse(struct isl_ctx *ctx,
 543         struct isl_mat *mat)
 544 {
 545         struct isl_mat *inv;
 546         int row;
 547         isl_int a, b;
 548
 549         mat = isl_mat_cow(ctx, mat);
 550         if (!mat)
 551                 return NULL;
 552
 553         inv = isl_mat_identity(ctx, mat->n_col);
 554         inv = isl_mat_cow(ctx, inv);
 555         if (!inv)
 556                 goto error;
 557
 558         isl_int_init(a);
 559         isl_int_init(b);
 560         for (row = 0; row < mat->n_row; ++row) {
 561                 int pivot, first, i, off;
 562                 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
 563                 if (pivot < 0) {
 564                         isl_int_clear(a);
 565                         isl_int_clear(b);
 566                         goto error;
 567                 }
 568                 pivot += row;
 569                 if (pivot != row)
 570                         exchange(ctx, mat, &inv, NULL, row, pivot, row);
 571                 if (isl_int_is_neg(mat->row[row][row]))
 572                         oppose(ctx, mat, &inv, NULL, row, row);
 573                 first = row+1;
 574                 while ((off = isl_seq_first_non_zero(mat->row[row]+first,
 575                                                     mat->n_col-first)) != -1) {
 576                         first += off;
 577                         isl_int_fdiv_q(a, mat->row[row][first],
 578                                                     mat->row[row][row]);
 579                         subtract(ctx, mat, &inv, NULL, row, row, first, a);
 580                         if (!isl_int_is_zero(mat->row[row][first]))
 581                                 exchange(ctx, mat, &inv, NULL, row, row, first);
 582                         else
 583                                 ++first;
 584                 }
 585                 for (i = 0; i < row; ++i) {
 586                         if (isl_int_is_zero(mat->row[row][i]))
 587                                 continue;
 588                         isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
 589                         isl_int_divexact(b, mat->row[row][i], a);
 590                         isl_int_divexact(a, mat->row[row][row], a);
 591                         isl_int_neg(a, a);
 592                         isl_mat_col_combine(mat, i, a, i, b, row);
 593                         isl_mat_col_combine(inv, i, a, i, b, row);
 594                 }
 595         }
 596         isl_int_clear(b);
 597
 598         isl_int_set(a, mat->row[0][0]);
 599         for (row = 1; row < mat->n_row; ++row)
 600                 isl_int_lcm(a, a, mat->row[row][row]);
 601         if (isl_int_is_zero(a)){
 602                 isl_int_clear(a);
 603                 goto error;
 604         }
 605         for (row = 0; row < mat->n_row; ++row) {
 606                 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
 607                 if (isl_int_is_one(mat->row[row][row]))
 608                         continue;
 609                 isl_mat_col_scale(inv, row, mat->row[row][row]);
 610         }
 611         isl_int_clear(a);
 612
 613         isl_mat_free(ctx, mat);
 614
 615         return inv;
 616 error:
 617         isl_mat_free(ctx, mat);
 618         return NULL;
 619 }
 620
 621 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
 622         struct isl_mat *mat, unsigned i, unsigned j)
 623 {
 624         isl_int *t;
 625
 626         if (!mat)
 627                 return NULL;
 628         mat = isl_mat_cow(ctx, mat);
 629         if (!mat)
 630                 return NULL;
 631         t = mat->row[i];
 632         mat->row[i] = mat->row[j];
 633         mat->row[j] = t;
 634         return mat;
 635 }
 636
 637 struct isl_mat *isl_mat_product(struct isl_ctx *ctx,
 638         struct isl_mat *left, struct isl_mat *right)
 639 {
 640         int i, j, k;
 641         struct isl_mat *prod;
 642
 643         if (!left || !right)
 644                 goto error;
 645         isl_assert(ctx, left->n_col == right->n_row, goto error);
 646         prod = isl_mat_alloc(ctx, left->n_row, right->n_col);
 647         if (!prod)
 648                 goto error;
 649         if (left->n_col == 0) {
 650                 for (i = 0; i < prod->n_row; ++i)
 651                         isl_seq_clr(prod->row[i], prod->n_col);
 652                 return prod;
 653         }
 654         for (i = 0; i < prod->n_row; ++i) {
 655                 for (j = 0; j < prod->n_col; ++j) {
 656                         isl_int_mul(prod->row[i][j],
 657                                     left->row[i][0], right->row[0][j]);
 658                         for (k = 1; k < left->n_col; ++k)
 659                                 isl_int_addmul(prod->row[i][j],
 660                                             left->row[i][k], right->row[k][j]);
 661                 }
 662         }
 663         isl_mat_free(ctx, left);
 664         isl_mat_free(ctx, right);
 665         return prod;
 666 error:
 667         isl_mat_free(ctx, left);
 668         isl_mat_free(ctx, right);
 669         return NULL;
 670 }
 671
 672 struct isl_basic_set *isl_basic_set_preimage(struct isl_ctx *ctx,
 673         struct isl_basic_set *bset, struct isl_mat *mat)
 674 {
 675         struct isl_mat *t;
 676         int i;
 677
 678         if (!bset || !mat)
 679                 goto error;
 680
 681         bset = isl_basic_set_cow(ctx, bset);
 682         if (!bset)
 683                 goto error;
 684
 685         isl_assert(ctx, bset->nparam == 0, goto error);
 686         isl_assert(ctx, bset->n_div == 0, goto error);
 687         isl_assert(ctx, 1+bset->dim == mat->n_row, goto error);
 688
 689         if (mat->n_col > mat->n_row)
 690                 bset = isl_basic_set_extend(ctx, bset, 0, mat->n_col-1, 0,
 691                                                 0, 0);
 692         else {
 693                 bset->dim -= mat->n_row - mat->n_col;
 694                 bset->extra += mat->n_row - mat->n_col;
 695         }
 696
 697         t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
 698         t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
 699         if (!t)
 700                 goto error;
 701         for (i = 0; i < bset->n_eq; ++i) {
 702                 isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
 703                 isl_seq_clr(bset->eq[i]+t->n_col, bset->extra);
 704         }
 705         isl_mat_free(ctx, t);
 706
 707         t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
 708         t = isl_mat_product(ctx, t, mat);
 709         if (!t)
 710                 goto error2;
 711         for (i = 0; i < bset->n_ineq; ++i) {
 712                 isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
 713                 isl_seq_clr(bset->ineq[i]+t->n_col, bset->extra);
 714         }
 715         isl_mat_free(ctx, t);
 716
 717         bset = isl_basic_set_simplify(ctx, bset);
 718         bset = isl_basic_set_finalize(ctx, bset);
 719
 720         return bset;
 721 error:
 722         isl_mat_free(ctx, mat);
 723 error2:
 724         isl_basic_set_free(ctx, bset);
 725         return NULL;
 726 }
 727
 728 struct isl_set *isl_set_preimage(struct isl_ctx *ctx,
 729         struct isl_set *set, struct isl_mat *mat)
 730 {
 731         int i;
 732
 733         set = isl_set_cow(ctx, set);
 734         if (!set)
 735                 return NULL;
 736
 737         for (i = 0; i < set->n; ++i) {
 738                 set->p[i] = isl_basic_set_preimage(ctx, set->p[i],
 739                                                     isl_mat_copy(ctx, mat));
 740                 if (!set->p[i])
 741                         goto error;
 742         }
 743         set->dim += mat->n_col;
 744         set->dim -= mat->n_row;
 745         isl_mat_free(ctx, mat);
 746         return set;
 747 error:
 748         isl_set_free(ctx, set);
 749         isl_mat_free(ctx, mat);
 750         return NULL;
 751 }
 752
 753 void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
 754                                 FILE *out, int indent)
 755 {
 756         int i, j;
 757
 758         if (!mat) {
 759                 fprintf(out, "null mat\n");
 760                 return;
 761         }
 762
 763         if (mat->n_row == 0)
 764                 fprintf(out, "%*s[]\n", indent, "");
 765
 766         for (i = 0; i < mat->n_row; ++i) {
 767                 if (!i)
 768                         fprintf(out, "%*s[[", indent, "");
 769                 else
 770                         fprintf(out, "%*s[", indent+1, "");
 771                 for (j = 0; j < mat->n_col; ++j) {
 772                         if (j)
 773                             fprintf(out, ",");
 774                         isl_int_print(out, mat->row[i][j]);
 775                 }
 776                 if (i == mat->n_row-1)
 777                         fprintf(out, "]]\n");
 778                 else
 779                         fprintf(out, "]\n");
 780         }
 781 }
 782
 783 struct isl_mat *isl_mat_drop_col(struct isl_ctx *ctx, struct isl_mat *mat,
 784                                 unsigned col)
 785 {
 786         int r;
 787
 788         if (!mat)
 789                 return NULL;
 790
 791         if (col != mat->n_col-1) {
 792                 for (r = 0; r < mat->n_row; ++r)
 793                         isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+1,
 794                                         mat->n_col - col - 1);
 795         }
 796         mat->n_col--;
 797         return mat;
 798 }
 799
 800 struct isl_mat *isl_mat_drop_rows(struct isl_ctx *ctx, struct isl_mat *mat,
 801                                 unsigned row, unsigned n)
 802 {
 803         int r;
 804
 805         if (!mat)
 806                 return NULL;
 807
 808         for (r = row; r+n < mat->n_row; ++r)
 809                 mat->row[r] = mat->row[r+n];
 810
 811         mat->n_row -= n;
 812         return mat;
 813 }