1	/*
2	* Copyright 2008-2009 Katholieke Universiteit Leuven
3	* Copyright 2010 INRIA Saclay
4	* Copyright 2014 Ecole Normale Superieure
5	* Copyright 2017 Sven Verdoolaege
6	*
7	* Use of this software is governed by the MIT license
8	*
9	* Written by Sven Verdoolaege, K.U.Leuven, Departement
10	* Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
11	* and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
12	* ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
13	* and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
14	*/
15
16	#include <isl_ctx_private.h>
17	#include <isl_map_private.h>
18	#include <isl/space.h>
19	#include <isl_seq.h>
20	#include <isl_mat_private.h>
21	#include <isl_vec_private.h>
22	#include <isl_space_private.h>
23	#include <isl_val_private.h>
24
25	isl_ctx *isl_mat_get_ctx(__isl_keep isl_mat *mat)
26	{
27	return mat ? mat->ctx : NULL;
28	}
29
30	/* Return a hash value that digests "mat".
31	*/
32	uint32_t isl_mat_get_hash(__isl_keep isl_mat *mat)
33	{
34	int i;
35	uint32_t hash;
36
37	if (!mat)
38	return 0;
39
40	hash = isl_hash_init();
41	isl_hash_byte(hash, mat->n_row & 0xFF);
42	isl_hash_byte(hash, mat->n_col & 0xFF);
43	for (i = 0; i < mat->n_row; ++i) {
44	uint32_t row_hash;
45
46	row_hash = isl_seq_get_hash(p: mat->row[i], len: mat->n_col);
47	isl_hash_hash(hash, row_hash);
48	}
49
50	return hash;
51	}
52
53	__isl_give isl_mat *isl_mat_alloc(isl_ctx *ctx,
54	unsigned n_row, unsigned n_col)
55	{
56	int i;
57	struct isl_mat *mat;
58
59	mat = isl_alloc_type(ctx, struct isl_mat);
60	if (!mat)
61	return NULL;
62
63	mat->row = NULL;
64	mat->block = isl_blk_alloc(ctx, n: n_row * n_col);
65	if (isl_blk_is_error(block: mat->block))
66	goto error;
67	mat->row = isl_calloc_array(ctx, isl_int *, n_row);
68	if (n_row && !mat->row)
69	goto error;
70
71	if (n_col != 0) {
72	for (i = 0; i < n_row; ++i)
73	mat->row[i] = mat->block.data + i * n_col;
74	}
75
76	mat->ctx = ctx;
77	isl_ctx_ref(ctx);
78	mat->ref = 1;
79	mat->n_row = n_row;
80	mat->n_col = n_col;
81	mat->max_col = n_col;
82	mat->flags = 0;
83
84	return mat;
85	error:
86	isl_blk_free(ctx, block: mat->block);
87	free(ptr: mat);
88	return NULL;
89	}
90
91	__isl_give isl_mat *isl_mat_extend(__isl_take isl_mat *mat,
92	unsigned n_row, unsigned n_col)
93	{
94	int i;
95	isl_int *old;
96	isl_int **row;
97
98	if (!mat)
99	return NULL;
100
101	if (mat->max_col >= n_col && mat->n_row >= n_row) {
102	if (mat->n_col < n_col)
103	mat->n_col = n_col;
104	return mat;
105	}
106
107	if (mat->max_col < n_col) {
108	struct isl_mat *new_mat;
109
110	if (n_row < mat->n_row)
111	n_row = mat->n_row;
112	new_mat = isl_mat_alloc(ctx: mat->ctx, n_row, n_col);
113	if (!new_mat)
114	goto error;
115	for (i = 0; i < mat->n_row; ++i)
116	isl_seq_cpy(dst: new_mat->row[i], src: mat->row[i], len: mat->n_col);
117	isl_mat_free(mat);
118	return new_mat;
119	}
120
121	mat = isl_mat_cow(mat);
122	if (!mat)
123	goto error;
124
125	old = mat->block.data;
126	mat->block = isl_blk_extend(ctx: mat->ctx, block: mat->block, new_n: n_row * mat->max_col);
127	if (isl_blk_is_error(block: mat->block))
128	goto error;
129	row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);
130	if (n_row && !row)
131	goto error;
132	mat->row = row;
133
134	for (i = 0; i < mat->n_row; ++i)
135	mat->row[i] = mat->block.data + (mat->row[i] - old);
136	for (i = mat->n_row; i < n_row; ++i)
137	mat->row[i] = mat->block.data + i * mat->max_col;
138	mat->n_row = n_row;
139	if (mat->n_col < n_col)
140	mat->n_col = n_col;
141
142	return mat;
143	error:
144	isl_mat_free(mat);
145	return NULL;
146	}
147
148	__isl_give isl_mat *isl_mat_sub_alloc6(isl_ctx *ctx, isl_int **row,
149	unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
150	{
151	int i;
152	struct isl_mat *mat;
153
154	mat = isl_alloc_type(ctx, struct isl_mat);
155	if (!mat)
156	return NULL;
157	mat->row = isl_alloc_array(ctx, isl_int *, n_row);
158	if (n_row && !mat->row)
159	goto error;
160	for (i = 0; i < n_row; ++i)
161	mat->row[i] = row[first_row+i] + first_col;
162	mat->ctx = ctx;
163	isl_ctx_ref(ctx);
164	mat->ref = 1;
165	mat->n_row = n_row;
166	mat->n_col = n_col;
167	mat->block = isl_blk_empty();
168	mat->flags = ISL_MAT_BORROWED;
169	return mat;
170	error:
171	free(ptr: mat);
172	return NULL;
173	}
174
175	__isl_give isl_mat *isl_mat_sub_alloc(__isl_keep isl_mat *mat,
176	unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
177	{
178	if (!mat)
179	return NULL;
180	return isl_mat_sub_alloc6(ctx: mat->ctx, row: mat->row, first_row, n_row,
181	first_col, n_col);
182	}
183
184	void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
185	unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
186	{
187	int i;
188
189	for (i = 0; i < n_row; ++i)
190	isl_seq_cpy(dst: dst[i]+dst_col, src: src[i]+src_col, len: n_col);
191	}
192
193	void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
194	unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
195	{
196	int i;
197
198	for (i = 0; i < n_row; ++i)
199	isl_seq_neg(dst: dst[i]+dst_col, src: src[i]+src_col, len: n_col);
200	}
201
202	__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)
203	{
204	if (!mat)
205	return NULL;
206
207	mat->ref++;
208	return mat;
209	}
210
211	__isl_give isl_mat *isl_mat_dup(__isl_keep isl_mat *mat)
212	{
213	int i;
214	struct isl_mat *mat2;
215
216	if (!mat)
217	return NULL;
218	mat2 = isl_mat_alloc(ctx: mat->ctx, n_row: mat->n_row, n_col: mat->n_col);
219	if (!mat2)
220	return NULL;
221	for (i = 0; i < mat->n_row; ++i)
222	isl_seq_cpy(dst: mat2->row[i], src: mat->row[i], len: mat->n_col);
223	return mat2;
224	}
225
226	__isl_give isl_mat *isl_mat_cow(__isl_take isl_mat *mat)
227	{
228	struct isl_mat *mat2;
229	if (!mat)
230	return NULL;
231
232	if (mat->ref == 1 && !ISL_F_ISSET(mat, ISL_MAT_BORROWED))
233	return mat;
234
235	mat2 = isl_mat_dup(mat);
236	isl_mat_free(mat);
237	return mat2;
238	}
239
240	__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
241	{
242	if (!mat)
243	return NULL;
244
245	if (--mat->ref > 0)
246	return NULL;
247
248	if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
249	isl_blk_free(ctx: mat->ctx, block: mat->block);
250	isl_ctx_deref(ctx: mat->ctx);
251	free(ptr: mat->row);
252	free(ptr: mat);
253
254	return NULL;
255	}
256
257	isl_size isl_mat_rows(__isl_keep isl_mat *mat)
258	{
259	return mat ? mat->n_row : isl_size_error;
260	}
261
262	isl_size isl_mat_cols(__isl_keep isl_mat *mat)
263	{
264	return mat ? mat->n_col : isl_size_error;
265	}
266
267	/* Check that "col" is a valid column position for "mat".
268	*/
269	static isl_stat check_col(__isl_keep isl_mat *mat, int col)
270	{
271	if (!mat)
272	return isl_stat_error;
273	if (col < 0 || col >= mat->n_col)
274	isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
275	"column out of range", return isl_stat_error);
276	return isl_stat_ok;
277	}
278
279	/* Check that "row" is a valid row position for "mat".
280	*/
281	static isl_stat check_row(__isl_keep isl_mat *mat, int row)
282	{
283	if (!mat)
284	return isl_stat_error;
285	if (row < 0 || row >= mat->n_row)
286	isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
287	"row out of range", return isl_stat_error);
288	return isl_stat_ok;
289	}
290
291	/* Check that there are "n" columns starting at position "first" in "mat".
292	*/
293	static isl_stat check_col_range(__isl_keep isl_mat *mat, unsigned first,
294	unsigned n)
295	{
296	if (!mat)
297	return isl_stat_error;
298	if (first + n > mat->n_col || first + n < first)
299	isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
300	"column position or range out of bounds",
301	return isl_stat_error);
302	return isl_stat_ok;
303	}
304
305	/* Check that there are "n" rows starting at position "first" in "mat".
306	*/
307	static isl_stat check_row_range(__isl_keep isl_mat *mat, unsigned first,
308	unsigned n)
309	{
310	if (!mat)
311	return isl_stat_error;
312	if (first + n > mat->n_row || first + n < first)
313	isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
314	"row position or range out of bounds",
315	return isl_stat_error);
316	return isl_stat_ok;
317	}
318
319	int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
320	{
321	if (check_row(mat, row) < 0)
322	return -1;
323	if (check_col(mat, col) < 0)
324	return -1;
325	isl_int_set(*v, mat->row[row][col]);
326	return 0;
327	}
328
329	/* Extract the element at row "row", oolumn "col" of "mat".
330	*/
331	__isl_give isl_val *isl_mat_get_element_val(__isl_keep isl_mat *mat,
332	int row, int col)
333	{
334	isl_ctx *ctx;
335
336	if (check_row(mat, row) < 0)
337	return NULL;
338	if (check_col(mat, col) < 0)
339	return NULL;
340	ctx = isl_mat_get_ctx(mat);
341	return isl_val_int_from_isl_int(ctx, n: mat->row[row][col]);
342	}
343
344	__isl_give isl_mat *isl_mat_set_element(__isl_take isl_mat *mat,
345	int row, int col, isl_int v)
346	{
347	mat = isl_mat_cow(mat);
348	if (check_row(mat, row) < 0)
349	return isl_mat_free(mat);
350	if (check_col(mat, col) < 0)
351	return isl_mat_free(mat);
352	isl_int_set(mat->row[row][col], v);
353	return mat;
354	}
355
356	__isl_give isl_mat *isl_mat_set_element_si(__isl_take isl_mat *mat,
357	int row, int col, int v)
358	{
359	mat = isl_mat_cow(mat);
360	if (check_row(mat, row) < 0)
361	return isl_mat_free(mat);
362	if (check_col(mat, col) < 0)
363	return isl_mat_free(mat);
364	isl_int_set_si(mat->row[row][col], v);
365	return mat;
366	}
367
368	/* Replace the element at row "row", column "col" of "mat" by "v".
369	*/
370	__isl_give isl_mat *isl_mat_set_element_val(__isl_take isl_mat *mat,
371	int row, int col, __isl_take isl_val *v)
372	{
373	if (!v)
374	return isl_mat_free(mat);
375	if (!isl_val_is_int(v))
376	isl_die(isl_val_get_ctx(v), isl_error_invalid,
377	"expecting integer value", goto error);
378	mat = isl_mat_set_element(mat, row, col, v: v->n);
379	isl_val_free(v);
380	return mat;
381	error:
382	isl_val_free(v);
383	return isl_mat_free(mat);
384	}
385
386	__isl_give isl_mat *isl_mat_diag(isl_ctx *ctx, unsigned n_row, isl_int d)
387	{
388	int i;
389	struct isl_mat *mat;
390
391	mat = isl_mat_alloc(ctx, n_row, n_col: n_row);
392	if (!mat)
393	return NULL;
394	for (i = 0; i < n_row; ++i) {
395	isl_seq_clr(p: mat->row[i], len: i);
396	isl_int_set(mat->row[i][i], d);
397	isl_seq_clr(p: mat->row[i]+i+1, len: n_row-(i+1));
398	}
399
400	return mat;
401	}
402
403	/* Create an "n_row" by "n_col" matrix with zero elements.
404	*/
405	__isl_give isl_mat *isl_mat_zero(isl_ctx *ctx, unsigned n_row, unsigned n_col)
406	{
407	int i;
408	isl_mat *mat;
409
410	mat = isl_mat_alloc(ctx, n_row, n_col);
411	if (!mat)
412	return NULL;
413	for (i = 0; i < n_row; ++i)
414	isl_seq_clr(p: mat->row[i], len: n_col);
415
416	return mat;
417	}
418
419	__isl_give isl_mat *isl_mat_identity(isl_ctx *ctx, unsigned n_row)
420	{
421	if (!ctx)
422	return NULL;
423	return isl_mat_diag(ctx, n_row, d: ctx->one);
424	}
425
426	/* Is "mat" a (possibly scaled) identity matrix?
427	*/
428	isl_bool isl_mat_is_scaled_identity(__isl_keep isl_mat *mat)
429	{
430	int i;
431
432	if (!mat)
433	return isl_bool_error;
434	if (mat->n_row != mat->n_col)
435	return isl_bool_false;
436
437	for (i = 0; i < mat->n_row; ++i) {
438	if (isl_seq_first_non_zero(p: mat->row[i], len: i) != -1)
439	return isl_bool_false;
440	if (isl_int_ne(mat->row[0][0], mat->row[i][i]))
441	return isl_bool_false;
442	if (isl_seq_first_non_zero(p: mat->row[i] + i + 1,
443	len: mat->n_col - (i + 1)) != -1)
444	return isl_bool_false;
445	}
446
447	return isl_bool_true;
448	}
449
450	__isl_give isl_vec *isl_mat_vec_product(__isl_take isl_mat *mat,
451	__isl_take isl_vec *vec)
452	{
453	int i;
454	struct isl_vec *prod;
455
456	if (!mat || !vec)
457	goto error;
458
459	isl_assert(mat->ctx, mat->n_col == vec->size, goto error);
460
461	prod = isl_vec_alloc(ctx: mat->ctx, size: mat->n_row);
462	if (!prod)
463	goto error;
464
465	for (i = 0; i < prod->size; ++i)
466	isl_seq_inner_product(p1: mat->row[i], p2: vec->el, len: vec->size,
467	prod: &prod->block.data[i]);
468	isl_mat_free(mat);
469	isl_vec_free(vec);
470	return prod;
471	error:
472	isl_mat_free(mat);
473	isl_vec_free(vec);
474	return NULL;
475	}
476
477	__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
478	__isl_take isl_vec *vec)
479	{
480	struct isl_mat *vec_mat;
481	int i;
482
483	if (!mat || !vec)
484	goto error;
485	vec_mat = isl_mat_alloc(ctx: vec->ctx, n_row: vec->size, n_col: 1);
486	if (!vec_mat)
487	goto error;
488	for (i = 0; i < vec->size; ++i)
489	isl_int_set(vec_mat->row[i][0], vec->el[i]);
490	vec_mat = isl_mat_inverse_product(left: mat, right: vec_mat);
491	isl_vec_free(vec);
492	if (!vec_mat)
493	return NULL;
494	vec = isl_vec_alloc(ctx: vec_mat->ctx, size: vec_mat->n_row);
495	if (vec)
496	for (i = 0; i < vec->size; ++i)
497	isl_int_set(vec->el[i], vec_mat->row[i][0]);
498	isl_mat_free(mat: vec_mat);
499	return vec;
500	error:
501	isl_mat_free(mat);
502	isl_vec_free(vec);
503	return NULL;
504	}
505
506	__isl_give isl_vec *isl_vec_mat_product(__isl_take isl_vec *vec,
507	__isl_take isl_mat *mat)
508	{
509	int i, j;
510	struct isl_vec *prod;
511
512	if (!mat || !vec)
513	goto error;
514
515	isl_assert(mat->ctx, mat->n_row == vec->size, goto error);
516
517	prod = isl_vec_alloc(ctx: mat->ctx, size: mat->n_col);
518	if (!prod)
519	goto error;
520
521	for (i = 0; i < prod->size; ++i) {
522	isl_int_set_si(prod->el[i], 0);
523	for (j = 0; j < vec->size; ++j)
524	isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
525	}
526	isl_mat_free(mat);
527	isl_vec_free(vec);
528	return prod;
529	error:
530	isl_mat_free(mat);
531	isl_vec_free(vec);
532	return NULL;
533	}
534
535	__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,
536	__isl_take isl_mat *right)
537	{
538	int i;
539	struct isl_mat *sum;
540
541	if (!left || !right)
542	goto error;
543
544	isl_assert(left->ctx, left->n_row == right->n_row, goto error);
545	isl_assert(left->ctx, left->n_row >= 1, goto error);
546	isl_assert(left->ctx, left->n_col >= 1, goto error);
547	isl_assert(left->ctx, right->n_col >= 1, goto error);
548	isl_assert(left->ctx,
549	isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
550	goto error);
551	isl_assert(left->ctx,
552	isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
553	goto error);
554
555	sum = isl_mat_alloc(ctx: left->ctx, n_row: left->n_row, n_col: left->n_col + right->n_col - 1);
556	if (!sum)
557	goto error;
558	isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
559	isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
560	isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
561
562	isl_seq_clr(p: sum->row[0]+1, len: sum->n_col-1);
563	for (i = 1; i < sum->n_row; ++i) {
564	isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
565	isl_int_addmul(sum->row[i][0],
566	right->row[0][0], right->row[i][0]);
567	isl_seq_scale(dst: sum->row[i]+1, src: left->row[i]+1, f: left->row[0][0],
568	len: left->n_col-1);
569	isl_seq_scale(dst: sum->row[i]+left->n_col,
570	src: right->row[i]+1, f: right->row[0][0],
571	len: right->n_col-1);
572	}
573
574	isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
575	isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
576	isl_mat_free(mat: left);
577	isl_mat_free(mat: right);
578	return sum;
579	error:
580	isl_mat_free(mat: left);
581	isl_mat_free(mat: right);
582	return NULL;
583	}
584
585	static void exchange(__isl_keep isl_mat *M, __isl_keep isl_mat **U,
586	__isl_keep isl_mat **Q, unsigned row, unsigned i, unsigned j)
587	{
588	int r;
589	for (r = row; r < M->n_row; ++r)
590	isl_int_swap(M->row[r][i], M->row[r][j]);
591	if (U) {
592	for (r = 0; r < (*U)->n_row; ++r)
593	isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
594	}
595	if (Q)
596	isl_mat_swap_rows(mat: *Q, i, j);
597	}
598
599	static void subtract(__isl_keep isl_mat *M, __isl_keep isl_mat **U,
600	__isl_keep isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
601	{
602	int r;
603	for (r = row; r < M->n_row; ++r)
604	isl_int_submul(M->row[r][j], m, M->row[r][i]);
605	if (U) {
606	for (r = 0; r < (*U)->n_row; ++r)
607	isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
608	}
609	if (Q) {
610	for (r = 0; r < (*Q)->n_col; ++r)
611	isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
612	}
613	}
614
615	static void oppose(__isl_keep isl_mat *M, __isl_keep isl_mat **U,
616	__isl_keep isl_mat **Q, unsigned row, unsigned col)
617	{
618	int r;
619	for (r = row; r < M->n_row; ++r)
620	isl_int_neg(M->row[r][col], M->row[r][col]);
621	if (U) {
622	for (r = 0; r < (*U)->n_row; ++r)
623	isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
624	}
625	if (Q)
626	isl_seq_neg(dst: (*Q)->row[col], src: (*Q)->row[col], len: (*Q)->n_col);
627	}
628
629	/* Given matrix M, compute
630	*
631	* M U = H
632	* M = H Q
633	*
634	* with U and Q unimodular matrices and H a matrix in column echelon form
635	* such that on each echelon row the entries in the non-echelon column
636	* are non-negative (if neg == 0) or non-positive (if neg == 1)
637	* and strictly smaller (in absolute value) than the entries in the echelon
638	* column.
639	* If U or Q are NULL, then these matrices are not computed.
640	*/
641	__isl_give isl_mat *isl_mat_left_hermite(__isl_take isl_mat *M, int neg,
642	__isl_give isl_mat **U, __isl_give isl_mat **Q)
643	{
644	isl_int c;
645	int row, col;
646
647	if (U)
648	*U = NULL;
649	if (Q)
650	*Q = NULL;
651	if (!M)
652	goto error;
653	if (U) {
654	*U = isl_mat_identity(ctx: M->ctx, n_row: M->n_col);
655	if (!*U)
656	goto error;
657	}
658	if (Q) {
659	*Q = isl_mat_identity(ctx: M->ctx, n_row: M->n_col);
660	if (!*Q)
661	goto error;
662	}
663
664	if (M->n_col == 0)
665	return M;
666
667	M = isl_mat_cow(mat: M);
668	if (!M)
669	goto error;
670
671	col = 0;
672	isl_int_init(c);
673	for (row = 0; row < M->n_row; ++row) {
674	int first, i, off;
675	first = isl_seq_abs_min_non_zero(p: M->row[row]+col, len: M->n_col-col);
676	if (first == -1)
677	continue;
678	first += col;
679	if (first != col)
680	exchange(M, U, Q, row, i: first, j: col);
681	if (isl_int_is_neg(M->row[row][col]))
682	oppose(M, U, Q, row, col);
683	first = col+1;
684	while ((off = isl_seq_first_non_zero(p: M->row[row]+first,
685	len: M->n_col-first)) != -1) {
686	first += off;
687	isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
688	subtract(M, U, Q, row, i: col, j: first, m: c);
689	if (!isl_int_is_zero(M->row[row][first]))
690	exchange(M, U, Q, row, i: first, j: col);
691	else
692	++first;
693	}
694	for (i = 0; i < col; ++i) {
695	if (isl_int_is_zero(M->row[row][i]))
696	continue;
697	if (neg)
698	isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
699	else
700	isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
701	if (isl_int_is_zero(c))
702	continue;
703	subtract(M, U, Q, row, i: col, j: i, m: c);
704	}
705	++col;
706	}
707	isl_int_clear(c);
708
709	return M;
710	error:
711	if (Q) {
712	isl_mat_free(mat: *Q);
713	*Q = NULL;
714	}
715	if (U) {
716	isl_mat_free(mat: *U);
717	*U = NULL;
718	}
719	isl_mat_free(mat: M);
720	return NULL;
721	}
722
723	/* Use row "row" of "mat" to eliminate column "col" from all other rows.
724	*/
725	static __isl_give isl_mat *eliminate(__isl_take isl_mat *mat, int row, int col)
726	{
727	int k;
728	isl_size nr, nc;
729	isl_ctx *ctx;
730
731	nr = isl_mat_rows(mat);
732	nc = isl_mat_cols(mat);
733	if (nr < 0 || nc < 0)
734	return isl_mat_free(mat);
735
736	ctx = isl_mat_get_ctx(mat);
737
738	for (k = 0; k < nr; ++k) {
739	if (k == row)
740	continue;
741	if (isl_int_is_zero(mat->row[k][col]))
742	continue;
743	mat = isl_mat_cow(mat);
744	if (!mat)
745	return NULL;
746	isl_seq_elim(dst: mat->row[k], src: mat->row[row], pos: col, len: nc, NULL);
747	isl_seq_normalize(ctx, p: mat->row[k], len: nc);
748	}
749
750	return mat;
751	}
752
753	/* Perform Gaussian elimination on the rows of "mat", but start
754	* from the final row and the final column.
755	* Any zero rows that result from the elimination are removed.
756	*
757	* In particular, for each column from last to first,
758	* look for the last row with a non-zero coefficient in that column,
759	* move it last (but before other rows moved last in previous steps) and
760	* use it to eliminate the column from the other rows.
761	*/
762	__isl_give isl_mat *isl_mat_reverse_gauss(__isl_take isl_mat *mat)
763	{
764	int k, row, last;
765	isl_size nr, nc;
766
767	nr = isl_mat_rows(mat);
768	nc = isl_mat_cols(mat);
769	if (nr < 0 || nc < 0)
770	return isl_mat_free(mat);
771
772	last = nc - 1;
773	for (row = nr - 1; row >= 0; --row) {
774	for (; last >= 0; --last) {
775	for (k = row; k >= 0; --k)
776	if (!isl_int_is_zero(mat->row[k][last]))
777	break;
778	if (k >= 0)
779	break;
780	}
781	if (last < 0)
782	break;
783	if (k != row)
784	mat = isl_mat_swap_rows(mat, i: k, j: row);
785	if (!mat)
786	return NULL;
787	if (isl_int_is_neg(mat->row[row][last]))
788	mat = isl_mat_row_neg(mat, row);
789	mat = eliminate(mat, row, col: last);
790	if (!mat)
791	return NULL;
792	}
793	mat = isl_mat_drop_rows(mat, row: 0, n: row + 1);
794
795	return mat;
796	}
797
798	/* Negate the lexicographically negative rows of "mat" such that
799	* all rows in the result are lexicographically non-negative.
800	*/
801	__isl_give isl_mat *isl_mat_lexnonneg_rows(__isl_take isl_mat *mat)
802	{
803	int i;
804	isl_size nr, nc;
805
806	nr = isl_mat_rows(mat);
807	nc = isl_mat_cols(mat);
808	if (nr < 0 || nc < 0)
809	return isl_mat_free(mat);
810
811	for (i = 0; i < nr; ++i) {
812	int pos;
813
814	pos = isl_seq_first_non_zero(p: mat->row[i], len: nc);
815	if (pos < 0)
816	continue;
817	if (isl_int_is_nonneg(mat->row[i][pos]))
818	continue;
819	mat = isl_mat_row_neg(mat, row: i);
820	if (!mat)
821	return NULL;
822	}
823
824	return mat;
825	}
826
827	/* Given a matrix "H" is column echelon form, what is the first
828	* zero column? That is how many initial columns are non-zero?
829	* Start looking at column "first_col" and only consider
830	* the columns to be of size "n_row".
831	* "H" is assumed to be non-NULL.
832	*
833	* Since "H" is in column echelon form, the first non-zero entry
834	* in a column is always in a later position compared to the previous column.
835	*/
836	static int hermite_first_zero_col(__isl_keep isl_mat *H, int first_col,
837	int n_row)
838	{
839	int row, col;
840
841	for (col = first_col, row = 0; col < H->n_col; ++col) {
842	for (; row < n_row; ++row)
843	if (!isl_int_is_zero(H->row[row][col]))
844	break;
845	if (row == n_row)
846	return col;
847	}
848
849	return H->n_col;
850	}
851
852	/* Return the rank of "mat", or isl_size_error in case of error.
853	*/
854	isl_size isl_mat_rank(__isl_keep isl_mat *mat)
855	{
856	int rank;
857	isl_mat *H;
858
859	H = isl_mat_left_hermite(M: isl_mat_copy(mat), neg: 0, NULL, NULL);
860	if (!H)
861	return isl_size_error;
862
863	rank = hermite_first_zero_col(H, first_col: 0, n_row: H->n_row);
864	isl_mat_free(mat: H);
865
866	return rank;
867	}
868
869	__isl_give isl_mat *isl_mat_right_kernel(__isl_take isl_mat *mat)
870	{
871	int rank;
872	struct isl_mat *U = NULL;
873	struct isl_mat *K;
874
875	mat = isl_mat_left_hermite(M: mat, neg: 0, U: &U, NULL);
876	if (!mat || !U)
877	goto error;
878
879	rank = hermite_first_zero_col(H: mat, first_col: 0, n_row: mat->n_row);
880	K = isl_mat_alloc(ctx: U->ctx, n_row: U->n_row, n_col: U->n_col - rank);
881	if (!K)
882	goto error;
883	isl_mat_sub_copy(ctx: K->ctx, dst: K->row, src: U->row, n_row: U->n_row, dst_col: 0, src_col: rank, n_col: U->n_col-rank);
884	isl_mat_free(mat);
885	isl_mat_free(mat: U);
886	return K;
887	error:
888	isl_mat_free(mat);
889	isl_mat_free(mat: U);
890	return NULL;
891	}
892
893	__isl_give isl_mat *isl_mat_lin_to_aff(__isl_take isl_mat *mat)
894	{
895	int i;
896	struct isl_mat *mat2;
897
898	if (!mat)
899	return NULL;
900	mat2 = isl_mat_alloc(ctx: mat->ctx, n_row: 1+mat->n_row, n_col: 1+mat->n_col);
901	if (!mat2)
902	goto error;
903	isl_int_set_si(mat2->row[0][0], 1);
904	isl_seq_clr(p: mat2->row[0]+1, len: mat->n_col);
905	for (i = 0; i < mat->n_row; ++i) {
906	isl_int_set_si(mat2->row[1+i][0], 0);
907	isl_seq_cpy(dst: mat2->row[1+i]+1, src: mat->row[i], len: mat->n_col);
908	}
909	isl_mat_free(mat);
910	return mat2;
911	error:
912	isl_mat_free(mat);
913	return NULL;
914	}
915
916	/* Given two matrices M1 and M2, return the block matrix
917	*
918	* [ M1 0 ]
919	* [ 0 M2 ]
920	*/
921	__isl_give isl_mat *isl_mat_diagonal(__isl_take isl_mat *mat1,
922	__isl_take isl_mat *mat2)
923	{
924	int i;
925	isl_mat *mat;
926
927	if (!mat1 || !mat2)
928	goto error;
929
930	mat = isl_mat_alloc(ctx: mat1->ctx, n_row: mat1->n_row + mat2->n_row,
931	n_col: mat1->n_col + mat2->n_col);
932	if (!mat)
933	goto error;
934	for (i = 0; i < mat1->n_row; ++i) {
935	isl_seq_cpy(dst: mat->row[i], src: mat1->row[i], len: mat1->n_col);
936	isl_seq_clr(p: mat->row[i] + mat1->n_col, len: mat2->n_col);
937	}
938	for (i = 0; i < mat2->n_row; ++i) {
939	isl_seq_clr(p: mat->row[mat1->n_row + i], len: mat1->n_col);
940	isl_seq_cpy(dst: mat->row[mat1->n_row + i] + mat1->n_col,
941	src: mat2->row[i], len: mat2->n_col);
942	}
943	isl_mat_free(mat: mat1);
944	isl_mat_free(mat: mat2);
945	return mat;
946	error:
947	isl_mat_free(mat: mat1);
948	isl_mat_free(mat: mat2);
949	return NULL;
950	}
951
952	static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
953	{
954	int i;
955
956	for (i = 0; i < n_row; ++i)
957	if (!isl_int_is_zero(row[i][col]))
958	return i;
959	return -1;
960	}
961
962	static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
963	{
964	int i, min = row_first_non_zero(row, n_row, col);
965	if (min < 0)
966	return -1;
967	for (i = min + 1; i < n_row; ++i) {
968	if (isl_int_is_zero(row[i][col]))
969	continue;
970	if (isl_int_abs_lt(row[i][col], row[min][col]))
971	min = i;
972	}
973	return min;
974	}
975
976	static isl_stat inv_exchange(__isl_keep isl_mat **left,
977	__isl_keep isl_mat **right, unsigned i, unsigned j)
978	{
979	*left = isl_mat_swap_rows(mat: *left, i, j);
980	*right = isl_mat_swap_rows(mat: *right, i, j);
981
982	if (!*left || !*right)
983	return isl_stat_error;
984	return isl_stat_ok;
985	}
986
987	static void inv_oppose(
988	__isl_keep isl_mat *left, __isl_keep isl_mat *right, unsigned row)
989	{
990	isl_seq_neg(dst: left->row[row]+row, src: left->row[row]+row, len: left->n_col-row);
991	isl_seq_neg(dst: right->row[row], src: right->row[row], len: right->n_col);
992	}
993
994	static void inv_subtract(__isl_keep isl_mat *left, __isl_keep isl_mat *right,
995	unsigned row, unsigned i, isl_int m)
996	{
997	isl_int_neg(m, m);
998	isl_seq_combine(dst: left->row[i]+row,
999	m1: left->ctx->one, src1: left->row[i]+row,
1000	m2: m, src2: left->row[row]+row,
1001	len: left->n_col-row);
1002	isl_seq_combine(dst: right->row[i], m1: right->ctx->one, src1: right->row[i],
1003	m2: m, src2: right->row[row], len: right->n_col);
1004	}
1005
1006	/* Compute inv(left)*right
1007	*/
1008	__isl_give isl_mat *isl_mat_inverse_product(__isl_take isl_mat *left,
1009	__isl_take isl_mat *right)
1010	{
1011	int row;
1012	isl_int a, b;
1013
1014	if (!left || !right)
1015	goto error;
1016
1017	isl_assert(left->ctx, left->n_row == left->n_col, goto error);
1018	isl_assert(left->ctx, left->n_row == right->n_row, goto error);
1019
1020	if (left->n_row == 0) {
1021	isl_mat_free(mat: left);
1022	return right;
1023	}
1024
1025	left = isl_mat_cow(mat: left);
1026	right = isl_mat_cow(mat: right);
1027	if (!left || !right)
1028	goto error;
1029
1030	isl_int_init(a);
1031	isl_int_init(b);
1032	for (row = 0; row < left->n_row; ++row) {
1033	int pivot, first, i, off;
1034	pivot = row_abs_min_non_zero(row: left->row+row, n_row: left->n_row-row, col: row);
1035	if (pivot < 0) {
1036	isl_int_clear(a);
1037	isl_int_clear(b);
1038	isl_assert(left->ctx, pivot >= 0, goto error);
1039	}
1040	pivot += row;
1041	if (pivot != row)
1042	if (inv_exchange(left: &left, right: &right, i: pivot, j: row) < 0)
1043	goto error;
1044	if (isl_int_is_neg(left->row[row][row]))
1045	inv_oppose(left, right, row);
1046	first = row+1;
1047	while ((off = row_first_non_zero(row: left->row+first,
1048	n_row: left->n_row-first, col: row)) != -1) {
1049	first += off;
1050	isl_int_fdiv_q(a, left->row[first][row],
1051	left->row[row][row]);
1052	inv_subtract(left, right, row, i: first, m: a);
1053	if (!isl_int_is_zero(left->row[first][row])) {
1054	if (inv_exchange(left: &left, right: &right, i: row, j: first) < 0)
1055	goto error;
1056	} else {
1057	++first;
1058	}
1059	}
1060	for (i = 0; i < row; ++i) {
1061	if (isl_int_is_zero(left->row[i][row]))
1062	continue;
1063	isl_int_gcd(a, left->row[row][row], left->row[i][row]);
1064	isl_int_divexact(b, left->row[i][row], a);
1065	isl_int_divexact(a, left->row[row][row], a);
1066	isl_int_neg(b, b);
1067	isl_seq_combine(dst: left->row[i] + i,
1068	m1: a, src1: left->row[i] + i,
1069	m2: b, src2: left->row[row] + i,
1070	len: left->n_col - i);
1071	isl_seq_combine(dst: right->row[i], m1: a, src1: right->row[i],
1072	m2: b, src2: right->row[row], len: right->n_col);
1073	}
1074	}
1075	isl_int_clear(b);
1076
1077	isl_int_set(a, left->row[0][0]);
1078	for (row = 1; row < left->n_row; ++row)
1079	isl_int_lcm(a, a, left->row[row][row]);
1080	if (isl_int_is_zero(a)){
1081	isl_int_clear(a);
1082	isl_assert(left->ctx, 0, goto error);
1083	}
1084	for (row = 0; row < left->n_row; ++row) {
1085	isl_int_divexact(left->row[row][row], a, left->row[row][row]);
1086	if (isl_int_is_one(left->row[row][row]))
1087	continue;
1088	isl_seq_scale(dst: right->row[row], src: right->row[row],
1089	f: left->row[row][row], len: right->n_col);
1090	}
1091	isl_int_clear(a);
1092
1093	isl_mat_free(mat: left);
1094	return right;
1095	error:
1096	isl_mat_free(mat: left);
1097	isl_mat_free(mat: right);
1098	return NULL;
1099	}
1100
1101	void isl_mat_col_scale(__isl_keep isl_mat *mat, unsigned col, isl_int m)
1102	{
1103	int i;
1104
1105	for (i = 0; i < mat->n_row; ++i)
1106	isl_int_mul(mat->row[i][col], mat->row[i][col], m);
1107	}
1108
1109	void isl_mat_col_combine(__isl_keep isl_mat *mat, unsigned dst,
1110	isl_int m1, unsigned src1, isl_int m2, unsigned src2)
1111	{
1112	int i;
1113	isl_int tmp;
1114
1115	isl_int_init(tmp);
1116	for (i = 0; i < mat->n_row; ++i) {
1117	isl_int_mul(tmp, m1, mat->row[i][src1]);
1118	isl_int_addmul(tmp, m2, mat->row[i][src2]);
1119	isl_int_set(mat->row[i][dst], tmp);
1120	}
1121	isl_int_clear(tmp);
1122	}
1123
1124	__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)
1125	{
1126	struct isl_mat *inv;
1127	int row;
1128	isl_int a, b;
1129
1130	mat = isl_mat_cow(mat);
1131	if (!mat)
1132	return NULL;
1133
1134	inv = isl_mat_identity(ctx: mat->ctx, n_row: mat->n_col);
1135	inv = isl_mat_cow(mat: inv);
1136	if (!inv)
1137	goto error;
1138
1139	isl_int_init(a);
1140	isl_int_init(b);
1141	for (row = 0; row < mat->n_row; ++row) {
1142	int pivot, first, i, off;
1143	pivot = isl_seq_abs_min_non_zero(p: mat->row[row]+row, len: mat->n_col-row);
1144	if (pivot < 0) {
1145	isl_int_clear(a);
1146	isl_int_clear(b);
1147	isl_assert(mat->ctx, pivot >= 0, goto error);
1148	}
1149	pivot += row;
1150	if (pivot != row)
1151	exchange(M: mat, U: &inv, NULL, row, i: pivot, j: row);
1152	if (isl_int_is_neg(mat->row[row][row]))
1153	oppose(M: mat, U: &inv, NULL, row, col: row);
1154	first = row+1;
1155	while ((off = isl_seq_first_non_zero(p: mat->row[row]+first,
1156	len: mat->n_col-first)) != -1) {
1157	first += off;
1158	isl_int_fdiv_q(a, mat->row[row][first],
1159	mat->row[row][row]);
1160	subtract(M: mat, U: &inv, NULL, row, i: row, j: first, m: a);
1161	if (!isl_int_is_zero(mat->row[row][first]))
1162	exchange(M: mat, U: &inv, NULL, row, i: row, j: first);
1163	else
1164	++first;
1165	}
1166	for (i = 0; i < row; ++i) {
1167	if (isl_int_is_zero(mat->row[row][i]))
1168	continue;
1169	isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
1170	isl_int_divexact(b, mat->row[row][i], a);
1171	isl_int_divexact(a, mat->row[row][row], a);
1172	isl_int_neg(a, a);
1173	isl_mat_col_combine(mat, dst: i, m1: a, src1: i, m2: b, src2: row);
1174	isl_mat_col_combine(mat: inv, dst: i, m1: a, src1: i, m2: b, src2: row);
1175	}
1176	}
1177	isl_int_clear(b);
1178
1179	isl_int_set(a, mat->row[0][0]);
1180	for (row = 1; row < mat->n_row; ++row)
1181	isl_int_lcm(a, a, mat->row[row][row]);
1182	if (isl_int_is_zero(a)){
1183	isl_int_clear(a);
1184	goto error;
1185	}
1186	for (row = 0; row < mat->n_row; ++row) {
1187	isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1188	if (isl_int_is_one(mat->row[row][row]))
1189	continue;
1190	isl_mat_col_scale(mat: inv, col: row, m: mat->row[row][row]);
1191	}
1192	isl_int_clear(a);
1193
1194	isl_mat_free(mat);
1195
1196	return inv;
1197	error:
1198	isl_mat_free(mat);
1199	isl_mat_free(mat: inv);
1200	return NULL;
1201	}
1202
1203	__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)
1204	{
1205	struct isl_mat *transpose = NULL;
1206	int i, j;
1207
1208	if (!mat)
1209	return NULL;
1210
1211	if (mat->n_col == mat->n_row) {
1212	mat = isl_mat_cow(mat);
1213	if (!mat)
1214	return NULL;
1215	for (i = 0; i < mat->n_row; ++i)
1216	for (j = i + 1; j < mat->n_col; ++j)
1217	isl_int_swap(mat->row[i][j], mat->row[j][i]);
1218	return mat;
1219	}
1220	transpose = isl_mat_alloc(ctx: mat->ctx, n_row: mat->n_col, n_col: mat->n_row);
1221	if (!transpose)
1222	goto error;
1223	for (i = 0; i < mat->n_row; ++i)
1224	for (j = 0; j < mat->n_col; ++j)
1225	isl_int_set(transpose->row[j][i], mat->row[i][j]);
1226	isl_mat_free(mat);
1227	return transpose;
1228	error:
1229	isl_mat_free(mat);
1230	return NULL;
1231	}
1232
1233	__isl_give isl_mat *isl_mat_swap_cols(__isl_take isl_mat *mat,
1234	unsigned i, unsigned j)
1235	{
1236	int r;
1237
1238	mat = isl_mat_cow(mat);
1239	if (check_col_range(mat, first: i, n: 1) < 0 ||
1240	check_col_range(mat, first: j, n: 1) < 0)
1241	return isl_mat_free(mat);
1242
1243	for (r = 0; r < mat->n_row; ++r)
1244	isl_int_swap(mat->row[r][i], mat->row[r][j]);
1245	return mat;
1246	}
1247
1248	__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,
1249	unsigned i, unsigned j)
1250	{
1251	isl_int *t;
1252
1253	if (!mat)
1254	return NULL;
1255	mat = isl_mat_cow(mat);
1256	if (check_row_range(mat, first: i, n: 1) < 0 ||
1257	check_row_range(mat, first: j, n: 1) < 0)
1258	return isl_mat_free(mat);
1259
1260	t = mat->row[i];
1261	mat->row[i] = mat->row[j];
1262	mat->row[j] = t;
1263	return mat;
1264	}
1265
1266	/* Calculate the product of two matrices.
1267	*
1268	* This function is optimized for operand matrices that contain many zeros and
1269	* skips multiplications where we know one of the operands is zero.
1270	*/
1271	__isl_give isl_mat *isl_mat_product(__isl_take isl_mat *left,
1272	__isl_take isl_mat *right)
1273	{
1274	int i, j, k;
1275	struct isl_mat *prod;
1276
1277	if (!left || !right)
1278	goto error;
1279	isl_assert(left->ctx, left->n_col == right->n_row, goto error);
1280	prod = isl_mat_alloc(ctx: left->ctx, n_row: left->n_row, n_col: right->n_col);
1281	if (!prod)
1282	goto error;
1283	if (left->n_col == 0) {
1284	for (i = 0; i < prod->n_row; ++i)
1285	isl_seq_clr(p: prod->row[i], len: prod->n_col);
1286	isl_mat_free(mat: left);
1287	isl_mat_free(mat: right);
1288	return prod;
1289	}
1290	for (i = 0; i < prod->n_row; ++i) {
1291	for (j = 0; j < prod->n_col; ++j)
1292	isl_int_mul(prod->row[i][j],
1293	left->row[i][0], right->row[0][j]);
1294	for (k = 1; k < left->n_col; ++k) {
1295	if (isl_int_is_zero(left->row[i][k]))
1296	continue;
1297	for (j = 0; j < prod->n_col; ++j)
1298	isl_int_addmul(prod->row[i][j],
1299	left->row[i][k], right->row[k][j]);
1300	}
1301	}
1302	isl_mat_free(mat: left);
1303	isl_mat_free(mat: right);
1304	return prod;
1305	error:
1306	isl_mat_free(mat: left);
1307	isl_mat_free(mat: right);
1308	return NULL;
1309	}
1310
1311	/* Replace the variables x in the rows q by x' given by x = M x',
1312	* with M the matrix mat.
1313	*
1314	* If the number of new variables is greater than the original
1315	* number of variables, then the rows q have already been
1316	* preextended. If the new number is smaller, then the coefficients
1317	* of the divs, which are not changed, need to be shifted down.
1318	* The row q may be the equalities, the inequalities or the
1319	* div expressions. In the latter case, has_div is true and
1320	* we need to take into account the extra denominator column.
1321	*/
1322	static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,
1323	unsigned n_div, int has_div, struct isl_mat *mat)
1324	{
1325	int i;
1326	struct isl_mat *t;
1327	int e;
1328
1329	if (mat->n_col >= mat->n_row)
1330	e = 0;
1331	else
1332	e = mat->n_row - mat->n_col;
1333	if (has_div)
1334	for (i = 0; i < n; ++i)
1335	isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);
1336	t = isl_mat_sub_alloc6(ctx: mat->ctx, row: q, first_row: 0, n_row: n, first_col: has_div, n_col: mat->n_row);
1337	t = isl_mat_product(left: t, right: mat);
1338	if (!t)
1339	return -1;
1340	for (i = 0; i < n; ++i) {
1341	isl_seq_swp_or_cpy(dst: q[i] + has_div, src: t->row[i], len: t->n_col);
1342	isl_seq_cpy(dst: q[i] + has_div + t->n_col,
1343	src: q[i] + has_div + t->n_col + e, len: n_div);
1344	isl_seq_clr(p: q[i] + has_div + t->n_col + n_div, len: e);
1345	}
1346	isl_mat_free(mat: t);
1347	return 0;
1348	}
1349
1350	/* Replace the variables x in bset by x' given by x = M x', with
1351	* M the matrix mat.
1352	*
1353	* If there are fewer variables x' then there are x, then we perform
1354	* the transformation in place, which means that, in principle,
1355	* this frees up some extra variables as the number
1356	* of columns remains constant, but we would have to extend
1357	* the div array too as the number of rows in this array is assumed
1358	* to be equal to extra.
1359	*/
1360	__isl_give isl_basic_set *isl_basic_set_preimage(
1361	__isl_take isl_basic_set *bset, __isl_take isl_mat *mat)
1362	{
1363	struct isl_ctx *ctx;
1364
1365	if (!bset || !mat)
1366	goto error;
1367
1368	ctx = bset->ctx;
1369	bset = isl_basic_set_cow(bset);
1370	if (isl_basic_set_check_no_params(bset) < 0)
1371	goto error;
1372
1373	isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1374	isl_assert(ctx, mat->n_col > 0, goto error);
1375
1376	if (mat->n_col > mat->n_row) {
1377	bset = isl_basic_set_add_dims(bset, type: isl_dim_set,
1378	n: mat->n_col - mat->n_row);
1379	if (!bset)
1380	goto error;
1381	} else if (mat->n_col < mat->n_row) {
1382	bset->dim = isl_space_cow(space: bset->dim);
1383	if (!bset->dim)
1384	goto error;
1385	bset->dim->n_out -= mat->n_row - mat->n_col;
1386	}
1387
1388	if (preimage(ctx, q: bset->eq, n: bset->n_eq, n_div: bset->n_div, has_div: 0,
1389	mat: isl_mat_copy(mat)) < 0)
1390	goto error;
1391
1392	if (preimage(ctx, q: bset->ineq, n: bset->n_ineq, n_div: bset->n_div, has_div: 0,
1393	mat: isl_mat_copy(mat)) < 0)
1394	goto error;
1395
1396	if (preimage(ctx, q: bset->div, n: bset->n_div, n_div: bset->n_div, has_div: 1, mat) < 0)
1397	goto error2;
1398
1399	ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
1400	ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1401	ISL_F_CLR(bset, ISL_BASIC_SET_SORTED);
1402	ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1403	ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1404
1405	bset = isl_basic_set_simplify(bset);
1406	bset = isl_basic_set_finalize(bset);
1407
1408	return bset;
1409	error:
1410	isl_mat_free(mat);
1411	error2:
1412	isl_basic_set_free(bset);
1413	return NULL;
1414	}
1415
1416	__isl_give isl_set *isl_set_preimage(
1417	__isl_take isl_set *set, __isl_take isl_mat *mat)
1418	{
1419	int i;
1420
1421	set = isl_set_cow(set);
1422	if (!set)
1423	goto error;
1424
1425	for (i = 0; i < set->n; ++i) {
1426	set->p[i] = isl_basic_set_preimage(bset: set->p[i],
1427	mat: isl_mat_copy(mat));
1428	if (!set->p[i])
1429	goto error;
1430	}
1431	if (mat->n_col != mat->n_row) {
1432	set->dim = isl_space_cow(space: set->dim);
1433	if (!set->dim)
1434	goto error;
1435	set->dim->n_out += mat->n_col;
1436	set->dim->n_out -= mat->n_row;
1437	}
1438	isl_mat_free(mat);
1439	ISL_F_CLR(set, ISL_SET_NORMALIZED);
1440	return set;
1441	error:
1442	isl_set_free(set);
1443	isl_mat_free(mat);
1444	return NULL;
1445	}
1446
1447	/* Replace the variables x starting at "first_col" in the rows "rows"
1448	* of some coefficient matrix by x' with x = M x' with M the matrix mat.
1449	* That is, replace the corresponding coefficients c by c M.
1450	*/
1451	isl_stat isl_mat_sub_transform(isl_int **row, unsigned n_row,
1452	unsigned first_col, __isl_take isl_mat *mat)
1453	{
1454	int i;
1455	isl_ctx *ctx;
1456	isl_mat *t;
1457
1458	if (!mat)
1459	return isl_stat_error;
1460	ctx = isl_mat_get_ctx(mat);
1461	t = isl_mat_sub_alloc6(ctx, row, first_row: 0, n_row, first_col, n_col: mat->n_row);
1462	t = isl_mat_product(left: t, right: mat);
1463	if (!t)
1464	return isl_stat_error;
1465	for (i = 0; i < n_row; ++i)
1466	isl_seq_swp_or_cpy(dst: row[i] + first_col, src: t->row[i], len: t->n_col);
1467	isl_mat_free(mat: t);
1468	return isl_stat_ok;
1469	}
1470
1471	void isl_mat_print_internal(__isl_keep isl_mat *mat, FILE *out, int indent)
1472	{
1473	int i, j;
1474
1475	if (!mat) {
1476	fprintf(stream: out, format: "%*snull mat\n", indent, "");
1477	return;
1478	}
1479
1480	if (mat->n_row == 0)
1481	fprintf(stream: out, format: "%*s[]\n", indent, "");
1482
1483	for (i = 0; i < mat->n_row; ++i) {
1484	if (!i)
1485	fprintf(stream: out, format: "%*s[[", indent, "");
1486	else
1487	fprintf(stream: out, format: "%*s[", indent+1, "");
1488	for (j = 0; j < mat->n_col; ++j) {
1489	if (j)
1490	fprintf(stream: out, format: ",");
1491	isl_int_print(out, mat->row[i][j], 0);
1492	}
1493	if (i == mat->n_row-1)
1494	fprintf(stream: out, format: "]]\n");
1495	else
1496	fprintf(stream: out, format: "]\n");
1497	}
1498	}
1499
1500	void isl_mat_dump(__isl_keep isl_mat *mat)
1501	{
1502	isl_mat_print_internal(mat, stderr, indent: 0);
1503	}
1504
1505	__isl_give isl_mat *isl_mat_drop_cols(__isl_take isl_mat *mat,
1506	unsigned col, unsigned n)
1507	{
1508	int r;
1509
1510	if (n == 0)
1511	return mat;
1512
1513	mat = isl_mat_cow(mat);
1514	if (check_col_range(mat, first: col, n) < 0)
1515	return isl_mat_free(mat);
1516
1517	if (col != mat->n_col-n) {
1518	for (r = 0; r < mat->n_row; ++r)
1519	isl_seq_cpy(dst: mat->row[r]+col, src: mat->row[r]+col+n,
1520	len: mat->n_col - col - n);
1521	}
1522	mat->n_col -= n;
1523	return mat;
1524	}
1525
1526	__isl_give isl_mat *isl_mat_drop_rows(__isl_take isl_mat *mat,
1527	unsigned row, unsigned n)
1528	{
1529	int r;
1530
1531	mat = isl_mat_cow(mat);
1532	if (check_row_range(mat, first: row, n) < 0)
1533	return isl_mat_free(mat);
1534
1535	for (r = row; r+n < mat->n_row; ++r)
1536	mat->row[r] = mat->row[r+n];
1537
1538	mat->n_row -= n;
1539	return mat;
1540	}
1541
1542	__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1543	unsigned col, unsigned n)
1544	{
1545	isl_mat *ext;
1546
1547	if (check_col_range(mat, first: col, n: 0) < 0)
1548	return isl_mat_free(mat);
1549	if (n == 0)
1550	return mat;
1551
1552	ext = isl_mat_alloc(ctx: mat->ctx, n_row: mat->n_row, n_col: mat->n_col + n);
1553	if (!ext)
1554	goto error;
1555
1556	isl_mat_sub_copy(ctx: mat->ctx, dst: ext->row, src: mat->row, n_row: mat->n_row, dst_col: 0, src_col: 0, n_col: col);
1557	isl_mat_sub_copy(ctx: mat->ctx, dst: ext->row, src: mat->row, n_row: mat->n_row,
1558	dst_col: col + n, src_col: col, n_col: mat->n_col - col);
1559
1560	isl_mat_free(mat);
1561	return ext;
1562	error:
1563	isl_mat_free(mat);
1564	return NULL;
1565	}
1566
1567	__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1568	unsigned first, unsigned n)
1569	{
1570	int i;
1571
1572	if (!mat)
1573	return NULL;
1574	mat = isl_mat_insert_cols(mat, col: first, n);
1575	if (!mat)
1576	return NULL;
1577
1578	for (i = 0; i < mat->n_row; ++i)
1579	isl_seq_clr(p: mat->row[i] + first, len: n);
1580
1581	return mat;
1582	}
1583
1584	__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1585	{
1586	if (!mat)
1587	return NULL;
1588
1589	return isl_mat_insert_zero_cols(mat, first: mat->n_col, n);
1590	}
1591
1592	__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1593	unsigned row, unsigned n)
1594	{
1595	isl_mat *ext;
1596
1597	if (check_row_range(mat, first: row, n: 0) < 0)
1598	return isl_mat_free(mat);
1599	if (n == 0)
1600	return mat;
1601
1602	ext = isl_mat_alloc(ctx: mat->ctx, n_row: mat->n_row + n, n_col: mat->n_col);
1603	if (!ext)
1604	goto error;
1605
1606	isl_mat_sub_copy(ctx: mat->ctx, dst: ext->row, src: mat->row, n_row: row, dst_col: 0, src_col: 0, n_col: mat->n_col);
1607	isl_mat_sub_copy(ctx: mat->ctx, dst: ext->row + row + n, src: mat->row + row,
1608	n_row: mat->n_row - row, dst_col: 0, src_col: 0, n_col: mat->n_col);
1609
1610	isl_mat_free(mat);
1611	return ext;
1612	error:
1613	isl_mat_free(mat);
1614	return NULL;
1615	}
1616
1617	__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1618	{
1619	if (!mat)
1620	return NULL;
1621
1622	return isl_mat_insert_rows(mat, row: mat->n_row, n);
1623	}
1624
1625	__isl_give isl_mat *isl_mat_insert_zero_rows(__isl_take isl_mat *mat,
1626	unsigned row, unsigned n)
1627	{
1628	int i;
1629
1630	mat = isl_mat_insert_rows(mat, row, n);
1631	if (!mat)
1632	return NULL;
1633
1634	for (i = 0; i < n; ++i)
1635	isl_seq_clr(p: mat->row[row + i], len: mat->n_col);
1636
1637	return mat;
1638	}
1639
1640	__isl_give isl_mat *isl_mat_add_zero_rows(__isl_take isl_mat *mat, unsigned n)
1641	{
1642	if (!mat)
1643	return NULL;
1644
1645	return isl_mat_insert_zero_rows(mat, row: mat->n_row, n);
1646	}
1647
1648	void isl_mat_col_submul(__isl_keep isl_mat *mat,
1649	int dst_col, isl_int f, int src_col)
1650	{
1651	int i;
1652
1653	for (i = 0; i < mat->n_row; ++i)
1654	isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1655	}
1656
1657	void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)
1658	{
1659	int i;
1660
1661	if (!mat)
1662	return;
1663
1664	for (i = 0; i < mat->n_row; ++i)
1665	isl_int_add(mat->row[i][dst_col],
1666	mat->row[i][dst_col], mat->row[i][src_col]);
1667	}
1668
1669	void isl_mat_col_mul(__isl_keep isl_mat *mat, int dst_col, isl_int f,
1670	int src_col)
1671	{
1672	int i;
1673
1674	for (i = 0; i < mat->n_row; ++i)
1675	isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1676	}
1677
1678	/* Add "f" times column "src_col" to column "dst_col" of "mat" and
1679	* return the result.
1680	*/
1681	__isl_give isl_mat *isl_mat_col_addmul(__isl_take isl_mat *mat, int dst_col,
1682	isl_int f, int src_col)
1683	{
1684	int i;
1685
1686	if (check_col(mat, col: dst_col) < 0 || check_col(mat, col: src_col) < 0)
1687	return isl_mat_free(mat);
1688
1689	for (i = 0; i < mat->n_row; ++i) {
1690	if (isl_int_is_zero(mat->row[i][src_col]))
1691	continue;
1692	mat = isl_mat_cow(mat);
1693	if (!mat)
1694	return NULL;
1695	isl_int_addmul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1696	}
1697
1698	return mat;
1699	}
1700
1701	/* Negate column "col" of "mat" and return the result.
1702	*/
1703	__isl_give isl_mat *isl_mat_col_neg(__isl_take isl_mat *mat, int col)
1704	{
1705	int i;
1706
1707	if (check_col(mat, col) < 0)
1708	return isl_mat_free(mat);
1709
1710	for (i = 0; i < mat->n_row; ++i) {
1711	if (isl_int_is_zero(mat->row[i][col]))
1712	continue;
1713	mat = isl_mat_cow(mat);
1714	if (!mat)
1715	return NULL;
1716	isl_int_neg(mat->row[i][col], mat->row[i][col]);
1717	}
1718
1719	return mat;
1720	}
1721
1722	/* Negate row "row" of "mat" and return the result.
1723	*/
1724	__isl_give isl_mat *isl_mat_row_neg(__isl_take isl_mat *mat, int row)
1725	{
1726	if (check_row(mat, row) < 0)
1727	return isl_mat_free(mat);
1728	if (isl_seq_first_non_zero(p: mat->row[row], len: mat->n_col) == -1)
1729	return mat;
1730	mat = isl_mat_cow(mat);
1731	if (!mat)
1732	return NULL;
1733	isl_seq_neg(dst: mat->row[row], src: mat->row[row], len: mat->n_col);
1734	return mat;
1735	}
1736
1737	__isl_give isl_mat *isl_mat_unimodular_complete(__isl_take isl_mat *M, int row)
1738	{
1739	int r;
1740	struct isl_mat *H = NULL, *Q = NULL;
1741
1742	if (!M)
1743	return NULL;
1744
1745	isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1746	M->n_row = row;
1747	H = isl_mat_left_hermite(M: isl_mat_copy(mat: M), neg: 0, NULL, Q: &Q);
1748	M->n_row = M->n_col;
1749	if (!H)
1750	goto error;
1751	for (r = 0; r < row; ++r)
1752	isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1753	for (r = row; r < M->n_row; ++r)
1754	isl_seq_cpy(dst: M->row[r], src: Q->row[r], len: M->n_col);
1755	isl_mat_free(mat: H);
1756	isl_mat_free(mat: Q);
1757	return M;
1758	error:
1759	isl_mat_free(mat: H);
1760	isl_mat_free(mat: Q);
1761	isl_mat_free(mat: M);
1762	return NULL;
1763	}
1764
1765	__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1766	__isl_take isl_mat *bot)
1767	{
1768	struct isl_mat *mat;
1769
1770	if (!top || !bot)
1771	goto error;
1772
1773	isl_assert(top->ctx, top->n_col == bot->n_col, goto error);
1774	if (top->n_row == 0) {
1775	isl_mat_free(mat: top);
1776	return bot;
1777	}
1778	if (bot->n_row == 0) {
1779	isl_mat_free(mat: bot);
1780	return top;
1781	}
1782
1783	mat = isl_mat_alloc(ctx: top->ctx, n_row: top->n_row + bot->n_row, n_col: top->n_col);
1784	if (!mat)
1785	goto error;
1786	isl_mat_sub_copy(ctx: mat->ctx, dst: mat->row, src: top->row, n_row: top->n_row,
1787	dst_col: 0, src_col: 0, n_col: mat->n_col);
1788	isl_mat_sub_copy(ctx: mat->ctx, dst: mat->row + top->n_row, src: bot->row, n_row: bot->n_row,
1789	dst_col: 0, src_col: 0, n_col: mat->n_col);
1790	isl_mat_free(mat: top);
1791	isl_mat_free(mat: bot);
1792	return mat;
1793	error:
1794	isl_mat_free(mat: top);
1795	isl_mat_free(mat: bot);
1796	return NULL;
1797	}
1798
1799	isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1800	{
1801	int i;
1802
1803	if (!mat1 || !mat2)
1804	return isl_bool_error;
1805
1806	if (mat1->n_row != mat2->n_row)
1807	return isl_bool_false;
1808
1809	if (mat1->n_col != mat2->n_col)
1810	return isl_bool_false;
1811
1812	for (i = 0; i < mat1->n_row; ++i)
1813	if (!isl_seq_eq(p1: mat1->row[i], p2: mat2->row[i], len: mat1->n_col))
1814	return isl_bool_false;
1815
1816	return isl_bool_true;
1817	}
1818
1819	__isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)
1820	{
1821	struct isl_mat *mat;
1822
1823	if (!vec)
1824	return NULL;
1825	mat = isl_mat_alloc(ctx: vec->ctx, n_row: 1, n_col: vec->size);
1826	if (!mat)
1827	goto error;
1828
1829	isl_seq_cpy(dst: mat->row[0], src: vec->el, len: vec->size);
1830
1831	isl_vec_free(vec);
1832	return mat;
1833	error:
1834	isl_vec_free(vec);
1835	return NULL;
1836	}
1837
1838	/* Return a copy of row "row" of "mat" as an isl_vec.
1839	*/
1840	__isl_give isl_vec *isl_mat_get_row(__isl_keep isl_mat *mat, unsigned row)
1841	{
1842	isl_vec *v;
1843
1844	if (!mat)
1845	return NULL;
1846	if (row >= mat->n_row)
1847	isl_die(mat->ctx, isl_error_invalid, "row out of range",
1848	return NULL);
1849
1850	v = isl_vec_alloc(ctx: isl_mat_get_ctx(mat), size: mat->n_col);
1851	if (!v)
1852	return NULL;
1853	isl_seq_cpy(dst: v->el, src: mat->row[row], len: mat->n_col);
1854
1855	return v;
1856	}
1857
1858	__isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,
1859	__isl_take isl_vec *bot)
1860	{
1861	return isl_mat_concat(top, bot: isl_mat_from_row_vec(vec: bot));
1862	}
1863
1864	__isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,
1865	unsigned dst_col, unsigned src_col, unsigned n)
1866	{
1867	isl_mat *res;
1868
1869	if (!mat)
1870	return NULL;
1871	if (n == 0 || dst_col == src_col)
1872	return mat;
1873
1874	res = isl_mat_alloc(ctx: mat->ctx, n_row: mat->n_row, n_col: mat->n_col);
1875	if (!res)
1876	goto error;
1877
1878	if (dst_col < src_col) {
1879	isl_mat_sub_copy(ctx: res->ctx, dst: res->row, src: mat->row, n_row: mat->n_row,
1880	dst_col: 0, src_col: 0, n_col: dst_col);
1881	isl_mat_sub_copy(ctx: res->ctx, dst: res->row, src: mat->row, n_row: mat->n_row,
1882	dst_col, src_col, n_col: n);
1883	isl_mat_sub_copy(ctx: res->ctx, dst: res->row, src: mat->row, n_row: mat->n_row,
1884	dst_col: dst_col + n, src_col: dst_col, n_col: src_col - dst_col);
1885	isl_mat_sub_copy(ctx: res->ctx, dst: res->row, src: mat->row, n_row: mat->n_row,
1886	dst_col: src_col + n, src_col: src_col + n,
1887	n_col: res->n_col - src_col - n);
1888	} else {
1889	isl_mat_sub_copy(ctx: res->ctx, dst: res->row, src: mat->row, n_row: mat->n_row,
1890	dst_col: 0, src_col: 0, n_col: src_col);
1891	isl_mat_sub_copy(ctx: res->ctx, dst: res->row, src: mat->row, n_row: mat->n_row,
1892	dst_col: src_col, src_col: src_col + n, n_col: dst_col - src_col);
1893	isl_mat_sub_copy(ctx: res->ctx, dst: res->row, src: mat->row, n_row: mat->n_row,
1894	dst_col, src_col, n_col: n);
1895	isl_mat_sub_copy(ctx: res->ctx, dst: res->row, src: mat->row, n_row: mat->n_row,
1896	dst_col: dst_col + n, src_col: dst_col + n,
1897	n_col: res->n_col - dst_col - n);
1898	}
1899	isl_mat_free(mat);
1900
1901	return res;
1902	error:
1903	isl_mat_free(mat);
1904	return NULL;
1905	}
1906
1907	/* Return the gcd of the elements in row "row" of "mat" in *gcd.
1908	* Return isl_stat_ok on success and isl_stat_error on failure.
1909	*/
1910	isl_stat isl_mat_row_gcd(__isl_keep isl_mat *mat, int row, isl_int *gcd)
1911	{
1912	if (check_row(mat, row) < 0)
1913	return isl_stat_error;
1914
1915	isl_seq_gcd(p: mat->row[row], len: mat->n_col, gcd);
1916
1917	return isl_stat_ok;
1918	}
1919
1920	void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1921	{
1922	int i;
1923	isl_int g;
1924
1925	isl_int_set_si(*gcd, 0);
1926	if (!mat)
1927	return;
1928
1929	isl_int_init(g);
1930	for (i = 0; i < mat->n_row; ++i) {
1931	isl_seq_gcd(p: mat->row[i], len: mat->n_col, gcd: &g);
1932	isl_int_gcd(*gcd, *gcd, g);
1933	}
1934	isl_int_clear(g);
1935	}
1936
1937	/* Return the result of scaling "mat" by a factor of "m".
1938	*/
1939	__isl_give isl_mat *isl_mat_scale(__isl_take isl_mat *mat, isl_int m)
1940	{
1941	int i;
1942
1943	if (isl_int_is_one(m))
1944	return mat;
1945
1946	mat = isl_mat_cow(mat);
1947	if (!mat)
1948	return NULL;
1949
1950	for (i = 0; i < mat->n_row; ++i)
1951	isl_seq_scale(dst: mat->row[i], src: mat->row[i], f: m, len: mat->n_col);
1952
1953	return mat;
1954	}
1955
1956	__isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1957	{
1958	int i;
1959
1960	if (isl_int_is_one(m))
1961	return mat;
1962
1963	mat = isl_mat_cow(mat);
1964	if (!mat)
1965	return NULL;
1966
1967	for (i = 0; i < mat->n_row; ++i)
1968	isl_seq_scale_down(dst: mat->row[i], src: mat->row[i], f: m, len: mat->n_col);
1969
1970	return mat;
1971	}
1972
1973	__isl_give isl_mat *isl_mat_scale_down_row(__isl_take isl_mat *mat, int row,
1974	isl_int m)
1975	{
1976	if (isl_int_is_one(m))
1977	return mat;
1978
1979	mat = isl_mat_cow(mat);
1980	if (!mat)
1981	return NULL;
1982
1983	isl_seq_scale_down(dst: mat->row[row], src: mat->row[row], f: m, len: mat->n_col);
1984
1985	return mat;
1986	}
1987
1988	__isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1989	{
1990	isl_int gcd;
1991
1992	if (!mat)
1993	return NULL;
1994
1995	isl_int_init(gcd);
1996	isl_mat_gcd(mat, gcd: &gcd);
1997	mat = isl_mat_scale_down(mat, m: gcd);
1998	isl_int_clear(gcd);
1999
2000	return mat;
2001	}
2002
2003	__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
2004	{
2005	mat = isl_mat_cow(mat);
2006	if (!mat)
2007	return NULL;
2008
2009	isl_seq_normalize(ctx: mat->ctx, p: mat->row[row], len: mat->n_col);
2010
2011	return mat;
2012	}
2013
2014	/* Number of initial non-zero columns.
2015	*/
2016	int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
2017	{
2018	int i;
2019
2020	if (!mat)
2021	return -1;
2022
2023	for (i = 0; i < mat->n_col; ++i)
2024	if (row_first_non_zero(row: mat->row, n_row: mat->n_row, col: i) < 0)
2025	break;
2026
2027	return i;
2028	}
2029
2030	/* Return a basis for the space spanned by the rows of "mat".
2031	* Any basis will do, so simply perform Gaussian elimination and
2032	* remove the empty rows.
2033	*/
2034	__isl_give isl_mat *isl_mat_row_basis(__isl_take isl_mat *mat)
2035	{
2036	return isl_mat_reverse_gauss(mat);
2037	}
2038
2039	/* Return rows that extend a basis of "mat1" to one
2040	* that covers both "mat1" and "mat2".
2041	* The Hermite normal form of the concatenation of the two matrices is
2042	*
2043	* [ Q1 ]
2044	* [ M1 ] = [ H1 0 0 ] [ Q2 ]
2045	* [ M2 ] = [ H2 H3 0 ] [ Q3 ]
2046	*
2047	* The number of columns in H1 and H3 determine the number of rows
2048	* in Q1 and Q2. Q1 is a basis for M1, while Q2 extends this basis
2049	* to also cover M2.
2050	*/
2051	__isl_give isl_mat *isl_mat_row_basis_extension(
2052	__isl_take isl_mat *mat1, __isl_take isl_mat *mat2)
2053	{
2054	isl_size n_row;
2055	int r1, r;
2056	isl_size n1;
2057	isl_mat *H, *Q;
2058
2059	n1 = isl_mat_rows(mat: mat1);
2060	H = isl_mat_concat(top: mat1, bot: mat2);
2061	H = isl_mat_left_hermite(M: H, neg: 0, NULL, Q: &Q);
2062	if (n1 < 0 || !H || !Q)
2063	goto error;
2064
2065	r1 = hermite_first_zero_col(H, first_col: 0, n_row: n1);
2066	r = hermite_first_zero_col(H, first_col: r1, n_row: H->n_row);
2067	n_row = isl_mat_rows(mat: Q);
2068	if (n_row < 0)
2069	goto error;
2070	Q = isl_mat_drop_rows(mat: Q, row: r, n: n_row - r);
2071	Q = isl_mat_drop_rows(mat: Q, row: 0, n: r1);
2072
2073	isl_mat_free(mat: H);
2074	return Q;
2075	error:
2076	isl_mat_free(mat: H);
2077	isl_mat_free(mat: Q);
2078	return NULL;
2079	}
2080
2081	/* Are the rows of "mat1" linearly independent of those of "mat2"?
2082	* That is, is there no linear dependence among the combined rows
2083	* that is not already present in either "mat1" or "mat2"?
2084	* In other words, is the rank of "mat1" and "mat2" combined equal
2085	* to the sum of the ranks of "mat1" and "mat2"?
2086	*/
2087	isl_bool isl_mat_has_linearly_independent_rows(__isl_keep isl_mat *mat1,
2088	__isl_keep isl_mat *mat2)
2089	{
2090	isl_size r1, r2, r;
2091	isl_mat *mat;
2092
2093	r1 = isl_mat_rank(mat: mat1);
2094	if (r1 < 0)
2095	return isl_bool_error;
2096	if (r1 == 0)
2097	return isl_bool_true;
2098	r2 = isl_mat_rank(mat: mat2);
2099	if (r2 < 0)
2100	return isl_bool_error;
2101	if (r2 == 0)
2102	return isl_bool_true;
2103
2104	mat = isl_mat_concat(top: isl_mat_copy(mat: mat1), bot: isl_mat_copy(mat: mat2));
2105	r = isl_mat_rank(mat);
2106	isl_mat_free(mat);
2107	if (r < 0)
2108	return isl_bool_error;
2109	return isl_bool_ok(b: r == r1 + r2);
2110	}
2111