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