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