---

dinic.cpp

#include "dinic.hpp"
#include "PG_algorithm.hpp"

Graph * g;
Graph::VertexSetHandle grayset1, grayset2, set_s_vertex, set_t_vertex;

VertexPropertyMap< vertex_level_tag, Graph > level;
VertexPropertyMap< vertex_potential_tag, Graph > potential;
EdgeMap< edgelist1_tag, Graph > network_edges;
EdgeMap< edgelist_residual_tag, Graph > residual_edges;
EdgePropertyMap< edge_flow_tag, EdgeMap< edgelist1_tag, Graph > >
       network_flow;
EdgePropertyMap< edge_flow_tag, EdgeMap< edgelist_residual_tag, Graph > >
       residual_flow;
EdgePropertyMap< edge_capacity_tag, EdgeMap< edgelist1_tag, Graph > >
       network_capacity;
EdgePropertyMap< edge_capacity_tag, EdgeMap< edgelist_residual_tag, Graph > > 
       residual_capacity;
EdgePropertyMap< edge_flag1_tag, EdgeMap< edgelist_residual_tag, Graph > >
       edgeflag;
EdgePropertyMap< edge_reference_tag, EdgeMap< edgelist_residual_tag, Graph > >
       edgereference;

size_t iteration = 0;

#define LOCALTIME
double crg_time = 0, clg_time = 0, pf_time = 0;

bool is_edge_marked( ResidualEdge e ) {
  return get( edgeflag, e );
}

struct create_residual_edge
{
  AllToAllBuffer buf;

  void operator()( NetworkEdge e ) {
    unsigned int flow = get( network_flow, e );
    unsigned int cap = get( network_capacity, e );
    Vertex v1 = source(e, *g);
    Vertex v2 = target(e, *g);

    if( flow < cap )
    {
      residual_edges_property_type prop( cap - flow, 
      residual_edges_property_type::next_type( true,
      residual_edges_property_type::next_type::next_type( get_edge_id(e),
      residual_edges_property_type::next_type::next_type::next_type( 0 ))));

      add_edge( residual_edges, v1, v2, prop, *g, buf );
    }
    if( flow > 0 )
    {
      residual_edges_property_type prop( flow,
      residual_edges_property_type::next_type( true,
      residual_edges_property_type::next_type::next_type( get_edge_id(e),
      residual_edges_property_type::next_type::next_type::next_type( 0 ))));

      add_edge( residual_edges, v2, v1, prop, *g, buf );
    }
  }
};

struct flow_tag {};
typedef property< flow_tag, unsigned int > first_pull_type;

struct pull_discover_vertex
{
  typedef on_discover_vertex Category;

  void operator() ( Vertex v )
  {
    unsigned int pot = 0;
    if( is_inside( v, set_t_vertex ) )
    {
      residual_edgelist_type::in_iterator eit, eend;
      boost::tie( eit, eend ) = in_edges( *g, v, residual_edges );
      for(; eit != eend ; eit++ ) pot += get( residual_capacity, *eit );
    }
    put( potential, v, *g, pot );
  }
};

struct pull_choose_edges
{
  typedef on_examine_vertex Category;

  void operator() ( Vertex v )
  {
    unsigned int pot = get( potential, v, *g );
    residual_edgelist_type::in_iterator eit, eend;
    boost::tie( eit, eend ) = in_edges( *g, v, residual_edges );
    // Variante 1
    for(; eit != eend && pot != 0 ; eit++ )
    {
      unsigned int m, cap = get( residual_capacity, *eit );
      m = std::min( cap, pot );
      put( residual_flow, *eit, m );
      put( edgeflag, *eit, true );
      pot -= m;
    }

    // Variante 2
/*    unsigned int cap = 0;
    residual_edgelist_type::in_iterator eit2;
    eit2 = eit;
    for(; eit != eend ; eit++ )
    {
      cap += get( residual_capacity, eit );
    }
    if( cap <= pot )
    {
      pot -= cap;
      for(; eit2 != eend ; eit2++ )
      {
        cap = get( residual_capacity, eit2 );
        put( residual_flow, eit2, cap );
        put( edgeflag, eit2, true );
      }
    }
    else
    {
      double f = ( ((double) pot) / ((double) cap) );
      for(; eit2 != eend && pot != 0 ; eit2++ )
      {
        cap = get( residual_capacity, eit2 );
        unsigned int m = (unsigned int)( ((double) cap) * f );
        if( m == 0 ) m = 1;
        put( residual_flow, eit2, m );
        put( edgeflag, eit2, true );
        pot -= m;
      }
      }*/
  }
};

