From 9ba53ea693fc8b68672354d7447e5358b262a5ec Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Fri, 29 Jun 2007 16:11:58 +0200 Subject: [PATCH] Polyhedron_Sample: use lp solver for finding range of a variable --- Makefile.am | 18 +++++----- polysign.c | 14 ++++++++ polysign.h | 18 ++++++++++ polysign_cdd.cc | 2 ++ polysign_cdd_template.cc | 70 ++++++++++++++++++++++++++++++++++----- polysign_cddf.cc | 5 +++ polysign_polylib.c | 28 ++++++++++++++++ sample.c | 86 +++++++++++++++++++++++++++--------------------- 8 files changed, 186 insertions(+), 55 deletions(-) diff --git a/Makefile.am b/Makefile.am index ed4c19b..5d61c86 100644 --- a/Makefile.am +++ b/Makefile.am @@ -46,6 +46,9 @@ endif if HAVE_GLPK BR_GLPK = basis_reduction_glpk.c endif +if HAVE_CDDLIB + POLYSIGN = polysign_cddf.cc polysign_cdd.cc +endif libbarvinok_la_SOURCES = \ barvinok/evalue.h \ barvinok/genfun.h \ @@ -76,6 +79,10 @@ libbarvinok_la_SOURCES = \ options.c \ param_util.c \ param_util.h \ + $(POLYSIGN) \ + polysign.c \ + polysign_polylib.c \ + polysign.h \ reduce_domain.c \ reduce_domain.h \ reducer.cc \ @@ -99,6 +106,9 @@ EXTRA_libbarvinok_la_SOURCES = \ basis_reduction_glpk.c \ bernstein.cc \ piputil.h \ + polysign_cddf.cc \ + polysign_cdd.cc \ + polysign.c \ initcdd.cc \ initcdd.h if HAVE_GINAC @@ -169,9 +179,6 @@ barvinok_summate_SOURCES = \ verify.c EXTRA_barvinok_summate_SOURCES = fdstream.cc fdstream.h 4coins_SOURCES = 4coins.cc -if HAVE_CDDLIB -POLYSIGN = polysign_cddf.cc polysign_cdd.cc -endif lexmin_SOURCES = \ lexmin.h \ lexmin.cc \ @@ -181,13 +188,8 @@ lexmin_SOURCES = \ edomain.h \ evalue_util.cc \ evalue_util.h \ - $(POLYSIGN) \ - polysign.c \ - polysign_polylib.c \ - polysign.h \ verify.h \ verify.c -EXTRA_lexmin_SOURCES = polysign_cddf.cc polysign_cdd.cc polysign.c test_approx_SOURCES = \ test_approx.c \ verify.h \ diff --git a/polysign.c b/polysign.c index ed2cc7f..4a5e85f 100644 --- a/polysign.c +++ b/polysign.c @@ -29,3 +29,17 @@ enum order_sign polyhedron_affine_sign(Polyhedron *D, Matrix *T, else assert(0); } + +enum lp_result polyhedron_range(Polyhedron *D, Value *obj, Value denom, + Value *min, Value *max, + struct barvinok_options *options) +{ + if (options->lp_solver == BV_LP_POLYLIB) + return PL_polyhedron_range(D, obj, denom, min, max, options); + else if (options->lp_solver == BV_LP_CDD) + return cdd_polyhedron_range(D, obj, denom, min, max, options); + else if (options->lp_solver == BV_LP_CDDF) + return cddf_polyhedron_range(D, obj, denom, min, max, options); + else + assert(0); +} diff --git a/polysign.h b/polysign.h index 7af60e4..37bd430 100644 --- a/polysign.h +++ b/polysign.h @@ -19,6 +19,24 @@ enum order_sign cdd_polyhedron_affine_sign(Polyhedron *D, Matrix *T, enum order_sign cddf_polyhedron_affine_sign(Polyhedron *D, Matrix *T, struct barvinok_options *options); +enum lp_result { lp_ok = 0, lp_unbounded, lp_empty }; + +/* Computes min and max of obj/denom over D, with min rounded up + * and max rounded down to an integer. obj has D->Dimension+1 + * elements, with the constant term in position D->Dimension. + */ +enum lp_result polyhedron_range(Polyhedron *D, Value *obj, Value denom, + Value *min, Value *max, + struct barvinok_options *options); +enum lp_result PL_polyhedron_range(Polyhedron *D, Value *obj, Value denom, + Value *min, Value *max, + struct barvinok_options *options); +enum lp_result cdd_polyhedron_range(Polyhedron *D, Value *obj, Value denom, + Value *min, Value *max, + struct barvinok_options *options); +enum lp_result cddf_polyhedron_range(Polyhedron *D, Value *obj, Value denom, + Value *min, Value *max, + struct barvinok_options *options); #if defined(__cplusplus) } #endif diff --git a/polysign_cdd.cc b/polysign_cdd.cc index 6abeb69..dc529fd 100644 --- a/polysign_cdd.cc +++ b/polysign_cdd.cc @@ -38,5 +38,7 @@ #define DD_rat_sign(sign,obj,val) sign = dd_sgn(val) #define DD_set_z(a,b) mpq_set_z(a,b) +#define DD_floor(a,b) mpz_fdiv_q(a,mpq_numref(b),mpq_denref(b)) +#define DD_ceil(a,b) mpz_cdiv_q(a,mpq_numref(b),mpq_denref(b)) #include "polysign_cdd_template.cc" diff --git a/polysign_cdd_template.cc b/polysign_cdd_template.cc index 7fc8beb..92b3dcc 100644 --- a/polysign_cdd_template.cc +++ b/polysign_cdd_template.cc @@ -8,12 +8,10 @@ #define EMPTY_DOMAIN -2 -static int polyhedron_affine_minmax(DD_LPObjectiveType obj, Polyhedron *P, - Matrix *T, bool rational) +static DD_LPType *solve_lp(DD_LPObjectiveType obj, Polyhedron *P, + Value *f) { - DD_LPType *lp; - assert(P->Dimension == T->NbColumns-1); - assert(T->NbRows == 2); + DD_LPType *lp; DD_rowrange irev = P->NbConstraints; DD_rowrange rows = irev + P->NbEq + 1; DD_colrange cols = 1 + P->Dimension; @@ -37,13 +35,50 @@ static int polyhedron_affine_minmax(DD_LPObjectiveType obj, Polyhedron *P, } /* objective function */ for (DD_colrange k = 0; k < P->Dimension; ++k) - DD_set_z(lp->A[rows-1][1+k], T->p[0][k]); - DD_set_z(lp->A[rows-1][0], T->p[0][P->Dimension]); + DD_set_z(lp->A[rows-1][1+k], f[k]); + DD_set_z(lp->A[rows-1][0], f[P->Dimension]); DD_ErrorType err = DD_NoError; DD_LPSolve(lp, DD_DualSimplex, &err); assert(err == DD_NoError); + return lp; +} + +static lp_result polyhedron_affine_minmax(DD_LPObjectiveType obj, Polyhedron *P, + Value *f, Value denom, Value *opt) +{ + lp_result res = lp_ok; + DD_LPType *lp = solve_lp(obj, P, f); + assert(value_one_p(denom)); + + switch(lp->LPS) { + case DD_Optimal: + if (obj == DD_LPmin) + DD_ceil(*opt, lp->optvalue); + else + DD_floor(*opt, lp->optvalue); + break; + case DD_DualInconsistent: + res = lp_unbounded; + break; + case DD_Inconsistent: + res = lp_empty; + break; + default: + assert(0); + } + DD_FreeLPData(lp); + return res; +} + +static int polyhedron_affine_minmax_sign(DD_LPObjectiveType obj, Polyhedron *P, + Matrix *T, bool rational) +{ + assert(P->Dimension == T->NbColumns-1); + assert(T->NbRows == 2); + DD_LPType *lp = solve_lp(obj, P, T->p[0]); + int sign; if (lp->LPS == DD_Optimal) { if (rational) @@ -75,12 +110,12 @@ enum order_sign cdd_polyhedron_affine_sign(Polyhedron *D, Matrix *T, INIT_CDD; bool rational = !POL_ISSET(options->MaxRays, POL_INTEGER); - int min = polyhedron_affine_minmax(DD_LPmin, D, T, rational); + int min = polyhedron_affine_minmax_sign(DD_LPmin, D, T, rational); if (min == EMPTY_DOMAIN) return order_undefined; if (min > 0) return order_gt; - int max = polyhedron_affine_minmax(DD_LPmax, D, T, rational); + int max = polyhedron_affine_minmax_sign(DD_LPmax, D, T, rational); assert(max != EMPTY_DOMAIN); if (max < 0) return order_lt; @@ -92,3 +127,20 @@ enum order_sign cdd_polyhedron_affine_sign(Polyhedron *D, Matrix *T, return order_ge; return order_unknown; } + +enum lp_result cdd_polyhedron_range(Polyhedron *D, Value *obj, Value denom, + Value *min, Value *max, + struct barvinok_options *options) +{ + lp_result res; + + if (emptyQ2(D)) + return lp_empty; + + INIT_CDD; + res = polyhedron_affine_minmax(DD_LPmin, D, obj, denom, min); + if (res != lp_ok) + return res; + res = polyhedron_affine_minmax(DD_LPmax, D, obj, denom, max); + return res; +} diff --git a/polysign_cddf.cc b/polysign_cddf.cc index 7304c3e..612623d 100644 --- a/polysign_cddf.cc +++ b/polysign_cddf.cc @@ -1,3 +1,5 @@ +#include + #define DD_LPType ddf_LPType #define DD_CreateLPData ddf_CreateLPData #define DD_LPObjectiveType ddf_LPObjectiveType @@ -27,6 +29,8 @@ } while (0) #define DD_set_z(a,b) a[0] = VALUE_TO_DOUBLE(b) +#define DD_floor(a,b) value_set_si(a,(int)floor(b[0]+ddf_almostzero)) +#define DD_ceil(a,b) value_set_si(a,(int)ceil(b[0]-ddf_almostzero)) #define cdd_polyhedron_affine_sign cddf_polyhedron_affine_sign #define DD_rat_sign(sign,obj,val) do { \ @@ -36,5 +40,6 @@ val[0] += ddf_almostzero; \ sign = ddf_sgn(val); \ } while (0) +#define cdd_polyhedron_range cddf_polyhedron_range #include "polysign_cdd_template.cc" diff --git a/polysign_polylib.c b/polysign_polylib.c index 19e512b..2af14dd 100644 --- a/polysign_polylib.c +++ b/polysign_polylib.c @@ -50,3 +50,31 @@ enum order_sign PL_polyhedron_affine_sign(Polyhedron *D, Matrix *T, Polyhedron_Free(I); return sign; } + +enum lp_result PL_polyhedron_range(Polyhedron *D, Value *obj, Value denom, + Value *min, Value *max, + struct barvinok_options *options) +{ + Polyhedron *I; + int bounded; + + Matrix *T = Matrix_Alloc(2, D->Dimension+1); + Vector_Copy(obj, T->p[0], D->Dimension+1); + value_assign(T->p[1][D->Dimension], denom); + + I = Polyhedron_Image(D, T, options->MaxRays); + Matrix_Free(T); + POL_ENSURE_VERTICES(I); + + if (emptyQ(I)) { + Polyhedron_Free(I); + return lp_empty; + } + bounded = line_minmax(I, min, max); + + if (!bounded) + return lp_unbounded; + if (value_gt(*min, *max)) + return lp_empty; + return lp_ok; +} diff --git a/sample.c b/sample.c index a149c68..4cbf459 100644 --- a/sample.c +++ b/sample.c @@ -2,6 +2,7 @@ #include #include #include +#include "polysign.h" #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type)) @@ -242,12 +243,13 @@ static Polyhedron *Polyhedron_RemoveFixedColumns(Polyhedron *P, Matrix **T) Vector *Polyhedron_Sample(Polyhedron *P, struct barvinok_options *options) { int i, j; - Vector *sample = NULL; + Vector *sample = NULL, *obj = NULL; Polyhedron *Q; - Matrix *T, *inv, *M; + Matrix *inv = NULL; Value min, max, tmp; Vector *v; int ok; + enum lp_result res; POL_ENSURE_VERTICES(P); if (emptyQ(P)) @@ -273,41 +275,41 @@ Vector *Polyhedron_Sample(Polyhedron *P, struct barvinok_options *options) return sample; } - Matrix *basis = Polyhedron_Reduced_Basis(P, options); - - T = Matrix_Alloc(P->Dimension+1, P->Dimension+1); - inv = Matrix_Alloc(P->Dimension+1, P->Dimension+1); - for (i = 0; i < P->Dimension; ++i) - for (j = 0; j < P->Dimension; ++j) - value_assign(T->p[i][j], basis->p[i][j]); - value_set_si(T->p[P->Dimension][P->Dimension], 1); - Matrix_Free(basis); - - M = Matrix_Copy(T); - ok = Matrix_Inverse(M, inv); - assert(ok); - Matrix_Free(M); - - Q = Polyhedron_Image(P, T, options->MaxRays); - - POL_ENSURE_VERTICES(Q); - value_init(min); value_init(max); - value_init(tmp); - mpz_cdiv_q(min, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]); - mpz_fdiv_q(max, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]); + obj = Vector_Alloc(P->Dimension+1); + value_set_si(obj->p[0], 1); + res = polyhedron_range(P, obj->p, obj->p[0], &min, &max, options); + assert(res == lp_ok); - for (j = 1; j < Q->NbRays; ++j) { - mpz_cdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]); - if (value_lt(tmp, min)) - value_assign(min, tmp); - mpz_fdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]); - if (value_gt(tmp, max)) - value_assign(max, tmp); + if (value_eq(min, max)) { + Q = P; + } else { + Matrix *basis = Polyhedron_Reduced_Basis(P, options); + Matrix *M; + Matrix *T = Matrix_Alloc(P->Dimension+1, P->Dimension+1); + inv = Matrix_Alloc(P->Dimension+1, P->Dimension+1); + for (i = 0; i < P->Dimension; ++i) + for (j = 0; j < P->Dimension; ++j) + value_assign(T->p[i][j], basis->p[i][j]); + value_set_si(T->p[P->Dimension][P->Dimension], 1); + Matrix_Free(basis); + + M = Matrix_Copy(T); + ok = Matrix_Inverse(M, inv); + assert(ok); + Matrix_Free(M); + + Q = Polyhedron_Image(P, T, options->MaxRays); + res = polyhedron_range(Q, obj->p, obj->p[0], &min, &max, options); + assert(res == lp_ok); + + Matrix_Free(T); } + value_init(tmp); + v = Vector_Alloc(1+Q->Dimension+1); value_set_si(v->p[1], -1); @@ -329,21 +331,29 @@ Vector *Polyhedron_Sample(Polyhedron *P, struct barvinok_options *options) S_sample = Polyhedron_Sample(S, options); Polyhedron_Free(S); if (S_sample) { - Vector *Q_sample = Vector_Alloc(Q->Dimension + 1); + Vector *Q_sample = obj; + obj = NULL; Matrix_Vector_Product(T, S_sample->p, Q_sample->p); Matrix_Free(T); Vector_Free(S_sample); - sample = Vector_Alloc(P->Dimension + 1); - Matrix_Vector_Product(inv, Q_sample->p, sample->p); - Vector_Free(Q_sample); + if (!inv) + sample = Q_sample; + else { + sample = Vector_Alloc(P->Dimension + 1); + Matrix_Vector_Product(inv, Q_sample->p, sample->p); + Vector_Free(Q_sample); + } break; } Matrix_Free(T); } - Matrix_Free(T); - Matrix_Free(inv); - Polyhedron_Free(Q); + if (obj) + Vector_Free(obj); + if (inv) + Matrix_Free(inv); + if (P != Q) + Polyhedron_Free(Q); Vector_Free(v); value_clear(min); -- 2.11.4.GIT