211 IntVector _thread; |
212 IntVector _thread; |
212 IntVector _rev_thread; |
213 IntVector _rev_thread; |
213 IntVector _succ_num; |
214 IntVector _succ_num; |
214 IntVector _last_succ; |
215 IntVector _last_succ; |
215 IntVector _dirty_revs; |
216 IntVector _dirty_revs; |
216 CharVector _forward; |
217 BoolVector _forward; |
217 CharVector _state; |
218 BoolVector _state; |
218 int _root; |
219 int _root; |
219 |
220 |
220 // Temporary data used in the current pivot iteration |
221 // Temporary data used in the current pivot iteration |
221 int in_arc, join, u_in, v_in, u_out, v_out; |
222 int in_arc, join, u_in, v_in, u_out, v_out; |
222 int first, second, right, last; |
223 int first, second, right, last; |
353 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
354 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
354 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
355 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
355 _next_arc(0) |
356 _next_arc(0) |
356 { |
357 { |
357 // The main parameters of the pivot rule |
358 // The main parameters of the pivot rule |
358 const double BLOCK_SIZE_FACTOR = 0.5; |
359 const double BLOCK_SIZE_FACTOR = 1.0; |
359 const int MIN_BLOCK_SIZE = 10; |
360 const int MIN_BLOCK_SIZE = 10; |
360 |
361 |
361 _block_size = std::max( int(BLOCK_SIZE_FACTOR * |
362 _block_size = std::max( int(BLOCK_SIZE_FACTOR * |
362 std::sqrt(double(_search_arc_num))), |
363 std::sqrt(double(_search_arc_num))), |
363 MIN_BLOCK_SIZE ); |
364 MIN_BLOCK_SIZE ); |
366 // Find next entering arc |
367 // Find next entering arc |
367 bool findEnteringArc() { |
368 bool findEnteringArc() { |
368 Cost c, min = 0; |
369 Cost c, min = 0; |
369 int cnt = _block_size; |
370 int cnt = _block_size; |
370 int e; |
371 int e; |
371 for (e = _next_arc; e < _search_arc_num; ++e) { |
372 for (e = _next_arc; e != _search_arc_num; ++e) { |
372 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
373 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
373 if (c < min) { |
374 if (c < min) { |
374 min = c; |
375 min = c; |
375 _in_arc = e; |
376 _in_arc = e; |
376 } |
377 } |
377 if (--cnt == 0) { |
378 if (--cnt == 0) { |
378 if (min < 0) goto search_end; |
379 if (min < 0) goto search_end; |
379 cnt = _block_size; |
380 cnt = _block_size; |
380 } |
381 } |
381 } |
382 } |
382 for (e = 0; e < _next_arc; ++e) { |
383 for (e = 0; e != _next_arc; ++e) { |
383 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
384 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
384 if (c < min) { |
385 if (c < min) { |
385 min = c; |
386 min = c; |
386 _in_arc = e; |
387 _in_arc = e; |
387 } |
388 } |
468 } |
469 } |
469 |
470 |
470 // Major iteration: build a new candidate list |
471 // Major iteration: build a new candidate list |
471 min = 0; |
472 min = 0; |
472 _curr_length = 0; |
473 _curr_length = 0; |
473 for (e = _next_arc; e < _search_arc_num; ++e) { |
474 for (e = _next_arc; e != _search_arc_num; ++e) { |
474 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
475 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
475 if (c < 0) { |
476 if (c < 0) { |
476 _candidates[_curr_length++] = e; |
477 _candidates[_curr_length++] = e; |
477 if (c < min) { |
478 if (c < min) { |
478 min = c; |
479 min = c; |
479 _in_arc = e; |
480 _in_arc = e; |
480 } |
481 } |
481 if (_curr_length == _list_length) goto search_end; |
482 if (_curr_length == _list_length) goto search_end; |
482 } |
483 } |
483 } |
484 } |
484 for (e = 0; e < _next_arc; ++e) { |
485 for (e = 0; e != _next_arc; ++e) { |
485 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
486 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
486 if (c < 0) { |
487 if (c < 0) { |
487 _candidates[_curr_length++] = e; |
488 _candidates[_curr_length++] = e; |
488 if (c < min) { |
489 if (c < min) { |
489 min = c; |
490 min = c; |
563 |
564 |
564 // Find next entering arc |
565 // Find next entering arc |
565 bool findEnteringArc() { |
566 bool findEnteringArc() { |
566 // Check the current candidate list |
567 // Check the current candidate list |
567 int e; |
568 int e; |
568 for (int i = 0; i < _curr_length; ++i) { |
569 for (int i = 0; i != _curr_length; ++i) { |
569 e = _candidates[i]; |
570 e = _candidates[i]; |
570 _cand_cost[e] = _state[e] * |
571 _cand_cost[e] = _state[e] * |
571 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
572 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
572 if (_cand_cost[e] >= 0) { |
573 if (_cand_cost[e] >= 0) { |
573 _candidates[i--] = _candidates[--_curr_length]; |
574 _candidates[i--] = _candidates[--_curr_length]; |
576 |
577 |
577 // Extend the list |
578 // Extend the list |
578 int cnt = _block_size; |
579 int cnt = _block_size; |
579 int limit = _head_length; |
580 int limit = _head_length; |
580 |
581 |
581 for (e = _next_arc; e < _search_arc_num; ++e) { |
582 for (e = _next_arc; e != _search_arc_num; ++e) { |
582 _cand_cost[e] = _state[e] * |
583 _cand_cost[e] = _state[e] * |
583 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
584 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
584 if (_cand_cost[e] < 0) { |
585 if (_cand_cost[e] < 0) { |
585 _candidates[_curr_length++] = e; |
586 _candidates[_curr_length++] = e; |
586 } |
587 } |
1358 _thread[old_rev_thread] = right; |
1359 _thread[old_rev_thread] = right; |
1359 _rev_thread[right] = old_rev_thread; |
1360 _rev_thread[right] = old_rev_thread; |
1360 } |
1361 } |
1361 |
1362 |
1362 // Update _rev_thread using the new _thread values |
1363 // Update _rev_thread using the new _thread values |
1363 for (int i = 0; i < int(_dirty_revs.size()); ++i) { |
1364 for (int i = 0; i != int(_dirty_revs.size()); ++i) { |
1364 u = _dirty_revs[i]; |
1365 u = _dirty_revs[i]; |
1365 _rev_thread[_thread[u]] = u; |
1366 _rev_thread[_thread[u]] = u; |
1366 } |
1367 } |
1367 |
1368 |
1368 // Update _pred, _forward, _last_succ and _succ_num for the |
1369 // Update _pred, _forward, _last_succ and _succ_num for the |
1428 // Update potentials in the subtree, which has been moved |
1429 // Update potentials in the subtree, which has been moved |
1429 int end = _thread[_last_succ[u_in]]; |
1430 int end = _thread[_last_succ[u_in]]; |
1430 for (int u = u_in; u != end; u = _thread[u]) { |
1431 for (int u = u_in; u != end; u = _thread[u]) { |
1431 _pi[u] += sigma; |
1432 _pi[u] += sigma; |
1432 } |
1433 } |
|
1434 } |
|
1435 |
|
1436 // Heuristic initial pivots |
|
1437 bool initialPivots() { |
|
1438 Value curr, total = 0; |
|
1439 std::vector<Node> supply_nodes, demand_nodes; |
|
1440 for (NodeIt u(_graph); u != INVALID; ++u) { |
|
1441 curr = _supply[_node_id[u]]; |
|
1442 if (curr > 0) { |
|
1443 total += curr; |
|
1444 supply_nodes.push_back(u); |
|
1445 } |
|
1446 else if (curr < 0) { |
|
1447 demand_nodes.push_back(u); |
|
1448 } |
|
1449 } |
|
1450 if (_sum_supply > 0) total -= _sum_supply; |
|
1451 if (total <= 0) return true; |
|
1452 |
|
1453 IntVector arc_vector; |
|
1454 if (_sum_supply >= 0) { |
|
1455 if (supply_nodes.size() == 1 && demand_nodes.size() == 1) { |
|
1456 // Perform a reverse graph search from the sink to the source |
|
1457 typename GR::template NodeMap<bool> reached(_graph, false); |
|
1458 Node s = supply_nodes[0], t = demand_nodes[0]; |
|
1459 std::vector<Node> stack; |
|
1460 reached[t] = true; |
|
1461 stack.push_back(t); |
|
1462 while (!stack.empty()) { |
|
1463 Node u, v = stack.back(); |
|
1464 stack.pop_back(); |
|
1465 if (v == s) break; |
|
1466 for (InArcIt a(_graph, v); a != INVALID; ++a) { |
|
1467 if (reached[u = _graph.source(a)]) continue; |
|
1468 int j = _arc_id[a]; |
|
1469 if (_cap[j] >= total) { |
|
1470 arc_vector.push_back(j); |
|
1471 reached[u] = true; |
|
1472 stack.push_back(u); |
|
1473 } |
|
1474 } |
|
1475 } |
|
1476 } else { |
|
1477 // Find the min. cost incomming arc for each demand node |
|
1478 for (int i = 0; i != int(demand_nodes.size()); ++i) { |
|
1479 Node v = demand_nodes[i]; |
|
1480 Cost c, min_cost = std::numeric_limits<Cost>::max(); |
|
1481 Arc min_arc = INVALID; |
|
1482 for (InArcIt a(_graph, v); a != INVALID; ++a) { |
|
1483 c = _cost[_arc_id[a]]; |
|
1484 if (c < min_cost) { |
|
1485 min_cost = c; |
|
1486 min_arc = a; |
|
1487 } |
|
1488 } |
|
1489 if (min_arc != INVALID) { |
|
1490 arc_vector.push_back(_arc_id[min_arc]); |
|
1491 } |
|
1492 } |
|
1493 } |
|
1494 } else { |
|
1495 // Find the min. cost outgoing arc for each supply node |
|
1496 for (int i = 0; i != int(supply_nodes.size()); ++i) { |
|
1497 Node u = supply_nodes[i]; |
|
1498 Cost c, min_cost = std::numeric_limits<Cost>::max(); |
|
1499 Arc min_arc = INVALID; |
|
1500 for (OutArcIt a(_graph, u); a != INVALID; ++a) { |
|
1501 c = _cost[_arc_id[a]]; |
|
1502 if (c < min_cost) { |
|
1503 min_cost = c; |
|
1504 min_arc = a; |
|
1505 } |
|
1506 } |
|
1507 if (min_arc != INVALID) { |
|
1508 arc_vector.push_back(_arc_id[min_arc]); |
|
1509 } |
|
1510 } |
|
1511 } |
|
1512 |
|
1513 // Perform heuristic initial pivots |
|
1514 for (int i = 0; i != int(arc_vector.size()); ++i) { |
|
1515 in_arc = arc_vector[i]; |
|
1516 if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] - |
|
1517 _pi[_target[in_arc]]) >= 0) continue; |
|
1518 findJoinNode(); |
|
1519 bool change = findLeavingArc(); |
|
1520 if (delta >= MAX) return false; |
|
1521 changeFlow(change); |
|
1522 if (change) { |
|
1523 updateTreeStructure(); |
|
1524 updatePotential(); |
|
1525 } |
|
1526 } |
|
1527 return true; |
1433 } |
1528 } |
1434 |
1529 |
1435 // Execute the algorithm |
1530 // Execute the algorithm |
1436 ProblemType start(PivotRule pivot_rule) { |
1531 ProblemType start(PivotRule pivot_rule) { |
1437 // Select the pivot rule implementation |
1532 // Select the pivot rule implementation |
1452 |
1547 |
1453 template <typename PivotRuleImpl> |
1548 template <typename PivotRuleImpl> |
1454 ProblemType start() { |
1549 ProblemType start() { |
1455 PivotRuleImpl pivot(*this); |
1550 PivotRuleImpl pivot(*this); |
1456 |
1551 |
|
1552 // Perform heuristic initial pivots |
|
1553 if (!initialPivots()) return UNBOUNDED; |
|
1554 |
1457 // Execute the Network Simplex algorithm |
1555 // Execute the Network Simplex algorithm |
1458 while (pivot.findEnteringArc()) { |
1556 while (pivot.findEnteringArc()) { |
1459 findJoinNode(); |
1557 findJoinNode(); |
1460 bool change = findLeavingArc(); |
1558 bool change = findLeavingArc(); |
1461 if (delta >= MAX) return UNBOUNDED; |
1559 if (delta >= MAX) return UNBOUNDED; |