diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/max_flow_expression.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/max_flow_expression.cc	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,109 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <lemon/graph_utils.h>
+#include <lemon/smart_graph.h>
+#include <lemon/list_graph.h>
+#include <lemon/dimacs.h>
+#include <lemon/time_measure.h>
+#include <lp_solver_glpk.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::Col Col;
+  typedef LPSolver::Row Row;
+  typedef Graph::EdgeMap<Col> EdgeIndexMap;
+  typedef Graph::NodeMap<Row> NodeIndexMap;
+  EdgeIndexMap edge_index_map(g);
+  NodeIndexMap node_index_map(g);
+  PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);
+
+  // nonnegativity of flow and capacity function
+  for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
+    Col col=lp.addCol();
+    edge_index_map.set(e, col);
+    // interesting property in GLPK:
+    // if you change the order of the following two lines, the 
+    // two runs of GLPK are extremely different
+      lp.setColLowerBound(col, 0);
+      lp.setColUpperBound(col, 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 constraints
+    if ((n!=s) && (n!=t)) {
+      node_index_map.set(n, lp.addRow(expr == 0.0));
+    }
+  }
+  lp.solveSimplex();
+  cout << "elapsed time: " << ts << endl;
+//   cout << "rows:" << endl;
+//   for (
+//        LPSolver::Rows::ClassIt i(lp.row_iter_map, 0);
+//        i!=INVALID;
+//        ++i) { 
+//     cout << i << " ";
+//   }
+//   cout << endl;
+//   cout << "cols:" << endl;
+//   for (
+//        LPSolver::Cols::ClassIt i(lp.col_iter_map, 0);
+//        i!=INVALID;
+//        ++i) { 
+//     cout << i << " ";
+//   }
+//   cout << endl;
+  lp.setMIP();
+  cout << "elapsed time: " << ts << endl;
+  for (LPSolver::Cols::ClassIt it(lp.col_iter_map ,1); it!=INVALID; ++it) {
+    lp.setColInt(it);
+  }
+  cout << "elapsed time: " << ts << endl;
+  lp.solveBandB();
+  cout << "elapsed time: " << ts << endl;
+}