Module: ppl/ppl Branch: pip Commit: 9f62791f13a0de4decf334ca5b3d604207bac6f7 URL: http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=9f62791f13a0d... Author: François Galea <francois.galea@uvsq.fr> Date: Tue Oct 20 18:18:19 2009 +0200 Parameter compatibility check now applies a revised dual simplex method. --- src/PIP_Tree.cc | 48 +++++++++++++++++++++++++++++++++++++++++------- src/PIP_Tree.defs.hh | 8 ++++---- 2 files changed, 45 insertions(+), 11 deletions(-) diff --git a/src/PIP_Tree.cc b/src/PIP_Tree.cc index fff5599..355f8f1 100644 --- a/src/PIP_Tree.cc +++ b/src/PIP_Tree.cc @@ -688,14 +688,23 @@ bool PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { Matrix s(ctx); s.add_row(cnst); - dimension_type i, i_, j, k, j_, j__; + dimension_type i, i_, j, k, j_, j__, var_i, var_j; dimension_type num_rows = s.num_rows(); dimension_type num_cols = s.num_columns(); + dimension_type num_vars = num_cols-1; + const dimension_type n_a_d = not_a_dimension(); std::vector<Coefficient> scaling(num_rows, 1); PPL_DIRTY_TEMP_COEFFICIENT(sij); PPL_DIRTY_TEMP_COEFFICIENT(mult); PPL_DIRTY_TEMP_COEFFICIENT(gcd); PPL_DIRTY_TEMP_COEFFICIENT(scale_factor); + PPL_DIRTY_TEMP_COEFFICIENT(scaling_i); + std::vector<dimension_type> mapping; + std::vector<bool> basis(num_vars, true); + for (j=1; j<=num_vars; ++j) + mapping.push_back(j); + Row p(num_cols, compute_capacity(num_cols, Matrix::max_num_columns()), + Row::Flags()); /* Perform simplex pivots on the context until we find an empty solution * or an optimum */ @@ -726,9 +735,37 @@ PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { // No negative RHS: optimum found return true; } + // Now we have a positive s[i][j] pivot - const Row& p = s[i]; + var_j = n_a_d; + var_i = n_a_d; + for (j_=0; j_<num_vars; ++j_) { + if (basis[j_]) { + if (mapping[j_] == j) + var_j = j_; + } else { + if (mapping[j_] == i) + var_i = j_; + } + } + /* update basis */ + if (var_j != n_a_d) { + basis[var_j] = false; + mapping[var_j] = i; + } + if (var_i != n_a_d) { + basis[var_i] = true; + mapping[var_i] = j; + } + + /* create the identity matrix row corresponding to basic variable j */ + for (j_=0; j_<num_cols; ++j_) + p[j_]=0; + p[j] = 1; + p.swap(s[i]); sij = p[j]; + scaling_i = scaling[i]; + scaling[i] = 1; // Perform a pivot operation on the matrix for (j_=0; j_<num_cols; ++j_) { @@ -739,8 +776,6 @@ PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { // 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; Row& row = s[k]; mult = row[j] * sij_; if (mult % sij != 0) { @@ -754,14 +789,13 @@ PIP_Solution_Node::compatibility_check(const Matrix &ctx, const Row &cnst) { } row[j_] -= mult / sij; } - s[i][j_] = 0; } - if (sij != scaling[i]) { + 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]; - mult = skj*scaling[i]; + mult = skj*scaling_i; if (mult % sij != 0) { // as above, we must perform row scaling gcd_assign(gcd, mult, sij); diff --git a/src/PIP_Tree.defs.hh b/src/PIP_Tree.defs.hh index 4141818..c2a8186 100644 --- a/src/PIP_Tree.defs.hh +++ b/src/PIP_Tree.defs.hh @@ -331,10 +331,10 @@ private: Checks whether a constraint is compatible with a context, ie. does not make the context empty. - The algorithm consists in performing simplex pivots on a Matrix consisting - in the original matrix with the constraint inserted. If the simplex - terminates with a solution, then the restrained context is not empty. - Otherwise, it is. + The method consists in applying the revised dual simplex algorithm on a + Matrix consisting in the original matrix with the constraint inserted. If + the simplex terminates with a feasible (optimal) solution, then the + restrained context is not empty. Otherwise, it is. */ static bool compatibility_check(const Matrix &ctx, const Row &cnst);