struct pull_on_edge_prepare
{
  typedef on_examine_edge Category;

  void operator() ( Vertex v1, ResidualEdge & e, Vertex v2, first_pull_type & prop )
  {
    get_property_value( prop, flow_tag() ) = get( residual_flow, e );
  }
  void operator() ( Vertex v1, ResidualEdge e, Vertex v2 ) {}
};

template< class Cat >
struct pull_on_edge
{
  typedef Cat Category;

  void operator() ( VertexID_t v1, EdgeID_t eid, VertexID_t v2id, first_pull_type & prop, no_property & )
  {
    Vertex v2 = g->get_known_vertex( v2id );
    ResidualEdge e = out_edge( *g, v2, eid, residual_edges );
    unsigned int m = get_property_value( prop, flow_tag() );
    put( potential, v2, *g, get( potential, v2, *g ) + m );
    put( residual_flow, e, m );
    put( edgeflag, e, true );
  }
  void operator() ( Vertex v1, ResidualEdge e, Vertex v2 )
  {
    unsigned int m = get( residual_flow, e );
    put( potential, v2, *g, get( potential, v2, *g ) + m );
  }
};

struct base_edge_pos_tag {};
struct remove_edge_tag {};
typedef property< base_edge_pos_tag, bool,
        property< remove_edge_tag, bool,
        property< flow_tag, unsigned int > > > first_transfer_type;

struct tf_new_vertex {
  typedef on_discover_vertex Category;
  void operator() ( Vertex v1 )
  {
    if( !is_inside( v1, set_s_vertex ) )
      put( potential, v1, *g, 0 );
  }
};

struct transfer_flow
{
  typedef on_examine_edge Category;

  void operator() ( Vertex v1, ResidualEdge & e, Vertex v2, first_transfer_type & prop )
  {
    unsigned int pot = get( potential, v1, *g ),
                 flow = get( residual_flow, e ),
                 cap = get( residual_capacity, e ), m;
    m = std::min( flow, pot );
    put( potential, v1, *g, pot - m );
    get_property_value( prop, flow_tag() ) = m;

    pair< NetworkEdge, bool > base_edge =
        out_edge_test( *g, v1, get( edgereference, e), network_edges );
    if( base_edge.second )
    {
      put( network_flow, base_edge.first, get( network_flow, base_edge.first ) + m );
      get_property_value( prop, base_edge_pos_tag() ) = false;
    }
    else
    {
      get_property_value( prop, base_edge_pos_tag() ) = true;
    }

    if( m == cap || (!is_inside( v1, set_s_vertex ) &&
                     in_degree( *g, v1, residual_edges ) == 0 ) )
    {
      remove_edge( e, *g, residual_edges );
      get_property_value( prop, remove_edge_tag() ) = true;
    }
    else
    {
      put( residual_capacity, e, cap - m );
      put( edgeflag, e, false );
      get_property_value( prop, remove_edge_tag() ) = false;
    }
  }
  void operator() ( Vertex v1, ResidualEdge e, Vertex v2 ){}
};

void if_remove_residual_edge( ResidualEdge e, Vertex v2, 
                              unsigned int cap_rest, bool b )
{
  if( b )
  {
    remove_edge( e, *g, residual_edges );
    if( in_degree( *g, v2, residual_edges ) == 0 )
    {
      residual_edgelist_type::out_iterator outit, outend;
      boost::tie( outit, outend ) = out_edges( *g, v2, residual_edges );
      for(; outit != outend ; outit++ )
      {
        if( get( edgeflag, *outit ) == false )
        {
          put( edgeflag, *outit, true );
          put( residual_flow, *outit, 0 );
        }
      }
    }
  }
  else
  {
    put( residual_capacity, e, cap_rest );
    put( edgeflag, e, false );
  }
}

