October 2009 - PPL-devel - BUGSENG Mailman

[GIT] ppl/ppl(pip): Better use of scaling and normalization to keep low coefficient values
by François Galea 16 Oct '09

16 Oct '09

Module: ppl/ppl Branch: pip Commit: 65b630d43bcbe16e2fb0bbde64a73f9c29b781d6 URL: http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=65b630d43bcb… Author: François Galea <francois.galea(a)uvsq.fr> Date: Fri Oct 16 09:38:10 2009 +0200 Better use of scaling and normalization to keep low coefficient values where possible. --- src/PIP_Tree.cc | 48 +++++++++++++++++++++++++++++++++++++++--------- 1 files changed, 39 insertions(+), 9 deletions(-) diff --git a/src/PIP_Tree.cc b/src/PIP_Tree.cc index 5199b41..9f09443 100644 --- a/src/PIP_Tree.cc +++ b/src/PIP_Tree.cc @@ -74,7 +74,9 @@ negate_assign(Row& x, const Row& y, const Coefficient& sc) { PPL_ASSERT(x.size() == y.size()); for (dimension_type i = x.size(); i-- > 0; ) x[i] = -y[i]; - x[0] -= sc; + PPL_DIRTY_TEMP_COEFFICIENT(mod); + mod_assign(mod, x[0], sc); + x[0] -= ((mod == 0) ? sc : mod); } // Update given context matrix using local artificials @@ -99,6 +101,29 @@ update_context(Variables_Set &params, Matrix &context, } } +void +row_normalize(Row& x) { + dimension_type size = x.size(); + if (size == 0) + return; + dimension_type j; + PPL_DIRTY_TEMP_COEFFICIENT(gcd); + gcd = x[0]; + for (j=1; j<size; ++j) { + const Coefficient& c = x[j]; + if (c != 0) { + gcd_assign(gcd, c, gcd); + if (gcd == 1) + return; + } + } + // Divide the coefficients by the GCD. + for (j=0; j<size; ++j) { + Coefficient& c = x[j]; + exact_div_assign(c, c, gcd); + } +} + } // namespace namespace IO_Operators { @@ -936,7 +961,10 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, if (!found) continue; Row row(tableau.t[i]); - row[0] -= tableau.get_denominator(); + const Coefficient& denom = tableau.get_denominator(); + PPL_DIRTY_TEMP_COEFFICIENT(mod); + mod_assign(mod, row[0], denom); + row[0] -= ((mod == 0) ? denom : mod); if (compatibility_check(context, row)) { if (i__ == n_a_d) i__ = i; @@ -1159,7 +1187,8 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, "variable coefficients: " << i << std::endl; #endif - Row &r = tableau.t[i]; + Row r(tableau.t[i]); + row_normalize(r); context.add_row(r); add_constraint(r, parameters); sign[i] = POSITIVE; @@ -1191,6 +1220,8 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, } i__ = best_i; + Row test(tableau.t[i__]); + row_normalize(test); #ifdef NOISY_PIP { using namespace IO_Operators; @@ -1198,14 +1229,13 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, Variables_Set::const_iterator p; dimension_type j; for (p = parameters.begin(), j=1; p != parameters.end(); ++p, ++j) - e += tableau.t[i__][j] * Variable(*p); - e += tableau.t[i__][0]; + e += test[j] * Variable(*p); + e += test[0]; std::cout << "Found row with mixed parameter sign: " << i__ << "\nSolution depends on the sign of parameter " << e << std::endl; } #endif - Row test(tableau.t[i__]); /* Create a solution Node to become "true" version of current Node */ PIP_Tree_Node *tru = new PIP_Solution_Node(*this, true); @@ -1392,19 +1422,19 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, #endif context.add_row(ctx1); context.add_row(ctx2); - cut_t[num_params] = d*d; + cut_t[num_params] = d; } // Generate new cut const Row& row_s = tableau.s[i]; for (j=0; j<num_vars; ++j) { mod_assign(mod, row_s[j], d); - cut_s[j] = d*mod; + cut_s[j] = mod; } for (j=0; j<num_params; ++j) { mod_assign(mod, row_t[j], d); if (mod != 0) - cut_t[j] = d*(mod - d); + cut_t[j] = mod - d; else cut_t[j] = 0; }

