From 51dd6b2b022ade7a1fc4ec8c404d9b81c7e961f5 Mon Sep 17 00:00:00 2001 From: "Zancanaro; Carlo" Date: Tue, 27 Nov 2012 14:56:56 +1100 Subject: Updated solver stuff. This really should have already been in here... --- clang/include/clang/Analysis/Analyses/Interval.h | 18 +- .../Analysis/Analyses/IntervalSolver/Complete.hpp | 5 +- .../Analysis/Analyses/IntervalSolver/MCFClass.h | 2438 ++++++++++++++++++++ .../Analysis/Analyses/IntervalSolver/MCFSimplex.h | 1087 +++++++++ .../Analyses/IntervalSolver/MaxStrategy.hpp | 2 +- .../Analysis/Analyses/IntervalSolver/OPTUtils.h | 475 ++++ .../Analysis/Analyses/IntervalSolver/Operator.hpp | 72 +- .../Analyses/IntervalSolver/VariableAssignment.hpp | 5 +- 8 files changed, 4074 insertions(+), 28 deletions(-) create mode 100755 clang/include/clang/Analysis/Analyses/IntervalSolver/MCFClass.h create mode 100644 clang/include/clang/Analysis/Analyses/IntervalSolver/MCFSimplex.h create mode 100755 clang/include/clang/Analysis/Analyses/IntervalSolver/OPTUtils.h (limited to 'clang/include') diff --git a/clang/include/clang/Analysis/Analyses/Interval.h b/clang/include/clang/Analysis/Analyses/Interval.h index 8bb8d0f..2b37cba 100644 --- a/clang/include/clang/Analysis/Analyses/Interval.h +++ b/clang/include/clang/Analysis/Analyses/Interval.h @@ -19,6 +19,19 @@ #include "llvm/ADT/DenseMap.h" #include "llvm/ADT/ImmutableSet.h" +#include "clang/Analysis/Analyses/IntervalSolver/Complete.hpp" + +#include +#include +#include + + +typedef Complete ZBar; +template<> +inline ZBar infinity() { + return ZBar(1, true); +} + namespace clang { class CFG; @@ -26,14 +39,17 @@ class CFGBlock; class Stmt; class DeclRefExpr; class SourceManager; + +typedef std::vector > ConstraintMatrix; + class IntervalAnalysis : public ManagedAnalysis { public: IntervalAnalysis(AnalysisDeclContext &analysisContext); virtual ~IntervalAnalysis(); - void runOnAllBlocks(const Decl&); + std::map > runOnAllBlocks(const Decl&, const ConstraintMatrix&); static IntervalAnalysis *create(AnalysisDeclContext &analysisContext) { return new IntervalAnalysis(analysisContext); diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp index 73c8f23..2aa657c 100644 --- a/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp +++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/Complete.hpp @@ -10,10 +10,11 @@ template T infinity(); template<> -double infinity() { +inline double infinity() { return std::numeric_limits::infinity(); } + template struct Complete { Complete() @@ -139,7 +140,7 @@ std::ostream& operator<<(std::ostream& cout, const Complete& num) { } template<> -Complete infinity() { +inline Complete infinity() { return Complete(1, true); } diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/MCFClass.h b/clang/include/clang/Analysis/Analyses/IntervalSolver/MCFClass.h new file mode 100755 index 0000000..9534776 --- /dev/null +++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/MCFClass.h @@ -0,0 +1,2438 @@ +/*--------------------------------------------------------------------------*/ +/*-------------------------- File MCFClass.h -------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @file + * Header file for the abstract base class MCFClass, which defines a standard + * interface for (linear or convex quadratic separable) Min Cost Flow Problem + * solvers, to be implemented as derived classes. + * + * \version 3.01 + * + * \date 30 - 09 - 2011 + * + * \author Alessandro Bertolini \n + * Operations Research Group \n + * Dipartimento di Informatica \n + * Universita' di Pisa \n + * + * \author Antonio Frangioni \n + * Operations Research Group \n + * Dipartimento di Informatica \n + * Universita' di Pisa \n + * + * \author Claudio Gentile \n + * Istituto di Analisi di Sistemi e Informatica \n + * Consiglio Nazionale delle Ricerche \n + * + * Copyright © 1996 - 2011 by Antonio Frangioni, Claudio Gentile + */ + +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*----------------------------- DEFINITIONS --------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#ifndef __MCFClass + #define __MCFClass /* self-identification: #endif at the end of the file */ + +/*--------------------------------------------------------------------------*/ +/*--------------------------------- MACROS ---------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @defgroup MCFCLASS_MACROS Compile-time switches in MCFClass.h + These macros control some important details of the class interface. + Although using macros for activating features of the interface is not + very C++, switching off some unused features may allow some + implementation to be more efficient in running time or memory. + @{ */ + +/*-------------------------------- USENAME0 --------------------------------*/ + +#define USENAME0 1 + +/**< Decides if 0 or 1 is the "name" of the first node. + If USENAME0 == 1, (warning: it has to be *exactly* 1), then the node + names go from 0 to n - 1, otherwise from 1 to n. Note that this does not + affect the position of the deficit in the deficit vectors, i.e., the + deficit of the i-th node - be its "name" `i' or `i - 1' - is always in + the i-th position of the vector. */ + +/*@} end( group( MCFCLASS_MACROS ) ) */ +/*--------------------------------------------------------------------------*/ +/*------------------------------ INCLUDES ----------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#include "OPTUtils.h" + +/* OPTtypes.h defines standard interfaces for timing and random routines, as + well as the namespace OPTtypes_di_unipi_it and the macro + OPT_USE_NAMESPACES, useful for switching off all namespaces in one blow + for those strange cases where they create problems. */ + +#include +#include +#include +#include + +// 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(). */ + + template + class Inf { + public: + Inf() {} + operator T() { return( std::numeric_limits::max() ); } + }; + +/*--------------------------------------------------------------------------*/ +/** Very small class to simplify extracting the "machine epsilon" for a + basic type (FNumber, CNumber); just use Eps(). */ + + template + class Eps { + public: + Eps() {} + operator T() { return( std::numeric_limits::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() * 100; + EpsCst = Eps() * 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 + + Then the node definition lines must be found, in the form + + n + + 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 + + There must be exactly 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, . 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() ) = 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()-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() ) = 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()-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() ) = 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()-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() if MCFGetStatus() == kUnfeasible and + - Inf() 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() ); + 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(); + 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()-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 , Index_Set ) + { + return( Inf() ); + } + +/**< 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() 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 ) {} + +/**< 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() ) = 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()-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() ) = 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()-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() ) + +/**< 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()-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 ) + +/**< 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() ) = 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()-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() ) = 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()-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() ) = 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()-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 = NULL , + cIndex = 0 , Index = Inf() ) //HERE + +/**< 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()-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 , 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() ) = 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()-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() ) = 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()-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() 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() , + cIndex nEN = Inf() ) = 0; + +/**< Change the starting and/or ending node of arc `name' to nSN and nEN. + Each parameter being Inf() 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() 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 + inline bool ETZ( T x , const T eps ) + { + if( numeric_limits::is_integer ) + return( x == 0 ); + else + return( (x <= eps ) && ( x >= -eps ) ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is greater than zero (possibly considering tolerances). */ + + template + inline bool GTZ( T x , const T eps ) + { + if( numeric_limits::is_integer ) + return( x > 0 ); + else + return( x > eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is greater than or equal to zero (possibly considering + tolerances). */ + + template + inline bool GEZ( T x , const T eps ) + { + if( numeric_limits::is_integer ) + return( x >= 0 ); + else + return( x >= - eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is less than zero (possibly considering tolerances). */ + + template + inline bool LTZ( T x , const T eps ) + { + if( numeric_limits::is_integer ) + return( x < 0 ); + else + return( x < - eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is less than or equal to zero (possibly considering + tolerances). */ + + template + inline bool LEZ( T x , const T eps ) + { + if( numeric_limits::is_integer ) + return( x <= 0 ); + else + return( x <= eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is greater than flow y (possibly considering tolerances). + */ + + template + inline bool GT( T x , T y , const T eps ) + { + if( numeric_limits::is_integer ) + return( x > y ); + else + return( x > y + eps ); + } + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ +/** true if flow x is less than flow y (possibly considering tolerances). */ + + template + inline bool LT( T x , T y , const T eps ) + { + if( numeric_limits::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::is_integer ) || + ( ! numeric_limits::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 + +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#undef throw + +#endif /* MCFClass.h included */ + +/*--------------------------------------------------------------------------*/ +/*------------------------- End File MCFClass.h ----------------------------*/ +/*--------------------------------------------------------------------------*/ diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/MCFSimplex.h b/clang/include/clang/Analysis/Analyses/IntervalSolver/MCFSimplex.h new file mode 100644 index 0000000..9174337 --- /dev/null +++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/MCFSimplex.h @@ -0,0 +1,1087 @@ +/*--------------------------------------------------------------------------*/ +/*------------------------- File MCFSimplex.h ------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @file + * Linear and Quadratic Min Cost Flow problems solver based on the (primal and + * dual) simplex algorithm. Conforms to the standard MCF interface defined in + * MCFClass.h. + * + * \Version 1.00 + * + * \date 29 - 08 - 2011 + * + * \author Alessandro Bertolini \n + * Antonio Frangioni \n + * Operations Research Group \n + * Dipartimento di Informatica \n + * Universita' di Pisa \n + * + * Copyright © 2008 - 2011 by Alessandro Bertolini, Antonio Frangioni + */ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*----------------------------- DEFINITIONS --------------------------------*/ +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#ifndef __MCFSimplex + #define __MCFSimplex /* self-identification: #endif at the end of the file */ + +/*--------------------------------------------------------------------------*/ +/*------------------------------ INCLUDES ----------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#include "MCFClass.h" + +/*@} -----------------------------------------------------------------------*/ +/*--------------------------- NAMESPACE ------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#if( OPT_USE_NAMESPACES ) +namespace MCFClass_di_unipi_it +{ +#endif + +/*--------------------------------------------------------------------------*/ +/*-------------------------------- MACROS ----------------------------------*/ +/*--------------------------------------------------------------------------*/ +/** @defgroup MCFSimplex_MACROS Compile-time switches in MCFSimplex.h + There is only one macro in MCFSimplex, but it is very important! + @{ */ + +#define QUADRATICCOST 0 + +/**< Setting QUADRATICCOST == 1 the solver can solve problems with linear and + quadratic costs too (but the latter only with the Primal Simplex). + The reason for having a macro is that when quadratic costs are present the + "arcType" struct has the additional field "quadraticCost" to hold it. + Furthermore, the field "ident" is not created because the solver doesn't + use the classical TLU tripartition. Instead, closed arcs and deleted arcs + are characterized as follows: + - closed arcs have the field "cost" to INFINITY (Inf()); + - 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() ); + +/*--------------------------------------------------------------------------*/ + + void MCFGetRC( CRow CR , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf() ) ; + + CNumber MCFGetRC( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFGetPi( CRow P , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf() ); + +/*--------------------------------------------------------------------------*/ + + 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() ); + + inline Index MCFSNde( cIndex i ); + + inline Index MCFENde( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFCosts( CRow Costv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf() ); + + inline CNumber MCFCost( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFQCoef( CRow Qv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf() ); + + inline CNumber MCFQCoef( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFUCaps( FRow UCapv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf() ); + + inline FNumber MCFUCap( cIndex i ); + +/*--------------------------------------------------------------------------*/ + + void MCFDfcts( FRow Dfctv , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf() ); + + 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() ); + + void ChgCost( Index arc , cCNumber NCost ); + +/*--------------------------------------------------------------------------*/ + + void ChgQCoef( cCRow NQCoef = NULL , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf() ); + + void ChgQCoef( Index arc , cCNumber NQCoef ); + +/*--------------------------------------------------------------------------*/ + + void ChgDfcts( cFRow NDfct , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf() ); + + void ChgDfct( Index nod , cFNumber NDfct ); + +/*--------------------------------------------------------------------------*/ + + void ChgUCaps( cFRow NCap , cIndex_Set nms = NULL , + cIndex strt = 0 , Index stp = Inf() ); + + 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() , + cIndex nEN = Inf() ); + + 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 + 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 + 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 + 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 + 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 + 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 + inline N* Father( N *n, A *a ); + +/**< Method to find the predecessor of the node in the tree. */ + +/*--------------------------------------------------------------------------*/ + + #if(QUADRATICCOST) + template + inline FONumber ReductCost( A *a ); + #else + template + 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 ) +{ + #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 +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 +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 + +//#endif /* MCFSimplex.h included */ +/*-------------------------------------------------------------------------*/ +/*-------------------------------------------------------------------------*/ + +/*-------------------------------------------------------------------------*/ +/*---------------------- End File MCFSimplex.h ----------------------------*/ +/*-------------------------------------------------------------------------*/ diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/MaxStrategy.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/MaxStrategy.hpp index da05dd8..5aca4b2 100644 --- a/clang/include/clang/Analysis/Analyses/IntervalSolver/MaxStrategy.hpp +++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/MaxStrategy.hpp @@ -186,7 +186,7 @@ struct Solver { } Domain solve(Variable& var) { - MaxExpression& rhs = *_system[var]; + MaxExpression& rhs = *_system[var]; // this will automatically work sufficiently to get the final // strategy for this variable's RHS _strategy.get(rhs); diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/OPTUtils.h b/clang/include/clang/Analysis/Analyses/IntervalSolver/OPTUtils.h new file mode 100755 index 0000000..66ad4af --- /dev/null +++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/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 © 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 + + /* For random routines, see OPTrand() below. */ +#endif + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ + +#if( OPT_TIMERS <= 4 ) + #if( ( OPT_TIMERS == 1 ) || ( OPT_TIMERS == 3 ) ) + #include + #else + #include + #endif + #include +#elif( OPT_TIMERS == 5 ) + #include +#elif( OPT_TIMERS == 6 ) + #include +#endif + +/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ + +#include +#include + +/* 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 +inline void DfltdSfInpt( istream *iStrm , T &Param , const T Dflt , + const char cmntc = '#' ) +{ + string comm( 1 , cmntc ); + + // the "! ! stream" trick is there to force the compiler to apply the + // stream -> bool conversion, in case it is too dumb to do it by itself + if( iStrm && ( ! ( ! (*iStrm) ) ) ) { + string buf; + + for( ;; ) { + if( ! ( (*iStrm) >> ws ) ) // skip whitespace + break; + + if( ! (*iStrm).good() ) // check stream is OK + break; + + getline( *iStrm , buf ); + + if( ! buf.empty() ) + if( buf.substr( 0 , 1 ).compare( comm ) ) { + stringstream ss; + ss << buf; + if( ! ( ss >> Param ) ) + break; + + return; + } + } + } + + Param = Dflt; + + } // end( DfltdSfInpt ) + +/* @} end( group( OPTtypes_FUNCTIONS ) ) */ +/*--------------------------------------------------------------------------*/ + +#if( OPT_USE_NAMESPACES ) + }; // end( namespace OPTtypes_di_unipi_it ) +#endif + +/*--------------------------------------------------------------------------*/ +/*--------------------------------------------------------------------------*/ + +#endif /* OPTutils.h included */ + +/*--------------------------------------------------------------------------*/ +/*---------------------- End File OPTutils.h -------------------------------*/ +/*--------------------------------------------------------------------------*/ diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/Operator.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/Operator.hpp index 8dbdc39..98ac6fd 100644 --- a/clang/include/clang/Analysis/Analyses/IntervalSolver/Operator.hpp +++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/Operator.hpp @@ -3,6 +3,7 @@ #include #include +#include template struct Operator { @@ -141,9 +142,9 @@ template struct MinCostFlow : public Operator { MinCostFlow(const std::vector& supplies, const std::vector > arcs) : _supplies(supplies), _arcs(arcs), _solver(0,0) { - MCFClass::FNumber* deficits = new MCFClass::FNumber[supplies.size()]; - MCFClass::Index* starts = new MCFClass::Index[arcs.size()]; - MCFClass::Index* ends = new MCFClass::Index[arcs.size()]; + deficits = new MCFClass::FNumber[supplies.size()]; + starts = new MCFClass::Index[arcs.size()]; + ends = new MCFClass::Index[arcs.size()]; for (int i = 0, size = supplies.size(); i < size; ++i) { deficits[i] = -supplies[i].template as(); @@ -152,13 +153,8 @@ struct MinCostFlow : public Operator { starts[i] = arcs[i].first; ends[i] = arcs[i].second; } - - _solver.LoadNet(supplies.size(), arcs.size(), // max nodes/arcs - supplies.size(), arcs.size(), // current nodes/arcs - NULL, NULL, // arcs have inf cap, zero cost (to begin) - deficits, // deficits for each node - starts, ends); // start/end of each edge - + } + ~MinCostFlow() { delete[] deficits; delete[] starts; delete[] ends; @@ -167,23 +163,50 @@ struct MinCostFlow : public Operator { assert(costs.size() == _arcs.size()); if (costs.size() == 0) return 0; + + _solver.LoadNet(_supplies.size(), _arcs.size(), // max nodes/arcs + _supplies.size(), _arcs.size(), // current nodes/arcs + NULL, NULL, // arcs have inf cap, zero cost (to begin) + deficits, // deficits for each node + starts, ends); // start/end of each edge + for (int i = 0, size = costs.size(); i < size; ++i) { _solver.ChgCost(i, costs[i].template as()); } - _solver.SolveMCF(); - if (_solver.MCFGetStatus() != MCFClass::kOK) { - return -infinity(); - } else { - MCFClass::FONumber num = _solver.MCFGetFO(); - if (num == MCFClass::Inf()) { - return infinity(); - } else if (num == -MCFClass::Inf()) { - return -infinity(); - } else { - if (((int)num) == -2147483648) + do { + _solver.SolveMCF(); + } while (_solver.MCFGetStatus() == MCFClass::kStopped); + + MCFClass::FONumber num; + switch (_solver.MCFGetStatus()) { + case MCFClass::kOK: + num = _solver.MCFGetFO(); + /*std::cout << "MCF solved for: " << num << std::endl; + std::cout << " with arguments ["; + for (unsigned int i = 0; i < costs.size(); ++i) { + if (i > 0) + std::cout << ", "; + std::cout << costs[i]; + } + std::cout << "]" << std::endl;*/ + if (num == MCFClass::Inf()) { + return infinity(); + } else if (num == -MCFClass::Inf()) { return -infinity(); - return num; - } + } else { + if ((long)num == std::numeric_limits::min()) + return -infinity(); + return num; + } + + case MCFClass::kUnbounded: + return -infinity(); + + case MCFClass::kUnfeasible: + return infinity(); + + default: + assert(false); } } void print(std::ostream& cout) const { @@ -220,6 +243,9 @@ private: std::vector _supplies; std::vector > _arcs; mutable MCFSimplex _solver; + MCFClass::FNumber* deficits; + MCFClass::Index* starts; + MCFClass::Index* ends; }; template diff --git a/clang/include/clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp b/clang/include/clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp index 4c4519a..04c4ea8 100644 --- a/clang/include/clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp +++ b/clang/include/clang/Analysis/Analyses/IntervalSolver/VariableAssignment.hpp @@ -64,11 +64,14 @@ struct DynamicVariableAssignment : public VariableAssignment { it = _touched.begin(), ei = _touched.end(); it != ei; - ++it) { + ) { Variable& var = _system.variable(*it); if (!_unstable.contains(var) && _old_values[var] != _values[var]) { changed.insert(var); + it++; _touched.remove(var); + } else { + it++; } } //_touched.clear(); -- cgit v1.2.3