X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=isl_mat.c;h=e6638385d624481e844f9e80c6cb159aaf4881b9;hb=84037b07483c1970b915b7e829598a3cbd2c5d7f;hp=15c6b49af14db635572fadda39e614819620e7a7;hpb=177d99db8865eca40bf8baf52da173f287b8d3c5;p=platform%2Fupstream%2Fisl.git

diff --git a/isl_mat.c b/isl_mat.c
index 15c6b49..e663838 100644
--- a/isl_mat.c
+++ b/isl_mat.c
@@ -1,3 +1,4 @@
+#include "isl_dim.h"
 #include "isl_seq.h"
 #include "isl_mat.h"
 #include "isl_map_private.h"
@@ -12,10 +13,9 @@ struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
 	if (!mat)
 		return NULL;
 
-	mat->block.size = 0;
 	mat->row = NULL;
 	mat->block = isl_blk_alloc(ctx, n_row * n_col);
-	if (!mat->block.data)
+	if (isl_blk_is_error(mat->block))
 		goto error;
 	mat->row = isl_alloc_array(ctx, isl_int *, n_row);
 	if (!mat->row)
@@ -239,10 +239,10 @@ static void exchange(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
 		isl_mat_swap_rows(ctx, *Q, i, j);
 }
 
-static void subtract(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
+static void subtract(struct isl_mat *M, struct isl_mat **U,
 	struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
 {
-	int r, c;
+	int r;
 	for (r = row; r < M->n_row; ++r)
 		isl_int_submul(M->row[r][j], m, M->row[r][i]);
 	if (U) {
@@ -276,11 +276,13 @@ static void oppose(struct isl_ctx *ctx, struct isl_mat *M, struct isl_mat **U,
  *
  * with U and Q unimodular matrices and H a matrix in column echelon form
  * such that on each echelon row the entries in the non-echelon column
- * are non-negative and stricly smaller than the entries in the echelon
- * column.  If U or Q are NULL, then these matrices are not computed.
+ * are non-negative (if neg == 0) or non-positive (if neg == 1)
+ * and stricly smaller (in absolute value) than the entries in the echelon
+ * column.
+ * If U or Q are NULL, then these matrices are not computed.
  */
 struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
-	struct isl_mat *M, struct isl_mat **U, struct isl_mat **Q)
+	struct isl_mat *M, int neg, struct isl_mat **U, struct isl_mat **Q)
 {
 	isl_int c;
 	int row, col;
@@ -322,7 +324,7 @@ struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
 						       M->n_col-first)) != -1) {
 			first += off;
 			isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
-			subtract(ctx, M, U, Q, row, col, first, c);
+			subtract(M, U, Q, row, col, first, c);
 			if (!isl_int_is_zero(M->row[row][first]))
 				exchange(ctx, M, U, Q, row, first, col);
 			else
@@ -331,10 +333,13 @@ struct isl_mat *isl_mat_left_hermite(struct isl_ctx *ctx,
 		for (i = 0; i < col; ++i) {
 			if (isl_int_is_zero(M->row[row][i]))
 				continue;
-			isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
+			if (neg)
+				isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
+			else
+				isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
 			if (isl_int_is_zero(c))
 				continue;
-			subtract(ctx, M, U, Q, row, col, i, c);
+			subtract(M, U, Q, row, col, i, c);
 		}
 		++col;
 	}
@@ -517,6 +522,108 @@ error:
 	return NULL;
 }
 
+void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
+{
+	int i;
+
+	for (i = 0; i < mat->n_row; ++i)
+		isl_int_mul(mat->row[i][col], mat->row[i][col], m);
+}
+
+void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
+	isl_int m1, unsigned src1, isl_int m2, unsigned src2)
+{
+	int i;
+	isl_int tmp;
+
+	isl_int_init(tmp);
+	for (i = 0; i < mat->n_row; ++i) {
+		isl_int_mul(tmp, m1, mat->row[i][src1]);
+		isl_int_addmul(tmp, m2, mat->row[i][src2]);
+		isl_int_set(mat->row[i][dst], tmp);
+	}
+	isl_int_clear(tmp);
+}
+
+struct isl_mat *isl_mat_right_inverse(struct isl_ctx *ctx,
+	struct isl_mat *mat)
+{
+	struct isl_mat *inv;
+	int row;
+	isl_int a, b;
+
+	mat = isl_mat_cow(ctx, mat);
+	if (!mat)
+		return NULL;
+
+	inv = isl_mat_identity(ctx, mat->n_col);
+	inv = isl_mat_cow(ctx, inv);
+	if (!inv)
+		goto error;
+
+	isl_int_init(a);
+	isl_int_init(b);
+	for (row = 0; row < mat->n_row; ++row) {
+		int pivot, first, i, off;
+		pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
+		if (pivot < 0) {
+			isl_int_clear(a);
+			isl_int_clear(b);
+			goto error;
+		}
+		pivot += row;
+		if (pivot != row)
+			exchange(ctx, mat, &inv, NULL, row, pivot, row);
+		if (isl_int_is_neg(mat->row[row][row]))
+			oppose(ctx, mat, &inv, NULL, row, row);
+		first = row+1;
+		while ((off = isl_seq_first_non_zero(mat->row[row]+first,
+						    mat->n_col-first)) != -1) {
+			first += off;
+			isl_int_fdiv_q(a, mat->row[row][first],
+						    mat->row[row][row]);
+			subtract(mat, &inv, NULL, row, row, first, a);
+			if (!isl_int_is_zero(mat->row[row][first]))
+				exchange(ctx, mat, &inv, NULL, row, row, first);
+			else
+				++first;
+		}
+		for (i = 0; i < row; ++i) {
+			if (isl_int_is_zero(mat->row[row][i]))
+				continue;
+			isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
+			isl_int_divexact(b, mat->row[row][i], a);
+			isl_int_divexact(a, mat->row[row][row], a);
+			isl_int_neg(a, a);
+			isl_mat_col_combine(mat, i, a, i, b, row);
+			isl_mat_col_combine(inv, i, a, i, b, row);
+		}
+	}
+	isl_int_clear(b);
+
+	isl_int_set(a, mat->row[0][0]);
+	for (row = 1; row < mat->n_row; ++row)
+		isl_int_lcm(a, a, mat->row[row][row]);
+	if (isl_int_is_zero(a)){
+		isl_int_clear(a);
+		goto error;
+	}
+	for (row = 0; row < mat->n_row; ++row) {
+		isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
+		if (isl_int_is_one(mat->row[row][row]))
+			continue;
+		isl_mat_col_scale(inv, row, mat->row[row][row]);
+	}
+	isl_int_clear(a);
+
+	isl_mat_free(ctx, mat);
+
+	return inv;
+error:
+	isl_mat_free(ctx, mat);
+	return NULL;
+}
+
 struct isl_mat *isl_mat_swap_rows(struct isl_ctx *ctx,
 	struct isl_mat *mat, unsigned i, unsigned j)
 {
@@ -577,39 +684,84 @@ struct isl_basic_set *isl_basic_set_preimage(struct isl_ctx *ctx,
 	if (!bset || !mat)
 		goto error;
 
-	bset = isl_basic_set_cow(ctx, bset);
+	bset = isl_basic_set_cow(bset);
 	if (!bset)
 		goto error;
 
-	isl_assert(ctx, bset->nparam == 0, goto error);
+	isl_assert(ctx, bset->dim->nparam == 0, goto error);
 	isl_assert(ctx, bset->n_div == 0, goto error);
-	isl_assert(ctx, 1+bset->dim == mat->n_row, goto error);
+	isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
+
+	if (mat->n_col > mat->n_row)
+		bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0,
+						0, 0);
+	else if (mat->n_col < mat->n_row) {
+		bset->dim = isl_dim_cow(bset->dim);
+		if (!bset->dim)
+			goto error;
+		bset->dim->n_out -= mat->n_row - mat->n_col;
+		bset->extra += mat->n_row - mat->n_col;
+	}
 
 	t = isl_mat_sub_alloc(ctx, bset->eq, 0, bset->n_eq, 0, mat->n_row);
 	t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
 	if (!t)
 		goto error;
-	for (i = 0; i < bset->n_eq; ++i)
+	for (i = 0; i < bset->n_eq; ++i) {
 		isl_seq_swp_or_cpy(bset->eq[i], t->row[i], t->n_col);
+		isl_seq_clr(bset->eq[i]+t->n_col, bset->extra);
+	}
 	isl_mat_free(ctx, t);
 
 	t = isl_mat_sub_alloc(ctx, bset->ineq, 0, bset->n_ineq, 0, mat->n_row);
-	t = isl_mat_product(ctx, t, isl_mat_copy(ctx, mat));
+	t = isl_mat_product(ctx, t, mat);
 	if (!t)
-		goto error;
-	for (i = 0; i < bset->n_ineq; ++i)
+		goto error2;
+	for (i = 0; i < bset->n_ineq; ++i) {
 		isl_seq_swp_or_cpy(bset->ineq[i], t->row[i], t->n_col);
+		isl_seq_clr(bset->ineq[i]+t->n_col, bset->extra);
+	}
 	isl_mat_free(ctx, t);
 
-	bset->dim -= mat->n_row - mat->n_col;
-	bset = isl_basic_set_simplify(ctx, bset);
-	bset = isl_basic_set_finalize(ctx, bset);
+	bset = isl_basic_set_simplify(bset);
+	bset = isl_basic_set_finalize(bset);
 
-	isl_mat_free(ctx, mat);
 	return bset;
 error:
 	isl_mat_free(ctx, mat);
-	isl_basic_set_free(ctx, bset);
+error2:
+	isl_basic_set_free(bset);
+	return NULL;
+}
+
+struct isl_set *isl_set_preimage(struct isl_ctx *ctx,
+	struct isl_set *set, struct isl_mat *mat)
+{
+	int i;
+
+	set = isl_set_cow(set);
+	if (!set)
+		return NULL;
+
+	ctx = set->ctx;
+	for (i = 0; i < set->n; ++i) {
+		set->p[i] = isl_basic_set_preimage(ctx, set->p[i],
+						    isl_mat_copy(ctx, mat));
+		if (!set->p[i])
+			goto error;
+	}
+	if (mat->n_col != mat->n_row) {
+		set->dim = isl_dim_cow(set->dim);
+		if (!set->dim)
+			goto error;
+		set->dim->n_out += mat->n_col;
+		set->dim->n_out -= mat->n_row;
+	}
+	isl_mat_free(ctx, mat);
+	return set;
+error:
+	isl_set_free(set);
+	isl_mat_free(ctx, mat);
 	return NULL;
 }
 