1 0

[GIT] ppl/ppl(pip): Corrected invalid compatibility check algorithm.
by François Galea 16 Oct '09

16 Oct '09

Module: ppl/ppl Branch: pip Commit: cd38635a96c022e5705aa7fdcbef65c7dc5e8295 URL: http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=cd38635a96c0… Author: François Galea <francois.galea(a)uvsq.fr> Date: Fri Oct 16 09:32:29 2009 +0200 Corrected invalid compatibility check algorithm. --- src/PIP_Tree.cc | 32 +++++++++++++++++++------------- tests/PIP_Problem/pipproblem1.cc | 2 +- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/PIP_Tree.cc b/src/PIP_Tree.cc index 833eb5e..5199b41 100644 --- a/src/PIP_Tree.cc +++ b/src/PIP_Tree.cc @@ -666,7 +666,12 @@ PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { dimension_type i, j, k, j_, j__; dimension_type num_rows = s.num_rows(); dimension_type num_cols = s.num_columns(); + std::vector<Coefficient> scaling(num_rows, 1); bool result = false; + PPL_DIRTY_TEMP_COEFFICIENT(sij); + PPL_DIRTY_TEMP_COEFFICIENT(mult); + PPL_DIRTY_TEMP_COEFFICIENT(gcd); + PPL_DIRTY_TEMP_COEFFICIENT(scale_factor); /* Perform simplex pivots on the context until we find an empty solution * or an optimum */ @@ -682,7 +687,7 @@ PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { } // Find a positive m[i][j] pivot j = 1; - Row &p = s[i]; + const Row& p = s[i]; while (j<num_cols && p[j] <= 0) j++; if (j == num_cols) { @@ -691,11 +696,11 @@ PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { break; } // Perform a pivot operation on the matrix - const Coefficient sij = p[j]; + sij = p[j]; for (j_=0; j_<num_cols; ++j_) { if (j_ == j) continue; - const Coefficient sij_ = p[j_]; + const Coefficient& sij_ = p[j_]; if (sij_ == 0) // if element j of pivot row is zero, nothing to do for this column continue; @@ -703,35 +708,36 @@ PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { if (k == i) continue; Row& row = s[k]; - Coefficient mult = row[j] * sij_; + mult = row[j] * sij_; if (mult % sij != 0) { // Must scale row to stay in integer case - Coefficient gcd; gcd_assign(gcd, mult, sij); - Coefficient scale_factor = sij/gcd; + exact_div_assign(scale_factor, sij, gcd); for (j__=0; j__<num_cols; ++j__) row[j__] *= scale_factor; mult *= scale_factor; + scaling[k] *= scale_factor; } row[j_] -= mult / sij; } s[i][j_] = 0; } - if (sij != 1) { + if (sij != scaling[i]) { // Update column only if pivot != 1 for (k=0; k<num_rows; ++k) { Row& row = s[k]; Coefficient& skj = row[j]; - if (skj % sij != 0) { + mult = skj*scaling[i]; + if (mult % sij != 0) { // as above, we must perform row scaling - Coefficient gcd; - gcd_assign(gcd, skj, sij); - Coefficient scale_factor = sij/gcd; + gcd_assign(gcd, mult, sij); + exact_div_assign(scale_factor, sij, gcd); for (j__=0; j__<num_cols; ++j__) row[j__] *= scale_factor; - skj *= scale_factor; + scaling[k] *= scale_factor; + mult *= scale_factor; } - skj /= sij; + exact_div_assign(skj, mult, sij); } } } diff --git a/tests/PIP_Problem/pipproblem1.cc b/tests/PIP_Problem/pipproblem1.cc index be6f68f..5b436fe 100644 --- a/tests/PIP_Problem/pipproblem1.cc +++ b/tests/PIP_Problem/pipproblem1.cc @@ -236,5 +236,5 @@ BEGIN_MAIN DO_TEST(test02); DO_TEST(test03); DO_TEST(test04); - //DO_TEST(test05); + DO_TEST(test05); END_MAIN

1 0

[GIT] ppl/ppl(floating_point): Added float_ibm_single declaration into class Float<T>.
by Fabio Biselli 15 Oct '09

15 Oct '09

Module: ppl/ppl Branch: floating_point Commit: 9c6033597683d6b3ccc16aa443a6c8e20354d872 URL: http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=9c6033597683… Author: Fabio Biselli <fabio.biselli(a)studenti.unipr.it> Date: Thu Oct 15 11:54:43 2009 +0200 Added float_ibm_single declaration into class Float<T>. --- src/Float.defs.hh | 6 ++++++ 1 files changed, 6 insertions(+), 0 deletions(-) diff --git a/src/Float.defs.hh b/src/Float.defs.hh index 303e988..f0eb5f8 100644 --- a/src/Float.defs.hh +++ b/src/Float.defs.hh @@ -275,6 +275,8 @@ public: typedef float_ieee754_single Binary; #elif PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_IEEE754_DOUBLE typedef float_ieee754_double Binary; +#elif PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_IBM_SINGLE + typedef float_ibm_single Binary; #elif PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_IEEE754_QUAD typedef float_ieee754_quad Binary; #elif PPL_CXX_FLOAT_BINARY_FORMAT == PPL_FLOAT_INTEL_DOUBLE_EXTENDED @@ -302,6 +304,8 @@ public: typedef float_ieee754_single Binary; #elif PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_DOUBLE typedef float_ieee754_double Binary; +#elif PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IBM_SINGLE + typedef float_ibm_single Binary; #elif PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_QUAD typedef float_ieee754_quad Binary; #elif PPL_CXX_DOUBLE_BINARY_FORMAT == PPL_FLOAT_INTEL_DOUBLE_EXTENDED @@ -329,6 +333,8 @@ public: typedef float_ieee754_single Binary; #elif PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_DOUBLE typedef float_ieee754_double Binary; +#elif PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IBM_SINGLE + typedef float_ibm_single Binary; #elif PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_IEEE754_QUAD typedef float_ieee754_quad Binary; #elif PPL_CXX_LONG_DOUBLE_BINARY_FORMAT == PPL_FLOAT_INTEL_DOUBLE_EXTENDED

1 0

[Fwd: Re: warning GMP -fexception message on configure]
by Roberto Bagnara 14 Oct '09

14 Oct '09

-------- Original Message -------- Subject: Re: [PPL-devel] warning GMP -fexception message on configure Date: Wed, 14 Oct 2009 08:15:30 -0700 From: Reimundo Heluani <rheluani(a)gmail.com> To: Roberto Bagnara <bagnara(a)cs.unipr.it> References: <1255460135.18699.37.camel(a)campari.homeip.net> <4AD596A2.4090800(a)cs.unipr.it> On Wed, 2009-10-14 at 11:15 +0200, Roberto Bagnara wrote: > Reimundo Heluani wrote: > > Hi there, trying to compile PPL-0.10.2 I get the message > > > > > > WARNING: CANNOT PROPAGATE EXCEPTIONS BACK FROM GMP: > > *** MEMORY EXHAUSTION MAY RESULT IN ABRUPT TERMINATION > > *** This is OK, if you do not plan to use the bounded memory > > capabilities > > *** offered by the PPL. Otherwise, if you are using GCC or the Intel; > > C/C++ compiler, please make sure you use a version of GMP compiled with > > the -fexceptions compiler option. > > ... > > .... > > > > And the end of the configure round. > > > > I just compiled (twice to make sure) GMP 4.3.1 with > > CPPFLAGS=-fexceptions --enable-cxx > > > > I'm trying to compile this on a CFLS x86_64 system. Is is normal? > > Hi Reimundo, > > from your message I gather you used, for GMP 4.3.1, a configuration > command of the form > > ./configure CPPFLAGS=-fexceptions --enable-cxx > > I would use a configuration command of this form instead > > CPPFLAGS=-fexceptions ./configure --enable-cxx Indeed the command I have used is of this form. I'm following CLFS and in that list someone told me that the compilation of GMP goes fine, just that the test on PPL is failing. Indeed in the final version of the system I do not get that warning. Just in the cross compiled one. Cheers, R. > > Please let us know how it goes. > All the best, > > Roberto > -- Prof. Roberto Bagnara Computer Science Group Department of Mathematics, University of Parma, Italy http://www.cs.unipr.it/~bagnara/ mailto:bagnara@cs.unipr.it

1 0

[GIT] ppl/ppl(floating_point): Added Floating_Point_Expression.templates. hh into documentation.
by Fabio Biselli 14 Oct '09

14 Oct '09

Module: ppl/ppl Branch: floating_point Commit: cacb4330a7f4af12093b5191ff6f453eba0f589a URL: http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=cacb4330a7f4… Author: Fabio Biselli <fabio.biselli(a)studenti.unipr.it> Date: Wed Oct 14 18:06:40 2009 +0200 Added Floating_Point_Expression.templates.hh into documentation. --- doc/devref.doxyconf-html.in | 1 + 1 files changed, 1 insertions(+), 0 deletions(-) diff --git a/doc/devref.doxyconf-html.in b/doc/devref.doxyconf-html.in index 7ecb038..17e02d1 100644 --- a/doc/devref.doxyconf-html.in +++ b/doc/devref.doxyconf-html.in @@ -317,6 +317,7 @@ INPUT = @srcdir@/definitions.dox \ @srcdir@/../src/Floating_Point_Expression.types.hh \ @srcdir@/../src/Floating_Point_Expression.defs.hh \ @srcdir@/../src/Floating_Point_Expression.inlines.hh \ + @srcdir@/../src/Floating_Point_Expression.templates.hh \ @srcdir@/../src/Constant_Floating_Point_Expression.types.hh \ @srcdir@/../src/Constant_Floating_Point_Expression.defs.hh \ @srcdir@/../src/Constant_Floating_Point_Expression.inlines.hh \

1 0

[GIT] ppl/ppl(pip): Added a test case illustrating an instability in the PIP solver.
by François Galea 14 Oct '09

14 Oct '09

Module: ppl/ppl Branch: pip Commit: 051c82f85918f442a031213156a288628c2be2c0 URL: http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=051c82f85918… Author: François Galea <francois.galea(a)uvsq.fr> Date: Wed Oct 14 12:28:29 2009 +0200 Added a test case illustrating an instability in the PIP solver. --- tests/PIP_Problem/pipproblem1.cc | 28 ++++++++++++++++++++++++++++ 1 files changed, 28 insertions(+), 0 deletions(-) diff --git a/tests/PIP_Problem/pipproblem1.cc b/tests/PIP_Problem/pipproblem1.cc index 191727c..be6f68f 100644 --- a/tests/PIP_Problem/pipproblem1.cc +++ b/tests/PIP_Problem/pipproblem1.cc @@ -202,6 +202,33 @@ test04() { return ok; } +bool +test05() { + Variable i(0); + Variable j(1); + Variable m(2); + Variable n(3); + Variables_Set params(m, n); + + Constraint_System cs; + cs.insert(3*j >= -2*i+8); + cs.insert(j <= 4*i - 4); + cs.insert(i <= n); + cs.insert(j <= m); + cs.insert(n >= 3); + + PIP_Problem pip(cs.space_dimension(), cs.begin(), cs.end(), params); + + bool ok = (pip.solve() == OPTIMIZED_PIP_PROBLEM); + if (ok) { + const PIP_Tree solution = pip.solution(); + display_solution(solution, params, Variables_Set(i, j), + cs.space_dimension()); + } + + return ok; +} + } // namespace BEGIN_MAIN @@ -209,4 +236,5 @@ BEGIN_MAIN DO_TEST(test02); DO_TEST(test03); DO_TEST(test04); + //DO_TEST(test05); END_MAIN

