add isl_mat_col_submul
[platform/upstream/isl.git] / 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 }