?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
isl 0.05.1
[isl.git] / isl_mat.c
blob18b996684bed32470c0d3e119bce2b28ea1987c7
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 isl_mat_free(inv);
860 return NULL;
861 }
862
863 struct isl_mat *isl_mat_transpose(struct isl_mat *mat)
864 {
865 struct isl_mat *transpose = NULL;
866 int i, j;
867
868 if (mat->n_col == mat->n_row) {
869 mat = isl_mat_cow(mat);
870 if (!mat)
871 return NULL;
872 for (i = 0; i < mat->n_row; ++i)
873 for (j = i + 1; j < mat->n_col; ++j)
874 isl_int_swap(mat->row[i][j], mat->row[j][i]);
875 return mat;
876 }
877 transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);
878 if (!transpose)
879 goto error;
880 for (i = 0; i < mat->n_row; ++i)
881 for (j = 0; j < mat->n_col; ++j)
882 isl_int_set(transpose->row[j][i], mat->row[i][j]);
883 isl_mat_free(mat);
884 return transpose;
885 error:
886 isl_mat_free(mat);
887 return NULL;
888 }
889
890 struct isl_mat *isl_mat_swap_cols(struct isl_mat *mat, unsigned i, unsigned j)
891 {
892 int r;
893
894 mat = isl_mat_cow(mat);
895 if (!mat)
896 return NULL;
897 isl_assert(mat->ctx, i < mat->n_col, goto error);
898 isl_assert(mat->ctx, j < mat->n_col, goto error);
899
900 for (r = 0; r < mat->n_row; ++r)
901 isl_int_swap(mat->row[r][i], mat->row[r][j]);
902 return mat;
903 error:
904 isl_mat_free(mat);
905 return NULL;
906 }
907
908 struct isl_mat *isl_mat_swap_rows(struct isl_mat *mat, unsigned i, unsigned j)
909 {
910 isl_int *t;
911
912 if (!mat)
913 return NULL;
914 mat = isl_mat_cow(mat);
915 if (!mat)
916 return NULL;
917 t = mat->row[i];
918 mat->row[i] = mat->row[j];
919 mat->row[j] = t;
920 return mat;
921 }
922
923 struct isl_mat *isl_mat_product(struct isl_mat *left, struct isl_mat *right)
924 {
925 int i, j, k;
926 struct isl_mat *prod;
927
928 if (!left || !right)
929 goto error;
930 isl_assert(left->ctx, left->n_col == right->n_row, goto error);
931 prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
932 if (!prod)
933 goto error;
934 if (left->n_col == 0) {
935 for (i = 0; i < prod->n_row; ++i)
936 isl_seq_clr(prod->row[i], prod->n_col);
937 return prod;
938 }
939 for (i = 0; i < prod->n_row; ++i) {
940 for (j = 0; j < prod->n_col; ++j) {
941 isl_int_mul(prod->row[i][j],
942 left->row[i][0], right->row[0][j]);
943 for (k = 1; k < left->n_col; ++k)
944 isl_int_addmul(prod->row[i][j],
945 left->row[i][k], right->row[k][j]);
946 }
947 }
948 isl_mat_free(left);
949 isl_mat_free(right);
950 return prod;
951 error:
952 isl_mat_free(left);
953 isl_mat_free(right);
954 return NULL;
955 }
956
957 /* Replace the variables x in the rows q by x' given by x = M x',
958 * with M the matrix mat.
959 *
960 * If the number of new variables is greater than the original
961 * number of variables, then the rows q have already been
962 * preextended. If the new number is smaller, then the coefficients
963 * of the divs, which are not changed, need to be shifted down.
964 * The row q may be the equalities, the inequalities or the
965 * div expressions. In the latter case, has_div is true and
966 * we need to take into account the extra denominator column.
967 */
968 static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,
969 unsigned n_div, int has_div, struct isl_mat *mat)
970 {
971 int i;
972 struct isl_mat *t;
973 int e;
974
975 if (mat->n_col >= mat->n_row)
976 e = 0;
977 else
978 e = mat->n_row - mat->n_col;
979 if (has_div)
980 for (i = 0; i < n; ++i)
981 isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);
982 t = isl_mat_sub_alloc(mat->ctx, q, 0, n, has_div, mat->n_row);
983 t = isl_mat_product(t, mat);
984 if (!t)
985 return -1;
986 for (i = 0; i < n; ++i) {
987 isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
988 isl_seq_cpy(q[i] + has_div + t->n_col,
989 q[i] + has_div + t->n_col + e, n_div);
990 isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
991 }
992 isl_mat_free(t);
993 return 0;
994 }
995
996 /* Replace the variables x in bset by x' given by x = M x', with
997 * M the matrix mat.
998 *
999 * If there are fewer variables x' then there are x, then we perform
1000 * the transformation in place, which that, in principle,
1001 * this frees up some extra variables as the number
1002 * of columns remains constant, but we would have to extend
1003 * the div array too as the number of rows in this array is assumed
1004 * to be equal to extra.
1005 */
1006 struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset,
1007 struct isl_mat *mat)
1008 {
1009 struct isl_ctx *ctx;
1010
1011 if (!bset || !mat)
1012 goto error;
1013
1014 ctx = bset->ctx;
1015 bset = isl_basic_set_cow(bset);
1016 if (!bset)
1017 goto error;
1018
1019 isl_assert(ctx, bset->dim->nparam == 0, goto error);
1020 isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1021 isl_assert(ctx, mat->n_col > 0, goto error);
1022
1023 if (mat->n_col > mat->n_row) {
1024 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1025 if (!bset)
1026 goto error;
1027 } else if (mat->n_col < mat->n_row) {
1028 bset->dim = isl_dim_cow(bset->dim);
1029 if (!bset->dim)
1030 goto error;
1031 bset->dim->n_out -= mat->n_row - mat->n_col;
1032 }
1033
1034 if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
1035 isl_mat_copy(mat)) < 0)
1036 goto error;
1037
1038 if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
1039 isl_mat_copy(mat)) < 0)
1040 goto error;
1041
1042 if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)
1043 goto error2;
1044
1045 ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
1046 ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1047 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1048 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1049 ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1050
1051 bset = isl_basic_set_simplify(bset);
1052 bset = isl_basic_set_finalize(bset);
1053
1054 return bset;
1055 error:
1056 isl_mat_free(mat);
1057 error2:
1058 isl_basic_set_free(bset);
1059 return NULL;
1060 }
1061
1062 struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat)
1063 {
1064 struct isl_ctx *ctx;
1065 int i;
1066
1067 set = isl_set_cow(set);
1068 if (!set)
1069 return NULL;
1070
1071 ctx = set->ctx;
1072 for (i = 0; i < set->n; ++i) {
1073 set->p[i] = isl_basic_set_preimage(set->p[i],
1074 isl_mat_copy(mat));
1075 if (!set->p[i])
1076 goto error;
1077 }
1078 if (mat->n_col != mat->n_row) {
1079 set->dim = isl_dim_cow(set->dim);
1080 if (!set->dim)
1081 goto error;
1082 set->dim->n_out += mat->n_col;
1083 set->dim->n_out -= mat->n_row;
1084 }
1085 isl_mat_free(mat);
1086 ISL_F_CLR(set, ISL_SET_NORMALIZED);
1087 return set;
1088 error:
1089 isl_set_free(set);
1090 isl_mat_free(mat);
1091 return NULL;
1092 }
1093
1094 void isl_mat_dump(struct isl_mat *mat, FILE *out, int indent)
1095 {
1096 int i, j;
1097
1098 if (!mat) {
1099 fprintf(out, "%*snull mat\n", indent, "");
1100 return;
1101 }
1102
1103 if (mat->n_row == 0)
1104 fprintf(out, "%*s[]\n", indent, "");
1105
1106 for (i = 0; i < mat->n_row; ++i) {
1107 if (!i)
1108 fprintf(out, "%*s[[", indent, "");
1109 else
1110 fprintf(out, "%*s[", indent+1, "");
1111 for (j = 0; j < mat->n_col; ++j) {
1112 if (j)
1113 fprintf(out, ",");
1114 isl_int_print(out, mat->row[i][j], 0);
1115 }
1116 if (i == mat->n_row-1)
1117 fprintf(out, "]]\n");
1118 else
1119 fprintf(out, "]\n");
1120 }
1121 }
1122
1123 struct isl_mat *isl_mat_drop_cols(struct isl_mat *mat, unsigned col, unsigned n)
1124 {
1125 int r;
1126
1127 mat = isl_mat_cow(mat);
1128 if (!mat)
1129 return NULL;
1130
1131 if (col != mat->n_col-n) {
1132 for (r = 0; r < mat->n_row; ++r)
1133 isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1134 mat->n_col - col - n);
1135 }
1136 mat->n_col -= n;
1137 return mat;
1138 }
1139
1140 struct isl_mat *isl_mat_drop_rows(struct isl_mat *mat, unsigned row, unsigned n)
1141 {
1142 int r;
1143
1144 mat = isl_mat_cow(mat);
1145 if (!mat)
1146 return NULL;
1147
1148 for (r = row; r+n < mat->n_row; ++r)
1149 mat->row[r] = mat->row[r+n];
1150
1151 mat->n_row -= n;
1152 return mat;
1153 }
1154
1155 __isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1156 unsigned col, unsigned n)
1157 {
1158 isl_mat *ext;
1159
1160 if (!mat)
1161 return NULL;
1162 if (n == 0)
1163 return mat;
1164
1165 ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1166 if (!ext)
1167 goto error;
1168
1169 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1170 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1171 col + n, col, mat->n_col - col);
1172
1173 isl_mat_free(mat);
1174 return ext;
1175 error:
1176 isl_mat_free(mat);
1177 return NULL;
1178 }
1179
1180 __isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1181 unsigned first, unsigned n)
1182 {
1183 int i;
1184
1185 if (!mat)
1186 return NULL;
1187 mat = isl_mat_insert_cols(mat, first, n);
1188 if (!mat)
1189 return NULL;
1190
1191 for (i = 0; i < mat->n_row; ++i)
1192 isl_seq_clr(mat->row[i] + first, n);
1193
1194 return mat;
1195 }
1196
1197 __isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1198 {
1199 if (!mat)
1200 return NULL;
1201
1202 return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1203 }
1204
1205 __isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1206 unsigned row, unsigned n)
1207 {
1208 isl_mat *ext;
1209
1210 if (!mat)
1211 return NULL;
1212 if (n == 0)
1213 return mat;
1214
1215 ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1216 if (!ext)
1217 goto error;
1218
1219 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1220 isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1221 mat->n_row - row, 0, 0, mat->n_col);
1222
1223 isl_mat_free(mat);
1224 return ext;
1225 error:
1226 isl_mat_free(mat);
1227 return NULL;
1228 }
1229
1230 __isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1231 {
1232 if (!mat)
1233 return NULL;
1234
1235 return isl_mat_insert_rows(mat, mat->n_row, n);
1236 }
1237
1238 void isl_mat_col_submul(struct isl_mat *mat,
1239 int dst_col, isl_int f, int src_col)
1240 {
1241 int i;
1242
1243 for (i = 0; i < mat->n_row; ++i)
1244 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1245 }
1246
1247 void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)
1248 {
1249 int i;
1250
1251 if (!mat)
1252 return;
1253
1254 for (i = 0; i < mat->n_row; ++i)
1255 isl_int_add(mat->row[i][dst_col],
1256 mat->row[i][dst_col], mat->row[i][src_col]);
1257 }
1258
1259 void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1260 {
1261 int i;
1262
1263 for (i = 0; i < mat->n_row; ++i)
1264 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1265 }
1266
1267 struct isl_mat *isl_mat_unimodular_complete(struct isl_mat *M, int row)
1268 {
1269 int r;
1270 struct isl_mat *H = NULL, *Q = NULL;
1271
1272 if (!M)
1273 return NULL;
1274
1275 isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1276 M->n_row = row;
1277 H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1278 M->n_row = M->n_col;
1279 if (!H)
1280 goto error;
1281 for (r = 0; r < row; ++r)
1282 isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1283 for (r = row; r < M->n_row; ++r)
1284 isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1285 isl_mat_free(H);
1286 isl_mat_free(Q);
1287 return M;
1288 error:
1289 isl_mat_free(H);
1290 isl_mat_free(Q);
1291 isl_mat_free(M);
1292 return NULL;
1293 }
1294
1295 __isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1296 __isl_take isl_mat *bot)
1297 {
1298 struct isl_mat *mat;
1299
1300 if (!top || !bot)
1301 goto error;
1302
1303 isl_assert(top->ctx, top->n_col == bot->n_col, goto error);
1304 if (top->n_row == 0) {
1305 isl_mat_free(top);
1306 return bot;
1307 }
1308 if (bot->n_row == 0) {
1309 isl_mat_free(bot);
1310 return top;
1311 }
1312
1313 mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1314 if (!mat)
1315 goto error;
1316 isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1317 0, 0, mat->n_col);
1318 isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1319 0, 0, mat->n_col);
1320 isl_mat_free(top);
1321 isl_mat_free(bot);
1322 return mat;
1323 error:
1324 isl_mat_free(top);
1325 isl_mat_free(bot);
1326 return NULL;
1327 }
1328
1329 int isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1330 {
1331 int i;
1332
1333 if (!mat1 || !mat2)
1334 return -1;
1335
1336 if (mat1->n_row != mat2->n_row)
1337 return 0;
1338
1339 if (mat1->n_col != mat2->n_col)
1340 return 0;
1341
1342 for (i = 0; i < mat1->n_row; ++i)
1343 if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col))
1344 return 0;
1345
1346 return 1;
1347 }
1348
1349 __isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)
1350 {
1351 struct isl_mat *mat;
1352
1353 if (!vec)
1354 return NULL;
1355 mat = isl_mat_alloc(vec->ctx, 1, vec->size);
1356 if (!mat)
1357 goto error;
1358
1359 isl_seq_cpy(mat->row[0], vec->el, vec->size);
1360
1361 isl_vec_free(vec);
1362 return mat;
1363 error:
1364 isl_vec_free(vec);
1365 return NULL;
1366 }
1367
1368 __isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,
1369 __isl_take isl_vec *bot)
1370 {
1371 return isl_mat_concat(top, isl_mat_from_row_vec(bot));
1372 }
1373
1374 __isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,
1375 unsigned dst_col, unsigned src_col, unsigned n)
1376 {
1377 isl_mat *res;
1378
1379 if (!mat)
1380 return NULL;
1381 if (n == 0 || dst_col == src_col)
1382 return mat;
1383
1384 res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1385 if (!res)
1386 goto error;
1387
1388 if (dst_col < src_col) {
1389 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1390 0, 0, dst_col);
1391 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1392 dst_col, src_col, n);
1393 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1394 dst_col + n, dst_col, src_col - dst_col);
1395 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1396 src_col + n, src_col + n,
1397 res->n_col - src_col - n);
1398 } else {
1399 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1400 0, 0, src_col);
1401 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1402 src_col, src_col + n, dst_col - src_col);
1403 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1404 dst_col, src_col, n);
1405 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1406 dst_col + n, dst_col + n,
1407 res->n_col - dst_col - n);
1408 }
1409 isl_mat_free(mat);
1410
1411 return res;
1412 error:
1413 isl_mat_free(mat);
1414 return NULL;
1415 }
1416
1417 void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1418 {
1419 int i;
1420 isl_int g;
1421
1422 isl_int_set_si(*gcd, 0);
1423 if (!mat)
1424 return;
1425
1426 isl_int_init(g);
1427 for (i = 0; i < mat->n_row; ++i) {
1428 isl_seq_gcd(mat->row[i], mat->n_col, &g);
1429 isl_int_gcd(*gcd, *gcd, g);
1430 }
1431 isl_int_clear(g);
1432 }
1433
1434 __isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1435 {
1436 int i;
1437
1438 if (!mat)
1439 return NULL;
1440
1441 for (i = 0; i < mat->n_row; ++i)
1442 isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1443
1444 return mat;
1445 }
1446
1447 __isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1448 {
1449 isl_int gcd;
1450
1451 if (!mat)
1452 return NULL;
1453
1454 isl_int_init(gcd);
1455 isl_mat_gcd(mat, &gcd);
1456 mat = isl_mat_scale_down(mat, gcd);
1457 isl_int_clear(gcd);
1458
1459 return mat;
1460 }
Integer set library