template< class Cat >
struct transfer_flow_other_side
{
  typedef Cat Category;

  void operator() ( VertexID_t v1, EdgeID_t eid, VertexID_t v2id, first_transfer_type & prop, no_property & )
  {
    Vertex v2 = g->get_known_vertex( v2id );
    ResidualEdge e = in_edge( *g, v2, eid, residual_edges );
    unsigned int m = get_property_value( prop, flow_tag() ),
                 cap = get( residual_capacity, e ),
                 pot = get( potential, v2, *g );
    put( potential, v2, *g, pot + m );

    if( get_property_value( prop, base_edge_pos_tag() ) == true )
    {
      NetworkEdge base_edge = out_edge( *g, v2, get( edgereference, e ), network_edges );
      put( network_flow, base_edge, get( network_flow, base_edge ) - m );
    }

    if_remove_residual_edge( e, v2, cap - m, get_property_value( prop, remove_edge_tag() ) );
  }
  void operator() ( Vertex v1, ResidualEdge e, Vertex v2 )
  {
    unsigned int pot1 = get( potential, v1, *g ), m, 
                 pot2 = get( potential, v2, *g ),
                 flow = get( residual_flow, e ),
                 cap = get( residual_capacity, e );
    m = std::min( flow, pot1 );
    put( potential, v1, *g, pot1 - m );
    put( potential, v2, *g, pot2 + m );

    pair< NetworkEdge, bool > base_edge_test =
        out_edge_test( *g, v1, get( edgereference, e), network_edges );
    if( base_edge_test.second )
    {
      put( network_flow, base_edge_test.first, get( network_flow, base_edge_test.first ) + m );
    }
    else
    {
      NetworkEdge base_edge = out_edge( *g, v2, get( edgereference, e ), network_edges );
      put( network_flow, base_edge, get( network_flow, base_edge ) - m );       
    }

    if_remove_residual_edge( e, v2, cap - m, m == cap || (!is_inside( v1, set_s_vertex ) &&
           in_degree( *g, v1, residual_edges ) == 0 ) );
  }
};

void create_residual_graph()
{
    if( _libbase.rank() == 0 ) std::cout<<"Erzeuge Residualgraph..."<<std::endl;
#ifdef LOCALTIME
  double local_time1, local_time2;
  local_time1 = MPI_Wtime();
#endif
  create_residual_edge cre;
  for_each_out_edge( *g, network_edges, cre );
  add_delayed_edges( residual_edges, *g, cre.buf );
#ifdef LOCALTIME
  local_time2 = MPI_Wtime();
  crg_time += (local_time2 - local_time1);
#endif
}

bool create_level_graph()
{
    if( _libbase.rank() == 0 ) std::cout<<"Erzeuge Levelgraph..."<<std::endl;
#ifdef LOCALTIME
  double local_time1, local_time2;
  local_time1 = MPI_Wtime();
#endif
  bool found = false;
  copy( set_s_vertex, grayset1 );

  bfs_is_vertex_inside_targetset< Graph::VertexSet >
      vis1( *set_t_vertex, &found );
  bfs_set_edgeproperty_both_sides< on_target_new_vertex, Graph,
      EdgeMap< edgelist_residual_tag, Graph >,
      EdgePropertyMap< edge_flag1_tag, EdgeMap< edgelist_residual_tag, Graph > > > vis2( *g, false );
  bfs_set_edgeproperty_both_sides< on_target_next_level, Graph,
      EdgeMap< edgelist_residual_tag, Graph >,
      EdgePropertyMap< edge_flag1_tag, EdgeMap< edgelist_residual_tag, Graph > > > vis3( *g, false );

  breadth_first_search( *g, *grayset1, *grayset2, &found, level, residual_edges,
                        make_bfs_visitor( bfs_distribution_types< no_property, no_property >(),
                           make_pair( vis3.get_functors(), make_pair( vis2.get_functors(), vis1 ))) );
 
  remove_edges_if( *g, residual_edges, is_edge_marked );
#ifdef LOCALTIME
  local_time2 = MPI_Wtime();
  clg_time += (local_time2 - local_time1);
#endif

  return PG_ReduceAnd( !found );
}


