239 |
239 |
240 /// Constructor |
240 /// Constructor |
241 BlockSearchPivotRule(NetworkSimplex &ns) : |
241 BlockSearchPivotRule(NetworkSimplex &ns) : |
242 _edge(ns._edge), _source(ns._source), _target(ns._target), |
242 _edge(ns._edge), _source(ns._source), _target(ns._target), |
243 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
243 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
244 _in_edge(ns._in_edge), _edge_num(ns._edge_num + ns._node_num), _next_edge(0) |
244 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0) |
245 { |
245 { |
246 // The main parameters of the pivot rule |
246 // The main parameters of the pivot rule |
247 const double BLOCK_SIZE_FACTOR = 2.0; |
247 const double BLOCK_SIZE_FACTOR = 0.5; |
248 const int MIN_BLOCK_SIZE = 10; |
248 const int MIN_BLOCK_SIZE = 10; |
249 |
249 |
250 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)), |
250 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)), |
251 MIN_BLOCK_SIZE ); |
251 MIN_BLOCK_SIZE ); |
252 } |
252 } |
253 |
253 |
254 /// Find next entering edge |
254 /// Find next entering edge |
255 bool findEnteringEdge() { |
255 bool findEnteringEdge() { |
256 Cost c, min = 0; |
256 Cost c, min = 0; |
257 int cnt = _block_size; |
257 int cnt = _block_size; |
258 int e, min_edge = _next_edge; |
258 int e; |
259 for (e = _next_edge; e < _edge_num; ++e) { |
259 for (e = _next_edge; e < _edge_num; ++e) { |
260 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
260 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
261 if (c < min) { |
261 if (c < min) { |
262 min = c; |
262 min = c; |
263 min_edge = e; |
263 _in_edge = e; |
264 } |
264 } |
265 if (--cnt == 0) { |
265 if (--cnt == 0) { |
266 if (min < 0) break; |
266 if (min < 0) goto search_end; |
267 cnt = _block_size; |
267 cnt = _block_size; |
268 } |
268 } |
269 } |
269 } |
270 if (min == 0 || cnt > 0) { |
270 for (e = 0; e < _next_edge; ++e) { |
271 for (e = 0; e < _next_edge; ++e) { |
271 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
272 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
272 if (c < min) { |
273 if (c < min) { |
273 min = c; |
274 min = c; |
274 _in_edge = e; |
275 min_edge = e; |
275 } |
276 } |
276 if (--cnt == 0) { |
277 if (--cnt == 0) { |
277 if (min < 0) goto search_end; |
278 if (min < 0) break; |
278 cnt = _block_size; |
279 cnt = _block_size; |
|
280 } |
|
281 } |
279 } |
282 } |
280 } |
283 if (min >= 0) return false; |
281 if (min >= 0) return false; |
284 _in_edge = min_edge; |
282 |
|
283 search_end: |
285 _next_edge = e; |
284 _next_edge = e; |
286 return true; |
285 return true; |
287 } |
286 } |
288 |
287 |
289 }; //class BlockSearchPivotRule |
288 }; //class BlockSearchPivotRule |
323 _edge(ns._edge), _source(ns._source), _target(ns._target), |
322 _edge(ns._edge), _source(ns._source), _target(ns._target), |
324 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
323 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
325 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0) |
324 _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0) |
326 { |
325 { |
327 // The main parameters of the pivot rule |
326 // The main parameters of the pivot rule |
328 const double LIST_LENGTH_FACTOR = 1.0; |
327 const double LIST_LENGTH_FACTOR = 0.25; |
329 const int MIN_LIST_LENGTH = 10; |
328 const int MIN_LIST_LENGTH = 10; |
330 const double MINOR_LIMIT_FACTOR = 0.1; |
329 const double MINOR_LIMIT_FACTOR = 0.1; |
331 const int MIN_MINOR_LIMIT = 3; |
330 const int MIN_MINOR_LIMIT = 3; |
332 |
331 |
333 _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edge_num)), |
332 _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edge_num)), |
339 } |
338 } |
340 |
339 |
341 /// Find next entering edge |
340 /// Find next entering edge |
342 bool findEnteringEdge() { |
341 bool findEnteringEdge() { |
343 Cost min, c; |
342 Cost min, c; |
344 int e, min_edge = _next_edge; |
343 int e; |
345 if (_curr_length > 0 && _minor_count < _minor_limit) { |
344 if (_curr_length > 0 && _minor_count < _minor_limit) { |
346 // Minor iteration: select the best eligible edge from the |
345 // Minor iteration: select the best eligible edge from the |
347 // current candidate list |
346 // current candidate list |
348 ++_minor_count; |
347 ++_minor_count; |
349 min = 0; |
348 min = 0; |
350 for (int i = 0; i < _curr_length; ++i) { |
349 for (int i = 0; i < _curr_length; ++i) { |
351 e = _candidates[i]; |
350 e = _candidates[i]; |
352 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
351 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
353 if (c < min) { |
352 if (c < min) { |
354 min = c; |
353 min = c; |
355 min_edge = e; |
354 _in_edge = e; |
356 } |
355 } |
357 if (c >= 0) { |
356 else if (c >= 0) { |
358 _candidates[i--] = _candidates[--_curr_length]; |
357 _candidates[i--] = _candidates[--_curr_length]; |
359 } |
358 } |
360 } |
359 } |
361 if (min < 0) { |
360 if (min < 0) return true; |
362 _in_edge = min_edge; |
|
363 return true; |
|
364 } |
|
365 } |
361 } |
366 |
362 |
367 // Major iteration: build a new candidate list |
363 // Major iteration: build a new candidate list |
368 min = 0; |
364 min = 0; |
369 _curr_length = 0; |
365 _curr_length = 0; |
371 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
367 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
372 if (c < 0) { |
368 if (c < 0) { |
373 _candidates[_curr_length++] = e; |
369 _candidates[_curr_length++] = e; |
374 if (c < min) { |
370 if (c < min) { |
375 min = c; |
371 min = c; |
376 min_edge = e; |
372 _in_edge = e; |
377 } |
373 } |
378 if (_curr_length == _list_length) break; |
374 if (_curr_length == _list_length) goto search_end; |
379 } |
375 } |
380 } |
376 } |
381 if (_curr_length < _list_length) { |
377 for (e = 0; e < _next_edge; ++e) { |
382 for (e = 0; e < _next_edge; ++e) { |
378 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
383 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
379 if (c < 0) { |
384 if (c < 0) { |
380 _candidates[_curr_length++] = e; |
385 _candidates[_curr_length++] = e; |
381 if (c < min) { |
386 if (c < min) { |
382 min = c; |
387 min = c; |
383 _in_edge = e; |
388 min_edge = e; |
|
389 } |
|
390 if (_curr_length == _list_length) break; |
|
391 } |
384 } |
|
385 if (_curr_length == _list_length) goto search_end; |
392 } |
386 } |
393 } |
387 } |
394 if (_curr_length == 0) return false; |
388 if (_curr_length == 0) return false; |
|
389 |
|
390 search_end: |
395 _minor_count = 1; |
391 _minor_count = 1; |
396 _in_edge = min_edge; |
|
397 _next_edge = e; |
392 _next_edge = e; |
398 return true; |
393 return true; |
399 } |
394 } |
400 |
395 |
401 }; //class CandidateListPivotRule |
396 }; //class CandidateListPivotRule |
449 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
444 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
450 _in_edge(ns._in_edge), _edge_num(ns._edge_num), |
445 _in_edge(ns._in_edge), _edge_num(ns._edge_num), |
451 _next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost) |
446 _next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost) |
452 { |
447 { |
453 // The main parameters of the pivot rule |
448 // The main parameters of the pivot rule |
454 const double BLOCK_SIZE_FACTOR = 1.5; |
449 const double BLOCK_SIZE_FACTOR = 1.0; |
455 const int MIN_BLOCK_SIZE = 10; |
450 const int MIN_BLOCK_SIZE = 10; |
456 const double HEAD_LENGTH_FACTOR = 0.1; |
451 const double HEAD_LENGTH_FACTOR = 0.1; |
457 const int MIN_HEAD_LENGTH = 3; |
452 const int MIN_HEAD_LENGTH = 3; |
458 |
453 |
459 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)), |
454 _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)), |
477 } |
472 } |
478 } |
473 } |
479 |
474 |
480 // Extend the list |
475 // Extend the list |
481 int cnt = _block_size; |
476 int cnt = _block_size; |
482 int last_edge = 0; |
|
483 int limit = _head_length; |
477 int limit = _head_length; |
484 |
478 |
485 for (int e = _next_edge; e < _edge_num; ++e) { |
479 for (e = _next_edge; e < _edge_num; ++e) { |
486 _cand_cost[e] = _state[e] * |
480 _cand_cost[e] = _state[e] * |
487 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
481 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
488 if (_cand_cost[e] < 0) { |
482 if (_cand_cost[e] < 0) { |
489 _candidates[_curr_length++] = e; |
483 _candidates[_curr_length++] = e; |
490 last_edge = e; |
|
491 } |
484 } |
492 if (--cnt == 0) { |
485 if (--cnt == 0) { |
493 if (_curr_length > limit) break; |
486 if (_curr_length > limit) goto search_end; |
494 limit = 0; |
487 limit = 0; |
495 cnt = _block_size; |
488 cnt = _block_size; |
496 } |
489 } |
497 } |
490 } |
498 if (_curr_length <= limit) { |
491 for (e = 0; e < _next_edge; ++e) { |
499 for (int e = 0; e < _next_edge; ++e) { |
492 _cand_cost[e] = _state[e] * |
500 _cand_cost[e] = _state[e] * |
493 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
501 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
494 if (_cand_cost[e] < 0) { |
502 if (_cand_cost[e] < 0) { |
495 _candidates[_curr_length++] = e; |
503 _candidates[_curr_length++] = e; |
496 } |
504 last_edge = e; |
497 if (--cnt == 0) { |
505 } |
498 if (_curr_length > limit) goto search_end; |
506 if (--cnt == 0) { |
499 limit = 0; |
507 if (_curr_length > limit) break; |
500 cnt = _block_size; |
508 limit = 0; |
|
509 cnt = _block_size; |
|
510 } |
|
511 } |
501 } |
512 } |
502 } |
513 if (_curr_length == 0) return false; |
503 if (_curr_length == 0) return false; |
514 _next_edge = last_edge + 1; |
504 |
|
505 search_end: |
515 |
506 |
516 // Make heap of the candidate list (approximating a partial sort) |
507 // Make heap of the candidate list (approximating a partial sort) |
517 make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
508 make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
518 _sort_func ); |
509 _sort_func ); |
519 |
510 |
520 // Pop the first element of the heap |
511 // Pop the first element of the heap |
521 _in_edge = _candidates[0]; |
512 _in_edge = _candidates[0]; |
|
513 _next_edge = e; |
522 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
514 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
523 _sort_func ); |
515 _sort_func ); |
524 _curr_length = std::min(_head_length, _curr_length - 1); |
516 _curr_length = std::min(_head_length, _curr_length - 1); |
525 return true; |
517 return true; |
526 } |
518 } |
929 _succ_num[_root] = all_node_num; |
921 _succ_num[_root] = all_node_num; |
930 _last_succ[_root] = _root - 1; |
922 _last_succ[_root] = _root - 1; |
931 _supply[_root] = 0; |
923 _supply[_root] = 0; |
932 _pi[_root] = 0; |
924 _pi[_root] = 0; |
933 |
925 |
934 // Store the edges in a mixed order |
926 // Store the edges |
935 int k = std::max(int(sqrt(_edge_num)), 10); |
|
936 int i = 0; |
927 int i = 0; |
937 for (EdgeIt e(_orig_graph); e != INVALID; ++e) { |
928 for (EdgeIt e(_orig_graph); e != INVALID; ++e) { |
938 _edge[i] = e; |
929 _edge[i++] = e; |
939 if ((i += k) >= _edge_num) i = (i % k) + 1; |
|
940 } |
930 } |
941 |
931 |
942 // Initialize edge maps |
932 // Initialize edge maps |
943 for (int i = 0; i != _edge_num; ++i) { |
933 for (int i = 0; i != _edge_num; ++i) { |
944 Edge e = _edge[i]; |
934 Edge e = _edge[i]; |
959 } |
949 } |
960 } |
950 } |
961 } |
951 } |
962 |
952 |
963 // Add artificial edges and initialize the spanning tree data structure |
953 // Add artificial edges and initialize the spanning tree data structure |
964 Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
954 Cost max_cost = std::numeric_limits<Cost>::max() / 2 + 1; |
965 Capacity max_cap = std::numeric_limits<Capacity>::max(); |
955 Capacity max_cap = std::numeric_limits<Capacity>::max(); |
966 for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) { |
956 for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) { |
967 _parent[u] = _root; |
957 _parent[u] = _root; |
968 _pred[u] = e; |
958 _pred[u] = e; |
969 _thread[u] = u + 1; |
959 _thread[u] = u + 1; |