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