isl_mat.c

   1 #include "isl_dim.h"
   2 #include "isl_seq.h"
   3 #include "isl_mat.h"
   4 #include "isl_map_private.h"
   5
   6 struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
   7         unsigned n_row, unsigned n_col)
   8 {
   9         int i;
  10         struct isl_mat *mat;
  11
  12         mat = isl_alloc_type(ctx, struct isl_mat);
  13         if (!mat)
  14                 return NULL;
  15
  16         mat->row = NULL;
  17         mat->block = isl_blk_alloc(ctx, n_row * n_col);
  18         if (isl_blk_is_error(mat->block))
  19                 goto error;
  20         mat->row = isl_alloc_array(ctx, isl_int *, n_row);
  21         if (!mat->row)
  22                 goto error;
  23
  24         for (i = 0; i < n_row; ++i)
  25                 mat->row[i] = mat->block.data + i * n_col;
  26
  27         mat->ref = 1;
  28         mat->n_row = n_row;
  29         mat->n_col = n_col;
  30         mat->flags = 0;
  31
  32         return mat;
  33 error:
  34         isl_blk_free(ctx, mat->block);
  35         free(mat);
  36         return NULL;
  37 }
  38
  39 struct isl_mat *isl_mat_sub_alloc(struct isl_ctx *ctx, isl_int **row,
  40         unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
  41 {
  42         int i;
  43         struct isl_mat *mat;
  44
  45         mat = isl_alloc_type(ctx, struct isl_mat);
  46         if (!mat)
  47                 return NULL;
  48         mat->row = isl_alloc_array(ctx, isl_int *, n_row);
  49         if (!mat->row)
  50                 goto error;
  51         for (i = 0; i < n_row; ++i)
  52                 mat->row[i] = row[first_row+i] + first_col;
  53         mat->ref = 1;
  54         mat->n_row = n_row;
  55         mat->n_col = n_col;
  56         mat->block = isl_blk_empty();
  57         mat->flags = ISL_MAT_BORROWED;
  58         return mat;
  59 error:
  60         free(mat);
  61         return NULL;
  62 }
  63
  64 void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
  65         unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
  66 {
  67         int i;
  68
  69         for (i = 0; i < n_row; ++i)
  70                 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
  71 }
  72
  73 void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
  74         unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
  75 {
  76         int i;
  77
  78         for (i = 0; i < n_row; ++i)
  79                 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
  80 }
  81
  82 struct isl_mat *isl_mat_copy(struct isl_ctx *ctx, struct isl_mat *mat)
  83 {
  84         if (!mat)
  85                 return NULL;
  86
  87         mat->ref++;
  88         return mat;
  89 }
  90
  91 struct isl_mat *isl_mat_dup(struct isl_ctx *ctx, struct isl_mat *mat)
  92 {
  93         int i;
  94         struct isl_mat *mat2;
  95
  96         if (!mat)
  97                 return NULL;
  98         mat2 = isl_mat_alloc(ctx, mat->n_row, mat->n_col);
  99         if (!mat2)
 100                 return NULL;
 101         for (i = 0; i < mat->n_row; ++i)
 102                 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
 103         return mat2;
 104 }
 105
 106 struct isl_mat *isl_mat_cow(struct isl_ctx *ctx, struct isl_mat *mat)
 107 {
 108         struct isl_mat *mat2;
 109         if (!mat)
 110                 return NULL;
 111
 112         if (mat->ref == 1 && !F_ISSET(mat, ISL_MAT_BORROWED))
 113                 return mat;
 114
 115         mat2 = isl_mat_dup(ctx, mat);
 116         isl_mat_free(ctx, mat);
 117         return mat2;
 118 }
 119
 120 void isl_mat_free(struct isl_ctx *ctx, struct isl_mat *mat)
 121 {
 122         if (!mat)
 123                 return;
 124
 125         if (--mat->ref > 0)
 126                 return;
 127
 128         if (!F_ISSET(mat, ISL_MAT_BORROWED))
 129                 isl_blk_free(ctx, mat->block);
 130         free(mat->row);
 131         free(mat);
 132 }
 133
 134 struct isl_mat *isl_mat_identity(struct isl_ctx *ctx, unsigned n_row)
 135 {
 136         int i;
 137         struct isl_mat *mat;
 138
 139         mat = isl_mat_alloc(ctx, n_row, n_row);
 140         if (!mat)
 141                 return NULL;
 142         for (i = 0; i < n_row; ++i) {
 143                 isl_seq_clr(mat->row[i], i);
 144                 isl_int_set_si(mat->row[i][i], 1);
 145                 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
 146         }
 147
 148         return mat;
 149 }
 150
 151 struct isl_vec *isl_mat_vec_product(struct isl_ctx *ctx,
 152         struct isl_mat *mat, struct isl_vec *vec)
 153 {
 154         int i;
 155         struct isl_vec *prod;
 156
 157         if (!mat || !vec)
 158                 goto error;
 159
 160         isl_assert(ctx, mat->n_col == vec->size, goto error);
 161
 162         prod = isl_vec_alloc(ctx, mat->n_row);
 163         if (!prod)
 164                 goto error;
 165
 166         for (i = 0; i < prod->size; ++i)
 167                 isl_seq_inner_product(mat->row[i], vec->block.data, vec->size,
 168                                         &prod->block.data[i]);
 169         isl_mat_free(ctx, mat);
 170         isl_vec_free(ctx, vec);
 171         return prod;
 172 error:
 173         isl_mat_free(ctx, mat);
 174         isl_vec_free(ctx, vec);
 175         return NULL;
 176 }
 177
 178 struct isl_mat *isl_mat_aff_direct_sum(struct isl_ctx *ctx,
 179         struct isl_mat *left, struct isl_mat *right)
 180 {
 181         int i;
 182         struct isl_mat *sum;
 183
 184         if (!left || !right)
 185                 goto error;
 186
 187         isl_assert(ctx, left->n_row == right->n_row, goto error);
 188         isl_assert(ctx, left->n_row >= 1, goto error);
 189         isl_assert(ctx, left->n_col >= 1, goto error);
 190         isl_assert(ctx, right->n_col >= 1, goto error);
 191         isl_assert(ctx,
 192             isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
 193             goto error);
 194         isl_assert(ctx,
 195             isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
 196             goto error);
 197
 198         sum = isl_mat_alloc(ctx, left->n_row, left->n_col + right->n_col - 1);
 199         if (!sum)
 200                 goto error;
 201         isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
 202         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
 203         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
 204
 205         isl_seq_clr(sum->row[0]+1, sum->n_col-1);
 206         for (i = 1; i < sum->n_row; ++i) {
 207                 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
 208                 isl_int_addmul(sum->row[i][0],
 209                                 right->row[0][0], right->row[i][0]);
 210                 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
 211                                 left->n_col-1);
 212                 isl_seq_scale(sum->row[i]+left->n_col,
 213                                 right->row[i]+1, right->row[0][0],
 214                                 right->n_col-1);
 215         }
 216
 217         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
 218         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
 219         isl_mat_free(ctx, left);
 220         isl_mat_free(ctx, right);
 221         return sum;
 222 error:
 223         isl_mat_free(ctx, left);
 224         isl_mat_free(ctx, right);
 225         return NULL;
 226 }
 227
 228 static void exchange(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
 229         struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
 230 {
 231         int r;
 232         for (r = row; r < M->n_row; ++r)
 233                 isl_int_swap(M->row[r][i], M->row[r][j]);
 234         if (U) {
 235                 for (r = 0; r < (*U)->n_row; ++r)
 236                         isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
 237         }
 238         if (Q)
 239                 isl_mat_swap_rows(ctx, *Q, i, j);
 240 }
 241
 242 static void subtract(struct isl_mat *M, struct isl_mat **U,
 243         struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
 244 {
 245         int r;
 246         for (r = row; r < M->n_row; ++r)
 247                 isl_int_submul(M->row[r][j], m, M->row[r][i]);
 248         if (U) {
 249                 for (r = 0; r < (*U)->n_row; ++r)
 250                         isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
 251         }
 252         if (Q) {
 253                 for (r = 0; r < (*Q)->n_col; ++r)
 254                         isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
 255         }
 256 }
 257
 258 static void oppose(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
 259         struct isl_mat **Q, unsigned row, unsigned col)
 260 {
 261         int r;
 262         for (r = row; r < M->n_row; ++r)
 263                 isl_int_neg(M->row[r][col], M->row[r][col]);
 264         if (U) {
 265                 for (r = 0; r < (*U)->n_row; ++r)
 266                         isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
 267         }
 268         if (Q)
 269                 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
 270 }
 271
 272 /* Given matrix M, compute
 273  *
 274  *              M U = H
 275  *              M   = H Q
 276  *
 277  * with U and Q unimodular matrices and H a matrix in column echelon form
 278  * such that on each echelon row the entries in the non-echelon column
 279  * are non-negative (if neg == 0) or non-positive (if neg == 1)
 280  * and stricly smaller (in absolute value) than the entries in the echelon
 281  * column.
 282  * If U or Q are NULL, then these matrices are not computed.
 283  */
 284 struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
 285         struct isl_mat *M, int neg, struct isl_mat **U, struct isl_mat **Q)
 286 {
 287         isl_int c;
 288         int row, col;
 289
 290         if (U)
 291                 *U = NULL;
 292         if (Q)
 293                 *Q = NULL;
 294         if (!M)
 295                 goto error;
 296         M = isl_mat_cow(ctx, M);
 297         if (!M)
 298                 goto error;
 299         if (U) {
 300                 *U = isl_mat_identity(ctx, M->n_col);
 301                 if (!*U)
 302                         goto error;
 303         }
 304         if (Q) {
 305                 *Q = isl_mat_identity(ctx, M->n_col);
 306                 if (!*Q)
 307                         goto error;
 308         }
 309
 310         col = 0;
 311         isl_int_init(c);
 312         for (row = 0; row < M->n_row; ++row) {
 313                 int first, i, off;
 314                 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
 315                 if (first == -1)
 316                         continue;
 317                 first += col;
 318                 if (first != col)
 319                         exchange(ctx, M, U, Q, row, first, col);
 320                 if (isl_int_is_neg(M->row[row][col]))
 321                         oppose(ctx, M, U, Q, row, col);
 322                 first = col+1;
 323                 while ((off = isl_seq_first_non_zero(M->row[row]+first,
 324                                                        M->n_col-first)) != -1) {
 325                         first += off;
 326                         isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
 327                         subtract(M, U, Q, row, col, first, c);
 328                         if (!isl_int_is_zero(M->row[row][first]))
 329                                 exchange(ctx, M, U, Q, row, first, col);
 330                         else
 331                                 ++first;
 332                 }
 333                 for (i = 0; i < col; ++i) {
 334                         if (isl_int_is_zero(M->row[row][i]))
 335                                 continue;
 336                         if (neg)
 337                                 isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
 338                         else
 339                                 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
 340                         if (isl_int_is_zero(c))
 341                                 continue;
 342                         subtract(M, U, Q, row, col, i, c);
 343                 }
 344                 ++col;
 345         }
 346         isl_int_clear(c);
 347
 348         return M;
 349 error:
 350         if (Q) {
 351                 isl_mat_free(ctx, *Q);
 352                 *Q = NULL;
 353         }
 354         if (U) {
 355                 isl_mat_free(ctx, *U);
 356                 *U = NULL;
 357         }
 358         return NULL;
 359 }
 360
 361 struct isl_mat *isl_mat_lin_to_aff(struct isl_ctx *ctx, struct isl_mat *mat)
 362 {
 363         int i;
 364         struct isl_mat *mat2;
 365
 366         if (!mat)
 367                 return NULL;
 368         mat2 = isl_mat_alloc(ctx, 1+mat->n_row, 1+mat->n_col);
 369         if (!mat2)
 370                 return NULL;
 371         isl_int_set_si(mat2->row[0][0], 1);
 372         isl_seq_clr(mat2->row[0]+1, mat->n_col);
 373         for (i = 0; i < mat->n_row; ++i) {
 374                 isl_int_set_si(mat2->row[1+i][0], 0);
 375                 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
 376         }
 377         isl_mat_free(ctx, mat);
 378         return mat2;
 379 }
 380
 381 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
 382 {
 383         int i;
 384
 385         for (i = 0; i < n_row; ++i)
 386                 if (!isl_int_is_zero(row[i][col]))
 387                         return i;
 388         return -1;
 389 }
 390
 391 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
 392 {
 393         int i, min = row_first_non_zero(row, n_row, col);
 394         if (min < 0)
 395                 return -1;
 396         for (i = min + 1; i < n_row; ++i) {
 397                 if (isl_int_is_zero(row[i][col]))
 398                         continue;
 399                 if (isl_int_abs_lt(row[i][col], row[min][col]))
 400                         min = i;
 401         }
 402         return min;
 403 }
 404
 405 static void inv_exchange(struct isl_ctx *ctx,
 406         struct isl_mat *left, struct isl_mat *right,
 407         unsigned i, unsigned j)
 408 {
 409         left = isl_mat_swap_rows(ctx, left, i, j);
 410         right = isl_mat_swap_rows(ctx, right, i, j);
 411 }
 412
 413 static void inv_oppose(struct isl_ctx *ctx,
 414         struct isl_mat *left, struct isl_mat *right, unsigned row)
 415 {
 416         isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
 417         isl_seq_neg(right->row[row], right->row[row], right->n_col);
 418 }
 419
 420 static void inv_subtract(struct isl_ctx *ctx,
 421         struct isl_mat *left, struct isl_mat *right,
 422         unsigned row, unsigned i, isl_int m)
 423 {
 424         isl_int_neg(m, m);
 425         isl_seq_combine(left->row[i]+row,
 426                         ctx->one, left->row[i]+row,
 427                         m, left->row[row]+row,
 428                         left->n_col-row);
 429         isl_seq_combine(right->row[i], ctx->one, right->row[i],
 430                         m, right->row[row], right->n_col);
 431 }
 432
 433 /* Compute inv(left)*right
 434  */
 435 struct isl_mat *isl_mat_inverse_product(struct isl_ctx *ctx,
 436         struct isl_mat *left, struct isl_mat *right)
 437 {
 438         int row;
 439         isl_int a, b;
 440
 441         if (!left || !right)
 442                 goto error;
 443
 444         isl_assert(ctx, left->n_row == left->n_col, goto error);
 445         isl_assert(ctx, left->n_row == right->n_row, goto error);
 446
 447         if (left->n_row == 0) {
 448                 isl_mat_free(ctx, left);
 449                 return right;
 450         }
 451
 452         left = isl_mat_cow(ctx, left);
 453         right = isl_mat_cow(ctx, right);
 454         if (!left || !right)
 455                 goto error;
 456
 457         isl_int_init(a);
 458         isl_int_init(b);
 459         for (row = 0; row < left->n_row; ++row) {
 460                 int pivot, first, i, off;
 461                 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
 462                 if (pivot < 0) {
 463                         isl_int_clear(a);
 464                         isl_int_clear(b);
 465                         isl_assert(ctx, pivot >= 0, goto error);
 466                 }
 467                 pivot += row;
 468                 if (pivot != row)
 469                         inv_exchange(ctx, left, right, pivot, row);
 470                 if (isl_int_is_neg(left->row[row][row]))
 471                         inv_oppose(ctx, left, right, row);
 472                 first = row+1;
 473                 while ((off = row_first_non_zero(left->row+first,
 474                                         left->n_row-first, row)) != -1) {
 475                         first += off;
 476                         isl_int_fdiv_q(a, left->row[first][row],
 477                                         left->row[row][row]);
 478                         inv_subtract(ctx, left, right, row, first, a);
 479                         if (!isl_int_is_zero(left->row[first][row]))
 480                                 inv_exchange(ctx, left, right, row, first);
 481                         else
 482                                 ++first;
 483                 }
 484                 for (i = 0; i < row; ++i) {
 485                         if (isl_int_is_zero(left->row[i][row]))
 486                                 continue;
 487                         isl_int_gcd(a, left->row[row][row], left->row[i][row]);
 488                         isl_int_divexact(b, left->row[i][row], a);
 489                         isl_int_divexact(a, left->row[row][row], a);
 490                         isl_int_neg(a, a);
 491                         isl_seq_combine(left->row[i]+row,
 492                                         a, left->row[i]+row,
 493                                         b, left->row[row]+row,
 494                                         left->n_col-row);
 495                         isl_seq_combine(right->row[i], a, right->row[i],
 496                                         b, right->row[row], right->n_col);
 497                 }
 498         }
 499         isl_int_clear(b);
 500
 501         isl_int_set(a, left->row[0][0]);
 502         for (row = 1; row < left->n_row; ++row)
 503                 isl_int_lcm(a, a, left->row[row][row]);
 504         if (isl_int_is_zero(a)){
 505                 isl_int_clear(a);
 506                 isl_assert(ctx, 0, goto error);
 507         }
 508         for (row = 0; row < left->n_row; ++row) {
 509                 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
 510                 if (isl_int_is_one(left->row[row][row]))
 511                         continue;
 512                 isl_seq_scale(right->row[row], right->row[row],
 513                                 left->row[row][row], right->n_col);
 514         }
 515         isl_int_clear(a);
 516
 517         isl_mat_free(ctx, left);
 518         return right;
 519 error:
 520         isl_mat_free(ctx, left);
 521         isl_mat_free(ctx, right);
 522         return NULL;
 523 }
 524
 525 void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
 526 {
 527         int i;
 528
 529         for (i = 0; i < mat->n_row; ++i)
 530                 isl_int_mul(mat->row[i][col], mat->row[i][col], m);
 531 }
 532
 533 void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
 534         isl_int m1, unsigned src1, isl_int m2, unsigned src2)
 535 {
 536         int i;
 537         isl_int tmp;
 538
 539         isl_int_init(tmp);
 540         for (i = 0; i < mat->n_row; ++i) {
 541                 isl_int_mul(tmp, m1, mat->row[i][src1]);
 542                 isl_int_addmul(tmp, m2, mat->row[i][src2]);
 543                 isl_int_set(mat->row[i][dst], tmp);
 544         }
 545         isl_int_clear(tmp);
 546 }
 547
 548 struct isl_mat *isl_mat_right_inverse(struct isl_ctx *ctx,
 549         struct isl_mat *mat)
 550 {
 551         struct isl_mat *inv;
 552         int row;
 553         isl_int a, b;
 554
 555         mat = isl_mat_cow(ctx, mat);
 556         if (!mat)
 557                 return NULL;
 558
 559         inv = isl_mat_identity(ctx, mat->n_col);
 560         inv = isl_mat_cow(ctx, inv);
 561         if (!inv)
 562                 goto error;
 563
 564         isl_int_init(a);
 565         isl_int_init(b);
 566         for (row = 0; row < mat->n_row; ++row) {
 567                 int pivot, first, i, off;
 568                 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
 569                 if (pivot < 0) {
 570                         isl_int_clear(a);
 571                         isl_int_clear(b);
 572                         goto error;
 573                 }
 574                 pivot += row;
 575                 if (pivot != row)
 576                         exchange(ctx, mat, &inv, NULL, row, pivot, row);
 577                 if (isl_int_is_neg(mat->row[row][row]))
 578                         oppose(ctx, mat, &inv, NULL, row, row);
 579                 first = row+1;
 580                 while ((off = isl_seq_first_non_zero(mat->row[row]+first,
 581                                                     mat->n_col-first)) != -1) {
 582                         first += off;
 583                         isl_int_fdiv_q(a, mat->row[row][first],
 584                                                     mat->row[row][row]);
 585                         subtract(mat, &inv, NULL, row, row, first, a);
 586                         if (!isl_int_is_zero(mat->row[row][first]))
 587                                 exchange(ctx, mat, &inv, NULL, row, row, first);
 588                         else
 589                                 ++first;
 590                 }
 591                 for (i = 0; i < row; ++i) {
 592                         if (isl_int_is_zero(mat->row[row][i]))
 593                                 continue;
 594                         isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
 595                         isl_int_divexact(b, mat->row[row][i], a);
 596                         isl_int_divexact(a, mat->row[row][row], a);
 597                         isl_int_neg(a, a);
 598                         isl_mat_col_combine(mat, i, a, i, b, row);
 599                         isl_mat_col_combine(inv, i, a, i, b, row);
 600                 }
 601         }
 602         isl_int_clear(b);
 603
 604         isl_int_set(a, mat->row[0][0]);
 605         for (row = 1; row < mat->n_row; ++row)
 606                 isl_int_lcm(a, a, mat->row[row][row]);
 607         if (isl_int_is_zero(a)){
 608                 isl_int_clear(a);
 609                 goto error;
 610         }
 611         for (row = 0; row < mat->n_row; ++row) {
 612                 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
 613                 if (isl_int_is_one(mat->row[row][row]))
 614                         continue;
 615                 isl_mat_col_scale(inv, row, mat->row[row][row]);
 616         }
 617         isl_int_clear(a);
 618
 619         isl_mat_free(ctx, mat);
 620
 621         return inv;
 622 error:
 623         isl_mat_free(ctx, mat);
 624         return NULL;
 625 }
 626
 627 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
 628         struct isl_mat *mat, unsigned i, unsigned j)
 629 {
 630         isl_int *t;
 631
 632         if (!mat)
 633                 return NULL;
 634         mat = isl_mat_cow(ctx, mat);
 635         if (!mat)
 636                 return NULL;
 637         t = mat->row[i];
 638         mat->row[i] = mat->row[j];
 639         mat->row[j] = t;
 640         return mat;
 641 }
 642
 643 struct isl_mat *isl_mat_product(struct isl_ctx *ctx,
 644         struct isl_mat *left, struct isl_mat *right)
 645 {
 646         int i, j, k;
 647         struct isl_mat *prod;
 648
 649         if (!left || !right)
 650                 goto error;
 651         isl_assert(ctx, left->n_col == right->n_row, goto error);
 652         prod = isl_mat_alloc(ctx, left->n_row, right->n_col);
 653         if (!prod)
 654                 goto error;
 655         if (left->n_col == 0) {
 656                 for (i = 0; i < prod->n_row; ++i)
 657                         isl_seq_clr(prod->row[i], prod->n_col);
 658                 return prod;
 659         }
 660         for (i = 0; i < prod->n_row; ++i) {
 661                 for (j = 0; j < prod->n_col; ++j) {
 662                         isl_int_mul(prod->row[i][j],
 663                                     left->row[i][0], right->row[0][j]);
 664                         for (k = 1; k < left->n_col; ++k)
 665                                 isl_int_addmul(prod->row[i][j],
 666                                             left->row[i][k], right->row[k][j]);
 667                 }
 668         }
 669         isl_mat_free(ctx, left);
 670         isl_mat_free(ctx, right);
 671         return prod;
 672 error:
 673         isl_mat_free(ctx, left);
 674         isl_mat_free(ctx, right);
 675         return NULL;
 676 }
 677
 678 /* Replace the variables x in bset by x' given by x = M x', with
 679  * M the matrix mat.
 680  *
 681  * If there are fewer variables x' then there are x, then we perform
 682  * the transformation in place, which that, in principle,
 683  * this frees up some extra variables as the number
 684  * of columns remains constant, but we would have to extend
 685  * the div array too as the number of rows in this array is assumed
 686  * to be equal to extra.
 687  */
 688 struct isl_basic_set *isl_basic_set_preimage(struct isl_ctx *ctx,
 689         struct isl_basic_set *bset, struct isl_mat *mat)
 690 {
 691         struct isl_mat *t;
 692         int i;
 693
 694         if (!bset || !mat)
 695                 goto error;
 696
 697         bset = isl_basic_set_cow(bset);
 698         if (!bset)
 699                 goto error;
 700
 701         isl_assert(ctx, bset->dim->nparam == 0, goto error);
 702         isl_assert(ctx, bset->n_div == 0, goto error);
 703         isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
 704
 705         if (mat->n_col > mat->n_row)
 706                 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0,
 707                                                 0, 0);
 708         else if (mat->n_col < mat->n_row) {
 709                 bset->dim = isl_dim_cow(bset->dim);
 710                 if (!bset->dim)
 711                         goto error;
 712                 bset->dim->n_out -= mat->n_row - mat->n_col;
 713         }
 714
 715         t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
 716         t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
 717         if (!t)
 718                 goto error;
 719         for (i = 0; i < bset->n_eq; ++i) {
 720                 isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
 721                 isl_seq_clr(bset->eq[i]+t->n_col, bset->extra);
 722         }
 723         isl_mat_free(ctx, t);
 724
 725         t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
 726         t = isl_mat_product(ctx, t, mat);
 727         if (!t)
 728                 goto error2;
 729         for (i = 0; i < bset->n_ineq; ++i) {
 730                 isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
 731                 isl_seq_clr(bset->ineq[i]+t->n_col, bset->extra);
 732         }
 733         isl_mat_free(ctx, t);
 734
 735         bset = isl_basic_set_simplify(bset);
 736         bset = isl_basic_set_finalize(bset);
 737
 738         return bset;
 739 error:
 740         isl_mat_free(ctx, mat);
 741 error2:
 742         isl_basic_set_free(bset);
 743         return NULL;
 744 }
 745
 746 struct isl_set *isl_set_preimage(struct isl_ctx *ctx,
 747         struct isl_set *set, struct isl_mat *mat)
 748 {
 749         int i;
 750
 751         set = isl_set_cow(set);
 752         if (!set)
 753                 return NULL;
 754
 755         ctx = set->ctx;
 756         for (i = 0; i < set->n; ++i) {
 757                 set->p[i] = isl_basic_set_preimage(ctx, set->p[i],
 758                                                     isl_mat_copy(ctx, mat));
 759                 if (!set->p[i])
 760                         goto error;
 761         }
 762         if (mat->n_col != mat->n_row) {
 763                 set->dim = isl_dim_cow(set->dim);
 764                 if (!set->dim)
 765                         goto error;
 766                 set->dim->n_out += mat->n_col;
 767                 set->dim->n_out -= mat->n_row;
 768         }
 769         isl_mat_free(ctx, mat);
 770         return set;
 771 error:
 772         isl_set_free(set);
 773         isl_mat_free(ctx, mat);
 774         return NULL;
 775 }
 776
 777 void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
 778                                 FILE *out, int indent)
 779 {
 780         int i, j;
 781
 782         if (!mat) {
 783                 fprintf(out, "null mat\n");
 784                 return;
 785         }
 786
 787         if (mat->n_row == 0)
 788                 fprintf(out, "%*s[]\n", indent, "");
 789
 790         for (i = 0; i < mat->n_row; ++i) {
 791                 if (!i)
 792                         fprintf(out, "%*s[[", indent, "");
 793                 else
 794                         fprintf(out, "%*s[", indent+1, "");
 795                 for (j = 0; j < mat->n_col; ++j) {
 796                         if (j)
 797                             fprintf(out, ",");
 798                         isl_int_print(out, mat->row[i][j], 0);
 799                 }
 800                 if (i == mat->n_row-1)
 801                         fprintf(out, "]]\n");
 802                 else
 803                         fprintf(out, "]\n");
 804         }
 805 }
 806
 807 struct isl_mat *isl_mat_drop_cols(struct isl_ctx *ctx, struct isl_mat *mat,
 808                                 unsigned col, unsigned n)
 809 {
 810         int r;
 811
 812         if (!mat)
 813                 return NULL;
 814
 815         if (col != mat->n_col-n) {
 816                 for (r = 0; r < mat->n_row; ++r)
 817                         isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
 818                                         mat->n_col - col - n);
 819         }
 820         mat->n_col -= n;
 821         return mat;
 822 }
 823
 824 struct isl_mat *isl_mat_drop_rows(struct isl_ctx *ctx, struct isl_mat *mat,
 825                                 unsigned row, unsigned n)
 826 {
 827         int r;
 828
 829         if (!mat)
 830                 return NULL;
 831
 832         for (r = row; r+n < mat->n_row; ++r)
 833                 mat->row[r] = mat->row[r+n];
 834
 835         mat->n_row -= n;
 836         return mat;
 837 }
 838
 839 void isl_mat_col_submul(struct isl_mat *mat,
 840                         int dst_col, isl_int f, int src_col)
 841 {
 842         int i;
 843
 844         for (i = 0; i < mat->n_row; ++i)
 845                 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
 846 }