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