add isl_vec_mat_product
[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->ctx = ctx;
28         isl_ctx_ref(ctx);
29         mat->ref = 1;
30         mat->n_row = n_row;
31         mat->n_col = n_col;
32         mat->max_col = n_col;
33         mat->flags = 0;
34
35         return mat;
36 error:
37         isl_blk_free(ctx, mat->block);
38         free(mat);
39         return NULL;
40 }
41
42 struct isl_mat *isl_mat_extend(struct isl_mat *mat,
43         unsigned n_row, unsigned n_col)
44 {
45         int i;
46         isl_int *old;
47
48         if (!mat)
49                 return NULL;
50
51         if (mat->max_col >= n_col && mat->n_row >= n_row) {
52                 if (mat->n_col < n_col)
53                         mat->n_col = n_col;
54                 return mat;
55         }
56
57         if (mat->max_col < n_col) {
58                 struct isl_mat *new_mat;
59
60                 if (n_row < mat->n_row)
61                         n_row = mat->n_row;
62                 new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
63                 if (!new_mat)
64                         goto error;
65                 for (i = 0; i < mat->n_row; ++i)
66                         isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
67                 isl_mat_free(mat);
68                 return new_mat;
69         }
70
71         mat = isl_mat_cow(mat);
72         if (!mat)
73                 goto error;
74
75         assert(mat->ref == 1);
76         old = mat->block.data;
77         mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
78         if (isl_blk_is_error(mat->block))
79                 goto error;
80         mat->row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);
81         if (!mat->row)
82                 goto error;
83
84         for (i = 0; i < mat->n_row; ++i)
85                 mat->row[i] = mat->block.data + (mat->row[i] - old);
86         for (i = mat->n_row; i < n_row; ++i)
87                 mat->row[i] = mat->block.data + i * mat->max_col;
88         mat->n_row = n_row;
89         if (mat->n_col < n_col)
90                 mat->n_col = n_col;
91
92         return mat;
93 error:
94         isl_mat_free(mat);
95         return NULL;
96 }
97
98 struct isl_mat *isl_mat_sub_alloc(struct isl_ctx *ctx, isl_int **row,
99         unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
100 {
101         int i;
102         struct isl_mat *mat;
103
104         mat = isl_alloc_type(ctx, struct isl_mat);
105         if (!mat)
106                 return NULL;
107         mat->row = isl_alloc_array(ctx, isl_int *, n_row);
108         if (!mat->row)
109                 goto error;
110         for (i = 0; i < n_row; ++i)
111                 mat->row[i] = row[first_row+i] + first_col;
112         mat->ctx = ctx;
113         isl_ctx_ref(ctx);
114         mat->ref = 1;
115         mat->n_row = n_row;
116         mat->n_col = n_col;
117         mat->block = isl_blk_empty();
118         mat->flags = ISL_MAT_BORROWED;
119         return mat;
120 error:
121         free(mat);
122         return NULL;
123 }
124
125 void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
126         unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
127 {
128         int i;
129
130         for (i = 0; i < n_row; ++i)
131                 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
132 }
133
134 void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
135         unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
136 {
137         int i;
138
139         for (i = 0; i < n_row; ++i)
140                 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
141 }
142
143 struct isl_mat *isl_mat_copy(struct isl_mat *mat)
144 {
145         if (!mat)
146                 return NULL;
147
148         mat->ref++;
149         return mat;
150 }
151
152 struct isl_mat *isl_mat_dup(struct isl_mat *mat)
153 {
154         int i;
155         struct isl_mat *mat2;
156
157         if (!mat)
158                 return NULL;
159         mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
160         if (!mat2)
161                 return NULL;
162         for (i = 0; i < mat->n_row; ++i)
163                 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
164         return mat2;
165 }
166
167 struct isl_mat *isl_mat_cow(struct isl_mat *mat)
168 {
169         struct isl_mat *mat2;
170         if (!mat)
171                 return NULL;
172
173         if (mat->ref == 1 && !ISL_F_ISSET(mat, ISL_MAT_BORROWED))
174                 return mat;
175
176         mat2 = isl_mat_dup(mat);
177         isl_mat_free(mat);
178         return mat2;
179 }
180
181 void isl_mat_free(struct isl_mat *mat)
182 {
183         if (!mat)
184                 return;
185
186         if (--mat->ref > 0)
187                 return;
188
189         if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
190                 isl_blk_free(mat->ctx, mat->block);
191         isl_ctx_deref(mat->ctx);
192         free(mat->row);
193         free(mat);
194 }
195
196 struct isl_mat *isl_mat_identity(struct isl_ctx *ctx, unsigned n_row)
197 {
198         int i;
199         struct isl_mat *mat;
200
201         mat = isl_mat_alloc(ctx, n_row, n_row);
202         if (!mat)
203                 return NULL;
204         for (i = 0; i < n_row; ++i) {
205                 isl_seq_clr(mat->row[i], i);
206                 isl_int_set_si(mat->row[i][i], 1);
207                 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
208         }
209
210         return mat;
211 }
212
213 struct isl_vec *isl_mat_vec_product(struct isl_mat *mat, struct isl_vec *vec)
214 {
215         int i;
216         struct isl_vec *prod;
217
218         if (!mat || !vec)
219                 goto error;
220
221         isl_assert(ctx, mat->n_col == vec->size, goto error);
222
223         prod = isl_vec_alloc(mat->ctx, mat->n_row);
224         if (!prod)
225                 goto error;
226
227         for (i = 0; i < prod->size; ++i)
228                 isl_seq_inner_product(mat->row[i], vec->el, vec->size,
229                                         &prod->block.data[i]);
230         isl_mat_free(mat);
231         isl_vec_free(vec);
232         return prod;
233 error:
234         isl_mat_free(mat);
235         isl_vec_free(vec);
236         return NULL;
237 }
238
239 struct isl_vec *isl_vec_mat_product(struct isl_vec *vec, struct isl_mat *mat)
240 {
241         int i, j;
242         struct isl_vec *prod;
243
244         if (!mat || !vec)
245                 goto error;
246
247         isl_assert(ctx, mat->n_row == vec->size, goto error);
248
249         prod = isl_vec_alloc(mat->ctx, mat->n_col);
250         if (!prod)
251                 goto error;
252
253         for (i = 0; i < prod->size; ++i) {
254                 isl_int_set_si(prod->el[i], 0);
255                 for (j = 0; j < vec->size; ++j)
256                         isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
257         }
258         isl_mat_free(mat);
259         isl_vec_free(vec);
260         return prod;
261 error:
262         isl_mat_free(mat);
263         isl_vec_free(vec);
264         return NULL;
265 }
266
267 struct isl_mat *isl_mat_aff_direct_sum(struct isl_mat *left,
268         struct isl_mat *right)
269 {
270         int i;
271         struct isl_mat *sum;
272
273         if (!left || !right)
274                 goto error;
275
276         isl_assert(ctx, left->n_row == right->n_row, goto error);
277         isl_assert(ctx, left->n_row >= 1, goto error);
278         isl_assert(ctx, left->n_col >= 1, goto error);
279         isl_assert(ctx, right->n_col >= 1, goto error);
280         isl_assert(ctx,
281             isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
282             goto error);
283         isl_assert(ctx,
284             isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
285             goto error);
286
287         sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
288         if (!sum)
289                 goto error;
290         isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
291         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
292         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
293
294         isl_seq_clr(sum->row[0]+1, sum->n_col-1);
295         for (i = 1; i < sum->n_row; ++i) {
296                 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
297                 isl_int_addmul(sum->row[i][0],
298                                 right->row[0][0], right->row[i][0]);
299                 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
300                                 left->n_col-1);
301                 isl_seq_scale(sum->row[i]+left->n_col,
302                                 right->row[i]+1, right->row[0][0],
303                                 right->n_col-1);
304         }
305
306         isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
307         isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
308         isl_mat_free(left);
309         isl_mat_free(right);
310         return sum;
311 error:
312         isl_mat_free(left);
313         isl_mat_free(right);
314         return NULL;
315 }
316
317 static void exchange(struct isl_mat *M, struct isl_mat **U,
318         struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
319 {
320         int r;
321         for (r = row; r < M->n_row; ++r)
322                 isl_int_swap(M->row[r][i], M->row[r][j]);
323         if (U) {
324                 for (r = 0; r < (*U)->n_row; ++r)
325                         isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
326         }
327         if (Q)
328                 isl_mat_swap_rows(*Q, i, j);
329 }
330
331 static void subtract(struct isl_mat *M, struct isl_mat **U,
332         struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
333 {
334         int r;
335         for (r = row; r < M->n_row; ++r)
336                 isl_int_submul(M->row[r][j], m, M->row[r][i]);
337         if (U) {
338                 for (r = 0; r < (*U)->n_row; ++r)
339                         isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
340         }
341         if (Q) {
342                 for (r = 0; r < (*Q)->n_col; ++r)
343                         isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
344         }
345 }
346
347 static void oppose(struct isl_mat *M, struct isl_mat **U,
348         struct isl_mat **Q, unsigned row, unsigned col)
349 {
350         int r;
351         for (r = row; r < M->n_row; ++r)
352                 isl_int_neg(M->row[r][col], M->row[r][col]);
353         if (U) {
354                 for (r = 0; r < (*U)->n_row; ++r)
355                         isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
356         }
357         if (Q)
358                 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
359 }
360
361 /* Given matrix M, compute
362  *
363  *              M U = H
364  *              M   = H Q
365  *
366  * with U and Q unimodular matrices and H a matrix in column echelon form
367  * such that on each echelon row the entries in the non-echelon column
368  * are non-negative (if neg == 0) or non-positive (if neg == 1)
369  * and stricly smaller (in absolute value) than the entries in the echelon
370  * column.
371  * If U or Q are NULL, then these matrices are not computed.
372  */
373 struct isl_mat *isl_mat_left_hermite(struct isl_mat *M, int neg,
374         struct isl_mat **U, struct isl_mat **Q)
375 {
376         isl_int c;
377         int row, col;
378
379         if (U)
380                 *U = NULL;
381         if (Q)
382                 *Q = NULL;
383         if (!M)
384                 goto error;
385         M = isl_mat_cow(M);
386         if (!M)
387                 goto error;
388         if (U) {
389                 *U = isl_mat_identity(M->ctx, M->n_col);
390                 if (!*U)
391                         goto error;
392         }
393         if (Q) {
394                 *Q = isl_mat_identity(M->ctx, M->n_col);
395                 if (!*Q)
396                         goto error;
397         }
398
399         col = 0;
400         isl_int_init(c);
401         for (row = 0; row < M->n_row; ++row) {
402                 int first, i, off;
403                 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
404                 if (first == -1)
405                         continue;
406                 first += col;
407                 if (first != col)
408                         exchange(M, U, Q, row, first, col);
409                 if (isl_int_is_neg(M->row[row][col]))
410                         oppose(M, U, Q, row, col);
411                 first = col+1;
412                 while ((off = isl_seq_first_non_zero(M->row[row]+first,
413                                                        M->n_col-first)) != -1) {
414                         first += off;
415                         isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
416                         subtract(M, U, Q, row, col, first, c);
417                         if (!isl_int_is_zero(M->row[row][first]))
418                                 exchange(M, U, Q, row, first, col);
419                         else
420                                 ++first;
421                 }
422                 for (i = 0; i < col; ++i) {
423                         if (isl_int_is_zero(M->row[row][i]))
424                                 continue;
425                         if (neg)
426                                 isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
427                         else
428                                 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
429                         if (isl_int_is_zero(c))
430                                 continue;
431                         subtract(M, U, Q, row, col, i, c);
432                 }
433                 ++col;
434         }
435         isl_int_clear(c);
436
437         return M;
438 error:
439         if (Q) {
440                 isl_mat_free(*Q);
441                 *Q = NULL;
442         }
443         if (U) {
444                 isl_mat_free(*U);
445                 *U = NULL;
446         }
447         return NULL;
448 }
449
450 struct isl_mat *isl_mat_right_kernel(struct isl_mat *mat)
451 {
452         int i, rank;
453         struct isl_mat *U = NULL;
454         struct isl_mat *K;
455
456         mat = isl_mat_left_hermite(mat, 0, &U, NULL);
457         if (!mat || !U)
458                 goto error;
459
460         for (i = 0, rank = 0; rank < mat->n_col; ++rank) {
461                 while (i < mat->n_row && isl_int_is_zero(mat->row[i][rank]))
462                         ++i;
463                 if (i >= mat->n_row)
464                         break;
465         }
466         K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank);
467         if (!K)
468                 goto error;
469         isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
470         isl_mat_free(mat);
471         isl_mat_free(U);
472         return K;
473 error:
474         isl_mat_free(mat);
475         isl_mat_free(U);
476         return NULL;
477 }
478
479 struct isl_mat *isl_mat_lin_to_aff(struct isl_mat *mat)
480 {
481         int i;
482         struct isl_mat *mat2;
483
484         if (!mat)
485                 return NULL;
486         mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
487         if (!mat2)
488                 return NULL;
489         isl_int_set_si(mat2->row[0][0], 1);
490         isl_seq_clr(mat2->row[0]+1, mat->n_col);
491         for (i = 0; i < mat->n_row; ++i) {
492                 isl_int_set_si(mat2->row[1+i][0], 0);
493                 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
494         }
495         isl_mat_free(mat);
496         return mat2;
497 }
498
499 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
500 {
501         int i;
502
503         for (i = 0; i < n_row; ++i)
504                 if (!isl_int_is_zero(row[i][col]))
505                         return i;
506         return -1;
507 }
508
509 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
510 {
511         int i, min = row_first_non_zero(row, n_row, col);
512         if (min < 0)
513                 return -1;
514         for (i = min + 1; i < n_row; ++i) {
515                 if (isl_int_is_zero(row[i][col]))
516                         continue;
517                 if (isl_int_abs_lt(row[i][col], row[min][col]))
518                         min = i;
519         }
520         return min;
521 }
522
523 static void inv_exchange(struct isl_mat *left, struct isl_mat *right,
524         unsigned i, unsigned j)
525 {
526         left = isl_mat_swap_rows(left, i, j);
527         right = isl_mat_swap_rows(right, i, j);
528 }
529
530 static void inv_oppose(
531         struct isl_mat *left, struct isl_mat *right, unsigned row)
532 {
533         isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
534         isl_seq_neg(right->row[row], right->row[row], right->n_col);
535 }
536
537 static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
538         unsigned row, unsigned i, isl_int m)
539 {
540         isl_int_neg(m, m);
541         isl_seq_combine(left->row[i]+row,
542                         left->ctx->one, left->row[i]+row,
543                         m, left->row[row]+row,
544                         left->n_col-row);
545         isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
546                         m, right->row[row], right->n_col);
547 }
548
549 /* Compute inv(left)*right
550  */
551 struct isl_mat *isl_mat_inverse_product(struct isl_mat *left,
552         struct isl_mat *right)
553 {
554         int row;
555         isl_int a, b;
556
557         if (!left || !right)
558                 goto error;
559
560         isl_assert(left->ctx, left->n_row == left->n_col, goto error);
561         isl_assert(left->ctx, left->n_row == right->n_row, goto error);
562
563         if (left->n_row == 0) {
564                 isl_mat_free(left);
565                 return right;
566         }
567
568         left = isl_mat_cow(left);
569         right = isl_mat_cow(right);
570         if (!left || !right)
571                 goto error;
572
573         isl_int_init(a);
574         isl_int_init(b);
575         for (row = 0; row < left->n_row; ++row) {
576                 int pivot, first, i, off;
577                 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
578                 if (pivot < 0) {
579                         isl_int_clear(a);
580                         isl_int_clear(b);
581                         isl_assert(ctx, pivot >= 0, goto error);
582                 }
583                 pivot += row;
584                 if (pivot != row)
585                         inv_exchange(left, right, pivot, row);
586                 if (isl_int_is_neg(left->row[row][row]))
587                         inv_oppose(left, right, row);
588                 first = row+1;
589                 while ((off = row_first_non_zero(left->row+first,
590                                         left->n_row-first, row)) != -1) {
591                         first += off;
592                         isl_int_fdiv_q(a, left->row[first][row],
593                                         left->row[row][row]);
594                         inv_subtract(left, right, row, first, a);
595                         if (!isl_int_is_zero(left->row[first][row]))
596                                 inv_exchange(left, right, row, first);
597                         else
598                                 ++first;
599                 }
600                 for (i = 0; i < row; ++i) {
601                         if (isl_int_is_zero(left->row[i][row]))
602                                 continue;
603                         isl_int_gcd(a, left->row[row][row], left->row[i][row]);
604                         isl_int_divexact(b, left->row[i][row], a);
605                         isl_int_divexact(a, left->row[row][row], a);
606                         isl_int_neg(a, a);
607                         isl_seq_combine(left->row[i]+row,
608                                         a, left->row[i]+row,
609                                         b, left->row[row]+row,
610                                         left->n_col-row);
611                         isl_seq_combine(right->row[i], a, right->row[i],
612                                         b, right->row[row], right->n_col);
613                 }
614         }
615         isl_int_clear(b);
616
617         isl_int_set(a, left->row[0][0]);
618         for (row = 1; row < left->n_row; ++row)
619                 isl_int_lcm(a, a, left->row[row][row]);
620         if (isl_int_is_zero(a)){
621                 isl_int_clear(a);
622                 isl_assert(ctx, 0, goto error);
623         }
624         for (row = 0; row < left->n_row; ++row) {
625                 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
626                 if (isl_int_is_one(left->row[row][row]))
627                         continue;
628                 isl_seq_scale(right->row[row], right->row[row],
629                                 left->row[row][row], right->n_col);
630         }
631         isl_int_clear(a);
632
633         isl_mat_free(left);
634         return right;
635 error:
636         isl_mat_free(left);
637         isl_mat_free(right);
638         return NULL;
639 }
640
641 void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
642 {
643         int i;
644
645         for (i = 0; i < mat->n_row; ++i)
646                 isl_int_mul(mat->row[i][col], mat->row[i][col], m);
647 }
648
649 void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
650         isl_int m1, unsigned src1, isl_int m2, unsigned src2)
651 {
652         int i;
653         isl_int tmp;
654
655         isl_int_init(tmp);
656         for (i = 0; i < mat->n_row; ++i) {
657                 isl_int_mul(tmp, m1, mat->row[i][src1]);
658                 isl_int_addmul(tmp, m2, mat->row[i][src2]);
659                 isl_int_set(mat->row[i][dst], tmp);
660         }
661         isl_int_clear(tmp);
662 }
663
664 struct isl_mat *isl_mat_right_inverse(struct isl_mat *mat)
665 {
666         struct isl_mat *inv;
667         int row;
668         isl_int a, b;
669
670         mat = isl_mat_cow(mat);
671         if (!mat)
672                 return NULL;
673
674         inv = isl_mat_identity(mat->ctx, mat->n_col);
675         inv = isl_mat_cow(inv);
676         if (!inv)
677                 goto error;
678
679         isl_int_init(a);
680         isl_int_init(b);
681         for (row = 0; row < mat->n_row; ++row) {
682                 int pivot, first, i, off;
683                 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
684                 if (pivot < 0) {
685                         isl_int_clear(a);
686                         isl_int_clear(b);
687                         goto error;
688                 }
689                 pivot += row;
690                 if (pivot != row)
691                         exchange(mat, &inv, NULL, row, pivot, row);
692                 if (isl_int_is_neg(mat->row[row][row]))
693                         oppose(mat, &inv, NULL, row, row);
694                 first = row+1;
695                 while ((off = isl_seq_first_non_zero(mat->row[row]+first,
696                                                     mat->n_col-first)) != -1) {
697                         first += off;
698                         isl_int_fdiv_q(a, mat->row[row][first],
699                                                     mat->row[row][row]);
700                         subtract(mat, &inv, NULL, row, row, first, a);
701                         if (!isl_int_is_zero(mat->row[row][first]))
702                                 exchange(mat, &inv, NULL, row, row, first);
703                         else
704                                 ++first;
705                 }
706                 for (i = 0; i < row; ++i) {
707                         if (isl_int_is_zero(mat->row[row][i]))
708                                 continue;
709                         isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
710                         isl_int_divexact(b, mat->row[row][i], a);
711                         isl_int_divexact(a, mat->row[row][row], a);
712                         isl_int_neg(a, a);
713                         isl_mat_col_combine(mat, i, a, i, b, row);
714                         isl_mat_col_combine(inv, i, a, i, b, row);
715                 }
716         }
717         isl_int_clear(b);
718
719         isl_int_set(a, mat->row[0][0]);
720         for (row = 1; row < mat->n_row; ++row)
721                 isl_int_lcm(a, a, mat->row[row][row]);
722         if (isl_int_is_zero(a)){
723                 isl_int_clear(a);
724                 goto error;
725         }
726         for (row = 0; row < mat->n_row; ++row) {
727                 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
728                 if (isl_int_is_one(mat->row[row][row]))
729                         continue;
730                 isl_mat_col_scale(inv, row, mat->row[row][row]);
731         }
732         isl_int_clear(a);
733
734         isl_mat_free(mat);
735
736         return inv;
737 error:
738         isl_mat_free(mat);
739         return NULL;
740 }
741
742 struct isl_mat *isl_mat_transpose(struct isl_mat *mat)
743 {
744         struct isl_mat *transpose = NULL;
745         int i, j;
746
747         if (mat->n_col == mat->n_row) {
748                 mat = isl_mat_cow(mat);
749                 if (!mat)
750                         return NULL;
751                 for (i = 0; i < mat->n_row; ++i)
752                         for (j = i + 1; j < mat->n_col; ++j)
753                                 isl_int_swap(mat->row[i][j], mat->row[j][i]);
754                 return mat;
755         }
756         transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);
757         if (!transpose)
758                 goto error;
759         for (i = 0; i < mat->n_row; ++i)
760                 for (j = 0; j < mat->n_col; ++j)
761                         isl_int_set(transpose->row[j][i], mat->row[i][j]);
762         isl_mat_free(mat);
763         return transpose;
764 error:
765         isl_mat_free(mat);
766         return NULL;
767 }
768
769 struct isl_mat *isl_mat_swap_cols(struct isl_mat *mat, unsigned i, unsigned j)
770 {
771         int r;
772
773         mat = isl_mat_cow(mat);
774         if (!mat)
775                 return NULL;
776         isl_assert(ctx, i < mat->n_col, goto error);
777         isl_assert(ctx, j < mat->n_col, goto error);
778
779         for (r = 0; r < mat->n_row; ++r)
780                 isl_int_swap(mat->row[r][i], mat->row[r][j]);
781         return mat;
782 error:
783         isl_mat_free(mat);
784         return NULL;
785 }
786
787 struct isl_mat *isl_mat_swap_rows(struct isl_mat *mat, unsigned i, unsigned j)
788 {
789         isl_int *t;
790
791         if (!mat)
792                 return NULL;
793         mat = isl_mat_cow(mat);
794         if (!mat)
795                 return NULL;
796         t = mat->row[i];
797         mat->row[i] = mat->row[j];
798         mat->row[j] = t;
799         return mat;
800 }
801
802 struct isl_mat *isl_mat_product(struct isl_mat *left, struct isl_mat *right)
803 {
804         int i, j, k;
805         struct isl_mat *prod;
806
807         if (!left || !right)
808                 goto error;
809         isl_assert(ctx, left->n_col == right->n_row, goto error);
810         prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
811         if (!prod)
812                 goto error;
813         if (left->n_col == 0) {
814                 for (i = 0; i < prod->n_row; ++i)
815                         isl_seq_clr(prod->row[i], prod->n_col);
816                 return prod;
817         }
818         for (i = 0; i < prod->n_row; ++i) {
819                 for (j = 0; j < prod->n_col; ++j) {
820                         isl_int_mul(prod->row[i][j],
821                                     left->row[i][0], right->row[0][j]);
822                         for (k = 1; k < left->n_col; ++k)
823                                 isl_int_addmul(prod->row[i][j],
824                                             left->row[i][k], right->row[k][j]);
825                 }
826         }
827         isl_mat_free(left);
828         isl_mat_free(right);
829         return prod;
830 error:
831         isl_mat_free(left);
832         isl_mat_free(right);
833         return NULL;
834 }
835
836 /* Replace the variables x in the rows q by x' given by x = M x',
837  * with M the matrix mat.
838  *
839  * If the number of new variables is greater than the original
840  * number of variables, then the rows q have already been
841  * preextended.  If the new number is smaller, then the coefficients
842  * of the divs, which are not changed, need to be shifted down.
843  * The row q may be the equalities, the inequalities or the
844  * div expressions.  In the latter case, has_div is true and
845  * we need to take into account the extra denominator column.
846  */
847 static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,
848         unsigned n_div, int has_div, struct isl_mat *mat)
849 {
850         int i;
851         struct isl_mat *t;
852         int e;
853
854         if (mat->n_col >= mat->n_row)
855                 e = 0;
856         else
857                 e = mat->n_row - mat->n_col;
858         if (has_div)
859                 for (i = 0; i < n; ++i)
860                         isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);
861         t = isl_mat_sub_alloc(mat->ctx, q, 0, n, has_div, mat->n_row);
862         t = isl_mat_product(t, mat);
863         if (!t)
864                 return -1;
865         for (i = 0; i < n; ++i) {
866                 isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
867                 isl_seq_cpy(q[i] + has_div + t->n_col,
868                             q[i] + has_div + t->n_col + e, n_div);
869                 isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
870         }
871         isl_mat_free(t);
872         return 0;
873 }
874
875 /* Replace the variables x in bset by x' given by x = M x', with
876  * M the matrix mat.
877  *
878  * If there are fewer variables x' then there are x, then we perform
879  * the transformation in place, which that, in principle,
880  * this frees up some extra variables as the number
881  * of columns remains constant, but we would have to extend
882  * the div array too as the number of rows in this array is assumed
883  * to be equal to extra.
884  */
885 struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset,
886         struct isl_mat *mat)
887 {
888         struct isl_ctx *ctx;
889
890         if (!bset || !mat)
891                 goto error;
892
893         ctx = bset->ctx;
894         bset = isl_basic_set_cow(bset);
895         if (!bset)
896                 goto error;
897
898         isl_assert(ctx, bset->dim->nparam == 0, goto error);
899         isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
900
901         if (mat->n_col > mat->n_row)
902                 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0,
903                                                 0, 0);
904         else if (mat->n_col < mat->n_row) {
905                 bset->dim = isl_dim_cow(bset->dim);
906                 if (!bset->dim)
907                         goto error;
908                 bset->dim->n_out -= mat->n_row - mat->n_col;
909         }
910
911         if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
912                         isl_mat_copy(mat)) < 0)
913                 goto error;
914
915         if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
916                         isl_mat_copy(mat)) < 0)
917                 goto error;
918
919         if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)
920                 goto error2;
921
922         ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
923         ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
924         ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
925         ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
926         ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
927
928         bset = isl_basic_set_simplify(bset);
929         bset = isl_basic_set_finalize(bset);
930
931         return bset;
932 error:
933         isl_mat_free(mat);
934 error2:
935         isl_basic_set_free(bset);
936         return NULL;
937 }
938
939 struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat)
940 {
941         struct isl_ctx *ctx;
942         int i;
943
944         set = isl_set_cow(set);
945         if (!set)
946                 return NULL;
947
948         ctx = set->ctx;
949         for (i = 0; i < set->n; ++i) {
950                 set->p[i] = isl_basic_set_preimage(set->p[i],
951                                                     isl_mat_copy(mat));
952                 if (!set->p[i])
953                         goto error;
954         }
955         if (mat->n_col != mat->n_row) {
956                 set->dim = isl_dim_cow(set->dim);
957                 if (!set->dim)
958                         goto error;
959                 set->dim->n_out += mat->n_col;
960                 set->dim->n_out -= mat->n_row;
961         }
962         isl_mat_free(mat);
963         ISL_F_CLR(set, ISL_SET_NORMALIZED);
964         return set;
965 error:
966         isl_set_free(set);
967         isl_mat_free(mat);
968         return NULL;
969 }
970
971 void isl_mat_dump(struct isl_mat *mat, FILE *out, int indent)
972 {
973         int i, j;
974
975         if (!mat) {
976                 fprintf(out, "%*snull mat\n", indent, "");
977                 return;
978         }
979
980         if (mat->n_row == 0)
981                 fprintf(out, "%*s[]\n", indent, "");
982
983         for (i = 0; i < mat->n_row; ++i) {
984                 if (!i)
985                         fprintf(out, "%*s[[", indent, "");
986                 else
987                         fprintf(out, "%*s[", indent+1, "");
988                 for (j = 0; j < mat->n_col; ++j) {
989                         if (j)
990                             fprintf(out, ",");
991                         isl_int_print(out, mat->row[i][j], 0);
992                 }
993                 if (i == mat->n_row-1)
994                         fprintf(out, "]]\n");
995                 else
996                         fprintf(out, "]\n");
997         }
998 }
999
1000 struct isl_mat *isl_mat_drop_cols(struct isl_mat *mat, unsigned col, unsigned n)
1001 {
1002         int r;
1003
1004         mat = isl_mat_cow(mat);
1005         if (!mat)
1006                 return NULL;
1007
1008         if (col != mat->n_col-n) {
1009                 for (r = 0; r < mat->n_row; ++r)
1010                         isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1011                                         mat->n_col - col - n);
1012         }
1013         mat->n_col -= n;
1014         return mat;
1015 }
1016
1017 struct isl_mat *isl_mat_drop_rows(struct isl_mat *mat, unsigned row, unsigned n)
1018 {
1019         int r;
1020
1021         mat = isl_mat_cow(mat);
1022         if (!mat)
1023                 return NULL;
1024
1025         for (r = row; r+n < mat->n_row; ++r)
1026                 mat->row[r] = mat->row[r+n];
1027
1028         mat->n_row -= n;
1029         return mat;
1030 }
1031
1032 void isl_mat_col_submul(struct isl_mat *mat,
1033                         int dst_col, isl_int f, int src_col)
1034 {
1035         int i;
1036
1037         for (i = 0; i < mat->n_row; ++i)
1038                 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1039 }
1040
1041 void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1042 {
1043         int i;
1044
1045         for (i = 0; i < mat->n_row; ++i)
1046                 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1047 }
1048
1049 struct isl_mat *isl_mat_unimodular_complete(struct isl_mat *M, int row)
1050 {
1051         int r;
1052         struct isl_mat *H = NULL, *Q = NULL;
1053
1054         isl_assert(ctx, M->n_row == M->n_col, goto error);
1055         M->n_row = row;
1056         H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1057         M->n_row = M->n_col;
1058         if (!H)
1059                 goto error;
1060         for (r = 0; r < row; ++r)
1061                 isl_assert(ctx, isl_int_is_one(H->row[r][r]), goto error);
1062         for (r = row; r < M->n_row; ++r)
1063                 isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1064         isl_mat_free(H);
1065         isl_mat_free(Q);
1066         return M;
1067 error:
1068         isl_mat_free(H);
1069         isl_mat_free(Q);
1070         isl_mat_free(M);
1071         return NULL;
1072 }