Home | History | Annotate | Line # | Download | only in dist
      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(mat->row[i], 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_row * n_col);
     65 	if (isl_blk_is_error(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, mat->block);
     87 	free(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(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(new_mat->row[i], mat->row[i], 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(mat->ctx, mat->block, n_row * mat->max_col);
    127 	if (isl_blk_is_error(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(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(mat->ctx, 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[i]+dst_col, src[i]+src_col, 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[i]+dst_col, src[i]+src_col, 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(mat->ctx, mat->n_row, mat->n_col);
    219 	if (!mat2)
    220 		return NULL;
    221 	for (i = 0; i < mat->n_row; ++i)
    222 		isl_seq_cpy(mat2->row[i], mat->row[i], 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(mat->ctx, mat->block);
    250 	isl_ctx_deref(mat->ctx);
    251 	free(mat->row);
    252 	free(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, 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->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_row);
    392 	if (!mat)
    393 		return NULL;
    394 	for (i = 0; i < n_row; ++i) {
    395 		isl_seq_clr(mat->row[i], i);
    396 		isl_int_set(mat->row[i][i], d);
    397 		isl_seq_clr(mat->row[i]+i+1, 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(mat->row[i], 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, 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(mat->row[i], 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(mat->row[i] + i + 1,
    443 					    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(mat->ctx, mat->n_row);
    462 	if (!prod)
    463 		goto error;
    464 
    465 	for (i = 0; i < prod->size; ++i)
    466 		isl_seq_inner_product(mat->row[i], vec->el, vec->size,
    467 					&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(vec->ctx, vec->size, 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(mat, vec_mat);
    491 	isl_vec_free(vec);
    492 	if (!vec_mat)
    493 		return NULL;
    494 	vec = isl_vec_alloc(vec_mat->ctx, 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(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(mat->ctx, 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(left->ctx, left->n_row, 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(sum->row[0]+1, 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(sum->row[i]+1, left->row[i]+1, left->row[0][0],
    568 				left->n_col-1);
    569 		isl_seq_scale(sum->row[i]+left->n_col,
    570 				right->row[i]+1, right->row[0][0],
    571 				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(left);
    577 	isl_mat_free(right);
    578 	return sum;
    579 error:
    580 	isl_mat_free(left);
    581 	isl_mat_free(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(*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((*Q)->row[col], (*Q)->row[col], (*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(M->ctx, M->n_col);
    655 		if (!*U)
    656 			goto error;
    657 	}
    658 	if (Q) {
    659 		*Q = isl_mat_identity(M->ctx, 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(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(M->row[row]+col, M->n_col-col);
    676 		if (first == -1)
    677 			continue;
    678 		first += col;
    679 		if (first != col)
    680 			exchange(M, U, Q, row, first, 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(M->row[row]+first,
    685 						       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, col, first, c);
    689 			if (!isl_int_is_zero(M->row[row][first]))
    690 				exchange(M, U, Q, row, first, 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, col, i, c);
    704 		}
    705 		++col;
    706 	}
    707 	isl_int_clear(c);
    708 
    709 	return M;
    710 error:
    711 	if (Q) {
    712 		isl_mat_free(*Q);
    713 		*Q = NULL;
    714 	}
    715 	if (U) {
    716 		isl_mat_free(*U);
    717 		*U = NULL;
    718 	}
    719 	isl_mat_free(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(mat->row[k], mat->row[row], col, nc, NULL);
    747 		isl_seq_normalize(ctx, mat->row[k], 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, k, 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, last);
    790 		if (!mat)
    791 			return NULL;
    792 	}
    793 	mat = isl_mat_drop_rows(mat, 0, 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(mat->row[i], 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, 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(isl_mat_copy(mat), 0, NULL, NULL);
    860 	if (!H)
    861 		return isl_size_error;
    862 
    863 	rank = hermite_first_zero_col(H, 0, H->n_row);
    864 	isl_mat_free(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(mat, 0, &U, NULL);
    876 	if (!mat || !U)
    877 		goto error;
    878 
    879 	rank = hermite_first_zero_col(mat, 0, mat->n_row);
    880 	K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank);
    881 	if (!K)
    882 		goto error;
    883 	isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
    884 	isl_mat_free(mat);
    885 	isl_mat_free(U);
    886 	return K;
    887 error:
    888 	isl_mat_free(mat);
    889 	isl_mat_free(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(mat->ctx, 1+mat->n_row, 1+mat->n_col);
    901 	if (!mat2)
    902 		goto error;
    903 	isl_int_set_si(mat2->row[0][0], 1);
    904 	isl_seq_clr(mat2->row[0]+1, 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(mat2->row[1+i]+1, mat->row[i], 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(mat1->ctx, mat1->n_row + mat2->n_row,
    931 				       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(mat->row[i], mat1->row[i], mat1->n_col);
    936 		isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
    937 	}
    938 	for (i = 0; i < mat2->n_row; ++i) {
    939 		isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
    940 		isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
    941 						    mat2->row[i], mat2->n_col);
    942 	}
    943 	isl_mat_free(mat1);
    944 	isl_mat_free(mat2);
    945 	return mat;
    946 error:
    947 	isl_mat_free(mat1);
    948 	isl_mat_free(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(*left, i, j);
    980 	*right = isl_mat_swap_rows(*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(left->row[row]+row, left->row[row]+row, left->n_col-row);
    991 	isl_seq_neg(right->row[row], right->row[row], 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(left->row[i]+row,
    999 			left->ctx->one, left->row[i]+row,
   1000 			m, left->row[row]+row,
   1001 			left->n_col-row);
   1002 	isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
   1003 			m, right->row[row], 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(left);
   1022 		return right;
   1023 	}
   1024 
   1025 	left = isl_mat_cow(left);
   1026 	right = isl_mat_cow(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(left->row+row, left->n_row-row, 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, &right, pivot, 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(left->row+first,
   1048 					left->n_row-first, 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, first, a);
   1053 			if (!isl_int_is_zero(left->row[first][row])) {
   1054 				if (inv_exchange(&left, &right, row, 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(left->row[i] + i,
   1068 					a, left->row[i] + i,
   1069 					b, left->row[row] + i,
   1070 					left->n_col - i);
   1071 			isl_seq_combine(right->row[i], a, right->row[i],
   1072 					b, right->row[row], 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(right->row[row], right->row[row],
   1089 				left->row[row][row], right->n_col);
   1090 	}
   1091 	isl_int_clear(a);
   1092 
   1093 	isl_mat_free(left);
   1094 	return right;
   1095 error:
   1096 	isl_mat_free(left);
   1097 	isl_mat_free(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(mat->ctx, mat->n_col);
   1135 	inv = isl_mat_cow(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(mat->row[row]+row, 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(mat, &inv, NULL, row, pivot, row);
   1152 		if (isl_int_is_neg(mat->row[row][row]))
   1153 			oppose(mat, &inv, NULL, row, row);
   1154 		first = row+1;
   1155 		while ((off = isl_seq_first_non_zero(mat->row[row]+first,
   1156 						    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(mat, &inv, NULL, row, row, first, a);
   1161 			if (!isl_int_is_zero(mat->row[row][first]))
   1162 				exchange(mat, &inv, NULL, row, row, 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, i, a, i, b, row);
   1174 			isl_mat_col_combine(inv, i, a, i, b, 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(inv, row, 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(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(mat->ctx, mat->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, i, 1) < 0 ||
   1240 	    check_col_range(mat, j, 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, i, 1) < 0 ||
   1257 	    check_row_range(mat, j, 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(left->ctx, left->n_row, 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(prod->row[i], prod->n_col);
   1286 		isl_mat_free(left);
   1287 		isl_mat_free(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(left);
   1303 	isl_mat_free(right);
   1304 	return prod;
   1305 error:
   1306 	isl_mat_free(left);
   1307 	isl_mat_free(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(mat->ctx, q, 0, n, has_div, mat->n_row);
   1337 	t = isl_mat_product(t, mat);
   1338 	if (!t)
   1339 		return -1;
   1340 	for (i = 0; i < n; ++i) {
   1341 		isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
   1342 		isl_seq_cpy(q[i] + has_div + t->n_col,
   1343 			    q[i] + has_div + t->n_col + e, n_div);
   1344 		isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
   1345 	}
   1346 	isl_mat_free(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, isl_dim_set,
   1378 						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(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, bset->eq, bset->n_eq, bset->n_div, 0,
   1389 			isl_mat_copy(mat)) < 0)
   1390 		goto error;
   1391 
   1392 	if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
   1393 			isl_mat_copy(mat)) < 0)
   1394 		goto error;
   1395 
   1396 	if (preimage(ctx, bset->div, bset->n_div, bset->n_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(set->p[i],
   1427 						    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(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, 0, n_row, first_col, mat->n_row);
   1462 	t = isl_mat_product(t, mat);
   1463 	if (!t)
   1464 		return isl_stat_error;
   1465 	for (i = 0; i < n_row; ++i)
   1466 		isl_seq_swp_or_cpy(row[i] + first_col, t->row[i], t->n_col);
   1467 	isl_mat_free(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(out, "%*snull mat\n", indent, "");
   1477 		return;
   1478 	}
   1479 
   1480 	if (mat->n_row == 0)
   1481 		fprintf(out, "%*s[]\n", indent, "");
   1482 
   1483 	for (i = 0; i < mat->n_row; ++i) {
   1484 		if (!i)
   1485 			fprintf(out, "%*s[[", indent, "");
   1486 		else
   1487 			fprintf(out, "%*s[", indent+1, "");
   1488 		for (j = 0; j < mat->n_col; ++j) {
   1489 			if (j)
   1490 			    fprintf(out, ",");
   1491 			isl_int_print(out, mat->row[i][j], 0);
   1492 		}
   1493 		if (i == mat->n_row-1)
   1494 			fprintf(out, "]]\n");
   1495 		else
   1496 			fprintf(out, "]\n");
   1497 	}
   1498 }
   1499 
   1500 void isl_mat_dump(__isl_keep isl_mat *mat)
   1501 {
   1502 	isl_mat_print_internal(mat, stderr, 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, 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(mat->row[r]+col, mat->row[r]+col+n,
   1520 					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, 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, col, 0) < 0)
   1548 		return isl_mat_free(mat);
   1549 	if (n == 0)
   1550 		return mat;
   1551 
   1552 	ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
   1553 	if (!ext)
   1554 		goto error;
   1555 
   1556 	isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
   1557 	isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
   1558 				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, first, n);
   1575 	if (!mat)
   1576 		return NULL;
   1577 
   1578 	for (i = 0; i < mat->n_row; ++i)
   1579 		isl_seq_clr(mat->row[i] + first, 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, 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, row, 0) < 0)
   1598 		return isl_mat_free(mat);
   1599 	if (n == 0)
   1600 		return mat;
   1601 
   1602 	ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
   1603 	if (!ext)
   1604 		goto error;
   1605 
   1606 	isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
   1607 	isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
   1608 				mat->n_row - row, 0, 0, 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, 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(mat->row[row + i], 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, 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, dst_col) < 0 || check_col(mat, 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(mat->row[row], mat->n_col) == -1)
   1729 		return mat;
   1730 	mat = isl_mat_cow(mat);
   1731 	if (!mat)
   1732 		return NULL;
   1733 	isl_seq_neg(mat->row[row], mat->row[row], 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(isl_mat_copy(M), 0, NULL, &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(M->row[r], Q->row[r], M->n_col);
   1755 	isl_mat_free(H);
   1756 	isl_mat_free(Q);
   1757 	return M;
   1758 error:
   1759 	isl_mat_free(H);
   1760 	isl_mat_free(Q);
   1761 	isl_mat_free(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(top);
   1776 		return bot;
   1777 	}
   1778 	if (bot->n_row == 0) {
   1779 		isl_mat_free(bot);
   1780 		return top;
   1781 	}
   1782 
   1783 	mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
   1784 	if (!mat)
   1785 		goto error;
   1786 	isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
   1787 			 0, 0, mat->n_col);
   1788 	isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
   1789 			 0, 0, mat->n_col);
   1790 	isl_mat_free(top);
   1791 	isl_mat_free(bot);
   1792 	return mat;
   1793 error:
   1794 	isl_mat_free(top);
   1795 	isl_mat_free(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(mat1->row[i], mat2->row[i], 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(vec->ctx, 1, vec->size);
   1826 	if (!mat)
   1827 		goto error;
   1828 
   1829 	isl_seq_cpy(mat->row[0], vec->el, 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(isl_mat_get_ctx(mat), mat->n_col);
   1851 	if (!v)
   1852 		return NULL;
   1853 	isl_seq_cpy(v->el, mat->row[row], 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, isl_mat_from_row_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(mat->ctx, mat->n_row, mat->n_col);
   1875 	if (!res)
   1876 		goto error;
   1877 
   1878 	if (dst_col < src_col) {
   1879 		isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
   1880 				 0, 0, dst_col);
   1881 		isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
   1882 				 dst_col, src_col, n);
   1883 		isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
   1884 				 dst_col + n, dst_col, src_col - dst_col);
   1885 		isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
   1886 				 src_col + n, src_col + n,
   1887 				 res->n_col - src_col - n);
   1888 	} else {
   1889 		isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
   1890 				 0, 0, src_col);
   1891 		isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
   1892 				 src_col, src_col + n, dst_col - src_col);
   1893 		isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
   1894 				 dst_col, src_col, n);
   1895 		isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
   1896 				 dst_col + n, dst_col + n,
   1897 				 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(mat->row[row], 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(mat->row[i], mat->n_col, &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(mat->row[i], mat->row[i], m, 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(mat->row[i], mat->row[i], m, 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(mat->row[row], mat->row[row], m, 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);
   1997 	mat = isl_mat_scale_down(mat, 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(mat->ctx, mat->row[row], 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(mat->row, mat->n_row, 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(mat1);
   2060 	H = isl_mat_concat(mat1, mat2);
   2061 	H = isl_mat_left_hermite(H, 0, NULL, &Q);
   2062 	if (n1 < 0 || !H || !Q)
   2063 		goto error;
   2064 
   2065 	r1 = hermite_first_zero_col(H, 0, n1);
   2066 	r = hermite_first_zero_col(H, r1, H->n_row);
   2067 	n_row = isl_mat_rows(Q);
   2068 	if (n_row < 0)
   2069 		goto error;
   2070 	Q = isl_mat_drop_rows(Q, r, n_row - r);
   2071 	Q = isl_mat_drop_rows(Q, 0, r1);
   2072 
   2073 	isl_mat_free(H);
   2074 	return Q;
   2075 error:
   2076 	isl_mat_free(H);
   2077 	isl_mat_free(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(mat1);
   2094 	if (r1 < 0)
   2095 		return isl_bool_error;
   2096 	if (r1 == 0)
   2097 		return isl_bool_true;
   2098 	r2 = isl_mat_rank(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(isl_mat_copy(mat1), isl_mat_copy(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(r == r1 + r2);
   2110 }
   2111