Ticket #415: network_simplex.patch

.h

@@ -243 +243 @@
 
   private:
 
+    // Handle multiplication of costs with _state and _pred_dir
+    // elements which are in {-1, 0, 1} and {-1, 1} respectively. This
+    // means the Cost type '*' operator is never used and need not
+    // exist.
+    static inline Cost unitMul(signed char state, Cost cost) {
+      // fastest for integral costs on x86:
+      // return cost * (int)state;
+
+      // fastest for integral costs on machines with predicated execution
+      // and/or high latency integer multiplies (e.g. arm, ia64):
+      if(state == 1)
+        return cost;
+      else if(state == -1)
+        return -cost;
+      else
+        return 0;
+    }
+
     // Implementation of the First Eligible pivot rule
     class FirstEligiblePivotRule
     {
@@ -274 +292 @@
       bool findEnteringArc() {
         Cost c;
         for (int e = _next_arc; e != _search_arc_num; ++e) {
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (c < 0) {
             _in_arc = e;
             _next_arc = e + 1;
@@ -282 +300 @@
           }
         }
         for (int e = 0; e != _next_arc; ++e) {
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (c < 0) {
             _in_arc = e;
             _next_arc = e + 1;
@@ -322 +340 @@
       bool findEnteringArc() {
         Cost c, min = 0;
         for (int e = 0; e != _search_arc_num; ++e) {
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (c < min) {
             min = c;
             _in_arc = e;
@@ -376 +394 @@
         int cnt = _block_size;
         int e;
         for (e = _next_arc; e != _search_arc_num; ++e) {
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (c < min) {
             min = c;
             _in_arc = e;
@@ -387 +405 @@
           }
         }
         for (e = 0; e != _next_arc; ++e) {
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (c < min) {
             min = c;
             _in_arc = e;
@@ -462 +480 @@
           min = 0;
           for (int i = 0; i < _curr_length; ++i) {
             e = _candidates[i];
-            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+            c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
             if (c < min) {
               min = c;
               _in_arc = e;
@@ -478 +496 @@
         min = 0;
         _curr_length = 0;
         for (e = _next_arc; e != _search_arc_num; ++e) {
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (c < 0) {
             _candidates[_curr_length++] = e;
             if (c < min) {
@@ -489 +507 @@
           }
         }
         for (e = 0; e != _next_arc; ++e) {
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (c < 0) {
             _candidates[_curr_length++] = e;
             if (c < min) {
@@ -575 +593 @@
         Cost c;
         for (int i = 0; i != _curr_length; ++i) {
           e = _candidates[i];
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (c < 0) {
             _cand_cost[e] = c;
           } else {
@@ -588 +606 @@
         int limit = _head_length;
 
         for (e = _next_arc; e != _search_arc_num; ++e) {
-          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          c = unitMul(_state[e], (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (c < 0) {
             _cand_cost[e] = c;
             _candidates[_curr_length++] = e;
@@ -600 +618 @@
           }
         }
         for (e = 0; e != _next_arc; ++e) {
-          _cand_cost[e] = _state[e] *
-            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          _cand_cost[e] = unitMul(_state[e], 
+            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]));
           if (_cand_cost[e] < 0) {
             _candidates[_curr_length++] = e;
           }
@@ -986 +1004 @@
       Number c = 0;
       for (ArcIt a(_graph); a != INVALID; ++a) {
         int i = _arc_id[a];
-        c += Number(_flow[i]) * Number(_cost[i]);
+        c += _cost[i] * _flow[i];
       }
       return c;
     }
@@ -1449 +1467 @@
     // Update potentials in the subtree that has been moved
     void updatePotential() {
       Cost sigma = _pi[v_in] - _pi[u_in] -
-                   _pred_dir[u_in] * _cost[in_arc];
+                   unitMul(_pred_dir[u_in], _cost[in_arc]);
       int end = _thread[_last_succ[u_in]];
       for (int u = u_in; u != end; u = _thread[u]) {
         _pi[u] += sigma;
@@ -1536 +1554 @@
       // Perform heuristic initial pivots
       for (int i = 0; i != int(arc_vector.size()); ++i) {
         in_arc = arc_vector[i];
-        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
-            _pi[_target[in_arc]]) >= 0) continue;
+        if (unitMul(_state[in_arc], (_cost[in_arc] + _pi[_source[in_arc]] -
+            _pi[_target[in_arc]])) >= 0) continue;
         findJoinNode();
         bool change = findLeavingArc();
         if (delta >= MAX) return false;