?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
isl_map_convex_hull: check for boundedness using recession cone
[isl.git] / isl_mat.c
blobe15967f0caa45af843943d62f5fcc68cf09033eb
1 #include "isl_dim.h"
2 #include "isl_seq.h"
3 #include "isl_mat.h"
4 #include "isl_map_private.h"
5
6 struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
7 unsigned n_row, unsigned n_col)
8 {
9 int i;
10 struct isl_mat *mat;
11
12 mat = isl_alloc_type(ctx, struct isl_mat);
13 if (!mat)
14 return NULL;
15
16 mat->row = NULL;
17 mat->block = isl_blk_alloc(ctx, n_row * n_col);
18 if (isl_blk_is_error(mat->block))
19 goto error;
20 mat->row = isl_alloc_array(ctx, isl_int *, n_row);
21 if (!mat->row)
22 goto error;
23
24 for (i = 0; i < n_row; ++i)
25 mat->row[i] = mat->block.data + i * n_col;
26
27 mat->ref = 1;
28 mat->n_row = n_row;
29 mat->n_col = n_col;
30 mat->flags = 0;
31
32 return mat;
33 error:
34 isl_blk_free(ctx, mat->block);
35 free(mat);
36 return NULL;
37 }
38
39 struct isl_mat *isl_mat_extend(struct isl_ctx *ctx, struct isl_mat *mat,
40 unsigned n_row, unsigned n_col)
41 {
42 int i;
43
44 if (!mat)
45 return NULL;
46
47 if (mat->n_col >= n_col && mat->n_row >= n_row)
48 return mat;
49
50 if (mat->n_col < n_col) {
51 struct isl_mat *new_mat;
52
53 new_mat = isl_mat_alloc(ctx, n_row, n_col);
54 if (!new_mat)
55 goto error;
56 for (i = 0; i < mat->n_row; ++i)
57 isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
58 isl_mat_free(ctx, mat);
59 return new_mat;
60 }
61
62 mat = isl_mat_cow(ctx, mat);
63 if (!mat)
64 goto error;
65
66 assert(mat->ref == 1);
67 mat->block = isl_blk_extend(ctx, mat->block, n_row * mat->n_col);
68 if (isl_blk_is_error(mat->block))
69 goto error;
70 mat->row = isl_realloc_array(ctx, mat->row, isl_int *, n_row);
71 if (!mat->row)
72 goto error;
73
74 for (i = 0; i < n_row; ++i)
75 mat->row[i] = mat->block.data + i * mat->n_col;
76 mat->n_row = n_row;
77
78 return mat;
79 error:
80 isl_mat_free(ctx, mat);
81 return NULL;
82 }
83
84 struct isl_mat *isl_mat_sub_alloc(struct isl_ctx *ctx, isl_int **row,
85 unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
86 {
87 int i;
88 struct isl_mat *mat;
89
90 mat = isl_alloc_type(ctx, struct isl_mat);
91 if (!mat)
92 return NULL;
93 mat->row = isl_alloc_array(ctx, isl_int *, n_row);
94 if (!mat->row)
95 goto error;
96 for (i = 0; i < n_row; ++i)
97 mat->row[i] = row[first_row+i] + first_col;
98 mat->ref = 1;
99 mat->n_row = n_row;
100 mat->n_col = n_col;
101 mat->block = isl_blk_empty();
102 mat->flags = ISL_MAT_BORROWED;
103 return mat;
104 error:
105 free(mat);
106 return NULL;
107 }
108
109 void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
110 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
111 {
112 int i;
113
114 for (i = 0; i < n_row; ++i)
115 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
116 }
117
118 void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
119 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
120 {
121 int i;
122
123 for (i = 0; i < n_row; ++i)
124 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
125 }
126
127 struct isl_mat *isl_mat_copy(struct isl_ctx *ctx, struct isl_mat *mat)
128 {
129 if (!mat)
130 return NULL;
131
132 mat->ref++;
133 return mat;
134 }
135
136 struct isl_mat *isl_mat_dup(struct isl_ctx *ctx, struct isl_mat *mat)
137 {
138 int i;
139 struct isl_mat *mat2;
140
141 if (!mat)
142 return NULL;
143 mat2 = isl_mat_alloc(ctx, mat->n_row, mat->n_col);
144 if (!mat2)
145 return NULL;
146 for (i = 0; i < mat->n_row; ++i)
147 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
148 return mat2;
149 }
150
151 struct isl_mat *isl_mat_cow(struct isl_ctx *ctx, struct isl_mat *mat)
152 {
153 struct isl_mat *mat2;
154 if (!mat)
155 return NULL;
156
157 if (mat->ref == 1 && !ISL_F_ISSET(mat, ISL_MAT_BORROWED))
158 return mat;
159
160 mat2 = isl_mat_dup(ctx, mat);
161 isl_mat_free(ctx, mat);
162 return mat2;
163 }
164
165 void isl_mat_free(struct isl_ctx *ctx, struct isl_mat *mat)
166 {
167 if (!mat)
168 return;
169
170 if (--mat->ref > 0)
171 return;
172
173 if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
174 isl_blk_free(ctx, mat->block);
175 free(mat->row);
176 free(mat);
177 }
178
179 struct isl_mat *isl_mat_identity(struct isl_ctx *ctx, unsigned n_row)
180 {
181 int i;
182 struct isl_mat *mat;
183
184 mat = isl_mat_alloc(ctx, n_row, n_row);
185 if (!mat)
186 return NULL;
187 for (i = 0; i < n_row; ++i) {
188 isl_seq_clr(mat->row[i], i);
189 isl_int_set_si(mat->row[i][i], 1);
190 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
191 }
192
193 return mat;
194 }
195
196 struct isl_vec *isl_mat_vec_product(struct isl_ctx *ctx,
197 struct isl_mat *mat, struct isl_vec *vec)
198 {
199 int i;
200 struct isl_vec *prod;
201
202 if (!mat || !vec)
203 goto error;
204
205 isl_assert(ctx, mat->n_col == vec->size, goto error);
206
207 prod = isl_vec_alloc(ctx, mat->n_row);
208 if (!prod)
209 goto error;
210
211 for (i = 0; i < prod->size; ++i)
212 isl_seq_inner_product(mat->row[i], vec->block.data, vec->size,
213 &prod->block.data[i]);
214 isl_mat_free(ctx, mat);
215 isl_vec_free(ctx, vec);
216 return prod;
217 error:
218 isl_mat_free(ctx, mat);
219 isl_vec_free(ctx, vec);
220 return NULL;
221 }
222
223 struct isl_mat *isl_mat_aff_direct_sum(struct isl_ctx *ctx,
224 struct isl_mat *left, struct isl_mat *right)
225 {
226 int i;
227 struct isl_mat *sum;
228
229 if (!left || !right)
230 goto error;
231
232 isl_assert(ctx, left->n_row == right->n_row, goto error);
233 isl_assert(ctx, left->n_row >= 1, goto error);
234 isl_assert(ctx, left->n_col >= 1, goto error);
235 isl_assert(ctx, right->n_col >= 1, goto error);
236 isl_assert(ctx,
237 isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
238 goto error);
239 isl_assert(ctx,
240 isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
241 goto error);
242
243 sum = isl_mat_alloc(ctx, left->n_row, left->n_col + right->n_col - 1);
244 if (!sum)
245 goto error;
246 isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
247 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
248 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
249
250 isl_seq_clr(sum->row[0]+1, sum->n_col-1);
251 for (i = 1; i < sum->n_row; ++i) {
252 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
253 isl_int_addmul(sum->row[i][0],
254 right->row[0][0], right->row[i][0]);
255 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
256 left->n_col-1);
257 isl_seq_scale(sum->row[i]+left->n_col,
258 right->row[i]+1, right->row[0][0],
259 right->n_col-1);
260 }
261
262 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
263 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
264 isl_mat_free(ctx, left);
265 isl_mat_free(ctx, right);
266 return sum;
267 error:
268 isl_mat_free(ctx, left);
269 isl_mat_free(ctx, right);
270 return NULL;
271 }
272
273 static void exchange(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
274 struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
275 {
276 int r;
277 for (r = row; r < M->n_row; ++r)
278 isl_int_swap(M->row[r][i], M->row[r][j]);
279 if (U) {
280 for (r = 0; r < (*U)->n_row; ++r)
281 isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
282 }
283 if (Q)
284 isl_mat_swap_rows(ctx, *Q, i, j);
285 }
286
287 static void subtract(struct isl_mat *M, struct isl_mat **U,
288 struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
289 {
290 int r;
291 for (r = row; r < M->n_row; ++r)
292 isl_int_submul(M->row[r][j], m, M->row[r][i]);
293 if (U) {
294 for (r = 0; r < (*U)->n_row; ++r)
295 isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
296 }
297 if (Q) {
298 for (r = 0; r < (*Q)->n_col; ++r)
299 isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
300 }
301 }
302
303 static void oppose(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
304 struct isl_mat **Q, unsigned row, unsigned col)
305 {
306 int r;
307 for (r = row; r < M->n_row; ++r)
308 isl_int_neg(M->row[r][col], M->row[r][col]);
309 if (U) {
310 for (r = 0; r < (*U)->n_row; ++r)
311 isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
312 }
313 if (Q)
314 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
315 }
316
317 /* Given matrix M, compute
318 *
319 * M U = H
320 * M = H Q
321 *
322 * with U and Q unimodular matrices and H a matrix in column echelon form
323 * such that on each echelon row the entries in the non-echelon column
324 * are non-negative (if neg == 0) or non-positive (if neg == 1)
325 * and stricly smaller (in absolute value) than the entries in the echelon
326 * column.
327 * If U or Q are NULL, then these matrices are not computed.
328 */
329 struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
330 struct isl_mat *M, int neg, struct isl_mat **U, struct isl_mat **Q)
331 {
332 isl_int c;
333 int row, col;
334
335 if (U)
336 *U = NULL;
337 if (Q)
338 *Q = NULL;
339 if (!M)
340 goto error;
341 M = isl_mat_cow(ctx, M);
342 if (!M)
343 goto error;
344 if (U) {
345 *U = isl_mat_identity(ctx, M->n_col);
346 if (!*U)
347 goto error;
348 }
349 if (Q) {
350 *Q = isl_mat_identity(ctx, M->n_col);
351 if (!*Q)
352 goto error;
353 }
354
355 col = 0;
356 isl_int_init(c);
357 for (row = 0; row < M->n_row; ++row) {
358 int first, i, off;
359 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
360 if (first == -1)
361 continue;
362 first += col;
363 if (first != col)
364 exchange(ctx, M, U, Q, row, first, col);
365 if (isl_int_is_neg(M->row[row][col]))
366 oppose(ctx, M, U, Q, row, col);
367 first = col+1;
368 while ((off = isl_seq_first_non_zero(M->row[row]+first,
369 M->n_col-first)) != -1) {
370 first += off;
371 isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
372 subtract(M, U, Q, row, col, first, c);
373 if (!isl_int_is_zero(M->row[row][first]))
374 exchange(ctx, M, U, Q, row, first, col);
375 else
376 ++first;
377 }
378 for (i = 0; i < col; ++i) {
379 if (isl_int_is_zero(M->row[row][i]))
380 continue;
381 if (neg)
382 isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
383 else
384 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
385 if (isl_int_is_zero(c))
386 continue;
387 subtract(M, U, Q, row, col, i, c);
388 }
389 ++col;
390 }
391 isl_int_clear(c);
392
393 return M;
394 error:
395 if (Q) {
396 isl_mat_free(ctx, *Q);
397 *Q = NULL;
398 }
399 if (U) {
400 isl_mat_free(ctx, *U);
401 *U = NULL;
402 }
403 return NULL;
404 }
405
406 struct isl_mat *isl_mat_right_kernel(struct isl_ctx *ctx, struct isl_mat *mat)
407 {
408 int i, rank;
409 struct isl_mat *U = NULL;
410 struct isl_mat *K;
411
412 mat = isl_mat_left_hermite(ctx, mat, 0, &U, NULL);
413 if (!mat || !U)
414 goto error;
415
416 for (i = 0, rank = 0; rank < mat->n_col; ++rank) {
417 while (i < mat->n_row && isl_int_is_zero(mat->row[i][rank]))
418 ++i;
419 if (i >= mat->n_row)
420 break;
421 }
422 K = isl_mat_alloc(ctx, U->n_row, U->n_col - rank);
423 if (!K)
424 goto error;
425 isl_mat_sub_copy(ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
426 isl_mat_free(ctx, mat);
427 isl_mat_free(ctx, U);
428 return K;
429 error:
430 isl_mat_free(ctx, mat);
431 isl_mat_free(ctx, U);
432 return NULL;
433 }
434
435 struct isl_mat *isl_mat_lin_to_aff(struct isl_ctx *ctx, struct isl_mat *mat)
436 {
437 int i;
438 struct isl_mat *mat2;
439
440 if (!mat)
441 return NULL;
442 mat2 = isl_mat_alloc(ctx, 1+mat->n_row, 1+mat->n_col);
443 if (!mat2)
444 return NULL;
445 isl_int_set_si(mat2->row[0][0], 1);
446 isl_seq_clr(mat2->row[0]+1, mat->n_col);
447 for (i = 0; i < mat->n_row; ++i) {
448 isl_int_set_si(mat2->row[1+i][0], 0);
449 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
450 }
451 isl_mat_free(ctx, mat);
452 return mat2;
453 }
454
455 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
456 {
457 int i;
458
459 for (i = 0; i < n_row; ++i)
460 if (!isl_int_is_zero(row[i][col]))
461 return i;
462 return -1;
463 }
464
465 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
466 {
467 int i, min = row_first_non_zero(row, n_row, col);
468 if (min < 0)
469 return -1;
470 for (i = min + 1; i < n_row; ++i) {
471 if (isl_int_is_zero(row[i][col]))
472 continue;
473 if (isl_int_abs_lt(row[i][col], row[min][col]))
474 min = i;
475 }
476 return min;
477 }
478
479 static void inv_exchange(struct isl_ctx *ctx,
480 struct isl_mat *left, struct isl_mat *right,
481 unsigned i, unsigned j)
482 {
483 left = isl_mat_swap_rows(ctx, left, i, j);
484 right = isl_mat_swap_rows(ctx, right, i, j);
485 }
486
487 static void inv_oppose(struct isl_ctx *ctx,
488 struct isl_mat *left, struct isl_mat *right, unsigned row)
489 {
490 isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
491 isl_seq_neg(right->row[row], right->row[row], right->n_col);
492 }
493
494 static void inv_subtract(struct isl_ctx *ctx,
495 struct isl_mat *left, struct isl_mat *right,
496 unsigned row, unsigned i, isl_int m)
497 {
498 isl_int_neg(m, m);
499 isl_seq_combine(left->row[i]+row,
500 ctx->one, left->row[i]+row,
501 m, left->row[row]+row,
502 left->n_col-row);
503 isl_seq_combine(right->row[i], ctx->one, right->row[i],
504 m, right->row[row], right->n_col);
505 }
506
507 /* Compute inv(left)*right
508 */
509 struct isl_mat *isl_mat_inverse_product(struct isl_ctx *ctx,
510 struct isl_mat *left, struct isl_mat *right)
511 {
512 int row;
513 isl_int a, b;
514
515 if (!left || !right)
516 goto error;
517
518 isl_assert(ctx, left->n_row == left->n_col, goto error);
519 isl_assert(ctx, left->n_row == right->n_row, goto error);
520
521 if (left->n_row == 0) {
522 isl_mat_free(ctx, left);
523 return right;
524 }
525
526 left = isl_mat_cow(ctx, left);
527 right = isl_mat_cow(ctx, right);
528 if (!left || !right)
529 goto error;
530
531 isl_int_init(a);
532 isl_int_init(b);
533 for (row = 0; row < left->n_row; ++row) {
534 int pivot, first, i, off;
535 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
536 if (pivot < 0) {
537 isl_int_clear(a);
538 isl_int_clear(b);
539 isl_assert(ctx, pivot >= 0, goto error);
540 }
541 pivot += row;
542 if (pivot != row)
543 inv_exchange(ctx, left, right, pivot, row);
544 if (isl_int_is_neg(left->row[row][row]))
545 inv_oppose(ctx, left, right, row);
546 first = row+1;
547 while ((off = row_first_non_zero(left->row+first,
548 left->n_row-first, row)) != -1) {
549 first += off;
550 isl_int_fdiv_q(a, left->row[first][row],
551 left->row[row][row]);
552 inv_subtract(ctx, left, right, row, first, a);
553 if (!isl_int_is_zero(left->row[first][row]))
554 inv_exchange(ctx, left, right, row, first);
555 else
556 ++first;
557 }
558 for (i = 0; i < row; ++i) {
559 if (isl_int_is_zero(left->row[i][row]))
560 continue;
561 isl_int_gcd(a, left->row[row][row], left->row[i][row]);
562 isl_int_divexact(b, left->row[i][row], a);
563 isl_int_divexact(a, left->row[row][row], a);
564 isl_int_neg(a, a);
565 isl_seq_combine(left->row[i]+row,
566 a, left->row[i]+row,
567 b, left->row[row]+row,
568 left->n_col-row);
569 isl_seq_combine(right->row[i], a, right->row[i],
570 b, right->row[row], right->n_col);
571 }
572 }
573 isl_int_clear(b);
574
575 isl_int_set(a, left->row[0][0]);
576 for (row = 1; row < left->n_row; ++row)
577 isl_int_lcm(a, a, left->row[row][row]);
578 if (isl_int_is_zero(a)){
579 isl_int_clear(a);
580 isl_assert(ctx, 0, goto error);
581 }
582 for (row = 0; row < left->n_row; ++row) {
583 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
584 if (isl_int_is_one(left->row[row][row]))
585 continue;
586 isl_seq_scale(right->row[row], right->row[row],
587 left->row[row][row], right->n_col);
588 }
589 isl_int_clear(a);
590
591 isl_mat_free(ctx, left);
592 return right;
593 error:
594 isl_mat_free(ctx, left);
595 isl_mat_free(ctx, right);
596 return NULL;
597 }
598
599 void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
600 {
601 int i;
602
603 for (i = 0; i < mat->n_row; ++i)
604 isl_int_mul(mat->row[i][col], mat->row[i][col], m);
605 }
606
607 void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
608 isl_int m1, unsigned src1, isl_int m2, unsigned src2)
609 {
610 int i;
611 isl_int tmp;
612
613 isl_int_init(tmp);
614 for (i = 0; i < mat->n_row; ++i) {
615 isl_int_mul(tmp, m1, mat->row[i][src1]);
616 isl_int_addmul(tmp, m2, mat->row[i][src2]);
617 isl_int_set(mat->row[i][dst], tmp);
618 }
619 isl_int_clear(tmp);
620 }
621
622 struct isl_mat *isl_mat_right_inverse(struct isl_ctx *ctx,
623 struct isl_mat *mat)
624 {
625 struct isl_mat *inv;
626 int row;
627 isl_int a, b;
628
629 mat = isl_mat_cow(ctx, mat);
630 if (!mat)
631 return NULL;
632
633 inv = isl_mat_identity(ctx, mat->n_col);
634 inv = isl_mat_cow(ctx, inv);
635 if (!inv)
636 goto error;
637
638 isl_int_init(a);
639 isl_int_init(b);
640 for (row = 0; row < mat->n_row; ++row) {
641 int pivot, first, i, off;
642 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
643 if (pivot < 0) {
644 isl_int_clear(a);
645 isl_int_clear(b);
646 goto error;
647 }
648 pivot += row;
649 if (pivot != row)
650 exchange(ctx, mat, &inv, NULL, row, pivot, row);
651 if (isl_int_is_neg(mat->row[row][row]))
652 oppose(ctx, mat, &inv, NULL, row, row);
653 first = row+1;
654 while ((off = isl_seq_first_non_zero(mat->row[row]+first,
655 mat->n_col-first)) != -1) {
656 first += off;
657 isl_int_fdiv_q(a, mat->row[row][first],
658 mat->row[row][row]);
659 subtract(mat, &inv, NULL, row, row, first, a);
660 if (!isl_int_is_zero(mat->row[row][first]))
661 exchange(ctx, mat, &inv, NULL, row, row, first);
662 else
663 ++first;
664 }
665 for (i = 0; i < row; ++i) {
666 if (isl_int_is_zero(mat->row[row][i]))
667 continue;
668 isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
669 isl_int_divexact(b, mat->row[row][i], a);
670 isl_int_divexact(a, mat->row[row][row], a);
671 isl_int_neg(a, a);
672 isl_mat_col_combine(mat, i, a, i, b, row);
673 isl_mat_col_combine(inv, i, a, i, b, row);
674 }
675 }
676 isl_int_clear(b);
677
678 isl_int_set(a, mat->row[0][0]);
679 for (row = 1; row < mat->n_row; ++row)
680 isl_int_lcm(a, a, mat->row[row][row]);
681 if (isl_int_is_zero(a)){
682 isl_int_clear(a);
683 goto error;
684 }
685 for (row = 0; row < mat->n_row; ++row) {
686 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
687 if (isl_int_is_one(mat->row[row][row]))
688 continue;
689 isl_mat_col_scale(inv, row, mat->row[row][row]);
690 }
691 isl_int_clear(a);
692
693 isl_mat_free(ctx, mat);
694
695 return inv;
696 error:
697 isl_mat_free(ctx, mat);
698 return NULL;
699 }
700
701 struct isl_mat *isl_mat_transpose(struct isl_ctx *ctx, struct isl_mat *mat)
702 {
703 int i, j;
704
705 mat = isl_mat_cow(ctx, mat);
706 if (!mat)
707 return NULL;
708 isl_assert(ctx, mat->n_col == mat->n_row, goto error);
709 for (i = 0; i < mat->n_row; ++i)
710 for (j = i + 1; j < mat->n_col; ++j)
711 isl_int_swap(mat->row[i][j], mat->row[j][i]);
712 return mat;
713 error:
714 isl_mat_free(ctx, mat);
715 return NULL;
716 }
717
718 struct isl_mat *isl_mat_swap_cols(struct isl_ctx *ctx,
719 struct isl_mat *mat, unsigned i, unsigned j)
720 {
721 int r;
722
723 mat = isl_mat_cow(ctx, mat);
724 if (!mat)
725 return NULL;
726 isl_assert(ctx, i < mat->n_col, goto error);
727 isl_assert(ctx, j < mat->n_col, goto error);
728
729 for (r = 0; r < mat->n_row; ++r)
730 isl_int_swap(mat->row[r][i], mat->row[r][j]);
731 return mat;
732 error:
733 isl_mat_free(ctx, mat);
734 return NULL;
735 }
736
737 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
738 struct isl_mat *mat, unsigned i, unsigned j)
739 {
740 isl_int *t;
741
742 if (!mat)
743 return NULL;
744 mat = isl_mat_cow(ctx, mat);
745 if (!mat)
746 return NULL;
747 t = mat->row[i];
748 mat->row[i] = mat->row[j];
749 mat->row[j] = t;
750 return mat;
751 }
752
753 struct isl_mat *isl_mat_product(struct isl_ctx *ctx,
754 struct isl_mat *left, struct isl_mat *right)
755 {
756 int i, j, k;
757 struct isl_mat *prod;
758
759 if (!left || !right)
760 goto error;
761 isl_assert(ctx, left->n_col == right->n_row, goto error);
762 prod = isl_mat_alloc(ctx, left->n_row, right->n_col);
763 if (!prod)
764 goto error;
765 if (left->n_col == 0) {
766 for (i = 0; i < prod->n_row; ++i)
767 isl_seq_clr(prod->row[i], prod->n_col);
768 return prod;
769 }
770 for (i = 0; i < prod->n_row; ++i) {
771 for (j = 0; j < prod->n_col; ++j) {
772 isl_int_mul(prod->row[i][j],
773 left->row[i][0], right->row[0][j]);
774 for (k = 1; k < left->n_col; ++k)
775 isl_int_addmul(prod->row[i][j],
776 left->row[i][k], right->row[k][j]);
777 }
778 }
779 isl_mat_free(ctx, left);
780 isl_mat_free(ctx, right);
781 return prod;
782 error:
783 isl_mat_free(ctx, left);
784 isl_mat_free(ctx, right);
785 return NULL;
786 }
787
788 /* Replace the variables x in bset by x' given by x = M x', with
789 * M the matrix mat.
790 *
791 * If there are fewer variables x' then there are x, then we perform
792 * the transformation in place, which that, in principle,
793 * this frees up some extra variables as the number
794 * of columns remains constant, but we would have to extend
795 * the div array too as the number of rows in this array is assumed
796 * to be equal to extra.
797 */
798 struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset,
799 struct isl_mat *mat)
800 {
801 struct isl_ctx *ctx;
802 struct isl_mat *t;
803 int i;
804
805 if (!bset || !mat)
806 goto error;
807
808 ctx = bset->ctx;
809 bset = isl_basic_set_cow(bset);
810 if (!bset)
811 goto error;
812
813 isl_assert(ctx, bset->dim->nparam == 0, goto error);
814 isl_assert(ctx, bset->n_div == 0, goto error);
815 isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
816
817 if (mat->n_col > mat->n_row)
818 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0,
819 0, 0);
820 else if (mat->n_col < mat->n_row) {
821 bset->dim = isl_dim_cow(bset->dim);
822 if (!bset->dim)
823 goto error;
824 bset->dim->n_out -= mat->n_row - mat->n_col;
825 }
826
827 t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
828 t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
829 if (!t)
830 goto error;
831 for (i = 0; i < bset->n_eq; ++i) {
832 isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
833 isl_seq_clr(bset->eq[i]+t->n_col, bset->extra);
834 }
835 isl_mat_free(ctx, t);
836
837 t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
838 t = isl_mat_product(ctx, t, mat);
839 if (!t)
840 goto error2;
841 for (i = 0; i < bset->n_ineq; ++i) {
842 isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
843 isl_seq_clr(bset->ineq[i]+t->n_col, bset->extra);
844 }
845 isl_mat_free(ctx, t);
846
847 bset = isl_basic_set_simplify(bset);
848 bset = isl_basic_set_finalize(bset);
849
850 return bset;
851 error:
852 isl_mat_free(ctx, mat);
853 error2:
854 isl_basic_set_free(bset);
855 return NULL;
856 }
857
858 struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat)
859 {
860 struct isl_ctx *ctx;
861 int i;
862
863 set = isl_set_cow(set);
864 if (!set)
865 return NULL;
866
867 ctx = set->ctx;
868 for (i = 0; i < set->n; ++i) {
869 set->p[i] = isl_basic_set_preimage(set->p[i],
870 isl_mat_copy(ctx, mat));
871 if (!set->p[i])
872 goto error;
873 }
874 if (mat->n_col != mat->n_row) {
875 set->dim = isl_dim_cow(set->dim);
876 if (!set->dim)
877 goto error;
878 set->dim->n_out += mat->n_col;
879 set->dim->n_out -= mat->n_row;
880 }
881 isl_mat_free(ctx, mat);
882 ISL_F_CLR(set, ISL_SET_NORMALIZED);
883 return set;
884 error:
885 isl_set_free(set);
886 isl_mat_free(ctx, mat);
887 return NULL;
888 }
889
890 void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
891 FILE *out, int indent)
892 {
893 int i, j;
894
895 if (!mat) {
896 fprintf(out, "%*snull mat\n", indent, "");
897 return;
898 }
899
900 if (mat->n_row == 0)
901 fprintf(out, "%*s[]\n", indent, "");
902
903 for (i = 0; i < mat->n_row; ++i) {
904 if (!i)
905 fprintf(out, "%*s[[", indent, "");
906 else
907 fprintf(out, "%*s[", indent+1, "");
908 for (j = 0; j < mat->n_col; ++j) {
909 if (j)
910 fprintf(out, ",");
911 isl_int_print(out, mat->row[i][j], 0);
912 }
913 if (i == mat->n_row-1)
914 fprintf(out, "]]\n");
915 else
916 fprintf(out, "]\n");
917 }
918 }
919
920 struct isl_mat *isl_mat_drop_cols(struct isl_ctx *ctx, struct isl_mat *mat,
921 unsigned col, unsigned n)
922 {
923 int r;
924
925 mat = isl_mat_cow(ctx, mat);
926 if (!mat)
927 return NULL;
928
929 if (col != mat->n_col-n) {
930 for (r = 0; r < mat->n_row; ++r)
931 isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
932 mat->n_col - col - n);
933 }
934 mat->n_col -= n;
935 return mat;
936 }
937
938 struct isl_mat *isl_mat_drop_rows(struct isl_ctx *ctx, struct isl_mat *mat,
939 unsigned row, unsigned n)
940 {
941 int r;
942
943 mat = isl_mat_cow(ctx, mat);
944 if (!mat)
945 return NULL;
946
947 for (r = row; r+n < mat->n_row; ++r)
948 mat->row[r] = mat->row[r+n];
949
950 mat->n_row -= n;
951 return mat;
952 }
953
954 void isl_mat_col_submul(struct isl_mat *mat,
955 int dst_col, isl_int f, int src_col)
956 {
957 int i;
958
959 for (i = 0; i < mat->n_row; ++i)
960 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
961 }
962
963 void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
964 {
965 int i;
966
967 for (i = 0; i < mat->n_row; ++i)
968 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
969 }
Integer set library