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