marci@1104: // -*- c++ -*-
marci@1104: #include <iostream>
marci@1104: #include <fstream>
marci@1104: 
marci@1143: #include <lemon/graph_utils.h>
marci@1104: #include <lemon/smart_graph.h>
marci@1104: #include <lemon/list_graph.h>
marci@1104: #include <lemon/dimacs.h>
marci@1104: #include <lemon/time_measure.h>
marci@1111: #include <lp_solver_base.h>
marci@1104: 
marci@1104: using std::cout;
marci@1104: using std::endl;
marci@1104: using namespace lemon;
marci@1104: 
marci@1104: template<typename Edge, typename EdgeIndexMap> 
marci@1104: class PrimalMap {
marci@1104: protected:
marci@1104:   LPGLPK* lp;
marci@1104:   EdgeIndexMap* edge_index_map;
marci@1104: public:
marci@1104:   PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : 
marci@1104:     lp(&_lp), edge_index_map(&_edge_index_map) { }
marci@1104:   double operator[](Edge e) const { 
marci@1104:     return lp->getPrimal((*edge_index_map)[e]);
marci@1104:   }
marci@1104: };
marci@1104: 
marci@1104: // Use a DIMACS max flow file as stdin.
marci@1104: // max_flow_expresion < dimacs_max_flow_file
marci@1104: 
marci@1104: int main(int, char **) {
marci@1104: 
marci@1104:   typedef ListGraph Graph;
marci@1104:   typedef Graph::Node Node;
marci@1104:   typedef Graph::Edge Edge;
marci@1104:   typedef Graph::EdgeIt EdgeIt;
marci@1104: 
marci@1104:   Graph g;
marci@1104:   
marci@1104:   Node s, t;
marci@1104:   Graph::EdgeMap<int> cap(g);
marci@1104:   readDimacs(std::cin, g, cap, s, t);
marci@1104:   Timer ts;
marci@1104:   
marci@1104:   typedef LPGLPK LPSolver;
marci@1104:   LPSolver lp;
marci@1104:   lp.setMaximize();
marci@1144:   typedef LPSolver::Col Col;
marci@1144:   typedef LPSolver::Row Row;
marci@1144:   typedef Graph::EdgeMap<Col> EdgeIndexMap;
marci@1144:   typedef Graph::NodeMap<Row> NodeIndexMap;
marci@1104:   EdgeIndexMap edge_index_map(g);
marci@1143:   NodeIndexMap node_index_map(g);
marci@1110:   PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);
marci@1104: 
marci@1111:   // nonnegativity of flow and capacity function
marci@1104:   for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
marci@1144:     Col col=lp.addCol();
marci@1144:     edge_index_map.set(e, col);
marci@1110:     // interesting property in GLPK:
marci@1110:     // if you change the order of the following two lines, the 
marci@1110:     // two runs of GLPK are extremely different
marci@1144:       lp.setColLowerBound(col, 0);
marci@1144:       lp.setColUpperBound(col, cap[e]);
marci@1104:   }
marci@1104:   
marci@1104:   for (Graph::NodeIt n(g); n!=INVALID; ++n) {
marci@1104:     LPSolver::Expression expr;
marci@1104:     for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
marci@1104:       expr+=edge_index_map[e];
marci@1104:     for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
marci@1104:       expr-=edge_index_map[e];
marci@1104:     // cost function
marci@1104:     if (n==s) {
marci@1104:       lp.setObjCoeffs(expr);      
marci@1104:     }
marci@1111:     // flow conservation constraints
marci@1104:     if ((n!=s) && (n!=t)) {
marci@1144:       node_index_map.set(n, lp.addRow(expr == 0.0));
marci@1104:     }
marci@1104:   }
marci@1104:   lp.solveSimplex();
marci@1110:   cout << "elapsed time: " << ts << endl;
marci@1104: }