Index: lemon/lp_base.h
===================================================================
--- lemon/lp_base.h	(revision 2363)
+++ lemon/lp_base.h	(revision 2364)
@@ -50,5 +50,5 @@
     
   public:
-
+    
     ///Possible outcomes of an LP solving procedure
     enum SolveExitStatus {
@@ -123,4 +123,5 @@
       friend class LpSolverBase;
       friend class MipSolverBase;
+      explicit Col(int _id) : id(_id) {}
     public:
       typedef Value ExprValue;
@@ -166,4 +167,5 @@
       int id;
       friend class LpSolverBase;
+      explicit Row(int _id) : id(_id) {}
     public:
       typedef Value ExprValue;
@@ -178,4 +180,20 @@
     };
 
+    class RowIt : public Row {
+      LpSolverBase *_lp;
+    public:
+      RowIt() {}
+      RowIt(LpSolverBase &lp) : _lp(&lp)
+      {
+        _lp->rows.firstFix(id);
+      }
+      RowIt(const Invalid&) : Row(INVALID) {}
+      RowIt &operator++() 
+      {
+        _lp->rows.nextFix(id);
+	return *this;
+      }
+    };
+
     static int id(const Row& row) { return row.id; }
 
@@ -188,4 +206,12 @@
     int _lpId(const Row& row) const {
       return rows.floatingId(id(row));
+    }
+
+    Col _item(int id, Col) const {
+      return Col(cols.fixId(id));
+    }
+
+    Row _item(int id, Row) const {
+      return Row(rows.fixId(id));
     }
 
@@ -582,9 +608,58 @@
   private:
 
-    template <typename _Base>
-    class MappedIterator {
+    template <typename _Expr>
+    class MappedOutputIterator {
     public:
 
-      typedef _Base Base;
+      typedef std::insert_iterator<_Expr> Base;
+
+      typedef std::output_iterator_tag iterator_category;
+      typedef void difference_type;
+      typedef void value_type;
+      typedef void reference;
+      typedef void pointer;
+      
+      MappedOutputIterator(const Base& _base, const LpSolverBase& _lp) 
+        : base(_base), lp(_lp) {}
+
+      MappedOutputIterator& operator*() {
+        return *this;
+      }
+
+      MappedOutputIterator& operator=(const std::pair<int, Value>& value) {
+        *base = std::make_pair(lp._item(value.first, typename _Expr::Key()), 
+                               value.second);
+        return *this;
+      }
+
+      MappedOutputIterator& operator++() {
+        ++base;
+        return *this;
+      }
+
+      MappedOutputIterator operator++(int) {
+        MappedOutputIterator tmp(*this);
+        ++base;
+        return tmp;
+      }
+
+      bool operator==(const MappedOutputIterator& it) const {
+        return base == it.base;
+      }
+
+      bool operator!=(const MappedOutputIterator& it) const {
+        return base != it.base;
+      }
+
+    private:
+      Base base;
+      const LpSolverBase& lp;
+    };
+
+    template <typename Expr>
+    class MappedInputIterator {
+    public:
+
+      typedef typename Expr::const_iterator Base;
 
       typedef typename Base::iterator_category iterator_category;
@@ -600,5 +675,5 @@
       };
 
-      MappedIterator(const Base& _base, const LpSolverBase& _lp) 
+      MappedInputIterator(const Base& _base, const LpSolverBase& _lp) 
         : base(_base), lp(_lp) {}
 
@@ -611,20 +686,20 @@
       }
 
