diff options
| -rw-r--r-- | clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp | 8 | ||||
| -rw-r--r-- | clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp | 12 | ||||
| -rw-r--r-- | clang/include/clang/Analysis/Analyses/IntervalSolver/Log.hpp | 6 | ||||
| -rw-r--r-- | clang/include/clang/Analysis/Analyses/IntervalSolver/MaxStrategy.hpp | 82 | ||||
| -rw-r--r-- | clang/include/clang/Analysis/Analyses/IntervalSolver/Operator.hpp | 82 | ||||
| -rw-r--r-- | clang/include/clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp | 34 | ||||
| -rw-r--r-- | clang/lib/Analysis/Interval.cpp | 496 | ||||
| -rwxr-xr-x | clang/lib/Analysis/MCFClass.h | 2436 | ||||
| -rw-r--r-- | clang/lib/Analysis/MCFSimplex.cpp | 68 | ||||
| -rw-r--r-- | clang/lib/Analysis/MCFSimplex.h | 1086 | ||||
| -rwxr-xr-x | clang/lib/Analysis/OPTUtils.h | 475 | ||||
| -rw-r--r-- | impl/Operator.hpp | 4 |
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 = ρ } - 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 © 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 © 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 © 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 |