1.1 --- a/lemon/network_simplex.h Fri Feb 06 21:52:34 2009 +0000
1.2 +++ b/lemon/network_simplex.h Mon Jun 01 15:35:27 2009 +0000
1.3 @@ -241,7 +241,7 @@
1.4 BlockSearchPivotRule(NetworkSimplex &ns) :
1.5 _edge(ns._edge), _source(ns._source), _target(ns._target),
1.6 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.7 - _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
1.8 + _in_edge(ns._in_edge), _edge_num(ns._edge_num + ns._node_num), _next_edge(0)
1.9 {
1.10 // The main parameters of the pivot rule
1.11 const double BLOCK_SIZE_FACTOR = 2.0;
1.12 @@ -591,14 +591,19 @@
1.13 // The number of edges in the original graph
1.14 int _edge_num;
1.15
1.16 - // The depth vector of the spanning tree structure
1.17 - IntVector _depth;
1.18 // The parent vector of the spanning tree structure
1.19 IntVector _parent;
1.20 // The pred_edge vector of the spanning tree structure
1.21 IntVector _pred;
1.22 // The thread vector of the spanning tree structure
1.23 IntVector _thread;
1.24 +
1.25 + IntVector _rev_thread;
1.26 + IntVector _succ_num;
1.27 + IntVector _last_succ;
1.28 +
1.29 + IntVector _dirty_revs;
1.30 +
1.31 // The forward vector of the spanning tree structure
1.32 BoolVector _forward;
1.33 // The state vector
1.34 @@ -882,11 +887,13 @@
1.35 _flow.resize(all_edge_num, 0);
1.36 _pi.resize(all_node_num, 0);
1.37
1.38 - _depth.resize(all_node_num);
1.39 _parent.resize(all_node_num);
1.40 _pred.resize(all_node_num);
1.41 + _forward.resize(all_node_num);
1.42 _thread.resize(all_node_num);
1.43 - _forward.resize(all_node_num);
1.44 + _rev_thread.resize(all_node_num);
1.45 + _succ_num.resize(all_node_num);
1.46 + _last_succ.resize(all_node_num);
1.47 _state.resize(all_edge_num, STATE_LOWER);
1.48
1.49 // Initialize node related data
1.50 @@ -915,10 +922,12 @@
1.51
1.52 // Set data for the artificial root node
1.53 _root = _node_num;
1.54 - _depth[_root] = 0;
1.55 _parent[_root] = -1;
1.56 _pred[_root] = -1;
1.57 _thread[_root] = 0;
1.58 + _rev_thread[0] = _root;
1.59 + _succ_num[_root] = all_node_num;
1.60 + _last_succ[_root] = _root - 1;
1.61 _supply[_root] = 0;
1.62 _pi[_root] = 0;
1.63
1.64 @@ -955,22 +964,30 @@
1.65 Cost max_cost = std::numeric_limits<Cost>::max() / 4;
1.66 Capacity max_cap = std::numeric_limits<Capacity>::max();
1.67 for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
1.68 - _thread[u] = u + 1;
1.69 - _depth[u] = 1;
1.70 _parent[u] = _root;
1.71 _pred[u] = e;
1.72 + _thread[u] = u + 1;
1.73 + _rev_thread[u + 1] = u;
1.74 + _succ_num[u] = 1;
1.75 + _last_succ[u] = u;
1.76 + _cap[e] = max_cap;
1.77 + _state[e] = STATE_TREE;
1.78 if (_supply[u] >= 0) {
1.79 + _forward[u] = true;
1.80 + _pi[u] = 0;
1.81 + _source[e] = u;
1.82 + _target[e] = _root;
1.83 _flow[e] = _supply[u];
1.84 - _forward[u] = true;
1.85 - _pi[u] = -max_cost;
1.86 - } else {
1.87 - _flow[e] = -_supply[u];
1.88 + _cost[e] = 0;
1.89 + }
1.90 + else {
1.91 _forward[u] = false;
1.92 _pi[u] = max_cost;
1.93 + _source[e] = _root;
1.94 + _target[e] = u;
1.95 + _flow[e] = -_supply[u];
1.96 + _cost[e] = max_cost;
1.97 }
1.98 - _cost[e] = max_cost;
1.99 - _cap[e] = max_cap;
1.100 - _state[e] = STATE_TREE;
1.101 }
1.102
1.103 return true;
1.104 @@ -980,11 +997,12 @@
1.105 void findJoinNode() {
1.106 int u = _source[_in_edge];
1.107 int v = _target[_in_edge];
1.108 - while (_depth[u] > _depth[v]) u = _parent[u];
1.109 - while (_depth[v] > _depth[u]) v = _parent[v];
1.110 while (u != v) {
1.111 - u = _parent[u];
1.112 - v = _parent[v];
1.113 + if (_succ_num[u] < _succ_num[v]) {
1.114 + u = _parent[u];
1.115 + } else {
1.116 + v = _parent[v];
1.117 + }
1.118 }
1.119 join = u;
1.120 }
1.121 @@ -1059,88 +1077,144 @@
1.122 _state[_in_edge] = -_state[_in_edge];
1.123 }
1.124 }
1.125 -
1.126 - // Update _thread and _parent vectors
1.127 - void updateThreadParent() {
1.128 - int u;
1.129 +
1.130 + // Update the tree structure
1.131 + void updateTreeStructure() {
1.132 + int u, w;
1.133 + int old_rev_thread = _rev_thread[u_out];
1.134 + int old_succ_num = _succ_num[u_out];
1.135 + int old_last_succ = _last_succ[u_out];
1.136 v_out = _parent[u_out];
1.137
1.138 - // Handle the case when join and v_out coincide
1.139 - bool par_first = false;
1.140 - if (join == v_out) {
1.141 - for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
1.142 - if (u == v_in) {
1.143 - par_first = true;
1.144 - while (_thread[u] != u_out) u = _thread[u];
1.145 - first = u;
1.146 + u = _last_succ[u_in]; // the last successor of u_in
1.147 + right = _thread[u]; // the node after it
1.148 +
1.149 + // Handle the case when old_rev_thread equals to v_in
1.150 + // (it also means that join and v_out coincide)
1.151 + if (old_rev_thread == v_in) {
1.152 + last = _thread[_last_succ[u_out]];
1.153 + } else {
1.154 + last = _thread[v_in];
1.155 + }
1.156 +
1.157 + // Update _thread and _parent along the stem nodes (i.e. the nodes
1.158 + // between u_in and u_out, whose parent have to be changed)
1.159 + _thread[v_in] = stem = u_in;
1.160 + _dirty_revs.clear();
1.161 + _dirty_revs.push_back(v_in);
1.162 + par_stem = v_in;
1.163 + while (stem != u_out) {
1.164 + // Insert the next stem node into the thread list
1.165 + new_stem = _parent[stem];
1.166 + _thread[u] = new_stem;
1.167 + _dirty_revs.push_back(u);
1.168 +
1.169 + // Remove the subtree of stem from the thread list
1.170 + w = _rev_thread[stem];
1.171 + _thread[w] = right;
1.172 + _rev_thread[right] = w;
1.173 +
1.174 + // Change the parent node and shift stem nodes
1.175 + _parent[stem] = par_stem;
1.176 + par_stem = stem;
1.177 + stem = new_stem;
1.178 +
1.179 + // Update u and right nodes
1.180 + u = _last_succ[stem] == _last_succ[par_stem] ?
1.181 + _rev_thread[par_stem] : _last_succ[stem];
1.182 + right = _thread[u];
1.183 + }
1.184 + _parent[u_out] = par_stem;
1.185 + _last_succ[u_out] = u;
1.186 + _thread[u] = last;
1.187 + _rev_thread[last] = u;
1.188 +
1.189 + // Remove the subtree of u_out from the thread list except for
1.190 + // the case when old_rev_thread equals to v_in
1.191 + // (it also means that join and v_out coincide)
1.192 + if (old_rev_thread != v_in) {
1.193 + _thread[old_rev_thread] = right;
1.194 + _rev_thread[right] = old_rev_thread;
1.195 + }
1.196 +
1.197 + // Update _rev_thread using the new _thread values
1.198 + for (int i = 0; i < int(_dirty_revs.size()); ++i) {
1.199 + u = _dirty_revs[i];
1.200 + _rev_thread[_thread[u]] = u;
1.201 + }
1.202 +
1.203 + // Update _last_succ for the stem nodes from u_out to u_in
1.204 + for (u = u_out; u != u_in; u = _parent[u]) {
1.205 + _last_succ[_parent[u]] = _last_succ[u];
1.206 + }
1.207 +
1.208 + // Set limits for updating _last_succ form v_in and v_out
1.209 + // towards the root
1.210 + int up_limit_in = -1;
1.211 + int up_limit_out = -1;
1.212 + if (_last_succ[join] == v_in) {
1.213 + up_limit_out = join;
1.214 + } else {
1.215 + up_limit_in = join;
1.216 + }
1.217 +
1.218 + // Update _last_succ from v_in towards the root
1.219 + for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
1.220 + u = _parent[u]) {
1.221 + _last_succ[u] = _last_succ[u_out];
1.222 + }
1.223 + // Update _last_succ from v_out towards the root
1.224 + if (join != old_rev_thread && v_in != old_rev_thread) {
1.225 + for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.226 + u = _parent[u]) {
1.227 + _last_succ[u] = old_rev_thread;
1.228 + }
1.229 + } else {
1.230 + for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
1.231 + u = _parent[u]) {
1.232 + _last_succ[u] = _last_succ[u_out];
1.233 }
1.234 }
1.235
1.236 - // Find the last successor of u_in (u) and the node after it (right)
1.237 - // according to the thread index
1.238 - for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
1.239 - right = _thread[u];
1.240 - if (_thread[v_in] == u_out) {
1.241 - for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
1.242 - if (last == u_out) last = _thread[last];
1.243 - }
1.244 - else last = _thread[v_in];
1.245 -
1.246 - // Update stem nodes
1.247 - _thread[v_in] = stem = u_in;
1.248 - par_stem = v_in;
1.249 - while (stem != u_out) {
1.250 - _thread[u] = new_stem = _parent[stem];
1.251 -
1.252 - // Find the node just before the stem node (u) according to
1.253 - // the original thread index
1.254 - for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
1.255 - _thread[u] = right;
1.256 -
1.257 - // Change the parent node of stem and shift stem and par_stem nodes
1.258 - _parent[stem] = par_stem;
1.259 - par_stem = stem;
1.260 - stem = new_stem;
1.261 -
1.262 - // Find the last successor of stem (u) and the node after it (right)
1.263 - // according to the thread index
1.264 - for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
1.265 - right = _thread[u];
1.266 - }
1.267 - _parent[u_out] = par_stem;
1.268 - _thread[u] = last;
1.269 -
1.270 - if (join == v_out && par_first) {
1.271 - if (first != v_in) _thread[first] = right;
1.272 - } else {
1.273 - for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
1.274 - _thread[u] = right;
1.275 - }
1.276 - }
1.277 -
1.278 - // Update _pred and _forward vectors
1.279 - void updatePredEdge() {
1.280 - int u = u_out, v;
1.281 + // Update _pred and _forward for the stem nodes from u_out to u_in
1.282 + u = u_out;
1.283 while (u != u_in) {
1.284 - v = _parent[u];
1.285 - _pred[u] = _pred[v];
1.286 - _forward[u] = !_forward[v];
1.287 - u = v;
1.288 + w = _parent[u];
1.289 + _pred[u] = _pred[w];
1.290 + _forward[u] = !_forward[w];
1.291 + u = w;
1.292 }
1.293 _pred[u_in] = _in_edge;
1.294 _forward[u_in] = (u_in == _source[_in_edge]);
1.295 +
1.296 + // Update _succ_num from v_in to join
1.297 + for (u = v_in; u != join; u = _parent[u]) {
1.298 + _succ_num[u] += old_succ_num;
1.299 + }
1.300 + // Update _succ_num from v_out to join
1.301 + for (u = v_out; u != join; u = _parent[u]) {
1.302 + _succ_num[u] -= old_succ_num;
1.303 + }
1.304 + // Update _succ_num for the stem nodes from u_out to u_in
1.305 + int tmp = 0;
1.306 + u = u_out;
1.307 + while (u != u_in) {
1.308 + w = _parent[u];
1.309 + tmp = _succ_num[u] - _succ_num[w] + tmp;
1.310 + _succ_num[u] = tmp;
1.311 + u = w;
1.312 + }
1.313 + _succ_num[u_in] = old_succ_num;
1.314 }
1.315
1.316 - // Update _depth and _potential vectors
1.317 - void updateDepthPotential() {
1.318 - _depth[u_in] = _depth[v_in] + 1;
1.319 + // Update potentials
1.320 + void updatePotential() {
1.321 Cost sigma = _forward[u_in] ?
1.322 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1.323 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1.324 - _pi[u_in] += sigma;
1.325 - for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
1.326 - _depth[u] = _depth[_parent[u]] + 1;
1.327 - if (_depth[u] <= _depth[u_in]) break;
1.328 + // Update in the lower subtree (which has been moved)
1.329 + int end = _thread[_last_succ[u_in]];
1.330 + for (int u = u_in; u != end; u = _thread[u]) {
1.331 _pi[u] += sigma;
1.332 }
1.333 }
1.334 @@ -1166,19 +1240,18 @@
1.335 template<class PivotRuleImplementation>
1.336 bool start() {
1.337 PivotRuleImplementation pivot(*this);
1.338 -
1.339 +
1.340 // Execute the network simplex algorithm
1.341 while (pivot.findEnteringEdge()) {
1.342 findJoinNode();
1.343 bool change = findLeavingEdge();
1.344 changeFlow(change);
1.345 if (change) {
1.346 - updateThreadParent();
1.347 - updatePredEdge();
1.348 - updateDepthPotential();
1.349 + updateTreeStructure();
1.350 + updatePotential();
1.351 }
1.352 }
1.353 -
1.354 +
1.355 // Check if the flow amount equals zero on all the artificial edges
1.356 for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
1.357 if (_flow[e] > 0) return false;
1.358 @@ -1199,7 +1272,7 @@
1.359 for (int i = 0; i != _node_num; ++i) {
1.360 (*_potential_result)[_node[i]] = _pi[i];
1.361 }
1.362 -
1.363 +
1.364 return true;
1.365 }
1.366