diff -r f196dc4f1b31 -r 23a54f889272 src/work/marci/lp/max_flow_expression.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/marci/lp/max_flow_expression.cc	Fri Jan 28 14:33:32 2005 +0000
@@ -0,0 +1,87 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <lemon/smart_graph.h>
+#include <lemon/list_graph.h>
+#include <lemon/dimacs.h>
+#include <lemon/time_measure.h>
+#include <lp_solver_wrapper_3.h>
+
+using std::cout;
+using std::endl;
+using namespace lemon;
+
+template<typename Edge, typename EdgeIndexMap> 
+class PrimalMap {
+protected:
+  LPGLPK* lp;
+  EdgeIndexMap* edge_index_map;
+public:
+  PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : 
+    lp(&_lp), edge_index_map(&_edge_index_map) { }
+  double operator[](Edge e) const { 
+    return lp->getPrimal((*edge_index_map)[e]);
+  }
+};
+
+// Use a DIMACS max flow file as stdin.
+// max_flow_expresion < dimacs_max_flow_file
+
+int main(int, char **) {
+
+  typedef ListGraph Graph;
+  typedef Graph::Node Node;
+  typedef Graph::Edge Edge;
+  typedef Graph::EdgeIt EdgeIt;
+
+  Graph g;
+  
+  Node s, t;
+  Graph::EdgeMap<int> cap(g);
+  readDimacs(std::cin, g, cap, s, t);
+  Timer ts;
+  
+  typedef LPGLPK LPSolver;
+  LPSolver lp;
+  lp.setMaximize();
+  typedef LPSolver::ColIt ColIt;
+  typedef LPSolver::RowIt RowIt;
+  typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
+  EdgeIndexMap edge_index_map(g);
+  PrimalMap<Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
+
+  // capacity function
+  for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
+    ColIt col_it=lp.addCol();
+    edge_index_map.set(e, col_it);
+    if (cap[e]==0)
+      lp.setColBounds(col_it, LPSolver::FIXED, 0, cap[e]);
+    else 
+      lp.setColBounds(col_it, LPSolver::DOUBLE, 0, cap[e]);
+  }
+  
+  for (Graph::NodeIt n(g); n!=INVALID; ++n) {
+    LPSolver::Expression expr;
+    for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
+      expr+=edge_index_map[e];
+    for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
+      expr-=edge_index_map[e];
+    // cost function
+    if (n==s) {
+      lp.setObjCoeffs(expr);      
+    }
+    // flow conservation
+    if ((n!=s) && (n!=t)) {
+      RowIt row_it=lp.addRow();
+      lp.setRowCoeffs(row_it, expr);
+      lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
+    }
+  }
+  lp.solveSimplex();
+  //std::cout << lp.colNum() << std::endl;
+  //std::cout << lp.rowNum() << std::endl;
+  //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
+  for (Graph::EdgeIt e(g); e!=INVALID; ++e) 
+    flow.set(e, lp_flow[e]);
+}