1 0

warning GMP -fexception message on configure
by Reimundo Heluani 14 Oct '09

14 Oct '09

Hi there, trying to compile PPL-0.10.2 I get the message WARNING: CANNOT PROPAGATE EXCEPTIONS BACK FROM GMP: *** MEMORY EXHAUSTION MAY RESULT IN ABRUPT TERMINATION *** This is OK, if you do not plan to use the bounded memory capabilities *** offered by the PPL. Otherwise, if you are using GCC or the Intel; C/C++ compiler, please make sure you use a version of GMP compiled with the -fexceptions compiler option. ... .... And the end of the configure round. I just compiled (twice to make sure) GMP 4.3.1 with CPPFLAGS=-fexceptions --enable-cxx I'm trying to compile this on a CFLS x86_64 system. Is is normal? Thanks R.

2 1

[GIT] ppl/ppl(floating_point): Modified widening in test05 and test06.
by Roberto Amadini 13 Oct '09

13 Oct '09

Module: ppl/ppl Branch: floating_point Commit: 286faa55e60189534de6fb214ec8518aa45ebb63 URL: http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=286faa55e601… Author: Roberto Amadini <r.amadini(a)virgilio.it> Date: Tue Oct 13 15:42:32 2009 +0200 Modified widening in test05 and test06. --- tests/Floating_Point_Expression/digitalfilters1.cc | 44 ++++++++++++-------- 1 files changed, 26 insertions(+), 18 deletions(-) diff --git a/tests/Floating_Point_Expression/digitalfilters1.cc b/tests/Floating_Point_Expression/digitalfilters1.cc index ae586bb..4dabd47 100644 --- a/tests/Floating_Point_Expression/digitalfilters1.cc +++ b/tests/Floating_Point_Expression/digitalfilters1.cc @@ -129,7 +129,7 @@ test01() { return !y.is_bounded(); } -// tests rate limiter using bounded differences abstract domain +// Tests rate limiter using bounded differences abstract domain // and ignoring rounding errors. bool test02() { @@ -219,7 +219,7 @@ test02() { return !y.is_bounded(); } -// tests rate limiter using octagons abstract domain +// Tests rate limiter using octagons abstract domain // and ignoring rounding errors. bool test03() { @@ -315,7 +315,7 @@ test03() { return y.is_bounded(); } -// tests rate limiter using polyhedra abstract domain +// Tests rate limiter using polyhedra abstract domain // and ignoring rounding errors. bool test04() { @@ -410,7 +410,7 @@ test04() { return y.is_bounded(); } -// tests rate limiter using octagons abstract domain and +// Tests rate limiter using octagons abstract domain and // linearization of floating point expressions. bool test05() { @@ -503,6 +503,7 @@ test05() { oc.upper_bound_assign(oc_begin); Constraint_System cs; PPL_DIRTY_TEMP_COEFFICIENT(M); + PPL_DIRTY_TEMP_COEFFICIENT(N); ANALYZER_FP_FORMAT max_analyzed = (2 - pow(2, -static_cast<ANALYZER_FP_FORMAT>(ANALYZED_FP_FORMAT::MANTISSA_BITS))) @@ -512,21 +513,26 @@ test05() { std::numeric_limits<ANALYZER_FP_FORMAT>::max(); max_analyzer = std::min(max_analyzer, max_analyzed); assign_r(M, max_analyzer, ROUND_DOWN); - cs.insert(Y <= M); - cs.insert(Y >= -M); - oc.limited_BHMZ05_extrapolation_assign(oc_begin, cs); - Box<FP_Interval> box(oc); + div_2exp_assign_r(N, M, 1, ROUND_DOWN); + // FIXME: By inserting the constraints Y <= M and Y >= -M, we obtain + // Y + Y <= 2 * M = +inf and -Y - Y <= -2 * M = +inf. + // For a more precise analysis, it is better to insert the + // constraints Y <= N and Y >= -N, where N = M / 2. + cs.insert(Y <= N); + cs.insert(Y >= -N); + FP_Octagonal_Shape oc_wid(oc); + oc_wid.limited_BHMZ05_extrapolation_assign(oc_begin, cs); + Box<FP_Interval> box(oc_wid); print_constraints(box, "*** after widening ***"); y = box.get_interval(Y); } nout << "*** Y in " << y << " ***" << endl; - return !y.is_bounded(); + return y.is_bounded(); } -// tests rate limiter using polyhedra abstract domain and +// Tests rate limiter using polyhedra abstract domain and // linearization of floating point expressions. -// FIXME: not pass at the moment. bool test06() { // Input signal. @@ -557,11 +563,12 @@ test06() { con_y.linearize(abstract_store, lf_abstract_store, lk); ph.affine_image(Y, lk); - for(unsigned short n = 0; y_begin != y && y.is_bounded(); ++n) { + for(unsigned short n = 0; y_begin != y; ++n) { nout << "*** n = " << n << " ***" << endl; ph_begin = ph; y_begin = y; + // X = [-128, 128]; D = [0, 16]; S = Y; R = X - S; Y = X; Con_FP_Expression con_x(-128, 128); con_x.linearize(abstract_store, lf_abstract_store, lk); @@ -633,23 +640,24 @@ test06() { assign_r(M, max_analyzer, ROUND_DOWN); cs.insert(Y <= M); cs.insert(Y >= -M); - ph.limited_BHRZ03_extrapolation_assign(ph_begin, cs); - Box<FP_Interval> box(ph); + NNC_Polyhedron ph_wid(ph); + ph_wid.limited_BHRZ03_extrapolation_assign(ph_begin, cs); + Box<FP_Interval> box(ph_wid); print_constraints(box, "*** after widening ***"); y = box.get_interval(Y); } nout << "*** Y in " << y << " ***" << endl; - return !y.is_bounded(); + return y.is_bounded(); } } // namespace -BEGIN_MAIN +BEGIN_MAIN /* DO_TEST(test01); DO_TEST(test02); DO_TEST(test03); - DO_TEST(test04); + DO_TEST(test04); */ DO_TEST(test05); - //DO_TEST(test06); + DO_TEST(test06); END_MAIN

