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