Updated solver stuff. This really should have already been in here...

11 files changed, 4391 insertions, 171 deletions

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