summaryrefslogtreecommitdiff
path: root/clang
diff options
context:
space:
mode:
Diffstat (limited to 'clang')
-rw-r--r--clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp10
-rw-r--r--clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp12
-rw-r--r--clang/lib/Analysis/CMakeLists.txt6
-rw-r--r--clang/lib/Analysis/Interval.cpp163
-rwxr-xr-xclang/lib/Analysis/MCFClass.h2436
-rw-r--r--clang/lib/Analysis/MCFSimplex.cpp4197
-rw-r--r--clang/lib/Analysis/MCFSimplex.h1086
-rwxr-xr-xclang/lib/Analysis/OPTUtils.h475
8 files changed, 8352 insertions, 33 deletions
diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp
index e3ec15a..664d71f 100644
--- a/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp
+++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp
@@ -99,6 +99,16 @@ struct Complete {
template<typename Z>
friend std::ostream& operator<<(std::ostream&, const Complete<Z>&);
+ template<typename S>
+ S as() const {
+ if (_infinity) {
+ if (_value > 0)
+ return infinity<S>();
+ return -infinity<S>();
+ }
+ return (S) _value;
+ }
+
private:
T _value;
bool _infinity;
diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp
index 3342cc7..5ee5405 100644
--- a/clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp
+++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/EquationSystem.hpp
@@ -76,9 +76,14 @@ struct EquationSystem {
}
Constant<Domain>& constant(const Domain& value) {
- Constant<Domain>* constant = new Constant<Domain>(value);
- _expressions.insert(constant);
- return *constant;
+ 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];
+ }
}
MaxExpression<Domain>* operator[](const Variable<Domain>& var) const {
@@ -127,6 +132,7 @@ 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/lib/Analysis/CMakeLists.txt b/clang/lib/Analysis/CMakeLists.txt
index 3d2251e..42f4e11 100644
--- a/clang/lib/Analysis/CMakeLists.txt
+++ b/clang/lib/Analysis/CMakeLists.txt
@@ -10,6 +10,7 @@ add_clang_library(clangAnalysis
Dominators.cpp
FormatString.cpp
Interval.cpp
+ MCFSimplex.cpp
LiveVariables.cpp
PostOrderCFGView.cpp
PrintfFormatString.cpp
@@ -21,8 +22,7 @@ add_clang_library(clangAnalysis
UninitializedValues.cpp
)
-add_library(lib_lemon IMPORTED emon)
+add_library(IMPORTED emon)
add_dependencies(clangAnalysis ClangAttrClasses ClangAttrList
- ClangDiagnosticAnalysis ClangDeclNodes ClangStmtNodes
- lib_lemon)
+ ClangDiagnosticAnalysis ClangDeclNodes ClangStmtNodes)
diff --git a/clang/lib/Analysis/Interval.cpp b/clang/lib/Analysis/Interval.cpp
index 7096d75..e643f5e 100644
--- a/clang/lib/Analysis/Interval.cpp
+++ b/clang/lib/Analysis/Interval.cpp
@@ -8,6 +8,8 @@
#include "clang/Analysis/Analyses/IntervalSolver/Complete.hpp"
#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"
@@ -41,6 +43,11 @@ std::ostream& operator<<(std::ostream& cout, const std::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) {
+ cout << "(" << v.first->name() << ", " << v.second << ")";
+ return cout;
+}
template<typename V>
std::ostream& operator<<(std::ostream& cout, const std::vector<V>& v) {
@@ -485,6 +492,86 @@ Condition fromComparison(const BinaryOperator* op, bool negate) {
return cond;
}
+
+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;
@@ -551,28 +638,41 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
if (negative) {
result = negate_result(result);
}
- EqnExpr* expression = &system.constant(result.second);
- 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");
- std::vector<EqnExpr*> additionArgs;
- additionArgs.push_back(expression);
-
- if (it->second == 1) {
- additionArgs.push_back(vars[it->first]);
- } else {
- std::vector<EqnExpr*> multiplicationArgs;
- multiplicationArgs.push_back(vars[it->first]);
- multiplicationArgs.push_back(&system.constant(it->second));
- additionArgs.push_back(&system.expression(new Multiplication<ZBar>(), multiplicationArgs));
- }
- expression = &system.expression(new Addition<ZBar>(), additionArgs);
+ 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(),
+ 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);
+ }
+ EqnExpr* minCostExpr = &system.expression(new MinCostFlow(varCosts), minCostArgs);
+
+ // add the constant factor to the min-cost bit
+ std::vector<EqnExpr*> additionArgs;
+ additionArgs.push_back(&system.constant(result.second));
+ additionArgs.push_back(minCostExpr);
+ expression = &system.expression(new Addition<ZBar>(), additionArgs);
+ } else {
+ expression = &system.constant(result.second);
}
-
+
+ // max(-inf, expr), so strategy iteration will work
std::vector<EqnExpr*> maxArgs;
maxArgs.push_back(&system.constant(-infinity<ZBar>()));
maxArgs.push_back(expression);
@@ -612,6 +712,7 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
EqnExpr* expr = NULL;
if (val == -infinity<ZBar>()) {
// don't do anything here: min(-inf, x) = -inf (for all x)
+ expr = &system.constant(-infinity<ZBar>());
} else if (val == infinity<ZBar>()) {
// no need to have a min here: min(inf, x) = x (for all x)
expr = jt->second;
@@ -623,7 +724,9 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
expr = &system.expression(new Minimum<ZBar>(), minArgs);
}
- if (expr) {
+ system[*var]->arguments().push_back(expr);
+
+ /*if (expr) {
std::set<std::string> ignore;
for (VarMap::iterator
variables = vars.begin(),
@@ -638,10 +741,13 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
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]));
@@ -656,7 +762,7 @@ void runOnBlock(const CFGBlock* block, EqnSys& system, BlockVars& block_vars) {
expr = &system.expression(new Guard<ZBar>(), guard_args);
}
system[*var]->arguments().push_back(expr);
- }
+ }*/
}
}
}
@@ -683,6 +789,7 @@ void IntervalAnalysis::runOnAllBlocks(const Decl& decl) {
BlockVars block_vars;
std::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());
@@ -690,11 +797,13 @@ void IntervalAnalysis::runOnAllBlocks(const Decl& decl) {
for (unsigned int i = func->getNumParams(); i > 0; i--) {
std::string name = func->getParamDecl(i-1)->getNameAsString();
- function_arguments.insert(name); // add it to the first block's vars
- function_arguments.insert(negate_string(name)); // add it to the first block's vars
-
- system[system.variable(name + '-' + block_id + "-pre")] = &system.maxExpression(infArg); // set the vars to infinity in the system here
- system[system.variable(negate_string(name) + '-' + block_id + "-pre")] = &system.maxExpression(infArg); // set the vars to infinity in the system here
+ // add the variables to the first block
+ function_arguments.insert(name);
+ function_arguments.insert(negate_string(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);
}
}
diff --git a/clang/lib/Analysis/MCFClass.h b/clang/lib/Analysis/MCFClass.h
new file mode 100755
index 0000000..41adb82
--- /dev/null
+++ b/clang/lib/Analysis/MCFClass.h
@@ -0,0 +1,2436 @@
+/*--------------------------------------------------------------------------*/
+/*-------------------------- File MCFClass.h -------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @file
+ * Header file for the abstract base class MCFClass, which defines a standard
+ * interface for (linear or convex quadratic separable) Min Cost Flow Problem
+ * solvers, to be implemented as derived classes.
+ *
+ * \version 3.01
+ *
+ * \date 30 - 09 - 2011
+ *
+ * \author Alessandro Bertolini \n
+ * Operations Research Group \n
+ * Dipartimento di Informatica \n
+ * Universita' di Pisa \n
+ *
+ * \author Antonio Frangioni \n
+ * Operations Research Group \n
+ * Dipartimento di Informatica \n
+ * Universita' di Pisa \n
+ *
+ * \author Claudio Gentile \n
+ * Istituto di Analisi di Sistemi e Informatica \n
+ * Consiglio Nazionale delle Ricerche \n
+ *
+ * Copyright &copy 1996 - 2011 by Antonio Frangioni, Claudio Gentile
+ */
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*----------------------------- DEFINITIONS --------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#ifndef __MCFClass
+ #define __MCFClass /* self-identification: #endif at the end of the file */
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------------- MACROS ---------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @defgroup MCFCLASS_MACROS Compile-time switches in MCFClass.h
+ These macros control some important details of the class interface.
+ Although using macros for activating features of the interface is not
+ very C++, switching off some unused features may allow some
+ implementation to be more efficient in running time or memory.
+ @{ */
+
+/*-------------------------------- USENAME0 --------------------------------*/
+
+#define USENAME0 1
+
+/**< Decides if 0 or 1 is the "name" of the first node.
+ If USENAME0 == 1, (warning: it has to be *exactly* 1), then the node
+ names go from 0 to n - 1, otherwise from 1 to n. Note that this does not
+ affect the position of the deficit in the deficit vectors, i.e., the
+ deficit of the i-th node - be its "name" `i' or `i - 1' - is always in
+ the i-th position of the vector. */
+
+/*@} end( group( MCFCLASS_MACROS ) ) */
+/*--------------------------------------------------------------------------*/
+/*------------------------------ INCLUDES ----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#include "OPTUtils.h"
+
+/* OPTtypes.h defines standard interfaces for timing and random routines, as
+ well as the namespace OPTtypes_di_unipi_it and the macro
+ OPT_USE_NAMESPACES, useful for switching off all namespaces in one blow
+ for those strange cases where they create problems. */
+
+#include <iomanip>
+#include <sstream>
+#include <limits>
+#include <cassert>
+
+// QUICK HACK
+#define throw(x) assert(false)
+
+/*--------------------------------------------------------------------------*/
+/*------------------------- NAMESPACE and USING ----------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#if( OPT_USE_NAMESPACES )
+namespace MCFClass_di_unipi_it
+{
+ /** @namespace MCFClass_di_unipi_it
+ The namespace MCFClass_di_unipi_it is defined to hold the MCFClass
+ class and all the relative stuff. It comprises the namespace
+ OPTtypes_di_unipi_it. */
+
+ using namespace OPTtypes_di_unipi_it;
+#endif
+
+/*@} end( group( MCFCLASS_CONSTANTS ) ) */
+/*--------------------------------------------------------------------------*/
+/*-------------------------- CLASS MCFClass --------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @defgroup MCFCLASS_CLASSES Classes in MCFClass.h
+ @{ */
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- GENERAL NOTES --------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** This abstract base class defines a standard interface for (linear or
+ convex quadartic separable) Min Cost Flow (MCF) problem solvers.
+
+ The data of the problem consist of a (directed) graph G = ( N , A ) with
+ n = |N| nodes and m = |A| (directed) arcs. Each node `i' has a deficit
+ b[ i ], i.e., the amount of flow that is produced/consumed by the node:
+ source nodes (which produce flow) have negative deficits and sink nodes
+ (which consume flow) have positive deficits. Each arc `(i, j)' has an
+ upper capacity U[ i , j ], a linear cost coefficient C[ i , j ] and a
+ (non negative) quadratic cost coefficient Q[ i , j ]. Flow variables
+ X[ i , j ] represents the amount of flow to be sent on arc (i, j).
+ Parallel arcs, i.e., multiple copies of the same arc `(i, j)' (with
+ possibily different costs and/or capacities) are in general allowed.
+ The formulation of the problem is therefore:
+ \f[
+ \min \sum_{ (i, j) \in A } C[ i , j ] X[ i, j ] +
+ Q[ i , j ] X[ i, j ]^2 / 2
+ \f]
+ \f[
+ (1) \sum_{ (j, i) \in A } X[ j , i ] -
+ \sum_{ (i, j) \in A } X[ i , j ] = b[ i ]
+ \hspace{1cm} i \in N
+ \f]
+ \f[
+ (2) 0 \leq X[ i , j ] \leq U[ i , j ]
+ \hspace{1cm} (i, j) \in A
+ \f]
+ The n equations (1) are the flow conservation constraints and the 2m
+ inequalities (2) are the flow nonnegativity and capacity constraints.
+ At least one of the flow conservation constraints is redundant, as the
+ demands must be balanced (\f$\sum_{ i \in N } b[ i ] = 0\f$); indeed,
+ exactly n - ConnectedComponents( G ) flow conservation constraints are
+ redundant, as demands must be balanced in each connected component of G.
+ Let us denote by QA and LA the disjoint subsets of A containing,
+ respectively, "quadratic" arcs (with Q[ i , j ] > 0) and "linear" arcs
+ (with Q[ i , j ] = 0); the (MCF) problem is linear if QA is empty, and
+ nonlinear (convex quadratic) if QA is nonempty.
+
+ The dual of the problem is:
+ \f[
+ \max \sum_{ i \in N } Pi[ i ] b[ i ] -
+ \sum_{ (i, j) \in A } W[ i , j ] U[ i , j ] -
+ \sum_{ (i, j) \in AQ } V[ i , j ]^2 / ( 2 * Q[ i , j ] )
+ \f]
+ \f[
+ (3.a) C[ i , j ] - Pi[ j ] + Pi[ i ] + W[ i , j ] - Z[ i , j ] = 0
+ \hspace{1cm} (i, j) \in AL
+ \f]
+ \f[
+ (3.b) C[ i , j ] - Pi[ j ] + Pi[ i ] + W[ i , j ] - Z[ i , j ] =
+ V[ i , j ]
+ \hspace{1cm} (i, j) \in AQ
+ \f]
+ \f[
+ (4.a) W[ i , j ] \geq 0 \hspace{1cm} (i, j) \in A
+ \f]
+ \f[
+ (4.b) Z[ i , j ] \geq 0 \hspace{1cm} (i, j) \in A
+ \f]
+
+ Pi[] is said the vector of node potentials for the problem, W[] are
+ bound variables and Z[] are slack variables. Given Pi[], the quantities
+ \f[
+ RC[ i , j ] = C[ i , j ] + Q[ i , j ] * X[ i , j ] - Pi[ j ] + Pi[ i ]
+ \f]
+ are said the "reduced costs" of arcs.
+
+ A primal and dual feasible solution pair is optimal if and only if the
+ complementary slackness conditions
+ \f[
+ RC[ i , j ] > 0 \Rightarrow X[ i , j ] = 0
+ \f]
+ \f[
+ RC[ i , j ] < 0 \Rightarrow X[ i , j ] = U[ i , j ]
+ \f]
+ are satisfied for all arcs (i, j) of A.
+
+ The MCFClass class provides an interface with methods for managing and
+ solving problems of this kind. Actually, the class can also be used as
+ an interface for more general NonLinear MCF problems, where the cost
+ function either nonseparable ( C( X ) ) or arc-separable
+ ( \f$\sum_{ (i, j) \in A } C_{i,j}( X[ i, j ] )\f$ ). However, solvers
+ for NonLinear MCF problems are typically objective-function-specific,
+ and there is no standard way for inputting a nonlinear function different
+ from a separable convex quadratic one, so only the simplest form is dealt
+ with in the interface, leaving more complex NonLinear parts to the
+ interface of derived classes. */
+
+class MCFClass {
+
+/*--------------------------------------------------------------------------*/
+/*----------------------- PUBLIC PART OF THE CLASS -------------------------*/
+/*--------------------------------------------------------------------------*/
+/*-- --*/
+/*-- The following methods and data are the actual interface of the --*/
+/*-- class: the standard user should use these methods and data only. --*/
+/*-- --*/
+/*--------------------------------------------------------------------------*/
+
+ public:
+
+/*--------------------------------------------------------------------------*/
+/*---------------------------- PUBLIC TYPES --------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @name Public types
+ The MCFClass defines four main public types:
+
+ - Index, the type of arc and node indices;
+
+ - FNumber, the type of flow variables, arc capacities, and node deficits;
+
+ - CNumber, the type of flow costs, node potentials, and arc reduced costs;
+
+ - FONumber, the type of objective function value.
+
+ By re-defining the types in this section, most MCFSolver should be made
+ to work with any reasonable choice of data type (= one that is capable of
+ properly representing the data of the instances to be solved). This may
+ be relevant due to an important property of MCF problems: *if all arc
+ capacities and node deficits are integer, then there exists an integral
+ optimal primal solution*, and *if all arc costs are integer, then there
+ exists an integral optimal dual solution*. Even more importantly, *many
+ solution algorithms will in fact produce an integral primal/dual
+ solution for free*, because *every primal/dual solution they generate
+ during the solution process is naturally integral*. Therefore, one can
+ use integer data types to represent everything connected with flows and/or
+ costs if the corresponding data is integer in all instances one needs to
+ solve. This directly translates in significant memory savings and/or speed
+ improvements.
+
+ *It is the user's responsibility to ensure that these types are set to
+ reasonable values*. So, the experienced user may want to experiment with
+ setting this data properly if memory footprint and/or speed is a primary
+ concern. Note, however, that *not all solution algorithms will happily
+ accept integer data*; one example are Interior-Point approaches, which
+ require both flow and cost variables to be continuous (float). So, the
+ viability of setting integer data (as well as its impact on performances)
+ is strictly related to the specific kind of algorithm used. Since these
+ types are common to all derived classes, they have to be set taking into
+ account the needs of all the solvers that are going to be used, and
+ adapting to the "worst case"; of course, FNumber == CNumber == double is
+ going to always be an acceptable "worst case" setting. MCFClass may in a
+ future be defined as a template class, with these as template parameters,
+ but this is currently deemed overkill and avoided.
+
+ Finally, note that the above integrality property only holds for *linear*
+ MCF problems. If any arc has a nonzero quadratic cost coefficient, optimal
+ flows and potentials may be fractional even if all the data of the problem
+ (comprised quadratic cost coefficients) is integer. Hence, for *quadratic*
+ MCF solvers, a setting like FNumber == CNumber == double is actually
+ *mandatory*, for any reasonable algorithm will typically misbehave
+ otherwise.
+ @{ */
+
+/*--------------------------------------------------------------------------*/
+
+ typedef unsigned int Index; ///< index of a node or arc ( >= 0 )
+ typedef Index *Index_Set; ///< set (array) of indices
+ typedef const Index cIndex; ///< a read-only index
+ typedef cIndex *cIndex_Set; ///< read-only index array
+
+/*--------------------------------------------------------------------------*/
+
+ typedef int SIndex; ///< index of a node or arc
+ typedef SIndex *SIndex_Set; ///< set (array) of indices
+ typedef const SIndex cSIndex; ///< a read-only index
+ typedef cSIndex *cSIndex_Set; ///< read-only index array
+
+/*--------------------------------------------------------------------------*/
+
+ typedef double FNumber; ///< type of arc flow
+ typedef FNumber *FRow; ///< vector of flows
+ typedef const FNumber cFNumber; ///< a read-only flow
+ typedef cFNumber *cFRow; ///< read-only flow array
+
+/*--------------------------------------------------------------------------*/
+
+ typedef double CNumber; ///< type of arc flow cost
+ typedef CNumber *CRow; ///< vector of costs
+ typedef const CNumber cCNumber; ///< a read-only cost
+ typedef cCNumber *cCRow; ///< read-only cost array
+
+/*--------------------------------------------------------------------------*/
+
+ typedef double FONumber;
+ /**< type of the objective function: has to hold sums of products of
+ FNumber(s) by CNumber(s) */
+
+ typedef const FONumber cFONumber; ///< a read-only o.f. value
+
+/*--------------------------------------------------------------------------*/
+/** Very small class to simplify extracting the "+ infinity" value for a
+ basic type (FNumber, CNumber, Index); just use Inf<type>(). */
+
+ template <typename T>
+ class Inf {
+ public:
+ Inf() {}
+ operator T() { return( std::numeric_limits<T>::max() ); }
+ };
+
+/*--------------------------------------------------------------------------*/
+/** Very small class to simplify extracting the "machine epsilon" for a
+ basic type (FNumber, CNumber); just use Eps<type>(). */
+
+ template <typename T>
+ class Eps {
+ public:
+ Eps() {}
+ operator T() { return( std::numeric_limits<T>::epsilon() ); }
+ };
+
+/*--------------------------------------------------------------------------*/
+/** Small class for exceptions. Derives from std::exception implementing the
+ virtual method what() -- and since what is virtual, remember to always
+ catch it by reference (catch exception &e) if you want the thing to work.
+ MCFException class are thought to be of the "fatal" type, i.e., problems
+ for which no solutions exists apart from aborting the program. Other kinds
+ of exceptions can be properly handled by defining derived classes with
+ more information. */
+
+ class MCFException : public exception {
+ public:
+ MCFException( const char *const msg = 0 ) { errmsg = msg; }
+
+ // wow, such a hack
+#undef throw
+ const char* what( void ) const throw () { return( errmsg ); }
+#define throw(x) assert(false)
+ private:
+ const char *errmsg;
+ };
+
+/*--------------------------------------------------------------------------*/
+/** Public enum describing the possible parameters of the MCF solver, to be
+ used with the methods SetPar() and GetPar(). */
+
+ enum MCFParam { kMaxTime = 0 , ///< max time
+ kMaxIter , ///< max number of iteration
+ kEpsFlw , ///< tolerance for flows
+ kEpsDfct , ///< tolerance for deficits
+ kEpsCst , ///< tolerance for costs
+ kReopt , ///< whether or not to reoptimize
+ kLastParam /**< dummy parameter: this is used to
+ allow derived classes to "extend"
+ the set of parameters. */
+ };
+
+/*--------------------------------------------------------------------------*/
+/** Public enum describing the possible status of the MCF solver. */
+
+ enum MCFStatus { kUnSolved = -1 , ///< no solution available
+ kOK = 0 , ///< optimal solution found
+
+ kStopped , ///< optimization stopped
+ kUnfeasible , ///< problem is unfeasible
+ kUnbounded , ///< problem is unbounded
+ kError ///< error in the solver
+ };
+
+/*--------------------------------------------------------------------------*/
+/** Public enum describing the possible reoptimization status of the MCF
+ solver. */
+
+ enum MCFAnswer { kNo = 0 , ///< no
+ kYes ///< yes
+ };
+
+/*--------------------------------------------------------------------------*/
+/** Public enum describing the possible file formats in WriteMCF(). */
+
+ enum MCFFlFrmt { kDimacs = 0 , ///< DIMACS file format for MCF
+ kQDimacs , ///< quadratic DIMACS file format for MCF
+ kMPS , ///< MPS file format for LP
+ kFWMPS ///< "Fixed Width" MPS format
+ };
+
+/*--------------------------------------------------------------------------*/
+/** Base class for representing the internal state of the MCF algorithm. */
+
+ class MCFState {
+ public:
+ MCFState( void ) {};
+ virtual ~MCFState() {};
+ };
+
+ typedef MCFState *MCFStatePtr; ///< pointer to a MCFState
+
+/*@} -----------------------------------------------------------------------*/
+/*--------------------------- PUBLIC METHODS -------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*---------------------------- CONSTRUCTOR ---------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @name Constructors
+ @{ */
+
+ MCFClass( cIndex nmx = 0 , cIndex mmx = 0 )
+ {
+ nmax = nmx;
+ mmax = mmx;
+ n = m = 0;
+
+ status = kUnSolved;
+ Senstv = true;
+
+ EpsFlw = Eps<FNumber>() * 100;
+ EpsCst = Eps<CNumber>() * 100;
+ EpsDfct = EpsFlw * ( nmax ? nmax : 100 );
+
+ MaxTime = 0;
+ MaxIter = 0;
+
+ MCFt = NULL;
+ }
+
+/**< Constructor of the class.
+
+ nmx and mmx, if provided, are taken to be respectively the maximum number
+ of nodes and arcs in the network. If nonzero values are passed, memory
+ allocation can be anticipated in the constructor, which is sometimes
+ desirable. The maximum values are stored in the protected fields nmax and
+ mmax, and can be changed with LoadNet() [see below]; however, changing
+ them typically requires memory allocation/deallocation, which is
+ sometimes undesirable outside the constructor.
+
+ After that an object has been constructed, no problem is loaded; this has
+ to be done with LoadNet() [see below]. Thus, it is an error to invoke any
+ method which requires the presence of a problem (typicall all except those
+ in the initializations part). The base class provides two protected fields
+ n and m for the current number of nodes and arcs, respectively, that are
+ set to 0 in the constructor precisely to indicate that no instance is
+ currently loaded. */
+
+/*@} -----------------------------------------------------------------------*/
+/*-------------------------- OTHER INITIALIZATIONS -------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @name Other initializations
+ @{ */
+
+ virtual void LoadNet( cIndex nmx = 0 , cIndex mmx = 0 , cIndex pn = 0 ,
+ cIndex pm = 0 , cFRow pU = NULL , cCRow pC = NULL ,
+ cFRow pDfct = NULL , cIndex_Set pSn = NULL ,
+ cIndex_Set pEn = NULL ) = 0;
+
+/**< Inputs a new network.
+
+ The parameters nmx and mmx are the new max number of nodes and arcs,
+ possibly overriding those set in the constructor [see above], altough at
+ the likely cost of memory allocation and deallocation. Passing nmx ==
+ mmx == 0 is intended as a signal to the solver to deallocate everything
+ and wait for new orders; in this case, all the other parameters are ignored.
+
+ Otherwise, in principle all the other parameters have to be provided.
+ Actually, some of them may not be needed for special classes of MCF
+ problems (e.g., costs in a MaxFlow problem, or start/end nodes in a
+ problem defined over a graph with fixed topology, such as a complete
+ graph). Also, passing NULL is allowed to set default values.
+
+ The meaning of the parameters is the following:
+
+ - pn is the current number of nodes of the network (<= nmax).
+ - pm is the number of arcs of the network (<= mmax).
+
+ - pU is the m-vector of the arc upper capacities; capacities must be
+ nonnegative, but can in principle be infinite (== F_INF); passing
+ pU == NULL means that all capacities are infinite;
+
+ - pC is the m-vector of the arc costs; costs must be finite (< C_INF);
+ passing pC == NULL means that all costs must be 0.
+
+ - pDfct is the n-vector of the node deficits; source nodes have negative
+ deficits and sink nodes have positive deficits; passing pDfct ==
+ NULL means that all deficits must be 0 (a circulation problem);
+
+ - pSn is the m-vector of the arc starting nodes; pSn == NULL is in
+ principle not allowed, unless the topology of the graph is fixed;
+ - pEn is the m-vector of the arc ending nodes; same comments as for pSn.
+
+ Note that node "names" in the arrays pSn and pEn must go from 1 to pn if
+ the macro USANAME0 [see above] is set to 0, while they must go from 0 to
+ pn - 1 if USANAME0 is set to 1. In both cases, however, the deficit of the
+ first node is read from the first (0-th) position of pDfct, that is if
+ USANAME0 == 0 then the deficit of the node with name `i' is read from
+ pDfct[ i - 1 ].
+
+ The data passed to LoadNet() can be used to specify that the arc `i' must
+ not "exist" in the problem. This is done by passing pC[ i ] == C_INF;
+ solvers which don't read costs are forced to read them in order to check
+ this, unless they provide alternative solver-specific ways to accomplish
+ the same tasks. These arcs are "closed", as for the effect of CloseArc()
+ [see below]. "invalid" costs (== C_INF) are set to 0 in order to being
+ subsequently capable of "opening" them back with OpenArc() [see below].
+ The way in which these non-existent arcs are phisically dealt with is
+ solver-specific; in some solvers, for instance, this could be obtained by
+ simply putting their capacity to zero. Details about these issues should
+ be found in the interface of derived classes.
+
+ Note that the quadratic part of the objective function, if any, is not
+ dealt with in LoadNet(); it can only be separately provided with
+ ChgQCoef() [see below]. By default, the problem is linear, i.e., all
+ coefficients of the second-order terms in the objective function are
+ assumed to be zero. */
+
+/*--------------------------------------------------------------------------*/
+
+ virtual inline void LoadDMX( istream &DMXs , bool IsQuad = false );
+
+/**< Read a MCF instance in DIMACS standard format from the istream. The
+ format is the following. The first line must be
+
+ p min <number of nodes> <number of arcs>
+
+ Then the node definition lines must be found, in the form
+
+ n <node number> <node supply>
+
+ Not all nodes need have a node definition line; these are given zero
+ supply, i.e., they are transhipment nodes (supplies are the inverse of
+ deficits, i.e., a node with positive supply is a source node). Finally,
+ the arc definition lines must be found, in the form
+
+ a <start node> <end node> <lower bound> <upper bound> <flow cost>
+
+ There must be exactly <number of arcs> arc definition lines in the file.
+
+ This method is *not* pure virtual because an implementation is provided by
+ the base class, using the LoadNet() method (which *is* pure virtual).
+ However, the method *is* virtual to allow derived classes to implement
+ more efficient versions, should they have any reason to do so.
+
+ \note Actually, the file format accepted by LoadDMX (at least in the
+ base class implementation) is more general than the DIMACS standard
+ format, in that it is allowed to mix node and arc definitions in
+ any order, while the DIMACS file requires all node information to
+ appear before all arc information.
+
+ \note Other than for the above, this method is assumed to allow for
+ *quadratic* Dimacs files, encoding for convex quadratic separable
+ Min Cost Flow instances. This is a simple extension where each arc
+ descriptor has a sixth field, <quadratic cost>. The provided
+ istream is assumed to be quadratic Dimacs file if IsQuad is true,
+ and a regular linear Dimacs file otherwise. */
+
+/*--------------------------------------------------------------------------*/
+
+ virtual void PreProcess( void ) {}
+
+/**< Extract a smaller/easier equivalent MCF problem. The data of the instance
+ is changed and the easier one is solved instead of the original one. In the
+ MCF case, preprocessing may involve reducing bounds, identifying
+ disconnected components of the graph etc. However, proprocessing is
+ solver-specific.
+
+ This method can be implemented by derived classes in their solver-specific
+ way. Preprocessing may reveal unboundedness or unfeasibility of the
+ problem; if that happens, PreProcess() should properly set the `status'
+ field, that can then be read with MCFGetStatus() [see below].
+
+ Note that preprocessing may destroy all the solution information. Also, it
+ may be allowed to change the data of the problem, such as costs/capacities
+ of the arcs.
+
+ A valid preprocessing is doing nothing, and that's what the default
+ implementation of this method (that is *not* pure virtual) does. */
+
+/*--------------------------------------------------------------------------*/
+
+ virtual inline void SetPar( int par , int val );
+
+/**< Set integer parameters of the algorithm.
+
+ @param par is the parameter to be set; the enum MCFParam can be used, but
+ 'par' is an int (every enum is an int) so that the method can
+ be extended by derived classes for the setting of their
+ parameters
+
+ @param value is the value to assign to the parameter.
+
+ The base class implementation handles these parameters:
+
+ - kMaxIter: the max number of iterations in which the MCF Solver can find
+ an optimal solution (default 0, which means no limit)
+
+ - kReopt: tells the solver if it has to reoptimize. The implementation in
+ the base class sets a flag, the protected \c bool field \c
+ Senstv; if true (default) this field instructs the MCF solver to
+ to try to exploit the information about the latest optimal
+ solution to speedup the optimization of the current problem,
+ while if the field is false the MCF solver should restart the
+ optimization "from scratch" discarding any previous information.
+ Usually reoptimization speeds up the computation considerably,
+ but this is not always true, especially if the data of the
+ problem changes a lot.*/
+
+/*--------------------------------------------------------------------------*/
+
+ virtual inline void SetPar( int par , double val );
+
+/**< Set float parameters of the algorithm.
+
+ @param par is the parameter to be set; the enum MCFParam can be used, but
+ 'par' is an int (every enum is an int) so that the method can
+ be extended by derived classes for the setting of their
+ parameters
+
+ @param value is the value to assign to the parameter.
+
+ The base class implementation handles these parameters:
+
+ - kEpsFlw: sets the tolerance for controlling if the flow on an arc is zero
+ to val. This also sets the tolerance for controlling if a node
+ deficit is zero (see kEpsDfct) to val * < max number of nodes >;
+ this value should be safe for graphs in which any node has less
+ than < max number of nodes > adjacent nodes, i.e., for all graphs
+ but for very dense ones with "parallel arcs"
+
+ - kEpsDfct: sets the tolerance for controlling if a node deficit is zero to
+ val, in case a better value than that autmatically set by
+ kEpsFlw (see above) is available (e.g., val * k would be good
+ if no node has more than k neighbours)
+
+ - kEpsCst: sets the tolerance for controlling if the reduced cost of an arc
+ is zero to val. A feasible solution satisfying eps-complementary
+ slackness, i.e., such that
+ \f[
+ RC[ i , j ] < - eps \Rightarrow X[ i , j ] = U[ ij ]
+ \f]
+ and
+ \f[
+ RC[ i , j ] > eps \Rightarrow X[ i , j ] == 0 ,
+ \f]
+ is known to be ( eps * n )-optimal.
+
+ - kMaxTime: sets the max time (in seconds) in which the MCF Solver can find
+ an optimal solution (default 0, which means no limit). */
+
+/*--------------------------------------------------------------------------*/
+
+ virtual inline void GetPar( int par , int &val );
+
+/**< This method returns one of the integer parameter of the algorithm.
+
+ @param par is the parameter to return [see SetPar( int ) for comments];
+
+ @param val upon return, it will contain the value of the parameter.
+
+ The base class implementation handles the parameters kMaxIter and kReopt.
+ */
+
+/*--------------------------------------------------------------------------*/
+
+ virtual inline void GetPar( int par , double &val );
+
+/**< This method returns one of the integer parameter of the algorithm.
+
+ @param par is the parameter to return [see SetPar( double ) for comments];
+
+ @param val upon return, it will contain the value of the parameter.
+
+ The base class implementation handles the parameters kEpsFlw, kEpsDfct,
+ kEpsCst, and kMaxTime. */
+
+/*--------------------------------------------------------------------------*/
+
+ virtual void SetMCFTime( bool TimeIt = true )
+ {
+ if( TimeIt )
+ if( MCFt )
+ MCFt->ReSet();
+ else
+ MCFt = new OPTtimers();
+ else {
+ delete MCFt;
+ MCFt = NULL;
+ }
+ }
+
+/**< Allocate an OPTtimers object [see OPTtypes.h] to be used for timing the
+ methods of the class. The time can be read with TimeMCF() [see below]. By
+ default, or if SetMCFTime( false ) is called, no timing is done. Note that,
+ since all the relevant methods ot the class are pure virtual, MCFClass can
+ only manage the OPTtimers object, but it is due to derived classes to
+ actually implement the timing.
+
+ @note time accumulates over the calls: calling SetMCFTime(), however,
+ resets the counters, allowing to time specific groups of calls.
+
+ @note of course, setting kMaxTime [see SetPar() above] to any nonzero
+ value has no effect unless SetMCFTime( true ) has been called. */
+
+/*@} -----------------------------------------------------------------------*/
+/*---------------------- METHODS FOR SOLVING THE PROBLEM -------------------*/
+/*--------------------------------------------------------------------------*/
+/** @name Solving the problem
+ @{ */
+
+ virtual void SolveMCF( void ) = 0;
+
+/**< Solver of the Min Cost Flow Problem. Attempts to solve the MCF instance
+ currently loaded in the object. */
+
+/*--------------------------------------------------------------------------*/
+
+ inline int MCFGetStatus( void )
+ {
+ return( status );
+ }
+
+/**< Returns an int describing the current status of the MCF solver. Possible
+ return values are:
+
+ - kUnSolved SolveMCF() has not been called yet, or the data of the
+ problem has been changed since the last call;
+
+ - kOK optimization has been carried out succesfully;
+
+ - kStopped optimization have been stopped before that the stopping
+ conditions of the solver applied, e.g. because of the
+ maximum allowed number of "iterations" [see SetPar( int )]
+ or the maximum allowed time [see SetPar( double )] has been
+ reached; this is not necessarily an error, as it might just
+ be required to re-call SolveMCF() giving it more "resources"
+ in order to solve the problem;
+
+ - kUnfeasible if the current MCF instance is (primal) unfeasible;
+
+ - kUnbounded if the current MCF instance is (primal) unbounded (this can
+ only happen if the solver actually allows F_INF capacities,
+ which is nonstandard in the interface);
+
+ - kError if there was an error during the optimization; this
+ typically indicates that computation cannot be resumed,
+ although solver-dependent ways of dealing with
+ solver-dependent errors may exist.
+
+ MCFClass has a protected \c int \c member \c status \c that can be used by
+ derived classes to hold status information and that is returned by the
+ standard implementation of this method. Note that \c status \c is an
+ \c int \c and not an \c enum \c, and that an \c int \c is returned by this
+ method, in order to allow the derived classes to extend the set of return
+ values if they need to do so. */
+
+/*@} -----------------------------------------------------------------------*/
+/*---------------------- METHODS FOR READING RESULTS -----------------------*/
+/*--------------------------------------------------------------------------*/
+/** @name Reading flow solution
+ @{ */
+
+ virtual void MCFGetX( FRow F , Index_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
+
+/**< Write the optimal flow solution in the vector F[]. If nms == NULL, F[]
+ will be in "dense" format, i.e., the flow relative to arc `i'
+ (i in 0 .. m - 1) is written in F[ i ]. If nms != NULL, F[] will be in
+ "sparse" format, i.e., the indices of the nonzero elements in the flow
+ solution are written in nms (that is then Inf<Index>()-terminated) and
+ the flow value of arc nms[ i ] is written in F[ i ]. Note that nms is
+ *not* guaranteed to be ordered. Also, note that, unlike MCFGetRC() and
+ MCFGetPi() [see below], nms is an *output* of the method.
+
+ The parameters `strt' and `stp' allow to restrict the output of the method
+ to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual cFRow MCFGetX( void )
+ {
+ return( NULL );
+ }
+
+/**< Return a read-only pointer to an internal data structure containing the
+ flow solution in "dense" format. Since this may *not always be available*,
+ depending on the implementation, this method can (uniformly) return NULL.
+ This is done by the base class already, so a derived class that does not
+ have the information ready does not need to implement the method. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual bool HaveNewX( void )
+ {
+ return( false );
+ }
+
+/**< Return true if a different (approximately) optimal primal solution is
+ available. If the method returns true, then any subsequent call to (any
+ form of) MCFGetX() will return a different primal solution w.r.t. the one
+ that was being returned *before* the call to HaveNewX(). This solution need
+ not be optimal (although, ideally, it has to be "good); this can be checked
+ by comparing its objective function value, that will be returned by a call
+ to MCFGetFO() [see below].
+
+ Any subsequent call of HaveNewX() that returns true produces a new
+ solution, until the first that returns false; from then on, no new
+ solutions will be generated until something changes in the problem's
+ data.
+
+ Note that a default implementation of HaveNewX() is provided which is
+ good for those solvers that only produce one optimal primal solution. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Reading potentials
+ @{ */
+
+ virtual void MCFGetPi( CRow P , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
+
+/**< Writes the optimal node potentials in the vector P[]. If nms == NULL,
+ the node potential of node `i' (i in 0 .. n - 1) is written in P[ i ]
+ (note that here node names always start from zero, regardless to the value
+ of USENAME0). If nms != NULL, it must point to a vector of indices in
+ 0 .. n - 1 (ordered in increasing sense and Inf<Index>()-terminated), and
+ the node potential of nms[ i ] is written in P[ i ]. Note that, unlike
+ MCFGetX() above, nms is an *input* of the method.
+
+ The parameters `strt' and `stp' allow to restrict the output of the method
+ to all and only the nodes `i' with strt <= i < min( MCFn() , stp ). `strt'
+ and `stp' work in "&&" with nms; that is, if nms != NULL then only the
+ values corresponding to nodes which are *both* in nms[] and whose index is
+ in the correct range are returned. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual cCRow MCFGetPi( void )
+ {
+ return( NULL );
+ }
+
+/**< Return a read-only pointer to an internal data structure containing the
+ node potentials. Since this may *not always be available*, depending on
+ the implementation, this method can (uniformly) return NULL.
+ This is done by the base class already, so a derived class that does not
+ have the information ready does not need to implement the method. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual bool HaveNewPi( void )
+ {
+ return( false );
+ }
+
+/**< Return true if a different (approximately) optimal dual solution is
+ available. If the method returns true, then any subsequent call to (any
+ form of) MCFGetPi() will return a different dual solution, and MCFGetRC()
+ [see below] will return the corresponding reduced costs. The new solution
+ need not be optimal (although, ideally, it has to be "good); this can be
+ checked by comparing its objective function value, that will be returned
+ by a call to MCFGetDFO() [see below].
+
+ Any subsequent call of HaveNewPi() that returns true produces a new
+ solution, until the first that returns false; from then on, no new
+ solutions will be generated until something changes in the problem's
+ data.
+
+ Note that a default implementation of HaveNewPi() is provided which is
+ good for those solvers that only produce one optimal dual solution. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Reading reduced costs
+ @{ */
+
+ virtual void MCFGetRC( CRow CR , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
+
+/**< Write the reduced costs corresponding to the current dual solution in
+ RC[]. If nms == NULL, the reduced cost of arc `i' (i in 0 .. m - 1) is
+ written in RC[ i ]; if nms != NULL, it must point to a vector of indices
+ in 0 .. m - 1 (ordered in increasing sense and Inf<Index>()-terminated),
+ and the reduced cost of arc nms[ i ] is written in RC[ i ]. Note that,
+ unlike MCFGetX() above, nms is an *input* of the method.
+
+ The parameters `strt' and `stp' allow to restrict the output of the method
+ to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
+ and `stp' work in "&&" with nms; that is, if nms != NULL then only the
+ values corresponding to arcs which are *both* in nms[] and whose index is
+ in the correct range are returned.
+
+ @note the output of MCFGetRC() will change after any call to HaveNewPi()
+ [see above] which returns true. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual cCRow MCFGetRC( void )
+ {
+ return( NULL );
+ }
+
+/**< Return a read-only pointer to an internal data structure containing the
+ reduced costs. Since this may *not always be available*, depending on the
+ implementation, this method can (uniformly) return NULL.
+ This is done by the base class already, so a derived class that does not
+ have the information ready does not need to implement the method.
+
+ @note the output of MCFGetRC() will change after any call to HaveNewPi()
+ [see above] which returns true. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual CNumber MCFGetRC( cIndex i ) = 0;
+
+/**< Return the reduced cost of the i-th arc. This information should be
+ cheapily available in most implementations.
+
+ @note the output of MCFGetRC() will change after any call to HaveNewPi()
+ [see above] which returns true. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Reading the objective function value
+ @{ */
+
+ virtual FONumber MCFGetFO( void ) = 0;
+
+/**< Return the objective function value of the primal solution currently
+ returned by MCFGetX().
+
+ If MCFGetStatus() == kOK, this is guaranteed to be the optimal objective
+ function value of the problem (to within the optimality tolerances), but
+ only prior to any call to HaveNewX() that returns true. MCFGetFO()
+ typically returns Inf<FONumber>() if MCFGetStatus() == kUnfeasible and
+ - Inf<FONumber>() if MCFGetStatus() == kUnbounded. If MCFGetStatus() ==
+ kStopped and MCFGetFO() returns a finite value, it must be an upper bound
+ on the optimal objective function value (typically, the objective function
+ value of one primal feasible solution). */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual FONumber MCFGetDFO( void )
+ {
+ switch( MCFGetStatus() ) {
+ case( kUnSolved ):
+ case( kStopped ):
+ case( kError ): return( -Inf<FONumber>() );
+ default: return( MCFGetFO() );
+ }
+ }
+
+/**< Return the objective function value of the dual solution currently
+ returned by MCFGetPi() / MCFGetRC(). This value (possibly) changes after
+ any call to HaveNewPi() that returns true. The relations between
+ MCFGetStatus() and MCFGetDFO() are analogous to these of MCFGetFO(),
+ except that a finite value corresponding to kStopped must be a lower
+ bound on the optimal objective function value (typically, the objective
+ function value one dual feasible solution).
+
+ A default implementation is provided for MCFGetDFO(), which is good for
+ MCF solvers where the primal and dual optimal solution values always
+ are identical (except if the problem is unfeasible/unbounded). */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Getting unfeasibility certificate
+ @{ */
+
+ virtual FNumber MCFGetUnfCut( Index_Set Cut )
+ {
+ *Cut = Inf<Index>();
+ return( 0 );
+ }
+
+/**< Return an unfeasibility certificate. In an unfeasible MCF problem,
+ unfeasibility can always be reduced to the existence of a cut (subset of
+ nodes of the graph) such as either:
+
+ - the inverse of the deficit of the cut (the sum of all the deficits of the
+ nodes in the cut) is larger than the forward capacity of the cut (sum of
+ the capacities of forward arcs in the cut); that is, the nodes in the
+ cut globally produce more flow than can be routed to sinks outside the
+ cut;
+
+ - the deficit of the cut is larger than the backward capacity of the cut
+ (sum of the capacities of backward arcs in the cut); that is, the nodes
+ in the cut globally require more flow than can be routed to them from
+ sources outside the cut.
+
+ When detecting unfeasibility, MCF solvers are typically capable to provide
+ one such cut. This information can be useful - typically, the only way to
+ make the problem feasible is to increase the capacity of at least one of
+ the forward/backward arcs of the cut -, and this method is provided for
+ getting it. It can be called only if MCFGetStatus() == kUnfeasible, and
+ should write in Cut the set of names of nodes in the unfeasible cut (note
+ that node names depend on USENAME0), Inf<Index>()-terminated, returning the
+ deficit of the cut (which allows to distinguish which of the two cases
+ above hold). In general, no special properties can be expected from the
+ returned cut, but solvers should be able to provide e.g. "small" cuts.
+
+ However, not all solvers may be (easily) capable of providing this
+ information; thus, returning 0 (no cut) is allowed, as in the base class
+ implementation, to signify that this information is not available. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Getting unboundedness certificate
+ @{ */
+
+ virtual Index MCFGetUnbCycl( Index_Set Pred , Index_Set ArcPred )
+ {
+ return( Inf<Index>() );
+ }
+
+/**< Return an unboundedness certificate. In an unbounded MCF problem,
+ unboundedness can always be reduced to the existence of a directed cycle
+ with negative cost and all arcs having infinite capacity.
+ When detecting unboundedness, MCF solvers are typically capable to provide
+ one such cycle. This information can be useful, and this method is provided
+ for getting it. It can be called only if MCFGetStatus() == kUnbounded, and
+ writes in Pred[] and ArcPred[], respectively, the node and arc predecessor
+ function of the cycle. That is, if node `i' belongs to the cycle then
+ `Pred[ i ]' contains the name of the predecessor of `j' of `i' in the cycle
+ (note that node names depend on USENAME0), and `ArcPred[ i ]' contains the
+ index of the arc joining the two (note that in general there may be
+ multiple copies of each arc). Entries of the vectors for nodes not
+ belonging to the cycle are in principle undefined, and the name of one node
+ belonging to the cycle is returned by the method.
+ Note that if there are multiple cycles with negative costs this method
+ will return just one of them (finding the cycle with most negative cost
+ is an NO-hard problem), although solvers should be able to produce cycles
+ with "large negative" cost.
+
+ However, not all solvers may be (easily) capable of providing this
+ information; thus, returning Inf<Index>() is allowed, as in the base class
+ implementation, to signify that this information is not available. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Saving/restoring the state of the solver
+ @{ */
+
+ virtual MCFStatePtr MCFGetState( void )
+ {
+ return( NULL );
+ }
+
+ /**< Save the state of the MCF solver. The MCFClass interface supports the
+ notion of saving and restoring the state of the MCF solver, such as the
+ current/optimal basis in a simplex solver. The "empty" class MCFState is
+ defined as a placeholder for state descriptions.
+
+ MCFGetState() creates and returns a pointer to an object of (a proper
+ derived class of) class MCFState which describes the current state of the
+ MCF solver. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void MCFPutState( MCFStatePtr S ) {}
+
+/**< Restore the solver to the state in which it was when the state `S' was
+ created with MCFGetState() [see above].
+
+ The typical use of this method is the following: a MCF problem is solved
+ and the "optimal state" is set aside. Then the problem changes and it is
+ re-solved. Then, the problem has to be changed again to a form that is
+ close to the original one but rather different from the second one (think
+ of a long backtracking in a Branch & Bound) to which the current "state"
+ referes. Then, the old optimal state can be expected to provide a better
+ starting point for reoptimization [see ReOptimize() below].
+
+ Note, however, that the state is only relative to the optimization process,
+ i.e., this operation is meaningless if the data of the problem has changed
+ in the meantime. So, if a state has to be used for speeding up
+ reoptimization, the following has to be done:
+
+ - first, the data of the solver is brought back to *exactly* the same as
+ it was at the moment where the state `S' was created (prior than this
+ operation a ReOptimize( false ) call is probably advisable);
+
+ - then, MCFPutState() is called (and ReOptimize( true ) is called);
+
+ - only afterwards the data of the problem is changed to the final state
+ and the problem is solved.
+
+ A "put state" operation does not "deplete" the state, which can therefore
+ be used more than once. Indeed, a state is constructed inside the solver
+ for each call to MCFGetState(), but the solver never deletes statuses;
+ this has to be done on the outside when they are no longer needed (the
+ solver must be "resistent" to deletion of the state at any moment).
+
+ Since not all the MCF solvers reoptimize (efficiently enough to make these
+ operations worth), an "empty" implementation that does nothing is provided
+ by the base class. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Time the code
+ @{ */
+
+ void TimeMCF( double &t_us , double &t_ss )
+ {
+ t_us = t_ss = 0;
+ if( MCFt )
+ MCFt->Read( t_us , t_ss );
+ }
+
+/**< Time the code. If called within any of the methods of the class that are
+ "actively timed" (this depends on the subclasses), this method returns the
+ user and sistem time (in seconds) since the start of that method. If
+ methods that are actively timed call other methods that are actively
+ timed, TimeMCF() returns the (...) time since the beginning of the
+ *outer* actively timed method. If called outside of any actively timed
+ method, this method returns the (...) time spent in all the previous
+ executions of all the actively timed methods of the class.
+
+ Implementing the proper calls to MCFt->Start() and MCFt->Stop() is due to
+ derived classes; these should at least be placed at the beginning and at
+ the end, respectively, of SolveMCF() and presumably the Chg***() methods,
+ that is, at least these methods should be "actively timed". */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ double TimeMCF( void )
+ {
+ return( MCFt ? MCFt->Read() : 0 );
+ }
+
+/**< Like TimeMCF(double,double) [see above], but returns the total time. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Check the solutions
+ @{ */
+
+ inline void CheckPSol( void );
+
+/**< Check that the primal solution returned by the solver is primal feasible.
+ (to within the tolerances set by SetPar(kEps****) [see above], if any). Also,
+ check that the objective function value is correct.
+
+ This method is implemented by the base class, using the above methods
+ for collecting the solutions and the methods of the next section for
+ reading the data of the problem; as such, they will work for any derived
+ class that properly implements all these methods. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ inline void CheckDSol( void );
+
+/**< Check that the dual solution returned by the solver is dual feasible.
+ (to within the tolerances set by SetPar(kEps****) [see above], if any). Also,
+ check that the objective function value is correct.
+
+ This method is implemented by the base class, using the above methods
+ for collecting the solutions and the methods of the next section for
+ reading the data of the problem; as such, they will work for any derived
+ class that properly implements all these methods. */
+
+/*@} -----------------------------------------------------------------------*/
+/*-------------- METHODS FOR READING THE DATA OF THE PROBLEM ---------------*/
+/*--------------------------------------------------------------------------*/
+/** @name Reading graph size
+ @{ */
+
+ inline Index MCFnmax( void )
+ {
+ return( nmax );
+ }
+
+/**< Return the maximum number of nodes for this instance of MCFClass.
+ The implementation of the method in the base class returns the protected
+ fields \c nmax, which is provided for derived classes to hold this
+ information. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ inline Index MCFmmax( void )
+ {
+ return( mmax );
+ }
+
+/**< Return the maximum number of arcs for this instance of MCFClass.
+ The implementation of the method in the base class returns the protected
+ fields \c mmax, which is provided for derived classes to hold this
+ information. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ inline Index MCFn( void )
+ {
+ return( n );
+ }
+
+/**< Return the number of nodes in the current graph.
+ The implementation of the method in the base class returns the protected
+ fields \c n, which is provided for derived classes to hold this
+ information. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ inline Index MCFm( void )
+ {
+ return( m );
+ }
+
+/**< Return the number of arcs in the current graph.
+ The implementation of the method in the base class returns the protected
+ fields \c m, which is provided for derived classes to hold this
+ information. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Reading graph topology
+ @{ */
+
+ virtual void MCFArcs( Index_Set Startv , Index_Set Endv ,
+ cIndex_Set nms = NULL , cIndex strt = 0 ,
+ Index stp = Inf<Index>() ) = 0;
+
+/**< Write the starting (tail) and ending (head) nodes of the arcs in Startv[]
+ and Endv[]. If nms == NULL, then the information relative to all arcs is
+ written into Startv[] and Endv[], otherwise Startv[ i ] and Endv[ i ]
+ contain the information relative to arc nms[ i ] (nms[] must be
+ Inf<Index>()-terminated).
+
+ The parameters `strt' and `stp' allow to restrict the output of the method
+ to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
+ and `stp' work in "&&" with nms; that is, if nms != NULL then only the
+ values corresponding to arcs which are *both* in nms and whose index is in
+ the correct range are returned.
+
+ Startv or Endv can be NULL, meaning that only the other information is
+ required.
+
+ @note If USENAME0 == 0 then the returned node names will be in the range
+ 1 .. n, while if USENAME0 == 1 the returned node names will be in
+ the range 0 .. n - 1.
+
+ @note If the graph is "dynamic", be careful to use MCFn() e MCFm() to
+ properly choose the dimension of nodes and arcs arrays. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual Index MCFSNde( cIndex i ) = 0;
+
+/**< Return the starting (tail) node of the arc `i'.
+
+ @note If USENAME0 == 0 then the returned node names will be in the range
+ 1 .. n, while if USENAME0 == 1 the returned node names will be in
+ the range 0 .. n - 1. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual Index MCFENde( cIndex i ) = 0;
+
+/**< Return the ending (head) node of the arc `i'.
+
+ @note If USENAME0 == 0 then the returned node names will be in the range
+ 1 .. n, while if USENAME0 == 1 the returned node names will be in
+ the range 0 .. n - 1. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual cIndex_Set MCFSNdes( void )
+ {
+ return( NULL );
+ }
+
+/**< Return a read-only pointer to an internal vector containing the starting
+ (tail) nodes for each arc. Since this may *not always be available*,
+ depending on the implementation, this method can (uniformly) return NULL.
+ This is done by the base class already, so a derived class that does not
+ have the information ready does not need to implement the method. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual cIndex_Set MCFENdes( void )
+ {
+ return( NULL );
+ }
+
+/**< Return a read-only pointer to an internal vector containing the ending
+ (head) nodes for each arc. Since this may *not always be available*,
+ depending on the implementation, this method can (uniformly) return NULL.
+ This is done by the base class already, so a derived class that does not
+ have the information ready does not need to implement the method. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Reading arc costs
+ @{ */
+
+ virtual void MCFCosts( CRow Costv , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
+
+/**< Write the arc costs into Costv[]. If nms == NULL, then all the costs are
+ written, otherwise Costv[ i ] contains the information relative to arc
+ nms[ i ] (nms[] must be Inf<Index>()-terminated).
+
+ The parameters `strt' and `stp' allow to restrict the output of the method
+ to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
+ and `stp' work in "&&" with nms; that is, if nms != NULL then only the
+ values corresponding to arcs which are *both* in nms and whose index is in
+ the correct range are returned. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual CNumber MCFCost( cIndex i ) = 0;
+
+/**< Return the cost of the i-th arc. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual cCRow MCFCosts( void )
+ {
+ return( NULL );
+ }
+
+/**< Return a read-only pointer to an internal vector containing the arc
+ costs. Since this may *not always be available*, depending on the
+ implementation, this method can (uniformly) return NULL.
+ This is done by the base class already, so a derived class that does not
+ have the information ready does not need to implement the method. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void MCFQCoef( CRow Qv , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() )
+
+/**< Write the quadratic coefficients of the arc costs into Qv[]. If
+ nms == NULL, then all the costs are written, otherwise Costv[ i ] contains
+ the information relative to arc nms[ i ] (nms[] must be
+ Inf<Index>()-terminated).
+
+ The parameters `strt' and `stp' allow to restrict the output of the method
+ to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
+ and `stp' work in "&&" with nms; that is, if nms != NULL then only the
+ values corresponding to arcs which are *both* in nms and whose index is in
+ the correct range are returned.
+
+ Note that the method is *not* pure virtual: an implementation is provided
+ for "pure linear" MCF solvers that only work with all zero quadratic
+ coefficients. */
+ {
+ if( nms ) {
+ while( *nms < strt )
+ nms++;
+
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(Qv++) = 0;
+ }
+ else {
+ if( stp > m )
+ stp = m;
+
+ while( stp-- > strt )
+ *(Qv++) = 0;
+ }
+ }
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual CNumber MCFQCoef( cIndex i )
+
+/**< Return the quadratic coefficients of the cost of the i-th arc. Note that
+ the method is *not* pure virtual: an implementation is provided for "pure
+ linear" MCF solvers that only work with all zero quadratic coefficients. */
+ {
+ return( 0 );
+ }
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual cCRow MCFQCoef( void )
+ {
+ return( NULL );
+ }
+
+/**< Return a read-only pointer to an internal vector containing the arc
+ costs. Since this may *not always be available*, depending on the
+ implementation, this method can (uniformly) return NULL.
+ This is done by the base class already, so a derived class that does not
+ have the information ready (such as "pure linear" MCF solvers that only
+ work with all zero quadratic coefficients) does not need to implement the
+ method. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Reading arc capacities
+ @{ */
+
+ virtual void MCFUCaps( FRow UCapv , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
+
+/**< Write the arc capacities into UCapv[]. If nms == NULL, then all the
+ capacities are written, otherwise UCapv[ i ] contains the information
+ relative to arc nms[ i ] (nms[] must be Inf<Index>()-terminated).
+
+ The parameters `strt' and `stp' allow to restrict the output of the method
+ to all and only the arcs `i' with strt <= i < min( MCFm() , stp ). `strt'
+ and `stp' work in "&&" with nms; that is, if nms != NULL then only the
+ values corresponding to arcs which are *both* in nms and whose index is in
+ the correct range are returned. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual FNumber MCFUCap( cIndex i ) = 0;
+
+/**< Return the capacity of the i-th arc. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual cFRow MCFUCaps( void )
+ {
+ return( NULL );
+ }
+
+/**< Return a read-only pointer to an internal vector containing the arc
+ capacities. Since this may *not always be available*, depending on the
+ implementation, this method can (uniformly) return NULL.
+ This is done by the base class already, so a derived class that does not
+ have the information ready does not need to implement the method. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Reading node deficits
+ @{ */
+
+ virtual void MCFDfcts( FRow Dfctv , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
+
+/**< Write the node deficits into Dfctv[]. If nms == NULL, then all the
+ defcits are written, otherwise Dfctvv[ i ] contains the information
+ relative to node nms[ i ] (nms[] must be Inf<Index>()-terminated).
+
+ The parameters `strt' and `stp' allow to restrict the output of the method
+ to all and only the nodes `i' with strt <= i < min( MCFn() , stp ). `strt'
+ and `stp' work in "&&" with nms; that is, if nms != NULL then only the
+ values corresponding to nodes which are *both* in nms and whose index is in
+ the correct range are returned.
+
+ @note Here node "names" (strt and stp, those contained in nms[] or `i'
+ in MCFDfct()) go from 0 to n - 1, regardless to the value of
+ USENAME0; hence, if USENAME0 == 0 then the first node is "named 1",
+ but its deficit is returned by MCFDfcts( Dfctv , NULL , 0 , 1 ). */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual FNumber MCFDfct( cIndex i ) = 0;
+
+/**< Return the deficit of the i-th node.
+
+ @note Here node "names" go from 0 to n - 1, regardless to the value of
+ USENAME0; hence, if USENAME0 == 0 then the first node is "named 1",
+ but its deficit is returned by MCFDfct( 0 ). */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual cFRow MCFDfcts( void )
+ {
+ return( NULL );
+ }
+
+/**< Return a read-only pointer to an internal vector containing the node
+ deficits. Since this may *not always be available*, depending on the
+ implementation, this method can (uniformly) return NULL.
+ This is done by the base class already, so a derived class that does not
+ have the information ready does not need to implement the method.
+
+ @note Here node "names" go from 0 to n - 1, regardless to the value of
+ USENAME0; hence, if USENAME0 == 0 then the first node is "named 1",
+ but its deficit is contained in MCFDfcts()[ 0 ]. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Write problem to file
+ @{ */
+
+ virtual inline void WriteMCF( ostream &oStrm , int frmt = 0 );
+
+/**< Write the current MCF problem to an ostream. This may be useful e.g. for
+ debugging purposes.
+
+ The base MCFClass class provides output in two different formats, depending
+ on the value of the parameter frmt:
+
+ - kDimacs the problem is written in DIMACS standard format, read by most
+ MCF codes available;
+
+ - kMPS the problem is written in the "modern version" (tab-separated)
+ of the MPS format, read by most LP/MIP solvers;
+
+ - kFWMPS the problem is written in the "old version" (fixed width
+ fields) of the MPS format; this is read by most LP/MIP
+ solvers, but some codes still require the old format.
+
+ The implementation of WriteMCF() in the base class uses all the above
+ methods for reading the data; as such it will work for any derived class
+ that properly implements this part of the interface, but it may not be
+ very efficient. Thus, the method is virtual to allow the derived classes
+ to either re-implement WriteMCF() for the above two formats in a more
+ efficient way, and/or to extend it to support other solver-specific
+ formats.
+
+ @note None of the above two formats supports quadratic MCFs, so if
+ nonzero quadratic coefficients are present, they are just ignored.
+ */
+
+/*@} -----------------------------------------------------------------------*/
+/*------------- METHODS FOR ADDING / REMOVING / CHANGING DATA --------------*/
+/*--------------------------------------------------------------------------*/
+/** @name Changing the costs
+ @{ */
+
+ virtual void ChgCosts( cCRow NCost , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
+
+/**< Change the arc costs. In particular, change the costs that are:
+
+ - listed in into the vector of indices `nms' (ordered in increasing
+ sense and Inf<Index>()-terminated),
+
+ - *and* whose name belongs to the interval [`strt', `stp').
+
+ That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th cost will be
+ changed to NCost[ i ]. If nms == NULL (as the default), *all* the
+ entries in the given range will be changed; if stp > MCFm(), then the
+ smaller bound is used.
+
+ @note changing the costs of arcs that *do not exist* is *not allowed*;
+ only arcs which have not been closed/deleted [see CloseArc() /
+ DelArc() below and LoadNet() above about C_INF costs] can be
+ touched with these methods. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void ChgCost( Index arc , cCNumber NCost ) = 0;
+
+/**< Change the cost of the i-th arc.
+
+ @note changing the costs of arcs that *do not exist* is *not allowed*;
+ only arcs which have not been closed/deleted [see CloseArc() /
+ DelArc() below and LoadNet() above about C_INF costs] can be
+ touched with these methods. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void ChgQCoef( cCRow NQCoef = NULL , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() )
+
+/**< Change the quadratic coefficients of the arc costs. In particular,
+ change the coefficients that are:
+
+ - listed in into the vector of indices `nms' (ordered in increasing
+ sense and Inf<Index>()-terminated),
+
+ - *and* whose name belongs to the interval [`strt', `stp').
+
+ That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th cost will be
+ changed to NCost[ i ]. If nms == NULL (as the default), *all* the
+ entries in the given range will be changed; if stp > MCFm(), then the
+ smaller bound is used. If NQCoef == NULL, all the specified coefficients
+ are set to zero.
+
+ Note that the method is *not* pure virtual: an implementation is provided
+ for "pure linear" MCF solvers that only work with all zero quadratic
+ coefficients.
+
+ @note changing the costs of arcs that *do not exist* is *not allowed*;
+ only arcs which have not been closed/deleted [see CloseArc() /
+ DelArc() below and LoadNet() above about C_INF costs] can be
+ touched with these methods. */
+ {
+ assert(NQCoef);
+ }
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void ChgQCoef( Index arc , cCNumber NQCoef )
+
+/**< Change the quadratic coefficient of the cost of the i-th arc.
+
+ Note that the method is *not* pure virtual: an implementation is provided
+ for "pure linear" MCF solvers that only work with all zero quadratic
+ coefficients.
+
+ @note changing the costs of arcs that *do not exist* is *not allowed*;
+ only arcs which have not been closed/deleted [see CloseArc() /
+ DelArc() below and LoadNet() above about C_INF costs] can be
+ touched with these methods. */
+ {
+ assert(NQCoef);
+ }
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Changing the capacities
+ @{ */
+
+ virtual void ChgUCaps( cFRow NCap , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
+
+/**< Change the arc capacities. In particular, change the capacities that are:
+
+ - listed in into the vector of indices `nms' (ordered in increasing sense
+ and Inf<Index>()-terminated),
+
+ - *and* whose name belongs to the interval [`strt', `stp').
+
+ That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th capacity will be
+ changed to NCap[ i ]. If nms == NULL (as the default), *all* the
+ entries in the given range will be changed; if stp > MCFm(), then the
+ smaller bound is used.
+
+ @note changing the capacities of arcs that *do not exist* is *not allowed*;
+ only arcs that have not been closed/deleted [see CloseArc() /
+ DelArc() below and LoadNet() above about C_INF costs] can be
+ touched with these methods. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void ChgUCap( Index arc , cFNumber NCap ) = 0;
+
+/**< Change the capacity of the i-th arc.
+
+ @note changing the capacities of arcs that *do not exist* is *not allowed*;
+ only arcs that have not been closed/deleted [see CloseArc() /
+ DelArc() below and LoadNet() above about C_INF costs] can be
+ touched with these methods. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Changing the deficits
+ @{ */
+
+ virtual void ChgDfcts( cFRow NDfct , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) = 0;
+
+/**< Change the node deficits. In particular, change the deficits that are:
+
+ - listed in into the vector of indices `nms' (ordered in increasing sense
+ and Inf<Index>()-terminated),
+
+ - *and* whose name belongs to the interval [`strt', `stp').
+
+ That is, if strt <= nms[ i ] < stp, then the nms[ i ]-th deficit will be
+ changed to NDfct[ i ]. If nms == NULL (as the default), *all* the
+ entries in the given range will be changed; if stp > MCFn(), then the
+ smaller bound is used.
+
+ Note that, in ChgDfcts(), node "names" (strt, stp or those contained
+ in nms[]) go from 0 to n - 1, regardless to the value of USENAME0; hence,
+ if USENAME0 == 0 then the first node is "named 1", but its deficit can be
+ changed e.g. with ChgDfcts( &new_deficit , NULL , 0 , 1 ).
+
+ @note changing the capacities of nodes that *do not exist* is *not
+ allowed*; only nodes that have not been deleted [see DelNode()
+ below] can be touched with these methods. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void ChgDfct( Index node , cFNumber NDfct ) = 0;
+
+/**< Change the deficit of the i-th node.
+
+ Note that, in ChgDfct[s](), node "names" (i, strt/ stp or those contained
+ in nms[]) go from 0 to n - 1, regardless to the value of USENAME0; hence,
+ if USENAME0 == 0 then the first node is "named 1", but its deficit can be
+ changed e.g. with ChgDfcts( 0 , new_deficit ).
+
+ @note changing the capacities of nodes that *do not exist* is *not
+ allowed*; only nodes that have not been deleted [see DelNode()
+ below] can be touched with these methods. */
+
+/*@} -----------------------------------------------------------------------*/
+/** @name Changing graph topology
+ @{ */
+
+ virtual void CloseArc( cIndex name ) = 0;
+
+/**< "Close" the arc `name'. Although all the associated information (name,
+ cost, capacity, end and start node) is kept, the arc is removed from the
+ problem until OpenArc( i ) [see below] is called.
+
+ "closed" arcs always have 0 flow, but are otherwise counted as any other
+ arc; for instance, MCFm() does *not* decrease as an effect of a call to
+ CloseArc(). How this closure is implemented is solver-specific. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual bool IsClosedArc( cIndex name ) = 0;
+
+/**< IsClosedArc() returns true if and only if the arc `name' is closed. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void DelNode( cIndex name ) = 0;
+
+/**< Delete the node `name'.
+
+ For any value of `name', all incident arcs to that node are closed [see
+ CloseArc() above] (*not* Deleted, see DelArc() below) and the deficit is
+ set to zero.
+
+ Il furthermore `name' is the last node, the number of nodes as reported by
+ MCFn() is reduced by at least one, until the n-th node is not a deleted
+ one. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void OpenArc( cIndex name ) = 0;
+
+/**< Restore the previously closed arc `name'. It is an error to open an arc
+ that has not been previously closed. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual Index AddNode( cFNumber aDfct ) = 0;
+
+/**< Add a new node with deficit aDfct, returning its name. Inf<Index>() is
+ returned if there is no room for a new node. Remember that the node names
+ are either { 0 .. nmax - 1 } or { 1 .. nmax }, depending on the value of
+ USENAME0. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void ChangeArc( cIndex name , cIndex nSN = Inf<Index>() ,
+ cIndex nEN = Inf<Index>() ) = 0;
+
+/**< Change the starting and/or ending node of arc `name' to nSN and nEN.
+ Each parameter being Inf<Index>() means to leave the previous starting or
+ ending node untouched. When this method is called `name' can be either the
+ name of a "normal" arc or that of a "closed" arc [see CloseArc() above]:
+ in the latter case, at the end of ChangeArc() the arc is *still closed*,
+ and it remains so until OpenArc( name ) [see above] is called. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual void DelArc( cIndex name ) = 0;
+
+/**< Delete the arc `name'. Unlike "closed" arcs, all the information
+ associated with a deleted arc is lost and `name' is made available as a
+ name for new arcs to be created with AddArc() [see below].
+
+ Il furthermore `name' is the last arc, the number of arcs as reported by
+ MCFm() is reduced by at least one, until the m-th arc is not a deleted
+ one. Otherwise, the flow on the arc is always ensured to be 0. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual bool IsDeletedArc( cIndex name ) = 0;
+
+/**< Return true if and only if the arc `name' is deleted. It should only
+ be called with name < MCFm(), as every other arc is deleted by
+ definition. */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ virtual Index AddArc( cIndex Start , cIndex End , cFNumber aU ,
+ cCNumber aC ) = 0;
+
+/**< Add the new arc ( Start , End ) with cost aC and capacity aU,
+ returning its name. Inf<Index>() is returned if there is no room for a new
+ arc. Remember that arc names go from 0 to mmax - 1. */
+
+/*@} -----------------------------------------------------------------------*/
+/*------------------------------ DESTRUCTOR --------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @name Destructor
+ @{ */
+
+ virtual ~MCFClass()
+ {
+ delete MCFt;
+ }
+
+/**< Destructor of the class. The implementation in the base class only
+ deletes the MCFt field. It is virtual, as it should be. */
+
+/*@} -----------------------------------------------------------------------*/
+/*-------------------- PROTECTED PART OF THE CLASS -------------------------*/
+/*--------------------------------------------------------------------------*/
+/*-- --*/
+/*-- The following fields are believed to be general enough to make it --*/
+/*-- worth adding them to the abstract base class. --*/
+/*-- --*/
+/*--------------------------------------------------------------------------*/
+
+ protected:
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- PROTECTED METHODS ----------------------------*/
+/*--------------------------------------------------------------------------*/
+/*-------------------------- MANAGING COMPARISONS --------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @name Managing comparisons.
+ The following methods are provided for making it easier to perform
+ comparisons, with and without tolerances.
+ @{ */
+
+/*--------------------------------------------------------------------------*/
+/** true if flow x is equal to zero (possibly considering tolerances). */
+
+ template<class T>
+ inline bool ETZ( T x , const T eps )
+ {
+ if( numeric_limits<T>::is_integer )
+ return( x == 0 );
+ else
+ return( (x <= eps ) && ( x >= -eps ) );
+ }
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
+/** true if flow x is greater than zero (possibly considering tolerances). */
+
+ template<class T>
+ inline bool GTZ( T x , const T eps )
+ {
+ if( numeric_limits<T>::is_integer )
+ return( x > 0 );
+ else
+ return( x > eps );
+ }
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
+/** true if flow x is greater than or equal to zero (possibly considering
+ tolerances). */
+
+ template<class T>
+ inline bool GEZ( T x , const T eps )
+ {
+ if( numeric_limits<T>::is_integer )
+ return( x >= 0 );
+ else
+ return( x >= - eps );
+ }
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
+/** true if flow x is less than zero (possibly considering tolerances). */
+
+ template<class T>
+ inline bool LTZ( T x , const T eps )
+ {
+ if( numeric_limits<T>::is_integer )
+ return( x < 0 );
+ else
+ return( x < - eps );
+ }
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
+/** true if flow x is less than or equal to zero (possibly considering
+ tolerances). */
+
+ template<class T>
+ inline bool LEZ( T x , const T eps )
+ {
+ if( numeric_limits<T>::is_integer )
+ return( x <= 0 );
+ else
+ return( x <= eps );
+ }
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
+/** true if flow x is greater than flow y (possibly considering tolerances).
+ */
+
+ template<class T>
+ inline bool GT( T x , T y , const T eps )
+ {
+ if( numeric_limits<T>::is_integer )
+ return( x > y );
+ else
+ return( x > y + eps );
+ }
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
+/** true if flow x is less than flow y (possibly considering tolerances). */
+
+ template<class T>
+ inline bool LT( T x , T y , const T eps )
+ {
+ if( numeric_limits<T>::is_integer )
+ return( x < y );
+ else
+ return( x < y - eps );
+ }
+
+/*@} -----------------------------------------------------------------------*/
+/*---------------------- PROTECTED DATA STRUCTURES -------------------------*/
+/*--------------------------------------------------------------------------*/
+
+ Index n; ///< total number of nodes
+ Index nmax; ///< maximum number of nodes
+
+ Index m; ///< total number of arcs
+ Index mmax; ///< maximum number of arcs
+
+ int status; ///< return status, see the comments to MCFGetStatus() above.
+ ///< Note that the variable is defined int to allow derived
+ ///< classes to return their own specialized status codes
+ bool Senstv; ///< true <=> the latest optimal solution should be exploited
+
+ OPTtimers *MCFt; ///< timer for performances evaluation
+
+ FNumber EpsFlw; ///< precision for comparing arc flows / capacities
+ FNumber EpsDfct; ///< precision for comparing node deficits
+ CNumber EpsCst; ///< precision for comparing arc costs
+
+ double MaxTime; ///< max time (in seconds) in which MCF Solver can find
+ ///< an optimal solution (0 = no limits)
+ int MaxIter; ///< max number of iterations in which MCF Solver can find
+ ///< an optimal solution (0 = no limits)
+
+/*--------------------------------------------------------------------------*/
+
+ }; // end( class MCFClass )
+
+/*@} -----------------------------------------------------------------------*/
+/*------------------- inline methods implementation ------------------------*/
+/*--------------------------------------------------------------------------*/
+
+inline void MCFClass::LoadDMX( istream &DMXs , bool IsQuad )
+{
+ // read first non-comment line - - - - - - - - - - - - - - - - - - - - - - -
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ char c;
+ for(;;) {
+ if( ! ( DMXs >> c ) )
+ throw( MCFException( "LoadDMX: error reading the input stream" ) );
+
+ if( c != 'c' ) // if it's not a comment
+ break;
+
+ do // skip the entire line
+ if( ! DMXs.get( c ) )
+ throw( MCFException( "LoadDMX: error reading the input stream" ) );
+ while( c != '\n' );
+ }
+
+ if( c != 'p' )
+ throw( MCFException( "LoadDMX: format error in the input stream" ) );
+
+ char buf[ 3 ]; // free space
+ DMXs >> buf; // skip (in has to be "min")
+
+ Index tn;
+ if( ! ( DMXs >> tn ) )
+ throw( MCFException( "LoadDMX: error reading number of nodes" ) );
+
+ Index tm;
+ if( ! ( DMXs >> tm ) )
+ throw( MCFException( "LoadDMX: error reading number of arcs" ) );
+
+ // allocate memory - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ FRow tU = new FNumber[ tm ]; // arc upper capacities
+ FRow tDfct = new FNumber[ tn ]; // node deficits
+ Index_Set tStartn = new Index[ tm ]; // arc start nodes
+ Index_Set tEndn = new Index[ tm ]; // arc end nodes
+ CRow tC = new CNumber[ tm ]; // arc costs
+ CRow tQ = NULL;
+ if( IsQuad )
+ tQ = new CNumber[ tm ]; // arc quadratic costs
+
+
+ for( Index i = 0 ; i < tn ; ) // all deficits are 0
+ tDfct[ i++ ] = 0; // unless otherwise stated
+
+ // read problem data - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ Index i = 0; // arc counter
+ for(;;) {
+ if( ! ( DMXs >> c ) )
+ break;
+
+ switch( c ) {
+ case( 'c' ): // comment line- - - - - - - - - - - - - - - - - - - - - - -
+ do
+ if( ! DMXs.get( c ) )
+ throw( MCFException( "LoadDMX: error reading the input stream" ) );
+ while( c != '\n' );
+ break;
+
+ case( 'n' ): // description of a node - - - - - - - - - - - - - - - - - -
+ Index j;
+ if( ! ( DMXs >> j ) )
+ throw( MCFException( "LoadDMX: error reading node name" ) );
+
+ if( ( j < 1 ) || ( j > tn ) )
+ throw( MCFException( "LoadDMX: invalid node name" ) );
+
+ FNumber Dfctj;
+ if( ! ( DMXs >> Dfctj ) )
+ throw( MCFException( "LoadDMX: error reading deficit" ) );
+
+ tDfct[ j - 1 ] -= Dfctj;
+ break;
+
+ case( 'a' ): // description of an arc - - - - - - - - - - - - - - - - - -
+ if( i == tm )
+ throw( MCFException( "LoadDMX: too many arc descriptors" ) );
+
+ if( ! ( DMXs >> tStartn[ i ] ) )
+ throw( MCFException( "LoadDMX: error reading start node" ) );
+
+ if( ( tStartn[ i ] < 1 ) || ( tStartn[ i ] > tn ) )
+ throw( MCFException( "LoadDMX: invalid start node" ) );
+
+ if( ! ( DMXs >> tEndn[ i ] ) )
+ throw( MCFException( "LoadDMX: error reading end node" ) );
+
+ if( ( tEndn[ i ] < 1 ) || ( tEndn[ i ] > tn ) )
+ throw( MCFException( "LoadDMX: invalid end node" ) );
+
+ if( tStartn[ i ] == tEndn[ i ] )
+ throw( MCFException( "LoadDMX: self-loops not permitted" ) );
+
+ FNumber LB;
+ if( ! ( DMXs >> LB ) )
+ throw( MCFException( "LoadDMX: error reading lower bound" ) );
+
+ if( ! ( DMXs >> tU[ i ] ) )
+ throw( MCFException( "LoadDMX: error reading upper bound" ) );
+
+ if( ! ( DMXs >> tC[ i ] ) )
+ throw( MCFException( "LoadDMX: error reading arc cost" ) );
+
+ if( tQ ) {
+ if( ! ( DMXs >> tQ[ i ] ) )
+ throw( MCFException( "LoadDMX: error reading arc quadratic cost" ) );
+
+ if( tQ[ i ] < 0 )
+ throw( MCFException( "LoadDMX: negative arc quadratic cost" ) );
+ }
+
+ if( tU[ i ] < LB )
+ throw( MCFException( "LoadDMX: lower bound > upper bound" ) );
+
+ tU[ i ] -= LB;
+ tDfct[ tStartn[ i ] - 1 ] += LB;
+ tDfct[ tEndn[ i ] - 1 ] -= LB;
+ #if( USENAME0 )
+ tStartn[ i ]--; // in the DIMACS format, node names start from 1
+ tEndn[ i ]--;
+ #endif
+ i++;
+ break;
+
+ default: // invalid code- - - - - - - - - - - - - - - - - - - - - - - - -
+ throw( MCFException( "LoadDMX: invalid code" ) );
+
+ } // end( switch( c ) )
+ } // end( for( ever ) )
+
+ if( i < tm )
+ throw( MCFException( "LoadDMX: too few arc descriptors" ) );
+
+ // call LoadNet- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ LoadNet( tn , tm , tn , tm , tU , tC , tDfct , tStartn , tEndn );
+
+ // then pass quadratic costs, if any
+
+ if( tQ )
+ ChgQCoef( tQ );
+
+ // delete the original data structures - - - - - - - - - - - - - - - - - - -
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ delete[] tQ;
+ delete[] tC;
+ delete[] tEndn;
+ delete[] tStartn;
+ delete[] tDfct;
+ delete[] tU;
+
+ } // end( MCFClass::LoadDMX )
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFClass::SetPar( int par , int val )
+{
+ switch( par ) {
+ case( kMaxIter ):
+ MaxIter = val;
+ break;
+
+ case( kReopt ):
+ Senstv = (val == kYes);
+ break;
+
+ default:
+ throw( MCFException( "Error using SetPar: unknown parameter" ) );
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFClass::SetPar( int par, double val )
+{
+ switch( par ) {
+ case( kEpsFlw ):
+ if( EpsFlw != FNumber( val ) ) {
+ EpsFlw = FNumber( val );
+ EpsDfct = EpsFlw * ( nmax ? nmax : 100 );
+ status = kUnSolved;
+ }
+ break;
+
+ case( kEpsDfct ):
+ if( EpsDfct != FNumber( val ) ) {
+ EpsDfct = FNumber( val );
+ status = kUnSolved;
+ }
+ break;
+
+ case( kEpsCst ):
+ if( EpsCst != CNumber( val ) ) {
+ EpsCst = CNumber( val );
+ status = kUnSolved;
+ }
+ break;
+
+ case( kMaxTime ):
+ MaxTime = val;
+ break;
+
+ default:
+ throw( MCFException( "Error using SetPar: unknown parameter" ) );
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFClass::GetPar( int par, int &val )
+{
+ switch( par ) {
+ case( kMaxIter ):
+ val = MaxIter;
+ break;
+
+ case( kReopt ):
+ if( Senstv ) val = kYes;
+ else val = kNo;
+ break;
+
+ default:
+ throw( MCFException( "Error using GetPar: unknown parameter" ) );
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFClass::GetPar( int par , double &val )
+{
+ switch( par ) {
+ case( kEpsFlw ):
+ val = double( EpsFlw );
+ break;
+
+ case( kEpsDfct ):
+ val = double( EpsDfct );
+ break;
+
+ case( kEpsCst ):
+ val = double( EpsCst );
+ break;
+
+ case( kMaxTime ):
+ val = MaxTime;
+ break;
+
+ default:
+ throw( MCFException( "Error using GetPar: unknown parameter" ) );
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFClass::CheckPSol( void )
+{
+ FRow tB = new FNumber[ MCFn() ];
+ MCFDfcts( tB );
+ #if( ! USENAME0 )
+ tB--;
+ #endif
+
+ FRow tX = new FNumber[ MCFm() ];
+ MCFGetX( tX );
+ FONumber CX = 0;
+ for( Index i = 0 ; i < MCFm() ; i++ )
+ if( GTZ( tX[ i ] , EpsFlw ) ) {
+ if( IsClosedArc( i ) )
+ throw( MCFException( "Closed arc with nonzero flow (CheckPSol)" ) );
+
+ if( GT( tX[ i ] , MCFUCap( i ) , EpsFlw ) )
+ throw( MCFException( "Arc flow exceeds capacity (CheckPSol)" ) );
+
+ if( LTZ( tX[ i ] , EpsFlw ) )
+ throw( MCFException( "Arc flow is negative (CheckPSol)" ) );
+
+ CX += ( FONumber( MCFCost( i ) ) +
+ FONumber( MCFQCoef( i ) ) * FONumber( tX[ i ] ) / 2 )
+ * FONumber( tX[ i ] );
+ tB[ MCFSNde( i ) ] += tX[ i ];
+ tB[ MCFENde( i ) ] -= tX[ i ];
+ }
+ delete[] tX;
+ #if( ! USENAME0 )
+ tB++;
+ #endif
+
+ for( Index i = 0 ; i < MCFn() ; i++ )
+ if( ! ETZ( tB[ i ] , EpsDfct ) )
+ throw( MCFException( "Node is not balanced (CheckPSol)" ) );
+
+ delete[] tB;
+ CX -= MCFGetFO();
+ if( ( CX >= 0 ? CX : - CX ) > EpsCst * MCFn() )
+ throw( MCFException( "Objective function value is wrong (CheckPSol)" ) );
+
+ } // end( CheckPSol )
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFClass::CheckDSol( void )
+{
+ CRow tPi = new CNumber[ MCFn() ];
+ MCFGetPi( tPi );
+
+ FONumber BY = 0;
+ for( Index i = 0 ; i < MCFn() ; i++ )
+ BY += FONumber( tPi[ i ] ) * FONumber( MCFDfct( i ) );
+
+ CRow tRC = new CNumber[ MCFm() ];
+ MCFGetRC( tRC );
+
+ FRow tX = new FNumber[ MCFm() ];
+ MCFGetX( tX );
+
+ #if( ! USENAME0 )
+ tPi--;
+ #endif
+
+ FONumber QdrtcCmpnt = 0; // the quadratic part of the objective
+ for( Index i = 0 ; i < MCFm() ; i++ ) {
+ if( IsClosedArc( i ) )
+ continue;
+
+ cFONumber tXi = FONumber( tX[ i ] );
+ cCNumber QCi = CNumber( MCFQCoef( i ) * tXi );
+ QdrtcCmpnt += QCi * tXi;
+ CNumber RCi = MCFCost( i ) + QCi
+ + tPi[ MCFSNde( i ) ] - tPi[ MCFENde( i ) ];
+ RCi -= tRC[ i ];
+
+ if( ! ETZ( RCi , EpsCst ) )
+ throw( MCFException( "Reduced Cost value is wrong (CheckDSol)" ) );
+
+ if( LTZ( tRC[ i ] , EpsCst ) ) {
+ BY += FONumber( tRC[ i ] ) * FONumber( MCFUCap( i ) );
+
+ if( LT( tX[ i ] , MCFUCap( i ) , EpsFlw ) )
+ throw( MCFException( "Csomplementary Slackness violated (CheckDSol)" ) );
+ }
+ else
+ if( GTZ( tRC[ i ] , EpsCst ) && GTZ( tX[ i ] , EpsFlw ) )
+ throw( MCFException( "Complementary Slackness violated (CheckDSol)" ) );
+
+ } // end( for( i ) )
+
+ delete[] tX;
+ delete[] tRC;
+ #if( ! USENAME0 )
+ tPi++;
+ #endif
+ delete[] tPi;
+
+ BY -= ( MCFGetFO() + QdrtcCmpnt / 2 );
+ if( ( BY >= 0 ? BY : - BY ) > EpsCst * MCFn() )
+ throw( MCFException( "Objective function value is wrong (CheckDSol)" ) );
+
+ } // end( CheckDSol )
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFClass::WriteMCF( ostream &oStrm , int frmt )
+{
+ if( ( ! numeric_limits<FNumber>::is_integer ) ||
+ ( ! numeric_limits<CNumber>::is_integer ) )
+ oStrm.precision( 12 );
+
+ switch( frmt ) {
+ case( kDimacs ): // DIMACS standard format - - - - - - - - - - - - - - - -
+ //- - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ case( kQDimacs ): // ... linear or quadratic
+
+ // compute and output preamble and size- - - - - - - - - - - - - - - - - -
+
+ oStrm << "p min " << MCFn() << " ";
+ {
+ Index tm = MCFm();
+ for( Index i = MCFm() ; i-- ; )
+ if( IsClosedArc( i ) || IsDeletedArc( i ) )
+ tm--;
+
+ oStrm << tm << endl;
+ }
+
+ // output arc information- - - - - - - - - - - - - - - - - - - - - - - - -
+
+ for( Index i = 0 ; i < MCFm() ; i++ )
+ if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) {
+ oStrm << "a\t";
+ #if( USENAME0 )
+ oStrm << MCFSNde( i ) + 1 << "\t" << MCFENde( i ) + 1 << "\t";
+ #else
+ oStrm << MCFSNde( i ) << "\t" << MCFENde( i ) << "\t";
+ #endif
+ oStrm << "0\t" << MCFUCap( i ) << "\t" << MCFCost( i );
+
+ if( frmt == kQDimacs )
+ oStrm << "\t" << MCFQCoef( i );
+
+ oStrm << endl;
+ }
+
+ // output node information - - - - - - - - - - - - - - - - - - - - - - - -
+
+ for( Index i = 0 ; i < MCFn() ; ) {
+ cFNumber Dfcti = MCFDfct( i++ );
+ if( Dfcti )
+ oStrm << "n\t" << i << "\t" << - Dfcti << endl;
+ }
+
+ break;
+
+ case( kMPS ): // MPS format(s)- - - - - - - - - - - - - - - - - - - - -
+ case( kFWMPS ): //- - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ // writing problem name- - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ oStrm << "NAME MCF" << endl;
+
+ // writing ROWS section: flow balance constraints- - - - - - - - - - - - -
+
+ oStrm << "ROWS" << endl << " N obj" << endl;
+
+ for( Index i = 0 ; i < MCFn() ; )
+ oStrm << " E c" << i++ << endl;
+
+ // writing COLUMNS section - - - - - - - - - - - - - - - - - - - - - - - -
+
+ oStrm << "COLUMNS" << endl;
+
+ if( frmt == kMPS ) { // "modern" MPS
+ for( Index i = 0 ; i < MCFm() ; i++ )
+ if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) {
+ oStrm << " x" << i << "\tobj\t" << MCFCost( i ) << "\tc";
+ #if( USENAME0 )
+ oStrm << MCFSNde( i ) << "\t-1" << endl << " x" << i << "\tc"
+ << MCFENde( i ) << "\t1" << endl;
+ #else
+ oStrm << MCFSNde( i ) - 1 << "\t-1" << endl << " x" << i << "\tc"
+ << MCFENde( i ) - 1 << "\t1" << endl;
+ #endif
+ }
+ }
+ else // "old" MPS
+ for( Index i = 0 ; i < MCFm() ; i++ ) {
+ ostringstream myField;
+ myField << "x" << i << ends;
+ oStrm << " " << setw( 8 ) << left << myField.str()
+ << " " << setw( 8 ) << left << "obj"
+ << " " << setw( 12 ) << left << MCFCost( i );
+
+ myField.seekp( 0 );
+ myField << "c" << MCFSNde( i ) - ( 1 - USENAME0 ) << ends;
+
+ oStrm << " " << setw( 8 ) << left << myField.str()
+ << " " << setw( 12 ) << left << -1 << endl;
+
+ myField.seekp( 0 );
+ myField << "x" << i << ends;
+
+ oStrm << " " << setw( 8 ) << left << myField.str();
+
+ myField.seekp( 0 );
+ myField << "c" << MCFENde( i ) - ( 1 - USENAME0 ) << ends;
+
+ oStrm << " " << setw( 8 ) << left << myField.str() << endl;
+ }
+
+ // writing RHS section: flow balance constraints - - - - - - - - - - - - -
+
+ oStrm << "RHS" << endl;
+
+ if( frmt == kMPS ) // "modern" MPS
+ for( Index i = 0 ; i < MCFn() ; i++ )
+ oStrm << " rhs\tc" << i << "\t" << MCFDfct( i ) << endl;
+ else // "old" MPS
+ for( Index i = 0 ; i < MCFn() ; i++ ) {
+ ostringstream myField;
+ oStrm << " " << setw( 8 ) << left << "rhs";
+ myField << "c" << i << ends;
+
+ oStrm << " " << setw( 8 ) << left << myField.str()
+ << " " << setw( 12 ) << left << MCFDfct( i ) << endl;
+ }
+
+ // writing BOUNDS section- - - - - - - - - - - - - - - - - - - - - - - - -
+
+ oStrm << "BOUNDS" << endl;
+
+ if( frmt == kMPS ) { // "modern" MPS
+ for( Index i = 0 ; i < MCFm() ; i++ )
+ if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) )
+ oStrm << " UP BOUND\tx" << i << "\t" << MCFUCap( i ) << endl;
+ }
+ else // "old" MPS
+ for( Index i = 0 ; i < MCFm() ; i++ )
+ if( ( ! IsClosedArc( i ) ) && ( ! IsDeletedArc( i ) ) ) {
+ ostringstream myField;
+ oStrm << " " << setw( 2 ) << left << "UP"
+ << " " << setw( 8 ) << left << "BOUND";
+
+ myField << "x" << i << ends;
+
+ oStrm << " " << setw( 8 ) << left << myField.str()
+ << " " << setw( 12 ) << left << MCFUCap( i ) << endl;
+ }
+
+ oStrm << "ENDATA" << endl;
+ break;
+
+ default: // unknown format - - - - - - - - - - - - - - - - - - - -
+ //- - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ oStrm << "Error: unknown format " << frmt << endl;
+ }
+ } // end( MCFClass::WriteMCF )
+
+/*--------------------------------------------------------------------------*/
+
+#if( OPT_USE_NAMESPACES )
+ }; // end( namespace MCFClass_di_unipi_it )
+#endif
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#endif /* MCFClass.h included */
+
+/*--------------------------------------------------------------------------*/
+/*------------------------- End File MCFClass.h ----------------------------*/
+/*--------------------------------------------------------------------------*/
diff --git a/clang/lib/Analysis/MCFSimplex.cpp b/clang/lib/Analysis/MCFSimplex.cpp
new file mode 100644
index 0000000..4b1b75a
--- /dev/null
+++ b/clang/lib/Analysis/MCFSimplex.cpp
@@ -0,0 +1,4197 @@
+/*--------------------------------------------------------------------------*/
+/*---------------------------- File MCFSimplex.C ---------------------------*/
+/*--------------------------------------------------------------------------*/
+/*-- --*/
+/*-- 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 --*/
+/*-- 29 - 08 - 2011 --*/
+/*-- --*/
+/*-- Implementation: --*/
+/*-- --*/
+/*-- Alessandro Bertolini --*/
+/*-- Antonio Frangioni --*/
+/*-- --*/
+/*-- Operations Research Group --*/
+/*-- Dipartimento di Informatica --*/
+/*-- Universita' di Pisa --*/
+/*-- --*/
+/*-- Copyright (C) 2008 - 2011 by Alessandro Bertolini, Antonio Frangioni --*/
+/*-- --*/
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*--------------------------- IMPLEMENTATION -------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+/*--------------------------------------------------------------------------*/
+/*------------------------------ INCLUDES ----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#include "MCFSimplex.h"
+#include <algorithm>
+#include <iostream>
+
+#include <cstdlib>
+#include <ctime>
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------------- USING -----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#if( OPT_USE_NAMESPACES )
+using namespace MCFClass_di_unipi_it;
+#endif
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------------- MACROS ----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#define LIMITATEPRECISION 1
+
+/* If LIMITATEPRECISION is 1, in the quadratic case the Primal Simplex accepts
+ entering arc in base only if the decrease of the o.f. value is bigger than
+ a fixed thresold (EpsOpt * oldFOValue / n). Otherwise, any strict decrease
+ in the o.f. value is accepted. */
+
+#define UNIPI_PRIMAL_INITIAL_SHOW 0
+
+/* If UNIPI_PRIMAL_INITIAL_SHOW = 1, Primal Simplex shows the initial condition
+ (arcs and nodes) of the network. */
+
+#define UNIPI_PRIMAL_ITER_SHOW 0
+
+/* If UNIPI_PRIMAL_FINAL_SHOW = x with x > 0, Primal Simplex shows the condition
+ (arcs and nodes) of the network every x iterations. */
+
+#define UNIPI_PRIMAL_FINAL_SHOW 0
+
+/* If UNIPI_PRIMAL_FINAL_SHOW = 1, Primal Simplex shows the final condition
+ (arcs and nodes) of the network. */
+
+#define UNIPI_DUAL_INITIAL_SHOW 0
+
+/* If UNIPI_DUAL_INITIAL_SHOW = 1, Dual Simplex shows the initial condition
+ (arcs and nodes) of the network. */
+
+#define UNIPI_DUAL_ITER_SHOW 0
+
+/* If UNIPI_DUAL_FINAL_SHOW = x with x > 0, Dual Simplex shows the condition
+ (arcs and nodes) of the network every x iterations. */
+
+#define UNIPI_DUAL_FINAL_SHOW 0
+
+/* If UNIPI_DUAL_FINAL_SHOW = 1, Dual Simplex shows the final condition
+ (arcs and nodes) of the network. */
+
+#define UNIPI_VIS_DUMMY_ARCS 1
+
+/* If UNIPI_VIS_DUMMY_ARCS = 1, Primal Simplex or Dual Simplex shows the
+ conditions of the dummy arcs. */
+
+#define UNIPI_VIS_ARC_UPPER 1
+#define UNIPI_VIS_ARC_COST 1
+#define UNIPI_VIS_ARC_Q_COST 1
+#define UNIPI_VIS_ARC_REDUCT_COST 1
+#define UNIPI_VIS_ARC_STATE 1
+#define UNIPI_VIS_NODE_BASIC_ARC 1
+
+/* When Primal Simplex or Dual Simplex shows the conditions of the network,
+ for every arcs the algorithm shows the flow, for every nodes it shows
+ balance and potential. These 6 flags decide if the algorithm shows a
+ particular value of the arcs/nodes; for example if
+ UNIPI_VIS_ARC_UPPER == 1, the algorithm shows the upper bounds of all
+ arcs. */
+
+#define FOSHOW 0
+
+/* If FOSHOW is 1, the algorithm shows the f.o. value every x iterations
+ (x = UNIPI_PRIMAL_ITER_SHOW or x = UNIPI_DUAL_ITER_SHOW). */
+
+#define OPTQUADRATIC 0
+
+/* If OPTQUADRATIC is 1 the Primal Simplex, in the quadratic case, tries to
+ optimize the update of the potential.
+ Unfortunately this doesn't work well: for this reason it is set to 0. */
+
+
+// QUICK HACK
+#define throw(x) assert(false);
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- FUNCTIONS ------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+template<class T>
+inline T ABS( const T x )
+{
+ return( x >= T( 0 ) ? x : - x );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class T>
+inline void Swap( T &v1 , T &v2 )
+{
+ T temp = v1;
+ v1 = v2;
+ v2 = temp;
+ }
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- CONSTANTS ------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#if( QUADRATICCOST == 0 )
+ static const char DELETED = -2; // ident for deleted arcs
+ static const char CLOSED = -1; // ident for closed arcs
+ static const char BASIC = 0; // ident for basis arcs
+ static const char AT_LOWER = 1; // ident for arcs in L
+ static const char AT_UPPER = 2; // ident for arcs in U
+#endif
+
+/* These macros will be used by method MemAllocCandidateList() to set the
+ values of numCandidateList and hotListSize. There are different macros,
+ according to:
+
+ - the used Simplex
+ - the size of the network
+ - (obviously) the different variables
+
+ This set of values tries to improve the performance of the two algorithms
+ according to diversified situations. */
+
+static const int PRIMAL_LOW_NUM_CANDIDATE_LIST = 30;
+static const int PRIMAL_MEDIUM_NUM_CANDIDATE_LIST = 50;
+static const int PRIMAL_HIGH_NUM_CANDIDATE_LIST = 200;
+static const int PRIMAL_LOW_HOT_LIST_SIZE = 5;
+static const int PRIMAL_MEDIUM_HOT_LIST_SIZE = 10;
+static const int PRIMAL_HIGH_HOT_LIST_SIZE = 20;
+static const int DUAL_LOW_NUM_CANDIDATE_LIST = 6;
+static const int DUAL_HIGH_NUM_CANDIDATE_LIST = 10;
+static const int DUAL_LOW_HOT_LIST_SIZE = 1;
+static const int DUAL_HIGH_HOT_LIST_SIZE = 2;
+
+
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- COSTRUCTOR -----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::MCFSimplex( cIndex nmx , cIndex mmx )
+ :
+ MCFClass( nmx , mmx )
+{
+ #if( QUADRATICCOST )
+ if( numeric_limits<FNumber>::is_integer )
+ throw( MCFException( "FNumber must be float if QUADRATICCOST == 1" ) );
+
+ if( numeric_limits<CNumber>::is_integer )
+ throw( MCFException( "CNumber must be float if QUADRATICCOST == 1" ) );
+ #endif
+
+ usePrimalSimplex = true;
+
+ newSession = true;
+ if( nmax && mmax )
+ MemAlloc();
+ else
+ nmax = mmax = 0;
+
+ #if( QUADRATICCOST )
+ recomputeFOLimits = 100;
+ // recomputeFOLimits represents the limit of the iteration in which
+ // quadratic Primal Simplex computes "manually" the f.o. value
+ EpsOpt = 1e-13;
+ // EpsOpt is the fixed precision of the quadratic Primal Simplex
+ #endif
+
+ pricingRule = kCandidateListPivot;
+ forcedNumCandidateList = 0;
+ forcedHotListSize = 0;
+ usePrimalSimplex = true;
+ nodesP = NULL;
+ nodesD = NULL;
+ arcsP = NULL;
+ arcsD = NULL;
+ candP = NULL;
+ candD = NULL;
+
+ if( numeric_limits<CNumber>::is_integer )
+ MAX_ART_COST = CNumber( 1e7 );
+ else
+ MAX_ART_COST = CNumber( 1e10 );
+
+ } // end( MCFSimplex )
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------- OTHER INITIALIZATIONS -------------------------*/
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::LoadNet( cIndex nmx , cIndex mmx , cIndex pn , cIndex pm ,
+ cFRow pU , cCRow pC , cFRow pDfct ,
+ cIndex_Set pSn , cIndex_Set pEn )
+{
+ MemDeAllocCandidateList();
+
+ if( ( nmx != nmax ) || ( mmx != mmax ) ) {
+ // if the size of the allocated memory changes
+
+ if( nmax && mmax ) { // if the memory was already allocated
+ MemDeAlloc(true); // deallocate the Primal
+ MemDeAlloc(false); // and the Dual data structures
+ nmax = mmax = 0;
+ }
+
+ if( nmx && mmx ) { // if the new dimensions of the memory are correct
+ nmax = nmx;
+ mmax = mmx;
+ MemAlloc();
+ }
+ }
+
+ if( ( ! nmax ) || ( ! mmax ) )
+ // if one of the two dimension of the memory isn't correct
+ nmax = mmax = 0;
+
+ if( nmax ) { // if the new dimensions of the memory are correct
+ n = pn;
+ m = pm;
+
+ if( usePrimalSimplex ) {
+ stopNodesP = nodesP + n;
+ dummyRootP = nodesP + nmax;
+ for( nodePType *node = nodesP ; node != stopNodesP ; node++ )
+ node->balance = pDfct[ node - nodesP ]; // initialize nodes
+
+ stopArcsP = arcsP + m;
+ dummyArcsP = arcsP + mmax;
+ stopDummyP = dummyArcsP + n;
+ for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++ ) {
+ // initialize real arcs
+ if (pC)
+ arc->cost = pC[ arc - arcsP ];
+ else
+ arc->cost = 0;
+ #if( QUADRATICCOST )
+ arc->quadraticCost = 0;
+ #endif
+ if (pU)
+ arc->upper = pU[ arc - arcsP ];
+ else
+ arc->upper = Inf<FNumber>();
+ arc->tail = nodesP + pSn[ arc - arcsP ] - 1;
+ arc->head = nodesP + pEn[ arc - arcsP ] - 1;
+ }
+ }
+ else {
+ stopNodesD = nodesD + n;
+ dummyRootD = nodesD + nmax;
+ for( nodeDType *node = nodesD ; node != stopNodesD ; node++ )
+ node->balance = pDfct[ node - nodesD ]; // initialize nodes
+
+ stopArcsD = arcsD + m;
+ dummyArcsD = arcsD + mmax;
+ stopDummyD = dummyArcsD + n;
+ for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) {
+ // initialize real arcs
+ arc->cost = pC[ arc - arcsD ];
+ #if( QUADRATICCOST )
+ arc->quadraticCost = 0;
+ #endif
+ arc->upper = pU[ arc - arcsD ];
+ arc->tail = nodesD + pSn[ arc - arcsD ] - 1;
+ arc->head = nodesD + pEn[ arc - arcsD ] - 1;
+ }
+
+ CreateAdditionalDualStructures();
+ }
+
+ if( pricingRule == kCandidateListPivot )
+ MemAllocCandidateList();
+
+ status = kUnSolved;
+ }
+ } // end( MCFSimplex::LoadNet )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::SetAlg( bool UsPrml , char WhchPrc )
+{
+ bool newUsePrimalSimplex = UsPrml;
+ bool oldUsePrimalSimplex = usePrimalSimplex;
+ char newPricingRule = WhchPrc;
+ char oldPricingRule = pricingRule;
+
+ usePrimalSimplex = newUsePrimalSimplex;
+ pricingRule = newPricingRule;
+
+ if( ( ! usePrimalSimplex ) && ( pricingRule == kDantzig) ) {
+ pricingRule = kFirstEligibleArc;
+ newPricingRule = pricingRule;
+ }
+
+ if( ( newUsePrimalSimplex != oldUsePrimalSimplex ) ||
+ ( newPricingRule != oldPricingRule ) ) {
+ MemDeAllocCandidateList();
+
+ if( newUsePrimalSimplex != oldUsePrimalSimplex ) {
+ #if( QUADRATICCOST )
+ throw(
+ MCFException( "Primal Simplex is the only option if QUADRATICCOST == 1"
+ ) );
+ }
+ #else
+ MemAlloc();
+ nodePType *nP = nodesP;
+ nodeDType *nD = nodesD;
+ arcPType *aP = arcsP;
+ arcDType *aD = arcsD;
+
+ if( newUsePrimalSimplex ) { // from Dual to Primal
+ if( nodesD == NULL )
+ return;
+
+ if( newSession )
+ CreateInitialDualBase();
+
+ dummyRootP = nodesP + nmax;
+ stopNodesP = nodesP + n;
+ dummyArcsP = arcsP + mmax;
+ stopArcsP = arcsP + m;
+ stopDummyP = dummyArcsP + n;
+ // Copy the old Dual data structure in a new Primal data structure
+ while( nD != stopNodesD ) {
+ nP->prevInT = nodesP + ( nD->prevInT - nodesD );
+ nP->nextInT = nodesP + ( nD->nextInT - nodesD );
+ nP->enteringTArc = arcsP + ( nD->enteringTArc - arcsD );
+ nP->balance = nD->balance;
+ nP->potential = nD->potential;
+ nP->subTreeLevel = nD->subTreeLevel;
+ nP++;
+ nD++;
+ }
+
+ dummyRootP->prevInT = NULL;
+ dummyRootP->nextInT = nodesP + ( dummyRootD->nextInT - nodesD );
+ dummyRootP->enteringTArc = arcsP + ( dummyRootD->enteringTArc - arcsD );
+ dummyRootP->balance = dummyRootD->balance;
+ dummyRootP->potential = dummyRootD->potential;
+ dummyRootP->subTreeLevel = dummyRootD->subTreeLevel;
+ while( aD != stopArcsD ) {
+ aP->tail = nodesP + ( aD->tail - nodesD );
+ aP->head = nodesP + ( aD->head - nodesD );
+ aP->flow = aD->flow;
+ aP->cost = aD->cost;
+ aP->ident = aD->ident;
+ aP->upper = aD->upper;
+ aP++;
+ aD++;
+ }
+
+ aP = dummyArcsP;
+ aD = dummyArcsD;
+ while( aD != stopDummyD ) {
+ aP->tail = nodesP + ( aD->tail - nodesD );
+ aP->head = nodesP + ( aD->head - nodesD );
+ aP->flow = aD->flow;
+ aP->cost = aD->cost;
+ aP->ident = aD->ident;
+ if( ( ETZ(aP->flow, EpsFlw) ) && (aP->ident == AT_UPPER) )
+ aP->ident = AT_LOWER;
+
+ aP->upper = Inf<FNumber>();
+ aP++;
+ aD++;
+ }
+
+ MemDeAlloc(false);
+ if( Senstv && ( status != kUnSolved ) ) {
+ nodePType *node = dummyRootP;
+ for( int i = 0 ; i < n ; i++ )
+ node = node->nextInT;
+
+ node->nextInT = NULL;
+ dummyRootP->prevInT = NULL;
+ dummyRootP->enteringTArc = NULL;
+ // balance the flow
+ CreateInitialPModifiedBalanceVector();
+ PostPVisit( dummyRootP );
+ // restore the primal admissibility
+ BalanceFlow( dummyRootP );
+ ComputePotential( dummyRootP );
+ }
+ else
+ status = kUnSolved;
+ }
+ else { // from Primal to Dual
+ if( nodesP == NULL )
+ return;
+
+ if( newSession )
+ CreateInitialPrimalBase();
+
+ dummyRootD = nodesD + nmax;
+ stopNodesD = nodesD + n;
+ dummyArcsD = arcsD + mmax;
+ stopArcsD = arcsD + m;
+ stopDummyD = dummyArcsD + n;
+ // Copy the old Primal data structure in a new Dual data structure
+ while( nP != stopNodesP ) {
+ nD->prevInT = nodesD + ( nP->prevInT - nodesP );
+ nD->nextInT = nodesD + ( nP->nextInT - nodesP );
+ nD->enteringTArc = arcsD + ( nP->enteringTArc - arcsP );
+ nD->balance = nP->balance;
+ nD->potential = nP->potential;
+ nD->subTreeLevel = nP->subTreeLevel;
+ nP++;
+ nD++;
+ }
+
+ dummyRootD->prevInT = NULL;
+ dummyRootD->nextInT = nodesD + ( dummyRootP->nextInT - nodesP );
+ dummyRootD->enteringTArc = NULL;
+ dummyRootD->balance = dummyRootP->balance;
+ dummyRootD->potential = dummyRootP->potential;
+ dummyRootD->subTreeLevel = dummyRootP->subTreeLevel;
+ while( aP != stopArcsP ) {
+ aD->tail = nodesD + ( aP->tail - nodesP );
+ aD->head = nodesD + ( aP->head - nodesP );
+ aD->flow = aP->flow;
+ aD->cost = aP->cost;
+ aD->ident = aP->ident;
+ aD->upper = aP->upper;
+ aP++;
+ aD++;
+ }
+
+ aP = dummyArcsP;
+ aD = dummyArcsD;
+ while( aP != stopDummyP ) {
+ aD->tail = nodesD + ( aP->tail - nodesP );
+ aD->head = nodesD + ( aP->head - nodesP );
+ aD->flow = aP->flow;
+ aD->cost = aP->cost;
+ aD->ident = aP->ident;
+ aD->upper = 0;
+ aP++;
+ aD++;
+ }
+
+ CreateAdditionalDualStructures();
+ MemDeAlloc(true);
+ nodeDType *node = dummyRootD;
+ for( int i = 0 ; i < n ; i++ )
+ node = node->nextInT;
+
+ node->nextInT = NULL;
+ dummyRootD->enteringTArc = NULL;
+ dummyRootD->prevInT = NULL;
+ if( Senstv && ( status != kUnSolved ) ) {
+ // fix every flow arc according to its reduct cost
+ for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) {
+ if( (arc->ident == AT_LOWER) && LTZ( ReductCost( arc ) , EpsCst ) ) {
+ arc->flow = arc->upper;
+ arc->ident = AT_UPPER;
+ }
+
+ if( ( arc->ident == AT_UPPER ) && GTZ( ReductCost( arc ) , EpsCst ) ) {
+ arc->flow = 0;
+ arc->ident = AT_LOWER;
+ }
+ }
+
+ // balance the flow
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ }
+ else
+ status = kUnSolved;
+ }
+ //#endif
+ }
+#endif
+ if( newPricingRule == kFirstEligibleArc )
+ if( newUsePrimalSimplex )
+ arcToStartP = arcsP;
+ else
+ arcToStartD = arcsD;
+
+ if( ( nmax && mmax ) && ( newPricingRule == kCandidateListPivot ) )
+ MemAllocCandidateList();
+ }
+ } // end( SetAlg )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::SetPar( int par, int val )
+{
+ switch( par ) {
+ case kAlgPrimal:
+ if( val == kYes )
+ SetAlg( true , pricingRule);
+
+ if( val == kNo )
+ SetAlg( false, pricingRule );
+
+ break;
+
+ case kAlgPricing:
+
+ if( ( val == kDantzig ) || ( val == kFirstEligibleArc ) ||
+ ( val == kCandidateListPivot ) )
+ SetAlg( usePrimalSimplex , val );
+
+ break;
+
+ case kNumCandList:
+
+ MemDeAllocCandidateList();
+ forcedNumCandidateList = val;
+ MemAllocCandidateList();
+ forcedNumCandidateList = 0;
+ forcedHotListSize = 0;
+ break;
+
+ case kHotListSize:
+
+ MemDeAllocCandidateList();
+ forcedHotListSize = val;
+ MemAllocCandidateList();
+ forcedNumCandidateList = 0;
+ forcedHotListSize = 0;
+ break;
+
+ case kRecomputeFOLimits:
+
+ recomputeFOLimits = val;
+ break;
+
+ default:
+
+ MCFClass::SetPar(par, val);
+ }
+ } // end( SetPar( int ) )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::SetPar( int par , double val )
+{
+ switch( par ) {
+ case kEpsOpt:
+
+ EpsOpt = val;
+ break;
+
+ default:
+ MCFClass::SetPar( par , val );
+ }
+ } // end( SetPar( double )
+
+/*-------------------------------------------------------------------------*/
+/*--------------- METHODS FOR SOLVING THE PROBLEM -------------------------*/
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::SolveMCF( void )
+{
+ if( MCFt )
+ MCFt->Start();
+
+ //if( status == kUnSolved )
+ if(usePrimalSimplex )
+ CreateInitialPrimalBase();
+ else
+ CreateInitialDualBase();
+
+ newSession = false;
+ if( usePrimalSimplex )
+ PrimalSimplex();
+ else
+ DualSimplex();
+
+ if( MCFt )
+ MCFt->Stop();
+
+ } // end( MCFSimplex::SolveMCF )
+
+/*--------------------------------------------------------------------------*/
+/*---------------------- METHODS FOR READING RESULTS -----------------------*/
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFGetX( FRow F , Index_Set nms , cIndex strt , Index stp )
+{
+ if( stp > m )
+ stp = m;
+
+ if( nms ) {
+ if( usePrimalSimplex )
+ for( Index i = strt ; i < stp ; i++ ) {
+ FNumber tXi = ( arcsP + i )->flow;
+ if( GTZ( tXi , EpsFlw ) ) {
+ *(F++) = tXi;
+ *(nms++) = i;
+ }
+ }
+ else
+ for( Index i = strt ; i < stp ; i++ ) {
+ FNumber tXi = ( arcsD + i )->flow;
+ if( GTZ( tXi , EpsFlw ) ) {
+ *(F++) = tXi;
+ *(nms++) = i;
+ }
+ }
+
+ *nms = Inf<Index>();
+ }
+ else
+ if( usePrimalSimplex )
+ for( Index i = strt; i < stp; i++ )
+ *(F++) = ( arcsP + i )->flow;
+ else
+ for( Index i = strt; i < stp; i++ )
+ *(F++) = ( arcsD + i )->flow;
+
+ } // end( MCFSimplex::MCFGetX( some ) )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFGetRC( CRow CR , cIndex_Set nms , cIndex strt , Index stp )
+{
+ if( nms ) {
+ while( *nms < strt )
+ nms++;
+
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(CR++) = CNumber( ReductCost( arcsP + h ) );
+ else
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(CR++) = ReductCost( arcsD + h );
+ }
+ else {
+ if( stp > m )
+ stp = m;
+
+ if( usePrimalSimplex )
+ for( arcPType* arc = arcsP + strt ; arc < arcsP + stp ; arc++ )
+ *(CR++) = CNumber( ReductCost( arc ) );
+ else
+ for( arcDType* arc = arcsD + strt ; arc < arcsD + stp ; arc++ )
+ *(CR++) = ReductCost( arc );
+ }
+ } // end( MCFSimplex::MCFGetRC( some ) )
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::CNumber MCFSimplex::MCFGetRC( cIndex i )
+{
+ if( usePrimalSimplex )
+ return CNumber( ReductCost( arcsP + i ) );
+ else
+ return( ReductCost( arcsD + i ) );
+
+ } // end( MCFSimplex::MCFGetRC( i ) )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFGetPi( CRow P , cIndex_Set nms , cIndex strt , Index stp )
+{
+ if( stp > n )
+ stp = n;
+
+ if( nms ) {
+ while( *nms < strt )
+ nms++;
+
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(P++) = CNumber( (nodesP + h)->potential );
+ else
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(P++) = (nodesD + h )->potential;
+ }
+ else
+ if( usePrimalSimplex )
+ for( nodePType *node = nodesP + strt ; node < ( nodesP + stp ) ; node++ )
+ *(P++) = CNumber( node->potential );
+ else
+ for( nodeDType *node = nodesD + strt ; node++ < ( nodesD + stp ) ; node++ )
+ *(P++) = node->potential;
+
+ } // end( MCFSimplex::MCFGetPi( some ) )
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::FONumber MCFSimplex::MCFGetFO( void )
+{
+ if( status == kOK )
+ return( (FONumber) GetFO() );
+ else
+ if( status == kUnbounded )
+ return( - Inf<FONumber>() );
+ else
+ return( Inf<FONumber>() );
+
+ } // end( MCFSimplex::MCFGetFO )
+
+/*-------------------------------------------------------------------------*/
+/*----------METHODS FOR READING THE DATA OF THE PROBLEM--------------------*/
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFArcs( Index_Set Startv , Index_Set Endv ,
+ cIndex_Set nms , cIndex strt , Index stp )
+{
+ if( stp > m )
+ stp = m;
+
+ if( nms ) {
+ while( *nms < strt )
+ nms++;
+
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(nms++) ) < stp ; ) {
+ if( Startv )
+ *(Startv++) = Index( (arcsP + h)->tail - nodesP) - USENAME0;
+ if( Endv )
+ *(Endv++) = Index( (arcsP + h)->head - nodesP ) - USENAME0;
+ }
+ else
+ for( Index h ; ( h = *(nms++) ) < stp ; ) {
+ if( Startv )
+ *(Startv++) = Index( (arcsD + h)->tail - nodesD) - USENAME0;
+ if( Endv )
+ *(Endv++) = Index( (arcsD + h)->head - nodesD ) - USENAME0;
+ }
+ }
+ else
+ if( usePrimalSimplex )
+ for( arcPType* arc = arcsP + strt ; arc < (arcsP + stp) ; arc++ ) {
+ if( Startv )
+ *(Startv++) = Index( arc->tail - nodesP ) - USENAME0;
+ if( Endv )
+ *(Endv++) = Index( arc->head - nodesP ) - USENAME0;
+ }
+ else
+ for( arcDType* arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ ) {
+ if( Startv )
+ *(Startv++) = Index( arc->tail - nodesD ) - USENAME0;
+ if( Endv )
+ *(Endv++) = Index( arc->head - nodesD ) - USENAME0;
+ }
+
+ } // end( MCFSimplex::MCFArcs )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFCosts( CRow Costv , cIndex_Set nms ,
+ cIndex strt , Index stp )
+{
+ if( stp > m )
+ stp = m;
+
+ if( nms ) {
+ while( *nms < strt )
+ nms++;
+
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(Costv++) = (arcsP + h)->cost;
+ else
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(Costv++) = (arcsD + h)->cost;
+ }
+ else
+ if( usePrimalSimplex )
+ for( arcPType* arc = arcsP + strt ; arc < (arcsP + stp) ; arc++ )
+ *(Costv++) = arc->cost;
+ else
+ for( arcDType* arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ )
+ *(Costv++) = arc->cost;
+
+ } // end( MCFSimplex::MCFCosts ( some ) )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFQCoef( CRow Qv , cIndex_Set nms ,
+ cIndex strt , Index stp )
+{
+ if( stp > m )
+ stp = m;
+
+ if( nms ) {
+ while( *nms < strt )
+ nms++;
+
+ #if( QUADRATICCOST )
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(Qv++) = (arcsP + h)->quadraticCost;
+ else
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(Qv++) = (arcsD + h)->quadraticCost;
+ #else
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(Qv++) = 0;
+ #endif
+ }
+ else
+ #if( QUADRATICCOST )
+ if( usePrimalSimplex )
+ for( arcPType* arc = arcsP + strt ; arc < ( arcsP + stp ) ; arc++ )
+ *(Qv++) = arc->quadraticCost;
+ else
+ for( arcDType* arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ )
+ *(Qv++) = arc->quadraticCost;
+ #else
+ for( Index h = strt ; h++ < stp ; )
+ *(Qv++) = 0;
+ #endif
+
+ } // end( MCFSimplex::MCFQCoef ( some ) )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFUCaps( FRow UCapv , cIndex_Set nms ,
+ cIndex strt , Index stp )
+{
+ if( stp > m )
+ stp = m;
+
+ if( nms ) {
+ while( *nms < strt )
+ nms++;
+
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(UCapv++) = (arcsP + h)->upper;
+ else
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(UCapv++) = (arcsD + h)->upper;
+ }
+ else
+ if( usePrimalSimplex )
+ for( arcPType* arc = arcsP + strt ; arc < (arcsP + stp ) ; arc++ )
+ *(UCapv++) = arc->upper;
+ else
+ for( arcDType* arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ )
+ *(UCapv++) = arc->upper;
+
+ } // end( MCFSimplex::MCFUCaps ( some ) )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFDfcts( FRow Dfctv , cIndex_Set nms ,
+ cIndex strt , Index stp )
+{
+ if( stp > n )
+ stp = n;
+
+ if( nms ) {
+ while( *nms < strt )
+ nms++;
+
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(Dfctv++) = ( nodesP + h )->balance;
+ else
+ for( Index h ; ( h = *(nms++) ) < stp ; )
+ *(Dfctv++) = (nodesD + h )->balance;
+ }
+ else
+ if( usePrimalSimplex )
+ for( nodePType* node = nodesP + strt ; node < ( nodesP + stp ) ; node++ )
+ *(Dfctv++) = node->balance;
+ else
+ for( nodeDType* node = nodesD + strt ; node < ( nodesD + stp ) ; node++ )
+ *(Dfctv++) = node->balance;
+
+ } // end( MCFSimplex::MCFDfcts )
+
+/*-------------------------------------------------------------------------*/
+/*--------- METHODS FOR ADDING / REMOVING / CHANGING DATA -----------------*/
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgCosts( cCRow NCost , cIndex_Set nms ,
+ cIndex strt , Index stp )
+{
+ if( stp > m )
+ stp = m;
+
+ if( nms ) {
+ while( *nms < strt ) {
+ nms++;
+ NCost++;
+ }
+
+ cIndex_Set tnms = nms; // nms may be needed below
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(tnms++) ) < stp ; )
+ arcsP[ h ].cost = *(NCost++);
+ else
+ for( Index h ; ( h = *(tnms++) ) < stp ; )
+ arcsD[ h ].cost = *(NCost++);
+ }
+ else
+ if( usePrimalSimplex )
+ for( arcPType *arc = arcsP + strt ; arc < (arcsP + stp) ; arc++ )
+ arc->cost = *(NCost++);
+ else
+ for( arcDType *arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ )
+ arc->cost = *(NCost++);
+
+ if( Senstv && ( status != kUnSolved ) )
+ if( usePrimalSimplex )
+ ComputePotential( dummyRootP );
+ else {
+ #if( QUADRATICCOST == 0 )
+ for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ )
+ if( arc->ident > BASIC )
+ if( GTZ( ReductCost( arc ) , EpsCst ) ) {
+ arc->flow = 0;
+ arc->ident = AT_LOWER;
+ }
+ else {
+ arc->flow = arc->upper;
+ arc->ident = AT_UPPER;
+ }
+
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ #endif
+ }
+ else
+ status = kUnSolved;
+
+ } // end( MCFSimplex::ChgCosts )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgCost( Index arc , cCNumber NCost )
+{
+ if( arc >= m )
+ return;
+
+ if( usePrimalSimplex )
+ ( arcsP + arc )->cost = NCost;
+ else
+ ( arcsD + arc )->cost = NCost;
+
+ if( Senstv && ( status != kUnSolved ) ) {
+ if( usePrimalSimplex ) {
+ nodePType *node = ( arcsP + arc )->tail;
+ if( ( ( arcsP + arc )->head)->subTreeLevel < node->subTreeLevel )
+ node = ( arcsP + arc )->head;
+
+ ComputePotential( dummyRootP );
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ nodeDType *node = ( arcsD + arc )->tail;
+ if( ( ( arcsD + arc )->head )->subTreeLevel < node->subTreeLevel )
+ node = ( arcsD + arc )->head;
+
+ ComputePotential( dummyRootD );
+ for( arcDType *a = arcsD ; a != stopArcsD ; a++)
+ if( a->ident > BASIC )
+ if( GTZ( ReductCost( a ) , EpsCst ) ) {
+ a->flow = 0;
+ a->ident = AT_LOWER;
+ }
+ else {
+ a->flow = a->upper;
+ a->ident = AT_UPPER;
+ }
+
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ #endif
+ }
+ }
+ else
+ status = kUnSolved;
+
+ } // end( MCFSimplex::ChgCost )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgQCoef( cCRow NQCoef , cIndex_Set nms ,
+ cIndex strt , Index stp )
+{
+ if( stp > m )
+ stp = m;
+
+ #if( QUADRATICCOST )
+ if( nms ) {
+ while( *nms < strt ) {
+ nms++;
+ NQCoef++;
+ }
+
+ cIndex_Set tnms = nms; // nms may be needed below
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(tnms++) ) < stp ; )
+ arcsP[ h ].quadraticCost = *(NQCoef++);
+ else
+ for( Index h ; ( h = *(tnms++) ) < stp ; )
+ arcsD[ h ].quadraticCost = *(NQCoef++);
+ }
+ else
+ if( usePrimalSimplex )
+ for( arcPType *arc = arcsP + strt ; arc < ( arcsP + stp ) ; arc++ )
+ arc->quadraticCost = *(NQCoef++);
+ else
+ for( arcDType *arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ )
+ arc->quadraticCost = *(NQCoef++);
+
+ if( Senstv && (status != kUnSolved ) )
+ ComputePotential( dummyRootP );
+ else
+ status = kUnSolved;
+ #else
+ if( NQCoef )
+ throw( MCFException( "ChgQCoef: nonzero coefficients not allowed" ) );
+ #endif
+} // end( MCFSimplex::ChgQCoef )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgQCoef( Index arc , cCNumber NQCoef )
+{
+ #if( QUADRATICCOST )
+ if( arc >= m )
+ return;
+
+ if( usePrimalSimplex )
+ ( arcsP + arc )->quadraticCost = NQCoef;
+ else
+ ( arcsD + arc )->quadraticCost = NQCoef;
+
+ if( Senstv && ( status != kUnSolved ) ) {
+ nodePType *node = ( arcsP + arc )->tail;
+ if( ( ( arcsP + arc )->head )->subTreeLevel < node->subTreeLevel )
+ node = ( arcsP + arc )->head;
+
+ ComputePotential( node );
+ }
+ else
+ status = kUnSolved;
+ #else
+ if( NQCoef )
+ throw( MCFException( "ChgQCoef: nonzero coefficients not allowed" ) );
+ #endif
+ } // end( MCFSimplex::ChgQCoef )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgDfcts( cFRow NDfct , cIndex_Set nms ,
+ cIndex strt , Index stp )
+{
+ if( stp > m )
+ stp = m;
+
+ if( nms ) {
+ while( *nms < strt ) {
+ nms++;
+ NDfct++;
+ }
+
+ cIndex_Set tnms = nms; // nms may be needed below
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(tnms++) ) < stp ; )
+ nodesP[ h ].balance = *(NDfct++);
+ else
+ for( Index h ; ( h = *(tnms++) ) < stp ; )
+ nodesD[ h ].balance = *(NDfct++);
+ }
+ else
+ if( usePrimalSimplex )
+ for( nodePType *node = nodesP + strt ; node < ( nodesP + stp ) ; node++ )
+ node->balance = *(NDfct++);
+ else
+ for( nodeDType *node = nodesD + strt ; node < ( nodesD + stp ) ; node++ )
+ node->balance = *(NDfct++);
+
+ if( Senstv && (status != kUnSolved ) )
+ if( usePrimalSimplex ) {
+ CreateInitialPModifiedBalanceVector();
+ PostPVisit( dummyRootP );
+ BalanceFlow( dummyRootP );
+ ComputePotential( dummyRootP );
+ }
+ else {
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ }
+ else
+ status = kUnSolved;
+
+ } // end( MCFSimplex::ChgDfcts )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgDfct( Index nod , cFNumber NDfct )
+{
+ if( nod > n )
+ return;
+
+ if( usePrimalSimplex )
+ ( nodesP + nod - 1 )->balance = NDfct;
+ else
+ ( nodesD + nod - 1 )->balance = NDfct;
+
+ if( Senstv && (status != kUnSolved ) )
+ if( usePrimalSimplex ) {
+ CreateInitialPModifiedBalanceVector();
+ PostPVisit( dummyRootP );
+ BalanceFlow( dummyRootP );
+ ComputePotential( dummyRootP );
+ }
+ else {
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ }
+ else
+ status = kUnSolved;
+
+ } // end( MCFSimplex::ChgDfct )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgUCaps( cFRow NCap , cIndex_Set nms ,
+ cIndex strt , Index stp )
+{
+ if( stp > m )
+ stp = m;
+
+ if( nms ) {
+ while( *nms < strt ) {
+ nms++;
+ NCap++;
+ }
+
+ cIndex_Set tnms = nms; // nms may be needed below
+ if( usePrimalSimplex )
+ for( Index h ; ( h = *(tnms++) ) < stp ; )
+ arcsP[ h ].upper = *(NCap++);
+ else
+ for( Index h ; ( h = *(tnms++) ) < stp ; )
+ arcsD[ h ].upper = *(NCap++);
+ }
+ else
+ if( usePrimalSimplex )
+ for( arcPType *arc = arcsP + strt ; arc < ( arcsP + stp ) ; arc++ )
+ arc->upper = *(NCap++);
+ else
+ for( arcDType *arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ )
+ arc->upper = *(NCap++);
+
+ if( Senstv && (status != kUnSolved ) ) {
+ if( usePrimalSimplex ) {
+ for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++)
+ #if( QUADRATICCOST )
+ if( GT( arc->flow , arc->upper , EpsFlw ) )
+ arc->flow = arc->upper;
+ #else
+ if( GT(arc->flow , arc->upper , EpsFlw ) ||
+ ( ( arc->ident == AT_UPPER ) &&
+ ( ! ETZ( arc->flow - arc->upper , EpsFlw ) ) ) )
+ arc->flow = arc->upper;
+ #endif
+
+ CreateInitialPModifiedBalanceVector();
+ PostPVisit( dummyRootP );
+ BalanceFlow( dummyRootP );
+ ComputePotential( dummyRootP );
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ )
+ if( ( GT( arc->flow , arc->upper , EpsFlw ) && ( arc->ident != BASIC ) ) ||
+ ( ( arc->ident == AT_UPPER ) &&
+ ( ! ETZ( arc->flow - arc->upper , EpsFlw ) ) ) )
+ arc->flow = arc->upper;
+
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ ComputePotential( dummyRootD );
+ #endif
+ }
+ }
+ else
+ status = kUnSolved;
+
+ } // end( MCFSimplex::ChgUCaps )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgUCap( Index arc , cFNumber NCap )
+{
+ if( arc >= m )
+ return;
+
+ if( usePrimalSimplex )
+ ( arcsP + arc )->upper = NCap;
+ else
+ ( arcsD + arc )->upper = NCap;
+
+ if( Senstv && (status != kUnSolved ) ) {
+ if( usePrimalSimplex ) {
+ #if( QUADRATICCOST )
+ if( GT( ( arcsP + arc )->flow , ( arcsP + arc )->upper , EpsFlw ) )
+ ( arcsP + arc )->flow = ( arcsP + arc )->upper;
+ #else
+ if( GT( ( arcsP + arc )->flow , ( arcsP + arc )->upper , EpsFlw ) ||
+ ( ( ( arcsP + arc )->ident == AT_UPPER ) &&
+ ( ! ETZ( ( arcsP + arc )->flow - ( arcsP + arc )->upper , EpsFlw ) ) ) )
+ ( arcsP + arc )->flow = ( arcsP + arc )->upper;
+ #endif
+
+ CreateInitialPModifiedBalanceVector();
+ PostPVisit( dummyRootP );
+ BalanceFlow( dummyRootP );
+ ComputePotential( dummyRootP );
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ if( ( GT( ( arcsD + arc )->flow , ( arcsD + arc )->upper , EpsFlw ) &&
+ ( ( ( arcsD + arc )->ident != BASIC ) ) ) ||
+ ( ( ( arcsD + arc )->ident == AT_UPPER ) &&
+ ( ! ETZ( ( arcsD + arc )->flow - ( arcsD + arc )->upper , EpsFlw ) ) ) ) {
+ ( arcsD + arc )->flow = ( arcsD + arc )->upper;
+ ( arcsD + arc )->ident = AT_UPPER;
+ }
+
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ ComputePotential( dummyRootD );
+ #endif
+ }
+ }
+ else
+ status = kUnSolved;
+
+ } // end( MCFSimplex::ChgUCap )
+
+/*-------------------------------------------------------------------------*/
+
+bool MCFSimplex::IsClosedArc( cIndex name )
+{
+ if( name >= m )
+ return( false );
+
+ #if( QUADRATICCOST )
+ return( ( arcsP + name )->cost == Inf<CNumber>() );
+ #else
+ if( usePrimalSimplex )
+ return( ( ( arcsP + name )->ident < BASIC) );
+ else
+ return( ( ( arcsD + name )->ident < BASIC) );
+ #endif
+ }
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::CloseArc( cIndex name )
+{
+ if( name >= m )
+ return;
+
+ if( usePrimalSimplex ) {
+ arcPType *arc = arcsP + name;
+ #if( QUADRATICCOST )
+ if( arc->cost == Inf<CNumber>() )
+ return;
+
+ arc->cost = Inf<CNumber>();
+ #else
+ if( arc->ident < BASIC )
+ return;
+
+ arc->ident = CLOSED;
+ #endif
+
+ arc->flow = 0;
+
+ if( Senstv && ( status != kUnSolved ) ) {
+ nodePType *node = NULL;
+ if( (arc->tail)->enteringTArc == arc )
+ node = arc->tail;
+
+ if( (arc->head)->enteringTArc == arc )
+ node = arc->head;
+
+ if( node ) {
+ node->enteringTArc = dummyArcsP + ( node - nodesP );
+ nodePType *last = CutAndUpdateSubtree( node , -node->subTreeLevel + 1 );
+ PasteSubtree( node , last , dummyRootP );
+ node->enteringTArc = dummyArcsP + ( node - nodesP );
+ }
+
+ CreateInitialPModifiedBalanceVector();
+ PostPVisit(dummyRootP);
+ BalanceFlow(dummyRootP);
+ ComputePotential(dummyRootP);
+ }
+ else
+ status = kUnSolved;
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ arcDType *arc = arcsD + name;
+ if( arc->ident < BASIC )
+ return;
+
+ arc->flow = 0;
+ arc->ident = CLOSED;
+
+ if( Senstv && ( status != kUnSolved ) ) {
+ nodeDType *node = NULL;
+ if( ( arc->tail )->enteringTArc == arc)
+ node = arc->tail;
+
+ if( ( arc->head )->enteringTArc == arc )
+ node = arc->head;
+
+ if( node ) {
+ node->enteringTArc = dummyArcsD + ( node - nodesD );
+ nodeDType *last = CutAndUpdateSubtree( node , -node->subTreeLevel + 1 );
+ PasteSubtree( node , last , dummyRootD );
+ node->enteringTArc = dummyArcsD + ( node - nodesD );
+ ComputePotential( dummyRootD );
+
+ for( arcDType *a = arcsD ; a != stopArcsD ; a++ )
+ if( a->ident > BASIC )
+ if( GTZ( ReductCost( a ) , EpsCst ) ) {
+ a->flow = 0;
+ a->ident = AT_LOWER;
+ }
+ else {
+ a->flow = a->upper;
+ a->ident = AT_UPPER;
+ }
+ }
+
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit(dummyRootD);
+ ComputePotential(dummyRootD);
+ }
+ else
+ status = kUnSolved;
+ #endif
+ }
+ } // end( MCFSimplex::CloseArc )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::DelNode( cIndex name )
+{
+ if( name >= n )
+ return;
+
+ if( usePrimalSimplex ) {
+ nodePType *node = nodesP + name;
+ nodePType *last = CutAndUpdateSubtree(node, -node->subTreeLevel);
+ nodePType *n = node->nextInT;
+ while( n ) {
+ if( n->subTreeLevel == 1 )
+ n->enteringTArc = dummyArcsP + ( n - nodesP );
+
+ n = n->nextInT;
+ }
+
+ PasteSubtree( node , last , dummyRootP );
+ n = node->nextInT;
+ dummyRootP->nextInT = n;
+ n->prevInT = dummyRootP;
+
+ for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++ ) {
+ if( ( ( arc->tail ) == node) || ( ( arc->head ) == node ) ) {
+ arc->flow = 0;
+ #if( QUADRATICCOST )
+ arc->cost = Inf<CNumber>();
+ #else
+ arc->ident = CLOSED;
+ #endif
+ }
+ }
+
+ for( arcPType *arc = dummyArcsP ; arc != stopDummyP ; arc++ ) {
+ if( ( ( arc->tail ) == node ) || ( ( arc->head ) == node ) ) {
+ arc->flow = 0;
+ #if( QUADRATICCOST )
+ arc->cost = Inf<CNumber>();
+ #else
+ arc->ident = CLOSED;
+ #endif
+ }
+ }
+
+ CreateInitialPModifiedBalanceVector();
+ PostPVisit( dummyRootP );
+ BalanceFlow( dummyRootP );
+ ComputePotential( dummyRootP );
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ nodeDType *node = nodesD + name;
+ nodeDType *last = CutAndUpdateSubtree( node , -node->subTreeLevel );
+ nodeDType *n = node->nextInT;
+ while( n ) {
+ if( n->subTreeLevel == 1 )
+ n->enteringTArc = dummyArcsD + ( n - nodesD );
+
+ n = n->nextInT;
+ }
+
+ PasteSubtree( node , last , dummyRootD );
+ n = node->nextInT;
+ dummyRootD->nextInT = n;
+ n->prevInT = dummyRootD;
+
+ for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ )
+ if( ( ( arc->tail ) == node) || ( ( arc->head ) == node ) ) {
+ arc->flow = 0;
+ arc->ident = CLOSED;
+ }
+
+ for( arcDType *arc = dummyArcsD ; arc != stopDummyD ; arc++ )
+ if( ( ( arc->tail ) == node ) || ( ( arc->head ) == node ) ) {
+ arc->flow = 0;
+ arc->ident = CLOSED;
+ }
+
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ ComputePotential( dummyRootD );
+ #endif
+ }
+ } // end( MCFSimplex::DelNode )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::OpenArc( cIndex name )
+{
+ if( name >= m )
+ return;
+
+ if( usePrimalSimplex ) {
+ /* Quadratic case is not implemented for a theory bug.
+ Infact a closed arc in the quadratic case has its cost fixed to infinity, and
+ it's impossible to restore the old value. */
+
+ #if( QUADRATICCOST == 0 )
+ arcPType *arc = arcsP + name;
+ if( arc->ident == CLOSED ) {
+ arc->ident = AT_LOWER;
+ arc->flow = 0;
+ }
+ #endif
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ arcDType *arc = arcsD + name;
+ if( arc->ident == CLOSED ) {
+ if( GTZ( ReductCost( arc ) , EpsCst ) )
+ arc->ident = AT_LOWER;
+ else {
+ arc->ident = AT_UPPER;
+ arc->flow = arc->upper;
+ if( Senstv && ( status != kUnSolved ) ) {
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ }
+ else
+ status = kUnSolved;
+ }
+ }
+ #endif
+ }
+ } // end( MCFSimplex:OpenArc )
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::Index MCFSimplex::AddNode( cFNumber aDfct )
+{
+ if( n >= nmax )
+ return( Inf<Index>() );
+
+ n++;
+ if( usePrimalSimplex ) {
+ nodePType *newNode = nodesP + n - 1;
+ stopArcsP->tail = newNode;
+ stopArcsP->head = dummyRootP;
+ stopArcsP->upper = Inf<FNumber>();
+ stopArcsP->flow = 0;
+ stopArcsP->cost = Inf<CNumber>();
+ #if( QUADRATICCOST )
+ stopArcsP->quadraticCost = 0;
+ #else
+ stopArcsP->ident = BASIC;
+ #endif
+ stopArcsP++;
+ newNode->balance = aDfct;
+ newNode->prevInT = dummyRootP;
+ newNode->nextInT = dummyRootP->nextInT;
+ (dummyRootP->nextInT)->prevInT = newNode;
+ dummyRootP->nextInT = newNode;
+ newNode->enteringTArc = stopArcsP--;
+ newNode->potential = 0;
+ #if(QUADRATICCOST)
+ newNode->sumQuadratic = 0;
+ #endif
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ nodeDType *newNode = nodesD + n - 1;
+ stopArcsD->tail = newNode;
+ stopArcsD->head = dummyRootD;
+ stopArcsD->upper = 0;
+ stopArcsD->flow = 0;
+ stopArcsD->cost = Inf<CNumber>();
+ stopArcsD->ident = BASIC;
+ newNode->balance = aDfct;
+ newNode->prevInT = dummyRootD;
+ newNode->nextInT = dummyRootD->nextInT;
+ (dummyRootD->nextInT)->prevInT = newNode;
+ dummyRootD->nextInT = newNode;
+ newNode->enteringTArc = stopArcsD;
+ newNode->potential = 0;
+ newNode->firstFs = stopArcsD;
+ newNode->firstBs = NULL;
+ stopArcsD->nextFs = NULL;
+ stopArcsD->nextBs = dummyRootD->firstBs;
+ dummyRootD->firstBs = stopArcsD;
+ stopArcsD++;
+ #endif
+ }
+
+ return( n );
+
+ } // end( MCFSimplex::AddNode )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::ChangeArc( cIndex name , cIndex nSN , cIndex nEN )
+{
+ if( name >= m )
+ return;
+
+ CloseArc( name );
+
+ if( usePrimalSimplex ) {
+ if( nSN <= n )
+ (arcsP + name)->tail = (nodesP + nSN + USENAME0 - 1);
+ if( nEN <= n )
+ (arcsP + name)->head = (nodesP + nEN + USENAME0 - 1);
+ }
+ else {
+ if( nSN <= n )
+ (arcsD + name)->tail = (nodesD + nSN + USENAME0 - 1);
+ if( nEN <= n )
+ (arcsD + name)->head = (nodesD + nEN + USENAME0 - 1);
+ }
+
+ OpenArc( name );
+
+ } // end( MCFSimplex::ChangeArc )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::DelArc( cIndex name )
+{
+ if( name >= m )
+ return;
+
+ if( usePrimalSimplex ) {
+ arcPType *arc = arcsP + name;
+ #if( QUADRATICCOST )
+ if( arc->upper == -Inf<FNumber>() )
+ return;
+
+ if( arc->cost < Inf<CNumber>() )
+ #else
+ if( arc->cost == DELETED )
+ return;
+
+ if( arc->cost >= BASIC )
+ #endif
+ CloseArc( name );
+
+ #if( QUADRATICCOST )
+ arc->upper = -Inf<FNumber>();
+ #else
+ arc->ident = DELETED;
+ #endif
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ arcDType *arc = arcsD + name;
+ if( arc->cost == DELETED )
+ return;
+
+ if( arc->cost >= BASIC )
+ CloseArc( name );
+
+ arc->ident = DELETED;
+ #endif
+ }
+ } // end( MCFSimplex::DelArc )
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::Index MCFSimplex::AddArc( cIndex Start , cIndex End ,
+ cFNumber aU , cCNumber aC )
+{
+ if( usePrimalSimplex ) {
+ arcPType *arc = arcsP;
+ #if( QUADRATICCOST )
+ while( ( arc->upper != -Inf<FNumber>() ) && ( arc != stopArcsP ) )
+ #else
+ while( ( arc->ident > DELETED ) && ( arc != stopArcsP ) )
+ #endif
+ arc++;
+
+ if( arc == stopArcsP ) {
+ if( m >= mmax )
+ return( Inf<Index>() );
+
+ m++;
+ stopArcsP++;
+ }
+
+ Index pos = ( arc - arcsP ) + 1;
+ arc->tail = nodesP + Start + USENAME0 - 1;
+ arc->head = nodesP + End + USENAME0 - 1;
+ arc->upper = aU;
+ arc->cost = aC;
+ arc->flow = 0;
+ #if( QUADRATICCOST )
+ arc->quadraticCost = 0;
+ #else
+ arc->ident = AT_LOWER;
+ #endif
+ ComputePotential( dummyRootP );
+ return( pos );
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ arcDType *arc = arcsD;
+ while( ( arc->ident > DELETED ) && ( arc != stopArcsD ) )
+ arc++;
+
+ if( arc == stopArcsD ) {
+ if( m >= mmax )
+ return( Inf<Index>() );
+
+ m++;
+ stopArcsD++;
+ }
+
+ Index pos = ( arc - arcsD ) + 1;
+ arc->tail = nodesD + Start + USENAME0 - 1;
+ arc->head = nodesD + End + USENAME0 - 1;
+ arc->upper = aU;
+ arc->cost = aC;
+ if( GEZ( ReductCost( arc ) , EpsCst ) ) {
+ arc->flow = 0;
+ arc->ident = AT_LOWER;
+ if( Senstv && ( status != kUnSolved ) ) {
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ ComputePotential( dummyRootD );
+ }
+ else
+ status = kUnSolved;
+ }
+ else {
+ arc->flow = arc->upper;
+ arc->ident = AT_UPPER;
+ if( Senstv && ( status != kUnSolved ) ) {
+ CreateInitialDModifiedBalanceVector();
+ PostDVisit( dummyRootD );
+ ComputePotential( dummyRootD );
+ }
+ else
+ status = kUnSolved;
+ }
+
+ ComputePotential( dummyRootD );
+ return( pos );
+ #endif
+ }
+ } // end( MCFSimplex::AddArc )
+
+/*--------------------------------------------------------------------------*/
+
+bool MCFSimplex::IsDeletedArc( cIndex name )
+{
+ if( name >= m )
+ return( false );
+
+ #if( QUADRATICCOST )
+ return( ( ( arcsP + name )->upper == -Inf<FNumber>() ) );
+ #else
+ if( usePrimalSimplex )
+ return( ( arcsP + name )->ident == DELETED );
+ else
+ return( ( arcsD + name )->ident == DELETED );
+ #endif
+ }
+
+/*--------------------------------------------------------------------------*/
+/*------------------------------ DESTRUCTOR --------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::~MCFSimplex()
+{
+ MemDeAllocCandidateList();
+ MemDeAlloc(true);
+ MemDeAlloc(false);
+ }
+
+/*--------------------------------------------------------------------------*/
+/*---------------------------- PRIVATE METHODS -----------------------------*/
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MemAlloc( void )
+{
+ if( usePrimalSimplex ) {
+ nodesP = new nodePType[ nmax + 1 ]; // array of nodes
+ arcsP = new arcPType[ mmax + nmax ]; // array of arcs
+ dummyArcsP = arcsP + mmax; // artificial arcs are in the last
+ // nmax positions of the array arcs[]
+ }
+ else {
+ nodesD = new nodeDType[ nmax + 1 ]; // array of nodes
+ arcsD = new arcDType[ mmax + nmax ]; // array of arcs
+ dummyArcsD = arcsD + mmax; // artificial arcs are in the last nmax
+ // positions of the array arcs[]
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MemDeAlloc( bool whatDeAlloc )
+{
+ if( whatDeAlloc ) {
+ delete[] nodesP;
+ delete[] arcsP;
+ nodesP = NULL;
+ arcsP = NULL;
+ }
+ else {
+ delete[] nodesD;
+ delete[] arcsD;
+ nodesD = NULL;
+ arcsD = NULL;
+ }
+ MemDeAllocCandidateList( );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MemAllocCandidateList( void )
+{
+ if( usePrimalSimplex ) {
+ if( m < 10000 ) {
+ numCandidateList = PRIMAL_LOW_NUM_CANDIDATE_LIST;
+ hotListSize = PRIMAL_LOW_HOT_LIST_SIZE;
+ }
+ else
+ if( m > 100000 ) {
+ numCandidateList = PRIMAL_HIGH_NUM_CANDIDATE_LIST;
+ hotListSize = PRIMAL_HIGH_HOT_LIST_SIZE ;
+ }
+ else {
+ numCandidateList = PRIMAL_MEDIUM_NUM_CANDIDATE_LIST;
+ hotListSize = PRIMAL_MEDIUM_HOT_LIST_SIZE;
+ }
+
+ #if( QUADRATICCOST )
+ int coef = 1;
+ // If the number of the arcs is more than 10000, numCandidateList and hotListSize
+ // are increased to improve the performance of the Quadratic Primal Simplex
+ if( m > 10000 )
+ coef = 10;
+
+ numCandidateList = numCandidateList * coef;
+ hotListSize = hotListSize * coef;
+ #endif
+
+ if( forcedNumCandidateList > 0 )
+ numCandidateList = forcedNumCandidateList;
+
+ if( forcedHotListSize > 0 )
+ hotListSize = forcedHotListSize;
+
+ candP = new primalCandidType[ hotListSize + numCandidateList + 1 ];
+ }
+ else {
+ if( n < 10000 ) {
+ numCandidateList = DUAL_LOW_NUM_CANDIDATE_LIST;
+ hotListSize = DUAL_LOW_HOT_LIST_SIZE;
+ }
+ else {
+ numCandidateList = DUAL_HIGH_NUM_CANDIDATE_LIST;
+ hotListSize = DUAL_HIGH_HOT_LIST_SIZE;
+ }
+
+ if( forcedNumCandidateList > 0 )
+ numCandidateList = forcedNumCandidateList;
+
+ if( forcedHotListSize > 0 )
+ hotListSize = forcedHotListSize;
+
+ candD = new dualCandidType[ hotListSize + numCandidateList + 1 ];
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MemDeAllocCandidateList( void )
+{
+ delete[] candP;
+ candP = NULL;
+ delete[] candD;
+ candD = NULL;
+}
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateInitialPrimalBase( void )
+{
+ arcPType *arc;
+ nodePType *node;
+ for( arc = arcsP ; arc != stopArcsP ; arc++ ) { // initialize real arcs
+ arc->flow = 0;
+ #if( QUADRATICCOST == 0 )
+ arc->ident = AT_LOWER;
+ #endif
+ }
+
+ for( arc = dummyArcsP ; arc != stopDummyP ; arc++ ) { // initialize dummy arcs
+ node = nodesP + ( arc - dummyArcsP );
+ if( node->balance > 0 ) { // sink nodes
+ arc->tail = dummyRootP;
+ arc->head = node;
+ arc->flow = node->balance;
+ }
+ else { // source nodes or transit node
+ arc->tail = node;
+ arc->head = dummyRootP;
+ arc->flow = -node->balance;
+ }
+
+ arc->cost = MAX_ART_COST;
+ #if( QUADRATICCOST )
+ arc->quadraticCost = 0;
+ #else
+ arc->ident = BASIC;
+ #endif
+ arc->upper = Inf<FNumber>();
+ }
+
+ dummyRootP->balance = 0;
+ dummyRootP->prevInT = NULL;
+ dummyRootP->nextInT = nodesP;
+ dummyRootP->enteringTArc = NULL;
+ #if( QUADRATICCOST )
+ dummyRootP->sumQuadratic = 0;
+ #endif
+ dummyRootP->potential = MAX_ART_COST;
+ dummyRootP->subTreeLevel = 0;
+ for( node = nodesP ; node != stopNodesP ; node++) { // initialize nodes
+ node->prevInT = node - 1;
+ node->nextInT = node + 1;
+ node->enteringTArc = dummyArcsP + (node - nodesP);
+ #if( QUADRATICCOST )
+ node->sumQuadratic = (node->enteringTArc)->quadraticCost;
+ #endif
+ if( node->balance > 0 ) // sink nodes
+ node->potential = 2 * MAX_ART_COST;
+ else // source nodes or transit node
+ node->potential = 0;
+
+ node->subTreeLevel = 1;
+ }
+
+ nodesP->prevInT = dummyRootP;
+ ( nodesP + n - 1 )->nextInT = NULL;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateInitialDualBase( void )
+{
+ arcDType *arc;
+ nodeDType *node;
+ for( arc = dummyArcsD ; arc != stopDummyD ; arc++ ) { // initialize dummy arcs
+ node = nodesD + ( arc - dummyArcsD );
+ arc->tail = node;
+ arc->head = dummyRootD;
+ arc->flow = -node->balance;
+ arc->cost = MAX_ART_COST;
+ #if( QUADRATICCOST )
+ arc->quadraticCost = 0;
+ #else
+ arc->ident = BASIC;
+ #endif
+ arc->upper = 0;
+ }
+
+ for( arc = arcsD ; arc != stopArcsD ; arc++ ) { // initialize real arcs
+ if( GTZ( arc->cost , EpsCst ) ) {
+ arc->flow = 0;
+ #if( QUADRATICCOST == 0 )
+ arc->ident = AT_LOWER;
+ #endif
+ }
+ else {
+ #if( QUADRATICCOST == 0 )
+ arc->ident = AT_UPPER;
+ #endif
+ arc->flow = arc->upper;
+ ( dummyArcsD + ( ( arc->tail ) - nodesD ) )->flow =
+ ( dummyArcsD + ( ( arc->tail ) - nodesD ) )->flow - arc->upper;
+
+ ( dummyArcsD + ( ( arc->head ) - nodesD ) )->flow =
+ ( dummyArcsD + ( ( arc->head ) - nodesD ) )->flow + arc->upper;
+ }
+ }
+
+ dummyRootD->balance = 0;
+ dummyRootD->prevInT = NULL;
+ dummyRootD->nextInT = nodesD;
+ dummyRootD->enteringTArc = NULL;
+ #if( QUADRATICCOST )
+ dummyRootD->sumQuadratic = 0;
+ #endif
+ dummyRootD->potential = MAX_ART_COST;
+ dummyRootD->subTreeLevel = 0;
+ for( node = nodesD ; node != stopNodesD ; node++ ) { // initialize nodes
+ node->prevInT = node - 1;
+ node->nextInT = node + 1;
+ node->enteringTArc = dummyArcsD + ( node - nodesD );
+ #if( QUADRATICCOST )
+ node->sumQuadratic = ( node->enteringTArc )->quadraticCost;
+ #endif
+ node->potential = 0;
+ node->subTreeLevel = 1;
+ node->whenInT2 = 0;
+ }
+
+ nodesD->prevInT = dummyRootD;
+ ( nodesD + n - 1 )->nextInT = NULL;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateAdditionalDualStructures( void )
+{
+ // this method creates, in a Dual context, the Backward Star and the
+ // Forward Star of every node
+
+ for( nodeDType *node = nodesD ; node != stopNodesD ; node++) {
+ // initialize nodes
+ node->firstBs = NULL;
+ node->firstFs = NULL;
+ node->numArcs = 0;
+ }
+
+ dummyRootD->firstBs = NULL;
+ dummyRootD->firstFs = NULL;
+ dummyRootD->numArcs = 0;
+ for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) {
+ // initialize real arcs
+ arc->nextFs = ( arc->tail )->firstFs;
+ ( arc->tail )->firstFs = arc;
+ arc->nextBs = ( arc->head )->firstBs;
+ ( arc->head )->firstBs = arc;
+ ( arc->tail )->numArcs++;
+ ( arc->head )->numArcs++;
+ }
+
+ ResetWhenInT2();
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PrimalSimplex( void )
+{
+ #if( UNIPI_PRIMAL_INITIAL_SHOW == 0 )
+ #if( UNIPI_PRIMAL_ITER_SHOW == 0 )
+ #if( UNIPI_PRIMAL_FINAL_SHOW == 0 )
+ //cout << endl;
+ #endif
+ #endif
+ #endif
+ #if( UNIPI_PRIMAL_INITIAL_SHOW )
+ cout << endl;
+ for( int t = 0; t < 3; t++ )
+ cout << "\t";
+ cout << "PRIMALE MCFSimplex: ARCHI E NODI ALL' INIZIO" << endl;
+ ShowSituation( 3 );
+ #endif
+ #if( QUADRATICCOST )
+ #if( LIMITATEPRECISION )
+ foValue = GetFO();
+ int cont = 0;
+ #endif
+ #endif
+
+ status = kUnSolved;
+ if( pricingRule != kCandidateListPivot )
+ arcToStartP = arcsP;
+
+ iterator = 0; // setting the initial arc for the Dantzig or First Elibigle Rule
+
+ arcPType *enteringArc;
+ arcPType *leavingArc;
+ if( pricingRule == kCandidateListPivot )
+ InitializePrimalCandidateList();
+
+ while( status == kUnSolved ) {
+ iterator++;
+ switch( pricingRule ) {
+ case( kDantzig ):
+ enteringArc = RuleDantzig();
+ break;
+ case( kFirstEligibleArc ):
+ enteringArc = PRuleFirstEligibleArc();
+ break;
+ case( kCandidateListPivot ):
+ enteringArc = RulePrimalCandidateListPivot();
+ break;
+ }
+
+ #if( QUADRATICCOST )
+ #if( LIMITATEPRECISION )
+ /* In the quadratic case with LIMITATEPRECISION == 1, the entering arcs
+ are selected according to a thresold.
+ This thresold is definited according to the old f.o. value.
+ If Primal Simplex doesn't find an entering arc, it calculates again
+ the f.o. value, and try again. */
+ if( enteringArc == NULL ) {
+ foValue = GetFO();
+ switch( pricingRule ) {
+ case( kDantzig ):
+ enteringArc = RuleDantzig();
+ break;
+ case( kFirstEligibleArc ):
+ enteringArc = PRuleFirstEligibleArc();
+ break;
+ case( kCandidateListPivot ):
+ enteringArc = RulePrimalCandidateListPivot();
+ break;
+ }
+ }
+ #endif
+ #endif
+
+ if( pricingRule != kCandidateListPivot ) {
+ // in every iteration the algorithm changes the initial arc for
+ // Dantzig and First Eligible Rule.
+ arcToStartP++;
+ if( arcToStartP == stopArcsP )
+ arcToStartP = arcsP;
+ }
+
+ if( enteringArc ) {
+ arcPType *arc;
+ nodePType *k1;
+ nodePType *k2;
+ /* If the reduced cost of the entering arc is > 0, the Primal Simplex
+ pushes flow in the cycle determinated by T and entering arc for decreases
+ flow in the entering arc: in the linear case entering arc's flow goes to 0,
+ in the quadratic case it decreases while it's possibile.
+ If the reduced cost of the entering arc is < 0, the Primal Simplex pushes
+ flow in the cycle determinated by T and entering arc for increases flow
+ in the entering arc: in the linear case entering arc's flow goes to upper
+ bound, in the quadratic case it increases while it's possibile. */
+
+ #if( QUADRATICCOST )
+ FONumber t;
+ FONumber theta;
+ FONumber deltaFO;
+ FNumber theta2;
+ CNumber Q = ( enteringArc->tail )->sumQuadratic +
+ ( enteringArc->head )->sumQuadratic + enteringArc->quadraticCost;
+ // Q is the sum of the quadratic coefficient in the cycle determinated by T
+ // and entering arc.
+ FONumber rc = ReductCost( enteringArc );
+ if( ETZ( Q, EpsCst ) )
+ theta = Inf<FNumber>(); // This value will be certainly decreased
+ else
+ theta = ABS( rc / Q );
+ // This is the best theta value (with best f.o. value decrease)
+
+ leavingArc = enteringArc;
+ nodePType *cycleRoot; // The radix of the cycle determinated by T and entering arc.
+ if( GTZ( rc , EpsCst ) ) {
+ #else
+ FNumber t;
+ FNumber theta;
+ if( enteringArc->ident == AT_UPPER ) {
+ #endif
+ /* Primal Simplex increases or decreases entering arc's flow.
+ "theta" is a positive value.
+ For this reason the algorithm uses two nodes ("k1" and "k2") to push
+ flow ("theta") from "k1" to "k2".
+ According to entering arc's reduct cost, the algorithm determinates
+ "k1" and "k2" */
+
+ k1 = enteringArc->head;
+ k2 = enteringArc->tail;
+ #if( QUADRATICCOST )
+ theta = min( theta , enteringArc->flow );
+ // The best value for theta is compared with the entering arc's bound
+ theta2 = - theta;
+ #else
+ theta = enteringArc->flow;
+ #endif
+ }
+ else {
+ k1 = enteringArc->tail;
+ k2 = enteringArc->head;
+ #if( QUADRATICCOST )
+ theta = min( theta , enteringArc->upper - enteringArc->flow );
+ // The best value for theta is compared with the entering arc's bound
+ theta2 = theta;
+ #else
+ theta = enteringArc->upper - enteringArc->flow;
+ #endif
+ }
+
+ nodePType *memK1 = k1;
+ nodePType *memK2 = k2;
+ leavingArc = NULL;
+ #if( QUADRATICCOST )
+ #if( LIMITATEPRECISION )
+ deltaFO = rc * theta2 + Q * theta2 * theta2 / 2;
+ #endif
+ bool leavingReducesFlow = GTZ( rc , EpsCst );
+ #else
+ bool leavingReducesFlow = GTZ( ReductCost( enteringArc ) , EpsCst );
+ #endif
+ // Compute "theta", find outgoing arc and "root" of the cycle
+ bool leave;
+ // Actual "theta" is compared with the bounds of the other cycle's arcs
+ while( k1 != k2 ) {
+ if( k1->subTreeLevel > k2->subTreeLevel ) {
+ arc = k1->enteringTArc;
+ if( arc->tail != k1 ) {
+ t = arc->upper - arc->flow;
+ leave = false;
+ }
+ else {
+ t = arc->flow;
+ leave = true;
+ }
+
+ if( t < theta ) {
+ // The algorithm has found a possible leaving arc
+ theta = t;
+ leavingArc = arc;
+ leavingReducesFlow = leave;
+ // If "leavingReducesFlow" == true, if this arc is selected to exit the base,
+ // it decreases its flow
+ }
+
+ k1 = Father( k1 , arc );
+ }
+ else {
+ arc = k2->enteringTArc;
+ if( arc->tail == k2 ) {
+ t = arc->upper - arc->flow;
+ leave = false;
+ }
+ else {
+ t = arc->flow;
+ leave = true;
+ }
+
+ if( t <= theta ) {
+ // The algorithm has found a possible leaving arc
+ theta = t;
+ leavingArc = arc;
+ leavingReducesFlow = leave;
+ // If "leavingReducesFlow" == true, if this arc is selected to exit the base,
+ // it decreases its flow
+ }
+
+ k2 = Father(k2, arc);
+ }
+ }
+
+ #if( QUADRATICCOST )
+ cycleRoot = k1;
+ #endif
+
+ if( leavingArc == NULL )
+ leavingArc = enteringArc;
+
+ if( theta >= Inf<FNumber>() ) {
+ status = kUnbounded;
+ break;
+ }
+
+ // Update flow with "theta"
+ k1 = memK1;
+ k2 = memK2;
+ #if( QUADRATICCOST )
+ if( enteringArc->tail == k1 )
+ theta2 = theta;
+ else
+ theta2 = -theta;
+
+ // "theta" is a positive value in every case.
+ // "theta2" is the real theta value according to the real
+ // direction of the entering arc
+ #if( LIMITATEPRECISION )
+ deltaFO = rc * theta2 + Q * theta2 * theta2 / 2;
+ // The decrease of the f.o. value in the quadratic case
+ #endif
+ #endif
+
+ if( ! ETZ(theta , EpsFlw ) ) {
+ if( enteringArc->tail == k1 )
+ enteringArc->flow = enteringArc->flow + theta;
+ else
+ enteringArc->flow = enteringArc->flow - theta;
+
+ while( k1 != k2 ) {
+ if( k1->subTreeLevel > k2->subTreeLevel ) {
+ arc = k1->enteringTArc;
+ if( arc->tail != k1 )
+ arc->flow = arc->flow + theta;
+ else
+ arc->flow = arc->flow - theta;
+
+ k1 = Father(k1, k1->enteringTArc);
+ }
+ else {
+ arc = k2->enteringTArc;
+ if( arc->tail == k2 )
+ arc->flow = arc->flow + theta;
+ else
+ arc->flow = arc->flow - theta;
+
+ k2 = Father(k2, k2->enteringTArc);
+ }
+ }
+ }
+
+ if( enteringArc != leavingArc ) {
+ bool leavingBringFlowInT2 = ( leavingReducesFlow ==
+ ( ( leavingArc->tail )->subTreeLevel > ( leavingArc->head )->subTreeLevel ) );
+ // "leavingBringFlowInT2" == true if leaving arc brings flow to the subtree T2
+ if( leavingBringFlowInT2 == ( memK1 == enteringArc->tail ) ) {
+ k2 = enteringArc->tail;
+ k1 = enteringArc->head;
+ }
+ else {
+ k2 = enteringArc->head;
+ k1 = enteringArc->tail;
+ }
+ }
+
+ #if( QUADRATICCOST == 0 )
+ if( leavingReducesFlow )
+ leavingArc->ident = AT_LOWER;
+ else
+ leavingArc->ident = AT_UPPER;
+
+ if( leavingArc != enteringArc ) {
+ enteringArc->ident = BASIC;
+ nodePType *h1;
+ nodePType *h2;
+ // "h1" is the node in the leaving arc with smaller tree's level
+ if( ( leavingArc->tail )->subTreeLevel < ( leavingArc->head )->subTreeLevel ) {
+ h1 = leavingArc->tail;
+ h2 = leavingArc->head;
+ }
+ else {
+ h1 = leavingArc->head;
+ h2 = leavingArc->tail;
+ }
+
+ UpdateT(leavingArc, enteringArc, h1, h2, k1, k2);
+ // Update potential of the subtree T2
+ k2 = enteringArc->head;
+ CNumber delta = ReductCost(enteringArc);
+ if( ( enteringArc->tail )->subTreeLevel > ( enteringArc->head )->subTreeLevel ) {
+ delta = -delta;
+ k2 = enteringArc->tail;
+ }
+
+ AddPotential( k2 , delta );
+ // In the linear case Primal Simplex only updates the potential of the nodes of
+ // subtree T2
+ }
+ #else
+ if( leavingArc != enteringArc ) {
+ nodePType *h1;
+ nodePType *h2;
+ // "h1" is the node in the leaving arc with smaller tree's level
+ if( ( leavingArc->tail )->subTreeLevel <
+ ( leavingArc->head )->subTreeLevel ) {
+ h1 = leavingArc->tail;
+ h2 = leavingArc->head;
+ }
+ else {
+ h1 = leavingArc->head;
+ h2 = leavingArc->tail;
+ }
+
+ // Update the basic tree
+ UpdateT( leavingArc , enteringArc , h1 , h2 , k1 , k2 );
+ }
+
+ #if( OPTQUADRATIC )
+ nodePType *h1;
+ nodePType *h2;
+ if( ( leavingArc->tail )->subTreeLevel <
+ ( leavingArc->head )->subTreeLevel ) {
+ h1 = leavingArc->tail;
+ h2 = leavingArc->head;
+ }
+ else {
+ h1 = leavingArc->head;
+ h2 = leavingArc->tail;
+ }
+
+ nodePType *node = h1;
+ nodePType *updateNode = h1;
+ if( h1 == cycleRoot )
+ ComputePotential( cycleRoot );
+ else {
+ while( node != cycleRoot ) {
+ arcPType *entArc = node->enteringTArc;
+ if( ! ETZ( entArc->quadraticCost , EpsCst ) )
+ updateNode = node;
+
+ node = Father( node , entArc );
+ }
+
+ ComputePotential( updateNode );
+ node = h2;
+ updateNode = h2;
+ while( node != cycleRoot ) {
+ arcPType *entArc = node->enteringTArc;
+ if( ! ETZ( entArc->quadraticCost , EpsCst ) )
+ updateNode = node;
+
+ node = Father( node , entArc );
+ }
+
+ ComputePotential( updateNode );
+ }
+ #else
+ // Update the potential of the node "manually"
+ ComputePotential( cycleRoot );
+ #endif
+
+ #if( LIMITATEPRECISION )
+ cont = cont + 1;
+ if( cont == recomputeFOLimits ) {
+ cont = 0;
+ foValue = GetFO(); // Calculate f.o. value manually
+ }
+ else
+ foValue = foValue + deltaFO;
+ // Calculate the f.o. value with the estimated decrease
+ #endif
+ #endif
+ }
+ else {
+ status = kOK;
+ // If one of dummy arcs has flow bigger than 0, the solution is unfeasible.
+ for( arcPType *arc = dummyArcsP ; arc != stopDummyP ; arc++ )
+ if( GTZ( arc->flow , EpsFlw ) )
+ status = kUnfeasible;
+ }
+
+ if( ( status == kUnSolved ) && MaxTime && MCFt ) {
+ double tu, ts;
+ TimeMCF( tu , ts );
+ if( MaxTime < tu + ts )
+ status = kStopped;
+ }
+
+ if( ( status == kUnSolved ) && MaxIter)
+ if( MaxIter < (int) iterator )
+ status = kStopped;
+
+ #if( UNIPI_PRIMAL_ITER_SHOW )
+ int it = (int) iterator;
+ if( it % UNIPI_PRIMAL_ITER_SHOW == 0 ) {
+ cout << endl;
+ for( int t = 0; t < 3; t++ )
+ cout << "\t";
+ cout << "PRIMALE MCFSimplex: ARCHI E NODI ALLA " << it << " ITERAZIONE" << endl;
+ ShowSituation( 3 );
+ #if( FOSHOW )
+ if( (int) iterator % FOSHOW == 0 )
+ clog << "Iteration = " << iterator << " of = "
+ #if( LIMITATEPRECISION && QUADRATICCOST )
+ << foValue
+ #else
+ << GetFO()
+ #endif
+ << endl;
+ #endif
+ }
+ #endif
+ }
+
+ #if( UNIPI_PRIMAL_FINAL_SHOW )
+ cout << endl;
+ for( int t = 0; t < 3; t++ )
+ cout << "\t";
+ cout << "PRIMALE UniPi: ARCHI E NODI ALLA FINE" << endl;
+ ShowSituation( 3 );
+ #endif
+
+ } // end( PrimalSimplex )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::DualSimplex( void )
+{
+ #if( UNIPI_PRIMAL_INITIAL_SHOW == 0 )
+ #if( UNIPI_PRIMAL_ITER_SHOW == 0 )
+ #if( UNIPI_PRIMAL_FINAL_SHOW == 0 )
+ cout << endl;
+ #endif
+ #endif
+ #endif
+ #if( UNIPI_DUAL_INITIAL_SHOW )
+ cout << endl;
+ for( int t = 0; t < 3; t++ )
+ cout << "\t";
+ cout << "DUALE MCFSimplex: ARCHI E NODI ALL' INIZIO" << endl;
+ ShowSituation( 3 );
+ #endif
+
+ if( pricingRule != kCandidateListPivot )
+ arcToStartD = arcsD;
+
+ iterator = 0;
+ arcDType *enteringArc;
+ arcDType *leavingArc;
+ if( pricingRule == kCandidateListPivot )
+ InitializeDualCandidateList();
+
+ status = kUnSolved;
+ while( status == kUnSolved ) {
+ iterator++;
+ if( iterator == Inf<iteratorType>() ) {
+ ResetWhenInT2(); // Restore to 0 every nodes' field "whenInT2"
+ iterator = 1;
+ }
+
+ switch( pricingRule ) {
+ case( kDantzig ):
+ leavingArc = DRuleFirstEligibleArc();
+ break; // si esegue cmq FEA
+ case( kFirstEligibleArc ):
+ leavingArc = DRuleFirstEligibleArc();
+ break;
+ case( kCandidateListPivot ):
+ leavingArc = RuleDualCandidateListPivot();
+ break;
+ }
+
+ if( pricingRule != kCandidateListPivot ) {
+ arcToStartD++;
+ if( arcToStartD == stopArcsD )
+ arcToStartD = dummyArcsD;
+
+ if( arcToStartD == stopDummyD )
+ arcToStartD = arcsD;
+ // Setting the initial arc for the Dantzig or First Elibigle Rule
+ }
+
+ if( leavingArc ) {
+ bool leavingArcInL = false;
+ bool leavingArcFromT1toT2;
+ if( LTZ( leavingArc->flow , EpsFlw ) )
+ leavingArcInL = true;
+
+ nodeDType *h1;
+ nodeDType *h2;
+ if( ( leavingArc->tail )->subTreeLevel <
+ ( leavingArc->head )->subTreeLevel ) {
+ h1 = leavingArc->tail;
+ h2 = leavingArc->head;
+ leavingArcFromT1toT2 = true;
+ }
+ else {
+ h1 = leavingArc->head;
+ h2 = leavingArc->tail;
+ leavingArcFromT1toT2 = false;
+ }
+
+ Index numOfT2Arcs = 0;
+ int level = h2->subTreeLevel;
+ nodeDType *node = h2;
+ node->whenInT2 = iterator;
+ nodeDType *lastNodeOfT2 = h2;
+ numOfT2Arcs = node->numArcs;
+ while( node->nextInT && ( ( node->nextInT )->subTreeLevel > level ) ) {
+ node = node->nextInT;
+ lastNodeOfT2 = node;
+ numOfT2Arcs = numOfT2Arcs + node->numArcs;
+ node->whenInT2 = iterator;
+ }
+
+ /* The Dual Simplex has determinated the leaving arc, and so the
+ subtrees T1 and T2. Dual Simplex scans T2 to fix the fields "whenInT2"
+ of T2's nodes to the iteration value, and counts the entering and
+ outgoing arcs from these nodes. According to this number, it decides
+ to scan the Backward and Forward of the subtree (T1 or T2) with the
+ minor number of entering/outgoing arcs from its nodes. */
+ enteringArc = NULL;
+ bool lv = ( leavingArcFromT1toT2 == leavingArcInL );
+ CNumber maxRc = Inf<CNumber>();
+ //Search arc in the Forward Star and Backward Star of nodes of T1
+ if( numOfT2Arcs > m ) {
+ // Dual Simplex starts from the node which follows the dummy root.
+ node = dummyRootD->nextInT;
+ bool fine = false;
+ 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( lastNodeOfT2->nextInT )
+ node = lastNodeOfT2->nextInT;
+ else
+ break;
+
+ // Search arc in the Backward Star of nodes of T1
+ arcDType *arc = node->firstBs;
+ while( arc ) {
+ if( ( arc->tail )->whenInT2 == iterator ) {
+ // Evaluate each arc from T2 to T1 which isn't in T
+ if( arc->ident == AT_LOWER ) {
+ if( lv ) {
+ CNumber rc = ABS( ReductCost( arc ) );
+ if( LT( rc , maxRc , EpsCst ) ) {
+ enteringArc = arc;
+ maxRc = rc;
+ /* If arc is appropriate to enter in T and its reduct cost is 0,
+ search is ended: this is the arc which enters in T */
+ if( ETZ( rc , EpsCst) ) {
+ fine = true;
+ break;
+ }
+ }
+ }
+ }
+
+ if( arc->ident == AT_UPPER ) {
+ if( ! lv ) {
+ CNumber rc = ABS( ReductCost( arc ) );
+ if( LT( rc , maxRc , EpsCst ) ) {
+ enteringArc = arc;
+ maxRc = rc;
+ /* If arc is appropriate to enter in T and its reduct cost is 0,
+ search is ended: this is the arc which enters in T */
+
+ if( ETZ( rc , EpsCst ) ) {
+ fine = true;
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ arc = arc->nextBs;
+ }
+
+ // Search arc in the Forward Star of nodes of T1
+ arc = node->firstFs;
+ while( arc ) {
+ if( ( arc->head )->whenInT2 == iterator ) {
+ // Evaluate each arc from T1 to T2 which isn't in T
+ if( arc->ident == AT_LOWER ) {
+ if( ! lv ) {
+ CNumber rc = ABS( ReductCost( arc ) );
+ if( LT( rc , maxRc , EpsCst ) ) {
+ enteringArc = arc;
+ maxRc = rc;
+ /* If arc is appropriate to enter in T and its reduct cost is 0,
+ search is ended: this is the arc which enters in T */
+ if( ETZ( rc , EpsCst ) ) {
+ fine = true;
+ break;
+ }
+ }
+ }
+ }
+
+ if( arc->ident == AT_UPPER ) {
+ if( lv ) {
+ CNumber rc = ABS( ReductCost( arc ) );
+ if( LT( rc , maxRc , EpsCst ) ) {
+ enteringArc = arc;
+ maxRc = rc;
+ /* If arc is appropriate to enter in T and its reduct cost is 0,
+ search is ended: this is the arc which enters in T */
+ if( ETZ( rc , EpsCst ) ) {
+ fine = true;
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ arc = arc->nextFs;
+ }
+
+ node = node->nextInT;
+ if( node == NULL )
+ fine = true;
+ }
+ }
+ // Search arc in the Forward Star and Backward Star of nodes of T2
+ else {
+ node = h2;
+ bool fine = false;
+ while( fine == false ) {
+ // Search arc in the Backward Star of nodes of T2
+ arcDType *arc = node->firstBs;
+ CNumber rc;
+ while( arc ) {
+ if( ( arc->tail )->whenInT2 != iterator ) {
+ // Evaluate each arc from T1 to T2 which isn't in T
+ if( arc->ident == AT_LOWER ) {
+ if( ! lv ) {
+ rc = ABS( ReductCost( arc ) );
+ if( LT( rc , maxRc , EpsCst ) ) {
+ enteringArc = arc;
+ maxRc = rc;
+ /* If arc is appropriate to enter in T and its reduct cost is 0,
+ search is ended: this is the arc which enters in T */
+ if( ETZ( rc , EpsCst ) ) {
+ fine = true;
+ break;
+ }
+ }
+ }
+ }
+
+ if( arc->ident == AT_UPPER ) {
+ if( lv ) {
+ rc = ABS( ReductCost( arc ) );
+ if( LT( rc , maxRc , EpsCst ) ) {
+ enteringArc = arc;
+ maxRc = rc;
+ /* If arc is appropriate to enter in T and its reduct cost is 0,
+ search is ended: this is the arc which enters in T */
+ if( ETZ( rc , EpsCst ) ) {
+ fine = true;
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ arc = arc->nextBs;
+ }
+
+ // Search arc in the Forward Star of nodes of T2
+ arc = node->firstFs;
+ while( arc ) {
+ if( ( arc->head )->whenInT2 != iterator ) {
+ // Evaluate each arc from T2 to T1 which isn't in T
+ if( arc->ident == AT_LOWER ) {
+ if( lv ) {
+ rc = ABS( ReductCost( arc ) );
+ if( LT( rc , maxRc , EpsCst ) ) {
+ enteringArc = arc;
+ maxRc = rc;
+ /* If arc is appropriate to enter in T and its reduct cost is 0,
+ search is ended: this is the arc which enters in T */
+ if( ETZ( rc , EpsCst ) ) {
+ fine = true;
+ break;
+ }
+ }
+ }
+ }
+
+ if( arc->ident == AT_UPPER ) {
+ if( ! lv ) {
+ rc = ABS( ReductCost( arc ) );
+ if( LT( rc , maxRc , EpsCst ) ) {
+ enteringArc = arc;
+ maxRc = rc;
+ /* If arc is appropriate to enter in T and its reduct cost is 0,
+ search is ended: this is the arc which enters in T */
+ if( ETZ( rc , EpsCst) ) {
+ fine = true;
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ arc = arc->nextFs;
+ }
+
+ if( node == lastNodeOfT2 )
+ fine = true;
+ else
+ node = node->nextInT;
+
+ }
+ }
+
+ if( enteringArc ) {
+ FNumber theta = -leavingArc->flow;
+ if( GTZ( leavingArc->flow , EpsFlw ) )
+ 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
+ determinated by T and entering arc for decreases flow in the entering arc:
+ if entering arc is in L, Dual Simplex pushs flow in the cycle
+ determinated by T and entering arc for increases flow in the entering arc:
+ Dual Simplex increases or decreases entering arc's flow.
+ theta is a positive value.
+ For this reason the algorithm uses two nodes (k1 and k2) to push flow
+ (theta) from k1 to k2. According to entering arc's reduct cost, the
+ algorithm determinates k1 and k2 */
+
+ if( enteringArc->ident == AT_UPPER ) {
+ k1 = enteringArc->head;
+ k2 = enteringArc->tail;
+ }
+ else {
+ k1 = enteringArc->tail;
+ k2 = enteringArc->head;
+ }
+
+ nodeDType *memK1 = k1;
+ nodeDType *memK2 = k2;
+ arcDType *arc;
+ k1 = memK1;
+ k2 = memK2;
+ // Update the flow
+ while( k1 != k2 ) {
+ if( k1->subTreeLevel > k2->subTreeLevel ) {
+ arc = k1->enteringTArc;
+ if( arc->tail != k1 )
+ arc->flow = arc->flow + theta;
+ else
+ arc->flow = arc->flow - theta;
+
+ k1 = Father(k1, k1->enteringTArc);
+ }
+ else {
+ arc = k2->enteringTArc;
+ if( arc->tail == k2 )
+ arc->flow = arc->flow + theta;
+ else
+ arc->flow = arc->flow - theta;
+
+ k2 = Father( k2 , k2->enteringTArc );
+ }
+ }
+
+ if(leavingArcInL )
+ leavingArc->ident = AT_LOWER;
+ else
+ leavingArc->ident = AT_UPPER;
+
+ bool leavingBringFlowInT2 = ( leavingArcInL ==
+ ( ( leavingArc->tail )->subTreeLevel >
+ ( leavingArc->head )->subTreeLevel ) );
+ // leavingBringFlowInT2 == true if leaving arc brings flow to the subtree T2
+ if( leavingBringFlowInT2 != ( memK1 == enteringArc->tail ) ) {
+ k2 = enteringArc->tail;
+ k1 = enteringArc->head;
+ }
+ else {
+ k2 = enteringArc->head;
+ k1 = enteringArc->tail;
+ }
+
+ if( enteringArc->ident == AT_LOWER )
+ enteringArc->flow = enteringArc->flow + theta;
+ else
+ enteringArc->flow = enteringArc->flow - theta;
+
+ enteringArc->ident = BASIC;
+ UpdateT( leavingArc , enteringArc , h1 , h2 , k1 , k2 );
+ // update potential of the subtree T2
+ k2 = enteringArc->head;
+ CNumber delta = ReductCost( enteringArc );
+ if( ( enteringArc->tail) ->subTreeLevel >
+ ( enteringArc->head )->subTreeLevel ) {
+ delta = -delta;
+ k2 = enteringArc->tail;
+ }
+
+ // Dual Simplex only updates the potential of the T2's nodes
+ AddPotential( k2 , delta );
+ }
+ else
+ status = kUnfeasible;
+ /* If Dual Simplex finds a leaving arc but it doesn't find an entering arc,
+ the algorithm stops. At this point Dual Simplex has an unfeasible primal
+ solution. */
+ }
+ else {
+ status = kOK;
+ // If one of dummy arcs has flow different than 0, the solution is unfeasible.
+ for( arcDType *arc = dummyArcsD ; arc != stopDummyD ; arc++ )
+ if( ! ETZ( arc->flow , EpsFlw ) ) {
+ status = kUnfeasible;
+ break;
+ }
+ }
+
+ if( ( status == kUnSolved ) && MaxTime ) {
+ double tu, ts;
+ TimeMCF( tu , ts );
+ if( MaxTime < tu + ts )
+ status = kStopped;
+ }
+
+ if( ( status == kUnSolved ) && MaxIter && MCFt )
+ if( MaxIter < (int) iterator )
+ status = kStopped;
+
+ #if( UNIPI_DUAL_ITER_SHOW )
+ if( (int) iterator % UNIPI_DUAL_ITER_SHOW == 0 ) {
+ cout << endl;
+ for( int t = 0; t < 3; t++ )
+ cout << "\t";
+ cout << "DUALE MCFSimplex: ARCHI E NODI ALLA " << iterator << " ITERAZIONE"
+ << endl;
+ ShowSituation( 3 );
+ #if( FOSHOW )
+ cout << "of = " << GetFO() << endl;
+ #endif
+ }
+ #endif
+ }
+
+ #if( UNIPI_DUAL_ITER_SHOW )
+ int it = (int) iterator;
+ if( it % UNIPI_DUAL_ITER_SHOW == 0 ) {
+ cout << endl;
+ for( int t = 0; t < 3; t++ )
+ cout << "\t";
+ cout << "DUALE MCFSimplex: ARCHI E NODI ALLA " << iterator << " ITERAZIONE"
+ << endl;
+ Showsituation( 3 );
+ }
+ #endif
+
+ } // end( DualSimplex )
+
+/*--------------------------------------------------------------------------*/
+
+template<class N, class A>
+void MCFSimplex::UpdateT( A *h , A *k , N *h1 , 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
+ coming root of T2). With the update of T, this path will be overturned:
+ node k2 will become the root of T2...
+ The subtree T2 must be reordered and the field "subTreeLevel", which
+ represents the depth in T of every node, of every T2's nodes is changed.
+ Variable delta represents the increase of "subTreeLevel" value for node
+ k2 and its descendants: probably this value is a negative value. */
+
+ int delta = (k1->subTreeLevel) + 1 - (k2->subTreeLevel);
+ N *root = k2;
+ N *dad;
+
+ /*To reorder T2, the method analyses every nodes of the path h2->k2, starting
+ from k2. For every node, it moves the node's subtree from its original
+ position to a new appropriate position. In particular k2 and its subtree
+ (k2 is the new root of T2, so the first nodes of the new T2) will be moved
+ next to node k1 (new father of k2), the next subtree will be moved beside
+ the last node of k2's subtree....
+ "previousNode" represents the node where the new subtree will be moved
+ beside in this iterative action. At the start "previousNode" is the node
+ k1 (T2 will be at the right of k1). */
+
+ N *previousNode = k1;
+ N *lastNode;
+ /* "arc1" is the entering arc in T (passed by parameters).
+ For every analysed node of path h2->k2, the method changes
+ "enteringTArc" but it must remember the old arc, which will be the
+ "enteringTArc" of the next analysed node.
+ At the start "arc1" is k (the new "enteringTArc" of k2). */
+
+ A *arc1 = k;
+ A *arc2;
+ bool fine = false;
+ while( fine == false ) {
+ // If analysed node is h2, this is the last iteration
+ if( root == h2 )
+ fine = true;
+
+ dad = Father( root , root->enteringTArc );
+ // Cut the root's subtree from T and update the "subLevelTree" of its nodes
+ lastNode = CutAndUpdateSubtree( root , delta );
+ // Paste the root's subtree in the right position;
+ PasteSubtree( root , lastNode , previousNode );
+ // In the next iteration the subtree will be moved beside the last node of
+ // the actual analysed subtree.
+ previousNode = lastNode;
+ /* The increase of the subtree in the next iteration is different from
+ the actual increase: in particual the increase increases itself (+2
+ at every iteration). */
+ delta = delta + 2;
+ /* A this point "enteringTArc" of actual root is stored in "arc2" and
+ changed; then "arc1" and "root" are changed. */
+ arc2 = root->enteringTArc;
+ root->enteringTArc = arc1;
+ arc1 = arc2;
+ root = dad;
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N>
+N* MCFSimplex::CutAndUpdateSubtree( N *root , int delta )
+{
+ int level = root->subTreeLevel;
+ N *node = root;
+ // The root of this subtree is passed by parameters, the last node is searched.
+ while ( ( node->nextInT ) && ( ( node->nextInT )->subTreeLevel > level ) ) {
+ node = node->nextInT;
+ // The "subTreeLevel" of every nodes of subtree is updated
+ node->subTreeLevel = node->subTreeLevel + delta;
+ }
+
+ root->subTreeLevel = root->subTreeLevel + delta;
+ /* The 2 neighbouring nodes of the subtree (the node at the left of the root
+ and the node at the right of the last node) is joined. */
+
+ if( root->prevInT )
+ ( root->prevInT )->nextInT = node->nextInT;
+ if( node->nextInT )
+ ( node->nextInT )->prevInT = root->prevInT;
+
+ return( node ); // the method returns the last node of the subtree
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N>
+void MCFSimplex::PasteSubtree( N *root , N *lastNode , N *previousNode )
+{
+ /* The method inserts subtree ("root" and "lastNode" are the extremity of the
+ subtree) after "previousNode". The method starts to identify the next node
+ of "previousNode" ("nextNode"), so it joins "root" with "previousNode" and
+ "lastNode" with "nextNode" (if exists). */
+
+ N *nextNode = previousNode->nextInT;
+ root->prevInT = previousNode;
+ previousNode->nextInT = root;
+ lastNode->nextInT = nextNode;
+ if( nextNode )
+ nextNode->prevInT = lastNode;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcPType* MCFSimplex::RuleDantzig( void )
+{
+ arcPType *arc = arcToStartP;
+ arcPType *enteringArc = NULL;
+ #if( QUADRATICCOST )
+ /* In the quadratic case used type for reduct cost is FONumber.
+ Value "lim" is the fixed thresold for the decrease of the f.o. value */
+ FONumber lim = EpsOpt * foValue / n;
+ FONumber RC;
+ FONumber maxValue = 0;
+ #else
+ CNumber RC;
+ CNumber maxValue = 0;
+ #endif
+
+ do {
+ // The method analyses every arc
+ #if( QUADRATICCOST )
+ RC = ReductCost( arc );
+ FNumber theta;
+ /* If reduct cost of arc is lower than 0, the flow of the arc must increase.
+ If reduct cost of arc is bigger than 0, the flow of the arc must decrease.
+ "theta" is the difference between lower (upper) bound and the actual flow.
+ */
+
+ if( LTZ( RC , EpsCst ) )
+ theta = arc->upper - arc->flow;
+
+ if( GTZ( RC , EpsCst ) )
+ theta = -arc->flow;
+
+ // If it's possible to increase (or decrease) the flow in this arc
+ if( ! ETZ( theta , EpsFlw ) ) {
+ /* "Q" is the sum of the quadratic coefficient of the arc belonging the
+ T path from tail's arc to head's arc
+ "Q" is always bigger than 0 or equals to 0.
+ If "Q" > 0, the value - RC / Q is the increase (decrease) of the flow
+ with the best decrease of f.o. value.
+ - RC/ Q must be compare with "theta" to avoid that the best increase
+ (decrease) of the flow violates the bounds of the arc.
+ This confront determines "theta". */
+
+ CNumber Q = ( arc->tail )->sumQuadratic + ( arc->head )->sumQuadratic +
+ arc->quadraticCost;
+
+ if( GTZ( Q , EpsCst ) )
+ if( GTZ( theta , EpsFlw ) )
+ theta = min( theta , - RC / Q );
+ else
+ theta = max( theta , - RC / Q );
+
+ /* Calculate the estimate decrease of the f.o. value if this arc is
+ selected and flow is increased (decreased) by "theta" */
+
+ CNumber deltaFO = RC * theta + Q * theta * theta / 2;
+ // if deltaFO < 0 this arc is appropriate; if deltaFO is lower than
+ // old decrease value, arc is the best arc.
+
+ if( deltaFO < maxValue ) {
+ maxValue = deltaFO;
+ enteringArc = arc;
+ }
+ }
+ #else
+ if( arc->ident > BASIC ) {
+ RC = ReductCost( arc );
+
+ if( ( LTZ( RC , EpsCst ) && ( arc->ident == AT_LOWER ) ) ||
+ ( GTZ( RC , EpsCst ) && ( arc->ident == AT_UPPER ) ) ) {
+
+ if( RC < 0 )
+ RC = -RC;
+
+ if( RC > maxValue ) {
+ maxValue = RC;
+ enteringArc = arc;
+ }
+ }
+ }
+ #endif
+
+ arc++;
+ if( arc == stopArcsP )
+ arc = arcsP;
+
+ } while( arc != arcToStartP );
+
+ #if( ( LIMITATEPRECISION ) && ( QUADRATICCOST ) )
+ if( -maxValue <= lim )
+ enteringArc = NULL;
+ #endif
+
+ return( enteringArc );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcPType* MCFSimplex::PRuleFirstEligibleArc( void )
+{
+ arcPType *arc = arcToStartP;
+ arcPType *enteringArc = NULL;
+ #if( QUADRATICCOST )
+ FONumber RC;
+ #else
+ CNumber RC;
+ #endif
+
+ do {
+ #if( QUADRATICCOST )
+ // In this method the "decrease f.o. value" criterion is not used.
+ RC = ReductCost( arc );
+ if( ( LTZ( RC , EpsCst ) && LT( arc->flow , arc->upper , EpsFlw ) ) ||
+ ( GTZ( RC , EpsCst ) && GTZ( arc->flow , EpsFlw ) ) )
+ enteringArc = arc;
+ #else
+ if( arc->ident > BASIC ) {
+ RC = ReductCost( arc );
+ if( ( LTZ( RC , EpsCst ) && ( arc->ident == AT_LOWER ) ) ||
+ ( GTZ( RC , EpsCst ) && ( arc->ident == AT_UPPER ) ) )
+ enteringArc = arc;
+ }
+ #endif
+
+ arc++;
+ if( arc == stopArcsP )
+ arc = dummyArcsP;
+ if( arc == stopDummyP )
+ arc = arcsP;
+
+ } while( ( enteringArc == NULL ) && ( arc != arcToStartP ) );
+
+ return( enteringArc );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcDType* MCFSimplex::DRuleFirstEligibleArc( void )
+{
+ arcDType *arc = arcToStartD;
+ arcDType *leavingArc = NULL;
+ do {
+ if( GT( arc->flow , arc->upper , EpsFlw ) || LTZ( arc->flow , EpsFlw ) )
+ leavingArc = arc;
+
+ arc++;
+ if( arc == stopArcsD )
+ arc = dummyArcsD;
+ if( arc == stopDummyD )
+ arc = arcsD;
+
+ } while( ( leavingArc == NULL ) && ( arc != arcToStartD ) );
+
+ return(leavingArc);
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcPType* MCFSimplex::RulePrimalCandidateListPivot( void )
+{
+ Index next = 0;
+ Index i;
+ Index minimeValue;
+ if( hotListSize < tempCandidateListSize )
+ minimeValue = hotListSize;
+ else
+ minimeValue = tempCandidateListSize;
+
+ #if( QUADRATICCOST )
+ // Check if the left arcs in the list continue to violate the dual condition
+ for( i = 2 ; i <= minimeValue ; i++ ) {
+ arcPType *arc = candP[ i ].arc;
+ FONumber red_cost = ReductCost( arc );
+ FNumber theta = 0;
+ /* If reduct cost of arc is lower than 0, the flow of the arc must increase.
+ If reduct cost of arc is bigger than 0, the flow of the arc must decrease.
+ "theta" is the difference between lower (upper) bound and the actual flow.
+ */
+
+ if( LTZ( red_cost , EpsCst ) )
+ theta = arc->upper - arc->flow;
+ else
+ if( GTZ( red_cost , EpsCst ) )
+ theta = - arc->flow;
+
+ // if it's possible to increase (or decrease) the flow in this arc
+ if( theta != 0 ) {
+ /* "Q" is the sum of the quadratic coefficient of the arc belonging the T
+ path from tail's arc to head's arc
+ "Q" is always bigger than 0 or equals to 0.
+ If "Q" > 0, the value - RC / Q is the increase (decrease) of the flow
+ with the best decrease of f.o. value.
+ - RC/ Q must be compare with "theta" to avoid that the best increase
+ (decrease) of the flow violates the bounds of the arc.
+ This confront determines "theta". */
+
+ CNumber Q = ( arc->tail )->sumQuadratic + ( arc->head )->sumQuadratic +
+ arc->quadraticCost;
+
+ if( GTZ( Q , EpsCst ) )
+ if( GTZ( theta , EpsFlw ) )
+ theta = min( theta , - red_cost / Q );
+ else
+ theta = max( theta , - red_cost / Q );
+
+ /* Calculate the estimate decrease of the f.o. value if this arc is
+ selected and flow is increased (decreased) by "theta" */
+
+ CNumber deltaFO = red_cost * theta + Q * theta * theta / 2;
+ #if( LIMITATEPRECISION )
+ // if deltaFO < 0 this arc is appropriate; if deltaFO is lower than
+ // old decrease value, arc is the best arc.
+ if( - deltaFO > ( EpsOpt * foValue / n ) ) {
+ #else
+ if( LTZ( deltaFO , EpsCst ) ) {
+ #endif
+ next++;
+ candP[ next ].arc = arc;
+ candP[ next ].absRC = -deltaFO;
+ }
+ }
+ }
+
+ tempCandidateListSize = next;
+ Index oldGroupPos = groupPos;
+ // Search other arcs to fill the list
+ do {
+ arcPType *arc;
+ for( arc = arcsP + groupPos ; arc < stopArcsP ; arc += numGroup ) {
+ FONumber red_cost = ReductCost( arc );
+ FNumber theta = 0;
+ /* If reduced cost is lower than 0, the flow of the arc must increase.
+ If reduced cost is larger than 0, the flow of the arc must decrease.
+ "theta" is the difference between lower (upper) bound and the actual
+ flow. */
+
+ if( LTZ( red_cost , EpsCst ) )
+ theta = arc->upper - arc->flow;
+ else
+ if( GTZ( red_cost , EpsCst ) )
+ theta = - arc->flow;
+
+ // if it's possible to increase (or decrease) the flow in this arc
+ if( theta != 0 ) {
+ /* "Q" is the sum of the quadratic coefficient of the arc belonging the T
+ path from tail's arc to head's arc
+ "Q" is always bigger than 0 or equals to 0.
+ If "Q" > 0, the value - RC / Q is the increase (decrease) of the flow
+ with the best decrease of f.o. value.
+ - RC/ Q must be compare with "theta" to avoid that the best increase
+ (decrease) of the flow violates the bounds of the arc.
+ This confront determines "theta". */
+ CNumber Q = ( arc->tail )->sumQuadratic + ( arc->head )->sumQuadratic +
+ arc->quadraticCost;
+
+ if( GTZ( Q , EpsCst ) )
+ if( GTZ( theta , EpsFlw ) )
+ theta = min( theta , - red_cost / Q );
+ else
+ theta = max( theta , - red_cost / Q );
+
+ /* Calculate the estimate decrease of the f.o. value if this arc is
+ selected and flow is increased (decreased) by "theta" */
+
+ CNumber deltaFO = red_cost * theta + Q * theta * theta / 2;
+ #if( LIMITATEPRECISION )
+ // if deltaFO < 0 this arc is appropriate; if deltaFO is lower than
+ // old decrease value, arc is the best arc.
+ if( -deltaFO > ( EpsOpt * foValue / n ) ) {
+ #else
+ if( LTZ( deltaFO , EpsCst ) ) {
+ #endif
+ tempCandidateListSize++;
+ candP[ tempCandidateListSize ].arc = arc;
+ candP[ tempCandidateListSize ].absRC = -deltaFO;
+ }
+ }
+ }
+
+ groupPos++;
+ if( groupPos == numGroup )
+ groupPos = 0;
+
+ } while( ( tempCandidateListSize < hotListSize ) &&
+ ( groupPos != oldGroupPos ) );
+ #else
+ // Check if the left arcs in the list continue to violate the dual condition
+ for( i = 2 ; i <= minimeValue ; i++ ) {
+ arcPType *arc = candP[i].arc;
+ CNumber red_cost = ReductCost( arc );
+
+ if( ( LTZ( red_cost , EpsCst ) && ( arc->ident == AT_LOWER ) ) ||
+ ( GTZ( red_cost , EpsCst ) && ( arc->ident == AT_UPPER ) ) ) {
+ next++;
+ candP[ next ].arc = arc;
+ candP[ next ].absRC = ABS( red_cost );
+ }
+ }
+
+ tempCandidateListSize = next;
+ Index oldGroupPos = groupPos;
+ // Search other arcs to fill the list
+ do {
+ arcPType *arc;
+ for( arc = arcsP + groupPos ; arc < stopArcsP ; arc += numGroup ) {
+ if( arc->ident == AT_LOWER ) {
+ CNumber red_cost = ReductCost( arc );
+ if( LTZ( red_cost , EpsCst ) ) {
+ tempCandidateListSize++;
+ candP[ tempCandidateListSize ].arc = arc;
+ candP[ tempCandidateListSize ].absRC = ABS( red_cost );
+ }
+ }
+ else
+ if( arc->ident == AT_UPPER ) {
+ CNumber red_cost = ReductCost( arc );
+ if( GTZ( red_cost , EpsCst ) ) {
+ tempCandidateListSize++;
+ candP[ tempCandidateListSize ].arc = arc;
+ candP[ tempCandidateListSize ].absRC = ABS( red_cost );
+ }
+ }
+ }
+
+ groupPos++;
+ if( groupPos == numGroup )
+ groupPos = 0;
+
+ } while( ( tempCandidateListSize < hotListSize ) && ( groupPos != oldGroupPos ) );
+ #endif
+
+ if( tempCandidateListSize ) {
+ SortPrimalCandidateList( 1 , tempCandidateListSize );
+ return( candP[ 1 ].arc );
+ }
+ else
+ return( NULL );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::InitializePrimalCandidateList( void )
+{
+ numGroup = ( ( m - 1 ) / numCandidateList ) + 1;
+ groupPos = 0;
+ tempCandidateListSize = 0;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::SortPrimalCandidateList( Index min , Index max )
+{
+ Index left = min;
+ Index right = max;
+ #if( QUADRATICCOST )
+ FONumber cut = candP[ ( left + right ) / 2 ].absRC;
+ #else
+ CNumber cut = candP[ ( left + right ) / 2 ].absRC;
+ #endif
+ do {
+ while( candP[ left ].absRC > cut)
+ left++;
+ while( cut > candP[ right ].absRC)
+ right--;
+
+ if( left < right )
+ Swap( candP[ left ] , candP[ right ] );
+
+ if(left <= right) {
+ left++;
+ right--;
+ }
+ } while( left <= right );
+
+ if( min < right )
+ SortPrimalCandidateList( min , right );
+ if( ( left < max ) && ( left <= hotListSize ) )
+ SortPrimalCandidateList( left , max );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcDType* MCFSimplex::RuleDualCandidateListPivot( void )
+{
+ Index next = 0;
+ // Check if the left arcs in the list continue to violate the primal condition
+ for( Index i = 2 ; ( i <= hotListSize ) && ( i <= tempCandidateListSize ) ;
+ i++ ) {
+ nodeDType *node = candD[ i ].node;
+ arcDType *arc = node->enteringTArc;
+ cFNumber flow = arc->flow;
+ if( LTZ( flow , EpsFlw ) ) {
+ next++;
+ candD[ next ].node = node;
+ candD[ next ].absInfeas = ABS( flow );
+ }
+
+ if( GT( flow , arc->upper , EpsFlw ) ) {
+ next++;
+ candD[ next ].node = node;
+ candD[ next ].absInfeas = flow - arc->upper;
+ }
+ }
+
+ tempCandidateListSize = next;
+ Index oldGroupPos = groupPos;
+ // Search other arcs to fill the list
+ do {
+ nodeDType *node = nodesD + groupPos;
+ for( node ; node < stopNodesD ; node += numGroup ) {
+ arcDType *arc = node->enteringTArc;
+ cFNumber flow = arc->flow;
+ if( LTZ( flow , EpsFlw ) ) {
+ tempCandidateListSize++;
+ candD[ tempCandidateListSize ].node = node;
+ candD[ tempCandidateListSize ].absInfeas = ABS( flow );
+ }
+
+ if( GT( flow , arc->upper , EpsFlw) ) {
+ tempCandidateListSize++;
+ candD[ tempCandidateListSize ].node = node;
+ candD[ tempCandidateListSize ].absInfeas = flow - arc->upper;
+ }
+ }
+
+ groupPos++;
+ if( groupPos == numGroup )
+ groupPos = 0;
+
+ } while( ( tempCandidateListSize < hotListSize ) &&
+ ( groupPos != oldGroupPos ) );
+
+ if( tempCandidateListSize ) {
+ SortDualCandidateList( 1 , tempCandidateListSize );
+ return( (candD[ 1 ].node)->enteringTArc );
+ }
+ else
+ return( NULL );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::InitializeDualCandidateList( void )
+{
+ numGroup = ( ( n - 1 ) / numCandidateList ) + 1;
+ groupPos = 0;
+ tempCandidateListSize = 0;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::SortDualCandidateList(Index min, Index max)
+{
+ Index left = min;
+ Index right = max;
+ FNumber cut = candD[ ( left + right ) / 2 ].absInfeas;
+ do {
+ while( candD[ left ].absInfeas > cut )
+ left++;
+ while( cut > candD[ right ].absInfeas )
+ right--;
+ if( left < right )
+ Swap( candD[left ] , candD[ right ] );
+ if( left <= right) {
+ left++;
+ right--;
+ }
+ } while( left <= right );
+
+ if( min < right )
+ SortDualCandidateList( min , right );
+ if( (left < max) && ( left <= hotListSize ) )
+ SortDualCandidateList( left , max );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N, class RCT>
+inline void MCFSimplex::AddPotential( N *r , RCT delta )
+{
+ int level = r->subTreeLevel;
+ N *n = r;
+
+ do {
+ n->potential = n->potential + delta;
+ n = n->nextInT;
+ } while ( ( n ) && ( n->subTreeLevel > level ) );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N>
+inline void MCFSimplex::ComputePotential( N *r )
+{
+ N *n = r;
+ int level = r->subTreeLevel;
+ FONumber cost;
+ // If "n" is not the dummy root, the potential of "r" is computed.
+ // If "n" is the dummy root, the potential of dummy root is a constant.
+
+ do {
+ if( n->enteringTArc ) {
+ cost = ( n->enteringTArc )->cost;
+ #if (QUADRATICCOST)
+ // Also field "sumQuadratic" is updated
+ n->sumQuadratic = ( Father( n , n->enteringTArc ) )->sumQuadratic +
+ ( n->enteringTArc )->quadraticCost;
+
+ if( ! ETZ( ( n->enteringTArc )->flow , EpsFlw ) )
+ cost = cost + ( ( n->enteringTArc )->quadraticCost * ( n->enteringTArc )->flow );
+ #endif
+
+ if( n == ( n->enteringTArc )->head )
+ n->potential = ( Father( n , n->enteringTArc ) )->potential + cost;
+ else
+ n->potential = ( Father( n , n->enteringTArc ) )->potential - cost;
+ }
+ n = n->nextInT;
+ } while( ( n ) && ( n->subTreeLevel > level ) );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateInitialPModifiedBalanceVector( void )
+{
+ int i = 0;
+ delete[] modifiedBalance;
+ modifiedBalance = new FNumber[ n ];
+ // Initialited every node's modifiedBalance to his balance
+ for ( nodePType *node = nodesP ; node != stopNodesP ; node++ ) {
+ modifiedBalance[i] = node->balance;
+ i++;
+ }
+
+ // Modify the vector according to the arcs out of base with flow non zero
+ // Scan the real arcs
+
+ for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++ ) {
+ #if( QUADRATICCOST )
+ if( ( ! ETZ( arc->flow , EpsFlw ) ) &&
+ ( ( arc->tail )->enteringTArc != arc ) &&
+ ( ( arc->head )->enteringTArc != arc ) ) {
+ i = (arc->tail) - nodesP;
+ modifiedBalance[ i ] += arc->flow;
+ i = (arc->head) - nodesP;
+ modifiedBalance[ i ] -= arc->flow;
+ }
+ #else
+ if( arc->ident == AT_UPPER ) {
+ i = (arc->tail) - nodesP;
+ modifiedBalance[ i ] += arc->upper;
+ i = (arc->head) - nodesP;
+ modifiedBalance[ i ] -= arc->upper;
+ }
+ #endif
+ }
+
+ // Scan the dummy arcs
+ for( arcPType *arc = dummyArcsP ; arc != stopDummyP ; arc++ ) {
+ #if( QUADRATICCOST )
+ if ( ( ! ETZ( arc->flow , EpsFlw ) ) &&
+ ( ( arc->tail )->enteringTArc != arc ) &&
+ ( ( arc->head )->enteringTArc != arc ) ) {
+ i = (arc->tail) - nodesP;
+ modifiedBalance[ i ] += arc->flow;
+ i = (arc->head) - nodesP;
+ modifiedBalance[ i ] -= arc->flow;
+ }
+ #else
+ if (arc->ident == AT_UPPER) {
+ i = (arc->tail) - nodesP;
+ modifiedBalance[ i ] += arc->upper;
+ i = (arc->head) - nodesP;
+ modifiedBalance[ i ] -= arc->upper;
+ }
+ #endif
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+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 )->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"
+ ( 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
+ nodePType *desc = r->nextInT;
+ // Call PostPVisit for every child of "r"
+ while( ( desc ) && ( desc->subTreeLevel > r->subTreeLevel ) ) {
+ if( desc->subTreeLevel - 1 == r->subTreeLevel ) { // desc is a son of r
+ PostPVisit( desc );
+
+ if( ( desc->enteringTArc )->head == r ) // enteringTArc of desc goes in r
+ modifiedBalance[ i ] -= ( desc->enteringTArc )->flow;
+ else // If enteringTArc of "desc" goes out "r"
+ modifiedBalance[ i ] += ( desc->enteringTArc )->flow;
+ }
+ desc = desc->nextInT;
+ }
+
+ 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 ];
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::BalanceFlow( nodePType *r )
+{
+ // used only by Primal Simplex to restore a primal feasible solution.
+ if( r == dummyRootP ) {
+ nodePType *node = dummyRootP->nextInT;
+ while( node ) {
+ // call this function recursively for every son of dummy root
+ if( node->subTreeLevel == 1 )
+ BalanceFlow( node );
+
+ node = node->nextInT;
+ }
+ }
+ else {
+ // The method controls if "r" is a leaf in T
+ bool rLeaf = false;
+ 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
+ // not feasible; in case adjust its flow
+ else { // If "r" isn't a leaf
+ nodePType *node = r->nextInT;
+ // Balance the flow of every child of "r"
+ while ( ( node ) && ( node->subTreeLevel > r->subTreeLevel ) ) {
+ if( node->subTreeLevel == r->subTreeLevel + 1 )
+ BalanceFlow( node );
+ node = node->nextInT;
+ }
+
+ // The method controls if entering basic arc in "r" is not feasible;
+ //in case adjust its flow
+ AdjustFlow( r );
+ }
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::AdjustFlow( nodePType *r )
+{
+ arcPType *arc = r->enteringTArc;
+ 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;
+ arc->tail = arc->head;
+ arc->head = temp;
+ arc->flow = -arc->flow;
+ }
+ else { // If entering arc of "r" is not a dummy arc
+ bool orientationDown = ( arc->head == r );
+ FNumber delta = 0;
+ if( LTZ( arc->flow , EpsFlw ) ) { // If flow is < 0
+ delta = -arc->flow;
+ arc->flow = 0;
+ #if( ! QUADRATICCOST )
+ arc->ident = AT_LOWER;
+ #endif
+ }
+
+ if( GT( arc->flow , arc->upper , EpsFlw ) ) {
+ // If flow goes over the capacity of the arc
+ delta = arc->upper - arc->flow;
+ arc->flow = arc->upper;
+ #if( ! QUADRATICCOST )
+ arc->ident = AT_UPPER;
+ #endif
+ }
+
+ /* This arc goes out from the basis, and the relative dummy arc goes in T.
+ Then the algorithm push flow in the cycle made by the arc and some arcs
+ of T to balance the flow. */
+
+ if( ! ETZ( delta , EpsFlw ) ) {
+ nodePType *node = Father( r , arc );
+ while( node != dummyRootP ) {
+ arc = node->enteringTArc;
+ if( ( arc->head == node ) == orientationDown )
+ arc->flow += delta;
+ else
+ arc->flow -= delta;
+
+ node = Father( node , arc );
+ }
+
+ arcPType *dummy = dummyArcsP + ( r - nodesP );
+ #if( ! QUADRATICCOST )
+ dummy->ident = BASIC;
+ #endif
+
+ /* Update the structure of the tree. If entering basic arc of "r" is
+ changed, subtree of "r"is moved next dummy root. */
+ r->enteringTArc = dummy;
+ int deltaLevel = 1 - r->subTreeLevel;
+ nodePType *lastNode = CutAndUpdateSubtree( r , deltaLevel );
+ PasteSubtree( r , lastNode , dummyRootP );
+ if( ( dummy->head == r ) != orientationDown )
+ dummy->flow += delta;
+ else
+ dummy->flow -= delta;
+
+ if( LTZ( dummy->flow , EpsFlw ) ) {
+ nodePType *temp = dummy->tail;
+ dummy->tail = dummy->head;
+ dummy->head = temp;
+ dummy->flow = -dummy->flow;
+ }
+ }
+ }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateInitialDModifiedBalanceVector( void )
+{
+ #if( ! QUADRATICCOST )
+ int i = 0;
+ modifiedBalance = new FNumber[ n ];
+ // Initialited every node's modifiedBalance to his balance
+ for( nodeDType *node = nodesD ; node != stopNodesD ; node++ ) {
+ modifiedBalance[ i ] = node->balance;
+ i++;
+ }
+
+ // Modify the vector according to the arcs out of base with flow non zero
+ // Scan the real arcs
+ for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ )
+ if( arc->ident == AT_UPPER ) {
+ i = (arc->tail) - nodesD;
+ modifiedBalance[ i ] += arc->upper;
+ i = (arc->head) - nodesD;
+ modifiedBalance[ i ] -= arc->upper;
+ }
+
+ // Scan the dummy arcs
+ for( arcDType *arc = dummyArcsD ; arc != stopDummyD ; arc++ )
+ if( arc->ident == AT_UPPER ) {
+ i = (arc->tail) - nodesD;
+ modifiedBalance[ i ] += arc->upper;
+ i = (arc->head) - nodesD;
+ modifiedBalance[ i ] -= arc->upper;
+ }
+ #endif
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PostDVisit( nodeDType *r )
+{
+ #if( ! QUADRATICCOST )
+ // The method controls if "r" is a leaf in T
+ bool rLeaf = false;
+ int i = r - nodesD;
+ 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"
+ ( 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
+ nodeDType *desc = r->nextInT;
+ // Call PostDVisit for every child of "r"
+ while( ( desc ) && ( desc->subTreeLevel > r->subTreeLevel ) ) {
+ if( desc->subTreeLevel -1 == r->subTreeLevel ) { // desc is a son of r
+ PostDVisit( desc );
+
+ if( ( desc->enteringTArc )->head == r ) // enteringTArc of desc goes in r
+ modifiedBalance[ i ] -= ( desc->enteringTArc )->flow;
+ else // If enteringTArc of "desc" goes out "r"
+ modifiedBalance[ i ] += ( desc->enteringTArc )->flow;
+ }
+
+ desc = desc->nextInT;
+ }
+
+ 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
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::ResetWhenInT2( void )
+{
+ for( nodeDType *n = nodesD ; n != stopNodesD ; n++)
+ n->whenInT2 = 0;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N, class A>
+inline N* MCFSimplex::Father( N *n , A *a )
+{
+ if( a == NULL )
+ return NULL;
+
+ if( a->tail == n )
+ return( a->head );
+ else
+ return( a->tail );
+ }
+
+/*-------------------------------------------------------------------------*/
+
+inline MCFSimplex::FONumber MCFSimplex::GetFO( void )
+{
+ FONumber fo = 0;
+ if( usePrimalSimplex ) {
+ arcPType *arco;
+ for( arco = arcsP ; arco != stopArcsP ; arco++ ) {
+ #if( QUADRATICCOST )
+ if( ! ETZ( arco->flow , EpsFlw ) )
+ fo += arco->flow * ( arco->cost + arco->flow * arco->quadraticCost / 2 );
+ #else
+ if( ( arco->ident == BASIC ) || ( arco->ident == AT_UPPER ) )
+ fo += arco->cost * arco->flow;
+ #endif
+ }
+
+ for( arco = dummyArcsP ; arco != stopDummyP ; arco++ ) {
+ #if( QUADRATICCOST )
+ if( ! ETZ( arco->flow , EpsFlw ) )
+ fo += arco->flow * ( arco->cost + arco->flow * arco->quadraticCost / 2 );
+ #else
+ if( ( arco->ident == BASIC ) || ( arco->ident == AT_UPPER ) )
+ fo += arco->cost * arco->flow;
+ #endif
+ }
+ }
+ else {
+ arcDType *a;
+ for( a = arcsD ; a != stopArcsD ; a++ ) {
+ #if (QUADRATICCOST)
+ fo += ( a->cost * a->flow ) + ( a->quadraticCost * a->flow * a->flow ) / 2;
+ #else
+ if( ( a->ident == BASIC ) || (a->ident == AT_UPPER ) )
+ fo += a->cost * a->flow;
+ #endif
+ }
+
+ for( a = dummyArcsD ; a != stopDummyD ; a++) {
+ #if (QUADRATICCOST)
+ fo += ( a->cost * a->flow ) + ( a->quadraticCost * a->flow * a->flow ) / 2;
+ #else
+ if( ( a->ident == BASIC ) || ( a->ident == AT_UPPER ) )
+ fo += a->cost * a->flow;
+ #endif
+ }
+ }
+
+ return( fo );
+ }
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::PrintPNode( nodePType *nodo )
+{
+ if( nodo )
+ if( nodo != dummyRootP )
+ cout << ( nodo - nodesP + 1 );
+ else
+ cout << "r";
+ else
+ cout << "..";
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PrintPArc( arcPType *arc )
+{
+ if( arc ) {
+ cout << "(";
+ PrintPNode( arc->tail );
+ cout << ", ";
+ PrintPNode( arc->head );
+ cout << ")";
+ }
+ else
+ cout << "..";
+}
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PrintDNode( nodeDType *nodo )
+{
+ if( nodo )
+ if( nodo != dummyRootD )
+ cout << ( nodo - nodesD + 1 );
+ else
+ cout << "r";
+ else
+ cout << "..";
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PrintDArc( arcDType *arc )
+{
+ if( arc ) {
+ cout << "(";
+ PrintDNode( arc->tail );
+ cout << ", ";
+ PrintDNode( arc->head );
+ cout << ")";
+ }
+ else
+ cout << "..";
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::nodePType* MCFSimplex::RecoverPNode( Index ind )
+{
+ if( ( ind < 0 ) || ( ind > n ) )
+ return( NULL );
+ if( ind )
+ return( nodesP + ind - 1 );
+ else
+ return( dummyRootP );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcPType* MCFSimplex::RecoverPArc( nodePType *tail ,
+ nodePType *head )
+{
+ if( ( tail == NULL ) || ( head == NULL ) )
+ return( NULL );
+
+ arcPType *arc = arcsP;
+ while( ( arc->tail != tail ) || ( arc->head != head ) ) {
+ arc++;
+ if( arc == stopArcsP )
+ arc = dummyArcsP;
+ if( arc == stopDummyP )
+ return( NULL );
+ }
+
+ return( arc );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::nodeDType* MCFSimplex::RecoverDNode( Index ind )
+{
+ if( ( ind < 0 ) || ( ind > n ) )
+ return( NULL );
+
+ if( ind )
+ return( nodesD + ind - 1 );
+ else
+ return( dummyRootD );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcDType* MCFSimplex::RecoverDArc( nodeDType *tail ,
+ nodeDType *head )
+{
+ if( ( tail == NULL ) || ( head == NULL ) )
+ return( NULL );
+
+ arcDType *arc = arcsD;
+ while( ( arc->tail != tail ) || ( arc->head != head ) ) {
+ arc++;
+ if( arc == stopArcsD )
+ arc = dummyArcsD;
+ if( arc == stopDummyD )
+ return( NULL );
+ }
+
+ return( arc );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::infoPNode( nodePType *node , int tab )
+{
+ for( int t = 0 ; t < tab ; t++ )
+ cout << "\t";
+ cout << "Nodo ";
+ PrintPNode( node );
+ cout << ": b = " << node->balance << " y = " << node->potential << endl;
+ #if( UNIPI_VIS_NODE_BASIC_ARC )
+ cout << ": TArc=";
+ PrintPArc( node->enteringTArc );
+ cout << endl;
+ #endif
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::infoPArc( arcPType *arc , int ind , int tab )
+{
+ for( int t = 0 ; t < tab ; t++ )
+ cout << "\t";
+ cout << "Arco ";
+ PrintPArc( arc );
+ cout << ": x = " << arc->flow;
+ #if( UNIPI_VIS_ARC_UPPER )
+ cout << " u = " << arc->upper;
+ #endif
+ #if( UNIPI_VIS_ARC_COST )
+ cout << " c = " << arc->cost;
+ #endif
+ #if( QUADRATICCOST )
+ #if( UNIPI_VIS_ARC_Q_COST )
+ cout << " q = " << arc->quadraticCost;
+ #endif
+ cout << endl;
+ for( int t = 0 ; t < tab ; t++ )
+ cout << "\t";
+ #if( UNIPI_VIS_ARC_REDUCT_COST )
+ cout << " rc = " << MCFGetRC( ind );
+ #endif
+ #else
+ cout << endl;
+ for( int t = 0 ; t < tab ; t++ )
+ cout << "\t";
+ #if( UNIPI_VIS_ARC_REDUCT_COST )
+ cout << " rc = " << MCFGetRC( ind );
+ #endif
+ #if( UNIPI_VIS_ARC_STATE )
+ switch( arc->ident ) {
+ case( BASIC ): cout << " in T"; break;
+ case( AT_LOWER ): cout << " in L"; break;
+ case( AT_UPPER ): cout << " in U"; break;
+ case( DELETED ): cout << " canceled"; break;
+ case( CLOSED ): cout << " closed";
+ }
+ #endif
+ #endif
+ cout << endl;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::infoDNode( nodeDType *node , int tab )
+{
+ for( int t = 0 ; t < tab; t++ )
+ cout << "\t";
+ cout << "Nodo ";
+ PrintDNode( node );
+ cout << ": b = " << node->balance << " y = " << node->potential;
+ #if( UNIPI_VIS_NODE_BASIC_ARC )
+ cout << ": TArc=";
+ PrintDArc( node->enteringTArc );
+ cout << endl;
+ #endif
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::infoDArc( arcDType *arc , int ind , int tab )
+{
+ for( int t = 0 ; t < tab ; t++ )
+ cout << "\t";
+ cout << "Arco ";
+ PrintDArc( arc );
+ cout << " x = " << arc->flow;
+ #if( UNIPI_VIS_ARC_UPPER )
+ cout << " u = " << arc->upper;
+ #endif
+ #if( UNIPI_VIS_ARC_COST )
+ cout << " c = " << arc->cost;
+ #endif
+ #if( QUADRATICCOST )
+ #if( UNIPI_VIS_ARC_Q_COST )
+ cout << " q = " << arc->quadraticCost;
+ #endif
+ cout << endl;
+ for( int t = 0 ; t < tab ; t++ )
+ cout << "\t";
+ #if( UNIPI_VIS_ARC_REDUCT_COST )
+ cout << " rc = " << MCFGetRC( ind );
+ #endif
+ #else
+ cout << endl;
+ for( int t = 0 ; t < tab ; t++ )
+ cout << "\t";
+ #if( UNIPI_VIS_ARC_REDUCT_COST )
+ cout << " rc = " << MCFGetRC( ind );
+ #endif
+ #if (UNIPI_VIS_ARC_STATE)
+ switch( arc->ident ) {
+ case( BASIC ): cout << " in T"; break;
+ case( AT_LOWER ): cout << " in L"; break;
+ case( AT_UPPER ): cout << " in U"; break;
+ case( DELETED ): cout << " canceled"; break;
+ case( CLOSED ): cout << " closed";
+ }
+ #endif
+ #endif
+ cout << endl;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::ShowSituation( int tab )
+{
+ if( usePrimalSimplex ) {
+ arcPType *arc;
+ nodePType *node;
+ int i = 0;
+ for( arc = arcsP ; arc != stopArcsP ; arc++ ) {
+ infoPArc( arc , i , tab );
+ i++;
+ }
+ cout << endl;
+ #if( UNIPI_VIS_DUMMY_ARCS )
+ i = 0;
+ for( arc = dummyArcsP ; arc != stopDummyP ; arc++ ) {
+ infoPArc( arc , i , tab );
+ i++;
+ }
+ cout << endl;
+ #endif
+ infoPNode( dummyRootP , tab );
+ for( node = nodesP ; node != stopNodesP ; node++ )
+ infoPNode( node , tab );
+ }
+ else {
+ arcDType *arc;
+ nodeDType *node;
+ int i = 0;
+ for( arc = arcsD ; arc != stopArcsD ; arc++ ) {
+ infoDArc( arc , i , tab );
+ i++;
+ }
+ cout << endl;
+ #if( UNIPI_VIS_DUMMY_ARCS )
+ i = 0;
+ for( arc = dummyArcsD ; arc != stopDummyD ; arc++) {
+ infoDArc( arc , i , tab );
+ i++;
+ }
+ cout << endl;
+ #endif
+ infoDNode( dummyRootD , tab );
+ for( node = nodesD ; node != stopNodesD ; node++ )
+ infoDNode( node , tab );
+ }
+ }
+
+/*-------------------------------------------------------------------------*/
+/*---------------------- End File MCFSimplex.C ----------------------------*/
+/*-------------------------------------------------------------------------*/
diff --git a/clang/lib/Analysis/MCFSimplex.h b/clang/lib/Analysis/MCFSimplex.h
new file mode 100644
index 0000000..a8d84f6
--- /dev/null
+++ b/clang/lib/Analysis/MCFSimplex.h
@@ -0,0 +1,1086 @@
+/*--------------------------------------------------------------------------*/
+/*------------------------- File MCFSimplex.h ------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @file
+ * Linear and Quadratic Min Cost Flow problems solver based on the (primal and
+ * dual) simplex algorithm. Conforms to the standard MCF interface defined in
+ * MCFClass.h.
+ *
+ * \Version 1.00
+ *
+ * \date 29 - 08 - 2011
+ *
+ * \author Alessandro Bertolini \n
+ * Antonio Frangioni \n
+ * Operations Research Group \n
+ * Dipartimento di Informatica \n
+ * Universita' di Pisa \n
+ *
+ * Copyright &copy 2008 - 2011 by Alessandro Bertolini, Antonio Frangioni
+ */
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*----------------------------- DEFINITIONS --------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#ifndef __MCFSimplex
+ #define __MCFSimplex /* self-identification: #endif at the end of the file */
+
+/*--------------------------------------------------------------------------*/
+/*------------------------------ INCLUDES ----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#include "MCFClass.h"
+
+/*@} -----------------------------------------------------------------------*/
+/*--------------------------- NAMESPACE ------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#if( OPT_USE_NAMESPACES )
+namespace MCFClass_di_unipi_it
+{
+#endif
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------------- MACROS ----------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @defgroup MCFSimplex_MACROS Compile-time switches in MCFSimplex.h
+ There is only one macro in MCFSimplex, but it is very important!
+ @{ */
+
+#define QUADRATICCOST 0
+
+/**< Setting QUADRATICCOST == 1 the solver can solve problems with linear and
+ quadratic costs too (but the latter only with the Primal Simplex).
+ The reason for having a macro is that when quadratic costs are present the
+ "arcType" struct has the additional field "quadraticCost" to hold it.
+ Furthermore, the field "ident" is not created because the solver doesn't
+ use the classical TLU tripartition. Instead, closed arcs and deleted arcs
+ are characterized as follows:
+ - closed arcs have the field "cost" to INFINITY (Inf<FNumber>());
+ - deleted arcs have the field "upper" to INFINITY and the "tail" and
+ "head" field are NULL.
+ Furthermore, the solver needs the variables "ignoredEnteringArc" and
+ "firstIgnoredEnteringArc", used to avoid nasty loops during the execution
+ of the Quadratic Primal Simplex algorithm.
+ If, instead, QUADRATICCOST == 0 then the solver can solve only problems
+ with linear costs. Hence, the field "quadraticCost" is useless and it
+ isn't created. Furthermore, Primal Simplex and Dual Simplex use the
+ tripartition TLU to divide the arcs, so the solver creates the field
+ "ident", which differentiates the set of the arcs in among deleted arcs,
+ closed arcs, arcs in T, arcs in L, arcs in U.
+ Thus, with QUADRATICCOST == 0 the solver cannot solve problems with
+ quadratic costs, but it does solve problems with linear costs faster. */
+
+/*@} end( group( MCFCLASS_MACROS ) ) */
+
+/*--------------------------------------------------------------------------*/
+/*---------------------------- CLASSES -------------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @defgroup MCFSimplex_CLASSES Classes in MCFSimplex.h
+ @{ */
+
+/** The MCFSimplex class derives from the abstract base class MCFClass, thus
+ sharing its (standard) interface, and implements both the Primal and
+ Dual network simplex algorithms for solving (Linear and Quadratic)
+ Min Cost Flow problems */
+
+class MCFSimplex: public MCFClass
+{
+
+/*--------------------------------------------------------------------------*/
+/*----------------------- PUBLIC PART OF THE CLASS -------------------------*/
+/*--------------------------------------------------------------------------*/
+/*-- --*/
+/*-- The following methods and data are the actual interface of the --*/
+/*-- class: the standard user should use these methods and data only. --*/
+/*-- --*/
+/*--------------------------------------------------------------------------*/
+
+ public:
+
+/*--------------------------------------------------------------------------*/
+/*---------------------------- PUBLIC TYPES --------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+/** Public enum describing the parameters of MCFSimplex. */
+
+ enum SimplexParam
+ {
+ kAlgPrimal = kLastParam , ///< parameter to set algorithm (Primal/Dual):
+ kAlgPricing , ///< parameter to set algorithm of pricing
+ kNumCandList , /**< parameter to set the number of candidate
+ list for Candidate List Pivot method */
+ kHotListSize , /**< parameter to set the size of Hot List
+ for Candidate List Pivot method */
+ kRecomputeFOLimits , /**< parameter to set the number of iterations
+ in which quadratic Primal Simplex computes
+ "manually" the f.o. value */
+ kEpsOpt /**< parameter to set the precision of the
+ objective function value for the
+ quadratic Primal Simplex */
+ };
+
+/** Public enum describing the pricing rules in MCFSimplex::SetAlg(). */
+
+ enum enumPrcngRl
+ {
+ kDantzig = 0, ///< Dantzig's rule (most violated constraint)
+ kFirstEligibleArc , ///< First eligible arc in round-robin
+ kCandidateListPivot ///< Candidate List Pivot Rule
+ };
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- PUBLIC METHODS -------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*---------------------------- CONSTRUCTOR ---------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+ MCFSimplex( cIndex nmx = 0 , cIndex mmx = 0 );
+
+/**< Constructor of the class, as in MCFClass::MCFClass(). */
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------- OTHER INITIALIZATIONS -------------------------*/
+/*--------------------------------------------------------------------------*/
+
+ void LoadNet( cIndex nmx = 0 , cIndex mmx = 0 , cIndex pn = 0 ,
+ cIndex pm = 0 , cFRow pU = NULL , cCRow pC = NULL ,
+ cFRow pDfct = NULL , cIndex_Set pSn = NULL ,
+ cIndex_Set pEn = NULL );
+
+/**< Inputs a new network, as in MCFClass::LoadNet(). */
+
+/*--------------------------------------------------------------------------*/
+
+ void SetAlg( bool UsPrml , char WhchPrc );
+
+/**< Selects which algorithm (Primal vs Dual Network Simplex), and which
+ pricing rule within the algorithm, is used.
+
+ If UsPrml == TRUE then the Primal Network Simplex algorithm is used,
+ otherwise the Dual Network Simplex is used.
+
+ The allowed values for WhchPrc are:
+
+ - kDantzig Dantzig's pricing rule, i.e., most violated dual
+ constraint; this can only be used with the Primal
+ Network Simplex
+
+ - kFirstEligibleArcA "dumb" rule, first eligible arc in round-robin;
+
+ - kCandidateListPivot Candidate List Pivot Rule
+
+ If this method is *not* called, the Primal Network Simplex with the
+ Candidate List Pivot Rule (the best setting on most instances) is
+ used. */
+
+/*--------------------------------------------------------------------------*/
+
+ void SetPar( int par , int val );
+
+/**< Set general integer parameters.
+
+ @param par is the parameter to set; since this method accepts an int
+ value, the enum SimplexParam can be used in addition to the
+ enum MCFParam to specify the integer parameters (every enum
+ is an int).
+
+ @param value is the value to assign to the parameter. */
+
+/*-------------------------------------------------------------------------*/
+
+ void SetPar( int par , double val );
+
+/**< Set general float parameters.
+
+ @param par is the parameter to set; since this method accepts an int
+ value, the enum SimplexParam can be used in addition to the
+ enum MCFParam to specify the float parameters (every enum
+ is an int).
+
+ @param value is the value to assign to the parameter. */
+
+/*--------------------------------------------------------------------------*/
+/*-------------------- METHODS FOR SOLVING THE PROBLEM ---------------------*/
+/*--------------------------------------------------------------------------*/
+
+ void SolveMCF( void );
+
+/*--------------------------------------------------------------------------*/
+/*---------------------- METHODS FOR READING RESULTS -----------------------*/
+/*--------------------------------------------------------------------------*/
+
+ void MCFGetX( FRow F , Index_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+/*--------------------------------------------------------------------------*/
+
+ void MCFGetRC( CRow CR , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() ) ;
+
+ CNumber MCFGetRC( cIndex i );
+
+/*--------------------------------------------------------------------------*/
+
+ void MCFGetPi( CRow P , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+/*--------------------------------------------------------------------------*/
+
+ FONumber MCFGetFO( void );
+
+/*--------------------------------------------------------------------------*/
+/*-------------- METHODS FOR READING THE DATA OF THE PROBLEM ---------------*/
+/*--------------------------------------------------------------------------*/
+
+ virtual void MCFArcs( Index_Set Startv , Index_Set Endv ,
+ cIndex_Set nms = NULL , cIndex strt = 0 ,
+ Index stp = Inf<Index>() );
+
+ inline Index MCFSNde( cIndex i );
+
+ inline Index MCFENde( cIndex i );
+
+/*--------------------------------------------------------------------------*/
+
+ void MCFCosts( CRow Costv , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+ inline CNumber MCFCost( cIndex i );
+
+/*--------------------------------------------------------------------------*/
+
+ void MCFQCoef( CRow Qv , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+ inline CNumber MCFQCoef( cIndex i );
+
+/*--------------------------------------------------------------------------*/
+
+ void MCFUCaps( FRow UCapv , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+ inline FNumber MCFUCap( cIndex i );
+
+/*--------------------------------------------------------------------------*/
+
+ void MCFDfcts( FRow Dfctv , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+ inline FNumber MCFDfct( cIndex i );
+
+/*--------------------------------------------------------------------------*/
+/*------------- METHODS FOR ADDING / REMOVING / CHANGING DATA --------------*/
+/*--------------------------------------------------------------------------*/
+/*------- Changing the costs, QCoef, deficits and upper capacities ---------*/
+/*--------------------------------------------------------------------------*/
+
+ void ChgCosts( cCRow NCost , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+ void ChgCost( Index arc , cCNumber NCost );
+
+/*--------------------------------------------------------------------------*/
+
+ void ChgQCoef( cCRow NQCoef = NULL , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+ void ChgQCoef( Index arc , cCNumber NQCoef );
+
+/*--------------------------------------------------------------------------*/
+
+ void ChgDfcts( cFRow NDfct , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+ void ChgDfct( Index nod , cFNumber NDfct );
+
+/*--------------------------------------------------------------------------*/
+
+ void ChgUCaps( cFRow NCap , cIndex_Set nms = NULL ,
+ cIndex strt = 0 , Index stp = Inf<Index>() );
+
+ void ChgUCap( Index arc , cFNumber NCap );
+
+/*--------------------------------------------------------------------------*/
+/*--------------- Modifying the structure of the graph ---------------------*/
+/*--------------------------------------------------------------------------*/
+
+ void CloseArc( cIndex name );
+
+ void DelNode( cIndex name );
+
+ bool IsClosedArc( cIndex name );
+
+ void OpenArc( cIndex name ) ;
+
+ Index AddNode( cFNumber aDfct );
+
+ void ChangeArc( cIndex name , cIndex nSS = Inf<Index>() ,
+ cIndex nEN = Inf<Index>() );
+
+ void DelArc( cIndex name );
+
+ Index AddArc( cIndex Start , cIndex End , cFNumber aU , cCNumber aC );
+
+ bool IsDeletedArc( cIndex name );
+
+/*--------------------------------------------------------------------------*/
+/*------------------------------ DESTRUCTOR --------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+ ~MCFSimplex();
+
+/*--------------------------------------------------------------------------*/
+/*--------------------- PRIVATE PART OF THE CLASS --------------------------*/
+/*--------------------------------------------------------------------------*/
+/*-- --*/
+/*-- Nobody should ever look at this part: everything that is under this --*/
+/*-- advice may be changed without notice in any new release of the code. --*/
+/*-- --*/
+/*--------------------------------------------------------------------------*/
+
+ private:
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------- PRIVATE DATA TYPES ----------------------------*/
+/*--------------------------------------------------------------------------*/
+/*-- --*/
+/*-- Let T \subseteq A be a spanning tree, and consider some node v \in V --*/
+/*-- \setminus { 0 }. There is an unique (undirected) path, denoted by --*/
+/*-- P(v), defined by T from v to the root node 0. The arc in P(v), which --*/
+/*-- is incident to v, is called the *basic arc* of v. The other terminal --*/
+/*-- node u of this basic arc is called the *father node* of v. The --*/
+/*-- spanning tree T is represented saving the basic arc of every node, --*/
+/*-- and maintaining the order of the nodes and the depth as to T root --*/
+/*-- after a Post-Visit of T. This order is saved in a bidirectional list --*/
+/*-- written in the node. --*/
+/*-- --*/
+/*-- The Primal Simplex uses a different data structure than the Dual --*/
+/*-- Simplex, because the Dual Simplex needs additional data (mainly the --*/
+/*-- Backward Star and Forward Star. Furthermore, the Primal Simplex uses --*/
+/*-- different data structures in the quadratic case. --*/
+/*-- --*/
+/*--------------------------------------------------------------------------*/
+
+ struct arcPType; // pre-declaration of the arc structure (pointers to arcs
+ // are contained in the node structure) for Primal Simplex
+
+ struct arcDType; // pre-declaration of the arc structure (pointers to arcs
+ // are contained in the node structure) for Dual Simplex
+
+ typedef double iteratorType; // type for the iteration counter and array
+ // "whenInT2"
+
+ struct nodePType { // node structure for Primal Simplex - - - - - - - -
+ nodePType *prevInT; // previous node in the order of the Post-Visit on T
+
+ nodePType *nextInT; // next node in the order of the Post-Visit on T
+
+ arcPType *enteringTArc; // entering basic arc of this node
+
+ FNumber balance; // supply/demand of this node; a node is called a
+ // supply node, a demand node, or a transshipment
+ // node depending upon whether balance is larger
+ // than, smaller than, or equal to zero
+ #if( QUADRATICCOST )
+ CNumber sumQuadratic; // the sum of the quadratic coefficients of the tree's arcs
+ // from root of T to the node
+
+ FONumber potential; // the node potential corresponding with the flow
+ // conservation constrait of this node
+ #else
+ CNumber potential; // the node potential corresponding with the flow
+ // conservation constrait of this node
+ #endif
+
+ int subTreeLevel; // the depth of the node in T as to T root
+
+ }; // end( struct( nodePType ) )
+
+ struct nodeDType { // node structure for Dual Simplex - - - - - - - - -
+ nodeDType *prevInT; // previous node in the order of the Post-Visit on T
+
+ nodeDType *nextInT; // next node in the order of the Post-Visit on T
+
+ arcDType *enteringTArc; // entering basic arc of this node
+
+ FNumber balance; // supply/demand of this node; a node is called a
+ // supply node, a demand node, or a transshipment
+ // node depending upon whether balance is larger
+ // than, smaller than, or equal to zero
+
+ #if( QUADRATICCOST )
+ CNumber sumQuadratic; // the sum of the quadratic coefficients of the tree's arcs
+ // from root of T to the node
+
+ FONumber potential; // the node potential corresponding with the flow
+ // conservation constrait of this node
+ #else
+ CNumber potential; // the node potential corresponding with the flow
+ // conservation constrait of this node
+ #endif
+
+ int subTreeLevel; // the depth of the node in T as to T root
+ iteratorType whenInT2; // the last iteration where a node is in subtree T2
+
+ Index numArcs; // the number of the arcs which enter/exit from node
+ arcDType *firstBs; // the first arc in the node's Backward Star
+ arcDType *firstFs; // the first arc in the node's Forward Star
+
+ }; // end( struct( nodeDType ) )
+
+ struct arcPType { // arc structure for Primal Simplex - - - - - - - -
+ nodePType *tail; // tail node
+ nodePType *head; // head node
+ FNumber flow; // arc flow
+ CNumber cost; // arc linear cost
+
+ #if( QUADRATICCOST )
+ CNumber quadraticCost; // arc quadratic cost
+ #else
+ char ident; // tells if arc is deleted, closed, in T, L, or U
+ #endif
+
+ FNumber upper; // arc upper bound
+ }; // end( struct( arcPType ) )
+
+ struct arcDType { // arc structure for Primal Simplex - - - - - - - -
+ nodeDType *tail; // tail node
+ nodeDType *head; // head node
+ FNumber flow; // arc flow
+ CNumber cost; // arc linear cost
+
+ #if( QUADRATICCOST )
+ CNumber quadraticCost; // arc quadratic cost
+ #else
+ char ident; // indicates if arc is deleted, closed, in T, in L, or in U
+ #endif
+
+ FNumber upper; // arc upper bound
+ arcDType *nextBs; // the next arc in the Backward Star of the arc's head
+ arcDType *nextFs; // the next arc in the Forward Star of the arc's tail
+
+ }; // end( struct( arcDType ) )
+
+ struct primalCandidType { // Primal Candidate List- - - - - - - - - - - - -
+ arcPType *arc; // pointer to the violating primal bound arc
+
+ #if(QUADRATICCOST )
+ FONumber absRC; // absolute value of the arc's reduced cost
+ #else
+ CNumber absRC; // absolute value of the arc's reduced cost
+ #endif
+ }; // end( struct( primalCandidateType ) )
+
+ struct dualCandidType { // Dual Candidate List- - - - - - - - - - - - - -
+ nodeDType *node; // deepest node violating the dual bound arc
+ FNumber absInfeas; // absolute value of the arc's flow infeasibility
+ }; // end( struct( dualCandidateType ) )
+
+/*--------------------------------------------------------------------------*/
+/*----------------------- PRIVATE DATA STRUCTURES -------------------------*/
+/*--------------------------------------------------------------------------*/
+
+ bool usePrimalSimplex; // TRUE if the Primal Network Simplex is used
+
+ char pricingRule; // which pricing rule is used
+
+ nodePType *nodesP; // vector of nodes: points to the n + 1 node
+ // structs (including the dummy root node)
+ // where the first node is indexed by zero
+ // and the last node is the dummy root node
+
+ nodePType *dummyRootP; // the dummy root node
+
+ nodePType *stopNodesP; // first infeasible node address = nodes + n
+
+ arcPType *arcsP; // vector of arcs; this variable points to
+ // the m arc structs.
+
+ arcPType *dummyArcsP; // vector of artificial dummy arcs: points
+ // to the artificial dummy arc variables and
+ // contains n arc structs. The artificial
+ // arcs are used to build artificial feasible
+ // starting bases. For each node i, there is
+ // exactly one dummy arc defined to connect
+ // the node i with the dummy root node.
+
+ arcPType *stopArcsP; // first infeasible arc address = arcs + m
+
+ arcPType *stopDummyP; // first infeasible dummy arc address
+ // = arcs + m + n
+
+ arcPType *arcToStartP; // Dantzig Rule and First Eligible Arc Rule
+ // start their search from this arc
+
+ nodeDType *nodesD; // vector of nodes: points to the n + 1 node
+ // structs (including the dummy root node)
+ // where the first node is indexed by zero
+ // and the last node is the dummy root node
+
+ nodeDType *dummyRootD; // the dummy root node
+
+ nodeDType *stopNodesD; // first infeasible node address = nodes + n
+
+ arcDType *arcsD; // vector of arcs; this variable points to
+ // the m arc structs.
+
+ arcDType *dummyArcsD; // vector of artificial dummy arcs: points to
+ // to the artificial dummy arc variables and
+ // contains n arc structs. The artificial
+ // arcs are used to build artificial feasible
+ // starting bases. For each node i, there is
+ // exactly one dummy arc defined to connect
+ // the node i with the dummy root node.
+
+ arcDType *stopArcsD; // first infeasible arc address = arcs + m
+
+ arcDType *stopDummyD; // first infeasible dummy arc address
+ // = arcs + m + n
+
+ arcDType *arcToStartD; // Dantzig Rule and First Eligible Arc Rule
+ // start their search from this arc
+
+ iteratorType iterator; // the current number of iterations
+
+ primalCandidType *candP; // every element points to an element of the
+ // arcs vector which contains an arc violating
+ // dual bound
+
+ dualCandidType *candD; // every element points to an element of the
+ // arcs vector which contains an arc violating
+ // primal bond
+
+ Index numGroup; // number of the candidate lists
+
+ Index tempCandidateListSize; // hot list dimension (it is variable)
+
+ Index groupPos; // contains the actual candidate list
+
+ Index numCandidateList; // number of candidate lists
+
+ Index hotListSize; // number of candidate lists and hot list dimension
+
+ Index forcedNumCandidateList; // value used to force the number of candidate list
+
+ Index forcedHotListSize; // value used to force the number of candidate list
+ // and hot list dimension
+
+ bool newSession; // true if algorithm is just started
+
+ CNumber MAX_ART_COST; // large cost for artificial arcs
+
+ FNumber *modifiedBalance; // vector of balance used by the PostVisit
+
+ FONumber EpsOpt; // the precision of the objective function value
+ // for the quadratic case of the Primal Simplex
+
+ int recomputeFOLimits; // after how many iterations the quadratic Primal
+ // Simplex computes "manually" the o.f. value
+
+ #if( QUADRATICCOST )
+ FONumber foValue; // the temporary objective function value
+ #endif
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------- PRIVATE METHODS -------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+ void MemAlloc( void );
+
+/**< Method to allocate memory for the main data structures. It creates the
+ vector of nodes, the vector of arcs (real and dummy). If the algorithm is
+ using the Dual Simplex, it creates also the vectors whenInT2, firstIn and
+ nextIn, usefull to identify the next entering. */
+
+/*--------------------------------------------------------------------------*/
+
+ void MemDeAlloc( bool whatDeAlloc );
+
+/**< Method to deallocate memory for the main data structures created in
+ MemAlloc(). */
+
+/*--------------------------------------------------------------------------*/
+
+ void MemAllocCandidateList( void );
+
+/**< Method to allocate memory for the data structure used by the Candidate
+ List Pivot Rule. It creates the vector of candP (or candD in the Dual
+ Simplex case), determining the size of this vector on the basis of number
+ of arcs or nodes. */
+
+/*--------------------------------------------------------------------------*/
+
+ void MemDeAllocCandidateList( void );
+
+/**< Method to deallocate memory for the data structures created in
+ MemAllocCandidateList(). */
+
+/*--------------------------------------------------------------------------*/
+
+ void CreateInitialPrimalBase( void );
+
+/**< Method to create an initial feasible primal base. Add one node (dummyRoot)
+ to the network and connect this dummy root to each real node with n dummy
+ arcs. For each source node i, a dummy arc (i, r) exists. For each sink or
+ transit node j, a dummy arc (r, j) exists. The source nodes send their flows
+ to the dummy root through their dummy arcs, so the dummy root send all to
+ the sink nodes. The dummy root balance is 0, the costs of dummy arcs are
+ fixed to "infinity". */
+
+/*--------------------------------------------------------------------------*/
+
+ void CreateInitialDualBase( void );
+
+/**< Method to create an initial feasible dual base. Add one node (dummyRoot)
+ to the network and connect this dummy root to each real node with n dummy
+ arcs. The source nodes send their flows to the dummy root through their
+ dummy arcs, so the dummy root send all to the sink nodes. The dummy root
+ balance is 0, the costs of dummy arcs are fixed to "infinity". */
+
+/*--------------------------------------------------------------------------*/
+
+ void CreateAdditionalDualStructures( void );
+
+/**< The Dual Simplex needs nodes' Backward and Forward Star to work. So when
+ the Primal Simplex runs, these structures don't exist. When the Dual
+ Simplex starts, these structure are created in this method. */
+
+/*--------------------------------------------------------------------------*/
+
+ void PrimalSimplex( void );
+
+/**< Main method to implement the Primal Simplex algorithm. */
+
+/*--------------------------------------------------------------------------*/
+
+ void DualSimplex( void );
+
+/**< Main method to implement the Dual Simplex algorithm. */
+
+/*--------------------------------------------------------------------------*/
+
+ template<class N, class A>
+ void UpdateT( A *h , A *k , N *h1 , N *h2 , N *k1 , N *k2 );
+
+ /**< Method to update the spanning tree T every iteration.
+ The spanning tree T is implemented with a bidirectional list stored
+ in the node structure, which represents the nodes' order after a
+ Post-Visit. The parameter of the method are the outgoing arc "h", the
+ incoming arc "k", and the four nodes on the extremity of these arcs
+ (for example h2 is the deepest node of the outgoing arc). Removing the
+ arc "h" splits T in two subtrees: T1 (which contains the root of T) and
+ T2, which will be re-connected by the incoming arc "k".
+ T2 will be reordered; in fact, the node "k2" becomes the root of T2
+ instead of "h2" and the hierarchy of T2 will be overturned. Then T2 will
+ be moved; the root of T2 is changed, therefore the predecessor of the
+ root will become the node "k1".
+
+ This method uses the methods cutAndUpdateSubtree() and pasteSubtree().
+ First it cuts the node "k2" and its subtree from T using
+ cutAndUpdateSubtree(). Moreover the method cutAndUpdateSubtree()
+ updates the field "subTreeLevel" of every subtree's nodes, since k2's
+ subtree will be moved from the bottom to the top of T2. Then the method
+ pasteSubtree() puts this subtree in the bidirectional list after the node
+ "k1". The same operations will be applied to the old precedessor of "k2"
+ (which will become one of the childs of "k2"). This second subtree will
+ be cut, the subTreeLevel fields will be updated, and it will be inserted
+ in the bidirectional list after the k2's subtree. This is iterated until
+ the node "h2" is reached. */
+
+/*--------------------------------------------------------------------------*/
+
+ template<class N>
+ N* CutAndUpdateSubtree( N *root, int delta );
+
+/**< This method cuts a generic subtree from the spanning tree T. Then it
+ updates the field "subTreeLevel" of every subtree's nodes adding the value
+ "delta". This method returns the last node of the subtree. */
+
+/*--------------------------------------------------------------------------*/
+
+ template<class N>
+ void PasteSubtree( N *root , N *lastNode , N *previousNode );
+
+/**< This method inserts a generic subtree with root passed by parameter into
+ the spanning tree T, between the nodes "previousNode" and "lastNode". */
+
+/*--------------------------------------------------------------------------*/
+
+ arcPType* RuleDantzig( void );
+
+/**< This method returns an arc which violates the dual conditions. It searchs
+ the arc with most violation of dual conditions in the entire set of real
+ arcs. It can be used only by the Primal Simplex in the case of networks
+ with linear costs. */
+
+/*--------------------------------------------------------------------------*/
+
+ arcPType* PRuleFirstEligibleArc( void );
+
+/**< This method returns the first found arc which violates the dual conditions
+ in the case of Primal Simplex, the primal condition in the case of Dual
+ Simplex. It can be used only in the case of networks with linear costs. */
+
+/*--------------------------------------------------------------------------*/
+
+ arcDType* DRuleFirstEligibleArc( void );
+
+/**< This method returns the first found arc which violates the dual conditions
+ in the case of Primal Simplex, the primal condition in the case of Dual
+ Simplex. It can be used only in the case of networks with linear costs. */
+
+/*--------------------------------------------------------------------------*/
+
+ arcPType* RulePrimalCandidateListPivot( void );
+
+/**< This method returns an arc which violates the dual conditions. It searches
+ the arc with most violation of dual conditions in a small set of candidate
+ arcs. In every iteration the method rebuilds this set of arcs executing three
+ phases:
+
+ - in the first phase it analyzes the remaining arcs and delete the arcs which
+ don't violate the dual condition any more;
+
+ - in the second phase it tries to fill the set, so it searchs other arcs
+ which violate the dual condition: the set of arcs is divided into "buckets"
+ which are searched sequentially until the candidate list is full; the
+ last visited bucket is retained, and the search is restarted from that
+ one at later iterations
+
+ - in the third phase the small set of candidate arcs is ordered according
+ to the violation of dual condition by the method SortPrimalCandidateList()
+ using an implementation of the algorithm "quicksort".
+
+ At last the method returns the first arc in the ordered small set. If the
+ arc doesn't exist (the set is empty), it returns NULL. */
+
+/*--------------------------------------------------------------------------*/
+
+ inline void InitializePrimalCandidateList( void );
+
+/**< Method to initialize some important details for Primal Candidate List
+ Rule. */
+
+/*--------------------------------------------------------------------------*/
+
+ inline void SortPrimalCandidateList( Index min , Index max );
+
+/**< Method to order the little set of candidate arcs according to
+ infeasibility of dual conditions. It implements the "quicksort"
+ algorithm. */
+
+/*--------------------------------------------------------------------------*/
+
+ arcDType* RuleDualCandidateListPivot( void );
+
+/**< Similar to RulePrimalCandidateListPivot() for the Dual Simplex. */
+
+/*--------------------------------------------------------------------------*/
+
+ inline void InitializeDualCandidateList( void );
+
+/**< Method to initialize some important details for Dual Candidate List Rule.
+ */
+
+/*--------------------------------------------------------------------------*/
+
+ inline void SortDualCandidateList( Index min , Index max );
+
+/**< Similar to SortPrimalCandidateList() for the Dual Simplex. */
+
+/*--------------------------------------------------------------------------*/
+
+ template<class N, class RCT>
+ inline void AddPotential( N *r , RCT delta );
+
+/**< Method to quickly update the dual solutions. During the change of the
+ base, the potential of node "k2" (deepest node in T of incoming arc "k",
+ and new root of T2) changes according to the new structure of T. In fact,
+ the precedessor of "k2" changes: now the predecessor of "k2"is "k1" (the
+ other node of the incoming arc "k"). This change of predecessor causes a
+ change of potential of "k2". The change of potential of "k2" causes the
+ changes of potential of all nodes of T2. This method computes the change
+ of potential of "k2" and applies it on all the nodes of T2. */
+
+/*--------------------------------------------------------------------------*/
+
+ template<class N>
+ inline void ComputePotential( N *r );
+
+/**< Method to update the dual solutions. It computes all the potential of the
+ nodes of the subtree which has r as root. */
+
+/*--------------------------------------------------------------------------*/
+
+ inline void ResetWhenInT2( void );
+
+/**< Method to order the small set of candidate arcs according to dual
+ infeasibility. It implements the algorithm "quicksort". */
+
+/*--------------------------------------------------------------------------*/
+
+ void CreateInitialPModifiedBalanceVector( void );
+
+/**< Method to initialize the vector of modified balance for the postvisit on
+ T in the Primal Simplex data structure. */
+
+/*--------------------------------------------------------------------------*/
+
+ void PostPVisit( nodePType *r );
+
+/**< Method to calculate the flow on the basic arcs with the Primal Simplex's
+ data structure. It uses the set of the upper bound arcs, the construction
+ of a modified balance vector and the postvisit on T. */
+
+/*--------------------------------------------------------------------------*/
+
+ void BalanceFlow( nodePType *r );
+
+/**< This method works after the method PostPVisit (in a Primal context). It
+ restores primal admissimibility on the r's subtree. It starts from the leaf
+ of the subtree and goes up to the root, using the method AdjustFlow. */
+
+/*--------------------------------------------------------------------------*/
+
+ void AdjustFlow( nodePType *r );
+
+/**< This method checks the primal admissibility of the r's basic entering arc.
+ If it is out of bounds, the method removes it from the tree (and keeps the
+ relative dummy arc) and push flow in the cycle (some tree's arc and the old
+ entering arc) to restores the right balances of the node. */
+
+/*--------------------------------------------------------------------------*/
+
+ void CreateInitialDModifiedBalanceVector( void );
+
+/**< Method to initialize the vector of modified balance for the postvisit on
+ T in the Dual Simplex data structure. */
+
+/*--------------------------------------------------------------------------*/
+
+ void PostDVisit( nodeDType *r );
+
+/**< Method to calculate the flow on the basic arcs with the Dual Simplex's data
+ structure, using the set of the upper bound arcs, the construction of a
+ modified balance vector and the postvisit on T. */
+
+/*--------------------------------------------------------------------------*/
+
+ template<class N, class A>
+ inline N* Father( N *n, A *a );
+
+/**< Method to find the predecessor of the node in the tree. */
+
+/*--------------------------------------------------------------------------*/
+
+ #if(QUADRATICCOST)
+ template<class A>
+ inline FONumber ReductCost( A *a );
+ #else
+ template<class A>
+ inline CNumber ReductCost( A *a );
+ #endif
+
+/**< Method to calculate the reduct cost of the arc. */
+
+/*--------------------------------------------------------------------------*/
+
+ inline FONumber GetFO( void );
+
+/**< Method to calculate the temporary (or the final) objective function
+ value. */
+
+/*--------------------------------------------------------------------------*/
+
+ void PrintPNode( nodePType *nodo );
+
+/**< Method to print the "name" of the node in the Primal Simplex. */
+
+/*--------------------------------------------------------------------------*/
+
+ void PrintPArc( arcPType *arc );
+
+/**< Method to print the "name" of the arc in the Primal Simplex. */
+
+/*--------------------------------------------------------------------------*/
+
+ void PrintDNode( nodeDType *nodo );
+
+/**< Method to print the "name" of the node in the Dual Simplex. */
+
+/*--------------------------------------------------------------------------*/
+
+ void PrintDArc( arcDType *arc );
+
+/**< Method to print the "name" of the arc in the Dual Simplex. */
+
+/*--------------------------------------------------------------------------*/
+
+ nodePType* RecoverPNode( Index ind );
+
+/**< Method to find a node (in the Primal Simplex) using its index. */
+
+/*--------------------------------------------------------------------------*/
+
+ arcPType* RecoverPArc( nodePType *tail , nodePType *head );
+
+/**< Method to find an arc (in the Primal Simplex) using 2 pointers to tail
+ node and head node. */
+
+/*--------------------------------------------------------------------------*/
+
+ nodeDType* RecoverDNode( Index ind );
+
+/**< Method to find a node (in the Dual Simplex) using its index. */
+
+/*--------------------------------------------------------------------------*/
+
+ arcDType* RecoverDArc( nodeDType *tail , nodeDType *head );
+
+/**< Method to find an arc (in the Dual Simplex) using 2 pointers to tail
+ node and head node. */
+
+/*--------------------------------------------------------------------------*/
+
+ void infoPNode( nodePType *node , int tab );
+
+/**< Method to print some information of the node (in the Primal Simplex). */
+
+/*--------------------------------------------------------------------------*/
+
+ void infoPArc( arcPType *arc , int ind , int tab );
+
+/**< Method to print some information of the arc (in the Primal Simplex). */
+
+/*--------------------------------------------------------------------------*/
+
+ void infoDNode( nodeDType *node , int tab );
+
+/**< Method to print some information of the node (in the Dual Simplex). */
+
+/*--------------------------------------------------------------------------*/
+
+ void infoDArc( arcDType *arc , int ind , int tab );
+
+/**< Method to print some information of the arc (in the Dual Simplex). */
+
+/*--------------------------------------------------------------------------*/
+
+ void ShowSituation( int tab );
+
+/**< Method to show the actual complete situation. */
+
+/*--------------------------------------------------------------------------*/
+
+ }; // end( class MCFSimplex )
+
+/* @} end( group( MCFSimplex_CLASSES ) ) */
+
+#endif /* MCFSimplex.h included */
+
+/*-------------------------------------------------------------------------*/
+/*-------------------inline methods implementation-------------------------*/
+/*-------------------------------------------------------------------------*/
+
+inline MCFClass::Index MCFSimplex::MCFSNde( MCFClass::cIndex i )
+{
+ if( usePrimalSimplex )
+ return( Index( ( (arcsP + i)->tail - nodesP + 1 ) - USENAME0 ) );
+ else
+ return( Index( ( (arcsD + i)->tail - nodesD + 1 ) - USENAME0 ) );
+ }
+
+/*-------------------------------------------------------------------------*/
+
+inline MCFClass::Index MCFSimplex::MCFENde( MCFClass::cIndex i )
+{
+ if( usePrimalSimplex )
+ return( Index( ( (arcsP + i)->head - nodesP + 1 ) - USENAME0 ) );
+ else
+ return( Index( ( (arcsD + i)->head - nodesD + 1 ) - USENAME0 ) );
+ }
+
+/*-------------------------------------------------------------------------*/
+
+inline MCFClass::CNumber MCFSimplex::MCFCost( MCFClass::cIndex i )
+{
+ if( usePrimalSimplex )
+ return( (arcsP + i)->cost );
+ else
+ return( (arcsD + i)->cost );
+ }
+
+/*-------------------------------------------------------------------------*/
+
+inline MCFClass::CNumber MCFSimplex::MCFQCoef( MCFClass::cIndex i )
+{
+ #if( QUADRATICCOST )
+ if( usePrimalSimplex )
+ return( (arcsP + i)->quadraticCost );
+ else
+ return( (arcsD + i)->quadraticCost );
+ #else
+ return( 0 );
+ #endif
+ }
+
+/*-------------------------------------------------------------------------*/
+
+inline MCFClass::FNumber MCFSimplex::MCFUCap( MCFClass::cIndex i )
+{
+ if( usePrimalSimplex )
+ return( (arcsP + i)->upper );
+ else
+ return( (arcsD + i)->upper );
+ }
+
+/*-------------------------------------------------------------------------*/
+
+inline MCFClass::FNumber MCFSimplex::MCFDfct( MCFClass::cIndex i )
+{
+ if( usePrimalSimplex )
+ return( (nodesP + i)->balance );
+ else
+ return( (nodesD + i)->balance );
+ }
+
+/*-------------------------------------------------------------------------*/
+
+#if( QUADRATICCOST )
+
+template <class A>
+inline MCFSimplex::FONumber MCFSimplex::ReductCost( A *a )
+{
+ FONumber redc = (a->tail)->potential - (a->head)->potential;
+ redc = redc + a->cost;
+ redc = redc + a->quadraticCost * a->flow;
+ return( redc );
+ }
+
+#else
+
+template <class A>
+inline MCFSimplex::CNumber MCFSimplex::ReductCost( A *a )
+{
+ CNumber redc = (a->tail)->potential - (a->head)->potential;
+ redc = redc + a->cost;
+ return( redc );
+ }
+
+#endif
+
+/*-------------------------------------------------------------------------*/
+
+#if( OPT_USE_NAMESPACES )
+}; // end( namespace MCFClass_di_unipi_it )
+#endif
+
+/*-------------------------------------------------------------------------*/
+/*-------------------------------------------------------------------------*/
+
+/*-------------------------------------------------------------------------*/
+/*---------------------- End File MCFSimplex.h ----------------------------*/
+/*-------------------------------------------------------------------------*/
diff --git a/clang/lib/Analysis/OPTUtils.h b/clang/lib/Analysis/OPTUtils.h
new file mode 100755
index 0000000..66ad4af
--- /dev/null
+++ b/clang/lib/Analysis/OPTUtils.h
@@ -0,0 +1,475 @@
+/*--------------------------------------------------------------------------*/
+/*--------------------------- File OPTutils.h ------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @file
+ * Small classes are provided for:
+ * - reading the time of a code;
+ * - generating random numbers.
+ *
+ * The classes can be adapted to different environments setting a
+ * compile-time switch in this file.
+ *
+ * Additionally, a function is provided for safely reading parameters
+ * out of a stream.
+ *
+ * \version 1.00
+ *
+ * \date 04 - 10 - 2010
+ *
+ * \author Antonio Frangioni \n
+ * Operations Research Group \n
+ * Dipartimento di Informatica \n
+ * Universita' di Pisa \n
+ *
+ * Copyright &copy 1994 - 2010 by Antonio Frangioni
+ */
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#ifndef __OPTutils
+ #define __OPTutils /* self-identification: #endif at the end of the file */
+
+/*--------------------------------------------------------------------------*/
+/*----------------------------- MACROS -------------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @defgroup OPTUTILS_MACROS Compile-time switches in OPTutils.h
+ These macros control how the classes OPTTimers and OPTrand are
+ implemented; choose the appropriate value for your environment,
+ or program a new version if no value suits you.
+ Also, namespaces can be eliminated if they create problems.
+ @{ */
+
+
+/*----------------------- OPT_USE_NAMESPACES -------------------------------*/
+
+#define OPT_USE_NAMESPACES 0
+
+/**< Setting OPT_USE_NAMESPACES == 0 should instruct all codes that use
+ OPTutils stuff to avoid using namespaces; to start with, the common
+ namespace OPTutils_di_unipi_it, that contains all the types defined
+ herein, is *not* defined. */
+
+/*---------------------------- OPT_TIMERS ----------------------------------*/
+
+#define OPT_TIMERS 5
+
+/**< The class OPTtimers is defined below to give an abstract interface to the
+ different timing routines that are used in different platforms. This is
+ needed since time-related functions are one of the less standard parts of
+ the C[++] library. The value of the OPT_TIMERS constant selects among the
+ different timing routines:
+
+ - 1 = Use the Unix times() routine in sys/times.h
+
+ - 2 = As 1 but uses sys/timeb.h (typical for Microsoft(TM) compilers)
+
+ - 3 = Still use times() of sys/times.h, but return wallclock time
+ rather than CPU time
+
+ - 4 = As 3 but uses sys/timeb.h (typical for Microsoft(TM) compilers)
+
+ - 5 = return the user time obtained with ANSI C clock() function; this
+ may result in more accurate running times w.r.t. but may be limited
+ to ~ 72 hours on systems where ints are 32bits.
+
+ - 6 = Use the Unix gettimeofday() routine of sys/time.h.
+
+ Any unsupported value would simply make the class to report constant
+ zero as the time.
+
+ The values 1 .. 4 rely on the constant CLK_TCK for converting between
+ clock ticks and seconds; for the case where the constant is not defined by
+ the compiler -- should not happen, but it does -- or it is defined in a
+ wrong way, the constant is re-defined below. */
+
+/*---------------------------- OPT_RANDOM ---------------------------------*/
+
+#define OPT_RANDOM 1
+
+/**< The class OPTrand is defined below to give an abstract interface to the
+ different random generators that are used in different platforms. This is
+ needed since random generators are one of the less standard parts of the
+ C[++] library. The value of the OPT_RANDOM constant selects among the
+ different timing routines:
+
+ - 0 = an hand-made implementation of a rather good random number generator
+ is used; note that this assumes that long ints >= 32 bits
+
+ - 1 = standard rand() / srand() pair, common to all C libreries but not
+ very sophisticated
+
+ - 2 = drand48() / srand48(), common on Unix architectures and pretty good.
+
+ Any unsupported value would simply make the functions to report constant
+ zero, which is not nice but useful to quickly fix problems if you don't
+ use random numbers at all. */
+
+/*@} -----------------------------------------------------------------------*/
+/*------------------------------ INCLUDES ----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#if( OPT_RANDOM )
+ #include <stdlib.h>
+
+ /* For random routines, see OPTrand() below. */
+#endif
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
+
+#if( OPT_TIMERS <= 4 )
+ #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 3 ) )
+ #include <sys/times.h>
+ #else
+ #include <sys/timeb.h>
+ #endif
+ #include <limits.h>
+#elif( OPT_TIMERS == 5 )
+ #include <time.h>
+#elif( OPT_TIMERS == 6 )
+ #include <sys/time.h>
+#endif
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
+
+#include <iostream>
+#include <sstream>
+
+/* For istream and the >> operator, used in DfltdSfInpt(). */
+
+/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- NAMESPACE ------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#if( OPT_USE_NAMESPACES )
+namespace OPTtypes_di_unipi_it
+{
+ /** @namespace OPTtypes_di_unipi_it
+ The namespace OPTtypes_di_unipi_it is defined to hold all the data
+ types, constants, classes and functions defined here. It also
+ comprises the namespace std. */
+#endif
+
+ using namespace std; // I know it's not elegant, but ...
+
+/*--------------------------------------------------------------------------*/
+/*----------------------------- NULL ---------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#ifndef NULL
+ #define NULL 0
+#endif
+
+/* Curiously enough, some compilers do not always define NULL; for instance,
+ sometimes it is defined in stdlib.h. */
+
+/*--------------------------------------------------------------------------*/
+/*---------------------------- CLK_TCK -------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#ifndef CLK_TCK
+ #define CLK_TCK 100
+#endif
+
+/* CLK_TCK is the constant used to transform (integer) time values in seconds
+ for times()-based routines. Its normal value is 100. */
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- OPT_TIMERS -----------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @defgroup OPTtypes_CLASSES Classes in OPTutils.h
+ @{ */
+
+#if( OPT_TIMERS )
+
+/** Provides a common interface to the different timing routines that are
+ available in different platforms. */
+
+class OPTtimers {
+
+ public: //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ /// constructor of the class
+ OPTtimers( void )
+ {
+ ReSet();
+ }
+
+ /// start the timer
+ void Start( void )
+ {
+ if( ! ticking ) {
+ #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 2 ) )
+ times( &buff );
+ t_u = buff.tms_utime;
+ t_s = buff.tms_stime;
+ #elif( ( OPT_TIMERS == 3 ) || ( OPT_TIMERS == 4 ) )
+ t_u = times( &buff );
+ #elif( OPT_TIMERS == 5 )
+ t_u = clock();
+ #elif( OPT_TIMERS == 6 )
+ struct timeval t;
+ gettimeofday( &t , NULL );
+ t_u = double( t.tv_sec + t.tv_usec * 1e-6 );
+ #endif
+
+ ticking = 1;
+ }
+ }
+
+ /// stop the timer
+ void Stop( void )
+ {
+ if( ticking ) {
+ Read( u , s );
+ ticking = 0;
+ }
+ }
+
+ /** Return the elapsed time. If the clock is ticking, return the *total*
+ time since the last Start() without stopping the clock; otherwise,
+ return the total elapsed time of all the past runs of the clock since
+ the last ReSet() [see below]. */
+
+ double Read( void )
+ {
+ double tu = 0;
+ double ts = 0;
+ Read( tu , ts );
+ return( tu + ts );
+ }
+
+ /// As Read( void ) but *adds* user and system time to tu and ts.
+ void Read( double &tu , double &ts )
+ {
+ if( ticking ) {
+ #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 2 ) )
+ times( &buff );
+ tu += ( double( buff.tms_utime - t_u ) ) / double( CLK_TCK );
+ ts += ( double( buff.tms_stime - t_s ) ) / double( CLK_TCK );
+ #elif( ( OPT_TIMERS == 3 ) || ( OPT_TIMERS == 4 ) )
+ tu += ( double( times( &buff ) - t_u ) ) / double( CLK_TCK );
+ #elif( OPT_TIMERS == 5 )
+ tu += ( double( clock() - t_u ) ) / double( CLOCKS_PER_SEC );
+ #elif( OPT_TIMERS == 6 )
+ struct timeval t;
+ gettimeofday( &t , NULL );
+ tu += double( t.tv_sec + t.tv_usec * 1e-6 ) - t_u;
+ #endif
+ }
+ else { tu += u; ts += s; }
+ }
+
+ /// reset the timer
+ void ReSet( void )
+ {
+ u = s = 0; ticking = 0;
+ }
+
+ private: //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ double u; // elapsed *user* time, in seconds
+ double s; // elapsed *system* time, in seconds
+ char ticking; // 1 if the clock is ticking
+
+ #if( ( OPT_TIMERS > 0 ) && ( OPT_TIMERS <= 5 ) )
+ clock_t t_u;
+
+ #if( OPT_TIMERS <= 4 )
+ struct tms buff;
+
+ #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 3 ) )
+ clock_t t_s;
+ #endif
+ #endif
+ #elif( OPT_TIMERS == 6 )
+ double t_u;
+ #endif
+
+ }; // end( class OPTtimers );
+
+#endif
+
+/*--------------------------------------------------------------------------*/
+/*------------------------------ OPTrand() ---------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+/** Provide a common interface to the different random generators that are
+ available in different platforms. */
+
+class OPTrand {
+
+ public: //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ /// constructor of the class
+ OPTrand( void )
+ {
+ #if( OPT_RANDOM == 0 )
+ A[ 0 ] = -1;
+ srand( long( 123456789 ) );
+ #else
+ OPTrand::srand( long( 1 ) );
+ #endif
+ }
+
+ /** Returns a random number uniformly distributed in [0, 1).
+ \note each object of class OPTrand has its own sequence, so that
+ multiple OPTrand objects being used within the same program do not
+ interfere with each other (as opposed to what C random routines
+ would do). */
+
+ double rand( void )
+ {
+ #if( OPT_RANDOM == 0 )
+ long nmbr = *(gb_fptr--);
+ if( nmbr < 0 )
+ nmbr = gb_flip_cycle();
+ return( double( nmbr ) / double( (unsigned long)0x80000000 ) );
+ #elif( OPT_RANDOM == 1 )
+ ::srand( myseed );
+ myseed = ::rand();
+ return( double( myseed ) / double( RAND_MAX ) );
+ #elif( OPT_RANDOM == 2 )
+ return( erand48( myseed ) );
+ #else
+ return( 0 ); // random, eh?
+ #endif
+ }
+
+ /// Seeds the random generator for this instance of OPTrand.
+ void srand( long seed )
+ {
+ #if( OPT_RANDOM == 0 )
+ long prev = seed , next = 1;
+ seed = prev = mod_diff( prev , 0 );
+ A[ 55 ] = prev;
+
+ for( long i = 21 ; i ; i = ( i + 21 ) % 55 ) {
+ A[ i ] = next;
+ next = mod_diff( prev , next );
+ if( seed & 1 )
+ seed = 0x40000000 + ( seed >> 1 );
+ else
+ seed >>= 1;
+
+ next = mod_diff( next , seed );
+ prev = A[ i ];
+ }
+
+ gb_flip_cycle();
+ gb_flip_cycle();
+ gb_flip_cycle();
+ gb_flip_cycle();
+ gb_flip_cycle();
+ #elif( OPT_RANDOM == 1 )
+ myseed = int( seed );
+ #elif( OPT_RANDOM == 2 )
+ long *sp = (long*)( &myseed );
+ *sp = seed; // copy higher 32 bits
+ myseed[ 2 ] = 0x330E; // initialize lower 16 bits
+ #endif
+ }
+
+ private: //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+ #if( OPT_RANDOM == 0 )
+ long A[ 56 ];
+ long *gb_fptr;
+
+ long mod_diff( long x , long y )
+ {
+ return( ( ( x ) - ( y ) ) & 0x7fffffff );
+ }
+
+ long gb_flip_cycle( void )
+ {
+ long *ii, *jj;
+ for( ii = &A[ 1 ] , jj = &A[ 32 ] ; jj <= &A[ 55 ] ; ii++ , jj++ )
+ *ii = mod_diff( *ii , *jj );
+
+ for( jj = &A[ 1 ] ; ii <= &A[ 55 ] ; ii++ , jj++ )
+ *ii = mod_diff( *ii , *jj );
+
+ gb_fptr = &A[ 54 ];
+
+ return A[ 55 ];
+ }
+ #elif( OPT_RANDOM == 1 )
+ int myseed;
+ #elif( OPT_RANDOM == 2 )
+ unsigned short int myseed[ 3 ];
+ #endif
+
+ }; // end( class( OPTrand ) )
+
+/* @} end( group( OPTtypes_CLASSES ) ) */
+/*--------------------------------------------------------------------------*/
+/*--------------------------- DfltdSfInpt() --------------------------------*/
+/*--------------------------------------------------------------------------*/
+/** @defgroup OPTtypes_FUNCTIONS Functions in OPTutils.h
+ @{ */
+
+/** Template function for reading parameters from a istream. The function is
+ "safe" because it works also if the istream is not given, is not be long
+ enough or contains erroneous things.
+
+ Given a &istream (possibly NULL), DfltdSfInpt() attempts to read Param out
+ of it, skipping any line that begins with the comment carachter (defaulted
+ to '#'), any blank line and any line starting with anything that can not
+ be interpreted as a `T'. If, for any reason, the read operation fails,
+ then the parameter is given the default value `Dflt'. Otherwise, all the
+ rest of the line up to the nearest newline ('\n') carachter is flushed. */
+
+template<class T>
+inline void DfltdSfInpt( istream *iStrm , T &Param , const T Dflt ,
+ const char cmntc = '#' )
+{
+ string comm( 1 , cmntc );
+
+ // the "! ! stream" trick is there to force the compiler to apply the
+ // stream -> bool conversion, in case it is too dumb to do it by itself
+ if( iStrm && ( ! ( ! (*iStrm) ) ) ) {
+ string buf;
+
+ for( ;; ) {
+ if( ! ( (*iStrm) >> ws ) ) // skip whitespace
+ break;
+
+ if( ! (*iStrm).good() ) // check stream is OK
+ break;
+
+ getline( *iStrm , buf );
+
+ if( ! buf.empty() )
+ if( buf.substr( 0 , 1 ).compare( comm ) ) {
+ stringstream ss;
+ ss << buf;
+ if( ! ( ss >> Param ) )
+ break;
+
+ return;
+ }
+ }
+ }
+
+ Param = Dflt;
+
+ } // end( DfltdSfInpt )
+
+/* @} end( group( OPTtypes_FUNCTIONS ) ) */
+/*--------------------------------------------------------------------------*/
+
+#if( OPT_USE_NAMESPACES )
+ }; // end( namespace OPTtypes_di_unipi_it )
+#endif
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#endif /* OPTutils.h included */
+
+/*--------------------------------------------------------------------------*/
+/*---------------------- End File OPTutils.h -------------------------------*/
+/*--------------------------------------------------------------------------*/