?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
add isl_vec_sort
[isl.git] / isl_mat.c
blob9ccf01ecf057a5bd295910bc416816a4d63da98d
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 *
4 * Use of this software is governed by the GNU LGPLv2.1 license
5 *
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8 */
9
10 #include <isl_ctx_private.h>
11 #include <isl/dim.h>
12 #include <isl/seq.h>
13 #include <isl_mat_private.h>
14 #include "isl_map_private.h"
15 #include <isl_dim_private.h>
16
17 struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
18 unsigned n_row, unsigned n_col)
19 {
20 int i;
21 struct isl_mat *mat;
22
23 mat = isl_alloc_type(ctx, struct isl_mat);
24 if (!mat)
25 return NULL;
26
27 mat->row = NULL;
28 mat->block = isl_blk_alloc(ctx, n_row * n_col);
29 if (isl_blk_is_error(mat->block))
30 goto error;
31 mat->row = isl_alloc_array(ctx, isl_int *, n_row);
32 if (!mat->row)
33 goto error;
34
35 for (i = 0; i < n_row; ++i)
36 mat->row[i] = mat->block.data + i * n_col;
37
38 mat->ctx = ctx;
39 isl_ctx_ref(ctx);
40 mat->ref = 1;
41 mat->n_row = n_row;
42 mat->n_col = n_col;
43 mat->max_col = n_col;
44 mat->flags = 0;
45
46 return mat;
47 error:
48 isl_blk_free(ctx, mat->block);
49 free(mat);
50 return NULL;
51 }
52
53 struct isl_mat *isl_mat_extend(struct isl_mat *mat,
54 unsigned n_row, unsigned n_col)
55 {
56 int i;
57 isl_int *old;
58 isl_int **row;
59
60 if (!mat)
61 return NULL;
62
63 if (mat->max_col >= n_col && mat->n_row >= n_row) {
64 if (mat->n_col < n_col)
65 mat->n_col = n_col;
66 return mat;
67 }
68
69 if (mat->max_col < n_col) {
70 struct isl_mat *new_mat;
71
72 if (n_row < mat->n_row)
73 n_row = mat->n_row;
74 new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
75 if (!new_mat)
76 goto error;
77 for (i = 0; i < mat->n_row; ++i)
78 isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
79 isl_mat_free(mat);
80 return new_mat;
81 }
82
83 mat = isl_mat_cow(mat);
84 if (!mat)
85 goto error;
86
87 old = mat->block.data;
88 mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
89 if (isl_blk_is_error(mat->block))
90 goto error;
91 row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);
92 if (!row)
93 goto error;
94 mat->row = row;
95
96 for (i = 0; i < mat->n_row; ++i)
97 mat->row[i] = mat->block.data + (mat->row[i] - old);
98 for (i = mat->n_row; i < n_row; ++i)
99 mat->row[i] = mat->block.data + i * mat->max_col;
100 mat->n_row = n_row;
101 if (mat->n_col < n_col)
102 mat->n_col = n_col;
103
104 return mat;
105 error:
106 isl_mat_free(mat);
107 return NULL;
108 }
109
110 __isl_give isl_mat *isl_mat_sub_alloc6(isl_ctx *ctx, isl_int **row,
111 unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
112 {
113 int i;
114 struct isl_mat *mat;
115
116 mat = isl_alloc_type(ctx, struct isl_mat);
117 if (!mat)
118 return NULL;
119 mat->row = isl_alloc_array(ctx, isl_int *, n_row);
120 if (!mat->row)
121 goto error;
122 for (i = 0; i < n_row; ++i)
123 mat->row[i] = row[first_row+i] + first_col;
124 mat->ctx = ctx;
125 isl_ctx_ref(ctx);
126 mat->ref = 1;
127 mat->n_row = n_row;
128 mat->n_col = n_col;
129 mat->block = isl_blk_empty();
130 mat->flags = ISL_MAT_BORROWED;
131 return mat;
132 error:
133 free(mat);
134 return NULL;
135 }
136
137 __isl_give isl_mat *isl_mat_sub_alloc(__isl_keep isl_mat *mat,
138 unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
139 {
140 if (!mat)
141 return NULL;
142 return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
143 first_col, n_col);
144 }
145
146 void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
147 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
148 {
149 int i;
150
151 for (i = 0; i < n_row; ++i)
152 isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
153 }
154
155 void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
156 unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
157 {
158 int i;
159
160 for (i = 0; i < n_row; ++i)
161 isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
162 }
163
164 struct isl_mat *isl_mat_copy(struct isl_mat *mat)
165 {
166 if (!mat)
167 return NULL;
168
169 mat->ref++;
170 return mat;
171 }
172
173 struct isl_mat *isl_mat_dup(struct isl_mat *mat)
174 {
175 int i;
176 struct isl_mat *mat2;
177
178 if (!mat)
179 return NULL;
180 mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
181 if (!mat2)
182 return NULL;
183 for (i = 0; i < mat->n_row; ++i)
184 isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
185 return mat2;
186 }
187
188 struct isl_mat *isl_mat_cow(struct isl_mat *mat)
189 {
190 struct isl_mat *mat2;
191 if (!mat)
192 return NULL;
193
194 if (mat->ref == 1 && !ISL_F_ISSET(mat, ISL_MAT_BORROWED))
195 return mat;
196
197 mat2 = isl_mat_dup(mat);
198 isl_mat_free(mat);
199 return mat2;
200 }
201
202 void isl_mat_free(struct isl_mat *mat)
203 {
204 if (!mat)
205 return;
206
207 if (--mat->ref > 0)
208 return;
209
210 if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
211 isl_blk_free(mat->ctx, mat->block);
212 isl_ctx_deref(mat->ctx);
213 free(mat->row);
214 free(mat);
215 }
216
217 int isl_mat_rows(__isl_keep isl_mat *mat)
218 {
219 return mat ? mat->n_row : -1;
220 }
221
222 int isl_mat_cols(__isl_keep isl_mat *mat)
223 {
224 return mat ? mat->n_col : -1;
225 }
226
227 int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
228 {
229 if (!mat)
230 return -1;
231 if (row < 0 || row >= mat->n_row)
232 isl_die(mat->ctx, isl_error_invalid, "row out of range",
233 return -1);
234 if (col < 0 || col >= mat->n_col)
235 isl_die(mat->ctx, isl_error_invalid, "column out of range",
236 return -1);
237 isl_int_set(*v, mat->row[row][col]);
238 return 0;
239 }
240
241 __isl_give isl_mat *isl_mat_set_element(__isl_take isl_mat *mat,
242 int row, int col, isl_int v)
243 {
244 mat = isl_mat_cow(mat);
245 if (!mat)
246 return NULL;
247 if (row < 0 || row >= mat->n_row)
248 isl_die(mat->ctx, isl_error_invalid, "row out of range",
249 goto error);
250 if (col < 0 || col >= mat->n_col)
251 isl_die(mat->ctx, isl_error_invalid, "column out of range",
252 goto error);
253 isl_int_set(mat->row[row][col], v);
254 return mat;
255 error:
256 isl_mat_free(mat);
257 return NULL;
258 }
259
260 __isl_give isl_mat *isl_mat_set_element_si(__isl_take isl_mat *mat,
261 int row, int col, int v)
262 {
263 mat = isl_mat_cow(mat);
264 if (!mat)
265 return NULL;
266 if (row < 0 || row >= mat->n_row)
267 isl_die(mat->ctx, isl_error_invalid, "row out of range",
268 goto error);
269 if (col < 0 || col >= mat->n_col)
270 isl_die(mat->ctx, isl_error_invalid, "column out of range",
271 goto error);
272 isl_int_set_si(mat->row[row][col], v);
273 return mat;
274 error:
275 isl_mat_free(mat);
276 return NULL;
277 }
278
279 struct isl_mat *isl_mat_identity(struct isl_ctx *ctx, unsigned n_row)
280 {
281 int i;
282 struct isl_mat *mat;
283
284 mat = isl_mat_alloc(ctx, n_row, n_row);
285 if (!mat)
286 return NULL;
287 for (i = 0; i < n_row; ++i) {
288 isl_seq_clr(mat->row[i], i);
289 isl_int_set_si(mat->row[i][i], 1);
290 isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
291 }
292
293 return mat;
294 }
295
296 struct isl_vec *isl_mat_vec_product(struct isl_mat *mat, struct isl_vec *vec)
297 {
298 int i;
299 struct isl_vec *prod;
300
301 if (!mat || !vec)
302 goto error;
303
304 isl_assert(mat->ctx, mat->n_col == vec->size, goto error);
305
306 prod = isl_vec_alloc(mat->ctx, mat->n_row);
307 if (!prod)
308 goto error;
309
310 for (i = 0; i < prod->size; ++i)
311 isl_seq_inner_product(mat->row[i], vec->el, vec->size,
312 &prod->block.data[i]);
313 isl_mat_free(mat);
314 isl_vec_free(vec);
315 return prod;
316 error:
317 isl_mat_free(mat);
318 isl_vec_free(vec);
319 return NULL;
320 }
321
322 __isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
323 __isl_take isl_vec *vec)
324 {
325 struct isl_mat *vec_mat;
326 int i;
327
328 if (!mat || !vec)
329 goto error;
330 vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
331 if (!vec_mat)
332 goto error;
333 for (i = 0; i < vec->size; ++i)
334 isl_int_set(vec_mat->row[i][0], vec->el[i]);
335 vec_mat = isl_mat_inverse_product(mat, vec_mat);
336 isl_vec_free(vec);
337 if (!vec_mat)
338 return NULL;
339 vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
340 if (vec)
341 for (i = 0; i < vec->size; ++i)
342 isl_int_set(vec->el[i], vec_mat->row[i][0]);
343 isl_mat_free(vec_mat);
344 return vec;
345 error:
346 isl_mat_free(mat);
347 isl_vec_free(vec);
348 return NULL;
349 }
350
351 struct isl_vec *isl_vec_mat_product(struct isl_vec *vec, struct isl_mat *mat)
352 {
353 int i, j;
354 struct isl_vec *prod;
355
356 if (!mat || !vec)
357 goto error;
358
359 isl_assert(mat->ctx, mat->n_row == vec->size, goto error);
360
361 prod = isl_vec_alloc(mat->ctx, mat->n_col);
362 if (!prod)
363 goto error;
364
365 for (i = 0; i < prod->size; ++i) {
366 isl_int_set_si(prod->el[i], 0);
367 for (j = 0; j < vec->size; ++j)
368 isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
369 }
370 isl_mat_free(mat);
371 isl_vec_free(vec);
372 return prod;
373 error:
374 isl_mat_free(mat);
375 isl_vec_free(vec);
376 return NULL;
377 }
378
379 struct isl_mat *isl_mat_aff_direct_sum(struct isl_mat *left,
380 struct isl_mat *right)
381 {
382 int i;
383 struct isl_mat *sum;
384
385 if (!left || !right)
386 goto error;
387
388 isl_assert(left->ctx, left->n_row == right->n_row, goto error);
389 isl_assert(left->ctx, left->n_row >= 1, goto error);
390 isl_assert(left->ctx, left->n_col >= 1, goto error);
391 isl_assert(left->ctx, right->n_col >= 1, goto error);
392 isl_assert(left->ctx,
393 isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
394 goto error);
395 isl_assert(left->ctx,
396 isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
397 goto error);
398
399 sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
400 if (!sum)
401 goto error;
402 isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
403 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
404 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
405
406 isl_seq_clr(sum->row[0]+1, sum->n_col-1);
407 for (i = 1; i < sum->n_row; ++i) {
408 isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
409 isl_int_addmul(sum->row[i][0],
410 right->row[0][0], right->row[i][0]);
411 isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
412 left->n_col-1);
413 isl_seq_scale(sum->row[i]+left->n_col,
414 right->row[i]+1, right->row[0][0],
415 right->n_col-1);
416 }
417
418 isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
419 isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
420 isl_mat_free(left);
421 isl_mat_free(right);
422 return sum;
423 error:
424 isl_mat_free(left);
425 isl_mat_free(right);
426 return NULL;
427 }
428
429 static void exchange(struct isl_mat *M, struct isl_mat **U,
430 struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
431 {
432 int r;
433 for (r = row; r < M->n_row; ++r)
434 isl_int_swap(M->row[r][i], M->row[r][j]);
435 if (U) {
436 for (r = 0; r < (*U)->n_row; ++r)
437 isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
438 }
439 if (Q)
440 isl_mat_swap_rows(*Q, i, j);
441 }
442
443 static void subtract(struct isl_mat *M, struct isl_mat **U,
444 struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
445 {
446 int r;
447 for (r = row; r < M->n_row; ++r)
448 isl_int_submul(M->row[r][j], m, M->row[r][i]);
449 if (U) {
450 for (r = 0; r < (*U)->n_row; ++r)
451 isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
452 }
453 if (Q) {
454 for (r = 0; r < (*Q)->n_col; ++r)
455 isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
456 }
457 }
458
459 static void oppose(struct isl_mat *M, struct isl_mat **U,
460 struct isl_mat **Q, unsigned row, unsigned col)
461 {
462 int r;
463 for (r = row; r < M->n_row; ++r)
464 isl_int_neg(M->row[r][col], M->row[r][col]);
465 if (U) {
466 for (r = 0; r < (*U)->n_row; ++r)
467 isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
468 }
469 if (Q)
470 isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
471 }
472
473 /* Given matrix M, compute
474 *
475 * M U = H
476 * M = H Q
477 *
478 * with U and Q unimodular matrices and H a matrix in column echelon form
479 * such that on each echelon row the entries in the non-echelon column
480 * are non-negative (if neg == 0) or non-positive (if neg == 1)
481 * and stricly smaller (in absolute value) than the entries in the echelon
482 * column.
483 * If U or Q are NULL, then these matrices are not computed.
484 */
485 struct isl_mat *isl_mat_left_hermite(struct isl_mat *M, int neg,
486 struct isl_mat **U, struct isl_mat **Q)
487 {
488 isl_int c;
489 int row, col;
490
491 if (U)
492 *U = NULL;
493 if (Q)
494 *Q = NULL;
495 if (!M)
496 goto error;
497 M = isl_mat_cow(M);
498 if (!M)
499 goto error;
500 if (U) {
501 *U = isl_mat_identity(M->ctx, M->n_col);
502 if (!*U)
503 goto error;
504 }
505 if (Q) {
506 *Q = isl_mat_identity(M->ctx, M->n_col);
507 if (!*Q)
508 goto error;
509 }
510
511 col = 0;
512 isl_int_init(c);
513 for (row = 0; row < M->n_row; ++row) {
514 int first, i, off;
515 first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
516 if (first == -1)
517 continue;
518 first += col;
519 if (first != col)
520 exchange(M, U, Q, row, first, col);
521 if (isl_int_is_neg(M->row[row][col]))
522 oppose(M, U, Q, row, col);
523 first = col+1;
524 while ((off = isl_seq_first_non_zero(M->row[row]+first,
525 M->n_col-first)) != -1) {
526 first += off;
527 isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
528 subtract(M, U, Q, row, col, first, c);
529 if (!isl_int_is_zero(M->row[row][first]))
530 exchange(M, U, Q, row, first, col);
531 else
532 ++first;
533 }
534 for (i = 0; i < col; ++i) {
535 if (isl_int_is_zero(M->row[row][i]))
536 continue;
537 if (neg)
538 isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
539 else
540 isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
541 if (isl_int_is_zero(c))
542 continue;
543 subtract(M, U, Q, row, col, i, c);
544 }
545 ++col;
546 }
547 isl_int_clear(c);
548
549 return M;
550 error:
551 if (Q) {
552 isl_mat_free(*Q);
553 *Q = NULL;
554 }
555 if (U) {
556 isl_mat_free(*U);
557 *U = NULL;
558 }
559 isl_mat_free(M);
560 return NULL;
561 }
562
563 struct isl_mat *isl_mat_right_kernel(struct isl_mat *mat)
564 {
565 int i, rank;
566 struct isl_mat *U = NULL;
567 struct isl_mat *K;
568
569 mat = isl_mat_left_hermite(mat, 0, &U, NULL);
570 if (!mat || !U)
571 goto error;
572
573 for (i = 0, rank = 0; rank < mat->n_col; ++rank) {
574 while (i < mat->n_row && isl_int_is_zero(mat->row[i][rank]))
575 ++i;
576 if (i >= mat->n_row)
577 break;
578 }
579 K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank);
580 if (!K)
581 goto error;
582 isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
583 isl_mat_free(mat);
584 isl_mat_free(U);
585 return K;
586 error:
587 isl_mat_free(mat);
588 isl_mat_free(U);
589 return NULL;
590 }
591
592 struct isl_mat *isl_mat_lin_to_aff(struct isl_mat *mat)
593 {
594 int i;
595 struct isl_mat *mat2;
596
597 if (!mat)
598 return NULL;
599 mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
600 if (!mat2)
601 goto error;
602 isl_int_set_si(mat2->row[0][0], 1);
603 isl_seq_clr(mat2->row[0]+1, mat->n_col);
604 for (i = 0; i < mat->n_row; ++i) {
605 isl_int_set_si(mat2->row[1+i][0], 0);
606 isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
607 }
608 isl_mat_free(mat);
609 return mat2;
610 error:
611 isl_mat_free(mat);
612 return NULL;
613 }
614
615 /* Given two matrices M1 and M2, return the block matrix
616 *
617 * [ M1 0 ]
618 * [ 0 M2 ]
619 */
620 __isl_give isl_mat *isl_mat_diagonal(__isl_take isl_mat *mat1,
621 __isl_take isl_mat *mat2)
622 {
623 int i;
624 isl_mat *mat;
625
626 if (!mat1 || !mat2)
627 goto error;
628
629 mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
630 mat1->n_col + mat2->n_col);
631 if (!mat)
632 goto error;
633 for (i = 0; i < mat1->n_row; ++i) {
634 isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
635 isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
636 }
637 for (i = 0; i < mat2->n_row; ++i) {
638 isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
639 isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
640 mat2->row[i], mat2->n_col);
641 }
642 isl_mat_free(mat1);
643 isl_mat_free(mat2);
644 return mat;
645 error:
646 isl_mat_free(mat1);
647 isl_mat_free(mat2);
648 return NULL;
649 }
650
651 static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
652 {
653 int i;
654
655 for (i = 0; i < n_row; ++i)
656 if (!isl_int_is_zero(row[i][col]))
657 return i;
658 return -1;
659 }
660
661 static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
662 {
663 int i, min = row_first_non_zero(row, n_row, col);
664 if (min < 0)
665 return -1;
666 for (i = min + 1; i < n_row; ++i) {
667 if (isl_int_is_zero(row[i][col]))
668 continue;
669 if (isl_int_abs_lt(row[i][col], row[min][col]))
670 min = i;
671 }
672 return min;
673 }
674
675 static void inv_exchange(struct isl_mat *left, struct isl_mat *right,
676 unsigned i, unsigned j)
677 {
678 left = isl_mat_swap_rows(left, i, j);
679 right = isl_mat_swap_rows(right, i, j);
680 }
681
682 static void inv_oppose(
683 struct isl_mat *left, struct isl_mat *right, unsigned row)
684 {
685 isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
686 isl_seq_neg(right->row[row], right->row[row], right->n_col);
687 }
688
689 static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
690 unsigned row, unsigned i, isl_int m)
691 {
692 isl_int_neg(m, m);
693 isl_seq_combine(left->row[i]+row,
694 left->ctx->one, left->row[i]+row,
695 m, left->row[row]+row,
696 left->n_col-row);
697 isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
698 m, right->row[row], right->n_col);
699 }
700
701 /* Compute inv(left)*right
702 */
703 struct isl_mat *isl_mat_inverse_product(struct isl_mat *left,
704 struct isl_mat *right)
705 {
706 int row;
707 isl_int a, b;
708
709 if (!left || !right)
710 goto error;
711
712 isl_assert(left->ctx, left->n_row == left->n_col, goto error);
713 isl_assert(left->ctx, left->n_row == right->n_row, goto error);
714
715 if (left->n_row == 0) {
716 isl_mat_free(left);
717 return right;
718 }
719
720 left = isl_mat_cow(left);
721 right = isl_mat_cow(right);
722 if (!left || !right)
723 goto error;
724
725 isl_int_init(a);
726 isl_int_init(b);
727 for (row = 0; row < left->n_row; ++row) {
728 int pivot, first, i, off;
729 pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
730 if (pivot < 0) {
731 isl_int_clear(a);
732 isl_int_clear(b);
733 isl_assert(left->ctx, pivot >= 0, goto error);
734 }
735 pivot += row;
736 if (pivot != row)
737 inv_exchange(left, right, pivot, row);
738 if (isl_int_is_neg(left->row[row][row]))
739 inv_oppose(left, right, row);
740 first = row+1;
741 while ((off = row_first_non_zero(left->row+first,
742 left->n_row-first, row)) != -1) {
743 first += off;
744 isl_int_fdiv_q(a, left->row[first][row],
745 left->row[row][row]);
746 inv_subtract(left, right, row, first, a);
747 if (!isl_int_is_zero(left->row[first][row]))
748 inv_exchange(left, right, row, first);
749 else
750 ++first;
751 }
752 for (i = 0; i < row; ++i) {
753 if (isl_int_is_zero(left->row[i][row]))
754 continue;
755 isl_int_gcd(a, left->row[row][row], left->row[i][row]);
756 isl_int_divexact(b, left->row[i][row], a);
757 isl_int_divexact(a, left->row[row][row], a);
758 isl_int_neg(b, b);
759 isl_seq_combine(left->row[i] + i,
760 a, left->row[i] + i,
761 b, left->row[row] + i,
762 left->n_col - i);
763 isl_seq_combine(right->row[i], a, right->row[i],
764 b, right->row[row], right->n_col);
765 }
766 }
767 isl_int_clear(b);
768
769 isl_int_set(a, left->row[0][0]);
770 for (row = 1; row < left->n_row; ++row)
771 isl_int_lcm(a, a, left->row[row][row]);
772 if (isl_int_is_zero(a)){
773 isl_int_clear(a);
774 isl_assert(left->ctx, 0, goto error);
775 }
776 for (row = 0; row < left->n_row; ++row) {
777 isl_int_divexact(left->row[row][row], a, left->row[row][row]);
778 if (isl_int_is_one(left->row[row][row]))
779 continue;
780 isl_seq_scale(right->row[row], right->row[row],
781 left->row[row][row], right->n_col);
782 }
783 isl_int_clear(a);
784
785 isl_mat_free(left);
786 return right;
787 error:
788 isl_mat_free(left);
789 isl_mat_free(right);
790 return NULL;
791 }
792
793 void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
794 {
795 int i;
796
797 for (i = 0; i < mat->n_row; ++i)
798 isl_int_mul(mat->row[i][col], mat->row[i][col], m);
799 }
800
801 void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
802 isl_int m1, unsigned src1, isl_int m2, unsigned src2)
803 {
804 int i;
805 isl_int tmp;
806
807 isl_int_init(tmp);
808 for (i = 0; i < mat->n_row; ++i) {
809 isl_int_mul(tmp, m1, mat->row[i][src1]);
810 isl_int_addmul(tmp, m2, mat->row[i][src2]);
811 isl_int_set(mat->row[i][dst], tmp);
812 }
813 isl_int_clear(tmp);
814 }
815
816 struct isl_mat *isl_mat_right_inverse(struct isl_mat *mat)
817 {
818 struct isl_mat *inv;
819 int row;
820 isl_int a, b;
821
822 mat = isl_mat_cow(mat);
823 if (!mat)
824 return NULL;
825
826 inv = isl_mat_identity(mat->ctx, mat->n_col);
827 inv = isl_mat_cow(inv);
828 if (!inv)
829 goto error;
830
831 isl_int_init(a);
832 isl_int_init(b);
833 for (row = 0; row < mat->n_row; ++row) {
834 int pivot, first, i, off;
835 pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
836 if (pivot < 0) {
837 isl_int_clear(a);
838 isl_int_clear(b);
839 isl_assert(mat->ctx, pivot >= 0, goto error);
840 }
841 pivot += row;
842 if (pivot != row)
843 exchange(mat, &inv, NULL, row, pivot, row);
844 if (isl_int_is_neg(mat->row[row][row]))
845 oppose(mat, &inv, NULL, row, row);
846 first = row+1;
847 while ((off = isl_seq_first_non_zero(mat->row[row]+first,
848 mat->n_col-first)) != -1) {
849 first += off;
850 isl_int_fdiv_q(a, mat->row[row][first],
851 mat->row[row][row]);
852 subtract(mat, &inv, NULL, row, row, first, a);
853 if (!isl_int_is_zero(mat->row[row][first]))
854 exchange(mat, &inv, NULL, row, row, first);
855 else
856 ++first;
857 }
858 for (i = 0; i < row; ++i) {
859 if (isl_int_is_zero(mat->row[row][i]))
860 continue;
861 isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
862 isl_int_divexact(b, mat->row[row][i], a);
863 isl_int_divexact(a, mat->row[row][row], a);
864 isl_int_neg(a, a);
865 isl_mat_col_combine(mat, i, a, i, b, row);
866 isl_mat_col_combine(inv, i, a, i, b, row);
867 }
868 }
869 isl_int_clear(b);
870
871 isl_int_set(a, mat->row[0][0]);
872 for (row = 1; row < mat->n_row; ++row)
873 isl_int_lcm(a, a, mat->row[row][row]);
874 if (isl_int_is_zero(a)){
875 isl_int_clear(a);
876 goto error;
877 }
878 for (row = 0; row < mat->n_row; ++row) {
879 isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
880 if (isl_int_is_one(mat->row[row][row]))
881 continue;
882 isl_mat_col_scale(inv, row, mat->row[row][row]);
883 }
884 isl_int_clear(a);
885
886 isl_mat_free(mat);
887
888 return inv;
889 error:
890 isl_mat_free(mat);
891 isl_mat_free(inv);
892 return NULL;
893 }
894
895 struct isl_mat *isl_mat_transpose(struct isl_mat *mat)
896 {
897 struct isl_mat *transpose = NULL;
898 int i, j;
899
900 if (mat->n_col == mat->n_row) {
901 mat = isl_mat_cow(mat);
902 if (!mat)
903 return NULL;
904 for (i = 0; i < mat->n_row; ++i)
905 for (j = i + 1; j < mat->n_col; ++j)
906 isl_int_swap(mat->row[i][j], mat->row[j][i]);
907 return mat;
908 }
909 transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);
910 if (!transpose)
911 goto error;
912 for (i = 0; i < mat->n_row; ++i)
913 for (j = 0; j < mat->n_col; ++j)
914 isl_int_set(transpose->row[j][i], mat->row[i][j]);
915 isl_mat_free(mat);
916 return transpose;
917 error:
918 isl_mat_free(mat);
919 return NULL;
920 }
921
922 struct isl_mat *isl_mat_swap_cols(struct isl_mat *mat, unsigned i, unsigned j)
923 {
924 int r;
925
926 mat = isl_mat_cow(mat);
927 if (!mat)
928 return NULL;
929 isl_assert(mat->ctx, i < mat->n_col, goto error);
930 isl_assert(mat->ctx, j < mat->n_col, goto error);
931
932 for (r = 0; r < mat->n_row; ++r)
933 isl_int_swap(mat->row[r][i], mat->row[r][j]);
934 return mat;
935 error:
936 isl_mat_free(mat);
937 return NULL;
938 }
939
940 struct isl_mat *isl_mat_swap_rows(struct isl_mat *mat, unsigned i, unsigned j)
941 {
942 isl_int *t;
943
944 if (!mat)
945 return NULL;
946 mat = isl_mat_cow(mat);
947 if (!mat)
948 return NULL;
949 t = mat->row[i];
950 mat->row[i] = mat->row[j];
951 mat->row[j] = t;
952 return mat;
953 }
954
955 struct isl_mat *isl_mat_product(struct isl_mat *left, struct isl_mat *right)
956 {
957 int i, j, k;
958 struct isl_mat *prod;
959
960 if (!left || !right)
961 goto error;
962 isl_assert(left->ctx, left->n_col == right->n_row, goto error);
963 prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
964 if (!prod)
965 goto error;
966 if (left->n_col == 0) {
967 for (i = 0; i < prod->n_row; ++i)
968 isl_seq_clr(prod->row[i], prod->n_col);
969 isl_mat_free(left);
970 isl_mat_free(right);
971 return prod;
972 }
973 for (i = 0; i < prod->n_row; ++i) {
974 for (j = 0; j < prod->n_col; ++j) {
975 isl_int_mul(prod->row[i][j],
976 left->row[i][0], right->row[0][j]);
977 for (k = 1; k < left->n_col; ++k)
978 isl_int_addmul(prod->row[i][j],
979 left->row[i][k], right->row[k][j]);
980 }
981 }
982 isl_mat_free(left);
983 isl_mat_free(right);
984 return prod;
985 error:
986 isl_mat_free(left);
987 isl_mat_free(right);
988 return NULL;
989 }
990
991 /* Replace the variables x in the rows q by x' given by x = M x',
992 * with M the matrix mat.
993 *
994 * If the number of new variables is greater than the original
995 * number of variables, then the rows q have already been
996 * preextended. If the new number is smaller, then the coefficients
997 * of the divs, which are not changed, need to be shifted down.
998 * The row q may be the equalities, the inequalities or the
999 * div expressions. In the latter case, has_div is true and
1000 * we need to take into account the extra denominator column.
1001 */
1002 static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,
1003 unsigned n_div, int has_div, struct isl_mat *mat)
1004 {
1005 int i;
1006 struct isl_mat *t;
1007 int e;
1008
1009 if (mat->n_col >= mat->n_row)
1010 e = 0;
1011 else
1012 e = mat->n_row - mat->n_col;
1013 if (has_div)
1014 for (i = 0; i < n; ++i)
1015 isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);
1016 t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1017 t = isl_mat_product(t, mat);
1018 if (!t)
1019 return -1;
1020 for (i = 0; i < n; ++i) {
1021 isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1022 isl_seq_cpy(q[i] + has_div + t->n_col,
1023 q[i] + has_div + t->n_col + e, n_div);
1024 isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1025 }
1026 isl_mat_free(t);
1027 return 0;
1028 }
1029
1030 /* Replace the variables x in bset by x' given by x = M x', with
1031 * M the matrix mat.
1032 *
1033 * If there are fewer variables x' then there are x, then we perform
1034 * the transformation in place, which that, in principle,
1035 * this frees up some extra variables as the number
1036 * of columns remains constant, but we would have to extend
1037 * the div array too as the number of rows in this array is assumed
1038 * to be equal to extra.
1039 */
1040 struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset,
1041 struct isl_mat *mat)
1042 {
1043 struct isl_ctx *ctx;
1044
1045 if (!bset || !mat)
1046 goto error;
1047
1048 ctx = bset->ctx;
1049 bset = isl_basic_set_cow(bset);
1050 if (!bset)
1051 goto error;
1052
1053 isl_assert(ctx, bset->dim->nparam == 0, goto error);
1054 isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1055 isl_assert(ctx, mat->n_col > 0, goto error);
1056
1057 if (mat->n_col > mat->n_row) {
1058 bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1059 if (!bset)
1060 goto error;
1061 } else if (mat->n_col < mat->n_row) {
1062 bset->dim = isl_dim_cow(bset->dim);
1063 if (!bset->dim)
1064 goto error;
1065 bset->dim->n_out -= mat->n_row - mat->n_col;
1066 }
1067
1068 if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
1069 isl_mat_copy(mat)) < 0)
1070 goto error;
1071
1072 if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
1073 isl_mat_copy(mat)) < 0)
1074 goto error;
1075
1076 if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)
1077 goto error2;
1078
1079 ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
1080 ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1081 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1082 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1083 ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1084
1085 bset = isl_basic_set_simplify(bset);
1086 bset = isl_basic_set_finalize(bset);
1087
1088 return bset;
1089 error:
1090 isl_mat_free(mat);
1091 error2:
1092 isl_basic_set_free(bset);
1093 return NULL;
1094 }
1095
1096 struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat)
1097 {
1098 struct isl_ctx *ctx;
1099 int i;
1100
1101 set = isl_set_cow(set);
1102 if (!set)
1103 return NULL;
1104
1105 ctx = set->ctx;
1106 for (i = 0; i < set->n; ++i) {
1107 set->p[i] = isl_basic_set_preimage(set->p[i],
1108 isl_mat_copy(mat));
1109 if (!set->p[i])
1110 goto error;
1111 }
1112 if (mat->n_col != mat->n_row) {
1113 set->dim = isl_dim_cow(set->dim);
1114 if (!set->dim)
1115 goto error;
1116 set->dim->n_out += mat->n_col;
1117 set->dim->n_out -= mat->n_row;
1118 }
1119 isl_mat_free(mat);
1120 ISL_F_CLR(set, ISL_SET_NORMALIZED);
1121 return set;
1122 error:
1123 isl_set_free(set);
1124 isl_mat_free(mat);
1125 return NULL;
1126 }
1127
1128 /* Replace the variables x starting at pos in the rows q
1129 * by x' with x = M x' with M the matrix mat.
1130 * That is, replace the corresponding coefficients c by c M.
1131 */
1132 static int transform(isl_ctx *ctx, isl_int **q, unsigned n,
1133 unsigned pos, __isl_take isl_mat *mat)
1134 {
1135 int i;
1136 isl_mat *t;
1137
1138 t = isl_mat_sub_alloc6(ctx, q, 0, n, pos, mat->n_row);
1139 t = isl_mat_product(t, mat);
1140 if (!t)
1141 return -1;
1142 for (i = 0; i < n; ++i)
1143 isl_seq_swp_or_cpy(q[i] + pos, t->row[i], t->n_col);
1144 isl_mat_free(t);
1145 return 0;
1146 }
1147
1148 /* Replace the variables x of type "type" starting at "first" in "bset"
1149 * by x' with x = M x' with M the matrix trans.
1150 * That is, replace the corresponding coefficients c by c M.
1151 *
1152 * The transformation matrix should be a square matrix.
1153 */
1154 __isl_give isl_basic_set *isl_basic_set_transform_dims(
1155 __isl_take isl_basic_set *bset, enum isl_dim_type type, unsigned first,
1156 __isl_take isl_mat *trans)
1157 {
1158 isl_ctx *ctx;
1159 unsigned pos;
1160
1161 bset = isl_basic_set_cow(bset);
1162 if (!bset || !trans)
1163 goto error;
1164
1165 ctx = isl_basic_set_get_ctx(bset);
1166 if (trans->n_row != trans->n_col)
1167 isl_die(trans->ctx, isl_error_invalid,
1168 "expecting square transformation matrix", goto error);
1169 if (first + trans->n_row > isl_basic_set_dim(bset, type))
1170 isl_die(trans->ctx, isl_error_invalid,
1171 "oversized transformation matrix", goto error);
1172
1173 pos = isl_basic_set_offset(bset, type) + first;
1174
1175 if (transform(ctx, bset->eq, bset->n_eq, pos, isl_mat_copy(trans)) < 0)
1176 goto error;
1177 if (transform(ctx, bset->ineq, bset->n_ineq, pos,
1178 isl_mat_copy(trans)) < 0)
1179 goto error;
1180 if (transform(ctx, bset->div, bset->n_div, 1 + pos,
1181 isl_mat_copy(trans)) < 0)
1182 goto error;
1183
1184 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1185 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1186
1187 isl_mat_free(trans);
1188 return bset;
1189 error:
1190 isl_mat_free(trans);
1191 isl_basic_set_free(bset);
1192 return NULL;
1193 }
1194
1195 void isl_mat_print_internal(__isl_keep isl_mat *mat, FILE *out, int indent)
1196 {
1197 int i, j;
1198
1199 if (!mat) {
1200 fprintf(out, "%*snull mat\n", indent, "");
1201 return;
1202 }
1203
1204 if (mat->n_row == 0)
1205 fprintf(out, "%*s[]\n", indent, "");
1206
1207 for (i = 0; i < mat->n_row; ++i) {
1208 if (!i)
1209 fprintf(out, "%*s[[", indent, "");
1210 else
1211 fprintf(out, "%*s[", indent+1, "");
1212 for (j = 0; j < mat->n_col; ++j) {
1213 if (j)
1214 fprintf(out, ",");
1215 isl_int_print(out, mat->row[i][j], 0);
1216 }
1217 if (i == mat->n_row-1)
1218 fprintf(out, "]]\n");
1219 else
1220 fprintf(out, "]\n");
1221 }
1222 }
1223
1224 void isl_mat_dump(__isl_keep isl_mat *mat)
1225 {
1226 isl_mat_print_internal(mat, stderr, 0);
1227 }
1228
1229 struct isl_mat *isl_mat_drop_cols(struct isl_mat *mat, unsigned col, unsigned n)
1230 {
1231 int r;
1232
1233 mat = isl_mat_cow(mat);
1234 if (!mat)
1235 return NULL;
1236
1237 if (col != mat->n_col-n) {
1238 for (r = 0; r < mat->n_row; ++r)
1239 isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1240 mat->n_col - col - n);
1241 }
1242 mat->n_col -= n;
1243 return mat;
1244 }
1245
1246 struct isl_mat *isl_mat_drop_rows(struct isl_mat *mat, unsigned row, unsigned n)
1247 {
1248 int r;
1249
1250 mat = isl_mat_cow(mat);
1251 if (!mat)
1252 return NULL;
1253
1254 for (r = row; r+n < mat->n_row; ++r)
1255 mat->row[r] = mat->row[r+n];
1256
1257 mat->n_row -= n;
1258 return mat;
1259 }
1260
1261 __isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1262 unsigned col, unsigned n)
1263 {
1264 isl_mat *ext;
1265
1266 if (!mat)
1267 return NULL;
1268 if (n == 0)
1269 return mat;
1270
1271 ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1272 if (!ext)
1273 goto error;
1274
1275 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1276 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1277 col + n, col, mat->n_col - col);
1278
1279 isl_mat_free(mat);
1280 return ext;
1281 error:
1282 isl_mat_free(mat);
1283 return NULL;
1284 }
1285
1286 __isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1287 unsigned first, unsigned n)
1288 {
1289 int i;
1290
1291 if (!mat)
1292 return NULL;
1293 mat = isl_mat_insert_cols(mat, first, n);
1294 if (!mat)
1295 return NULL;
1296
1297 for (i = 0; i < mat->n_row; ++i)
1298 isl_seq_clr(mat->row[i] + first, n);
1299
1300 return mat;
1301 }
1302
1303 __isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1304 {
1305 if (!mat)
1306 return NULL;
1307
1308 return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1309 }
1310
1311 __isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1312 unsigned row, unsigned n)
1313 {
1314 isl_mat *ext;
1315
1316 if (!mat)
1317 return NULL;
1318 if (n == 0)
1319 return mat;
1320
1321 ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1322 if (!ext)
1323 goto error;
1324
1325 isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1326 isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1327 mat->n_row - row, 0, 0, mat->n_col);
1328
1329 isl_mat_free(mat);
1330 return ext;
1331 error:
1332 isl_mat_free(mat);
1333 return NULL;
1334 }
1335
1336 __isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1337 {
1338 if (!mat)
1339 return NULL;
1340
1341 return isl_mat_insert_rows(mat, mat->n_row, n);
1342 }
1343
1344 __isl_give isl_mat *isl_mat_insert_zero_rows(__isl_take isl_mat *mat,
1345 unsigned row, unsigned n)
1346 {
1347 int i;
1348
1349 mat = isl_mat_insert_rows(mat, row, n);
1350 if (!mat)
1351 return NULL;
1352
1353 for (i = 0; i < n; ++i)
1354 isl_seq_clr(mat->row[row + i], mat->n_col);
1355
1356 return mat;
1357 }
1358
1359 __isl_give isl_mat *isl_mat_add_zero_rows(__isl_take isl_mat *mat, unsigned n)
1360 {
1361 if (!mat)
1362 return NULL;
1363
1364 return isl_mat_insert_zero_rows(mat, mat->n_row, n);
1365 }
1366
1367 void isl_mat_col_submul(struct isl_mat *mat,
1368 int dst_col, isl_int f, int src_col)
1369 {
1370 int i;
1371
1372 for (i = 0; i < mat->n_row; ++i)
1373 isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1374 }
1375
1376 void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)
1377 {
1378 int i;
1379
1380 if (!mat)
1381 return;
1382
1383 for (i = 0; i < mat->n_row; ++i)
1384 isl_int_add(mat->row[i][dst_col],
1385 mat->row[i][dst_col], mat->row[i][src_col]);
1386 }
1387
1388 void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1389 {
1390 int i;
1391
1392 for (i = 0; i < mat->n_row; ++i)
1393 isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1394 }
1395
1396 struct isl_mat *isl_mat_unimodular_complete(struct isl_mat *M, int row)
1397 {
1398 int r;
1399 struct isl_mat *H = NULL, *Q = NULL;
1400
1401 if (!M)
1402 return NULL;
1403
1404 isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1405 M->n_row = row;
1406 H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1407 M->n_row = M->n_col;
1408 if (!H)
1409 goto error;
1410 for (r = 0; r < row; ++r)
1411 isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1412 for (r = row; r < M->n_row; ++r)
1413 isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1414 isl_mat_free(H);
1415 isl_mat_free(Q);
1416 return M;
1417 error:
1418 isl_mat_free(H);
1419 isl_mat_free(Q);
1420 isl_mat_free(M);
1421 return NULL;
1422 }
1423
1424 __isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1425 __isl_take isl_mat *bot)
1426 {
1427 struct isl_mat *mat;
1428
1429 if (!top || !bot)
1430 goto error;
1431
1432 isl_assert(top->ctx, top->n_col == bot->n_col, goto error);
1433 if (top->n_row == 0) {
1434 isl_mat_free(top);
1435 return bot;
1436 }
1437 if (bot->n_row == 0) {
1438 isl_mat_free(bot);
1439 return top;
1440 }
1441
1442 mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1443 if (!mat)
1444 goto error;
1445 isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1446 0, 0, mat->n_col);
1447 isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1448 0, 0, mat->n_col);
1449 isl_mat_free(top);
1450 isl_mat_free(bot);
1451 return mat;
1452 error:
1453 isl_mat_free(top);
1454 isl_mat_free(bot);
1455 return NULL;
1456 }
1457
1458 int isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1459 {
1460 int i;
1461
1462 if (!mat1 || !mat2)
1463 return -1;
1464
1465 if (mat1->n_row != mat2->n_row)
1466 return 0;
1467
1468 if (mat1->n_col != mat2->n_col)
1469 return 0;
1470
1471 for (i = 0; i < mat1->n_row; ++i)
1472 if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col))
1473 return 0;
1474
1475 return 1;
1476 }
1477
1478 __isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)
1479 {
1480 struct isl_mat *mat;
1481
1482 if (!vec)
1483 return NULL;
1484 mat = isl_mat_alloc(vec->ctx, 1, vec->size);
1485 if (!mat)
1486 goto error;
1487
1488 isl_seq_cpy(mat->row[0], vec->el, vec->size);
1489
1490 isl_vec_free(vec);
1491 return mat;
1492 error:
1493 isl_vec_free(vec);
1494 return NULL;
1495 }
1496
1497 __isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,
1498 __isl_take isl_vec *bot)
1499 {
1500 return isl_mat_concat(top, isl_mat_from_row_vec(bot));
1501 }
1502
1503 __isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,
1504 unsigned dst_col, unsigned src_col, unsigned n)
1505 {
1506 isl_mat *res;
1507
1508 if (!mat)
1509 return NULL;
1510 if (n == 0 || dst_col == src_col)
1511 return mat;
1512
1513 res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1514 if (!res)
1515 goto error;
1516
1517 if (dst_col < src_col) {
1518 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1519 0, 0, dst_col);
1520 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1521 dst_col, src_col, n);
1522 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1523 dst_col + n, dst_col, src_col - dst_col);
1524 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1525 src_col + n, src_col + n,
1526 res->n_col - src_col - n);
1527 } else {
1528 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1529 0, 0, src_col);
1530 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1531 src_col, src_col + n, dst_col - src_col);
1532 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1533 dst_col, src_col, n);
1534 isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1535 dst_col + n, dst_col + n,
1536 res->n_col - dst_col - n);
1537 }
1538 isl_mat_free(mat);
1539
1540 return res;
1541 error:
1542 isl_mat_free(mat);
1543 return NULL;
1544 }
1545
1546 void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1547 {
1548 int i;
1549 isl_int g;
1550
1551 isl_int_set_si(*gcd, 0);
1552 if (!mat)
1553 return;
1554
1555 isl_int_init(g);
1556 for (i = 0; i < mat->n_row; ++i) {
1557 isl_seq_gcd(mat->row[i], mat->n_col, &g);
1558 isl_int_gcd(*gcd, *gcd, g);
1559 }
1560 isl_int_clear(g);
1561 }
1562
1563 __isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1564 {
1565 int i;
1566
1567 if (!mat)
1568 return NULL;
1569
1570 for (i = 0; i < mat->n_row; ++i)
1571 isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1572
1573 return mat;
1574 }
1575
1576 __isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1577 {
1578 isl_int gcd;
1579
1580 if (!mat)
1581 return NULL;
1582
1583 isl_int_init(gcd);
1584 isl_mat_gcd(mat, &gcd);
1585 mat = isl_mat_scale_down(mat, gcd);
1586 isl_int_clear(gcd);
1587
1588 return mat;
1589 }
1590
1591 /* Number of initial non-zero columns.
1592 */
1593 int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
1594 {
1595 int i;
1596
1597 if (!mat)
1598 return -1;
1599
1600 for (i = 0; i < mat->n_col; ++i)
1601 if (row_first_non_zero(mat->row, mat->n_row, i) < 0)
1602 break;
1603
1604 return i;
1605 }
Integer set library