?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
isl_hash_table_init: take minimal size instead of number of bits needed
[isl.git] / isl_mat.c
blob1230a18717dd297c4918275d7a613abb36995224
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 && !ISL_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 (!ISL_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_right_kernel(struct isl_ctx *ctx, struct isl_mat *mat)
362 {
363 int i, rank;
364 struct isl_mat *U = NULL;
365 struct isl_mat *K;
366
367 mat = isl_mat_left_hermite(ctx, mat, 0, &U, NULL);
368 if (!mat || !U)
369 goto error;
370
371 for (i = 0, rank = 0; rank < mat->n_col; ++rank) {
372 while (i < mat->n_row && isl_int_is_zero(mat->row[i][rank]))
373 ++i;
374 if (i >= mat->n_row)
375 break;
376 }
377 K = isl_mat_alloc(ctx, U->n_row, U->n_col - rank);
378 if (!K)
379 goto error;
380 isl_mat_sub_copy(ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
381 isl_mat_free(ctx, mat);
382 isl_mat_free(ctx, U);
383 return K;
384 error:
385 isl_mat_free(ctx, mat);
386 isl_mat_free(ctx, U);
387 return NULL;
388 }
389
390 struct isl_mat *isl_mat_lin_to_aff(struct isl_ctx *ctx, struct isl_mat *mat)
391 {
392 int i;
393 struct isl_mat *mat2;
394
395 if (!mat)
396 return NULL;
397 mat2 = isl_mat_alloc(ctx, 1+mat->n_row, 1+mat->n_col);
398 if (!mat2)
399 return NULL;
400 isl_int_set_si(mat2->row[0][0], 1);
401 isl_seq_clr(mat2->row[0]+1, mat->n_col);
402 for (i = 0; i < mat->n_row; ++i) {
403 isl_int_set_si(mat2->row[1+i][0], 0);
404 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
405 }
406 isl_mat_free(ctx, mat);
407 return mat2;
408 }
409
410 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
411 {
412 int i;
413
414 for (i = 0; i < n_row; ++i)
415 if (!isl_int_is_zero(row[i][col]))
416 return i;
417 return -1;
418 }
419
420 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
421 {
422 int i, min = row_first_non_zero(row, n_row, col);
423 if (min < 0)
424 return -1;
425 for (i = min + 1; i < n_row; ++i) {
426 if (isl_int_is_zero(row[i][col]))
427 continue;
428 if (isl_int_abs_lt(row[i][col], row[min][col]))
429 min = i;
430 }
431 return min;
432 }
433
434 static void inv_exchange(struct isl_ctx *ctx,
435 struct isl_mat *left, struct isl_mat *right,
436 unsigned i, unsigned j)
437 {
438 left = isl_mat_swap_rows(ctx, left, i, j);
439 right = isl_mat_swap_rows(ctx, right, i, j);
440 }
441
442 static void inv_oppose(struct isl_ctx *ctx,
443 struct isl_mat *left, struct isl_mat *right, unsigned row)
444 {
445 isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
446 isl_seq_neg(right->row[row], right->row[row], right->n_col);
447 }
448
449 static void inv_subtract(struct isl_ctx *ctx,
450 struct isl_mat *left, struct isl_mat *right,
451 unsigned row, unsigned i, isl_int m)
452 {
453 isl_int_neg(m, m);
454 isl_seq_combine(left->row[i]+row,
455 ctx->one, left->row[i]+row,
456 m, left->row[row]+row,
457 left->n_col-row);
458 isl_seq_combine(right->row[i], ctx->one, right->row[i],
459 m, right->row[row], right->n_col);
460 }
461
462 /* Compute inv(left)*right
463 */
464 struct isl_mat *isl_mat_inverse_product(struct isl_ctx *ctx,
465 struct isl_mat *left, struct isl_mat *right)
466 {
467 int row;
468 isl_int a, b;
469
470 if (!left || !right)
471 goto error;
472
473 isl_assert(ctx, left->n_row == left->n_col, goto error);
474 isl_assert(ctx, left->n_row == right->n_row, goto error);
475
476 if (left->n_row == 0) {
477 isl_mat_free(ctx, left);
478 return right;
479 }
480
481 left = isl_mat_cow(ctx, left);
482 right = isl_mat_cow(ctx, right);
483 if (!left || !right)
484 goto error;
485
486 isl_int_init(a);
487 isl_int_init(b);
488 for (row = 0; row < left->n_row; ++row) {
489 int pivot, first, i, off;
490 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
491 if (pivot < 0) {
492 isl_int_clear(a);
493 isl_int_clear(b);
494 isl_assert(ctx, pivot >= 0, goto error);
495 }
496 pivot += row;
497 if (pivot != row)
498 inv_exchange(ctx, left, right, pivot, row);
499 if (isl_int_is_neg(left->row[row][row]))
500 inv_oppose(ctx, left, right, row);
501 first = row+1;
502 while ((off = row_first_non_zero(left->row+first,
503 left->n_row-first, row)) != -1) {
504 first += off;
505 isl_int_fdiv_q(a, left->row[first][row],
506 left->row[row][row]);
507 inv_subtract(ctx, left, right, row, first, a);
508 if (!isl_int_is_zero(left->row[first][row]))
509 inv_exchange(ctx, left, right, row, first);
510 else
511 ++first;
512 }
513 for (i = 0; i < row; ++i) {
514 if (isl_int_is_zero(left->row[i][row]))
515 continue;
516 isl_int_gcd(a, left->row[row][row], left->row[i][row]);
517 isl_int_divexact(b, left->row[i][row], a);
518 isl_int_divexact(a, left->row[row][row], a);
519 isl_int_neg(a, a);
520 isl_seq_combine(left->row[i]+row,
521 a, left->row[i]+row,
522 b, left->row[row]+row,
523 left->n_col-row);
524 isl_seq_combine(right->row[i], a, right->row[i],
525 b, right->row[row], right->n_col);
526 }
527 }
528 isl_int_clear(b);
529
530 isl_int_set(a, left->row[0][0]);
531 for (row = 1; row < left->n_row; ++row)
532 isl_int_lcm(a, a, left->row[row][row]);
533 if (isl_int_is_zero(a)){
534 isl_int_clear(a);
535 isl_assert(ctx, 0, goto error);
536 }
537 for (row = 0; row < left->n_row; ++row) {
538 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
539 if (isl_int_is_one(left->row[row][row]))
540 continue;
541 isl_seq_scale(right->row[row], right->row[row],
542 left->row[row][row], right->n_col);
543 }
544 isl_int_clear(a);
545
546 isl_mat_free(ctx, left);
547 return right;
548 error:
549 isl_mat_free(ctx, left);
550 isl_mat_free(ctx, right);
551 return NULL;
552 }
553
554 void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
555 {
556 int i;
557
558 for (i = 0; i < mat->n_row; ++i)
559 isl_int_mul(mat->row[i][col], mat->row[i][col], m);
560 }
561
562 void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
563 isl_int m1, unsigned src1, isl_int m2, unsigned src2)
564 {
565 int i;
566 isl_int tmp;
567
568 isl_int_init(tmp);
569 for (i = 0; i < mat->n_row; ++i) {
570 isl_int_mul(tmp, m1, mat->row[i][src1]);
571 isl_int_addmul(tmp, m2, mat->row[i][src2]);
572 isl_int_set(mat->row[i][dst], tmp);
573 }
574 isl_int_clear(tmp);
575 }
576
577 struct isl_mat *isl_mat_right_inverse(struct isl_ctx *ctx,
578 struct isl_mat *mat)
579 {
580 struct isl_mat *inv;
581 int row;
582 isl_int a, b;
583
584 mat = isl_mat_cow(ctx, mat);
585 if (!mat)
586 return NULL;
587
588 inv = isl_mat_identity(ctx, mat->n_col);
589 inv = isl_mat_cow(ctx, inv);
590 if (!inv)
591 goto error;
592
593 isl_int_init(a);
594 isl_int_init(b);
595 for (row = 0; row < mat->n_row; ++row) {
596 int pivot, first, i, off;
597 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
598 if (pivot < 0) {
599 isl_int_clear(a);
600 isl_int_clear(b);
601 goto error;
602 }
603 pivot += row;
604 if (pivot != row)
605 exchange(ctx, mat, &inv, NULL, row, pivot, row);
606 if (isl_int_is_neg(mat->row[row][row]))
607 oppose(ctx, mat, &inv, NULL, row, row);
608 first = row+1;
609 while ((off = isl_seq_first_non_zero(mat->row[row]+first,
610 mat->n_col-first)) != -1) {
611 first += off;
612 isl_int_fdiv_q(a, mat->row[row][first],
613 mat->row[row][row]);
614 subtract(mat, &inv, NULL, row, row, first, a);
615 if (!isl_int_is_zero(mat->row[row][first]))
616 exchange(ctx, mat, &inv, NULL, row, row, first);
617 else
618 ++first;
619 }
620 for (i = 0; i < row; ++i) {
621 if (isl_int_is_zero(mat->row[row][i]))
622 continue;
623 isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
624 isl_int_divexact(b, mat->row[row][i], a);
625 isl_int_divexact(a, mat->row[row][row], a);
626 isl_int_neg(a, a);
627 isl_mat_col_combine(mat, i, a, i, b, row);
628 isl_mat_col_combine(inv, i, a, i, b, row);
629 }
630 }
631 isl_int_clear(b);
632
633 isl_int_set(a, mat->row[0][0]);
634 for (row = 1; row < mat->n_row; ++row)
635 isl_int_lcm(a, a, mat->row[row][row]);
636 if (isl_int_is_zero(a)){
637 isl_int_clear(a);
638 goto error;
639 }
640 for (row = 0; row < mat->n_row; ++row) {
641 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
642 if (isl_int_is_one(mat->row[row][row]))
643 continue;
644 isl_mat_col_scale(inv, row, mat->row[row][row]);
645 }
646 isl_int_clear(a);
647
648 isl_mat_free(ctx, mat);
649
650 return inv;
651 error:
652 isl_mat_free(ctx, mat);
653 return NULL;
654 }
655
656 struct isl_mat *isl_mat_transpose(struct isl_ctx *ctx, struct isl_mat *mat)
657 {
658 int i, j;
659
660 mat = isl_mat_cow(ctx, mat);
661 if (!mat)
662 return NULL;
663 isl_assert(ctx, mat->n_col == mat->n_row, goto error);
664 for (i = 0; i < mat->n_row; ++i)
665 for (j = i + 1; j < mat->n_col; ++j)
666 isl_int_swap(mat->row[i][j], mat->row[j][i]);
667 return mat;
668 error:
669 isl_mat_free(ctx, mat);
670 return NULL;
671 }
672
673 struct isl_mat *isl_mat_swap_cols(struct isl_ctx *ctx,
674 struct isl_mat *mat, unsigned i, unsigned j)
675 {
676 int r;
677
678 mat = isl_mat_cow(ctx, mat);
679 if (!mat)
680 return NULL;
681 isl_assert(ctx, i < mat->n_col, goto error);
682 isl_assert(ctx, j < mat->n_col, goto error);
683
684 for (r = 0; r < mat->n_row; ++r)
685 isl_int_swap(mat->row[r][i], mat->row[r][j]);
686 return mat;
687 error:
688 isl_mat_free(ctx, mat);
689 return NULL;
690 }
691
692 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
693 struct isl_mat *mat, unsigned i, unsigned j)
694 {
695 isl_int *t;
696
697 if (!mat)
698 return NULL;
699 mat = isl_mat_cow(ctx, mat);
700 if (!mat)
701 return NULL;
702 t = mat->row[i];
703 mat->row[i] = mat->row[j];
704 mat->row[j] = t;
705 return mat;
706 }
707
708 struct isl_mat *isl_mat_product(struct isl_ctx *ctx,
709 struct isl_mat *left, struct isl_mat *right)
710 {
711 int i, j, k;
712 struct isl_mat *prod;
713
714 if (!left || !right)
715 goto error;
716 isl_assert(ctx, left->n_col == right->n_row, goto error);
717 prod = isl_mat_alloc(ctx, left->n_row, right->n_col);
718 if (!prod)
719 goto error;
720 if (left->n_col == 0) {
721 for (i = 0; i < prod->n_row; ++i)
722 isl_seq_clr(prod->row[i], prod->n_col);
723 return prod;
724 }
725 for (i = 0; i < prod->n_row; ++i) {
726 for (j = 0; j < prod->n_col; ++j) {
727 isl_int_mul(prod->row[i][j],
728 left->row[i][0], right->row[0][j]);
729 for (k = 1; k < left->n_col; ++k)
730 isl_int_addmul(prod->row[i][j],
731 left->row[i][k], right->row[k][j]);
732 }
733 }
734 isl_mat_free(ctx, left);
735 isl_mat_free(ctx, right);
736 return prod;
737 error:
738 isl_mat_free(ctx, left);
739 isl_mat_free(ctx, right);
740 return NULL;
741 }
742
743 /* Replace the variables x in bset by x' given by x = M x', with
744 * M the matrix mat.
745 *
746 * If there are fewer variables x' then there are x, then we perform
747 * the transformation in place, which that, in principle,
748 * this frees up some extra variables as the number
749 * of columns remains constant, but we would have to extend
750 * the div array too as the number of rows in this array is assumed
751 * to be equal to extra.
752 */
753 struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset,
754 struct isl_mat *mat)
755 {
756 struct isl_ctx *ctx;
757 struct isl_mat *t;
758 int i;
759
760 if (!bset || !mat)
761 goto error;
762
763 ctx = bset->ctx;
764 bset = isl_basic_set_cow(bset);
765 if (!bset)
766 goto error;
767
768 isl_assert(ctx, bset->dim->nparam == 0, goto error);
769 isl_assert(ctx, bset->n_div == 0, goto error);
770 isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
771
772 if (mat->n_col > mat->n_row)
773 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0,
774 0, 0);
775 else if (mat->n_col < mat->n_row) {
776 bset->dim = isl_dim_cow(bset->dim);
777 if (!bset->dim)
778 goto error;
779 bset->dim->n_out -= mat->n_row - mat->n_col;
780 }
781
782 t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
783 t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
784 if (!t)
785 goto error;
786 for (i = 0; i < bset->n_eq; ++i) {
787 isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
788 isl_seq_clr(bset->eq[i]+t->n_col, bset->extra);
789 }
790 isl_mat_free(ctx, t);
791
792 t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
793 t = isl_mat_product(ctx, t, mat);
794 if (!t)
795 goto error2;
796 for (i = 0; i < bset->n_ineq; ++i) {
797 isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
798 isl_seq_clr(bset->ineq[i]+t->n_col, bset->extra);
799 }
800 isl_mat_free(ctx, t);
801
802 bset = isl_basic_set_simplify(bset);
803 bset = isl_basic_set_finalize(bset);
804
805 return bset;
806 error:
807 isl_mat_free(ctx, mat);
808 error2:
809 isl_basic_set_free(bset);
810 return NULL;
811 }
812
813 struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat)
814 {
815 struct isl_ctx *ctx;
816 int i;
817
818 set = isl_set_cow(set);
819 if (!set)
820 return NULL;
821
822 ctx = set->ctx;
823 for (i = 0; i < set->n; ++i) {
824 set->p[i] = isl_basic_set_preimage(set->p[i],
825 isl_mat_copy(ctx, mat));
826 if (!set->p[i])
827 goto error;
828 }
829 if (mat->n_col != mat->n_row) {
830 set->dim = isl_dim_cow(set->dim);
831 if (!set->dim)
832 goto error;
833 set->dim->n_out += mat->n_col;
834 set->dim->n_out -= mat->n_row;
835 }
836 isl_mat_free(ctx, mat);
837 ISL_F_CLR(set, ISL_SET_NORMALIZED);
838 return set;
839 error:
840 isl_set_free(set);
841 isl_mat_free(ctx, mat);
842 return NULL;
843 }
844
845 void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
846 FILE *out, int indent)
847 {
848 int i, j;
849
850 if (!mat) {
851 fprintf(out, "%*snull mat\n", indent, "");
852 return;
853 }
854
855 if (mat->n_row == 0)
856 fprintf(out, "%*s[]\n", indent, "");
857
858 for (i = 0; i < mat->n_row; ++i) {
859 if (!i)
860 fprintf(out, "%*s[[", indent, "");
861 else
862 fprintf(out, "%*s[", indent+1, "");
863 for (j = 0; j < mat->n_col; ++j) {
864 if (j)
865 fprintf(out, ",");
866 isl_int_print(out, mat->row[i][j], 0);
867 }
868 if (i == mat->n_row-1)
869 fprintf(out, "]]\n");
870 else
871 fprintf(out, "]\n");
872 }
873 }
874
875 struct isl_mat *isl_mat_drop_cols(struct isl_ctx *ctx, struct isl_mat *mat,
876 unsigned col, unsigned n)
877 {
878 int r;
879
880 mat = isl_mat_cow(ctx, mat);
881 if (!mat)
882 return NULL;
883
884 if (col != mat->n_col-n) {
885 for (r = 0; r < mat->n_row; ++r)
886 isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
887 mat->n_col - col - n);
888 }
889 mat->n_col -= n;
890 return mat;
891 }
892
893 struct isl_mat *isl_mat_drop_rows(struct isl_ctx *ctx, struct isl_mat *mat,
894 unsigned row, unsigned n)
895 {
896 int r;
897
898 mat = isl_mat_cow(ctx, mat);
899 if (!mat)
900 return NULL;
901
902 for (r = row; r+n < mat->n_row; ++r)
903 mat->row[r] = mat->row[r+n];
904
905 mat->n_row -= n;
906 return mat;
907 }
908
909 void isl_mat_col_submul(struct isl_mat *mat,
910 int dst_col, isl_int f, int src_col)
911 {
912 int i;
913
914 for (i = 0; i < mat->n_row; ++i)
915 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
916 }
917
918 void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
919 {
920 int i;
921
922 for (i = 0; i < mat->n_row; ++i)
923 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
924 }
Integer set library