lemon-0.x: comparison lemon/network

equal deleted inserted replaced

-:13ab2cd16858
+:bd106abd230b
 /// Constructor
 BlockSearchPivotRule(NetworkSimplex &ns) :
 _edge(ns._edge), _source(ns._source), _target(ns._target),
 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
-_in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
+_in_edge(ns._in_edge), _edge_num(ns._edge_num + ns._node_num), _next_edge(0)
 {
 // The main parameters of the pivot rule
 const double BLOCK_SIZE_FACTOR = 2.0;
 const int MIN_BLOCK_SIZE = 10;
 // The number of nodes in the original graph
 int _node_num;
 // The number of edges in the original graph
 int _edge_num;
-// The depth vector of the spanning tree structure
-IntVector _depth;
 // The parent vector of the spanning tree structure
 IntVector _parent;
 // The pred_edge vector of the spanning tree structure
 IntVector _pred;
 // The thread vector of the spanning tree structure
 IntVector _thread;
+IntVector _rev_thread;
+IntVector _succ_num;
+IntVector _last_succ;
+IntVector _dirty_revs;
 // The forward vector of the spanning tree structure
 BoolVector _forward;
 // The state vector
 IntVector _state;
 // The root node
 _cost.resize(all_edge_num);
 _supply.resize(all_node_num);
 _flow.resize(all_edge_num, 0);
 _pi.resize(all_node_num, 0);
-_depth.resize(all_node_num);
 _parent.resize(all_node_num);
 _pred.resize(all_node_num);
+_forward.resize(all_node_num);
 _thread.resize(all_node_num);
-_forward.resize(all_node_num);
+_rev_thread.resize(all_node_num);
+_succ_num.resize(all_node_num);
+_last_succ.resize(all_node_num);
 _state.resize(all_edge_num, STATE_LOWER);
 // Initialize node related data
 bool valid_supply = true;
 if (_orig_supply) {
 }
 if (!valid_supply) return false;
 // Set data for the artificial root node
 _root = _node_num;
-_depth[_root] = 0;
 _parent[_root] = -1;
 _pred[_root] = -1;
 _thread[_root] = 0;
+_rev_thread[0] = _root;
+_succ_num[_root] = all_node_num;
+_last_succ[_root] = _root - 1;
 _supply[_root] = 0;
 _pi[_root] = 0;
 // Store the edges in a mixed order
 int k = std::max(int(sqrt(_edge_num)), 10);
 // Add artificial edges and initialize the spanning tree data structure
 Cost max_cost = std::numeric_limits<Cost>::max() / 4;
 Capacity max_cap = std::numeric_limits<Capacity>::max();
 for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
-_thread[u] = u + 1;
-_depth[u] = 1;
 _parent[u] = _root;
 _pred[u] = e;
+_thread[u] = u + 1;
+_rev_thread[u + 1] = u;
+_succ_num[u] = 1;
+_last_succ[u] = u;
+_cap[e] = max_cap;
+_state[e] = STATE_TREE;
 if (_supply[u] >= 0) {
+_forward[u] = true;
+_pi[u] = 0;
+_source[e] = u;
+_target[e] = _root;
 _flow[e] = _supply[u];
-_forward[u] = true;
+_cost[e] = 0;
-_pi[u] = -max_cost;
+}
-} else {
+else {
-_flow[e] = -_supply[u];
 _forward[u] = false;
 _pi[u] = max_cost;
-}
+_source[e] = _root;
-_cost[e] = max_cost;
+_target[e] = u;
-_cap[e] = max_cap;
+_flow[e] = -_supply[u];
-_state[e] = STATE_TREE;
+_cost[e] = max_cost;
+}
 }
 return true;
 }
 // Find the join node
 void findJoinNode() {
 int u = _source[_in_edge];
 int v = _target[_in_edge];
-while (_depth[u] > _depth[v]) u = _parent[u];
-while (_depth[v] > _depth[u]) v = _parent[v];
 while (u != v) {
-u = _parent[u];
+if (_succ_num[u] < _succ_num[v]) {
-v = _parent[v];
+u = _parent[u];
+} else {
+v = _parent[v];
+}
 }
 join = u;
 }
 // Find the leaving edge of the cycle and returns true if the
 (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
 } else {
 _state[_in_edge] = -_state[_in_edge];
 }
 }
