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