1.1 --- a/lemon/capacity_scaling.h Thu Nov 12 23:34:35 2009 +0100
1.2 +++ b/lemon/capacity_scaling.h Thu Nov 12 23:45:15 2009 +0100
1.3 @@ -133,7 +133,7 @@
1.4 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.5
1.6 typedef std::vector<int> IntVector;
1.7 - typedef std::vector<bool> BoolVector;
1.8 + typedef std::vector<char> BoolVector;
1.9 typedef std::vector<Value> ValueVector;
1.10 typedef std::vector<Cost> CostVector;
1.11
1.12 @@ -196,6 +196,7 @@
1.13 private:
1.14
1.15 int _node_num;
1.16 + bool _geq;
1.17 const IntVector &_first_out;
1.18 const IntVector &_target;
1.19 const CostVector &_cost;
1.20 @@ -210,10 +211,10 @@
1.21 public:
1.22
1.23 ResidualDijkstra(CapacityScaling& cs) :
1.24 - _node_num(cs._node_num), _first_out(cs._first_out),
1.25 - _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap),
1.26 - _excess(cs._excess), _pi(cs._pi), _pred(cs._pred),
1.27 - _dist(cs._node_num)
1.28 + _node_num(cs._node_num), _geq(cs._sum_supply < 0),
1.29 + _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
1.30 + _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
1.31 + _pred(cs._pred), _dist(cs._node_num)
1.32 {}
1.33
1.34 int run(int s, Value delta = 1) {
1.35 @@ -232,7 +233,8 @@
1.36 heap.pop();
1.37
1.38 // Traverse outgoing residual arcs
1.39 - for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
1.40 + int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
1.41 + for (int a = _first_out[u]; a != last_out; ++a) {
1.42 if (_res_cap[a] < delta) continue;
1.43 v = _target[a];
1.44 switch (heap.state(v)) {
1.45 @@ -687,22 +689,25 @@
1.46 }
1.47 if (_sum_supply > 0) return INFEASIBLE;
1.48
1.49 - // Initialize maps
1.50 + // Initialize vectors
1.51 for (int i = 0; i != _root; ++i) {
1.52 _pi[i] = 0;
1.53 _excess[i] = _supply[i];
1.54 }
1.55
1.56 // Remove non-zero lower bounds
1.57 + const Value MAX = std::numeric_limits<Value>::max();
1.58 + int last_out;
1.59 if (_have_lower) {
1.60 for (int i = 0; i != _root; ++i) {
1.61 - for (int j = _first_out[i]; j != _first_out[i+1]; ++j) {
1.62 + last_out = _first_out[i+1];
1.63 + for (int j = _first_out[i]; j != last_out; ++j) {
1.64 if (_forward[j]) {
1.65 Value c = _lower[j];
1.66 if (c >= 0) {
1.67 - _res_cap[j] = _upper[j] < INF ? _upper[j] - c : INF;
1.68 + _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
1.69 } else {
1.70 - _res_cap[j] = _upper[j] < INF + c ? _upper[j] - c : INF;
1.71 + _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
1.72 }
1.73 _excess[i] -= c;
1.74 _excess[_target[j]] += c;
1.75 @@ -718,15 +723,16 @@
1.76 }
1.77
1.78 // Handle negative costs
1.79 - for (int u = 0; u != _root; ++u) {
1.80 - for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
1.81 - Value rc = _res_cap[a];
1.82 - if (_cost[a] < 0 && rc > 0) {
1.83 - if (rc == INF) return UNBOUNDED;
1.84 - _excess[u] -= rc;
1.85 - _excess[_target[a]] += rc;
1.86 - _res_cap[a] = 0;
1.87 - _res_cap[_reverse[a]] += rc;
1.88 + for (int i = 0; i != _root; ++i) {
1.89 + last_out = _first_out[i+1] - 1;
1.90 + for (int j = _first_out[i]; j != last_out; ++j) {
1.91 + Value rc = _res_cap[j];
1.92 + if (_cost[j] < 0 && rc > 0) {
1.93 + if (rc >= MAX) return UNBOUNDED;
1.94 + _excess[i] -= rc;
1.95 + _excess[_target[j]] += rc;
1.96 + _res_cap[j] = 0;
1.97 + _res_cap[_reverse[j]] += rc;
1.98 }
1.99 }
1.100 }
1.101 @@ -736,24 +742,21 @@
1.102 _pi[_root] = 0;
1.103 _excess[_root] = -_sum_supply;
1.104 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.105 - int u = _target[a];
1.106 - if (_excess[u] < 0) {
1.107 - _res_cap[a] = -_excess[u] + 1;
1.108 - } else {
1.109 - _res_cap[a] = 1;
1.110 - }
1.111 - _res_cap[_reverse[a]] = 0;
1.112 + int ra = _reverse[a];
1.113 + _res_cap[a] = -_sum_supply + 1;
1.114 + _res_cap[ra] = 0;
1.115 _cost[a] = 0;
1.116 - _cost[_reverse[a]] = 0;
1.117 + _cost[ra] = 0;
1.118 }
1.119 } else {
1.120 _pi[_root] = 0;
1.121 _excess[_root] = 0;
1.122 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.123 + int ra = _reverse[a];
1.124 _res_cap[a] = 1;
1.125 - _res_cap[_reverse[a]] = 0;
1.126 + _res_cap[ra] = 0;
1.127 _cost[a] = 0;
1.128 - _cost[_reverse[a]] = 0;
1.129 + _cost[ra] = 0;
1.130 }
1.131 }
1.132
1.133 @@ -762,8 +765,9 @@
1.134 // With scaling
1.135 Value max_sup = 0, max_dem = 0;
1.136 for (int i = 0; i != _node_num; ++i) {
1.137 - if ( _excess[i] > max_sup) max_sup = _excess[i];
1.138 - if (-_excess[i] > max_dem) max_dem = -_excess[i];
1.139 + Value ex = _excess[i];
1.140 + if ( ex > max_sup) max_sup = ex;
1.141 + if (-ex > max_dem) max_dem = -ex;
1.142 }
1.143 Value max_cap = 0;
1.144 for (int j = 0; j != _res_arc_num; ++j) {
1.145 @@ -789,7 +793,8 @@
1.146
1.147 // Handle non-zero lower bounds
1.148 if (_have_lower) {
1.149 - for (int j = 0; j != _res_arc_num - _node_num + 1; ++j) {
1.150 + int limit = _first_out[_root];
1.151 + for (int j = 0; j != limit; ++j) {
1.152 if (!_forward[j]) _res_cap[j] += _lower[j];
1.153 }
1.154 }
1.155 @@ -812,8 +817,11 @@
1.156 ResidualDijkstra _dijkstra(*this);
1.157 while (true) {
1.158 // Saturate all arcs not satisfying the optimality condition
1.159 + int last_out;
1.160 for (int u = 0; u != _node_num; ++u) {
1.161 - for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
1.162 + last_out = _sum_supply < 0 ?
1.163 + _first_out[u+1] : _first_out[u+1] - 1;
1.164 + for (int a = _first_out[u]; a != last_out; ++a) {
1.165 int v = _target[a];
1.166 Cost c = _cost[a] + _pi[u] - _pi[v];
1.167 Value rc = _res_cap[a];
1.168 @@ -830,8 +838,9 @@
1.169 _excess_nodes.clear();
1.170 _deficit_nodes.clear();
1.171 for (int u = 0; u != _node_num; ++u) {
1.172 - if (_excess[u] >= _delta) _excess_nodes.push_back(u);
1.173 - if (_excess[u] <= -_delta) _deficit_nodes.push_back(u);
1.174 + Value ex = _excess[u];
1.175 + if (ex >= _delta) _excess_nodes.push_back(u);
1.176 + if (ex <= -_delta) _deficit_nodes.push_back(u);
1.177 }
1.178 int next_node = 0, next_def_node = 0;
1.179
2.1 --- a/lemon/network_simplex.h Thu Nov 12 23:34:35 2009 +0100
2.2 +++ b/lemon/network_simplex.h Thu Nov 12 23:45:15 2009 +0100
2.3 @@ -164,7 +164,7 @@
2.4 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
2.5
2.6 typedef std::vector<int> IntVector;
2.7 - typedef std::vector<bool> BoolVector;
2.8 + typedef std::vector<char> CharVector;
2.9 typedef std::vector<Value> ValueVector;
2.10 typedef std::vector<Cost> CostVector;
2.11
2.12 @@ -212,8 +212,8 @@
2.13 IntVector _succ_num;
2.14 IntVector _last_succ;
2.15 IntVector _dirty_revs;
2.16 - BoolVector _forward;
2.17 - IntVector _state;
2.18 + CharVector _forward;
2.19 + CharVector _state;
2.20 int _root;
2.21
2.22 // Temporary data used in the current pivot iteration
2.23 @@ -221,6 +221,8 @@
2.24 int first, second, right, last;
2.25 int stem, par_stem, new_stem;
2.26 Value delta;
2.27 +
2.28 + const Value MAX;
2.29
2.30 public:
2.31
2.32 @@ -242,7 +244,7 @@
2.33 const IntVector &_source;
2.34 const IntVector &_target;
2.35 const CostVector &_cost;
2.36 - const IntVector &_state;
2.37 + const CharVector &_state;
2.38 const CostVector &_pi;
2.39 int &_in_arc;
2.40 int _search_arc_num;
2.41 @@ -294,7 +296,7 @@
2.42 const IntVector &_source;
2.43 const IntVector &_target;
2.44 const CostVector &_cost;
2.45 - const IntVector &_state;
2.46 + const CharVector &_state;
2.47 const CostVector &_pi;
2.48 int &_in_arc;
2.49 int _search_arc_num;
2.50 @@ -333,7 +335,7 @@
2.51 const IntVector &_source;
2.52 const IntVector &_target;
2.53 const CostVector &_cost;
2.54 - const IntVector &_state;
2.55 + const CharVector &_state;
2.56 const CostVector &_pi;
2.57 int &_in_arc;
2.58 int _search_arc_num;
2.59 @@ -406,7 +408,7 @@
2.60 const IntVector &_source;
2.61 const IntVector &_target;
2.62 const CostVector &_cost;
2.63 - const IntVector &_state;
2.64 + const CharVector &_state;
2.65 const CostVector &_pi;
2.66 int &_in_arc;
2.67 int _search_arc_num;
2.68 @@ -509,7 +511,7 @@
2.69 const IntVector &_source;
2.70 const IntVector &_target;
2.71 const CostVector &_cost;
2.72 - const IntVector &_state;
2.73 + const CharVector &_state;
2.74 const CostVector &_pi;
2.75 int &_in_arc;
2.76 int _search_arc_num;
2.77 @@ -631,9 +633,9 @@
2.78 /// but it is usually slower. Therefore it is disabled by default.
2.79 NetworkSimplex(const GR& graph, bool arc_mixing = false) :
2.80 _graph(graph), _node_id(graph), _arc_id(graph),
2.81 + MAX(std::numeric_limits<Value>::max()),
2.82 INF(std::numeric_limits<Value>::has_infinity ?
2.83 - std::numeric_limits<Value>::infinity() :
2.84 - std::numeric_limits<Value>::max())
2.85 + std::numeric_limits<Value>::infinity() : MAX)
2.86 {
2.87 // Check the value types
2.88 LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
2.89 @@ -1020,9 +1022,9 @@
2.90 for (int i = 0; i != _arc_num; ++i) {
2.91 Value c = _lower[i];
2.92 if (c >= 0) {
2.93 - _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
2.94 + _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
2.95 } else {
2.96 - _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
2.97 + _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
2.98 }
2.99 _supply[_source[i]] -= c;
2.100 _supply[_target[i]] += c;
2.101 @@ -1214,7 +1216,7 @@
2.102 for (int u = first; u != join; u = _parent[u]) {
2.103 e = _pred[u];
2.104 d = _forward[u] ?
2.105 - _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
2.106 + _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
2.107 if (d < delta) {
2.108 delta = d;
2.109 u_out = u;
2.110 @@ -1225,7 +1227,7 @@
2.111 for (int u = second; u != join; u = _parent[u]) {
2.112 e = _pred[u];
2.113 d = _forward[u] ?
2.114 - (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
2.115 + (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
2.116 if (d <= delta) {
2.117 delta = d;
2.118 u_out = u;
2.119 @@ -1424,7 +1426,7 @@
2.120 while (pivot.findEnteringArc()) {
2.121 findJoinNode();
2.122 bool change = findLeavingArc();
2.123 - if (delta >= INF) return UNBOUNDED;
2.124 + if (delta >= MAX) return UNBOUNDED;
2.125 changeFlow(change);
2.126 if (change) {
2.127 updateTreeStructure();