From b32aee8ebb8d0ad4204dc975d27d2b474cf98dcc Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Tue, 30 Oct 2007 21:45:35 +0100 Subject: [PATCH] Use zsolve to compute Hilbert basis of a cone --- Makefile.am | 9 +++- cone_hilbert_basis.c | 24 +++++++++ hilbert.c | 146 +++++++++++++++++++++++++++++++++++++++++++++++++++ hilbert.h | 11 ++++ testlib.cc | 18 +++++++ 5 files changed, 206 insertions(+), 2 deletions(-) create mode 100644 cone_hilbert_basis.c create mode 100644 hilbert.c create mode 100644 hilbert.h diff --git a/Makefile.am b/Makefile.am index 1a96520..854c2b8 100644 --- a/Makefile.am +++ b/Makefile.am @@ -18,6 +18,8 @@ polylib/libpolylibgmp.la: FORCE $(MAKE) $(AM_MAKEFLAGS) -C polylib libpolylibgmp.la piplib/libpiplibMP.la: FORCE $(MAKE) $(AM_MAKEFLAGS) -C piplib libpiplibMP.la +zsolve/libzsolve.la: FORCE + $(MAKE) $(AM_MAKEFLAGS) -C zsolve libzsolve.la AM_CPPFLAGS = -I$(srcdir)/bernstein/include -I$(srcdir)/lib \ @POLYLIB_CPPFLAGS@ @PIPLIB_CPPFLAGS@ @GINACLIB_CPPFLAGS@ @@ -30,7 +32,8 @@ noinst_PROGRAMS = test testlib randomtest \ remove_redundant_equalities \ barvinok_union polytope_volume test_approx \ barvinok_summate verify_lexsmaller piptest \ - polyhedron_sample 4coins lexmin @bv_barvinok_maximize@ + polyhedron_sample 4coins lexmin @bv_barvinok_maximize@ \ + cone_hilbert_basis EXTRA_PROGRAMS = barvinok_maximize pkginclude_HEADERS = \ barvinok/NTL_QQ.h \ @@ -91,6 +94,8 @@ libbarvinok_la_SOURCES = \ euler.h \ genfun_constructor.cc \ genfun_constructor.h \ + hilbert.c \ + hilbert.h \ lattice_point.cc \ lattice_point.h \ options.c \ @@ -148,7 +153,7 @@ if BUNDLED_PIPLIB endif libbarvinok_la_LIBADD = @LTLIBOBJS@ $(POLYLIB_LA) @POLYLIB_LIBS@ \ $(PIPLIB_LA) @PIPLIB_LIBS@ $(BERNSTEIN_LA) \ - lib/libgnu.la + zsolve/libzsolve.la lib/libgnu.la libbarvinok_la_LDFLAGS = @BV_LDFLAGS@ -version-info @versioninfo@ LDADD = libbarvinok.la @POLYLIB_LIBS@ @PIPLIB_LIBS@ diff --git a/cone_hilbert_basis.c b/cone_hilbert_basis.c new file mode 100644 index 0000000..71a4a64 --- /dev/null +++ b/cone_hilbert_basis.c @@ -0,0 +1,24 @@ +#include +#include +#include "hilbert.h" + +int main(int argc, char **argv) +{ + Matrix *M; + Polyhedron *C; + struct barvinok_options *options = barvinok_options_new_with_defaults(); + + M = Matrix_Read(); + C = Constraints2Polyhedron(M, options->MaxRays); + Matrix_Free(M); + + M = Cone_Hilbert_Basis(C, options->MaxRays); + + Polyhedron_Free(C); + + Matrix_Print(stdout, P_VALUE_FMT, M); + Matrix_Free(M); + + barvinok_options_free(options); + return 0; +} diff --git a/hilbert.c b/hilbert.c new file mode 100644 index 0000000..a25eb81 --- /dev/null +++ b/hilbert.c @@ -0,0 +1,146 @@ +#define Vector ZSolveVector +#define Matrix ZSolveMatrix +#include "zsolve/libzsolve.h" +#undef Vector +#undef Matrix +#include "hilbert.h" + +static ZSolveMatrix Matrix2zsolve(Matrix *M) +{ + int i, j; + ZSolveMatrix zmatrix; + + zmatrix = createMatrix(M->NbColumns-2, M->NbRows); + for (i = 0; i < M->NbRows; ++i) + for (j = 0; j < M->NbColumns-2; ++j) { + assert(mpz_cmp_si(M->p[i][1+j], -MAXINT) > 0); + assert(mpz_cmp_si(M->p[i][1+j], MAXINT) < 0); + zmatrix->Data[i*zmatrix->Width+j] = mpz_get_si(M->p[i][1+j]); + } + + return zmatrix; +} + +static Matrix *VectorArray2Matrix(VectorArray array, unsigned cols) +{ + int i, j; + Matrix *M = Matrix_Alloc(array->Size, cols+1); + + for (i = 0; i < array->Size; ++i) { + for (j = 0; j < cols; ++j) + value_set_si(M->p[i][j], array->Data[i][j]); + value_set_si(M->p[i][cols], 1); + } + return M; +} + +static void Polyhedron_Remove_Positivity_Constraint(Polyhedron *P) +{ + int i; + + for (i = 0; i < P->NbConstraints; ++i) { + if (First_Non_Zero(P->Constraint[i]+1, P->Dimension) != -1) + continue; + if (i < P->NbConstraints-1) + Vector_Exchange(P->Constraint[i], + P->Constraint[P->NbConstraints-1], + P->Dimension+2); + P->NbConstraints--; + --i; + } +} + +static Matrix *Polyhedron2standard_form(Polyhedron *P, Matrix **T) +{ + int i; + unsigned dim = P->Dimension; + Matrix *M2 = Matrix_Alloc(P->NbConstraints+1, dim+1); + Matrix *H, *Q; + + Polyhedron_Remove_Positivity_Constraint(P); + for (i = 0; i < P->NbConstraints; ++i) { + assert(value_zero_p(P->Constraint[i][1+dim])); + Vector_Copy(P->Constraint[i]+1, M2->p[i], dim); + } + value_set_si(M2->p[P->NbConstraints][dim], 1); + neg_left_hermite(M2, &H, &Q, T); + Matrix_Free(Q); + Matrix_Free(M2); + + M2 = Matrix_Alloc(P->NbConstraints, 2+dim+(P->NbConstraints-P->NbEq)); + + for (i = 0; i < P->NbEq; ++i) + Vector_Copy(H->p[i], M2->p[i]+1, dim); + for (i = P->NbEq; i < P->NbConstraints; ++i) { + Vector_Copy(H->p[i], M2->p[i]+1, dim); + value_set_si(M2->p[i][1+dim+i-P->NbEq], -1); + } + Matrix_Free(H); + return M2; +} + +/* Assumes C is a linear cone (i.e. with apex zero). + * All equalities are removed first to speed up the computation + * in zsolve. + */ +Matrix *Cone_Hilbert_Basis(Polyhedron *C, unsigned MaxRays) +{ + unsigned dim; + int i; + Matrix *M2, *M3, *T; + Matrix *CV = NULL; + LinearSystem initialsystem; + ZSolveMatrix matrix; + ZSolveVector rhs; + ZSolveContext ctx; + + remove_all_equalities(&C, NULL, NULL, &CV, 0, MaxRays); + dim = C->Dimension; + + for (i = 0; i < C->NbConstraints; ++i) + assert(value_zero_p(C->Constraint[i][1+dim]) || + First_Non_Zero(C->Constraint[i]+1, dim) == -1); + + M2 = Polyhedron2standard_form(C, &T); + matrix = Matrix2zsolve(M2); + Matrix_Free(M2); + + rhs = createVector(matrix->Height); + for (i = 0; i < matrix->Height; i++) + rhs[i] = 0; + + initialsystem = createLinearSystem(); + setLinearSystemMatrix(initialsystem, matrix); + deleteMatrix(matrix); + + setLinearSystemRHS(initialsystem, rhs); + deleteVector(rhs); + + setLinearSystemLimit(initialsystem, -1, 0, MAXINT, 0); + setLinearSystemEquationType(initialsystem, -1, EQUATION_EQUAL, 0); + + ctx = createZSolveContextFromSystem(initialsystem, NULL, 0, 0, NULL, NULL); + zsolveSystem(ctx, 0); + + M2 = VectorArray2Matrix(ctx->Homs, C->Dimension); + deleteZSolveContext(ctx, 1); + Matrix_Transposition(T); + M3 = Matrix_Alloc(M2->NbRows, M2->NbColumns); + Matrix_Product(M2, T, M3); + Matrix_Free(M2); + Matrix_Free(T); + + if (CV) { + Matrix *T, *M; + T = Transpose(CV); + M = Matrix_Alloc(M3->NbRows, T->NbColumns); + Matrix_Product(M3, T, M); + Matrix_Free(M3); + Matrix_Free(CV); + Matrix_Free(T); + Polyhedron_Free(C); + M3 = M; + } + + return M3; +} diff --git a/hilbert.h b/hilbert.h new file mode 100644 index 0000000..aff1571 --- /dev/null +++ b/hilbert.h @@ -0,0 +1,11 @@ +#include + +#if defined(__cplusplus) +extern "C" { +#endif + +Matrix *Cone_Hilbert_Basis(Polyhedron *C, unsigned MaxRays); + +#if defined(__cplusplus) +} +#endif diff --git a/testlib.cc b/testlib.cc index bf3067a..4b19838 100644 --- a/testlib.cc +++ b/testlib.cc @@ -11,6 +11,7 @@ #include "lattice_point.h" #include "counter.h" #include "bernoulli.h" +#include "hilbert.h" #include "matrix_read.h" using std::cout; @@ -469,6 +470,22 @@ int test_bernoulli_sum(struct barvinok_options *options) evalue_free(sum2); } +int test_hilbert(struct barvinok_options *options) +{ + Matrix *M = matrix_read_from_str( + "2 4\n" + " 1 4 -3 0 \n" + " 1 3 2 0 \n"); + Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays); + Matrix_Free(M); + + M = Cone_Hilbert_Basis(P, options->MaxRays); + Polyhedron_Free(P); + assert(M->NbRows = 5); + assert(M->NbColumns = 3); + Matrix_Free(M); +} + int main(int argc, char **argv) { struct barvinok_options *options = barvinok_options_new_with_defaults(); @@ -482,5 +499,6 @@ int main(int argc, char **argv) test_todd(options); test_bernoulli(options); test_bernoulli_sum(options); + test_hilbert(options); barvinok_options_free(options); } -- 2.11.4.GIT