1.1 --- a/lemon/network_simplex.h Wed Sep 30 08:36:43 2009 +0200
1.2 +++ b/lemon/network_simplex.h Wed Sep 30 08:41:06 2009 +0200
1.3 @@ -161,8 +161,6 @@
1.4
1.5 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.6
1.7 - typedef std::vector<Arc> ArcVector;
1.8 - typedef std::vector<Node> NodeVector;
1.9 typedef std::vector<int> IntVector;
1.10 typedef std::vector<bool> BoolVector;
1.11 typedef std::vector<Value> ValueVector;
1.12 @@ -364,33 +362,32 @@
1.13 bool findEnteringArc() {
1.14 Cost c, min = 0;
1.15 int cnt = _block_size;
1.16 - int e, min_arc = _next_arc;
1.17 + int e;
1.18 for (e = _next_arc; e < _search_arc_num; ++e) {
1.19 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.20 if (c < min) {
1.21 min = c;
1.22 - min_arc = e;
1.23 + _in_arc = e;
1.24 }
1.25 if (--cnt == 0) {
1.26 - if (min < 0) break;
1.27 + if (min < 0) goto search_end;
1.28 cnt = _block_size;
1.29 }
1.30 }
1.31 - if (min == 0 || cnt > 0) {
1.32 - for (e = 0; e < _next_arc; ++e) {
1.33 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.34 - if (c < min) {
1.35 - min = c;
1.36 - min_arc = e;
1.37 - }
1.38 - if (--cnt == 0) {
1.39 - if (min < 0) break;
1.40 - cnt = _block_size;
1.41 - }
1.42 + for (e = 0; e < _next_arc; ++e) {
1.43 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.44 + if (c < min) {
1.45 + min = c;
1.46 + _in_arc = e;
1.47 + }
1.48 + if (--cnt == 0) {
1.49 + if (min < 0) goto search_end;
1.50 + cnt = _block_size;
1.51 }
1.52 }
1.53 if (min >= 0) return false;
1.54 - _in_arc = min_arc;
1.55 +
1.56 + search_end:
1.57 _next_arc = e;
1.58 return true;
1.59 }
1.60 @@ -428,7 +425,7 @@
1.61 _next_arc(0)
1.62 {
1.63 // The main parameters of the pivot rule
1.64 - const double LIST_LENGTH_FACTOR = 1.0;
1.65 + const double LIST_LENGTH_FACTOR = 0.25;
1.66 const int MIN_LIST_LENGTH = 10;
1.67 const double MINOR_LIMIT_FACTOR = 0.1;
1.68 const int MIN_MINOR_LIMIT = 3;
1.69 @@ -445,7 +442,7 @@
1.70 /// Find next entering arc
1.71 bool findEnteringArc() {
1.72 Cost min, c;
1.73 - int e, min_arc = _next_arc;
1.74 + int e;
1.75 if (_curr_length > 0 && _minor_count < _minor_limit) {
1.76 // Minor iteration: select the best eligible arc from the
1.77 // current candidate list
1.78 @@ -456,16 +453,13 @@
1.79 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.80 if (c < min) {
1.81 min = c;
1.82 - min_arc = e;
1.83 + _in_arc = e;
1.84 }
1.85 - if (c >= 0) {
1.86 + else if (c >= 0) {
1.87 _candidates[i--] = _candidates[--_curr_length];
1.88 }
1.89 }
1.90 - if (min < 0) {
1.91 - _in_arc = min_arc;
1.92 - return true;
1.93 - }
1.94 + if (min < 0) return true;
1.95 }
1.96
1.97 // Major iteration: build a new candidate list
1.98 @@ -477,27 +471,26 @@
1.99 _candidates[_curr_length++] = e;
1.100 if (c < min) {
1.101 min = c;
1.102 - min_arc = e;
1.103 + _in_arc = e;
1.104 }
1.105 - if (_curr_length == _list_length) break;
1.106 + if (_curr_length == _list_length) goto search_end;
1.107 }
1.108 }
1.109 - if (_curr_length < _list_length) {
1.110 - for (e = 0; e < _next_arc; ++e) {
1.111 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.112 - if (c < 0) {
1.113 - _candidates[_curr_length++] = e;
1.114 - if (c < min) {
1.115 - min = c;
1.116 - min_arc = e;
1.117 - }
1.118 - if (_curr_length == _list_length) break;
1.119 + for (e = 0; e < _next_arc; ++e) {
1.120 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.121 + if (c < 0) {
1.122 + _candidates[_curr_length++] = e;
1.123 + if (c < min) {
1.124 + min = c;
1.125 + _in_arc = e;
1.126 }
1.127 + if (_curr_length == _list_length) goto search_end;
1.128 }
1.129 }
1.130 if (_curr_length == 0) return false;
1.131 +
1.132 + search_end:
1.133 _minor_count = 1;
1.134 - _in_arc = min_arc;
1.135 _next_arc = e;
1.136 return true;
1.137 }
1.138 @@ -549,7 +542,7 @@
1.139 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
1.140 {
1.141 // The main parameters of the pivot rule
1.142 - const double BLOCK_SIZE_FACTOR = 1.5;
1.143 + const double BLOCK_SIZE_FACTOR = 1.0;
1.144 const int MIN_BLOCK_SIZE = 10;
1.145 const double HEAD_LENGTH_FACTOR = 0.1;
1.146 const int MIN_HEAD_LENGTH = 3;
1.147 @@ -578,39 +571,35 @@
1.148
1.149 // Extend the list
1.150 int cnt = _block_size;
1.151 - int last_arc = 0;
1.152 int limit = _head_length;
1.153
1.154 - for (int e = _next_arc; e < _search_arc_num; ++e) {
1.155 + for (e = _next_arc; e < _search_arc_num; ++e) {
1.156 _cand_cost[e] = _state[e] *
1.157 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.158 if (_cand_cost[e] < 0) {
1.159 _candidates[_curr_length++] = e;
1.160 - last_arc = e;
1.161 }
1.162 if (--cnt == 0) {
1.163 - if (_curr_length > limit) break;
1.164 + if (_curr_length > limit) goto search_end;
1.165 limit = 0;
1.166 cnt = _block_size;
1.167 }
1.168 }
1.169 - if (_curr_length <= limit) {
1.170 - for (int e = 0; e < _next_arc; ++e) {
1.171 - _cand_cost[e] = _state[e] *
1.172 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.173 - if (_cand_cost[e] < 0) {
1.174 - _candidates[_curr_length++] = e;
1.175 - last_arc = e;
1.176 - }
1.177 - if (--cnt == 0) {
1.178 - if (_curr_length > limit) break;
1.179 - limit = 0;
1.180 - cnt = _block_size;
1.181 - }
1.182 + for (e = 0; e < _next_arc; ++e) {
1.183 + _cand_cost[e] = _state[e] *
1.184 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.185 + if (_cand_cost[e] < 0) {
1.186 + _candidates[_curr_length++] = e;
1.187 + }
1.188 + if (--cnt == 0) {
1.189 + if (_curr_length > limit) goto search_end;
1.190 + limit = 0;
1.191 + cnt = _block_size;
1.192 }
1.193 }
1.194 if (_curr_length == 0) return false;
1.195 - _next_arc = last_arc + 1;
1.196 +
1.197 + search_end:
1.198
1.199 // Make heap of the candidate list (approximating a partial sort)
1.200 make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.201 @@ -618,6 +607,7 @@
1.202
1.203 // Pop the first element of the heap
1.204 _in_arc = _candidates[0];
1.205 + _next_arc = e;
1.206 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.207 _sort_func );
1.208 _curr_length = std::min(_head_length, _curr_length - 1);
1.209 @@ -633,7 +623,11 @@
1.210 /// The constructor of the class.
1.211 ///
1.212 /// \param graph The digraph the algorithm runs on.
1.213 - NetworkSimplex(const GR& graph) :
1.214 + /// \param arc_mixing Indicate if the arcs have to be stored in a
1.215 + /// mixed order in the internal data structure.
1.216 + /// In special cases, it could lead to better overall performance,
1.217 + /// but it is usually slower. Therefore it is disabled by default.
1.218 + NetworkSimplex(const GR& graph, bool arc_mixing = false) :
1.219 _graph(graph), _node_id(graph), _arc_id(graph),
1.220 INF(std::numeric_limits<Value>::has_infinity ?
1.221 std::numeric_limits<Value>::infinity() :
1.222 @@ -671,31 +665,33 @@
1.223 _last_succ.resize(all_node_num);
1.224 _state.resize(max_arc_num);
1.225
1.226 - // Copy the graph (store the arcs in a mixed order)
1.227 + // Copy the graph
1.228 int i = 0;
1.229 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.230 _node_id[n] = i;
1.231 }
1.232 - int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.233 - i = 0;
1.234 - for (ArcIt a(_graph); a != INVALID; ++a) {
1.235 - _arc_id[a] = i;
1.236 - _source[i] = _node_id[_graph.source(a)];
1.237 - _target[i] = _node_id[_graph.target(a)];
1.238 - if ((i += k) >= _arc_num) i = (i % k) + 1;
1.239 + if (arc_mixing) {
1.240 + // Store the arcs in a mixed order
1.241 + int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.242 + int i = 0, j = 0;
1.243 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.244 + _arc_id[a] = i;
1.245 + _source[i] = _node_id[_graph.source(a)];
1.246 + _target[i] = _node_id[_graph.target(a)];
1.247 + if ((i += k) >= _arc_num) i = ++j;
1.248 + }
1.249 + } else {
1.250 + // Store the arcs in the original order
1.251 + int i = 0;
1.252 + for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
1.253 + _arc_id[a] = i;
1.254 + _source[i] = _node_id[_graph.source(a)];
1.255 + _target[i] = _node_id[_graph.target(a)];
1.256 + }
1.257 }
1.258
1.259 - // Initialize maps
1.260 - for (int i = 0; i != _node_num; ++i) {
1.261 - _supply[i] = 0;
1.262 - }
1.263 - for (int i = 0; i != _arc_num; ++i) {
1.264 - _lower[i] = 0;
1.265 - _upper[i] = INF;
1.266 - _cost[i] = 1;
1.267 - }
1.268 - _have_lower = false;
1.269 - _stype = GEQ;
1.270 + // Reset parameters
1.271 + reset();
1.272 }
1.273
1.274 /// \name Parameters
1.275 @@ -768,7 +764,6 @@
1.276 /// This function sets the supply values of the nodes.
1.277 /// If neither this function nor \ref stSupply() is used before
1.278 /// calling \ref run(), the supply of each node will be set to zero.
1.279 - /// (It makes sense only if non-zero lower bounds are given.)
1.280 ///
1.281 /// \param map A node map storing the supply values.
1.282 /// Its \c Value type must be convertible to the \c Value type
1.283 @@ -789,7 +784,6 @@
1.284 /// and the required flow value.
1.285 /// If neither this function nor \ref supplyMap() is used before
1.286 /// calling \ref run(), the supply of each node will be set to zero.
1.287 - /// (It makes sense only if non-zero lower bounds are given.)
1.288 ///
1.289 /// Using this function has the same effect as using \ref supplyMap()
1.290 /// with such a map in which \c k is assigned to \c s, \c -k is