1 0

[GIT] ppl/ppl(pip): Integrality check only on the original variables ( not slack ones).
by François Galea 13 Oct '09

13 Oct '09

Module: ppl/ppl Branch: pip Commit: 8e8d79ad7cef965e68277a902d51a9dcfe8b7d8a URL: http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=8e8d79ad7cef… Author: François Galea <francois.galea(a)uvsq.fr> Date: Tue Oct 13 10:50:48 2009 +0200 Integrality check only on the original variables (not slack ones). --- src/PIP_Tree.cc | 13 ++++++++----- 1 files changed, 8 insertions(+), 5 deletions(-) diff --git a/src/PIP_Tree.cc b/src/PIP_Tree.cc index 91406f7..833eb5e 100644 --- a/src/PIP_Tree.cc +++ b/src/PIP_Tree.cc @@ -1265,8 +1265,11 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, // Look for a row with non integer parameter coefficients (first is okay) const Coefficient& d = tableau.get_denominator(); - for (i=0; i<num_rows; ++i) { - const Row& row = tableau.t[i]; + for (i=0; i<num_vars; ++i) { + if (basis[i]) + // basic variable = 0 -> integer + continue; + const Row& row = tableau.t[mapping[i]]; for (j=0; j<num_params; ++j) { if (row[j] % d != 0) goto endsearch; @@ -1274,7 +1277,7 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, } endsearch: - if (i == num_rows) { + if (i == num_vars) { /* The solution is integer */ #ifdef NOISY_PIP std::cout << "Solution found for problem in current node." @@ -1288,8 +1291,8 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, /* Look for row which will generate the "deepest" cut */ Coefficient score; Coefficient best = 0; - dimension_type best_i = i; - for (i_ = i; i_ < num_rows; ++i_) { + dimension_type best_i = 0; + for (i_ = 0; i_ < num_rows; ++i_) { const Row& row = tableau.t[i_]; score = 0; for (j = 0; j < num_params; ++j) {

