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