[PPL-devel] [GIT] ppl/ppl(pip): Implemented simplex pivot on rational matrices.

[François Galea]

Module: ppl/ppl
Branch: pip
Commit: b8e1cbd1782c3fa553331b6199668f4bb26581b1
URL:    http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=b8e1cbd1782c3fa553331b6199668f4bb26581b1

Author: François Galea <francois.galea at uvsq.fr>
Date:   Tue Sep 15 15:09:08 2009 +0200

Implemented simplex pivot on rational matrices.

---

 src/PIP_Tree.cc      |  111 +++++++++++++++++++++++++++++++++++++++++++++++++-
 src/PIP_Tree.defs.hh |    3 +
 2 files changed, 112 insertions(+), 2 deletions(-)

diff --git a/src/PIP_Tree.cc b/src/PIP_Tree.cc
index 173d39f..c9f576a 100644
--- a/src/PIP_Tree.cc
+++ b/src/PIP_Tree.cc
@@ -253,6 +253,16 @@ PIP_Solution_Node::Rational_Matrix::normalize() {
 }
 
 void
+PIP_Solution_Node::Rational_Matrix::scale(const Coefficient &ratio) {
+  dimension_type i, j;
+  dimension_type i_max = num_rows();
+  dimension_type j_max = num_columns();
+  for (i=0; i<i_max; ++i)
+    for (j=0; j<j_max; ++j)
+      rows[i][j] *= ratio;
+}
+
+void
 PIP_Solution_Node::Rational_Matrix::ascii_dump(std::ostream& s) const {
   s << "denominator " << denominator << "\n";
   Matrix::ascii_dump(s);
@@ -472,6 +482,7 @@ PIP_Solution_Node::solve(PIP_Tree_Node** parent_ref,
   Matrix context(ctx);
   merge_assign(context, constraints_);
   const dimension_type n_a_d = not_a_dimension();
+  Coefficient gcd;
 
   // Main loop of the simplex algorithm
   for(;;) {
@@ -539,7 +550,7 @@ PIP_Solution_Node::solve(PIP_Tree_Node** parent_ref,
       std::cout << "Found row with negative parameters: " << i_
                 << std::endl;
 #endif
-      dimension_type j;
+      dimension_type j, k;
       Coefficient z(0);
       Coefficient sij, cij, cij_;
       Coefficient c;
@@ -590,11 +601,107 @@ PIP_Solution_Node::solve(PIP_Tree_Node** parent_ref,
                 << std::endl;
 #endif
 
+      /* ** Perform pivot operation ** */
+
+      /* Check if column j_ or row i_ correspond to a problem variable */
+      dimension_type var_j = n_a_d;
+      dimension_type 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_ */
+      Row prs(num_vars, tableau.s.capacity(), Row::Flags());
+      Row prt(num_params, tableau.t.capacity(), Row::Flags());
+      prs[j_] = tableau.s.get_denominator();
+      /* swap it with pivot row which would become identity after pivoting */
+      prs.swap(tableau.s[i_]);
+      prt.swap(tableau.t[i_]);
+      sign[i_] = ZERO;
+      for (j=0; j<num_vars; ++j) {
+        if (j==j_)
+          continue;
+        for (k=0; k<num_rows; ++k) {
+          Coefficient mult = prs[j] * tableau.s[k][j_];
+          if (mult % sij != 0) {
+            // Must scale matrix to stay in integer case
+            gcd_assign(gcd, mult, sij);
+            Coefficient scale_factor = sij/gcd;
+            tableau.s.scale(scale_factor);
+            mult *= scale_factor;
+          }
+          tableau.s[k][j] -= mult / sij;
+        }
+      }
+
+      /* create the identity matrix row corresponding to basic variable j_ */
+      for (j=0; j<num_params; ++j) {
+        for (k=0; k<num_rows; ++k) {
+          Coefficient c = prt[j] * tableau.s[k][j_];
+          if (c % sij != 0) {
+            // Must scale matrix to stay in integer case
+            gcd_assign(gcd, c, sij);
+            Coefficient scale_factor = sij/gcd;
+            tableau.t.scale(scale_factor);
+            c *= scale_factor;
+          }
+          c /= sij;
+          tableau.t[k][j] -= c;
+
+          s = sign[k];
+          if (s != MIXED) {
+            switch (s) {
+               case ZERO:
+                if (c > 0)
+                  sign[k] = NEGATIVE;
+                else if (c < 0)
+                  sign[k] = POSITIVE;
+                break;
+              case POSITIVE:
+                if (c > 0)
+                  sign[k] = MIXED;
+                break;
+              case NEGATIVE:
+                if (c < 0)
+                  sign[k] = MIXED;
+                break;
+              default:
+                break;
+            }
+          }
+        }
+      }
+
+      for (k=0; k<num_rows; ++k) {
+        Coefficient &c = tableau.s[k][j_];
+        if (c % sij != 0) {
+          Coefficient gcd;
+          gcd_assign(gcd, c, sij);
+          Coefficient scale_factor = sij/gcd;
+          tableau.s.scale(scale_factor);
+        }
+        c /= sij;
+      }
     }
+    //FIXME: to be finished
 
   } // Main loop of the simplex algorithm
 
-  //FIXME
   return OPTIMIZED_PIP_PROBLEM;
 }
 
diff --git a/src/PIP_Tree.defs.hh b/src/PIP_Tree.defs.hh
index c8b2a16..e88aa0f 100644
--- a/src/PIP_Tree.defs.hh
+++ b/src/PIP_Tree.defs.hh
@@ -180,6 +180,9 @@ private:
     //! Copy constructor.
     Rational_Matrix(const Rational_Matrix& y);
 
+    //! Multiplies all coefficients and denominator with ratio.
+    void scale(const Coefficient &ratio);
+
     //! Normalizes the modulo of coefficients so that they are mutually prime.
     /*!
       Computes the Greatest Common Divisor (GCD) among the elements of