From 6eb17965a17c5055d7c79996069d3a54ace84062 Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Thu, 10 Aug 2006 10:53:07 +0200 Subject: [PATCH] util.c: compress_variables: extracted from lexmin.cc --- barvinok/util.h | 2 ++ lexmin.cc | 83 +++++++-------------------------------------------------- util.c | 78 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 89 insertions(+), 74 deletions(-) diff --git a/barvinok/util.h b/barvinok/util.h index 108bbbc..ee518dd 100644 --- a/barvinok/util.h +++ b/barvinok/util.h @@ -86,6 +86,8 @@ void Polyhedron_pprint(FILE *out, Polyhedron *P, int dim, int nparam, Polyhedron *Cone_over_Polyhedron(Polyhedron *P); +Matrix *compress_variables(Matrix *Equalities, unsigned nparam); + const char *barvinok_version(); #if defined(__cplusplus) diff --git a/lexmin.cc b/lexmin.cc index d6d1a9c..4a73fd8 100644 --- a/lexmin.cc +++ b/lexmin.cc @@ -1902,82 +1902,17 @@ static Matrix *compress_parameters(Polyhedron **P, unsigned nparam, unsigned Max static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays) { - unsigned dim = (*P)->Dimension - nparam; - Matrix *M, *H, *Q, *U, *C, *ratH, *invH, *Ul, *T1, *T2, *T; - Value mone; - int n; - - for (n = 0; n < (*P)->NbEq; ++n) - if (First_Non_Zero((*P)->Constraint[n]+1, dim) == -1) - break; - if (n == 0) - return NULL; - value_init(mone); - value_set_si(mone, -1); - M = Matrix_Alloc(n, dim); - C = Matrix_Alloc(n+1, nparam+1); - for (int i = 0; i < n; ++i) { - Vector_Copy((*P)->Constraint[i]+1, M->p[i], dim); - Vector_Scale((*P)->Constraint[i]+1+dim, C->p[i], mone, nparam+1); - } - value_set_si(C->p[n][nparam], 1); - left_hermite(M, &H, &Q, &U); - Matrix_Free(M); - Matrix_Free(Q); - value_clear(mone); + /* Matrix "view" of equalities */ + Matrix M; + M.NbRows = (*P)->NbEq; + M.NbColumns = (*P)->Dimension+2; + M.p_Init = (*P)->p_Init; + M.p = (*P)->Constraint; - /* we will need to treat scalings later */ - if (nparam > 0) - for (int i = 0; i < n; ++i) - assert(value_one_p(H->p[i][i])); - - ratH = Matrix_Alloc(n+1, n+1); - invH = Matrix_Alloc(n+1, n+1); - for (int i = 0; i < n; ++i) - Vector_Copy(H->p[i], ratH->p[i], n); - value_set_si(ratH->p[n][n], 1); - int ok = Matrix_Inverse(ratH, invH); - Matrix_Free(H); - Matrix_Free(ratH); - assert(ok); - T1 = Matrix_Alloc(n+1, nparam+1); - Matrix_Product(invH, C, T1); - if (nparam == 0 && value_notone_p(T1->p[n][nparam])) { - for (int i = 0; i < n; ++i) { - if (!mpz_divisible_p(T1->p[i][nparam], T1->p[n][nparam])) { - Matrix_Free(T1); - Matrix_Free(C); - Matrix_Free(invH); - Matrix_Free(U); - return NULL; - } - value_division(T1->p[i][nparam], T1->p[i][nparam], T1->p[n][nparam]); - } - value_set_si(T1->p[n][nparam], 1); - } - Matrix_Free(C); - Matrix_Free(invH); - Ul = Matrix_Alloc(dim+1, n+1); - for (int i = 0; i < dim; ++i) - Vector_Copy(U->p[i], Ul->p[i], n); - value_set_si(Ul->p[dim][n], 1); - T2 = Matrix_Alloc(dim+1, nparam+1); - Matrix_Product(Ul, T1, T2); - Matrix_Free(Ul); - Matrix_Free(T1); - - T = Matrix_Alloc(dim+nparam+1, (dim-n)+nparam+1); - for (int i = 0; i < dim; ++i) { - Vector_Copy(U->p[i]+n, T->p[i], dim-n); - Vector_Copy(T2->p[i], T->p[i]+dim-n, nparam+1); - } - for (int i = 0; i < nparam+1; ++i) - value_set_si(T->p[dim+i][(dim-n)+i], 1); - assert(value_one_p(T2->p[dim][nparam])); - Matrix_Free(U); - Matrix_Free(T2); + Matrix *T = compress_variables(&M, nparam); - *P = Polyhedron_Preimage(*P, T, MaxRays); + if (T) + *P = Polyhedron_Preimage(*P, T, MaxRays); return T; } diff --git a/util.c b/util.c index 302b4f9..c84e679 100644 --- a/util.c +++ b/util.c @@ -1400,6 +1400,84 @@ Polyhedron *Cone_over_Polyhedron(Polyhedron *P) return C; } +Matrix *compress_variables(Matrix *Equalities, unsigned nparam) +{ + unsigned dim = (Equalities->NbColumns-2) - nparam; + Matrix *M, *H, *Q, *U, *C, *ratH, *invH, *Ul, *T1, *T2, *T; + Value mone; + int n, i; + + for (n = 0; n < Equalities->NbRows; ++n) + if (First_Non_Zero(Equalities->p[n]+1, dim) == -1) + break; + if (n == 0) + return Identity(dim+nparam+1); + value_init(mone); + value_set_si(mone, -1); + M = Matrix_Alloc(n, dim); + C = Matrix_Alloc(n+1, nparam+1); + for (i = 0; i < n; ++i) { + Vector_Copy(Equalities->p[i]+1, M->p[i], dim); + Vector_Scale(Equalities->p[i]+1+dim, C->p[i], mone, nparam+1); + } + value_set_si(C->p[n][nparam], 1); + left_hermite(M, &H, &Q, &U); + Matrix_Free(M); + Matrix_Free(Q); + value_clear(mone); + + /* we will need to treat scalings later */ + if (nparam > 0) + for (i = 0; i < n; ++i) + assert(value_one_p(H->p[i][i])); + + ratH = Matrix_Alloc(n+1, n+1); + invH = Matrix_Alloc(n+1, n+1); + for (i = 0; i < n; ++i) + Vector_Copy(H->p[i], ratH->p[i], n); + value_set_si(ratH->p[n][n], 1); + int ok = Matrix_Inverse(ratH, invH); + Matrix_Free(H); + Matrix_Free(ratH); + assert(ok); + T1 = Matrix_Alloc(n+1, nparam+1); + Matrix_Product(invH, C, T1); + Matrix_Free(C); + Matrix_Free(invH); + if (nparam == 0 && value_notone_p(T1->p[n][nparam])) { + for (i = 0; i < n; ++i) { + if (!mpz_divisible_p(T1->p[i][nparam], T1->p[n][nparam])) { + Matrix_Free(T1); + Matrix_Free(U); + return NULL; + } + value_division(T1->p[i][nparam], T1->p[i][nparam], T1->p[n][nparam]); + } + value_set_si(T1->p[n][nparam], 1); + } + Ul = Matrix_Alloc(dim+1, n+1); + for (i = 0; i < dim; ++i) + Vector_Copy(U->p[i], Ul->p[i], n); + value_set_si(Ul->p[dim][n], 1); + T2 = Matrix_Alloc(dim+1, nparam+1); + Matrix_Product(Ul, T1, T2); + Matrix_Free(Ul); + Matrix_Free(T1); + + T = Matrix_Alloc(dim+nparam+1, (dim-n)+nparam+1); + for (i = 0; i < dim; ++i) { + Vector_Copy(U->p[i]+n, T->p[i], dim-n); + Vector_Copy(T2->p[i], T->p[i]+dim-n, nparam+1); + } + for (i = 0; i < nparam+1; ++i) + value_set_si(T->p[dim+i][(dim-n)+i], 1); + assert(value_one_p(T2->p[dim][nparam])); + Matrix_Free(U); + Matrix_Free(T2); + + return T; +} + const char *barvinok_version(void) { return -- 2.11.4.GIT