1.1 --- a/lemon/network_simplex.h Mon Jun 01 16:53:59 2009 +0000
1.2 +++ b/lemon/network_simplex.h Thu Jun 04 01:19:06 2009 +0000
1.3 @@ -241,10 +241,10 @@
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 + ns._node_num), _next_edge(0)
1.8 + _in_edge(ns._in_edge), _edge_num(ns._edge_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 + const double BLOCK_SIZE_FACTOR = 0.5;
1.13 const int MIN_BLOCK_SIZE = 10;
1.14
1.15 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
1.16 @@ -255,33 +255,32 @@
1.17 bool findEnteringEdge() {
1.18 Cost c, min = 0;
1.19 int cnt = _block_size;
1.20 - int e, min_edge = _next_edge;
1.21 + int e;
1.22 for (e = _next_edge; e < _edge_num; ++e) {
1.23 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.24 if (c < min) {
1.25 min = c;
1.26 - min_edge = e;
1.27 + _in_edge = e;
1.28 }
1.29 if (--cnt == 0) {
1.30 - if (min < 0) break;
1.31 + if (min < 0) goto search_end;
1.32 cnt = _block_size;
1.33 }
1.34 }
1.35 - if (min == 0 || cnt > 0) {
1.36 - for (e = 0; e < _next_edge; ++e) {
1.37 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.38 - if (c < min) {
1.39 - min = c;
1.40 - min_edge = e;
1.41 - }
1.42 - if (--cnt == 0) {
1.43 - if (min < 0) break;
1.44 - cnt = _block_size;
1.45 - }
1.46 + for (e = 0; e < _next_edge; ++e) {
1.47 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.48 + if (c < min) {
1.49 + min = c;
1.50 + _in_edge = e;
1.51 + }
1.52 + if (--cnt == 0) {
1.53 + if (min < 0) goto search_end;
1.54 + cnt = _block_size;
1.55 }
1.56 }
1.57 if (min >= 0) return false;
1.58 - _in_edge = min_edge;
1.59 +
1.60 + search_end:
1.61 _next_edge = e;
1.62 return true;
1.63 }
1.64 @@ -325,7 +324,7 @@
1.65 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
1.66 {
1.67 // The main parameters of the pivot rule
1.68 - const double LIST_LENGTH_FACTOR = 1.0;
1.69 + const double LIST_LENGTH_FACTOR = 0.25;
1.70 const int MIN_LIST_LENGTH = 10;
1.71 const double MINOR_LIMIT_FACTOR = 0.1;
1.72 const int MIN_MINOR_LIMIT = 3;
1.73 @@ -341,7 +340,7 @@
1.74 /// Find next entering edge
1.75 bool findEnteringEdge() {
1.76 Cost min, c;
1.77 - int e, min_edge = _next_edge;
1.78 + int e;
1.79 if (_curr_length > 0 && _minor_count < _minor_limit) {
1.80 // Minor iteration: select the best eligible edge from the
1.81 // current candidate list
1.82 @@ -352,16 +351,13 @@
1.83 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.84 if (c < min) {
1.85 min = c;
1.86 - min_edge = e;
1.87 + _in_edge = e;
1.88 }
1.89 - if (c >= 0) {
1.90 + else if (c >= 0) {
1.91 _candidates[i--] = _candidates[--_curr_length];
1.92 }
1.93 }
1.94 - if (min < 0) {
1.95 - _in_edge = min_edge;
1.96 - return true;
1.97 - }
1.98 + if (min < 0) return true;
1.99 }
1.100
1.101 // Major iteration: build a new candidate list
1.102 @@ -373,27 +369,26 @@
1.103 _candidates[_curr_length++] = e;
1.104 if (c < min) {
1.105 min = c;
1.106 - min_edge = e;
1.107 + _in_edge = e;
1.108 }
1.109 - if (_curr_length == _list_length) break;
1.110 + if (_curr_length == _list_length) goto search_end;
1.111 }
1.112 }
1.113 - if (_curr_length < _list_length) {
1.114 - for (e = 0; e < _next_edge; ++e) {
1.115 - c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.116 - if (c < 0) {
1.117 - _candidates[_curr_length++] = e;
1.118 - if (c < min) {
1.119 - min = c;
1.120 - min_edge = e;
1.121 - }
1.122 - if (_curr_length == _list_length) break;
1.123 + for (e = 0; e < _next_edge; ++e) {
1.124 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.125 + if (c < 0) {
1.126 + _candidates[_curr_length++] = e;
1.127 + if (c < min) {
1.128 + min = c;
1.129 + _in_edge = e;
1.130 }
1.131 + if (_curr_length == _list_length) goto search_end;
1.132 }
1.133 }
1.134 if (_curr_length == 0) return false;
1.135 +
1.136 + search_end:
1.137 _minor_count = 1;
1.138 - _in_edge = min_edge;
1.139 _next_edge = e;
1.140 return true;
1.141 }
1.142 @@ -451,7 +446,7 @@
1.143 _next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost)
1.144 {
1.145 // The main parameters of the pivot rule
1.146 - const double BLOCK_SIZE_FACTOR = 1.5;
1.147 + const double BLOCK_SIZE_FACTOR = 1.0;
1.148 const int MIN_BLOCK_SIZE = 10;
1.149 const double HEAD_LENGTH_FACTOR = 0.1;
1.150 const int MIN_HEAD_LENGTH = 3;
1.151 @@ -479,39 +474,35 @@
1.152
1.153 // Extend the list
1.154 int cnt = _block_size;
1.155 - int last_edge = 0;
1.156 int limit = _head_length;
1.157
1.158 - for (int e = _next_edge; e < _edge_num; ++e) {
1.159 + for (e = _next_edge; e < _edge_num; ++e) {
1.160 _cand_cost[e] = _state[e] *
1.161 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.162 if (_cand_cost[e] < 0) {
1.163 _candidates[_curr_length++] = e;
1.164 - last_edge = e;
1.165 }
1.166 if (--cnt == 0) {
1.167 - if (_curr_length > limit) break;
1.168 + if (_curr_length > limit) goto search_end;
1.169 limit = 0;
1.170 cnt = _block_size;
1.171 }
1.172 }
1.173 - if (_curr_length <= limit) {
1.174 - for (int e = 0; e < _next_edge; ++e) {
1.175 - _cand_cost[e] = _state[e] *
1.176 - (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.177 - if (_cand_cost[e] < 0) {
1.178 - _candidates[_curr_length++] = e;
1.179 - last_edge = e;
1.180 - }
1.181 - if (--cnt == 0) {
1.182 - if (_curr_length > limit) break;
1.183 - limit = 0;
1.184 - cnt = _block_size;
1.185 - }
1.186 + for (e = 0; e < _next_edge; ++e) {
1.187 + _cand_cost[e] = _state[e] *
1.188 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.189 + if (_cand_cost[e] < 0) {
1.190 + _candidates[_curr_length++] = e;
1.191 + }
1.192 + if (--cnt == 0) {
1.193 + if (_curr_length > limit) goto search_end;
1.194 + limit = 0;
1.195 + cnt = _block_size;
1.196 }
1.197 }
1.198 if (_curr_length == 0) return false;
1.199 - _next_edge = last_edge + 1;
1.200 +
1.201 + search_end:
1.202
1.203 // Make heap of the candidate list (approximating a partial sort)
1.204 make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.205 @@ -519,6 +510,7 @@
1.206
1.207 // Pop the first element of the heap
1.208 _in_edge = _candidates[0];
1.209 + _next_edge = e;
1.210 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.211 _sort_func );
1.212 _curr_length = std::min(_head_length, _curr_length - 1);
1.213 @@ -931,12 +923,10 @@
1.214 _supply[_root] = 0;
1.215 _pi[_root] = 0;
1.216
1.217 - // Store the edges in a mixed order
1.218 - int k = std::max(int(sqrt(_edge_num)), 10);
1.219 + // Store the edges
1.220 int i = 0;
1.221 for (EdgeIt e(_orig_graph); e != INVALID; ++e) {
1.222 - _edge[i] = e;
1.223 - if ((i += k) >= _edge_num) i = (i % k) + 1;
1.224 + _edge[i++] = e;
1.225 }
1.226
1.227 // Initialize edge maps
1.228 @@ -961,7 +951,7 @@
1.229 }
1.230
1.231 // Add artificial edges and initialize the spanning tree data structure
1.232 - Cost max_cost = std::numeric_limits<Cost>::max() / 4;
1.233 + Cost max_cost = std::numeric_limits<Cost>::max() / 2 + 1;
1.234 Capacity max_cap = std::numeric_limits<Capacity>::max();
1.235 for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
1.236 _parent[u] = _root;