210 IntVector _thread; |
211 IntVector _thread; |
211 IntVector _rev_thread; |
212 IntVector _rev_thread; |
212 IntVector _succ_num; |
213 IntVector _succ_num; |
213 IntVector _last_succ; |
214 IntVector _last_succ; |
214 IntVector _dirty_revs; |
215 IntVector _dirty_revs; |
215 CharVector _forward; |
216 BoolVector _forward; |
216 CharVector _state; |
217 BoolVector _state; |
217 int _root; |
218 int _root; |
218 |
219 |
219 // Temporary data used in the current pivot iteration |
220 // Temporary data used in the current pivot iteration |
220 int in_arc, join, u_in, v_in, u_out, v_out; |
221 int in_arc, join, u_in, v_in, u_out, v_out; |
221 int first, second, right, last; |
222 int first, second, right, last; |
352 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
353 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
353 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
354 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
354 _next_arc(0) |
355 _next_arc(0) |
355 { |
356 { |
356 // The main parameters of the pivot rule |
357 // The main parameters of the pivot rule |
357 const double BLOCK_SIZE_FACTOR = 0.5; |
358 const double BLOCK_SIZE_FACTOR = 1.0; |
358 const int MIN_BLOCK_SIZE = 10; |
359 const int MIN_BLOCK_SIZE = 10; |
359 |
360 |
360 _block_size = std::max( int(BLOCK_SIZE_FACTOR * |
361 _block_size = std::max( int(BLOCK_SIZE_FACTOR * |
361 std::sqrt(double(_search_arc_num))), |
362 std::sqrt(double(_search_arc_num))), |
362 MIN_BLOCK_SIZE ); |
363 MIN_BLOCK_SIZE ); |
365 // Find next entering arc |
366 // Find next entering arc |
366 bool findEnteringArc() { |
367 bool findEnteringArc() { |
367 Cost c, min = 0; |
368 Cost c, min = 0; |
368 int cnt = _block_size; |
369 int cnt = _block_size; |
369 int e; |
370 int e; |
370 for (e = _next_arc; e < _search_arc_num; ++e) { |
371 for (e = _next_arc; e != _search_arc_num; ++e) { |
371 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
372 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
372 if (c < min) { |
373 if (c < min) { |
373 min = c; |
374 min = c; |
374 _in_arc = e; |
375 _in_arc = e; |
375 } |
376 } |
376 if (--cnt == 0) { |
377 if (--cnt == 0) { |
377 if (min < 0) goto search_end; |
378 if (min < 0) goto search_end; |
378 cnt = _block_size; |
379 cnt = _block_size; |
379 } |
380 } |
380 } |
381 } |
381 for (e = 0; e < _next_arc; ++e) { |
382 for (e = 0; e != _next_arc; ++e) { |
382 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
383 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
383 if (c < min) { |
384 if (c < min) { |
384 min = c; |
385 min = c; |
385 _in_arc = e; |
386 _in_arc = e; |
386 } |
387 } |
467 } |
468 } |
468 |
469 |
469 // Major iteration: build a new candidate list |
470 // Major iteration: build a new candidate list |
470 min = 0; |
471 min = 0; |
471 _curr_length = 0; |
472 _curr_length = 0; |
472 for (e = _next_arc; e < _search_arc_num; ++e) { |
473 for (e = _next_arc; e != _search_arc_num; ++e) { |
473 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
474 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
474 if (c < 0) { |
475 if (c < 0) { |
475 _candidates[_curr_length++] = e; |
476 _candidates[_curr_length++] = e; |
476 if (c < min) { |
477 if (c < min) { |
477 min = c; |
478 min = c; |
478 _in_arc = e; |
479 _in_arc = e; |
479 } |
480 } |
480 if (_curr_length == _list_length) goto search_end; |
481 if (_curr_length == _list_length) goto search_end; |
481 } |
482 } |
482 } |
483 } |
483 for (e = 0; e < _next_arc; ++e) { |
484 for (e = 0; e != _next_arc; ++e) { |
484 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
485 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
485 if (c < 0) { |
486 if (c < 0) { |
486 _candidates[_curr_length++] = e; |
487 _candidates[_curr_length++] = e; |
487 if (c < min) { |
488 if (c < min) { |
488 min = c; |
489 min = c; |
562 |
563 |
563 // Find next entering arc |
564 // Find next entering arc |
564 bool findEnteringArc() { |
565 bool findEnteringArc() { |
565 // Check the current candidate list |
566 // Check the current candidate list |
566 int e; |
567 int e; |
567 for (int i = 0; i < _curr_length; ++i) { |
568 for (int i = 0; i != _curr_length; ++i) { |
568 e = _candidates[i]; |
569 e = _candidates[i]; |
569 _cand_cost[e] = _state[e] * |
570 _cand_cost[e] = _state[e] * |
570 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
571 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
571 if (_cand_cost[e] >= 0) { |
572 if (_cand_cost[e] >= 0) { |
572 _candidates[i--] = _candidates[--_curr_length]; |
573 _candidates[i--] = _candidates[--_curr_length]; |
575 |
576 |
576 // Extend the list |
577 // Extend the list |
577 int cnt = _block_size; |
578 int cnt = _block_size; |
578 int limit = _head_length; |
579 int limit = _head_length; |
579 |
580 |
580 for (e = _next_arc; e < _search_arc_num; ++e) { |
581 for (e = _next_arc; e != _search_arc_num; ++e) { |
581 _cand_cost[e] = _state[e] * |
582 _cand_cost[e] = _state[e] * |
582 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
583 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
583 if (_cand_cost[e] < 0) { |
584 if (_cand_cost[e] < 0) { |
584 _candidates[_curr_length++] = e; |
585 _candidates[_curr_length++] = e; |
585 } |
586 } |
1326 _thread[old_rev_thread] = right; |
1327 _thread[old_rev_thread] = right; |
1327 _rev_thread[right] = old_rev_thread; |
1328 _rev_thread[right] = old_rev_thread; |
1328 } |
1329 } |
1329 |
1330 |
1330 // Update _rev_thread using the new _thread values |
1331 // Update _rev_thread using the new _thread values |
1331 for (int i = 0; i < int(_dirty_revs.size()); ++i) { |
1332 for (int i = 0; i != int(_dirty_revs.size()); ++i) { |
1332 u = _dirty_revs[i]; |
1333 u = _dirty_revs[i]; |
1333 _rev_thread[_thread[u]] = u; |
1334 _rev_thread[_thread[u]] = u; |
1334 } |
1335 } |
1335 |
1336 |
1336 // Update _pred, _forward, _last_succ and _succ_num for the |
1337 // Update _pred, _forward, _last_succ and _succ_num for the |
1396 // Update potentials in the subtree, which has been moved |
1397 // Update potentials in the subtree, which has been moved |
1397 int end = _thread[_last_succ[u_in]]; |
1398 int end = _thread[_last_succ[u_in]]; |
1398 for (int u = u_in; u != end; u = _thread[u]) { |
1399 for (int u = u_in; u != end; u = _thread[u]) { |
1399 _pi[u] += sigma; |
1400 _pi[u] += sigma; |
1400 } |
1401 } |
|
1402 } |
|
1403 |
|
1404 // Heuristic initial pivots |
|
1405 bool initialPivots() { |
|
1406 Value curr, total = 0; |
|
1407 std::vector<Node> supply_nodes, demand_nodes; |
|
1408 for (NodeIt u(_graph); u != INVALID; ++u) { |
|
1409 curr = _supply[_node_id[u]]; |
|
1410 if (curr > 0) { |
|
1411 total += curr; |
|
1412 supply_nodes.push_back(u); |
|
1413 } |
|
1414 else if (curr < 0) { |
|
1415 demand_nodes.push_back(u); |
|
1416 } |
|
1417 } |
|
1418 if (_sum_supply > 0) total -= _sum_supply; |
|
1419 if (total <= 0) return true; |
|
1420 |
|
1421 IntVector arc_vector; |
|
1422 if (_sum_supply >= 0) { |
|
1423 if (supply_nodes.size() == 1 && demand_nodes.size() == 1) { |
|
1424 // Perform a reverse graph search from the sink to the source |
|
1425 typename GR::template NodeMap<bool> reached(_graph, false); |
|
1426 Node s = supply_nodes[0], t = demand_nodes[0]; |
|
1427 std::vector<Node> stack; |
|
1428 reached[t] = true; |
|
1429 stack.push_back(t); |
|
1430 while (!stack.empty()) { |
|
1431 Node u, v = stack.back(); |
|
1432 stack.pop_back(); |
|
1433 if (v == s) break; |
|
1434 for (InArcIt a(_graph, v); a != INVALID; ++a) { |
|
1435 if (reached[u = _graph.source(a)]) continue; |
|
1436 int j = _arc_id[a]; |
|
1437 if (_cap[j] >= total) { |
|
1438 arc_vector.push_back(j); |
|
1439 reached[u] = true; |
|
1440 stack.push_back(u); |
|
1441 } |
|
1442 } |
|
1443 } |
|
1444 } else { |
|
1445 // Find the min. cost incomming arc for each demand node |
|
1446 for (int i = 0; i != int(demand_nodes.size()); ++i) { |
|
1447 Node v = demand_nodes[i]; |
|
1448 Cost c, min_cost = std::numeric_limits<Cost>::max(); |
|
1449 Arc min_arc = INVALID; |
|
1450 for (InArcIt a(_graph, v); a != INVALID; ++a) { |
|
1451 c = _cost[_arc_id[a]]; |
|
1452 if (c < min_cost) { |
|
1453 min_cost = c; |
|
1454 min_arc = a; |
|
1455 } |
|
1456 } |
|
1457 if (min_arc != INVALID) { |
|
1458 arc_vector.push_back(_arc_id[min_arc]); |
|
1459 } |
|
1460 } |
|
1461 } |
|
1462 } else { |
|
1463 // Find the min. cost outgoing arc for each supply node |
|
1464 for (int i = 0; i != int(supply_nodes.size()); ++i) { |
|
1465 Node u = supply_nodes[i]; |
|
1466 Cost c, min_cost = std::numeric_limits<Cost>::max(); |
|
1467 Arc min_arc = INVALID; |
|
1468 for (OutArcIt a(_graph, u); a != INVALID; ++a) { |
|
1469 c = _cost[_arc_id[a]]; |
|
1470 if (c < min_cost) { |
|
1471 min_cost = c; |
|
1472 min_arc = a; |
|
1473 } |
|
1474 } |
|
1475 if (min_arc != INVALID) { |
|
1476 arc_vector.push_back(_arc_id[min_arc]); |
|
1477 } |
|
1478 } |
|
1479 } |
|
1480 |
|
1481 // Perform heuristic initial pivots |
|
1482 for (int i = 0; i != int(arc_vector.size()); ++i) { |
|
1483 in_arc = arc_vector[i]; |
|
1484 if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] - |
|
1485 _pi[_target[in_arc]]) >= 0) continue; |
|
1486 findJoinNode(); |
|
1487 bool change = findLeavingArc(); |
|
1488 if (delta >= MAX) return false; |
|
1489 changeFlow(change); |
|
1490 if (change) { |
|
1491 updateTreeStructure(); |
|
1492 updatePotential(); |
|
1493 } |
|
1494 } |
|
1495 return true; |
1401 } |
1496 } |
1402 |
1497 |
1403 // Execute the algorithm |
1498 // Execute the algorithm |
1404 ProblemType start(PivotRule pivot_rule) { |
1499 ProblemType start(PivotRule pivot_rule) { |
1405 // Select the pivot rule implementation |
1500 // Select the pivot rule implementation |
1420 |
1515 |
1421 template <typename PivotRuleImpl> |
1516 template <typename PivotRuleImpl> |
1422 ProblemType start() { |
1517 ProblemType start() { |
1423 PivotRuleImpl pivot(*this); |
1518 PivotRuleImpl pivot(*this); |
1424 |
1519 |
|
1520 // Perform heuristic initial pivots |
|
1521 if (!initialPivots()) return UNBOUNDED; |
|
1522 |
1425 // Execute the Network Simplex algorithm |
1523 // Execute the Network Simplex algorithm |
1426 while (pivot.findEnteringArc()) { |
1524 while (pivot.findEnteringArc()) { |
1427 findJoinNode(); |
1525 findJoinNode(); |
1428 bool change = findLeavingArc(); |
1526 bool change = findLeavingArc(); |
1429 if (delta >= MAX) return UNBOUNDED; |
1527 if (delta >= MAX) return UNBOUNDED; |