COIN-OR::LEMON - Graph Library

Changeset 655:6ac5d9ae1d3d in lemon


Ignore:
Timestamp:
04/03/09 18:59:15 (15 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Message:

Support real types + numerical stability fix in NS (#254)

  • Real types are supported by appropriate inicialization.
  • A feature of the XTI spanning tree structure is removed to ensure numerical stability (could cause problems using integer types). The node potentials are updated always on the lower subtree, in order to prevent overflow problems. The former method isn't notably faster during to our tests.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r654 r655  
    5555  /// algorithm. By default it is the same as \c F.
    5656  ///
    57   /// \warning Both value types must be signed integer types.
     57  /// \warning Both value types must be signed and all input data must
     58  /// be integer.
    5859  ///
    5960  /// \note %NetworkSimplex provides five different pivot rule
     
    10451046
    10461047      // Initialize arc maps
    1047       Flow max_cap = std::numeric_limits<Flow>::max();
    1048       Cost max_cost = std::numeric_limits<Cost>::max() / 4;
     1048      Flow inf_cap =
     1049        std::numeric_limits<Flow>::has_infinity ?
     1050        std::numeric_limits<Flow>::infinity() :
     1051        std::numeric_limits<Flow>::max();
    10491052      if (_pupper && _pcost) {
    10501053        for (int i = 0; i != _arc_num; ++i) {
     
    10701073        } else {
    10711074          for (int i = 0; i != _arc_num; ++i)
    1072             _cap[i] = max_cap;
     1075            _cap[i] = inf_cap;
    10731076        }
    10741077        if (_pcost) {
     
    10791082            _cost[i] = 1;
    10801083        }
     1084      }
     1085     
     1086      // Initialize artifical cost
     1087      Cost art_cost;
     1088      if (std::numeric_limits<Cost>::is_exact) {
     1089        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
     1090      } else {
     1091        art_cost = std::numeric_limits<Cost>::min();
     1092        for (int i = 0; i != _arc_num; ++i) {
     1093          if (_cost[i] > art_cost) art_cost = _cost[i];
     1094        }
     1095        art_cost = (art_cost + 1) * _node_num;
    10811096      }
    10821097
     
    11011116        _parent[u] = _root;
    11021117        _pred[u] = e;
    1103         _cost[e] = max_cost;
    1104         _cap[e] = max_cap;
     1118        _cost[e] = art_cost;
     1119        _cap[e] = inf_cap;
    11051120        _state[e] = STATE_TREE;
    11061121        if (_supply[u] >= 0) {
    11071122          _flow[e] = _supply[u];
    11081123          _forward[u] = true;
    1109           _pi[u] = -max_cost;
     1124          _pi[u] = -art_cost;
    11101125        } else {
    11111126          _flow[e] = -_supply[u];
    11121127          _forward[u] = false;
    1113           _pi[u] = max_cost;
     1128          _pi[u] = art_cost;
    11141129        }
    11151130      }
     
    13281343        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
    13291344        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
    1330       if (_succ_num[u_in] > _node_num / 2) {
    1331         // Update in the upper subtree (which contains the root)
    1332         int before = _rev_thread[u_in];
    1333         int after = _thread[_last_succ[u_in]];
    1334         _thread[before] = after;
    1335         _pi[_root] -= sigma;
    1336         for (int u = _thread[_root]; u != _root; u = _thread[u]) {
    1337           _pi[u] -= sigma;
    1338         }
    1339         _thread[before] = u_in;
    1340       } else {
    1341         // Update in the lower subtree (which has been moved)
    1342         int end = _thread[_last_succ[u_in]];
    1343         for (int u = u_in; u != end; u = _thread[u]) {
    1344           _pi[u] += sigma;
    1345         }
     1345      // Update potentials in the subtree, which has been moved
     1346      int end = _thread[_last_succ[u_in]];
     1347      for (int u = u_in; u != end; u = _thread[u]) {
     1348        _pi[u] += sigma;
    13461349      }
    13471350    }
Note: See TracChangeset for help on using the changeset viewer.