-      MappedIterator& operator++() {
+      MappedInputIterator& operator++() {
         ++base;
         return *this;
       }
 
-      MappedIterator& operator++(int) {
-        MappedIterator tmp(*this);
+      MappedInputIterator operator++(int) {
+        MappedInputIterator tmp(*this);
         ++base;
         return tmp;
       }
 
-      bool operator==(const MappedIterator& it) const {
+      bool operator==(const MappedInputIterator& it) const {
         return base == it.base;
       }
 
-      bool operator!=(const MappedIterator& it) const {
+      bool operator!=(const MappedInputIterator& it) const {
         return base != it.base;
       }
@@ -638,7 +713,12 @@
 
     /// STL compatible iterator for lp col
-    typedef MappedIterator<Expr::const_iterator> LpRowIterator;
+    typedef MappedInputIterator<Expr> ConstRowIterator;
     /// STL compatible iterator for lp row
-    typedef MappedIterator<DualExpr::const_iterator> LpColIterator;
+    typedef MappedInputIterator<DualExpr> ConstColIterator;
+
+    /// STL compatible iterator for lp col
+    typedef MappedOutputIterator<Expr> RowIterator;
+    /// STL compatible iterator for lp row
+    typedef MappedOutputIterator<DualExpr> ColIterator;
 
     //Abstract virtual functions
@@ -660,9 +740,12 @@
     virtual void _getColName(int col, std::string & name) = 0;
     virtual void _setColName(int col, const std::string & name) = 0;
-    virtual void _setRowCoeffs(int i, LpRowIterator b, LpRowIterator e) = 0;
-    virtual void _setColCoeffs(int i, LpColIterator b, LpColIterator e) = 0;
+    virtual void _setRowCoeffs(int i, ConstRowIterator b, 
+                               ConstRowIterator e) = 0;
+    virtual void _getRowCoeffs(int i, RowIterator b) = 0;
+    virtual void _setColCoeffs(int i, ConstColIterator b, 
+                               ConstColIterator e) = 0;
+    virtual void _getColCoeffs(int i, ColIterator b) = 0;
     virtual void _setCoeff(int row, int col, Value value) = 0;
     virtual Value _getCoeff(int row, int col) = 0;
-
     virtual void _setColLowerBound(int i, Value value) = 0;
     virtual Value _getColLowerBound(int i) = 0;
@@ -787,6 +870,16 @@
     void col(Col c,const DualExpr &e) {
       e.simplify();
-      _setColCoeffs(_lpId(c), LpColIterator(e.begin(), *this), 
-                    LpColIterator(e.end(), *this));
+      _setColCoeffs(_lpId(c), ConstColIterator(e.begin(), *this), 
+                    ConstColIterator(e.end(), *this));
+    }
+
+    ///Get a column (i.e a dual constraint) of the LP
+
+    ///\param r is the column to get
+    ///\return the dual expression associated to the column
+    DualExpr col(Col c) {
+      DualExpr e;
+      _getColCoeffs(_lpId(c), ColIterator(std::inserter(e, e.end()), *this));
+      return e;
     }
 
@@ -884,7 +977,7 @@
     void row(Row r, Value l,const Expr &e, Value u) {
       e.simplify();
-      _setRowCoeffs(_lpId(r), LpRowIterator(e.begin(), *this),
-                    LpRowIterator(e.end(), *this));
-       _setRowBounds(_lpId(r),l-e.constComp(),u-e.constComp());
+      _setRowCoeffs(_lpId(r), ConstRowIterator(e.begin(), *this),
+                    ConstRowIterator(e.end(), *this));
+      _setRowBounds(_lpId(r),l-e.constComp(),u-e.constComp());
     }
 
@@ -896,4 +989,15 @@
       row(r, c.lowerBounded()?c.lowerBound():-INF,
           c.expr(), c.upperBounded()?c.upperBound():INF);
+    }
+
+    
+    ///Get a row (i.e a constraint) of the LP
+
+    ///\param r is the row to get
+    ///\return the expression associated to the row
+    Expr row(Row r) {
+      Expr e;
+      _getRowCoeffs(_lpId(r), RowIterator(std::inserter(e, e.end()), *this));
+      return e;
     }
 
@@ -1178,4 +1282,19 @@
     }
 
+    ///Get the objective function
+
+    ///\return the objective function as a linear expression of type \ref Expr.
+    Expr obj() {
+      Expr e;
+      for (ColIt it(*this); it != INVALID; ++it) {
+        double c = objCoeff(it);
+        if (c != 0.0) {
+          e.insert(std::make_pair(it, c));
+        }
+      }
+      return e;
+    }
+    
+
     ///Maximize
     void max() { _setMax(); }
Index: lemon/lp_cplex.cc
===================================================================
--- lemon/lp_cplex.cc	(revision 2363)
+++ lemon/lp_cplex.cc	(revision 2364)
@@ -110,5 +110,5 @@
   
   ///\warning Data at index 0 is ignored in the arrays.
-  void LpCplex::_setRowCoeffs(int i, LpRowIterator b, LpRowIterator e)
+  void LpCplex::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e)
   {
     std::vector<int> indices;
@@ -116,5 +116,5 @@
     std::vector<Value> values;
 
-    for(LpRowIterator it=b; it!=e; ++it) {
+    for(ConstRowIterator it=b; it!=e; ++it) {
       indices.push_back(it->first);
       values.push_back(it->second);
@@ -125,6 +125,10 @@
 			    &rowlist[0], &indices[0], &values[0]); 
   }
-  
-  void LpCplex::_setColCoeffs(int i, LpColIterator b, LpColIterator e)
+
+  void LpSoplex::_getRowCoeffs(int i, RowIterator b) {
+    /// \todo implement
+  }
+  
+  void LpCplex::_setColCoeffs(int i, ConstColIterator b, ConstColIterator e)
   {
     std::vector<int> indices;
@@ -132,5 +136,5 @@
     std::vector<Value> values;
 
-    for(LpColIterator it=b; it!=e; ++it) {
+    for(ConstColIterator it=b; it!=e; ++it) {
       indices.push_back(it->first);
       values.push_back(it->second);
@@ -140,4 +144,8 @@
     status = CPXchgcoeflist(env, lp, values.size(), 
 			    &indices[0], &collist[0], &values[0]); 
+  }
+
+  void LpSoplex::_getColCoeffs(int i, ColIterator b) {
+    /// \todo implement
   }
   
Index: lemon/lp_cplex.h
===================================================================
--- lemon/lp_cplex.h	(revision 2361)
+++ lemon/lp_cplex.h	(revision 2364)
@@ -63,6 +63,8 @@
     virtual void _getColName(int col,       std::string & name);
     virtual void _setColName(int col, const std::string & name);
-    virtual void _setRowCoeffs(int i, LpRowIterator b, LpRowIterator e);
-    virtual void _setColCoeffs(int i, LpColIterator b, LpColIterator e);
+    virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e);
+    virtual void _getRowCoeffs(int i, RowIterator b);
+    virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e);
+    virtual void _getColCoeffs(int i, ColIterator b);
     virtual void _setCoeff(int row, int col, Value value);
     virtual Value _getCoeff(int row, int col);
Index: lemon/lp_glpk.cc
===================================================================
--- lemon/lp_glpk.cc	(revision 2363)
+++ lemon/lp_glpk.cc	(revision 2364)
@@ -127,5 +127,5 @@
   }
   
-  void LpGlpk::_setRowCoeffs(int i, LpRowIterator b, LpRowIterator e) 
+  void LpGlpk::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e) 
   {
     std::vector<int> indices;
@@ -135,5 +135,5 @@
     values.push_back(0);
 
-    for(LpRowIterator it=b; it!=e; ++it) {
+    for(ConstRowIterator it=b; it!=e; ++it) {
       indices.push_back(it->first);
       values.push_back(it->second);
@@ -142,6 +142,21 @@
     lpx_set_mat_row(lp, i, values.size() - 1, &indices[0], &values[0]);
   }
-  
-  void LpGlpk::_setColCoeffs(int i, LpColIterator b, LpColIterator e) {
+
+  void LpGlpk::_getRowCoeffs(int i, RowIterator b) 
+  {
+    int length = lpx_get_mat_row(lp, i, 0, 0);
+    
+    std::vector<int> indices(length + 1);
+    std::vector<Value> values(length + 1);
+    
+    lpx_get_mat_row(lp, i, &indices[0], &values[0]);
+    
+    for (int i = 1; i <= length; ++i) {
+      *b = std::make_pair(indices[i], values[i]);
+      ++b;
+    }
+  }
+  
+  void LpGlpk::_setColCoeffs(int i, ConstColIterator b, ConstColIterator e) {
 
     std::vector<int> indices;
@@ -151,5 +166,5 @@
     values.push_back(0);
 
-    for(LpColIterator it=b; it!=e; ++it) {
+    for(ConstColIterator it=b; it!=e; ++it) {
       indices.push_back(it->first);
       values.push_back(it->second);
@@ -159,4 +174,18 @@
   }
 
+  void LpGlpk::_getColCoeffs(int i, ColIterator b) 
+  {
+    int length = lpx_get_mat_col(lp, i, 0, 0);
+    
+    std::vector<int> indices(length + 1);
+    std::vector<Value> values(length + 1);
+    
+    lpx_get_mat_col(lp, i, &indices[0], &values[0]);
+    
+    for (int i = 1; i <= length; ++i) {
+      *b = std::make_pair(indices[i], values[i]);
+      ++b;
+    }
+  }
 
   void LpGlpk::_setCoeff(int row, int col, Value value) 
@@ -224,6 +253,6 @@
     int length=lpx_get_mat_row(lp, row, 0, 0);
     
-    std::vector<int> indices(length + 2);
-    std::vector<Value> values(length + 2);
+    std::vector<int> indices(length + 1);
+    std::vector<Value> values(length + 1);
     
     lpx_get_mat_row(lp, row, &indices[0], &values[0]);
Index: lemon/lp_glpk.h
===================================================================
--- lemon/lp_glpk.h	(revision 2328)
+++ lemon/lp_glpk.h	(revision 2364)
@@ -58,6 +58,8 @@
     virtual void _getColName(int col,       std::string & name);
     virtual void _setColName(int col, const std::string & name);
-    virtual void _setRowCoeffs(int i, LpRowIterator b, LpRowIterator e);
-    virtual void _setColCoeffs(int i, LpColIterator b, LpColIterator e);
+    virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e);
+    virtual void _getRowCoeffs(int i, RowIterator b);
+    virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e);
+    virtual void _getColCoeffs(int i, ColIterator b);
     virtual void _setCoeff(int row, int col, Value value);
     virtual Value _getCoeff(int row, int col);
Index: lemon/lp_skeleton.cc
===================================================================
--- lemon/lp_skeleton.cc	(revision 2328)
+++ lemon/lp_skeleton.cc	(revision 2364)
@@ -59,8 +59,14 @@
   
   
-  void LpSkeleton::_setRowCoeffs(int, LpRowIterator, LpRowIterator) {
+  void LpSkeleton::_setRowCoeffs(int, ConstRowIterator, ConstRowIterator) {
+  }
+
+  void LpSkeleton::_getRowCoeffs(int, RowIterator) {
   }
   
-  void LpSkeleton::_setColCoeffs(int, LpColIterator, LpColIterator) {
+  void LpSkeleton::_setColCoeffs(int, ConstColIterator, ConstColIterator) {
+  }
+
+  void LpSkeleton::_getColCoeffs(int, ColIterator) {
   }
 
Index: lemon/lp_skeleton.h
===================================================================
--- lemon/lp_skeleton.h	(revision 2363)
+++ lemon/lp_skeleton.h	(revision 2364)
@@ -50,7 +50,11 @@
 
     /// \e
-    virtual void _setRowCoeffs(int i, LpRowIterator b, LpRowIterator e);
+    virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e);
     /// \e
-    virtual void _setColCoeffs(int i, LpColIterator b, LpColIterator e);
+    virtual void _getRowCoeffs(int i, RowIterator b);
+    /// \e
+    virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e);
+    /// \e
+    virtual void _getColCoeffs(int i, ColIterator b);
     
     /// Set one element of the coefficient matrix
Index: lemon/lp_soplex.cc
===================================================================
--- lemon/lp_soplex.cc	(revision 2363)
+++ lemon/lp_soplex.cc	(revision 2364)
@@ -100,22 +100,38 @@
   
 
-  void LpSoplex::_setRowCoeffs(int i, LpRowIterator b, LpRowIterator e) {
+  void LpSoplex::_setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e) {
     for (int j = 0; j < soplex->nCols(); ++j) {
       soplex->changeElement(i, j, 0.0);
     }
-    for(LpRowIterator it = b; it != e; ++it) {
+    for(ConstRowIterator it = b; it != e; ++it) {
       soplex->changeElement(i, it->first, it->second);
     }
     solved = false;
   }
-  
-  void LpSoplex::_setColCoeffs(int j, LpColIterator b, LpColIterator e) {
+
+  void LpSoplex::_getRowCoeffs(int i, RowIterator b) {
+    const soplex::SVector& vec = soplex->rowVector(i);
+    for (int k = 0; k < vec.size(); ++k) {
+      *b = std::make_pair(vec.index(k), vec.value(k));
+      ++b;
+    }
+  }
+  
+  void LpSoplex::_setColCoeffs(int j, ConstColIterator b, ConstColIterator e) {
     for (int i = 0; i < soplex->nRows(); ++i) {
       soplex->changeElement(i, j, 0.0);
     }
-    for(LpColIterator it = b; it != e; ++it) {
+    for(ConstColIterator it = b; it != e; ++it) {
       soplex->changeElement(it->first, j, it->second);
     }
     solved = false;
+  }
+
+  void LpSoplex::_getColCoeffs(int i, ColIterator b) {
+    const soplex::SVector& vec = soplex->colVector(i);
+    for (int k = 0; k < vec.size(); ++k) {
+      *b = std::make_pair(vec.index(k), vec.value(k));
+      ++b;
+    }
   }
   
Index: lemon/lp_soplex.h
===================================================================
--- lemon/lp_soplex.h	(revision 2363)
+++ lemon/lp_soplex.h	(revision 2364)
@@ -73,6 +73,8 @@
     virtual void _getColName(int col, std::string & name);
     virtual void _setColName(int col, const std::string & name);
-    virtual void _setRowCoeffs(int i, LpRowIterator b, LpRowIterator e);
-    virtual void _setColCoeffs(int i, LpColIterator b, LpColIterator e);
+    virtual void _setRowCoeffs(int i, ConstRowIterator b, ConstRowIterator e);
+    virtual void _getRowCoeffs(int i, RowIterator b);
+    virtual void _setColCoeffs(int i, ConstColIterator b, ConstColIterator e);
+    virtual void _getColCoeffs(int i, ColIterator b);
     virtual void _setCoeff(int row, int col, Value value);
     virtual Value _getCoeff(int row, int col);
Index: lemon/lp_utils.h
===================================================================
--- lemon/lp_utils.h	(revision 2363)
+++ lemon/lp_utils.h	(revision 2364)
@@ -69,9 +69,14 @@
   private:
 
+    enum SubSection {
+      CONSTRAINTS, BOUNDS, OBJECTIVE
+    };
+
     enum Relation {
       LE, EQ, GE
     };
 
-    std::istream& readConstraint(std::istream& is, LpSolverBase::Constr& c) {
+    std::istream& readConstraint(std::istream& is) {
+      LpSolverBase::Constr c;
       char x;
       LpSolverBase::Expr e1, e2;
@@ -83,11 +88,25 @@
       is >> std::ws;
       readExpression(is, e2);
+      is >> std::ws;
       if (!is.get(x)) {
         if (op1 == LE) {
+          if (e1.size() == 0) {
+            c = e2 >= e1;
+          } else {
+            c = e1 <= e2;
+          }
           c = e1 <= e2;
         } else if (op1 == GE) {
-          c = e1 >= e2;
-        } else {
-          c = e1 == e2;
+          if (e1.size() == 0) {
+            c = e2 <= e1;
+          } else {
+            c = e1 >= e2;
+          }
+        } else {
+          if (e1.size() == 0) {
+            c = e2 == e1;
+          } else {
+            c = e1 == e2;
+          }
         }
       } else {
@@ -109,5 +128,102 @@
           c = e1.constComp() >= e2 >= e3.constComp();
         }
-      }
+        is >> std::ws;
+        if (is.get(x)) {
+          throw DataFormatError("Wrong variable bounds format");
+        }
+      }
+      lp.addRow(c);
+      return is;
+    }
+
+    std::istream& readObjective(std::istream& is) {
+      is >> std::ws;
+      std::string sense;
+      is >> sense;
+      if (sense != "min" && sense != "max") {
+        throw DataFormatError("Wrong objective function format");
+      }
+      LpSolverBase::Expr expr;
+      is >> std::ws;
+      readExpression(is, expr);
+      lp.setObj(expr);
+      if (sense == "min") {
+        lp.min();
+      } else {
+        lp.max();
+      }
+      return is;
+    }
+
+    std::istream& readBounds(std::istream& is) {
+      LpSolverBase::Col c;
+      char x;
+      is >> std::ws;
+      if (!is.get(x)) {
+        throw DataFormatError("Wrong variable bounds format");
+      }
+      if (varFirstChar(x)) {
+        is.putback(x);
+        readCol(is, c);
+        is >> std::ws;
+        Relation op1;
+        readRelation(is, op1);
+        double v;
+        readNum(is, v);
+        if (op1 == EQ) {
+          lp.colLowerBound(c, v);
+          lp.colUpperBound(c, v);
+        } else if (op1 == LE) {
+          lp.colUpperBound(c, v);
+        } else {
+          lp.colLowerBound(c, v);
+        }
+        is >> std::ws;
+        if (is.get(x)) {
+          throw DataFormatError("Wrong variable bounds format");
+        }
+      } else {
+        is.putback(x);
+        double v;
+        readNum(is, v);
+        is >> std::ws;
+        Relation op1;
+        readRelation(is, op1);
+        is >> std::ws;
+        readCol(is, c);
+        is >> std::ws;
+        if (is.get(x)) {
+          is.putback(x);
+          is >> std::ws;
+          Relation op2;
+          readRelation(is, op2);
+          double w;
+          is >> std::ws;
+          readNum(is, w);
+          if (op2 != op1 || op1 == EQ) {
+            throw DataFormatError("Wrong range format");
+          }
+          if (op1 == LE) {
+            lp.colLowerBound(c, v);
+            lp.colUpperBound(c, w);
+          } else {
+            lp.colLowerBound(c, w);
+            lp.colUpperBound(c, v);
+          }
+          if (is.get(x)) {
+            throw DataFormatError("Wrong variable bounds format");
+          }
+        } else {
+          if (op1 == EQ) {
+            lp.colLowerBound(c, v);
+            lp.colUpperBound(c, v);
+          } else if (op1 == LE) {
+            lp.colLowerBound(c, v);
+          } else {
+            lp.colUpperBound(c, v);
+          }
+        }
+      }
+      return is;
     }
 
@@ -278,38 +394,7 @@
     static bool varChar(char c) {
       return (c >= '0' && c <= '9') || 
-        (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z');
-    }
-
-
-    void addConstraint(const LpSolverBase::Constr& constr) {
-      if (constr.expr().size() != 1) {
-        lp.addRow(constr);
-      } else {
-        Lp::Expr e = constr.expr();
-        LpSolverBase::Col col = e.begin()->first;
-        double coeff = e.begin()->second;
-        double lb = LpSolverBase::NaN;
-        double ub = LpSolverBase::NaN;
-        if (coeff > 0) {
-          if (constr.upperBounded()) {
-            lb = (constr.lowerBound() - e.constComp()) / coeff;
-          }
-          if (constr.lowerBounded()) {
-            ub = (constr.upperBound() - e.constComp()) / coeff;
-          }
-        } else if (coeff < 0) {
-          if (constr.upperBounded()) {
-            lb = (constr.upperBound() - e.constComp()) / coeff;
-          }
-          if (constr.lowerBounded()) {
-            ub = (constr.lowerBound() - e.constComp()) / coeff;
-          }
-        } else {
-          lb = -LpSolverBase::INF;
-          ub = LpSolverBase::INF;
-        }
-        lp.colLowerBound(col, lb);
-        lp.colUpperBound(col, ub);
-      }
+        (c >= 'a' && c <= 'z') || 
+        (c >= 'A' && c <= 'Z') ||
+        c == '[' || c == ']';
     }
     
@@ -321,24 +406,28 @@
     virtual void read(std::istream& is) {
       std::string line;
+      SubSection sub = CONSTRAINTS;
       while (getline(is, line)) {
         std::istringstream ls(line);
-        std::string sense;
-        ls >> sense;
-        if (sense == "min" || sense == "max") {
-          LpSolverBase::Expr expr;
-          ls >> std::ws;
-          readExpression(ls, expr);
-          lp.setObj(expr);
-          if (sense == "min") {
-            lp.min();
-          } else {
-            lp.max();
-          }
+        std::string type;
+        ls >> type;
+        if (type == "constraints") {
+          sub = CONSTRAINTS;
+        } else if (type == "bounds") {
+          sub = BOUNDS;
+        } else if (type == "objective") {
+          sub = OBJECTIVE;
         } else {
           ls.str(line);
-          LpSolverBase::Constr constr;
-          ls >> std::ws;
-          readConstraint(ls, constr);
-          addConstraint(constr);
+          switch (sub) {
+          case CONSTRAINTS:
+            readConstraint(ls);
+            break;
+          case BOUNDS:
+            readBounds(ls);
+            break;
+          case OBJECTIVE:
+            readObjective(ls);
+            break;
+          }
         }        
       }
@@ -361,332 +450,175 @@
 
 
-//   class LpWriter : public LemonWriter::SectionWriter {
-//     typedef LemonWriter::SectionWriter Parent;
-//   public:
-
-
-//     /// \brief Constructor.
-//     ///
-//     /// Constructor for LpWriter. It creates the LpWriter and attach
-//     /// it into the given LemonWriter. The lp writer will add
-//     /// variables, constraints and objective function to the
-//     /// given lp solver.
-//     LpWriter(LemonWriter& _writer, LpSolverBase& _lp, 
-//              const std::string& _name = std::string())
-//       : Parent(_writer), lp(_lp), name(_name) {} 
-
-
-//     /// \brief Destructor.
-//     ///
-//     /// Destructor for NodeSetWriter.
-//     virtual ~LpWriter() {}
-
-//   private:
-//     LpWriter(const LpWriter&);
-//     void operator=(const LpWriter&);
+  class LpWriter : public LemonWriter::SectionWriter {
+    typedef LemonWriter::SectionWriter Parent;
+  public:
+
+
+    /// \brief Constructor.
+    ///
+    /// Constructor for LpWriter. It creates the LpWriter and attach
+    /// it into the given LemonWriter. The lp writer will add
+    /// variables, constraints and objective function to the
+    /// given lp solver.
+    LpWriter(LemonWriter& _writer, LpSolverBase& _lp, 
+             const std::string& _name = std::string())
+      : Parent(_writer), lp(_lp), name(_name) {} 
+
+
+    /// \brief Destructor.
+    ///
+    /// Destructor for NodeSetWriter.
+    virtual ~LpWriter() {}
+
+  private:
+    LpWriter(const LpWriter&);
+    void operator=(const LpWriter&);
   
-//   protected:
-
-//     /// \brief Gives back true when the SectionWriter can process 
-//     /// the section with the given header line.
-//     ///
-//     /// It gives back the header line of the \c \@lp section.
-//     virtual std::string header() {
-//       std::ostringstream ls(line);
-//       ls << "@lp " << name;
-//       return ls.str();
-//     }
-
-//   private:
-
-//     std::ostream& writeConstraint(std::ostream& is, LpSolverBase::Constr& c) {
-//       char x;
-
+  protected:
+
+    /// \brief Gives back true when the SectionWriter can process 
+    /// the section with the given header line.
+    ///
+    /// It gives back the header line of the \c \@lp section.
+    virtual std::string header() {
+      std::ostringstream ls;
+      ls << "@lp " << name;
+      return ls.str();
+    }
+
+  private:
+
+    void createCols() {
+      int num = 1;
+
+      std::map<std::string, LpSolverBase::Col> inverse;
+
+      for (LpSolverBase::ColIt it(lp); it != INVALID; ++it) {
+        std::string name = lp.colName(it);
+        if (!name.empty() && inverse.find(name) == inverse.end()) {
+          cols.insert(std::make_pair(it, name));
+          inverse.insert(std::make_pair(name, it));
+        } else {
+          std::ostringstream ls;
+          ls << "var" << num;
+          ++num;
+          cols.insert(std::make_pair(it, ls.str()));
+          inverse.insert(std::make_pair(ls.str(), it));
+        }
+      }
+    }
+    
+    void writeExpression(std::ostream& os, const LpSolverBase::Expr& e) {
+      bool first = true;
+      for (LpSolverBase::Expr::const_iterator jt = e.begin(); 
+           jt != e.end(); ++jt) {
+        if (jt->second < 0.0) {
+          os << "- ";
+          if (jt->second != -1.0) {
+            os << -jt->second << " * ";
+          }
+        } else if (jt->second > 0.0) {
+          if (!first) {
+            os << "+ ";
+          }
+          if (jt->second != 1.0) {
+            os << jt->second << " * ";
+          }          
+        }
+        first = false;
+        os << cols[jt->first] << " ";
+      }
+      if (e.constComp() < 0.0) {
+        os << "- " << -e.constComp() << " ";
+      } else if (e.constComp() > 0.0) {
+        if (!first) {
+          os << "+ ";
+        }
+        os << e.constComp() << " ";
+      } 
+    }
+
+  protected:
+
+    /// \brief Writer function of the section.
+    ///
+    /// It writes the content of the section.
+    virtual void write(std::ostream& os) {
+      createCols();
+
+      os << "constraints" << std::endl;
       
-//       LpSolverBase::Expr e1, e2;
-//       Relation op1;
-//       is >> std::ws;
-//       writexpression(is, e1);
-//       is >> std::ws;
-//       writeRelation(is, op1);
-//       is >> std::ws;
-//       writexpression(is, e2);
-//       if (!is.get(x)) {
-//         if (op1 == LE) {
-//           c = e1 <= e2;
-//         } else if (op1 == GE) {
-//           c = e1 >= e2;
-//         } else {
-//           c = e1 == e2;
-//         }
-//       } else {
-//         is.putback(x);
-//         LpSolverBase::Expr e3;
-//         Relation op2;
-//         writeRelation(is, op2);
-//         is >> std::ws;
-//         writexpression(is, e3);
-//         if (!e1.empty() || !e3.empty()) {
-//           throw DataFormatError("Wrong range format");
-//         }
-//         if (op2 != op1 || op1 == EQ) {
-//           throw DataFormatError("Wrong range format");
-//         }
-//         if (op1 == LE) {
-//           c = e1.constComp() <= e2 <= e3.constComp();
-//         } else {
-//           c = e1.constComp() >= e2 >= e3.constComp();
-//         }
-//       }
-//     }
-
-//     std::ostream& writexpression(std::ostream& is, LpSolverBase::Expr& e) {
-//       LpSolverBase::Col c;
-//       double d;
-//       char x;
-//       writelement(is, c, d);
-//       if (c != INVALID) {
-//         e += d * c;
-//       } else {
-//         e += d;
-//       }
-//       is >> std::ws;
-//       while (is.get(x) && (x == '+' || x == '-')) {
-//         is >> std::ws;
-//         writelement(is, c, d);
-//         if (c != INVALID) {
-//           e += (x == '+' ? d : -d) * c;
-//         } else {
-//           e += (x == '+' ? d : -d);
-//         }
-//         is >> std::ws;
-//       }
-//       if (!is) {
-//         is.clear();
-//       } else {
-//         is.putback(x);
-//       }
-//       return is;
-//     }
-
-//     std::ostream& writelement(std::ostream& is, 
-//                               LpSolverBase::Col& c, double& d) { 
-//       d = 1.0;
-//       c = INVALID;
-//       char x, y;
-//       if (!is.get(x)) throw DataFormatError("Cannot find lp element");
-//       if (x == '+' || x == '-') {
-//         is >> std::ws;
-//         d *= x == '-' ? -1 : 1;
-//         while (is.get(x) && (x == '+' || x == '-')) {
-//           d *= x == '-' ? -1 : 1;
-//           is >> std::ws;
-//         }
-//         if (!is) throw DataFormatError("Cannot find lp element");
-//       }
-//       if (numFirstChar(x)) {
-//         is.putback(x);
-//         double e;
-//         writeNum(is, e);
-//         d *= e;
-//       } else if (varFirstChar(x)) {
-//         is.putback(x);
-//         LpSolverBase::Col f;
-//         writeCol(is, f);
-//         c = f;
-//       } else {
-//         throw DataFormatError("Invalid expression format");          
-//       }
-//       is >> std::ws;
-//       while (is.get(y) && (y == '*' || y == '/')) {
-//         is >> std::ws;
-//         if (!is.get(x)) throw DataFormatError("Cannot find lp element");
-//         if (x == '+' || x == '-') {
-//           is >> std::ws;
-//           d *= x == '-' ? -1 : 1;
-//           while (is.get(x) && (x == '+' || x == '-')) {
-//             d *= x == '-' ? -1 : 1;
-//             is >> std::ws;
-//           }
-//           if (!is) throw DataFormatError("Cannot find lp element");
-//         }
-//         if (numFirstChar(x)) {
-//           is.putback(x);
-//           double e;
-//           writeNum(is, e);
-//           if (y == '*') {
-//             d *= e;
-//           } else {
-//             d /= e;
-//           }
-//         } else if (varFirstChar(x)) {
-//           is.putback(x);
-//           LpSolverBase::Col f;
-//           writeCol(is, f);
-//           if (y == '*') {
-//             if (c == INVALID) {
-//               c = f;
-//             } else {
-//               throw DataFormatError("Quadratic element in expression");
-//             }
-//           } else {
-//             throw DataFormatError("Division by variable");
-//           }
-//         } else {
-//           throw DataFormatError("Invalid expression format");          
-//         }
-//         is >> std::ws;
-//       }
-//       if (!is) {
-//         is.clear();
-//       } else {
-//         is.putback(y);
-//       }
-//       return is;
-//     }
-
-//     std::ostream& writeCol(std::ostream& is, LpSolverBase::Col& c) {
-//       char x;
-//       std::string var;
-//       while (is.get(x) && varChar(x)) {
-//         var += x;
-//       }
-//       if (!is) {
-//         is.clear();
-//       } else {
-//         is.putback(x);
-//       }
-//       ColMap::const_iterator it = cols.find(var);
-//       if (cols.find(var) != cols.end()) {
-//         c = it->second;
-//       } else {
-//         c = lp.addCol();
-//         cols.insert(std::make_pair(var, c));
-//         lp.colName(c, var);
-//       }
-//       return is;
-//     }
-
-//     std::ostream& writeNum(std::ostream& is, double& d) {
-//       is >> d;
-//       if (!is) throw DataFormatError("Wrong number format");
-//       return is;
-//     }
-
-//     std::ostream& writeRelation(std::ostream& is, Relation& op) {
-//       char x, y;
-//       if (!is.get(x) || !(x == '<' || x == '=' || x == '>')) {
-//         throw DataFormatError("Wrong relation operator");
-//       }
-//       if (!is.get(y) || y != '=') {
-//         throw DataFormatError("Wrong relation operator");
-//       }
-//       switch (x) {
-//       case '<': op = LE; 
-//         break;
-//       case '=': op = EQ; 
-//         break;
-//       case '>': op = GE; 
-//         break;
-//       }
-//       return is;
-//     }
-
-//     static bool relationFirstChar(char c) {
-//       return c == '<' || c == '=' || c == '>';
-//     }
-
-//     static bool varFirstChar(char c) {
-//       return (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z');
-//     }
-
-//     static bool numFirstChar(char c) {
-//       return (c >= '0' && c <= '9') || c == '.';
-//     }
-
-//     static bool varChar(char c) {
-//       return (c >= '0' && c <= '9') || 
-//         (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z');
-//     }
-
-
-//     void addConstraint(const LpSolverBase::Constr& constr) {
-//       if (constr.expr().size() != 1) {
-//         lp.addRow(constr);
-//       } else {
-//         Lp::Expr e = constr.expr();
-//         LpSolverBase::Col col = e.begin()->first;
-//         double coeff = e.begin()->second;
-//         double lb = LpSolverBase::NaN;
-//         double ub = LpSolverBase::NaN;
-//         if (coeff > 0) {
-//           if (constr.upperBounded()) {
-//             lb = (constr.lowerBound() - e.constComp()) / coeff;
-//           }
-//           if (constr.lowerBounded()) {
-//             ub = (constr.upperBound() - e.constComp()) / coeff;
-//           }
-//         } else if (coeff < 0) {
-//           if (constr.upperBounded()) {
-//             lb = (constr.upperBound() - e.constComp()) / coeff;
-//           }
-//           if (constr.lowerBounded()) {
-//             ub = (constr.lowerBound() - e.constComp()) / coeff;
-//           }
-//         } else {
-//           lb = -LpSolverBase::INF;
-//           ub = LpSolverBase::INF;
-//         }
-//         lp.colLowerBound(col, lb);
-//         lp.colUpperBound(col, ub);
-//       }
-//     }
-    
-//   protected:
-
-//     /// \brief Writer function of the section.
-//     ///
-//     /// It writes the content of the section.
-//     virtual void write(std::ostream& is) {
-//       std::string line;
-//       std::map<std::string, LpSolverBase::Col> vars;
-//       while (getline(is, line)) {
-//         std::istringstream ls(line);
-//         std::string sense;
-//         ls >> sense;
-//         if (sense == "min" || sense == "max") {
-//           LpSolverBase::Expr expr;
-//           ls >> std::ws;
-//           writeExpression(ls, expr);
-//           lp.setObj(expr);
-//           if (sense == "min") {
-//             lp.min();
-//           } else {
-//             lp.max();
-//           }
-//         } else {
-//           ls.str(line);
-//           LpSolverBase::Constr constr;
-//           ls >> std::ws;
-//           writeConstraint(ls, constr);
-//           addConstraint(constr);
-//         }        
-//       }
-//     }
+      for (LpSolverBase::RowIt it(lp); it != INVALID; ++it) {
+        double lower, upper;
+        lp.getRowBounds(it, lower, upper);
+        if (lower == -LpSolverBase::INF && upper == LpSolverBase::INF) 
+          continue;
+        os << "   ";
+        
+        if (lower != upper) {
+          if (lower != -LpSolverBase::INF) {
+            os << lower << " <= ";
+          }
+          
+          writeExpression(os, lp.row(it));
+          
+          if (upper != LpSolverBase::INF) {
+            os << "<= " << upper << " ";
+          }
+        } else {
+
+          writeExpression(os, lp.row(it));
+          os << "== " << upper << " ";
+          
+        }
+
+        os << std::endl;
+      }
+
+      os << "bounds" << std::endl;
       
-//     virtual void missing() {
-//       ErrorMessage msg;
-//       msg << "Lp section not found in file: @lp " << name;
-//       throw IoParameterError(msg.message());
-//     }
-
-//   private:
-
-//     typedef std::map<std::string, LpSolverBase::Col> ColMap;
+      for (LpSolverBase::ColIt it(lp); it != INVALID; ++it) {
+        double lower = lp.colLowerBound(it);
+        double upper = lp.colUpperBound(it);
+        if (lower == -LpSolverBase::INF && upper == LpSolverBase::INF) 
+          continue;
+        os << "   ";
+
+        if (lower != upper) {
+          if (lower != -LpSolverBase::INF) {
+            os << lower << " <= ";
+          }
+          
+          os << cols[it] << " ";
+          
+          if (upper != LpSolverBase::INF) {
+            os << "<= " << upper << " ";
+          }
+        } else {
+          os << cols[it] << " == " << upper << " ";
+        }
+        os << std::endl;
+      }
+
+      os << "objective" << std::endl;
+      os << "   ";
+      if (lp.is_min()) {
+        os << "min ";
+      } else {
+        os << "max ";
+      }
+      writeExpression(os, lp.obj());
+      os << std::endl;
+    }
       
-//     LpSolverBase& lp;
-//     std::string name;
-//     ColMap cols;
-//   };
+
+  private:
+
+    typedef std::map<LpSolverBase::Col, std::string> ColMap;
+      
+    LpSolverBase& lp;
+    std::string name;
+    ColMap cols;
+  };
 
 }