@@ -618,6 +770,11 @@ void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
 {
 	int i, j;
 
+	if (!mat) {
+		fprintf(out, "null mat\n");
+		return;
+	}
+
 	if (mat->n_row == 0)
 		fprintf(out, "%*s[]\n", indent, "");
 
@@ -629,7 +786,7 @@ void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
 		for (j = 0; j < mat->n_col; ++j) {
 			if (j)
 			    fprintf(out, ",");
-			isl_int_print(out, mat->row[i][j]);
+			isl_int_print(out, mat->row[i][j], 0);
 		}
 		if (i == mat->n_row-1)
 			fprintf(out, "]]\n");
@@ -638,20 +795,20 @@ void isl_mat_dump(struct isl_ctx *ctx, struct isl_mat *mat,
 	}
 }
 
-struct isl_mat *isl_mat_drop_col(struct isl_ctx *ctx, struct isl_mat *mat,
-				unsigned col)
+struct isl_mat *isl_mat_drop_cols(struct isl_ctx *ctx, struct isl_mat *mat,
+				unsigned col, unsigned n)
 {
 	int r;
 
 	if (!mat)
 		return NULL;
 
-	if (col != mat->n_col-1) {
+	if (col != mat->n_col-n) {
 		for (r = 0; r < mat->n_row; ++r)
-			isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+1,
-					mat->n_col - col - 1);
+			isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
+					mat->n_col - col - n);
 	}
-	mat->n_col--;
+	mat->n_col -= n;
 	return mat;
 }