1 0

[GIT] ppl/ppl(pip): Optimized pivot operation.
by François Galea 13 Oct '09

13 Oct '09

Module: ppl/ppl Branch: pip Commit: 386c5e6645338fbfb50271b82e9c6f621ec27280 URL: http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=386c5e664533… Author: François Galea <francois.galea(a)uvsq.fr> Date: Tue Oct 13 10:49:20 2009 +0200 Optimized pivot operation. --- src/PIP_Tree.cc | 65 ++++++++++++++++++++++++++++++++++-------------------- 1 files changed, 41 insertions(+), 24 deletions(-) diff --git a/src/PIP_Tree.cc b/src/PIP_Tree.cc index 7a470dd..91406f7 100644 --- a/src/PIP_Tree.cc +++ b/src/PIP_Tree.cc @@ -696,6 +696,9 @@ PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { if (j_ == j) continue; const Coefficient sij_ = p[j_]; + if (sij_ == 0) + // if element j of pivot row is zero, nothing to do for this column + continue; for (k=0; k<num_rows; ++k) { if (k == i) continue; @@ -714,19 +717,22 @@ PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { } s[i][j_] = 0; } - for (k=0; k<num_rows; ++k) { - Row& row = s[k]; - Coefficient& skj = row[j]; - if (skj % sij != 0) { - // as above, we must perform row scaling - Coefficient gcd; - gcd_assign(gcd, skj, sij); - Coefficient scale_factor = sij/gcd; - for (j__=0; j__<num_cols; ++j__) - row[j__] *= scale_factor; - skj *= scale_factor; + if (sij != 1) { + // Update column only if pivot != 1 + for (k=0; k<num_rows; ++k) { + Row& row = s[k]; + Coefficient& skj = row[j]; + if (skj % sij != 0) { + // as above, we must perform row scaling + Coefficient gcd; + gcd_assign(gcd, skj, sij); + Coefficient scale_factor = sij/gcd; + for (j__=0; j__<num_cols; ++j__) + row[j__] *= scale_factor; + skj *= scale_factor; + } + skj /= sij; } - skj /= sij; } } @@ -1036,8 +1042,12 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, for (j=0; j<num_vars; ++j) { if (j==j_) continue; + const Coefficient& prsj = prs[j]; + if (prsj == 0) + // if element j of pivot row is zero, nothing to do for this column + continue; for (k=0; k<num_rows; ++k) { - Coefficient mult = prs[j] * tableau.s[k][j_]; + Coefficient mult = prsj * tableau.s[k][j_]; if (mult % sij != 0) { // Must scale matrix to stay in integer case gcd_assign(gcd, mult, sij); @@ -1051,8 +1061,12 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, /* Compute columns t[*][j] : t[k][j] -= t[k][j_] * prt[j] / sij */ for (j=0; j<num_params; ++j) { + const Coefficient& prtj = prt[j]; + if (prtj == 0) + // if element j of pivot row is zero, nothing to do for this column + continue; for (k=0; k<num_rows; ++k) { - Coefficient c = prt[j] * tableau.s[k][j_]; + Coefficient c = prtj * tableau.s[k][j_]; if (c % sij != 0) { // Must scale matrix to stay in integer case gcd_assign(gcd, c, sij); @@ -1088,17 +1102,20 @@ PIP_Solution_Node::solve(PIP_Tree_Node*& parent_ref, const Matrix& ctx, } /* compute column s[*][j_] : s[k][j_] /= sij */ - for (k=0; k<num_rows; ++k) { - Coefficient& c = tableau.s[k][j_]; - Coefficient numerator = c * sij_denom; - if (numerator % sij != 0) { - Coefficient gcd; - gcd_assign(gcd, numerator, sij); - Coefficient scale_factor = sij/gcd; - tableau.scale(scale_factor); - numerator *= scale_factor; + if (sij != sij_denom) { + // Update column only if pivot != 1 + for (k=0; k<num_rows; ++k) { + Coefficient& c = tableau.s[k][j_]; + Coefficient numerator = c * sij_denom; + if (numerator % sij != 0) { + Coefficient gcd; + gcd_assign(gcd, numerator, sij); + Coefficient scale_factor = sij/gcd; + tableau.scale(scale_factor); + numerator *= scale_factor; + } + c = numerator / sij; } - c = numerator / sij; } solution_valid = false; }

1 0