Implementation stuff.

12 files changed, 326 insertions, 4463 deletions

diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp
index 664d71f..73c8f23 100644
--- a/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp
+++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp
@@ -4,9 +4,15 @@
 #include <cassert>
 #include <ostream>
 #include <istream>
+#include <limits>
 
 template<typename T>
-T infinity() { }
+T infinity();
+
+template<>
+double infinity() {
+  return std::numeric_limits<double>::infinity();
+}
 
 template<typename T>
 struct Complete {
diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp
index 5ee5405..3342cc7 100644
--- a/clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp
+++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp
@@ -76,14 +76,9 @@ struct EquationSystem {
   }
 
   Constant<Domain>& constant(const Domain& value) {
-    if (_constants.find(value) == _constants.end()) {
-      Constant<Domain>* constant = new Constant<Domain>(value);
-      _expressions.insert(constant);
-      _constants[value] = constant;
-      return *constant;
-    } else {
-      return *_constants[value];
-    }
+    Constant<Domain>* constant = new Constant<Domain>(value);
+    _expressions.insert(constant);
+    return *constant;
   }
 
   MaxExpression<Domain>* operator[](const Variable<Domain>& var) const {
@@ -132,7 +127,6 @@ struct EquationSystem {
   private:
   std::set<Operator<Domain>*> _operators;
   std::set<Expression<Domain>*> _expressions;
-  std::map<Domain,Constant<Domain>*> _constants;
   std::vector<Variable<Domain>*> _variables;
   std::map<std::string, Variable<Domain>*> _variable_names;
   IdMap<MaxExpression<Domain>, Variable<Domain>*>* _expr_to_var;
diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/Log.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/Log.hpp
index 90ab7e5..7502ac8 100644
--- a/clang/include/clang/Analysis/Analyses/IntervalSolver/Log.hpp
+++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/Log.hpp
@@ -1,17 +1,13 @@
 #ifndef LOG_HPP
 #define LOG_HPP
 
-// could not be hackier, but C++ is annoying
-#define protected public
-#include <streambuf>
-#undef protected
-
 #include <string>
 #include <iostream>
 #include <map>
 #include <cstdio>
 
 namespace solver_log {
+
   struct Logger : public std::ostream {
     Logger(std::streambuf* buffer, const std::string& name)
       : std::ostream(NULL) { }
diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/MaxStrategy.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/MaxStrategy.hpp
index 9ae7394..da05dd8 100644
--- a/clang/include/clang/Analysis/Analyses/IntervalSolver/MaxStrategy.hpp
+++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/MaxStrategy.hpp
@@ -39,22 +39,13 @@ struct DynamicMaxStrategy : public MaxStrategy<Domain> {
       _values(system.maxExpressionCount(), 0),
       _stable(system.maxExpressionCount()),
       _influence(system.maxExpressionCount(),
-                 IdSet<MaxExpression<Domain> >(system.maxExpressionCount())),
-      _changed(false)
+                 IdSet<MaxExpression<Domain> >(system.maxExpressionCount()))
   {}
 
   void setRho(DynamicVariableAssignment<Domain>& rho) {
     _rho = &rho;
   }
 
-  void hasChanged(bool c) {
-    _changed = c;
-  }
-
-  bool hasChanged() const {
-    return _changed;
-  }
-
   unsigned int get(const MaxExpression<Domain>& e) const {
     solver_log::strategy << indent() << "Look up " << e << std::endl;
     return _values[e];
@@ -66,19 +57,6 @@ struct DynamicMaxStrategy : public MaxStrategy<Domain> {
     return _values[e];
   }
 
-  void invalidate(const Variable<Domain>& v) {
-    solver_log::strategy << indent() << "Invalidating " << v << " - " << *_system[v] << std::endl; 
-
-    IdSet<MaxExpression<Domain> > infl = _influence[*_system[v]];
-    for (typename IdSet<MaxExpression<Domain> >::iterator
-           it = infl.begin(),
-           end = infl.end();
-         it != end;
-         ++it) {
-      invalidate(_system.maxExpression(*it));
-    }
-  }
-
   void invalidate(const MaxExpression<Domain>& v) {
     if (_stable.contains(v)) {
       solver_log::strategy << indent() << "Invalidating " << v << std::endl; 
@@ -111,20 +89,18 @@ private:
       if (val != _values[x]) {
         solver_log::strategy << x << " => " << *x.arguments()[val] << std::endl;
 
-        IdSet<MaxExpression<Domain> > oldInfluence = _influence[x];
-        //_influence[x].clear();
         _values[x] = val;
-        _changed = true;
 
         _rho->invalidate(*_system.varFromExpr(x));
         
 	//_rho->stabilise();
 
-        _stable.filter(oldInfluence);
-
+        IdSet<MaxExpression<Domain> > infl = _influence[x];
+        _stable.filter(infl); 
+        
         for (typename IdSet<MaxExpression<Domain> >::iterator
-               it = oldInfluence.begin(),
-               end = oldInfluence.end();
+               it = infl.begin(),
+               end = infl.end();
              it != end;
              ++it) {
           solve(_system.maxExpression(*it));
@@ -141,7 +117,7 @@ private:
 
   struct DependencyAssignment : public VariableAssignment<Domain>{
     DependencyAssignment(DynamicMaxStrategy& strat,
-                         VariableAssignment<Domain>& rho,
+                         DynamicVariableAssignment<Domain>& rho,
                          const MaxExpression<Domain>& expr)
       : _strat(strat),
         _rho(rho),
@@ -149,14 +125,25 @@ private:
 	_current(strat._system.variableCount()) {
     }
     const Domain operator[](const Variable<Domain>& var) {
-      // solve the strategy for this variable, too
-      _strat.solve(*_strat._system[var]);
+      // force evaluation to get touched variables
+      Domain value = _rho[var];
+
       _strat._influence[*_strat._system[var]].insert(_expr);
-      return _rho[var];
+
+      // invalidate touched variables
+      IdSet<Variable<Domain> > changed = _rho.get_changed();
+      for (typename IdSet<Variable<Domain> >::iterator
+             it = changed.begin(),
+             ei = changed.end();
+           it != ei;
+           ++it) {
+        _strat.invalidate(*_strat._system[_strat._system.variable(*it)]);
+      }
+      return value;
     }
   private:
     DynamicMaxStrategy& _strat;
-    VariableAssignment<Domain>& _rho;
+    DynamicVariableAssignment<Domain>& _rho;
     const MaxExpression<Domain>& _expr;
     IdSet<Variable<Domain> > _current;
   };
@@ -186,7 +173,6 @@ private:
   IdMap<MaxExpression<Domain>,unsigned int> _values;
   IdSet<MaxExpression<Domain> > _stable;
   IdMap<MaxExpression<Domain>,IdSet<MaxExpression<Domain> > > _influence;
-  bool _changed;
 };
 
 
@@ -200,28 +186,12 @@ struct Solver {
   }
 
   Domain solve(Variable<Domain>& var) {
-    MaxExpression<Domain>& rhs = *_system[var];
-
-    // _rho automatically keeps up to date with changes made to the
-    // strategy, so you don't need to stabilise it
-
+      MaxExpression<Domain>& rhs = *_system[var];
+    // this will automatically work sufficiently to get the final
+    // strategy for this variable's RHS
     _strategy.get(rhs);
 
-
-    /*
-    do {
-      _strategy.hasChanged(false);
-      
-      solver_log::debug << "Stabilise assignment (toplevel)" << std::endl;
-      _rho.stabilise();
-      
-      solver_log::debug << "Improve strategy (toplevel)" << std::endl;
-      // improve strategy
-      _strategy.get(rhs);
-      solver_log::debug << (_strategy.hasChanged() ? "Strategy has changed - loop again" : "Strategy has not changed - terminate") << std::endl;
-    } while (_strategy.hasChanged());
-    */
-
+    // this will automatically solve for the value of the RHS, if required
     return _rho[var];
   }
 private:
diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/Operator.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/Operator.hpp
index 64ef096..0b08866 100644
--- a/clang/include/clang/Analysis/Analyses/IntervalSolver/Operator.hpp
+++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/Operator.hpp
@@ -135,6 +135,88 @@ struct Guard : public Operator<Domain> {
   }
 };
 
+#include "MCFSimplex.h"
+
+template<typename Domain>
+struct MinCostFlow : public Operator<Domain> {
+  MinCostFlow(const std::vector<Domain>& supplies, const std::vector<std::pair<int,int> > arcs)
+    : _supplies(supplies), _arcs(arcs), _solver(0,0) {
+    MCFClass::FNumber* deficits = new MCFClass::FNumber[supplies.size()];
+    MCFClass::Index* starts = new MCFClass::Index[arcs.size()];
+    MCFClass::Index* ends = new MCFClass::Index[arcs.size()];
+
+    for (int i = 0, size = supplies.size(); i < size; ++i) {
+      deficits[i] = -supplies[i].template as<MCFClass::FNumber>();
+    }
+    for (int i = 0, size = arcs.size(); i < size; ++i) {
+      starts[i] = arcs[i].first;
+      ends[i] = arcs[i].second;
+    }
+
+    _solver.LoadNet(supplies.size(), arcs.size(), // max nodes/arcs
+                    supplies.size(), arcs.size(), // current nodes/arcs
+                    NULL, NULL, // arcs have inf cap, zero cost (to begin)
+                    deficits, // deficits for each node
+                    starts, ends); // start/end of each edge
+
+    delete[] deficits;
+    delete[] starts;
+    delete[] ends;
+  }
+  Domain eval (const std::vector<Domain>& costs) const {
+    assert(costs.size() == _arcs.size());
+    if (costs.size() == 0)
+      return 0;
+    for (int i = 0, size = costs.size(); i < size; ++i) {
+      _solver.ChgCost(i, costs[i].template as<MCFClass::CNumber>());
+    }
+    _solver.SolveMCF();
+    if (_solver.MCFGetStatus() == MCFClass::kUnfeasible){
+      return -infinity<Domain>();
+    } else if (_solver.MCFGetFO() == MCFClass::Inf<MCFClass::FONumber>()) {
+      return infinity<Domain>();
+    } else if (_solver.MCFGetFO() == -MCFClass::Inf<MCFClass::FONumber>()) {
+      return -infinity<Domain>();        
+    } else {
+      return _solver.MCFGetFO();
+    }
+  }
+  void print(std::ostream& cout) const {
+    std::string supplyString = "[";
+    for (int i = 0, size = _supplies.size(); i < size; ++i) {
+      if (i > 0)
+        supplyString += ",";
+      std::stringstream stream;
+      stream << _supplies[i];
+      supplyString += stream.str();
+    }
+    supplyString += ']';
+
+    std::string arcString = "[";
+    for (int i = 0, size = _arcs.size(); i < size; ++i) {
+      if (i > 0)
+        arcString += ",";
+      {
+        std::stringstream stream;
+        stream << _arcs[i].first;
+        arcString += stream.str() + ":";
+      }
+      {
+        std::stringstream stream;
+        stream << _arcs[i].second;
+        arcString += stream.str();
+      }
+    }
+    arcString += ']';
+
+    cout << "MCF<" << supplyString << "," << arcString << ">";
+  }
+private:
+  std::vector<Domain> _supplies;
+  std::vector<std::pair<int,int> > _arcs;
+  mutable MCFSimplex _solver;
+};
+
 /*#include "TemplateConstraintMatrix.hpp"
 
 template<typename Domain>
diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp
index 746e5ef..4c4519a 100644
--- a/clang/include/clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp
+++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp
@@ -24,9 +24,11 @@ struct DynamicVariableAssignment : public VariableAssignment<Domain> {
   ) : _system(system),
       _strategy(strat),
       _values(system.variableCount(), value),
+      _old_values(system.variableCount(), value),
       _unstable(system.variableCount()),
       _influence(system.variableCount(),
-                 IdSet<Variable<Domain> >(system.variableCount()))
+                 IdSet<Variable<Domain> >(system.variableCount())),
+      _touched(system.variableCount())
   { }
 
   const Domain operator[](const Variable<Domain>& var) {
@@ -34,18 +36,13 @@ struct DynamicVariableAssignment : public VariableAssignment<Domain> {
     return _values[var];
   }
 
-  /*void stabilise() {
-    if (!_unstable.empty()) {
-      Variable<Domain>& var = _system.variable(*_unstable.begin());
-      solve(var);
-    }
-    }*/
-
   void invalidate(const Variable<Domain>& x) {
     if (!_unstable.contains(x)) {
       solver_log::fixpoint << indent() << "Invalidating " << x << std::endl;
  
       _unstable.insert(x);
+      _touched.insert(x);
+      _old_values[x] = _values[x];
       _values[x] = infinity<Domain>();
 
       IdSet<Variable<Domain> > infl = _influence[x];
@@ -61,6 +58,23 @@ struct DynamicVariableAssignment : public VariableAssignment<Domain> {
     }
   }
 
+  IdSet<Variable<Domain> > get_changed() {
+    IdSet<Variable<Domain> > changed;
+    for (typename IdSet<Variable<Domain> >::iterator
+           it = _touched.begin(),
+           ei = _touched.end();
+         it != ei;
+         ++it) {
+      Variable<Domain>& var = _system.variable(*it);
+      if (!_unstable.contains(var) && _old_values[var] != _values[var]) {
+        changed.insert(var);
+        _touched.remove(var);
+      }
+    }
+    //_touched.clear();
+    return changed;
+  }
+
 private:
 
   void solve(const Variable<Domain>& x) {
@@ -79,8 +93,6 @@ private:
       if (val != _values[x]) {
         solver_log::fixpoint << x << " = " << val << std::endl;
  
-        _strategy.invalidate(x);
-
         IdSet<Variable<Domain> > oldInfluence = _influence[x];
         _influence[x].clear();
         _values[x] = val;
@@ -118,8 +130,10 @@ private:
   const EquationSystem<Domain>& _system;
   DynamicMaxStrategy<Domain>& _strategy;
   IdMap<Variable<Domain>, Domain> _values;
+  IdMap<Variable<Domain>, Domain> _old_values;
   IdSet<Variable<Domain> > _unstable;
   IdMap<Variable<Domain>,IdSet<Variable<Domain> > > _influence;
+  IdSet<Variable<Domain> > _touched;
 };
 
 #endif
diff --git a/clang/lib/Analysis/Interval.cpp b/clang/lib/Analysis/Interval.cpp
index e643f5e..4b274ed 100644
--- a/clang/lib/Analysis/Interval.cpp
+++ b/clang/lib/Analysis/Interval.cpp
@@ -9,8 +9,6 @@
 #include "clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp"
 #include "clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp"
 
-#include "MCFSimplex.h"
-  
 #include "llvm/ADT/PostOrderIterator.h"
 #include "llvm/ADT/DenseMap.h"
 #include "llvm/Support/Process.h"
@@ -22,8 +20,14 @@
 #include <set>
 
 using namespace clang;
+using namespace std;
+
+template<typename T>
+T neg(const T& t) {
+  return -t;
+}
 
-std::string negate_string(const std::string& str) {
+string operator-(const string& str) {
   if (str[0] == '-')
     return str.substr(1);
   return '-' + str;
@@ -31,28 +35,28 @@ std::string negate_string(const std::string& str) {
 
 #include <sstream>
 template<typename T>
-std::string toString(const T& obj) {
-  std::stringstream stream;
+string toString(const T& obj) {
+  stringstream stream;
   stream << obj;
   return stream.str();
 }
 
 #include <ostream>
 template<typename K,typename V>
-std::ostream& operator<<(std::ostream& cout, const std::pair<K,V>& v) {
+ostream& operator<<(ostream& cout, const pair<K,V>& v) {
   cout << "(" << v.first << ", " << v.second << ")";
   return cout;
 }
 template<typename V>
-std::ostream& operator<<(std::ostream& cout, const std::pair<Variable<Complete<int64_t> >*, V>& v) {
+ostream& operator<<(ostream& cout, const pair<Variable<Complete<int64_t> >*, V>& v) {
   cout << "(" << v.first->name() << ", " << v.second << ")";
   return cout;
 }
 
 template<typename V>
-std::ostream& operator<<(std::ostream& cout, const std::vector<V>& v) {
+ostream& operator<<(ostream& cout, const vector<V>& v) {
   cout << "[";
-  for(typename std::vector<V>::const_iterator it = v.begin(), ei = v.end();
+  for(typename vector<V>::const_iterator it = v.begin(), ei = v.end();
       it != ei;
       ++it) {
     if (it != v.begin())
@@ -64,9 +68,9 @@ std::ostream& operator<<(std::ostream& cout, const std::vector<V>& v) {
 }
 
 template<typename K,typename V>
-std::ostream& operator<<(std::ostream& cout, const std::map<K,V>& v) {
+ostream& operator<<(ostream& cout, const map<K,V>& v) {
   cout << "{";
-  for (typename std::map<K,V>::const_iterator it = v.begin(), ei = v.end();
+  for (typename map<K,V>::const_iterator it = v.begin(), ei = v.end();
        it != ei;
        ++it) {
     if (it != v.begin())
@@ -91,219 +95,89 @@ template<>
 ZBar infinity() {
   return ZBar(1, true);
 }
-//typedef std::map<std::string, ZBar> Vector;
 
-struct Vector : public std::map<std::string, ZBar> {
+
+
+
+
+struct Vector : public map<string, ZBar> {
   Vector(const ZBar& val=0)
     : _val(val) { }
-  ZBar operator[](const std::string& key) const {
-    if (this->find(key) != this->end())
-      return this->find(key)->second;
+  ZBar operator[](const string& key) const {
+    const_iterator it = this->find(key);
+    if (it != this->end())
+      return it->second;
     return _val;
   }
-  ZBar& operator[](const std::string& key) {
-    if (this->find(key) != this->end())
-      return this->find(key)->second;
-    std::pair<iterator,bool> p = this->insert(std::pair<const std::string, ZBar>(key, _val));
+  ZBar& operator[](const string& key) {
+    iterator it = this->find(key);
+    if (it != this->end())
+      return it->second;
+    pair<iterator,bool> p = this->insert(pair<const string, ZBar>(key, _val));
     return p.first->second;
   }
   ZBar _val;
 };
 
-Vector negate_vector(const Vector& v) {
+Vector operator-(const Vector& v) {
   Vector result;
   for (Vector::const_iterator it = v.begin(),
          ei = v.end();
        it != ei;
        ++it) {
-    result[negate_string(it->first)] = it->second;
+    result[-it->first] = it->second;
   }
   return result;
 }
 
-typedef std::pair<Vector, ZBar> Result; // a "slice" of an equation
-
-Result negate_result(const Result& r) {
-  return Result(negate_vector(r.first), -r.second);
-}
-
-//typedef std::map<std::string, Result> LinearEquation; // one `Result` per variable
-struct LinearEquation : public std::map<std::string, Result> {
-  Result operator[](const std::string& key) const {
-    if (this->find(key) != this->end())
-      return this->find(key)->second;
-    Result r;
-    r.first[key] = 1;
-    r.second = 0;
-    return r;
-  }
-  Result& operator[](const std::string& key) {
-    if (this->find(key) != this->end())
-      return this->find(key)->second;
-    Result r;
-    r.first[key] = 1;
-    r.second = 0;
-    std::pair<iterator,bool> p = this->insert(std::pair<const std::string, Result>(key, r));
-    return p.first->second;
-  }
-};
-
-typedef Vector Condition;
-
-typedef EquationSystem<ZBar> EqnSys;
-typedef Expression<ZBar> EqnExpr;
-typedef Variable<ZBar> EqnVar;
-
-
-struct LinearOperator : public Operator<Vector> {
-  LinearOperator(const LinearEquation* result)
-    : _values(result) {}
-
-  Vector eval(const std::vector<Vector>& vector) const {
-    assert(vector.size() == 1);
-    const Vector& v = vector[0];
-    Vector result = v;
-    for (LinearEquation::const_iterator it = _values->begin(),
-           ei = _values->end();
-         it != ei;
-         ++it) {
-      ZBar subresult = 0;
-      for (Vector::const_iterator jt = it->second.first.begin(),
-             ej = it->second.first.end();
-           jt != ej;
-           ++jt) {
-        subresult += jt->second * v[jt->first];
+Vector operator+(const Vector& left, const Vector& right) {
+  Vector::const_iterator
+    left_iter = left.begin(),
+    left_end = left.end(),
+    right_iter = right.begin(),
+    right_end = right.end();
+
+  Vector result(left._val + right._val);
+  while (left_iter != left_end && right_iter != right_end) {
+    if (left_iter->first == right_iter->first) {
+      result[left_iter->first] = left_iter->second + right_iter->second;
+      left_iter++;
+      right_iter++;
+    } else {
+      if (left_iter->first < right_iter->first) {
+        result[left_iter->first] = left_iter->second;
+        left_iter++;
+      } else {
+        result[right_iter->first] = right_iter->second;
+        right_iter++;
       }
-      subresult += it->second.second;
-      result[it->first] = subresult;
     }
-    return result;
   }
+  Vector::const_iterator it = (right_iter == right_end ? left_iter : right_iter);
+  Vector::const_iterator end = (right_iter == right_end ? left_end : right_end);
+  for (; it != end; ++it)
+    result[it->first] = it->second;
 
-  void print(std::ostream& cout) const {
-    cout << "linear[" << *_values << "]";
-  }
-
-  const LinearEquation* _values;
-};
-
-
-
-template<class F, class M>
-void transform_values(const F& f, M& map) {
-  for (typename M::iterator it = map.begin(),
-         ei = map.end();
-       it != ei;
-       ++it) {
-    it->second = f(it->second);
-  }
-}
-
-template<class M, class F> 
-M merge_maps_with(const F& f, const M& left, const M& right) {
-  M result;
-  typename M::const_iterator first1 = left.begin(), last1 = left.end(),
-    first2 = right.begin(), last2 = right.end();
-  for (; first1 != last1 && first2 != last2;) {
-    if (first2->first < first1->first) {
-      result[first2->first] = first2->second;
-      ++first2;
-    } else if (first1->first == first2->first) {
-      result[first1->first] = f(first1->second, first2->second);
-      ++first1;
-      ++first2;
-    } else {
-      result[first1->first] = first1->second;
-      ++first1;
-    }
-  }
-  while (first1 != last1) {
-    result[first1->first] = first1->second;
-    ++first1;
-  }
-  while (first2 != last2) {
-    result[first2->first] = first2->second;
-    ++first2;
-  }
   return result;
 }
+Vector operator*(const ZBar& left, const Vector& right) {
+  Vector result(left * right._val);
 
-template<class T>
-T max(const T& l, const T& r) {
-  return (l < r ? l : r);
-}
-template<class T>
-T negate(const T& v) {
-  return -v;
-}
-template<class T>
-T addValues(const T& l, const T& r) {
-  return l + r;
-}
-
-Vector operator-(const Vector& vector) {
-  Vector result(-vector._val);
-  for (Vector::const_iterator it = vector.begin(),
-         ei = vector.end();
-       it != ei;
+  for (Vector::const_iterator
+         it = right.begin(),
+         end = right.end();
+       it != end;
        ++it) {
-    result[it->first] = -it->second;
+    result[it->first] = left * it->second;
   }
-  return result;
-}
-
-Vector operator+(const Vector& left, const Vector& right) {
-  return merge_maps_with(addValues<ZBar>, left, right);
-}
-
-Vector operator-(const Vector& left, const Vector& right) {
-  return merge_maps_with(addValues<ZBar>, left, -right);
-}
 
-Vector operator*(const Vector& left, const ZBar& right) {
-  Vector result;
-  for (Vector::const_iterator it = left.begin(),
-         ei = left.end();
-       it != ei;
-       ++it) {
-    result[it->first] = (it->second * right);
-  }
   return result;
 }
-Vector operator*(const ZBar& left, const Vector& right) {
+Vector operator*(const Vector& left, const ZBar& right) {
   return right * left;
 }
-bool operator<(const Vector& left, const Vector& right) {
-  bool equal = true;
-  for (Vector::const_iterator it = left.begin(),
-         ei = left.end();
-       it != ei;
-       ++it) {
-    if (it->second < right[it->first]) {
-      equal = false;
-    } else if (it->second > right[it->first]) {
-      return false;
-    }
-  }
-  for (Vector::const_iterator it = right.begin(),
-         ei = right.end();
-       it != ei;
-       ++it) {
-    if (left[it->first] < it->second) {
-      equal = false;
-    } else if (left[it->first] > it->second) {
-      return false;
-    }
-  }
-  return equal ? left._val < right._val : true;
-}
-
-template<>
-Vector infinity<Vector>() {
-  return Vector(infinity<ZBar>());
-}
 
-std::ostream& operator<<(std::ostream& cout, const Vector& v) {
+ostream& operator<<(ostream& cout, const Vector& v) {
   cout << "{";
   for (Vector::const_iterator it = v.begin(), ei = v.end();
        it != ei;
@@ -317,6 +191,22 @@ std::ostream& operator<<(std::ostream& cout, const Vector& v) {
 
 
 
+typedef pair<Vector, ZBar> Result; // a "slice" of an equation
+Result operator-(const Result& r) {
+  return Result(-r.first, -r.second);
+}
+
+typedef Vector Condition;
+
+typedef EquationSystem<ZBar> EqnSys;
+typedef Expression<ZBar> EqnExpr;
+typedef Variable<ZBar> EqnVar;
+
+
+
+
+
+
 
 
 
@@ -329,7 +219,7 @@ Result fromInteger(const IntegerLiteral* expr) {
 }
 
 Result fromDeclExpr(const DeclRefExpr* expr) {
-  Vector val;
+  Vector val(0);
   val[expr->getNameInfo().getAsString()] = 1;
   return Result(val, 0);
 }
@@ -356,14 +246,13 @@ Result fromBinary(const BinaryOperator* op) {
   case BO_Assign:
     return right;
   case BO_Sub:
-    right = negate_result(right);
+    right = -right;
     //transform_values(negate<ZBar>, right.first);
     //right.second *= -1;
   case BO_Add:
     {
       Result result;
-      result.first = merge_maps_with(addValues<ZBar>,
-                                     left.first, right.first);
+      result.first = left.first + right.first;
       result.second = left.second + right.second;
       return result;
     }
@@ -384,7 +273,7 @@ Result fromBinary(const BinaryOperator* op) {
       if (scalar >= 0) {
         return scalar * value;
       } else {
-        return -scalar * negate_result(value);
+        return scalar * -value;
       }
     }
   case BO_LT:
@@ -405,15 +294,13 @@ Result fromExpr(const Expr* stmt) {
     return fromInteger(static_cast<const IntegerLiteral*>(stmt));
   case Stmt::DeclRefExprClass:
     return fromDeclExpr(static_cast<const DeclRefExpr*>(stmt));
-  case Stmt::UnaryOperatorClass:
-    return fromUnary(static_cast<const UnaryOperator*>(stmt));
   case Stmt::BinaryOperatorClass:
     return fromBinary(static_cast<const BinaryOperator*>(stmt));
   }
   const Expr* expr = stmt->IgnoreParenCasts();
   if (stmt != expr)
     return fromExpr(expr);
-  llvm::errs() << "we shouldn't get here...\n";
+  assert(false);
   return Result();
 }
 
@@ -433,7 +320,7 @@ Condition fromComparison(const BinaryOperator* op, bool negate) {
     const Expr* right = op->getRHS()->IgnoreParenCasts();
 
     bool flip = false;
-    std::string name;
+    string name;
     int64_t value;
     if (left->getStmtClass() == Stmt::DeclRefExprClass) {
       name = static_cast<const DeclRefExpr*>(left)->getNameInfo().getAsString();
@@ -465,13 +352,13 @@ Condition fromComparison(const BinaryOperator* op, bool negate) {
     switch (operation) {
     case BO_LT:
       if (negate)
-	cond[negate_string(name)] = -value;
+	cond[-name] = -value;
       else
 	cond[name] = value - 1;
       break;
     case BO_LE:
       if (negate)
-	cond[negate_string(name)] = -(value + 1);
+	cond[-name] = -(value + 1);
       else
 	cond[name] = value;
       break;
@@ -479,13 +366,13 @@ Condition fromComparison(const BinaryOperator* op, bool negate) {
       if (negate)
 	cond[name] = value - 1;
       else
-	cond[negate_string(name)] = -value;
+	cond[-name] = -value;
       break;
     case BO_GT:
       if (negate)
 	cond[name] = value;
       else
-	cond['-' + name] = -(value + 1);
+	cond[-name] = -(value + 1);
       break;
     }
   }
@@ -493,97 +380,20 @@ Condition fromComparison(const BinaryOperator* op, bool negate) {
 }
 
 
-struct MinCostFlow : public Operator<ZBar> {
-  MinCostFlow(const std::vector<std::pair<EqnVar*, ZBar> >& supplies)
-    : solver(0, 0) {
-    assert(supplies.size() % 2 == 0);
-
-    if (supplies.size() == 0)
-      return;
-
-    MCFClass::FNumber* deficits = new MCFClass::FNumber[supplies.size()];
-    MCFClass::Index* starts = new MCFClass::Index[supplies.size()];
-    MCFClass::Index* ends = new MCFClass::Index[supplies.size()];
-
-    for (int i = 0, size = supplies.size(); i < size; i += 2) {
-      deficits[i] = -supplies[i].second.as<MCFClass::FNumber>();
-      deficits[i+1] = -supplies[i+1].second.as<MCFClass::FNumber>();
-
-      starts[i] = i+1;
-      ends[i] = i+2;
-
-      starts[i+1] = i+2;
-      ends[i+1] = i+1;
-
-      printableSupply.push_back(supplies[i].second);
-      printableSupply.push_back(supplies[i+1].second);
-    }
-
-    solver.LoadNet(supplies.size(), supplies.size(), // max nodes/arcs
-                   supplies.size(), supplies.size(), // current nodes/arcs
-                   NULL, NULL, // arcs have inf cap, zero cost (to begin)
-                   deficits, // deficits for each node
-                   starts, ends); // start/end of each edge
-
-    delete[] deficits;
-    delete[] starts;
-    delete[] ends;
-  }
-
-  ZBar eval(const std::vector<ZBar>& args) const {
-    if (args.size() == 0)
-      return 0;
-    for (int i = 0, size = args.size(); i < size; ++i) {
-      //if (args[i] == infinity<ZBar>()) {
-      //  if (!solver.IsClosedArc(i))
-      //      solver.CloseArc(i);
-      //} else {
-      //  if (solver.IsClosedArc(i))
-      //    solver.OpenArc(i);
-      solver.ChgCost(i, args[i].as<MCFClass::CNumber>());
-      //}
-    }
-    std::cerr << "MCF" << this << "<" << printableSupply << ">(" << args << ")" << " = ";
-    solver.SolveMCF();
-    if (solver.MCFGetStatus() == MCFClass::kUnfeasible){
-      std::cerr << -infinity<ZBar>() << std::endl;
-      return -infinity<ZBar>();
-    } else {
-      if (solver.MCFGetFO() == MCFClass::Inf<MCFClass::FONumber>()) {
-        std::cerr << infinity<ZBar>() << std::endl;
-        return infinity<ZBar>();
-      } else if (solver.MCFGetFO() == -MCFClass::Inf<MCFClass::FONumber>()) {
-        std::cerr << -infinity<ZBar>() << std::endl;
-        return -infinity<ZBar>();        
-      } else {
-        ZBar result = solver.MCFGetFO();
-        std::cerr << result << std::endl;
-        return result;
-      }
-    }
-  }
-
-  virtual void print(std::ostream& cout) const {
-    cout << "MCF" << this << "<" << printableSupply << ">";
-  }
-  
-  mutable MCFSimplex solver;
-  std::vector<ZBar> printableSupply;
-};
 
 
 /* Blocks */
 
-typedef std::map<std::string, unsigned int> Counters;
-typedef std::map<std::string, EqnVar*> VarMap;
-typedef std::map<const CFGBlock*, std::set<std::string> > BlockVars;
+typedef map<string, unsigned int> Counters;
+typedef map<string, EqnVar*> VarMap;
+typedef map<const CFGBlock*, set<string> > BlockVars;
 
 void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
   Counters counters;
-  std::string block_id = toString(block->getBlockID());
+  string block_id = toString(block->getBlockID());
   VarMap vars;
 
-  for (std::set<std::string>::iterator it = block_vars[block].begin(),
+  for (set<string>::iterator it = block_vars[block].begin(),
 	 ei = block_vars[block].end();
        it != ei;
        ++it) {
@@ -597,7 +407,7 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
     const CFGStmt* cfg_stmt = it->getAs<CFGStmt>();
     const Stmt* stmt = cfg_stmt->getStmt();
 
-    std::string name = "";
+    string name = "";
     Result result;
     if (stmt->getStmtClass() == Stmt::BinaryOperatorClass) {
       const BinaryOperator* binop = static_cast<const BinaryOperator*>(stmt);
@@ -630,41 +440,55 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
     if (name == "")
       continue;
 
-    std::string count = toString(counters[name]);
+    string count = toString(counters[name]);
     EqnVar* var = &system.variable(name + '-' + block_id + '[' + count + ']');
-    EqnVar* negative_var = &system.variable(negate_string(name) + '-' + block_id + '[' + count + ']');
+    EqnVar* negative_var = &system.variable(-name + '-' + block_id + '[' + count + ']');
     counters[name]++;
     for (int negative = 0; negative < 2; ++negative) { // one loop for positive, the other for negative
       if (negative) {
-	result = negate_result(result);
+	result = -result;
       }
       
       EqnExpr* expression;
 
       if (result.first.size() > 0) {
-        // set up the min-cost stuff
-        std::vector<EqnExpr*> minCostArgs;
-        std::vector<std::pair<EqnVar*, ZBar> > varCosts;
-        for (Vector::iterator it = result.first.begin(),
+        // make sure all our variables exist in vars
+        for (Vector::iterator
+               it = result.first.begin(),
                ei = result.first.end();
              it != ei;
              ++it) {
           if (!vars[it->first])
             vars[it->first] = &system.variable(it->first + '-' + block_id + "-pre");
+        }
 
-          EqnVar& var = *vars[it->first];
-          EqnVar& negVar = *vars[negate_string(it->first)];
-
-          varCosts.push_back(std::pair<EqnVar*, ZBar>(&var, it->second));
-          varCosts.push_back(std::pair<EqnVar*, ZBar>(&negVar, -it->second));
-
-          minCostArgs.push_back(&var);
-          minCostArgs.push_back(&negVar);
+        // set up the min-cost-flow operator
+        vector<ZBar> supplies;
+        vector<pair<int,int> > arcs;
+        vector<EqnExpr*> minCostArgs;
+        ZBar dummy_value = 0;
+        supplies.push_back(dummy_value); // dummy node to suck up flow
+        int index = 1; // the solver uses 1-indexing, for some reason
+        for (map<std::string,EqnVar*>::iterator
+               it = vars.begin(),
+               ei = vars.end();
+             it != ei;
+             it++) {
+          index++;
+          supplies.push_back(result.first[it->first]);
+          dummy_value -= result.first[it->first];
+          if (it->first[0] == '-')
+            arcs.push_back(pair<int,int>(1,index));
+          else 
+            arcs.push_back(pair<int,int>(index,1));
+          minCostArgs.push_back(vars[it->first]);
         }
-        EqnExpr* minCostExpr = &system.expression(new MinCostFlow(varCosts), minCostArgs);
+        supplies[0] = dummy_value;
+
+        EqnExpr* minCostExpr = &system.expression(new MinCostFlow<ZBar>(supplies, arcs), minCostArgs);
         
         // add the constant factor to the min-cost bit
-        std::vector<EqnExpr*> additionArgs;
+        vector<EqnExpr*> additionArgs;
         additionArgs.push_back(&system.constant(result.second));
         additionArgs.push_back(minCostExpr);
         expression = &system.expression(new Addition<ZBar>(), additionArgs);
@@ -673,7 +497,7 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
       }
 
       // max(-inf, expr), so strategy iteration will work
-      std::vector<EqnExpr*> maxArgs;
+      vector<EqnExpr*> maxArgs;
       maxArgs.push_back(&system.constant(-infinity<ZBar>()));
       maxArgs.push_back(expression);
       if (negative)
@@ -682,9 +506,9 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
 	system[*var] = &system.maxExpression(maxArgs);
     }
     vars[name] = var;
-    vars[negate_string(name)] = negative_var;
+    vars[-name] = negative_var;
     block_vars[block].insert(name);
-    block_vars[block].insert(negate_string(name));
+    block_vars[block].insert(-name);
   }
 
   // add to our successor entry values
@@ -704,7 +528,7 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
       ZBar val = cond[jt->first];
       EqnVar* var = &system.variable(jt->first + '-' + toString((*it)->getBlockID()) + "-pre");
       if (system[*var] == NULL) {
-	std::vector<EqnExpr*> maxArgs;
+	vector<EqnExpr*> maxArgs;
 	maxArgs.push_back(&system.constant(-infinity<ZBar>()));
 	system[*var] = &system.maxExpression(maxArgs);
       }
@@ -718,51 +542,13 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
 	expr = jt->second;
       } else {
 	// need a min here
-	std::vector<EqnExpr*> minArgs;
+	vector<EqnExpr*> minArgs;
 	minArgs.push_back(&system.constant(val));
 	minArgs.push_back(jt->second);
 	expr = &system.expression(new Minimum<ZBar>(), minArgs);
       }
 
       system[*var]->arguments().push_back(expr);
-
-      /*if (expr) {
-	std::set<std::string> ignore;
-	for (VarMap::iterator
-	       variables = vars.begin(),
-	       variables_end = vars.end();
-	     variables != variables_end;
-	     ++variables) {
-	  if (ignore.find(variables->first) != ignore.end())
-	    continue;
-	  ignore.insert(negate_string(variables->first));
-
-	  std::vector<EqnExpr*> plusArgs;
-	  for (int negate = 0; negate < 2; ++negate) {
-	    std::string var_name = negate ? negate_string(variables->first) : variables->first;
-            if (cond[var_name] == infinity<ZBar>()) {
-              // min(x, inf) = x
-              plusArgs.push_back(vars[var_name]);
-            } else if (cond[var_name] == -infinity<ZBar>()) {
-              // min(x, -inf) = -inf
-              plusArgs.push_back(&system.constant(-infinity<ZBar>()));
-            } else {
-              // min(x, y) = ?
-              std::vector<EqnExpr*> minArgs;
-              minArgs.push_back(vars[var_name]);
-              minArgs.push_back(&system.constant(cond[var_name]));
-              plusArgs.push_back(&system.expression(new Minimum<ZBar>(), minArgs));
-            }
-	  }
-
-	  std::vector<EqnExpr*> guard_args;
-	  guard_args.push_back(&system.expression(new Addition<ZBar>(), plusArgs)); // value
-	  guard_args.push_back(&system.constant(0)); // lower-bound (so value must be >= this)
-	  guard_args.push_back(expr); // result
-	  expr = &system.expression(new Guard<ZBar>(), guard_args);
-	}
-	system[*var]->arguments().push_back(expr);
-        }*/
     }
   }
 }
@@ -788,27 +574,27 @@ void IntervalAnalysis::runOnAllBlocks(const Decl& decl) {
   EqnSys system; 
   BlockVars block_vars;
 
-  std::vector<EqnExpr*> infArg; 
+  vector<EqnExpr*> infArg; 
   infArg.push_back(&system.constant(-infinity<ZBar>())); // left-most argument has to be -infinity
   infArg.push_back(&system.constant(infinity<ZBar>()));
-  std::set<std::string>& function_arguments = block_vars[&cfg->getEntry()];
-  std::string block_id = toString(cfg->getEntry().getBlockID());
+  set<string>& function_arguments = block_vars[&cfg->getEntry()];
+  string block_id = toString(cfg->getEntry().getBlockID());
   if (const FunctionDecl* func = dyn_cast<const FunctionDecl>(&decl)) {
     for (unsigned int i = func->getNumParams(); i > 0; i--) {
-      std::string name = func->getParamDecl(i-1)->getNameAsString();
+      string name = func->getParamDecl(i-1)->getNameAsString();
 
       // add the variables to the first block
       function_arguments.insert(name);
-      function_arguments.insert(negate_string(name));
+      function_arguments.insert(neg(name));
       
       // set the vars to infinity (unconstrained)
       system[system.variable(name + '-' + block_id + "-pre")] = &system.maxExpression(infArg);
-      system[system.variable(negate_string(name) + '-' + block_id + "-pre")] = &system.maxExpression(infArg);
+      system[system.variable(neg(name) + '-' + block_id + "-pre")] = &system.maxExpression(infArg);
     }
   }
 
-  std::set<const CFGBlock*> seen;
-  std::deque<const CFGBlock*> todo;
+  set<const CFGBlock*> seen;
+  deque<const CFGBlock*> todo;
   todo.push_back(&cfg->getEntry());
 
   while (!todo.empty()) {
@@ -835,7 +621,7 @@ void IntervalAnalysis::runOnAllBlocks(const Decl& decl) {
   Solver<ZBar> solver(system);
 
   for (unsigned int i = 0, size = system.variableCount(); i < size; ++i) {
-    EqnVar& var = system.variable(size - i - 1);
+    EqnVar& var = system.variable(i);
     llvm::errs() << toString(var.name()) << " = " << toString(solver.solve(var)) << "\n";
   }
 }
diff --git a/clang/lib/Analysis/MCFClass.h b/clang/lib/Analysis/MCFClass.h
deleted file mode 100755
index 41adb82..0000000
--- a/clang/lib/Analysis/MCFClass.h
+++ /dev/null
@@ -1,2436 +0,0 @@
-/*--------------------------------------------------------------------------*/
-/*-------------------------- File MCFClass.h -------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @file
- * Header file for the abstract base class MCFClass, which defines a standard
- * interface for (linear or convex quadratic separable) Min Cost Flow Problem
- * solvers, to be implemented as derived classes.
- *
- * \version 3.01
- *
- * \date 30 - 09 - 2011
- *
- * \author Alessandro Bertolini \n
- *         Operations Research Group \n
- *         Dipartimento di Informatica \n
- *         Universita' di Pisa \n
- *
- * \author Antonio Frangioni \n
- *         Operations Research Group \n
- *         Dipartimento di Informatica \n
- *         Universita' di Pisa \n
- *
- * \author Claudio Gentile \n
- *         Istituto di Analisi di Sistemi e Informatica \n
- *         Consiglio Nazionale delle Ricerche \n
- *
- * Copyright &copy 1996 - 2011 by Antonio Frangioni, Claudio Gentile
- */
-
-/*--------------------------------------------------------------------------*/
-/*--------------------------------------------------------------------------*/
-/*----------------------------- DEFINITIONS --------------------------------*/
-/*--------------------------------------------------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#ifndef __MCFClass
- #define __MCFClass  /* self-identification: #endif at the end of the file */
-
-/*--------------------------------------------------------------------------*/
-/*--------------------------------- MACROS ---------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @defgroup MCFCLASS_MACROS Compile-time switches in MCFClass.h
-    These macros control some important details of the class interface.
-    Although using macros for activating features of the interface is not
-    very C++, switching off some unused features may allow some
-    implementation to be more efficient in running time or memory.
-    @{ */
-
-/*-------------------------------- USENAME0 --------------------------------*/
-
-#define USENAME0 1
-
-/**< Decides if 0 or 1 is the "name" of the first node.
-   If USENAME0 == 1, (warning: it has to be *exactly* 1), then the node
-   names go from 0 to n - 1, otherwise from 1 to n. Note that this does not
-   affect the position of the deficit in the deficit vectors, i.e., the
-   deficit of the i-th node - be its "name" `i' or `i - 1' - is always in
-   the i-th position of the vector. */
-
-/*@}  end( group( MCFCLASS_MACROS ) ) */ 
-/*--------------------------------------------------------------------------*/
-/*------------------------------ INCLUDES ----------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#include "OPTUtils.h"
-
-/* OPTtypes.h defines standard interfaces for timing and random routines, as
-   well as the namespace OPTtypes_di_unipi_it and the macro
-   OPT_USE_NAMESPACES, useful for switching off all namespaces in one blow
-   for those strange cases where they create problems. */
-
-#include <iomanip>
-#include <sstream>
-#include <limits>
-#include <cassert>
-
-// QUICK HACK
-#define throw(x) assert(false)
-
-/*--------------------------------------------------------------------------*/
-/*------------------------- NAMESPACE and USING ----------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#if( OPT_USE_NAMESPACES )
-namespace MCFClass_di_unipi_it
-{
- /** @namespace MCFClass_di_unipi_it
-     The namespace MCFClass_di_unipi_it is defined to hold the MCFClass
-     class and all the relative stuff. It comprises the namespace
-     OPTtypes_di_unipi_it. */
-
- using namespace OPTtypes_di_unipi_it;
-#endif
-
-/*@}  end( group( MCFCLASS_CONSTANTS ) ) */
-/*--------------------------------------------------------------------------*/
-/*-------------------------- CLASS MCFClass --------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @defgroup MCFCLASS_CLASSES Classes in MCFClass.h
-    @{ */
-
-/*--------------------------------------------------------------------------*/
-/*--------------------------- GENERAL NOTES --------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** This abstract base class defines a standard interface for (linear or
-    convex quadartic separable) Min Cost Flow (MCF) problem solvers.
-
-    The data of the problem consist of a (directed) graph G = ( N , A ) with
-    n = |N| nodes and m = |A| (directed) arcs. Each node `i' has a deficit
-    b[ i ], i.e., the amount of flow that is produced/consumed by the node:
-    source nodes (which produce flow) have negative deficits and sink nodes
-    (which consume flow) have positive deficits. Each arc `(i, j)' has an
-    upper capacity U[ i , j ], a linear cost coefficient C[ i , j ] and a
-    (non negative) quadratic cost coefficient Q[ i , j ]. Flow variables
-    X[ i , j ] represents the amount of flow  to be sent on arc (i, j).
-    Parallel arcs, i.e., multiple copies of the same arc `(i, j)' (with
-    possibily different costs and/or capacities) are in general allowed.
-    The formulation of the problem is therefore:
-    \f[
-     \min \sum_{ (i, j) \in A } C[ i , j ] X[ i, j ] +
-                                Q[ i , j ] X[ i, j ]^2 / 2
-    \f]
-    \f[
-     (1) \sum_{ (j, i) \in A } X[ j , i ] -
-         \sum_{ (i, j) \in A } X[ i , j ] = b[ i ]
-         \hspace{1cm} i \in N
-    \f]
-    \f[
-     (2) 0 \leq X[ i , j ] \leq U[ i , j ]
-         \hspace{1cm} (i, j) \in A
-    \f]
-    The n equations (1) are the flow conservation constraints and the 2m
-    inequalities (2) are the flow nonnegativity and capacity constraints.
-    At least one of the flow conservation constraints is redundant, as the
-    demands must be balanced (\f$\sum_{ i \in N } b[ i ] = 0\f$); indeed,
-    exactly n - ConnectedComponents( G ) flow conservation constraints are
-    redundant, as demands must be balanced in each connected component of G.
-    Let us denote by QA and LA the disjoint subsets of A containing,
-    respectively, "quadratic" arcs (with Q[ i , j ] > 0) and "linear" arcs
-    (with Q[ i , j ] = 0); the (MCF) problem is linear if QA is empty, and
-    nonlinear (convex quadratic) if QA is nonempty.
-
-    The dual of the problem is:
-    \f[
-     \max \sum_{ i \in N } Pi[ i ] b[ i ] -
-          \sum_{ (i, j) \in A } W[ i , j ] U[ i , j ] -
-          \sum_{ (i, j) \in AQ } V[ i , j ]^2 / ( 2 * Q[ i , j ] )
-    \f]
-    \f[
-     (3.a) C[ i , j ] - Pi[ j ] + Pi[ i ] + W[ i , j ] - Z[ i , j ] = 0
-           \hspace{1cm} (i, j) \in AL
-    \f]
-    \f[
-     (3.b) C[ i , j ] - Pi[ j ] + Pi[ i ] + W[ i , j ] - Z[ i , j ] =
-           V[ i , j ]
-           \hspace{1cm} (i, j) \in AQ
-    \f]
-    \f[
-     (4.a) W[ i , j ] \geq 0 \hspace{1cm} (i, j) \in A
-    \f]
-    \f[
-     (4.b) Z[ i , j ] \geq 0 \hspace{1cm} (i, j) \in A
-    \f]
-
-    Pi[] is said the vector of node potentials for the problem, W[] are
-    bound variables and Z[] are slack variables. Given Pi[], the quantities
-    \f[
-     RC[ i , j ] =  C[ i , j ] + Q[ i , j ] * X[ i , j ] - Pi[ j ] + Pi[ i ]
-    \f]
-    are said the "reduced costs" of arcs.
-
-    A primal and dual feasible solution pair is optimal if and only if the
-    complementary slackness conditions
-    \f[
-     RC[ i , j ] > 0 \Rightarrow X[ i , j ] = 0
-    \f]
-    \f[
-     RC[ i , j ] < 0 \Rightarrow X[ i , j ] = U[ i , j ]
-    \f]
-    are satisfied for all arcs (i, j) of A.
-
-    The MCFClass class provides an interface with methods for managing and
-    solving problems of this kind. Actually, the class can also be used as
-    an interface for more general NonLinear MCF problems, where the cost
-    function either nonseparable ( C( X ) ) or arc-separable
-    ( \f$\sum_{ (i, j) \in A } C_{i,j}( X[ i, j ] )\f$ ). However, solvers
-    for NonLinear MCF problems are typically objective-function-specific,
-    and there is no standard way for inputting a nonlinear function different
-    from a separable convex quadratic one, so only the simplest form is dealt
-    with in the interface, leaving more complex NonLinear parts to the
-    interface of derived classes. */
-
-class MCFClass {
-
-/*--------------------------------------------------------------------------*/
-/*----------------------- PUBLIC PART OF THE CLASS -------------------------*/
-/*--------------------------------------------------------------------------*/
-/*--                                                                      --*/
-/*--  The following methods and data are the actual interface of the      --*/
-/*--  class: the standard user should use these methods and data only.    --*/
-/*--                                                                      --*/
-/*--------------------------------------------------------------------------*/
-
- public:
-
-/*--------------------------------------------------------------------------*/
-/*---------------------------- PUBLIC TYPES --------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @name Public types
-    The MCFClass defines four main public types:
-
-    - Index, the type of arc and node indices;
-
-    - FNumber, the type of flow variables, arc capacities, and node deficits;
-
-    - CNumber, the type of flow costs, node potentials, and arc reduced costs;
-
-    - FONumber, the type of objective function value.
-
-    By re-defining the types in this section, most MCFSolver should be made
-    to work with any reasonable choice of data type (= one that is capable of
-    properly representing the data of the instances to be solved). This may
-    be relevant due to an important property of MCF problems: *if all arc
-    capacities and node deficits are integer, then there exists an integral
-    optimal primal solution*, and  *if all arc costs are integer, then there
-    exists an integral optimal dual solution*. Even more importantly, *many
-    solution algorithms will in fact produce an integral primal/dual
-    solution for free*, because *every primal/dual solution they generate
-    during the solution process is naturally integral*. Therefore, one can
-    use integer data types to represent everything connected with flows and/or
-    costs if the corresponding data is integer in all instances one needs to
-    solve. This directly translates in significant memory savings and/or speed
-    improvements.
-
-    *It is the user's responsibility to ensure that these types are set to
-    reasonable values*. So, the experienced user may want to experiment with
-    setting this data properly if memory footprint and/or speed is a primary
-    concern. Note, however, that *not all solution algorithms will happily
-    accept integer data*; one example are Interior-Point approaches, which
-    require both flow and cost variables to be continuous (float). So, the
-    viability of setting integer data (as well as its impact on performances)
-    is strictly related to the specific kind of algorithm used. Since these
-    types are common to all derived classes, they have to be set taking into
-    account the needs of all the solvers that are going to be used, and
-    adapting to the "worst case"; of course, FNumber == CNumber == double is
-    going to always be an acceptable "worst case" setting. MCFClass may in a
-    future be defined as a template class, with these as template parameters,
-    but this is currently deemed overkill and avoided.
-
-    Finally, note that the above integrality property only holds for *linear*
-    MCF problems. If any arc has a nonzero quadratic cost coefficient, optimal
-    flows and potentials may be fractional even if all the data of the problem
-    (comprised quadratic cost coefficients) is integer. Hence, for *quadratic*
-    MCF solvers, a setting like FNumber == CNumber == double is actually
-    *mandatory*, for any reasonable algorithm will typically misbehave
-    otherwise.
-    @{ */
-
-/*--------------------------------------------------------------------------*/
-
- typedef unsigned int    Index;           ///< index of a node or arc ( >= 0 )
- typedef Index          *Index_Set;       ///< set (array) of indices
- typedef const Index    cIndex;           ///< a read-only index
- typedef cIndex        *cIndex_Set;       ///< read-only index array
-
-/*--------------------------------------------------------------------------*/
-
- typedef int             SIndex;           ///< index of a node or arc 
- typedef SIndex         *SIndex_Set;       ///< set (array) of indices
- typedef const SIndex   cSIndex;           ///< a read-only index
- typedef cSIndex       *cSIndex_Set;       ///< read-only index array
-
-/*--------------------------------------------------------------------------*/
-
- typedef double          FNumber;        ///< type of arc flow
- typedef FNumber        *FRow;           ///< vector of flows
- typedef const FNumber  cFNumber;        ///< a read-only flow
- typedef cFNumber      *cFRow;           ///< read-only flow array
-
-/*--------------------------------------------------------------------------*/
-
- typedef double          CNumber;        ///< type of arc flow cost
- typedef CNumber        *CRow;           ///< vector of costs
- typedef const CNumber  cCNumber;        ///< a read-only cost
- typedef cCNumber      *cCRow;           ///< read-only cost array
-
-/*--------------------------------------------------------------------------*/
-
- typedef double          FONumber; 
- /**< type of the objective function: has to hold sums of products of
-    FNumber(s) by CNumber(s) */
-
- typedef const FONumber cFONumber;       ///< a read-only o.f. value
-
-/*--------------------------------------------------------------------------*/
-/** Very small class to simplify extracting the "+ infinity" value for a
-    basic type (FNumber, CNumber, Index); just use Inf<type>(). */
-
-   template <typename T>
-   class Inf {
-    public:
-     Inf() {}
-     operator T() { return( std::numeric_limits<T>::max() ); }
-    };
-
-/*--------------------------------------------------------------------------*/
-/** Very small class to simplify extracting the "machine epsilon" for a
-    basic type (FNumber, CNumber); just use Eps<type>(). */
-
-   template <typename T>
-   class Eps {
-    public:
-     Eps() {}
-     operator T() { return( std::numeric_limits<T>::epsilon() ); }
-    };
-
-/*--------------------------------------------------------------------------*/
-/** Small class for exceptions. Derives from std::exception implementing the
-    virtual method what() -- and since what is virtual, remember to always
-    catch it by reference (catch exception &e) if you want the thing to work.
-    MCFException class are thought to be of the "fatal" type, i.e., problems
-    for which no solutions exists apart from aborting the program. Other kinds
-    of exceptions can be properly handled by defining derived classes with
-    more information. */
-
- class MCFException : public exception {
- public:
-  MCFException( const char *const msg = 0 ) { errmsg = msg; }
-
-  // wow, such a hack
-#undef throw
-  const char* what( void ) const throw () { return( errmsg ); }
-#define throw(x) assert(false)
- private:
-  const char *errmsg;
-  };
-
-/*--------------------------------------------------------------------------*/
-/** Public enum describing the possible parameters of the MCF solver, to be
-    used with the methods SetPar() and GetPar(). */
-
-  enum MCFParam { kMaxTime = 0 ,     ///< max time 
-                  kMaxIter ,         ///< max number of iteration
-                  kEpsFlw ,          ///< tolerance for flows
-                  kEpsDfct ,         ///< tolerance for deficits
-                  kEpsCst ,          ///< tolerance for costs
-                  kReopt ,           ///< whether or not to reoptimize
-                  kLastParam         /**< dummy parameter: this is used to
-                                        allow derived classes to "extend"
-                                        the set of parameters. */
-                  };
-
-/*--------------------------------------------------------------------------*/
-/** Public enum describing the possible status of the MCF solver. */
-
-  enum MCFStatus { kUnSolved = -1 , ///< no solution available
-                   kOK = 0 ,        ///< optimal solution found
-
-                   kStopped ,       ///< optimization stopped
-                   kUnfeasible ,    ///< problem is unfeasible
-                   kUnbounded ,     ///< problem is unbounded
-                   kError           ///< error in the solver
-                   };
-
-/*--------------------------------------------------------------------------*/
-/** Public enum describing the possible reoptimization status of the MCF
-    solver. */
-
-  enum MCFAnswer { kNo = 0 ,  ///< no 
-                   kYes       ///< yes
-                   };
-
-/*--------------------------------------------------------------------------*/
-/** Public enum describing the possible file formats in WriteMCF(). */
-
-  enum MCFFlFrmt { kDimacs = 0 ,    ///< DIMACS file format for MCF
-                   kQDimacs ,       ///< quadratic DIMACS file format for MCF
-                   kMPS ,           ///< MPS file format for LP
-                   kFWMPS           ///< "Fixed Width" MPS format
-                   };
-
-/*--------------------------------------------------------------------------*/
-/** Base class for representing the internal state of the MCF algorithm. */
-
- class MCFState {
- public:
-   MCFState( void ) {};
-   virtual ~MCFState() {};
- };
-
- typedef MCFState *MCFStatePtr;  ///< pointer to a MCFState
-
-/*@} -----------------------------------------------------------------------*/
-/*--------------------------- PUBLIC METHODS -------------------------------*/
-/*--------------------------------------------------------------------------*/
-/*---------------------------- CONSTRUCTOR ---------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @name Constructors
-    @{ */
-
-   MCFClass( cIndex nmx = 0 , cIndex mmx = 0 )
-   {
-    nmax = nmx;
-    mmax = mmx;
-    n = m = 0;
-
-    status = kUnSolved;
-    Senstv = true;
-
-    EpsFlw = Eps<FNumber>() * 100;
-    EpsCst = Eps<CNumber>() * 100;
-    EpsDfct = EpsFlw * ( nmax ? nmax : 100 );
-
-    MaxTime = 0;
-    MaxIter = 0;
-
-    MCFt = NULL;
-    }
-
-/**< Constructor of the class.
-
-   nmx and mmx, if provided, are taken to be respectively the maximum number 
-   of nodes and arcs in the network. If nonzero values are passed, memory
-   allocation can be anticipated in the constructor, which is sometimes
-   desirable. The maximum values are stored in the protected fields nmax and
-   mmax, and can be changed with LoadNet() [see below]; however, changing
-   them typically requires memory allocation/deallocation, which is
-   sometimes undesirable outside the constructor.
-
-   After that an object has been constructed, no problem is loaded; this has
-   to be done with LoadNet() [see below]. Thus, it is an error to invoke any
-   method which requires the presence of a problem (typicall all except those
-   in the initializations part). The base class provides two protected fields
-   n and m for the current number of nodes and arcs, respectively, that are
-   set to 0 in the constructor precisely to indicate that no instance is
-   currently loaded. */
-
-/*@} -----------------------------------------------------------------------*/
-/*-------------------------- OTHER INITIALIZATIONS -------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @name Other initializations
-    @{ */
-
-   virtual void LoadNet( cIndex nmx = 0 , cIndex mmx = 0 , cIndex pn = 0 ,
-                         cIndex pm = 0 , cFRow pU = NULL , cCRow pC = NULL ,
-                         cFRow pDfct = NULL , cIndex_Set pSn = NULL ,
-                         cIndex_Set pEn = NULL ) = 0;
-
-/**< Inputs a new network.
- 
-   The parameters nmx and mmx are the new max number of nodes and arcs,
-   possibly overriding those set in the constructor [see above], altough at
-   the likely cost of memory allocation and deallocation. Passing nmx == 
-   mmx == 0 is intended as a signal to the solver to deallocate everything
-   and wait for new orders; in this case, all the other parameters are ignored.
-
-   Otherwise, in principle all the other parameters have to be provided.
-   Actually, some of them may not be needed for special classes of MCF
-   problems (e.g., costs in a MaxFlow problem, or start/end nodes in a
-   problem defined over a graph with fixed topology, such as a complete
-   graph). Also, passing NULL is allowed to set default values.
-
-   The meaning of the parameters is the following:
-
-   - pn     is the current number of nodes of the network (<= nmax).
-   - pm     is the number of arcs of the network (<= mmax).
-
-   - pU     is the m-vector of the arc upper capacities; capacities must be
-            nonnegative, but can in principle be infinite (== F_INF); passing
-            pU == NULL means that all capacities are infinite;
-
-   - pC     is the m-vector of the arc costs; costs must be finite (< C_INF);
-            passing pC == NULL means that all costs must be 0.
-
-   - pDfct  is the n-vector of the node deficits; source nodes have negative
-            deficits and sink nodes have positive deficits; passing pDfct ==
-            NULL means that all deficits must be 0 (a circulation problem);
-
-   - pSn    is the m-vector of the arc starting nodes; pSn == NULL is in
-            principle not allowed, unless the topology of the graph is fixed;
-   - pEn    is the m-vector of the arc ending nodes; same comments as for pSn.
-
-   Note that node "names" in the arrays pSn and pEn must go from 1 to pn if
-   the macro USANAME0 [see above] is set to 0, while they must go from 0 to
-   pn - 1 if USANAME0 is set to 1. In both cases, however, the deficit of the
-   first node is read from the first (0-th) position of pDfct, that is if
-   USANAME0 == 0 then the deficit of the node with name `i' is read from
-   pDfct[ i - 1 ].
-
-   The data passed to LoadNet() can be used to specify that the arc `i' must
-   not "exist" in the problem. This is done by passing pC[ i ] == C_INF;
-   solvers which don't read costs are forced to read them in order to check
-   this, unless they provide alternative solver-specific ways to accomplish
-   the same tasks. These arcs are "closed", as for the effect of CloseArc()
-   [see below]. "invalid" costs (== C_INF) are set to 0 in order to being
-   subsequently capable of "opening" them back with OpenArc()  [see below].
-   The way in which these non-existent arcs are phisically dealt with is
-   solver-specific; in some solvers, for instance, this could be obtained by
-   simply putting their capacity to zero. Details about these issues should
-   be found in the interface of derived classes.
-
-   Note that the quadratic part of the objective function, if any, is not
-   dealt with in LoadNet(); it can only be separately provided with
-   ChgQCoef() [see below]. By default, the problem is linear, i.e., all
-   coefficients of the second-order terms in the objective function are
-   assumed to be zero. */
-
-/*--------------------------------------------------------------------------*/
-
-   virtual inline void LoadDMX( istream &DMXs , bool IsQuad = false );
-
-/**< Read a MCF instance in DIMACS standard format from the istream. The
-   format is the following. The first line must be
-
-      p min <number of nodes> <number of arcs>
-
-   Then the node definition lines must be found, in the form
-
-      n <node number> <node supply>
-
-   Not all nodes need have a node definition line; these are given zero
-   supply, i.e., they are transhipment nodes (supplies are the inverse of
-   deficits, i.e., a node with positive supply is a source node). Finally,
-   the arc definition lines must be found, in the form
-
-      a <start node> <end node> <lower bound> <upper bound> <flow cost>
-
-   There must be exactly <number of arcs> arc definition lines in the file.
-
-   This method is *not* pure virtual because an implementation is provided by
-   the base class, using the LoadNet() method (which *is* pure virtual).
-   However, the method *is* virtual to allow derived classes to implement
-   more efficient versions, should they have any reason to do so.
-
-   \note Actually, the file format accepted by LoadDMX (at least in the
-         base class implementation) is more general than the DIMACS standard
-         format, in that it is allowed to mix node and arc definitions in
-         any order, while the DIMACS file requires all node information to
-         appear before all arc information.
-
-   \note Other than for the above, this method is assumed to allow for
-         *quadratic* Dimacs files, encoding for convex quadratic separable
-         Min Cost Flow instances. This is a simple extension where each arc
-         descriptor has a sixth field, <quadratic cost>. The provided
-         istream is assumed to be quadratic Dimacs file if IsQuad is true,
-         and a regular linear Dimacs file otherwise. */
-
-/*--------------------------------------------------------------------------*/
-
-   virtual void PreProcess( void ) {}
-
-/**< Extract a smaller/easier equivalent MCF problem. The data of the instance
-   is changed and the easier one is solved instead of the original one. In the
-   MCF case, preprocessing may involve reducing bounds, identifying
-   disconnected components of the graph etc. However, proprocessing is
-   solver-specific.
-
-   This method can be implemented by derived classes in their solver-specific
-   way. Preprocessing may reveal unboundedness or unfeasibility of the
-   problem; if that happens, PreProcess() should properly set the `status'
-   field, that can then be read with MCFGetStatus() [see below].
- 
-   Note that preprocessing may destroy all the solution information. Also, it
-   may be allowed to change the data of the problem, such as costs/capacities
-   of the arcs.
-
-   A valid preprocessing is doing nothing, and that's what the default
-   implementation of this method (that is *not* pure virtual) does. */
-
-/*--------------------------------------------------------------------------*/
-
-   virtual inline void SetPar( int par , int val );
-
-/**< Set integer parameters of the algorithm.
-
-   @param par   is the parameter to be set; the enum MCFParam can be used, but
-                'par' is an int (every enum is an int) so that the method can
-                be extended by derived classes for the setting of their
-                parameters
-
-   @param value is the value to assign to the parameter.  
-
-   The base class implementation handles these parameters: 
-
-   - kMaxIter: the max number of iterations in which the MCF Solver can find
-               an optimal solution (default 0, which means no limit)
-
-   - kReopt:   tells the solver if it has to reoptimize. The implementation in
-               the base class sets a flag, the protected \c bool field \c
-               Senstv; if true (default) this field instructs the MCF solver to
-               to try to exploit the information about the latest optimal
-	       solution to speedup the optimization of the current problem,
-	       while if the field is false the MCF solver should restart the
-	       optimization "from scratch" discarding any previous information.
-	       Usually reoptimization speeds up the computation considerably,
-	       but this is not always true, especially if the data of the
-	       problem changes a lot.*/
-
-/*--------------------------------------------------------------------------*/
-
-   virtual inline void SetPar( int par , double val );
-
-/**< Set float parameters of the algorithm.
-
-   @param par   is the parameter to be set; the enum MCFParam can be used, but
-                'par' is an int (every enum is an int) so that the method can
-                be extended by derived classes for the setting of their
-                parameters
-
-   @param value is the value to assign to the parameter.  
-
-   The base class implementation handles these parameters: 
-
-   - kEpsFlw:  sets the tolerance for controlling if the flow on an arc is zero 
-               to val. This also sets the tolerance for controlling if a node
-	       deficit is zero (see kEpsDfct) to val * < max number of nodes >;
-	       this value should be safe for graphs in which any node has less
-	       than < max number of nodes > adjacent nodes, i.e., for all graphs
-	       but for very dense ones with "parallel arcs"
-
-   - kEpsDfct: sets the tolerance for controlling if a node deficit is zero to 
-               val, in case a better value than that autmatically set by
-               kEpsFlw (see above) is available (e.g., val * k would be good
-	       if no node has more than k neighbours)
-
-   - kEpsCst:  sets the tolerance for controlling if the reduced cost of an arc
-               is zero to val. A feasible solution satisfying eps-complementary
-               slackness, i.e., such that
-               \f[
-                RC[ i , j ] < - eps \Rightarrow X[ i , j ] = U[ ij ]
-               \f]
-               and
-               \f[
-                RC[ i , j ] > eps \Rightarrow X[ i , j ] == 0 ,
-               \f]
-               is known to be ( eps * n )-optimal.
-
-   - kMaxTime: sets the max time (in seconds) in which the MCF Solver can find
-               an optimal solution  (default 0, which means no limit). */
-
-/*--------------------------------------------------------------------------*/
-
-   virtual inline void GetPar( int par , int &val );
-
-/**< This method returns one of the integer parameter of the algorithm.
-
-   @param par  is the parameter to return [see SetPar( int ) for comments];
-
-   @param val  upon return, it will contain the value of the parameter.
-
-   The base class implementation handles the parameters kMaxIter and kReopt.
-   */
-
-/*--------------------------------------------------------------------------*/
-
-   virtual inline void GetPar( int par , double &val );
-
-/**< This method returns one of the integer parameter of the algorithm.
-
-   @param par  is the parameter to return [see SetPar( double ) for comments];
-
-   @param val  upon return, it will contain the value of the parameter.
-
-   The base class implementation handles the parameters kEpsFlw, kEpsDfct,
-   kEpsCst, and kMaxTime. */
-
-/*--------------------------------------------------------------------------*/
-
-   virtual void SetMCFTime( bool TimeIt = true )
-   {
-    if( TimeIt )
-     if( MCFt )
-      MCFt->ReSet();
-     else
-      MCFt = new OPTtimers();
-    else {
-     delete MCFt;
-     MCFt = NULL;
-     }
-    }
-
-/**< Allocate an OPTtimers object [see OPTtypes.h] to be used for timing the
-   methods of the class. The time can be read with TimeMCF() [see below]. By
-   default, or if SetMCFTime( false ) is called, no timing is done. Note that,
-   since all the relevant methods ot the class are pure virtual, MCFClass can
-   only manage the OPTtimers object, but it is due to derived classes to
-   actually implement the timing.
-
-   @note time accumulates over the calls: calling SetMCFTime(), however,
-         resets the counters, allowing to time specific groups of calls.
-
-   @note of course, setting kMaxTime [see SetPar() above] to any nonzero
-         value has no effect unless SetMCFTime( true ) has been called. */
-
-/*@} -----------------------------------------------------------------------*/
-/*---------------------- METHODS FOR SOLVING THE PROBLEM -------------------*/
-/*--------------------------------------------------------------------------*/
-/** @name Solving the problem
-    @{ */
-
-   virtual void SolveMCF( void ) = 0;
-
-/**< Solver of the Min Cost Flow Problem. Attempts to solve the MCF instance
-   currently loaded in the object. */
-
-/*--------------------------------------------------------------------------*/
-
-   inline int MCFGetStatus( void )
-   {
-    return( status );
-    }
-
-/**< Returns an int describing the current status of the MCF solver. Possible
-   return values are:
-
-   - kUnSolved    SolveMCF() has not been called yet, or the data of the
-                  problem has been changed since the last call;
-
-   - kOK          optimization has been carried out succesfully;
-
-   - kStopped     optimization have been stopped before that the stopping
-                  conditions of the solver applied, e.g. because of the
-                  maximum allowed number of "iterations" [see SetPar( int )]
-                  or the maximum allowed time [see SetPar( double )] has been
-		  reached; this is not necessarily an error, as it might just
-                  be required to re-call SolveMCF() giving it more "resources"
-                  in order to solve the problem;
-
-   - kUnfeasible  if the current MCF instance is (primal) unfeasible;
-
-   - kUnbounded   if the current MCF instance is (primal) unbounded (this can
-                  only happen if the solver actually allows F_INF capacities,
-                  which is nonstandard in the interface);
-
-   - kError       if there was an error during the optimization; this
-                  typically indicates that computation cannot be resumed,
-                  although solver-dependent ways of dealing with
-                  solver-dependent errors may exist.
-
-   MCFClass has a protected \c int \c member \c status \c that can be used by
-   derived classes to hold status information and that is returned by the
-   standard implementation of this method. Note that \c status \c is an
-   \c int \c and not an \c enum \c, and that an \c int \c is returned by this
-   method, in order to allow the derived classes to extend the set of return
-   values if they need to do so. */
-
-/*@} -----------------------------------------------------------------------*/
-/*---------------------- METHODS FOR READING RESULTS -----------------------*/
-/*--------------------------------------------------------------------------*/
-/** @name Reading flow solution
-    @{ */
-
-   virtual void MCFGetX( FRow F ,  Index_Set nms = NULL ,
-                         cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
-
-/**< Write the optimal flow solution in the vector F[]. If nms == NULL, F[]
-   will be in "dense" format, i.e., the flow relative to arc `i'
-   (i in 0 .. m - 1) is written in F[ i ]. If nms != NULL, F[] will be in 
-   "sparse" format, i.e., the indices of the nonzero elements in the flow
-   solution are written in nms (that is then Inf<Index>()-terminated) and
-   the flow value of arc nms[ i ] is written in F[ i ]. Note that nms is
-   *not* guaranteed to be ordered. Also, note that, unlike MCFGetRC() and
-   MCFGetPi() [see below], nms is an *output* of the method.
-
-   The parameters `strt' and `stp' allow to restrict the output of the method
-   to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual cFRow MCFGetX( void )
-   {
-    return( NULL );
-    }
-
-/**< Return a read-only pointer to an internal data structure containing the
-   flow solution in "dense" format. Since this may *not always be available*,
-   depending on the implementation, this method can (uniformly) return NULL.
-   This is done by the base class already, so a derived class that does not
-   have the information ready does not need to implement the method. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual bool HaveNewX( void )
-   {
-    return( false );
-    }
-
-/**< Return true if a different (approximately) optimal primal solution is
-   available. If the method returns true, then any subsequent call to (any
-   form of) MCFGetX() will return a different primal solution w.r.t. the one
-   that was being returned *before* the call to HaveNewX(). This solution need
-   not be optimal (although, ideally, it has to be "good); this can be checked
-   by comparing its objective function value, that will be returned by a call
-   to MCFGetFO() [see below].
-
-   Any subsequent call of HaveNewX() that returns true produces a new
-   solution, until the first that returns false; from then on, no new
-   solutions will be generated until something changes in the problem's
-   data.
-
-   Note that a default implementation of HaveNewX() is provided which is
-   good for those solvers that only produce one optimal primal solution. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Reading potentials
-    @{ */
-
-   virtual void MCFGetPi(  CRow P ,  cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
-
-/**< Writes the optimal node potentials in the vector P[]. If nms == NULL,
-   the node potential of node `i' (i in 0 .. n - 1) is written in P[ i ]
-   (note that here node names always start from zero, regardless to the value
-   of USENAME0). If nms != NULL, it must point to a vector of indices in
-   0 .. n - 1 (ordered in increasing sense and Inf<Index>()-terminated), and
-   the node potential of nms[ i ] is written in P[ i ]. Note that, unlike
-   MCFGetX() above, nms is an *input* of the method.
-
-   The parameters `strt' and `stp' allow to restrict the output of the method
-   to all and only the nodes `i' with strt <= i < min( MCFn() , stp ). `strt'
-   and `stp' work in "&&" with nms; that is, if nms != NULL then only the
-   values corresponding to nodes which are *both* in nms[] and whose index is
-   in the correct range are returned. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual cCRow MCFGetPi( void )
-   {
-    return( NULL );
-    }
-
-/**< Return a read-only pointer to an internal data structure containing the
-   node potentials. Since this may *not always be available*, depending on
-   the implementation, this method can (uniformly) return NULL.
-   This is done by the base class already, so a derived class that does not
-   have the information ready does not need to implement the method. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual bool HaveNewPi( void )
-   {
-    return( false );
-    }
-
-/**< Return true if a different (approximately) optimal dual solution is
-   available. If the method returns true, then any subsequent call to (any
-   form of) MCFGetPi() will return a different dual solution, and MCFGetRC()
-   [see below] will return the corresponding reduced costs. The new solution
-   need not be optimal (although, ideally, it has to be "good); this can be
-   checked by comparing its objective function value, that will be returned
-   by a call to MCFGetDFO() [see below].
-
-   Any subsequent call of HaveNewPi() that returns true produces a new
-   solution, until the first that returns false; from then on, no new
-   solutions will be generated until something changes in the problem's
-   data.
-
-   Note that a default implementation of HaveNewPi() is provided which is
-   good for those solvers that only produce one optimal dual solution. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Reading reduced costs
-    @{ */
-
-   virtual void MCFGetRC(  CRow CR ,  cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
-
-/**< Write the reduced costs corresponding to the current dual solution in
-   RC[]. If nms == NULL, the reduced cost of arc `i' (i in 0 .. m - 1) is
-   written in RC[ i ]; if nms != NULL, it must point to a vector of indices
-   in 0 .. m - 1 (ordered in increasing sense and Inf<Index>()-terminated),
-   and the reduced cost of arc nms[ i ] is written in RC[ i ]. Note that,
-   unlike MCFGetX() above, nms is an *input* of the method.
-
-   The parameters `strt' and `stp' allow to restrict the output of the method
-   to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
-   and `stp' work in "&&" with nms; that is, if nms != NULL then only the
-   values corresponding to arcs which are *both* in nms[] and whose index is
-   in the correct range are returned.
-
-   @note the output of MCFGetRC() will change after any call to HaveNewPi()
-         [see above] which returns true. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual cCRow MCFGetRC( void )
-   {
-    return( NULL );
-    }
-
-/**< Return a read-only pointer to an internal data structure containing the
-   reduced costs. Since this may *not always be available*, depending on the
-   implementation, this method can (uniformly) return NULL.
-   This is done by the base class already, so a derived class that does not
-   have the information ready does not need to implement the method.
-
-   @note the output of MCFGetRC() will change after any call to HaveNewPi()
-         [see above] which returns true. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual CNumber MCFGetRC( cIndex i ) = 0;
-
-/**< Return the reduced cost of the i-th arc. This information should be
-   cheapily available in most implementations.
-
-   @note the output of MCFGetRC() will change after any call to HaveNewPi()
-         [see above] which returns true. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Reading the objective function value
-   @{ */
-
-   virtual FONumber MCFGetFO( void ) = 0;
-
-/**< Return the objective function value of the primal solution currently
-   returned by MCFGetX().
-   
-   If MCFGetStatus() == kOK, this is guaranteed to be the optimal objective
-   function value of the problem (to within the optimality tolerances), but
-   only prior to any call to HaveNewX() that returns true. MCFGetFO()
-   typically returns Inf<FONumber>() if MCFGetStatus() == kUnfeasible and
-   - Inf<FONumber>() if MCFGetStatus() == kUnbounded. If MCFGetStatus() ==
-   kStopped and MCFGetFO() returns a finite value, it must be an upper bound
-   on the optimal objective function value (typically, the objective function
-   value of one primal feasible solution). */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual FONumber MCFGetDFO( void )
-   {
-    switch( MCFGetStatus() ) {
-     case( kUnSolved ):
-     case( kStopped ):
-     case( kError ):    return( -Inf<FONumber>() );
-     default:           return( MCFGetFO() );
-     }
-    }
-
-/**< Return the objective function value of the dual solution currently
-   returned by MCFGetPi() / MCFGetRC(). This value (possibly) changes after
-   any call to HaveNewPi() that returns true. The relations between
-   MCFGetStatus() and MCFGetDFO() are analogous to these of MCFGetFO(),
-   except that a finite value corresponding to kStopped must be a lower
-   bound on the optimal objective function value (typically, the objective
-   function value one dual feasible solution).
-
-   A default implementation is provided for MCFGetDFO(), which is good for
-   MCF solvers where the primal and dual optimal solution values always
-   are identical (except if the problem is unfeasible/unbounded). */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Getting unfeasibility certificate
-   @{ */
-
-   virtual FNumber MCFGetUnfCut( Index_Set Cut )
-   {
-    *Cut = Inf<Index>();
-    return( 0 );
-    }
-
-/**< Return an unfeasibility certificate. In an unfeasible MCF problem,
-   unfeasibility can always be reduced to the existence of a cut (subset of
-   nodes of the graph) such as either:
-
-   - the inverse of the deficit of the cut (the sum of all the deficits of the
-     nodes in the cut) is larger than the forward capacity of the cut (sum of
-     the capacities of forward arcs in the cut); that is, the nodes in the
-     cut globally produce more flow than can be routed to sinks outside the
-     cut;
-
-   - the deficit of the cut is larger than the backward capacity of the cut
-     (sum of the capacities of backward arcs in the cut); that is, the nodes
-     in the cut globally require more flow than can be routed to them from
-     sources outside the cut.
-
-   When detecting unfeasibility, MCF solvers are typically capable to provide
-   one such cut. This information can be useful - typically, the only way to
-   make the problem feasible is to increase the capacity of at least one of
-   the forward/backward arcs of the cut -, and this method is provided for
-   getting it. It can be called only if MCFGetStatus() == kUnfeasible, and
-   should write in Cut the set of names of nodes in the unfeasible cut (note
-   that node names depend on USENAME0), Inf<Index>()-terminated, returning the
-   deficit of the cut (which allows to distinguish which of the two cases
-   above hold). In general, no special properties can be expected from the
-   returned cut, but solvers should be able to provide e.g. "small" cuts.
-
-   However, not all solvers may be (easily) capable of providing this
-   information; thus, returning 0 (no cut) is allowed, as in the base class
-   implementation, to signify that this information is not available. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Getting unboundedness certificate
-   @{ */
-
-   virtual Index MCFGetUnbCycl( Index_Set Pred , Index_Set ArcPred )
-   {
-    return( Inf<Index>() );
-    }
-
-/**< Return an unboundedness certificate. In an unbounded MCF problem,
-   unboundedness can always be reduced to the existence of a directed cycle
-   with negative cost and all arcs having infinite capacity.
-   When detecting unboundedness, MCF solvers are typically capable to provide
-   one such cycle. This information can be useful, and this method is provided
-   for getting it. It can be called only if MCFGetStatus() == kUnbounded, and
-   writes in Pred[] and ArcPred[], respectively, the node and arc predecessor
-   function of the cycle. That is, if node `i' belongs to the cycle then
-   `Pred[ i ]' contains the name of the predecessor of `j' of `i' in the cycle
-   (note that node names depend on USENAME0), and `ArcPred[ i ]' contains the
-   index of the arc joining the two (note that in general there may be
-   multiple copies of each arc). Entries of the vectors for nodes not
-   belonging to the cycle are in principle undefined, and the name of one node
-   belonging to the cycle is returned by the method. 
-   Note that if there are multiple cycles with negative costs this method
-   will return just one of them (finding the cycle with most negative cost
-   is an NO-hard problem), although solvers should be able to produce cycles
-   with "large negative" cost.
-
-   However, not all solvers may be (easily) capable of providing this
-   information; thus, returning Inf<Index>() is allowed, as in the base class
-   implementation, to signify that this information is not available. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Saving/restoring the state of the solver
-    @{ */
-
-   virtual MCFStatePtr MCFGetState( void )
-   {
-    return( NULL );
-    }
-
- /**< Save the state of the MCF solver. The MCFClass interface supports the
-   notion of saving and restoring the state of the MCF solver, such as the
-   current/optimal basis in a simplex solver. The "empty" class MCFState is
-   defined as a placeholder for state descriptions.
-
-   MCFGetState() creates and returns a pointer to an object of (a proper
-   derived class of) class MCFState which describes the current state of the
-   MCF solver. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void MCFPutState( MCFStatePtr S ) {}
-
-/**< Restore the solver to the state in which it was when the state `S' was
-   created with MCFGetState() [see above].
-
-   The typical use of this method is the following: a MCF problem is solved
-   and the "optimal state" is set aside. Then the problem changes and it is
-   re-solved. Then, the problem has to be changed again to a form that is
-   close to the original one but rather different from the second one (think
-   of a long backtracking in a Branch & Bound) to which the current "state"
-   referes. Then, the old optimal state can be expected to provide a better
-   starting point for reoptimization [see ReOptimize() below].
-
-   Note, however, that the state is only relative to the optimization process,
-   i.e., this operation is meaningless if the data of the problem has changed
-   in the meantime. So, if a state has to be used for speeding up
-   reoptimization, the following has to be done:
-
-   - first, the data of the solver is brought back to *exactly* the same as
-     it was at the moment where the state `S' was created (prior than this
-     operation a ReOptimize( false ) call is probably advisable);
-
-   - then, MCFPutState() is called (and ReOptimize( true ) is called);
-
-   - only afterwards the data of the problem is changed to the final state
-     and the problem is solved.
-
-   A "put state" operation does not "deplete" the state, which can therefore
-   be used more than once. Indeed, a state is constructed inside the solver
-   for each call to MCFGetState(), but the solver never deletes statuses;
-   this has to be done on the outside when they are no longer needed (the
-   solver must be "resistent" to deletion of the state at any moment).
-
-   Since not all the MCF solvers reoptimize (efficiently enough to make these
-   operations worth), an "empty" implementation that does nothing is provided
-   by the base class. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Time the code
-    @{ */
-
-   void TimeMCF( double &t_us , double &t_ss )
-   {
-    t_us = t_ss = 0;
-    if( MCFt )
-     MCFt->Read( t_us , t_ss ); 
-    }
-
-/**< Time the code. If called within any of the methods of the class that are
-   "actively timed" (this depends on the subclasses), this method returns the
-   user and sistem time (in seconds) since the start of that method. If
-   methods that are actively timed call other methods that are actively
-   timed, TimeMCF() returns the (...) time since the beginning of the
-   *outer* actively timed method. If called outside of any actively timed
-   method, this method returns the (...) time spent in all the previous
-   executions of all the actively timed methods of the class.
-
-   Implementing the proper calls to MCFt->Start() and MCFt->Stop() is due to
-   derived classes; these should at least be placed at the beginning and at
-   the end, respectively, of SolveMCF() and presumably the Chg***() methods,
-   that is, at least these methods should be "actively timed". */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   double TimeMCF( void )
-   {
-    return( MCFt ? MCFt->Read() : 0 );
-    }
-
-/**< Like TimeMCF(double,double) [see above], but returns the total time. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Check the solutions
-   @{ */
-
-   inline void CheckPSol( void );
-
-/**< Check that the primal solution returned by the solver is primal feasible.
-   (to within the tolerances set by SetPar(kEps****) [see above], if any). Also,
-   check that the objective function value is correct.
-
-   This method is implemented by the base class, using the above methods
-   for collecting the solutions and the methods of the next section for
-   reading the data of the problem; as such, they will work for any derived
-   class that properly implements all these methods. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   inline void CheckDSol( void );
-
-/**< Check that the dual solution returned by the solver is dual feasible.
-   (to within the tolerances set by SetPar(kEps****) [see above], if any). Also,
-   check that the objective function value is correct.
-
-   This method is implemented by the base class, using the above methods
-   for collecting the solutions and the methods of the next section for
-   reading the data of the problem; as such, they will work for any derived
-   class that properly implements all these methods. */
-
-/*@} -----------------------------------------------------------------------*/
-/*-------------- METHODS FOR READING THE DATA OF THE PROBLEM ---------------*/
-/*--------------------------------------------------------------------------*/
-/** @name Reading graph size
-    @{ */
-
-   inline Index MCFnmax( void )
-   {
-    return( nmax );
-    }
-
-/**< Return the maximum number of nodes for this instance of MCFClass.
-   The implementation of the method in the base class returns the protected
-   fields \c nmax, which is provided for derived classes to hold this
-   information. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   inline Index MCFmmax( void )
-   {
-    return( mmax );
-    }
-
-/**< Return the maximum number of arcs for this instance of MCFClass.
-   The implementation of the method in the base class returns the protected
-   fields \c mmax, which is provided for derived classes to hold this
-   information. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   inline Index MCFn( void )
-   {
-    return( n );
-    }
-
-/**< Return the number of nodes in the current graph.
-   The implementation of the method in the base class returns the protected
-   fields \c n, which is provided for derived classes to hold this
-   information. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   inline Index MCFm( void )
-   {
-    return( m );
-    }
-
-/**< Return the number of arcs in the current graph.
-   The implementation of the method in the base class returns the protected
-   fields \c m, which is provided for derived classes to hold this
-   information. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Reading graph topology
-   @{ */
-
-   virtual void MCFArcs( Index_Set Startv ,  Index_Set Endv ,
-                         cIndex_Set nms = NULL , cIndex strt = 0 ,
-                         Index stp = Inf<Index>() ) = 0;
-
-/**< Write the starting (tail) and ending (head) nodes of the arcs in Startv[]
-   and Endv[]. If nms == NULL, then the information relative to all arcs is
-   written into Startv[] and Endv[], otherwise Startv[ i ] and Endv[ i ]
-   contain the information relative to arc nms[ i ] (nms[] must be
-   Inf<Index>()-terminated).
-
-   The parameters `strt' and `stp' allow to restrict the output of the method
-   to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
-   and `stp' work in "&&" with nms; that is, if nms != NULL then only the
-   values corresponding to arcs which are *both* in nms and whose index is in
-   the correct range are returned.
-
-   Startv or Endv can be NULL, meaning that only the other information is
-   required.
-
-   @note If USENAME0 == 0 then the returned node names will be in the range
-         1 .. n, while if USENAME0 == 1 the returned node names will be in
-         the range 0 .. n - 1.
-
-   @note If the graph is "dynamic", be careful to use MCFn() e MCFm() to
-         properly choose the dimension of nodes and arcs arrays. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual Index MCFSNde( cIndex i ) = 0;
-
-/**< Return the starting (tail) node of the arc `i'.
-
-   @note If USENAME0 == 0 then the returned node names will be in the range
-         1 .. n, while if USENAME0 == 1 the returned node names will be in
-         the range 0 .. n - 1. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual Index MCFENde( cIndex i ) = 0;
-
-/**< Return the ending (head) node of the arc `i'.
-
-   @note If USENAME0 == 0 then the returned node names will be in the range
-         1 .. n, while if USENAME0 == 1 the returned node names will be in
-         the range 0 .. n - 1. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual cIndex_Set MCFSNdes( void )
-   {
-    return( NULL );
-    }
-
-/**< Return a read-only pointer to an internal vector containing the starting
-   (tail) nodes for each arc. Since this may *not always be available*,
-   depending on the implementation, this method can (uniformly) return NULL.
-   This is done by the base class already, so a derived class that does not
-   have the information ready does not need to implement the method. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual cIndex_Set MCFENdes( void )
-   {
-    return( NULL );
-    }
-
-/**< Return a read-only pointer to an internal vector containing the ending
-   (head) nodes for each arc. Since this may *not always be available*,
-   depending on the implementation, this method can (uniformly) return NULL.
-   This is done by the base class already, so a derived class that does not
-   have the information ready does not need to implement the method. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Reading arc costs
-   @{ */
-
-   virtual void MCFCosts( CRow Costv , cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
-
-/**< Write the arc costs into Costv[]. If nms == NULL, then all the costs are
-   written, otherwise Costv[ i ] contains the information relative to arc
-   nms[ i ] (nms[] must be Inf<Index>()-terminated).
-
-   The parameters `strt' and `stp' allow to restrict the output of the method
-   to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
-   and `stp' work in "&&" with nms; that is, if nms != NULL then only the
-   values corresponding to arcs which are *both* in nms and whose index is in
-   the correct range are returned. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual CNumber MCFCost( cIndex i ) = 0;
-
-/**< Return the cost of the i-th arc. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual cCRow MCFCosts( void )
-   {
-    return( NULL );
-    }
-
-/**< Return a read-only pointer to an internal vector containing the arc
-   costs. Since this may *not always be available*, depending on the
-   implementation, this method can (uniformly) return NULL.
-   This is done by the base class already, so a derived class that does not
-   have the information ready does not need to implement the method. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void MCFQCoef( CRow Qv , cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() )
-
-/**< Write the quadratic coefficients of the arc costs into Qv[]. If
-   nms == NULL, then all the costs are written, otherwise Costv[ i ] contains
-   the information relative to arc nms[ i ] (nms[] must be
-   Inf<Index>()-terminated).
-
-   The parameters `strt' and `stp' allow to restrict the output of the method
-   to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
-   and `stp' work in "&&" with nms; that is, if nms != NULL then only the
-   values corresponding to arcs which are *both* in nms and whose index is in
-   the correct range are returned.
-
-   Note that the method is *not* pure virtual: an implementation is provided
-   for "pure linear" MCF solvers that only work with all zero quadratic
-   coefficients. */
-   {
-    if( nms ) {
-     while( *nms < strt )
-      nms++;
-
-     for( Index h ; ( h = *(nms++) ) < stp ; )
-      *(Qv++) = 0;
-     }
-    else {
-     if( stp > m )
-      stp = m;
-
-     while( stp-- > strt )
-      *(Qv++) = 0;
-     }
-    }
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual CNumber MCFQCoef( cIndex i )
-
-/**< Return the quadratic coefficients of the cost of the i-th arc. Note that
-   the method is *not* pure virtual: an implementation is provided for "pure
-   linear" MCF solvers that only work with all zero quadratic coefficients. */
-   {
-    return( 0 );
-    }
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual cCRow MCFQCoef( void )
-   {
-    return( NULL );
-    }
-
-/**< Return a read-only pointer to an internal vector containing the arc
-   costs. Since this may *not always be available*, depending on the
-   implementation, this method can (uniformly) return NULL.
-   This is done by the base class already, so a derived class that does not
-   have the information ready (such as "pure linear" MCF solvers that only
-   work with all zero quadratic coefficients) does not need to implement the
-   method. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Reading arc capacities
-    @{ */
-
-   virtual void MCFUCaps( FRow UCapv , cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
-
-/**< Write the arc capacities into UCapv[]. If nms == NULL, then all the
-   capacities are written, otherwise UCapv[ i ] contains the information
-   relative to arc nms[ i ] (nms[] must be Inf<Index>()-terminated).
-
-   The parameters `strt' and `stp' allow to restrict the output of the method
-   to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
-   and `stp' work in "&&" with nms; that is, if nms != NULL then only the
-   values corresponding to arcs which are *both* in nms and whose index is in
-   the correct range are returned. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual FNumber MCFUCap( cIndex i ) = 0;
-
-/**< Return the capacity of the i-th arc. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual cFRow MCFUCaps( void )
-   {
-    return( NULL );
-    }
-
-/**< Return a read-only pointer to an internal vector containing the arc
-   capacities. Since this may *not always be available*, depending on the
-   implementation, this method can (uniformly) return NULL.
-   This is done by the base class already, so a derived class that does not
-   have the information ready does not need to implement the method. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Reading node deficits
-    @{ */
-
-   virtual void MCFDfcts( FRow Dfctv , cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
-
-/**< Write the node deficits into Dfctv[]. If nms == NULL, then all the
-   defcits are written, otherwise Dfctvv[ i ] contains the information
-   relative to node nms[ i ] (nms[] must be Inf<Index>()-terminated).
-
-   The parameters `strt' and `stp' allow to restrict the output of the method
-   to all and only the nodes `i' with strt <= i < min( MCFn() , stp ). `strt'
-   and `stp' work in "&&" with nms; that is, if nms != NULL then only the
-   values corresponding to nodes which are *both* in nms and whose index is in
-   the correct range are returned.
-
-   @note Here node "names" (strt and stp, those contained in nms[] or `i'
-         in MCFDfct()) go from 0 to n - 1, regardless to the value of
-         USENAME0; hence, if USENAME0 == 0 then the first node is "named 1",
-         but its deficit is returned by MCFDfcts( Dfctv , NULL , 0 , 1 ). */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual FNumber MCFDfct( cIndex i ) = 0;
-
-/**< Return the deficit of the i-th node.
-
-   @note Here node "names" go from 0 to n - 1, regardless to the value of
-         USENAME0; hence, if USENAME0 == 0 then the first node is "named 1",
-         but its deficit is returned by MCFDfct( 0 ). */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual cFRow MCFDfcts( void )
-   {
-    return( NULL );
-    }
-
-/**< Return a read-only pointer to an internal vector containing the node
-   deficits.  Since this may *not always be available*, depending on the
-   implementation, this method can (uniformly) return NULL.
-   This is done by the base class already, so a derived class that does not
-   have the information ready does not need to implement the method.
-
-   @note Here node "names" go from 0 to n - 1, regardless to the value of
-         USENAME0; hence, if USENAME0 == 0 then the first node is "named 1",
-         but its deficit is contained in MCFDfcts()[ 0 ]. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Write problem to file
-    @{ */
-
-   virtual inline void WriteMCF( ostream &oStrm , int frmt = 0 );
-
-/**< Write the current MCF problem to an ostream. This may be useful e.g. for
-   debugging purposes.
-
-   The base MCFClass class provides output in two different formats, depending
-   on the value of the parameter frmt:
-
-   - kDimacs  the problem is written in DIMACS standard format, read by most
-              MCF codes available;
-
-   - kMPS     the problem is written in the "modern version" (tab-separated)
-              of the MPS format, read by most LP/MIP solvers;
-
-    - kFWMPS  the problem is written in the "old version" (fixed width
-              fields) of the MPS format; this is read by most LP/MIP
-              solvers, but some codes still require the old format.
-
-   The implementation of WriteMCF() in the base class uses all the above
-   methods for reading the data; as such it will work for any derived class
-   that properly implements this part of the interface, but it may not be
-   very efficient. Thus, the method is virtual to allow the derived classes
-   to either re-implement WriteMCF() for the above two formats in a more
-   efficient way, and/or to extend it to support other solver-specific
-   formats.
-
-   @note None of the above two formats supports quadratic MCFs, so if
-         nonzero quadratic coefficients are present, they are just ignored.
-   */
-
-/*@} -----------------------------------------------------------------------*/
-/*------------- METHODS FOR ADDING / REMOVING / CHANGING DATA --------------*/
-/*--------------------------------------------------------------------------*/
-/** @name Changing the costs
-    @{ */
-
-   virtual void ChgCosts( cCRow NCost , cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
-
-/**< Change the arc costs. In particular, change the costs that are:
-
-   - listed in into the vector of indices `nms' (ordered in increasing
-     sense and Inf<Index>()-terminated),
-
-   - *and* whose name belongs to the interval [`strt', `stp').
-
-   That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th cost will be
-   changed to NCost[ i ]. If nms == NULL (as the default), *all* the
-   entries in the given range will be changed; if stp > MCFm(), then the
-   smaller bound is used.
-
-   @note changing the costs of arcs that *do not exist* is *not allowed*;
-         only arcs which have not been closed/deleted [see CloseArc() /
-         DelArc() below and LoadNet() above about C_INF costs] can be
-         touched with these methods. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void ChgCost( Index arc , cCNumber NCost ) = 0;
-
-/**< Change the cost of the i-th arc.
-
-   @note changing the costs of arcs that *do not exist* is *not allowed*;
-         only arcs which have not been closed/deleted [see CloseArc() /
-         DelArc() below and LoadNet() above about C_INF costs] can be
-         touched with these methods. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void ChgQCoef( cCRow NQCoef = NULL , cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() )
-
-/**< Change the quadratic coefficients of the arc costs. In particular,
-   change the coefficients that are:
-
-   - listed in into the vector of indices `nms' (ordered in increasing
-     sense and Inf<Index>()-terminated),
-
-   - *and* whose name belongs to the interval [`strt', `stp').
-
-   That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th cost will be
-   changed to NCost[ i ]. If nms == NULL (as the default), *all* the
-   entries in the given range will be changed; if stp > MCFm(), then the
-   smaller bound is used. If NQCoef == NULL, all the specified coefficients
-   are set to zero.
-
-   Note that the method is *not* pure virtual: an implementation is provided
-   for "pure linear" MCF solvers that only work with all zero quadratic
-   coefficients.
-
-   @note changing the costs of arcs that *do not exist* is *not allowed*;
-         only arcs which have not been closed/deleted [see CloseArc() /
-         DelArc() below and LoadNet() above about C_INF costs] can be
-         touched with these methods. */
-   {
-     assert(NQCoef);
-   }
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void ChgQCoef(  Index arc , cCNumber NQCoef )
-
-/**< Change the quadratic coefficient of the cost of the i-th arc.
-
-   Note that the method is *not* pure virtual: an implementation is provided
-   for "pure linear" MCF solvers that only work with all zero quadratic
-   coefficients.
-
-   @note changing the costs of arcs that *do not exist* is *not allowed*;
-         only arcs which have not been closed/deleted [see CloseArc() /
-         DelArc() below and LoadNet() above about C_INF costs] can be
-         touched with these methods. */
-   {
-     assert(NQCoef);
-   }
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Changing the capacities
-    @{ */
-
-   virtual void ChgUCaps( cFRow NCap , cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
-
-/**< Change the arc capacities. In particular, change the capacities that are:
-
-   - listed in into the vector of indices `nms' (ordered in increasing sense
-     and Inf<Index>()-terminated),
-
-   - *and* whose name belongs to the interval [`strt', `stp').
-
-   That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th capacity will be
-   changed to NCap[ i ]. If nms == NULL (as the default), *all* the
-   entries in the given range will be changed; if stp > MCFm(), then the
-   smaller bound is used.
-
-   @note changing the capacities of arcs that *do not exist* is *not allowed*;
-         only arcs that have not been closed/deleted [see CloseArc() /
-         DelArc() below and LoadNet() above about C_INF costs] can be
-         touched with these methods. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void ChgUCap( Index arc , cFNumber NCap ) = 0;
-
-/**< Change the capacity of the i-th arc.
-
-   @note changing the capacities of arcs that *do not exist* is *not allowed*;
-         only arcs that have not been closed/deleted [see CloseArc() /
-         DelArc() below and LoadNet() above about C_INF costs] can be
-         touched with these methods. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Changing the deficits
-    @{ */
-
-   virtual void ChgDfcts( cFRow NDfct , cIndex_Set nms = NULL ,
-                          cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
-
-/**< Change the node deficits. In particular, change the deficits that are:
-
-   - listed in into the vector of indices `nms' (ordered in increasing sense
-     and Inf<Index>()-terminated),
-
-   - *and* whose name belongs to the interval [`strt', `stp').
-
-   That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th deficit will be
-   changed to NDfct[ i ]. If nms == NULL (as the default), *all* the
-   entries in the given range will be changed; if stp > MCFn(), then the
-   smaller bound is used.
-
-   Note that, in ChgDfcts(), node "names" (strt, stp or those contained
-   in nms[]) go from 0 to n - 1, regardless to the value of USENAME0; hence,
-   if USENAME0 == 0 then the first node is "named 1", but its deficit can be
-   changed e.g. with ChgDfcts( &new_deficit , NULL , 0 , 1 ).
-
-   @note changing the capacities of nodes that *do not exist* is *not
-         allowed*; only nodes that have not been deleted [see DelNode()
-         below] can be touched with these methods. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void ChgDfct( Index node , cFNumber NDfct ) = 0;
-
-/**< Change the deficit of the i-th node.
-
-   Note that, in ChgDfct[s](), node "names" (i, strt/ stp or those contained
-   in nms[]) go from 0 to n - 1, regardless to the value of USENAME0; hence,
-   if USENAME0 == 0 then the first node is "named 1", but its deficit can be
-   changed e.g. with ChgDfcts( 0 , new_deficit ).
-
-   @note changing the capacities of nodes that *do not exist* is *not
-         allowed*; only nodes that have not been deleted [see DelNode()
-         below] can be touched with these methods. */
-
-/*@} -----------------------------------------------------------------------*/
-/** @name Changing graph topology
-    @{ */
-
-   virtual void CloseArc( cIndex name ) = 0;
-
-/**< "Close" the arc `name'. Although all the associated information (name,
-   cost, capacity, end and start node) is kept, the arc is removed from the
-   problem until OpenArc( i ) [see below] is called.
-
-   "closed" arcs always have 0 flow, but are otherwise counted as any other
-   arc; for instance, MCFm() does *not* decrease as an effect of a call to
-   CloseArc(). How this closure is implemented is solver-specific. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual bool IsClosedArc( cIndex name ) = 0;
-
-/**< IsClosedArc() returns true if and only if the arc `name' is closed. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void DelNode( cIndex name ) = 0;
-
-/**< Delete the node `name'.
-
-   For any value of `name', all incident arcs to that node are closed [see
-   CloseArc() above] (*not* Deleted, see DelArc() below) and the deficit is
-   set to zero.
-
-   Il furthermore `name' is the last node, the number of nodes as reported by
-   MCFn() is reduced by at least one, until the n-th node is not a deleted
-   one. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void OpenArc( cIndex name ) = 0;
-
-/**< Restore the previously closed arc `name'. It is an error to open an arc
-   that has not been previously closed. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual Index AddNode( cFNumber aDfct ) = 0;
-
-/**< Add a new node with deficit aDfct, returning its name. Inf<Index>() is
-   returned if there is no room for a new node. Remember that the node names
-   are either { 0 .. nmax - 1 } or { 1 .. nmax }, depending on the value of
-   USENAME0. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void ChangeArc( cIndex name , cIndex nSN = Inf<Index>() ,
-                           cIndex nEN = Inf<Index>() ) = 0;
-
-/**< Change the starting and/or ending node of arc `name' to nSN and nEN.
-   Each parameter being Inf<Index>() means to leave the previous starting or
-   ending node untouched. When this method is called `name' can be either the
-   name of a "normal" arc or that of a "closed" arc [see CloseArc() above]:
-   in the latter case, at the end of ChangeArc() the arc is *still closed*,
-   and it remains so until OpenArc( name ) [see above] is called. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual void DelArc( cIndex name ) = 0;
-
-/**< Delete the arc `name'. Unlike "closed" arcs, all the information
-   associated with a deleted arc is lost and `name' is made available as a
-   name for new arcs to be created with AddArc() [see below].
-
-   Il furthermore `name' is the last arc, the number of arcs as reported by
-   MCFm() is reduced by at least one, until the m-th arc is not a deleted
-   one. Otherwise, the flow on the arc is always ensured to be 0. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual bool IsDeletedArc( cIndex name ) = 0;
-
-/**< Return true if and only if the arc `name' is deleted. It should only
-   be called with name < MCFm(), as every other arc is deleted by
-   definition. */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-   virtual Index AddArc( cIndex Start , cIndex End , cFNumber aU ,
-                         cCNumber aC ) = 0; 
-
-/**< Add the new arc ( Start , End ) with cost aC and capacity aU,
-   returning its name. Inf<Index>() is returned if there is no room for a new
-   arc. Remember that arc names go from 0 to mmax - 1. */
-
-/*@} -----------------------------------------------------------------------*/
-/*------------------------------ DESTRUCTOR --------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @name Destructor
-    @{ */
-
-   virtual ~MCFClass()
-   {
-    delete MCFt;
-    }
-
-/**< Destructor of the class. The implementation in the base class only
-   deletes the MCFt field. It is virtual, as it should be. */
-
-/*@} -----------------------------------------------------------------------*/
-/*-------------------- PROTECTED PART OF THE CLASS -------------------------*/
-/*--------------------------------------------------------------------------*/
-/*--                                                                      --*/
-/*--  The following fields are believed to be general enough to make it   --*/
-/*--  worth adding them to the abstract base class.                       --*/
-/*--                                                                      --*/
-/*--------------------------------------------------------------------------*/
-
- protected:
-
-/*--------------------------------------------------------------------------*/
-/*--------------------------- PROTECTED METHODS ----------------------------*/
-/*--------------------------------------------------------------------------*/
-/*-------------------------- MANAGING COMPARISONS --------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @name Managing comparisons.
-    The following methods are provided for making it easier to perform
-    comparisons, with and without tolerances.
-    @{ */
-
-/*--------------------------------------------------------------------------*/
-/** true if flow x is equal to zero (possibly considering tolerances). */
-        
-   template<class T>
-   inline bool ETZ( T x , const T eps )
-   {
-    if( numeric_limits<T>::is_integer )
-     return( x == 0 );
-    else 
-     return( (x <= eps ) && ( x >= -eps ) );
-    }
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-/** true if flow x is greater than zero (possibly considering tolerances). */
-
-   template<class T>
-   inline bool GTZ( T x , const T eps )
-   {
-    if( numeric_limits<T>::is_integer )
-     return( x > 0 );
-    else
-     return( x > eps );
-    }
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-/** true if flow x is greater than or equal to zero (possibly considering
-    tolerances). */
-
-   template<class T>
-   inline bool GEZ( T x , const T eps )
-   {
-    if( numeric_limits<T>::is_integer )
-     return( x >= 0 );
-    else
-     return( x >= - eps );
-    }
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-/** true if flow x is less than zero (possibly considering tolerances). */
-
-   template<class T>
-   inline bool LTZ( T x , const T eps )
-   {
-    if( numeric_limits<T>::is_integer )
-     return( x < 0 );
-    else
-     return( x < - eps );
-    }
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-/** true if flow x is less than or equal to zero (possibly considering
-    tolerances). */
-
-   template<class T>
-   inline bool LEZ( T x , const T eps )
-   {
-    if( numeric_limits<T>::is_integer )
-     return( x <= 0 );
-    else
-     return( x <= eps );
-    }
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-/** true if flow x is greater than flow y (possibly considering tolerances).
- */
-
-   template<class T>
-   inline bool GT( T x , T y , const T eps )
-   {
-    if( numeric_limits<T>::is_integer )
-     return( x > y );
-    else
-     return( x > y + eps );
-    }
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-/** true if flow x is less than flow y (possibly considering tolerances). */
-
-   template<class T>
-   inline bool LT( T x , T y , const T eps )
-   {
-    if( numeric_limits<T>::is_integer )
-     return( x < y );
-    else
-     return( x < y - eps );
-    }
-
-/*@} -----------------------------------------------------------------------*/
-/*---------------------- PROTECTED DATA STRUCTURES -------------------------*/
-/*--------------------------------------------------------------------------*/
-
- Index n;      ///< total number of nodes
- Index nmax;   ///< maximum number of nodes
-
- Index m;      ///< total number of arcs
- Index mmax;   ///< maximum number of arcs
-
- int status;   ///< return status, see the comments to MCFGetStatus() above.
-               ///< Note that the variable is defined int to allow derived
-               ///< classes to return their own specialized status codes
- bool Senstv;  ///< true <=> the latest optimal solution should be exploited
-
- OPTtimers *MCFt;   ///< timer for performances evaluation
-
- FNumber EpsFlw;   ///< precision for comparing arc flows / capacities
- FNumber EpsDfct;  ///< precision for comparing node deficits
- CNumber EpsCst;   ///< precision for comparing arc costs
-
- double MaxTime;   ///< max time (in seconds) in which MCF Solver can find 
-                   ///< an optimal solution (0 = no limits)
- int MaxIter;      ///< max number of iterations in which MCF Solver can find 
-                   ///< an optimal solution (0 = no limits)
-
-/*--------------------------------------------------------------------------*/
-
- };   // end( class MCFClass )
-
-/*@} -----------------------------------------------------------------------*/
-/*------------------- inline methods implementation ------------------------*/
-/*--------------------------------------------------------------------------*/
-
-inline void MCFClass::LoadDMX( istream &DMXs , bool IsQuad )
-{
- // read first non-comment line - - - - - - - - - - - - - - - - - - - - - - -
- // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
- char c;
- for(;;) {
-  if( ! ( DMXs >> c ) )
-   throw( MCFException( "LoadDMX: error reading the input stream" ) );
-
-  if( c != 'c' )  // if it's not a comment
-   break;
-
-  do              // skip the entire line
-   if( ! DMXs.get( c ) )
-    throw( MCFException( "LoadDMX: error reading the input stream" ) );
-  while( c != '\n' );
-  }
-
- if( c != 'p' )
-  throw( MCFException( "LoadDMX: format error in the input stream" ) );
-
- char buf[ 3 ];    // free space
- DMXs >> buf;     // skip (in has to be "min")
-
- Index tn;
- if( ! ( DMXs >> tn ) )
-  throw( MCFException( "LoadDMX: error reading number of nodes" ) );
-
- Index tm;
- if( ! ( DMXs >> tm ) )
-  throw( MCFException( "LoadDMX: error reading number of arcs" ) );
-
- // allocate memory - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
- FRow      tU      = new FNumber[ tm ];  // arc upper capacities
- FRow      tDfct   = new FNumber[ tn ];  // node deficits
- Index_Set tStartn = new Index[ tm ];    // arc start nodes
- Index_Set tEndn   = new Index[ tm ];    // arc end nodes
- CRow      tC      = new CNumber[ tm ];  // arc costs
- CRow      tQ      = NULL;
- if( IsQuad )
-  tQ = new CNumber[ tm ];                // arc quadratic costs
-
-
- for( Index i = 0 ; i < tn ; )           // all deficits are 0
-  tDfct[ i++ ] = 0;                      // unless otherwise stated
-
- // read problem data - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
- Index i = 0;  // arc counter
- for(;;) {
-  if( ! ( DMXs >> c ) )
-   break;
-
-  switch( c ) {
-   case( 'c' ):  // comment line- - - - - - - - - - - - - - - - - - - - - - -
-    do
-     if( ! DMXs.get( c ) )
-      throw( MCFException( "LoadDMX: error reading the input stream" ) );
-    while( c != '\n' );
-    break;
-
-   case( 'n' ):  // description of a node - - - - - - - - - - - - - - - - - -
-    Index j;
-    if( ! ( DMXs >> j ) )
-     throw( MCFException( "LoadDMX: error reading node name" ) );
-
-    if( ( j < 1 ) || ( j > tn ) )
-     throw( MCFException( "LoadDMX: invalid node name" ) );
-
-    FNumber Dfctj;
-    if( ! ( DMXs >> Dfctj ) )
-     throw( MCFException( "LoadDMX: error reading deficit" ) );
-
-    tDfct[ j - 1 ] -= Dfctj;
-    break;
-
-   case( 'a' ):  // description of an arc - - - - - - - - - - - - - - - - - -
-    if( i == tm )
-     throw( MCFException( "LoadDMX: too many arc descriptors" ) );
-
-    if( ! ( DMXs >> tStartn[ i ] ) )
-     throw( MCFException( "LoadDMX: error reading start node" ) );
-
-    if( ( tStartn[ i ] < 1 ) || ( tStartn[ i ] > tn ) )
-     throw( MCFException( "LoadDMX: invalid start node" ) );
-
-    if( ! ( DMXs >> tEndn[ i ] ) )
-     throw( MCFException( "LoadDMX: error reading end node" ) );
-
-    if( ( tEndn[ i ] < 1 ) || ( tEndn[ i ] > tn ) )
-     throw( MCFException( "LoadDMX: invalid end node" ) );
-
-    if( tStartn[ i ] == tEndn[ i ] )
-     throw( MCFException( "LoadDMX: self-loops not permitted" ) );
-
-    FNumber LB;
-    if( ! ( DMXs >> LB ) )
-     throw( MCFException( "LoadDMX: error reading lower bound" ) );
-
-    if( ! ( DMXs >> tU[ i ] ) )
-     throw( MCFException( "LoadDMX: error reading upper bound" ) );
-
-    if( ! ( DMXs >> tC[ i ] ) )
-     throw( MCFException( "LoadDMX: error reading arc cost" ) );
-
-    if( tQ ) {
-     if( ! ( DMXs >> tQ[ i ] ) )
-      throw( MCFException( "LoadDMX: error reading arc quadratic cost" ) );
-
-     if( tQ[ i ] < 0 )
-      throw( MCFException( "LoadDMX: negative arc quadratic cost" ) );
-     }
-
-    if( tU[ i ] < LB )
-     throw( MCFException( "LoadDMX: lower bound > upper bound" ) );
-
-    tU[ i ] -= LB;
-    tDfct[ tStartn[ i ] - 1 ] += LB;
-    tDfct[ tEndn[ i ] - 1 ] -= LB;
-    #if( USENAME0 )
-     tStartn[ i ]--;  // in the DIMACS format, node names start from 1
-     tEndn[ i ]--;
-    #endif
-    i++;
-    break; 
-
-   default:  // invalid code- - - - - - - - - - - - - - - - - - - - - - - - -
-    throw( MCFException( "LoadDMX: invalid code" ) );
-
-   }  // end( switch( c ) )
-  }  // end( for( ever ) )
-
- if( i < tm )
-  throw( MCFException( "LoadDMX: too few arc descriptors" ) );
-
- // call LoadNet- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
- LoadNet( tn , tm , tn , tm , tU , tC , tDfct , tStartn , tEndn );
-
- // then pass quadratic costs, if any
-
- if( tQ )
-  ChgQCoef( tQ );
-
- // delete the original data structures - - - - - - - - - - - - - - - - - - -
- // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
- delete[] tQ;
- delete[] tC;
- delete[] tEndn;
- delete[] tStartn;
- delete[] tDfct;
- delete[] tU;
-
- }  // end( MCFClass::LoadDMX )
-
-/*--------------------------------------------------------------------------*/
-
-inline void MCFClass::SetPar( int par , int val )
-{
- switch( par ) {
- case( kMaxIter ):
-  MaxIter = val;
-  break;
-
- case( kReopt ):
-  Senstv = (val == kYes);
-  break;
-  
- default:
-  throw( MCFException( "Error using SetPar: unknown parameter" ) );
-  }
- }
-
-/*--------------------------------------------------------------------------*/
-
-inline void MCFClass::SetPar( int par, double val )
-{
- switch( par ) {
- case( kEpsFlw ):
-  if( EpsFlw != FNumber( val ) ) {
-   EpsFlw = FNumber( val );
-   EpsDfct = EpsFlw * ( nmax ? nmax : 100 );
-   status = kUnSolved;
-   }
-  break;
-
- case( kEpsDfct ):
-  if( EpsDfct != FNumber( val ) ) {
-   EpsDfct = FNumber( val );
-   status = kUnSolved;
-   }
-  break;
-
- case( kEpsCst ):
-  if( EpsCst != CNumber( val ) ) {
-   EpsCst = CNumber( val );
-   status = kUnSolved;
-   }
-  break;
-  
- case( kMaxTime ):
-  MaxTime = val;
-  break;
-
- default:
-  throw( MCFException( "Error using SetPar: unknown parameter" ) );
-  }
- }
-
-/*--------------------------------------------------------------------------*/
-
-inline void MCFClass::GetPar( int par, int &val )
-{
- switch( par ) {
- case( kMaxIter ):
-  val = MaxIter;
-  break;
-
- case( kReopt ):
-  if( Senstv ) val = kYes;
-  else         val = kNo;
-  break;
-  
- default:
-  throw( MCFException( "Error using GetPar: unknown parameter" ) );
-  }
- }
-
-/*--------------------------------------------------------------------------*/
-
-inline void MCFClass::GetPar( int par , double &val )
-{
- switch( par ) {
- case( kEpsFlw ):
-  val = double( EpsFlw );
-  break;
-
- case( kEpsDfct ):
-  val = double( EpsDfct );
-  break;
-
- case( kEpsCst ):
-  val = double( EpsCst );
-  break;
-  
- case( kMaxTime ):
-  val = MaxTime;
-  break;
-
- default:
-  throw( MCFException( "Error using GetPar: unknown parameter" ) );
-  }
- }
-
-/*--------------------------------------------------------------------------*/
-
-inline void MCFClass::CheckPSol( void )
-{
- FRow tB = new FNumber[ MCFn() ];
- MCFDfcts( tB );
- #if( ! USENAME0 )
-  tB--;
- #endif
-
- FRow tX = new FNumber[ MCFm() ];
- MCFGetX( tX );
- FONumber CX = 0;
- for( Index i = 0 ; i < MCFm() ; i++ )
-  if( GTZ( tX[ i ] , EpsFlw ) ) {
-   if( IsClosedArc( i ) )
-    throw( MCFException( "Closed arc with nonzero flow (CheckPSol)" ) );
-
-   if( GT( tX[ i ] , MCFUCap( i ) , EpsFlw ) )
-    throw( MCFException( "Arc flow exceeds capacity (CheckPSol)" ) );
-
-   if( LTZ( tX[ i ] , EpsFlw ) )
-    throw( MCFException( "Arc flow is negative (CheckPSol)" ) );
-
-   CX += ( FONumber( MCFCost( i ) ) +
-           FONumber( MCFQCoef( i ) ) * FONumber( tX[ i ] ) / 2 )
-         * FONumber( tX[ i ] );
-   tB[ MCFSNde( i ) ] += tX[ i ];
-   tB[ MCFENde( i ) ] -= tX[ i ];
-   }
- delete[] tX;
- #if( ! USENAME0 )
-  tB++;
- #endif
-
- for( Index i = 0 ; i < MCFn() ; i++ )
-  if( ! ETZ( tB[ i ] , EpsDfct ) )
-   throw( MCFException( "Node is not balanced (CheckPSol)" ) );
-
- delete[] tB;
- CX -= MCFGetFO(); 
- if( ( CX >= 0 ? CX : - CX ) > EpsCst * MCFn() )
-  throw( MCFException( "Objective function value is wrong (CheckPSol)" ) );
-
- }  // end( CheckPSol )
-
-/*--------------------------------------------------------------------------*/
-
-inline void MCFClass::CheckDSol( void )
-{
- CRow tPi = new CNumber[ MCFn() ];
- MCFGetPi( tPi );
-
- FONumber BY = 0;
- for( Index i = 0 ; i < MCFn() ; i++ )
-  BY += FONumber( tPi[ i ] ) * FONumber( MCFDfct( i ) );
-
- CRow tRC = new CNumber[ MCFm() ];
- MCFGetRC( tRC );
-
- FRow tX = new FNumber[ MCFm() ];
- MCFGetX( tX );
-
- #if( ! USENAME0 )
-  tPi--;
- #endif
-
- FONumber QdrtcCmpnt = 0;  // the quadratic part of the objective
- for( Index i = 0 ; i < MCFm() ; i++ ) {
-  if( IsClosedArc( i ) )
-   continue;
-
-  cFONumber tXi = FONumber( tX[ i ] );
-  cCNumber QCi = CNumber( MCFQCoef( i ) * tXi );
-  QdrtcCmpnt += QCi * tXi;
-  CNumber RCi = MCFCost( i ) + QCi
-                + tPi[ MCFSNde( i ) ] - tPi[ MCFENde( i ) ];
-  RCi -= tRC[ i ];
-
-  if( ! ETZ( RCi , EpsCst ) )
-   throw( MCFException( "Reduced Cost value is wrong (CheckDSol)" ) );
-
-  if( LTZ( tRC[ i ] , EpsCst ) ) {
-   BY += FONumber( tRC[ i ] ) * FONumber( MCFUCap( i ) );
-
-   if( LT( tX[ i ] , MCFUCap( i ) , EpsFlw ) )
-    throw( MCFException( "Csomplementary Slackness violated (CheckDSol)" ) );
-   }
-  else
-   if( GTZ( tRC[ i ] , EpsCst ) && GTZ( tX[ i ] , EpsFlw ) )
-    throw( MCFException( "Complementary Slackness violated (CheckDSol)" ) );
-
-  }  // end( for( i ) )
-
- delete[] tX;
- delete[] tRC;
- #if( ! USENAME0 )
-  tPi++;
- #endif
- delete[] tPi;
-
- BY -= ( MCFGetFO() + QdrtcCmpnt / 2 );
- if( ( BY >= 0 ? BY : - BY ) > EpsCst * MCFn() )
-  throw( MCFException( "Objective function value is wrong (CheckDSol)" ) );
-
- }  // end( CheckDSol )
-
-/*--------------------------------------------------------------------------*/
-
-inline void MCFClass::WriteMCF( ostream &oStrm , int frmt )
-{
- if( ( ! numeric_limits<FNumber>::is_integer ) ||
-     ( ! numeric_limits<CNumber>::is_integer ) )
-  oStrm.precision( 12 );
-
- switch( frmt ) {
-  case( kDimacs ):  // DIMACS standard format - - - - - - - - - - - - - - - -
-                    //- - - - - - - - - - - - - - - - - - - - - - - - - - - -
-  case( kQDimacs ):  // ... linear or quadratic
-
-   // compute and output preamble and size- - - - - - - - - - - - - - - - - -
-
-   oStrm << "p min " << MCFn() << " ";
-   {
-    Index tm = MCFm();
-    for(  Index i = MCFm() ; i-- ; )
-     if( IsClosedArc( i ) || IsDeletedArc( i ) )
-      tm--;
-
-    oStrm << tm << endl;
-    }
-
-   // output arc information- - - - - - - - - - - - - - - - - - - - - - - - -
-
-   for(  Index i = 0 ; i < MCFm() ; i++ )
-    if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) {
-     oStrm << "a\t";
-     #if( USENAME0 )
-      oStrm << MCFSNde( i ) + 1 << "\t" << MCFENde( i ) + 1 << "\t";
-     #else
-      oStrm << MCFSNde( i ) << "\t" << MCFENde( i ) << "\t";
-     #endif
-     oStrm << "0\t" << MCFUCap( i ) << "\t" << MCFCost( i );
-
-     if( frmt == kQDimacs )
-      oStrm << "\t" << MCFQCoef( i );
-
-     oStrm << endl;
-     }
-
-   // output node information - - - - - - - - - - - - - - - - - - - - - - - -
-
-   for(  Index i = 0 ; i < MCFn() ; ) {
-    cFNumber Dfcti = MCFDfct( i++ );
-    if( Dfcti )
-     oStrm << "n\t" << i << "\t" << - Dfcti << endl;
-    }
-
-   break;
-
-  case( kMPS ):     // MPS format(s)- - - - - - - - - - - - - - - - - - - - -
-  case( kFWMPS ):   //- - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
-   // writing problem name- - - - - - - - - - - - - - - - - - - - - - - - - -
-
-   oStrm << "NAME      MCF" << endl;
-
-   // writing ROWS section: flow balance constraints- - - - - - - - - - - - -
-
-   oStrm << "ROWS" << endl << " N  obj" << endl;
-
-   for(  Index i = 0 ; i < MCFn() ; )
-    oStrm << " E  c" << i++ << endl;
-
-   // writing COLUMNS section - - - - - - - - - - - - - - - - - - - - - - - -
-
-   oStrm << "COLUMNS" << endl;
-
-   if( frmt == kMPS ) {  // "modern" MPS
-    for(  Index i = 0 ; i < MCFm() ; i++ )
-     if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) {
-      oStrm << " x" << i << "\tobj\t" << MCFCost( i ) << "\tc";
-      #if( USENAME0 )
-       oStrm << MCFSNde( i ) << "\t-1" << endl << " x" << i << "\tc"
-             << MCFENde( i ) << "\t1" << endl;
-      #else
-       oStrm << MCFSNde( i ) - 1 << "\t-1" << endl << " x" << i << "\tc"
-             << MCFENde( i ) - 1 << "\t1" << endl;
-      #endif
-      }
-     }
-   else              // "old" MPS
-    for(  Index i = 0 ; i < MCFm() ; i++ ) {
-     ostringstream myField;
-     myField << "x" << i << ends;
-     oStrm << "    " << setw( 8 )  << left << myField.str()
-           << "  "   << setw( 8 )  << left << "obj"
-           << "  "   << setw( 12 ) << left << MCFCost( i );
-
-     myField.seekp( 0 );
-     myField << "c" << MCFSNde( i ) - ( 1 - USENAME0 ) << ends;
-
-     oStrm << "   " << setw( 8 )  << left << myField.str()
-           << "  "  << setw( 12 ) << left << -1 << endl;
-
-     myField.seekp( 0 );
-     myField << "x" << i << ends;
-
-     oStrm << "    " << setw( 8 ) << left << myField.str();
-
-     myField.seekp( 0 );
-     myField << "c" << MCFENde( i ) - ( 1 - USENAME0 ) << ends;
-
-     oStrm << "  " << setw( 8 ) << left << myField.str() << endl;
-     }
-
-   // writing RHS section: flow balance constraints - - - - - - - - - - - - -
-
-   oStrm << "RHS" << endl;
-
-   if( frmt == kMPS )  // "modern" MPS
-    for( Index i = 0 ; i < MCFn() ; i++ )
-     oStrm << " rhs\tc" << i << "\t" << MCFDfct( i ) << endl;
-   else              // "old" MPS
-    for( Index i = 0 ; i < MCFn() ; i++ ) {
-     ostringstream myField;
-     oStrm << "    " << setw( 8 ) << left << "rhs";
-     myField << "c" << i << ends;
-
-     oStrm << "  " << setw( 8 )  << left << myField.str()
-           << "  " << setw( 12 ) << left << MCFDfct( i ) << endl;
-     }
-
-   // writing BOUNDS section- - - - - - - - - - - - - - - - - - - - - - - - -
-
-   oStrm << "BOUNDS" << endl;
-
-   if( frmt == kMPS ) {  // "modern" MPS
-    for( Index i = 0 ; i < MCFm() ; i++ )
-     if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) )
-      oStrm << " UP BOUND\tx" << i << "\t" << MCFUCap( i ) << endl;
-    }
-   else              // "old" MPS
-    for(  Index i = 0 ; i < MCFm() ; i++ )
-     if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) {
-      ostringstream myField;
-      oStrm << " " << setw( 2 ) << left << "UP"
-            << " " << setw( 8 ) << left << "BOUND";
-
-      myField << "x" << i << ends;
-
-      oStrm << "  " << setw( 8 )  << left << myField.str()
-            << "  " << setw( 12 ) << left << MCFUCap( i ) << endl;
-      }
-
-   oStrm << "ENDATA" << endl;
-   break;
-
-  default:          // unknown format - - - - - - - - - - - - - - - - - - - -
-                    //- - - - - - - - - - - - - - - - - - - - - - - - - - - -
-   oStrm << "Error: unknown format " << frmt << endl;
-  }
- }  // end( MCFClass::WriteMCF )
-
-/*--------------------------------------------------------------------------*/
-
-#if( OPT_USE_NAMESPACES )
- };  // end( namespace MCFClass_di_unipi_it )
-#endif
-
-/*--------------------------------------------------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#endif  /* MCFClass.h included */
-
-/*--------------------------------------------------------------------------*/
-/*------------------------- End File MCFClass.h ----------------------------*/
-/*--------------------------------------------------------------------------*/
diff --git a/clang/lib/Analysis/MCFSimplex.cpp b/clang/lib/Analysis/MCFSimplex.cpp
index 4b1b75a..3da685d 100644
--- a/clang/lib/Analysis/MCFSimplex.cpp
+++ b/clang/lib/Analysis/MCFSimplex.cpp
@@ -30,7 +30,7 @@
 /*------------------------------ INCLUDES ----------------------------------*/
 /*--------------------------------------------------------------------------*/
 
-#include "MCFSimplex.h"
+#include "clang/Analysis/Analyses/IntervalSolver/MCFSimplex.h"
 #include <algorithm>
 #include <iostream>
 
@@ -116,9 +116,10 @@ using namespace MCFClass_di_unipi_it;
    optimize the update of the potential.
    Unfortunately this doesn't work well: for this reason it is set to 0. */
 
-
+#ifndef throw
 // QUICK HACK
 #define throw(x) assert(false);
+#endif
 
 /*--------------------------------------------------------------------------*/
 /*--------------------------- FUNCTIONS ------------------------------------*/
@@ -410,7 +411,7 @@ void MCFSimplex::SetAlg( bool UsPrml , char WhchPrc )
      MemDeAlloc(false);
      if( Senstv && ( status != kUnSolved ) ) {
       nodePType *node = dummyRootP;
-      for( int i = 0 ; i < n ; i++ )
+      for( unsigned int i = 0 ; i < n ; i++ )
        node = node->nextInT;
 
       node->nextInT = NULL;
@@ -483,7 +484,7 @@ void MCFSimplex::SetAlg( bool UsPrml , char WhchPrc )
     CreateAdditionalDualStructures();
     MemDeAlloc(true);
     nodeDType *node = dummyRootD;
-    for( int i = 0 ; i < n ; i++ )
+    for( unsigned int i = 0 ; i < n ; i++ )
      node = node->nextInT;
 
     node->nextInT = NULL;
@@ -513,11 +514,12 @@ void MCFSimplex::SetAlg( bool UsPrml , char WhchPrc )
    //#endif
    }
 #endif
- if( newPricingRule == kFirstEligibleArc )
+ if( newPricingRule == kFirstEligibleArc ) {
   if( newUsePrimalSimplex )
    arcToStartP = arcsP;
   else
    arcToStartD = arcsD;
+ }
 
  if( ( nmax && mmax ) && ( newPricingRule == kCandidateListPivot ) )
   MemAllocCandidateList();
@@ -941,13 +943,13 @@ void MCFSimplex::ChgCosts( cCRow NCost , cIndex_Set nms ,
    for( arcDType *arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ )
     arc->cost = *(NCost++); 
 
- if( Senstv && ( status != kUnSolved ) )
+ if( Senstv && ( status != kUnSolved ) ) {
   if( usePrimalSimplex )
    ComputePotential( dummyRootP );
   else {
    #if( QUADRATICCOST == 0 )
     for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ )
-     if( arc->ident > BASIC )
+     if( arc->ident > BASIC ) {
       if( GTZ( ReductCost( arc ) , EpsCst ) ) {
        arc->flow = 0;
        arc->ident = AT_LOWER;
@@ -956,13 +958,15 @@ void MCFSimplex::ChgCosts( cCRow NCost , cIndex_Set nms ,
        arc->flow = arc->upper; 
        arc->ident = AT_UPPER;
        }
+     }
 
     CreateInitialDModifiedBalanceVector();
     PostDVisit( dummyRootD );
    #endif
    }
- else
+ } else {
   status = kUnSolved;
+ }
 
  }  // end( MCFSimplex::ChgCosts )
 
@@ -993,8 +997,8 @@ void MCFSimplex::ChgCost( Index arc , cCNumber NCost )
      node = ( arcsD + arc )->head;
 
     ComputePotential( dummyRootD );
-    for( arcDType *a = arcsD ; a != stopArcsD ; a++)
-     if( a->ident > BASIC )
+    for( arcDType *a = arcsD ; a != stopArcsD ; a++) 
+     if( a->ident > BASIC ) {
       if( GTZ( ReductCost( a ) , EpsCst ) ) {
        a->flow = 0;
        a->ident = AT_LOWER;
@@ -1003,6 +1007,7 @@ void MCFSimplex::ChgCost( Index arc , cCNumber NCost )
        a->flow = a->upper; 
        a->ident = AT_UPPER;
        }
+     }
 
     CreateInitialDModifiedBalanceVector();
     PostDVisit( dummyRootD );
@@ -1016,8 +1021,8 @@ void MCFSimplex::ChgCost( Index arc , cCNumber NCost )
 
 /*-------------------------------------------------------------------------*/
 
-void MCFSimplex::ChgQCoef( cCRow NQCoef , cIndex_Set nms ,
-			   cIndex strt , Index stp )
+void MCFSimplex::ChgQCoef( cCRow NQCoef , cIndex_Set ,
+			   cIndex , Index stp )
 {
  if( stp > m )
   stp = m;
@@ -1057,7 +1062,7 @@ void MCFSimplex::ChgQCoef( cCRow NQCoef , cIndex_Set nms ,
 
 /*-------------------------------------------------------------------------*/
 
-void MCFSimplex::ChgQCoef( Index arc , cCNumber NQCoef )
+void MCFSimplex::ChgQCoef( Index , cCNumber NQCoef )
 {
  #if( QUADRATICCOST )
   if( arc >= m )
@@ -1362,7 +1367,7 @@ void MCFSimplex::CloseArc( cIndex name )
      ComputePotential( dummyRootD );
 
      for( arcDType *a = arcsD ; a != stopArcsD ; a++ )
-      if( a->ident > BASIC )
+      if( a->ident > BASIC ) {
        if( GTZ( ReductCost( a ) , EpsCst ) ) {
 	a->flow = 0;
 	a->ident = AT_LOWER;
@@ -1371,6 +1376,7 @@ void MCFSimplex::CloseArc( cIndex name )
 	a->flow = a->upper; 
 	a->ident = AT_UPPER;
         }
+      }
      }
 
     CreateInitialDModifiedBalanceVector();
@@ -2564,11 +2570,12 @@ void MCFSimplex::DualSimplex( void )
     while( fine == false ) {
      /* If node is the root of subtree T2, Dual Simplex jumps to the node
 	(if exists) which follows the last node of T2 */
-     if( node == h2 )
+     if( node == h2 ) {
       if( lastNodeOfT2->nextInT )
        node = lastNodeOfT2->nextInT;
       else
        break;
+     }
 
      // Search arc in the Backward Star of nodes of T1
      arcDType *arc = node->firstBs;
@@ -2760,7 +2767,6 @@ void MCFSimplex::DualSimplex( void )
      theta = leavingArc->flow - leavingArc->upper;
     // Initial value of theta is the infeasibility of the leaving arc
 
-    FNumber t;
     nodeDType *k1;
     nodeDType *k2;
     /* if entering arc is in U, Dual Simplex pushs flow in the cycle 
@@ -2905,7 +2911,7 @@ void MCFSimplex::DualSimplex( void )
 /*--------------------------------------------------------------------------*/
 
 template<class N, class A>
-void MCFSimplex::UpdateT( A *h , A *k , N *h1 , N *h2 , N *k1 , N *k2 )
+void MCFSimplex::UpdateT( A * , A *k , N * , N *h2 , N *k1 , N *k2 )
 {
  /* In subtree T2 there is a path from node h2 (deepest node of the leaving
     arc h and root of T2) to node k2 (deepest node of the leaving arc h and
@@ -3421,7 +3427,7 @@ MCFSimplex::arcDType* MCFSimplex::RuleDualCandidateListPivot( void )
  // Search other arcs to fill the list
  do {
   nodeDType *node = nodesD + groupPos;
-  for( node ; node < stopNodesD ; node += numGroup ) {
+  for( ; node < stopNodesD ; node += numGroup ) {
    arcDType *arc = node->enteringTArc;
    cFNumber flow = arc->flow;
    if( LTZ( flow , EpsFlw ) ) {
@@ -3598,18 +3604,19 @@ void MCFSimplex::PostPVisit( nodePType *r )
  // The method controls if "r" is a leaf in T
  bool rLeaf = false;
  int i = r - nodesP;
- if( r->nextInT ) 
+ if( r->nextInT ) {
   if( ( r->nextInT )->subTreeLevel <= r->subTreeLevel )
    rLeaf = true;
   else
    rLeaf = true;
+ }
 
- if( rLeaf )  // If "r" is a leaf
+ if( rLeaf ) { // If "r" is a leaf
   if( ( r->enteringTArc)->head == r ) // If enteringTArc of "r" goes in "r"
    ( r->enteringTArc )->flow = modifiedBalance[ i ];
   else // If enteringTArc of "r" goes out "r"
    ( r->enteringTArc )->flow = - modifiedBalance[ i ];
- else { // If "r" isn't a leaf
+ } else { // If "r" isn't a leaf
   nodePType *desc = r->nextInT;
   // Call PostPVisit for every child of "r"
   while( ( desc ) && ( desc->subTreeLevel > r->subTreeLevel ) ) {
@@ -3624,12 +3631,13 @@ void MCFSimplex::PostPVisit( nodePType *r )
    desc = desc->nextInT;
    }
 
-  if( r != dummyRootP )
+  if( r != dummyRootP ) {
    if( ( r->enteringTArc )->head == r ) // If enteringTArc of "r" goes in "r"
     ( r->enteringTArc )->flow = modifiedBalance[ i ];
    else // If enteringTArc of "r" goes out "r"
     ( r->enteringTArc )->flow = - modifiedBalance[ i ];
   }
+  }
  }
 
 /*--------------------------------------------------------------------------*/
@@ -3650,11 +3658,12 @@ void MCFSimplex::BalanceFlow( nodePType *r )
  else {
   // The method controls if "r" is a leaf in T
   bool rLeaf = false;
-  if( r->nextInT )
+  if( r->nextInT ) {
    if( ( r->nextInT )->subTreeLevel <= r->subTreeLevel )
     rLeaf = true;
    else
     rLeaf = true;
+  }
 
   if( rLeaf )  // If "r" is a leaf
    AdjustFlow( r );  // The method controls if entering basic arc in "r" is
@@ -3680,7 +3689,7 @@ void MCFSimplex::BalanceFlow( nodePType *r )
 void MCFSimplex::AdjustFlow( nodePType *r )
 {
  arcPType *arc = r->enteringTArc;
- if( arc >= dummyArcsP ) // If entering arc of "r" is a dummy arc
+ if( arc >= dummyArcsP ) { // If entering arc of "r" is a dummy arc
   if( LTZ( arc->flow , EpsFlw ) ) {
    // If this dummy arc has flow < 0, the algorithm overturns the arc
    nodePType *temp = arc->tail;
@@ -3749,6 +3758,7 @@ void MCFSimplex::AdjustFlow( nodePType *r )
    }
   }
  }
+ }
 
 /*--------------------------------------------------------------------------*/
 
@@ -3792,11 +3802,12 @@ void MCFSimplex::PostDVisit( nodeDType *r )
   // The method controls if "r" is a leaf in T
   bool rLeaf = false;
   int i = r - nodesD;
-  if( r->nextInT )
+  if( r->nextInT ) {
    if( ( r->nextInT )->subTreeLevel <= r->subTreeLevel )
     rLeaf = true;
    else
     rLeaf = true;
+  }
 
   if( rLeaf ) // If "r" is a leaf
    if( ( r->enteringTArc)->head == r ) // If enteringTArc of "r" goes in "r"
@@ -3819,12 +3830,13 @@ void MCFSimplex::PostDVisit( nodeDType *r )
     desc = desc->nextInT;
     }
 
-   if( r != dummyRootD )
+   if( r != dummyRootD ) {
     if( ( r->enteringTArc )->head == r ) // If enteringTArc of "r" goes in "r"
      ( r->enteringTArc )->flow = modifiedBalance[ i ];
     else // If enteringTArc of "r" goes out "r"
      ( r->enteringTArc )->flow = - modifiedBalance[ i ];
    }
+   }
  #endif
  }
 
@@ -3961,7 +3973,7 @@ void MCFSimplex::PrintDArc( arcDType *arc )
 
 MCFSimplex::nodePType* MCFSimplex::RecoverPNode( Index ind ) 
 {
- if( ( ind < 0 ) || ( ind > n ) )
+ if( ( ind > n ) )
   return( NULL );
  if( ind )
   return( nodesP + ind - 1 );
@@ -3993,7 +4005,7 @@ MCFSimplex::arcPType* MCFSimplex::RecoverPArc( nodePType *tail ,
 
 MCFSimplex::nodeDType* MCFSimplex::RecoverDNode( Index ind )
 {
- if( ( ind < 0 ) || ( ind > n ) ) 
+ if( ( ind > n ) ) 
   return( NULL );
 
  if( ind )
diff --git a/clang/lib/Analysis/MCFSimplex.h b/clang/lib/Analysis/MCFSimplex.h
deleted file mode 100644
index a8d84f6..0000000
--- a/clang/lib/Analysis/MCFSimplex.h
+++ /dev/null
@@ -1,1086 +0,0 @@
-/*--------------------------------------------------------------------------*/
-/*------------------------- File MCFSimplex.h ------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @file
- * Linear and Quadratic Min Cost Flow problems solver based on the (primal and
- * dual) simplex algorithm. Conforms to the standard MCF interface defined in
- * MCFClass.h.
- *
- * \Version 1.00
- *
- * \date 29 - 08 - 2011
- *
- * \author Alessandro Bertolini \n
- *         Antonio Frangioni \n
- *         Operations Research Group \n
- *         Dipartimento di Informatica \n
- *         Universita' di Pisa \n
- *
- * Copyright &copy 2008 - 2011 by Alessandro Bertolini, Antonio Frangioni
- */
-/*--------------------------------------------------------------------------*/
-/*--------------------------------------------------------------------------*/
-/*----------------------------- DEFINITIONS --------------------------------*/
-/*--------------------------------------------------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#ifndef __MCFSimplex
- #define __MCFSimplex  /* self-identification: #endif at the end of the file */
-
-/*--------------------------------------------------------------------------*/
-/*------------------------------ INCLUDES ----------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#include "MCFClass.h"
-
-/*@} -----------------------------------------------------------------------*/
-/*--------------------------- NAMESPACE ------------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#if( OPT_USE_NAMESPACES )
-namespace MCFClass_di_unipi_it
-{
-#endif
-
-/*--------------------------------------------------------------------------*/
-/*-------------------------------- MACROS ----------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @defgroup MCFSimplex_MACROS Compile-time switches in MCFSimplex.h
-    There is only one macro in MCFSimplex, but it is very important!
-    @{ */
-
-#define QUADRATICCOST 0
-
-/**< Setting QUADRATICCOST == 1 the solver can solve problems with linear and 
-   quadratic costs too (but the latter only with the Primal Simplex).
-   The reason for having a macro is that when quadratic costs are present the
-   "arcType" struct has the additional field "quadraticCost" to hold it.
-   Furthermore, the field "ident" is not created because the solver doesn't
-   use the classical TLU tripartition. Instead, closed arcs and deleted arcs
-   are characterized as follows:
-   - closed arcs have the field "cost" to INFINITY (Inf<FNumber>());
-   - deleted arcs have the field "upper" to INFINITY and the "tail" and
-     "head" field are NULL.
-   Furthermore, the solver needs the variables "ignoredEnteringArc" and
-   "firstIgnoredEnteringArc", used to avoid nasty loops during the execution
-   of the Quadratic Primal Simplex algorithm.
-   If, instead, QUADRATICCOST == 0 then the solver can solve only problems
-   with linear costs. Hence, the field "quadraticCost" is useless and it
-   isn't created. Furthermore, Primal Simplex and Dual Simplex use the
-   tripartition TLU to divide the arcs, so the solver creates the field
-   "ident", which differentiates the set of the arcs in among deleted arcs,
-   closed arcs, arcs in T, arcs in L, arcs in U.
-   Thus, with QUADRATICCOST == 0 the solver cannot solve problems with
-   quadratic costs, but it does solve problems with linear costs faster. */
-
-/*@}  end( group( MCFCLASS_MACROS ) ) */ 
-
-/*--------------------------------------------------------------------------*/
-/*---------------------------- CLASSES -------------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @defgroup MCFSimplex_CLASSES Classes in MCFSimplex.h
-    @{ */
-
-/** The MCFSimplex class derives from the abstract base class MCFClass, thus
-    sharing its (standard) interface, and implements both the Primal and
-    Dual network simplex algorithms for solving (Linear and Quadratic) 
-    Min Cost Flow problems */
-
-class MCFSimplex: public MCFClass 
-{
-
-/*--------------------------------------------------------------------------*/
-/*----------------------- PUBLIC PART OF THE CLASS -------------------------*/
-/*--------------------------------------------------------------------------*/
-/*--                                                                      --*/
-/*--  The following methods and data are the actual interface of the      --*/
-/*--  class: the standard user should use these methods and data only.    --*/
-/*--                                                                      --*/
-/*--------------------------------------------------------------------------*/
-
- public:
-
-/*--------------------------------------------------------------------------*/
-/*---------------------------- PUBLIC TYPES --------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-/** Public enum describing the parameters of MCFSimplex. */
-
- enum SimplexParam 
- { 
-  kAlgPrimal = kLastParam , ///< parameter to set algorithm (Primal/Dual):
-  kAlgPricing ,             ///< parameter to set algorithm of pricing
-  kNumCandList ,            /**< parameter to set the number of candidate
-			         list for Candidate List Pivot method */
-  kHotListSize ,            /**< parameter to set the size of Hot List
-			         for Candidate List Pivot method */
-  kRecomputeFOLimits ,      /**< parameter to set the number of iterations
-                                 in which quadratic Primal Simplex computes 
-                                 "manually" the f.o. value */
-  kEpsOpt                   /**< parameter to set the precision of the 
-                                 objective function value for the 
-				 quadratic Primal Simplex */
-  };
-    
-/** Public enum describing the pricing rules in MCFSimplex::SetAlg(). */
-
- enum enumPrcngRl 
- { 
-  kDantzig = 0,        ///< Dantzig's rule (most violated constraint)
-  kFirstEligibleArc ,  ///< First eligible arc in round-robin
-  kCandidateListPivot  ///< Candidate List Pivot Rule
-  };
-
-/*--------------------------------------------------------------------------*/
-/*--------------------------- PUBLIC METHODS -------------------------------*/
-/*--------------------------------------------------------------------------*/
-/*---------------------------- CONSTRUCTOR ---------------------------------*/
-/*--------------------------------------------------------------------------*/
-
- MCFSimplex( cIndex nmx = 0 , cIndex mmx = 0 );
-
-/**< Constructor of the class, as in MCFClass::MCFClass(). */
-
-/*--------------------------------------------------------------------------*/
-/*-------------------------- OTHER INITIALIZATIONS -------------------------*/
-/*--------------------------------------------------------------------------*/
-
- void LoadNet( cIndex nmx = 0 , cIndex mmx = 0 , cIndex pn = 0 ,
-	       cIndex pm = 0 , cFRow pU = NULL , cCRow pC = NULL ,
-	       cFRow pDfct = NULL , cIndex_Set pSn = NULL ,
-	       cIndex_Set pEn = NULL );
-
-/**< Inputs a new network, as in MCFClass::LoadNet(). */
-
-/*--------------------------------------------------------------------------*/
-
- void SetAlg( bool UsPrml , char WhchPrc );
-
-/**< Selects which algorithm (Primal vs Dual Network Simplex), and which
-   pricing rule within the algorithm, is used.
-
-   If UsPrml == TRUE then the Primal Network Simplex algorithm is used,
-   otherwise the Dual Network Simplex is used.
-
-   The allowed values for WhchPrc are:
-
-   - kDantzig            Dantzig's pricing rule, i.e., most violated dual 
-                         constraint; this can only be used with the Primal 
-                         Network Simplex
-
-   - kFirstEligibleArcA  "dumb" rule, first eligible arc in round-robin;
-
-   - kCandidateListPivot Candidate List Pivot Rule
-
-   If this method is *not* called, the Primal Network Simplex with the
-   Candidate List Pivot Rule (the best setting on most instances) is
-   used. */
-
-/*--------------------------------------------------------------------------*/
-
- void SetPar( int par , int val );
-
-/**< Set general integer parameters.
-
-   @param par   is the parameter to set; since this method accepts an int
-                value, the enum SimplexParam can be used in addition to the
-                enum MCFParam to specify the integer parameters (every enum
-		is an int).
-
-   @param value is the value to assign to the parameter. */
-
-/*-------------------------------------------------------------------------*/
-
- void SetPar( int par , double val );
-
-/**< Set general float parameters.
-
-   @param par   is the parameter to set; since this method accepts an int
-                value, the enum SimplexParam can be used in addition to the
-                enum MCFParam to specify the float parameters (every enum
-		is an int).
-
-   @param value is the value to assign to the parameter. */
-
-/*--------------------------------------------------------------------------*/
-/*-------------------- METHODS FOR SOLVING THE PROBLEM ---------------------*/
-/*--------------------------------------------------------------------------*/
-
- void SolveMCF( void );
-
-/*--------------------------------------------------------------------------*/
-/*---------------------- METHODS FOR READING RESULTS -----------------------*/
-/*--------------------------------------------------------------------------*/
-
- void MCFGetX( FRow F , Index_Set nms = NULL ,
-	       cIndex strt = 0 , Index stp = Inf<Index>() );
-
-/*--------------------------------------------------------------------------*/
-
- void MCFGetRC( CRow CR , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() ) ;
-
- CNumber MCFGetRC( cIndex i );
-
-/*--------------------------------------------------------------------------*/
-
- void MCFGetPi( CRow P , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() );
-
-/*--------------------------------------------------------------------------*/
-
- FONumber MCFGetFO( void );
-
-/*--------------------------------------------------------------------------*/
-/*-------------- METHODS FOR READING THE DATA OF THE PROBLEM ---------------*/
-/*--------------------------------------------------------------------------*/
-
- virtual void MCFArcs( Index_Set Startv , Index_Set Endv ,
-		       cIndex_Set nms = NULL , cIndex strt = 0 ,
-		       Index stp = Inf<Index>() );
-
- inline Index MCFSNde( cIndex i );
-
- inline Index MCFENde( cIndex i );
-
-/*--------------------------------------------------------------------------*/
-  
- void MCFCosts( CRow Costv , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() );
-
- inline CNumber MCFCost( cIndex i );
-
-/*--------------------------------------------------------------------------*/
-
- void MCFQCoef( CRow Qv , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() );
-
- inline CNumber MCFQCoef( cIndex i );
-
-/*--------------------------------------------------------------------------*/
-
- void MCFUCaps( FRow UCapv , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() );
-
- inline FNumber MCFUCap( cIndex i );
-
-/*--------------------------------------------------------------------------*/
-
- void MCFDfcts( FRow Dfctv , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() );
-
- inline FNumber MCFDfct( cIndex i );
-
-/*--------------------------------------------------------------------------*/
-/*------------- METHODS FOR ADDING / REMOVING / CHANGING DATA --------------*/
-/*--------------------------------------------------------------------------*/
-/*------- Changing the costs, QCoef, deficits and upper capacities ---------*/
-/*--------------------------------------------------------------------------*/
-
- void ChgCosts( cCRow NCost , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() );
-
- void ChgCost( Index arc , cCNumber NCost );
-
-/*--------------------------------------------------------------------------*/
-
- void ChgQCoef( cCRow NQCoef = NULL , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() );
-
- void ChgQCoef( Index arc , cCNumber NQCoef );
-
-/*--------------------------------------------------------------------------*/
-
- void ChgDfcts( cFRow NDfct , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() );
-
- void ChgDfct( Index nod , cFNumber NDfct );
-
-/*--------------------------------------------------------------------------*/
-
- void ChgUCaps( cFRow NCap , cIndex_Set nms = NULL ,
-		cIndex strt = 0 , Index stp = Inf<Index>() );
-
- void ChgUCap( Index arc , cFNumber NCap );
-
-/*--------------------------------------------------------------------------*/
-/*--------------- Modifying the structure of the graph ---------------------*/
-/*--------------------------------------------------------------------------*/
-
- void CloseArc( cIndex name );
-
- void DelNode( cIndex name );
-
- bool IsClosedArc( cIndex name );
-
- void OpenArc( cIndex name ) ;
-
- Index AddNode( cFNumber aDfct );
-
- void ChangeArc( cIndex name , cIndex nSS = Inf<Index>() ,
-		 cIndex nEN = Inf<Index>() );
-
- void DelArc( cIndex name );
-
- Index AddArc( cIndex Start , cIndex End , cFNumber aU , cCNumber aC );
-
- bool IsDeletedArc( cIndex name );
-
-/*--------------------------------------------------------------------------*/
-/*------------------------------ DESTRUCTOR --------------------------------*/
-/*--------------------------------------------------------------------------*/
-
- ~MCFSimplex();
-
-/*--------------------------------------------------------------------------*/
-/*--------------------- PRIVATE PART OF THE CLASS --------------------------*/
-/*--------------------------------------------------------------------------*/
-/*--                                                                      --*/
-/*-- Nobody should ever look at this part: everything that is under this  --*/
-/*-- advice may be changed without notice in any new release of the code. --*/
-/*--                                                                      --*/
-/*--------------------------------------------------------------------------*/
-
- private:
-
-/*--------------------------------------------------------------------------*/
-/*-------------------------- PRIVATE DATA TYPES ----------------------------*/
-/*--------------------------------------------------------------------------*/
-/*--                                                                      --*/
-/*-- Let T \subseteq A be a spanning tree, and consider some node v \in V --*/
-/*-- \setminus { 0 }. There is an unique (undirected) path, denoted by    --*/
-/*-- P(v), defined by T from v to the root node 0. The arc in P(v), which --*/
-/*-- is incident to v, is called the *basic arc* of v. The other terminal --*/
-/*-- node u of this basic arc is called the *father node* of v. The       --*/
-/*-- spanning tree T is represented saving the basic arc of every node,   --*/
-/*-- and maintaining the order of the nodes and the depth as to T root    --*/
-/*-- after a Post-Visit of T. This order is saved in a bidirectional list --*/
-/*-- written in the node.                                                 --*/
-/*--                                                                      --*/
-/*-- The Primal Simplex uses a different data structure than the Dual     --*/
-/*-- Simplex, because the Dual Simplex needs additional data (mainly the  --*/
-/*-- Backward Star and Forward Star. Furthermore, the Primal Simplex uses --*/
-/*-- different data structures in the quadratic case.                     --*/
-/*--                                                                      --*/
-/*--------------------------------------------------------------------------*/
-
- struct arcPType;   // pre-declaration of the arc structure (pointers to arcs
-                    // are contained in the node structure) for Primal Simplex
-
- struct arcDType;   // pre-declaration of the arc structure (pointers to arcs
-                    // are contained in the node structure) for Dual Simplex
-
- typedef double iteratorType;  // type for the iteration counter and array
-                               // "whenInT2"
-
- struct nodePType {       // node structure for Primal Simplex - - - - - - - -
-  nodePType *prevInT;     // previous node in the order of the Post-Visit on T
-
-  nodePType *nextInT;     // next node in the order of the Post-Visit on T
-
-  arcPType *enteringTArc; // entering basic arc of this node     
-
-  FNumber balance;        // supply/demand of this node; a node is called a
-                          // supply node, a demand node, or a transshipment
-                          // node depending upon whether balance is larger
-                          // than, smaller than, or equal to zero
-  #if( QUADRATICCOST )
-   CNumber sumQuadratic;  // the sum of the quadratic coefficients of the tree's arcs
-                          // from root of T to the node
-
-   FONumber potential;    // the node potential corresponding with the flow
-                          // conservation constrait of this node
-  #else
-   CNumber potential;      // the node potential corresponding with the flow
-                           // conservation constrait of this node
-  #endif
-
-  int subTreeLevel;        // the depth of the node in T as to T root
-
-  };                       // end( struct( nodePType ) )
-
- struct nodeDType {       // node structure for Dual Simplex - - - - - - - - -
-  nodeDType *prevInT;     // previous node in the order of the Post-Visit on T
-
-  nodeDType *nextInT;     // next node in the order of the Post-Visit on T
-
-  arcDType *enteringTArc; // entering basic arc of this node     
-
-  FNumber balance;        // supply/demand of this node; a node is called a
-                          // supply node, a demand node, or a transshipment
-                          // node depending upon whether balance is larger
-                          // than, smaller than, or equal to zero
-
- #if( QUADRATICCOST )
-  CNumber sumQuadratic;   // the sum of the quadratic coefficients of the tree's arcs
-                          // from root of T to the node
-        
-  FONumber potential;     // the node potential corresponding with the flow
-                          // conservation constrait of this node
- #else
-  CNumber potential;      // the node potential corresponding with the flow
-                          // conservation constrait of this node
- #endif
-
-  int subTreeLevel;       // the depth of the node in T as to T root
-  iteratorType whenInT2;  // the last iteration where a node is in subtree T2
-
-  Index numArcs;          // the number of the arcs which enter/exit from node
-  arcDType *firstBs;      // the first arc in the node's Backward Star
-  arcDType *firstFs;      // the first arc in the node's Forward Star
-
-  };                      // end( struct( nodeDType ) )
-
- struct arcPType {        // arc structure for Primal Simplex - - - - - - - -
-  nodePType *tail;        // tail node
-  nodePType *head;        // head node
-  FNumber flow;           // arc flow
-  CNumber cost;           // arc linear cost
-
-  #if( QUADRATICCOST )
-   CNumber quadraticCost; // arc quadratic cost
-  #else
-   char ident;            // tells if arc is deleted, closed, in T, L, or U
-  #endif
-
-  FNumber upper;          // arc upper bound
-  };                      // end( struct( arcPType ) )
-
- struct arcDType {        // arc structure for Primal Simplex - - - - - - - -
-  nodeDType *tail;        // tail node
-  nodeDType *head;        // head node
-  FNumber flow;           // arc flow
-  CNumber cost;           // arc linear cost
-
-  #if( QUADRATICCOST )
-   CNumber quadraticCost; // arc quadratic cost
-  #else
-   char ident;            // indicates if arc is deleted, closed, in T, in L, or in U
-  #endif
-
-  FNumber upper;          // arc upper bound
-  arcDType *nextBs;       // the next arc in the Backward Star of the arc's head
-  arcDType *nextFs;       // the next arc in the Forward Star of the arc's tail
-
-  };                      // end( struct( arcDType ) )
-
- struct primalCandidType {  // Primal Candidate List- - - - - - - - - - - - -
-  arcPType *arc;            // pointer to the violating primal bound arc
-
-  #if(QUADRATICCOST )
-   FONumber absRC;          // absolute value of the arc's reduced cost
-  #else
-   CNumber absRC;           // absolute value of the arc's reduced cost
-  #endif
-  };                        // end( struct( primalCandidateType ) )
-
- struct dualCandidType {    // Dual Candidate List- - - - - - - - - - - - - -
-  nodeDType *node;          //  deepest node violating the dual bound arc
-  FNumber absInfeas;        // absolute value of the arc's flow infeasibility
-  };                        // end( struct( dualCandidateType ) )
-
-/*--------------------------------------------------------------------------*/
-/*----------------------- PRIVATE DATA STRUCTURES  -------------------------*/
-/*--------------------------------------------------------------------------*/
-
- bool usePrimalSimplex;         // TRUE if the Primal Network Simplex is used
-
- char pricingRule;              // which pricing rule is used
-
- nodePType *nodesP;             // vector of nodes: points to the n + 1 node
-                                // structs (including the dummy root node)
-                                // where the first node is indexed by zero
-                                // and the last node is the dummy root node
-
- nodePType *dummyRootP;         // the dummy root node
-
- nodePType *stopNodesP;         // first infeasible node address = nodes + n 
-  
- arcPType *arcsP;               // vector of arcs; this variable points to
-                                // the m arc structs.
-
- arcPType *dummyArcsP;          // vector of artificial dummy arcs: points
-                                // to the artificial dummy arc variables and
-                                // contains n arc structs. The artificial
-                                // arcs are used to build artificial feasible
-                                // starting bases. For each node i, there is
-                                // exactly one dummy arc defined to connect
-                                // the node i with the dummy root node.
-    
- arcPType *stopArcsP;           // first infeasible arc address = arcs + m 
-
- arcPType *stopDummyP;          // first infeasible dummy arc address
-                                // = arcs + m + n
-
- arcPType *arcToStartP;         // Dantzig Rule and First Eligible Arc Rule
-                                // start their search from this arc
-
- nodeDType *nodesD;             // vector of nodes: points to the n + 1 node
-                                // structs (including the dummy root node)
-                                // where the first node is indexed by zero
-                                // and the last node is the dummy root node
-
- nodeDType *dummyRootD;         // the dummy root node
-
- nodeDType *stopNodesD;         // first infeasible node address = nodes + n 
-  
- arcDType *arcsD;               // vector of arcs; this variable points to
-                                // the m arc structs.
-
- arcDType *dummyArcsD;          // vector of artificial dummy arcs: points to
-                                // to the artificial dummy arc variables and
-                                // contains n arc structs. The artificial
-                                // arcs are used to build artificial feasible
-                                // starting bases. For each node i, there is
-                                // exactly one dummy arc defined to connect
-                                // the node i with the dummy root node.
-
- arcDType *stopArcsD;           // first infeasible arc address = arcs + m 
-
- arcDType *stopDummyD;          // first infeasible dummy arc address
-                                // = arcs + m + n
-
- arcDType *arcToStartD;         // Dantzig Rule and First Eligible Arc Rule
-                                // start their search from this arc
-
- iteratorType iterator;         // the current number of iterations
-    
- primalCandidType *candP;       // every element points to an element of the
-                                // arcs vector which contains an arc violating 
-                                // dual bound
-
- dualCandidType *candD;         // every element points to an element of the
-                                // arcs vector which contains an arc violating 
-                                // primal bond
-
- Index numGroup;                // number of the candidate lists
-    
- Index tempCandidateListSize;   // hot list dimension (it is variable)
-    
- Index groupPos;                // contains the actual candidate list
-    
- Index numCandidateList;        // number of candidate lists
-    
- Index hotListSize;             // number of candidate lists and hot list dimension
-
- Index forcedNumCandidateList;  // value used to force the number of candidate list
-    
- Index forcedHotListSize;       // value used to force the number of candidate list
-                                // and hot list dimension
-
- bool newSession;               // true if algorithm is just started
-  
- CNumber MAX_ART_COST;          // large cost for artificial arcs
-
- FNumber *modifiedBalance;      // vector of balance used by the PostVisit
-
- FONumber EpsOpt;               // the precision of the objective function value
-                                // for the quadratic case of the Primal Simplex
-
- int recomputeFOLimits;         // after how many iterations the quadratic Primal
-                                // Simplex computes "manually" the o.f. value
-
- #if( QUADRATICCOST )
-  FONumber foValue;             // the temporary objective function value
- #endif
-
-/*--------------------------------------------------------------------------*/
-/*-------------------------- PRIVATE METHODS -------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-  void MemAlloc( void );
-
-/**< Method to allocate memory for the main data structures. It creates the
-   vector of nodes, the vector of arcs (real and dummy). If the algorithm is
-   using the Dual Simplex, it creates also the vectors whenInT2, firstIn and
-   nextIn, usefull to identify the next entering. */
-
-/*--------------------------------------------------------------------------*/
-
-  void MemDeAlloc( bool whatDeAlloc );
-
-/**< Method to deallocate memory for the main data structures created in
-   MemAlloc(). */
-
-/*--------------------------------------------------------------------------*/
-
-  void MemAllocCandidateList( void );
-
-/**< Method to allocate memory for the data structure used by the Candidate
-   List Pivot Rule. It creates the vector of candP (or candD in the Dual
-   Simplex case), determining the size of this vector on the basis of number
-   of arcs or nodes. */
-
-/*--------------------------------------------------------------------------*/
-
-  void MemDeAllocCandidateList( void );
-
-/**< Method to deallocate memory for the data structures created in
-   MemAllocCandidateList(). */
-
-/*--------------------------------------------------------------------------*/
-
-  void CreateInitialPrimalBase( void );
-
-/**< Method to create an initial feasible primal base. Add one node (dummyRoot)
-   to the network and connect this dummy root to each real node with n dummy
-   arcs. For each source node i, a dummy arc (i, r) exists. For each sink or
-   transit node j, a dummy arc (r, j) exists. The source nodes send their flows
-   to the dummy root through their dummy arcs, so the dummy root send all to
-   the sink nodes. The dummy root balance is 0, the costs of dummy arcs are
-   fixed to "infinity". */
-
-/*--------------------------------------------------------------------------*/
-
-  void CreateInitialDualBase( void );
-
-/**< Method to create an initial feasible dual base. Add one node (dummyRoot)
-   to the network and connect this dummy root to each real node with n dummy
-   arcs. The source nodes send their flows to the dummy root through their
-   dummy arcs, so the dummy root send all to the sink nodes. The dummy root
-   balance is 0, the costs of dummy arcs are fixed to "infinity". */
-
-/*--------------------------------------------------------------------------*/
-
-  void CreateAdditionalDualStructures( void );
-
-/**< The Dual Simplex needs nodes' Backward and Forward Star to work. So when
-   the Primal Simplex runs, these structures don't exist. When the Dual
-   Simplex starts, these structure are created in this method. */
-
-/*--------------------------------------------------------------------------*/
-
-  void PrimalSimplex( void );
-
-/**< Main method to implement the Primal Simplex algorithm. */
-
-/*--------------------------------------------------------------------------*/
-
-  void DualSimplex( void );
-
-/**< Main method to implement the Dual Simplex algorithm. */
-
-/*--------------------------------------------------------------------------*/
-
-  template<class N, class A>
-  void UpdateT( A *h , A *k , N *h1 , N *h2 , N *k1 , N *k2 );
-
-  /**< Method to update the spanning tree T every iteration.
-     The spanning tree T is implemented with a bidirectional list stored
-     in the node structure, which represents the nodes' order after a
-     Post-Visit. The parameter of the method are the outgoing arc "h", the
-     incoming arc "k", and the four nodes on the extremity of these arcs
-     (for example h2 is the deepest node of the outgoing arc). Removing the
-     arc "h" splits T in two subtrees: T1 (which contains the root of T) and
-     T2, which will be re-connected by the incoming arc "k".
-     T2 will be reordered; in fact, the node "k2" becomes the root of T2 
-     instead of "h2" and the hierarchy of T2 will be overturned. Then T2 will
-     be moved; the root of T2 is changed, therefore the predecessor of the
-     root will become the node "k1".
-
-     This method uses the methods cutAndUpdateSubtree() and pasteSubtree().
-     First it cuts the node "k2" and its subtree from T using
-     cutAndUpdateSubtree(). Moreover the method cutAndUpdateSubtree()
-     updates the field "subTreeLevel" of every subtree's nodes, since k2's
-     subtree will be moved from the bottom to the top of T2. Then the method
-     pasteSubtree() puts this subtree in the bidirectional list after the node
-     "k1". The same operations will be applied to the old precedessor of "k2" 
-     (which will become one of the childs of "k2"). This second subtree will
-     be cut, the subTreeLevel fields will be updated, and it will be inserted
-     in the bidirectional list after the k2's subtree. This is iterated until
-     the node "h2" is reached. */
-
-/*--------------------------------------------------------------------------*/
-
-  template<class N>
-  N* CutAndUpdateSubtree( N *root, int delta );
-
-/**< This method cuts a generic subtree from the spanning tree T. Then it
-   updates the field "subTreeLevel" of every subtree's nodes adding the value
-   "delta". This method returns the last node of the subtree. */
-
-/*--------------------------------------------------------------------------*/
-
-  template<class N>
-  void PasteSubtree( N *root , N *lastNode , N *previousNode );
-
-/**< This method inserts a generic subtree with root passed by parameter into
-   the spanning tree T, between the nodes "previousNode" and "lastNode". */
-
-/*--------------------------------------------------------------------------*/
-
-  arcPType* RuleDantzig( void );
-
-/**< This method returns an arc which violates the dual conditions. It searchs
-   the arc with most violation of dual conditions in the entire set of real
-   arcs. It can be used only by the Primal Simplex in the case of networks
-   with linear costs. */
-
-/*--------------------------------------------------------------------------*/
-
-  arcPType* PRuleFirstEligibleArc( void );
-
-/**< This method returns the first found arc which violates the dual conditions
-   in the case of Primal Simplex, the primal condition in the case of Dual
-   Simplex. It can be used only in the case of networks with linear costs. */
-
-/*--------------------------------------------------------------------------*/
-
-  arcDType* DRuleFirstEligibleArc( void );
-
-/**< This method returns the first found arc which violates the dual conditions
-   in the case of Primal Simplex, the primal condition in the case of Dual
-   Simplex. It can be used only in the case of networks with linear costs. */
-
-/*--------------------------------------------------------------------------*/
-
-  arcPType* RulePrimalCandidateListPivot( void );
-
-/**< This method returns an arc which violates the dual conditions. It searches
-   the arc with most violation of dual conditions in a small set of candidate
-   arcs. In every iteration the method rebuilds this set of arcs executing three
-   phases:
-
-   - in the first phase it analyzes the remaining arcs and delete the arcs which
-     don't violate the dual condition any more;
-
-   - in the second phase it tries to fill the set, so it searchs other arcs 
-     which violate the dual condition: the set of arcs is divided into "buckets"
-     which are searched sequentially until the candidate list is full; the
-     last visited bucket is retained, and the search is restarted from that
-     one at later iterations
-
-   - in the third phase the small set of candidate arcs is ordered according
-     to the violation of dual condition by the method SortPrimalCandidateList() 
-     using an implementation of the algorithm "quicksort".
-
-    At last the method returns the first arc in the ordered small set. If the
-    arc doesn't exist (the set is empty), it returns NULL. */
-
-/*--------------------------------------------------------------------------*/
-
-  inline void InitializePrimalCandidateList( void );
-
-/**< Method to initialize some important details for Primal Candidate List
-   Rule. */
-
-/*--------------------------------------------------------------------------*/
-
-  inline void SortPrimalCandidateList( Index min , Index max );
-
-/**< Method to order the little set of candidate arcs according to
-   infeasibility of dual conditions. It implements the "quicksort"
-   algorithm.  */
-
-/*--------------------------------------------------------------------------*/
-
-  arcDType* RuleDualCandidateListPivot( void );
-
-/**< Similar to RulePrimalCandidateListPivot() for the Dual Simplex. */
-
-/*--------------------------------------------------------------------------*/
-
-  inline void InitializeDualCandidateList( void );
-
-/**< Method to initialize some important details for Dual Candidate List Rule.
-    */
-
-/*--------------------------------------------------------------------------*/
-
-  inline void SortDualCandidateList( Index min , Index max );
-
-/**< Similar to SortPrimalCandidateList() for the Dual Simplex. */
-
-/*--------------------------------------------------------------------------*/
-
-  template<class N, class RCT>
-  inline void AddPotential( N *r , RCT delta );
-
-/**< Method to quickly update the dual solutions. During the change of the
-   base, the potential of node "k2" (deepest node in T of incoming arc "k",
-   and new root of T2) changes according to the new structure of T. In fact,
-   the precedessor of "k2" changes: now the predecessor of "k2"is "k1" (the
-   other node of the incoming arc "k"). This change of predecessor causes a
-   change of potential of "k2". The change of potential of "k2" causes the
-   changes of potential of all nodes of T2. This method computes the change
-   of potential of "k2" and applies it on all the nodes of T2. */
-
-/*--------------------------------------------------------------------------*/
-
-  template<class N>
-  inline void ComputePotential( N *r );
-
-/**< Method to update the dual solutions. It computes all the potential of the
-   nodes of the subtree which has r as root. */
-
-/*--------------------------------------------------------------------------*/
-
-  inline void ResetWhenInT2( void );
-
-/**< Method to order the small set of candidate arcs according to dual
-   infeasibility. It implements the algorithm "quicksort". */
-
-/*--------------------------------------------------------------------------*/
-
-  void CreateInitialPModifiedBalanceVector( void );
-
-/**< Method to initialize the vector of modified balance for the postvisit on
-   T in the Primal Simplex data structure. */
-
-/*--------------------------------------------------------------------------*/
-    
-  void PostPVisit( nodePType *r );
-
-/**< Method to calculate the flow on the basic arcs with the Primal Simplex's
-   data structure. It uses the set of the upper bound arcs, the construction
-   of a modified balance vector and the postvisit on T. */
-
-/*--------------------------------------------------------------------------*/
-
-  void BalanceFlow( nodePType *r );
-
-/**< This method works after the method PostPVisit (in a Primal context). It
-   restores primal admissimibility on the r's subtree. It starts from the leaf
-   of the subtree and goes up to the root, using the method AdjustFlow. */
-
-/*--------------------------------------------------------------------------*/
-
-  void AdjustFlow( nodePType *r );
-
-/**< This method checks the primal admissibility of the r's basic entering arc.
-    If it is out of bounds, the method removes it from the tree (and keeps the
-    relative dummy arc) and push flow in the cycle (some tree's arc and the old
-    entering arc) to restores the right balances of the node. */
-
-/*--------------------------------------------------------------------------*/
-
-  void CreateInitialDModifiedBalanceVector( void );
-
-/**< Method to initialize the vector of modified balance for the postvisit on
-   T in the Dual Simplex data structure. */
-
-/*--------------------------------------------------------------------------*/
-    
-  void PostDVisit( nodeDType *r );
-
-/**< Method to calculate the flow on the basic arcs with the Dual Simplex's data
-   structure, using the set of the upper bound arcs, the construction of a
-   modified balance vector and the postvisit on T. */
-
-/*--------------------------------------------------------------------------*/
-
-  template<class N, class A>
-  inline N* Father( N *n, A *a );
-
-/**< Method to find the predecessor of the node in the tree. */
-
-/*--------------------------------------------------------------------------*/
-
- #if(QUADRATICCOST)
-  template<class A>
-  inline FONumber ReductCost( A *a );
- #else
-  template<class A>
-  inline CNumber ReductCost( A *a );
- #endif
-
-/**< Method to calculate the reduct cost of the arc. */
-
-/*--------------------------------------------------------------------------*/
-
-  inline FONumber GetFO( void );
-
-/**< Method to calculate the temporary (or the final) objective function
-   value. */
-
-/*--------------------------------------------------------------------------*/
-    
-  void PrintPNode( nodePType *nodo );
-
-/**< Method to print the "name" of the node in the Primal Simplex. */
-
-/*--------------------------------------------------------------------------*/
-
-  void PrintPArc( arcPType *arc );
-
-/**< Method to print the "name" of the arc in the Primal Simplex. */
-
-/*--------------------------------------------------------------------------*/
-
-  void PrintDNode( nodeDType *nodo );
-
-/**< Method to print the "name" of the node in the Dual Simplex. */
-
-/*--------------------------------------------------------------------------*/
-
-  void PrintDArc( arcDType *arc );
-
-/**< Method to print the "name" of the arc in the Dual Simplex. */
-
-/*--------------------------------------------------------------------------*/
-
-  nodePType* RecoverPNode( Index ind );
-
-/**< Method to find a node (in the Primal Simplex) using its index. */
-
-/*--------------------------------------------------------------------------*/
-
-  arcPType* RecoverPArc( nodePType *tail , nodePType *head );
-
-/**< Method to find an arc (in the Primal Simplex) using 2 pointers to tail
-   node and head node. */
-
-/*--------------------------------------------------------------------------*/
-
-  nodeDType* RecoverDNode( Index ind );
-
-/**< Method to find a node (in the Dual Simplex) using its index. */
-
-/*--------------------------------------------------------------------------*/
-
-  arcDType* RecoverDArc( nodeDType *tail , nodeDType *head ); 
-
-/**< Method to find an arc (in the Dual Simplex) using 2 pointers to tail
-   node and head node. */
-
-/*--------------------------------------------------------------------------*/
-
-  void infoPNode( nodePType *node , int tab );
-
-/**< Method to print some information of the node (in the Primal Simplex). */
-
-/*--------------------------------------------------------------------------*/
-
-  void infoPArc( arcPType *arc , int ind , int tab ); 
-
-/**< Method to print some information of the arc (in the Primal Simplex). */
-
-/*--------------------------------------------------------------------------*/
-
-  void infoDNode( nodeDType *node , int tab );
-
-/**< Method to print some information of the node (in the Dual Simplex). */
-
-/*--------------------------------------------------------------------------*/
-
-  void infoDArc( arcDType *arc , int ind , int tab );
-
-/**< Method to print some information of the arc (in the Dual Simplex). */
-
-/*--------------------------------------------------------------------------*/
-
-  void ShowSituation( int tab );
-
-/**< Method to show the actual complete situation. */
-
-/*--------------------------------------------------------------------------*/
-    
-  };  // end( class MCFSimplex )
-
-/* @} end( group( MCFSimplex_CLASSES ) ) */
-
-#endif  /* MCFSimplex.h included */
-
-/*-------------------------------------------------------------------------*/
-/*-------------------inline methods implementation-------------------------*/
-/*-------------------------------------------------------------------------*/
-
-inline MCFClass::Index MCFSimplex::MCFSNde( MCFClass::cIndex i )
-{
- if( usePrimalSimplex )
-  return( Index( ( (arcsP + i)->tail - nodesP + 1 ) - USENAME0 ) );
- else
-  return( Index( ( (arcsD + i)->tail - nodesD + 1 ) - USENAME0 ) );
- }
-
-/*-------------------------------------------------------------------------*/
-
-inline MCFClass::Index MCFSimplex::MCFENde( MCFClass::cIndex i )
-{
- if( usePrimalSimplex )
-  return( Index( ( (arcsP + i)->head - nodesP + 1 ) - USENAME0 ) );
- else
-  return( Index( ( (arcsD + i)->head - nodesD + 1 ) - USENAME0 ) );
- }
-
-/*-------------------------------------------------------------------------*/
-
-inline MCFClass::CNumber MCFSimplex::MCFCost( MCFClass::cIndex i )
-{
- if( usePrimalSimplex )
-  return( (arcsP + i)->cost );
- else
-  return( (arcsD + i)->cost );
- }
-
-/*-------------------------------------------------------------------------*/
-
-inline MCFClass::CNumber MCFSimplex::MCFQCoef( MCFClass::cIndex i )
-{
- #if( QUADRATICCOST )
-  if( usePrimalSimplex )
-   return( (arcsP + i)->quadraticCost );
-  else
-   return( (arcsD + i)->quadraticCost );
- #else
-  return( 0 );
- #endif
- }
-
-/*-------------------------------------------------------------------------*/
-
-inline MCFClass::FNumber MCFSimplex::MCFUCap( MCFClass::cIndex i )
-{
- if( usePrimalSimplex )
-  return( (arcsP + i)->upper );
- else
-  return( (arcsD + i)->upper );
- }
-
-/*-------------------------------------------------------------------------*/
-
-inline MCFClass::FNumber MCFSimplex::MCFDfct( MCFClass::cIndex i )
-{
- if( usePrimalSimplex )
-  return( (nodesP + i)->balance );
- else
-  return( (nodesD + i)->balance );
- }
-
-/*-------------------------------------------------------------------------*/
-
-#if( QUADRATICCOST )
-
-template <class A>
-inline MCFSimplex::FONumber MCFSimplex::ReductCost( A *a )
-{
- FONumber redc = (a->tail)->potential - (a->head)->potential;
- redc = redc + a->cost;
- redc = redc + a->quadraticCost * a->flow;
- return( redc );
- }
-
-#else
-
-template <class A>
-inline MCFSimplex::CNumber MCFSimplex::ReductCost( A *a )
-{
- CNumber redc = (a->tail)->potential - (a->head)->potential;
- redc = redc + a->cost;
- return( redc );
- }
-
-#endif
-
-/*-------------------------------------------------------------------------*/
- 
-#if( OPT_USE_NAMESPACES )
-};  // end( namespace MCFClass_di_unipi_it )
-#endif
-
-/*-------------------------------------------------------------------------*/
-/*-------------------------------------------------------------------------*/
-
-/*-------------------------------------------------------------------------*/
-/*---------------------- End File MCFSimplex.h ----------------------------*/
-/*-------------------------------------------------------------------------*/
diff --git a/clang/lib/Analysis/OPTUtils.h b/clang/lib/Analysis/OPTUtils.h
deleted file mode 100755
index 66ad4af..0000000
--- a/clang/lib/Analysis/OPTUtils.h
+++ /dev/null
@@ -1,475 +0,0 @@
-/*--------------------------------------------------------------------------*/
-/*--------------------------- File OPTutils.h ------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @file
- * Small classes are provided for:
- * - reading the time of a code;
- * - generating random numbers.
- *
- * The classes can be adapted to different environments setting a
- * compile-time switch in this file.
- *
- * Additionally, a function is provided for safely reading parameters
- * out of a stream.
- *
- * \version 1.00
- *
- * \date 04 - 10 - 2010
- *
- * \author Antonio Frangioni \n
- *         Operations Research Group \n
- *         Dipartimento di Informatica \n
- *         Universita' di Pisa \n
- *
- * Copyright &copy 1994 - 2010 by Antonio Frangioni
- */
-/*--------------------------------------------------------------------------*/
-/*--------------------------------------------------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#ifndef __OPTutils
- #define __OPTutils   /* self-identification: #endif at the end of the file */
-
-/*--------------------------------------------------------------------------*/
-/*----------------------------- MACROS -------------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @defgroup OPTUTILS_MACROS Compile-time switches in OPTutils.h
-    These macros control how the classes OPTTimers and OPTrand are
-    implemented; choose the appropriate value for your environment,
-    or program a new version if no value suits you.
-    Also, namespaces can be eliminated if they create problems.
-    @{ */
-
-
-/*----------------------- OPT_USE_NAMESPACES -------------------------------*/
-
-#define OPT_USE_NAMESPACES 0
-
-/**< Setting OPT_USE_NAMESPACES == 0 should instruct all codes that use
-   OPTutils stuff to avoid using namespaces; to start with, the common
-   namespace OPTutils_di_unipi_it, that contains all the types defined
-   herein, is *not* defined. */
-
-/*---------------------------- OPT_TIMERS ----------------------------------*/
-
-#define OPT_TIMERS 5
-
-/**< The class OPTtimers is defined below to give an abstract interface to the
-   different timing routines that are used in different platforms. This is
-   needed since time-related functions are one of the less standard parts of
-   the C[++] library. The value of the OPT_TIMERS constant selects among the
-   different timing routines:
-
-   - 1 = Use the Unix times() routine in sys/times.h
-
-   - 2 = As 1 but uses sys/timeb.h (typical for Microsoft(TM) compilers)
-
-   - 3 = Still use times() of sys/times.h, but return wallclock time
-         rather than CPU time
-
-   - 4 = As 3 but uses sys/timeb.h (typical for Microsoft(TM) compilers)
-
-   - 5 = return the user time obtained with ANSI C clock() function; this
-         may result in more accurate running times w.r.t. but may be limited
-         to ~ 72 hours on systems where ints are 32bits.
-
-   - 6 = Use the Unix gettimeofday() routine of sys/time.h.
-
-   Any unsupported value would simply make the class to report constant
-   zero as the time.
-
-   The values 1 .. 4 rely on the constant CLK_TCK for converting between
-   clock ticks and seconds; for the case where the constant is not defined by
-   the compiler -- should not happen, but it does -- or it is defined in a
-   wrong way, the constant is re-defined below. */
-
-/*---------------------------- OPT_RANDOM ---------------------------------*/
-
-#define OPT_RANDOM 1
-
-/**< The class OPTrand is defined below to give an abstract interface to the
-   different random generators that are used in different platforms. This is
-   needed since random generators are one of the less standard parts of the
-   C[++] library. The value of the OPT_RANDOM constant selects among the
-   different timing routines:
-
-   - 0 = an hand-made implementation of a rather good random number generator
-         is used; note that this assumes that long ints >= 32 bits
-
-   - 1 = standard rand() / srand() pair, common to all C libreries but not
-         very sophisticated
-
-   - 2 = drand48() / srand48(), common on Unix architectures and pretty good.
-
-   Any unsupported value would simply make the functions to report constant
-   zero, which is not nice but useful to quickly fix problems if you don't
-   use random numbers at all. */
-
-/*@} -----------------------------------------------------------------------*/
-/*------------------------------ INCLUDES ----------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#if( OPT_RANDOM )
- #include <stdlib.h>
-
- /* For random routines, see OPTrand() below. */
-#endif
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-
-#if( OPT_TIMERS <= 4 )
- #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 3 ) )
-  #include <sys/times.h>
- #else
-  #include <sys/timeb.h>
- #endif
- #include <limits.h>
-#elif( OPT_TIMERS == 5 )
- #include <time.h>
-#elif( OPT_TIMERS == 6 )
- #include <sys/time.h>
-#endif
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-
-#include <iostream>
-#include <sstream>
-
-/* For istream and the >> operator, used in DfltdSfInpt(). */
-
-/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-
-/*--------------------------------------------------------------------------*/
-/*--------------------------- NAMESPACE ------------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#if( OPT_USE_NAMESPACES )
-namespace OPTtypes_di_unipi_it
-{
- /** @namespace OPTtypes_di_unipi_it
-     The namespace OPTtypes_di_unipi_it is defined to hold all the data
-     types, constants, classes and functions defined here. It also
-     comprises the namespace std. */
-#endif
-
- using namespace std;  // I know it's not elegant, but ...
-
-/*--------------------------------------------------------------------------*/
-/*----------------------------- NULL ---------------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#ifndef NULL
- #define NULL 0
-#endif
-
-/* Curiously enough, some compilers do not always define NULL; for instance,
-   sometimes it is defined in stdlib.h. */
-
-/*--------------------------------------------------------------------------*/
-/*---------------------------- CLK_TCK -------------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#ifndef CLK_TCK
- #define CLK_TCK 100
-#endif
-
-/* CLK_TCK is the constant used to transform (integer) time values in seconds
-   for times()-based routines. Its normal value is 100. */
-
-/*--------------------------------------------------------------------------*/
-/*--------------------------- OPT_TIMERS -----------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @defgroup OPTtypes_CLASSES Classes in OPTutils.h
-    @{ */
-
-#if( OPT_TIMERS )
-
-/** Provides a common interface to the different timing routines that are
-    available in different platforms. */
-
-class OPTtimers {
-
- public:  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
-  /// constructor of the class
-  OPTtimers( void )
-  {
-   ReSet();
-   }
-
-  /// start the timer
-  void Start( void )
-  {
-   if( ! ticking ) {
-    #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 2 ) )
-     times( &buff );
-     t_u = buff.tms_utime;
-     t_s = buff.tms_stime;
-    #elif( ( OPT_TIMERS == 3 ) || ( OPT_TIMERS == 4 ) )
-     t_u = times( &buff );
-    #elif( OPT_TIMERS == 5 )
-     t_u = clock();
-    #elif( OPT_TIMERS == 6 )
-     struct timeval t;
-     gettimeofday( &t , NULL );
-     t_u = double( t.tv_sec + t.tv_usec * 1e-6 );
-    #endif
-
-    ticking = 1;
-    }
-   }
-
-  /// stop the timer
-  void Stop( void )
-  {
-   if( ticking ) {
-    Read( u , s );
-    ticking = 0;
-    }
-   }
-
-  /** Return the elapsed time. If the clock is ticking, return the *total*
-    time since the last Start() without stopping the clock; otherwise,
-    return the total elapsed time of all the past runs of the clock since
-    the last ReSet() [see below]. */
-
-  double Read( void )
-  {
-   double tu = 0;
-   double ts = 0;
-   Read( tu , ts );
-   return( tu + ts );
-   }
-
-  /// As Read( void ) but *adds* user and system time to tu and ts.
-  void Read( double &tu , double &ts )
-  {
-   if( ticking ) {
-    #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 2 ) )
-     times( &buff );
-     tu += ( double( buff.tms_utime - t_u ) ) / double( CLK_TCK );
-     ts += ( double( buff.tms_stime - t_s ) ) / double( CLK_TCK );
-    #elif( ( OPT_TIMERS == 3 ) || ( OPT_TIMERS == 4 ) )
-     tu += ( double( times( &buff ) - t_u ) ) / double( CLK_TCK );
-    #elif( OPT_TIMERS == 5 )
-     tu += ( double( clock() - t_u ) ) / double( CLOCKS_PER_SEC );
-    #elif( OPT_TIMERS == 6 )
-     struct timeval t;
-     gettimeofday( &t , NULL );
-     tu += double( t.tv_sec + t.tv_usec * 1e-6 ) - t_u;
-    #endif
-    }
-   else { tu += u; ts += s; }
-   }
-
-  /// reset the timer
-  void ReSet( void )
-  {
-   u = s = 0; ticking = 0;
-   }
-
- private:  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
- double u;      // elapsed *user* time, in seconds
- double s;      // elapsed *system* time, in seconds
- char ticking;  // 1 if the clock is ticking
-
- #if( ( OPT_TIMERS > 0 ) && ( OPT_TIMERS <= 5 ) )
-  clock_t t_u;
-
-  #if( OPT_TIMERS <= 4 )
-   struct tms buff;
-
-   #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 3 ) )
-    clock_t t_s;
-   #endif
-  #endif
- #elif( OPT_TIMERS == 6 )
-  double t_u;
- #endif
-
- };  // end( class OPTtimers );
-
-#endif
-
-/*--------------------------------------------------------------------------*/
-/*------------------------------ OPTrand() ---------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-/** Provide a common interface to the different random generators that are
-    available in different platforms. */
-
-class OPTrand {
-
- public:  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
-  /// constructor of the class
-  OPTrand( void )
-  {
-   #if( OPT_RANDOM == 0 )
-    A[ 0 ] = -1;
-    srand( long( 123456789 ) );
-   #else
-    OPTrand::srand( long( 1 ) );
-   #endif
-   }
-
-  /** Returns a random number uniformly distributed in [0, 1).
-      \note each object of class OPTrand has its own sequence, so that
-      multiple OPTrand objects being used within the same program do not
-      interfere with each other (as opposed to what C random routines
-      would do). */
-
-  double rand( void )
-  {
-   #if( OPT_RANDOM == 0 )
-    long nmbr = *(gb_fptr--);
-    if( nmbr < 0 )
-     nmbr = gb_flip_cycle();
-    return( double( nmbr ) / double( (unsigned long)0x80000000 ) );
-   #elif( OPT_RANDOM == 1 )
-    ::srand( myseed );
-    myseed = ::rand();
-    return( double( myseed ) / double( RAND_MAX ) );
-   #elif( OPT_RANDOM == 2 )
-    return( erand48( myseed ) );
-   #else
-    return( 0 );  // random, eh?
-   #endif
-   }
-
-  /// Seeds the random generator for this instance of OPTrand.
-  void srand( long seed )
-  {
-   #if( OPT_RANDOM == 0 )
-    long prev = seed , next = 1;
-    seed = prev = mod_diff( prev , 0 );
-    A[ 55 ] = prev;
-
-    for( long i = 21 ; i ; i = ( i + 21 ) % 55 ) {
-     A[ i ] = next;
-     next = mod_diff( prev , next );
-     if( seed & 1 )
-      seed = 0x40000000 + ( seed >> 1 );
-     else
-      seed >>= 1;
-
-     next = mod_diff( next , seed );
-     prev = A[ i ];
-     }
-
-    gb_flip_cycle();
-    gb_flip_cycle();
-    gb_flip_cycle();
-    gb_flip_cycle();
-    gb_flip_cycle();
-   #elif( OPT_RANDOM == 1 )
-    myseed = int( seed );
-   #elif( OPT_RANDOM == 2 )
-    long *sp = (long*)( &myseed );
-    *sp = seed;                // copy higher 32 bits
-    myseed[ 2 ] = 0x330E;      // initialize lower 16 bits
-   #endif
-   }
-
- private:  //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-
- #if( OPT_RANDOM == 0 )
-  long A[ 56 ];
-  long *gb_fptr;
-
-  long mod_diff( long x , long y )
-  {
-   return( ( ( x ) - ( y ) ) & 0x7fffffff );
-   }
-
-  long gb_flip_cycle( void )
-  {
-   long *ii, *jj;
-   for( ii = &A[ 1 ] , jj = &A[ 32 ] ; jj <= &A[ 55 ] ; ii++ , jj++ )
-    *ii = mod_diff( *ii , *jj );
-
-   for( jj = &A[ 1 ] ; ii <= &A[ 55 ] ; ii++ , jj++ )
-    *ii = mod_diff( *ii , *jj );
-
-   gb_fptr = &A[ 54 ];
-
-   return A[ 55 ];
-   }
- #elif( OPT_RANDOM == 1 )
-  int myseed;
- #elif( OPT_RANDOM == 2 )
-  unsigned short int myseed[ 3 ];
- #endif
-
- };  // end( class( OPTrand ) )
-
-/* @} end( group( OPTtypes_CLASSES ) ) */
-/*--------------------------------------------------------------------------*/
-/*--------------------------- DfltdSfInpt() --------------------------------*/
-/*--------------------------------------------------------------------------*/
-/** @defgroup OPTtypes_FUNCTIONS Functions in OPTutils.h
-    @{ */
-
-/** Template function for reading parameters from a istream. The function is
-   "safe" because it works also if the istream is not given, is not be long
-   enough or contains erroneous things.
-
-   Given a &istream (possibly NULL), DfltdSfInpt() attempts to read Param out
-   of it, skipping any line that begins with the comment carachter (defaulted
-   to '#'), any blank line and any line starting with anything that can not
-   be interpreted as a `T'. If, for any reason, the read operation fails,
-   then the parameter is given the default value `Dflt'. Otherwise, all the
-   rest of the line up to the nearest newline ('\n') carachter is flushed. */
-
-template<class T>
-inline void DfltdSfInpt( istream *iStrm , T &Param , const T Dflt ,
-                         const char cmntc = '#' )
-{
- string comm( 1 , cmntc );
-
- // the "! ! stream" trick is there to force the compiler to apply the
- // stream -> bool conversion, in case it is too dumb to do it by itself
- if( iStrm && ( ! ( ! (*iStrm) ) ) ) {
-  string buf;
-
-  for( ;; ) {
-   if( ! ( (*iStrm) >> ws ) )   // skip whitespace
-    break;
-
-   if( ! (*iStrm).good() )      // check stream is OK
-    break;
-
-   getline( *iStrm , buf );
-
-   if( ! buf.empty() )
-    if( buf.substr( 0 , 1 ).compare( comm ) ) {
-     stringstream ss;
-     ss << buf;
-     if( ! ( ss >> Param ) )
-      break;
- 
-     return;
-     }
-   }
-  }
-
- Param = Dflt; 
-
- }  // end( DfltdSfInpt )
-
-/* @} end( group( OPTtypes_FUNCTIONS ) ) */
-/*--------------------------------------------------------------------------*/
-
-#if( OPT_USE_NAMESPACES )
- };  // end( namespace OPTtypes_di_unipi_it )
-#endif
-
-/*--------------------------------------------------------------------------*/
-/*--------------------------------------------------------------------------*/
-
-#endif  /* OPTutils.h included */
-
-/*--------------------------------------------------------------------------*/
-/*---------------------- End File OPTutils.h -------------------------------*/
-/*--------------------------------------------------------------------------*/
diff --git a/impl/Operator.hpp b/impl/Operator.hpp
index 4a573b8..0b08866 100644
--- a/impl/Operator.hpp
+++ b/impl/Operator.hpp
@@ -153,8 +153,8 @@ struct MinCostFlow : public Operator<Domain> {
       ends[i] = arcs[i].second;
     }
 
-    _solver.LoadNet(supplies.size(), supplies.size(), // max nodes/arcs
-                    supplies.size(), supplies.size(), // current nodes/arcs
+    _solver.LoadNet(supplies.size(), arcs.size(), // max nodes/arcs
+                    supplies.size(), arcs.size(), // current nodes/arcs
                     NULL, NULL, // arcs have inf cap, zero cost (to begin)
                     deficits, // deficits for each node
                     starts, ends); // start/end of each edge