1.1 --- a/lemon/network_simplex.h Fri Sep 25 11:58:34 2009 +0200
1.2 +++ b/lemon/network_simplex.h Sat Sep 26 07:16:22 2009 +0200
1.3 @@ -362,33 +362,32 @@
1.4 bool findEnteringArc() {
1.5 Cost c, min = 0;
1.6 int cnt = _block_size;
1.7 - int e, min_arc = _next_arc;
1.8 + int e;
1.9 for (e = _next_arc; e < _search_arc_num; ++e) {
1.10 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.11 if (c < min) {
1.12 min = c;
1.13 - min_arc = e;
1.14 + _in_arc = e;
1.15 }
1.16 if (--cnt == 0) {
1.17 - if (min < 0) break;
1.18 + if (min < 0) goto search_end;
1.19 cnt = _block_size;
1.20 }
1.21 }
1.22 - if (min == 0 || cnt > 0) {
1.23 - for (e = 0; e < _next_arc; ++e) {
1.24 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.25 - if (c < min) {
1.26 - min = c;
1.27 - min_arc = e;
1.28 - }
1.29 - if (--cnt == 0) {
1.30 - if (min < 0) break;
1.31 - cnt = _block_size;
1.32 - }
1.33 + for (e = 0; e < _next_arc; ++e) {
1.34 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.35 + if (c < min) {
1.36 + min = c;
1.37 + _in_arc = e;
1.38 + }
1.39 + if (--cnt == 0) {
1.40 + if (min < 0) goto search_end;
1.41 + cnt = _block_size;
1.42 }
1.43 }
1.44 if (min >= 0) return false;
1.45 - _in_arc = min_arc;
1.46 +
1.47 + search_end:
1.48 _next_arc = e;
1.49 return true;
1.50 }
1.51 @@ -426,7 +425,7 @@
1.52 _next_arc(0)
1.53 {
1.54 // The main parameters of the pivot rule
1.55 - const double LIST_LENGTH_FACTOR = 1.0;
1.56 + const double LIST_LENGTH_FACTOR = 0.25;
1.57 const int MIN_LIST_LENGTH = 10;
1.58 const double MINOR_LIMIT_FACTOR = 0.1;
1.59 const int MIN_MINOR_LIMIT = 3;
1.60 @@ -443,7 +442,7 @@
1.61 /// Find next entering arc
1.62 bool findEnteringArc() {
1.63 Cost min, c;
1.64 - int e, min_arc = _next_arc;
1.65 + int e;
1.66 if (_curr_length > 0 && _minor_count < _minor_limit) {
1.67 // Minor iteration: select the best eligible arc from the
1.68 // current candidate list
1.69 @@ -454,16 +453,13 @@
1.70 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.71 if (c < min) {
1.72 min = c;
1.73 - min_arc = e;
1.74 + _in_arc = e;
1.75 }
1.76 - if (c >= 0) {
1.77 + else if (c >= 0) {
1.78 _candidates[i--] = _candidates[--_curr_length];
1.79 }
1.80 }
1.81 - if (min < 0) {
1.82 - _in_arc = min_arc;
1.83 - return true;
1.84 - }
1.85 + if (min < 0) return true;
1.86 }
1.87
1.88 // Major iteration: build a new candidate list
1.89 @@ -475,27 +471,26 @@
1.90 _candidates[_curr_length++] = e;
1.91 if (c < min) {
1.92 min = c;
1.93 - min_arc = e;
1.94 + _in_arc = e;
1.95 }
1.96 - if (_curr_length == _list_length) break;
1.97 + if (_curr_length == _list_length) goto search_end;
1.98 }
1.99 }
1.100 - if (_curr_length < _list_length) {
1.101 - for (e = 0; e < _next_arc; ++e) {
1.102 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.103 - if (c < 0) {
1.104 - _candidates[_curr_length++] = e;
1.105 - if (c < min) {
1.106 - min = c;
1.107 - min_arc = e;
1.108 - }
1.109 - if (_curr_length == _list_length) break;
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 + _in_arc = e;
1.117 }
1.118 + if (_curr_length == _list_length) goto search_end;
1.119 }
1.120 }
1.121 if (_curr_length == 0) return false;
1.122 +
1.123 + search_end:
1.124 _minor_count = 1;
1.125 - _in_arc = min_arc;
1.126 _next_arc = e;
1.127 return true;
1.128 }
1.129 @@ -547,7 +542,7 @@
1.130 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
1.131 {
1.132 // The main parameters of the pivot rule
1.133 - const double BLOCK_SIZE_FACTOR = 1.5;
1.134 + const double BLOCK_SIZE_FACTOR = 1.0;
1.135 const int MIN_BLOCK_SIZE = 10;
1.136 const double HEAD_LENGTH_FACTOR = 0.1;
1.137 const int MIN_HEAD_LENGTH = 3;
1.138 @@ -576,39 +571,35 @@
1.139
1.140 // Extend the list
1.141 int cnt = _block_size;
1.142 - int last_arc = 0;
1.143 int limit = _head_length;
1.144
1.145 - for (int e = _next_arc; e < _search_arc_num; ++e) {
1.146 + for (e = _next_arc; e < _search_arc_num; ++e) {
1.147 _cand_cost[e] = _state[e] *
1.148 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.149 if (_cand_cost[e] < 0) {
1.150 _candidates[_curr_length++] = e;
1.151 - last_arc = e;
1.152 }
1.153 if (--cnt == 0) {
1.154 - if (_curr_length > limit) break;
1.155 + if (_curr_length > limit) goto search_end;
1.156 limit = 0;
1.157 cnt = _block_size;
1.158 }
1.159 }
1.160 - if (_curr_length <= limit) {
1.161 - for (int e = 0; e < _next_arc; ++e) {
1.162 - _cand_cost[e] = _state[e] *
1.163 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.164 - if (_cand_cost[e] < 0) {
1.165 - _candidates[_curr_length++] = e;
1.166 - last_arc = e;
1.167 - }
1.168 - if (--cnt == 0) {
1.169 - if (_curr_length > limit) break;
1.170 - limit = 0;
1.171 - cnt = _block_size;
1.172 - }
1.173 + for (e = 0; e < _next_arc; ++e) {
1.174 + _cand_cost[e] = _state[e] *
1.175 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.176 + if (_cand_cost[e] < 0) {
1.177 + _candidates[_curr_length++] = e;
1.178 + }
1.179 + if (--cnt == 0) {
1.180 + if (_curr_length > limit) goto search_end;
1.181 + limit = 0;
1.182 + cnt = _block_size;
1.183 }
1.184 }
1.185 if (_curr_length == 0) return false;
1.186 - _next_arc = last_arc + 1;
1.187 +
1.188 + search_end:
1.189
1.190 // Make heap of the candidate list (approximating a partial sort)
1.191 make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.192 @@ -616,6 +607,7 @@
1.193
1.194 // Pop the first element of the heap
1.195 _in_arc = _candidates[0];
1.196 + _next_arc = e;
1.197 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.198 _sort_func );
1.199 _curr_length = std::min(_head_length, _curr_length - 1);
1.200 @@ -631,7 +623,11 @@
1.201 /// The constructor of the class.
1.202 ///
1.203 /// \param graph The digraph the algorithm runs on.
1.204 - NetworkSimplex(const GR& graph) :
1.205 + /// \param arc_mixing Indicate if the arcs have to be stored in a
1.206 + /// mixed order in the internal data structure.
1.207 + /// In special cases, it could lead to better overall performance,
1.208 + /// but it is usually slower. Therefore it is disabled by default.
1.209 + NetworkSimplex(const GR& graph, bool arc_mixing = false) :
1.210 _graph(graph), _node_id(graph), _arc_id(graph),
1.211 INF(std::numeric_limits<Value>::has_infinity ?
1.212 std::numeric_limits<Value>::infinity() :
1.213 @@ -669,18 +665,29 @@
1.214 _last_succ.resize(all_node_num);
1.215 _state.resize(max_arc_num);
1.216
1.217 - // Copy the graph (store the arcs in a mixed order)
1.218 + // Copy the graph
1.219 int i = 0;
1.220 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.221 _node_id[n] = i;
1.222 }
1.223 - int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.224 - i = 0;
1.225 - for (ArcIt a(_graph); a != INVALID; ++a) {
1.226 - _arc_id[a] = i;
1.227 - _source[i] = _node_id[_graph.source(a)];
1.228 - _target[i] = _node_id[_graph.target(a)];
1.229 - if ((i += k) >= _arc_num) i = (i % k) + 1;
1.230 + if (arc_mixing) {
1.231 + // Store the arcs in a mixed order
1.232 + int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.233 + int i = 0, j = 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 = ++j;
1.239 + }
1.240 + } else {
1.241 + // Store the arcs in the original order
1.242 + int i = 0;
1.243 + for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
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 + }
1.248 }
1.249
1.250 // Reset parameters