Add an MCF operator to the separate solver

For the solver utility it'd be good to have MCF problems, so here they
are!

Format is:
  MCF<supplies, arcs>(cost*)

Supplies is a [int,int,int,...], where each int represents a new node
Arcs is [int:int, int:int, int:int, ...] where each int:int pair
represents an edge from the first to the second (1 indexed from the
"supplies" array).

Costs is the argument to the function. There must be as many costs as
arcs, and they are set from left to right, in order.

1 files changed, 4198 insertions, 0 deletions

diff --git a/impl/MCFSimplex.cpp b/impl/MCFSimplex.cpp
new file mode 100644
index 0000000..60c70d1
--- /dev/null
+++ b/impl/MCFSimplex.cpp
@@ -0,0 +1,4198 @@
+/*--------------------------------------------------------------------------*/
+/*---------------------------- File MCFSimplex.C ---------------------------*/
+/*--------------------------------------------------------------------------*/
+/*--                                                                      --*/
+/*-- Linear and Quadratic Min Cost Flow problems solver based on the      --*/
+/*-- (primal and dual) simplex algorithm. Conforms to the standard MCF    --*/
+/*-- interface defined in MCFClass.h.                                     --*/
+/*--                                                                      --*/
+/*--                            VERSION 1.00                              --*/
+/*--                           29 - 08 - 2011                             --*/
+/*--                                                                      --*/
+/*--                           Implementation:                            --*/
+/*--                                                                      --*/
+/*--                         Alessandro Bertolini                         --*/
+/*--                          Antonio Frangioni                           --*/
+/*--                                                                      --*/
+/*--                       Operations Research Group                      --*/
+/*--                      Dipartimento di Informatica                     --*/
+/*--                         Universita' di Pisa                          --*/
+/*--                                                                      --*/
+/*-- Copyright (C) 2008 - 2011 by Alessandro Bertolini, Antonio Frangioni --*/
+/*--                                                                      --*/
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*--------------------------- IMPLEMENTATION -------------------------------*/
+/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+/*--------------------------------------------------------------------------*/
+/*------------------------------ INCLUDES ----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#include "MCFSimplex.h"
+#include <algorithm>
+#include <iostream>
+
+#include <cstdlib>
+#include <ctime>
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------------- USING -----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#if( OPT_USE_NAMESPACES )
+using namespace MCFClass_di_unipi_it;
+#endif
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------------- MACROS ----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#define LIMITATEPRECISION 1
+
+/* If LIMITATEPRECISION is 1, in the quadratic case the Primal Simplex accepts
+   entering arc in base only if the decrease of the o.f. value is bigger than 
+   a fixed thresold (EpsOpt * oldFOValue / n). Otherwise, any strict decrease
+   in the o.f. value is accepted. */
+   
+#define UNIPI_PRIMAL_INITIAL_SHOW 0
+
+/* If UNIPI_PRIMAL_INITIAL_SHOW = 1, Primal Simplex shows the initial condition
+   (arcs and nodes) of the network. */
+
+#define UNIPI_PRIMAL_ITER_SHOW 0
+
+/* If UNIPI_PRIMAL_FINAL_SHOW = x with x > 0, Primal Simplex shows the condition
+   (arcs and nodes) of the network every x iterations. */
+
+#define UNIPI_PRIMAL_FINAL_SHOW 0
+
+/* If UNIPI_PRIMAL_FINAL_SHOW = 1, Primal Simplex shows the final condition
+   (arcs and nodes) of the network. */
+
+#define UNIPI_DUAL_INITIAL_SHOW 0
+
+/* If UNIPI_DUAL_INITIAL_SHOW = 1, Dual Simplex shows the initial condition
+   (arcs and nodes) of the network. */
+
+#define UNIPI_DUAL_ITER_SHOW 0
+
+/* If UNIPI_DUAL_FINAL_SHOW = x with x > 0, Dual Simplex shows the condition
+   (arcs and nodes) of the network every x iterations. */
+
+#define UNIPI_DUAL_FINAL_SHOW 0
+
+/* If UNIPI_DUAL_FINAL_SHOW = 1, Dual Simplex shows the final condition
+   (arcs and nodes) of the network. */
+
+#define UNIPI_VIS_DUMMY_ARCS 1
+
+/* If UNIPI_VIS_DUMMY_ARCS = 1, Primal Simplex or Dual Simplex shows the
+   conditions of the dummy arcs. */
+
+#define UNIPI_VIS_ARC_UPPER 1
+#define UNIPI_VIS_ARC_COST 1
+#define UNIPI_VIS_ARC_Q_COST 1
+#define UNIPI_VIS_ARC_REDUCT_COST 1
+#define UNIPI_VIS_ARC_STATE 1
+#define UNIPI_VIS_NODE_BASIC_ARC 1
+
+/* When Primal Simplex or Dual Simplex shows the conditions of the network, 
+   for every arcs the algorithm shows the flow, for every nodes it shows
+   balance and potential. These 6 flags decide if the algorithm shows a
+   particular value of the arcs/nodes;  for example if
+   UNIPI_VIS_ARC_UPPER == 1, the algorithm shows the upper bounds of all
+   arcs. */
+
+#define FOSHOW 0
+
+/* If FOSHOW is 1, the algorithm shows the f.o. value every x iterations 
+   (x = UNIPI_PRIMAL_ITER_SHOW or x = UNIPI_DUAL_ITER_SHOW). */
+
+#define OPTQUADRATIC 0
+
+/* If OPTQUADRATIC is 1 the Primal Simplex, in the quadratic case, tries to
+   optimize the update of the potential.
+   Unfortunately this doesn't work well: for this reason it is set to 0. */
+
+#ifndef throw
+// QUICK HACK
+#define throw(x) assert(false);
+#endif
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- FUNCTIONS ------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+template<class T>
+inline T ABS( const T x )
+{
+ return( x >= T( 0 ) ? x : - x );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class T>
+inline void Swap( T &v1 , T &v2 )
+{
+ T temp = v1;
+ v1 = v2;
+ v2 = temp;
+ }
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- CONSTANTS ------------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+#if( QUADRATICCOST == 0 )
+ static const char DELETED  = -2;  // ident for deleted arcs
+ static const char CLOSED   = -1;  // ident for closed arcs
+ static const char BASIC    =  0;  // ident for basis arcs
+ static const char AT_LOWER =  1;  // ident for arcs in L
+ static const char AT_UPPER =  2;  // ident for arcs in U
+#endif
+
+/* These macros will be used by method MemAllocCandidateList() to set the
+   values of numCandidateList and hotListSize. There are different macros,
+   according to:
+
+   - the used Simplex
+   - the size of the network 
+   - (obviously) the different variables
+
+   This set of values tries to improve the performance of the two algorithms
+   according to diversified situations. */
+
+static const int PRIMAL_LOW_NUM_CANDIDATE_LIST =  30;
+static const int PRIMAL_MEDIUM_NUM_CANDIDATE_LIST =  50;
+static const int PRIMAL_HIGH_NUM_CANDIDATE_LIST =  200;
+static const int PRIMAL_LOW_HOT_LIST_SIZE =  5;
+static const int PRIMAL_MEDIUM_HOT_LIST_SIZE =  10;
+static const int PRIMAL_HIGH_HOT_LIST_SIZE =  20;
+static const int DUAL_LOW_NUM_CANDIDATE_LIST =  6;
+static const int DUAL_HIGH_NUM_CANDIDATE_LIST =  10;
+static const int DUAL_LOW_HOT_LIST_SIZE =  1;
+static const int DUAL_HIGH_HOT_LIST_SIZE =  2;
+
+
+
+/*--------------------------------------------------------------------------*/
+/*--------------------------- COSTRUCTOR -----------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::MCFSimplex( cIndex nmx , cIndex mmx )
+            :
+            MCFClass( nmx , mmx )
+{
+ #if( QUADRATICCOST )
+  if( numeric_limits<FNumber>::is_integer )
+   throw( MCFException( "FNumber must be float if QUADRATICCOST == 1" ) );
+
+  if( numeric_limits<CNumber>::is_integer )
+   throw( MCFException( "CNumber must be float if QUADRATICCOST == 1" ) );
+ #endif
+
+ usePrimalSimplex = true;
+
+ newSession = true;
+ if( nmax && mmax )
+  MemAlloc();
+ else
+  nmax = mmax = 0;
+
+ #if( QUADRATICCOST )
+  recomputeFOLimits = 100;
+  // recomputeFOLimits represents the limit of the iteration in which 
+  // quadratic Primal Simplex computes "manually" the f.o. value
+  EpsOpt = 1e-13;
+  // EpsOpt is the fixed precision of the quadratic Primal Simplex
+ #endif
+
+ pricingRule = kCandidateListPivot;
+ forcedNumCandidateList = 0;
+ forcedHotListSize = 0;
+ usePrimalSimplex = true;
+ nodesP = NULL;
+ nodesD = NULL;
+ arcsP = NULL;
+ arcsD = NULL;
+ candP = NULL;
+ candD = NULL;
+
+ if( numeric_limits<CNumber>::is_integer )
+  MAX_ART_COST = CNumber( 1e7 );
+ else
+  MAX_ART_COST = CNumber( 1e10 );
+
+ }  // end( MCFSimplex )
+
+/*--------------------------------------------------------------------------*/
+/*-------------------------- OTHER INITIALIZATIONS -------------------------*/
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::LoadNet( cIndex nmx , cIndex mmx , cIndex pn , cIndex pm ,
+                          cFRow pU , cCRow pC , cFRow pDfct ,
+                          cIndex_Set pSn , cIndex_Set pEn )
+{
+ MemDeAllocCandidateList();
+
+ if( ( nmx != nmax ) || ( mmx != mmax ) ) {
+  // if the size of the allocated memory changes
+
+  if( nmax && mmax )  {  // if the memory was already allocated
+   MemDeAlloc(true);         // deallocate the Primal 
+   MemDeAlloc(false);        // and the Dual data structures
+   nmax = mmax = 0;
+   }
+
+  if( nmx && mmx ) {   // if the new dimensions of the memory are correct
+   nmax = nmx;
+   mmax = mmx;
+   MemAlloc();
+   }
+  }
+
+ if( ( ! nmax ) || ( ! mmax ) )
+  // if one of the two dimension of the memory isn't correct
+  nmax = mmax = 0;
+
+ if( nmax ) {  // if the new dimensions of the memory are correct
+  n = pn;
+  m = pm;
+
+  if( usePrimalSimplex ) {
+   stopNodesP = nodesP + n;
+   dummyRootP = nodesP + nmax;
+   for( nodePType *node = nodesP ; node != stopNodesP ; node++ )
+    node->balance = pDfct[ node - nodesP ];  // initialize nodes
+
+   stopArcsP = arcsP + m;
+   dummyArcsP = arcsP + mmax;
+   stopDummyP = dummyArcsP + n;
+   for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++ ) {
+    // initialize real arcs
+    if (pC)
+      arc->cost = pC[ arc - arcsP ];
+    else
+      arc->cost = 0;
+    #if( QUADRATICCOST )
+     arc->quadraticCost = 0; 
+    #endif
+    if (pU)
+      arc->upper = pU[ arc - arcsP ];
+    else
+      arc->upper = Inf<FNumber>();
+    arc->tail = nodesP + pSn[ arc - arcsP ] - 1;
+    arc->head = nodesP + pEn[ arc - arcsP ] - 1;
+    }
+   }
+  else {
+   stopNodesD = nodesD + n;
+   dummyRootD = nodesD + nmax;
+   for( nodeDType *node = nodesD ; node != stopNodesD ; node++ )
+    node->balance = pDfct[ node - nodesD ];  // initialize nodes
+
+   stopArcsD = arcsD + m;
+   dummyArcsD = arcsD + mmax;
+   stopDummyD = dummyArcsD + n;
+   for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) {
+    // initialize real arcs
+    arc->cost = pC[ arc - arcsD ];
+    #if( QUADRATICCOST )
+     arc->quadraticCost = 0; 
+    #endif
+    arc->upper = pU[ arc - arcsD ];
+    arc->tail = nodesD + pSn[ arc - arcsD ] - 1;
+    arc->head = nodesD + pEn[ arc - arcsD ] - 1;
+    }
+
+   CreateAdditionalDualStructures();
+   }
+
+  if( pricingRule == kCandidateListPivot )
+   MemAllocCandidateList();
+
+  status = kUnSolved;
+  }
+ }  // end( MCFSimplex::LoadNet )
+                
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::SetAlg( bool UsPrml , char WhchPrc )
+{
+ bool newUsePrimalSimplex = UsPrml;
+ bool oldUsePrimalSimplex = usePrimalSimplex;
+ char newPricingRule = WhchPrc;
+ char oldPricingRule = pricingRule;
+
+ usePrimalSimplex = newUsePrimalSimplex;
+ pricingRule = newPricingRule;
+
+ if( ( ! usePrimalSimplex ) && ( pricingRule == kDantzig) ) {
+  pricingRule = kFirstEligibleArc;
+  newPricingRule = pricingRule;
+  }
+
+ if( ( newUsePrimalSimplex != oldUsePrimalSimplex ) ||
+     ( newPricingRule != oldPricingRule ) ) {
+  MemDeAllocCandidateList();
+
+  if( newUsePrimalSimplex != oldUsePrimalSimplex ) {
+   #if( QUADRATICCOST )
+     throw(
+      MCFException( "Primal Simplex is the only option if QUADRATICCOST == 1"
+                    ) );
+   }
+   #else
+    MemAlloc();
+    nodePType *nP = nodesP;
+    nodeDType *nD = nodesD;
+    arcPType *aP = arcsP;
+    arcDType *aD = arcsD;
+
+    if( newUsePrimalSimplex ) { // from Dual to Primal
+     if( nodesD == NULL )
+      return;
+
+     if( newSession )
+      CreateInitialDualBase();
+
+     dummyRootP = nodesP + nmax;
+     stopNodesP = nodesP + n;
+     dummyArcsP = arcsP + mmax;
+     stopArcsP = arcsP + m;
+     stopDummyP = dummyArcsP + n;
+     // Copy the old Dual data structure in a new Primal data structure
+     while( nD != stopNodesD ) {
+      nP->prevInT = nodesP + ( nD->prevInT - nodesD );
+      nP->nextInT = nodesP + ( nD->nextInT - nodesD );
+      nP->enteringTArc = arcsP + ( nD->enteringTArc - arcsD );
+      nP->balance = nD->balance;
+      nP->potential = nD->potential;
+      nP->subTreeLevel = nD->subTreeLevel;
+      nP++;
+      nD++;
+      }
+
+     dummyRootP->prevInT = NULL;
+     dummyRootP->nextInT = nodesP + ( dummyRootD->nextInT - nodesD );
+     dummyRootP->enteringTArc = arcsP + ( dummyRootD->enteringTArc - arcsD );
+     dummyRootP->balance = dummyRootD->balance;
+     dummyRootP->potential = dummyRootD->potential;
+     dummyRootP->subTreeLevel = dummyRootD->subTreeLevel;
+     while( aD != stopArcsD ) {
+      aP->tail = nodesP + ( aD->tail - nodesD );
+      aP->head = nodesP + ( aD->head - nodesD );
+      aP->flow = aD->flow;
+      aP->cost = aD->cost;
+      aP->ident = aD->ident;
+      aP->upper = aD->upper;
+      aP++;
+      aD++;
+      }
+
+     aP = dummyArcsP;
+     aD = dummyArcsD;
+     while( aD != stopDummyD ) {
+      aP->tail = nodesP + ( aD->tail - nodesD );
+      aP->head = nodesP + ( aD->head - nodesD );
+      aP->flow = aD->flow;
+      aP->cost = aD->cost;
+      aP->ident = aD->ident;
+      if( ( ETZ(aP->flow, EpsFlw) ) && (aP->ident == AT_UPPER) )
+       aP->ident = AT_LOWER;
+
+      aP->upper = Inf<FNumber>();
+      aP++;
+      aD++;
+      }
+
+     MemDeAlloc(false);
+     if( Senstv && ( status != kUnSolved ) ) {
+      nodePType *node = dummyRootP;
+      for( int i = 0 ; i < n ; i++ )
+       node = node->nextInT;
+
+      node->nextInT = NULL;
+      dummyRootP->prevInT = NULL;
+      dummyRootP->enteringTArc = NULL;
+      // balance the flow
+      CreateInitialPModifiedBalanceVector();
+      PostPVisit( dummyRootP );
+      // restore the primal admissibility
+      BalanceFlow( dummyRootP );
+      ComputePotential( dummyRootP );
+      }
+     else
+      status = kUnSolved;
+     }
+   else {  // from Primal to Dual
+    if( nodesP == NULL )
+     return;
+
+    if( newSession )
+     CreateInitialPrimalBase();
+
+    dummyRootD = nodesD + nmax;
+    stopNodesD = nodesD + n;
+    dummyArcsD = arcsD + mmax;
+    stopArcsD = arcsD + m;
+    stopDummyD = dummyArcsD + n;
+    // Copy the old Primal data structure in a new Dual data structure
+    while( nP != stopNodesP ) {
+     nD->prevInT = nodesD + ( nP->prevInT - nodesP );
+     nD->nextInT = nodesD + ( nP->nextInT - nodesP );
+     nD->enteringTArc = arcsD + ( nP->enteringTArc - arcsP );
+     nD->balance = nP->balance;
+     nD->potential = nP->potential;
+     nD->subTreeLevel = nP->subTreeLevel;
+     nP++;
+     nD++;
+     }
+
+    dummyRootD->prevInT = NULL;
+    dummyRootD->nextInT = nodesD + ( dummyRootP->nextInT - nodesP );
+    dummyRootD->enteringTArc = NULL;
+    dummyRootD->balance = dummyRootP->balance;
+    dummyRootD->potential = dummyRootP->potential;
+    dummyRootD->subTreeLevel = dummyRootP->subTreeLevel;
+    while( aP != stopArcsP ) {
+     aD->tail = nodesD + ( aP->tail - nodesP );
+     aD->head = nodesD + ( aP->head - nodesP );
+     aD->flow = aP->flow;
+     aD->cost = aP->cost;
+     aD->ident = aP->ident;
+     aD->upper = aP->upper;
+     aP++;
+     aD++;
+     }
+
+    aP = dummyArcsP;
+    aD = dummyArcsD;
+    while( aP != stopDummyP ) {
+     aD->tail = nodesD + ( aP->tail - nodesP );
+     aD->head = nodesD + ( aP->head - nodesP );
+     aD->flow = aP->flow;
+     aD->cost = aP->cost;
+     aD->ident = aP->ident;
+     aD->upper = 0;
+     aP++;
+     aD++;
+     }
+
+    CreateAdditionalDualStructures();
+    MemDeAlloc(true);
+    nodeDType *node = dummyRootD;
+    for( int i = 0 ; i < n ; i++ )
+     node = node->nextInT;
+
+    node->nextInT = NULL;
+    dummyRootD->enteringTArc = NULL;
+    dummyRootD->prevInT = NULL;
+    if( Senstv && ( status != kUnSolved ) ) {
+     // fix every flow arc according to its reduct cost
+     for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) {
+      if( (arc->ident == AT_LOWER) && LTZ( ReductCost( arc ) , EpsCst ) ) {
+       arc->flow = arc->upper;
+       arc->ident = AT_UPPER;
+       }
+
+      if( ( arc->ident == AT_UPPER ) && GTZ( ReductCost( arc ) , EpsCst ) ) {
+       arc->flow = 0;
+       arc->ident = AT_LOWER;
+       }
+      }
+
+     // balance the flow
+     CreateInitialDModifiedBalanceVector();
+     PostDVisit( dummyRootD );
+     }
+    else
+     status = kUnSolved;
+    }
+   //#endif
+   }
+#endif
+ if( newPricingRule == kFirstEligibleArc )
+  if( newUsePrimalSimplex )
+   arcToStartP = arcsP;
+  else
+   arcToStartD = arcsD;
+
+ if( ( nmax && mmax ) && ( newPricingRule == kCandidateListPivot ) )
+  MemAllocCandidateList();
+  }
+ }  // end( SetAlg )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::SetPar( int par, int val )
+{
+ switch( par ) {
+ case kAlgPrimal:
+  if( val == kYes )
+   SetAlg( true , pricingRule);
+
+  if( val == kNo )
+   SetAlg( false, pricingRule );
+
+  break;
+
+ case kAlgPricing:
+ 
+  if( ( val == kDantzig ) || ( val == kFirstEligibleArc ) ||
+      ( val == kCandidateListPivot ) )
+   SetAlg( usePrimalSimplex , val );
+
+  break;
+
+ case kNumCandList:
+
+  MemDeAllocCandidateList();
+  forcedNumCandidateList = val;
+  MemAllocCandidateList();
+  forcedNumCandidateList = 0;
+  forcedHotListSize = 0;
+  break;
+
+ case kHotListSize:
+
+  MemDeAllocCandidateList();
+  forcedHotListSize = val;
+  MemAllocCandidateList();
+  forcedNumCandidateList = 0;
+  forcedHotListSize = 0;
+  break;
+
+ case kRecomputeFOLimits:
+
+  recomputeFOLimits = val;
+  break;
+
+ default:
+
+  MCFClass::SetPar(par, val);
+ }
+ }  // end( SetPar( int ) )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::SetPar( int par , double val )
+{
+ switch( par ) {
+ case kEpsOpt:
+
+  EpsOpt = val;
+  break;
+
+ default:
+  MCFClass::SetPar( par , val );
+ }
+ }  // end( SetPar( double )
+
+/*-------------------------------------------------------------------------*/
+/*--------------- METHODS FOR SOLVING THE PROBLEM -------------------------*/
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::SolveMCF( void )
+{
+ if( MCFt )
+  MCFt->Start();
+
+ //if( status == kUnSolved )
+  if(usePrimalSimplex )
+   CreateInitialPrimalBase();
+  else
+   CreateInitialDualBase();
+
+ newSession = false;
+ if( usePrimalSimplex )
+  PrimalSimplex();
+ else
+  DualSimplex();
+
+ if( MCFt )
+  MCFt->Stop();
+
+ }  // end( MCFSimplex::SolveMCF )
+
+/*--------------------------------------------------------------------------*/
+/*---------------------- METHODS FOR READING RESULTS -----------------------*/
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFGetX( FRow F , Index_Set nms , cIndex strt , Index stp )
+{
+ if( stp > m )
+  stp = m;
+
+ if( nms ) {
+  if( usePrimalSimplex )
+   for( Index i = strt ; i < stp ; i++ ) {
+    FNumber tXi = ( arcsP + i )->flow;
+    if( GTZ( tXi , EpsFlw ) ) {
+     *(F++) = tXi;
+     *(nms++) = i;
+     }
+    }
+  else
+   for( Index i = strt ; i < stp ; i++ ) {
+    FNumber tXi = ( arcsD + i )->flow;
+    if( GTZ( tXi , EpsFlw ) ) {
+     *(F++) = tXi;
+     *(nms++) = i;
+     }
+    }
+
+  *nms = Inf<Index>();
+  }        
+ else
+  if( usePrimalSimplex )
+   for( Index i = strt; i < stp; i++ )
+    *(F++) = ( arcsP + i )->flow;
+  else
+   for( Index i = strt; i < stp; i++ )
+    *(F++) = ( arcsD + i )->flow;
+
+ }  // end( MCFSimplex::MCFGetX( some ) )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFGetRC( CRow CR , cIndex_Set nms , cIndex strt , Index stp )
+{
+ if( nms ) {
+  while( *nms < strt )
+   nms++;
+
+  if( usePrimalSimplex )
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+   *(CR++) = CNumber( ReductCost( arcsP + h ) );
+  else
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+    *(CR++) = ReductCost( arcsD + h );
+  }
+ else {
+  if( stp > m )
+   stp = m;
+
+  if( usePrimalSimplex )
+   for( arcPType* arc = arcsP + strt ; arc < arcsP + stp ; arc++ )
+    *(CR++) = CNumber( ReductCost( arc ) );
+  else
+   for( arcDType* arc = arcsD + strt ; arc < arcsD + stp ; arc++ )
+    *(CR++) = ReductCost( arc );
+  }
+ }  // end( MCFSimplex::MCFGetRC( some ) )
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::CNumber MCFSimplex::MCFGetRC( cIndex i )
+{
+ if( usePrimalSimplex )
+  return CNumber( ReductCost( arcsP + i ) );
+ else
+  return( ReductCost( arcsD + i ) );
+
+ }  // end( MCFSimplex::MCFGetRC( i ) )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFGetPi( CRow P , cIndex_Set nms , cIndex strt , Index stp )
+{
+ if( stp > n )
+  stp = n;
+
+ if( nms ) {
+  while( *nms < strt )
+   nms++;
+
+  if( usePrimalSimplex )
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+    *(P++) = CNumber( (nodesP + h)->potential );
+  else
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+    *(P++) = (nodesD + h )->potential;
+  }
+ else
+  if( usePrimalSimplex )
+   for( nodePType *node = nodesP + strt ; node < ( nodesP + stp ) ; node++ )
+    *(P++) = CNumber( node->potential );
+  else
+   for( nodeDType *node = nodesD + strt ; node++ < ( nodesD + stp ) ; node++ )
+    *(P++) = node->potential;
+
+ }  // end(  MCFSimplex::MCFGetPi( some ) )
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::FONumber MCFSimplex::MCFGetFO( void )
+{
+ if( status == kOK )
+  return( (FONumber) GetFO() );
+ else
+  if( status == kUnbounded ) 
+   return( - Inf<FONumber>() );
+  else
+   return( Inf<FONumber>() );
+
+ }  // end( MCFSimplex::MCFGetFO )
+
+/*-------------------------------------------------------------------------*/
+/*----------METHODS FOR READING THE DATA OF THE PROBLEM--------------------*/
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFArcs( Index_Set Startv , Index_Set Endv ,
+			  cIndex_Set nms , cIndex strt , Index stp )
+{
+ if( stp > m )
+  stp = m;
+
+ if( nms ) {
+  while( *nms < strt )
+   nms++;
+
+  if( usePrimalSimplex )
+   for( Index h ; ( h = *(nms++) ) < stp ; ) {
+    if( Startv )
+     *(Startv++) = Index( (arcsP + h)->tail - nodesP) - USENAME0;
+    if( Endv )
+     *(Endv++) = Index( (arcsP + h)->head - nodesP ) - USENAME0;
+    }
+  else
+   for( Index h ; ( h = *(nms++) ) < stp ; ) {
+    if( Startv )
+     *(Startv++) = Index( (arcsD + h)->tail - nodesD) - USENAME0;
+    if( Endv )
+     *(Endv++) = Index( (arcsD + h)->head - nodesD ) - USENAME0;
+    }
+  }
+ else
+  if( usePrimalSimplex )
+   for( arcPType* arc = arcsP + strt ; arc < (arcsP + stp) ; arc++ ) {
+    if( Startv )
+     *(Startv++) = Index( arc->tail - nodesP ) - USENAME0;
+    if( Endv )
+     *(Endv++) = Index( arc->head - nodesP ) - USENAME0;
+    }
+  else
+   for( arcDType* arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ ) {
+    if( Startv )
+     *(Startv++) = Index( arc->tail - nodesD ) - USENAME0;
+    if( Endv )
+     *(Endv++) = Index( arc->head - nodesD ) - USENAME0;
+    }
+
+ }  // end( MCFSimplex::MCFArcs )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFCosts( CRow Costv , cIndex_Set nms ,
+			   cIndex strt , Index stp )
+{
+ if( stp > m )
+  stp = m;
+
+ if( nms ) {
+  while( *nms < strt )
+   nms++;
+
+  if( usePrimalSimplex )
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+    *(Costv++) = (arcsP + h)->cost;
+  else
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+    *(Costv++) = (arcsD + h)->cost;
+  }
+ else
+  if( usePrimalSimplex )
+   for( arcPType* arc = arcsP + strt ; arc < (arcsP + stp) ; arc++ )
+    *(Costv++) = arc->cost;           
+  else
+   for( arcDType* arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ )
+    *(Costv++) = arc->cost;           
+
+ }  // end( MCFSimplex::MCFCosts ( some ) )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFQCoef( CRow Qv , cIndex_Set nms  ,
+			   cIndex strt , Index stp )
+{
+ if( stp > m )
+  stp = m;
+
+ if( nms ) {
+  while( *nms < strt )
+   nms++;
+
+  #if( QUADRATICCOST )
+   if( usePrimalSimplex )
+    for( Index h ; ( h = *(nms++) ) < stp ; )
+     *(Qv++) = (arcsP + h)->quadraticCost;
+   else
+    for( Index h ; ( h = *(nms++) ) < stp ; )
+     *(Qv++) = (arcsD + h)->quadraticCost;
+  #else
+    for( Index h ; ( h = *(nms++) ) < stp ; )
+     *(Qv++) = 0;
+  #endif
+  }
+ else
+  #if( QUADRATICCOST )
+   if( usePrimalSimplex )
+    for( arcPType* arc = arcsP + strt ; arc < ( arcsP + stp ) ; arc++ )
+     *(Qv++) = arc->quadraticCost;
+   else
+    for( arcDType* arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ )
+     *(Qv++) = arc->quadraticCost;
+  #else
+   for( Index h = strt ; h++ < stp ; )
+    *(Qv++) = 0;           
+  #endif
+
+ }  // end( MCFSimplex::MCFQCoef ( some ) )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFUCaps( FRow UCapv , cIndex_Set nms ,
+			   cIndex strt , Index stp ) 
+{
+ if( stp > m )
+  stp = m;
+
+ if( nms ) {
+  while( *nms < strt )
+   nms++;
+
+  if( usePrimalSimplex )
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+    *(UCapv++) = (arcsP + h)->upper;
+  else
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+    *(UCapv++) = (arcsD + h)->upper;
+  }
+ else
+  if( usePrimalSimplex )
+   for( arcPType* arc = arcsP + strt ; arc <  (arcsP + stp ) ; arc++ )
+    *(UCapv++) = arc->upper;
+  else
+   for( arcDType* arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ )
+    *(UCapv++) = arc->upper;
+
+ }  // end( MCFSimplex::MCFUCaps ( some ) )
+ 
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::MCFDfcts( FRow Dfctv , cIndex_Set nms ,
+			   cIndex strt , Index stp )
+{
+ if( stp > n )
+  stp = n;
+
+ if( nms ) {
+  while( *nms < strt )
+   nms++;
+
+  if( usePrimalSimplex )
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+    *(Dfctv++) = ( nodesP + h )->balance;
+  else
+   for( Index h ; ( h = *(nms++) ) < stp ; )
+    *(Dfctv++) = (nodesD + h )->balance;
+  }
+ else
+  if( usePrimalSimplex )
+   for( nodePType* node = nodesP + strt ; node < ( nodesP + stp ) ; node++ )
+    *(Dfctv++) = node->balance;
+  else
+   for( nodeDType* node = nodesD + strt ; node < ( nodesD + stp ) ; node++ )
+    *(Dfctv++) = node->balance;
+
+ }  // end( MCFSimplex::MCFDfcts )
+
+/*-------------------------------------------------------------------------*/
+/*--------- METHODS FOR ADDING / REMOVING / CHANGING DATA -----------------*/
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgCosts( cCRow NCost , cIndex_Set nms ,
+			   cIndex strt , Index stp )
+{
+ if( stp > m )
+  stp = m;
+
+ if( nms ) {
+  while( *nms < strt ) {
+   nms++;
+   NCost++;
+   }
+
+  cIndex_Set tnms = nms;  // nms may be needed below
+  if( usePrimalSimplex )
+   for( Index h ; ( h = *(tnms++) ) < stp ; )
+    arcsP[ h ].cost = *(NCost++);
+  else
+   for( Index h ; ( h = *(tnms++) ) < stp ; )
+    arcsD[ h ].cost = *(NCost++);
+  }
+ else
+  if( usePrimalSimplex )
+   for( arcPType *arc = arcsP + strt ; arc < (arcsP + stp) ; arc++ )
+    arc->cost = *(NCost++); 
+  else
+   for( arcDType *arc = arcsD + strt ; arc < (arcsD + stp) ; arc++ )
+    arc->cost = *(NCost++); 
+
+ if( Senstv && ( status != kUnSolved ) )
+  if( usePrimalSimplex )
+   ComputePotential( dummyRootP );
+  else {
+   #if( QUADRATICCOST == 0 )
+    for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ )
+     if( arc->ident > BASIC )
+      if( GTZ( ReductCost( arc ) , EpsCst ) ) {
+       arc->flow = 0;
+       arc->ident = AT_LOWER;
+       }
+      else {
+       arc->flow = arc->upper; 
+       arc->ident = AT_UPPER;
+       }
+
+    CreateInitialDModifiedBalanceVector();
+    PostDVisit( dummyRootD );
+   #endif
+   }
+ else
+  status = kUnSolved;
+
+ }  // end( MCFSimplex::ChgCosts )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgCost( Index arc , cCNumber NCost )
+{
+ if( arc >= m )
+  return;
+
+ if( usePrimalSimplex )
+  ( arcsP + arc )->cost = NCost;
+ else
+  ( arcsD + arc )->cost = NCost;
+
+ if( Senstv && ( status != kUnSolved ) ) {
+  if( usePrimalSimplex ) {
+   nodePType *node = ( arcsP + arc )->tail;
+   if( ( ( arcsP + arc )->head)->subTreeLevel < node->subTreeLevel )
+    node = ( arcsP + arc )->head;
+
+   ComputePotential( dummyRootP );
+   }
+  else {
+   #if( QUADRATICCOST == 0 )
+    nodeDType *node = ( arcsD + arc )->tail;
+    if( ( ( arcsD + arc )->head )->subTreeLevel < node->subTreeLevel )
+     node = ( arcsD + arc )->head;
+
+    ComputePotential( dummyRootD );
+    for( arcDType *a = arcsD ; a != stopArcsD ; a++)
+     if( a->ident > BASIC )
+      if( GTZ( ReductCost( a ) , EpsCst ) ) {
+       a->flow = 0;
+       a->ident = AT_LOWER;
+       }
+      else {
+       a->flow = a->upper; 
+       a->ident = AT_UPPER;
+       }
+
+    CreateInitialDModifiedBalanceVector();
+    PostDVisit( dummyRootD );
+   #endif
+   }
+  }
+ else
+  status = kUnSolved;
+
+ }  // end( MCFSimplex::ChgCost )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgQCoef( cCRow NQCoef , cIndex_Set nms ,
+			   cIndex strt , Index stp )
+{
+ if( stp > m )
+  stp = m;
+
+ #if( QUADRATICCOST )
+  if( nms ) {
+   while( *nms < strt ) {
+    nms++;
+    NQCoef++;
+    }
+
+   cIndex_Set tnms = nms;  // nms may be needed below
+   if( usePrimalSimplex )
+    for( Index h ; ( h = *(tnms++) ) < stp ; )
+     arcsP[ h ].quadraticCost = *(NQCoef++);
+   else
+    for( Index h ; ( h = *(tnms++) ) < stp ; )
+     arcsD[ h ].quadraticCost = *(NQCoef++);
+   }
+  else
+   if( usePrimalSimplex )
+    for( arcPType *arc = arcsP + strt ; arc < ( arcsP + stp ) ; arc++ )
+     arc->quadraticCost = *(NQCoef++); 
+   else
+    for( arcDType *arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ )
+     arc->quadraticCost = *(NQCoef++); 
+
+  if( Senstv && (status != kUnSolved ) )
+   ComputePotential( dummyRootP );
+  else
+   status = kUnSolved;
+ #else
+  if( NQCoef )
+   throw( MCFException( "ChgQCoef: nonzero coefficients not allowed" ) );
+ #endif
+}  // end( MCFSimplex::ChgQCoef )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgQCoef( Index arc , cCNumber NQCoef )
+{
+ #if( QUADRATICCOST )
+  if( arc >= m )
+   return;
+
+  if( usePrimalSimplex )
+   ( arcsP + arc )->quadraticCost = NQCoef;
+  else
+   ( arcsD + arc )->quadraticCost = NQCoef;
+
+  if( Senstv && ( status != kUnSolved ) ) {
+   nodePType *node = ( arcsP + arc )->tail;
+   if( ( ( arcsP + arc )->head )->subTreeLevel < node->subTreeLevel )
+    node = ( arcsP + arc )->head;
+
+   ComputePotential( node );
+   }
+  else
+   status = kUnSolved;
+ #else
+  if( NQCoef )
+   throw( MCFException( "ChgQCoef: nonzero coefficients not allowed" ) );
+ #endif
+ }  // end( MCFSimplex::ChgQCoef )
+
+/*-------------------------------------------------------------------------*/
+    
+void MCFSimplex::ChgDfcts( cFRow NDfct , cIndex_Set nms ,
+			   cIndex strt , Index stp )
+{
+ if( stp > m )
+  stp = m;
+
+ if( nms ) {
+  while( *nms < strt ) {
+   nms++;
+   NDfct++;
+   }
+
+  cIndex_Set tnms = nms;  // nms may be needed below
+  if( usePrimalSimplex )
+   for( Index h ; ( h = *(tnms++) ) < stp ; )
+    nodesP[ h ].balance = *(NDfct++);
+  else
+   for( Index h ; ( h = *(tnms++) ) < stp ; )
+    nodesD[ h ].balance = *(NDfct++);
+  }
+ else
+  if( usePrimalSimplex )
+   for( nodePType *node = nodesP + strt ; node < ( nodesP + stp ) ; node++ )
+    node->balance = *(NDfct++);
+  else
+   for( nodeDType *node = nodesD + strt ; node < ( nodesD + stp ) ; node++ )
+    node->balance = *(NDfct++);
+
+ if( Senstv && (status != kUnSolved ) )
+  if( usePrimalSimplex ) {
+   CreateInitialPModifiedBalanceVector();
+   PostPVisit( dummyRootP );
+   BalanceFlow( dummyRootP );
+   ComputePotential( dummyRootP );
+   }
+  else {
+   CreateInitialDModifiedBalanceVector();
+   PostDVisit( dummyRootD );
+   }
+ else
+  status = kUnSolved;
+
+ }  // end( MCFSimplex::ChgDfcts )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgDfct( Index nod , cFNumber NDfct )
+{ 
+ if( nod > n )
+  return;
+
+ if( usePrimalSimplex )
+  ( nodesP + nod - 1 )->balance = NDfct;
+ else
+  ( nodesD + nod - 1 )->balance = NDfct;
+
+ if( Senstv && (status != kUnSolved ) )
+  if( usePrimalSimplex ) {
+   CreateInitialPModifiedBalanceVector();
+   PostPVisit( dummyRootP );
+   BalanceFlow( dummyRootP );
+   ComputePotential( dummyRootP );
+   }
+  else {
+   CreateInitialDModifiedBalanceVector();
+   PostDVisit( dummyRootD );
+   }
+ else
+  status = kUnSolved;
+
+ }  // end( MCFSimplex::ChgDfct )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgUCaps( cFRow NCap , cIndex_Set nms ,
+			   cIndex strt , Index stp )
+{
+ if( stp > m )
+  stp = m;
+
+ if( nms ) {
+  while( *nms < strt ) {
+   nms++;
+   NCap++;
+   }
+
+  cIndex_Set tnms = nms;  // nms may be needed below
+  if( usePrimalSimplex )
+   for( Index h ; ( h = *(tnms++) ) < stp ; )
+    arcsP[ h ].upper = *(NCap++);
+  else
+   for( Index h ; ( h = *(tnms++) ) < stp ; )
+    arcsD[ h ].upper = *(NCap++);
+  }
+ else
+  if( usePrimalSimplex )
+   for( arcPType *arc = arcsP + strt ; arc < ( arcsP + stp ) ; arc++ )
+    arc->upper = *(NCap++);
+  else
+   for( arcDType *arc = arcsD + strt ; arc < ( arcsD + stp ) ; arc++ )
+    arc->upper = *(NCap++); 
+
+ if( Senstv && (status != kUnSolved ) ) {
+  if( usePrimalSimplex ) {
+   for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++)
+    #if( QUADRATICCOST )
+     if( GT( arc->flow , arc->upper , EpsFlw ) ) 
+      arc->flow = arc->upper;
+    #else
+     if( GT(arc->flow , arc->upper , EpsFlw ) || 
+	 ( ( arc->ident == AT_UPPER ) &&
+	   ( ! ETZ( arc->flow - arc->upper , EpsFlw ) ) ) )
+      arc->flow = arc->upper;
+    #endif
+
+   CreateInitialPModifiedBalanceVector();
+   PostPVisit( dummyRootP );
+   BalanceFlow( dummyRootP );
+   ComputePotential( dummyRootP );
+   }
+  else {
+   #if( QUADRATICCOST == 0 )
+    for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ )
+     if( ( GT( arc->flow , arc->upper , EpsFlw ) && ( arc->ident != BASIC ) ) ||
+	 ( ( arc->ident == AT_UPPER ) &&
+	   ( ! ETZ( arc->flow - arc->upper , EpsFlw ) ) ) )
+      arc->flow = arc->upper;
+
+    CreateInitialDModifiedBalanceVector();
+    PostDVisit( dummyRootD );
+    ComputePotential( dummyRootD );
+   #endif
+   }
+  }
+ else
+  status = kUnSolved;
+
+ }  // end( MCFSimplex::ChgUCaps )
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::ChgUCap( Index arc , cFNumber NCap )
+{
+ if( arc >= m )
+  return;
+
+ if( usePrimalSimplex )
+  ( arcsP + arc )->upper = NCap;
+ else
+  ( arcsD + arc )->upper = NCap;
+
+ if( Senstv && (status != kUnSolved ) ) {
+  if( usePrimalSimplex ) {
+   #if( QUADRATICCOST )
+    if( GT( ( arcsP + arc )->flow , ( arcsP + arc )->upper , EpsFlw ) ) 
+     ( arcsP + arc )->flow = ( arcsP + arc )->upper;
+   #else
+    if( GT( ( arcsP + arc )->flow , ( arcsP + arc )->upper , EpsFlw ) ||
+	( ( ( arcsP + arc )->ident == AT_UPPER ) &&
+	  ( ! ETZ( ( arcsP + arc )->flow - ( arcsP + arc )->upper , EpsFlw ) ) ) )
+     ( arcsP + arc )->flow = ( arcsP + arc )->upper;
+   #endif
+
+   CreateInitialPModifiedBalanceVector();
+   PostPVisit( dummyRootP );
+   BalanceFlow( dummyRootP );
+   ComputePotential( dummyRootP );
+   }
+  else {
+   #if( QUADRATICCOST == 0 )
+    if( ( GT( ( arcsD + arc )->flow , ( arcsD + arc )->upper , EpsFlw ) &&
+	  ( ( ( arcsD + arc )->ident != BASIC ) ) ) || 
+	( ( ( arcsD + arc )->ident == AT_UPPER ) &&
+	  ( ! ETZ( ( arcsD + arc )->flow - ( arcsD + arc )->upper , EpsFlw ) ) ) ) {
+     ( arcsD + arc )->flow = ( arcsD + arc )->upper;
+     ( arcsD + arc )->ident = AT_UPPER;
+     }
+
+    CreateInitialDModifiedBalanceVector();
+    PostDVisit( dummyRootD );
+    ComputePotential( dummyRootD );
+   #endif
+   }
+  }
+ else
+  status = kUnSolved;
+
+ }  // end( MCFSimplex::ChgUCap )
+
+/*-------------------------------------------------------------------------*/
+
+bool MCFSimplex::IsClosedArc( cIndex name )
+{
+ if( name >= m )
+  return( false );
+
+ #if( QUADRATICCOST )
+  return( ( arcsP + name )->cost == Inf<CNumber>() );
+ #else
+  if( usePrimalSimplex )
+   return( ( ( arcsP + name )->ident < BASIC) );
+  else
+   return( ( ( arcsD + name )->ident < BASIC) );
+ #endif
+ }
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::CloseArc( cIndex name )
+{
+ if( name >= m )
+  return;
+
+ if( usePrimalSimplex ) {
+  arcPType *arc = arcsP + name;
+  #if( QUADRATICCOST )
+   if( arc->cost == Inf<CNumber>() )
+    return;
+
+   arc->cost = Inf<CNumber>();
+  #else
+   if( arc->ident < BASIC )
+    return;
+
+   arc->ident = CLOSED;
+  #endif
+
+  arc->flow = 0;
+
+  if( Senstv && ( status != kUnSolved ) ) {
+   nodePType *node = NULL;
+   if( (arc->tail)->enteringTArc == arc )
+    node = arc->tail;
+
+   if( (arc->head)->enteringTArc == arc )
+    node = arc->head;
+
+   if( node ) {
+    node->enteringTArc = dummyArcsP + ( node - nodesP );
+    nodePType *last = CutAndUpdateSubtree( node , -node->subTreeLevel + 1 );
+    PasteSubtree( node , last , dummyRootP );
+    node->enteringTArc = dummyArcsP + ( node - nodesP );
+    }
+
+   CreateInitialPModifiedBalanceVector();
+   PostPVisit(dummyRootP);
+   BalanceFlow(dummyRootP);                
+   ComputePotential(dummyRootP);
+   }
+  else
+   status = kUnSolved;
+  }
+ else {
+  #if( QUADRATICCOST == 0 )
+   arcDType *arc = arcsD + name;
+   if( arc->ident < BASIC )
+    return;
+
+   arc->flow = 0;
+   arc->ident = CLOSED;
+
+   if( Senstv && ( status != kUnSolved ) ) {
+    nodeDType *node = NULL;
+    if( ( arc->tail )->enteringTArc == arc)
+     node = arc->tail;
+
+    if( ( arc->head )->enteringTArc == arc )
+     node = arc->head;
+
+    if( node ) {
+     node->enteringTArc = dummyArcsD + ( node - nodesD );
+     nodeDType *last = CutAndUpdateSubtree( node , -node->subTreeLevel + 1 );
+     PasteSubtree( node , last , dummyRootD );
+     node->enteringTArc = dummyArcsD + ( node - nodesD );
+     ComputePotential( dummyRootD );
+
+     for( arcDType *a = arcsD ; a != stopArcsD ; a++ )
+      if( a->ident > BASIC )
+       if( GTZ( ReductCost( a ) , EpsCst ) ) {
+	a->flow = 0;
+	a->ident = AT_LOWER;
+        }
+       else {
+	a->flow = a->upper; 
+	a->ident = AT_UPPER;
+        }
+     }
+
+    CreateInitialDModifiedBalanceVector();
+    PostDVisit(dummyRootD);
+    ComputePotential(dummyRootD);
+    }
+   else
+    status = kUnSolved;
+  #endif
+  }
+ }  // end( MCFSimplex::CloseArc )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::DelNode( cIndex name )
+{
+ if( name >= n )
+  return;
+
+ if( usePrimalSimplex ) {
+  nodePType *node = nodesP + name;
+  nodePType *last = CutAndUpdateSubtree(node, -node->subTreeLevel);
+  nodePType *n = node->nextInT;
+  while( n ) {
+   if( n->subTreeLevel == 1 ) 
+    n->enteringTArc = dummyArcsP + ( n - nodesP );
+
+   n = n->nextInT;
+   }
+
+  PasteSubtree( node , last , dummyRootP );
+  n = node->nextInT;
+  dummyRootP->nextInT = n;
+  n->prevInT = dummyRootP;
+  
+  for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++ ) {
+   if( ( ( arc->tail ) == node) || ( ( arc->head ) == node ) ) {
+    arc->flow = 0;
+    #if( QUADRATICCOST )
+     arc->cost = Inf<CNumber>();
+    #else
+     arc->ident = CLOSED;
+    #endif
+    }
+   }
+
+  for( arcPType *arc = dummyArcsP ; arc != stopDummyP ; arc++ ) {
+   if( ( ( arc->tail ) == node ) || ( ( arc->head ) == node ) ) {
+    arc->flow = 0;
+    #if( QUADRATICCOST )
+     arc->cost = Inf<CNumber>();
+    #else
+     arc->ident = CLOSED;
+    #endif
+    }
+   }
+
+  CreateInitialPModifiedBalanceVector();
+  PostPVisit( dummyRootP );
+  BalanceFlow( dummyRootP );
+  ComputePotential( dummyRootP );
+  }
+ else {
+  #if( QUADRATICCOST == 0 )
+   nodeDType *node = nodesD + name;
+   nodeDType *last = CutAndUpdateSubtree( node , -node->subTreeLevel );
+   nodeDType *n = node->nextInT;
+   while( n ) {
+    if( n->subTreeLevel == 1 )
+     n->enteringTArc = dummyArcsD + ( n - nodesD );
+
+    n = n->nextInT;
+    }
+
+   PasteSubtree( node , last , dummyRootD );
+   n = node->nextInT;
+   dummyRootD->nextInT = n;
+   n->prevInT = dummyRootD;
+
+   for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ )
+    if( ( ( arc->tail ) == node) || ( ( arc->head ) == node ) ) {
+     arc->flow = 0;
+     arc->ident = CLOSED;
+     }
+
+   for( arcDType *arc = dummyArcsD ; arc != stopDummyD ; arc++ )
+    if( ( ( arc->tail ) == node ) || ( ( arc->head ) == node ) ) {
+     arc->flow = 0;
+     arc->ident = CLOSED;
+     }
+
+   CreateInitialDModifiedBalanceVector();
+   PostDVisit( dummyRootD );
+   ComputePotential( dummyRootD );
+  #endif
+  }
+ }  // end( MCFSimplex::DelNode )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::OpenArc( cIndex name )
+{
+ if( name >= m )
+  return;
+
+ if( usePrimalSimplex ) {
+  /* Quadratic case is not implemented for a theory bug.
+     Infact a closed arc in the quadratic case has its cost fixed to infinity, and
+     it's impossible to restore the old value. */
+
+  #if( QUADRATICCOST == 0 )
+   arcPType *arc = arcsP + name;
+   if( arc->ident == CLOSED ) {
+    arc->ident = AT_LOWER;
+    arc->flow = 0; 
+    }
+  #endif
+  }
+ else {
+  #if( QUADRATICCOST == 0 )
+   arcDType *arc = arcsD + name;
+   if( arc->ident == CLOSED ) {
+    if( GTZ( ReductCost( arc ) , EpsCst ) )
+     arc->ident = AT_LOWER;
+    else {
+     arc->ident = AT_UPPER;
+     arc->flow = arc->upper;
+     if( Senstv && ( status != kUnSolved ) ) {
+      CreateInitialDModifiedBalanceVector();
+      PostDVisit( dummyRootD );
+      }
+     else
+      status = kUnSolved;
+     }
+    }
+  #endif
+  }
+ }  // end( MCFSimplex:OpenArc )
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::Index MCFSimplex::AddNode( cFNumber aDfct )
+{
+ if( n >= nmax )
+  return( Inf<Index>() );        
+
+ n++;
+ if( usePrimalSimplex ) {
+  nodePType *newNode = nodesP + n - 1;
+  stopArcsP->tail = newNode;
+  stopArcsP->head = dummyRootP;
+  stopArcsP->upper = Inf<FNumber>();
+  stopArcsP->flow = 0;
+  stopArcsP->cost = Inf<CNumber>();
+  #if( QUADRATICCOST )
+   stopArcsP->quadraticCost = 0;
+  #else
+   stopArcsP->ident = BASIC;
+  #endif
+  stopArcsP++;
+  newNode->balance = aDfct;
+  newNode->prevInT = dummyRootP;
+  newNode->nextInT = dummyRootP->nextInT;
+  (dummyRootP->nextInT)->prevInT = newNode;
+  dummyRootP->nextInT = newNode;
+  newNode->enteringTArc = stopArcsP--;
+  newNode->potential = 0;
+  #if(QUADRATICCOST)
+   newNode->sumQuadratic = 0;
+  #endif
+  }
+ else {
+  #if( QUADRATICCOST == 0 )
+   nodeDType *newNode = nodesD + n - 1;
+   stopArcsD->tail = newNode;
+   stopArcsD->head = dummyRootD;
+   stopArcsD->upper = 0;
+   stopArcsD->flow = 0;
+   stopArcsD->cost = Inf<CNumber>();
+   stopArcsD->ident = BASIC;
+   newNode->balance = aDfct;
+   newNode->prevInT = dummyRootD;
+   newNode->nextInT = dummyRootD->nextInT;
+   (dummyRootD->nextInT)->prevInT = newNode;
+   dummyRootD->nextInT = newNode;
+   newNode->enteringTArc = stopArcsD;
+   newNode->potential = 0;
+   newNode->firstFs = stopArcsD;
+   newNode->firstBs = NULL;
+   stopArcsD->nextFs = NULL;
+   stopArcsD->nextBs = dummyRootD->firstBs;
+   dummyRootD->firstBs = stopArcsD;
+   stopArcsD++;
+  #endif
+  }
+
+ return( n );
+
+ }  // end( MCFSimplex::AddNode )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::ChangeArc( cIndex name , cIndex nSN , cIndex nEN )
+{
+ if( name >= m )
+  return;
+
+ CloseArc( name );
+
+ if( usePrimalSimplex ) {
+  if( nSN <= n )
+   (arcsP + name)->tail = (nodesP + nSN + USENAME0 - 1);
+  if( nEN <= n )
+   (arcsP + name)->head = (nodesP + nEN + USENAME0 - 1);
+  }
+ else {
+  if( nSN <= n )
+   (arcsD + name)->tail = (nodesD + nSN + USENAME0 - 1);
+  if( nEN <= n )
+   (arcsD + name)->head = (nodesD + nEN + USENAME0 - 1);
+  }
+
+ OpenArc( name );
+
+ }  // end( MCFSimplex::ChangeArc )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::DelArc( cIndex name )
+{
+ if( name >= m )
+  return;
+
+ if( usePrimalSimplex ) {
+  arcPType *arc = arcsP + name;
+  #if( QUADRATICCOST )
+   if( arc->upper == -Inf<FNumber>() )
+    return;
+
+   if( arc->cost < Inf<CNumber>() )
+  #else
+   if( arc->cost == DELETED )
+    return;
+
+   if( arc->cost >= BASIC )
+  #endif
+    CloseArc( name );
+
+  #if( QUADRATICCOST )
+   arc->upper = -Inf<FNumber>();
+  #else
+   arc->ident = DELETED;
+  #endif
+  }
+ else {
+  #if( QUADRATICCOST == 0 )
+   arcDType *arc = arcsD + name;
+   if( arc->cost == DELETED )
+    return;
+
+   if( arc->cost >= BASIC )
+    CloseArc( name );
+
+   arc->ident = DELETED;
+  #endif
+  }
+ }  // end( MCFSimplex::DelArc )
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::Index MCFSimplex::AddArc( cIndex Start , cIndex End ,
+				      cFNumber aU , cCNumber aC ) 
+{
+ if( usePrimalSimplex ) {
+  arcPType *arc = arcsP;
+  #if( QUADRATICCOST )
+   while( ( arc->upper != -Inf<FNumber>() ) && ( arc != stopArcsP ) )
+  #else
+   while( ( arc->ident > DELETED ) && ( arc != stopArcsP ) )
+  #endif
+    arc++;
+
+  if( arc == stopArcsP ) {
+   if( m >= mmax )
+    return( Inf<Index>() );
+
+   m++;
+   stopArcsP++;
+   }
+
+  Index pos = ( arc - arcsP ) + 1;
+  arc->tail = nodesP + Start + USENAME0 - 1;
+  arc->head = nodesP + End + USENAME0 - 1;
+  arc->upper = aU;
+  arc->cost = aC;
+  arc->flow = 0;
+  #if( QUADRATICCOST )
+   arc->quadraticCost = 0;
+  #else
+   arc->ident = AT_LOWER;
+  #endif
+  ComputePotential( dummyRootP );
+  return( pos );
+  }
+ else {
+  #if( QUADRATICCOST == 0 )
+   arcDType *arc = arcsD;
+   while( ( arc->ident > DELETED ) && ( arc != stopArcsD ) )
+    arc++;
+
+   if( arc == stopArcsD ) {
+    if( m >= mmax )
+     return( Inf<Index>() );
+
+    m++;
+    stopArcsD++;
+    }
+
+   Index pos = ( arc - arcsD ) + 1;
+   arc->tail = nodesD + Start + USENAME0 - 1;
+   arc->head = nodesD + End + USENAME0 - 1;
+   arc->upper = aU;
+   arc->cost = aC;
+   if( GEZ( ReductCost( arc ) , EpsCst ) ) {
+    arc->flow = 0;
+    arc->ident = AT_LOWER;
+    if( Senstv && ( status != kUnSolved ) ) {
+     CreateInitialDModifiedBalanceVector();
+     PostDVisit( dummyRootD );
+     ComputePotential( dummyRootD );
+     }
+    else
+     status = kUnSolved;
+    }
+   else {
+    arc->flow = arc->upper;
+    arc->ident = AT_UPPER;
+    if( Senstv && ( status != kUnSolved ) ) {
+     CreateInitialDModifiedBalanceVector();
+     PostDVisit( dummyRootD );
+     ComputePotential( dummyRootD );
+     }
+    else
+     status = kUnSolved;
+    }
+
+   ComputePotential( dummyRootD );
+   return( pos );
+  #endif
+  }
+ }  // end( MCFSimplex::AddArc )
+
+/*--------------------------------------------------------------------------*/
+
+bool MCFSimplex::IsDeletedArc( cIndex name )
+{
+ if( name >= m )
+  return( false );
+
+ #if( QUADRATICCOST )
+  return( ( ( arcsP + name )->upper == -Inf<FNumber>() ) );
+ #else
+  if( usePrimalSimplex )
+   return( ( arcsP + name )->ident == DELETED );
+  else
+   return( ( arcsD + name )->ident == DELETED );
+ #endif
+ }
+
+/*--------------------------------------------------------------------------*/
+/*------------------------------ DESTRUCTOR --------------------------------*/
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::~MCFSimplex()
+{
+ MemDeAllocCandidateList();
+ MemDeAlloc(true);
+ MemDeAlloc(false);
+ }
+
+/*--------------------------------------------------------------------------*/
+/*---------------------------- PRIVATE METHODS -----------------------------*/
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MemAlloc( void )
+{
+ if( usePrimalSimplex )        {
+  nodesP = new nodePType[ nmax + 1 ];   // array of nodes
+  arcsP = new arcPType[ mmax + nmax ];  // array of arcs
+  dummyArcsP = arcsP + mmax;            // artificial arcs are in the last
+                                        // nmax positions of the array arcs[]
+  }
+ else {
+  nodesD = new nodeDType[ nmax + 1 ];   // array of nodes
+  arcsD = new arcDType[ mmax + nmax ];  // array of arcs
+  dummyArcsD = arcsD + mmax;            // artificial arcs are in the last nmax
+                                        // positions of the array arcs[]
+  }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MemDeAlloc( bool whatDeAlloc )
+{
+ if( whatDeAlloc ) {
+  delete[] nodesP;
+  delete[] arcsP;
+  nodesP = NULL;
+  arcsP = NULL;
+  }
+ else {
+  delete[] nodesD;
+  delete[] arcsD;
+  nodesD = NULL;
+  arcsD = NULL;
+ }
+ MemDeAllocCandidateList( );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MemAllocCandidateList( void )
+{
+ if( usePrimalSimplex ) {
+  if( m < 10000 ) {
+   numCandidateList = PRIMAL_LOW_NUM_CANDIDATE_LIST; 
+   hotListSize = PRIMAL_LOW_HOT_LIST_SIZE;
+   }
+  else
+   if( m > 100000 ) {
+    numCandidateList = PRIMAL_HIGH_NUM_CANDIDATE_LIST; 
+    hotListSize = PRIMAL_HIGH_HOT_LIST_SIZE ;
+    }
+   else {
+    numCandidateList = PRIMAL_MEDIUM_NUM_CANDIDATE_LIST; 
+    hotListSize = PRIMAL_MEDIUM_HOT_LIST_SIZE;
+    }
+
+  #if( QUADRATICCOST )
+   int coef = 1;
+   // If the number of the arcs is more than 10000, numCandidateList and hotListSize 
+   // are increased to improve the performance of the Quadratic Primal Simplex
+   if( m > 10000 )
+    coef = 10;
+
+   numCandidateList = numCandidateList * coef;
+   hotListSize = hotListSize * coef;
+  #endif
+
+  if( forcedNumCandidateList > 0 )
+   numCandidateList = forcedNumCandidateList;
+
+  if( forcedHotListSize > 0 )
+   hotListSize = forcedHotListSize;
+
+  candP = new primalCandidType[ hotListSize + numCandidateList + 1 ];
+  }
+ else {
+  if( n < 10000 ) {
+   numCandidateList = DUAL_LOW_NUM_CANDIDATE_LIST;
+   hotListSize = DUAL_LOW_HOT_LIST_SIZE;
+   }
+  else {
+   numCandidateList = DUAL_HIGH_NUM_CANDIDATE_LIST;
+   hotListSize = DUAL_HIGH_HOT_LIST_SIZE;
+   }
+
+  if( forcedNumCandidateList > 0 )
+   numCandidateList = forcedNumCandidateList;
+
+  if( forcedHotListSize > 0 )
+   hotListSize = forcedHotListSize;
+
+  candD = new dualCandidType[ hotListSize + numCandidateList + 1 ];
+  }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::MemDeAllocCandidateList( void )
+{
+ delete[] candP;
+ candP = NULL;
+ delete[] candD; 
+ candD = NULL;
+}
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateInitialPrimalBase( void )
+{
+ arcPType *arc;
+ nodePType *node;
+ for( arc = arcsP ; arc != stopArcsP ; arc++ ) {  // initialize real arcs
+  arc->flow = 0;
+  #if( QUADRATICCOST == 0 )
+   arc->ident = AT_LOWER;
+  #endif
+  }
+
+ for( arc = dummyArcsP ; arc != stopDummyP ; arc++ ) {  // initialize dummy arcs
+  node = nodesP + ( arc - dummyArcsP );
+  if( node->balance > 0 ) {  // sink nodes 
+   arc->tail = dummyRootP;
+   arc->head = node;
+   arc->flow = node->balance;
+   }
+  else {  // source nodes or transit node
+   arc->tail = node;
+   arc->head = dummyRootP;
+   arc->flow = -node->balance;
+   }
+
+  arc->cost = MAX_ART_COST;
+  #if( QUADRATICCOST )
+   arc->quadraticCost = 0; 
+  #else
+   arc->ident = BASIC;
+  #endif
+  arc->upper = Inf<FNumber>();
+  }
+
+ dummyRootP->balance = 0;
+ dummyRootP->prevInT = NULL;
+ dummyRootP->nextInT = nodesP;
+ dummyRootP->enteringTArc = NULL;
+ #if( QUADRATICCOST )
+  dummyRootP->sumQuadratic = 0;
+ #endif
+ dummyRootP->potential = MAX_ART_COST;
+ dummyRootP->subTreeLevel = 0;
+ for( node = nodesP ; node != stopNodesP ; node++) {  // initialize nodes
+  node->prevInT = node - 1;
+  node->nextInT = node + 1;
+  node->enteringTArc = dummyArcsP + (node - nodesP);
+  #if( QUADRATICCOST )
+   node->sumQuadratic = (node->enteringTArc)->quadraticCost;
+  #endif
+  if( node->balance > 0 )  // sink nodes
+   node->potential = 2 * MAX_ART_COST;
+  else                     // source nodes or transit node
+   node->potential = 0;
+
+  node->subTreeLevel = 1;
+  }
+
+ nodesP->prevInT = dummyRootP;
+ ( nodesP + n - 1 )->nextInT = NULL;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateInitialDualBase( void )
+{
+ arcDType *arc;
+ nodeDType *node;
+ for( arc = dummyArcsD ; arc != stopDummyD ; arc++ ) {  // initialize dummy arcs
+  node = nodesD + ( arc - dummyArcsD );
+  arc->tail = node;
+  arc->head = dummyRootD;
+  arc->flow = -node->balance;
+  arc->cost = MAX_ART_COST;
+  #if( QUADRATICCOST )
+   arc->quadraticCost = 0;
+  #else
+   arc->ident = BASIC;
+  #endif
+  arc->upper = 0;
+  }
+
+ for( arc = arcsD ; arc != stopArcsD ; arc++ ) {  // initialize real arcs
+  if( GTZ( arc->cost , EpsCst ) ) {
+   arc->flow = 0;
+   #if( QUADRATICCOST == 0 )
+    arc->ident = AT_LOWER;
+   #endif
+   }
+  else {
+   #if( QUADRATICCOST == 0 )
+    arc->ident = AT_UPPER;
+   #endif
+   arc->flow = arc->upper;
+   ( dummyArcsD + ( ( arc->tail ) - nodesD ) )->flow =
+     ( dummyArcsD + ( ( arc->tail ) - nodesD ) )->flow - arc->upper;
+
+   ( dummyArcsD + ( ( arc->head ) - nodesD ) )->flow =
+    ( dummyArcsD + ( ( arc->head ) - nodesD ) )->flow + arc->upper;
+   }
+  }
+
+ dummyRootD->balance = 0;
+ dummyRootD->prevInT = NULL;
+ dummyRootD->nextInT = nodesD;
+ dummyRootD->enteringTArc = NULL;
+ #if( QUADRATICCOST )
+  dummyRootD->sumQuadratic = 0;
+ #endif
+ dummyRootD->potential = MAX_ART_COST;
+ dummyRootD->subTreeLevel = 0;
+ for( node = nodesD ; node != stopNodesD ; node++ ) {  // initialize nodes
+  node->prevInT = node - 1;
+  node->nextInT = node + 1;
+  node->enteringTArc = dummyArcsD + ( node - nodesD );
+  #if( QUADRATICCOST )
+   node->sumQuadratic = ( node->enteringTArc )->quadraticCost;
+  #endif
+  node->potential = 0;
+  node->subTreeLevel = 1;
+  node->whenInT2 = 0;
+  }
+
+ nodesD->prevInT = dummyRootD;
+ ( nodesD + n - 1 )->nextInT = NULL;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateAdditionalDualStructures( void )
+{
+ // this method creates, in a Dual context, the Backward Star and the
+ // Forward Star of every node
+
+ for( nodeDType *node = nodesD ; node != stopNodesD ; node++) {
+  // initialize nodes
+  node->firstBs = NULL;
+  node->firstFs = NULL;
+  node->numArcs = 0;
+  }
+
+ dummyRootD->firstBs = NULL;
+ dummyRootD->firstFs = NULL;
+ dummyRootD->numArcs = 0;
+ for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ ) {
+  // initialize real arcs
+  arc->nextFs = ( arc->tail )->firstFs;
+  ( arc->tail )->firstFs = arc;
+  arc->nextBs = ( arc->head )->firstBs;
+  ( arc->head )->firstBs = arc;
+  ( arc->tail )->numArcs++;
+  ( arc->head )->numArcs++;
+  }
+
+ ResetWhenInT2();
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PrimalSimplex( void )
+{
+ #if( UNIPI_PRIMAL_INITIAL_SHOW == 0 )
+  #if( UNIPI_PRIMAL_ITER_SHOW == 0 )
+   #if( UNIPI_PRIMAL_FINAL_SHOW == 0 ) 
+  //cout << endl;
+   #endif
+  #endif
+ #endif
+ #if( UNIPI_PRIMAL_INITIAL_SHOW )
+  cout << endl;
+  for( int t = 0; t < 3; t++ )
+   cout << "\t";
+  cout << "PRIMALE MCFSimplex: ARCHI E NODI ALL' INIZIO" << endl;
+  ShowSituation( 3 );
+ #endif
+ #if( QUADRATICCOST )
+  #if( LIMITATEPRECISION )
+   foValue = GetFO();
+   int cont = 0;
+  #endif
+ #endif
+
+ status = kUnSolved;
+ if( pricingRule != kCandidateListPivot )
+  arcToStartP = arcsP;
+
+ iterator = 0;  // setting the initial arc for the Dantzig or First Elibigle Rule
+
+ arcPType *enteringArc;
+ arcPType *leavingArc;
+ if( pricingRule == kCandidateListPivot )
+  InitializePrimalCandidateList();
+
+ while( status == kUnSolved ) {
+  iterator++;
+  switch( pricingRule ) {
+  case( kDantzig ):
+   enteringArc = RuleDantzig();
+   break;
+  case( kFirstEligibleArc ):
+   enteringArc = PRuleFirstEligibleArc();
+   break;
+  case( kCandidateListPivot ):
+   enteringArc = RulePrimalCandidateListPivot();
+   break;
+  }
+
+ #if( QUADRATICCOST )
+  #if( LIMITATEPRECISION )
+   /* In the quadratic case with LIMITATEPRECISION == 1, the entering arcs
+      are selected according to a thresold.
+      This thresold is definited according to the old f.o. value.
+      If Primal Simplex doesn't find an entering arc, it calculates again
+      the f.o. value, and try again. */
+   if( enteringArc == NULL ) {
+    foValue = GetFO();
+    switch( pricingRule ) {
+    case( kDantzig ):
+     enteringArc = RuleDantzig();
+     break;
+    case( kFirstEligibleArc ):
+     enteringArc = PRuleFirstEligibleArc();
+     break;
+    case( kCandidateListPivot ):
+     enteringArc = RulePrimalCandidateListPivot();
+     break;
+     }
+    }
+  #endif
+ #endif
+
+ if( pricingRule != kCandidateListPivot ) {
+  // in every iteration the algorithm changes the initial arc for
+  // Dantzig and First Eligible Rule.
+  arcToStartP++;
+  if( arcToStartP == stopArcsP )
+   arcToStartP = arcsP;
+  }
+
+ if( enteringArc ) {
+  arcPType *arc;
+  nodePType *k1;
+  nodePType *k2;
+  /* If the reduced cost of the entering arc is > 0, the Primal Simplex
+     pushes flow in the cycle determinated by T and entering arc for decreases 
+     flow in the entering arc: in the linear case entering arc's flow goes to 0,
+     in the quadratic case it decreases while it's possibile.
+     If the reduced cost of the entering arc is < 0, the Primal Simplex pushes
+     flow in the cycle  determinated by T and entering arc for increases flow
+     in the entering arc: in the linear case entering arc's flow goes to upper
+     bound, in the quadratic case it increases while it's possibile.  */
+
+  #if( QUADRATICCOST )
+   FONumber t;
+   FONumber theta; 
+   FONumber deltaFO;
+   FNumber theta2;
+   CNumber Q = ( enteringArc->tail )->sumQuadratic +
+               ( enteringArc->head )->sumQuadratic + enteringArc->quadraticCost;
+   // Q is the sum of the quadratic coefficient in the cycle determinated by T
+   // and entering arc.
+   FONumber rc = ReductCost( enteringArc );
+   if( ETZ( Q, EpsCst ) )
+    theta = Inf<FNumber>();  // This value will be certainly decreased 
+   else
+    theta = ABS( rc / Q );
+    // This is the best theta value (with best f.o. value decrease) 
+
+   leavingArc = enteringArc;
+   nodePType *cycleRoot; // The radix of the cycle determinated by T and entering arc.
+   if( GTZ( rc , EpsCst ) ) {
+  #else
+   FNumber t;
+   FNumber theta;
+   if( enteringArc->ident == AT_UPPER ) {
+  #endif
+    /*  Primal Simplex increases or decreases entering arc's flow.
+	"theta" is a positive value.
+	For this reason the algorithm uses two nodes ("k1" and "k2") to push
+	flow ("theta") from "k1" to "k2".
+	According to entering arc's reduct cost, the algorithm determinates
+	"k1" and "k2" */
+
+    k1 = enteringArc->head;
+    k2 = enteringArc->tail;
+    #if( QUADRATICCOST )
+     theta = min( theta , enteringArc->flow );
+     // The best value for theta is compared with the entering arc's bound
+     theta2 = - theta;
+    #else
+     theta = enteringArc->flow;
+    #endif
+    }
+   else {
+    k1 = enteringArc->tail;
+    k2 = enteringArc->head;        
+    #if( QUADRATICCOST )
+     theta = min( theta , enteringArc->upper - enteringArc->flow );
+     // The best value for theta is compared with the entering arc's bound
+     theta2 = theta;
+    #else
+     theta = enteringArc->upper - enteringArc->flow;
+    #endif
+    }
+
+   nodePType *memK1 = k1;
+   nodePType *memK2 = k2;
+   leavingArc = NULL;
+   #if( QUADRATICCOST )
+    #if( LIMITATEPRECISION )
+     deltaFO = rc * theta2 + Q * theta2 * theta2 / 2;
+    #endif
+    bool leavingReducesFlow = GTZ( rc , EpsCst );
+   #else
+    bool leavingReducesFlow = GTZ( ReductCost( enteringArc ) , EpsCst );
+   #endif
+   // Compute "theta", find outgoing arc and "root" of the cycle
+   bool leave;
+   // Actual "theta" is compared with the bounds of the other cycle's arcs
+   while( k1 != k2 ) {
+    if( k1->subTreeLevel > k2->subTreeLevel ) {
+     arc = k1->enteringTArc;
+     if( arc->tail != k1 ) {
+      t = arc->upper - arc->flow;
+      leave = false;
+      }
+     else {
+      t = arc->flow;
+      leave = true;
+      }
+
+     if( t < theta ) {
+      // The algorithm has found a possible leaving arc
+      theta = t;
+      leavingArc = arc;
+      leavingReducesFlow = leave;
+      // If "leavingReducesFlow" == true, if this arc is selected to exit the base, 
+      // it decreases its flow
+      }
+
+     k1 = Father( k1 , arc );
+     }
+    else {
+     arc = k2->enteringTArc;
+     if( arc->tail == k2 ) {
+      t = arc->upper - arc->flow;
+      leave = false;
+      }
+     else {
+      t = arc->flow;
+      leave = true;
+      }
+
+     if( t <= theta ) {
+      // The algorithm has found a possible leaving arc
+      theta = t;
+      leavingArc = arc;
+      leavingReducesFlow = leave;
+      // If "leavingReducesFlow" == true, if this arc is selected to exit the base, 
+      // it decreases its flow
+      }
+
+     k2 = Father(k2, arc);
+     }
+    }
+
+   #if( QUADRATICCOST )
+    cycleRoot = k1;
+   #endif
+
+   if( leavingArc == NULL )
+    leavingArc = enteringArc;
+
+   if( theta >= Inf<FNumber>() ) {
+    status = kUnbounded;
+    break;
+    }
+
+   // Update flow with "theta"
+   k1 = memK1;
+   k2 = memK2;
+   #if( QUADRATICCOST )
+    if( enteringArc->tail == k1 )
+     theta2 = theta;
+    else
+     theta2 = -theta;
+
+    // "theta" is a positive value in every case.
+    // "theta2" is the real theta value according to the real
+    // direction of the entering arc
+    #if( LIMITATEPRECISION )
+     deltaFO = rc * theta2 + Q * theta2 * theta2 / 2;
+     // The decrease of the f.o. value in the quadratic case
+    #endif
+   #endif
+
+     if( ! ETZ(theta , EpsFlw ) ) {
+      if( enteringArc->tail == k1 )
+       enteringArc->flow = enteringArc->flow + theta;
+      else
+       enteringArc->flow = enteringArc->flow - theta;
+
+      while( k1 != k2 ) {
+       if( k1->subTreeLevel > k2->subTreeLevel ) {
+	arc = k1->enteringTArc;
+	if( arc->tail != k1 )
+	 arc->flow = arc->flow + theta;
+	else
+	 arc->flow = arc->flow - theta;
+
+	k1 = Father(k1, k1->enteringTArc);
+        }
+       else {
+	arc = k2->enteringTArc;
+	if( arc->tail == k2 )
+	 arc->flow = arc->flow + theta;
+	else
+	 arc->flow = arc->flow - theta;
+
+	k2 = Father(k2, k2->enteringTArc);
+        }
+       }
+      }
+
+     if( enteringArc != leavingArc ) {
+      bool leavingBringFlowInT2 = ( leavingReducesFlow == 
+	( ( leavingArc->tail )->subTreeLevel > ( leavingArc->head )->subTreeLevel ) );
+      // "leavingBringFlowInT2" == true if leaving arc brings flow to the subtree T2
+      if( leavingBringFlowInT2 == ( memK1 == enteringArc->tail ) ) {
+       k2 = enteringArc->tail;
+       k1 = enteringArc->head;
+       }
+      else {
+       k2 = enteringArc->head; 
+       k1 = enteringArc->tail;
+       }
+      }
+
+     #if( QUADRATICCOST == 0 )
+      if( leavingReducesFlow )
+       leavingArc->ident = AT_LOWER;
+      else
+       leavingArc->ident = AT_UPPER;
+
+      if( leavingArc != enteringArc ) {
+       enteringArc->ident = BASIC;
+       nodePType *h1;
+       nodePType *h2;
+       // "h1" is the node in the leaving arc with smaller tree's level 
+       if( ( leavingArc->tail )->subTreeLevel < ( leavingArc->head )->subTreeLevel ) {
+	h1 = leavingArc->tail;
+	h2 = leavingArc->head;
+        }
+       else {
+	h1 = leavingArc->head;
+	h2 = leavingArc->tail;
+        }
+
+       UpdateT(leavingArc, enteringArc, h1, h2, k1, k2);
+       // Update potential of the subtree T2
+       k2 = enteringArc->head;
+       CNumber delta = ReductCost(enteringArc);
+       if( ( enteringArc->tail )->subTreeLevel > ( enteringArc->head )->subTreeLevel ) {
+	delta = -delta;
+	k2 = enteringArc->tail;
+        }
+
+       AddPotential( k2 , delta );
+       // In the linear case Primal Simplex only updates the potential of the nodes of
+       // subtree T2
+       }
+     #else
+      if( leavingArc != enteringArc ) {
+       nodePType *h1;
+       nodePType *h2;
+       // "h1" is the node in the leaving arc with smaller tree's level 
+       if( ( leavingArc->tail )->subTreeLevel <
+	   ( leavingArc->head )->subTreeLevel ) {
+	h1 = leavingArc->tail;
+	h2 = leavingArc->head;
+        }
+       else {
+	h1 = leavingArc->head;
+	h2 = leavingArc->tail;
+        } 
+
+       // Update the basic tree
+       UpdateT( leavingArc , enteringArc , h1 , h2 , k1 , k2 );
+       }
+
+      #if( OPTQUADRATIC )
+       nodePType *h1;
+       nodePType *h2;
+       if( ( leavingArc->tail )->subTreeLevel <
+	   ( leavingArc->head )->subTreeLevel ) {
+	h1 = leavingArc->tail;
+	h2 = leavingArc->head;
+        }
+       else {
+	h1 = leavingArc->head;
+	h2 = leavingArc->tail;
+        }
+
+       nodePType *node = h1;
+       nodePType *updateNode = h1;
+       if( h1 == cycleRoot )
+	ComputePotential( cycleRoot );
+       else {
+	while( node != cycleRoot ) {
+	 arcPType *entArc = node->enteringTArc;
+	 if( ! ETZ( entArc->quadraticCost , EpsCst ) )
+	  updateNode = node;
+
+	 node = Father( node , entArc );
+	 }
+
+	ComputePotential( updateNode );
+	node = h2;
+	updateNode = h2;
+	while( node != cycleRoot ) {
+	 arcPType *entArc = node->enteringTArc;
+	 if( ! ETZ( entArc->quadraticCost , EpsCst ) )
+	  updateNode = node;
+
+	 node = Father( node , entArc );
+	 }
+
+	ComputePotential( updateNode );
+        }
+      #else
+       // Update the potential of the node "manually"
+       ComputePotential( cycleRoot );
+      #endif
+
+      #if( LIMITATEPRECISION )
+       cont = cont + 1;
+       if( cont == recomputeFOLimits ) {
+	cont = 0;
+	foValue = GetFO();  // Calculate f.o. value manually
+        }
+       else
+	foValue = foValue + deltaFO;
+        // Calculate the f.o. value with the estimated decrease
+      #endif
+     #endif
+    }
+   else {
+    status = kOK;
+    // If one of dummy arcs has flow bigger than 0, the solution is unfeasible.
+    for( arcPType *arc = dummyArcsP ; arc != stopDummyP ; arc++ )
+     if( GTZ( arc->flow , EpsFlw ) ) 
+      status = kUnfeasible;
+    }
+
+   if( ( status == kUnSolved ) && MaxTime && MCFt ) {
+    double tu, ts;
+    TimeMCF( tu , ts );
+    if( MaxTime < tu + ts )
+     status = kStopped;
+    }
+
+   if( ( status == kUnSolved ) && MaxIter)
+    if( MaxIter < (int) iterator )
+     status = kStopped;
+
+   #if( UNIPI_PRIMAL_ITER_SHOW )
+    int it = (int) iterator;
+    if( it % UNIPI_PRIMAL_ITER_SHOW == 0 ) {        
+     cout << endl;
+     for( int t = 0; t < 3; t++ )
+      cout << "\t";
+     cout << "PRIMALE MCFSimplex: ARCHI E NODI ALLA " << it << " ITERAZIONE" << endl;
+     ShowSituation( 3 );
+     #if( FOSHOW )
+      if( (int) iterator % FOSHOW == 0 ) 
+       clog << "Iteration = " << iterator << " of = "
+        #if( LIMITATEPRECISION && QUADRATICCOST )
+	    << foValue
+        #else
+	    << GetFO()
+        #endif
+            << endl;
+     #endif
+     }
+   #endif
+  }
+
+ #if( UNIPI_PRIMAL_FINAL_SHOW )
+  cout << endl;
+  for( int t = 0; t < 3; t++ )
+   cout << "\t";
+  cout << "PRIMALE UniPi: ARCHI E NODI ALLA FINE" << endl;
+  ShowSituation( 3 );
+ #endif
+
+ }  // end( PrimalSimplex )
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::DualSimplex( void )
+{
+ #if( UNIPI_PRIMAL_INITIAL_SHOW == 0 )
+  #if( UNIPI_PRIMAL_ITER_SHOW == 0 )
+   #if( UNIPI_PRIMAL_FINAL_SHOW == 0 ) 
+    cout << endl;
+   #endif
+  #endif
+ #endif
+ #if( UNIPI_DUAL_INITIAL_SHOW )
+  cout << endl;
+  for( int t = 0; t < 3; t++ )
+   cout << "\t";
+  cout << "DUALE MCFSimplex: ARCHI E NODI ALL' INIZIO" << endl;
+  ShowSituation( 3 );
+ #endif
+
+ if( pricingRule != kCandidateListPivot )
+  arcToStartD = arcsD;
+
+ iterator = 0;
+ arcDType *enteringArc;
+ arcDType *leavingArc;
+ if( pricingRule == kCandidateListPivot )
+  InitializeDualCandidateList();
+
+ status = kUnSolved;
+ while( status == kUnSolved ) {
+  iterator++;
+  if( iterator == Inf<iteratorType>() ) {
+   ResetWhenInT2(); // Restore to 0 every nodes' field "whenInT2"
+   iterator = 1;
+   }
+
+  switch( pricingRule ) {
+  case( kDantzig ):
+   leavingArc = DRuleFirstEligibleArc();
+   break; // si esegue cmq FEA
+  case( kFirstEligibleArc ):
+   leavingArc = DRuleFirstEligibleArc();
+   break;
+  case( kCandidateListPivot ):
+   leavingArc = RuleDualCandidateListPivot();
+   break;
+   }
+
+  if( pricingRule != kCandidateListPivot ) {
+   arcToStartD++;
+   if( arcToStartD == stopArcsD )
+    arcToStartD = dummyArcsD;
+
+   if( arcToStartD == stopDummyD )
+    arcToStartD = arcsD;
+   // Setting the initial arc for the Dantzig or First Elibigle Rule
+   }
+
+  if( leavingArc ) {
+   bool leavingArcInL = false;
+   bool leavingArcFromT1toT2;
+   if( LTZ( leavingArc->flow , EpsFlw ) )
+    leavingArcInL = true;
+
+   nodeDType *h1;
+   nodeDType *h2;
+   if( ( leavingArc->tail )->subTreeLevel <
+       ( leavingArc->head )->subTreeLevel ) {
+    h1 = leavingArc->tail;
+    h2 = leavingArc->head;
+    leavingArcFromT1toT2 = true;
+    }
+   else {
+    h1 = leavingArc->head;
+    h2 = leavingArc->tail;
+    leavingArcFromT1toT2 = false;
+    }
+
+   Index numOfT2Arcs = 0;
+   int level = h2->subTreeLevel;
+   nodeDType *node = h2;
+   node->whenInT2 = iterator;
+   nodeDType *lastNodeOfT2 = h2;
+   numOfT2Arcs = node->numArcs;
+   while( node->nextInT && ( ( node->nextInT )->subTreeLevel > level ) ) {
+    node = node->nextInT;
+    lastNodeOfT2 = node;
+    numOfT2Arcs = numOfT2Arcs + node->numArcs;
+    node->whenInT2 = iterator;
+    }
+
+   /* The Dual Simplex has determinated the leaving arc, and so the
+      subtrees T1 and T2. Dual Simplex scans T2 to fix the fields "whenInT2"
+      of T2's nodes to the iteration value, and counts the entering and
+      outgoing arcs from these nodes. According to this number, it decides
+      to scan the Backward and Forward of the subtree (T1 or T2) with the
+      minor number of entering/outgoing arcs from its nodes. */
+   enteringArc = NULL;
+   bool lv = ( leavingArcFromT1toT2 == leavingArcInL );
+   CNumber maxRc = Inf<CNumber>();
+   //Search arc in the Forward Star and Backward Star of nodes of T1
+   if( numOfT2Arcs > m ) {
+    // Dual Simplex starts from the node which follows the dummy root.
+    node = dummyRootD->nextInT;
+    bool fine = false;
+    while( fine == false ) {
+     /* If node is the root of subtree T2, Dual Simplex jumps to the node
+	(if exists) which follows the last node of T2 */
+     if( node == h2 )
+      if( lastNodeOfT2->nextInT )
+       node = lastNodeOfT2->nextInT;
+      else
+       break;
+
+     // Search arc in the Backward Star of nodes of T1
+     arcDType *arc = node->firstBs;
+     while( arc ) {
+      if( ( arc->tail )->whenInT2 == iterator ) {
+       // Evaluate each arc from T2 to T1 which isn't in T
+       if( arc->ident == AT_LOWER ) {
+	if( lv ) {
+	 CNumber rc = ABS( ReductCost( arc ) );
+	 if( LT( rc , maxRc , EpsCst ) ) {
+	  enteringArc = arc;
+	  maxRc = rc;
+	  /* If arc is appropriate to enter in T and its reduct cost is 0, 
+             search is ended: this is the arc which enters in T */
+	  if( ETZ( rc , EpsCst) ) {
+	   fine = true;
+	   break;
+	   }
+	  }
+	 }
+        }
+
+       if( arc->ident == AT_UPPER ) {
+	if( ! lv ) {
+	 CNumber rc = ABS( ReductCost( arc ) );
+	 if( LT( rc , maxRc , EpsCst ) ) {
+	  enteringArc = arc;
+	  maxRc = rc;
+	  /* If arc is appropriate to enter in T and its reduct cost is 0, 
+             search is ended: this is the arc which enters in T */
+
+	  if( ETZ( rc , EpsCst ) ) {
+	   fine = true;
+	   break;
+	   }
+	  }
+	 }
+        }
+       }
+
+      arc = arc->nextBs;
+      }
+
+     // Search arc in the Forward Star of nodes of T1
+     arc = node->firstFs;
+     while( arc ) {
+      if( ( arc->head )->whenInT2 == iterator ) {
+       // Evaluate each arc from T1 to T2 which isn't in T
+       if( arc->ident == AT_LOWER ) {
+	if( ! lv ) {
+	 CNumber rc = ABS( ReductCost( arc ) );
+	 if( LT( rc , maxRc , EpsCst ) ) {
+	  enteringArc = arc;
+	  maxRc = rc;
+	  /* If arc is appropriate to enter in T and its reduct cost is 0, 
+	     search is ended: this is the arc which enters in T */
+	  if( ETZ( rc , EpsCst ) ) {
+	   fine = true;
+	   break;
+	   }
+	  }
+	 }
+        }
+
+       if( arc->ident == AT_UPPER ) {
+	if( lv ) {
+	 CNumber rc = ABS( ReductCost( arc ) );
+	 if( LT( rc , maxRc , EpsCst ) ) {
+	  enteringArc = arc;
+	  maxRc = rc;
+	  /* If arc is appropriate to enter in T and its reduct cost is 0, 
+	     search is ended: this is the arc which enters in T */
+	  if( ETZ( rc , EpsCst ) ) {
+	   fine = true;
+	   break;
+	   }
+	  }
+	 }
+        }
+       }
+
+      arc = arc->nextFs;
+      }
+
+     node = node->nextInT;
+     if( node == NULL )
+      fine = true;
+     }
+    }
+   // Search arc in the Forward Star and Backward Star of nodes of T2
+   else {
+    node = h2;
+    bool fine = false;
+    while( fine == false ) {
+     // Search arc in the Backward Star of nodes of T2
+     arcDType *arc = node->firstBs;
+     CNumber rc;
+     while( arc ) {
+      if( ( arc->tail )->whenInT2 != iterator ) {
+       // Evaluate each arc from T1 to T2 which isn't in T
+       if( arc->ident == AT_LOWER ) {
+	if( ! lv ) {
+	 rc = ABS( ReductCost( arc ) );
+	 if( LT( rc , maxRc , EpsCst ) ) {
+	  enteringArc = arc;
+	  maxRc = rc;
+	  /* If arc is appropriate to enter in T and its reduct cost is 0, 
+	     search is ended: this is the arc which enters in T */
+	  if( ETZ( rc , EpsCst ) ) {
+	   fine = true;
+	   break;
+	   }
+	  }
+	 }
+        }
+
+       if( arc->ident == AT_UPPER ) {
+	if( lv ) {
+	 rc = ABS( ReductCost( arc ) );
+	 if( LT( rc , maxRc , EpsCst ) ) {
+	  enteringArc = arc;
+	  maxRc = rc;
+	  /* If arc is appropriate to enter in T and its reduct cost is 0, 
+	     search is ended: this is the arc which enters in T */
+	  if( ETZ( rc , EpsCst ) ) {
+	   fine = true;
+	   break;
+	   }
+	  }
+	 }
+        }
+       }
+
+      arc = arc->nextBs;
+      }
+
+     // Search arc in the Forward Star of nodes of T2
+     arc = node->firstFs;
+     while( arc ) {
+      if( ( arc->head )->whenInT2 != iterator ) {
+       // Evaluate each arc from T2 to T1 which isn't in T
+       if( arc->ident == AT_LOWER ) {
+	if( lv ) {
+	 rc = ABS( ReductCost( arc ) );
+	 if( LT( rc , maxRc , EpsCst ) ) {
+	  enteringArc = arc;
+	  maxRc = rc;
+	  /* If arc is appropriate to enter in T and its reduct cost is 0, 
+	     search is ended: this is the arc which enters in T */
+	  if( ETZ( rc , EpsCst ) ) {
+	   fine = true;
+	   break;
+	   }
+	  }
+	 }
+        }
+
+       if( arc->ident == AT_UPPER ) {
+	if( ! lv ) {
+	 rc = ABS( ReductCost( arc ) );
+	 if( LT( rc , maxRc , EpsCst ) ) {
+	  enteringArc = arc;
+	  maxRc = rc;
+	  /* If arc is appropriate to enter in T and its reduct cost is 0, 
+	     search is ended: this is the arc which enters in T */
+	  if( ETZ( rc , EpsCst) ) {
+	   fine = true;
+	   break;
+	   }
+	  }
+	 }
+        }
+       }
+
+      arc = arc->nextFs;
+      }
+
+     if( node == lastNodeOfT2 )
+      fine = true;
+     else
+      node = node->nextInT;
+
+     }
+    }
+
+   if( enteringArc ) {
+    FNumber theta = -leavingArc->flow;
+    if( GTZ( leavingArc->flow , EpsFlw ) )
+     theta = leavingArc->flow - leavingArc->upper;
+    // Initial value of theta is the infeasibility of the leaving arc
+
+    FNumber t;
+    nodeDType *k1;
+    nodeDType *k2;
+    /* if entering arc is in U, Dual Simplex pushs flow in the cycle 
+       determinated by T and entering arc for decreases flow in the entering arc:
+       if entering arc is in L, Dual Simplex pushs flow in the cycle 
+       determinated by T and entering arc for increases flow in the entering arc:
+       Dual Simplex increases or decreases entering arc's flow.
+       theta is a positive value.
+       For this reason the algorithm uses two nodes (k1 and k2) to push flow
+       (theta) from k1 to k2. According to entering arc's reduct cost, the
+       algorithm determinates k1 and k2 */
+
+    if( enteringArc->ident == AT_UPPER ) {
+     k1 = enteringArc->head;
+     k2 = enteringArc->tail;
+     }
+    else {
+     k1 = enteringArc->tail;
+     k2 = enteringArc->head;
+     }
+
+    nodeDType *memK1 = k1;
+    nodeDType *memK2 = k2;
+    arcDType *arc;                                
+    k1 = memK1;
+    k2 = memK2;
+    // Update the flow
+    while( k1 != k2 ) {
+     if( k1->subTreeLevel > k2->subTreeLevel ) {
+      arc = k1->enteringTArc;
+      if( arc->tail != k1 )
+       arc->flow = arc->flow + theta;
+      else
+       arc->flow = arc->flow - theta;
+
+      k1 = Father(k1, k1->enteringTArc);
+      }
+     else {
+      arc = k2->enteringTArc;
+      if( arc->tail == k2 )
+       arc->flow = arc->flow + theta;
+      else
+       arc->flow = arc->flow - theta;
+
+      k2 = Father( k2 , k2->enteringTArc );
+      }
+     }
+
+    if(leavingArcInL )
+     leavingArc->ident = AT_LOWER;
+    else
+     leavingArc->ident = AT_UPPER;
+
+    bool leavingBringFlowInT2 = ( leavingArcInL == 
+	 ( ( leavingArc->tail )->subTreeLevel >
+	   ( leavingArc->head )->subTreeLevel ) );
+    // leavingBringFlowInT2 == true if leaving arc brings flow to the subtree T2
+    if( leavingBringFlowInT2 != ( memK1 == enteringArc->tail ) ) {
+     k2 = enteringArc->tail;
+     k1 = enteringArc->head;
+     }
+    else {
+     k2 = enteringArc->head; 
+     k1 = enteringArc->tail;
+     }
+
+    if( enteringArc->ident == AT_LOWER )
+     enteringArc->flow = enteringArc->flow + theta;
+    else
+     enteringArc->flow = enteringArc->flow - theta;
+
+    enteringArc->ident = BASIC;
+    UpdateT( leavingArc , enteringArc , h1 , h2 , k1 , k2 );
+    // update potential of the subtree T2
+    k2 = enteringArc->head;
+    CNumber delta = ReductCost( enteringArc );
+    if( ( enteringArc->tail) ->subTreeLevel >
+	( enteringArc->head )->subTreeLevel ) {
+     delta = -delta;
+     k2 = enteringArc->tail;
+     }
+
+    // Dual Simplex only updates the potential of the T2's nodes
+    AddPotential( k2 , delta );
+    }
+   else
+    status = kUnfeasible;
+    /* If Dual Simplex finds a leaving arc but it doesn't find an entering arc,
+       the algorithm stops. At this point Dual Simplex has an unfeasible primal
+       solution. */
+   }
+  else {
+   status = kOK;
+   // If one of dummy arcs has flow different than 0, the solution is unfeasible.
+   for( arcDType *arc = dummyArcsD ; arc != stopDummyD ; arc++ )
+    if( ! ETZ( arc->flow , EpsFlw ) ) {
+     status = kUnfeasible;
+     break;
+     }
+   }
+
+  if( ( status == kUnSolved ) && MaxTime ) {
+   double tu, ts;
+   TimeMCF( tu , ts );
+   if( MaxTime < tu + ts )
+    status = kStopped;
+   }
+
+  if( ( status == kUnSolved ) && MaxIter && MCFt )
+   if( MaxIter < (int) iterator )
+    status = kStopped;
+
+  #if( UNIPI_DUAL_ITER_SHOW )
+   if( (int) iterator % UNIPI_DUAL_ITER_SHOW == 0 ) {
+    cout << endl;
+    for( int t = 0; t < 3; t++ )
+     cout << "\t";
+    cout << "DUALE MCFSimplex: ARCHI E NODI ALLA " << iterator << " ITERAZIONE"
+	 << endl;
+    ShowSituation( 3 );
+    #if( FOSHOW )
+     cout << "of = " << GetFO() << endl;
+    #endif
+    }
+  #endif
+  }
+
+ #if( UNIPI_DUAL_ITER_SHOW )
+  int it = (int) iterator;
+  if( it % UNIPI_DUAL_ITER_SHOW == 0 ) {        
+   cout << endl;
+   for( int t = 0; t < 3; t++ )
+    cout << "\t";
+   cout << "DUALE MCFSimplex: ARCHI E NODI ALLA " << iterator << " ITERAZIONE"
+	<< endl;
+   Showsituation( 3 );
+   }
+ #endif
+
+ }  // end( DualSimplex )
+
+/*--------------------------------------------------------------------------*/
+
+template<class N, class A>
+void MCFSimplex::UpdateT( A *h , A *k , N *h1 , N *h2 , N *k1 , N *k2 )
+{
+ /* In subtree T2 there is a path from node h2 (deepest node of the leaving
+    arc h and root of T2) to node k2 (deepest node of the leaving arc h and
+    coming root of T2). With the update of T, this path will be overturned:
+    node k2 will become the root of T2...
+    The subtree T2 must be reordered and the field "subTreeLevel", which
+    represents the depth in T of every node, of every T2's nodes is changed.
+    Variable delta represents the increase of "subTreeLevel" value for node
+    k2 and its descendants: probably this value is a negative value. */
+
+ int delta = (k1->subTreeLevel) + 1 - (k2->subTreeLevel);
+ N *root = k2;
+ N *dad;
+ 
+ /*To reorder T2, the method analyses every nodes of the path h2->k2, starting
+   from k2. For every node, it moves the node's subtree from its original
+   position to a new appropriate position. In particular k2 and its subtree
+   (k2 is the new root of T2, so the first nodes of the new T2) will be moved
+   next to node k1 (new father of k2), the next subtree will be moved beside
+   the last node of k2's subtree....
+   "previousNode" represents the node where the new subtree will be moved
+   beside in this iterative action. At the start "previousNode" is the node
+   k1 (T2 will be at the right of k1). */
+
+ N *previousNode = k1;
+ N *lastNode;
+ /* "arc1" is the entering arc in T (passed by parameters).
+    For every analysed node of path h2->k2, the method changes
+    "enteringTArc" but it must remember the old arc, which will be the
+    "enteringTArc" of the next analysed node.
+    At the start "arc1" is k (the new "enteringTArc" of k2). */
+
+ A *arc1 = k;
+ A *arc2;
+ bool fine = false;
+ while( fine == false ) {
+  // If analysed node is h2, this is the last iteration
+  if( root == h2 )
+   fine = true;
+
+  dad = Father( root , root->enteringTArc );
+  // Cut the root's subtree from T and update the "subLevelTree" of its nodes
+  lastNode = CutAndUpdateSubtree( root , delta );
+  // Paste the root's subtree in the right position;
+  PasteSubtree( root , lastNode , previousNode );
+  // In the next iteration the subtree will be moved beside the last node of
+  // the actual analysed subtree.
+  previousNode = lastNode;
+  /* The increase of the subtree in the next iteration is different from
+     the actual increase: in particual the increase increases itself (+2 
+     at every iteration). */
+  delta = delta + 2; 
+  /* A this point "enteringTArc" of actual root is stored in "arc2" and
+     changed; then "arc1" and "root" are changed. */
+  arc2 = root->enteringTArc;
+  root->enteringTArc = arc1;
+  arc1 = arc2;
+  root = dad;
+  } 
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N>
+N* MCFSimplex::CutAndUpdateSubtree( N *root , int delta )
+{
+ int level = root->subTreeLevel;
+ N *node = root;
+ // The root of this subtree is passed by parameters, the last node is searched.
+ while ( ( node->nextInT ) && ( ( node->nextInT )->subTreeLevel > level ) ) {
+  node = node->nextInT;
+  // The "subTreeLevel" of every nodes of subtree is updated
+  node->subTreeLevel = node->subTreeLevel + delta;
+  }
+
+ root->subTreeLevel = root->subTreeLevel + delta;
+ /* The 2 neighbouring nodes of the subtree (the node at the left of the root
+    and the node at the right of the last node) is joined. */
+
+ if( root->prevInT )
+  ( root->prevInT )->nextInT = node->nextInT;
+ if( node->nextInT )
+  ( node->nextInT )->prevInT = root->prevInT;
+
+ return( node );  // the method returns the last node of the subtree
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N>
+void MCFSimplex::PasteSubtree( N *root , N *lastNode , N *previousNode )
+{
+ /* The method inserts subtree ("root" and "lastNode" are the extremity of the
+    subtree) after "previousNode". The method starts to identify the next node
+    of "previousNode" ("nextNode"), so it joins "root" with "previousNode" and
+    "lastNode" with "nextNode" (if exists). */
+
+ N *nextNode = previousNode->nextInT;
+ root->prevInT = previousNode;
+ previousNode->nextInT = root;
+ lastNode->nextInT = nextNode;
+ if( nextNode )
+  nextNode->prevInT = lastNode;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcPType* MCFSimplex::RuleDantzig( void )
+{
+ arcPType *arc = arcToStartP;
+ arcPType *enteringArc = NULL;
+ #if( QUADRATICCOST )
+  /* In the quadratic case used type for reduct cost is FONumber.
+     Value "lim" is the fixed thresold for the decrease of the f.o. value */
+  FONumber lim = EpsOpt * foValue / n;
+  FONumber RC;
+  FONumber maxValue = 0;
+ #else
+  CNumber RC;
+  CNumber maxValue = 0;
+ #endif
+
+ do {
+  // The method analyses every arc 
+  #if( QUADRATICCOST )
+   RC = ReductCost( arc );
+   FNumber theta;
+   /* If reduct cost of arc is lower than 0, the flow of the arc must increase.
+      If reduct cost of arc is bigger than 0, the flow of the arc must decrease.
+      "theta" is the difference between lower (upper) bound and the actual flow.
+      */
+
+   if( LTZ( RC , EpsCst ) )
+    theta = arc->upper - arc->flow;
+
+   if( GTZ( RC , EpsCst ) ) 
+    theta = -arc->flow;
+
+   // If it's possible to increase (or decrease) the flow in this arc
+   if( ! ETZ( theta , EpsFlw ) ) {
+    /* "Q" is the sum of the quadratic coefficient of the arc belonging the
+       T path from tail's arc to head's arc
+       "Q" is always bigger than 0 or equals to 0.
+       If "Q" > 0, the value - RC / Q is the increase (decrease) of the flow
+       with the best decrease of f.o. value.
+       - RC/ Q must be compare with "theta" to avoid that the best increase
+       (decrease) of the flow violates the bounds of the arc.
+       This confront determines "theta". */
+
+    CNumber Q = ( arc->tail )->sumQuadratic + ( arc->head )->sumQuadratic +
+                arc->quadraticCost;
+
+    if( GTZ( Q , EpsCst ) )
+     if( GTZ( theta , EpsFlw ) )
+      theta = min( theta , - RC / Q );        
+     else
+      theta = max( theta , - RC / Q );        
+
+    /* Calculate the estimate decrease of the f.o. value if this arc is
+       selected and flow is increased (decreased) by "theta" */
+
+    CNumber deltaFO = RC * theta + Q * theta * theta / 2;
+    // if deltaFO < 0 this arc is appropriate; if deltaFO is lower than
+    // old decrease value, arc is the best arc.
+
+    if( deltaFO < maxValue ) {
+     maxValue = deltaFO;
+     enteringArc = arc;
+     }
+    }
+  #else
+   if( arc->ident > BASIC ) {
+    RC = ReductCost( arc );
+
+    if( ( LTZ( RC , EpsCst ) && ( arc->ident == AT_LOWER ) ) || 
+	( GTZ( RC , EpsCst ) && ( arc->ident == AT_UPPER ) ) ) {
+
+     if( RC < 0 )
+      RC = -RC;
+
+     if( RC > maxValue ) {
+      maxValue = RC;
+      enteringArc = arc;
+      }
+     }
+    }
+  #endif
+
+  arc++;
+  if( arc == stopArcsP )
+   arc = arcsP;
+
+  } while( arc != arcToStartP );
+
+ #if( ( LIMITATEPRECISION ) && ( QUADRATICCOST ) )
+  if( -maxValue <= lim )
+   enteringArc = NULL;
+ #endif
+
+ return( enteringArc );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcPType* MCFSimplex::PRuleFirstEligibleArc( void )
+{
+ arcPType *arc = arcToStartP;
+ arcPType *enteringArc = NULL;
+ #if( QUADRATICCOST )
+  FONumber RC;
+ #else
+  CNumber RC;
+ #endif
+
+ do {
+  #if( QUADRATICCOST )
+   // In this method the "decrease f.o. value" criterion is not used.
+   RC = ReductCost( arc );
+   if( ( LTZ( RC , EpsCst ) && LT( arc->flow , arc->upper , EpsFlw ) ) || 
+       ( GTZ( RC , EpsCst ) && GTZ( arc->flow , EpsFlw ) ) )
+    enteringArc = arc;
+  #else
+   if( arc->ident > BASIC ) {
+    RC = ReductCost( arc );
+    if( ( LTZ( RC , EpsCst ) && ( arc->ident == AT_LOWER ) ) ||
+	( GTZ( RC , EpsCst ) && ( arc->ident == AT_UPPER ) ) )
+     enteringArc = arc;
+    }
+  #endif
+
+   arc++;
+   if( arc == stopArcsP )
+    arc = dummyArcsP;
+   if( arc == stopDummyP )
+    arc = arcsP;
+
+  } while( ( enteringArc == NULL ) && ( arc != arcToStartP ) );
+
+ return( enteringArc );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcDType* MCFSimplex::DRuleFirstEligibleArc( void )
+{
+ arcDType *arc = arcToStartD;
+ arcDType *leavingArc = NULL;
+ do {
+  if( GT( arc->flow , arc->upper , EpsFlw ) || LTZ( arc->flow , EpsFlw ) )
+   leavingArc = arc;
+
+  arc++;
+  if( arc == stopArcsD )
+   arc = dummyArcsD;
+  if( arc == stopDummyD )
+   arc = arcsD;
+
+  } while( ( leavingArc == NULL ) && ( arc != arcToStartD ) );
+
+ return(leavingArc);
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcPType* MCFSimplex::RulePrimalCandidateListPivot( void )
+{
+ Index next = 0;
+ Index i;
+ Index minimeValue;
+ if( hotListSize < tempCandidateListSize )
+  minimeValue = hotListSize;
+ else
+  minimeValue = tempCandidateListSize;
+
+ #if( QUADRATICCOST )
+  // Check if the left arcs in the list continue to violate the dual condition
+  for( i = 2 ; i <= minimeValue ; i++ ) {
+   arcPType *arc = candP[ i ].arc;
+   FONumber red_cost = ReductCost( arc );
+   FNumber theta = 0;
+   /* If reduct cost of arc is lower than 0, the flow of the arc must increase.
+      If reduct cost of arc is bigger than 0, the flow of the arc must decrease.
+      "theta" is the difference between lower (upper) bound and the actual flow.
+      */
+
+   if( LTZ( red_cost , EpsCst ) )
+    theta = arc->upper - arc->flow;
+   else
+    if( GTZ( red_cost , EpsCst ) )
+     theta = - arc->flow;
+
+   // if it's possible to increase (or decrease) the flow in this arc
+   if( theta != 0 ) {
+    /* "Q" is the sum of the quadratic coefficient of the arc belonging the T
+       path from tail's arc to head's arc
+       "Q" is always bigger than 0 or equals to 0.
+       If "Q" > 0, the value - RC / Q is the increase (decrease) of the flow
+       with the best decrease of f.o. value.
+       - RC/ Q must be compare with "theta" to avoid that the best increase
+       (decrease) of the flow violates the bounds of the arc.
+       This confront determines "theta". */
+
+    CNumber Q = ( arc->tail )->sumQuadratic + ( arc->head )->sumQuadratic +
+                arc->quadraticCost;
+
+    if( GTZ( Q , EpsCst  ) )
+     if( GTZ( theta , EpsFlw ) )
+      theta = min( theta , - red_cost / Q );
+     else
+      theta = max( theta , - red_cost / Q );
+
+    /* Calculate the estimate decrease of the f.o. value if this arc is
+       selected and flow is increased (decreased) by "theta" */
+
+    CNumber deltaFO = red_cost * theta + Q * theta * theta / 2;
+    #if( LIMITATEPRECISION )
+     // if deltaFO < 0 this arc is appropriate; if deltaFO is lower than
+     // old decrease value, arc is the best arc.
+     if( - deltaFO > ( EpsOpt * foValue / n ) ) {
+    #else
+     if( LTZ( deltaFO , EpsCst ) ) {
+    #endif
+      next++;
+      candP[ next ].arc = arc;
+      candP[ next ].absRC = -deltaFO;
+      }
+     }
+    }
+
+   tempCandidateListSize = next;
+   Index oldGroupPos = groupPos;
+   // Search other arcs to fill the list
+   do {
+    arcPType *arc;
+    for( arc = arcsP + groupPos ; arc < stopArcsP ; arc += numGroup ) {
+     FONumber red_cost = ReductCost( arc );
+     FNumber theta = 0;
+     /* If reduced cost is lower than 0, the flow of the arc must increase.
+	If reduced cost is larger than 0, the flow of the arc must decrease.
+	"theta" is the difference between lower (upper) bound and the actual
+	flow. */
+
+     if( LTZ( red_cost , EpsCst ) ) 
+      theta = arc->upper - arc->flow;
+     else
+      if( GTZ( red_cost , EpsCst ) ) 
+       theta = - arc->flow;
+
+     // if it's possible to increase (or decrease) the flow in this arc
+     if( theta != 0 ) {
+      /* "Q" is the sum of the quadratic coefficient of the arc belonging the T
+	 path from tail's arc to head's arc
+	 "Q" is always bigger than 0 or equals to 0.
+	 If "Q" > 0, the value - RC / Q is the increase (decrease) of the flow
+	 with the best decrease of f.o. value.
+	 - RC/ Q must be compare with "theta" to avoid that the best increase
+	 (decrease) of the flow violates the bounds of the arc.
+	 This confront determines "theta". */
+      CNumber Q = ( arc->tail )->sumQuadratic + ( arc->head )->sumQuadratic +
+                  arc->quadraticCost;
+
+      if( GTZ( Q , EpsCst  ) )
+       if( GTZ( theta , EpsFlw ) )
+	theta = min( theta , - red_cost / Q );        
+       else
+	theta = max( theta , - red_cost / Q );        
+
+      /* Calculate the estimate decrease of the f.o. value if this arc is
+	 selected and flow is increased (decreased) by "theta" */
+
+      CNumber deltaFO = red_cost * theta + Q * theta * theta / 2;
+      #if( LIMITATEPRECISION )
+       // if deltaFO < 0 this arc is appropriate; if deltaFO is lower than
+       // old decrease value, arc is the best arc.
+       if( -deltaFO > ( EpsOpt * foValue / n ) ) {
+      #else
+       if( LTZ( deltaFO , EpsCst ) ) {
+      #endif
+	tempCandidateListSize++;
+	candP[ tempCandidateListSize ].arc = arc;
+	candP[ tempCandidateListSize ].absRC = -deltaFO;
+        }
+       }
+      }
+
+     groupPos++;
+     if( groupPos == numGroup )
+      groupPos = 0;
+
+     } while( ( tempCandidateListSize < hotListSize ) &&
+	      ( groupPos != oldGroupPos ) );
+ #else
+  // Check if the left arcs in the list continue to violate the dual condition
+  for( i = 2 ; i <= minimeValue ; i++ ) {
+   arcPType *arc = candP[i].arc;
+   CNumber red_cost = ReductCost( arc );
+
+   if( ( LTZ( red_cost , EpsCst ) && ( arc->ident == AT_LOWER ) ) ||
+       ( GTZ( red_cost , EpsCst ) && ( arc->ident == AT_UPPER ) ) ) {
+    next++;
+    candP[ next ].arc = arc;
+    candP[ next ].absRC = ABS( red_cost );
+    }
+   }
+
+  tempCandidateListSize = next;
+  Index oldGroupPos = groupPos;
+  // Search other arcs to fill the list
+  do {
+   arcPType *arc;
+   for( arc = arcsP + groupPos ; arc < stopArcsP ; arc += numGroup ) {
+    if( arc->ident == AT_LOWER ) {
+     CNumber red_cost = ReductCost( arc );
+     if( LTZ( red_cost , EpsCst ) ) {
+      tempCandidateListSize++;
+      candP[ tempCandidateListSize ].arc = arc;
+      candP[ tempCandidateListSize ].absRC = ABS( red_cost );
+      }
+     }
+    else
+     if( arc->ident == AT_UPPER ) {
+      CNumber red_cost = ReductCost( arc );
+      if( GTZ( red_cost , EpsCst ) ) {
+       tempCandidateListSize++;
+       candP[ tempCandidateListSize ].arc = arc;
+       candP[ tempCandidateListSize ].absRC = ABS( red_cost );
+       }
+      }
+    }
+
+   groupPos++;
+   if( groupPos == numGroup )
+    groupPos = 0;
+
+   } while( ( tempCandidateListSize < hotListSize ) && ( groupPos != oldGroupPos ) );
+ #endif
+
+ if( tempCandidateListSize ) {
+  SortPrimalCandidateList( 1 , tempCandidateListSize );
+  return( candP[ 1 ].arc );
+  }
+ else
+  return( NULL );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::InitializePrimalCandidateList( void )
+{
+ numGroup = ( ( m - 1 ) / numCandidateList ) + 1;
+ groupPos = 0;
+ tempCandidateListSize = 0;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::SortPrimalCandidateList( Index min , Index max )
+{
+ Index left = min;
+ Index right = max;
+ #if( QUADRATICCOST )
+  FONumber cut = candP[ ( left + right ) / 2 ].absRC;
+ #else
+  CNumber cut = candP[ ( left + right ) / 2 ].absRC;
+ #endif
+ do {
+  while( candP[ left ].absRC > cut) 
+   left++;
+  while( cut > candP[ right ].absRC) 
+   right--;
+
+  if( left < right )
+   Swap( candP[ left ] , candP[ right ] );
+
+  if(left <= right) {
+   left++;
+   right--;
+   }
+  } while( left <= right );
+
+ if( min < right ) 
+  SortPrimalCandidateList( min , right );
+ if( ( left < max ) && ( left <= hotListSize ) ) 
+  SortPrimalCandidateList( left , max );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcDType* MCFSimplex::RuleDualCandidateListPivot( void )
+{
+ Index next = 0;
+ // Check if the left arcs in the list continue to violate the primal condition
+ for( Index i = 2 ; ( i <= hotListSize ) && ( i <= tempCandidateListSize ) ;
+      i++ ) {
+  nodeDType *node = candD[ i ].node;
+  arcDType *arc = node->enteringTArc;
+  cFNumber flow = arc->flow;
+  if( LTZ( flow , EpsFlw ) ) {
+   next++;
+   candD[ next ].node = node;
+   candD[ next ].absInfeas = ABS( flow );
+   }
+
+  if( GT( flow , arc->upper , EpsFlw ) ) {
+   next++;
+   candD[ next ].node = node;
+   candD[ next ].absInfeas = flow - arc->upper;
+   }
+  }
+
+ tempCandidateListSize = next;
+ Index oldGroupPos = groupPos;
+ // Search other arcs to fill the list
+ do {
+  nodeDType *node = nodesD + groupPos;
+  for( node ; node < stopNodesD ; node += numGroup ) {
+   arcDType *arc = node->enteringTArc;
+   cFNumber flow = arc->flow;
+   if( LTZ( flow , EpsFlw ) ) {
+    tempCandidateListSize++;
+    candD[ tempCandidateListSize ].node = node;
+    candD[ tempCandidateListSize ].absInfeas = ABS( flow );
+    }
+
+   if( GT( flow , arc->upper , EpsFlw) ) {
+    tempCandidateListSize++;
+    candD[ tempCandidateListSize ].node = node;
+    candD[ tempCandidateListSize ].absInfeas = flow - arc->upper;
+    }
+   }
+
+  groupPos++;
+  if( groupPos == numGroup )
+   groupPos = 0;
+
+  } while( ( tempCandidateListSize < hotListSize ) &&
+	   ( groupPos != oldGroupPos ) );
+
+ if( tempCandidateListSize ) {
+  SortDualCandidateList( 1 , tempCandidateListSize );
+  return( (candD[ 1 ].node)->enteringTArc );
+  }
+ else
+  return( NULL );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::InitializeDualCandidateList( void )
+{
+ numGroup = ( ( n - 1 ) / numCandidateList ) + 1;
+ groupPos = 0;
+ tempCandidateListSize = 0;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::SortDualCandidateList(Index min, Index max)
+{
+ Index left = min;
+ Index right = max;
+ FNumber cut = candD[ ( left + right ) / 2 ].absInfeas;
+ do {
+  while( candD[ left ].absInfeas > cut )
+   left++;
+  while( cut > candD[ right ].absInfeas )
+   right--;
+  if( left < right ) 
+   Swap( candD[left ] , candD[ right ] );
+  if( left <= right) {
+   left++;
+   right--;
+   }
+  } while( left <= right );
+
+ if( min < right )
+  SortDualCandidateList( min , right );
+ if( (left < max) && ( left <= hotListSize ) ) 
+  SortDualCandidateList( left , max );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N, class RCT>
+inline void MCFSimplex::AddPotential( N *r , RCT delta )
+{
+ int level = r->subTreeLevel;
+ N *n = r;
+ 
+ do {
+  n->potential = n->potential + delta;
+  n = n->nextInT;
+  } while ( ( n ) && ( n->subTreeLevel > level ) );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N>
+inline void MCFSimplex::ComputePotential( N *r )
+{
+ N *n = r;
+ int level = r->subTreeLevel;
+ FONumber cost;
+ // If "n" is not the dummy root, the potential of "r" is computed.
+ // If "n" is the dummy root, the potential of dummy root is a constant.
+ 
+ do {
+  if( n->enteringTArc ) {
+   cost = ( n->enteringTArc )->cost;
+   #if (QUADRATICCOST)
+    // Also field "sumQuadratic" is updated
+    n->sumQuadratic = ( Father( n , n->enteringTArc ) )->sumQuadratic +
+                      ( n->enteringTArc )->quadraticCost;
+
+    if( ! ETZ( ( n->enteringTArc )->flow , EpsFlw ) )
+     cost = cost + ( ( n->enteringTArc )->quadraticCost * ( n->enteringTArc )->flow );
+   #endif
+
+   if( n == ( n->enteringTArc )->head ) 
+    n->potential = ( Father( n , n->enteringTArc ) )->potential + cost;
+   else
+    n->potential = ( Father( n , n->enteringTArc ) )->potential - cost;
+   }
+  n = n->nextInT;
+  } while( ( n ) && ( n->subTreeLevel > level ) );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateInitialPModifiedBalanceVector( void )
+{
+ int i = 0;
+ delete[] modifiedBalance;
+ modifiedBalance = new FNumber[ n ];
+ // Initialited every node's modifiedBalance to his balance
+ for ( nodePType *node = nodesP ; node != stopNodesP ; node++ ) {
+  modifiedBalance[i] = node->balance;
+  i++;
+  }
+
+ // Modify the vector according to the arcs out of base with flow non zero
+ // Scan the real arcs
+
+ for( arcPType *arc = arcsP ; arc != stopArcsP ; arc++ ) {
+  #if( QUADRATICCOST )
+   if( ( ! ETZ( arc->flow , EpsFlw ) ) &&
+       ( ( arc->tail )->enteringTArc != arc ) && 
+       ( ( arc->head )->enteringTArc != arc ) ) {
+    i = (arc->tail) - nodesP;
+    modifiedBalance[ i ] += arc->flow;
+    i = (arc->head) - nodesP;
+    modifiedBalance[ i ] -= arc->flow;
+    }
+  #else
+   if( arc->ident == AT_UPPER ) {
+    i = (arc->tail) - nodesP;
+    modifiedBalance[ i ] += arc->upper;
+    i = (arc->head) - nodesP;
+    modifiedBalance[ i ] -= arc->upper;
+    }
+  #endif
+  }
+
+ // Scan the dummy arcs
+ for( arcPType *arc = dummyArcsP ; arc != stopDummyP ; arc++ ) {
+  #if( QUADRATICCOST )
+   if ( ( ! ETZ( arc->flow , EpsFlw ) ) &&
+	( ( arc->tail )->enteringTArc != arc ) && 
+	( ( arc->head )->enteringTArc != arc ) ) {
+    i = (arc->tail) - nodesP;
+    modifiedBalance[ i ] += arc->flow;
+    i = (arc->head) - nodesP;
+    modifiedBalance[ i ] -= arc->flow;
+    }
+  #else
+   if (arc->ident == AT_UPPER) {
+    i = (arc->tail) - nodesP;
+    modifiedBalance[ i ] += arc->upper;
+    i = (arc->head) - nodesP;
+    modifiedBalance[ i ] -= arc->upper;
+   }
+  #endif
+  }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PostPVisit( nodePType *r )
+{
+ // The method controls if "r" is a leaf in T
+ bool rLeaf = false;
+ int i = r - nodesP;
+ if( r->nextInT ) 
+  if( ( r->nextInT )->subTreeLevel <= r->subTreeLevel )
+   rLeaf = true;
+  else
+   rLeaf = true;
+
+ if( rLeaf )  // If "r" is a leaf
+  if( ( r->enteringTArc)->head == r ) // If enteringTArc of "r" goes in "r"
+   ( r->enteringTArc )->flow = modifiedBalance[ i ];
+  else // If enteringTArc of "r" goes out "r"
+   ( r->enteringTArc )->flow = - modifiedBalance[ i ];
+ else { // If "r" isn't a leaf
+  nodePType *desc = r->nextInT;
+  // Call PostPVisit for every child of "r"
+  while( ( desc ) && ( desc->subTreeLevel > r->subTreeLevel ) ) {
+   if( desc->subTreeLevel - 1 == r->subTreeLevel ) { // desc is a son of r
+    PostPVisit( desc );
+
+    if( ( desc->enteringTArc )->head == r ) // enteringTArc of desc goes in r
+     modifiedBalance[ i ] -= ( desc->enteringTArc )->flow;
+    else // If enteringTArc of "desc" goes out "r"
+     modifiedBalance[ i ] += ( desc->enteringTArc )->flow;
+    }
+   desc = desc->nextInT;
+   }
+
+  if( r != dummyRootP )
+   if( ( r->enteringTArc )->head == r ) // If enteringTArc of "r" goes in "r"
+    ( r->enteringTArc )->flow = modifiedBalance[ i ];
+   else // If enteringTArc of "r" goes out "r"
+    ( r->enteringTArc )->flow = - modifiedBalance[ i ];
+  }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::BalanceFlow( nodePType *r )
+{
+ // used only by Primal Simplex to restore a primal feasible solution.
+ if( r == dummyRootP ) {
+  nodePType *node = dummyRootP->nextInT;
+  while( node ) {
+   // call this function recursively for every son of dummy root
+   if( node->subTreeLevel == 1 ) 
+    BalanceFlow( node );
+
+   node = node->nextInT;
+   }
+  }
+ else {
+  // The method controls if "r" is a leaf in T
+  bool rLeaf = false;
+  if( r->nextInT )
+   if( ( r->nextInT )->subTreeLevel <= r->subTreeLevel )
+    rLeaf = true;
+   else
+    rLeaf = true;
+
+  if( rLeaf )  // If "r" is a leaf
+   AdjustFlow( r );  // The method controls if entering basic arc in "r" is
+                     // not feasible; in case adjust its flow
+  else { // If "r" isn't a leaf
+   nodePType *node = r->nextInT;
+   // Balance the flow of every child of "r"
+   while ( ( node ) && ( node->subTreeLevel > r->subTreeLevel ) ) {
+    if( node->subTreeLevel == r->subTreeLevel + 1 ) 
+     BalanceFlow( node );
+    node = node->nextInT;
+    }
+
+   // The method controls if entering basic arc in "r" is not feasible;
+   //in case adjust its flow
+   AdjustFlow( r );
+   }
+  }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::AdjustFlow( nodePType *r )
+{
+ arcPType *arc = r->enteringTArc;
+ if( arc >= dummyArcsP ) // If entering arc of "r" is a dummy arc
+  if( LTZ( arc->flow , EpsFlw ) ) {
+   // If this dummy arc has flow < 0, the algorithm overturns the arc
+   nodePType *temp = arc->tail;
+   arc->tail = arc->head;
+   arc->head = temp;
+   arc->flow = -arc->flow;
+   }
+ else {  // If entering arc of "r" is not a dummy arc
+  bool orientationDown = ( arc->head == r );
+  FNumber delta = 0;
+  if( LTZ( arc->flow , EpsFlw ) ) { // If flow is < 0
+   delta = -arc->flow;
+   arc->flow = 0;
+   #if( ! QUADRATICCOST )
+    arc->ident = AT_LOWER;
+   #endif
+   }
+
+  if( GT( arc->flow , arc->upper , EpsFlw ) ) {
+   // If flow goes over the capacity of the arc
+   delta = arc->upper - arc->flow;
+   arc->flow = arc->upper;
+   #if( ! QUADRATICCOST )
+    arc->ident = AT_UPPER;
+   #endif
+   }
+
+  /* This arc goes out from the basis, and the relative dummy arc goes in T.
+     Then the algorithm push flow in the cycle made by the arc and some arcs
+     of T to balance the flow. */
+
+  if( ! ETZ( delta , EpsFlw ) ) {
+   nodePType *node = Father( r , arc );
+   while( node != dummyRootP ) {
+    arc = node->enteringTArc;
+    if( ( arc->head == node ) == orientationDown )
+     arc->flow += delta;
+    else
+     arc->flow -= delta;
+
+    node = Father( node , arc );
+    }
+
+   arcPType *dummy = dummyArcsP + ( r - nodesP );
+   #if( ! QUADRATICCOST )
+    dummy->ident = BASIC;
+   #endif
+
+   /* Update the structure of the tree. If entering basic arc of "r" is
+      changed, subtree of "r"is moved next dummy root. */
+   r->enteringTArc = dummy;
+   int deltaLevel = 1 - r->subTreeLevel;
+   nodePType *lastNode = CutAndUpdateSubtree( r , deltaLevel ); 
+   PasteSubtree( r , lastNode , dummyRootP );
+   if( ( dummy->head == r ) != orientationDown )
+    dummy->flow += delta;
+   else
+    dummy->flow -= delta;
+
+   if( LTZ( dummy->flow , EpsFlw ) ) {
+    nodePType *temp = dummy->tail;
+    dummy->tail = dummy->head;
+    dummy->head = temp;
+    dummy->flow = -dummy->flow;
+    }
+   }
+  }
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::CreateInitialDModifiedBalanceVector( void )
+{
+ #if( ! QUADRATICCOST )
+  int i = 0;
+  modifiedBalance = new FNumber[ n ];
+  // Initialited every node's modifiedBalance to his balance
+  for( nodeDType *node = nodesD ; node != stopNodesD ; node++ ) {
+   modifiedBalance[ i ] = node->balance;
+   i++;
+   }
+
+  // Modify the vector according to the arcs out of base with flow non zero
+  // Scan the real arcs
+  for( arcDType *arc = arcsD ; arc != stopArcsD ; arc++ )
+   if( arc->ident == AT_UPPER ) {
+    i = (arc->tail) - nodesD;
+    modifiedBalance[ i ] += arc->upper;
+    i = (arc->head) - nodesD;
+    modifiedBalance[ i ] -= arc->upper;
+    }
+
+  // Scan the dummy arcs
+  for( arcDType *arc = dummyArcsD ; arc != stopDummyD ; arc++ )
+   if( arc->ident == AT_UPPER ) {
+    i = (arc->tail) - nodesD;
+    modifiedBalance[ i ] += arc->upper;
+    i = (arc->head) - nodesD;
+    modifiedBalance[ i ] -= arc->upper;
+    }
+ #endif
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PostDVisit( nodeDType *r )
+{
+ #if( ! QUADRATICCOST )
+  // The method controls if "r" is a leaf in T
+  bool rLeaf = false;
+  int i = r - nodesD;
+  if( r->nextInT )
+   if( ( r->nextInT )->subTreeLevel <= r->subTreeLevel )
+    rLeaf = true;
+   else
+    rLeaf = true;
+
+  if( rLeaf ) // If "r" is a leaf
+   if( ( r->enteringTArc)->head == r ) // If enteringTArc of "r" goes in "r"
+    ( r->enteringTArc )->flow = modifiedBalance[ i ];
+   else // If enteringTArc of "r" goes out "r"
+    ( r->enteringTArc )->flow = - modifiedBalance[ i ];
+  else { // If "r" isn't a leaf
+   nodeDType *desc = r->nextInT;
+   // Call PostDVisit for every child of "r"
+   while( ( desc ) && ( desc->subTreeLevel > r->subTreeLevel ) ) {
+    if( desc->subTreeLevel -1 == r->subTreeLevel ) { // desc is a son of r
+     PostDVisit( desc );
+
+     if( ( desc->enteringTArc )->head == r ) // enteringTArc of desc goes in r
+      modifiedBalance[ i ] -= ( desc->enteringTArc )->flow;
+     else // If enteringTArc of "desc" goes out "r"
+      modifiedBalance[ i ] += ( desc->enteringTArc )->flow;
+     }
+
+    desc = desc->nextInT;
+    }
+
+   if( r != dummyRootD )
+    if( ( r->enteringTArc )->head == r ) // If enteringTArc of "r" goes in "r"
+     ( r->enteringTArc )->flow = modifiedBalance[ i ];
+    else // If enteringTArc of "r" goes out "r"
+     ( r->enteringTArc )->flow = - modifiedBalance[ i ];
+   }
+ #endif
+ }
+
+/*--------------------------------------------------------------------------*/
+
+inline void MCFSimplex::ResetWhenInT2( void )
+{
+ for( nodeDType *n = nodesD ; n != stopNodesD ; n++)
+  n->whenInT2 = 0;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+template<class N, class A>
+inline N* MCFSimplex::Father( N *n , A *a )
+{
+ if( a == NULL )
+  return NULL;
+
+ if( a->tail == n )
+  return( a->head );
+ else
+  return( a->tail );
+ }
+
+/*-------------------------------------------------------------------------*/
+
+inline MCFSimplex::FONumber MCFSimplex::GetFO( void )
+{
+ FONumber fo = 0;
+ if( usePrimalSimplex ) {
+  arcPType *arco;
+  for( arco = arcsP ; arco != stopArcsP ; arco++ ) {
+   #if( QUADRATICCOST ) 
+    if( ! ETZ( arco->flow , EpsFlw ) ) 
+     fo += arco->flow * ( arco->cost + arco->flow * arco->quadraticCost / 2 );
+   #else
+    if( ( arco->ident == BASIC ) || ( arco->ident == AT_UPPER ) )
+     fo += arco->cost * arco->flow;
+   #endif
+   }
+
+  for( arco = dummyArcsP ; arco != stopDummyP ; arco++ ) {
+   #if( QUADRATICCOST ) 
+    if( ! ETZ( arco->flow , EpsFlw ) )
+     fo += arco->flow * ( arco->cost + arco->flow * arco->quadraticCost / 2 );
+   #else
+    if( ( arco->ident == BASIC ) || ( arco->ident == AT_UPPER ) ) 
+     fo += arco->cost * arco->flow;
+   #endif
+   }
+  }
+ else {
+  arcDType *a;
+  for( a = arcsD ; a != stopArcsD ; a++ ) {
+   #if (QUADRATICCOST) 
+    fo += ( a->cost * a->flow ) + ( a->quadraticCost * a->flow * a->flow ) / 2;
+   #else
+    if( ( a->ident == BASIC ) || (a->ident == AT_UPPER ) ) 
+     fo += a->cost * a->flow;
+   #endif
+   }
+
+  for( a = dummyArcsD ; a != stopDummyD ; a++) {
+   #if (QUADRATICCOST) 
+    fo += ( a->cost * a->flow ) + ( a->quadraticCost * a->flow * a->flow ) / 2;
+   #else
+    if( ( a->ident == BASIC ) || ( a->ident == AT_UPPER ) ) 
+     fo += a->cost * a->flow;
+   #endif
+   }
+  }
+
+ return( fo );
+ }
+
+/*-------------------------------------------------------------------------*/
+
+void MCFSimplex::PrintPNode( nodePType *nodo )
+{
+ if( nodo )
+  if( nodo != dummyRootP )
+   cout << ( nodo - nodesP + 1 );
+  else
+   cout << "r";
+ else
+  cout << "..";
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PrintPArc( arcPType *arc )
+{
+ if( arc ) {
+  cout << "(";
+  PrintPNode( arc->tail );
+  cout << ", ";
+  PrintPNode( arc->head );
+  cout << ")";
+  }
+ else
+  cout << "..";
+}
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PrintDNode( nodeDType *nodo )
+{
+ if( nodo )
+  if( nodo != dummyRootD )
+   cout << ( nodo - nodesD + 1 );
+  else
+   cout << "r";
+ else
+  cout << "..";
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::PrintDArc( arcDType *arc )
+{
+ if( arc ) {
+  cout << "(";
+  PrintDNode( arc->tail );
+  cout << ", ";
+  PrintDNode( arc->head );
+  cout << ")";
+  }
+ else
+  cout << "..";
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::nodePType* MCFSimplex::RecoverPNode( Index ind ) 
+{
+ if( ( ind < 0 ) || ( ind > n ) )
+  return( NULL );
+ if( ind )
+  return( nodesP + ind - 1 );
+ else
+  return( dummyRootP );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcPType* MCFSimplex::RecoverPArc( nodePType *tail ,
+					       nodePType *head )
+{
+ if( ( tail == NULL ) || ( head == NULL ) )
+  return( NULL );
+
+ arcPType *arc = arcsP;
+ while( ( arc->tail != tail ) || ( arc->head != head ) ) {
+  arc++;
+  if( arc == stopArcsP )
+   arc = dummyArcsP;
+   if( arc == stopDummyP )
+    return( NULL );
+  }
+
+ return( arc );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::nodeDType* MCFSimplex::RecoverDNode( Index ind )
+{
+ if( ( ind < 0 ) || ( ind > n ) ) 
+  return( NULL );
+
+ if( ind )
+  return( nodesD + ind - 1 );
+ else
+  return( dummyRootD );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+MCFSimplex::arcDType* MCFSimplex::RecoverDArc( nodeDType *tail ,
+					       nodeDType *head )
+{
+ if( ( tail == NULL ) || ( head == NULL ) ) 
+  return( NULL );
+
+ arcDType *arc = arcsD;
+ while( ( arc->tail != tail ) || ( arc->head != head ) ) {
+  arc++;
+  if( arc == stopArcsD )
+   arc = dummyArcsD;
+  if( arc == stopDummyD )
+   return( NULL );
+  }
+
+ return( arc );
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::infoPNode( nodePType *node , int tab )
+{
+ for( int t = 0 ; t < tab ; t++ )
+  cout << "\t";
+ cout << "Nodo ";
+ PrintPNode( node );
+ cout << ": b = " << node->balance << " y = " << node->potential << endl;
+ #if( UNIPI_VIS_NODE_BASIC_ARC )
+  cout << ": TArc=";
+  PrintPArc( node->enteringTArc );
+  cout << endl;
+ #endif
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::infoPArc( arcPType *arc , int ind , int tab )
+{
+ for( int t = 0 ; t < tab ; t++ )
+  cout << "\t";
+ cout << "Arco ";
+ PrintPArc( arc );
+ cout << ": x = " << arc->flow;
+ #if( UNIPI_VIS_ARC_UPPER )
+  cout << " u = " << arc->upper;
+ #endif
+ #if( UNIPI_VIS_ARC_COST )
+  cout << " c = " << arc->cost;
+ #endif
+ #if( QUADRATICCOST )
+  #if( UNIPI_VIS_ARC_Q_COST )
+   cout << " q = " << arc->quadraticCost;
+  #endif
+  cout << endl;
+  for( int t = 0 ; t < tab ; t++ )
+   cout << "\t";
+  #if( UNIPI_VIS_ARC_REDUCT_COST )
+   cout << " rc = " << MCFGetRC( ind );
+  #endif
+ #else
+  cout << endl;
+  for( int t = 0 ; t < tab ; t++ )
+   cout << "\t";
+  #if( UNIPI_VIS_ARC_REDUCT_COST )
+   cout << " rc = " << MCFGetRC( ind );
+  #endif
+  #if( UNIPI_VIS_ARC_STATE )
+   switch( arc->ident ) {
+   case( BASIC ):    cout << " in T"; break;
+   case( AT_LOWER ): cout << " in L"; break;
+   case( AT_UPPER ): cout << " in U"; break;
+   case( DELETED ):  cout << " canceled"; break;
+   case( CLOSED ):   cout << " closed";
+   }
+  #endif
+ #endif
+ cout << endl;
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::infoDNode( nodeDType *node , int tab )
+{
+ for( int t = 0 ; t < tab; t++ )
+  cout << "\t";
+ cout << "Nodo ";
+ PrintDNode( node );
+ cout << ": b = " << node->balance << " y = " << node->potential;
+ #if( UNIPI_VIS_NODE_BASIC_ARC )    
+  cout << ": TArc=";
+  PrintDArc( node->enteringTArc );
+  cout << endl;
+ #endif
+ }
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::infoDArc( arcDType *arc , int ind , int tab )
+{
+ for( int t = 0 ; t < tab ; t++ )
+  cout << "\t";
+ cout << "Arco ";
+ PrintDArc( arc );
+ cout << " x = " << arc->flow;
+ #if( UNIPI_VIS_ARC_UPPER )
+  cout << " u = " << arc->upper;
+ #endif
+ #if( UNIPI_VIS_ARC_COST )
+  cout << " c = " << arc->cost;
+ #endif
+ #if( QUADRATICCOST )
+  #if( UNIPI_VIS_ARC_Q_COST )
+   cout << " q = " << arc->quadraticCost;
+  #endif
+  cout << endl;
+  for( int t = 0 ; t < tab ; t++ )
+   cout << "\t";
+  #if( UNIPI_VIS_ARC_REDUCT_COST )
+   cout << " rc = " << MCFGetRC( ind );
+  #endif
+ #else
+  cout << endl;
+  for( int t = 0 ; t < tab ; t++ )
+   cout << "\t";
+  #if( UNIPI_VIS_ARC_REDUCT_COST )
+   cout << " rc = " << MCFGetRC( ind );
+  #endif
+  #if (UNIPI_VIS_ARC_STATE)
+  switch( arc->ident ) {
+  case( BASIC ):    cout << " in T"; break;
+  case( AT_LOWER ): cout << " in L"; break;
+  case( AT_UPPER ): cout << " in U"; break;
+  case( DELETED ):  cout << " canceled"; break;
+  case( CLOSED ):   cout << " closed";
+  }
+  #endif
+ #endif
+ cout << endl;
+ } 
+
+/*--------------------------------------------------------------------------*/
+
+void MCFSimplex::ShowSituation( int tab )
+{
+ if( usePrimalSimplex ) {
+  arcPType *arc;
+  nodePType *node;
+  int i = 0;
+  for( arc = arcsP ; arc != stopArcsP ; arc++ ) {
+   infoPArc( arc , i , tab );
+   i++;
+   }
+  cout << endl;
+  #if( UNIPI_VIS_DUMMY_ARCS )
+   i = 0;
+   for( arc = dummyArcsP ; arc != stopDummyP ; arc++ ) {
+    infoPArc( arc , i , tab );                
+    i++;
+    }
+   cout << endl;
+  #endif
+  infoPNode( dummyRootP , tab );
+  for( node = nodesP ; node != stopNodesP ; node++ )
+   infoPNode( node , tab );
+  }
+ else {
+  arcDType *arc;
+  nodeDType *node;
+  int i = 0;
+  for( arc = arcsD ; arc != stopArcsD ; arc++ ) {
+   infoDArc( arc , i , tab );
+   i++;
+   }
+  cout << endl;
+  #if( UNIPI_VIS_DUMMY_ARCS )
+   i = 0;
+   for( arc = dummyArcsD ; arc != stopDummyD ; arc++) {
+    infoDArc( arc , i , tab );
+    i++;
+    }
+   cout << endl;
+  #endif
+  infoDNode( dummyRootD , tab );
+  for( node = nodesD ; node != stopNodesD ; node++ )
+   infoDNode( node , tab );
+  }
+ }
+
+/*-------------------------------------------------------------------------*/
+/*---------------------- End File MCFSimplex.C ----------------------------*/
+/*-------------------------------------------------------------------------*/