bool pull_flow()
{
    if( _libbase.rank() == 0 ) std::cout<<"Ermittle Fluss im Levelgraph..."<<std::endl;
#ifdef LOCALTIME
  double local_time1, local_time2;
  local_time1 = MPI_Wtime();
#endif
  bool found;
  FilteredEdgeMap< bool (*)(residual_edgelist_type::Edge), 
      ReverseEdgeMap< edgelist_residual_tag, Graph > > emap(is_edge_marked);
  Graph::VertexSet::set_type gset1, gset2;
  copy( set_t_vertex, gset1 );

  pull_discover_vertex vis1;
  pull_choose_edges vis2;
  pull_on_edge_prepare vis3;
  pull_on_edge< on_target_new_vertex > vis4;
  pull_on_edge< on_target_next_level > vis5;

  breadth_first_search( *g, gset1, gset2, &found, level, emap,
                       make_bfs_visitor( bfs_distribution_types< first_pull_type, no_property >(),
                                         make_pair( vis5, make_pair( vis4, make_pair( vis3, make_pair( vis1, vis2)))) ) );

  copy( set_s_vertex, gset1 );
  FilteredEdgeMap< bool (*)(residual_edgelist_type::Edge),
      EdgeMap< edgelist_residual_tag, Graph > > emap2(is_edge_marked);
  transfer_flow tvis1;
  tf_new_vertex tvis2;
  transfer_flow_other_side< on_target_new_vertex > tvis3;
  transfer_flow_other_side< on_target_next_level > tvis4;

  breadth_first_search( *g, gset1, gset2, &found, level, emap2,
                       make_bfs_visitor( bfs_distribution_types< first_transfer_type, no_property >(),
                          make_pair( tvis1, make_pair( tvis2, make_pair( tvis3, tvis4 ))) ) );
#ifdef LOCALTIME
  local_time2 = MPI_Wtime();
  pf_time += (local_time2 - local_time1);
#endif

  bool has_blocking_flow = false;
  if( ! (*set_t_vertex).empty() )
  {
    if( in_degree( *g, *((*set_t_vertex).begin()), residual_edges ) == 0 )
      has_blocking_flow = true;
  }
  return PG_ReduceOr( has_blocking_flow );
}

void dinic()
{
  double  starttime, endtime;
  starttime = MPI_Wtime();
  copy( set_s_vertex, grayset1 );
  bfs_do_distribution( *g, *grayset1, *grayset2, level, network_edges );

  for(;;)
  {
    create_residual_graph();
    if( create_level_graph() ) break;
    while( !pull_flow() );

    remove_edges( *g, residual_edges );
    iteration++;
  }
  endtime = MPI_Wtime();

  if( ! (*set_s_vertex).empty() )
  {
    unsigned int flow = 0;
    network_edgelist_type::out_iterator eit, eend;
    boost::tie( eit, eend ) = out_edges( *g, *(*set_s_vertex).begin(), network_edges );
    for(; eit != eend ; eit++ )
    {
      flow += get( network_flow, *eit );
    }
    std::cout << "\nAusgehender Fluss von s: " << flow << std::endl;
    std::cout << "Laufzeit: " << (endtime-starttime) << "s" << std::endl;
  }
#ifdef LOCALTIME
  for(size_t i = 0 ; i<_libbase.size() ; i++ )
  {
    MPI_Barrier( _libbase.comm() );
    if( i == _libbase.rank() )
    {
      std::cout << "Residualgraph: " << crg_time << "s" << std::endl;
      std::cout << "Levelgraph: " << clg_time << "s" << std::endl;
      std::cout << "Pull Flow: " << pf_time << "s" << std::endl;
   }
  }
#endif
}

/*
   Hauptroutine
*/
int main( int argc, char *argv[] )
{
  ParGraph_init( argc, argv );
  g = new Graph();
  grayset1 = add_vertex_set( *g );
  grayset2 = add_vertex_set( *g );
  set_s_vertex = add_vertex_set( *g );
  set_t_vertex = add_vertex_set( *g );

  dinic_read_graph( argc, argv );
  dinic();

  delete g;
  ParGraph_finalize();
  return 0;
}

---