-// Update _thread and _parent vectors
+// Update the tree structure
-void updateThreadParent() {
+void updateTreeStructure() {
-int u;
+int u, w;
+int old_rev_thread = _rev_thread[u_out];
+int old_succ_num = _succ_num[u_out];
+int old_last_succ = _last_succ[u_out];
 v_out = _parent[u_out];
-// Handle the case when join and v_out coincide
+u = _last_succ[u_in];  // the last successor of u_in
-bool par_first = false;
+right = _thread[u];    // the node after it
-if (join == v_out) {
-for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
+// Handle the case when old_rev_thread equals to v_in
-if (u == v_in) {
+// (it also means that join and v_out coincide)
-par_first = true;
+if (old_rev_thread == v_in) {
-while (_thread[u] != u_out) u = _thread[u];
+last = _thread[_last_succ[u_out]];
-first = u;
+} else {
-}
+last = _thread[v_in];
 }
-// Find the last successor of u_in (u) and the node after it (right)
+// Update _thread and _parent along the stem nodes (i.e. the nodes
-// according to the thread index
+// between u_in and u_out, whose parent have to be changed)
-for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
-right = _thread[u];
-if (_thread[v_in] == u_out) {
-for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
-if (last == u_out) last = _thread[last];
-}
-else last = _thread[v_in];
-// Update stem nodes
 _thread[v_in] = stem = u_in;
+_dirty_revs.clear();
+_dirty_revs.push_back(v_in);
 par_stem = v_in;
 while (stem != u_out) {
-_thread[u] = new_stem = _parent[stem];
+// Insert the next stem node into the thread list
+new_stem = _parent[stem];
-// Find the node just before the stem node (u) according to
+_thread[u] = new_stem;
-// the original thread index
+_dirty_revs.push_back(u);
-for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
-_thread[u] = right;
+// Remove the subtree of stem from the thread list
+w = _rev_thread[stem];
-// Change the parent node of stem and shift stem and par_stem nodes
+_thread[w] = right;
-_parent[stem] = par_stem;
+_rev_thread[right] = w;
+// Change the parent node and shift stem nodes
+_parent[stem] = par_stem;
 par_stem = stem;
 stem = new_stem;
-// Find the last successor of stem (u) and the node after it (right)
+// Update u and right nodes
-// according to the thread index
+u = _last_succ[stem] == _last_succ[par_stem] ?
-for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
+_rev_thread[par_stem] : _last_succ[stem];
 right = _thread[u];
 }
 _parent[u_out] = par_stem;
+_last_succ[u_out] = u;
 _thread[u] = last;
+_rev_thread[last] = u;
-if (join == v_out && par_first) {
-if (first != v_in) _thread[first] = right;
+// Remove the subtree of u_out from the thread list except for
+// the case when old_rev_thread equals to v_in
+// (it also means that join and v_out coincide)
+if (old_rev_thread != v_in) {
+_thread[old_rev_thread] = right;
+_rev_thread[right] = old_rev_thread;
+}
+// Update _rev_thread using the new _thread values
+for (int i = 0; i < int(_dirty_revs.size()); ++i) {
+u = _dirty_revs[i];
+_rev_thread[_thread[u]] = u;
+}
+// Update _last_succ for the stem nodes from u_out to u_in
+for (u = u_out; u != u_in; u = _parent[u]) {
+_last_succ[_parent[u]] = _last_succ[u];
+}
+// Set limits for updating _last_succ form v_in and v_out
+// towards the root
+int up_limit_in = -1;
+int up_limit_out = -1;
+if (_last_succ[join] == v_in) {
+up_limit_out = join;
 } else {
-for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
+up_limit_in = join;
-_thread[u] = right;
+}
-}
-}
+// Update _last_succ from v_in towards the root
+for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
-// Update _pred and _forward vectors
+u = _parent[u]) {
-void updatePredEdge() {
+_last_succ[u] = _last_succ[u_out];
-int u = u_out, v;
+}
+// Update _last_succ from v_out towards the root
+if (join != old_rev_thread && v_in != old_rev_thread) {
+for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
+u = _parent[u]) {
+_last_succ[u] = old_rev_thread;
+}
+} else {
+for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
+u = _parent[u]) {
+_last_succ[u] = _last_succ[u_out];
+}
+}
+// Update _pred and _forward for the stem nodes from u_out to u_in
+u = u_out;
 while (u != u_in) {
-v = _parent[u];
+w = _parent[u];
-_pred[u] = _pred[v];
+_pred[u] = _pred[w];
-_forward[u] = !_forward[v];
+_forward[u] = !_forward[w];
-u = v;
+u = w;
 }
 _pred[u_in] = _in_edge;
 _forward[u_in] = (u_in == _source[_in_edge]);
-}
+// Update _succ_num from v_in to join
-// Update _depth and _potential vectors
+for (u = v_in; u != join; u = _parent[u]) {
-void updateDepthPotential() {
+_succ_num[u] += old_succ_num;
-_depth[u_in] = _depth[v_in] + 1;
+}
+// Update _succ_num from v_out to join
+for (u = v_out; u != join; u = _parent[u]) {
+_succ_num[u] -= old_succ_num;
+}
+// Update _succ_num for the stem nodes from u_out to u_in
+int tmp = 0;
+u = u_out;
+while (u != u_in) {
+w = _parent[u];
+tmp = _succ_num[u] - _succ_num[w] + tmp;
+_succ_num[u] = tmp;
+u = w;
+}
+_succ_num[u_in] = old_succ_num;
+}
+// Update potentials
+void updatePotential() {
 Cost sigma = _forward[u_in] ?
 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
-_pi[u_in] += sigma;
+// Update in the lower subtree (which has been moved)
-for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
+int end = _thread[_last_succ[u_in]];
-_depth[u] = _depth[_parent[u]] + 1;
+for (int u = u_in; u != end; u = _thread[u]) {
-if (_depth[u] <= _depth[u_in]) break;
 _pi[u] += sigma;
 }
 }
 // Execute the algorithm
 }
 template<class PivotRuleImplementation>
 bool start() {
 PivotRuleImplementation pivot(*this);
 // Execute the network simplex algorithm
 while (pivot.findEnteringEdge()) {
 findJoinNode();
 bool change = findLeavingEdge();
 changeFlow(change);
 if (change) {
-updateThreadParent();
+updateTreeStructure();
-updatePredEdge();
+updatePotential();
-updateDepthPotential();
+}
 }
-}
 // Check if the flow amount equals zero on all the artificial edges
 for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
 if (_flow[e] > 0) return false;
 }
 }
 // Copy potential values to _potential_result
 for (int i = 0; i != _node_num; ++i) {
 (*_potential_result)[_node[i]] = _pi[i];
 }
 return true;
 }
 }; //class NetworkSimplex

changeset 2637	bafe730864da
parent 2634	e98bbe64cca4
child 2638	61bf01404ede