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