Merge
authorAlpar Juttner <alpar@cs.elte.hu>
Wed, 29 Apr 2009 19:22:14 +0100
changeset 693e01957e96c67
parent 692 cb8270a98660
parent 691 8d289c89d43e
child 694 dcba640438c7
child 697 a8dfe89b7719
child 698 3adf5e2d1e62
child 842 dc9316203edf
Merge
     1.1 --- a/CMakeLists.txt	Wed Apr 29 16:15:29 2009 +0100
     1.2 +++ b/CMakeLists.txt	Wed Apr 29 19:22:14 2009 +0100
     1.3 @@ -17,8 +17,6 @@
     1.4  FIND_PACKAGE(CPLEX)
     1.5  FIND_PACKAGE(COIN)
     1.6  
     1.7 -ADD_DEFINITIONS(-DHAVE_CONFIG_H)
     1.8 -
     1.9  IF(MSVC)
    1.10    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4250 /wd4355 /wd4800 /wd4996")
    1.11  # Suppressed warnings:
    1.12 @@ -28,8 +26,6 @@
    1.13  # C4996: 'function': was declared deprecated
    1.14  ENDIF(MSVC)
    1.15  
    1.16 -ADD_DEFINITIONS(-DHAVE_CONFIG_H)
    1.17 -
    1.18  INCLUDE(CheckTypeSize)
    1.19  CHECK_TYPE_SIZE("long long" LEMON_LONG_LONG)
    1.20  
     2.1 --- a/Makefile.am	Wed Apr 29 16:15:29 2009 +0100
     2.2 +++ b/Makefile.am	Wed Apr 29 19:22:14 2009 +0100
     2.3 @@ -11,11 +11,12 @@
     2.4  	m4/lx_check_cplex.m4 \
     2.5  	m4/lx_check_glpk.m4 \
     2.6  	m4/lx_check_soplex.m4 \
     2.7 -	m4/lx_check_clp.m4 \
     2.8 -	m4/lx_check_cbc.m4 \
     2.9 +	m4/lx_check_coin.m4 \
    2.10  	CMakeLists.txt \
    2.11  	cmake/FindGhostscript.cmake \
    2.12 +	cmake/FindCPLEX.cmake \
    2.13  	cmake/FindGLPK.cmake \
    2.14 +	cmake/FindCOIN.cmake \
    2.15  	cmake/version.cmake.in \
    2.16  	cmake/version.cmake \
    2.17  	cmake/nsis/lemon.ico \
     3.1 --- a/cmake/FindCOIN.cmake	Wed Apr 29 16:15:29 2009 +0100
     3.2 +++ b/cmake/FindCOIN.cmake	Wed Apr 29 19:22:14 2009 +0100
     3.3 @@ -1,28 +1,48 @@
     3.4  SET(COIN_ROOT_DIR "" CACHE PATH "COIN root directory")
     3.5  
     3.6  FIND_PATH(COIN_INCLUDE_DIR coin/CoinUtilsConfig.h
     3.7 -  PATHS ${COIN_ROOT_DIR}/include)
     3.8 -
     3.9 -FIND_LIBRARY(COIN_CBC_LIBRARY libCbc
    3.10 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.11 -FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY libCbcSolver
    3.12 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.13 -FIND_LIBRARY(COIN_CGL_LIBRARY libCgl
    3.14 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.15 -FIND_LIBRARY(COIN_CLP_LIBRARY libClp
    3.16 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.17 -FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY libCoinUtils
    3.18 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.19 -FIND_LIBRARY(COIN_OSI_LIBRARY libOsi
    3.20 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.21 -FIND_LIBRARY(COIN_OSI_CBC_LIBRARY libOsiCbc
    3.22 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.23 -FIND_LIBRARY(COIN_OSI_CLP_LIBRARY libOsiClp
    3.24 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.25 -FIND_LIBRARY(COIN_OSI_VOL_LIBRARY libOsiVol
    3.26 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.27 -FIND_LIBRARY(COIN_VOL_LIBRARY libVol
    3.28 -  PATHS ${COIN_ROOT_DIR}/lib)
    3.29 +  HINTS ${COIN_ROOT_DIR}/include
    3.30 +)
    3.31 +FIND_LIBRARY(COIN_CBC_LIBRARY
    3.32 +  NAMES Cbc libCbc
    3.33 +  HINTS ${COIN_ROOT_DIR}/lib
    3.34 +)
    3.35 +FIND_LIBRARY(COIN_CBC_SOLVER_LIBRARY
    3.36 +  NAMES CbcSolver libCbcSolver
    3.37 +  HINTS ${COIN_ROOT_DIR}/lib
    3.38 +)
    3.39 +FIND_LIBRARY(COIN_CGL_LIBRARY
    3.40 +  NAMES Cgl libCgl
    3.41 +  HINTS ${COIN_ROOT_DIR}/lib
    3.42 +)
    3.43 +FIND_LIBRARY(COIN_CLP_LIBRARY
    3.44 +  NAMES Clp libClp
    3.45 +  HINTS ${COIN_ROOT_DIR}/lib
    3.46 +)
    3.47 +FIND_LIBRARY(COIN_COIN_UTILS_LIBRARY
    3.48 +  NAMES CoinUtils libCoinUtils
    3.49 +  HINTS ${COIN_ROOT_DIR}/lib
    3.50 +)
    3.51 +FIND_LIBRARY(COIN_OSI_LIBRARY
    3.52 +  NAMES Osi libOsi
    3.53 +  HINTS ${COIN_ROOT_DIR}/lib
    3.54 +)
    3.55 +FIND_LIBRARY(COIN_OSI_CBC_LIBRARY
    3.56 +  NAMES OsiCbc libOsiCbc
    3.57 +  HINTS ${COIN_ROOT_DIR}/lib
    3.58 +)
    3.59 +FIND_LIBRARY(COIN_OSI_CLP_LIBRARY
    3.60 +  NAMES OsiClp libOsiClp
    3.61 +  HINTS ${COIN_ROOT_DIR}/lib
    3.62 +)
    3.63 +FIND_LIBRARY(COIN_OSI_VOL_LIBRARY
    3.64 +  NAMES OsiVol libOsiVol
    3.65 +  HINTS ${COIN_ROOT_DIR}/lib
    3.66 +)
    3.67 +FIND_LIBRARY(COIN_VOL_LIBRARY
    3.68 +  NAMES Vol libVol
    3.69 +  HINTS ${COIN_ROOT_DIR}/lib
    3.70 +)
    3.71  
    3.72  INCLUDE(FindPackageHandleStandardArgs)
    3.73  FIND_PACKAGE_HANDLE_STANDARD_ARGS(COIN DEFAULT_MSG
     4.1 --- a/cmake/FindCPLEX.cmake	Wed Apr 29 16:15:29 2009 +0100
     4.2 +++ b/cmake/FindCPLEX.cmake	Wed Apr 29 19:22:14 2009 +0100
     4.3 @@ -1,21 +1,32 @@
     4.4 +SET(CPLEX_ROOT_DIR "" CACHE PATH "CPLEX root directory")
     4.5 +
     4.6  FIND_PATH(CPLEX_INCLUDE_DIR
     4.7    ilcplex/cplex.h
     4.8 -  PATHS "C:/ILOG/CPLEX91/include")
     4.9 -
    4.10 +  PATHS "C:/ILOG/CPLEX91/include"
    4.11 +  PATHS "/opt/ilog/cplex91/include"
    4.12 +  HINTS ${CPLEX_ROOT_DIR}/include
    4.13 +)
    4.14  FIND_LIBRARY(CPLEX_LIBRARY
    4.15 -  NAMES cplex91
    4.16 -  PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda")
    4.17 +  cplex91
    4.18 +  PATHS "C:/ILOG/CPLEX91/lib/msvc7/stat_mda"
    4.19 +  PATHS "/opt/ilog/cplex91/bin"
    4.20 +  HINTS ${CPLEX_ROOT_DIR}/bin
    4.21 +)
    4.22  
    4.23  INCLUDE(FindPackageHandleStandardArgs)
    4.24  FIND_PACKAGE_HANDLE_STANDARD_ARGS(CPLEX DEFAULT_MSG CPLEX_LIBRARY CPLEX_INCLUDE_DIR)
    4.25  
    4.26  FIND_PATH(CPLEX_BIN_DIR
    4.27    cplex91.dll
    4.28 -  PATHS "C:/ILOG/CPLEX91/bin/x86_win32")
    4.29 +  PATHS "C:/ILOG/CPLEX91/bin/x86_win32"
    4.30 +)
    4.31  
    4.32  IF(CPLEX_FOUND)
    4.33    SET(CPLEX_INCLUDE_DIRS ${CPLEX_INCLUDE_DIR})
    4.34    SET(CPLEX_LIBRARIES ${CPLEX_LIBRARY})
    4.35 +  IF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
    4.36 +    SET(CPLEX_LIBRARIES "${CPLEX_LIBRARIES};m;pthread")
    4.37 +  ENDIF(CMAKE_SYSTEM_NAME STREQUAL "Linux")
    4.38  ENDIF(CPLEX_FOUND)
    4.39  
    4.40  MARK_AS_ADVANCED(CPLEX_LIBRARY CPLEX_INCLUDE_DIR CPLEX_BIN_DIR)
     5.1 --- a/cmake/FindGLPK.cmake	Wed Apr 29 16:15:29 2009 +0100
     5.2 +++ b/cmake/FindGLPK.cmake	Wed Apr 29 19:22:14 2009 +0100
     5.3 @@ -1,16 +1,50 @@
     5.4 +SET(GLPK_ROOT_DIR "" CACHE PATH "GLPK root directory")
     5.5 +
     5.6  SET(GLPK_REGKEY "[HKEY_LOCAL_MACHINE\\SOFTWARE\\GnuWin32\\Glpk;InstallPath]")
     5.7  GET_FILENAME_COMPONENT(GLPK_ROOT_PATH ${GLPK_REGKEY} ABSOLUTE)
     5.8  
     5.9  FIND_PATH(GLPK_INCLUDE_DIR
    5.10    glpk.h
    5.11 -  PATHS ${GLPK_REGKEY}/include)
    5.12 +  PATHS ${GLPK_REGKEY}/include
    5.13 +  HINTS ${GLPK_ROOT_DIR}/include
    5.14 +)
    5.15 +FIND_LIBRARY(GLPK_LIBRARY
    5.16 +  glpk
    5.17 +  PATHS ${GLPK_REGKEY}/lib
    5.18 +  HINTS ${GLPK_ROOT_DIR}/lib
    5.19 +)
    5.20  
    5.21 -FIND_LIBRARY(GLPK_LIBRARY
    5.22 -  NAMES glpk
    5.23 -  PATHS ${GLPK_REGKEY}/lib)
    5.24 +IF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY)
    5.25 +  FILE(READ ${GLPK_INCLUDE_DIR}/glpk.h GLPK_GLPK_H)
    5.26 +
    5.27 +  STRING(REGEX MATCH "define[ ]+GLP_MAJOR_VERSION[ ]+[0-9]+" GLPK_MAJOR_VERSION_LINE "${GLPK_GLPK_H}")
    5.28 +  STRING(REGEX REPLACE "define[ ]+GLP_MAJOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MAJOR "${GLPK_MAJOR_VERSION_LINE}")
    5.29 +
    5.30 +  STRING(REGEX MATCH "define[ ]+GLP_MINOR_VERSION[ ]+[0-9]+" GLPK_MINOR_VERSION_LINE "${GLPK_GLPK_H}")
    5.31 +  STRING(REGEX REPLACE "define[ ]+GLP_MINOR_VERSION[ ]+([0-9]+)" "\\1" GLPK_VERSION_MINOR "${GLPK_MINOR_VERSION_LINE}")
    5.32 +
    5.33 +  SET(GLPK_VERSION_STRING "${GLPK_VERSION_MAJOR}.${GLPK_VERSION_MINOR}")
    5.34 +
    5.35 +  IF(GLPK_FIND_VERSION)
    5.36 +    IF(GLPK_FIND_VERSION_COUNT GREATER 2)
    5.37 +      MESSAGE(SEND_ERROR "unexpected version string")
    5.38 +    ENDIF(GLPK_FIND_VERSION_COUNT GREATER 2)
    5.39 +
    5.40 +    MATH(EXPR GLPK_REQUESTED_VERSION "${GLPK_FIND_VERSION_MAJOR}*100 + ${GLPK_FIND_VERSION_MINOR}")
    5.41 +    MATH(EXPR GLPK_FOUND_VERSION "${GLPK_VERSION_MAJOR}*100 + ${GLPK_VERSION_MINOR}")
    5.42 +
    5.43 +    IF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
    5.44 +      SET(GLPK_PROPER_VERSION_FOUND FALSE)
    5.45 +    ELSE(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
    5.46 +      SET(GLPK_PROPER_VERSION_FOUND TRUE)
    5.47 +    ENDIF(GLPK_FOUND_VERSION LESS GLPK_REQUESTED_VERSION)
    5.48 +  ELSE(GLPK_FIND_VERSION)
    5.49 +    SET(GLPK_PROPER_VERSION_FOUND TRUE)
    5.50 +  ENDIF(GLPK_FIND_VERSION)
    5.51 +ENDIF(GLPK_INCLUDE_DIR AND GLPK_LIBRARY)
    5.52  
    5.53  INCLUDE(FindPackageHandleStandardArgs)
    5.54 -FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR)
    5.55 +FIND_PACKAGE_HANDLE_STANDARD_ARGS(GLPK DEFAULT_MSG GLPK_LIBRARY GLPK_INCLUDE_DIR GLPK_PROPER_VERSION_FOUND)
    5.56  
    5.57  IF(GLPK_FOUND)
    5.58    SET(GLPK_INCLUDE_DIRS ${GLPK_INCLUDE_DIR})
     6.1 --- a/doc/groups.dox	Wed Apr 29 16:15:29 2009 +0100
     6.2 +++ b/doc/groups.dox	Wed Apr 29 19:22:14 2009 +0100
     6.3 @@ -352,17 +352,17 @@
     6.4  minimum total cost from a set of supply nodes to a set of demand nodes
     6.5  in a network with capacity constraints (lower and upper bounds)
     6.6  and arc costs.
     6.7 -Formally, let \f$G=(V,A)\f$ be a digraph,
     6.8 -\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
     6.9 +Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
    6.10 +\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
    6.11  upper bounds for the flow values on the arcs, for which
    6.12 -\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
    6.13 -\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
    6.14 -on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
    6.15 +\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
    6.16 +\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
    6.17 +on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
    6.18  signed supply values of the nodes.
    6.19  If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
    6.20  supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
    6.21  \f$-sup(u)\f$ demand.
    6.22 -A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
    6.23 +A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
    6.24  of the following optimization problem.
    6.25  
    6.26  \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    6.27 @@ -404,7 +404,7 @@
    6.28  
    6.29  The dual solution of the minimum cost flow problem is represented by node 
    6.30  potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
    6.31 -An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
    6.32 +An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
    6.33  is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
    6.34  node potentials the following \e complementary \e slackness optimality
    6.35  conditions hold.
    6.36 @@ -413,15 +413,15 @@
    6.37     - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
    6.38     - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
    6.39     - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
    6.40 - - For all \f$u\in V\f$:
    6.41 + - For all \f$u\in V\f$ nodes:
    6.42     - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
    6.43       then \f$\pi(u)=0\f$.
    6.44   
    6.45  Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
    6.46 -\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
    6.47 +\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
    6.48  \f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
    6.49  
    6.50 -All algorithms provide dual solution (node potentials) as well
    6.51 +All algorithms provide dual solution (node potentials) as well,
    6.52  if an optimal flow is found.
    6.53  
    6.54  LEMON contains several algorithms for solving minimum cost flow problems.
     7.1 --- a/lemon/Makefile.am	Wed Apr 29 16:15:29 2009 +0100
     7.2 +++ b/lemon/Makefile.am	Wed Apr 29 19:22:14 2009 +0100
     7.3 @@ -15,7 +15,8 @@
     7.4  	lemon/random.cc \
     7.5  	lemon/bits/windows.cc
     7.6  
     7.7 -
     7.8 +nodist_lemon_HEADERS = lemon/config.h	
     7.9 +	
    7.10  lemon_libemon_la_CXXFLAGS = \
    7.11  	$(AM_CXXFLAGS) \
    7.12  	$(GLPK_CFLAGS) \
    7.13 @@ -57,11 +58,11 @@
    7.14  	lemon/assert.h \
    7.15  	lemon/bfs.h \
    7.16  	lemon/bin_heap.h \
    7.17 +	lemon/cbc.h \
    7.18  	lemon/circulation.h \
    7.19  	lemon/clp.h \
    7.20  	lemon/color.h \
    7.21  	lemon/concept_check.h \
    7.22 -	lemon/config.h \
    7.23  	lemon/connectivity.h \
    7.24  	lemon/counter.h \
    7.25  	lemon/core.h \
     8.1 --- a/lemon/circulation.h	Wed Apr 29 16:15:29 2009 +0100
     8.2 +++ b/lemon/circulation.h	Wed Apr 29 19:22:14 2009 +0100
     8.3 @@ -64,15 +64,15 @@
     8.4      /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
     8.5      typedef SM SupplyMap;
     8.6  
     8.7 -    /// \brief The type of the flow values.
     8.8 -    typedef typename SupplyMap::Value Flow;
     8.9 +    /// \brief The type of the flow and supply values.
    8.10 +    typedef typename SupplyMap::Value Value;
    8.11  
    8.12      /// \brief The type of the map that stores the flow values.
    8.13      ///
    8.14      /// The type of the map that stores the flow values.
    8.15      /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
    8.16      /// concept.
    8.17 -    typedef typename Digraph::template ArcMap<Flow> FlowMap;
    8.18 +    typedef typename Digraph::template ArcMap<Value> FlowMap;
    8.19  
    8.20      /// \brief Instantiates a FlowMap.
    8.21      ///
    8.22 @@ -104,7 +104,7 @@
    8.23      /// \brief The tolerance used by the algorithm
    8.24      ///
    8.25      /// The tolerance used by the algorithm to handle inexact computation.
    8.26 -    typedef lemon::Tolerance<Flow> Tolerance;
    8.27 +    typedef lemon::Tolerance<Value> Tolerance;
    8.28  
    8.29    };
    8.30  
    8.31 @@ -187,8 +187,8 @@
    8.32      typedef TR Traits;
    8.33      ///The type of the digraph the algorithm runs on.
    8.34      typedef typename Traits::Digraph Digraph;
    8.35 -    ///The type of the flow values.
    8.36 -    typedef typename Traits::Flow Flow;
    8.37 +    ///The type of the flow and supply values.
    8.38 +    typedef typename Traits::Value Value;
    8.39  
    8.40      ///The type of the lower bound map.
    8.41      typedef typename Traits::LowerMap LowerMap;
    8.42 @@ -221,7 +221,7 @@
    8.43      Elevator* _level;
    8.44      bool _local_level;
    8.45  
    8.46 -    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
    8.47 +    typedef typename Digraph::template NodeMap<Value> ExcessMap;
    8.48      ExcessMap* _excess;
    8.49  
    8.50      Tolerance _tol;
    8.51 @@ -530,7 +530,7 @@
    8.52            (*_excess)[_g.target(e)] += (*_lo)[e];
    8.53            (*_excess)[_g.source(e)] -= (*_lo)[e];
    8.54          } else {
    8.55 -          Flow fc = -(*_excess)[_g.target(e)];
    8.56 +          Value fc = -(*_excess)[_g.target(e)];
    8.57            _flow->set(e, fc);
    8.58            (*_excess)[_g.target(e)] = 0;
    8.59            (*_excess)[_g.source(e)] -= fc;
    8.60 @@ -563,11 +563,11 @@
    8.61        while((act=_level->highestActive())!=INVALID) {
    8.62          int actlevel=(*_level)[act];
    8.63          int mlevel=_node_num;
    8.64 -        Flow exc=(*_excess)[act];
    8.65 +        Value exc=(*_excess)[act];
    8.66  
    8.67          for(OutArcIt e(_g,act);e!=INVALID; ++e) {
    8.68            Node v = _g.target(e);
    8.69 -          Flow fc=(*_up)[e]-(*_flow)[e];
    8.70 +          Value fc=(*_up)[e]-(*_flow)[e];
    8.71            if(!_tol.positive(fc)) continue;
    8.72            if((*_level)[v]<actlevel) {
    8.73              if(!_tol.less(fc, exc)) {
    8.74 @@ -591,7 +591,7 @@
    8.75          }
    8.76          for(InArcIt e(_g,act);e!=INVALID; ++e) {
    8.77            Node v = _g.source(e);
    8.78 -          Flow fc=(*_flow)[e]-(*_lo)[e];
    8.79 +          Value fc=(*_flow)[e]-(*_lo)[e];
    8.80            if(!_tol.positive(fc)) continue;
    8.81            if((*_level)[v]<actlevel) {
    8.82              if(!_tol.less(fc, exc)) {
    8.83 @@ -661,13 +661,13 @@
    8.84  
    8.85      ///@{
    8.86  
    8.87 -    /// \brief Returns the flow on the given arc.
    8.88 +    /// \brief Returns the flow value on the given arc.
    8.89      ///
    8.90 -    /// Returns the flow on the given arc.
    8.91 +    /// Returns the flow value on the given arc.
    8.92      ///
    8.93      /// \pre Either \ref run() or \ref init() must be called before
    8.94      /// using this function.
    8.95 -    Flow flow(const Arc& arc) const {
    8.96 +    Value flow(const Arc& arc) const {
    8.97        return (*_flow)[arc];
    8.98      }
    8.99  
   8.100 @@ -750,7 +750,7 @@
   8.101          if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
   8.102        for(NodeIt n(_g);n!=INVALID;++n)
   8.103          {
   8.104 -          Flow dif=-(*_supply)[n];
   8.105 +          Value dif=-(*_supply)[n];
   8.106            for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
   8.107            for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
   8.108            if(_tol.negative(dif)) return false;
   8.109 @@ -765,10 +765,10 @@
   8.110      ///\sa barrierMap()
   8.111      bool checkBarrier() const
   8.112      {
   8.113 -      Flow delta=0;
   8.114 -      Flow inf_cap = std::numeric_limits<Flow>::has_infinity ?
   8.115 -        std::numeric_limits<Flow>::infinity() :
   8.116 -        std::numeric_limits<Flow>::max();
   8.117 +      Value delta=0;
   8.118 +      Value inf_cap = std::numeric_limits<Value>::has_infinity ?
   8.119 +        std::numeric_limits<Value>::infinity() :
   8.120 +        std::numeric_limits<Value>::max();
   8.121        for(NodeIt n(_g);n!=INVALID;++n)
   8.122          if(barrier(n))
   8.123            delta-=(*_supply)[n];
     9.1 --- a/lemon/core.h	Wed Apr 29 16:15:29 2009 +0100
     9.2 +++ b/lemon/core.h	Wed Apr 29 19:22:14 2009 +0100
     9.3 @@ -22,7 +22,7 @@
     9.4  #include <vector>
     9.5  #include <algorithm>
     9.6  
     9.7 -#include <lemon/core.h>
     9.8 +#include <lemon/config.h>
     9.9  #include <lemon/bits/enable_if.h>
    9.10  #include <lemon/bits/traits.h>
    9.11  #include <lemon/assert.h>
    10.1 --- a/lemon/network_simplex.h	Wed Apr 29 16:15:29 2009 +0100
    10.2 +++ b/lemon/network_simplex.h	Wed Apr 29 19:22:14 2009 +0100
    10.3 @@ -30,9 +30,6 @@
    10.4  
    10.5  #include <lemon/core.h>
    10.6  #include <lemon/math.h>
    10.7 -#include <lemon/maps.h>
    10.8 -#include <lemon/circulation.h>
    10.9 -#include <lemon/adaptors.h>
   10.10  
   10.11  namespace lemon {
   10.12  
   10.13 @@ -50,14 +47,19 @@
   10.14    ///
   10.15    /// In general this class is the fastest implementation available
   10.16    /// in LEMON for the minimum cost flow problem.
   10.17 -  /// Moreover it supports both direction of the supply/demand inequality
   10.18 -  /// constraints. For more information see \ref ProblemType.
   10.19 +  /// Moreover it supports both directions of the supply/demand inequality
   10.20 +  /// constraints. For more information see \ref SupplyType.
   10.21 +  ///
   10.22 +  /// Most of the parameters of the problem (except for the digraph)
   10.23 +  /// can be given using separate functions, and the algorithm can be
   10.24 +  /// executed using the \ref run() function. If some parameters are not
   10.25 +  /// specified, then default values will be used.
   10.26    ///
   10.27    /// \tparam GR The digraph type the algorithm runs on.
   10.28 -  /// \tparam F The value type used for flow amounts, capacity bounds
   10.29 +  /// \tparam V The value type used for flow amounts, capacity bounds
   10.30    /// and supply values in the algorithm. By default it is \c int.
   10.31    /// \tparam C The value type used for costs and potentials in the
   10.32 -  /// algorithm. By default it is the same as \c F.
   10.33 +  /// algorithm. By default it is the same as \c V.
   10.34    ///
   10.35    /// \warning Both value types must be signed and all input data must
   10.36    /// be integer.
   10.37 @@ -65,34 +67,92 @@
   10.38    /// \note %NetworkSimplex provides five different pivot rule
   10.39    /// implementations, from which the most efficient one is used
   10.40    /// by default. For more information see \ref PivotRule.
   10.41 -  template <typename GR, typename F = int, typename C = F>
   10.42 +  template <typename GR, typename V = int, typename C = V>
   10.43    class NetworkSimplex
   10.44    {
   10.45    public:
   10.46  
   10.47 -    /// The flow type of the algorithm
   10.48 -    typedef F Flow;
   10.49 -    /// The cost type of the algorithm
   10.50 +    /// The type of the flow amounts, capacity bounds and supply values
   10.51 +    typedef V Value;
   10.52 +    /// The type of the arc costs
   10.53      typedef C Cost;
   10.54 -#ifdef DOXYGEN
   10.55 -    /// The type of the flow map
   10.56 -    typedef GR::ArcMap<Flow> FlowMap;
   10.57 -    /// The type of the potential map
   10.58 -    typedef GR::NodeMap<Cost> PotentialMap;
   10.59 -#else
   10.60 -    /// The type of the flow map
   10.61 -    typedef typename GR::template ArcMap<Flow> FlowMap;
   10.62 -    /// The type of the potential map
   10.63 -    typedef typename GR::template NodeMap<Cost> PotentialMap;
   10.64 -#endif
   10.65  
   10.66    public:
   10.67  
   10.68 -    /// \brief Enum type for selecting the pivot rule.
   10.69 +    /// \brief Problem type constants for the \c run() function.
   10.70      ///
   10.71 -    /// Enum type for selecting the pivot rule for the \ref run()
   10.72 +    /// Enum type containing the problem type constants that can be
   10.73 +    /// returned by the \ref run() function of the algorithm.
   10.74 +    enum ProblemType {
   10.75 +      /// The problem has no feasible solution (flow).
   10.76 +      INFEASIBLE,
   10.77 +      /// The problem has optimal solution (i.e. it is feasible and
   10.78 +      /// bounded), and the algorithm has found optimal flow and node
   10.79 +      /// potentials (primal and dual solutions).
   10.80 +      OPTIMAL,
   10.81 +      /// The objective function of the problem is unbounded, i.e.
   10.82 +      /// there is a directed cycle having negative total cost and
   10.83 +      /// infinite upper bound.
   10.84 +      UNBOUNDED
   10.85 +    };
   10.86 +    
   10.87 +    /// \brief Constants for selecting the type of the supply constraints.
   10.88 +    ///
   10.89 +    /// Enum type containing constants for selecting the supply type,
   10.90 +    /// i.e. the direction of the inequalities in the supply/demand
   10.91 +    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
   10.92 +    ///
   10.93 +    /// The default supply type is \c GEQ, since this form is supported
   10.94 +    /// by other minimum cost flow algorithms and the \ref Circulation
   10.95 +    /// algorithm, as well.
   10.96 +    /// The \c LEQ problem type can be selected using the \ref supplyType()
   10.97      /// function.
   10.98      ///
   10.99 +    /// Note that the equality form is a special case of both supply types.
  10.100 +    enum SupplyType {
  10.101 +
  10.102 +      /// This option means that there are <em>"greater or equal"</em>
  10.103 +      /// supply/demand constraints in the definition, i.e. the exact
  10.104 +      /// formulation of the problem is the following.
  10.105 +      /**
  10.106 +          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
  10.107 +          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
  10.108 +              sup(u) \quad \forall u\in V \f]
  10.109 +          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
  10.110 +      */
  10.111 +      /// It means that the total demand must be greater or equal to the 
  10.112 +      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
  10.113 +      /// negative) and all the supplies have to be carried out from 
  10.114 +      /// the supply nodes, but there could be demands that are not 
  10.115 +      /// satisfied.
  10.116 +      GEQ,
  10.117 +      /// It is just an alias for the \c GEQ option.
  10.118 +      CARRY_SUPPLIES = GEQ,
  10.119 +
  10.120 +      /// This option means that there are <em>"less or equal"</em>
  10.121 +      /// supply/demand constraints in the definition, i.e. the exact
  10.122 +      /// formulation of the problem is the following.
  10.123 +      /**
  10.124 +          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
  10.125 +          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
  10.126 +              sup(u) \quad \forall u\in V \f]
  10.127 +          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
  10.128 +      */
  10.129 +      /// It means that the total demand must be less or equal to the 
  10.130 +      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
  10.131 +      /// positive) and all the demands have to be satisfied, but there
  10.132 +      /// could be supplies that are not carried out from the supply
  10.133 +      /// nodes.
  10.134 +      LEQ,
  10.135 +      /// It is just an alias for the \c LEQ option.
  10.136 +      SATISFY_DEMANDS = LEQ
  10.137 +    };
  10.138 +    
  10.139 +    /// \brief Constants for selecting the pivot rule.
  10.140 +    ///
  10.141 +    /// Enum type containing constants for selecting the pivot rule for
  10.142 +    /// the \ref run() function.
  10.143 +    ///
  10.144      /// \ref NetworkSimplex provides five different pivot rule
  10.145      /// implementations that significantly affect the running time
  10.146      /// of the algorithm.
  10.147 @@ -131,71 +191,15 @@
  10.148        ALTERING_LIST
  10.149      };
  10.150      
  10.151 -    /// \brief Enum type for selecting the problem type.
  10.152 -    ///
  10.153 -    /// Enum type for selecting the problem type, i.e. the direction of
  10.154 -    /// the inequalities in the supply/demand constraints of the
  10.155 -    /// \ref min_cost_flow "minimum cost flow problem".
  10.156 -    ///
  10.157 -    /// The default problem type is \c GEQ, since this form is supported
  10.158 -    /// by other minimum cost flow algorithms and the \ref Circulation
  10.159 -    /// algorithm as well.
  10.160 -    /// The \c LEQ problem type can be selected using the \ref problemType()
  10.161 -    /// function.
  10.162 -    ///
  10.163 -    /// Note that the equality form is a special case of both problem type.
  10.164 -    enum ProblemType {
  10.165 -
  10.166 -      /// This option means that there are "<em>greater or equal</em>"
  10.167 -      /// constraints in the defintion, i.e. the exact formulation of the
  10.168 -      /// problem is the following.
  10.169 -      /**
  10.170 -          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
  10.171 -          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
  10.172 -              sup(u) \quad \forall u\in V \f]
  10.173 -          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
  10.174 -      */
  10.175 -      /// It means that the total demand must be greater or equal to the 
  10.176 -      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
  10.177 -      /// negative) and all the supplies have to be carried out from 
  10.178 -      /// the supply nodes, but there could be demands that are not 
  10.179 -      /// satisfied.
  10.180 -      GEQ,
  10.181 -      /// It is just an alias for the \c GEQ option.
  10.182 -      CARRY_SUPPLIES = GEQ,
  10.183 -
  10.184 -      /// This option means that there are "<em>less or equal</em>"
  10.185 -      /// constraints in the defintion, i.e. the exact formulation of the
  10.186 -      /// problem is the following.
  10.187 -      /**
  10.188 -          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
  10.189 -          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
  10.190 -              sup(u) \quad \forall u\in V \f]
  10.191 -          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
  10.192 -      */
  10.193 -      /// It means that the total demand must be less or equal to the 
  10.194 -      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
  10.195 -      /// positive) and all the demands have to be satisfied, but there
  10.196 -      /// could be supplies that are not carried out from the supply
  10.197 -      /// nodes.
  10.198 -      LEQ,
  10.199 -      /// It is just an alias for the \c LEQ option.
  10.200 -      SATISFY_DEMANDS = LEQ
  10.201 -    };
  10.202 -
  10.203    private:
  10.204  
  10.205      TEMPLATE_DIGRAPH_TYPEDEFS(GR);
  10.206  
  10.207 -    typedef typename GR::template ArcMap<Flow> FlowArcMap;
  10.208 -    typedef typename GR::template ArcMap<Cost> CostArcMap;
  10.209 -    typedef typename GR::template NodeMap<Flow> FlowNodeMap;
  10.210 -
  10.211      typedef std::vector<Arc> ArcVector;
  10.212      typedef std::vector<Node> NodeVector;
  10.213      typedef std::vector<int> IntVector;
  10.214      typedef std::vector<bool> BoolVector;
  10.215 -    typedef std::vector<Flow> FlowVector;
  10.216 +    typedef std::vector<Value> ValueVector;
  10.217      typedef std::vector<Cost> CostVector;
  10.218  
  10.219      // State constants for arcs
  10.220 @@ -213,32 +217,23 @@
  10.221      int _arc_num;
  10.222  
  10.223      // Parameters of the problem
  10.224 -    FlowArcMap *_plower;
  10.225 -    FlowArcMap *_pupper;
  10.226 -    CostArcMap *_pcost;
  10.227 -    FlowNodeMap *_psupply;
  10.228 -    bool _pstsup;
  10.229 -    Node _psource, _ptarget;
  10.230 -    Flow _pstflow;
  10.231 -    ProblemType _ptype;
  10.232 -
  10.233 -    // Result maps
  10.234 -    FlowMap *_flow_map;
  10.235 -    PotentialMap *_potential_map;
  10.236 -    bool _local_flow;
  10.237 -    bool _local_potential;
  10.238 +    bool _have_lower;
  10.239 +    SupplyType _stype;
  10.240 +    Value _sum_supply;
  10.241  
  10.242      // Data structures for storing the digraph
  10.243      IntNodeMap _node_id;
  10.244 -    ArcVector _arc_ref;
  10.245 +    IntArcMap _arc_id;
  10.246      IntVector _source;
  10.247      IntVector _target;
  10.248  
  10.249      // Node and arc data
  10.250 -    FlowVector _cap;
  10.251 +    ValueVector _lower;
  10.252 +    ValueVector _upper;
  10.253 +    ValueVector _cap;
  10.254      CostVector _cost;
  10.255 -    FlowVector _supply;
  10.256 -    FlowVector _flow;
  10.257 +    ValueVector _supply;
  10.258 +    ValueVector _flow;
  10.259      CostVector _pi;
  10.260  
  10.261      // Data for storing the spanning tree structure
  10.262 @@ -257,7 +252,16 @@
  10.263      int in_arc, join, u_in, v_in, u_out, v_out;
  10.264      int first, second, right, last;
  10.265      int stem, par_stem, new_stem;
  10.266 -    Flow delta;
  10.267 +    Value delta;
  10.268 +
  10.269 +  public:
  10.270 +  
  10.271 +    /// \brief Constant for infinite upper bounds (capacities).
  10.272 +    ///
  10.273 +    /// Constant for infinite upper bounds (capacities).
  10.274 +    /// It is \c std::numeric_limits<Value>::infinity() if available,
  10.275 +    /// \c std::numeric_limits<Value>::max() otherwise.
  10.276 +    const Value INF;
  10.277  
  10.278    private:
  10.279  
  10.280 @@ -659,25 +663,68 @@
  10.281      ///
  10.282      /// \param graph The digraph the algorithm runs on.
  10.283      NetworkSimplex(const GR& graph) :
  10.284 -      _graph(graph),
  10.285 -      _plower(NULL), _pupper(NULL), _pcost(NULL),
  10.286 -      _psupply(NULL), _pstsup(false), _ptype(GEQ),
  10.287 -      _flow_map(NULL), _potential_map(NULL),
  10.288 -      _local_flow(false), _local_potential(false),
  10.289 -      _node_id(graph)
  10.290 +      _graph(graph), _node_id(graph), _arc_id(graph),
  10.291 +      INF(std::numeric_limits<Value>::has_infinity ?
  10.292 +          std::numeric_limits<Value>::infinity() :
  10.293 +          std::numeric_limits<Value>::max())
  10.294      {
  10.295 -      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
  10.296 -                   std::numeric_limits<Flow>::is_signed,
  10.297 -        "The flow type of NetworkSimplex must be signed integer");
  10.298 -      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
  10.299 -                   std::numeric_limits<Cost>::is_signed,
  10.300 -        "The cost type of NetworkSimplex must be signed integer");
  10.301 -    }
  10.302 +      // Check the value types
  10.303 +      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
  10.304 +        "The flow type of NetworkSimplex must be signed");
  10.305 +      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
  10.306 +        "The cost type of NetworkSimplex must be signed");
  10.307 +        
  10.308 +      // Resize vectors
  10.309 +      _node_num = countNodes(_graph);
  10.310 +      _arc_num = countArcs(_graph);
  10.311 +      int all_node_num = _node_num + 1;
  10.312 +      int all_arc_num = _arc_num + _node_num;
  10.313  
  10.314 -    /// Destructor.
  10.315 -    ~NetworkSimplex() {
  10.316 -      if (_local_flow) delete _flow_map;
  10.317 -      if (_local_potential) delete _potential_map;
  10.318 +      _source.resize(all_arc_num);
  10.319 +      _target.resize(all_arc_num);
  10.320 +
  10.321 +      _lower.resize(all_arc_num);
  10.322 +      _upper.resize(all_arc_num);
  10.323 +      _cap.resize(all_arc_num);
  10.324 +      _cost.resize(all_arc_num);
  10.325 +      _supply.resize(all_node_num);
  10.326 +      _flow.resize(all_arc_num);
  10.327 +      _pi.resize(all_node_num);
  10.328 +
  10.329 +      _parent.resize(all_node_num);
  10.330 +      _pred.resize(all_node_num);
  10.331 +      _forward.resize(all_node_num);
  10.332 +      _thread.resize(all_node_num);
  10.333 +      _rev_thread.resize(all_node_num);
  10.334 +      _succ_num.resize(all_node_num);
  10.335 +      _last_succ.resize(all_node_num);
  10.336 +      _state.resize(all_arc_num);
  10.337 +
  10.338 +      // Copy the graph (store the arcs in a mixed order)
  10.339 +      int i = 0;
  10.340 +      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  10.341 +        _node_id[n] = i;
  10.342 +      }
  10.343 +      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
  10.344 +      i = 0;
  10.345 +      for (ArcIt a(_graph); a != INVALID; ++a) {
  10.346 +        _arc_id[a] = i;
  10.347 +        _source[i] = _node_id[_graph.source(a)];
  10.348 +        _target[i] = _node_id[_graph.target(a)];
  10.349 +        if ((i += k) >= _arc_num) i = (i % k) + 1;
  10.350 +      }
  10.351 +      
  10.352 +      // Initialize maps
  10.353 +      for (int i = 0; i != _node_num; ++i) {
  10.354 +        _supply[i] = 0;
  10.355 +      }
  10.356 +      for (int i = 0; i != _arc_num; ++i) {
  10.357 +        _lower[i] = 0;
  10.358 +        _upper[i] = INF;
  10.359 +        _cost[i] = 1;
  10.360 +      }
  10.361 +      _have_lower = false;
  10.362 +      _stype = GEQ;
  10.363      }
  10.364  
  10.365      /// \name Parameters
  10.366 @@ -689,21 +736,19 @@
  10.367      /// \brief Set the lower bounds on the arcs.
  10.368      ///
  10.369      /// This function sets the lower bounds on the arcs.
  10.370 -    /// If neither this function nor \ref boundMaps() is used before
  10.371 -    /// calling \ref run(), the lower bounds will be set to zero
  10.372 -    /// on all arcs.
  10.373 +    /// If it is not used before calling \ref run(), the lower bounds
  10.374 +    /// will be set to zero on all arcs.
  10.375      ///
  10.376      /// \param map An arc map storing the lower bounds.
  10.377 -    /// Its \c Value type must be convertible to the \c Flow type
  10.378 +    /// Its \c Value type must be convertible to the \c Value type
  10.379      /// of the algorithm.
  10.380      ///
  10.381      /// \return <tt>(*this)</tt>
  10.382 -    template <typename LOWER>
  10.383 -    NetworkSimplex& lowerMap(const LOWER& map) {
  10.384 -      delete _plower;
  10.385 -      _plower = new FlowArcMap(_graph);
  10.386 +    template <typename LowerMap>
  10.387 +    NetworkSimplex& lowerMap(const LowerMap& map) {
  10.388 +      _have_lower = true;
  10.389        for (ArcIt a(_graph); a != INVALID; ++a) {
  10.390 -        (*_plower)[a] = map[a];
  10.391 +        _lower[_arc_id[a]] = map[a];
  10.392        }
  10.393        return *this;
  10.394      }
  10.395 @@ -711,63 +756,23 @@
  10.396      /// \brief Set the upper bounds (capacities) on the arcs.
  10.397      ///
  10.398      /// This function sets the upper bounds (capacities) on the arcs.
  10.399 -    /// If none of the functions \ref upperMap(), \ref capacityMap()
  10.400 -    /// and \ref boundMaps() is used before calling \ref run(),
  10.401 -    /// the upper bounds (capacities) will be set to
  10.402 -    /// \c std::numeric_limits<Flow>::max() on all arcs.
  10.403 +    /// If it is not used before calling \ref run(), the upper bounds
  10.404 +    /// will be set to \ref INF on all arcs (i.e. the flow value will be
  10.405 +    /// unbounded from above on each arc).
  10.406      ///
  10.407      /// \param map An arc map storing the upper bounds.
  10.408 -    /// Its \c Value type must be convertible to the \c Flow type
  10.409 +    /// Its \c Value type must be convertible to the \c Value type
  10.410      /// of the algorithm.
  10.411      ///
  10.412      /// \return <tt>(*this)</tt>
  10.413 -    template<typename UPPER>
  10.414 -    NetworkSimplex& upperMap(const UPPER& map) {
  10.415 -      delete _pupper;
  10.416 -      _pupper = new FlowArcMap(_graph);
  10.417 +    template<typename UpperMap>
  10.418 +    NetworkSimplex& upperMap(const UpperMap& map) {
  10.419        for (ArcIt a(_graph); a != INVALID; ++a) {
  10.420 -        (*_pupper)[a] = map[a];
  10.421 +        _upper[_arc_id[a]] = map[a];
  10.422        }
  10.423        return *this;
  10.424      }
  10.425  
  10.426 -    /// \brief Set the upper bounds (capacities) on the arcs.
  10.427 -    ///
  10.428 -    /// This function sets the upper bounds (capacities) on the arcs.
  10.429 -    /// It is just an alias for \ref upperMap().
  10.430 -    ///
  10.431 -    /// \return <tt>(*this)</tt>
  10.432 -    template<typename CAP>
  10.433 -    NetworkSimplex& capacityMap(const CAP& map) {
  10.434 -      return upperMap(map);
  10.435 -    }
  10.436 -
  10.437 -    /// \brief Set the lower and upper bounds on the arcs.
  10.438 -    ///
  10.439 -    /// This function sets the lower and upper bounds on the arcs.
  10.440 -    /// If neither this function nor \ref lowerMap() is used before
  10.441 -    /// calling \ref run(), the lower bounds will be set to zero
  10.442 -    /// on all arcs.
  10.443 -    /// If none of the functions \ref upperMap(), \ref capacityMap()
  10.444 -    /// and \ref boundMaps() is used before calling \ref run(),
  10.445 -    /// the upper bounds (capacities) will be set to
  10.446 -    /// \c std::numeric_limits<Flow>::max() on all arcs.
  10.447 -    ///
  10.448 -    /// \param lower An arc map storing the lower bounds.
  10.449 -    /// \param upper An arc map storing the upper bounds.
  10.450 -    ///
  10.451 -    /// The \c Value type of the maps must be convertible to the
  10.452 -    /// \c Flow type of the algorithm.
  10.453 -    ///
  10.454 -    /// \note This function is just a shortcut of calling \ref lowerMap()
  10.455 -    /// and \ref upperMap() separately.
  10.456 -    ///
  10.457 -    /// \return <tt>(*this)</tt>
  10.458 -    template <typename LOWER, typename UPPER>
  10.459 -    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
  10.460 -      return lowerMap(lower).upperMap(upper);
  10.461 -    }
  10.462 -
  10.463      /// \brief Set the costs of the arcs.
  10.464      ///
  10.465      /// This function sets the costs of the arcs.
  10.466 @@ -779,12 +784,10 @@
  10.467      /// of the algorithm.
  10.468      ///
  10.469      /// \return <tt>(*this)</tt>
  10.470 -    template<typename COST>
  10.471 -    NetworkSimplex& costMap(const COST& map) {
  10.472 -      delete _pcost;
  10.473 -      _pcost = new CostArcMap(_graph);
  10.474 +    template<typename CostMap>
  10.475 +    NetworkSimplex& costMap(const CostMap& map) {
  10.476        for (ArcIt a(_graph); a != INVALID; ++a) {
  10.477 -        (*_pcost)[a] = map[a];
  10.478 +        _cost[_arc_id[a]] = map[a];
  10.479        }
  10.480        return *this;
  10.481      }
  10.482 @@ -797,17 +800,14 @@
  10.483      /// (It makes sense only if non-zero lower bounds are given.)
  10.484      ///
  10.485      /// \param map A node map storing the supply values.
  10.486 -    /// Its \c Value type must be convertible to the \c Flow type
  10.487 +    /// Its \c Value type must be convertible to the \c Value type
  10.488      /// of the algorithm.
  10.489      ///
  10.490      /// \return <tt>(*this)</tt>
  10.491 -    template<typename SUP>
  10.492 -    NetworkSimplex& supplyMap(const SUP& map) {
  10.493 -      delete _psupply;
  10.494 -      _pstsup = false;
  10.495 -      _psupply = new FlowNodeMap(_graph);
  10.496 +    template<typename SupplyMap>
  10.497 +    NetworkSimplex& supplyMap(const SupplyMap& map) {
  10.498        for (NodeIt n(_graph); n != INVALID; ++n) {
  10.499 -        (*_psupply)[n] = map[n];
  10.500 +        _supply[_node_id[n]] = map[n];
  10.501        }
  10.502        return *this;
  10.503      }
  10.504 @@ -820,71 +820,39 @@
  10.505      /// calling \ref run(), the supply of each node will be set to zero.
  10.506      /// (It makes sense only if non-zero lower bounds are given.)
  10.507      ///
  10.508 +    /// Using this function has the same effect as using \ref supplyMap()
  10.509 +    /// with such a map in which \c k is assigned to \c s, \c -k is
  10.510 +    /// assigned to \c t and all other nodes have zero supply value.
  10.511 +    ///
  10.512      /// \param s The source node.
  10.513      /// \param t The target node.
  10.514      /// \param k The required amount of flow from node \c s to node \c t
  10.515      /// (i.e. the supply of \c s and the demand of \c t).
  10.516      ///
  10.517      /// \return <tt>(*this)</tt>
  10.518 -    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
  10.519 -      delete _psupply;
  10.520 -      _psupply = NULL;
  10.521 -      _pstsup = true;
  10.522 -      _psource = s;
  10.523 -      _ptarget = t;
  10.524 -      _pstflow = k;
  10.525 +    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
  10.526 +      for (int i = 0; i != _node_num; ++i) {
  10.527 +        _supply[i] = 0;
  10.528 +      }
  10.529 +      _supply[_node_id[s]] =  k;
  10.530 +      _supply[_node_id[t]] = -k;
  10.531        return *this;
  10.532      }
  10.533      
  10.534 -    /// \brief Set the problem type.
  10.535 +    /// \brief Set the type of the supply constraints.
  10.536      ///
  10.537 -    /// This function sets the problem type for the algorithm.
  10.538 -    /// If it is not used before calling \ref run(), the \ref GEQ problem
  10.539 +    /// This function sets the type of the supply/demand constraints.
  10.540 +    /// If it is not used before calling \ref run(), the \ref GEQ supply
  10.541      /// type will be used.
  10.542      ///
  10.543 -    /// For more information see \ref ProblemType.
  10.544 +    /// For more information see \ref SupplyType.
  10.545      ///
  10.546      /// \return <tt>(*this)</tt>
  10.547 -    NetworkSimplex& problemType(ProblemType problem_type) {
  10.548 -      _ptype = problem_type;
  10.549 +    NetworkSimplex& supplyType(SupplyType supply_type) {
  10.550 +      _stype = supply_type;
  10.551        return *this;
  10.552      }
  10.553  
  10.554 -    /// \brief Set the flow map.
  10.555 -    ///
  10.556 -    /// This function sets the flow map.
  10.557 -    /// If it is not used before calling \ref run(), an instance will
  10.558 -    /// be allocated automatically. The destructor deallocates this
  10.559 -    /// automatically allocated map, of course.
  10.560 -    ///
  10.561 -    /// \return <tt>(*this)</tt>
  10.562 -    NetworkSimplex& flowMap(FlowMap& map) {
  10.563 -      if (_local_flow) {
  10.564 -        delete _flow_map;
  10.565 -        _local_flow = false;
  10.566 -      }
  10.567 -      _flow_map = &map;
  10.568 -      return *this;
  10.569 -    }
  10.570 -
  10.571 -    /// \brief Set the potential map.
  10.572 -    ///
  10.573 -    /// This function sets the potential map, which is used for storing
  10.574 -    /// the dual solution.
  10.575 -    /// If it is not used before calling \ref run(), an instance will
  10.576 -    /// be allocated automatically. The destructor deallocates this
  10.577 -    /// automatically allocated map, of course.
  10.578 -    ///
  10.579 -    /// \return <tt>(*this)</tt>
  10.580 -    NetworkSimplex& potentialMap(PotentialMap& map) {
  10.581 -      if (_local_potential) {
  10.582 -        delete _potential_map;
  10.583 -        _local_potential = false;
  10.584 -      }
  10.585 -      _potential_map = &map;
  10.586 -      return *this;
  10.587 -    }
  10.588 -    
  10.589      /// @}
  10.590  
  10.591      /// \name Execution Control
  10.592 @@ -896,13 +864,12 @@
  10.593      ///
  10.594      /// This function runs the algorithm.
  10.595      /// The paramters can be specified using functions \ref lowerMap(),
  10.596 -    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
  10.597 -    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
  10.598 -    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
  10.599 +    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
  10.600 +    /// \ref supplyType().
  10.601      /// For example,
  10.602      /// \code
  10.603      ///   NetworkSimplex<ListDigraph> ns(graph);
  10.604 -    ///   ns.boundMaps(lower, upper).costMap(cost)
  10.605 +    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
  10.606      ///     .supplyMap(sup).run();
  10.607      /// \endcode
  10.608      ///
  10.609 @@ -910,33 +877,44 @@
  10.610      /// that have been given are kept for the next call, unless
  10.611      /// \ref reset() is called, thus only the modified parameters
  10.612      /// have to be set again. See \ref reset() for examples.
  10.613 +    /// However the underlying digraph must not be modified after this
  10.614 +    /// class have been constructed, since it copies and extends the graph.
  10.615      ///
  10.616      /// \param pivot_rule The pivot rule that will be used during the
  10.617      /// algorithm. For more information see \ref PivotRule.
  10.618      ///
  10.619 -    /// \return \c true if a feasible flow can be found.
  10.620 -    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
  10.621 -      return init() && start(pivot_rule);
  10.622 +    /// \return \c INFEASIBLE if no feasible flow exists,
  10.623 +    /// \n \c OPTIMAL if the problem has optimal solution
  10.624 +    /// (i.e. it is feasible and bounded), and the algorithm has found
  10.625 +    /// optimal flow and node potentials (primal and dual solutions),
  10.626 +    /// \n \c UNBOUNDED if the objective function of the problem is
  10.627 +    /// unbounded, i.e. there is a directed cycle having negative total
  10.628 +    /// cost and infinite upper bound.
  10.629 +    ///
  10.630 +    /// \see ProblemType, PivotRule
  10.631 +    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
  10.632 +      if (!init()) return INFEASIBLE;
  10.633 +      return start(pivot_rule);
  10.634      }
  10.635  
  10.636      /// \brief Reset all the parameters that have been given before.
  10.637      ///
  10.638      /// This function resets all the paramaters that have been given
  10.639      /// before using functions \ref lowerMap(), \ref upperMap(),
  10.640 -    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
  10.641 -    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
  10.642 -    /// \ref flowMap() and \ref potentialMap().
  10.643 +    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
  10.644      ///
  10.645      /// It is useful for multiple run() calls. If this function is not
  10.646      /// used, all the parameters given before are kept for the next
  10.647      /// \ref run() call.
  10.648 +    /// However the underlying digraph must not be modified after this
  10.649 +    /// class have been constructed, since it copies and extends the graph.
  10.650      ///
  10.651      /// For example,
  10.652      /// \code
  10.653      ///   NetworkSimplex<ListDigraph> ns(graph);
  10.654      ///
  10.655      ///   // First run
  10.656 -    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
  10.657 +    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
  10.658      ///     .supplyMap(sup).run();
  10.659      ///
  10.660      ///   // Run again with modified cost map (reset() is not called,
  10.661 @@ -947,29 +925,22 @@
  10.662      ///   // Run again from scratch using reset()
  10.663      ///   // (the lower bounds will be set to zero on all arcs)
  10.664      ///   ns.reset();
  10.665 -    ///   ns.capacityMap(cap).costMap(cost)
  10.666 +    ///   ns.upperMap(capacity).costMap(cost)
  10.667      ///     .supplyMap(sup).run();
  10.668      /// \endcode
  10.669      ///
  10.670      /// \return <tt>(*this)</tt>
  10.671      NetworkSimplex& reset() {
  10.672 -      delete _plower;
  10.673 -      delete _pupper;
  10.674 -      delete _pcost;
  10.675 -      delete _psupply;
  10.676 -      _plower = NULL;
  10.677 -      _pupper = NULL;
  10.678 -      _pcost = NULL;
  10.679 -      _psupply = NULL;
  10.680 -      _pstsup = false;
  10.681 -      _ptype = GEQ;
  10.682 -      if (_local_flow) delete _flow_map;
  10.683 -      if (_local_potential) delete _potential_map;
  10.684 -      _flow_map = NULL;
  10.685 -      _potential_map = NULL;
  10.686 -      _local_flow = false;
  10.687 -      _local_potential = false;
  10.688 -
  10.689 +      for (int i = 0; i != _node_num; ++i) {
  10.690 +        _supply[i] = 0;
  10.691 +      }
  10.692 +      for (int i = 0; i != _arc_num; ++i) {
  10.693 +        _lower[i] = 0;
  10.694 +        _upper[i] = INF;
  10.695 +        _cost[i] = 1;
  10.696 +      }
  10.697 +      _have_lower = false;
  10.698 +      _stype = GEQ;
  10.699        return *this;
  10.700      }
  10.701  
  10.702 @@ -985,7 +956,7 @@
  10.703      /// \brief Return the total cost of the found flow.
  10.704      ///
  10.705      /// This function returns the total cost of the found flow.
  10.706 -    /// The complexity of the function is O(e).
  10.707 +    /// Its complexity is O(e).
  10.708      ///
  10.709      /// \note The return type of the function can be specified as a
  10.710      /// template parameter. For example,
  10.711 @@ -997,15 +968,12 @@
  10.712      /// function.
  10.713      ///
  10.714      /// \pre \ref run() must be called before using this function.
  10.715 -    template <typename Num>
  10.716 -    Num totalCost() const {
  10.717 -      Num c = 0;
  10.718 -      if (_pcost) {
  10.719 -        for (ArcIt e(_graph); e != INVALID; ++e)
  10.720 -          c += (*_flow_map)[e] * (*_pcost)[e];
  10.721 -      } else {
  10.722 -        for (ArcIt e(_graph); e != INVALID; ++e)
  10.723 -          c += (*_flow_map)[e];
  10.724 +    template <typename Number>
  10.725 +    Number totalCost() const {
  10.726 +      Number c = 0;
  10.727 +      for (ArcIt a(_graph); a != INVALID; ++a) {
  10.728 +        int i = _arc_id[a];
  10.729 +        c += Number(_flow[i]) * Number(_cost[i]);
  10.730        }
  10.731        return c;
  10.732      }
  10.733 @@ -1021,18 +989,22 @@
  10.734      /// This function returns the flow on the given arc.
  10.735      ///
  10.736      /// \pre \ref run() must be called before using this function.
  10.737 -    Flow flow(const Arc& a) const {
  10.738 -      return (*_flow_map)[a];
  10.739 +    Value flow(const Arc& a) const {
  10.740 +      return _flow[_arc_id[a]];
  10.741      }
  10.742  
  10.743 -    /// \brief Return a const reference to the flow map.
  10.744 +    /// \brief Return the flow map (the primal solution).
  10.745      ///
  10.746 -    /// This function returns a const reference to an arc map storing
  10.747 -    /// the found flow.
  10.748 +    /// This function copies the flow value on each arc into the given
  10.749 +    /// map. The \c Value type of the algorithm must be convertible to
  10.750 +    /// the \c Value type of the map.
  10.751      ///
  10.752      /// \pre \ref run() must be called before using this function.
  10.753 -    const FlowMap& flowMap() const {
  10.754 -      return *_flow_map;
  10.755 +    template <typename FlowMap>
  10.756 +    void flowMap(FlowMap &map) const {
  10.757 +      for (ArcIt a(_graph); a != INVALID; ++a) {
  10.758 +        map.set(a, _flow[_arc_id[a]]);
  10.759 +      }
  10.760      }
  10.761  
  10.762      /// \brief Return the potential (dual value) of the given node.
  10.763 @@ -1042,19 +1014,22 @@
  10.764      ///
  10.765      /// \pre \ref run() must be called before using this function.
  10.766      Cost potential(const Node& n) const {
  10.767 -      return (*_potential_map)[n];
  10.768 +      return _pi[_node_id[n]];
  10.769      }
  10.770  
  10.771 -    /// \brief Return a const reference to the potential map
  10.772 -    /// (the dual solution).
  10.773 +    /// \brief Return the potential map (the dual solution).
  10.774      ///
  10.775 -    /// This function returns a const reference to a node map storing
  10.776 -    /// the found potentials, which form the dual solution of the
  10.777 -    /// \ref min_cost_flow "minimum cost flow" problem.
  10.778 +    /// This function copies the potential (dual value) of each node
  10.779 +    /// into the given map.
  10.780 +    /// The \c Cost type of the algorithm must be convertible to the
  10.781 +    /// \c Value type of the map.
  10.782      ///
  10.783      /// \pre \ref run() must be called before using this function.
  10.784 -    const PotentialMap& potentialMap() const {
  10.785 -      return *_potential_map;
  10.786 +    template <typename PotentialMap>
  10.787 +    void potentialMap(PotentialMap &map) const {
  10.788 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  10.789 +        map.set(n, _pi[_node_id[n]]);
  10.790 +      }
  10.791      }
  10.792  
  10.793      /// @}
  10.794 @@ -1063,245 +1038,82 @@
  10.795  
  10.796      // Initialize internal data structures
  10.797      bool init() {
  10.798 -      // Initialize result maps
  10.799 -      if (!_flow_map) {
  10.800 -        _flow_map = new FlowMap(_graph);
  10.801 -        _local_flow = true;
  10.802 +      if (_node_num == 0) return false;
  10.803 +
  10.804 +      // Check the sum of supply values
  10.805 +      _sum_supply = 0;
  10.806 +      for (int i = 0; i != _node_num; ++i) {
  10.807 +        _sum_supply += _supply[i];
  10.808        }
  10.809 -      if (!_potential_map) {
  10.810 -        _potential_map = new PotentialMap(_graph);
  10.811 -        _local_potential = true;
  10.812 +      if ( !((_stype == GEQ && _sum_supply <= 0) ||
  10.813 +             (_stype == LEQ && _sum_supply >= 0)) ) return false;
  10.814 +
  10.815 +      // Remove non-zero lower bounds
  10.816 +      if (_have_lower) {
  10.817 +        for (int i = 0; i != _arc_num; ++i) {
  10.818 +          Value c = _lower[i];
  10.819 +          if (c >= 0) {
  10.820 +            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
  10.821 +          } else {
  10.822 +            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
  10.823 +          }
  10.824 +          _supply[_source[i]] -= c;
  10.825 +          _supply[_target[i]] += c;
  10.826 +        }
  10.827 +      } else {
  10.828 +        for (int i = 0; i != _arc_num; ++i) {
  10.829 +          _cap[i] = _upper[i];
  10.830 +        }
  10.831        }
  10.832  
  10.833 -      // Initialize vectors
  10.834 -      _node_num = countNodes(_graph);
  10.835 -      _arc_num = countArcs(_graph);
  10.836 -      int all_node_num = _node_num + 1;
  10.837 -      int all_arc_num = _arc_num + _node_num;
  10.838 -      if (_node_num == 0) return false;
  10.839 -
  10.840 -      _arc_ref.resize(_arc_num);
  10.841 -      _source.resize(all_arc_num);
  10.842 -      _target.resize(all_arc_num);
  10.843 -
  10.844 -      _cap.resize(all_arc_num);
  10.845 -      _cost.resize(all_arc_num);
  10.846 -      _supply.resize(all_node_num);
  10.847 -      _flow.resize(all_arc_num);
  10.848 -      _pi.resize(all_node_num);
  10.849 -
  10.850 -      _parent.resize(all_node_num);
  10.851 -      _pred.resize(all_node_num);
  10.852 -      _forward.resize(all_node_num);
  10.853 -      _thread.resize(all_node_num);
  10.854 -      _rev_thread.resize(all_node_num);
  10.855 -      _succ_num.resize(all_node_num);
  10.856 -      _last_succ.resize(all_node_num);
  10.857 -      _state.resize(all_arc_num);
  10.858 -
  10.859 -      // Initialize node related data
  10.860 -      bool valid_supply = true;
  10.861 -      Flow sum_supply = 0;
  10.862 -      if (!_pstsup && !_psupply) {
  10.863 -        _pstsup = true;
  10.864 -        _psource = _ptarget = NodeIt(_graph);
  10.865 -        _pstflow = 0;
  10.866 -      }
  10.867 -      if (_psupply) {
  10.868 -        int i = 0;
  10.869 -        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  10.870 -          _node_id[n] = i;
  10.871 -          _supply[i] = (*_psupply)[n];
  10.872 -          sum_supply += _supply[i];
  10.873 +      // Initialize artifical cost
  10.874 +      Cost ART_COST;
  10.875 +      if (std::numeric_limits<Cost>::is_exact) {
  10.876 +        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
  10.877 +      } else {
  10.878 +        ART_COST = std::numeric_limits<Cost>::min();
  10.879 +        for (int i = 0; i != _arc_num; ++i) {
  10.880 +          if (_cost[i] > ART_COST) ART_COST = _cost[i];
  10.881          }
  10.882 -        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
  10.883 -                       (_ptype == LEQ && sum_supply >= 0);
  10.884 -      } else {
  10.885 -        int i = 0;
  10.886 -        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  10.887 -          _node_id[n] = i;
  10.888 -          _supply[i] = 0;
  10.889 -        }
  10.890 -        _supply[_node_id[_psource]] =  _pstflow;
  10.891 -        _supply[_node_id[_ptarget]] = -_pstflow;
  10.892 -      }
  10.893 -      if (!valid_supply) return false;
  10.894 -
  10.895 -      // Infinite capacity value
  10.896 -      Flow inf_cap =
  10.897 -        std::numeric_limits<Flow>::has_infinity ?
  10.898 -        std::numeric_limits<Flow>::infinity() :
  10.899 -        std::numeric_limits<Flow>::max();
  10.900 -
  10.901 -      // Initialize artifical cost
  10.902 -      Cost art_cost;
  10.903 -      if (std::numeric_limits<Cost>::is_exact) {
  10.904 -        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
  10.905 -      } else {
  10.906 -        art_cost = std::numeric_limits<Cost>::min();
  10.907 -        for (int i = 0; i != _arc_num; ++i) {
  10.908 -          if (_cost[i] > art_cost) art_cost = _cost[i];
  10.909 -        }
  10.910 -        art_cost = (art_cost + 1) * _node_num;
  10.911 +        ART_COST = (ART_COST + 1) * _node_num;
  10.912        }
  10.913  
  10.914 -      // Run Circulation to check if a feasible solution exists
  10.915 -      typedef ConstMap<Arc, Flow> ConstArcMap;
  10.916 -      ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
  10.917 -      FlowNodeMap *csup = NULL;
  10.918 -      bool local_csup = false;
  10.919 -      if (_psupply) {
  10.920 -        csup = _psupply;
  10.921 -      } else {
  10.922 -        csup = new FlowNodeMap(_graph, 0);
  10.923 -        (*csup)[_psource] =  _pstflow;
  10.924 -        (*csup)[_ptarget] = -_pstflow;
  10.925 -        local_csup = true;
  10.926 +      // Initialize arc maps
  10.927 +      for (int i = 0; i != _arc_num; ++i) {
  10.928 +        _flow[i] = 0;
  10.929 +        _state[i] = STATE_LOWER;
  10.930        }
  10.931 -      bool circ_result = false;
  10.932 -      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
  10.933 -        // GEQ problem type
  10.934 -        if (_plower) {
  10.935 -          if (_pupper) {
  10.936 -            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
  10.937 -              circ(_graph, *_plower, *_pupper, *csup);
  10.938 -            circ_result = circ.run();
  10.939 -          } else {
  10.940 -            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
  10.941 -              circ(_graph, *_plower, inf_arc_map, *csup);
  10.942 -            circ_result = circ.run();
  10.943 -          }
  10.944 -        } else {
  10.945 -          if (_pupper) {
  10.946 -            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
  10.947 -              circ(_graph, zero_arc_map, *_pupper, *csup);
  10.948 -            circ_result = circ.run();
  10.949 -          } else {
  10.950 -            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
  10.951 -              circ(_graph, zero_arc_map, inf_arc_map, *csup);
  10.952 -            circ_result = circ.run();
  10.953 -          }
  10.954 -        }
  10.955 -      } else {
  10.956 -        // LEQ problem type
  10.957 -        typedef ReverseDigraph<const GR> RevGraph;
  10.958 -        typedef NegMap<FlowNodeMap> NegNodeMap;
  10.959 -        RevGraph rgraph(_graph);
  10.960 -        NegNodeMap neg_csup(*csup);
  10.961 -        if (_plower) {
  10.962 -          if (_pupper) {
  10.963 -            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
  10.964 -              circ(rgraph, *_plower, *_pupper, neg_csup);
  10.965 -            circ_result = circ.run();
  10.966 -          } else {
  10.967 -            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
  10.968 -              circ(rgraph, *_plower, inf_arc_map, neg_csup);
  10.969 -            circ_result = circ.run();
  10.970 -          }
  10.971 -        } else {
  10.972 -          if (_pupper) {
  10.973 -            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
  10.974 -              circ(rgraph, zero_arc_map, *_pupper, neg_csup);
  10.975 -            circ_result = circ.run();
  10.976 -          } else {
  10.977 -            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
  10.978 -              circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
  10.979 -            circ_result = circ.run();
  10.980 -          }
  10.981 -        }
  10.982 -      }
  10.983 -      if (local_csup) delete csup;
  10.984 -      if (!circ_result) return false;
  10.985 -
  10.986 +      
  10.987        // Set data for the artificial root node
  10.988        _root = _node_num;
  10.989        _parent[_root] = -1;
  10.990        _pred[_root] = -1;
  10.991        _thread[_root] = 0;
  10.992        _rev_thread[0] = _root;
  10.993 -      _succ_num[_root] = all_node_num;
  10.994 +      _succ_num[_root] = _node_num + 1;
  10.995        _last_succ[_root] = _root - 1;
  10.996 -      _supply[_root] = -sum_supply;
  10.997 -      if (sum_supply < 0) {
  10.998 -        _pi[_root] = -art_cost;
  10.999 -      } else {
 10.1000 -        _pi[_root] = art_cost;
 10.1001 -      }
 10.1002 -
 10.1003 -      // Store the arcs in a mixed order
 10.1004 -      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
 10.1005 -      int i = 0;
 10.1006 -      for (ArcIt e(_graph); e != INVALID; ++e) {
 10.1007 -        _arc_ref[i] = e;
 10.1008 -        if ((i += k) >= _arc_num) i = (i % k) + 1;
 10.1009 -      }
 10.1010 -
 10.1011 -      // Initialize arc maps
 10.1012 -      if (_pupper && _pcost) {
 10.1013 -        for (int i = 0; i != _arc_num; ++i) {
 10.1014 -          Arc e = _arc_ref[i];
 10.1015 -          _source[i] = _node_id[_graph.source(e)];
 10.1016 -          _target[i] = _node_id[_graph.target(e)];
 10.1017 -          _cap[i] = (*_pupper)[e];
 10.1018 -          _cost[i] = (*_pcost)[e];
 10.1019 -          _flow[i] = 0;
 10.1020 -          _state[i] = STATE_LOWER;
 10.1021 -        }
 10.1022 -      } else {
 10.1023 -        for (int i = 0; i != _arc_num; ++i) {
 10.1024 -          Arc e = _arc_ref[i];
 10.1025 -          _source[i] = _node_id[_graph.source(e)];
 10.1026 -          _target[i] = _node_id[_graph.target(e)];
 10.1027 -          _flow[i] = 0;
 10.1028 -          _state[i] = STATE_LOWER;
 10.1029 -        }
 10.1030 -        if (_pupper) {
 10.1031 -          for (int i = 0; i != _arc_num; ++i)
 10.1032 -            _cap[i] = (*_pupper)[_arc_ref[i]];
 10.1033 -        } else {
 10.1034 -          for (int i = 0; i != _arc_num; ++i)
 10.1035 -            _cap[i] = inf_cap;
 10.1036 -        }
 10.1037 -        if (_pcost) {
 10.1038 -          for (int i = 0; i != _arc_num; ++i)
 10.1039 -            _cost[i] = (*_pcost)[_arc_ref[i]];
 10.1040 -        } else {
 10.1041 -          for (int i = 0; i != _arc_num; ++i)
 10.1042 -            _cost[i] = 1;
 10.1043 -        }
 10.1044 -      }
 10.1045 -      
 10.1046 -      // Remove non-zero lower bounds
 10.1047 -      if (_plower) {
 10.1048 -        for (int i = 0; i != _arc_num; ++i) {
 10.1049 -          Flow c = (*_plower)[_arc_ref[i]];
 10.1050 -          if (c != 0) {
 10.1051 -            _cap[i] -= c;
 10.1052 -            _supply[_source[i]] -= c;
 10.1053 -            _supply[_target[i]] += c;
 10.1054 -          }
 10.1055 -        }
 10.1056 -      }
 10.1057 +      _supply[_root] = -_sum_supply;
 10.1058 +      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
 10.1059  
 10.1060        // Add artificial arcs and initialize the spanning tree data structure
 10.1061        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 10.1062 +        _parent[u] = _root;
 10.1063 +        _pred[u] = e;
 10.1064          _thread[u] = u + 1;
 10.1065          _rev_thread[u + 1] = u;
 10.1066          _succ_num[u] = 1;
 10.1067          _last_succ[u] = u;
 10.1068 -        _parent[u] = _root;
 10.1069 -        _pred[u] = e;
 10.1070 -        _cost[e] = art_cost;
 10.1071 -        _cap[e] = inf_cap;
 10.1072 +        _cost[e] = ART_COST;
 10.1073 +        _cap[e] = INF;
 10.1074          _state[e] = STATE_TREE;
 10.1075 -        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
 10.1076 +        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
 10.1077            _flow[e] = _supply[u];
 10.1078            _forward[u] = true;
 10.1079 -          _pi[u] = -art_cost + _pi[_root];
 10.1080 +          _pi[u] = -ART_COST + _pi[_root];
 10.1081          } else {
 10.1082            _flow[e] = -_supply[u];
 10.1083            _forward[u] = false;
 10.1084 -          _pi[u] = art_cost + _pi[_root];
 10.1085 +          _pi[u] = ART_COST + _pi[_root];
 10.1086          }
 10.1087        }
 10.1088  
 10.1089 @@ -1336,13 +1148,14 @@
 10.1090        }
 10.1091        delta = _cap[in_arc];
 10.1092        int result = 0;
 10.1093 -      Flow d;
 10.1094 +      Value d;
 10.1095        int e;
 10.1096  
 10.1097        // Search the cycle along the path form the first node to the root
 10.1098        for (int u = first; u != join; u = _parent[u]) {
 10.1099          e = _pred[u];
 10.1100 -        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
 10.1101 +        d = _forward[u] ?
 10.1102 +          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
 10.1103          if (d < delta) {
 10.1104            delta = d;
 10.1105            u_out = u;
 10.1106 @@ -1352,7 +1165,8 @@
 10.1107        // Search the cycle along the path form the second node to the root
 10.1108        for (int u = second; u != join; u = _parent[u]) {
 10.1109          e = _pred[u];
 10.1110 -        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
 10.1111 +        d = _forward[u] ? 
 10.1112 +          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
 10.1113          if (d <= delta) {
 10.1114            delta = d;
 10.1115            u_out = u;
 10.1116 @@ -1374,7 +1188,7 @@
 10.1117      void changeFlow(bool change) {
 10.1118        // Augment along the cycle
 10.1119        if (delta > 0) {
 10.1120 -        Flow val = _state[in_arc] * delta;
 10.1121 +        Value val = _state[in_arc] * delta;
 10.1122          _flow[in_arc] += val;
 10.1123          for (int u = _source[in_arc]; u != join; u = _parent[u]) {
 10.1124            _flow[_pred[u]] += _forward[u] ? -val : val;
 10.1125 @@ -1526,7 +1340,7 @@
 10.1126      }
 10.1127  
 10.1128      // Execute the algorithm
 10.1129 -    bool start(PivotRule pivot_rule) {
 10.1130 +    ProblemType start(PivotRule pivot_rule) {
 10.1131        // Select the pivot rule implementation
 10.1132        switch (pivot_rule) {
 10.1133          case FIRST_ELIGIBLE:
 10.1134 @@ -1540,41 +1354,55 @@
 10.1135          case ALTERING_LIST:
 10.1136            return start<AlteringListPivotRule>();
 10.1137        }
 10.1138 -      return false;
 10.1139 +      return INFEASIBLE; // avoid warning
 10.1140      }
 10.1141  
 10.1142      template <typename PivotRuleImpl>
 10.1143 -    bool start() {
 10.1144 +    ProblemType start() {
 10.1145        PivotRuleImpl pivot(*this);
 10.1146  
 10.1147        // Execute the Network Simplex algorithm
 10.1148        while (pivot.findEnteringArc()) {
 10.1149          findJoinNode();
 10.1150          bool change = findLeavingArc();
 10.1151 +        if (delta >= INF) return UNBOUNDED;
 10.1152          changeFlow(change);
 10.1153          if (change) {
 10.1154            updateTreeStructure();
 10.1155            updatePotential();
 10.1156          }
 10.1157        }
 10.1158 -
 10.1159 -      // Copy flow values to _flow_map
 10.1160 -      if (_plower) {
 10.1161 -        for (int i = 0; i != _arc_num; ++i) {
 10.1162 -          Arc e = _arc_ref[i];
 10.1163 -          _flow_map->set(e, (*_plower)[e] + _flow[i]);
 10.1164 -        }
 10.1165 -      } else {
 10.1166 -        for (int i = 0; i != _arc_num; ++i) {
 10.1167 -          _flow_map->set(_arc_ref[i], _flow[i]);
 10.1168 +      
 10.1169 +      // Check feasibility
 10.1170 +      if (_sum_supply < 0) {
 10.1171 +        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 10.1172 +          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
 10.1173          }
 10.1174        }
 10.1175 -      // Copy potential values to _potential_map
 10.1176 -      for (NodeIt n(_graph); n != INVALID; ++n) {
 10.1177 -        _potential_map->set(n, _pi[_node_id[n]]);
 10.1178 +      else if (_sum_supply > 0) {
 10.1179 +        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 10.1180 +          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
 10.1181 +        }
 10.1182 +      }
 10.1183 +      else {
 10.1184 +        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 10.1185 +          if (_flow[e] != 0) return INFEASIBLE;
 10.1186 +        }
 10.1187        }
 10.1188  
 10.1189 -      return true;
 10.1190 +      // Transform the solution and the supply map to the original form
 10.1191 +      if (_have_lower) {
 10.1192 +        for (int i = 0; i != _arc_num; ++i) {
 10.1193 +          Value c = _lower[i];
 10.1194 +          if (c != 0) {
 10.1195 +            _flow[i] += c;
 10.1196 +            _supply[_source[i]] += c;
 10.1197 +            _supply[_target[i]] -= c;
 10.1198 +          }
 10.1199 +        }
 10.1200 +      }
 10.1201 +
 10.1202 +      return OPTIMAL;
 10.1203      }
 10.1204  
 10.1205    }; //class NetworkSimplex
    11.1 --- a/lemon/preflow.h	Wed Apr 29 16:15:29 2009 +0100
    11.2 +++ b/lemon/preflow.h	Wed Apr 29 19:22:14 2009 +0100
    11.3 @@ -46,13 +46,13 @@
    11.4      typedef CAP CapacityMap;
    11.5  
    11.6      /// \brief The type of the flow values.
    11.7 -    typedef typename CapacityMap::Value Flow;
    11.8 +    typedef typename CapacityMap::Value Value;
    11.9  
   11.10      /// \brief The type of the map that stores the flow values.
   11.11      ///
   11.12      /// The type of the map that stores the flow values.
   11.13      /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
   11.14 -    typedef typename Digraph::template ArcMap<Flow> FlowMap;
   11.15 +    typedef typename Digraph::template ArcMap<Value> FlowMap;
   11.16  
   11.17      /// \brief Instantiates a FlowMap.
   11.18      ///
   11.19 @@ -84,7 +84,7 @@
   11.20      /// \brief The tolerance used by the algorithm
   11.21      ///
   11.22      /// The tolerance used by the algorithm to handle inexact computation.
   11.23 -    typedef lemon::Tolerance<Flow> Tolerance;
   11.24 +    typedef lemon::Tolerance<Value> Tolerance;
   11.25  
   11.26    };
   11.27  
   11.28 @@ -125,7 +125,7 @@
   11.29      ///The type of the capacity map.
   11.30      typedef typename Traits::CapacityMap CapacityMap;
   11.31      ///The type of the flow values.
   11.32 -    typedef typename Traits::Flow Flow;
   11.33 +    typedef typename Traits::Value Value;
   11.34  
   11.35      ///The type of the flow map.
   11.36      typedef typename Traits::FlowMap FlowMap;
   11.37 @@ -151,7 +151,7 @@
   11.38      Elevator* _level;
   11.39      bool _local_level;
   11.40  
   11.41 -    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
   11.42 +    typedef typename Digraph::template NodeMap<Value> ExcessMap;
   11.43      ExcessMap* _excess;
   11.44  
   11.45      Tolerance _tolerance;
   11.46 @@ -470,7 +470,7 @@
   11.47        }
   11.48  
   11.49        for (NodeIt n(_graph); n != INVALID; ++n) {
   11.50 -        Flow excess = 0;
   11.51 +        Value excess = 0;
   11.52          for (InArcIt e(_graph, n); e != INVALID; ++e) {
   11.53            excess += (*_flow)[e];
   11.54          }
   11.55 @@ -519,7 +519,7 @@
   11.56        _level->initFinish();
   11.57  
   11.58        for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
   11.59 -        Flow rem = (*_capacity)[e] - (*_flow)[e];
   11.60 +        Value rem = (*_capacity)[e] - (*_flow)[e];
   11.61          if (_tolerance.positive(rem)) {
   11.62            Node u = _graph.target(e);
   11.63            if ((*_level)[u] == _level->maxLevel()) continue;
   11.64 @@ -531,7 +531,7 @@
   11.65          }
   11.66        }
   11.67        for (InArcIt e(_graph, _source); e != INVALID; ++e) {
   11.68 -        Flow rem = (*_flow)[e];
   11.69 +        Value rem = (*_flow)[e];
   11.70          if (_tolerance.positive(rem)) {
   11.71            Node v = _graph.source(e);
   11.72            if ((*_level)[v] == _level->maxLevel()) continue;
   11.73 @@ -564,11 +564,11 @@
   11.74          int num = _node_num;
   11.75  
   11.76          while (num > 0 && n != INVALID) {
   11.77 -          Flow excess = (*_excess)[n];
   11.78 +          Value excess = (*_excess)[n];
   11.79            int new_level = _level->maxLevel();
   11.80  
   11.81            for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   11.82 -            Flow rem = (*_capacity)[e] - (*_flow)[e];
   11.83 +            Value rem = (*_capacity)[e] - (*_flow)[e];
   11.84              if (!_tolerance.positive(rem)) continue;
   11.85              Node v = _graph.target(e);
   11.86              if ((*_level)[v] < level) {
   11.87 @@ -591,7 +591,7 @@
   11.88            }
   11.89  
   11.90            for (InArcIt e(_graph, n); e != INVALID; ++e) {
   11.91 -            Flow rem = (*_flow)[e];
   11.92 +            Value rem = (*_flow)[e];
   11.93              if (!_tolerance.positive(rem)) continue;
   11.94              Node v = _graph.source(e);
   11.95              if ((*_level)[v] < level) {
   11.96 @@ -637,11 +637,11 @@
   11.97  
   11.98          num = _node_num * 20;
   11.99          while (num > 0 && n != INVALID) {
  11.100 -          Flow excess = (*_excess)[n];
  11.101 +          Value excess = (*_excess)[n];
  11.102            int new_level = _level->maxLevel();
  11.103  
  11.104            for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  11.105 -            Flow rem = (*_capacity)[e] - (*_flow)[e];
  11.106 +            Value rem = (*_capacity)[e] - (*_flow)[e];
  11.107              if (!_tolerance.positive(rem)) continue;
  11.108              Node v = _graph.target(e);
  11.109              if ((*_level)[v] < level) {
  11.110 @@ -664,7 +664,7 @@
  11.111            }
  11.112  
  11.113            for (InArcIt e(_graph, n); e != INVALID; ++e) {
  11.114 -            Flow rem = (*_flow)[e];
  11.115 +            Value rem = (*_flow)[e];
  11.116              if (!_tolerance.positive(rem)) continue;
  11.117              Node v = _graph.source(e);
  11.118              if ((*_level)[v] < level) {
  11.119 @@ -778,12 +778,12 @@
  11.120  
  11.121        Node n;
  11.122        while ((n = _level->highestActive()) != INVALID) {
  11.123 -        Flow excess = (*_excess)[n];
  11.124 +        Value excess = (*_excess)[n];
  11.125          int level = _level->highestActiveLevel();
  11.126          int new_level = _level->maxLevel();
  11.127  
  11.128          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  11.129 -          Flow rem = (*_capacity)[e] - (*_flow)[e];
  11.130 +          Value rem = (*_capacity)[e] - (*_flow)[e];
  11.131            if (!_tolerance.positive(rem)) continue;
  11.132            Node v = _graph.target(e);
  11.133            if ((*_level)[v] < level) {
  11.134 @@ -806,7 +806,7 @@
  11.135          }
  11.136  
  11.137          for (InArcIt e(_graph, n); e != INVALID; ++e) {
  11.138 -          Flow rem = (*_flow)[e];
  11.139 +          Value rem = (*_flow)[e];
  11.140            if (!_tolerance.positive(rem)) continue;
  11.141            Node v = _graph.source(e);
  11.142            if ((*_level)[v] < level) {
  11.143 @@ -897,18 +897,18 @@
  11.144      ///
  11.145      /// \pre Either \ref run() or \ref init() must be called before
  11.146      /// using this function.
  11.147 -    Flow flowValue() const {
  11.148 +    Value flowValue() const {
  11.149        return (*_excess)[_target];
  11.150      }
  11.151  
  11.152 -    /// \brief Returns the flow on the given arc.
  11.153 +    /// \brief Returns the flow value on the given arc.
  11.154      ///
  11.155 -    /// Returns the flow on the given arc. This method can
  11.156 +    /// Returns the flow value on the given arc. This method can
  11.157      /// be called after the second phase of the algorithm.
  11.158      ///
  11.159      /// \pre Either \ref run() or \ref init() must be called before
  11.160      /// using this function.
  11.161 -    Flow flow(const Arc& arc) const {
  11.162 +    Value flow(const Arc& arc) const {
  11.163        return (*_flow)[arc];
  11.164      }
  11.165  
    12.1 --- a/test/lp_test.cc	Wed Apr 29 16:15:29 2009 +0100
    12.2 +++ b/test/lp_test.cc	Wed Apr 29 19:22:14 2009 +0100
    12.3 @@ -21,9 +21,7 @@
    12.4  #include "test_tools.h"
    12.5  #include <lemon/tolerance.h>
    12.6  
    12.7 -#ifdef HAVE_CONFIG_H
    12.8  #include <lemon/config.h>
    12.9 -#endif
   12.10  
   12.11  #ifdef LEMON_HAVE_GLPK
   12.12  #include <lemon/glpk.h>
    13.1 --- a/test/min_cost_flow_test.cc	Wed Apr 29 16:15:29 2009 +0100
    13.2 +++ b/test/min_cost_flow_test.cc	Wed Apr 29 19:22:14 2009 +0100
    13.3 @@ -18,6 +18,7 @@
    13.4  
    13.5  #include <iostream>
    13.6  #include <fstream>
    13.7 +#include <limits>
    13.8  
    13.9  #include <lemon/list_graph.h>
   13.10  #include <lemon/lgf_reader.h>
   13.11 @@ -33,57 +34,57 @@
   13.12  
   13.13  char test_lgf[] =
   13.14    "@nodes\n"
   13.15 -  "label  sup1 sup2 sup3 sup4 sup5\n"
   13.16 -  "    1    20   27    0   20   30\n"
   13.17 -  "    2    -4    0    0   -8   -3\n"
   13.18 -  "    3     0    0    0    0    0\n"
   13.19 -  "    4     0    0    0    0    0\n"
   13.20 -  "    5     9    0    0    6   11\n"
   13.21 -  "    6    -6    0    0   -5   -6\n"
   13.22 -  "    7     0    0    0    0    0\n"
   13.23 -  "    8     0    0    0    0    3\n"
   13.24 -  "    9     3    0    0    0    0\n"
   13.25 -  "   10    -2    0    0   -7   -2\n"
   13.26 -  "   11     0    0    0  -10    0\n"
   13.27 -  "   12   -20  -27    0  -30  -20\n"
   13.28 -  "\n"
   13.29 +  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
   13.30 +  "    1    20   27    0   30   20   30\n"
   13.31 +  "    2    -4    0    0    0   -8   -3\n"
   13.32 +  "    3     0    0    0    0    0    0\n"
   13.33 +  "    4     0    0    0    0    0    0\n"
   13.34 +  "    5     9    0    0    0    6   11\n"
   13.35 +  "    6    -6    0    0    0   -5   -6\n"
   13.36 +  "    7     0    0    0    0    0    0\n"
   13.37 +  "    8     0    0    0    0    0    3\n"
   13.38 +  "    9     3    0    0    0    0    0\n"
   13.39 +  "   10    -2    0    0    0   -7   -2\n"
   13.40 +  "   11     0    0    0    0  -10    0\n"
   13.41 +  "   12   -20  -27    0  -30  -30  -20\n"
   13.42 +  "\n"                
   13.43    "@arcs\n"
   13.44 -  "       cost  cap low1 low2\n"
   13.45 -  " 1  2    70   11    0    8\n"
   13.46 -  " 1  3   150    3    0    1\n"
   13.47 -  " 1  4    80   15    0    2\n"
   13.48 -  " 2  8    80   12    0    0\n"
   13.49 -  " 3  5   140    5    0    3\n"
   13.50 -  " 4  6    60   10    0    1\n"
   13.51 -  " 4  7    80    2    0    0\n"
   13.52 -  " 4  8   110    3    0    0\n"
   13.53 -  " 5  7    60   14    0    0\n"
   13.54 -  " 5 11   120   12    0    0\n"
   13.55 -  " 6  3     0    3    0    0\n"
   13.56 -  " 6  9   140    4    0    0\n"
   13.57 -  " 6 10    90    8    0    0\n"
   13.58 -  " 7  1    30    5    0    0\n"
   13.59 -  " 8 12    60   16    0    4\n"
   13.60 -  " 9 12    50    6    0    0\n"
   13.61 -  "10 12    70   13    0    5\n"
   13.62 -  "10  2   100    7    0    0\n"
   13.63 -  "10  7    60   10    0    0\n"
   13.64 -  "11 10    20   14    0    6\n"
   13.65 -  "12 11    30   10    0    0\n"
   13.66 +  "       cost  cap low1 low2 low3\n"
   13.67 +  " 1  2    70   11    0    8    8\n"
   13.68 +  " 1  3   150    3    0    1    0\n"
   13.69 +  " 1  4    80   15    0    2    2\n"
   13.70 +  " 2  8    80   12    0    0    0\n"
   13.71 +  " 3  5   140    5    0    3    1\n"
   13.72 +  " 4  6    60   10    0    1    0\n"
   13.73 +  " 4  7    80    2    0    0    0\n"
   13.74 +  " 4  8   110    3    0    0    0\n"
   13.75 +  " 5  7    60   14    0    0    0\n"
   13.76 +  " 5 11   120   12    0    0    0\n"
   13.77 +  " 6  3     0    3    0    0    0\n"
   13.78 +  " 6  9   140    4    0    0    0\n"
   13.79 +  " 6 10    90    8    0    0    0\n"
   13.80 +  " 7  1    30    5    0    0   -5\n"
   13.81 +  " 8 12    60   16    0    4    3\n"
   13.82 +  " 9 12    50    6    0    0    0\n"
   13.83 +  "10 12    70   13    0    5    2\n"
   13.84 +  "10  2   100    7    0    0    0\n"
   13.85 +  "10  7    60   10    0    0   -3\n"
   13.86 +  "11 10    20   14    0    6  -20\n"
   13.87 +  "12 11    30   10    0    0  -10\n"
   13.88    "\n"
   13.89    "@attributes\n"
   13.90    "source 1\n"
   13.91    "target 12\n";
   13.92  
   13.93  
   13.94 -enum ProblemType {
   13.95 +enum SupplyType {
   13.96    EQ,
   13.97    GEQ,
   13.98    LEQ
   13.99  };
  13.100  
  13.101  // Check the interface of an MCF algorithm
  13.102 -template <typename GR, typename Flow, typename Cost>
  13.103 +template <typename GR, typename Value, typename Cost>
  13.104  class McfClassConcept
  13.105  {
  13.106  public:
  13.107 @@ -94,53 +95,46 @@
  13.108        checkConcept<concepts::Digraph, GR>();
  13.109  
  13.110        MCF mcf(g);
  13.111 +      const MCF& const_mcf = mcf;
  13.112  
  13.113        b = mcf.reset()
  13.114               .lowerMap(lower)
  13.115               .upperMap(upper)
  13.116 -             .capacityMap(upper)
  13.117 -             .boundMaps(lower, upper)
  13.118               .costMap(cost)
  13.119               .supplyMap(sup)
  13.120               .stSupply(n, n, k)
  13.121 -             .flowMap(flow)
  13.122 -             .potentialMap(pot)
  13.123               .run();
  13.124 -      
  13.125 -      const MCF& const_mcf = mcf;
  13.126  
  13.127 -      const typename MCF::FlowMap &fm = const_mcf.flowMap();
  13.128 -      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
  13.129 -
  13.130 -      v = const_mcf.totalCost();
  13.131 -      double x = const_mcf.template totalCost<double>();
  13.132 +      c = const_mcf.totalCost();
  13.133 +      x = const_mcf.template totalCost<double>();
  13.134        v = const_mcf.flow(a);
  13.135 -      v = const_mcf.potential(n);
  13.136 -
  13.137 -      ignore_unused_variable_warning(fm);
  13.138 -      ignore_unused_variable_warning(pm);
  13.139 -      ignore_unused_variable_warning(x);
  13.140 +      c = const_mcf.potential(n);
  13.141 +      const_mcf.flowMap(fm);
  13.142 +      const_mcf.potentialMap(pm);
  13.143      }
  13.144  
  13.145      typedef typename GR::Node Node;
  13.146      typedef typename GR::Arc Arc;
  13.147 -    typedef concepts::ReadMap<Node, Flow> NM;
  13.148 -    typedef concepts::ReadMap<Arc, Flow> FAM;
  13.149 +    typedef concepts::ReadMap<Node, Value> NM;
  13.150 +    typedef concepts::ReadMap<Arc, Value> VAM;
  13.151      typedef concepts::ReadMap<Arc, Cost> CAM;
  13.152 +    typedef concepts::WriteMap<Arc, Value> FlowMap;
  13.153 +    typedef concepts::WriteMap<Node, Cost> PotMap;
  13.154  
  13.155      const GR &g;
  13.156 -    const FAM &lower;
  13.157 -    const FAM &upper;
  13.158 +    const VAM &lower;
  13.159 +    const VAM &upper;
  13.160      const CAM &cost;
  13.161      const NM &sup;
  13.162      const Node &n;
  13.163      const Arc &a;
  13.164 -    const Flow &k;
  13.165 -    Flow v;
  13.166 +    const Value &k;
  13.167 +    FlowMap fm;
  13.168 +    PotMap pm;
  13.169      bool b;
  13.170 -
  13.171 -    typename MCF::FlowMap &flow;
  13.172 -    typename MCF::PotentialMap &pot;
  13.173 +    double x;
  13.174 +    typename MCF::Value v;
  13.175 +    typename MCF::Cost c;
  13.176    };
  13.177  
  13.178  };
  13.179 @@ -151,7 +145,7 @@
  13.180             typename SM, typename FM >
  13.181  bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
  13.182                  const SM& supply, const FM& flow,
  13.183 -                ProblemType type = EQ )
  13.184 +                SupplyType type = EQ )
  13.185  {
  13.186    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
  13.187  
  13.188 @@ -208,21 +202,25 @@
  13.189  // Run a minimum cost flow algorithm and check the results
  13.190  template < typename MCF, typename GR,
  13.191             typename LM, typename UM,
  13.192 -           typename CM, typename SM >
  13.193 -void checkMcf( const MCF& mcf, bool mcf_result,
  13.194 +           typename CM, typename SM,
  13.195 +           typename PT >
  13.196 +void checkMcf( const MCF& mcf, PT mcf_result,
  13.197                 const GR& gr, const LM& lower, const UM& upper,
  13.198                 const CM& cost, const SM& supply,
  13.199 -               bool result, typename CM::Value total,
  13.200 +               PT result, bool optimal, typename CM::Value total,
  13.201                 const std::string &test_id = "",
  13.202 -               ProblemType type = EQ )
  13.203 +               SupplyType type = EQ )
  13.204  {
  13.205    check(mcf_result == result, "Wrong result " + test_id);
  13.206 -  if (result) {
  13.207 -    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
  13.208 +  if (optimal) {
  13.209 +    typename GR::template ArcMap<typename SM::Value> flow(gr);
  13.210 +    typename GR::template NodeMap<typename CM::Value> pi(gr);
  13.211 +    mcf.flowMap(flow);
  13.212 +    mcf.potentialMap(pi);
  13.213 +    check(checkFlow(gr, lower, upper, supply, flow, type),
  13.214            "The flow is not feasible " + test_id);
  13.215      check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
  13.216 -    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
  13.217 -                         mcf.potentialMap()),
  13.218 +    check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
  13.219            "Wrong potentials " + test_id);
  13.220    }
  13.221  }
  13.222 @@ -231,11 +229,13 @@
  13.223  {
  13.224    // Check the interfaces
  13.225    {
  13.226 -    typedef int Flow;
  13.227 -    typedef int Cost;
  13.228      typedef concepts::Digraph GR;
  13.229 -    checkConcept< McfClassConcept<GR, Flow, Cost>,
  13.230 -                  NetworkSimplex<GR, Flow, Cost> >();
  13.231 +    checkConcept< McfClassConcept<GR, int, int>,
  13.232 +                  NetworkSimplex<GR> >();
  13.233 +    checkConcept< McfClassConcept<GR, double, double>,
  13.234 +                  NetworkSimplex<GR, double> >();
  13.235 +    checkConcept< McfClassConcept<GR, int, double>,
  13.236 +                  NetworkSimplex<GR, int, double> >();
  13.237    }
  13.238  
  13.239    // Run various MCF tests
  13.240 @@ -244,8 +244,8 @@
  13.241  
  13.242    // Read the test digraph
  13.243    Digraph gr;
  13.244 -  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
  13.245 -  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
  13.246 +  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
  13.247 +  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
  13.248    ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
  13.249    Node v, w;
  13.250  
  13.251 @@ -255,14 +255,56 @@
  13.252      .arcMap("cap", u)
  13.253      .arcMap("low1", l1)
  13.254      .arcMap("low2", l2)
  13.255 +    .arcMap("low3", l3)
  13.256      .nodeMap("sup1", s1)
  13.257      .nodeMap("sup2", s2)
  13.258      .nodeMap("sup3", s3)
  13.259      .nodeMap("sup4", s4)
  13.260      .nodeMap("sup5", s5)
  13.261 +    .nodeMap("sup6", s6)
  13.262      .node("source", v)
  13.263      .node("target", w)
  13.264      .run();
  13.265 +  
  13.266 +  // Build a test digraph for testing negative costs
  13.267 +  Digraph ngr;
  13.268 +  Node n1 = ngr.addNode();
  13.269 +  Node n2 = ngr.addNode();
  13.270 +  Node n3 = ngr.addNode();
  13.271 +  Node n4 = ngr.addNode();
  13.272 +  Node n5 = ngr.addNode();
  13.273 +  Node n6 = ngr.addNode();
  13.274 +  Node n7 = ngr.addNode();
  13.275 +  
  13.276 +  Arc a1 = ngr.addArc(n1, n2);
  13.277 +  Arc a2 = ngr.addArc(n1, n3);
  13.278 +  Arc a3 = ngr.addArc(n2, n4);
  13.279 +  Arc a4 = ngr.addArc(n3, n4);
  13.280 +  Arc a5 = ngr.addArc(n3, n2);
  13.281 +  Arc a6 = ngr.addArc(n5, n3);
  13.282 +  Arc a7 = ngr.addArc(n5, n6);
  13.283 +  Arc a8 = ngr.addArc(n6, n7);
  13.284 +  Arc a9 = ngr.addArc(n7, n5);
  13.285 +  
  13.286 +  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
  13.287 +  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
  13.288 +  Digraph::NodeMap<int> ns(ngr, 0);
  13.289 +  
  13.290 +  nl2[a7] =  1000;
  13.291 +  nl2[a8] = -1000;
  13.292 +  
  13.293 +  ns[n1] =  100;
  13.294 +  ns[n4] = -100;
  13.295 +  
  13.296 +  nc[a1] =  100;
  13.297 +  nc[a2] =   30;
  13.298 +  nc[a3] =   20;
  13.299 +  nc[a4] =   80;
  13.300 +  nc[a5] =   50;
  13.301 +  nc[a6] =   10;
  13.302 +  nc[a7] =   80;
  13.303 +  nc[a8] =   30;
  13.304 +  nc[a9] = -120;
  13.305  
  13.306    // A. Test NetworkSimplex with the default pivot rule
  13.307    {
  13.308 @@ -271,63 +313,77 @@
  13.309      // Check the equality form
  13.310      mcf.upperMap(u).costMap(c);
  13.311      checkMcf(mcf, mcf.supplyMap(s1).run(),
  13.312 -             gr, l1, u, c, s1, true,  5240, "#A1");
  13.313 +             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
  13.314      checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
  13.315 -             gr, l1, u, c, s2, true,  7620, "#A2");
  13.316 +             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
  13.317      mcf.lowerMap(l2);
  13.318      checkMcf(mcf, mcf.supplyMap(s1).run(),
  13.319 -             gr, l2, u, c, s1, true,  5970, "#A3");
  13.320 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
  13.321      checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
  13.322 -             gr, l2, u, c, s2, true,  8010, "#A4");
  13.323 +             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
  13.324      mcf.reset();
  13.325      checkMcf(mcf, mcf.supplyMap(s1).run(),
  13.326 -             gr, l1, cu, cc, s1, true,  74, "#A5");
  13.327 +             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
  13.328      checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
  13.329 -             gr, l2, cu, cc, s2, true,  94, "#A6");
  13.330 +             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
  13.331      mcf.reset();
  13.332      checkMcf(mcf, mcf.run(),
  13.333 -             gr, l1, cu, cc, s3, true,   0, "#A7");
  13.334 -    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
  13.335 -             gr, l2, u, cc, s3, false,   0, "#A8");
  13.336 +             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
  13.337 +    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
  13.338 +             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
  13.339 +    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
  13.340 +    checkMcf(mcf, mcf.run(),
  13.341 +             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
  13.342  
  13.343      // Check the GEQ form
  13.344 -    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
  13.345 +    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
  13.346      checkMcf(mcf, mcf.run(),
  13.347 -             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
  13.348 -    mcf.problemType(mcf.GEQ);
  13.349 +             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
  13.350 +    mcf.supplyType(mcf.GEQ);
  13.351      checkMcf(mcf, mcf.lowerMap(l2).run(),
  13.352 -             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
  13.353 -    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
  13.354 +             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
  13.355 +    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
  13.356      checkMcf(mcf, mcf.run(),
  13.357 -             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
  13.358 +             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
  13.359  
  13.360      // Check the LEQ form
  13.361 -    mcf.reset().problemType(mcf.LEQ);
  13.362 -    mcf.upperMap(u).costMap(c).supplyMap(s5);
  13.363 +    mcf.reset().supplyType(mcf.LEQ);
  13.364 +    mcf.upperMap(u).costMap(c).supplyMap(s6);
  13.365      checkMcf(mcf, mcf.run(),
  13.366 -             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
  13.367 +             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
  13.368      checkMcf(mcf, mcf.lowerMap(l2).run(),
  13.369 -             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
  13.370 -    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
  13.371 +             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
  13.372 +    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
  13.373      checkMcf(mcf, mcf.run(),
  13.374 -             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
  13.375 +             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
  13.376 +
  13.377 +    // Check negative costs
  13.378 +    NetworkSimplex<Digraph> nmcf(ngr);
  13.379 +    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
  13.380 +    checkMcf(nmcf, nmcf.run(),
  13.381 +      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
  13.382 +    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
  13.383 +      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
  13.384 +    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
  13.385 +    checkMcf(nmcf, nmcf.run(),
  13.386 +      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
  13.387    }
  13.388  
  13.389    // B. Test NetworkSimplex with each pivot rule
  13.390    {
  13.391      NetworkSimplex<Digraph> mcf(gr);
  13.392 -    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
  13.393 +    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
  13.394  
  13.395      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
  13.396 -             gr, l2, u, c, s1, true,  5970, "#B1");
  13.397 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
  13.398      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
  13.399 -             gr, l2, u, c, s1, true,  5970, "#B2");
  13.400 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
  13.401      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
  13.402 -             gr, l2, u, c, s1, true,  5970, "#B3");
  13.403 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
  13.404      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
  13.405 -             gr, l2, u, c, s1, true,  5970, "#B4");
  13.406 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
  13.407      checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
  13.408 -             gr, l2, u, c, s1, true,  5970, "#B5");
  13.409 +             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
  13.410    }
  13.411  
  13.412    return 0;
    14.1 --- a/test/mip_test.cc	Wed Apr 29 16:15:29 2009 +0100
    14.2 +++ b/test/mip_test.cc	Wed Apr 29 19:22:14 2009 +0100
    14.3 @@ -18,9 +18,7 @@
    14.4  
    14.5  #include "test_tools.h"
    14.6  
    14.7 -#ifdef HAVE_CONFIG_H
    14.8  #include <lemon/config.h>
    14.9 -#endif
   14.10  
   14.11  #ifdef LEMON_HAVE_CPLEX
   14.12  #include <lemon/cplex.h>
    15.1 --- a/tools/dimacs-solver.cc	Wed Apr 29 16:15:29 2009 +0100
    15.2 +++ b/tools/dimacs-solver.cc	Wed Apr 29 19:22:14 2009 +0100
    15.3 @@ -119,8 +119,8 @@
    15.4  
    15.5    ti.restart();
    15.6    NetworkSimplex<Digraph, Value> ns(g);
    15.7 -  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
    15.8 -  if (sum_sup > 0) ns.problemType(ns.LEQ);
    15.9 +  ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
   15.10 +  if (sum_sup > 0) ns.supplyType(ns.LEQ);
   15.11    if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
   15.12    ti.restart();
   15.13    bool res = ns.run();